Numerical Solution of Flow Field
From ThermalFluidsPedia
Line 2:  Line 2:  
Related:  Related:  
  +  ==Special Difficulties==  
  [[Staggered grid]]<  +  For a twodimensional incompressible flow problem without a body force, the continuity and momentum equations in the Cartesian coordinate system are 
  +  
  [[The SIMPLE Algorithm  +  { class="wikitable" border="0" 
+    
+   width="100%" <center>  
+  <math>\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}=0</math>  
+  </center>  
+  {{EquationRef(1)}}  
+  }  
+  
+  { class="wikitable" border="0"  
+    
+   width="100%" <center>  
+  <math>\rho u\frac{\partial u}{\partial x}+\rho v\frac{\partial u}{\partial y}=\frac{\partial }{\partial x}\left( \mu \frac{\partial u}{\partial x} \right)+\frac{\partial }{\partial y}\left( \mu \frac{\partial u}{\partial y} \right)\frac{\partial p}{\partial x}</math>  
+  </center>  
+  {{EquationRef(1)}}  
+  }  
+  
+  { class="wikitable" border="0"  
+    
+   width="100%" <center>  
+  <math>\rho u\frac{\partial v}{\partial x}+\rho v\frac{\partial v}{\partial y}=\frac{\partial }{\partial x}\left( \mu \frac{\partial v}{\partial x} \right)+\frac{\partial }{\partial y}\left( \mu \frac{\partial v}{\partial y} \right)\frac{\partial p}{\partial y}</math>  
+  </center>  
+  {{EquationRef(1)}}  
+  }  
+  It is observed that the momentum equations in each direction could be expressed in the same format as eq. (4.200) by replacing the general variable <math>\varphi </math> by the components of velocity in that particular (x, or y) direction. The source term for the momentum equations in the x and y directions can be expressed as <math>S_{u}=\partial p/\partial x</math> and <math>S_{v}=\partial p/\partial y,</math> respectively. It seems that we can directly apply the discretization schemes discussed in the preceding subsection in order to obtain the solution of the flow field. The only additional work necessary is to include the pressure gradient in each momentum equation.  
+  With the aid of Fig. 4.17, the pressure gradient in the xdirection can be expressed as  
+  
+  { class="wikitable" border="0"  
+    
+   width="100%" <center>  
+  <math>\left( \frac{\partial p}{\partial x} \right)_{P}=\frac{p_{e}p_{w}}{\Delta x}</math>  
+  </center>  
+  {{EquationRef(1)}}  
+  }  
+  If central difference is employed, the pressures at the faces of the control volume become  
+  
+  <math>p_{e}=\frac{p_{E}+p_{P}}{2},\text{ }p_{w}=\frac{p_{P}+p_{W}}{2}</math>  
+  
+  Substituting the above expressions into eq. (4.289), the pressure gradient in the xdirection becomes  
+  
+  { class="wikitable" border="0"  
+    
+   width="100%" <center>  
+  <math>\left( \frac{\partial p}{\partial x} \right)_{P}=\frac{p_{E}p_{W}}{2\Delta x}</math>  
+  </center>  
+  {{EquationRef(1)}}  
+  }  
+  Similarly, the pressure gradient in the ydirection can be expressed as  
+  
+  { class="wikitable" border="0"  
+    
+   width="100%" <center>  
+  <math>\left( \frac{\partial p}{\partial y} \right)_{P}=\frac{p_{N}p_{S}}{2\Delta y}</math>  
+  </center>  
+  {{EquationRef(1)}}  
+  }  
+  [[Image:Fig4.22.pngthumb400 pxalt=Checkerboard pressure field Figure 1: Checkerboard pressure field.]]  
+  It can be seen that the pressure gradient at point P is related to the pressures of the neighbor grid points and is not related to its own pressure. Thus, the pressure gradient is obtained by using the pressures of two alternative points, so that the accuracy of the pressure gradient for grid size Δx (or Δy) is equivalent to that obtained for coarse grid size 2Δx (or 2Δy). Another consequence of eqs. (4.290) and (4.291) is that if we have a pressure distribution as shown in Fig. 1 – referred to as a checkerboard pressure field – the discretization scheme represented by eqs. (4.290) and (4.291) will obtain <math>\partial p/\partial x=0</math> and <math>\partial p/\partial y=0</math> throughout the computational domain, i.e., eqs. (4.290) and (4.291) cannot recognize the difference between a checkerboard pressure field and a uniform pressure field. In other words, if a checkerboard pressure field is supposed to be the real pressure field for any reason, the discretization scheme based on eqs. (4.290) and (4.291) can neither detect nor eliminate such an unrealistic effect.  
+  
+  While it is straightforward to use the momentum equations (4.287) and (4.288) to solve for the velocity components in the x and ydirections, the pressure only appears as a source term in the momentum equations and there is not a separate equation for pressure itself. Since the continuity equation, (4.286), is not used, it would be logical to use the continuity equation to get the pressure field. Although a correct pressure field will yield a velocity field that satisfies the continuity equation, an algorithm must be developed to obtain the correct pressure field. The checkerboard pressure field problem and the lack of an equation for pressure are associated with the coupling between the pressure and velocity. An unrealistic (e.g., checkerboard) pressure field can be obtained if such a coupling relationship is not reflected in the algorithm (e.g., using central difference scheme to discretize the pressure gradient in the momentum equations).  
+  
+  ==Staggered grid==  
+  {{Comp Method for Forced Convection Category}}  
+  [[Image:Fig4.23.pngthumb400 pxalt=Staggered grid: (a) control volume for all other variables, (b) control volume for u, and (c) control volume for v. Figure 1: Staggered grid: (a) control volume for all other variables, (b) control volume for u, and (c) control volume for v.]]  
+  
+  To overcome the checkerboard pressure field and develop an effective algorithm for the pressure field, one can use a staggered grid system (Patankar and Spalding, 1972; Patankar, 1980) that stores the pressure and all other variables on the main grid but calculates the velocity at the face of the control volume (see Fig. 1). The control volume for the velocity component in the xdirection, u, is staggered from the control volume for all other variables to the right direction by half a grid. Similarly, the control volume for v is staggered up by a half grid. The discretized equations for u and v will be obtained by integrating the momentum equation in their control volumes shown in Fig. 1 (b) and (c), respectively.  
+  The discretization equations in the staggered grid system for all variables except velocity are the same as those presented in the preceding subsection. Since velocity is defined at the face of the control volume, the grid Peclet numbers can be directly calculated from the velocity component, i.e., no assumption on the velocity profile between the grid points is needed. For the momentum equation in the xdirection, we will need to integrate eq. (4.287) for the control volumes of u shown in Fig. 1 (b). The result can be expressed as  
+  
+  { class="wikitable" border="0"  
+    
+   width="100%" <center>  
+  <math>a_{e}u_{e}=\sum{a_{nb}u_{nb}+b+A_{e}(p_{P}p_{E})}</math>  
+  </center>  
+  {{EquationRef(1)}}  
+  }  
+  where the first term on the right hand side represents the summation of all the neighbor points of e. The effect of pressure has been separated from the source term, and the pressures at P and E were used to calculate the velocity ue. For a twodimensional problem, the area on which the pressure difference acts on is <math>A_{e}=\Delta y</math>.  
+  Similarly, the discretization equation for v can be obtained by integrating eq.(4.287) for the control volume of vn shown in Fig. 1(c), i.e.,  
+  
+  { class="wikitable" border="0"  
+    
+   width="100%" <center>  
+  <math>a_{n}v_{n}=\sum{a_{nb}v_{nb}+b+A_{n}(p_{P}p_{N})}</math>  
+  </center>  
+  {{EquationRef(2)}}  
+  }  
+  where the area on which pressure difference acts on is <math>A_{n}=\Delta x</math> for a twodimensional problem.  
+  For a threedimensional problem, eqs(4.292) and (4.293) are still valid except the neighbor points from the top and bottom should be included in the first term on the righthand side. The areas on which the pressure difference act on should be modified to <math>A_{e}=\Delta y\Delta z</math> and <math>A_{n}=\Delta x\Delta z</math>. In addition to eqs. (4.292) and (4.293), another equation for the velocity component in the zdirection, wt, is also needed, i.e.,  
+  
+  { class="wikitable" border="0"  
+    
+   width="100%" <center>  
+  <math>a_{t}w_{t}=\sum{a_{nb}w_{nb}+b+A_{t}(p_{P}p_{T})}</math>  
+  </center>  
+  {{EquationRef(3)}}  
+  }  
+  which is obtained by integrating the momentum equation in the zdirection for the control volume of wt that are staggered to the positive zdirection by a half grid. The area on which the pressure difference, <math>p_{P}p_{T}</math>, acts on is <math>A_{t}=\Delta x\Delta y</math>.  
+  
+  ==Pressure Correction Equation==  
+  {{Comp Method for Forced Convection Category}}  
+  The ability (of one) to solve eqs. (4.292) – (4.294) to obtain the velocity field relies on the accuracy of the pressure field, which is unknown a priori. If we assume that the pressure field is p* – which can be obtained from the previous time step, iteration, or guessing – the corresponding velocity field can be obtained by solving eqs. (4.292) – (4.294).  
+  
+  { class="wikitable" border="0"  
+    
+   width="100%" <center>  
+  <math>a_{e}u_{e}^{*}=\sum{a_{nb}u_{nb}^{*}+b+A_{e}(p_{P}^{*}p_{E}^{*})}</math>  
+  </center>  
+  {{EquationRef(1)}}  
+  }  
+  
+  { class="wikitable" border="0"  
+    
+   width="100%" <center>  
+  <math>a_{n}v_{n}^{*}=\sum{a_{nb}v_{nb}^{*}+b+A_{n}(p_{P}^{*}p_{N}^{*})}</math>  
+  </center>  
+  {{EquationRef(2)}}  
+  }  
+  
+  { class="wikitable" border="0"  
+    
+   width="100%" <center>  
+  <math>a_{t}w_{t}^{*}=\sum{a_{nb}w_{nb}^{*}+b+A_{t}(p_{P}^{*}p_{T}^{*})}</math>  
+  </center>  
+  {{EquationRef(3)}}  
+  }  
+  Since the velocity field obtained from eqs. (4.295) – (4.297) is based on a guessed pressure field, it may not satisfy the continuity equation. An algorithm based on pressure correction (Patankar, 1980) can be used to obtain the correct pressure field and consequently, the velocity field. If the correct pressure can be obtained by  
+  
+  { class="wikitable" border="0"  
+    
+   width="100%" <center>  
+  <math>\begin{matrix}{}\\\end{matrix}p=p^{*}+{p}'</math>  
+  </center>  
+  {{EquationRef(4)}}  
+  }  
+  where <math>{p}'</math> is pressure correction, the new velocity field based on the corrected pressure field will be  
+  
+  { class="wikitable" border="0"  
+    
+   width="100%" <center>  
+  <math>\begin{matrix}{}\\\end{matrix}u=u^{*}+{u}',\text{ }v=v^{*}+{v}',\text{ }w=w^{*}+{w}'</math>  
+  </center>  
+  {{EquationRef(5)}}  
+  }  
+  where <math>{u}',\text{ }{v}',\text{ and }{w}'</math> are velocity corrections. The corrected velocity field satisfies eqs(4.292) – (4.294), i.e.,  
+  
+  <math>a_{e}(u_{e}^{*}+{u}'_{e})=\sum{a_{nb}(u_{nb}^{*}+{u}'_{nb})+b+A_{e}[(p_{P}^{*}+{p}'_{P})(p_{E}^{*}+{p}'_{E})}]</math>  
+  
+  
+  <math>a_{n}(v_{n}^{*}+{v}'_{n})=\sum{a_{nb}(v_{nb}^{*}+{v}'_{nb})+b+A_{n}[(p_{P}^{*}+{p}'_{P})(p_{N}^{*}+{p}'_{N})}]</math>  
+  
+  
+  <math>a_{t}(w_{t}^{*}+{w}'_{t})=\sum{a_{nb}(w_{nb}^{*}+{w}'_{nb})+b+A_{t}[(p_{P}^{*}+{p}'_{P})(p_{T}^{*}+{p}'_{T})}]</math>  
+  
+  Subtracting eqs. (4.295) – (4.297) from the above three equations yields the following equations for the velocity corrections:  
+  
+  { class="wikitable" border="0"  
+    
+   width="100%" <center>  
+  <math>a_{e}{u}'_{e}=\sum{a_{nb}{u}'_{nb}+A_{e}({p}'_{P}{p}'_{E})}</math>  
+  </center>  
+  {{EquationRef(6)}}  
+  }  
+  
+  { class="wikitable" border="0"  
+    
+   width="100%" <center>  
+  <math>a_{n}{v}'_{n}=\sum{a_{nb}{v}'_{nb}+A_{n}({p}'_{P}{p}'_{N})}</math>  
+  </center>  
+  {{EquationRef(7)}}  
+  }  
+  
+  { class="wikitable" border="0"  
+    
+   width="100%" <center>  
+  <math>a_{t}{w}'_{t}=\sum{a_{nb}{w}'_{nb}+A_{t}({p}'_{P}{p}'_{T})}</math>  
+  </center>  
+  {{EquationRef(8)}}  
+  }  
+  It can be seen from eqs. (4.300) – (4.302) that the correction of the velocity can be from (1) the difference on pressure correction at two neighbor points in the main grid (the second terms), and (2) velocity correction from all neighbor points (the first terms). The solution of the velocity correction from eqs. (4.300) – (4.302) will be very complicated because the velocity correction at neighbor points will be related to pressure corrections at more grid points, and eventually velocity correction at any point will be related to the pressure corrections in the entire computational domain. Assuming the contribution of the latter is negligible, which is equivalent to assuming <math>a_{nb}=0</math> in eqs. (4.300) – (4.302), the velocity correction becomes  
+  
+  { class="wikitable" border="0"  
+    
+   width="100%" <center>  
+  <math>\begin{matrix}{}\\\end{matrix}{u}'_{e}=d_{e}({p}'_{P}{p}'_{E}),\text{ }{v}'_{n}=d_{n}({p}'_{P}{p}'_{N}),\text{ }{w}'_{t}=d_{t}({p}'_{P}{p}'_{T})</math>  
+  </center>  
+  {{EquationRef(9)}}  
+  }  
+  where  
+  
+  { class="wikitable" border="0"  
+    
+   width="100%" <center>  
+  <math>\begin{matrix}{}\\\end{matrix}d_{e}=A_{e}/a_{e},\text{ }d_{n}=A_{n}/a_{n},\text{ }d_{t}=A_{t}/a_{t}</math>  
+  </center>  
+  {{EquationRef(10)}}  
+  }  
+  Substituting eq. (4.303) into eq. (4.299), we have  
+  
+  { class="wikitable" border="0"  
+    
+   width="100%" <center>  
+  <math>u_{e}=u_{e}^{*}+d_{e}({p}'_{P}{p}'_{E})</math>  
+  </center>  
+  {{EquationRef(11)}}  
+  }  
+  
+  { class="wikitable" border="0"  
+    
+   width="100%" <center>  
+  <math>v_{n}=v_{n}^{*}+d_{n}({p}'_{P}{p}'_{N})</math>  
+  </center>  
+  {{EquationRef(12)}}  
+  }  
+  
+  { class="wikitable" border="0"  
+    
+   width="100%" <center>  
+  <math>w_{t}=w_{t}^{*}+d_{t}({p}'_{P}{p}'_{T})</math>  
+  </center>  
+  {{EquationRef(13)}}  
+  }  
+  which can be used to obtain the corrected velocity field from the pressure correction.  
+  The assumption of neglecting the first terms in eqs. (4.300) – (4.302) significantly simplifies the equations in order to obtain the corrected velocity field. If, however, this simplification results in unacceptable error in the final velocity field, it may not be used. As the iteration progresses, <math>v^{*}</math>, gradually approaches the true velocity, the velocity correction at each point approaches zero. After the converged velocity is obtained, the velocity correction from all the neighbor points will be zero. Therefore, neglecting the first terms in eqs. (4.300) – (4.302) will not affect the ultimate velocity field. It will, however, affect how quickly a converged solution can be obtained.  
+  We now turn our attention to the algebraic equation for the pressure correction. Since the starred velocity field (u*, v*, and w*) does not satisfy the continuity equation, we hope that the corrected velocity will satisfy the continuity equation. For a threedimensional flow, the continuity equation is  
+  
+  { class="wikitable" border="0"  
+    
+   width="100%" <center>  
+  <math>\frac{\partial \rho }{\partial t}+\frac{\partial (\rho u)}{\partial x}+\frac{\partial (\rho v)}{\partial y}+\frac{\partial (\rho w)}{\partial z}=0</math>  
+  </center>  
+  {{EquationRef(14)}}  
+  }  
+  Integrating eq. (4.308) over the main control volume [see Fig. 4.23(a) for the twodimensional projection of the threedimensional control volume] and using a fullyimplicit scheme, we obtain  
+  
+  <math>\begin{align}  
+  & \frac{\rho _{P}\rho _{P}^{0}}{\Delta t}\Delta x\Delta y\Delta z+[(\rho u)_{e}(\rho u)_{w}]\Delta y\Delta z \\  
+  & +[(\rho v)_{n}(\rho v)_{s}]\Delta x\Delta z+[(\rho w)_{t}(\rho w)_{b}]\Delta x\Delta y=0 \\  
+  \end{align}</math>  
+  
+  Substituting eqs. (4.305) – (4.307) into the above equation, the following algebraic equation for pressure correction is obtained:  
+  
+  { class="wikitable" border="0"  
+    
+   width="100%" <center>  
+  <math>\begin{matrix}{}\\\end{matrix}a_{P}{p}'_{P}=a_{E}{p}'_{E}+a_{W}{p}'_{W}+a_{N}{p}'_{N}+a_{S}{p}'_{S}+a_{T}{p}'_{T}+a_{B}{p}'_{B}+b</math>  
+  </center>  
+  {{EquationRef(15)}}  
+  }  
+  where  
+  
+  { class="wikitable" border="0"  
+    
+   width="100%" <center>  
+  <math>\begin{matrix}{}\\\end{matrix}a_{E}=\rho _{e}d_{e}\Delta y\Delta z,\text{ }a_{W}=\rho _{w}d_{w}\Delta y\Delta z</math>  
+  </center>  
+  {{EquationRef(16)}}  
+  }  
+  
+  { class="wikitable" border="0"  
+    
+   width="100%" <center>  
+  <math>\begin{matrix}{}\\\end{matrix}a_{N}=\rho _{n}d_{n}\Delta z\Delta x,\text{ }a_{S}=\rho _{s}d_{s}\Delta z\Delta x</math>  
+  </center>  
+  {{EquationRef(17)}}  
+  }  
+  
+  { class="wikitable" border="0"  
+    
+   width="100%" <center>  
+  <math>\begin{matrix}{}\\\end{matrix}a_{T}=\rho _{t}d_{t}\Delta x\Delta y,\text{ }a_{B}=\rho _{b}d_{b}\Delta x\Delta y</math>  
+  </center>  
+  {{EquationRef(18)}}  
+  }  
+  
+  { class="wikitable" border="0"  
+    
+   width="100%" <center>  
+  <math>\begin{matrix}{}\\\end{matrix}a_{P}=a_{E}+a_{W}+a_{N}+a_{S}+a_{T}+a_{B}</math>  
+  </center>  
+  {{EquationRef(19)}}  
+  }  
+  { class="wikitable" border="0"  
+    
+   width="100%" <center>  
+  <math>\begin{align}  
+  & b=\frac{\rho _{P}^{0}\rho _{P}}{\Delta t}\Delta x\Delta y\Delta z+[(\rho u^{*})_{w}(\rho u^{*})_{e}]\Delta y\Delta z \\  
+  & +[(\rho v^{*})_{s}(\rho v^{*})_{n}]\Delta x\Delta z+[(\rho w^{*})_{b}(\rho w^{*})_{t}]\Delta x\Delta y \\  
+  \end{align}</math>  
+  </center>  
+  {{EquationRef(20)}}  
+  }  
+  
+  [[Image:Fig4.24.pngthumb400 pxalt=Control volume for continuity equation at boundary Figure 1: Control volume for continuity equation at boundary.]]  
+  
+  If the starred velocity terms (u*, v*, and w*) already satisfy the continuity equation, we have <math>b=0</math> from eq. (4.314). Therefore, b defined in eq. (4.314) can be viewed as the negative value of the “residual” mass for grid point P. If the continuity equation is satisfied, such residual mass should be zero. Practically, the maximum absolute value of b in the entire computational domain as well as the absolute value of the summation of b in the entire domain are used as criteria for the convergence.  
+  In order to apply eq(4.309) to obtain the pressure correction, the boundary conditions for pressure correction have to be specified. If the pressure at the boundary is given, no correction for the pressure at the boundary is needed so <math>{p}'=0</math> at the boundary. If the boundary velocity is given (see Fig. 1), the velocity correction for the boundary velocity (ue in Fig. 1) will be zero. It follows from eq. (4.303) that de for the control volume shown in Fig. 1 has to be zero. Therefore, one can set <math>a_{E}=0</math> [see eq. (4.310)] if ue in Fig. 1 is given.  
+  
+  ==The SIMPLE Algorithm==  
+  {{Comp Method for Forced Convection Category}}  
+  The above approach for the solution of the incompressible flow field was named SIMPLE, which stands for SemiImplicit Method for PressureLinked Equations. It was originally proposed by Patankar and Spalding (1972) and summarized in Patankar (1980). This algorithm is a semiimplicit method because the first terms on the righthand side of eqs. (4.300) – (4.302) are neglected. If these terms are retained, we will have to solve for the velocity corrections ( ) for the entire flow field simultaneously and the algorithm will become fullyimplicit. The procedure for SIMPLE algorithm can be summarized as following:<BR>  
+  1. Guess a pressure field, p*.<BR>  
+  2. Obtain the starred velocity field, (u*, v*, and w*) from eqs. (4.295) – (4.297). <BR>  
+  3. Solve the pressure correction equation (4.309) to get .<BR>  
+  4. Obtain a new pressure field from eq. (4.298).<BR>  
+  5. Improve the velocity field using eqs. (4.305) – (4.307). <BR>  
+  6. Obtain solutions for other ’s from eq. (4.274). If the flow field is not affected by a particular , it should be solved after a converged solution for flow field is obtained.<BR>  
+  7. Treat the corrected pressure in Step 4 as a new starred pressure and go back to step 2. The iteration procedure is repeated until a converged solution is obtained. <BR>  
+  Since the invention of the SIMPLE algorithm in 1972, it has evolved into the classic approach for computational fluid dynamics and heat transfer. Although there are several improved versions in the literature, SIMPLE still remains one of basic the most powerful algorithms. The objective of this section is to introduce the readers to the world of computational fluid dynamics and heat transfer without the overwhelming numerical details. Additional information can be found in the seminal books of Patankar (1980) and Tao (2001), as well as the Handbook of Numerical Heat Transfer (Minkowycz et al., 2006). 
Revision as of 06:32, 21 July 2010
Computational methodologies for forced convection 
Related:
Contents 
Special Difficulties
For a twodimensional incompressible flow problem without a body force, the continuity and momentum equations in the Cartesian coordinate system are



It is observed that the momentum equations in each direction could be expressed in the same format as eq. (4.200) by replacing the general variable by the components of velocity in that particular (x, or y) direction. The source term for the momentum equations in the x and y directions can be expressed as and respectively. It seems that we can directly apply the discretization schemes discussed in the preceding subsection in order to obtain the solution of the flow field. The only additional work necessary is to include the pressure gradient in each momentum equation. With the aid of Fig. 4.17, the pressure gradient in the xdirection can be expressed as

If central difference is employed, the pressures at the faces of the control volume become
Substituting the above expressions into eq. (4.289), the pressure gradient in the xdirection becomes

Similarly, the pressure gradient in the ydirection can be expressed as

It can be seen that the pressure gradient at point P is related to the pressures of the neighbor grid points and is not related to its own pressure. Thus, the pressure gradient is obtained by using the pressures of two alternative points, so that the accuracy of the pressure gradient for grid size Δx (or Δy) is equivalent to that obtained for coarse grid size 2Δx (or 2Δy). Another consequence of eqs. (4.290) and (4.291) is that if we have a pressure distribution as shown in Fig. 1 – referred to as a checkerboard pressure field – the discretization scheme represented by eqs. (4.290) and (4.291) will obtain and throughout the computational domain, i.e., eqs. (4.290) and (4.291) cannot recognize the difference between a checkerboard pressure field and a uniform pressure field. In other words, if a checkerboard pressure field is supposed to be the real pressure field for any reason, the discretization scheme based on eqs. (4.290) and (4.291) can neither detect nor eliminate such an unrealistic effect.
While it is straightforward to use the momentum equations (4.287) and (4.288) to solve for the velocity components in the x and ydirections, the pressure only appears as a source term in the momentum equations and there is not a separate equation for pressure itself. Since the continuity equation, (4.286), is not used, it would be logical to use the continuity equation to get the pressure field. Although a correct pressure field will yield a velocity field that satisfies the continuity equation, an algorithm must be developed to obtain the correct pressure field. The checkerboard pressure field problem and the lack of an equation for pressure are associated with the coupling between the pressure and velocity. An unrealistic (e.g., checkerboard) pressure field can be obtained if such a coupling relationship is not reflected in the algorithm (e.g., using central difference scheme to discretize the pressure gradient in the momentum equations).
Staggered grid
Computational methodologies for forced convection 
To overcome the checkerboard pressure field and develop an effective algorithm for the pressure field, one can use a staggered grid system (Patankar and Spalding, 1972; Patankar, 1980) that stores the pressure and all other variables on the main grid but calculates the velocity at the face of the control volume (see Fig. 1). The control volume for the velocity component in the xdirection, u, is staggered from the control volume for all other variables to the right direction by half a grid. Similarly, the control volume for v is staggered up by a half grid. The discretized equations for u and v will be obtained by integrating the momentum equation in their control volumes shown in Fig. 1 (b) and (c), respectively. The discretization equations in the staggered grid system for all variables except velocity are the same as those presented in the preceding subsection. Since velocity is defined at the face of the control volume, the grid Peclet numbers can be directly calculated from the velocity component, i.e., no assumption on the velocity profile between the grid points is needed. For the momentum equation in the xdirection, we will need to integrate eq. (4.287) for the control volumes of u shown in Fig. 1 (b). The result can be expressed as

where the first term on the right hand side represents the summation of all the neighbor points of e. The effect of pressure has been separated from the source term, and the pressures at P and E were used to calculate the velocity ue. For a twodimensional problem, the area on which the pressure difference acts on is A_{e} = Δy. Similarly, the discretization equation for v can be obtained by integrating eq.(4.287) for the control volume of vn shown in Fig. 1(c), i.e.,

where the area on which pressure difference acts on is A_{n} = Δx for a twodimensional problem. For a threedimensional problem, eqs(4.292) and (4.293) are still valid except the neighbor points from the top and bottom should be included in the first term on the righthand side. The areas on which the pressure difference act on should be modified to A_{e} = ΔyΔz and A_{n} = ΔxΔz. In addition to eqs. (4.292) and (4.293), another equation for the velocity component in the zdirection, wt, is also needed, i.e.,

which is obtained by integrating the momentum equation in the zdirection for the control volume of wt that are staggered to the positive zdirection by a half grid. The area on which the pressure difference, p_{P} − p_{T}, acts on is A_{t} = ΔxΔy.
Pressure Correction Equation
Computational methodologies for forced convection 
The ability (of one) to solve eqs. (4.292) – (4.294) to obtain the velocity field relies on the accuracy of the pressure field, which is unknown a priori. If we assume that the pressure field is p* – which can be obtained from the previous time step, iteration, or guessing – the corresponding velocity field can be obtained by solving eqs. (4.292) – (4.294).



Since the velocity field obtained from eqs. (4.295) – (4.297) is based on a guessed pressure field, it may not satisfy the continuity equation. An algorithm based on pressure correction (Patankar, 1980) can be used to obtain the correct pressure field and consequently, the velocity field. If the correct pressure can be obtained by

where p' is pressure correction, the new velocity field based on the corrected pressure field will be

where u', v', and w' are velocity corrections. The corrected velocity field satisfies eqs(4.292) – (4.294), i.e.,
Subtracting eqs. (4.295) – (4.297) from the above three equations yields the following equations for the velocity corrections:



It can be seen from eqs. (4.300) – (4.302) that the correction of the velocity can be from (1) the difference on pressure correction at two neighbor points in the main grid (the second terms), and (2) velocity correction from all neighbor points (the first terms). The solution of the velocity correction from eqs. (4.300) – (4.302) will be very complicated because the velocity correction at neighbor points will be related to pressure corrections at more grid points, and eventually velocity correction at any point will be related to the pressure corrections in the entire computational domain. Assuming the contribution of the latter is negligible, which is equivalent to assuming a_{nb} = 0 in eqs. (4.300) – (4.302), the velocity correction becomes

where

Substituting eq. (4.303) into eq. (4.299), we have



which can be used to obtain the corrected velocity field from the pressure correction. The assumption of neglecting the first terms in eqs. (4.300) – (4.302) significantly simplifies the equations in order to obtain the corrected velocity field. If, however, this simplification results in unacceptable error in the final velocity field, it may not be used. As the iteration progresses, v^{ * }, gradually approaches the true velocity, the velocity correction at each point approaches zero. After the converged velocity is obtained, the velocity correction from all the neighbor points will be zero. Therefore, neglecting the first terms in eqs. (4.300) – (4.302) will not affect the ultimate velocity field. It will, however, affect how quickly a converged solution can be obtained. We now turn our attention to the algebraic equation for the pressure correction. Since the starred velocity field (u*, v*, and w*) does not satisfy the continuity equation, we hope that the corrected velocity will satisfy the continuity equation. For a threedimensional flow, the continuity equation is

Integrating eq. (4.308) over the main control volume [see Fig. 4.23(a) for the twodimensional projection of the threedimensional control volume] and using a fullyimplicit scheme, we obtain
Substituting eqs. (4.305) – (4.307) into the above equation, the following algebraic equation for pressure correction is obtained:

where





If the starred velocity terms (u*, v*, and w*) already satisfy the continuity equation, we have b = 0 from eq. (4.314). Therefore, b defined in eq. (4.314) can be viewed as the negative value of the “residual” mass for grid point P. If the continuity equation is satisfied, such residual mass should be zero. Practically, the maximum absolute value of b in the entire computational domain as well as the absolute value of the summation of b in the entire domain are used as criteria for the convergence. In order to apply eq(4.309) to obtain the pressure correction, the boundary conditions for pressure correction have to be specified. If the pressure at the boundary is given, no correction for the pressure at the boundary is needed so p' = 0 at the boundary. If the boundary velocity is given (see Fig. 1), the velocity correction for the boundary velocity (ue in Fig. 1) will be zero. It follows from eq. (4.303) that de for the control volume shown in Fig. 1 has to be zero. Therefore, one can set a_{E} = 0 [see eq. (4.310)] if ue in Fig. 1 is given.
The SIMPLE Algorithm
Computational methodologies for forced convection 
The above approach for the solution of the incompressible flow field was named SIMPLE, which stands for SemiImplicit Method for PressureLinked Equations. It was originally proposed by Patankar and Spalding (1972) and summarized in Patankar (1980). This algorithm is a semiimplicit method because the first terms on the righthand side of eqs. (4.300) – (4.302) are neglected. If these terms are retained, we will have to solve for the velocity corrections ( ) for the entire flow field simultaneously and the algorithm will become fullyimplicit. The procedure for SIMPLE algorithm can be summarized as following:
1. Guess a pressure field, p*.
2. Obtain the starred velocity field, (u*, v*, and w*) from eqs. (4.295) – (4.297).
3. Solve the pressure correction equation (4.309) to get .
4. Obtain a new pressure field from eq. (4.298).
5. Improve the velocity field using eqs. (4.305) – (4.307).
6. Obtain solutions for other ’s from eq. (4.274). If the flow field is not affected by a particular , it should be solved after a converged solution for flow field is obtained.
7. Treat the corrected pressure in Step 4 as a new starred pressure and go back to step 2. The iteration procedure is repeated until a converged solution is obtained.
Since the invention of the SIMPLE algorithm in 1972, it has evolved into the classic approach for computational fluid dynamics and heat transfer. Although there are several improved versions in the literature, SIMPLE still remains one of basic the most powerful algorithms. The objective of this section is to introduce the readers to the world of computational fluid dynamics and heat transfer without the overwhelming numerical details. Additional information can be found in the seminal books of Patankar (1980) and Tao (2001), as well as the Handbook of Numerical Heat Transfer (Minkowycz et al., 2006).