# Numerical Solution of Flow Field

(Difference between revisions)
 Revision as of 06:37, 21 July 2010 (view source)← Older edit Current revision as of 19:45, 23 July 2010 (view source) (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 eq. (4.200) by replacing the general variable $\varphi$ 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 $S_{u}=-\partial p/\partial x$ and $S_{v}=-\partial p/\partial y,$ 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. + [[Image:Fig4.17.png|thumb|400 px|alt=Control volume for one-dimensional problem | Control volume for one-dimensional problem.]] - With the aid of Fig. 4.17, the pressure gradient in the x-direction can be expressed as + 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 $\varphi$ 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 $S_{u}=-\partial p/\partial x$ and $S_{v}=-\partial p/\partial y,$ 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 x-direction 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 - $p_{e}=\frac{p_{E}+p_{P}}{2},\text{ }p_{w}=\frac{p_{P}+p_{W}}{2}$ +
$p_{e}=\frac{p_{E}+p_{P}}{2},\text{ }p_{w}=\frac{p_{P}+p_{W}}{2} - Substituting the above expressions into eq. (4.289), the pressure gradient in the x-direction becomes + 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.png|thumb|400 px|alt=Checkerboard pressure field |Figure 1: Checkerboard pressure field.]] + [[Image:Fig4.22.png|thumb|400 px|alt=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 Δ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 [itex]\partial p/\partial x=0$ and $\partial p/\partial y=0$ 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. + + 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 $\partial p/\partial x=0$ and $\partial p/\partial y=0$ 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 (4.287) and (4.288) 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, (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). + 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.png|thumb|400 px|alt=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.]] + [[Image:Fig4.23.png|thumb|400 px|alt=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 (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 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 Fig. 1 (b) and (c), respectively. + To overcome the checkerboard pressure field and develop an effective algorithm for the pressure field, one can use a staggered grid system Patankar, S.V., and Spalding, D.B., 1972, “A Calculation Procedure for Heat, and Mass and Momentum Transfer in Three-Dimensional Parabolic Flows,” International Journal of Heat and Mass Transfer, Vol. 17, pp. 1787-1806. Patankar, S.V., 1980, Numerical Heat Transfer and Fluid Flow, Hemisphere, Washington, DC. 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 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. (4.287) for the control volumes of u shown in Fig. 1 (b).  The result can be expressed as + + 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 two-dimensional problem, the area on which the pressure difference acts on is $A_{e}=\Delta y$. 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 two-dimensional problem, the area on which the pressure difference acts on is $A_{e}=\Delta 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., + 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., {| class="wikitable" border="0" {| class="wikitable" border="0" Line 87: Line 90: |} |} where the area on which pressure difference acts on is $A_{n}=\Delta x$ for a two-dimensional problem. where the area on which pressure difference acts on is $A_{n}=\Delta x$ for a two-dimensional problem. - For a three-dimensional 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 right-hand side. The areas on which the pressure difference act on should be modified to $A_{e}=\Delta y\Delta z$ and $A_{n}=\Delta x\Delta z$. In addition to eqs. (4.292) and (4.293), another equation for the velocity component in the z-direction, wt, is also needed, i.e., + + For a three-dimensional 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 right-hand side. The areas on which the pressure difference act on should be modified to $A_{e}=\Delta y\Delta z$ and $A_{n}=\Delta x\Delta z$. In addition to eqs. (7) and (8), another equation for the velocity component in the ''z''-direction, ''wt'', 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 z-direction for the control volume of wt  that are staggered to the positive z-direction by a half grid. The area on which the pressure difference, $p_{P}-p_{T}$, acts on is $A_{t}=\Delta x\Delta y$. + which is obtained by integrating the momentum equation in the ''z''-direction for the control volume of ''wt'' that are staggered to the positive ''z''-direction by a half grid. The area on which the pressure difference, $p_{P}-p_{T}$, acts on is $A_{t}=\Delta x\Delta y$. ==Pressure Correction Equation== ==Pressure Correction Equation== - 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). + 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. (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 + 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 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 ${u}',\text{ }{v}',\text{ and }{w}'$ are velocity corrections. The corrected velocity field satisfies eqs(4.292) – (4.294), i.e., + where ${u}',\text{ }{v}',\text{ and }{w}'$ are velocity corrections. The corrected velocity field satisfies eqs. (7) – (9), i.e., $a_{e}(u_{e}^{*}+{u}'_{e})=\sum{a_{nb}(u_{nb}^{*}+{u}'_{nb})+b+A_{e}[(p_{P}^{*}+{p}'_{P})-(p_{E}^{*}+{p}'_{E})}]$ $a_{e}(u_{e}^{*}+{u}'_{e})=\sum{a_{nb}(u_{nb}^{*}+{u}'_{nb})+b+A_{e}[(p_{P}^{*}+{p}'_{P})-(p_{E}^{*}+{p}'_{E})}]$ Line 152: Line 156: $a_{t}(w_{t}^{*}+{w}'_{t})=\sum{a_{nb}(w_{nb}^{*}+{w}'_{nb})+b+A_{t}[(p_{P}^{*}+{p}'_{P})-(p_{T}^{*}+{p}'_{T})}]$ $a_{t}(w_{t}^{*}+{w}'_{t})=\sum{a_{nb}(w_{nb}^{*}+{w}'_{nb})+b+A_{t}[(p_{P}^{*}+{p}'_{P})-(p_{T}^{*}+{p}'_{T})}]$ - Subtracting eqs. (4.295) – (4.297) from the above three equations yields the following equations for the velocity corrections: + 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. (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 + 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 {| class="wikitable" border="0" {| class="wikitable" border="0" Line 195: Line 199: |{{EquationRef|(19)}} |{{EquationRef|(19)}} |} |} - Substituting eq. (4.303) into eq. (4.299), we have + 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. (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 three-dimensional 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, $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 three-dimensional flow, the continuity equation is {| class="wikitable" border="0" {| class="wikitable" border="0" Line 231: Line 236: |{{EquationRef|(23)}} |{{EquationRef|(23)}} |} |} - Integrating eq. (4.308) over the main control volume [see Fig. 4.23(a) for the two-dimensional projection of the three-dimensional control volume] and using a fully-implicit scheme, we obtain + Integrating eq. (23) over the main control volume [see (a) in above figure for the two-dimensional projection of the three-dimensional control volume] and using a fully-implicit scheme, we obtain - \begin{align} + [itex]\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} + \end{align}
- Substituting eqs. (4.305) – (4.307) into the above equation, the following algebraic equation for pressure correction is obtained: + 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.png|thumb|400 px|alt=Control volume for continuity equation at boundary |Figure 1: Control volume for continuity equation at boundary.]] + [[Image:Fig4.24.png|thumb|400 px|alt=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 $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. + 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(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. + + 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 (''ue'' 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 ''ue'' 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 Semi-Implicit Method for Pressure-Linked Equations. It was originally proposed by Patankar and Spalding (1972) and summarized in Patankar (1980). This algorithm is a semi-implicit method because the first terms on the right-hand 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 fully-implicit. The procedure for SIMPLE algorithm can be summarized as following:
+ The above approach for the solution of the incompressible flow field was named SIMPLE, which stands for Semi-Implicit Method for Pressure-Linked Equations. It was originally proposed by Patankar and Spalding and summarized in Patankar . This algorithm is a semi-implicit method because the first terms on the right-hand side of eqs. (15) – (17) are neglected. If these terms are retained, we will have to solve for the velocity corrections (${{{u}'}_{e}},\text{ }{{{v}'}_{n}},\text{ and }{{{w}'}_{t}}$) for the entire flow field simultaneously and the algorithm will become fully-implicit. 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).
+ 1.  Guess a pressure field, ''p''*.
- 3.  Solve the pressure correction equation (4.309) to get .
+ - 4.  Obtain a new pressure field from eq. (4.298).
+ 2.  Obtain the starred velocity field, (''u''*, ''v''*, and ''w''*) from eqs. (10) – (12).
- 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.
+ 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  $\phi$’s from the follow equation +
$\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}$ +
+ If the flow field is not affected by a particular $\phi$, 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.
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). + + 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 and Tao Tao 2001, W.Q., Numerical Heat Transfer, 2nd Ed., Xi’an Jiaotong University Press, Xi’an, China (in Chinese). , as well as the Handbook of Numerical Heat Transfer 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.. + + ==References== + {{Reflist}}

## Special Difficulties

For a two-dimensional incompressible flow problem without a body force, the continuity and momentum equations in the Cartesian coordinate system are $\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}=0$ (1) $\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}$ (2) $\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}$ (3)

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 $\varphi$ 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 $S_{u}=-\partial p/\partial x$ and $S_{v}=-\partial p/\partial y,$ 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 x-direction can be expressed as $\left( \frac{\partial p}{\partial x} \right)_{P}=\frac{p_{e}-p_{w}}{\Delta x}$ (4)

If central difference is employed, the pressures at the faces of the control volume become $p_{e}=\frac{p_{E}+p_{P}}{2},\text{ }p_{w}=\frac{p_{P}+p_{W}}{2}$

Substituting the above expressions into eq. (4), the pressure gradient in the x-direction becomes $\left( \frac{\partial p}{\partial x} \right)_{P}=\frac{p_{E}-p_{W}}{2\Delta x}$ (5)

Similarly, the pressure gradient in the y-direction can be expressed as $\left( \frac{\partial p}{\partial y} \right)_{P}=\frac{p_{N}-p_{S}}{2\Delta y}$ (6)

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 $\partial p/\partial x=0$ and $\partial p/\partial y=0$ 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 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

To overcome the checkerboard pressure field and develop an effective algorithm for the pressure field, one can use a staggered grid system   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 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 $a_{e}u_{e}=\sum{a_{nb}u_{nb}+b+A_{e}(p_{P}-p_{E})}$ (7)

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 two-dimensional problem, the area on which the pressure difference acts on is Ae = Δy. Similarly, the discretization equation for v can be obtained by integrating eq.(2) for the control volume of vn shown in (c) of the right figure, i.e., $a_{n}v_{n}=\sum{a_{nb}v_{nb}+b+A_{n}(p_{P}-p_{N})}$ (8)

where the area on which pressure difference acts on is An = Δx for a two-dimensional problem.

For a three-dimensional 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 right-hand side. The areas on which the pressure difference act on should be modified to Ae = ΔyΔz and An = ΔxΔz. In addition to eqs. (7) and (8), another equation for the velocity component in the z-direction, wt, is also needed, i.e., $a_{t}w_{t}=\sum{a_{nb}w_{nb}+b+A_{t}(p_{P}-p_{T})}$ (9)

which is obtained by integrating the momentum equation in the z-direction for the control volume of wt that are staggered to the positive z-direction by a half grid. The area on which the pressure difference, pPpT, acts on is At = Δ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). $a_{e}u_{e}^{*}=\sum{a_{nb}u_{nb}^{*}+b+A_{e}(p_{P}^{*}-p_{E}^{*})}$ (10) $a_{n}v_{n}^{*}=\sum{a_{nb}v_{nb}^{*}+b+A_{n}(p_{P}^{*}-p_{N}^{*})}$ (11) $a_{t}w_{t}^{*}=\sum{a_{nb}w_{nb}^{*}+b+A_{t}(p_{P}^{*}-p_{T}^{*})}$ (12)

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  can be used to obtain the correct pressure field and consequently, the velocity field. If the correct pressure can be obtained by $\begin{matrix}{}\\\end{matrix}p=p^{*}+{p}'$ (13)

where p' is pressure correction, the new velocity field based on the corrected pressure field will be $\begin{matrix}{}\\\end{matrix}u=u^{*}+{u}',\text{ }v=v^{*}+{v}',\text{ }w=w^{*}+{w}'$ (14)

where u', v', and w' are velocity corrections. The corrected velocity field satisfies eqs. (7) – (9), i.e., $a_{e}(u_{e}^{*}+{u}'_{e})=\sum{a_{nb}(u_{nb}^{*}+{u}'_{nb})+b+A_{e}[(p_{P}^{*}+{p}'_{P})-(p_{E}^{*}+{p}'_{E})}]$ $a_{n}(v_{n}^{*}+{v}'_{n})=\sum{a_{nb}(v_{nb}^{*}+{v}'_{nb})+b+A_{n}[(p_{P}^{*}+{p}'_{P})-(p_{N}^{*}+{p}'_{N})}]$ $a_{t}(w_{t}^{*}+{w}'_{t})=\sum{a_{nb}(w_{nb}^{*}+{w}'_{nb})+b+A_{t}[(p_{P}^{*}+{p}'_{P})-(p_{T}^{*}+{p}'_{T})}]$

Subtracting eqs. (10) – (12) from the above three equations yields the following equations for the velocity corrections: $a_{e}{u}'_{e}=\sum{a_{nb}{u}'_{nb}+A_{e}({p}'_{P}-{p}'_{E})}$ (15) $a_{n}{v}'_{n}=\sum{a_{nb}{v}'_{nb}+A_{n}({p}'_{P}-{p}'_{N})}$ (16) $a_{t}{w}'_{t}=\sum{a_{nb}{w}'_{nb}+A_{t}({p}'_{P}-{p}'_{T})}$ (17)

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 anb = 0 in eqs. (15) – (17), the velocity correction becomes $\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})$ (18)

where $\begin{matrix}{}\\\end{matrix}d_{e}=A_{e}/a_{e},\text{ }d_{n}=A_{n}/a_{n},\text{ }d_{t}=A_{t}/a_{t}$ (19)

Substituting eq. (18) into eq. (19), we have $u_{e}=u_{e}^{*}+d_{e}({p}'_{P}-{p}'_{E})$ (20) $v_{n}=v_{n}^{*}+d_{n}({p}'_{P}-{p}'_{N})$ (21) $w_{t}=w_{t}^{*}+d_{t}({p}'_{P}-{p}'_{T})$ (22)

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 three-dimensional flow, the continuity equation is $\frac{\partial \rho }{\partial t}+\frac{\partial (\rho u)}{\partial x}+\frac{\partial (\rho v)}{\partial y}+\frac{\partial (\rho w)}{\partial z}=0$ (23)

Integrating eq. (23) over the main control volume [see (a) in above figure for the two-dimensional projection of the three-dimensional control volume] and using a fully-implicit scheme, we obtain \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}

Substituting eqs. (20) – (22) into the above equation, the following algebraic equation for pressure correction is obtained: $\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$ (24)

where $\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$ (25) $\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$ (26) $\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$ (27) $\begin{matrix}{}\\\end{matrix}a_{P}=a_{E}+a_{W}+a_{N}+a_{S}+a_{T}+a_{B}$ (28) \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} (29)

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 (ue 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 aE = 0 [see eq. (25)] if ue 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 Semi-Implicit Method for Pressure-Linked Equations. It was originally proposed by Patankar and Spalding  and summarized in Patankar . This algorithm is a semi-implicit method because the first terms on the right-hand 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 fully-implicit. 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 $\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}$

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  and Tao , as well as the Handbook of Numerical Heat Transfer .

## References

1. 1.0 1.1 Patankar, S.V., and Spalding, D.B., 1972, “A Calculation Procedure for Heat, and Mass and Momentum Transfer in Three-Dimensional Parabolic Flows,” International Journal of Heat and Mass Transfer, Vol. 17, pp. 1787-1806.
2. 2.0 2.1 2.2 2.3 Patankar, S.V., 1980, Numerical Heat Transfer and Fluid Flow, Hemisphere, Washington, DC.
3. Tao 2001, W.Q., Numerical Heat Transfer, 2nd Ed., Xi’an Jiaotong University Press, Xi’an, China (in Chinese).
4. Minkowycz, M.J., Sparrow, E.M., and Murthy, J.Y., eds., 2006, Handbook of Numerical Heat Transfer, 2nd ed. John Wiley & Sons, Hoboken, NJ.
5. Faghri, A., Zhang, Y., and Howell, J. R., 2010, Advanced Heat and Mass Transfer, Global Digital Press, Columbia, MO.