Unsteady Engine Analysis with Moving Mesh in OpenFOAM®

By Prof. Federico Piscaglia
Dipartimento di Energia, Politecnico di Milano (Italy)

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 [1]. 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 [2]. 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:

Simulation of an engine-like geometry consisting of a flat-top cylinder head with a fixed, axis-centered valve and a low-speed piston [3]. Geometry specification as well as a list of main functioning parameters are represented in Fig. 1.
Simulation of an optical engine with a two-valve head with simple intake and exhaust port/runner geometries and a pancake-shape transparent combustion chamber (TCC) [4].

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 [10] 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 = ρ U_p D ⁄ μ = 2000

Figure 1: Top) Geometry of the experimental apparatus used by Morse et al. [10]. 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 U_p as the piston velocity

U_p(t) = S ⁄ 2 sin(ωt) i_c

where S is the piston stroke and i_c the unit vector parallel to cylinder axis, the point velocity up was calculated as:

u_p(x_p) = U_p ⋅ ( z_max − z_p ) ⁄ ( z_max − z_piston )

In the equation above, x_p = (x_p,y_p,z_p) is the point position, z_max is the z-coordinate of the farthest moving point from the cylinder and z_piston is the z-coordinate of the piston.

A compressible solver for compressible flows with dynamic meshes, developed by the authors [11] 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 [12]. 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).

Figure 2: Turbulent velocity field from LES calculation by OpenFOAM.

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 [3], 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. [10] (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 [3].

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 [10].

Figure 4: Comparison of simulated (left) and experimental (right) streamlines within the cylinder region at: a) CA =36° ATDC; b) 90° ATDC; c) 144° ATDC; d) 270° ATDC.

The spatial distribution of the mean turbulent kinetic energy during the intake stroke is plotted in Fig. 5. With increasing piston speed, k_max 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 k_max 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.

Figure 5: Mean turbulent kinetic energy at different crank angles during the intake stroke.

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 slidingInterface, the layerAdditionRemoval and the attachDetach of boundaries. In addition, a new class to handle the interaction with the Finite Volume (FV) mesh has been designed. The slidingInterface mesh modifier creates a seamless joint between two regions of the FV domain, even if the mesh is non-conformal at the interface. The standard implementation of slidingInterface has been deeply modified to allow for restoration of the original (split) mesh configuration after the initial coupling, without any loss of quality. The resulting implementation can be used to simulate partially overlapping interfaces with relative motion. Also, the merge algorithm has been strengthened to correctly handle unstructured meshes, non-manifold patches and shared points between the sides of the non-conformal interface. Thanks to these features, it is possible to stitch mesh regions in complex cases, such as unstructured grids, sharp corners and progressive refinement boxes. The layerAdditionRemoval class has been improved in its parallel operation, in order to deal with complex domain decompositions and avoid clashes with unstructured regions or mesh boundaries. An adaptive time-stepping procedure has been implemented in the dynamicFvMeshClass to ensure the mesh validity at any time and avoid geometrical inconsistencies (e.g. negative volume cells). The attachDetach class has been revised, too, by implementing a geometric point matching algorithm that makes it more robust when included in complex meshes. The construction, activation and interaction of the mesh modifiers with the FV mesh are managed at a low level; in this way, the use of topology modifiers transparent to the developer can be added, as new point motion strategies based on topological changes are developed. As a consequence, the user can set up a large number of different dynamic mesh cases very quickly with little or no programming effort. The present implementation is completely based on the mesh definition of the OpenFOAM version released by ESI-OpenCFD^® and works fully in parallel. A full description of the advanced algorithm developed for mesh motion with non-conformal grids is reported in [9, 11].

Example of the strategy to generate the engine mesh for the simulation of a transparent combustion engine [4] with canted valves using Pointwise.

Figure 7: Example of the strategy to generate the engine mesh for the simulation of a transparent combustion engine [4] 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 [11] 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.

Figure 9: LES of an IC engine with canted valves [11].

ACKNOWLEDGMENTS

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 [14].

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.

References

S. Pope, Turbulent flows. Cambridge University Press, 2000.
A. Amsden, P. O'Rourke, and T. Butler, KIVA-II: A Computer Program for Chemically Reactive Flows with Sprays. LA 11560-MS, Los Alamos National Laboratory, 1989.
A. Montorfano, F. Piscaglia, and A. Onorati, “A LES study on the evolution of turbulent structures in moving engine geometries by an open-source CFD code,” SAE Technical Paper 2014-01-1147, 2014.
Engine Combustion Network, http://www.sandia.gov/ecn/engines/engineFlows/TCCEngine.php, 2013.
F. Piscaglia, A. Montorfano, and A. Onorati, “Development of a non-reflecting boundary condition for multidimensional nonlinear duct acoustic computation,” Journal of Sound and Vibration, vol. 332, no. 4, pp. 922–935, 2013.
F. Piscaglia, A. Montorfano, and A. Onorati, “Improving the simulation of the acoustic performance of complex silencers for ICE by a Multi-Dimensional non-linear approach,” SAE Int. J. Engines, vol. 2, no. 5, pp. 633–648, 2012. doi:10.4271/2012-01-0828.
F. Piscaglia, A. Montorfano, A. Onorati, and F. Brusiani, “Boundary conditions and subgrid scale models for LES simulation of IC engines.” Accepted for publication on Oil & Gas Science and Technology – Revue d'IFP Energies nouvelles, 2013.
F. Piscaglia, A. Montorfano, and A. Onorati, “Towards the LES simulation of IC Engines with parallel topologically changing meshes,” SAE Int. J. Engines, vol. 6, no. 2, pp. 926–940, 2013.
F. Piscaglia, A. Montorfano, and A. Onorati, “Development of fully-automatic parallel algorithms for mesh handling in the openfoam-2.2.x technology,” SAE Technical Paper 2013-24-0027, 2013.
A. P. Morse, J. H. Whitelaw, and M. Yanneskis, “Turbulent flow measurements by laser-doppler anemometry in motored piston-cylinder assemblies.,” Journal of Fluids Engineering, vol. 101, pp. 208–216, 1979.
F. Piscaglia, A. Montorfano, and A. Onorati, “A moving mesh strategy to perform adaptive Large Eddy Simulation of IC engines in OpenFOAM.” International Multidimensional Engine Modeling User's Group Meeting 2014, The Detroit Downtown Courtyard by Marriott Hotel, Detroit, MI (USA), April 7, 2014. https://imem.cray.com/2014/Meeting-2014/9-Piscaglia-Milano-IMEM2014.pdf.
H.Weller, “Controlling the computational modes of the arbitrarily structured C grid,” Monthly Weather Review, vol. 140, pp. 3220–3234, 2012.
F. Piscaglia, A. Montorfano, and A. Onorati, “Large eddy simulation of multiple cycles in a valve/piston assembly in the OpenFOAM technology,” In proceedings of International Journal of Engine Research - special issue on Cyclic Dispersion in Combustion Engines, 2015.
Interview to F. Piscaglia, “Faster and open engine simulations for automotive industry,” Prace Digest 2013, Industrial Edition, 2013.