Internal combustion engine flows for motored operation

Application Challenge AC2-10 © copyright ERCOFTAC 2024

CFD Simulations

Introduction

The TU Darmstadt engine has been investigated by three different groups using LES (Large Eddy Simulation) and hybrid URANS (unsteady Reynolds-averaged Navier-Stokes)/LES. These research groups are located at Technische Universität Bergakademie Freiberg (TUBF), Universität Duisburg-Essen (UDE) and Technische Universität Darmstadt (TUD). In the follow, the different approaches will be presented, including general information about the code and the physical modelling, computational domain and mesh treatment, initial and boundary conditions and computational requirements. Additional and detailed information can be found in the scientific papers listed in Table 3.

Overview of simulation

Test data for the motored case is available at three different operating points. A test case considering an engine speed of 800 RPM is investigated in the numerical studies, which are summarized in Table 3.

**Table 3:** Summary of simulations carried out at 800 RPM
Group	Approach	CFD Code / Publication
TUBF	URANS / LES	Ansys CFX R16.0 / Buhl et al. [9, 10]
UDE	LES	OpenFoam-2.3.x / Janas et al. [21]
TUD	LES	KIVA-4mpi / He et al. [20]

Further, simulations were carried out by Nguyen et al. [31] and Baumann et al. [6], which are not part of this comparison.

CFD codes

In the following, the governing equations, turbulence modelling and solver are described for the three different CFD codes.

Ansys CFX R16.0

For this compressible case, the filtered governing equations for mass, momentum and total enthalpy ${(H=h(T,p)+0.5u_{j}^{2})}$ are solved. They are defined as:

${{\dfrac {\partial {\overline {\rho }}}{\partial t}}+{\dfrac {\partial \left({\overline {\rho }}{\widetilde {u}}_{j}\right)}{\partial x_{j}}}=0}$	${\text{(4.1)}}$
${{\dfrac {\partial \left({\overline {\rho }}{\widetilde {u}}_{i}\right)}{\partial t}}+{\dfrac {\partial \left({\overline {\rho }}{\widetilde {u}}_{i}{\widetilde {u}}_{j}\right)}{\partial x_{j}}}=-{\dfrac {\partial {\overline {p}}}{\partial x_{i}}}+{\dfrac {\partial }{\partial x_{j}}}\left(\left(\mu +\mu _{t}\right){\dfrac {\partial {\widetilde {u}}_{i}}{\partial x_{j}}}\right)}$	${\text{(4.2)}}$
${{\dfrac {\partial \left({\overline {\rho }}{\widetilde {H}}\right)}{\partial t}}-{\dfrac {\partial {\overline {p}}}{\partial t}}+{\dfrac {\partial \left({\overline {\rho }}{\widetilde {u}}_{j}{\widetilde {H}}\right)}{\partial x_{j}}}={\dfrac {\partial }{\partial x_{j}}}\left(\lambda {\dfrac {\partial {\widetilde {T}}}{\partial x_{j}}}+{\dfrac {\mu _{t}}{Pr_{t}}}{\dfrac {\partial {\widetilde {h}}}{\partial x_{j}}}\right)}$	${\text{(4.3)}}$

with the thermal conductivity ${\lambda }$ and the turbulent Prandtl number ${Pr_{t}}$ (set to 0.9 in the following). A non-density-weighted filtered variable is denoted by ${\overline {\Phi }}$ , while ${\widetilde {\Phi }}$ is the density-weighted filtered variable. This system is complemented by the filtered equation of state for perfect gases, defined as:

{{\overline {p}}={\overline {\rho }}{\dfrac {R}{W}}{\widetilde {T}}}

{\text{(4.4)}}

with the molecular weight ${W}$ and the universal gas constant ${R}$ . The eddy viscosity ( ${\mu _{t}}$ ) is calculated based on the scale-adaptive simulation (SAS-SST) turbulence model [27]. In case of insufficient spatial or temporal resolution (for a SRS), it reverts to the SST model [26] and maintains a valid base of modelling. Its key element is the von Kármàn length scale [38]. ( ${L_{vK}}$ ), introduced into the scale-determining ${\omega }$ equation. In unstable flows, ${L_{vK}}$ adjusts the eddy viscosity to a level which allows the formation of a turbulent spectrum [27, 45, 39, 25]. ${L_{vK}}$ is defined as

{L_{vK}=\kappa {\dfrac {\sqrt {2{\widetilde {S}}_{ij}{\widetilde {S}}_{ij}}}{{\widetilde {u}}''}}}

{\text{(4.5)}}

with

{{\widetilde {S}}_{ij}={\dfrac {1}{2}}\left({\dfrac {\partial {\widetilde {u}}_{i}}{\partial x_{j}}}+{\dfrac {\partial {\widetilde {u}}_{j}}{\partial x_{i}}}\right){\text{,}}\quad {\widetilde {u}}''={\sqrt {{\dfrac {\partial ^{2}{\widetilde {u}}_{i}}{\partial x_{k}^{2}}}{\dfrac {\partial ^{2}{\widetilde {u}}_{i}}{\partial x_{j}^{2}}}}}}

{\text{(4.6)}}

and the von Kármán constant ${\kappa }$ .

ANSYS CFX R16.0 is a node-based, conservative and time-implicit finite-volume code [36, 35, 46]. The mass flow is evaluated such that a pressure-velocity coupling is achieved [37]. The discrete systems of equations are solved using a coupled algebraic multi-grid method [36]. To minimize numerical diffusion, a second-order scheme in space [45] and a second-order backward scheme in time are used.

OpenFOAM-2.3.x

The following governing equations for mass, momentum and absolute internal energy are solved on a moving grid:

${{\dfrac {\partial {\overline {\rho }}}{\partial t}}+{\dfrac {\partial \left({\overline {\rho }}{\widetilde {u}}_{j}\right)}{\partial x_{j}}}=0}$	${\text{(4.7)}}$
${{\dfrac {\partial \left({\overline {\rho }}{\widetilde {u}}_{i}\right)}{\partial t}}+{\dfrac {\partial \left({\overline {\rho }}{\widetilde {u}}_{i}{\widetilde {u}}_{j}\right)}{\partial x_{j}}}={\dfrac {\partial {\overline {\tau }}_{ij}}{\partial x_{j}}}+{\dfrac {\partial \tau _{ij}^{sgs}}{\partial x_{j}}}-{\dfrac {\partial {\overline {p}}}{\partial x_{i}}}}$	${\text{(4.8)}}$
${{\dfrac {\partial \left({\overline {\rho }}{\widetilde {e}}\right)}{\partial t}}+{\dfrac {\partial \left({\overline {\rho }}{\widetilde {u}}_{j}{\widetilde {e}}\right)}{\partial x_{j}}}+{\dfrac {\partial \left({\overline {\rho }}{\widetilde {K}}\right)}{\partial t}}+{\dfrac {\partial \left({\overline {\rho }}{\widetilde {K}}{\widetilde {u}}_{j}\right)}{\partial x_{j}}}={\dfrac {\partial }{\partial x_{j}}}\left(\alpha _{eff}{\dfrac {\partial {\widetilde {e}}}{\partial x_{j}}}\right)-{\dfrac {\partial }{\partial x_{j}}}\left({\overline {p}}{\widetilde {u}}_{j}\right)}$	${\text{(4.9)}}$

Here, the overline denotes LES filtered and the tilde Favre filtered quantities. Further, ${{\overline {\tau }}_{ij}}$ denotes the viscous stress tensor. Detailed information regarding equations 4.7–4.9 and the utilized variables can be found in [21, 22].

The standard Smagorinsky model [40] with a model constant of 0.062 is used to account for the unresolved subgrid stresses τ_sgs and the Sutherland law is used to calculate the molecular viscosity. The pressure-velocity-density coupling for flows at arbitrary Mach-number, proposed by Demirdžić, is employed to calculate the pressure [14, 16]. Since the time step size is not limited according to the speed of sound, this approach allows bigger time steps, using an implicit second-order scheme. The momentum equation is solved using a central differencing scheme (CDS). In regions with Mach-numbers greater than 0.3 a TVD-scheme is used to avoid numerical problems [29] and unacceptable numerical dissipation. A total variation diminishing scheme (TVD-scheme) is utilized to discretize the convective scalar-fluxes, using the Sweby limiter [43].

The time step size employed is calculated from the maximum CFL number of 2. The smallest time step occurs when the valves are opening and closing, since the smallest mesh resolution can be found inside the valve gap and the highest velocities. This leads to time steps smaller than one micro second. After the valves are closed, the time step size increases to almost 50 micro seconds during the compression and expansion stroke.

KIVA-4mpi

The Favre-filtered governing equations for mass, momentum and internal energy are discretized on a moving mesh according to the Arbitrary Lagrangian-Eulerian (ALE) approach [2, 44]. The overline characterizes LES filtered and the tilde Favre filtered quantities. Further information regarding the different terms of equations 4.10–4.13 can be found in [2, 44, 20].

${{\dfrac {\partial {\overline {\rho }}}{\partial t}}+{\dfrac {\partial \left({\overline {\rho }}{\widetilde {u}}_{j}\right)}{\partial x_{j}}}=0}$	${\text{(4.10)}}$

${{\dfrac {\partial \left({\overline {\rho }}{\widetilde {u}}_{i}\right)}{\partial t}}+{\dfrac {\partial \left({\overline {\rho }}{\widetilde {u}}_{i}{\widetilde {u}}_{j}\right)}{\partial x_{j}}}}=-{\dfrac {\partial {\overline {P}}}{\partial x_{i}}}+{\dfrac {\partial }{\partial x_{j}}}\left(\left({\overline {\mu }}+\mu _{t}\right)\left({\dfrac {\partial {\widetilde {u}}_{i}}{\partial x_{j}}}+{\dfrac {\partial {\widetilde {u}}_{j}}{\partial x_{i}}}-{\dfrac {2}{3}}{\dfrac {\partial {\widetilde {u}}_{k}}{\partial x_{k}}}\delta _{ij}\right)\right)+{\overline {\rho }}g_{i}$	${\text{(4.11)}}$
${{\text{with}}\quad {\overline {P}}={\overline {p}}+{\dfrac {1}{3}}{\overline {\rho }}\tau _{kk}^{sgs}}$

${{\dfrac {\partial \left({\overline {\rho }}{\widetilde {e}}\right)}{\partial t}}+{\dfrac {\partial \left({\overline {\rho }}{\widetilde {u}}_{j}{\widetilde {e}}\right)}{\partial x_{j}}}=-{\overline {p}}{\dfrac {\partial {\widetilde {u}}_{j}}{\partial x_{j}}}+{\widetilde {\tau }}_{ij}{\dfrac {\partial {\widetilde {u}}_{i}}{\partial x_{j}}}+{\dfrac {\partial }{\partial x_{j}}}\left({\dfrac {\left({\overline {\mu }}+{\dfrac {c_{v}}{c_{p}}}\mu _{t}\right)c_{p}}{Pr}}{\dfrac {\partial {\widetilde {T}}}{\partial x_{j}}}\right)}$	${\text{(4.12)}}$

The internal energy ${e}$ is defined as follows:

{e=h-{\dfrac {p}{\rho }}=e_{0}+\int _{T_{0}}^{T}c_{v}dT-{\dfrac {p}{\rho }}\quad {\text{, with}}\quad h=h_{0}+\int _{T_{0}}^{T}c_{p}dT\quad {\text{and}}\quad c_{v}=\left.{\dfrac {\partial e}{\partial T}}\right|_{v}}

{\text{(4.13)}}

The classical Smagorinsky model [40] with a Smagorinsky constant set to 0.1 is utilized to model the turbulent kinematic viscosity. The quasi-second-order upwind (QSOU) scheme [2] is used as the spatial and the first-order implicit Euler as the temporal discretization scheme. A staggered arrangement is employed to store the variables and the time step size is adjusted during the computation according to accuracy requirements and to retain a CFL number smaller than one.