Numerical Solution of Flow Field
From ThermalFluidsPedia
(7 intermediate revisions not shown)  
Line 26:  Line 26:  
{{EquationRef(3)}}  {{EquationRef(3)}}  
}  }  
  It is observed that the momentum equations in each direction could be expressed in the same format as  +  [[Image:Fig4.17.pngthumb400 pxalt=Control volume for onedimensional problem  Control volume for onedimensional problem.]] 
  With the aid of  +  It is observed that the momentum equations in each direction could be expressed in the same format as the first equation in [[Computational methodologies for forced convection]] 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 the figure to the right, the pressure gradient in the xdirection can be expressed as  
{ class="wikitable" border="0"  { class="wikitable" border="0"  
Line 38:  Line 39:  
If central difference is employed, the pressures at the faces of the control volume become  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>  +  <center><math>p_{e}=\frac{p_{E}+p_{P}}{2},\text{ }p_{w}=\frac{p_{P}+p_{W}}{2}</math></center> 
  Substituting the above expressions into eq. (4  +  Substituting the above expressions into eq. (4), the pressure gradient in the ''x''direction becomes 
{ class="wikitable" border="0"  { class="wikitable" border="0"  
Line 58:  Line 59:  
{{EquationRef(6)}}  {{EquationRef(6)}}  
}  }  
  [[Image:Fig4.22.pngthumb400 pxalt=Checkerboard pressure field   +  [[Image:Fig4.22.pngthumb400 pxalt=Checkerboard pressure field  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  +  
+  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. (5) and (6) is that if we have a pressure distribution as shown in the figure to the right – referred to as a checkerboard pressure field – the discretization scheme represented by eqs. (5) and (6) will obtain <math>\partial p/\partial x=0</math> and <math>\partial p/\partial y=0</math> throughout the computational domain, i.e., eqs. (5) and (6) 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. (5) and (6) can neither detect nor eliminate such an unrealistic effect.  
  While it is straightforward to use the momentum equations (  +  While it is straightforward to use the momentum equations (2) and (3) to solve for the velocity components in the ''x'' and ''y''directions, 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, (1), 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==  ==Staggered grid==  
  [[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.   +  [[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.  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  +  To overcome the checkerboard pressure field and develop an effective algorithm for the pressure field, one can use a staggered grid system <ref name="PS1972">Patankar, S.V., and Spalding, D.B., 1972, “A Calculation Procedure for Heat, and Mass and Momentum Transfer in ThreeDimensional Parabolic Flows,” International Journal of Heat and Mass Transfer, Vol. 17, pp. 17871806.</ref> <ref name="P1980">Patankar, S.V., 1980, Numerical Heat Transfer and Fluid Flow, Hemisphere, Washington, DC.</ref> that stores the pressure and all other variables on the main grid but calculates the velocity at the face of the control volume (see figure to the right). The control volume for the velocity component in the ''x''direction, ''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 (b) and (c) of the right figure, 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  +  
+  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 ''x''direction, we will need to integrate eq. (2) for the control volumes of ''u'' shown in (b) of the figure to the right. The result can be expressed as  
{ class="wikitable" border="0"  { class="wikitable" border="0"  
Line 77:  Line 80:  
}  }  
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>.  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.(  +  Similarly, the discretization equation for v can be obtained by integrating eq.(2) for the control volume of ''v''<sub>n</sub> shown in (c) of the right figure, i.e., 
{ class="wikitable" border="0"  { class="wikitable" border="0"  
Line 87:  Line 90:  
}  }  
where the area on which pressure difference acts on is <math>A_{n}=\Delta x</math> for a twodimensional problem.  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(  +  
+  For a threedimensional problem, eqs(7) and (8) 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. (7) and (8), another equation for the velocity component in the ''z''direction, ''w<sub>t</sub>'', is also needed, i.e.,  
{ class="wikitable" border="0"  { class="wikitable" border="0"  
Line 96:  Line 100:  
{{EquationRef(9)}}  {{EquationRef(9)}}  
}  }  
  which is obtained by integrating the momentum equation in the zdirection for the control volume of  +  which is obtained by integrating the momentum equation in the ''z''direction for the control volume of ''w<sub>t</sub>'' that are staggered to the positive ''z''direction 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==  ==Pressure Correction Equation==  
  The ability (of one) to solve eqs. (  +  The ability (of one) to solve eqs. (7) – (9) 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. (7) – (9). 
{ class="wikitable" border="0"  { class="wikitable" border="0"  
Line 124:  Line 128:  
{{EquationRef(12)}}  {{EquationRef(12)}}  
}  }  
  Since the velocity field obtained from eqs. (  +  Since the velocity field obtained from eqs. (10) – (12) is based on a guessed pressure field, it may not satisfy the continuity equation. An algorithm based on pressure correction <ref name="P1980"/> 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"  { class="wikitable" border="0"  
Line 142:  Line 146:  
{{EquationRef(14)}}  {{EquationRef(14)}}  
}  }  
  where <math>{u}',\text{ }{v}',\text{ and }{w}'</math> are velocity corrections. The corrected velocity field satisfies eqs(  +  where <math>{u}',\text{ }{v}',\text{ and }{w}'</math> are velocity corrections. The corrected velocity field satisfies eqs. (7) – (9), 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_{e}(u_{e}^{*}+{u}'_{e})=\sum{a_{nb}(u_{nb}^{*}+{u}'_{nb})+b+A_{e}[(p_{P}^{*}+{p}'_{P})(p_{E}^{*}+{p}'_{E})}]</math>  
Line 152:  Line 156:  
<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>  <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. (  +  Subtracting eqs. (10) – (12) from the above three equations yields the following equations for the velocity corrections: 
{ class="wikitable" border="0"  { class="wikitable" border="0"  
Line 177:  Line 181:  
{{EquationRef(17)}}  {{EquationRef(17)}}  
}  }  
  It can be seen from eqs. (  +  It can be seen from eqs. (15) – (17) 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. (15) – (17) 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. (15) – (17), the velocity correction becomes 
{ class="wikitable" border="0"  { class="wikitable" border="0"  
Line 195:  Line 199:  
{{EquationRef(19)}}  {{EquationRef(19)}}  
}  }  
  Substituting eq. (  +  Substituting eq. (18) into eq. (19), we have 
{ class="wikitable" border="0"  { class="wikitable" border="0"  
Line 220:  Line 224:  
{{EquationRef(22)}}  {{EquationRef(22)}}  
}  }  
  which can be used to obtain the corrected velocity field from the pressure correction.  +  which can be used to obtain the corrected velocity field from the pressure correction. 
  The assumption of neglecting the first terms in eqs. (  +  
  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  +  The assumption of neglecting the first terms in eqs. (15) – (17) 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. (15) – (17) 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"  { class="wikitable" border="0"  
Line 231:  Line 236:  
{{EquationRef(23)}}  {{EquationRef(23)}}  
}  }  
  Integrating eq. (  +  Integrating eq. (23) over the main control volume [see (a) in above figure for the twodimensional projection of the threedimensional control volume] and using a fullyimplicit scheme, we obtain 
  <math>\begin{align}  +  <center><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 \\  & \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 \\  & +[(\rho v)_{n}(\rho v)_{s}]\Delta x\Delta z+[(\rho w)_{t}(\rho w)_{b}]\Delta x\Delta y=0 \\  
  \end{align}</math>  +  \end{align}</math></center> 
  Substituting eqs. (  +  Substituting eqs. (20) – (22) into the above equation, the following algebraic equation for pressure correction is obtained: 
{ class="wikitable" border="0"  { class="wikitable" border="0"  
Line 291:  Line 296:  
}  }  
  [[Image:Fig4.24.pngthumb400 pxalt=Control volume for continuity equation at boundary   +  [[Image:Fig4.24.pngthumb400 pxalt=Control volume for continuity equation at boundary  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. (  +  If the starred velocity terms (''u''*, ''v''*, and ''w''*) already satisfy the continuity equation, we have <math>b=0</math> from eq. (29). Therefore, b defined in eq. (29) 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(  +  
+  In order to apply eq(24) 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 figure to the right), the velocity correction for the boundary velocity (''u<sub>e</sub>'' in figure to the right) will be zero. It follows from eq. (18) that de for the control volume shown in the figure has to be zero. Therefore, one can set <math>a_{E}=0</math> [see eq. (25)] if ''u<sub>e</sub>'' in the figure is given.  
==The SIMPLE Algorithm==  ==The SIMPLE Algorithm==  
  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  +  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 <ref name="PS1972"/> and summarized in Patankar <ref name="P1980"/>. This algorithm is a semiimplicit method because the first terms on the righthand side of eqs. (15) – (17) are neglected. If these terms are retained, we will have to solve for the velocity corrections (<math>{{{u}'}_{e}},\text{ }{{{v}'}_{n}},\text{ and }{{{w}'}_{t}}</math>) 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. (  +  1. Guess a pressure field, ''p''*.<BR> 
  3. Solve the pressure correction equation (  +  
  4. Obtain a new pressure field from eq. (  +  2. Obtain the starred velocity field, (''u''*, ''v''*, and ''w''*) from eqs. (10) – (12). <BR> 
  5. Improve the velocity field using eqs. (  +  
  6. Obtain solutions for other ’s from  +  3. Solve the pressure correction equation (24) to get <math>{p}'</math>.<BR> 
+  
+  4. Obtain a new pressure field from eq. (13).<BR>  
+  
+  5. Improve the velocity field using eqs. (20) – (22). <BR>  
+  
+  6. Obtain solutions for other <math>\phi </math>’s from the follow equation  
+  <center><math>\begin{matrix}  
+  {} & {{a}_{P}}{{\phi }_{P}}={{a}_{E}}{{\phi }_{E}}+{{a}_{W}}{{\phi }_{W}}+{{a}_{N}}{{\phi }_{N}}+{{a}_{S}}{{\phi }_{S}}+{{a}_{T}}{{\phi }_{T}}+{{a}_{B}}{{\phi }_{B}}+b \\  
+  \end{matrix}</math>  
+  </center>  
+  If the flow field is not affected by a particular <math>\phi </math>, 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>  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  +  
+  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 <ref name="P1980"/> and Tao <ref name="T2001">Tao 2001, W.Q., Numerical Heat Transfer, 2nd Ed., Xi’an Jiaotong University Press, Xi’an, China (in Chinese). </ref>, as well as the Handbook of Numerical Heat Transfer <ref name="M2006">Minkowycz, M.J., Sparrow, E.M., and Murthy, J.Y., eds., 2006, Handbook of Numerical Heat Transfer, 2nd ed. John Wiley & Sons, Hoboken, NJ.</ref><ref>Faghri, A., Zhang, Y., and Howell, J. R., 2010, Advanced Heat and Mass Transfer, Global Digital Press, Columbia, MO.</ref>.  
+  
+  ==References==  
+  {{Reflist}} 
Current revision as of 19:45, 23 July 2010
Computational methodologies for forced convection 
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 the first equation in Computational methodologies for forced convection 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 the figure to the right, 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), 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. (5) and (6) is that if we have a pressure distribution as shown in the figure to the right – referred to as a checkerboard pressure field – the discretization scheme represented by eqs. (5) and (6) will obtain and throughout the computational domain, i.e., eqs. (5) and (6) 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. (5) and (6) can neither detect nor eliminate such an unrealistic effect.
While it is straightforward to use the momentum equations (2) and (3) 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, (1), 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
To overcome the checkerboard pressure field and develop an effective algorithm for the pressure field, one can use a staggered grid system ^{[1]} ^{[2]} that stores the pressure and all other variables on the main grid but calculates the velocity at the face of the control volume (see figure to the right). 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 (b) and (c) of the right figure, 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. (2) for the control volumes of u shown in (b) of the figure to the right. 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.(2) for the control volume of v_{n} shown in (c) of the right figure, i.e.,

where the area on which pressure difference acts on is A_{n} = Δx for a twodimensional problem.
For a threedimensional problem, eqs(7) and (8) 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. (7) and (8), another equation for the velocity component in the zdirection, w_{t}, is also needed, i.e.,

which is obtained by integrating the momentum equation in the zdirection for the control volume of w_{t} 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
The ability (of one) to solve eqs. (7) – (9) 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. (7) – (9).



Since the velocity field obtained from eqs. (10) – (12) is based on a guessed pressure field, it may not satisfy the continuity equation. An algorithm based on pressure correction ^{[2]} 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. (7) – (9), i.e.,
Subtracting eqs. (10) – (12) from the above three equations yields the following equations for the velocity corrections:



It can be seen from eqs. (15) – (17) 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. (15) – (17) 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. (15) – (17), the velocity correction becomes

where

Substituting eq. (18) into eq. (19), 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. (15) – (17) 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. (15) – (17) 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. (23) over the main control volume [see (a) in above figure for the twodimensional projection of the threedimensional control volume] and using a fullyimplicit scheme, we obtain
Substituting eqs. (20) – (22) 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. (29). Therefore, b defined in eq. (29) 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(24) 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 figure to the right), the velocity correction for the boundary velocity (u_{e} in figure to the right) will be zero. It follows from eq. (18) that de for the control volume shown in the figure has to be zero. Therefore, one can set a_{E} = 0 [see eq. (25)] if u_{e} in the figure is given.
The SIMPLE Algorithm
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 ^{[1]} and summarized in Patankar ^{[2]}. This algorithm is a semiimplicit method because the first terms on the righthand side of eqs. (15) – (17) are neglected. If these terms are retained, we will have to solve for the velocity corrections (u'_{e}, v'_{n}, and w'_{t}) 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. (10) – (12).
3. Solve the pressure correction equation (24) to get p'.
4. Obtain a new pressure field from eq. (13).
5. Improve the velocity field using eqs. (20) – (22).
6. Obtain solutions for other φ’s from the follow equation
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 ^{[2]} and Tao ^{[3]}, as well as the Handbook of Numerical Heat Transfer ^{[4]}^{[5]}.
References
 ↑ ^{1.0} ^{1.1} Patankar, S.V., and Spalding, D.B., 1972, “A Calculation Procedure for Heat, and Mass and Momentum Transfer in ThreeDimensional Parabolic Flows,” International Journal of Heat and Mass Transfer, Vol. 17, pp. 17871806.
 ↑ ^{2.0} ^{2.1} ^{2.2} ^{2.3} Patankar, S.V., 1980, Numerical Heat Transfer and Fluid Flow, Hemisphere, Washington, DC.
 ↑ Tao 2001, W.Q., Numerical Heat Transfer, 2nd Ed., Xi’an Jiaotong University Press, Xi’an, China (in Chinese).
 ↑ Minkowycz, M.J., Sparrow, E.M., and Murthy, J.Y., eds., 2006, Handbook of Numerical Heat Transfer, 2nd ed. John Wiley & Sons, Hoboken, NJ.
 ↑ Faghri, A., Zhang, Y., and Howell, J. R., 2010, Advanced Heat and Mass Transfer, Global Digital Press, Columbia, MO.