CFD Simulations AC2-10
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\"{a}t Bergakademie Freiberg (TUBF), Universit\"{a}t Duisburg-Essen (UDE) and Technische Universit\"{a}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 \ref{tab:simulations}.
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 \ref{tab:simulations}.
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. \cite{Nguyen2016} and Baumann et al. \cite{Baumann2014}, 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 are solved. They are defined as:
with the thermal conductivity and the turbulent Prandtl number (set to 0.9 in the following). A non-density-weighted filtered variable is denoted by , while is the density-weighted filtered variable. This system is complemented by the filtered equation of state for perfect gases, defined as:
with the molecular weight and the universal gas constant . The eddy viscosity () is calculated based on the scale-adaptive simulation (SAS-SST) turbulence model \cite{Menter2010}. In case of insufficient spatial or temporal resolution (for a SRS), it reverts to the SST model \cite{Menter1994} and maintains a valid base of modelling. Its key element is the von K\'arm\`an length scale \cite{Rotta2010} (), introduced into the scale-determining equation. In unstable flows, adjusts the eddy viscosity to a level which allows the formation of a turbulent spectrum \cite{Menter2010,Travin2002,Schaefer2010,Lucius2010}. is defined as
with
and the von Kármán constant .
ANSYS CFX R16.0 is a node-based, conservative and time-implicit finite-volume code \cite{Raw1996,Raithby1979,VanDoormaal1984}. The mass flow is evaluated such that a pressure-velocity coupling is achieved \cite{Rhie1983}. The discrete systems of equations are solved using a coupled algebraic multi-grid method \cite{Raw1996}. To minimize numerical diffusion, a second-order scheme in space \cite{Travin2002} 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:
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 \ref{eq:ofmass}-\ref{eq:ofenergy} and the utilized variables can be found in \cite{Janas2017,JanasDiss2017}.\\
The standard Smagorinsky model \cite{Smagorinsky1963} with a model constant of 0.062 is used to account for the unresolved subgrid stresses $\tau_{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\v{z}i\'{c}, is employed to calculate the pressure \cite{Demirdzic1993,Ferziger2002}. 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 \cite{Nguyen2014} and unacceptable numerical dissipation. A total variation diminishing scheme (TVD-scheme) is utilized to discretize the convective scalar-fluxes, using the Sweby limiter \cite{Sweby1984}.\\
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 \cite{Amsden1989,Torres2006}. The overline characterizes LES filtered and the tilde Favre filtered quantities. Further information regarding the different terms of equations \ref{eq:kivamass}-\ref{eq:kivaenergydef} can be found in \cite{Amsden1989,Torres2006,He2017}.
The internal energy is defined as follows:
The classical Smagorinsky model \cite{Smagorinsky1963} with a Smagorinsky constant set to 0.1 is utilized to model the turbulent kinematic viscosity. The quasi-second-order upwind (QSOU) scheme \cite{Amsden1989} 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.
Computational Domain and Mesh Treatment
As illustrated in Figure \ref{fig:fdomain}, the numerical domain is characterized by temporarily disabled ports for the time periods when the corresponding valves are closed. When the valves reopen, the ports are added to the numerical domain. Due to this disabling and enabling strategy, the flow field of each port has to be reinitialized within each individual operating cycle. In this study the reinitialization is done by one (previously calculated) flow field for each port, to minimize the computational costs.
![]() |
Figure 12: Ansys CFX - Numerical domain during the entire operation cycle. |
Figure \ref{fig:fmesh} illustrates the mesh topology within the combustion chamber, shown on the valve middle plane. The hybrid, and body-conformal mesh consists of tetrahedral, prism and hexahedral elements with a maximum cell size of 1\,mm. Mesh refinements of 0.25\,mm are defined at the intake jet region, the valve stem and the separation edge within the intake port. The boundary layer is resolved by 10 prism layers with a minimum height of 25\,$\mathrm{\mu m}$ next to the wall. From 310 to 36\degree bTDC and from 37 to 314\degree bTDC an extrusion volume is generated to reduce the number of cells. A key grid approach \cite{Buhl2017} is used, while the effects of mesh morphing and interpolation were investigated in a previous study \cite{Hasse2010}. The time step width is set to 0.1 CAD (≈ 2.08×10-5s) for all simulations.
![]() |
![]() |
Figure 13: Ansys CFX - General mesh topology within the combustion chamber. Mesh refinement at the intake jet region, valve stem and intake port separation edge are highlighted by orange box. |
OpenFOAM-2.3.x
The numerical grid is built up of equidistant hexahedral cells, locally refined at the spark plug and the valve seat region (see left of Figure 14. Inside the combustion chamber an average cell size of 0.5\,mm is employed, whereas inside the valve seat region, the cell size is 0.125 to 1\,mm. At the spark plug the cell size is 0.25\,mm and inside the inlet and outlet manifolds, at a distance of the valves, the mesh resolution is reduced to 2 mm. In total, the used grids contain 2.0 million cells at TDC and up to 5.7 million cells at BDC. To close the valves, so called 'curtains' are used, which are internal walls around the valve seat.
![]() |
![]() |
![]() |
Figure 14: OpenFoam - Left: Engine grid on a cut through the valves with opened intake valve. Middle: Top view of the cylinder head. Right: Full computational domain.} |
The mesh motion strategy involves 144 target grids for every 5 CAD for one full cycle, which are created prior to the simulation. Here, 5 CAD correspond to a piston motion of less than 4mm, where the mesh quality stays high. Every 5 CAD the results are mapped to a new grid. The mesh motion is realized by shifting each mesh point by a distance resulting from the multiplication of the time step size and deformation velocity, which is obtained by a Laplace equation for a grid velocity field. Further, all convective fluxes are calculated relative to the moving grid \cite{Demirdzic1988}. Further Information can be found in \cite{JanasDiss2017}.
Contributed by: Carl Philip Ding,Rene Honza, Elias Baum, Andreas Dreizler — Fachgebiet Reaktive Strömungen und Messtechnik (RSM),Technische Universität Darmstadt, Germany
Contributed by: Brian Peterson — School of Engineering, University of Edinburgh, Scotland UK
Contributed by: Chao He , Wibke Leudesdorff, Guido Kuenne, Benjamin Böhm, Amsini Sadiki, Johannes Janicka — Fachgebiet Energie und Kraftwerkstechnik (EKT), Technische Universität Darmstadt, Germany
Contributed by: Peter Janas, Andreas Kempf — Institut für Verbrennung und Gasdynamik (IVG), Lehrstuhl für Fluiddynamik, Universität Duisburg-Essen, Germany
Contributed by: Stefan Buhl, Christian Hasse — Fachgebiet Simulation reaktiver Thermo-Fluid Systeme (STFS), Technische Universität Darmstadt, Germany; former: Professur Numerische Thermofluiddynamik (NTFD), Technische Universität Bergakademie Freiberg, Germany
© copyright ERCOFTAC 2018