Impact of Grid Density on the Analysis of the In-Cylinder Flow of an Optical Engine

The evaluation of Internal Combustion Engine (ICE) flows by 3D-CFD strongly depends on a combination of mutually interacting factors, among which grid resolution, closure model, numerics. A careful choice should be made in order to limit the extremely high computational cost and numerical problems arising from the combination of refined grids, high-order numeric schemes and complex geometries typical of ICEs. The paper focuses on the comparison between different grid strategies: in particular, attention is focused firstly on near-wall grid through the comparison between multi-layer and single-layer grids, and secondly on core grid density. The performance of each grid strategy is assessed in terms of accuracy and computational efficiency. A detailed comparison is presented against PIV flow measurements of the Spray Guided Darmstadt Engine available at the Darmstadt University of Technology. As many research groups are simultaneously working on the Darmstadt engine using different CFD codes and meshing approaches, it constitutes a perfect environment for both method validation and scientific cooperation. A motored engine condition is chosen and the flow evolution throughout the engine cycle is evaluated on two different section planes. Pros and cons of each grid strategy are highlighted and motivated.


INTRODUCTION
CFD simulation of in-cylinder flows is now a fundamental tool during the design and development of new ICE concepts. Thermo-physical processes taking place in the cylinder are extremely complex and they strictly influence the engine efficiency, both in terms of pollution and performance [1]. Using a CFD tool it is possible to study a multitude of ICE phenomena, among which fuel/air entrainment and mixing [2][3][4][5][6][7][8][9][10], combustion [11][12][13], heat transfer [14][15][16][17][18] and pollutant formation [19][20][21][22][23]. The available literature acknowledges the remarkable potential of Large Eddy Simulation (LES) as opposite to the more standardized Reynolds Averaged Navier Stokes (RANS) to simulate the turbulent nature of engineering flows [24]. Such turbulent nature is primarily responsible for cycle-to-cycle variability of incylinder phenomena, which in turn severely affects efficiency. However, LES simulations require a high number of cycles to be run in order to perform statistical analysis of the outcomes, thus the computational cost is much higher than the one of RANS simulations, which require only few cycles to ensure cyclic repeatability. This is the main reason why RANS simulation is the most used approach for industrial applications. RANS simulations are based on a time/ensemble averaging operation which splits any flow variable in a timeindependent contribution and in a turbulent fluctuation; the latter involves the whole range of turbulent scales responsible for the deviation of the instantaneous signal from the averaged one. Nonetheless, RANS simulations have also been used to study phenomena that are extremely variable in time, such as knock [25][26][27][28].
The k-ε turbulence model is the most used for RANS simulations. It belongs to the twoequation eddy viscosity models class, in which transport equations are solved for two turbulent quantities, specifically k, the turbulent kinetic energy, and ε, the turbulent kinetic energy dissipation. The k-ε model is arguably the simplest complete turbulence model and it has an extremely wide range of applicability [29,30], moreover many different versions have been proposed since its introduction [31]. Gosman and Watkins in 1977 extended the original k-ε formulation [17] to the engine flows by including the effect of compressibility in the constitutive equation and by taking into account the density variations over time. Later, a socalled Realizable k-ε model was proposed by Shih et al. [32], containing a new transport equation for the turbulent dissipation rate based on the dynamic fluctuation of the vorticity. In the calculation of the turbulent eddy viscosity, the coefficient C_μ is expressed as function of mean flow and turbulence properties. These improvements make the Realizable k-ε more accurate than the standard k-ε in many applications [33]. For this reason, the Realizable k-ε model is chosen in the following paragraphs.
An equally impacting aspect in ICE CFD is the grid strategy and the mesh motion technology adopted. The authors carried out a detailed study on the impact of different mesh motion technologies in [34]. In the present paper, a single strategy, namely the morph and remapping technology is adopted, while focus is made on the impact of grid density. In RANS simulations, near-wall phenomena are modelled by a wall treatment. In commercial CFD software there are two main strategies. In the first one, the CFD tool spatially resolves the near-wall profiles of velocity and temperature by working with a refined computational near-wall grid and exploiting a low-Reynolds number or a two-layer wall treatment [27]. In the second one, near wall flow temperature in the proximity of the walls are reconstructed by using algebraic wall functions; the resulting high-Reynolds number wall treatment avoids the need of fine near-wall grids typical of the low-Reynolds number counterpart, reducing the computational effort [14,18,35]. In this work, a two-layer all y+ model is adopted; a multi-layer near-wall cell grid is compared to a single layer one to assess the performance in terms of quality of results vs. computational effort.
To this aim, the Darmstadt Spray-Guided Engine operated under motored condition is simulated and each grid strategy is evaluated in terms of computational efficiency and fidelity of the outcomes to the measured particle image velocimetry (PIV) flow fields.
The paper is organized as follows: initially, a description is provided of the experimental apparatus; this description is then followed by the numerical setup and an overview of each grid strategy. Subsequently, results are shown to compare calculated and experimental flow fields. Finally, conclusions are drawn.

EXPERIMENTAL APPARATUS
The so-called Darmstadt Engine is an optical access four-stroke single-cylinder engine at Technische Universität Darmstadt, illustrated in Fig. 1. It is a direct-injection spark-ignition engine featuring a twin cam, overhead-valve pent-roof cylinder head. Two different cylinder heads are available for the experimental tests, one is intended for a wall guided direct injection engine operation, the other for a spray guided one. In the spray-guided version, object of our study, the injector is centrally mounted close to the spark-plug, which is on the exhaust valve side.  Table 1. The design of the test bench allows conditioning pressure, temperature, and gas composition of the intake port, as well as port fuel injection of liquid and gaseous fuels. Two cylindrically shaped plenums are present both at the intake and at the exhaust to limit the effect of pressure perturbations due to singlecylinder operation. The intake port hosts two different pressure sensors, labeled P in,1 and P in,2 . In particular, P in,2 provides intake port pressure measurements close to the intake valves. Pressure signals inside the cylinder as well as in the intake ports were acquired every 0.5 Crank Angle degrees (CA), at the above-mentioned positions, also shown in Fig.2. The optical access to the cylinder is available thanks to an optical-grade fused-silica cylinder liner and flat piston window, the latter connected to a Bowditch piston extension for mirror access. This setup guarantees the visibility of the first 55mm of the stroke plus an additional 8mm in the pentroof area. PIV measurements were conducted on two different section planes, the symmetryplane, and the valve-plane (hereafter referred to as Y= 0 mm, Y= 19 mm). PIV images were acquired every 5 CA using a monochrome CMOS high-speed camera (Vision Research, Phantom v1610) and a 105mm F2.8 macro Sigma lens. More information about the experimental setup can be found in [36,37].

Numerical Setup
In this study, the RANS investigation of the Darmstadt spray-guided engine is performed by using a commercial CFD-3D software tool, SimCenter STAR-CCM+ v2020.1 licensed by Siemens PLM.
Two different core-grid densities, combined with two different near-wall grid strategies are adopted, using the same mesh motion technology based on morphing/remapping methodology. Besides different grid strategies, the simulations share the same physics models and settings. Time is solved by the Pressure Implicit with Splitting of Operator (PISO) algorithm. A monotone advection and reconstruction scheme (MARS) is used for temporal and spatial discretization. Turbulence is modelled using the Realizable K-Epsilon in combination with the All-y + wall treatment model. The Ideal Gas equation of state is employed to compute the density and its derivatives as a function of temperature and pressure. Sutherland's law is used for the Dynamic Viscosity. A calibrated 1D model of the engine in the GT-Power tool is used to provide time-dependent boundary conditions (BCs) for pressure and temperature. A detailed comparison between the adopted boundary conditions and the experiments is carried out in a companion paper.

Calculation Grid
The computational domain does not include all the components of the experimental apparatus. As a matter of fact, to limit the total number of cells, both plenums and the first part of the intake ports, until the P in,2 .sensor position, are not included in the CFD domain. Concerning the exhaust side, the computational domain is truncated at the location of the pressure sensor shown in Fig. 2. The meshing tool uses a cell trimming process in combination with a dedicated prismatic mesher for the near-wall grid. A mesh extrusion of 30mm with 10 layers was used both at inlet and outlet boundaries to ensure stability, without modifying the geometry. To compensate, the original ports geometry was shortened by the mentioned 30mm.

Mesh #1
Of the three different grids that were tested, Mesh #1 is the finest one. The in-cylinder core grid spacing is equal to 0.75mm, with fixed and moving control volumes (CVs) near the valves to reduce the grid size to 0.375mm. By using two different CVs the mesh in the intake port, as well as in the exhaust port, is refined to 0.75mm near the valves and then coarsened to 1.5mm moving further away from the valve stem region. A multi-layer approach is adopted for the near-wall grid over the entire computational domain. Near-wall grid is composed by four layers with a near-wall thickness equal to 0.05mm and a total thickness of 0.5mm.
A CV is added to accurately represent the complex geometry of the spark plug zone. Here the grid size is reduced to 0.1875mm, while the near-wall grid is reduced to two layers with a near-wall cell size equal to 0.04mm and a total thickness of 0.08mm. With this setup, the maximum cell count in the domain is around 5.1 million (of which, 3 million devoted to the in the cylinder region) at top dead center (TDC), and 3.2 million (1.3 million in the cylinder region) at bottom dead center (BDC).

Mesh #2
Mesh #2 uses the same setting as Mesh #1. However, a single layer replaces the previously adopted four in the near-wall region, with a thickness equal to 0.3mm, the only exceptions being the valve seat regions and the spark plug region. The former uses two layers, each 0.15mm thick, to ensure an adequate number of grid points in the valve curtain during the opening/closing valve periods. The latter retains the same setting as Mesh #1. In Mesh #2, the maximum cell count in the domain is around 3.4 million (2.2 million in the cylinder region) at TDC, decreasing to 1.8 million (0.9 million in the cylinder region) at BDC.

Mesh #3
Mesh #3 has the same near-wall prism-layer settings as Mesh #2; however, the core grid size is increased by scaling the spacing by a factor of nearly 60%. The same settings as those for Mesh #1 and #2 are adopted at the spark plug. The total number of cells is approximately 1.75 million (1 million in the cylinder region) at TDC, while the cells at the BDC are around 1.2 million (0.7 million in the cylinder region).

RESULTS
Three consecutive RANS cycles are carried out for each of the above-mentioned grids to discard dependency of the results on the initial conditions; results of the latest cycles are shown in the present paper for the sake of brevity. A quantitative comparison with the experimental measurements is carried out in order to evaluate the impact of the different grids on the accuracy of the solution. The comparison between high speed PIV measurements and RANS computations focuses on three different CA positions onto the two available PIV planes. The chosen CAs for the analysis represent significant engine events: • middle of intake stroke at 470 CA (maximum intake valve lift) • BDC (540 CA) • middle of compression stroke (630 CA) Temperature measurements are not available; therefore, the attention is focused on the incylinder pressure traces and the averaged intake port mass flow rate. As shown in Figure 3, Mesh #2 and #3 slightly overestimate the in-cylinder pressure at the end of the compression stroke, while Mesh #1 is well aligned to the experimental trace. Nevertheless, differences seem to be negligible for the sake of the present activity. As for the intake average mass flow rate, Mesh #1 matches the experimental value of 11.50 kg/h +/-2% with a value of 11.5012 kg/h, while Mesh #2 and #3 underestimate it with a value of 11.4276 kg/h and 11.4373 kg/h respectively. However, all the values are within the tolerance range.

CFD vs PIV
To make an in-depht comparison, CFD results (i.e. flow variables evaluated on the cell centroids) are processed in MATLAB and linearly interpolated on the same experimental grid as the PIV one. PIV data are stored on a 132x129 grid with a base size of 0.6023 mm on Y=0 plane and averaged over more than 220 cycles, while on Y=19 plane PIV data are stored on a 125x104 grid with a base size of 0.5036 mm and averaged over 250 cycles.
Even though the flow fields over the two different PIV planes are compared at aforementioned CAs, only those over Y=19, identified as the most representative plane, are shown for the sake of brevity. In Figure 4a) 4b) 4c) and 4d)the velocity components U and W are shown over the Y=19 plane at 470 CA. Both Mesh #2 and #3 match the penetration of the flow incoming from the intake port, while Mesh #1 slightly overestimates it. Mesh#1, however, seems to better predict the bending of the incoming jet in comparison with the other two grids. Furthermore, only the finest grid accurately resolves the secondary flow near the piston top in the lower right corner of the field of view (FOV). In Figure 4e) 4f) 4g) and 4h) the velocity components U and W are shown at 540 CA, while in figure 4i) 4l) 4m) and 4n) they are shown at 630 CA, over Y=19 plane. The velocity fields obtained using Mesh #1 grid are the most accurate during the full cycle length, since they closely reproduce the experimental flow structures.

Alignment Parameter
To perform a more quantitative analysis of the accuracy of calculated fields, the alignment parameter is employed to rate the local deviation of the calculated velocity field from the experimental one, assumed as a reference. The alignment parameter is defined as follow: being A k an experimental vector field while B k is the calculated counterpart, therefore is the angle between the two considered vectors. If the alignment parameter is equal to 1, there is a perfect match between the directions of the PIV and of the calculated field; conversely, the directions are opposite if the alignment parameter is equal to -1. Given its definition, the alignment parameter provides only information on the relative direction of the considered fields, while no information on their magnitude can be inferred. The local AP field are plane-averaged on Y = 0 and Y = 19 planes to obtain two global indicators, ���� (one for each section plane), as suggested by [38,39]. The two indicators ���� are obtained using the following equation (2), with N being the total number of sampling points in the FOV: (2) Plane-averaging is conditioned on the experimentally limited FOV, discarding out-ofview areas (marked with Not a Number (NaN) value in all the figures shown here); furthermore, it attributes an equal relevance to all areas. In order to avoid misleading observations an original variation to such indices was proposed in [40]. A weighted averaging (wa) operation is introduced, based on the sum of the 2D PIV and numerically predicted kinetic energies (KE) for the k−th point of the FOV (Eq. (3)). = The �������� is defined as in Eq. (4), and it is used to draw more engineering relevant observations, by attributing major relevance to the correct simulation of high kinetic energy regions. In Figure 5 the alignment parameter field is shown over the Y=0, figure 5a) 5b) and 5c), and over Y=19, figure 5d) 5e) and 5f) respectively for Mesh #1, #2 and #3, at the maximum intake valves lift (470 CA). The general agreement between measured and predicted fields is confirmed for all the three grids. However, some differences can be observed in low speed flow areas such as vortex cores and recirculation zones and close to the walls; the AP field emphasizes such differences and it is therefore a useful tool to perform detailed comparative analyses.
The visual comparison can be supported by a more quantitative evaluation using the aforementioned ���� parameter. This achieves values of 0.8524 on Y=0 and to 0.8874 on Y=19 for the refined grid (Mesh #1), Mesh #2 reaches 0.8544 and 0.8503 and Mesh #3 0.8893 and 0.8522 respectively. It appears therefore difficult to clearly identify a ranking for the three grids. A better understanding comes from the use of the KE-wa method: in this case Mesh #1, due to the multi-layer approach and the higher core grid density, reaches a higher value than the other two both on Y=0 and Y=19 planes. More in depth, Mesh #1 reaches a value of AP wa equal to 0.9812 and 0.9597 over Y=0 and Y=19 mm respectively, instead of 0.9766 and 0.9382 for Mesh #2 and 0.9805 and 0.9423 for Mesh #3 grid. In Figure 6 the alignment parameter field, at the mid-compression stroke, is shown over the Y=0, figure 6a) 6b) and 6c), and over Y=19, figure 6d) 6e) and 6f) respectively for Mesh #1, #2 and #3. Considering the KE-wa method Mesh #1 reaches a value of AP wa equal to 0.9828 and 0.9798 over Y=0 and Y=19 mm respectively, instead of the 0.9760 and 0.9730 of Mesh #2 and 0.9680 and 0.9814 achieved by using Mesh #3 grid. So, Mesh #1 achieves a higher KE-wa value over the Y=0 plane and a value comparable to Mesh #3 over the Y=19 plane.

Magnitude Index
The relative direction-only information of AP is enriched through the Magnitude Index (MI) indicator [38], which considers both relative orientation and magnitude of two vector fields A and B, as defined in Eq. (5).
A perfect agreement between the two evaluated fields is expressed by MI = 1, while an MI = 0 indicates the mismatch of the two. This is a more comprehensive indicator than AP, thanks to the combined magnitude and direction evaluation. Areas of low MI are the vortex cores and recirculation areas in the proximity of cylinder walls. As for the AP indicator also for the MI a plane averaging and a weighted averaging is performed, respectively defined as in Eq (6) and (7).
For the sake of brevity the results of the MI field over the two considered section planes are reported only for the maximum intake valves lift (470 CA), results at the other two considered CA follow the same trend and are omitted for the sake of brevity. In Figure 7 the magnitude index field at 470 CA is shown over the Y=0, figure 7a) 7b) and 7c), and over Y=19, figure 7d) 7e) and 7f) respectively for Mesh #1, #2 and #3. Once again, a visual comparison suggests that the highest values are achieved by the refined grid (Mesh #1). Considering the KE-wa method Mesh #1 reaches a value of MI wa equal to 0.8712 and 0.8222 over Y=0 and Y=19 mm respectively, instead of the 0.8314 and 0.7726 of Mesh #2 and 0.8394 and 0.7803 achieved by using the Mesh #3 grid. Averaged values confirm therefore that Mesh #1 outcomes are those closer to the experimental PIV data.

SUMMARY/CONCLUSIONS
Three different computational grids are used for the analysis of the Darmstadt engine under motored conditions. Results are analyzed through the comparison with available PIV experiments; all the tested grids allow the software to adequately represent the macroscopic characteristics of the flow structures inside the cylinder. Therefore, they can all be considered reliable and accurate. However, differences in terms of accuracy and computational efficiency can be analyzed more in depth using quantitative indices such as the Alignment Parameter and the Magnitude Index.
Mesh #1 is the finest computational grid. A multi-layer approach is used for the prism layer, in conjunction with the highest core grid density; results show the best overall level of accuracy over the three different grids. Despite it requires higher computational cost, the high scalability performance of the adopted tool leads to complete a cycle in approximately 40 hours on 64 CPUs.
Mesh #2 shares with Mesh #1 the same core grid density, but here a single near wall layer grid is used instead of a multi-layer approach, with a near-wall cells size 6 times thicker than the one of Mesh #1. Almost halving the total amount of cells, the time spent for a full cycle on equal CPU number is reduced to around 31 hours. However, the multi-layer approach provides better results over the entire cycle.
Mesh #3 shares the near-wall mesh with Mesh #2, but a coarser core grid size is adopted in order to limit the total amount of cells. The full cycle, by using this computational grid settings, is lowered to 21 hours approximately. The results are anyway reasonably good compared to the experimental data.
Aside from the computational cost each grid requires, this study shows how the core grid size and near-wall cell size notably impact on the results, thus confirming that having a highquality/high-density mesh is of paramount importance in order to obtain high quality results.
Mesh #1, the one with the highest grid density and a multi-layer near-wall approach, is the one that provides the best results, even if it has the highest computational costs. Since the adopted software shows very high scalability performance and stability, a higher number of CPUs can be used to reduce the time needed to cover a full engine cycle.