Unsteady Engine Analysis with Moving Mesh in OpenFOAM®
By Prof. Federico Piscaglia
The effect of turbulence on the thermodynamic efficiency, brake power and emissions is significant among all the physical phenomena that occur inside an engine cylinder because its influence extends from volumetric efficiency to air/fuel mixing, combustion and heat transfer.
Historically, turbulent flows have been simulated mainly by models based on Reynolds Averaged Navier-Stokes equations (RANS), either in their original version or in the unsteady formulation (URANS) for slowly-varying flows . In internal combustion (IC) engine simulation, URANS approaches have proved to provide very good predictions of phase-averaged flow fields: macroscopic features of charge motion – like swirl and tumble vortexes – can be estimated with good accuracy . On the other hand, most of the time-varying quantities characterizing in-cylinder flows cannot be resolved by a model based on implicit time- or ensemble-averaging methods like URANS. In this context, the research work described here aims to investigate how small-scale turbulence, cycle-to-cycle variability (CCV) and in-cycle evolution of three-dimensional structures (jets and vortices) in IC engines can be simulated only by a time-resolved (rather than time-averaged) approach like Large Eddy Simulation (LES). The computational tool used in this paper is the open-source computational fluid dynamics (CFD) software OpenFOAM® released by ESI-OpenCFD, which has been extended to perform LES simulation with topologically changing grids. Examples of application of the proposed methodology are the simulation of two- and four-stroke engines with moving valves, the non-uniform corrosion of solid propellants for aerospace applications, and the simulation of wind turbines.
In this paper, two different test cases related to the simulation of IC engines are reported:
The two cases were used to study the interaction of the LES models implemented in OpenFOAM [3, 5, 6, 7, 8, 9] with different dynamic mesh motion strategies. In particular, the two cases correspond to two different levels of complexity for the mesh motion strategy that was used to move the piston and the valves.
Case 1: a flat-top cylinder head with a fixed, axis-centered valve
In this first case, the geometry experimentally investigated in  has been modeled. The piston had a bore diameter of 75 mm and stroke was 60 mm. Clearance height was 30 mm from cylinder head, thus the geometric compression ratio is 3. The poppet valve was coaxial with respect to the piston axis and it remained at a fixed height through the whole engine cycle. Piston motion was purely harmonic with a frequency of 200 rpm; as a consequence, the engine Reynolds number (considering air as a working fluid) was
Re = ρ Up D ⁄ μ = 2000
Figure 1: Top) Geometry of the experimental apparatus used by Morse et al. . The same configuration has been modeled for the present paper, with the addition of an upper plenum (not shown). Bottom) Cut view of the finite volume grid used for the simulations. The whole mesh had about 4.8 million hexahedral elements, including the plenum (not shown).
The geometry has been discretized with a pure structured hexahedral mesh (in Fig. 1 bottom). The cylinder inlet section was connected to a plenum (not represented) the volume of which was about 15 times the cylinder volume. The purpose of the plenum was to avoid reflections from the boundaries and thus improve the stability of the simulation. The mesh was generated using Pointwise. To have the best possible control on the mesh quality around the valve, a structured surface mesh was manually created. From this, a boundary layer mesh with hexahedral cells was extruded around the valve. The grid was made by fully structured hexahedral elements. In order to capture the jet flow instability out of the valve seat an oriented mesh was used. The total number of elements was approximately 4.8 million at the Top Dead Center, with a minimum cell thickness of about 0.02 mm in the valve seat.
For this case, dynamic mesh handling has been based on a point motion strategy without topological changes. Mesh points have been divided into three sets according to the motion law to which they were subjected (see Fig. 1 bottom). All points above the cylinder head were fixed since there were no moving boundaries there; as a consequence, they belonged to the set named “static mesh” and their displacement velocity was zero everywhere. Similarly, points located below the cylinder head (“head points” in Fig. 1 bottom) were not moving to preserve the cell quality in a critical region like the interface between valve seat and cylinder. Finally, points belonging to the cylinder region (with the exception of the “head points” just mentioned) did move axially to account for the deforming geometry. The cylinder axis was aligned to the z axis of the global reference frame (Fig. 1); the following relations are therefore written accordingly. Given Up as the piston velocity
Up(t) = S ⁄ 2 sin(ωt) ic
where S is the piston stroke and ic the unit vector parallel to cylinder axis, the point velocity up was calculated as:
up(xp) = Up ⋅ ( zmax − zp ) ⁄ ( zmax − zpiston )
In the equation above, xp = (xp,yp,zp) is the point position, zmax is the z-coordinate of the farthest moving point from the cylinder and zpiston is the z-coordinate of the piston.
A compressible solver for compressible flows with dynamic meshes, developed by the authors  in OpenFOAM, was used. Second-order backward differencing scheme was applied for discretizing the temporal derivatives, whereas momentum convection was performed with the Linear-Upwind Stabilized Transport (LUST) scheme, a low-dissipation method specifically developed for LES . For the remaining differential terms, pure second-order differencing schemes were used, with the exception of energy, for which an upwind-biased method was employed for stability. Starting from a quiescent velocity field, ten cycles were simulated in total. A complex flow field consisting of a number of interacting large-scale coherent structures and a large number of small vortices was obtained (shown in Fig. 2).
Ten full engine cycles were simulated; the last eight were used for post-processing to calculate the ensemble-and azimuthally-averages and fluctuations of the velocity. Very good agreement with the experimentally-determined means and fluctuations at all reported measurement locations at 36°, 90°, 144° and 270° crank-angle degree (CA) was observed. As widely discussed in , LES could reproduce the measurements as reported exemplarily in Fig. 3 at 36° CA over a plane z = 20 mm below the cylinder head.
In Fig. 4 the streamlines stemming from simulation results (on the left) are combined with those obtained from experimental data by Morse et al.  (on the right). The purpose of this series of images is to provide a qualitative overview of the mean velocity field inside the cylinder. More details about the mesh motion strategy adopted and the results achieved on this test case are shown in .
Figure 3: Profiles of mean axial velocity (left) and axial RMS fluctuations (right) for CA= 36° ATDC (after top dead center), at different distances from the cylinder head (conventionally z = 0). — simulations; ■ experiments .
The spatial distribution of the mean turbulent kinetic energy during the intake stroke is plotted in Fig. 5. With increasing piston speed, kmax increases and it reaches its maximum value around 90° CA in the proximity of the outer shear layer of the incoming jet. During the second half of the intake stroke (from 90° to 144° CA), the incoming jet impinges directly on the cylinder wall. Turbulent kinetic energy kmax dissipates and the maximum value of the axial component of the stress tensor shifts close to the cylinder wall; higher values in the shear stress v′r v′θ are connected to an increased anisotropy level with higher piston speeds.
The most promising aspect of LES is the capability of predicting unsteady flow phenomena. Cyclic Combustion Variability (CCV) is usually detected and measured on the basis of some quantities obtained by an integration over the volume, that are examined as they evolve through the different engine cycles. The parameter that has been chosen to monitor CCV is the specific kinetic energy:
〈 E(t) 〉 = ∫ ∫ ∫ ½ ũ(x,t) ũ(x,t) dx dy dz
In Fig. 6a, the normalized value of the integrated kinetic energy E, as defined in the equation above, is plotted against global crank angle for the cylinder region. The quasi-periodic trend is clearly visible; cycle-to-cycle variation is present, and it is especially apparent at the local maxima, when θ ≈ 120° + k 360° . In particular, the difference in absolute value between the highest peak (which belongs to cycle n. 7) and the lowest one (cycle n. 5) amounts to nearly 12 percent. A more detailed discussion can be found in [3, 13].
Figure 6: a) Volume integral of resolved kinetic energy 〈 E(θ)〉 versus global crank-angle; b) Maximum per-cycle value of 〈 E(θ)〉: c) Cycle-to-Cycle Variability of 〈 E(θ)〉: versus engine crank angle: — simulations ▬ average over eight engine cycles (3 to 10); d) Variance of 〈 E(θ)〉: calculated over eight engine cycles (3 to 10).
Case 2: Two-valve head with simple intake and exhaust port/runner geometries and a pancake-shape combustion chamber
The second part of the presented research includes the simulation of IC
engines by an extension of the dynamic mesh handling developed in
OpenFOAM, consisting of automatic mesh motion algorithms based on
topological changes [9, 11]. Some topology modifiers originally
available in the official distribution of OpenFOAM have been revised to
enhance their operation: the
Figure 7: Example of the strategy to generate the engine mesh for the simulation of a transparent combustion engine  with canted valves using Pointwise. Mesh motion is based on sliding non-conformal grids as described in [9,11].
Figure 8: Top) During layer addition/removal, newly inserted faces are dynamically included into the face zone, preserving mesh topology; layer removal is automatically deactivated when only one layer of cells exists between the master face zone and a mesh boundary; Bottom) non-conformal interface between the spark plug and the cylinder mesh.
In this case, layers of hexahedra were grown into the fluid domain from duct surfaces with a smooth transition to tetrahedra in the interior. The mesh in the valve seat was fully hexahedral, while non-conformal interfaces were located between the valve and the cylinder mesh and between the cylinder mesh and the spark plug. Automatic dynamic addition/removal of cell layers was applied at run-time during the simulation, both on the piston surface and on the valves. Mesh size was about 30 million cells at the top dead center. For this work, Pointwise has proved to be a very helpful tool: the software is able to extract internal surfaces from the volume mesh and to export them in stereolithography (STL) format. Also, surface grids can be projected onto prescribed database geometries. These features are very important to facilitate a mesh motion technique based on a non-conformal interface. Some details about the mesh features are reported in Fig. 7-9. Results achieved on preliminary RANS simulation are very accurate and the compressible solver developed  shows very good performance across non-conformal interfaces. A detailed check of continuity and interpolation errors has been performed; also, no problem of dissipation or generation of spurious solution artifacts has been detected. The flow solver and the mesh motion solver developed in OpenFOAM proved to have very good scalability.
The implemented algorithms and methodologies and the first engine case presented in this document have been selected by the CINECA computing center (Bologna, Italy) as high scalable applications to test on the Blue-Gene/Q architecture within the PRACE FP7 European project .
Tests on the second engine case are currently running on the HPC cluster made available by the Argonne National Lab (Chicago, Illinois) within the PETSc-Foam project.
Visualizations presented were generated using ParaView® by Kitware. OpenFOAM in its original implementation is licensed and distributed by the OpenFOAM Foundation and developed by OpenCFD.