**4. Numerical method**

#### **4.1. Transformation of coordinates**

Discretizationand integration of the transport equations described previously are performed using a non-orthogonal generalized mesh as shown in **Figure 1**. The independent Cartesian coordinates (*x, z*), describing the physical domain, are thus replaced by a boundary-fitted coordinate system (*ξ*, *ζ*), defined by the mesh lines which may have, locally, any orientation and inclination. The computational domain (as opposed to the physical domain) is the space considered in terms of the boundary-fitted coordinate system(*ξ, ζ*) as depicted in **Figure 1**. In the computational domain, the mesh spacing is considered, for convenience, as unitary, i.e., *Δξ* =*Δζ* =*1*, and the mesh lines are always horizontal (*ξ* lines) or vertical (*ζ* lines). The mesh arrangement is of the collocated type (as opposed to the staggered mesh) with the two velocity components and scalar quantities (temperature, pressure, turbulence kinetic energy and its dissipation rate or frequency, as well as concentrations) located at the control volume (CV) centre.

**Figure 1.** The two domains. (a) Physical domain, showing the mesh lines; (b) Computational domain.

Transformation of the original equations is accomplished by replacing the independent variables, using the chain rule, which states that, generically:

$$\frac{\partial \, \phi}{\partial \mathbf{x}} = \frac{\partial \phi}{\partial \xi} \frac{\partial \, \xi}{\partial \mathbf{x}} + \frac{\partial \phi}{\partial \zeta} \frac{\partial \zeta}{\partial \mathbf{x}} = \xi\_x \frac{\partial \phi}{\partial \xi} + \zeta\_x \frac{\partial \phi}{\partial \zeta} \tag{35}$$

The Jacobian of the transformation

The constants are *α*1= 5/9; *β*1= 3/40; *σ<sup>k</sup>*<sup>1</sup> = 0.85; *σω*<sup>1</sup> = 0.5; *α*2= 0.44; *β*<sup>2</sup> = 0.0828; *σ<sup>k</sup>*<sup>2</sup> = 1; *σω*2 = 0.856;

The near wall treatment for momentum and turbulence equations follows the proposal described in [5]. The basic principle behind the automatic wall functions is to switch from a low-Reynolds number formulation to a wall function based on the grid nodes proximity to the wall. According to these authors, the automatic wall treatment avoids the deterioration of

The known solutions for *ω* in the viscous (linear) and in the logarithmic near wall region are

t

c

*0 075 y 03 y* (31)

*1 vis* (32)

<sup>+</sup> <sup>+</sup> = = *vis U U 1 1 u u <sup>y</sup> E y* (33)

= + (34)

log ; . .

 w

For the turbulence kinetic energy, a zero flux along the direction perpendicular to the wall is assigned. In turn, for the momentum equations, a similar reasoning applies, with expressions

( )

ln

 t

Discretizationand integration of the transport equations described previously are performed using a non-orthogonal generalized mesh as shown in **Figure 1**. The independent Cartesian coordinates (*x, z*), describing the physical domain, are thus replaced by a boundary-fitted coordinate system (*ξ*, *ζ*), defined by the mesh lines which may have, locally, any orientation

with *U*1 being the fluid velocity relative to the wall velocity. The wall shear stress is computed

c

*6 u*

results typical of low-Reynolds models when applied on too coarse meshes.

n

w

for the shear velocity in the viscous and in the logarithmic region:

t

*vis* = = *<sup>2</sup>*

 ww = + log *2 2*

log ;

 t

( ) ( ) log *4 4 <sup>4</sup> vis uu u*

tt

w

270 Numerical Simulation - From Brain Imaging to Turbulent Flows

The imposed value for *ω* at the first node close to a wall is

*β\** = 0.09.

as follows:

**4. Numerical method**

**4.1. Transformation of coordinates**

$$J = \begin{vmatrix} \boldsymbol{\chi}\_{\mathcal{\xi}} & \boldsymbol{\chi}\_{\mathcal{\zeta}} \\ \boldsymbol{\varepsilon}\_{\mathcal{\xi}} & \boldsymbol{\varepsilon}\_{\mathcal{\zeta}} \end{vmatrix} = \boldsymbol{\chi}\_{\mathcal{\xi}} \boldsymbol{\varepsilon}\_{\mathcal{\zeta}} - \boldsymbol{\chi}\_{\mathcal{\zeta}} \boldsymbol{\varepsilon}\_{\mathcal{\xi}} \tag{36}$$

represents the ratio between the physical size of the control volume and its computational size (unitary). The derivatives *ξx*, *ξz*, *ζ<sup>x</sup> ζz* are the contravariant metrics of the transformation. They are computed from the covariant metrics, *xξ*, *zξ*, *xζ*, *zζ* as follows:

$$\xi\_{\mathbf{x}} = \boldsymbol{\varpi}\_{\boldsymbol{\zeta}} \;/\ J ; \quad \xi\_{\boldsymbol{\varpi}} = -\mathbf{x}\_{\boldsymbol{\zeta}} \;/\ J ; \quad \boldsymbol{\zeta}\_{\mathbf{x}} = -\boldsymbol{\varpi}\_{\boldsymbol{\xi}} \;/\ J ; \quad \boldsymbol{\zeta}\_{\boldsymbol{\varpi}} = \mathbf{x}\_{\boldsymbol{\xi}} \;/\ J \tag{37}$$

To obtain the strong conservative form in the boundary-fitted coordinate system (*ξ, ζ*), the transport equations are transformed through the application of the chain rule (35) and multiplied by the Jacobian of the transformation. The metric identity

$$\frac{\partial}{\partial \xi} (J\xi\_x) + \frac{\partial}{\partial \zeta} (J\zeta\_x) = 0; \quad \frac{\partial}{\partial \xi} (J\xi\_z) + \frac{\partial}{\partial \zeta} (J\zeta\_z) = 0 \tag{38}$$

is then used to recast some terms. The result, for a generic variable *ϕ*, is

$$\begin{split} \int \frac{\hat{\mathcal{I}}\left(\rho\phi\right)}{\hat{\mathcal{O}}t} + \frac{\hat{\mathcal{O}}}{\hat{\mathcal{O}}\tilde{\xi}} \big(J\rho U\phi\big) + \frac{\hat{\mathcal{O}}}{\hat{\mathcal{O}}\tilde{\xi}} \big(J\rho W\phi\big) &=\\ \frac{\partial}{\partial\tilde{\xi}} \Big[ \Gamma J \text{g}^{IJ} \frac{\partial\tilde{\phi}}{\partial\tilde{\xi}}\Big] + \frac{\hat{\mathcal{O}}}{\partial\tilde{\xi}} \Big[ \Gamma J \text{g}^{13} \frac{\partial\tilde{\phi}}{\partial\tilde{\xi}}\Big] + \frac{\hat{\mathcal{O}}}{\partial\tilde{\xi}} \Big[ \Gamma J \text{g}^{13} \frac{\partial\tilde{\phi}}{\partial\tilde{\xi}}\Big] + \frac{\hat{\mathcal{O}}}{\partial\tilde{\xi}} \Big[ \Gamma J \text{g}^{13} \frac{\partial\tilde{\phi}}{\partial\tilde{\xi}}\Big] + J S\_{\Phi} \end{split} \tag{39}$$

The terms *g11*, *g13* and *g33* are the contravariant metric relations given by

$$\mathbf{g}^{11} = \boldsymbol{\xi}\_{\mathbf{x}}^{2} + \boldsymbol{\xi}\_{z}^{2}; \quad \mathbf{g}^{33} = \boldsymbol{\zeta}\_{\mathbf{x}}^{2} + \boldsymbol{\zeta}\_{z}^{2}; \quad \mathbf{g}^{13} = \boldsymbol{\xi}\_{\mathbf{x}}^{} \boldsymbol{\zeta}\_{\mathbf{x}} + \boldsymbol{\xi}\_{z} \boldsymbol{\zeta}\_{z} \tag{40}$$

The non-orthogonal term *g13* is null if the mesh is locally orthogonal. *U* and *W*, in Eq. (39), are the contravariant velocities. The terms *JρU* and *JρW* represent mass fluxes through the control volume faces along the computational directions *ξ* and *ζ*, respectively, and are computed as follows:

$$F\_{\tilde{\xi}} = J \, \rho U = \rho \left( \mathbf{z}\_{\zeta} \, \mu - \mathbf{x}\_{\zeta} \, \mathbf{w} \right) \tag{41}$$

$$F\_{\zeta'} = J \, \rho W = \rho \left(\mathbf{x}\_{\xi} \, \mathbf{w} - \mathbf{z}\_{\xi} \boldsymbol{\mu}\right) \tag{42}$$

Note that, in this case, the sub-index for the fluxes (such as in *Fξ* ) represents the flux direction (and not a derivative, such as for the metrics). Eq. (39) may be rewritten as

$$J\frac{\partial\left(\rho\phi\right)}{\partial t} + \frac{\partial}{\partial\xi}\left(F\_{\xi}\phi\right) + \frac{\partial}{\partial\zeta}\left(F\_{\zeta}\phi\right) = \frac{\partial}{\partial\xi}\left[\Gamma J \mathbf{g}^{IJ} \frac{\partial\phi}{\partial\xi}\right] + \frac{\partial}{\partial\zeta}\left[\Gamma J \mathbf{g}^{13} \frac{\partial\phi}{\partial\zeta}\right] + J\left(S\_{cross} + S\_{\phi}\right) \tag{43}$$

where the cross-derivatives were incorporated into the source term *Scross*. A similar procedure is applied to obtain the generalized form of remaining equations, leading to the following form for the Navier-Stokes and continuity equations:

$$\begin{aligned} \begin{aligned} \displaystyle J\frac{\overrightarrow{\mathcal{C}}(\rho u)}{\partial t} + \frac{\overrightarrow{\mathcal{C}}}{\partial \overrightarrow{\xi}}(F\_{\xi}u) + \frac{\overrightarrow{\mathcal{C}}}{\partial \overrightarrow{\zeta}}(F\_{\zeta}u) &= \\ \displaystyle \frac{\overrightarrow{\mathcal{C}}}{\partial \overleftarrow{\xi}} \bigg[ \Gamma J \text{g}^{1/} \frac{\partial u}{\partial \xi} \bigg] + \frac{\overrightarrow{\mathcal{C}}}{\partial \overleftarrow{\xi}} \bigg[ \Gamma J \text{g}^{13} \frac{\partial u}{\partial \zeta} \bigg] - z\_{\zeta} \frac{\overrightarrow{\mathcal{C}} p}{\overrightarrow{\zeta} \overline{\xi}} + z\_{\xi} \frac{\overrightarrow{\mathcal{C}} p}{\overrightarrow{\mathcal{C}} \zeta} + S\_{cross} \\\\ \displaystyle J\frac{\overrightarrow{\mathcal{C}}(\rho w)}{\partial t} + \frac{\overrightarrow{\mathcal{C}}}{\partial \overleftarrow{\xi}}(F\_{\xi}w) + \frac{\overrightarrow{\mathcal{C}}}{\partial \overleftarrow{\zeta}} \bigg(F\_{\xi}w \bigg) = \frac{\overrightarrow{\mathcal{C}}}{\overrightarrow{\mathcal{C}} \overline{\xi}} \bigg[ \Gamma J \text{g}^{11} \frac{\partial w}{\partial \overline{\xi}} \bigg] + \frac{\overrightarrow{\mathcal{C}}}{\partial \overleftarrow{\zeta}} \bigg[ \Gamma J \text{g}^{13} \frac{\partial w}{\partial \overline{\zeta}} \bigg] \\\\ \square \end{aligned} \tag{45}$$

$$\mathbf{x} + \mathbf{x}\_{\xi} \frac{\partial \mathbf{p}}{\partial \xi} - \mathbf{x}\_{\xi} \frac{\partial \mathbf{p}}{\partial \xi} + \mathbf{S}\_{crossw} + I \tag{43}$$

$$J\frac{\partial \boldsymbol{\hat{\rho}}}{\partial t} + \frac{\partial}{\partial \boldsymbol{\xi}} \Big[\boldsymbol{\rho}\left(\boldsymbol{z}\_{\boldsymbol{\xi}}\boldsymbol{\mu} - \mathbf{x}\_{\boldsymbol{\xi}}\boldsymbol{w}\right)\Big] + \frac{\partial}{\partial \boldsymbol{\xi}} \Big[\boldsymbol{\rho}\left(\mathbf{x}\_{\boldsymbol{\xi}}\boldsymbol{w} - \boldsymbol{z}\_{\boldsymbol{\xi}}\boldsymbol{u}\right)\Big] = J\frac{\partial \boldsymbol{\hat{\rho}}}{\partial \boldsymbol{t}} + \frac{\partial}{\partial \boldsymbol{\xi}} \Big(J\boldsymbol{\rho}U\Big) + \frac{\partial}{\partial \boldsymbol{\xi}} \Big(J\boldsymbol{\rho}W\Big) = 0\tag{46}$$

#### **4.2. Integration**

represents the ratio between the physical size of the control volume and its computational size (unitary). The derivatives *ξx*, *ξz*, *ζ<sup>x</sup> ζz* are the contravariant metrics of the transformation. They

/; /; /; /

To obtain the strong conservative form in the boundary-fitted coordinate system (*ξ, ζ*), the transport equations are transformed through the application of the chain rule (35) and

 z

( *J J0 J J0*

*11 33 13 13*

=+ =+ = +

*F J U zu xw*

 r

*F J W xw zu*

= = -

 r

r

r

(and not a derivative, such as for the metrics). Eq. (39) may be rewritten as

= = -

 zz

 xx

Note that, in this case, the sub-index for the fluxes (such as in *Fξ* ) represents the flux direction

x

z

The non-orthogonal term *g13* is null if the mesh is locally orthogonal. *U* and *W*, in Eq. (39), are the contravariant velocities. The terms *JρU* and *JρW* represent mass fluxes through the control volume faces along the computational directions *ξ* and *ζ*, respectively, and are computed as

 z z

 ¶

é ùé ùé ùé ù + +++ ê úê úê úê ú ë ûë ûë ûë û

 ¶x

*Jg Jg Jg Jg JS*

*xx zz* ) ( ) ; ( ) ( )

+= +=

xz

 ¶f

 ¶z  ¶

 ¶z

 xz

; ; *11 2 2 33 2 2 13 x z x z xx zz ggg* (40)

 G

 xz

xz

¶¶ ¶¶ (38)

xx

 z*xz x z* = =- =- = *zJ xJ zJ xJ* (37)

f

(39)

 ¶f

 ¶x

( ) (41)

( ) (42)

are computed from the covariant metrics, *xξ*, *zξ*, *xζ*, *zζ* as follows:

multiplied by the Jacobian of the transformation. The metric identity

¶¶ ¶¶

is then used to recast some terms. The result, for a generic variable *ϕ*, is

 r f

 ¶f

 ¶z

The terms *g11*, *g13* and *g33* are the contravariant metric relations given by

xz

xz

( ) ( ) ( )

rf

 ¶

x x  ¶

 ¶z

++ =

GGG

*J JU JW*

¶

 ¶x

 ¶f

 ¶x ¶z

*t*

¶

¶ rf

¶

¶x

follows:

 z

z

272 Numerical Simulation - From Brain Imaging to Turbulent Flows

xx

The integration and solution method for the transport equations are entirely based on the methodology in [6], with some of the suggestions described in [7] and incorporating the necessary modifications for the generalized mesh approach. The general Eq. (43) may be written as

$$J\frac{\partial\left(\rho\phi\right)}{\partial t} + \frac{\partial}{\partial\tilde{\xi}}\bigg(F\_{\xi}\phi - \Gamma J \text{g}^{II}\frac{\partial\phi}{\partial\tilde{\xi}}\bigg) + \frac{\partial}{\partial\tilde{\xi}}\bigg(F\_{\xi}\phi - \Gamma J \text{g}^{II}\frac{\partial\phi}{\partial\tilde{\xi}}\bigg) = J\left(S\_{\text{cross}} + S\_{\phi}\right) \tag{47}$$

The integration of the previous equation in its control volume leads to

$$\begin{aligned} \frac{\partial \left(\rho \phi\right)}{\partial t} + \left[F\_s \phi\_e - D\_e \left(\phi\_E - \phi\_P\right)\right] - \left[F\_o \phi\_o - D\_o \left(\phi\_P - \phi\_o\right)\right] + \\\ \left[F\_i \phi\_i - D\_i \left(\phi\_r - \phi\_P\right)\right] - \left[F\_b \phi\_b - D\_b \left(\phi\_P - \phi\_B\right)\right] = J\left(S\_{1\_p} + S\_{2\_p} \phi\_P\right) \end{aligned} \tag{48}$$

where the source term has been written as a linear combination involving *ϕP. F* and *D* represent the advective fluxes and the diffusive coefficients, respectively:

$$\begin{aligned} F\_e &= F\_{\tilde{\xi}e} = \left[ \rho \left( z\_{\tilde{\xi}} \boldsymbol{\mu} - \mathbf{x}\_{\tilde{\xi}} \boldsymbol{\mathcal{W}} \right) \right]\_{\mathcal{C}}; \; F\_o = F\_{\tilde{\xi}o} = \left[ \rho \left( z\_{\tilde{\xi}} \boldsymbol{\mu} - \mathbf{x}\_{\tilde{\xi}} \boldsymbol{\mathcal{W}} \right) \right]\_{\mathcal{O}} \\\\ F\_t &= F\_{\tilde{\xi}t} = \left[ \rho \left( \mathbf{x}\_{\tilde{\xi}} \boldsymbol{\mu} - z\_{\tilde{\xi}} \boldsymbol{\mu} \right) \right]\_t; \; \; F\_b = F\_{\tilde{\xi}b} = \left[ \rho \left( \mathbf{x}\_{\tilde{\xi}} \boldsymbol{\mu} - z\_{\tilde{\xi}} \boldsymbol{\mu} \right) \right]\_b \end{aligned} \tag{49}$$

$$D\_{\mathcal{E}} = \left(\Gamma J \mathbf{g}^{II}\right)\_{\mathcal{E}}; \quad D\_o = \left(\Gamma J \mathbf{g}^{II}\right)\_o; \quad D\_t = \left(\Gamma J \mathbf{g}^{I\mathcal{I}}\right)\_I; \quad D\_b = \left(\Gamma J \mathbf{g}^{I\mathcal{I}}\right)\_b \tag{50}$$

In the previous expressions, the subscripts indicate the location relative to the *CV* centre, in the computational domain, with the uppercase denoting neighbour nodes and the lowercase denoting neighbour faces: *E,e*: East; *O,o*: West; *T,t*: Top; *B,b*: Bottom. For simplicity, in Eq. (48), the subscripts *ξ* and *ζ* were dropped from the fluxes *F* (in fact, fluxes across the eastern and western faces are always along the *ξ* direction and fluxes across the top and bottom faces are always along the *ζ* direction).

For the solution of the equations, it is necessary to evaluate the values of *ϕ* in the *CV* faces (i.e., *ϕe*, *ϕo*, *ϕt*, *ϕb*). These values are computed as a function of both *ϕP* and the values in the neighbour nodes *ϕE*, *ϕO*, *ϕ<sup>T</sup>* , *ϕ<sup>B</sup>* according to the adopted advection scheme (to be descri‐ bed later on). Eq. (48) may, then, be written in the following standard form:

$$a\_P \phi\_P = a\_E \phi\_E + a\_O \phi\_O + a\_T \phi\_T + a\_B \phi\_B + b \tag{51}$$

or, in a more compact manner,

$$a\_p \phi\_p = \sum\_{nb} a\_{nb} \mu\_{nb} + b \tag{51a}$$

with "*nb*" indicating that the sum is to be performed for all the neighbouring locations.

It is necessary to compute the face values *ϕe*, *ϕo*, *ϕt*, *ϕb* as a function of *ϕP*, *ϕE*, *ϕO*, *ϕ<sup>T</sup>* , *ϕ<sup>B</sup>* in order to obtain the coefficients *aE*, *aO*, *aT* , *aB* of Eq. (51). It is known that a simple arithmetic average is not a physically plausible solution since, due to the presence of a flow, the property *ϕ*, being advected, tends to assume a value closer to the upwind value. Several advection schemes may be adopted, being the simplest one the *upwind* scheme. According to this scheme, the property *ϕ* in the CVface takes the upwind value, i.e.,for example, *Fe* >*0* ⇒*ϕ<sup>e</sup>* =*ϕP*; *Fe* <*0* ⇒*ϕ<sup>e</sup>* =*ϕE* . The complete mathematical formulation for the upwind scheme is

$$\begin{aligned} a\_E &= D\_e + \left[ -F\_e, 0 \right] \; ; \; a\_O = D\_o + \left[ F\_o, 0 \right] \\ a\_T &= D\_t + \left[ -F\_t, 0 \right] \; ; \; a\_B = D\_b + \left[ F\_b, 0 \right] \end{aligned} \tag{52}$$

Due to its first-order character, the upwind scheme is often not used due to associated numerical diffusion. The Quick scheme is third-order accurate. The deferred correction version of Hayase [8] combines the first-order upwind scheme with a third-order correction, *b*quick, added to the source term as follows:

$$\begin{aligned} \label{eq:1} \delta\_{quark} &= \\ & \frac{I}{8} \| \boldsymbol{F}\_{e}, \boldsymbol{0} \| \left( -\phi\_{O} - 2\phi\_{P} + 3\phi\_{E} \right) - \frac{I}{8} \| -\boldsymbol{F}\_{e}, \boldsymbol{0} \| \left( 3\phi\_{P} - 2\phi\_{E} - 3\phi\_{EE} \right) + \\ & \frac{I}{8} \| \boldsymbol{F}\_{o}, \boldsymbol{0} \| \left( 3\phi\_{O} - 2\phi\_{P} - \phi\_{E} \right) - \frac{I}{8} \| -\boldsymbol{F}\_{o}, \boldsymbol{0} \| \left( -\phi\_{OO} - 2\phi\_{O} + 3\phi\_{P} \right) + \\ & \frac{I}{8} \| \boldsymbol{F}\_{t}, \boldsymbol{0} \| \left( -\phi\_{B} - 2\phi\_{P} + 3\phi\_{T} \right) - \frac{I}{8} \| -\boldsymbol{F}\_{t}, \boldsymbol{0} \| \left( 3\phi\_{P} - 2\phi\_{T} - 3\phi\_{TT} \right) + \\ & \frac{I}{8} \| \boldsymbol{F}\_{b}, \boldsymbol{0} \| \left( 3\phi\_{B} - 2\phi\_{P} - \phi\_{T} \right) - \frac{I}{8} \| -\boldsymbol{F}\_{b}, \boldsymbol{0} \| \left( -\phi\_{BB} - 2\phi\_{B} + 3\phi\_{P} \right) + \end{aligned} \tag{53}$$

The Quick scheme, although third-order accurate, presents oscillations that may lead to some unrealistic behaviour. Total variation diminishing (TVD) schemes, also implemented in the present code, were developed to provide second-order accurate solutions that are free or nearlyfree from oscillations. For further information on this, please refer to, e.g., [9].

#### **4.3. Pressure-velocity coupling**

where the source term has been written as a linear combination involving *ϕP. F* and *D* represent

*e e* ( ) ; *o o* ( ) *e o F F zu xw F F zu xw* == - == -

éù éù

 ) ; ;; ( ) ( ) ( ) *<sup>11</sup> <sup>11</sup> <sup>33</sup> <sup>33</sup> eotb eot b D Jg D Jg D Jg D Jg* (50)

In the previous expressions, the subscripts indicate the location relative to the *CV* centre, in the computational domain, with the uppercase denoting neighbour nodes and the lowercase denoting neighbour faces: *E,e*: East; *O,o*: West; *T,t*: Top; *B,b*: Bottom. For simplicity, in Eq. (48), the subscripts *ξ* and *ζ* were dropped from the fluxes *F* (in fact, fluxes across the eastern and western faces are always along the *ξ* direction and fluxes across the top and bottom faces

For the solution of the equations, it is necessary to evaluate the values of *ϕ* in the *CV* faces (i.e., *ϕe*, *ϕo*, *ϕt*, *ϕb*). These values are computed as a function of both *ϕP* and the values in the neighbour nodes *ϕE*, *ϕO*, *ϕ<sup>T</sup>* , *ϕ<sup>B</sup>* according to the adopted advection scheme (to be descri‐

*t t* ( ) ; *b b* ( ) *t b F F xw zu F F xw zu* == - = = -

==== (GGGG

bed later on). Eq. (48) may, then, be written in the following standard form:

f

*a aa aab PP EE OO TT BB*

*P P nb nb nb a au b*

*Fe* <*0* ⇒*ϕ<sup>e</sup>* =*ϕE* . The complete mathematical formulation for the upwind scheme is

with "*nb*" indicating that the sum is to be performed for all the neighbouring locations.

It is necessary to compute the face values *ϕe*, *ϕo*, *ϕt*, *ϕb* as a function of *ϕP*, *ϕE*, *ϕO*, *ϕ<sup>T</sup>* , *ϕ<sup>B</sup>* in order to obtain the coefficients *aE*, *aO*, *aT* , *aB* of Eq. (51). It is known that a simple arithmetic average is not a physically plausible solution since, due to the presence of a flow, the property *ϕ*, being advected, tends to assume a value closer to the upwind value. Several advection schemes may be adopted, being the simplest one the *upwind* scheme. According to this scheme, the property *ϕ* in the CVface takes the upwind value, i.e.,for example, *Fe* >*0* ⇒*ϕ<sup>e</sup>* =*ϕP*;

f

ffff

 x

éùéù

ëûëû

 z  zz

 r

 xx

ëû ëû (49)

= + +++ (51)

=å + (51a)

 r

the advective fluxes and the diffusive coefficients, respectively:

r

 xx

r

 zz

x

z

274 Numerical Simulation - From Brain Imaging to Turbulent Flows

are always along the *ζ* direction).

or, in a more compact manner,

EasyCFD adopts the SIMPLEC algorithm (Semi-Implicit Method for Pressure-Linked Equa‐ tions-Consistent) [7], which is based on the original formulation *SIMPLE* [6]. Due to the nonstaggered mesh arrangement (collocated mesh), the Rie-Chow interpolation procedure [10], with the modifications proposed in [11] and [12], is implemented. Let us consider the *u* momentum conservation Eq. (44). After integration, the evaluation of this equation leads to

$$a\_p \mu\_p = \sum\_{nb} a\_{nb} \mu\_{nb} - z\_{\zeta} \frac{\partial \lrcorner p}{\partial \zeta} + z\_{\zeta} \frac{\partial \lrcorner p}{\partial \zeta} + S\_{cross} \tag{54}$$

During the iterative process, velocities *u*\* are computed from the available velocity field *um* and pressure field *p*\* obtained at the previous step as follows:

$$a\_P \left( I + \frac{I}{E} \right) \boldsymbol{\mu}\_p^\* = \sum\_{nb} a\_{nb} \boldsymbol{\mu}\_{nb}^\* + \frac{a\_P}{E} \boldsymbol{\mu}\_p^m - \boldsymbol{z}\_\zeta \frac{\boldsymbol{\mathcal{O}} \boldsymbol{p}^\*}{\boldsymbol{\mathcal{O}} \boldsymbol{\xi}} + \boldsymbol{z}\_\xi \frac{\boldsymbol{\mathcal{O}} \boldsymbol{p}^\*}{\boldsymbol{\mathcal{O}} \boldsymbol{\zeta}} + S\_{cross} \tag{55}$$

where *E* is the under-relaxation factor [7]. During the iterative process, unless convergence is reached, the velocity field *u\** obtained from the solution of the momentum equations does not satisfy the continuity requirement unless the correct pressure field *p* is employed (instead of an incorrect pressure field *p*\* ). This means that, if the correct pressure field was employed, a mass-conservative velocity field would be obtained:

$$a\_P \left( I + \frac{I}{E} \right) \mu\_P = \sum\_{nb} a\_{nb} \mu\_{nb} + \frac{a\_P}{E} u\_p^m - z\_{\zeta} \frac{\partial \, p}{\partial \xi} + z\_{\xi} \frac{\partial \, p}{\partial \zeta} + b \tag{56}$$

Thus, the pressure and velocity fields *p*\* and *u*\* should be corrected by a certain amount *p*' and *u*' :

$$p = p^\* + p'; \quad u\_P = u\_P^\* + u\_P' \tag{57}$$

Subtracting Eq. (56) from Eq. (55) and taking into account Eq. (57), one obtains

$$a\_p \left( I + \frac{I}{E} \right) \stackrel{\cdot}{\mu\_p} = \sum\_{nb} a\_{nb} \stackrel{\cdot}{\mu\_{nb}} - z\_{\xi} \frac{\stackrel{\cdot}{\sigma} \stackrel{\cdot}{p}}{\partial \xi} + z\_{\xi} \frac{\stackrel{\cdot}{\sigma} \stackrel{\cdot}{p}}{\partial \xi} \tag{58}$$

The keystone of the SIMPLECalgorithm consists on the subtraction of the term ∑ *nb anbuP* ' on both sides of the equation and subsequent dropping, leading to

$$
\left[a\_P \left(I + \frac{I}{E}\right) - \sum\_{nb} a\_{nb}\right] \stackrel{\cdot}{u}\_P = \sum\_{pb} a\_{nb} \left(\mu\_{b\mathfrak{F}}^{\cdot} - u\_P\right) - z\_{\zeta} \frac{\partial \not\!p}{\partial \xi} + z\_{\xi} \frac{\partial \not\!p}{\partial \zeta} \tag{59}
$$

or, for simplicity:

$$
\tilde{a}\_P \stackrel{\circ}{u}\_P = -\mathbb{Z}\_{\tilde{\zeta}} \frac{\stackrel{\circ}{\mathcal{P}} \stackrel{\circ}{p}}{\stackrel{\circ}{\mathcal{P}}} + \mathbb{Z}\_{\tilde{\xi}} \frac{\stackrel{\circ}{\mathcal{P}} \stackrel{\circ}{p}}{\stackrel{\circ}{\mathcal{Q}}} \tag{60}
$$

where

$$
\tilde{a}\_P = a\_P \left( I + \frac{I}{E} \right) - \sum\_{nb} a\_{nb} \tag{61}
$$

The equations for the velocity correction are obtained through the pressure correction field, recurring to the previous equation:

$$\dot{\tilde{u}}\_{P} = -\frac{\boldsymbol{z}\_{\zeta}}{\tilde{\tilde{a}}\_{P}} \frac{\partial \dot{\tilde{p}}}{\partial \tilde{\xi}} + \frac{\boldsymbol{z}\_{\tilde{\xi}}}{\tilde{\tilde{a}}\_{P}} \frac{\partial \dot{\tilde{p}}}{\partial \tilde{\xi}}; \quad \dot{\tilde{w}}\_{P} = -\frac{\boldsymbol{x}\_{\tilde{\xi}}}{\tilde{\tilde{a}}\_{P}} \frac{\partial \dot{\tilde{p}}}{\partial \tilde{\xi}} + \frac{\boldsymbol{x}\_{\zeta}}{\tilde{\tilde{a}}\_{P}} \frac{\partial \dot{\tilde{p}}}{\partial \tilde{\xi}}\tag{62}$$

leading to

where *E* is the under-relaxation factor [7]. During the iterative process, unless convergence is reached, the velocity field *u\** obtained from the solution of the momentum equations does not satisfy the continuity requirement unless the correct pressure field *p* is employed (instead of

z

ç ÷ += +- + + <sup>å</sup> è ø *<sup>P</sup>*

*1 a pp a 1 u au u z z b*

'

Subtracting Eq. (56) from Eq. (55) and taking into account Eq. (57), one obtains

' '

The keystone of the SIMPLECalgorithm consists on the subtraction of the term ∑

æ ö + =å - + ç ÷

*<sup>1</sup> p p a 1 u au z z*

z

' '

 ¶

> ¶z

 x

*nb*

¶

¶x

*<sup>P</sup> P nb nb nb*

( ) ' '' é ù æ ö ê ú ç ÷ +- = - å å ë û è ø *<sup>P</sup> nb P nb nb P nb nb <sup>1</sup> a1 a u a u u*

> '

> > z

¶

æ ö = +- ç ÷ <sup>å</sup> è ø %*P P nb*

*<sup>1</sup> a a1 a*

¶x

*E*

sides of the equation and subsequent dropping, leading to

*E*

or, for simplicity:

where

¶

¶x

). This means that, if the correct pressure field was employed, a

 x

*E E* (56)

 ¶

> ¶z

\* \* =+ =+ ¢ *PPP p p p; u u u* (57)

' '

 ¶

> ¶z

z

¶

¶x

%*P P* =- + *p p au z z* (60)

è ø (58)

 x

and *u*\* should be corrected by a certain amount *p*' and

' '

 ¶


 ¶z

*<sup>E</sup>* (61)

 x *nb anbuP* ' on both

an incorrect pressure field *p*\*

*u*' :

mass-conservative velocity field would be obtained:

*<sup>P</sup> <sup>m</sup> <sup>P</sup> P nb nb nb*

æ ö

276 Numerical Simulation - From Brain Imaging to Turbulent Flows

Thus, the pressure and velocity fields *p*\*

$$\boldsymbol{u}\_{P} = \stackrel{\bullet}{\boldsymbol{u}\_{P}} - \frac{\boldsymbol{\pi}\_{\boldsymbol{\zeta}}}{\tilde{\boldsymbol{a}}\_{P}} \frac{\boldsymbol{\mathcal{O}} \stackrel{\bullet}{\boldsymbol{p}}}{\tilde{\boldsymbol{\alpha}}\_{P}^{\sf F}} + \frac{\boldsymbol{\pi}\_{\boldsymbol{\xi}}}{\tilde{\boldsymbol{a}}\_{P}} \frac{\boldsymbol{\mathcal{O}} \stackrel{\bullet}{\boldsymbol{p}}}{\tilde{\boldsymbol{\alpha}}\_{P}^{\sf F}} \; ; \; \boldsymbol{w}\_{P} = \stackrel{\bullet}{\boldsymbol{w}\_{P}} + \frac{\boldsymbol{\chi}\_{\boldsymbol{\zeta}}}{\tilde{\boldsymbol{a}}\_{P}} \frac{\boldsymbol{\mathcal{O}} \stackrel{\bullet}{\boldsymbol{p}}}{\tilde{\boldsymbol{a}}\_{P}^{\sf F}} - \frac{\boldsymbol{\chi}\_{\boldsymbol{\xi}}}{\tilde{\boldsymbol{a}}\_{P}} \frac{\boldsymbol{\mathcal{O}} \stackrel{\bullet}{\boldsymbol{p}}}{\tilde{\boldsymbol{a}}\_{P}^{\sf F}} \tag{63}$$

#### **4.4. The pressure correction equation**

As previously stated, the objective of the pressure correction is to produce a pressure field such that the solution of the momentum equations is a mass-conservative velocity field. Conse‐ quently, the equations for solving the *p'* field must be obtained from the continuity equation. The discretized form of this equation is obtained directly from the integration of Eq. (2):

$$\frac{1}{2}J\frac{\rho\_p - \rho\_p^0}{\Delta t} + \left[\rho\left(z\_\xi \boldsymbol{\mu} - \mathbf{x}\_\xi \boldsymbol{\mu}\right)\right]\_\mathbf{e} - \left[\rho\left(z\_\xi \boldsymbol{\mu} - \mathbf{x}\_\xi \boldsymbol{\mu}\right)\right]\_\mathbf{e} + \left[\rho\left(\mathbf{x}\_\xi \boldsymbol{\mu} - z\_\xi \boldsymbol{\mu}\right)\right]\_\mathbf{e} - \left[\rho\left(\mathbf{x}\_\xi \boldsymbol{\mu} - z\_\xi \boldsymbol{\mu}\right)\right]\_\mathbf{b} = 0\tag{64}$$

As one may see, velocities are, now, needed at the control volume faces. Taking Eq. (63), for the *u* velocity component, at the "*e* " and "*o*" faces and the same equation for the *w* component at the "*t*" and "*b*" faces, substituting into Eq. (64) and rearranging the terms leads to

$$\begin{split} \frac{1}{\Delta t} \frac{\rho\_{P} - \rho\_{P}^{0}}{\Delta t} + F\_{e}^{\*} &- F\_{o}^{\*} - F\_{i}^{\*} - \left( \frac{\rho \mathbf{g}\_{13}}{\tilde{\mathbf{a}}\_{P}} \frac{\mathcal{C} \mathbf{p}^{\cdot}}{\partial \tilde{\mathbf{y}}\_{\sigma}^{\prime}} \right)\_{e} + \left( \frac{\rho \mathbf{g}\_{13}}{\tilde{\mathbf{a}}\_{P}} \frac{\mathcal{C} \mathbf{p}^{\cdot}}{\partial \tilde{\mathbf{y}}\_{\sigma}^{\prime}} \right)\_{o} - \left( \frac{\rho \mathbf{g}\_{11}}{\tilde{\mathbf{a}}\_{P}} \frac{\mathcal{C} \mathbf{p}^{\cdot}}{\partial \tilde{\mathbf{y}}\_{\sigma}^{\prime}} \right)\_{i} + \left( \frac{\rho \mathbf{g}\_{11}}{\tilde{\mathbf{a}}\_{P}} \frac{\mathcal{C} \mathbf{p}^{\cdot}}{\partial \tilde{\mathbf{y}}\_{\sigma}^{\prime}} \right)\_{b} + \\\left( \frac{\rho \mathbf{g}\_{13}}{\tilde{\mathbf{a}}\_{P}} \frac{\mathcal{C} \mathbf{p}^{\cdot}}{\partial \tilde{\mathbf{y}}\_{\sigma}^{\prime}} \right)\_{o} - \left( \frac{\rho \mathbf{g}\_{13}}{\tilde{\mathbf{a}}\_{P}} \frac{\mathcal{C} \mathbf{p}^{\cdot}}{\partial \tilde{\mathbf{y}}\_{\sigma}^{\prime}} \right)\_{i} - \left( \frac{\rho \mathbf{g}\_{13}}{\tilde{\mathbf{a}}\_{P}} \frac{\mathcal{C} \mathbf{p}^{\cdot}}{\partial \tilde{\mathbf{y}}\_{\sigma}^{\prime}} \right)\_{b} = 0 \end{split} \tag{65}$$

The terms *gij* are the covariant metric relations defined as

$$\mathbf{g}\_{11} = \mathbf{x}\_{\boldsymbol{\xi}}^{2} + \mathbf{z}\_{\boldsymbol{\xi}}^{2} = \mathbf{g}^{33}J^{2}; \ \mathbf{g}\_{13} = \mathbf{x}\_{\boldsymbol{\xi}}^{2} + \mathbf{z}\_{\boldsymbol{\xi}}^{2} = \mathbf{g}^{ll}J^{2}; \ \ \mathbf{g}\_{13} = \mathbf{x}\_{\boldsymbol{\xi}}\mathbf{x}\_{\boldsymbol{\xi}} + \mathbf{z}\_{\boldsymbol{\xi}}\mathbf{z}\_{\boldsymbol{\xi}} = -\mathbf{g}^{ll}J^{2} \tag{66}$$

The derivatives are evaluated as

$$
\left(\frac{\partial \dot{\boldsymbol{p'}}}{\partial \dot{\boldsymbol{\xi}}}\right)\_{\boldsymbol{\epsilon}} = \dot{\boldsymbol{p}\_{\boldsymbol{\varepsilon}}} - \dot{\boldsymbol{p}\_{\boldsymbol{p}}} \cdot \left(\frac{\partial \dot{\boldsymbol{p'}}}{\partial \dot{\boldsymbol{\xi}}}\right)\_{\boldsymbol{o}} = \dot{\boldsymbol{p}\_{\boldsymbol{p}}} - \dot{\boldsymbol{p}\_{o}} \cdot \left(\frac{\partial \dot{\boldsymbol{p'}}}{\partial \boldsymbol{\xi}}\right)\_{\boldsymbol{\epsilon}} = \dot{\boldsymbol{p}\_{\boldsymbol{r}}} - \dot{\boldsymbol{p}\_{\boldsymbol{p}}} \cdot \left(\frac{\partial \dot{\boldsymbol{p'}}}{\partial \dot{\boldsymbol{\xi}}}\right)\_{\boldsymbol{o}} = \dot{\boldsymbol{p}\_{\boldsymbol{r}}} - \dot{\boldsymbol{p}\_{\boldsymbol{o}}} \tag{67}
$$

and, for the cross-derivatives,

$$
\begin{split}
\left(\frac{\partial \hat{\boldsymbol{p}} \cdot \mathbf{\dot{p}}}{\partial \tilde{\xi}}\right)\_{t} &= 0.25 \left(\dot{\boldsymbol{p}}\_{E} + \dot{\boldsymbol{p}}\_{TE} - \dot{\boldsymbol{p}}\_{O} - \dot{\boldsymbol{p}}\_{TO}\right); \; \left(\frac{\partial \dot{\boldsymbol{p}} \cdot \mathbf{\dot{\zeta}}}{\partial \tilde{\xi}}\right)\_{b} = 0.25 \left(\dot{\boldsymbol{p}}\_{E} + \dot{\boldsymbol{p}}\_{BE} - \dot{\boldsymbol{p}}\_{O} - \dot{\boldsymbol{p}}\_{BO}\right) \\
\left(\frac{\partial \boldsymbol{p}}{\partial \tilde{\xi}}\right)\_{o} &= 0.25 \left(\dot{\boldsymbol{p}}\_{T} + \dot{\boldsymbol{p}}\_{TE} - \dot{\boldsymbol{p}}\_{B} - \dot{\boldsymbol{p}}\_{BE}\right); \; \left(\frac{\partial \boldsymbol{p}}{\partial \tilde{\xi}}\right)\_{o} = 0.25 \left(\dot{\boldsymbol{p}}\_{T} + \dot{\boldsymbol{p}}\_{TO} - \dot{\boldsymbol{p}}\_{B} - \dot{\boldsymbol{p}}\_{BO}\right)
\end{split}
\tag{68}
$$

Introducing the discretization expressed by Eqs. (68) and 67 into Eq. (65) allows us to obtain the pressure correction equation:

$$a\_p \dot{p\_p} = a\_E \dot{p\_E} + a\_O \dot{p\_O} + a\_T \dot{p\_T} + a\_B \dot{p\_B} + b \tag{69}$$

where

$$a\_E = \left(\frac{\mathbf{g}\_{13}}{\left\langle \tilde{a}\_p \,\,/\rho \right\rangle} \right)\_e; \quad a\_O = \left(\frac{\mathbf{g}\_{13}}{\left\langle \tilde{a}\_p \,\,/\rho \right\rangle} \right)\_o; \quad a\_T = \left(\frac{\mathbf{g}\_{II}}{\left\langle \tilde{a}\_p \,\,/\rho \right\rangle} \right)\_l; \quad a\_B = \left(\frac{\mathbf{g}\_{II}}{\left\langle \tilde{a}\_p \,\,/\rho \right\rangle} \right)\_b \tag{70}$$

$$a\_P = a\_E + a\_O + a\_T + a\_B \tag{71}$$

$$\begin{split} b &= J \frac{\rho\_p^0 - \rho\_p}{\Delta t} - F\_e^\* + F\_o^\* - F\_r^\* + F\_b^\* - \\ \left( \frac{\mathbf{g}\_{13}}{\{\widetilde{a}\_p \, /\rho\}} \frac{\partial \boldsymbol{p}^\cdot}{\partial \widetilde{\boldsymbol{\gamma}}^\cdot} \right)\_e &+ \left( \frac{\mathbf{g}\_{13}}{\{\widetilde{a}\_p \, /\rho\}} \frac{\partial \boldsymbol{p}^\cdot}{\partial \widetilde{\boldsymbol{\gamma}}^\cdot} \right)\_o - \left( \frac{\mathbf{g}\_{13}}{\{\widetilde{a}\_p \, /\rho\}} \frac{\partial \boldsymbol{p}^\cdot}{\partial \widetilde{\boldsymbol{\gamma}}^\cdot} \right)\_l + \left( \frac{\mathbf{g}\_{13}}{\{\widetilde{a}\_p \, /\rho\}} \frac{\partial \boldsymbol{p}^\cdot}{\partial \widetilde{\boldsymbol{\gamma}}^\cdot} \right)\_b \end{split} \tag{72}$$

The symbol denotes a linear interpolation from the control volume centre (where the momentum equations are defined) to the control volume faces (where the fluxes for the continuity equations are needed). The pressure correction field *p'* obtained from the solution of Eq. (69) is employed for correcting pressure and velocity corrections are obtained from Eq. (63). Note that, since these equations are defined at the control volume centre, *p'* derivatives are evaluated as

A Numerical Procedure for 2D Fluid Flow Simulation in Unstructured Meshes http://dx.doi.org/10.5772/63077 279

$$\left(\frac{\partial \vec{p'}}{\partial \vec{\xi}}\right)\_p = 0.\mathfrak{S}\left(\vec{p\_e} - \vec{p\_o}\right); \quad \left(\frac{\partial \vec{p'}}{\partial \vec{\xi}'}\right)\_p = 0.\mathfrak{S}\left(\vec{p\_r} - \vec{p\_s}\right) \tag{73}$$

One may note that pressure values at the control volume centre are not included in the evaluation of these derivatives (the same applies for the pressure derivatives in the momentum equations). This may lead to the well-known checkerboard pattern for the pressure field. To avoid this effect, the Chie-Row interpolation method proposes that the mass fluxes, to be evaluated at the control volume interfaces for all the transport equations (*Fe*, *Fo*, *Ft* and *Fb*, in Eq. (48)), be corrected using the pressure correction *p'* field (instead of being computed from the corrected control volume centre velocities). The correction equations for fluxes are obtained from the correction equations for velocities:

$$F\_{e,o} = F\_{e,o}" - \left[ \left< \frac{\rho}{\tilde{a}\_p} \right> \left( \mathbf{g}\_{13} \frac{\partial \, \dot{p}}{\partial \tilde{\xi}} - \mathbf{g}\_{13} \frac{\partial \, \dot{p}}{\partial \tilde{\xi}} \right) \right]\_{e,o} \tag{74}$$

$$F\_{\mathfrak{r},b} = F\_{\mathfrak{r},b} \left[ \left< \frac{\rho}{\tilde{a}\_p} \right> \left( \mathbf{g}\_{11} \frac{\partial \mathbf{p}^\cdot}{\partial \tilde{\xi}^\cdot} - \mathbf{g}\_{13} \frac{\partial \mathbf{p}^\cdot}{\partial \tilde{\xi}^\cdot} \right) \right]\_{\mathfrak{r},b} \tag{75}$$

The "starred" fluxes *F\** at the control volume interfaces are obtained by interpolating the momentum equation. The keystone of the method is that the pressure gradient is not interpo‐ lated, but, instead, is obtained directly from the pressure at contiguous control volume centres. Considering Eq. (55), evaluated in terms of fluxes, one may write

$$\left|F\_{e,o}\right.^{\*} = \frac{F\_{e,o}^{m}}{I+E} + \hat{F}\_{e,o} + \left<\frac{\rho}{\tilde{a}\_p}\right>\_{e,o} \left(\mathbf{g}\_{13}\frac{\partial\left.\tilde{p}^{\*}\right|}{\partial\tilde{\xi}} - \mathbf{g}\_{13}\frac{\partial\left.\tilde{p}^{\*}\right|}{\partial\tilde{\zeta}}\right)\_{e,o} \tag{76}$$

$$\left(F\_{t,b}\right)^{\*} = \frac{F\_{t,b}^{m}}{I+E} + \hat{F}\_{t,b} + \left\langle \frac{\rho}{\tilde{a}\_P} \right\rangle\_{t,b} \left(\mathbf{g}\_{II}\frac{\partial p}{\partial \zeta} - \mathbf{g}\_{I3}\frac{\partial p}{\partial \zeta} \right)\_{t,b} \tag{77}$$

The terms *F* ^ *<sup>e</sup>*,*<sup>o</sup>* and *F* ^ *<sup>t</sup>*,*<sup>b</sup>* represent the fluxes

'' ' '

æö æö æö æö ç÷ ç÷ ç÷ ç÷ =- =- =- = èø èø èø èø *E P P O T P P B eo t b*

 ¶

' ' ' ' ' ' ' '

( ) ( )

*E TE O TO E BE O BO*

'' '' '' ''

 ¶

*p p 0 25 p p p p 0 25 p p p p*

*p p 0 25 p p p p 0 25 p p p p*

' ' '''

//// rrrr

*P P PP eotb*

'' ''

 r ¶x

////

*aaaa*

*PP PP eotb*

The symbol denotes a linear interpolation from the control volume centre (where the momentum equations are defined) to the control volume faces (where the fluxes for the continuity equations are needed). The pressure correction field *p'* obtained from the solution of Eq. (69) is employed for correcting pressure and velocity corrections are obtained from Eq. (63). Note that, since these equations are defined at the control volume centre, *p'* derivatives

¶¶¶¶

*gp gp gp gp*

æ öæ öæ öæ ö ç ÷ç ÷ç ÷ç ÷ + -+ è øè øè øè ø %% %%

*13 13 13 13*

*gggg <sup>a</sup> ; a ; a ; a aaaa* (70)

æö æöæöæö ==== ç÷ ç÷ç÷ç÷ èø èøèøèø % %%% *33 33 11 11*

*EOTB*

\*\*\*\*

*eot b*

 r ¶z  ¶

¶x

¶z

Introducing the discretization expressed by Eqs. (68) and 67 into Eq. (65) allows us to obtain

( ) ( )

*T TE B BE T TO B BO*

*PP EE OO TT BB ap ap ap ap ap b* = + +++ (69)

*P EOT B aaaaa* =+++ (71)

 r ¶x

'' '' '' ''

*pp p p p p; p p; p p; p p* (67)

¶z  ¶

¶z

(68)

(72)

¶¶

278 Numerical Simulation - From Brain Imaging to Turbulent Flows

and, for the cross-derivatives,

the pressure correction equation:

*0 P P*

r ¶z

are evaluated as

r r

*t*

D

*bJ F F F F*


¶x

' '

. ; .

æ ö æ ö ç ÷ = + -- ç ÷ = + -- ç ÷ ç ÷ è ø è ø

æ ö æ ö ç ÷ = + -- ç ÷ = + -- ç ÷ ç ÷ è ø è ø

. ; .

' '

*t b*

*e o*

¶x

¶

¶

where

¶z

¶x

$$
\hat{F}\_{e,o} = \boldsymbol{z}\_{\boldsymbol{\zeta}} \left< \rho \hat{u} \right>\_{e,o} - \boldsymbol{x}\_{\boldsymbol{\zeta}} \left< \rho \hat{w} \right>\_{e,o} \tag{78}
$$

$$
\hat{F}\_{\iota,b} = \chi\_{\varepsilon} \left< \rho \hat{w} \right>\_{\iota,b} - z\_{\zeta} \left< \rho \hat{u} \right>\_{\iota,b} \tag{79}
$$

where the revised velocities *û* and *ŵ* are obtained from the momentum equations as follows:

$$
\hat{a}\_P \hat{u}\_P = \sum\_{nb} a\_{nb} \hat{u}\_{nb} + \underbrace{a\_P}\_{} \underbrace{a\_p}\_{} - z\_\zeta \underbrace{\hat{\mathcal{Q}} \not{p} \not{q}}\_{} + z\_\xi \underbrace{\hat{\mathcal{Q}} \not{p} \not{q} \not{q}}\_{} + b \Rightarrow \hat{u}\_P = \frac{\sum\_{nb} a\_{nb} \hat{u}\_{nb} + b}{\tilde{a}\_P} \tag{80}
$$

$$\tilde{a}\_P \hat{w}\_P = \sum\_{nb} a\_{nb} \hat{w}\_{nb} + \frac{a\_P}{\mathcal{K}} \bigvee\_p^m + \mathbf{x}\_\zeta \bigvee\_{\tilde{\mathcal{K}}}^p - \mathbf{x}\_\tilde{\zeta} \bigvee\_{\tilde{\mathcal{K}}}^p + b \implies \hat{w}\_P = \frac{\frac{\sum a\_{nb} \hat{w}\_{nb} + b}{nb}}{\tilde{a}\_P} \tag{81}$$

Note that the source term *b* includes all the contributions (e.g., transient term, cross-derivatives and buoyancy for the *w* velocity component), except the under-relaxation and pressure gradient.

#### **4.5. Solution of the equations**

The solution of the equations previously described is carried out with an iterative procedure. For accelerating the convergence rate, two relaxation factors (described next) are applied.

The solution of the equation is sub- or over-relaxed in the following manner:

$$
\phi = f \phi\_{component} + (I - f)\phi\_{previous} \tag{82}
$$

Values of *f*lower than unity lead to sub-relaxation, while values greater than unity over-relax the solution process. In the present code, the value *f =1.1* is employed for the transport equations, while, for the pressure correction equation, the Point Successive Over Relaxation (PSOR) method [13] combined with the additive correction multigrid method [14] is employed.

The whole flow field calculation is considered to be converged when all the normalized residuals are lower that a predefined value *Rconv* :

$$\|\mathbb{I}R\_u, R\_w, R\_m, R\_\mathbb{I}, R\_k, R\_\varepsilon\| < R\_{\text{conv}}\tag{83}$$

The total normalized residual for the transport equations of *ϕ* (*ϕ* =*u*, *w*, *T* , *k*, *ε*, *ω*) is deter‐ mined as follows:

$$R\_{\phi} = \Sigma \left( \left| \frac{a\_{\rho} \phi\_{p} - \Sigma a\_{nb} u\_{nb} + b \right| J}{a\_{\rho} \left( \phi\_{\text{max}} - \phi\_{\text{min}} \right)} \right) / \Sigma\_{\text{all}} J \tag{84}$$

where (*ϕ*max −*ϕ*min) quantifies the amplitude of the variable *ϕ* in the calculation domain.
