3D CFD simulation of a gaseous fuel injection in a hydrogen-fueled internal combustion engine

Nowadays, one of the hottest topic in the automotive engineering community is the reduction of fossil fuels. Hydrogen is an alternative energy source that is already providing clean, renewable, and efficient power being used in fuel cells. Despite being developed since a few decades, fuel cells are affected by several hurdles, the most impacting one being their cost per unit power. While waiting for their cost reduction and mass-market penetration, hydrogen-fueled internal combustion engines (H2ICEs) can be a rapidly applicable solution to reduce pollution caused by the combustion of fossil fuels. Such engines benefit from the advanced technology of modern internal combustion engines (ICEs) and the advantages related to hydrogen combustion, although some modifications are needed for conventional liquid-fueled engines to run on hydrogen. The gaseous injection of hydrogen directly into the combustion chamber is a challenge both for the designers and for the injection system suppliers. To reduce uncertainties, time, and development cost, computational fluid dynamics (CFD) tools appear extremely useful, since they can accurately predict mixture formation and combustion before the expensive production/testing phase. The high-pressure gaseous injection which takes place in Direct-Injected H2ICEs promotes a super-sonic flow with very high gradients in the zone between the bulk of the injected hydrogen and the flow already inside the combustion chamber. To develop a methodology for an accurate simulation of these phenomena, the SoPHy Engine of the Engine Combustion Network group (ECN) is used and presented. This engine is fed through a single nozzle H2-injector; planar laser-induced fluorescence (PLIF) data are available for comparison with the CFD outcomes.


Introduction
Owing to its physical and thermophysical properties [1,2], hydrogen has shown promising benefits for applications as fuel in ICEs or fuel cells. However, fuel cells, which, so far offer the highest well-to-wheels benefit in terms of efficiency and emissions [3][4][5], are struggling to emerge on the market because of their expensiveness and the lack of hydrogen-refueling infrastructures. Many experts think that it will take no less than 10 years before their mass production [3,6,7]. Conversely, H2ICEs can take advantage of the current ICEs manufacturing infrastructure to be mass-produced on budget resources. Furthermore, several studies show that advanced H2ICE powered vehicles can outperform gasoline ones [2,8]. For these reasons, H2ICEs are seen as a bridging technology that will help to accelerate the development of H2-based infrastructures. Advanced H2ICEs are based on direct injection (DI) fuel systems. In contrast to port-fuel injection (PFI) systems, H2-DI systems not only avoid the loss of power due to the H2 low volumetric efficiency, but also can reduce combustion anomalies such as knock, backfire, and pre-ignition, which have harmful effects on ICE performance and emissions [9]. DI can also reduce heat transfer losses, which are relatively high in H2ICEs, due to the low flame quenching distance [10]. Several works show that the extent of benefit from DI largely depends on the injection strategy and on the subsequent mixture formation process. The injection strategy must be carefully set to find the trade-off between efficiency and NOx emissions. However, this requires deep knowledge of the injection event, in-cylinder mixing phenomenon, and the resulting mixture distribution before combustion takes place. Therefore, understanding the details of the injection process is an essential step to improve the efficiency of H2ICEs. As in commercial diesel and modern gasoline direct injection (GDI) ICEs [11,12], H2 injectors feature multi-hole nozzles to aid the fuel/air mixing process. However, to gain a fundamental understanding of the physical processes involved in mixture preparation and as a benchmark for the CFD simulation validation, single-hole injectors have proven to be very valuable. As a step towards this goal, i.e. the deep knowledge of the gaseous injection and subsequent mixing process, the outcomes of the experimental work at Sandia National Laboratories in an optically accessible DI-H2ICE, e.g. the SoPHy Engine, is compared with the outcomes of our CFD simulations. In the simulations, the commercial code SimCenter STAR-CCM+ v2021.1, licensed by Siemens DISW, is used for 3D modeling of the in-cylinder processes based on the Reynoldsaveraged Navier Stokes (RANS) equations. Previous works in literature [13,14] used two computational grids covering the portion of the engine cycle before the injection event and the remaining part of the cycle respectively. In this work, a single computational grid is used for the entire engine cycle, with an in-house developed adaptive refinement methodology that allows a correct resolution of the grid during the H2 injection. A sensitivity study is presented for both the grid size during the injection and the time-step. As aforementioned, the outcomes of 3D-CFD simulations are compared with the experimental available one in terms of H2 mole-fraction over the vertical plane through the injector nozzle symmetry axis.

Experimental Apparatus
SoPHy engine is a passenger-car-sized optically accessible single-cylinder engine, which uses a four-valve pent-roof head designed by GM research and is adapted for H2 operation and to allow optical access. The main customization to the head consists of two flat windows that provide optical access to each side of the pent-roof. The intake system consists of two intake ports, one for each valve, both straight and parallel to each other, tilted at 40° with respect to the fire deck. A Detailed list of engine and operating point parameters is reported in Table 1, while the engine geometry is shown in Figure 1. In addition to the pent-roof windows, optical access to the combustion chamber is guaranteed by fused silica liner and a flat Bowditch-type extended piston, with a window covering the central 65 mm of the 92 mm bore. H2 from a high-pressure vessel is supplied directly into the combustion chamber by a solenoid-driven injector, developed by Westport Inc. The H2 injector is placed in the sparkplug hole by an adapter, thus it has an offset of 3.2 mm towards the exhaust side due to the asymmetry of the pent-roof. As a consequence, this setup does not allow firing operation of the engine. As shown in Figure 2 a), the injector has a single hole-nozzle with a hole diameter of 1.46 mm. The hole is tilted at 50° with respect to the cylinder axis. The engine operates in motored conditions at a constant speed of 1500 rpm. Nitrogen is supplied as a bulk gas to make the intake charge inert and to avoid spectroscopic complications. H2 is injected at 100 bar for 17.5°CA in a single injection per cycle, which results in a reconstructed equivalence ratio (φ) equal to 0.25. As a matter of fact, the experimental equivalence ratio is calculated from steady-state bulk-gas and H2 flow rates, with nitrogen as a replacement for air. The average fuel flow is experimentally measured.  PLIF of gaseous-acetone, used as a fuel tracer, is performed to obtain quantitative images of the H2 mole-fraction in the combustion chamber. The methodology used to measure the in-cylinder fuel distribution is described in detail by Halter F. et al. in [15] and is therefore briefly discussed here. For the sake of science, differential diffusion of acetone, i.e. the tracer, and H2, the fuel, is estimated to be negligible due to the large spatial structures and the high Reynolds number of the post-injection flow. The acetone fluorescence is excited by a quadrupled Nd:YAG at 266 nm laser. Images are separately acquired for the horizontal plane, as well as for the vertical one, shown in Figure 2 b), through the symmetry axis of the injector. Small regions in the vertical plane, especially for early crank angles, may have errors larger than 5% because of window fouling and signal reflections, critical for this imaging configuration.

Numerical Setup
In this study, the in-cylinder add-on of SimCenter STAR-CCM+ is used to carry out the 3D-CFD simulation of the engine cycle. A RANS approach to turbulence is adopted [16], which is preferred to other more refined techniques such as DES [17,18] or LES [19][20][21][22] since the focus is made at present on the average H2 gaseous-injection event, without claiming to cover phenomena such as CoV [23], or combustion anomalies [24,25]. The injector nozzle-hole is modeled by a mass flow inlet boundary, with specified experimental pressure and temperature as well as mass flow rate.

Computational Grid
The adopted software features a built-in meshing tool based on a cell trimming process coupled with a dedicated prismatic mesher for the near-wall grid. Mesh motion is handled by a morphing/remapping technique. Whenever cell quality drops below user-editable quality metrics, a new mesh of the entire computational domain is generated, and results are conservatively interpolated onto the new grid. The in-cylinder core grid spacing is equal to 1.0 mm, with fixed and moving control volumes (CVs) to reduce cell size near the valves down to 0.5mm at valve opening/closing events. Furthermore, crank angle-dependent control volumes, shown in Figure 4, are activated to handle the gaseous injection with a proper grid resolution. The near-wall grid of the intake and exhaust ports is composed of a single layer with a near-wall thickness equal to 0.5mm, while inside the combustion chamber and over the valve CVs the number of layers is increased to two maintaining the same total thickness. The maximum cell count is around 1.9 million cells at bottom dead center (BDC), and 1.2 at top dead center (TDC). The mesh is shown in Figure 3. For the sake of gris sensitivity analysis, a second computational grid, shown in Figure 5, is generated in which the grid resolution controlled by the injector CVs during the injection phase, is doubled. It has a lower number of cells across the nozzle-hole diameter, i.e. only 10 cells instead of 20, with the fourth CV having the same grid size as the in-cylinder region, thus leading to a maximum cell number of 1.7 million cells.

Physical Models
The Turbulence model adopted for the simulation is the Realizable Two-Layer [26] kepsilon in combination with the All-y+ wall treatment model, using a 2 nd order discretization scheme. A monotone advection and reconstruction scheme (MARS) is used for temporal, as well as spatial discretization. Time is solved by the Pressure Implicit with Splitting of Operator (PISO) algorithm. 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, while the Specific Heat is computed using a Mass-Weighted Mixture method. Wall heat transfer is modeled using the GruMo-UniMORE thermal wall function [27][28][29][30].

Results
CFD outcomes are firstly compared to their experimental counterpart in terms of H2 molefraction over the Y section plane across the injector axis. A very restricted dataset of detailed PLIF measurements is available, together with a very low-resolution PIV [31]. Thus, PIV measurements are discarded for the present comparison, although a qualitative analysis highlights a satisfactory prediction of the main flow structures. Focus is then made on PLIF. For the sake of understanding, a delta value is computed to quantitatively compare the CFD outcomes and the experimental counterparts, computed as in the following eq, Eq. 1: Where is the k-th point of the PLIF measurements, and is the computed counterpart, linearly interpolated onto the same grid as the PLIF data. As aforementioned, two different sensitivity analyses are performed, the former focusing on grid density and the latter on the computational timestep. Furtermore, based on previous works in literature [13,31], different injection mass flow rate ramps are analyzed.
Experimental SOI is -140CA before firing top dead center (BfTDC). However, during the experimental campaign, the ECN scientists highlighted a delay in the effective SOI of nearly 3CAs, thus leading to an official SOI of -137 CAs BfTDC. While instantaneous H2 mass flow rate is not measured, a mass flow rate profile can be reconstructed based on the injector needle lift profile and on H2 injected mass per cycle. For the sake of brevity, only , the most interesting PLIF fields are reported hereafter, i.e. those during and immediately after the injection event. Starting from the grid sensitivity analysis, the in-cylinder H2 mole fraction distributions resulting from the two CFD meshes are Figure 6 together with their experimental counterpart. A limited yet visible grid dependency is observed, with the coarses mesh in the rightmost columns exhibiting higher H2 diffusion in the combustion chamber than the finer grid and the experiments. A better agreement in terms of both jet penetration and jet diffusion is therefore found using the refined grid. To confirm this, H2 mole fraction Delta fields are reported in Figure 7 for the two different simulations. Flow structures are similar for the two computed fields, but the H2 jet computed by using the coarse grid exhibits a higher momentum than that of the refined grid, leading to an overestimation of the penetration with respect to the experimental counterpart at 590CA degrees. Such overestimation persists at 600 CA degrees, Figure 7 d) e) and f), and, to some extent, is still visible until 620 CA degrees. Moving to timestep sensitivity, results presented so far were obtained using a constant timestep of 0.001 CA degrees during the whole injection process. Thus, being the injection duration 17.5 CAs degrees, the total timesteps required for a full cycle increase from less than 14000 for the unfuelled case to more than 30000 for the fuelled one. This further emphasizes the huge increase of computational cost of the fuelled case, as a result of the relevant increase of both cell count and timesteps. Similarly to the previous analysis, an attempt to reduce the computational cost is hereafter presented using larger timesteps. Results are reported in Figure 8. The qualitative comparison of the H2 mole fraction field, shown in Figure 8, suggests that there is not a clear dependency of the simulation outcomes on the adopted timestep. Such observation is confirmed by previous works available on literature [13,14,31], which relied on timesteps 10 to 50 times larger than those used in this work. However, using coarser timesteps, the convective courant number increases up to critical values, which in turn make the number of PISO corrector stages hit the maximum allowable during most of the injection process. The observation of the H2 delta fields, shown in Figure 9, confirms that the computed H2 field is negligibly affected by the timestep size. Looking at the total time to cover a complete engine cycle, the adoption of coarser timesteps dramatically reduces the CPU cost of the simulation: while the refined timestep case requires almost 4 days, the coarse timestep case takes just 1 day and a half, on equal number of CPUs, i.e. 64. Conversely, coarsening the grids on equal timestep leads to less than halved elapsed time, with the penalty of an increased diffusion of the H2 jet. Therefore, it can be assessed that for the specific case under investigation timestep coarsening is a better choice than gid coarsening.
Moving to the analysis of the instantaneous mass flow rate profiles, and particularly to the slope if the opening/closing ramps, previous authors experienced a better reproduction of the Schlieren data with a ramp of 4 CA degrees, instead of the original 1 CA degree one. It is worth to remind that no experimental evidence is available on the instantaneous injection law.  Figure 10 shows once again the computed H2 mole fraction fields comparing the two assumed ramps to the experimental counterpart. Focusing on 590 CA, a slightly better accuracy, in terms of jet penetration and diffusion, seems to be achieved by the longer ramp case. Conversely, for the remaining CAs, and especially for the one close to end of injection (EOI), a relevantly lower concentration of H2 is detected in the jet region, driven by the anticipated reduction of the injected mass. Delta fields for this comparison are not be reported in this work since the differences between the computed fields are evident from the fields reported in Figure 10.

Summary/Conclusions
The H2 gaseous injection of the SoPHy engine was successfully computed in a RANS framework, using commercial 3D-CFD tools, i.e. SimCenter STAR-CCM+ In-Cylinder addon v2021.1, licensed by Siemens DISW. A methodology to control the grid size CAs dependent was developed and used, increasing the computational efficiency and the solution accuracy.
Two different sensitivities were carried out, and the outcomes compared with their experimental counterpart, show that instead of using a coarse computational grid during the injection phase is better to use a coarse time-step, similar to the one used for motored simulations. Thus, reducing the time requested for the simulation of a full engine cycle. It must be considered that, in this work, the scalability shown by the software [32], is not used thus it will possible to use the more refined setup (both in terms of timestep as well as of computational grid) increasing the CPUs number to reduce at minimum the elapsed time to cover a full engine cycle. The ramp of 4CAs for the mass flow rate seems to lead to a loss of accuracy in the computed fields with respect to the experiments. A different ramp setting could lead to different outcomes, an asymmetric ramp mass flow rate setup, 4CAs going up and only 1CA going down from the static mass flow rate, is currently on our cluster.