Numerical Solution of Flow Field

From Thermal-FluidsPedia

**Related Topics Catalog**
Computational methodologies for forced convection One-Dimensional Steady-State Convection and Diffusion Central Difference Scheme Upwind Scheme Hybrid Scheme Exponential and Power Law Schemes A Generalized Expression of Discretization Schemes Multidimensional Convection and Diffusion Problems Numerical Solution of Flow Field Special Difficulties Staggered grid Pressure Correction Equation The SIMPLE Algorithm Numerical Simulation of Interfaces and Free Surfaces Application of Computational Methods

1 Special Difficulties
2 Staggered grid
3 Pressure Correction Equation
4 The SIMPLE Algorithm
5 References

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)

Control volume for one-dimensional problem.

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)

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. (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

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 ^[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 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 $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.,

$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 $A n = Δ 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 $A e = Δ y Δ z$ and $A n = Δ x Δ z$ . In addition to eqs. (7) and (8), another equation for the velocity component in the z-direction, w_t, 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 w_t 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 = Δ 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 ^[2] 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 $a n b = 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)

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. (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 Semi-Implicit Method for Pressure-Linked Equations. It was originally proposed by Patankar and Spalding ^[1] and summarized in Patankar ^[2]. 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 ^[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 Three-Dimensional Parabolic Flows,” International Journal of Heat and Mass Transfer, Vol. 17, pp. 1787-1806.
↑ ^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.

Retrieved from "http://www.thermalfluidscentral.org/encyclopedia/index.php/Numerical_Solution_of_Flow_Field"

Numerical Solution of Flow Field

From Thermal-FluidsPedia

Contents

Special Difficulties

Staggered grid

Pressure Correction Equation

The SIMPLE Algorithm

References

Views

Personal tools

Navigation

Search