2012
Validation Guide
©2012 Next Limit Technologies
Contents
1 Theory
1
2 Lid-driven cavity flow
9
3 Natural convection in a cavity
13
4 NACA-0012 airfoil at Re = 500
19
5 S825 airfoil
27
6 Vortex cell
35
7 Automotive aerodynamics
39
8 Multi-phase flows
45
ii
1 Theory
In the literature there are several particle-based numerical approaches to solve the computational fluid dynamics. They can be classified in three main categories: (i) algorithms modeling the behavior of the fluid at molecular level (e.g. Direct Simulation Montecarlo); (ii) algorithms which solve the equations at a macroscopic level, such as Smoothed Particle Hydrodynamics (SPH) or Vortex Particle Method (VPM), and finally, (iii) methods based on a mesoscopic framework, such as the Lattice Gas Automata (LGA) and Lattice Boltzmann Method (LBM). The algorithms that work at molecular level have a limited application, and they are used in theoretical analysis. The methods that solve macroscopic continuum equations are employed most frequently, but they also present several problems. SPH-like schemes are computationally expensive and in their less sophisticated implementations show lack of consistency and have problems imposing accurate boundary conditions. VPM schemes have also a high computational cost and besides, they require additional solvers (e.g. schemes based on boundary element method) to solve the pressure field, since they only model the rotational part of the flow. Finally, LGA and LBM schemes have been intensively studied in the last years being their affinity to the computational calculation their main advantage. Their main disadvantage is the complexity to analyze theoretically the emergent behavior of the system at macroscopic level from the laws imposed at mesoscopic level. XFlow’s approach to the fluid physics takes the main ideas behind these schemes and extends them to overcome most of the limitations present on these schemes.
1.1
Lattice Gas Automata
LGA schemes are simple models that allow to solve the behavior of gases. The main idea is that the particles move discretely in a d-dimensional lattice in one of
2
Theory
Figure 1.1: HPP model.
the predetermined direction at discrete times t = 0, 1, 2, ... and with velocity ci , i = 0, ..., b, also predetermined. The simplest model is the HPP, introduced by Hardy, Pomeau and de Pazzis, in which the particles move in a two-dimensional square grid and in four directions, as shown in Figure 1.1. The state of an element of the lattice at instant t is given by the occupation number ni (r, t), with i = 0, ..., b, being ni = 1 presence and ni = 0 absence of particles moving in direction i. The equation that governs the evolution of the system is as follows: ni (r + ci ∆t, t + ∆t) = ni (r, t) + Ωi (n1 , ..., nb )
(1.1)
where Ωi is the collision operator, which for each previous state (n1 , ..., nb ) computes C a post-collision state (nC 1 , ..., nb ) conserving the mass, linear momentum and energy; r is a position in the lattice and ci a velocity. From a statistical point of view, a system is constituted by a large number of elements which are macroscopically equivalent to the system under study. The macroscopic density and linear momentum are: b
ρ=
1X ni , b i=1
1.2
b
ρv =
1X n i ci b
(1.2)
i=1
Boltzmann’s transport equation
Boltzmann´s transport equation is defined as follows: fi (r + ci ∆t, t + ∆t) = fi (r, t) + ΩB i (f1 , .., fb )
(1.3)
where fi is the distribution function in the direction i and ΩB i the collision operator. From this equation and by means of the Chapman-Enskog expansion, the compressible Navier-Stokes equations can be recovered [1]. The Chapman-Enskog
1.3 Lattice Boltzmann Method
3
expansion shows that it is possible to design LGA schemes that recover the hydrodynamic macroscopic behavior at low Mach numbers. The main advantage of these methods is their great affinity with computers. They are easily programmed and very efficient. Some schemes have isotropy problems (do not satisfy Galilean invariance) and produce very noisy results. The main contribution of LGA schemes is that they are precursor of the Lattice Boltzmann method.
1.3
Lattice Boltzmann Method
The origins of the Lattice Boltzmann Method (LBM) ([2, 3, 4]) lie in the LGA schemes. While the LGA schemes use discrete numbers to represent the state of the molecules, the LBM method makes use of statistical distribution functions with real variables, preserving by construction the conservation of mass, linear momentum and energy. It can be shown that if the collision operator is simplified under the Bhatnagar-Gross-Krook (BGK) approximation [5], the resulting scheme reproduces the hydrodynamic regime also for low Mach numbers. This operator is defined as follows: 1 (1.4) ΩBGK = (fieq − fi ) i τ eq where fi is the local equilibrium function and τ is the relaxation characteristic time (which is related to the macroscopic viscosity). Usually, the equilibrium distribution function adopts the following expression: ciα vα vα vβ ciα ciβ eq ( 2 − δαβ ) (1.5) fi (r, t) = ti ρ 1 + 2 + cs 2c2s cs where cs is the sound speed, v the macroscopic velocity, δ the Kronecker delta, and ti are built preserving the isotropy in space.
(a) D2Q7
(b) D2Q9
Figure 1.2: Most common LBM schemes in two-dimensions. LBM schemes are classified as a function of the spatial dimensions d and the number of distribution functions b, resulting the notation DdQb. The most common
4
Theory
(a) D3Q19
(b) D3Q27
Figure 1.3: Most common LBM schemes in three-dimensions.
schemes in two dimensions are the D2Q7 and D2Q9 represented in Figure 1.2, while in three dimensions the most used schemes are the D3Q13, D3Q15, D3Q19 and D3Q27 plotted in Figure 1.3. Finally, the multiscale Chapman-Enskog expansion gives us the relation between the macroscopic viscosity and the relaxation parameter: 1 ν = c2s (τ − ) 2
(1.6)
For a positive viscosity, the relaxation time must be greater than 0.5. The most interesting aspect is that these schemes are able to model a wide range of viscosities (0, ∞) in an efficient way even using explicit formulations. General references are the review by Chen and Doolen (1998) [6] and the book by Succi (2001) [7].
1.4
Turbulence modeling
Following the dimensional analysis proposed by Kolmogorov at high Reynolds numbers, the flow tends to break in smaller eddies to transform the kinetic energy into internal energy. This process is known as ”Kolmogorov cascade” and it explains the turbulence phenomenon. The time necessary to break an eddy in the flow is in the order of Tbreak ∼
Leddy νeddy
(1.7)
and the time to dissipate the kinetic energy through viscosity is expressed as Tdissip,visc ∼
L2eddy ν
(1.8)
1.4 Turbulence modeling
Problem
Re Nelements
5
Blood flow Bioengineering
Vehicle airflow at low speed
Vehicle airflow
Plane airflow
103 - 104 106 - 109
105 - 106 1011 - 1013
107 - 108 1015 - 1018
> 109 > 1020
Table 1.1: Number of elements required for DNS.
For large eddies and high Reynolds numbers, the break time is smaller than time employed to dissipate the energy and this produces the Kolmogorov cascade. The kinetic energy of a turbulent structure can be estimated by 2 Ec eddy ∼ Veddy
The specific kinetic energy dissipation ratio is as follows εbreak ∼
2 Veddy
Tbreak
=
3 Veddy
Leddy
The smallest eddies present in the flow (of size ∼ Lcritical ) have a break time equal to time necessary to transform their kinetic energy to viscous energy (Tbreak ∼ Tdissip,visc ). Then the kinetic energy dissipation ratio can be estimated by 3 3 Veddy Vcritical ν3 ∼ ∼ 4 ∼ εbreak Lcritical Leddy Lcritical
(1.9)
3 L3 3 Veddy Leddy L3critical 3 eddy Lcritical ∼ = Reeddy 3 Lcritical ν 3 Leddy L3eddy
(1.10)
Lcritical −3/4 ∼ Reeddy Leddy
(1.11)
εviscous ∼ and thus,
Finally,
Taking into account this relation, if we want to solve explicitly every eddy in a three-dimensional flow, the number of elements are in the order of Nelements ∼
Leddy Lcritical
3
9/4
∼ Reeddy
(1.12)
Table 1.1 summarizes the number of elements necessary to solve in a direct way the turbulence at different Reynolds numbers. Direct numerical simulation (DNS) of turbulence is currently possible only for low Reynolds numbers, while for typical industrial Reynolds numbers some modeling for the unresolved scales is required. Reynolds-Averaged Navier-Stokes (RANS) approach models the turbulence in a global way. This approach is nowadays the most widely adopted and calculates
6
Theory
Figure 1.4: Turbulence modeling.
values averaged in time, removing the time dependence of the solution. Although calculating averaged results is computationally less expensive, new terms appear in the Navier-Stokes equations that have to be modeled by new transport equations. Moreover there are several RANS models, each one suitable for a specific problem, and the parameters of each model need to be adjusted empirically. Another approach to the turbulence problem is the Large Eddy Simulation (LES). These schemes solve the turbulence in a local way, modeling only the smallest scales, and are closer to the physics (see Figure 1.4). The turbulence at smallest scales has been extensively studied and its behavior can be reproduced without using arbitrary parameters. These are the type of schemes employed in XFlow. The LES scheme adopted by default by XFlow is the Wall-Adapting Local Eddyviscosity (WALE), which has good properties both near and far of the wall and with laminar and turbulent flows. This model is formulated as follows: νef f ective = νmolecular + νturbulent νturbulent =
Sαβ
∆2f ilterscale
(Gdαβ Gdαβ )3/2
(Sαβ Sαβ )5/2 + (Gdαβ Gdαβ )5/4 1 ∂vα ∂vβ = + 2 ∂rβ ∂rα
1 2 1 2 2 Gdαβ = (gαβ + gβα ) − δαβ gγγ 2 3 gαβ =
∂vα ∂rβ
(1.13) (1.14)
(1.15) (1.16) (1.17)
REFERENCES
7
References [1] Y.H. Qian, D. D’Humieres, and P. Lallemand. Lattice BGK models for NavierStokes equation. Europhysics Letters, 17:479, 1992. [2] G. McNamara and G. Zanetti. Use of the Boltzmann equation to simulate latticegas automata. Physical Review Letters, 61:2332–2335, 1988. [3] F.J. Higuera and J. Jimenez. Boltzmann approach to lattice gas simulations. Europhysics Letters, 9:663–668, 1989. [4] H. Chen, S. Chen, and W. Matthaeus. Recovery of the Navier-Stokes equations using a lattice-gas Boltzmann method. Physical Review A, 45:5339, 1992. [5] P.L. Bhatnagar, E.P. Gross, and M. Krook. A model for collision processes in gases. Phys. Rev., 94:511, 1954. [6] S. Chen and G. Doolen. Lattice Bolzmann method for fluid flows. Annual Reviews of Fluid Mechanics, 30:329–64, 1998. [7] S. Succi. The Lattice Boltzmann Equation for Fluid Dynamics and Beyond. Clarendon Press, 2001.
8
REFERENCES
2 Lid-driven cavity flow
The lid-driven cavity is a classical benchmark problem for viscous incompressible fluid flow. It consists in a cavity where the upper boundary moves to the right, and causes a rotation in the cavity. Side and bottom walls of the cavity are considered no-slip, while the velocity at the upper wall is imposed to vx = 1 m s−1 (see Figure 2.1). Although there is a discontinuity of the boundary conditions at the two top corners where the side wall meet the lid, this corner singularity plays a minor role in the overall solution field.
Figure 2.1: Geometry and boundary conditions. The solution field will depend on the Reynolds number Re =
ρ vx L µ
(2.1)
This validation case computes the laminar incompressible flow for a 2D driven cavity at Re = 1000 and compares XFlow’s results with the ones presented in [1] and [2]. Fluid properties are ρ = 1 kg m−3 and µ = 0.001 Pa s, in a L = 1 m square cavity discretized with 128x128, 256x256 and 512x512 cells lattices.
10
Lid-driven cavity flow
y
128x128
256x256
512x512
Ref. [1]
Ref. [2]
0.9688 0.9531 0.7344 0.5 0.2813 0.1016 0.0625
0.58578 0.476 0.18963 -0.06255 -0.28175 -0.29968 -0.20118
0.57217 0.46869 0.18866 -0.06218 -0.28020 -0.30286 -0.20493
0.57970 0.47213 0.18908 -0.06216 -0.28064 -0.30078 -0.20248
0.57492 0.46604 0.18719 -0.0608 -0.27805 -0.2973 -0.20196
0.58031 0.47239 0.18861 -0.06205 -0.2804 -0.30029 -0.20227
Table 2.1: Horizontal velocity vx along the vertical centerline. Comparison of XFlow’s results for different resolutions with the reference solutions.
The solution computed with XFlow is plotted in Figure 2.2. The vx velocity
Figure 2.2: Velocity field.
distribution along vertical centerline is shown in Figure 2.3 and Table 2.1. Additional quantitative comparisons of pressure and vorticity are presented in Tables 2.2 and 2.3. For this problem it is essential to activate the high order boundary conditions (Project Tree > Engine > Advanced Options > High order boundary conditions). This imposes the velocity at the lid through a high order scheme. While using high order BC the velocity relative error is O(10−3 ), with the default BC the relative error is O(10−2 ).
11
Figure 2.3: Horizontal velocity vx distribution along the vertical centerline. y
128x128
256x256
512x512
Ref. [1]
0.9688 0.9531 0.7344 0.5 0.2813 0.1016 0.0625
0.055238 0.053995 0.013019 0.0 0.043511 0.112102 0.117450
0.054295 0.053059 0.012739 0.0 0.042518 0.110042 0.115454
0.054077 0.052855 0.012738 0.0 0.042447 0.109613 0.114904
0.051514 0.050329 0.012122 0.0 0.040377 0.104187 0.109200
Table 2.2: Pressure p along the vertical centerline for Re = 1000. Comparison of XFlow’s results for different resolutions with a reference solution. y
128x128
256x256
512x512
Ref. [2]
0.9688 0.9531 0.7344 0.5 0.2813 0.1016 0.0625
9.51527 5.04854 2.10078 2.07863 2.27638 1.64462 2.32787
9.11829 4.63587 2.09089 2.06865 2.26145 1.61598 2.30427
9.43159 4.81708 2.09418 2.07152 2.26792 1.63825 2.31797
9.47810 4.86280 2.09090 2.06690 2.26780 -1.63520 -2.31740
Table 2.3: Vorticity ω along the vertical centerline. Comparison of XFlow’s results for different resolutions with a reference solution.
12
REFERENCES
References [1] U. Ghia, K.N. Ghia, and C.T. Shin. High-Re solutions for incompressible flow using the Navier-Stokes equations and a multigrid method. Journal of Computational Physics, 48:387–411, 1982. [2] C.-H. Bruneau and M. Saad. The 2d lid-driven cavity problem revisited. Computers & Fluids, 35:326–348, 2006.
3 Natural convection in a cavity
Buoyancy-driven flow in a square cavity with vertical sides which are differentially heated allows to test XFlow in thermal problems where the buoyancy forces are modeled using the Boussinesq approximation. The aim is to calculate the flow and thermal field for Rayleigh numbers of 103 and 106 , and compare the results with the benchmark solutions published by De Vahl et al. [1].
Figure 3.1: Geometry and boundary conditions. The boundary conditions for the problem (see Figure 3.1) involve two vertical walls at different temperatures, leading to a thermal gradient across the domain. This thermal gradient produces varying buoyancy forces between the walls that drive the flow. Horizontal walls are considered adiabatic. The Boussinesq approximation assumes the linear variation of the density as a function of the temperature θ, ρ = ρ0 [1 − α(θ − θ0 )] and that the thermodynamic properties of incompressible fluids are constant except
14
Natural convection in a cavity
when considering the body force ρg in the momentum equation: ∇·v =0
(3.1)
dv = ∇ · σ + ρ0 [1 − α(θ − θ0 )]g dt dθ ρ0 cp = k∇ · (∇θ) + Φ dt ρ0
(3.2) (3.3)
where α is the thermal expansion coefficient, cp the specific heat at constant pressure, k the thermal diffusion coefficient, θ0 a reference temperature, ρ0 = ρ(θ0 ), and Φ is the viscous heat dissipation. Benchmark solutions for the natural convection problem were published by [1]. These solutions are compared against the numerical results obtained with XFlow. The non-dimensional distances and velocities used in the results are X= with κ =
x , L
Y =
y , L
Vx =
vx L , κ
and Vy =
vy L κ
k the thermal diffusivity. The Rayleigh number is defined as ρ cp Ra =
ρ g α ∆θ L3 µκ
(3.4)
∆θ = |θ1 − θ2 | is the temperature difference between the two vertical walls. For the simulations presented here, the solution domain consists of a L = 1 m square, using 20x20, 40x40, 80x80, and 160x160 cell lattices. The fluid properties for the Ra = 103 case are: ρ0 = 1 kg m−3 , g = −10 m s−2 , α = 0.1 K−1 , ∆θ = 1 K (e.g. θ1 = 293.65 K, θ2 = 292.65 K), µ = 0.0266458 Pa s, k = 0.0375293 kg m s−3 K−1 , cp = 1 m2 s−2 K−1 µ cp = The last three values are chosen to match the Prantl number of air, P r = k 6 0.71. For the Ra = 10 case, only change µ = 0.0008426 Pa s and k = 0.0011868 kg m s−3 K−1 . The viscous heat dissipation has been neglected as in [1]. The following quantities will be analyzed at the steady state: N u0 , average heat flux in the hot vertical wall
Z N u0 = 0
1
∂θ dy ∂x X=0
15
max Vx , maximum value of the non-dimensional horizontal velocity on the vertical centerline and its location Ymax max Vy , maximum value of the non-dimensional vertical velocity on the horizontal centerline and its location Xmax
Tables 3.1 and 3.2 show the comparison of XFlow’s results with the benchmark solution published in [1] for Ra = 103 and Ra = 106 . There is a good agreement in both cases although the accuracy decreases for the highest Rayleigh number. This deterioration is consistent with the findings of most contributors reported in [1]. In this test, the maximum and minimum locations are cell centered, none of the interpolation techniques suggested by [2] have been employed.
max Vx at Ymax max Vy at Xmax N u0
20x20
40x40
80x80
160x160
Ref.[1]
3.5428 0.825 3.6359 0.175 1.0916
3.6325 0.8125 3.6941 0.1875 1.1072
3.6422 0.8187 3.7023 0.1812 1.1129
3.6468 0.8156 3.7012 0.1781 1.1155
3.649 0.813 3.697 0.178 1.117
Table 3.1: Case Ra = 103 . Comparison of XFlow’s results for different resolutions with the benchmark solution.
max Vx at Ymax max Vy at Xmax N u0
20x20
40x40
80x80
160x160
Ref.[1]
54.448 0.875 127.259 0.075 5.924
57.750 0.863 193.479 0.038 7.975
63.805 0.856 215.178 0.044 8.688
64.902 0.847 219.817 0.041 8.823
64.630 0.850 219.360 0.038 8.817
Table 3.2: Case Ra = 106 . Comparison of XFlow’s results for different resolutions with the benchmark solution. For this case it is essential to deactivate the viscous heat dissipation Φ (in Project Tree > Engine > Advanced Options > Enable viscous term in Energy equation: Off), set a Courant ≈ 1, and do not use turbulence model neither wall models (Off-resolved). Too large Courant values introduce an excessive compressibility in the flow and the results do not correspond to the Boussinesq incompressible reference solution.
16
Natural convection in a cavity
(a)
(b)
(c)
(d)
Figure 3.2: Result fields for Ra = 103 : (a) temperature, (b) horizontal velocity, (c) vertical velocity, and (d) vorticity.
17
(a)
(b)
(c)
(d)
Figure 3.3: Result fields for Ra = 106 : (a) temperature, (b) horizontal velocity, (c) vertical velocity, and (d) vorticity.
18
REFERENCES
References [1] G. DeVahl Davis and I.P. Jones. Natural convection in a square cavity: a comparison exercise. International Journal for Numerical Methods in Fluids, 3:227–248, 1983. [2] G. DeVahl Davis. Natural convection of air in a square cavity: a benchmark numerical solution. International Journal for Numerical Methods in Fluids, 3:249–264, 1983.
4 NACA-0012 airfoil at Re = 500
The NACA-0012 airfoil is a widely-used wing section that has zero camber and a maximum thickness to chord ratio of 12 percent. Its two-dimensional profile is symmetric and has very smooth aerodynamic shape, as shows Figure 4.1. This validation case presents the results for the flow past a NACA-0012 at zero angleof-attack and Reynolds number 500 using XFlow, and compares the results with reference data. Due to lack of experimental data for such a low Reynolds number, the comparison will be based on the CFL3D code from the National Aeronautics and Space Administration (NASA) [1, 2].
Figure 4.1: NACA-0012 profile. This two-dimensional single phase external aerodynamics analysis has been run using a virtual wind tunnel of dimensions 60 × 40 m and a NACA-0012 profile of chord length L = 1 m. In order to reach a Reynolds number equal to 500 based on the chord length, the simulation parameters have been set according to Table 4.1. The spatial resolution chosen is 2.56 m for the far field, and 0.005 m around the airfoil profile and within the wake area, as shown in Figure 4.2. The spatial discretization has been defined through a region of refinement instead of using an adaptive refinement in order to ensure that the symmetry of the NACA-0012 is respected. The discretization ended up with 1.3 million elements in 9 levels of refinement. Since the flow is laminar, no wall functions have been used to model the boundary layer. Due to the fact that XFlow solver is inherently transient, the analysis has been
20
NACA-0012 airfoil at Re = 500
Free-stream velocity Density Dynamic viscosity Chord length Reynolds number Angle-of-attack
vref ρ µ L Re α
50 m s−1 1 kg m−3 0.1 Pa s 1m 500 0 degree
Table 4.1: Simulation conditions.
run until the aerodynamic coefficients stabilize in time. The time step used was 0.004 s, which corresponds to a Courant number of 1 with respect to the lattice size and the free-stream velocity.
Figure 4.2: Spatial resolution in XFlow. Profiles of normalized X and Y velocity components, as well as pressure coefficients, have been analyzed at five vertical sections: x/L = 0, 0.25, 0.5, 0.75 and 1 (see Figure 4.3). Figure 4.4 shows the X-component of the velocity field normalized by the reference velocity vref = 50 m s−1 . The results of XFlow and CFL3D are perfectly matching for the five sections. The profiles are as expected in both codes: they tend to zero in the airfoil thickness and to one (or slightly more) on the sides where the boundary layer is fully developed.
21
Figure 4.3: Data plot lines.
Cd Cl
XFlow
CFL3D
Objective
0.1705 10−13
0.1741 −0.538 × 10−5
0.1741 0
Table 4.2: Aerodynamic coefficients comparison.
For the Y-component of the velocity field (normalized by vref ), again XFlow results are almost perfectly matching with those of CFL3D, as shown in Figure 4.5. Nevertheless, one can observe some differences close to the airfoil wall, specially for section x/L = 0.5. This might be due to the size of the first element within the boundary layer which is not fine enough. However, the differences between the two codes are small. pstatic , where pstatic is the The pressure coefficient Cp is defined as Cp = 1 2 2 ρ vref gauge static pressure. Figure 4.6 shows the pressure coefficient distribution at the five sections. More differences are now noticeable, especially at x/L = 0.5, 0.75 and 1.0. In general, the Cp tends to be slightly over-estimated. The aerodynamic coefficients predicted by XFlow are quite similar to the ones from CFL3D, see Table 4.2. The drag coefficient has a relative error of -2.0678% with respect to CFL3D results, whereas the lift coefficient is actually even more accurate since it should be equal to zero due to the symmetry of the NACA profile at zero angle-of-attack. Aknowledgements: Validation data have been kindly provided by courtesy of NASA Langley Research Center and David P. Lockard.
22
NACA-0012 airfoil at Re = 500
Figure 4.4: X-component of velocity u(x, y) at x/L = 0.0, 0.25, 0.50, 0.75 and 1.0.
23
Figure 4.5: Y-component of velocity v(x, y) at x/L = 0.0, 0.25, 0.50, 0.75 and 1.0.
24
NACA-0012 airfoil at Re = 500
Figure 4.6: Pressure coefficient at x/L = 0.0, 0.25, 0.5, 0.75 and 1.0.
REFERENCES
25
References [1] C. Rumsey, R. Biedron, and J. Thomas. CFL3D: Its history and some recent applications. Technical report, NASA TM-112861, 1997. [2] D.P. Lockard, L.-S. Luo, S.D. Milder, and B.A. Singer. Evaluation of PowerFLOW for Aerodynamic Applications. Journal of Statistical Physics, 107(1/2):423–478, 2002.
26
REFERENCES
5 S825 airfoil
The S825 airfoil has been designed for horizontal-axis wind turbine applications by the National Renewable Energy Laboratory (Colorado, USA). The report of the design and experimentation of the S825 airfoil [1] exposes the different objectives and constraints set for the design, as well as the methodology of measurements which have been conducted in the NASA Langley Low-Turbulence Pressure Tunnel (LTPT) [2]. As explained in [1], the main objectives were, first, to reach a maximum lift coefficient of at least 1.40 at a Reynolds number of 2 × 106 . Second, a low profiledrag coefficients should be obtained between 0.40 and 1.20 of the lift coefficient. Two main constraints were to keep the zero-lift pitching-moment coefficient greater than -0.15, and also to have an airfoil thickness equal to 17% of the chord. The final twodimensional design is as shown in Figure 5.1, with a chord length equal to 0.45715 m.
Figure 5.1: S825 airfoil shape.
Experiments have been conducted at different Reynolds numbers based on the chord length, however this validation case will only treat the Reynolds number 2×106 since it has been used for most of the data provided by [1]. The Mach number is 0.1 and the experimentation has been done with transition free (smooth) and with transition fixed by roughness at specific locations. The objective of this case is to validate the pressure distribution and aerodynamic forces predicted by XFlow at low Mach number and different angles of attack (AoA). The calculations have been performed with XFlow for a range of angles of attack
28
S825 airfoil
between -4 and 10 degrees every two degrees. All the calculations are transient due to the nature of the XFlow solver and use the Wall-Adapting Local-Eddy turbulence model, which belongs to the Large Eddy Simulation (LES) approach. Wall function models in XFlow assume that the boundary layer is fully turbulent, therefore it is not possible to model transition or prescribe a transition location. Two-dimensional single phase analyses have been performed using a virtual wind tunnel of 60 m × 40 m and a velocity at the inlet of 43.7493 m s−1 . The angle of attack is varying by rotating the geometry instead of projecting the inlet velocity vector, since XFlow allows easy manipulation of the geometry. The fluid has a density of 1 kg m−3 and a dynamic viscosity of 10−5 Pa s accordingly to the Reynolds number based on the airfoil chord length (Re = 2×106 ). The simulations have been run for 1 second of physical time, with a time step of 0.002 s. The resolution scale at the far field is 1.28 m, using the refinement near walls and dynamically adapting to the wake available in XFlow. The walls and the wake are resolved with a scale of 0.0025 m, as shown in Figure 5.2.
Figure 5.2: Resolution refinement near the airfoil and the wake. The solution for the static pressure and velocity flow variables at final time for zero angle of attack can be observed in Figure 5.3. For each angle of attack, the curve of pressure coefficients (Cp ) has been extracted in XFlow using a cutting plane field distribution which projects the selected field on the upper and lower sides of the airfoil. The Cp has been computed as following, being Vref equal to 43.7493 m s−1 : pstatic Cp = 1 (5.1) 2 2 ρ Vref
29
(a) Static pressure
(b) Velocity Figure 5.3: Static pressure and velocity flow fields at final time for AoA = 0 degrees.
For the angles of attack -4, -2, 0, 2, 4, 6, 8 and 10 degrees, the pressure coefficient distribution along the airfoil has been compared with the transition free experimental data presented in [1]. The results for angle of attack between 0 and 6 degrees are in good agreement with the experiments, as shows Figure 5.4. On the upper side of the airfoil, the pressure coefficients are slightly under-estimated when the angle increases but still match reasonably with the experimental data. However, when angles are increased to 8 and 10 degrees of AoA then results are getting less accurate, as shown in Figure 5.5. The pressure coefficient tends to be more under-estimated near the leading edge of the upper part. This could be explained by the lack of transition model or the LES model, which is not fully consistent for 2D simulations. Another last series of angles of attack have been studied, this time for negative incidence. Again, XFlow predicts with accuracy the pressure coefficient distribution for -2 and -4 degrees, as shown in Figure 5.6. Finally, Figure 5.7 compares the angle of attack vs. lift coefficient for theoretical [1], experimental [1] and XFlow results. For positive angles, the lift coefficient is slightly over-predicted by XFlow but in good agreement with the theoretical results.
30
S825 airfoil
(a)
(b)
(c)
(d)
Figure 5.4: Airfoil pressure coefficient distribution for different AoA: a) 0 degrees, b) 2 degrees, c) 4 degrees, d) 6 degrees.
31
(a)
(b)
Figure 5.5: Airfoil pressure coefficient distribution for different AoA: a) 8 degrees, b) 10 degrees.
(a)
(b)
Figure 5.6: Airfoil pressure coefficient distribution for different AoA: a) -2 degrees, b) -4 degrees.
32
S825 airfoil
Figure 5.7: AoA (α) vs. lift coefficient (Cl ) for theoretical, experimental and XFlow results.
REFERENCES
33
References [1] D. Somers. Design and experimental results for the s825 airfoil. Technical report, National Renewable Energy Laboratory, 2005. [2] R. J. McGhee, W. D. Beasley, and J. M. Foster. Recent modifications and calibration of the langley low-turbulence pressure tunnel. Technical report, NASA TP-2328, 1984.
34
REFERENCES
6 Vortex cell
Trapping vortices is a technique that prevents vortex shedding in flows past bluff bodies. Vortices forming near bluff bodies tend to be shed downstream but, if the vortex is kept near the body at all times, it is called trapped. This validation case compares XFlow results with experimental data for the flow inside a vortex trapping cavity (vortex cell). The geometry consists in a rectangular channel of section 520 × 52 mm with a spherical vortex cell of 45 mm depth located at mid-length as shown in Figure 6.1. The boundary condition at the inlet is set to a constant fluid velocity of UL = 36 m/s and the gauge pressure to 0 Pa at the outlet. The fluid has been initialized to UL in the whole domain except inside the vortex cell where it is 0 m/s in order to reach quicker the pseudo-steady state.
Figure 6.1: Vortex cell geometry. The experimental data from [1, 2] provide the normalized X-component of the velocity measured along a vertical line going from the bottom of the sphere up to the upper wall of the channel (see Figure 6.2). The vertical coordinate along the line is normalized by the line length L = 52 mm. The two-dimensional studies led by [3] based on steady RANS turbulence models show how sensitive are the numerical results for the vortex cell flow depending on the turbulence model and the choice of the numerical scheme. XFlow uses
36
Vortex cell
Figure 6.2: Line for measurement.
Large Eddy Simulation turbulence models, which are inherently three-dimensional. Furthermore, although the vortex cell seems a two-dimensional flow, turbulence effects (important near and inside the cavity) need a three-dimensional analysis to be accurately modeled. Unfortunately three-dimensional analyses may involve a large number of elements and long simulation times. The refinement algorithms available in XFlow (near the walls and adaptive wake) allow to minimize the number of elements, but tend to introduce numerical dissipation when passing from one element size to another and has been found to be inaccurate especially in the boundary layer that detaches from the leading edge of the vortex cell. The following results were obtained using a uniform resolution of 1 mm in the whole domain and a time step of 10−6 s. This resolution leads to a total of 1.1 million elements. The total simulation time solved is 0.9 s at a frequency of 500 Hz. Averaged results are required in order to analyze the pseudo-steady state of the solution. Figure 6.3 shows that XFlow 3D results are globally in good agreement with the experimental data. The areas of less accuracy are at y/L around -0.7 and the peak around 0.1. Nevertheless XFlow is able to predict the experimental velocity profile at the cavity entry (−0.4 < y/L < 0) better than the RANS calculations and correctly predicts the vortex speed at y/L = −0.8. Finally, Figure 6.4 shows the averaged velocity field computed by XFlow. It is possible to observe the creation of vortices at the leading edge of the cavity, what is not evident in steady calculations.
37
Figure 6.3: Comparison of XFlow 3D results with experimental data and RANS results from [3].
Figure 6.4: XFlow averaged velocity field.
38
REFERENCES
References [1] P. A. Baranov, S. V. Guvernyuk, M. A. Zubin, and S. A. Isaev. Numerical and physical modeling of the circulation in a vortex cell in the wall of a rectilinear channel. Fluid Dynamics, 35:663–673, 2000. [2] S.A. Isaev, S.V. Guvernyuk, M.A. Zubin, and Y.S. Prigorodov. Numerical and physical modeling of a low-velocity air flow in a channel with a circular vortex cell. Journal of Engineering Physics and Thermophysics, 73:337–344, 2000. [3] R. Donelli, P. Iannelli, S. Chernyshenko, A. Iollo, and L. Zannetti. Flow models for a vortex cell. AIAA Journal, 47:451–467, 2009.
7 Automotive aerodynamics
In automotive aerodynamics it is usual to use reference geometries to validate CFD codes [1]. This section uses the well known ASMO model, which comprises a squareback rear, smooth surfaces, boat tailing, underbody diffuser and no pressure induced boundary layer separation. The geometry does not have a well defined separation line and is characterized by a low drag shape. For this model, experimental data from Daimler Benz and Volvo model scale wind tunnel are available. Mesh is one of the major issues in classic CFD approaches for aerodynamical problems. Mesh quality may be as important as mesh resolution when high accuracy in the calculations is aimed for. The particle-based approach of XFlow avoids the costly generation of a good mesh. For this validation the 1/5 wind tunnel test model has been adopted. Vehicle’s length, width and height are 0.81 m, 0.29 m and 0.27 m respectively, while wind tunnel dimensions are 9 x 1.5 x 3 m (see Figure 7.1). This corresponds to a blockage ratio of 1.38%. The wind tunnel domain type available in XFlow is used, with an inlet uniform velocity of 50 m/s. The fluid properties are density ρ = 1 kg m−3 and dynamic viscosity µ = 1.5 × 10−5 Pa s.
Figure 7.1: ASMO body and wind tunnel geometries.
40
Automotive aerodynamics
The Reynolds number for this case is 2.7 × 106 , taking the length of the vehicle as reference. In this problem it is essential to resolve the turbulent wake properly. Particle resolution in the far field is 0.1 m, while in the wake and on the model surface scales up to 2.5 mm are resolved. Dynamic wake refinement is applied, so that the specified particle resolution is automatically adopted in regions with high turbulence, while less turbulent regions are treated with fewer particles. In addition, XFlow uses a Large Eddy Simulation (LES) approach for modeling the turbulence, in particular the Wall-Adapting Local Eddy-viscosity (WALE) model (see Section 1). The turbulent wake structure can be observed in Figures 7.2 and 7.3, together with the instantaneous pressure field and the skin friction distribution in Figure 7.4.
Figure 7.2: Snapshot of isosurface of vorticity.
41
(a)
(b)
(c) Figure 7.3: Instantaneous velocity field: (a) on the vehicle surface, (b) in the Y = 0.07 m plane, and (c) in the symmetry plane.
42
Automotive aerodynamics
Figure 7.4: Instantaneous pressure field in the symmetry plane and skin friction distribution.
Figures 7.5 to 7.8 show validation results of surface pressure measurements. XFlow results are compared with experimental data obtained from Volvo and Daimler Benz. Data are available in the symmetry plane and are shown for roof, underbody, front and base region of the vehicle. It can be seen that the comparison with the measurements is good, although some deviations can be observed especially in the base pressure, which is slightly underpredicted. However the proper level of the base pressure is not known exactly, as there is a large difference between both experiments.
Figure 7.5: Front pressure distribution along the symmetry plane.
43
Figure 7.6: Roof pressure distribution along the symmetry plane.
Figure 7.7: Base pressure distribution along the symmetry plane.
Figure 7.8: Underbody pressure distribution along the symmetry plane.
44
REFERENCES
Typically, drag stabilizes in a characteristic time of the order of the flow traveling the vehicle length. In 0.1 seconds, the flow has traveled more than six times the whole body. The time averaged drag between 0.05 and 0.1 seconds of physical simulation is Cd = 0.151 (see Figure 7.9), in good agreement with the values measured in the experiments shown in Table 7.1. [2] showed that using transient CFD simulations, surface pressure values can be computed fairly accurate. The overall drag coefficient however, is not predicted satisfactorily. RANS turbulence models tend to overestimate the drag. LES transient calculations are becoming an integral part of the aerodynamic development process and affordable using XFlow.
Figure 7.9: Overall drag history. XFlow Experiments Volvo Experiments Daimler Benz
0.151 0.158 0.153
Table 7.1: Drag values for ASMO model.
References [1] G. Le Good and K. Garry. On the use of reference models in automotive aerodynamics. SAE paper, 2004-01-1308. [2] S. Perzon and L. Davidson. On transient modeling of the flow around vehicles using the reynolds equations. In ACFD 2000 Beijing, pages 720–727, 2000.
8 Multi-phase flows
When heavy fluid lies above lighter, the equilibrium in unstable and a small perturbation of the interface from the horizontal will grow with time, producing the phenomenon known as Rayleigh-Taylor instability [1]. This instability is a prototype problem for computational studies of multi-phase flows. The problem consists of two layers of fluid initially at rest in the rectangular domain Ω = (−d/2, d/2) × (−2d, 2d), see Figure 8.1. The flow is characterized by the density difference between the two fluids and their effective viscosity.
Figure 8.1: Initial configuration and physical properties of the fluids. The density difference is represented by the Atwood number At = (ρA − ρB )/(ρA + ρB ). The Reynolds number is defined as Re = ρA d3/2 g 1/2 /µ, where d is the reference length, g the gravity acceleration and µ the dynamic viscosity of the fluids (assumed uniform). The growth and evolution of Rayleigh–Taylor instability has been investigated among others by Tryggvason [2] for inviscid incompressible flows, and by Guermond & Quartapelle [3] and Ding et al. [4] for viscous flows. None of these studies has taken into account surface tension. We compare the XFlow results with those of [3] and [4] at At = 0.5 and Re = 3000. The initial position of the perturbed interface is y(x) = −0.1 d cos(2πx/d).
46
Figure 8.2: Vertical position of spike and bubble vs. time. XFlow solution with reference ones.
Multi-phase flows
Comparison of
Computations are carried out on a 200 × 800 grid and the time step is automatically set to 0.000144 s. Free-slip condition is enforced at all walls. The tracking of the interface is done using the marker-and-cell method. Results on the vertical position of the tip of the falling and rising fluid (spike and bubble, respectively) are shown in Figure 8.2. XFlow solution is in good agreement with the reference results [3, 4]. The evolution of the instability is shown in Figure 8.3 at dimensionless times √ t˜ =0, 1, 1.25, 1.5, 1.75, 2, 2.25, 2.5, where t˜ = t g At. Around t˜ = 1.5 the heavy fluid begins to roll up into two counter-rotating vortices (see also Figure 8.4). Later, around t˜ = 2, these two vortices become unstable and a pair of secondary vortices appear at the tails of the roll-ups. The roll-ups and vortices in the heavy fluid spike are due to the Kelvin–Helmholtz instability. The shapes of the fluid interface obtained with XFlow compare well with those of the reference results [3, 4].
47
t˜ = 0
t˜ = 1
t˜ = 1.25
t˜ = 1.5
t˜ = 1.75
t˜ = 2
t˜ = 2.25
t˜ = 2.5
Figure 8.3: Rayleigh-Taylor instability evolution.
48
REFERENCES
Figure 8.4: Velocity field at t˜ = 1.5 s.
References [1] D.H. Sharp. An overview of Rayleigh-Taylor instability. Physica D, 12:3–18, 1984. [2] G. Tryggvason. Numerical simulations of the Rayleigh-Taylor instability. Journal of Computational Physics, 75:253–282, 1988. [3] J.-L. Guermond and L. Quartapelle. A projection FEM for variable density incompressible flows. Journal of Computational Physics, 165:167–188, 2000. [4] H. Ding, P. Spelt, and C. Shu. Diffuse interface model for incompressible two-phase flows with large density ratios. Journal of Computational Physics, 226:2078–2095, 2007.
