CHEMICAL ENGINEERING DEPARTMENT
INDIAN INSTITUTE OF TECHNOLOGY MADRAS
CFDT PROJECT Steady Flow of an Incompressible Fluid in a Lid Driven Cavity Jason Ryan Picardo CH11D026 4/5/2012
Course: Computational Fluid Dynamics Tech Course No: CH 6020
Table of Contents Declaration 1
Introduction ........................................................................................................................ 4
2
The SIMPLE Method ......................................................................................................... 5
3
Discretization of the Governing Equations ........................................................................ 6 3.1
Momentum Equations ................................................................................................. 7
3.2
Pressure Correction Equation ...................................................................................... 9
3.3
Boundary Conditions................................................................................................. 11
3.4
Under-Relaxation of Pressure and Velocity Fields ................................................... 12
4
Results and Discussion ..................................................................................................... 13
5
Comparison of Results with Published Work................................................................... 16
6
Conclusions ...................................................................................................................... 17 Appendix .......................................................................................................................... 18 References ........................................................................................................................ 24
2
Declaration
The work presented in this report is the original work done by Jason Ryan Picardo, Roll No. CH11D026 of the Chemical Engineering Department, IIT Madras towards partial requirement for the completion of the course CH6020- Computational Fluid Dynamics Tech in the Jan-May semester of 2012.
Jason Ryan Picardo
3
1
Introduction
In this report a Computational Fluid Dynamics (CFD) solution to the 'Lid Driven Cavity' problem is presented. The steady flow of an incompressible fluid is considered. The SIMPLE method (Semi-Implicit Method for Pressure-Linked Equations), originally devised by Patankar and Spalding,1 to deal with incompressible steady flows, is adopted. The 'Lid Driven Cavity' problem is a standard problem that is often used to test new CFD codes and solution methods. The problem domain has a simple geometry with easily implementable boundary conditions. These characteristics have made it a popular benchmark problem. A standard two-dimensional case consists of a rectangular cavity with its top lid moving at a constant velocity (Vw). A schematic of the problem is depicted in Figure 1. The cavity is assumed to be infinite in the z-direction. The motion of the top lid creates a 2-D circulatory flow pattern inside the cavity. While the flow remains steady and laminar at moderate Reynolds numbers (Re), it does exhibit vortices, which grow with Re, at the lower sharp corners.
Figure 1. Schematic of the Lid Driven Cavity problem. Vw is the constant velocity at which the wall moves. The cavity has a height H in the y direction and a width W in the x-direction. The special case of H = W corresponds to a square cavity. The SIMPLE method is briefly described in Section 2. The discretized equations and the pressure correction equation is presented in Section 3. Results for various Re are presented and discussed in Section 4. In order to validate the present calculations, some results are compared with the work of Ghia and co-workers 2 in Section 0. Several obstacles were faced during the development of the code and the implementation of the SIMPLE method. This lead to new learning and insights which are discussed in the Conclusions.
4
2
The SIMPLE Method
The SIMPLE method, first formulated by Patankar and Spalding, 1 has proved to be an effective method for the solution of CFD problems involving incompressible fluids, especially when the steady solution is of interest. The method is described in detail in a text by Patankar. 3 It is essentially a pressure correction type method. For any given 'guessed' pressure field, the velocity fields are readily computed using the Navier-Stokes (N-S) Momentum Equations (density is a constant for incompressible fluids). However, the computed velocity fields will not, in general, satisfy continuity unless the guessed pressure field is 'correct'. In a pressure correction method, a pressure correction equation is developed which applies a correction to the initial guessed pressure such that the corresponding velocity fields satisfy the Continuity Equation. In the SIMPLE method, the pressure correction equation is only approximate, as shall be seen in Section 3.2. Moreover, the N-S Momentum equations are nonlinear and highly coupled. Thus an iterative method is necessary to arrive at the final solution. A flow chart of the major steps in the SIMPLE method is presented in Figure 2.
Figure 2. Flow Chart of the SIMPLE Method Before the SIMPLE method can be implemented, the N-S Momentum Equations must be suitably discretized. The exact form of the pressure correction equation will depend on the manner in which the Momentum equations are discretized. Suitable boundary conditions must also be specified with a means to implement them in the CFD code. This is dealt with in the next section.
5
3
Discretization of the Governing Equations
In this work a staggered grid is used, following Patankar. 3 This prevents the occurrence of checkerboard oscillations. The pressure is calculated at the (i , j) nodes while the x and y velocity components are calculated at (i+1/2 , j) and (i, j+1/2) respectively. Thus the grid for the x- velocity is staggered to the right (in the positive x direction) and the grid for the yvelocity is staggered to the top ( in the positive y direction) with respect to the pressure grid. A typical staggered grid is depicted in Figure 3.
Figure 3. A typical staggered grid. The pressure is calculated at the (i , j) nodes while the x and y velocity components are calculated at (i+1/2 , j) and (i, j+1/2) respectively.
The non-dimensional continuity equation for an incompressible fluid and the N-S Momentum Equations in conservative form are presented below.
u W v 0 x H y 2 u uu W vu P 1 2u W 2u t x x Re x 2 H y 2 H y
2 v uv W vv P 1 2 v W 2 v t x y Re x 2 H y 2 H y
(1)
(2)
(3)
6
Here, the following characteristic scales are used to non-dimensionalize the equations. u c vc Vw Pc u c2
xc W
Re
yc H
tc
W uc
(4)
u cW
The SIMPLE method elucidated by Patankar 3 is based on a Finite Volume Method (FVM) discretization of the equations. Here we adopt a Finite Difference (FDM) discretization scheme involving a fully implicit forward difference in time and central differences in space for both the convective and diffusive terms. A similar scheme is presented in the discussion of the SIMPLE method by Anderson. 4 The use of the conservative form of the N-S equations leads to a final set of discretized equations which are similar to those obtained using the FVM.4 In this problem we are interested in obtaining the steady flow solution. In spite of this we use the unsteady N-S Momentum Equations. It is much easier to converge to the solution from an arbitrary initial guess of the pressure and velocity fields if the unsteady equations are used. It must be noted that the SIMPLE method does not result in time accurate solutions. Here, the discrete time steps should be interpreted as sequential iteration steps which ultimately converge to the steady flow solution. The time step value (t) is simply a parameter which affects the speed of convergence and the stability of the iterations. High values can cause the iterations to become unstable. Anderson uses a fully explicit time differencing scheme which makes the calculations of the velocity fields simple. However, this makes the method susceptible to instability at relatively large values of the time step. In this work, a fully implicit scheme is adopted. While this requires the simultaneous solution of the equations at all grid points, it provides greater stability to the method.
3.1
Momentum Equations
Based on the previous discussion, a finite difference (FDM) based discretization of the XMomentum equation yields uin112, j ui,n1 2, j t
,n
n 1 W vn u H
u
,n
u n1
i 1 2 , j 1
i 3 2 , j
i 1 2 , j
2x ,n
vs u n1
2y
u ,n u n 1
i 1 2 , j 1
Pi 1,n, j Pi , j,n
(5)
x
2 n 1 n 1 n 1 n 1 n 1 n 1 1 ui 3 2, j 2ui 1 2, j ui 1 2, j W ui 1 2, j 1 2ui 1 2, j ui 1 2, j 1 Re x 2 y 2 H
7
After collecting terms Eq. (5) results in a linear equation for the velocity at (i+1/2 , j) in terms of the neighbouring nodes. u
n 1 i 1 2 , j
2 1 ui,n3 2, j 2 2 1 W n 1 u i 3 2 , j 2 2 Re x 2 H Re y 2x t Re x 2 vn,n W ui,1n 2, j 1 1 W n 1 ui 1 2, j 1 2 Re x 2 2x 2y H Re y H
u
n 1 i 1 2 , j
u
n 1 i 1 2 , j 1
2 ,n ,n vs ,n W ui,1n 2, j 1 W Pi 1, j Pi , j 2 x t 2y H Re y H
(6)
with ,n
vn
,n
vs
vi,1n, j 1 2 vi,,jn1 2 2 vi,,jn1 2
,n i 1, j 1 2
v
2
Similarly, the Y- Momentum equation may be discretized to give the following equation for the node (i, j+1/2)
n 1 i , j 1 2
v
2 2 1 vi,,jn3 2 W 2 2 1 W W n 1 vi , j 3 2 2 2 2 t Re x H Re y 2y H Re y H
n 1 i , j 1 2
v
,n 2 vi,,jn1 2 W 1 W n1 uf 1 v i 1, j 1 2 2 2 2x Re x 2y H Re y H
,n ,n ub,n vi,,jn1 2 1 Pi , j 1 Pi , j vin11, j 1 2 2 y t 2x Re x with
uf ub
,n
,n
(7)
ui,1n 2, j 1 ui,1n 2, j u
,n i 1 2, j 1
2 ui,1n 2, j 2
Clearly these equations must be solved simultaneously for all nodes. This may be done using a Gauss-Siedel iterative solution method with under-relaxation. In the calculations an under-relaxation factor of 0.7 is used. Thus given an initial pressure and velocity field (denoted by superscripted asterisk) new velocity fields may be calculated which satisfy the N-S momentum equations. This takes care of the first step in the iterative loop in Figure 2.
8
However, the computed velocity field will not, in general, satisfy the continuity equation since the pressure field was an initial guess value. The pressure correction equation is necessary to rectify this by incorporating suitable corrections to the pressure and velocity fields. 3.2
Pressure Correction Equation
In this section the discretized pressure correction equation is developed which enforces conservation of mass on the calculated velocity fields. The Equation of Continuity (Mass Conservation Equation) for an incompressible fluid is given in Eq. (1). The discretized equation for the node (i , j) is presented below. Central differences are used for both derivatives.
u~i1 2, j u~i1 2, j x
~ ~ W vi , j 1 2 vi , j 1 2 0 y H
(8)
u~ and v~ are velocity fields which satisfy continuity. The velocity fields computed from the momentum equations, in the previous step, will not satisfy the continuity equation. A velocity correction ( u and v ) is defined such that the corrected velocity field satisfies the continuity Eq. (8).
u~ u u v~ v v
(9)
The velocity correction ( u and v ) may be calculated from a corresponding correction to the pressure field ( P ). ~ P P P
(10)
If these two corrections are made to satisfy the N-S Momentum Equations (Eqs. (6(7) then ~ the continuity satisfying velocity field ( u~ and v~ ) and the corrected pressure field ( P ) will satisfy the governing equations of fluid motion for incompressible flow. This would then be the solution to the problem. However, to simplify calculations, only the direct effect of the pressure correction on the velocity at any node is accounted for i.e. the terms involving velocities at neighbouring nodes are dropped.
9
ui1 2, j vi, j 1 2
P
Pi, j
Ax Pi, j 1 Pi, j
i 1, j
Ay
(11)
with 2 1 2 2 W A 2 2 t Re x H Re y
This simplification will lead to a velocity and pressure field which satisfy mass conservation but not the N-S momentum equations. This is step two in the SIMPLE procedure (Figure 2). These fields are then taken as the new initial guesses and the process is repeated. Finally, the iterations will converge to a solution that satisfies the governing equations. The simplification in the pressure correction equation does not affect the final steady flow solution. The velocity correction is substituted into Eq. (8) along with the velocity field calculated in Step 1 of Figure 2 from Eqs. (6)) and ((7) to yield a discretized equation for the pressure correction.
2 2 Pi, j 2 Ay 2 Ax
W H
2
1 1 Pi1, j Pi1, j 2 2 Ax Ax
1 W 2 1 W 2 q( i , j ) Pi , j 1 P Ay 2 H i , j 1 Ay 2 H with q( i , j )
ui 1 2, j ui 1 2, j x
(12)
W vi , j 1 2 vi , j 1 2 y H
2 1 2 2 W A 2 2 t Re x H Re y
This is a pressure correction formula which has the form of a Poisson equation for the corrected pressure. q(i , j) is a source term which is in fact the residuals from the continuity equation of the incorrect velocity field at every node. When q(i , j) is identically zero, it implies that the field computed at Step 1 of Figure 2 satisfies continuity and is in fact the solution. This is reflected by the fact that the pressure corrections calculated from Eq. (12) will be zero. Hence the values of q(i , j) serve as a useful criteria to check for convergence. In this work the root mean square value ( ) of q(i , j) across all nodes is used as the convergence criterion (Eq. (13)).
10
q( i , j )
2
i,j
(13)
with q( i , j )
ui 1 2, j ui 1 2, j x
W vi , j 1 2 vi , j 1 2 y H
As discussed in Section 3, the solution at various 'time steps' does not correspond to any physically relevant transient solution due to the approximation in the pressure correction equation. They should be interpreted solely as iteration steps which converge to the steady solution. 4 Nevertheless, if the velocity correction is not relaxed, then the SIMPLE procedure converges to the final solution via a series of intermediate velocity fields which satisfy continuity. As Patankar 3 points out, this is advantageous in terms of numerical stability and for the calculation of scalar transport variables (heat, concentration etc.) as they can be calculated based on velocity fields which satisfy mass conservation.
3.3
Boundary Conditions
Boundary conditions must be specified for the momentum equations (Eqs. (6) and (7)) and the pressure correction equation (Eq. (12)). The velocity components are zero everywhere along the stationary walls from the no-slip condition and the kinematic condition. Along the top moving wall the y-velocity is zero while the x-velocity is taken to be that of the wall (1 in dimensionless variable form- see Eq.(4) and Figure 1). A homogeneous Neumann condition is applied for the derivative of the pressure correction field. This is justifiable since the normal components of velocity are zero at all the walls according to the kinematic condition.3,4
P 0 x P 0 y
at
x 0,1 (14)
at
y 0,1
11
3.4 Under-Relaxation of Pressure and Velocity Fields The iterations may become divergent for many problems. To ensure convergence an underrelaxation factor may be used while updating the pressure and velocity fields.
P ,n1 P ,n P P u ,n1 u n1 u u
(15)
v ,n1 v n1 v v The value must be between zero and one for under-relaxation. A value of zero corresponds to no correction. Patankar recommends values of P = 0.8 and u = v = 0.5. Anderson uses a more conservative estimate of P = 0.1 for solution of the Couette flow problem. In this work, a value of P = 0.1 , 0.05 is used depending on Re and the number of grid points. Smaller values of relaxation are required at higher Re and a larger number of grid points. u and v are both set to one. Thus the intermediate velocity field at various iterations satisfies mass conservation. The t value also acts as an additional relaxation factor and can be adjusted to affect convergence rate and stability. For the Couette problem Anderson uses a small value of t = 0.001 since an explicit time differencing scheme is adopted in that work. Here, a larger value of t = 0.01 is used as a fully implicit scheme of time discretization has been adopted (Section 3.1).
12
4
Results and Discussion
A code is developed in MATLAB based on the previously outlined solution method. A sample code for a particular case of Re = 100 and a 80 X 80 grid is given in the Appendix. Two simulations were carried out for the case of a square cavity. The first at a Re = 100 and the second at Re = 400. An 80 x 80 grid was used for the first case and a 81 X 81 grid was used for the second. These computations took about 16 and 24 hrs respectively. The convergence history of the computations are presented in Figure 4. Here, is the convergence criterion defined in Eq. (13).
Figure 4. Convergence history for the calculations for the cases of (a) Re = 100 on 80 X 80 grid and (b) Re = 400 on 81 X 81 grid. Streamline plots of the velocity field for both cases are presented in Figure 5. The circulatory flow pattern is clearly visible. The vortices formed at the lower edges of the cavity due to boundary layer separation grow as the Reynolds number is increased. This is due to the dominance of inertial forces over the viscous forces. In both cases, the top lid has a positive scaled velocity of 1 (in the positive x direction). In Figure 6, contours of the x and y velocity components are displayed for the case of Re = 400 on 81 X 81 grid. The velocity components switch direction across the cavity, as is consistent with a circulatory flow. Contours of pressure are displayed for the same case in Figure 7. The stagnation point is clearly visible at the top right corner of the cavity where the pressure is highest. The vortex regions at the two lower corners of the cavity are at a relatively higher pressure. This is characteristic of boundary layer separation situations.
13
Figure 5. Streamline plots for (a) Re = 100 on 81 X 81 grid and (b) Re = 400 on 81 X 81 grid.
Figure 6. Contours of the x and y components of velocity for the case of Re = 400 on 81 X 81 grid.
Figure 7. Contours of the pressure for the case of Re = 400 on 81 X 81 grid.
14
First, an interesting variant of the present problem is briefly considered. That is the case when both the top and bottom plates are moved at constant velocities, in the same direction. The resultant flow field will have some relation to the flow field inside slugs/drops moving within a second, immiscible, continuous fluid. Both fluids flow in a channel of diameter comparable to the slug/drop diameter. Such a scenario arises in two phase flow in microchannels and has several technological applications. In particular, mass transfer in liquid-liquid extraction in micro channels is enhanced due to the circulatory flow inside the slugs/drops. 5 Figure 8 displays the streamlines for a first computation at Re = 5 and a grid of 40 X 40. This low Re is typical for the case of a microchannel flow.
Figure 8. Streamlines for the two moving wall problem- reminiscent of circulation in drops/slugs. Re = 5 with a grid of 40 X 40. The flow field is symmetric as expected, with two counter rotating central vortices. No secondary vortices are formed since Re is very small. The velocity contours are presented in Figure 9 below. Here, both the plates are moved in the negative x direction uniformly.
Figure 9. Contours of the x and y components of velocity for the the two moving wall problem. Re = 5 with a grid of 40 X 40.
15
5
Comparison of Results with Published Work
In this section, the calculations done using the SIMPLE scheme are compared with the results of Ghia and co-workers. 2 The authors use a vorticity-stream function formulation of the 2D N-S Momentum Equations with a coupled strongly implicit multigrid method. Tabular results of velocity are presented which allows for convenient comparison. Figure 10 displays the comparison of results for Re = 100 and Re = 400. Ghia and co-workers 2 use a comparatively fine mesh of 129 X 129 grids. The results of the present work compare well for the case of Re = 100 , considering the coarse 81 X 81 grid used. However, at the higher Re = 400 a discrepancy is observed in the results. Clearly, a finer mesh must be used. This will require a longer computation time.
Figure 10. Comparison of results of the present work with the results of Ghia and coworkers.2 Comparisons are depicted for Re = 100 on 80 X 80 grid and Re = 400 on 81 X 81 grid with the corresponding results of Ghia et. al. on a 129 X 129 grid.
16
6
Conclusions
The SIMPLE method is successfully applied to solve the 'Lid Driven Cavity' problem. A finite difference scheme with a fully implicit time discretization is used. The method is implemented as a code in MATLAB. The calculations were found to converge monotonically. The results calculated for Re = 100, 400 on relatively coarse grids are compared with the published work of Ghia and co-workers 2. The comparison reveals a satisfactory agreement for lower Reynolds numbers (Re = 100). While a finer mesh must be used for better results, particularly at high Re, the SIMPLE approach is found to be adequate to solve this steady, incompressible fluid flow problem. An important point is the use of the unsteady N-S Momentum Equations as opposed to the steady state equations. In prior attempts at developing a SIMPLE based CFD solution, the steady state equations were also used without much success by the author. It is found that the SIMPLE method with the steady N-S Momentum equations does not converge for problems with Re > 50 and/or with a grid larger than 20 X 20. Even when the iterations converge, a very small relaxation factor of P = 0.001. This is to be contrasted with the value of 0.1/0.05 which is successfully used in the calculations with the unsteady equations. It seems that the SIMPLE method based on the unsteady equations is much more suited to the calculation of the final steady flow field via a series of iterations beginning with an arbitrary initial guess. The values of the variables at successive iterations do not correspond to a transient solution since the SIMPLE pressure correction equation is only approximate. Hence the time step, t is a factor similar to a relaxation factor which controls the rate of convergence and the stability of the iterations. It must be noted that large values of t result in a set of discretized equations which tend to the steady state discretized equations. Thus, based on the previous discussion, small values of t should be used to ensure convergence. An interesting variant of the problem occurs when both the top and bottom plates are moved at constant velocities, in the same direction. The resultant flow field is reminiscent of the flow field inside slugs/drops moving within a second, immiscible, continuous fluid. This flow field arises in microchannel flows and has much technological relevance. Since the computation time is quite large, this report presents only the first results of the application of the SIMPLE method to the Lid Driven Cavity problem. However, there are several aspects which entail further calculations. Firstly, simulations should be run for higher grid sizes in an attempt to establish grid independency. These results should then be compared with published work.
17
7
Appendix A sample MATLAB code for a particular case of Re = 100 and a 80 X 80 grid.
% This code solves the moving lid problem. The steady state solution is % computed for the case of incompressible viscous flow via the SIMPLE % algorithm of Patankar and Spalding. % Author Details % Jason Ryan Picardo- CH11D026 % Chemical Engg Dept- IIT Madras %TIME IMPLICIT clear all %Problem Parameters Uw=1; %velocity of moving lid (m/s) W=1; %Width of cavity (m) H=1; %Height of cavity (m) Re=100; %Reynolds No %Convergence array convergence=zeros(1,2000); k=0; %Grid Properties np=80; nu=np-1; nv=np-1; % x is alligned along the width (index i) origin at bottom left dx=1/(np-1); % y is alligned along the height (index i) origin at bottom left dy=1/(np-1); %artificial time dt=0.01; %Initialization Pstr=zeros(np,np); Ustr=zeros(np,nu); Vstr=zeros(nv,np); %A(j,i) Ustr(np,:)=Uw; qsum=1; %----------------------------------------------------------------------------------%SIMPLE LOOP while sqrt(qsum)>10^-3 k=k+1; U=zeros(np,nu); U(np,:)=Uw; V=zeros(nv,np); relax=0.7; %Calculation of U field from Pstr and Ustr and Vstr U(np,:)=Uw; % top boundary moving, all other boundaries are at 0 residual1=1; residual2=1;
18
while residual1>10^(-8) if(residual1>residual2&&residual1>100) display('U cal failed') display(residual1) break; end residual2=residual1; for i=2:nu-1 for j=2:np-1 Vnstr=(Vstr(j,i+1)+Vstr(j,i))/2; Vsstr=(Vstr(j-1,i+1)+Vstr(j-1,i))/2; U(j,i)=(1relax)*U(j,i)+relax*((1/(1/dt+2*(1/Re)/dx^2+2*(1/Re)/(dy*H/W)^2))*(Ustr(j,i )/dt+U(j,i+1)*(-Ustr(j,i+1)/(2*dx)+(1/Re)/dx^2)+... U(j,i-1)*(Ustr(j,i-1)/(2*dx)+(1/Re)/dx^2)+U(j+1,i)*(Vnstr/(2*(dy*H/W))+(1/Re)/(dy*H/W)^2)+U(j1,i)*(Vsstr/(2*(dy*H/W))+(1/Re)/(dy*H/W)^2)... -1/(dx)*(Pstr(j,i+1)-Pstr(j,i)))); end end delF=ones(1,(nu-2)*(np-2)); c=1; for i=2:nu-1 for j=2:np-1 Vnstr=(Vstr(j,i+1)+Vstr(j,i))/2; Vsstr=(Vstr(j-1,i+1)+Vstr(j-1,i))/2; delF(c)= U(j,i)(1/(1/dt+2*(1/Re)/dx^2+2*(1/Re)/(dy*H/W)^2))*(Ustr(j,i)/dt+U(j,i+1)*(Ustr(j,i+1)/(2*dx)+(1/Re)/dx^2)+... U(j,i-1)*(Ustr(j,i-1)/(2*dx)+(1/Re)/dx^2)+U(j+1,i)*(Vnstr/(2*(dy*H/W))+(1/Re)/(dy*H/W)^2)+U(j1,i)*(Vsstr/(2*(dy*H/W))+(1/Re)/(dy*H/W)^2)... -1/(dx)*(Pstr(j,i+1)-Pstr(j,i))); c=c+1; end end sum=0; for c=1:length(delF) sum=sum+delF(c)^2; end residual1=sqrt(sum/length(delF)); end %display('U calculation') %display(residual); %Calculation of V field from Pstr and Ustr and Vstr residual1=1; residual2=1; relax=0.7; while residual1>10^(-8) if(residual1>residual2&&residual1>100)
19
display('V cal failed') display(residual1) break; end residual2=residual1; for i=2:np-1 for j=2:nv-1 Ufstr=(Ustr(j+1,i)+Ustr(j,i))/2; Ubstr=(Ustr(j+1,i-1)+Ustr(j,i-1))/2; V(j,i)=(1relax)*V(j,i)+relax*((1/(1/dt+2*(1/Re)/dx^2+2*(1/Re)/(dy*H/W)^2))*(Vstr(j,i )/dt+V(j,i+1)*(-Ufstr/(2*dx)+(1/Re)/dx^2)+... V(j,i-1)*(Ubstr/(2*dx)+(1/Re)/dx^2)+V(j+1,i)*(Vstr(j+1,i)/(2*(dy*H/W))+(1/Re)/(dy*H/W)^2)+V(j-1,i)*(Vstr(j1,i)/(2*(dy*H/W))+(1/Re)/(dy*H/W)^2)... -1/((dy*H/W))*(Pstr(j+1,i)-Pstr(j,i)))); end end delF=ones(1,(nv-2)*(np-2)); c=1; for i=2:np-1 for j=2:nv-1 Ufstr=(Ustr(j+1,i)+Ustr(j,i))/2; Ubstr=(Ustr(j+1,i-1)+Ustr(j,i-1))/2; delF(c)= V(j,i)(1/(1/dt+2*(1/Re)/dx^2+2*(1/Re)/(dy*H/W)^2))*(Vstr(j,i)/dt+V(j,i+1)*(Ufstr/(2*dx)+(1/Re)/dx^2)+... V(j,i-1)*(Ubstr/(2*dx)+(1/Re)/dx^2)+V(j+1,i)*(Vstr(j+1,i)/(2*(dy*H/W))+(1/Re)/(dy*H/W)^2)+V(j-1,i)*(Vstr(j1,i)/(2*(dy*H/W))+(1/Re)/(dy*H/W)^2)... -1/((dy*H/W))*(Pstr(j+1,i)-Pstr(j,i))); c=c+1; end end sum=0; for c=1:length(delF) sum=sum+delF(c)^2; end residual1=sqrt(sum/length(delF)); end %display('V calculation') %display(residual) %Updating velocity fields Ustr=U; Vstr=V; %calculation of error in continuity from Ustr and Vstr Q=zeros(np,np); qsum=0; for i=2:np-1 for j=2:np-1 Q(j,i)=(Ustr(j,i)-Ustr(j,i-1))/dx+(Vstr(j,i)-Vstr(j-1,i))/(dy*H/W); qsum=qsum+Q(j,i)^2;
20
end end display(k) sqrt(qsum) convergence(k)=sqrt(qsum); if(k>2&&convergence(k)>1.5*convergence(k-1)) display('Iterations are diverging') break; end if(sqrt(qsum)<10^-3) display('convergence achieved to 10^-3 in continuity') break end save('Re_V4_100_np_80_con v_o3','Ustr','Vstr','Pstr','convergence','k','Re','np','dx','dy') %Calculation of P correction from Pstr and Ustr and Vstr Pcor=zeros(np,np); relax=0.7; residual1=1; residual2=1; while residual1>10^(-8) if(residual1>residual2&&residual1>100) display('Pcor cal failed') display(residual1) break; end residual2=residual1; A=-(1/(dx))*(1/(1/dt+2*(1/Re)/dx^2+2*(1/Re)/(dy*H/W)^2)); B=-(1/((dy*H/W)))*(1/(1/dt+2*(1/Re)/dx^2+2*(1/Re)/(dy*H/W)^2)); for i=2:np-1 for j=2:np-1 if(i==2&&j==2) Pcor(i,j)=(1relax)*Pcor(i,j)+(relax)*(1/(A/dx+B/(dy*H/W))*(A/dx*Pcor(i+1,j)+B/(dy*H/W)* Pcor(i,j+1)+Q(i,j))); elseif(i==2&&j==np-1) Pcor(i,j)=(1relax)*Pcor(i,j)+(relax)*(1/(A/dx+B/(dy*H/W))*(A/dx*Pcor(i+1,j)+B/(dy*H/W)* Pcor(i,j-1)+Q(i,j))); elseif(i==np-1&&j==2) Pcor(i,j)=(1relax)*Pcor(i,j)+(relax)*(1/(A/dx+B/(dy*H/W))*(A/dx*Pcor(i1,j)+B/(dy*H/W)*Pcor(i,j+1)+Q(i,j))); elseif(i==np-1&&j==np-1) Pcor(i,j)=(1relax)*Pcor(i,j)+(relax)*(1/(A/dx+B/(dy*H/W))*(+A/dx*Pcor(i1,j)+B/(dy*H/W)*Pcor(i,j-1)+Q(i,j))); elseif(i==2) Pcor(i,j)=(1relax)*Pcor(i,j)+(relax)*(1/(A/dx+2*B/(dy*H/W))*(A/dx*Pcor(i+1,j)... +B/(dy*H/W)*Pcor(i,j+1)+B/(dy*H/W)*Pcor(i,j-1)+Q(i,j))); elseif(i==np-1)
21
Pcor(i,j)=(1relax)*Pcor(i,j)+(relax)*(1/(A/dx+2*B/(dy*H/W))*(A/dx*Pcor(i-1,j)... +B/(dy*H/W)*Pcor(i,j+1)+B/(dy*H/W)*Pcor(i,j-1)+Q(i,j))); elseif(j==2) Pcor(i,j)=(1relax)*Pcor(i,j)+(relax)*(1/(2*A/dx+B/(dy*H/W))*(A/dx*Pcor(i+1,j)+A/dx*Pcor (i-1,j)... +B/(dy*H/W)*Pcor(i,j+1)+Q(i,j))); elseif(j==np-1) Pcor(i,j)=(1relax)*Pcor(i,j)+(relax)*(1/(2*A/dx+B/(dy*H/W))*(A/dx*Pcor(i+1,j)+A/dx*Pcor (i-1,j)... +B/(dy*H/W)*Pcor(i,j-1)+Q(i,j))); else Pcor(i,j)=(1relax)*Pcor(i,j)+(relax)*(1/(2*A/dx+2*B/(dy*H/W))*(A/dx*Pcor(i+1,j)+A/dx*Pc or(i-1,j)... +B/(dy*H/W)*Pcor(i,j+1)+B/(dy*H/W)*Pcor(i,j-1)+Q(i,j))); end end end delF=ones(1,(np-2)*(np-2)); c=1; for i=2:np-1 for j=2:np-1 if(i==2&&j==2) delF(c)=Pcor(i,j)1/(A/dx+B/(dy*H/W))*(A/dx*Pcor(i+1,j)+B/(dy*H/W)*Pcor(i,j+1)+Q(i,j)); elseif(i==2&&j==np-1) delF(c)=Pcor(i,j)1/(A/dx+B/(dy*H/W))*(A/dx*Pcor(i+1,j)+B/(dy*H/W)*Pcor(i,j-1)+Q(i,j)); elseif(i==np-1&&j==2) delF(c)=Pcor(i,j)-1/(A/dx+B/(dy*H/W))*(A/dx*Pcor(i1,j)+B/(dy*H/W)*Pcor(i,j+1)+Q(i,j)); elseif(i==np-1&&j==np-1) delF(c)=Pcor(i,j)-1/(A/dx+B/(dy*H/W))*(+A/dx*Pcor(i1,j)+B/(dy*H/W)*Pcor(i,j-1)+Q(i,j)); elseif(i==2) delF(c)=Pcor(i,j)1/(A/dx+2*B/(dy*H/W))*(A/dx*Pcor(i+1,j)... +B/(dy*H/W)*Pcor(i,j+1)+B/(dy*H/W)*Pcor(i,j-1)+Q(i,j)); elseif(i==np-1) delF(c)=Pcor(i,j)-1/(A/dx+2*B/(dy*H/W))*(A/dx*Pcor(i1,j)... +B/(dy*H/W)*Pcor(i,j+1)+B/(dy*H/W)*Pcor(i,j-1)+Q(i,j)); elseif(j==2) delF(c)=Pcor(i,j)1/(2*A/dx+B/(dy*H/W))*(A/dx*Pcor(i+1,j)+A/dx*Pcor(i-1,j)... +B/(dy*H/W)*Pcor(i,j+1)+Q(i,j)); elseif(j==np-1) delF(c)=Pcor(i,j)1/(2*A/dx+B/(dy*H/W))*(A/dx*Pcor(i+1,j)+A/dx*Pcor(i-1,j)... +B/(dy*H/W)*Pcor(i,j-1)+Q(i,j)); else delF(c)=Pcor(i,j)1/(2*A/dx+2*B/(dy*H/W))*(A/dx*Pcor(i+1,j)+A/dx*Pcor(i-1,j)... +B/(dy*H/W)*Pcor(i,j+1)+B/(dy*H/W)*Pcor(i,j-1)+Q(i,j)); end c=c+1;
22
end end sum=0; for c=1:length(delF) sum=sum+delF(c)^2; end residual1=sqrt(sum/length(delF)); end Pcor(1,:)=Pcor(2,:); Pcor(np,:)=Pcor(np-1,:); Pcor(:,np)=Pcor(:,np-1); Pcor(:,1)=Pcor(:,2); %display('Pcor calculation') %display(residual); Pcor; Ucor=zeros(np,nu); Vcor=zeros(nv,np); for i=2:nu-1 for j=2:np-1 Ucor(j,i)=(1/(1/dt+2*(1/Re)/dx^2+2*(1/Re)/(dy*H/W)^2))*... (-1/(dx)*(Pcor(j,i+1)-Pcor(j,i))); end end for i=2:np-1 for j=2:nv-1 Vcor(j,i)=(1/(1/dt+2*(1/Re)/dx^2+2*(1/Re)/(dy*H/W)^2))*... (-1/((dy*H/W))*(Pcor(j+1,i)-Pcor(j,i))); end end alp=0.1; al=1; Ustr=Ustr+al*Ucor; Vstr=Vstr+al*Vcor; Pstr=Pstr+alp*Pcor; %calculation of error in continuity from Ustr and Vstr-test % Q=zeros(np,np); % qsum=0; % for i=2:np-1 % for j=2:np-1 % Q(j,i)=(Ustr(j,i)-Ustr(j,i-1))/dx+(Vstr(j,i)-Vstr(j-1,i))/dy; % qsum=qsum+Q(j,i)^2; % end % end % qsum; end save('Re_V4_100_np_80_conv_o3','Ustr','Vstr','Pstr','convergence','k','Re', 'np','dx','dy')
23
8
References
1. Patankar, S. V.; Spalding, D. B. A calculation procedure for Heat, Mass and Momentum Transfer in Three-Dimensional Parabolic Flows. Int. J. Heat Mass Transfer 1972, 15, 1787-1806. 2. Ghia, U.; Ghia, K. N.; Shin, C. T. High-Re Solutions for Incompressible Flow using the Navier-Stokes Equations and a Multigrid Method. J. Comput. Phys. 1982, 48, 387-411. 3. Patankar, S. V. Numerical Heat Transfer and Fluid Flow; Hemisphere: New York, 1980. 4. Anderson, J. D. J. Computational Fluid Dynamics. The basics with applications.; McGraw-Hill: New York, 1995. 5. Fries, D. M.; Voitl, T.; Rohr, P. R. Liquid Extraction of Vanillin in Rectangular Microreactors. Chem. Eng. Technol. 2008, 31 (8), 1182-1187.
24
