1.3 Fundamentals of Computational Fluid Dynamics
1.3.1 Continuous gas phase
Any numerical simulation has to start with a suitable, mathematical description of the problem in adapted coordinates. In a next step, the problem is discretised: the consid-ered spatial domain is divided in a finite number of small volumes and the continuous equations are replaced by algebraic ones. The system of the latter is then solved by an efficient solver algorithm. Finally, the vast amount of resulting information has to be filtered and studied by the user with the help of a post-processor, which is not addressed further.
1.3.1.1 Mathematical description
Together with suitable initial and boundary conditions, conservation equations from continuum mechanics for mass, momentum and energy describe the gaseous flow mathe-matically. In the following, they are given in differential form and cartesian coordinates, see [19], [76]. The index summation according to Einstein is used1. ρ signifies the density, µ the dynamic viscosity and λ the thermal conductivity of the fluid. The velocity component in direction i is denoted asvi.
1If an indexioccurs twice in any term, the term is summed with the indexi taking values from1 to3.
∂(ρ·vi)
∂t +∂(ρ·vi·vj)
∂xj = ∂τij
∂xj − ∂p
∂xi +ρ·fV,i+Sp,i (i= 1,2,3), (1.3) where the components τij of the stress tensor matrix for a Newtonian fluid are defined as
τij =µ· ∂vi
∂xj +∂vj
∂xi −2 3 · ∂vk
∂xk ·δij
. (1.4)
The term fV,i stands for all volume forces, e.g. gravity, centrifugal, Coriolis or buoyancy forces, arising in the chosen reference system. The source termSp,irefers mainly to the drag force in coupled two-phase flows. δij denotes the Kronecker symbol.
• Energy conservation:
∂(ρ·h)
∂t +∂(ρ·vi·h)
∂xi =τij · ∂vi
∂xj +∂p
∂t +vi· ∂p
∂xi − ∂q˙i
∂xi +Sh, (1.5) where h=e+p/ρspecifies the enthalpy with ethe internal energy, both defined per mass unit. q˙i =−λ·∂x∂T
i (Fourier’s law) describes the heat flux in direction i due to a temperature gradient.
Additional energy fluxes e.g. due to a concentration gradient as well as thermal radiation are omitted in the equation, cf. [7] for further details.
The source term Sh is coupled to the vapour mass sources due to liquid drop evaporation.
Besides the velocity components vi and the enthalpy h, the five Equations 1.2-1.5 contain the unknown density, pressure and temperature fields: ρ = ρ(x1, x2, x3, t), p=p(x1, x2, x3, t) and T =T(x1, x2, x3, t). The system of equations must therefore be completed by state equations
p=p(ρ, T) and h=h(ρ, T), (1.6) which describe the thermodynamic properties of the fluid. For an ideal gas, they are given as p=ρ· R ·T and h=cP ·T with R the fluid-specific gas constant and cP the specific heat capacity per mass unit for constant pressure.
So far, the equations have been presented for a single, pure fluid. In case of a multi-component gaseous phase with varying composition, e.g. due to evaporation or chemical reactions (the latter being not considered in this work), a continuity equation has to be solved for every species n = 1. . . N:
∂ρn
∂t +∂(ρn·vn,i)
∂xi
=−∂jn,i
∂xi
+Sn . (1.7)
Sndescribes the production of the species. In comparison to Equation 1.2 which is still applicable for the overall mixture and can be used instead of one of the n equations for the individual species, diffusive fluxesjn,i appear in these equations. They describe the motion of a species relative to the mean fluid motion and comprise mass diffusion due to mechanical driving forces (Fick’s law), i.e. due to a concentration gradient, a pressure gradient or possibly due to an external force acting unequally on the various species (e.g. an electromagnetic force), as well as mass diffusion due to a temperature gradient (Soret effect). The contributions due to a concentration and a temperature gradient are often the most important ones.
Momentum and energy equation, cf. Equations 1.3 and 1.5, can be employed further in case of a multi-component fluid, if values, like the velocity, are considered as averages over all species. However, analogous to the diffusive mass fluxes in Equation 1.7, ad-ditional energy fluxes besides the Fourier term appear from the interdiffusion between the different species due to different mechanical driving forces (Dufour effect), cf. [7]
for further details.
1.3.1.2 Reynolds averaging for turbulent flows
The flows addressed in this work and in most industrial applications show a very com-plex and irregular behaviour in space and time: they are turbulent. To resolve the as-sociated high-frequency fluctuations, which can consist of many scales, in simulations, the spatial and temporal resolution have to capture the smallest vortex and fluctuation as well as the largest. Such a direct numerical simulation (DNS) is much too expensive in most engineering applications, see [25] and references therein.
In large-eddy simulations (LES), only the large scale eddies are resolved. These trans-port conserved quantities very effectively and comprise most of the turbulent energy.
The latter is passed through all length scales in a cascade of decaying eddies and finally dissipates into heat at the smallest scale. Small eddies are not resolved but modelled.
Yet, LES is still very demanding with respect to computational time and capacity.
Therefore, a statistical model is usually chosen with all turbulent fluctuations being modelled. In many applications and in the approach considered in this work, the re-sulting average quantities describe the respective flows sufficiently.
According to Reynolds, the instantaneous value of a quantity Φ(x1, x2, x3, t) (e.g. ve-locity components, pressure, temperature, density etc.) is divided in an averaged and a stochastic fluctuation term. The latter is high-frequent in comparison to the stationary or slowly changing averaged term:
Φ(x1, x2, x3, t) = ¯Φ(x1, x2, x3, t) + Φ0(x1, x2, x3, t), (1.8) where the ensemble average2 is defined as
Φ(x¯ 1, x2, x3, t) = lim
n→∞
1 n
n
X
k=1
Φ(k)(x1, x2, x3, t). (1.9) For statistically stationary flows this averaging can be replaced by a simpler time-averaging:
Φ(x¯ 1, x2, x3, t) = ¯Φ(x1, x2, x3) = lim
t∗→∞
1 t∗
Z t+t∗/2 t−t∗/2
Φ(x1, x2, x3, t0)dt0, (1.10)
2The flow is realisedntimes and every time, the value ofΦis noted at the same time and position.
∂t ∂xj ∂xj V,i p,i ∂xj An incompressible flow is volume-conserving, i.e. ∂ρ∂t +vi·∂x∂ρ
i = 0. In combination with the continuity Equation 1.2 this equals the condition ∂x∂vi
i = 0.
To solve for the statistical means, the new unknownsρ·v0i·vj0, which are called Reynolds stresses and which arise due to the non-linearity of the Navier-Stokes equations, have to be approximated. For the transport equation of any other scalar quantity, e.g. the energy, additional turbulent scalar fluxes, like ρ·vi0·e, appear and have to be treated, too.
This so-called turbulence modelling can be done with zero-equation models, where the stresses are modelled by algebraic equations. This is not very accurate and one- or two-equation models are preferable. These solve one or two additional differential two-equations for the turbulent kinetic energykand possibly the dissipation rateε. Another possibility are Reynolds stress models, where transport equations are solved for all stresses and the dissipation rate. These models promise better accuracy but are also more expensive and in general less robust.
Throughout this work, the standard two-equation modelsk-ε, see [36], and Shear-Stress-Transport (SST), cf. [3], [41], are considered. Both belong to the so-called eddy viscosity turbulence models: the effects of turbulent fluctuations in the flow are described by the introduction of a turbulent viscosityµt. This is done in analogy to laminar flows, where energy dissipation is mainly due to viscous forces. Yet, µt is not a fluid property like µ but dependent on the flow.
k-ε-model. Using the Boussinesq approximation, see [8], the Reynolds stresses are modelled in a form analogous to the components of the friction part of the stress tensor matrix for Newtonian fluids:
ρ·v0i·vj0 =−µt· ∂v¯i
∂x¯j +∂v¯j
∂x¯i
+ 2
3·ρ·δij ·k . (1.12) k = 0.5·v0i·vi0 defines the turbulent kinetic energy.
The problem is not closed yet as µt and k are unknown. The ansatz µt=Cµ·ρ· k2
ε (1.13)
and the definition of the dissipation rate ε= µ
ρ · ∂v0i
∂xj · ∂vi0
∂xj (1.14)
lead to the transport equations for k and ε, which are used for closure:
∂(ρ·k)
∂t + ∂
∂xi (ρ·¯vi·k) = ∂
∂xi
µ+ µt σk
· ∂k
∂xi +Gk−ρ·ε , (1.15)
∂(ρ·ε)
∂t + ∂
∂xi (ρ·v¯i·ε) =
µ+ µt σε
· ∂ε
∂xi + ε
k ·(Cε1·Gk−Cε2·ρ·ε) . The values of the empirical constants are commonly set to Cµ = 0.09, σk = 1.0, σε = 1.33, Cε1 = 1.44and Cε2 = 1.92. Gk is the production term of turbulent kinetic energy and given as
Gk=−ρ·vi0·vj0 · ∂v¯i
∂x¯j ≈µt· ∂v¯i
∂xj + ∂v¯j
∂xi
· ∂v¯i
∂x¯j . (1.16) Source terms due to two-phase interactions do not occur, cf. Section 1.3.2.2.
The k-ε-model has proved to be very robust and efficient in many applications. The underlying assumptions are valid in highly turbulent flows with isotropic turbulence (due tokandε being scalars). The latter is no longer fulfilled in the proximity of walls, i.e. in the viscous sublayer. The turbulent transport equations lose their validity there and the near-wall turbulence has to be modelled separately.
If the computational grid, see the following subsection, is fine enough to resolve the laminar sublayer, so-called “low Reynolds” models are to be used. Yet, the thickness of the viscous boundary layer scales with the Reynolds number, hl∼Re−0.5. For large values of Re, the first grid point often lies already in the fully developed, turbulent sublayer. There, the well-known logarithmic wall law can be used to describe the velocity profile:
v+ = vT
vτ = 1
κ ·ln(y+) +C . (1.17)
vT is the known velocity component tangential to the wall in a distance∆yfrom the wall, which is presumed in the already fully turbulent boundary layer. vτ = (τw/ρ)0.5signifies the friction velocity with τw the wall shear stress. κ ≈ 0.41 denotes the empirical von Kármán constant. y+ is the non-dimensional distance from the wall defined as y+ =ρ·vτ ·∆y/µ. C is a constant depending on the wall roughness.
For further details on the implementation of near-wall treatment, especially in CFX, cf. [3].
SST-model. Another two-equation turbulence model is the k-ω-model, see [91]. It is very similar to thek-ε-model but uses the turbulent frequency ω instead of the eddy dissipation rate ε. The turbulent viscosity µt is linked to k and ω via:
µt=ρ· k
ω. (1.18)
The transport equations to close the problem become:
∂(ρ·k)
∂t + ∂
∂xi(ρ·v¯i·k) = ∂
∂xi
µ+ µt
σk
· ∂k
∂xi +Gk−βk·ρ·k·ω , (1.19)
∂(ρ·ω)
∂t + ∂
∂xi (ρ·v¯i·ω) =
µ+ µt
σω
· ∂ω
∂xi +αω· ω
k ·Gk−βω·ρ·ω2. The empirical constants are usually given as βk = 0.09, σk = 2.0, σω = 2.0, αω = 5/9 andβω = 0.075. The turbulent production termGkis identical to that in thek-ε-model.
1.3.1.3 Discretisation
In the previous two subsections, the coupled basic equations which describe the gaseous phase have been presented shortly. Analytical solutions for these equations in combina-tion with appropriate initial and boundary condicombina-tions are only available for very simple flows under ideal conditions. Usually a numerical solution has to be found which re-quires a discretisation of the problem in space and time.
First, the considered domain has to be split in a finite number of small control volumes, see Figure 1.4. ANSYS Icem CFD is used for mesh generation in this work. Structured, regular grids are very efficient in numerical algorithms yet they should be oriented to the flow to prevent systematic errors, see [22]. For complex geometries such as a combustion chamber and for automatic grid generation unstructured grids are more flexible, cf. [56].
Therefore these are used for in-cylinder calculations at BMW and the meshes in this thesis are set up accordingly with tetrahedron cells. The solution variables and fluid properties are stored at the element nodes. To account for a better resolution of near-wall flow, flat prism layers are created on the surfaces.
Per default the maximal edge length of the tetrahedra and the total height of the prism layers is chosen as 1mm in the simulations carried out.
Figure 1.4: Definition of an exemplary two-dimensional control volume in CFX.
Second, the partial differential equations as well as the boundary conditions have to be discretised into a manageable, algebraic system. For that purpose, different math-ematical descriptions are well-known: the volume, the element, the finite-difference as well as the spectral method, see [25], [67], [85].
In CFX, the finite-volume method is implemented, which has the advantage to be con-servative in mass, momentum and energy in a global as well as in a discrete sense,
i.e. for every control volume. Its main idea can be summarised in looking at the general form of a transport equation:
∂(ρ·Φ)
∂t
| {z }
local change in time
+∂(ρ·vi·Φ)
∂xi
| {z }
convection
=
∂
ΓΦ· ∂x∂Φ
i
∂xi
| {z }
diffusion
+ SΦ
|{z}
sources, sinks
. (1.20)
All governing Equations 1.2-1.5 and 1.15 can be written in this form. The change of the respective transport quantityΦin time is given by the convective and diffuse transport and its sources and sinks. The diffusion coefficient Γφ and the source terms SΦ can be easily determined by comparing Equation 1.20 with the respective governing equation.
In the finite-volume method, the transport equation is considered in integral form for every control volume. Applying Gauss’ divergence theorem, it becomes:
Z
A
ρ·vi·Φ−ΓΦ· ∂Φ
∂xi
| {z }
convective and diffusive flux
d−→ A =
Z
V
−∂(ρ·Φ)
∂t +SΦ
dV . (1.21)
The terms of this equation are then discretised, i.e. approximated by discrete values:
for the convective flows, the values of Φat the surface border of a control volume have to be approximated. The upwind differencing scheme is a very common, first-order scheme, where the values are simply approximated by the values of Φ at the node in upstream direction. In contrast to higher order schemes it is bounded and does not yield oscillatory solutions in case of large gradients of the transport quantity. Yet, it is numerically diffusive, i.e. the leading error resembles a diffusive flux. More accurate but less robust is the second-order central differencing scheme which uses an interpolation between neighboured nodes to get the value ofΦ at the border of a control volume. In this work, the “high-resolution” scheme in CFX is chosen. It is a flux-blending technique, where a blending factor is calculated based on the gradients in the flow. This factor is used to interpolate between the mentioned first-order and second-order schemes in such a way that the discretisation is as close to second order as possible without local oscillations.
The diffusive fluxes are less problematic and always discretised second-order using the central differencing scheme.
For the volume integrals on the right-hand side of Equation 1.21, the midpoint rule is applied. The time derivate is finally approximated by an implicit second-order backward Euler approach.
1.3.1.4 Solving the system of equations
From discretisation an algebraic nonlinear system results. It has to be solved numeri-cally at the representative positions of every control volume and for a finite number of time values. To do so, it is linearised and assembled into a solution matrix:
A−→ Φ =−→
b . (1.22)
CFX uses a multigrid accelerated solver with an Incomplete Lower Upper factorization technique for the solution matrix, see [3], [67]. The velocity componentsvi (i= 1,2,3) and the densityρare treated as a single algebraic system first. After this hydrodynamic system, volume fractions, additional variables, energy, turbulence and finally the liquid
stopped. A target for the global balances, called conservation target, can additionally be set. In a transient calculation the solver then proceeds to the next timestep.