**5.1. Flow past an airfoil**

, , , <sup>ˆ</sup> ˆ ˆ *t b t b t b F xw zu* = -

\*

\*

 z

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

x

x

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

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.

( ) *computed previous*

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

e

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

,,,,, *RRRRRR R uwmT k*

 f

¶

¶z- *<sup>p</sup> <sup>x</sup>*

<sup>+</sup> *<sup>p</sup> <sup>z</sup>*

¶

¶z

 r

\*

\*

ˆ

ˆ

+Þ = %

= +- *f 1f* (82)

< *conv* § ¨ (83)

+Þ = %

*nb P*

> *nb P*

(79)

ˆ

*au b b u <sup>a</sup>* (80)

*P*

ˆ

*aw b b w <sup>a</sup>* (81)

*P*

*nb nb*

å +

*nb nb*

å +

x

z

z

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

ff

residuals are lower that a predefined value *Rconv* :

<sup>+</sup> *<sup>p</sup> <sup>x</sup>*

¶

¶x

¶

¶x- *<sup>p</sup> <sup>z</sup>*

% ˆ ˆ = + <sup>å</sup> *<sup>P</sup> <sup>P</sup> <sup>m</sup> P P nb nb*

280 Numerical Simulation - From Brain Imaging to Turbulent Flows

% ˆ ˆˆ = + <sup>å</sup> *<sup>P</sup> <sup>P</sup> <sup>m</sup> P P nb nb*

*nb <sup>a</sup> aw a w w <sup>E</sup>*

**4.5. Solution of the equations**

gradient.

mined as follows:

*nb <sup>a</sup> au a u u <sup>E</sup>* r

A calculation was performed to compute the aerodynamic coefficients of a NACA 0012 airfoil operating at a Reynolds number of 6×106 . The obtained values for the drag and lift coefficients are compared with existing experimental data. The first step is to define the calculation domain, which should be large enough in order to avoid numerical blockage effects. **Figure 2** represents the domain, for a 1 m airfoil cord length. Lateral boundaries are assigned a free slip condition and a uniform velocity profile with 5% turbulence intensity is imposed at the inlet. A mass conservative condition is applied at the outlet boundary. After a mesh independency study, a total of approximately 250,000 mesh nodes were employed, with three mesh refinement regions. Particular care was taken near the airfoil surface, were *y+* values ranging from, typically, 0.1 to 6, with an average value of 1.7 all around, were obtained. **Figure 3** depicts the mesh employed.

For the present simulations, both the first-order hybrid [6] and the Quick advection schemes were employed, along with the *k-ε* and *k-ω* SST turbulence models. Experimental data are reported by Abbott and Von Doenhoff [15] and Ladson [16].

**Figure 2.** Domain dimensions. Airfoil cord is 1 m.

**Figure 3.** Non-structured quadrilateral mesh.

Results for the dependence of the lift coefficient with the airfoil angle of attack *α* are shown in **Figure 4**. As expected, the lift coefficient presents a linear dependence with the angle of attack *α* up to the onset of separation, which occurs at *α* =*16* . The two advection schemes give similar results up to separation, after which the lift drop in the stall region is more pronounced with the hybrid scheme. Separation is completely established at *α* =*18* and for *α* ≥*20* the flow becomes unsteady. Both turbulence models perform very well in the linear region. After separation, data are more spread. The difficulty to obtain reliable experimental data in this circumstance is commonly recognized as the flow is unsteady and presents a 3D behaviour. Comparing the two advection schemes, the Quick scheme performs better, particularly after separation.

**Figure 4.** Lift coefficient vs. angle of attack. Influence of (a) turbulence modeland (b) advection scheme.

**Figure 5** depicts the relation between lift and drag coefficients. In this case, the *k-ω* SST turbulence model performs substantially better than the *k-ε* model. This is certainly due to the fact that the friction component plays an important role in drag, and thus correctly resolving the boundary layer in the very proximity of the wall is crucial for the drag computation. It is also interesting to note that the advection scheme plays a very important role, with the higher order scheme Quick showing much better than the hybrid scheme, when results are compared with the experimental data.

**Figure 5.** Drag coefficient vs. lift coefficient. Influence of (a) turbulence model and (b) advection scheme.

#### **5.2. Natural convection inside a cavity**

**Figure 3.** Non-structured quadrilateral mesh.

282 Numerical Simulation - From Brain Imaging to Turbulent Flows

separation.

Results for the dependence of the lift coefficient with the airfoil angle of attack *α* are shown in **Figure 4**. As expected, the lift coefficient presents a linear dependence with the angle of attack *α* up to the onset of separation, which occurs at *α* =*16* . The two advection schemes give similar results up to separation, after which the lift drop in the stall region is more pronounced with the hybrid scheme. Separation is completely established at *α* =*18* and for *α* ≥*20* the flow becomes unsteady. Both turbulence models perform very well in the linear region. After separation, data are more spread. The difficulty to obtain reliable experimental data in this circumstance is commonly recognized as the flow is unsteady and presents a 3D behaviour. Comparing the two advection schemes, the Quick scheme performs better, particularly after

**Figure 4.** Lift coefficient vs. angle of attack. Influence of (a) turbulence modeland (b) advection scheme.

**Figure 5** depicts the relation between lift and drag coefficients. In this case, the *k-ω* SST turbulence model performs substantially better than the *k-ε* model. This is certainly due to the fact that the friction component plays an important role in drag, and thus correctly resolving the boundary layer in the very proximity of the wall is crucial for the drag computation. It is also interesting to note that the advection scheme plays a very important role, with the higher The natural convection flow in a cavity is a classical test for numerical methods in fluid dynamics. The cavity is a square shape (cf. **Figure 6**) with adiabatic horizontal walls. A constant temperature is imposed in each the vertical wall.

**Figure 6.** Problem definition for the natural convection inside a cavity.

The problem is governed by the following non-dimensional parameters:

$$\text{Prandtl number, } Pr = \frac{\nu}{a} \tag{85}$$

where the thermal diffusivity is *α* =*λ* /(*ρcp*), *ν* is the kinematic viscosity, *λ* is the thermal conductivity, *ρ* is the density and *cp* is the specific heat:

$$\text{Rayleigh number}, Ra = \frac{\text{g}\,\beta\,\Delta T\,L^{\circ}\text{Pr}}{\nu^{\circ}}\tag{86}$$

where *β* is the thermal expansion coefficient, *g* is the gravity acceleration and *ΔT* =*Th* −*Tc* is the temperature difference between the vertical walls. The transition between laminar and turbulent flow takes place approximately for *Ra* =*10 <sup>7</sup>* . In the present work, simulations were conducted for *Ra* =*10 <sup>5</sup>* (laminar regime). Air is the operating fluid, for which Pr *= 0.71*. The hybrid advection scheme was used and the Boussinesq approximation was adopted. Compu‐ tations were performed in a non-uniform grid, with 82 × 82 = 6400 nodes. Reference results are reported in [17] and [18] for several Rayleigh numbers in laminar regime, comparing solutions given by several authors. Results for laminar and turbulent flow are also presented in [19].

**Figure 7(a)** and **(b)** displays isothermal lines generated using a constant value spacing between the minimum and the maximum verified within the domain. **Figure 8(a)** and **(b)** shows the flow streamlines. The flow, in the steady-state situation, is characterized by a large vortex filling the cavity, rotating in the clockwise direction. Two small vortices rotating in the same direction are located near the cavity centre. For this case, the minimum and the maximum streamline values used in the visualization do not correspond to the total amplitude of the stream function within the domain. These values were, instead, adjusted in EasyCFD to correspond to those employed in [19]. The agreement between the calculations and those reported in the literature is very good. Vahl Davis and Jones [18] present normalized maximum values for the *u* and for *w* components of velocity

$$
\mu^\* = \frac{\mu L}{a}; \quad \mathbf{w}^\* = \frac{\mathbf{w}L}{a} \tag{87}
$$

occurring in the vertical and horizontal symmetry lines, respectively. **Table 1** shows the results obtained with EasyCFD, the reference values in [17] and the range of variation for the *37* contributions reported in [18]. This range does not include the minimum and maximum reported values since these clearly fall outside the general trend of the remaining contributions.

**Figure 7.** Isothermal lines. Ra = 105 . (a) EasyCFD and (b) Dixit and Babu [19].

**Figure 8.** Streamlines. Ra = 105 . (a) EasyCFD and (b) Dixit and Babu [19].


**Table 1.** Normalized maximum values for the *u* and *w* velocity components.
