Flow Field Investigation of a Single Engine Valve Using PIV, POD, and LES

Abstract

Due to stringent emission regulations, it is of practical significance to understand cycle-to-cycle variations in the combustion of fossil or renewable fuels to reach future emission regulations. The present study aims to conduct a parametric investigation to analyse the influence of the valve lift and different mass flows of an inlet valve of the test engine “Flex-OeCoS” on the flow structures. To gain a deeper understanding of the flow behaviour, an optical test bench for 2D Particle Image Velocimetry (PIV) and a Large Eddy Simulation (LES) are used. Turbulence phenomena are investigated using Proper Orthogonal Decomposition (POD) with a quadruple decomposition and the Reynolds stress transport equation. The results show good agreement between the PIV and LES. Moreover, the main flow structures are primarily affected by valve lift while being unaffected by mass flow variation. The turbulent kinetic energy within the flow field increases quadratically to the mass flow and to the decreasing valve lift, where large high-energetic flow structures are observed in the vicinity of the jet and small low-energetic structures are homogeneously distributed within the flow field. Furthermore, the convective flux, the turbulent diffusive flux, the rate of change, and the production of specific Reynolds stress are the dominant terms within the specific Reynolds stress transport equation.

1. Introduction

1.1. Motivation

Although rising awareness of climate change and more stringent emission regulations lead to a rising share of electric vehicles, the internal combustion engine will continue to play a key role particularly in the cargo transport area [1]. As environmentally harmful gases must be reduced to avoid the impact on our planet, future emission regulations will lead to increasing challenges in engine development. These gases mainly originate in the combustion of fossil or renewable fuels, which will remain indispensable in the future of the transport and energy generation sectors due to growing global demand and the advantage of high energy density [2]. Therefore, greenhouse gases cannot be entirely avoided in the future, and countermeasures from many sides of the problem must be applied.

The challenge to reduce fuel consumption and harmful emissions ranks among the primary goals of combustion engine development. The highest potential for improvement lies in thermal efficiency, where many different combustion processes have been thoroughly investigated (lean combustion, external exhaust gas recirculation (EGR), stratified combustion, or homogeneous charge compression ignition) [3]. One important limiting factor regarding combustion optimisation is cycle-to-cycle fluctuations, which can lead to higher engine emissions, lower efficiencies, and abnormal combustion phenomena [4,5,6,7]. Therefore, a great potential for the development of future high efficiency combustion engines lies in understanding and reducing cycle-to-cycle variations. One specific influencing factor is the inflow jet oscillations, which strongly affect the flow structure in the cylinder. This is particularly important for the velocity field around the spark plug at the time of ignition, but of course also affects other processes such as mixture formation and combustion.

The root cause of cycle-to-cycle fluctuations is a complex problem which has been investigated for decades [8,9,10,11,12], and, since the engine combustion process is the result of multiple parallel and successive mechanisms, the underlying causes influence each other. This leads to a proposed cause–effect chain of cycle-to-cycle variations [13], which includes the inlet flow, the injection strategy and mixture formation, and the thermofluidic state at the spark plug at the time of ignition. The first, and probably most important factor, is the inlet inflow, which primarily dictates the degree of turbulence and, therefore, most of the subsequent processes of mixture formation and combustion [4,14]. Therefore, the air flow through the engine valves plays a paramount role in the mixing and combustion processes and has been widely studied [15,16,17,18,19,20,21,22,23,24,25,26,27].

In combustion engines, different components such as the intake ports [16,17,18,21,28], fuel injection [11,21,29], piston shape [30,31,32], tumble/swirl flaps [26], intake valves [33,34] and valve lift [35,36,37] are designed to influence air motion. Two main flow types can be found in the cylinder, namely swirl, which spins along the cylinder axis, and tumble, which spins perpendicular to the cylinder axis. In diesel engines, the intake pipes are usually designed to improve swirl, while in gasoline engines the tumble motion is prioritised [17,28,38,39].

Investigations regarding cycle-to-cycle variations show high fluctuations in the tumble motion when the inlet valve closes, further in its breakdown, which eventually leads to differences in turbulence for the combustion [9,10,11,24,35,40,41]. Tumble flow is especially present in modern highly efficient gasoline engines. Earlier results have shown that engines with strong tumble flows cause stronger turbulence than those with swirl flows but are more prone to cycle-to-cycle variations in the flow field in the combustion chamber up to the ignition point [12]. Moreover, a recent numerical investigation showed that the main varying secondary flows, caused by inflow jet oscillations, seem to superpose the main charge flow and create variations in spark plug velocity which influence flame kernel growth [42]. According to Hasse [13], flow field fluctuations in the vicinity of the spark plug at ignition timing seem to be a leading factor in cycle-to-cycle variations.

One way to study the in-cylinder flow is to use traditional non-intrusive velocimetry techniques such as Laser Doppler Anemometry (LDA) [43,44], Magnetic Resonance Velocimetry (MRV) [45,46], and Particle Image Velocimetry (PIV) [26,47,48,49,50], as they have supplied valuable knowledge about evolving in-cylinder flows. The (high-speed) PIV method is a well-suited measurement technique to study (time-resolved) in-cylinder flow fields as well as to obtain information on large-scale fluctuations and cycle-to-cycle-variations [51,52]. Moreover, intake inflow jet oscillations, especially detachment at the valve seat as well as secondary flow structures, should be resolvable. One common way to further process the results obtained from PIV and extend the information obtained from this technique is to apply the Proper Orthogonal Decomposition (POD) filtering approach [17,24,47,53,54,55,56]. With this method, the measured flow fields can be separated into dominant, coherent, and turbulent small-scale structures. The results obtained from POD are valuable reference data for the validation of numerical models of cycle-to-cycle variations. In addition, the POD technique has been used to identify the spatial structure and temporal evolution of an engine flow obtained through 3D high-speed PIV [54].

Experimental studies tend to focus on one or two planes, such as the symmetry or tumble plane [57,58], the parallel planes [59], the cross-tumble plane (orthogonal to the symmetry plane) [60], and horizontal planes [53]. Although there are many optically-accessible engine test facilities that allow the application of various optical techniques; sometimes, different boundary conditions and better accessibility are required, and, therefore, ad-hoc test rigs are designed and built which replicate the air flow [17,32,33,34,39,61,62,63,64,65,66]. However, a disadvantage of such test rigs is often reduced flexibility in terms of operating parameters (such as intake temperature and pressure) and a restriction to steady-state flows with fixed valve lift positions.

The use of 3D Computational Fluid Dynamics (CFD) simulation offers the possibility of investigating complex phenomena such as cycle-to-cycle variations in detail. Effects concerning the internal cylinder flow and combustion can be studied at a significantly lower monetary cost compared to experimental measurements. For engine development, commonly, the Reynolds Averaged Navier Stokes (RANS) simulation method is used [52]. However, this method is not suitable for the investigation of cyclical fluctuations due to its averaging nature; therefore, the Large Eddy Simulation (LES) or Detached Eddy Simulation (DES) methods are more commonly used. In LES, large vortex structures are calculated directly, while small vortices, which are computationally expensive to resolve, are modelled. LES has been proven to be in good agreement with measurements in terms of cycle-to-cycle variations [51,67]. Compared to RANS, LES has higher demands on the temporal and spatial resolution of the flow field. Furthermore, several cycles must be calculated to reproduce cycle-to-cycle variations, leading to even higher computing times. As an alternative to LES, DES provides a good compromise between RANS and LES. To reduce computational costs, the RANS approach is used for near-wall areas or in regions where the grid resolution is not sufficiently fine to resolve smaller turbulent structures. Therefore, the spatial resolution in those areas can be coarser than for LES. For free shear layers or separated flows far from walls, the LES approach is used. DES of cycle-to-cycle variations has been carried out in a study by Hasse et al. [68] for a simplified engine setup and showed good agreement with measurements. The most used turbulence models in CFD include the κ-ε turbulent model [9,19,20,22,23,25,28,31,32,36,62,64,65,66,69], as this model uses two transport variables and performs well far from the boundary walls. LES also enables the prediction of the turbulent flow [10,17,19,32,41,45,50,64,65,70,71,72,73,74]; this method appeared as a solution to the issues faced by Direct Numerical Simulations (DNS), which is challenging and computationally expensive.

1.2. State of the Art

Recent studies have focused on investigating the effects of different engine parameters on the in-cylinder flow through PIV.

With respect to the air movement of swirl and tumble, a study by Rabault et al. [33] investigated the admission flow of an engine cylinder through measurements in the swirl and tumble planes using 3D PIV in a steady flow test bench. It was found that there was a recirculation bubble in the tumble plane despite the absence of a piston, which was due to the entrainment of gas into the jet. In addition, complex jet-dominated vortices were found in the swirl plane near the cylinder top, as well as counter-rotating vortex pairs further down the cylinder. An innovative study by Bücker et al. [48] investigated the tumbling flow in an optical engine using multi-planar PIV to quantify the tumble vortex structure in 3D by extending the two-dimensional criterion previously presented by Graftieaux [53]. A recent study by Aljarf et al. [75] characterised the in-cylinder gas flow of a compression–ignition engine under motoring and pre-ignition conditions using high-speed PIV. The results propose that increases in the swirl ratio, the mean-velocity, and the turbulent kinetic energy are related to an increase in the engine speed. These are also influenced by the movement of the induced swirl centres across the piston cavity during expansion to exhaust strokes.

Regarding the in-cylinder vortices, a study by Graftieaux et al. [53] combines PIV and POD with vortex identification algorithms to study unsteady turbulent swirling flows, resulting in additional computed quantities such as angular momentum and the location of the swirl centre, while a study by Coletta et al. proposed a simple engineering method to identify and characterise vortices within a flow field using classic two-component PIV [76].

In-cylinder flow has been studied with respect to mixture and cycle-to-cycle variations. A study by Bode et al. [77] used high-speed 3D PIV to study the influence of in-cylinder flow on cycle-to-cycle variations in a direct-injection spark-ignition engine. The main finding is that the favourable flow moving the flame kernel towards the centre of the cylinder reduced combustion duration. A study by Kim et al. [78] combined PIV and CFD to analyse the characteristics of in-cylinder flow and mixture formation by calculating the velocity, turbulent kinetic energy, and tumble ratio and centre.

The use of test benches has been used to study the flow under steady conditions and fixed valve position, which is a common practice in industry simulations as it allows different boundary conditions from those in reciprocating engine test facilities and provides valuable data for the validation of numerical models.

A study by Kapitza et al. [17] investigated the flow fields generated by a four-valve DISI cylinder head in a steady-state test-bench by using high-speed stereoscopic PIV. They developed a procedure that relies on Singular Value Decomposition (SVD), and it was shown that 20 SVD spatial modes are enough to reconstruct the flow fields from a real in-cylinder flow. Furthermore, the analysis of the transient behaviour showed that increasing the mass flow or decreasing the valve lift led to higher fluctuations in the flow.

A research work by El-Adawy et al. [34] studied the in-cylinder flow of an IC engine by using two types of test-benches under steady-state conditions by applying PIV among other methods. It was concluded that, at a higher valve lift, a stronger tumble motion is generated with an increased magnitude of turbulent kinetic energy and vorticity. A later study by El-Adawy et al. [63] investigated the in-cylinder flow at different tumble planes under steady-state conditions using stereoscopic PIV. They showed that at a high valve lift the vortex centre did not change significantly among the tumble planes. Moreover, at a higher velocity magnitude, a rotational-like flow structure was seen at maximum valve lift. Nevertheless, the higher velocities improve the air motion in terms of turbulent kinetic energy in the mid-cylinder plane compared to the other, mid-injector and mid-valve, planes.

A novel work from Vester et al. [61] studied the in-cylinder flow during the intake stroke in a unique test rig that allowed them to obtain 3D by reconstructing the images from a stereoscopic setup. They presented the possibilities that this technique offers to evaluate flow fields and turbulent structures and to study the air movement of swirl and tumble in an engine. A follow-up study from Vester et al. [39] investigated the effects on tumble and swirl of two-cylinder heads with geometrically different inlet ports using multi-planar PIV. The twin-port was designed to create more swirl, while the single port was intended to create more tumble. They showed that the twin-port cylinder presented fewer cycle-to-cycle variations and gave rise to a more coherent flow field consisting of a strong swirling motion.

While all the previous studies presented very valuable information, a parametric investigation to analyse the effect of the valve lift under different mass flows has not been carried out. Therefore, the main objective of this investigation is to study stationary intake flows in a simplified engine geometry by means of low-speed PIV with the aim of identifying significant flow phenomena such as intake flow jet oscillation, flow detachment at valve seat, or secondary flow and tumble due to high valve lift. Furthermore, POD analysis will be carried out to distinguish dominant, coherent, and turbulent flow structures. Although the boundary conditions are not too close to a real engine with a moving piston and transient flow, the results with a simplified geometry will give an insight into the flow phenomena and will provide valuable data to validate simulations under stationary conditions that can later be used to explore different conditions.

2. Fundamentals

In the following section, the necessary fundamentals are explained with the associated formulas, starting with the basics of turbulent flow and continuing with the derivation of the Reynolds stress transport equations, including some models for the Reynolds stress transport via diffusion through turbulence. The method for Proper Orthogonal Decomposition (POD) is then presented, followed by the quadruple decomposition of the POD modes into mean value and dominant, coherent, and turbulent structures.

2.1. Turbulent Flows

Turbulent flows are generally three-dimensional and spatial with eddies within a wide range of length scales; the motion of those eddies leads to an effective exchange (diffusion) of mass, momentum, and heat compared to laminar flows. The large eddies within the flow are dominated by inertia effects and are inviscid since viscous effects are negligible. The largest and most energetic eddies are affected by the mean flow direction and, therefore, the vortices are stretched in order to fulfil momentum conservation. As a result, their structure is anisotropic. The rotational rate of the eddies increases and their radius decreases. In this way, kinetic energy is transferred from larger to smaller eddies, the so-called energy cascade. The smallest eddies, which have the lowest energy and are isotropic, are dominated by viscous effects, and their energy is dissipated into thermal energy. This is the reason for the increased energy losses in turbulent flows [79].

For a three-dimensional flow situation, three covariances and three variances of all the velocity components can be calculated which cause six additional stress components (normal and shear stresses). These covariances and variances can be rewritten by multiplication, with the negative value of the density $\rho$ as the so-called Reynolds stress tensor ${\overleftrightarrow{\tau }}_{RS}$ in Equation (1). ${\overleftrightarrow{\tau }}_{RS}$ can be found in the ensemble-averaged momentum equation and describes the momentum transfer that occurs due to turbulent eddies.

$${\overleftrightarrow{\tau }}_{RS}=\left(\begin{array}{ccc}-\overline{\rho {u^{\prime }}_x^2}& -\overline{\rho u_x^{\prime }{u}_y^{\prime }}& -\overline{\rho u_x^{\prime }{u}_z^{\prime }}\\ -\overline{\rho u_x^{\prime }{u}_y^{\prime }}& \overline{\rho {u^{\prime }}_y^2}& -\overline{\rho u_y^{\prime }{u}_z^{\prime }}\\ -\overline{\rho u_x^{\prime }{u}_z^{\prime }}& -\overline{\rho u_y^{\prime }{u}_z^{\prime }}& \overline{\rho {u^{\prime }}_z^2}\end{array}\right)$$

(1)

The turbulent kinetic energy $k$, which is defined in Equation (2), is basically calculated from the trace of the Reynolds stress tensor ${\overleftrightarrow{\tau }}_{RS}$ and differs by a multiplicative factor of $-\frac{\rho }{2}$. $k$ is the turbulent energy per unit mass of the turbulence at a given location [79].

$$k=\frac{1}{2}\cdot \overline{\left({u^{\prime }}_x^2+{u^{\prime }}_y^2+{u^{\prime }}_z^2\right)}$$

(2)

The rate of dissipation of turbulent kinetic energy $\epsilon$ is defined in Equation (3), where $\epsilon$ is calculated from the elements of the fluctuating rate-of-strain tensor [80]. However, since the Kolmogorov scales are not resolved with the Particle Image Velocimetry (PIV) method, Equation (3) is only conditionally applicable to the PIV measurement results. To calculate $\epsilon$ from the PIV results with a given spatial resolution $\mathsf{\Delta }$, Wang et al. 2021 [81] proposed the Smagorinsky sub-grid model (SSM). Using this model, the dissipation of turbulent kinetic energy is approximated as the subgrid scale dissipation rate at the PIV resolution ${\epsilon }_{\mathrm{SSM}}$ shown in the Equation (4), where $C_{\mathrm{SSM}}$ is the Smagorinsky constant. This constant can be chosen in the range of $C_{\mathrm{SSM}}\approx 0.11\dots 0.21$ according to [81,82,83,84,85,86]. The Smagorinsky constant should be chosen depending on the PIV measurement conditions such as the interrogation window overlap and if the data is two- or three-dimensional [82].

$$\begin{gathered} \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\epsilon =2\nu \left(\overline{{\left(\frac{\partial u_x^{\prime }}{\partial x}\right)}^2}+\overline{{\left(\frac{\partial u_y^{\prime }}{\partial y}\right)}^2}+\overline{{\left(\frac{\partial u_z^{\prime }}{\partial z}\right)}^2}+\frac{1}{2}\overline{{\left(\frac{\partial u_x^{\prime }}{\partial z}+\frac{\partial u_z^{\prime }}{\partial x}\right)}^2}+\frac{1}{2}\overline{{\left(\frac{\partial u_x^{\prime }}{\partial y}+\frac{\partial u_y^{\prime }}{\partial x}\right)}^2} \\ +\frac{1}{2}\overline{{\left(\frac{\partial u_y^{\prime }}{\partial z}+\frac{\partial u_z^{\prime }}{\partial y}\right)}^2}\right) \end{gathered}$$

(3)

$$\begin{gathered} {\epsilon }_{\mathrm{S}\mathrm{S}\mathrm{M}}={\left(C_{\mathrm{S}\mathrm{S}\mathrm{M}}\cdot \mathrm{\Delta }\right)}^2{\left(2\cdot \overline{{\left(\frac{\partial u_x^{\prime }}{\partial x}\right)}^2}+2\cdot \overline{{\left(\frac{\partial u_y^{\prime }}{\partial y}\right)}^2}+2\cdot \overline{{\left(\frac{\partial u_z^{\prime }}{\partial z}\right)}^2}+\overline{{\left(\frac{\partial u_x^{\prime }}{\partial y}\right)}^2}+\overline{{\left(\frac{\partial u_x^{\prime }}{\partial z}\right)}^2} \\ +\overline{{\left(\frac{\partial u_y^{\prime }}{\partial x}\right)}^2}+\overline{{\left(\frac{\partial u_y^{\prime }}{\partial z}\right)}^2}+\overline{{\left(\frac{\partial u_z^{\prime }}{\partial x}\right)}^2}+\overline{{\left(\frac{\partial u_z^{\prime }}{\partial y}\right)}^2}+2 \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\cdot \left(\overline{\left(\frac{\partial u_x^{\prime }}{\partial y}\right)\left(\frac{\partial u_y^{\prime }}{\partial x}\right)}+\overline{\left(\frac{\partial u_x^{\prime }}{\partial z}\right)\left(\frac{\partial u_z^{\prime }}{\partial x}\right)}+\overline{\left(\frac{\partial u_z^{\prime }}{\partial y}\right)\left(\frac{\partial u_y^{\prime }}{\partial z}\right)}\right)\right)}^{3/2} \end{gathered}$$

(4)

2.2. Reynolds Stress Transport Equation

Starting from the vectorial mass specific momentum balance shown in Equation (5), the Reynolds stress transport equation can be derived. The next step is to split the velocity vector $\overrightarrow{u}$ and the pressure $p$ into their mean values and the fluctuations $\overrightarrow{u}=\overline{\overrightarrow{u}}+{\overrightarrow{u}}^{\prime }$ and $p=\overline{p}+p^{\prime }$, respectively, and to insert them into the specific momentum balance. After that, the whole equation is averaged, and the equation is multiplied with a fluctuation component ${u_i}^{\prime }$ with $i=x,y,z$. The compressible terms are not considered since the flow is assumed to be incompressible and the kinematic viscosity $\nu$ is assumed to be constant. This results in transport equations for all six components of the specific Reynolds stresses shown in Equation (6).

$$\frac{\partial \overrightarrow{u}}{\partial t}+\left(\overrightarrow{u}\cdot \nabla \right)\cdot \overrightarrow{u}=-\frac{1}{\rho }\nabla p+\nu {\nabla }^2\overrightarrow{u}+K$$

(5)

$$\begin{array}{cc}\frac{\partial \overline{u_i^{\prime }\cdot u_j^{\prime }}}{\partial t}+\nabla \left(\overline{u_i^{\prime }\cdot u_j^{\prime }}\cdot \overline{\overrightarrow{u}}\right)\hfill & \\ & =\nabla \left(\nu \nabla \overline{u_i^{\prime }\cdot u_j^{\prime }}\right)-\nabla \left(\overline{u_i^{\prime }\cdot u_j^{\prime }\cdot {\overrightarrow{u}}^{\prime }}\right)-\frac{1}{\rho }\left(\frac{\partial \overline{p^{\prime }\cdot u_i^{\prime }}}{\partial x_j}+\frac{\partial \overline{p^{\prime }\cdot u_j^{\prime }}}{\partial x_i}\right)\hfill \\ & -\left(\left(\overline{u_i^{\prime }\cdot {\overrightarrow{u}}^{\prime }}\cdot \nabla \right){\overline{u}}_j+\left(\overline{u_j^{\prime }\cdot {\overrightarrow{u}}^{\prime }}\cdot \nabla \right){\overline{u}}_i\right)+\frac{1}{\rho }\overline{p^{\prime }\left(\frac{\partial \overline{u_i^{\prime }}}{\partial x_j}+\frac{\partial \overline{u_j^{\prime }}}{\partial x_i}\right)}\hfill \\ & -2\nu \overline{\left(\nabla u_i^{\prime }\right)\cdot \left(\nabla u_j^{\prime }\right)}\hfill \end{array}$$

(6)

The specific Reynolds stress transport Equation (6) can be rewritten in a simplified form as shown in Equation (7), where $R_{ij}$ describes the rate of change of the specific Reynolds stress component in the control volume and $C_{ij}$ is the convective flux of specific Reynolds stress over the surfaces of the control volume. $D_{MT,ij}$, $D_{TT,ij}$, and $D_{PVF,ij}$ describe the diffusive flux due to molecular transport, turbulent transport, and pressure/velocity fluctuations, respectively. $P_{ij}$ defines the production term of specific Reynolds stress, ${\mathsf{\Pi }}_{ij}$ is equal to the pressure–strain correlation, and ${{\rm E}}_{ij}$ is the dissipation rate of specific Reynolds stress.

$$R_{ij}+C_{ij}=D_{MT,ij}-D_{TT,ij}-D_{PVF,ij}-P_{ij}+{\mathsf{\Pi }}_{ij}-E_{ij}$$

(7)

In a Reynolds stress model, the three terms $D_{TT,ij}$, $D_{PVF,ij}$, and ${\mathsf{\Pi }}_{ij}$ are modelled because they contain products of three different fluctuation components. Three different models of $D_{TT,ij}$, namely the models by “Daly & Harlow” (DH) [87], “Hanjalic & Launder” (HL) [88], and “Mellor & Herring” (MH) [89], are presented in [90,91] and shown in the Equations (8)–(10) with ${\tilde{\tau }}_{ij}=\overline{u_i^{\prime }{u}_j^{\prime }}$ being the mass specific Reynolds stresses. It must be noted that the constant $C_S$ should be chosen as $C_S=0.22$ for the model “Daly & Harlow” (DH) [87] and as $C_S=0.11$ for the “Hanjalic & Launder” (HL) [88] and “Mellor & Herring” (MH) [89] models.

$$D_{TT,DH,ijk}=\frac{\partial }{\partial x_k}\left[-C_s\frac{k}{\epsilon }{\tilde{\tau }}_{km}\frac{\partial {\tilde{\tau }}_{ij}}{\partial x_m}\right]$$

(8)

$$D_{TT,HL,ijk}=\frac{\partial }{\partial x_k}\left[-C_s\frac{k}{\epsilon }\left({\tilde{\tau }}_{im}\frac{\partial {\tilde{\tau }}_{jk}}{\partial x_m}+{\tilde{\tau }}_{jm}\frac{\partial {\tilde{\tau }}_{ik}}{\partial x_m}+{\tilde{\tau }}_{km}\frac{\partial {\tilde{\tau }}_{ij}}{\partial x_m}\right)\right]$$

(9)

$$D_{TT,MH,ijk}=\frac{\partial }{\partial x_k}\left[-\frac{2}{3}{C}_s\frac{k^2}{\epsilon }\left(\frac{\partial {\tilde{\tau }}_{jk}}{\partial x_i}+\frac{\partial {\tilde{\tau }}_{ik}}{\partial x_i}+\frac{\partial {\tilde{\tau }}_{ij}}{\partial x_k}\right)\right]$$

(10)

2.3. Proper Orthogonal Decomposition (POD)

In the field of statistics, the POD method is called Principal Component Analysis (PCA) and corresponds to a reduction of the dimensionality of a large data set into a smaller one which contains most of the information from the initial data set. POD is often used in fluid-dynamics to evaluate turbulent flows from PIV measurements. The idea behind POD is to decompose the random vector fields capturing the turbulent fluctuations into a set of deterministic functions. Each of those functions captures some amount of the total fluctuating kinetic energy of the flow [92].

2.3.1. The Direct POD Method

The main idea behind POD is to decompose the random vector field $u^{\prime }\left(x,t\right)$ into a set of deterministic spatial functions ${\mathsf{\Phi }}_{\mathrm{k}}\left(x\right)$ (spatial modes) and their random time coefficients $a_k\left(t\right)$ as shown in Equation (11) [92].

$$u^{\prime }\left(x,t\right)={\displaystyle \sum }_{k=1}^{\infty }{a}_k\left(t\right){\mathsf{\Phi }}_k\left(x\right)$$

(11)

First, the snapshot matrix $U$ shown in Equation (12) with the fluctuations is built. Therefore, the averaged velocity is subtracted from the velocities of the data set (PIV). $U$ is an $m\times n$ matrix with n = 23,474 data points (total number of velocity vectors) for the velocity components in x- and y-directions and m = 1000 instant times. There are 13,671 measuring points for which the two velocity components in the x- and y-directions are measured [92].

$$U=\left(\begin{array}{cccc}{u}_{11}^{\prime }& u_{12}^{\prime }& \cdots & u_{1n}^{\prime }\\ u_{21}^{\prime }& u_{22}^{\prime }& \cdots & u_{2n}^{\prime }\\ ⋮& ⋮& \ddots & ⋮\\ u_{m1}^{\prime }& u_{m2}^{\prime }& \cdots & u_{mn}^{\prime }\end{array}\right)$$

(12)

Then, the covariance matrix $C$ is computed with the diagonal elements being the variances and the off-diagonal elements being the covariances shown in Equation (13). $C$ is a symmetric $n\times n$ matrix [92].

$$C=\frac{1}{m-1}{U}^{T}U=\frac{1}{m-1}\left(\begin{array}{cccc}{\displaystyle \sum _{i=1}^m}{u^{\prime }}_{i1}^2& {\displaystyle \sum _{i=1}^m}{u}_{i1}^{\prime }\cdot u_{i2}^{\prime }& \cdots & {\displaystyle \sum _{i=1}^m}{u}_{i1}^{\prime }\cdot u_{in}^{\prime }\\ {\displaystyle \sum _{i=1}^m}{u}_{i2}^{\prime }\cdot u_{i1}^{\prime }& {\displaystyle \sum _{i=1}^m}{u^{\prime }}_{i2}^2& \cdots & {\displaystyle \sum _{i=1}^m}{u}_{i2}^{\prime }\cdot u_{in}^{\prime }\\ ⋮& ⋮& \ddots & ⋮\\ {\displaystyle \sum _{i=1}^m}{u}_{in}^{\prime }\cdot u_{i1}^{\prime }& {\displaystyle \sum _{i=1}^m}{u}_{in}^{\prime }\cdot u_{i2}^{\prime }& \cdots & {\displaystyle \sum _{i=1}^m}{u^{\prime }}_{in}^2\end{array}\right)=\left(\begin{array}{cccc}{c}_{11}& c_{12}& \cdots & c_{1n}\\ c_{21}& c_{22}& \cdots & c_{2n}\\ ⋮& ⋮& \ddots & ⋮\\ c_{n1}& c_{n2}& \cdots & c_{nn}\end{array}\right)$$

(13)

Next, the eigenvalue problem $C\mathsf{\Phi }=\lambda \mathsf{\Phi }$ is solved. The eigenvectors $\mathsf{\Phi }$, the so-called POD modes, are orthonormal and show the correlation of the fluctuations and how they move together. The time coefficient $A=\left(a_1,a_2,\dots ,a_k\right)$can be computed from $A=U\mathsf{\Phi }$. The eigenvalues $\lambda$ correspond to the amount of energy in each mode. Finally, the eigenvalues $\lambda$ and modes $\mathsf{\Phi }$ are ranked in descending order, where the first few modes contain the most energy. [92]

The original data within the snapshot matrix $U$ can be reconstructed with the sum of all the time coefficients $a_k\left(t\right)$ multiplied by the corresponding spatial mode ${\mathsf{\Phi }}_k\left(x\right)$, as shown in Equation (14) [92].

$$u^{\prime }\left(x,t\right)={\displaystyle \sum }_{k=1}^{1000}{a}_k\left(t\right){\mathsf{\Phi }}_k\left(x\right)$$

(14)

2.3.2. The Snapshot POD Method

The so-called snapshot POD works in essentially the same way as the direct POD method. However, time and space are interchanged in the algorithm. Instead of the deterministic spatial modes $\mathsf{\Phi }\left(x\right)$ and the random time coefficients $a\left(t\right)$, the deterministic temporal modes $\mathsf{\Phi }\left(t\right)$ with the random spatial coefficients $a\left(x\right)$are computed [92].

The snapshot method is faster since $m<n$ and the last eigenvalues of the direct method are zero anyway. The snapshot matrix $U_s$ can be computed from $U_s=U^T$ and is an $m\times n$ matrix. The covariance matrix $C_s$ is an $m\times m$ matrix. The eigenvalues ${\lambda }_s=\lambda$ are equal to the ones obtained from the direct method. Since the coefficients and modes of the snapshot method differ by a multiplicative factor from those of the direct method, the spatial modes $\mathsf{\Phi }$ can be calculated from the spatial coefficients ${\mathsf{\Phi }}_s$ by normalisation [92].

2.4. POD Quadruple Decomposition

According to Roudnitzky et al. 2006 [93], Wu et al. 2022 [94], Rulli et al. 2021 [95], Vu et al. 2012 [24], and Druault et al. 2005 [96], the turbulent structures can be split into dominant, coherent, and turbulent structures using POD. The velocity $u\left(x,t\right)$ is split according to Equation (15) into a mean value $\overline{u}\left(x,t\right)$and the fluctuating part $u^{\prime }\left(x,t\right)$. According to Equation (16), the fluctuating part is further split into the dominant, coherent, and turbulent structures.

$$u\left(x,t\right)=\overline{u}\left(x,t\right)+u^{\prime }\left(x,t\right)$$

(15)

$$u^{\prime }\left(x,t\right)=u^{\mathrm{dominant}}\left(x,t\right)+u^{\mathrm{coherent}}\left(x,t\right)+u^{\mathrm{turbulent}}\left(x,t\right)$$

(16)

The relevance index $R_i$ according to Wu et al. 2022 [94] shown in Equation (17) measures the similarity of two vector fields. Using the relevance index $R_i$, the vector fields of two neighboring POD modes are compared regarding their similarity. The cut-off mode between the dominant and the coherent structures can be determined by a certain threshold of $R_i$ between 0.9 and 1 [94]. The relevance index $R_i$ was calculated for the x- and y-velocity components in the whole flow field and a target value of 0.9935 was used based on the present results to find the first cut-off mode, which splits the dominant and the coherent structures.

Roudnitzky et al. 2006 [93] proposed a criterion to identify the cut-off number between the coherent and the turbulent structures; it assumes that homogeneous isotropic turbulent flow follows Gaussian properties. To define a suitable cut-off number, the Skewness $S_i$and the Flatness $T_i$(kurtosis coefficient) coefficients defined in Equations (18) and (19) are used. The Skewness $S_i$ defines the third standardised central moment of a distribution. For turbulent structures, a value of zero is expected. The Flatness $T_i$ measures the tailedness of a probability distribution and describes the shape of this distribution. For a normal distribution, $T_i$ is equal to three [95].

$$R_i=\frac{\overline{\overrightarrow{a_i{\mathsf{\Phi }}_i}·\overrightarrow{a_{i+1}{\mathsf{\Phi }}_{i+1}}}}{\overline{\overrightarrow{a_i{\mathsf{\Phi }}_i}}\cdot \overline{\overrightarrow{a_{i+1}{\mathsf{\Phi }}_{i+1}}}}$$

(17)

$$S_i=\frac{\overline{\overrightarrow{{\left(a_i{\mathsf{\Phi }}_i\right)}^3}}}{{\overline{\overrightarrow{{\left(a_i{\mathsf{\Phi }}_i\right)}^2}}}^{3/2 }}$$

(18)

$$T_i=\frac{\overline{\overrightarrow{{\left(a_i{\mathsf{\Phi }}_i\right)}^4}}}{{\overline{\overrightarrow{{\left(a_i{\mathsf{\Phi }}_i\right)}^2}}}^{2 }}$$

(19)

3. Methodology

A novel optical test bench was designed and commissioned to analyse the influence of the mass flow rate and the valve lift on the flow phenomena around a single intake valve, which is described together with the optical technique and the image processing within the following chapter. Further, it has been used to provide a reference data set for the validation of the Large Eddy Simulation (LES) model, which is presented in the next chapter.

3.1. Experimental Facility

To gain an insight into flow phenomena and to prove the feasibility of measurement and modelling techniques, a simplified geometry of Flex-OeCoS [97] was designed as Optical Test Bench 1 (OTB1), shown in Figure 1 on the right-hand side. For OTB1, the inlet valves geometry of Flex-OeCoS, shown in Figure 1 on the left-hand side with some simplifications, was used. The single inlet valve is centred in a rectangular flow section. The geometry features a straight cylinder roof geometry, while the valve is slightly tilted in both planes for 7° and 8° respectively. The air reaches the valve at an angle of 60° and exits the inflow section shortly after redirecting. The maximum mass flow of about 60 kg/h per valve was derived from measurements at Flex-OeCoS; the maximum valve lift is 5.5 mm.

In Figure 2 on the left-hand side, a cross-section through the optical chamber is shown. The inlet section was rotated so that the inlet valve was vertical with the upward flow direction. A single valve was used to avoid interaction with other valves and the walls, thus reducing complexity. The width of the section was 130 mm, which corresponded to the piston diameter of Flex-OeCoS. The angle of inclination of the wall corresponded to the cylinder head of Flex-OeCoS, and the valve lift could be adjusted manually. In Figure 2 in the middle and on the right-hand side, the whole setup is shown. The optical chamber was placed on a table with an upwards flow direction. The straight, 700 mm long inlet section was necessary to assure that the flow at the valve inlet was fully developed. A temperature measurement was placed at the beginning of the inlet section and a static pressure measurement was placed 350 mm further downstream, so it was unaffected by the thermocouple. The outlet section was 1 m long, so the flow around the valve was unaffected by the outlet. Additional pressure and temperature measurements were placed at the end of the optical chamber and at the outlet. Furthermore, mass flow was measured by a Coriolis mass flow meter.

Figure 3 shows the piping and instrumentation diagram (P&ID) of OTB1. As with the original setup of Flex-OeCoS, pressurised air from the FHNW-internal network at around 8 bar was used as the air supply. The mass flow, which was regulated with a valve, was measured using a Coriolis DN8 mass flow meter. A 20 L air vessel dampened possible pressure oscillations originating from the air supply. Liquid seeding was used as tracer particles, which was generated using a bi-pass air flow. The fog was then diluted by the remaining air and entered the optical chamber through the valve under investigation.

3.2. Optical Techniques and Image Processing

The main measurement technique used in this project was Particle Image Velocimetry (PIV) and will be described. The optical setup is shown in Figure 4 (left). For the low-speed experiments, a Litron Nano L 135-15 PIV Nd:YAG laser was used as a light source, generating two circular beams with 5 mm diameter at 532 nm wavelength. The pulse length was 6–9 ns with an average power of 135 mJ. The beams were firstly aligned to the optical chamber with mirrors and then focused by the first lens with focal length f = 1000 mm. The cylindrical lens with f = −30 mm then diverged the beam to the desired height before was is again collimated by the second cylindrical lens with f = 300 mm. The generated light sheet was then again focused perpendicular to the plane with a cylindrical lens with f = 500 mm. The resulting light sheet was about 0.5 mm thick and 30 mm high. The scattered light was recorded by a 12-bit CCD camera (sensicam qe), commercially known as Imager Intense by LaVision with a resolution of 1376 × 1040 px. The camera was equipped with a 55 mm f/2.8 Micro-NIKKOR objective. The chosen field of view resulted in a resolution of approximately 11 px/mm. In double-frame mode, this camera can operate at a maximum frequency of 4 Hz.

The flow must contain tracer particles. For example, in gaseous flows, boron nitride particles are often used as solid particles or Di-Ethyl-Hexyl-Sebacate (DEHS) as liquid particles. As previously shown in Figure 3, the feed air was enriched with liquid seeding particles of about 1 µm generated by a Laskin nozzle aerosol generator PivPart45 by PIVTEC. As seeding medium, a SAFEX^® inside fog fluid flash/reflex B50, which consists of an aerosol solution of non-toxic highly pure polyols, was typically used in theatres, was established. The small particles are necessary so they can follow the flow while also not being so small that the scattered light is insufficient. The partial flow through the aerosol generator can be adjusted to modify particle density. The advantages of liquid seeding as opposed to solid seeding are easier aerosol generation, easier cleaning, and non-toxicity. On the downside, they scatter less light than solid particles [24].

A laser generated two laser-pulses with a fairly short time difference δt, depending on flow velocity. These pulses were spread into light sheets by suitable optical components to illuminate a 2D plane of particles in the flow, which scattered their light. This scattered light was recorded by a camera, as shown in Figure 5 on the left with “exposure 1” and “exposure 2”. The two images were then divided into interrogation spots, typically of 16 × 16 px, 32 × 32 px or 64 × 64 px size. An algorithm then compared the two exposures to find the dominant flow direction by calculating the cross-correlation between the interrogation spot of exposure 1 and the proximate interrogation spots of exposure 2. By doing so for the whole image, a flow field could be derived, as shown in Figure 5.

The raw double-images were processed with the commercially available software DaVis 10.1 by LaVision. An intensity normalization filter was used to eliminate local background noise before calculating the local flow velocities with a multi-pass, cross-correlation algorithm. One initial pass with a resolution of 32 × 32 px and 50% overlap was performed to find the highest velocities before applying four final passes with 16 × 16 px and 75% overlap. The overlap led to a higher resolution despite choosing a coarser interrogation window size, thus not losing dynamic range. The 2D vector field had a vector spacing of 4 mm. The resulting field of view is shown in Figure 4 (right); it was about 100 mm high and 120 mm wide. Since the laser sheet was only about 30 mm high, the field of view was cropped on the top and bottom since no scattered light could be observed.

4. LES Model Description

4.1. CFD Setup

The simulations were conducted with the CFD Solver Converge v3.0. By default, the near wall treatment was resolved by a high-y+ approach, which resolved the flow field outside of the viscous sublayer and the buffer zone and modelled the quantities inside the named zones. As an LES model the Dynamic Structure Model was used, signified by a one-equation term.

A density-based PISO solver was used for the Navier-Stokes-Equations. As a solver for the governing equations (Momentum, Pressure, Density, and Energy) CONVERGE offers the pointwise successive over-relaxation (SOR) algorithm, which is less computationally expensive and is recommended by CONVERGE for engine simulations. The iteration limitations and convergence tolerances from

Table 1. Overview of SOR settings for the governing equations regarding convergence tolerance, minimum and maximum iterations, and SOR Relaxation value.

	Convergence Tolerance	Min. Iterations	Max. Iterations	SOR Relaxation
Momentum	1 × 10−5	0	50	1.0
Pressure	1 × 10−8	2	500	1.1
Density	1 × 10−4	0	2	1.0
Energy	1 × 10−4	0	2	1.0

were used for the SOR algorithm.

To resolve the flow field both spatial and timewise accurately to the defined grid, a convective CFL-number of $CF{L}_{\mathrm{conv}}=1$ was chosen. This ensured that by the convective flow for each grid point the flow field quantities would be calculated.

As inflow boundary condition, a mass flow was chosen which was measured at the test bench. Inflow temperature was set to the measured ambient temperature value $T_{\mathrm{in}}=293 \mathrm{K}$, and air was chosen as medium. As outflow boundary condition, the measured ambient pressure $p_{\mathrm{out}}=1.013 \mathrm{bar}$ and ambient temperature $T_{\mathrm{out}}=293 \mathrm{K}$ were used. Initial conditions were a non-moving fluid with ambient temperature $T_{\mathrm{ini}}=293 \mathrm{K}$ and ambient pressure $p_{\mathrm{ini}}=1.013 \mathrm{bar}$.

For the determination of the quasi-stationary time domain, the maximum, minimum, and mean (spatial) pressure of the ‘Region of Interest’ as shown in Figure 6 (left) would be analysed. The maximum and minimum are single values out of the domain, while the mean value is a volume average of the whole domain. The ‘Region of Interest’ was evaluated because it is the domain where the intake valves are and, therefore, the scope of these investigations.

The results are shown in Figure 7 for one configuration, OP_2_39.6, which will be further explained in the following chapter. This configuration is representative of all other configurations regarding the definition of a transient time domain and a quasi-stationary domain.

As seen for the pressure curves in Figure 7 (left), the maximum and minimum values (continuous line) reached their extrema around 1 ms and then reached a plateau after approximately 10 ms, where they stayed until the end of the simulation at 40 ms. The mean pressure (continuous line) remained constant from the beginning.

The moving average weighted the pressure development of the minimum and maximum and reached a final value with proceeding physical time, which shows that already after 10 ms a quasi-stationary time domain was obtained. Therefore, time-averaging after this point was valid and the independence of the transient time domain was verified.

Figure 7 (right) shows the relative error between the actual pressure values and the moving average value for maximum, minimum, and mean pressure. After 10 ms the relative error of these values was within a range of +/− 1.5% tolerance, which is acceptable. Furthermore, it has to be said that due to the flow separations around the intake valves the pressure signal oscillated, which is reasonable for an LES which resolves the fluctuations spatially and temporally. Therefore, the fluctuations, for example, in the minimum pressure are no sign of a not-converged simulation, but rather a feature of LES.

4.2. Grid Study

To minimise the numerical error, a grid study with different grid sizes was conducted. The aim was to define a grid where the numerical error is small enough that it does not influence the overall result of the simulation. In this case, the discharge coefficient of the engine valve in Equation (20), which is a common non-dimensional quantity for valves, was chosen as an evaluation criterion.

$$C_d=\frac{\dot{m}}{A\sqrt{2\rho \mathsf{\Delta }p}}$$

(20)

The CFD-Domain was divided into two different mesh regions (see Figure 6), where the upper region close to the outflow had a coarser mesh, e.g., bigger cell size, because the gradients of pressure and velocity were not strong due to laminar flow behaviour and therefore no refinement was necessary. In the ‘Mesh Refinement Region’ upstream the valve gap, inside the valve gap, and in the wake of the flow separation, stronger gradients were expected; therefore, this region had to be refined to resolve the pressure and velocity fields appropriately. The cell size of the ‘Mesh Refinement Region’ was a quarter of that of the ‘Base Grid Size’ region.

In

Table 2. Test matrix for different grid sizes.

Case Name	Coarse	Medium	Medium-Fine	Fine	Extra-Fine
Base grid size	8 mm	4 mm	2.5 mm	2 mm	1.5 mm
Refinement region	2 mm	1 mm	0.625 mm	0.5 mm	0.375 mm
Total cell number	~440,000	~3,300,000	~13,900,000	~27,000,000	~64,500,000

, the investigated mesh sizes are shown. The ‘Fine’, ‘Medium’, and ‘Coarse’ meshes are signified by a duplication of the base grid size starting from the ‘Fine’ mesh. With these three simulations a final decision about the mesh quality could not be attained, because with a finer mesh an asymptotic behaviour from the discharge coefficient is expected but was not observed with these three simulations. The simulations of the ‘Extra-Fine’ and ‘Medium-Fine’ meshes were conducted to achieve a distinctive trend, which is shown by the trendline in Figure 8 on the left-hand side. The difference in the discharge coefficient between the ‘Fine’ and ‘Extra-Fine’ meshes was small, and with further refinement the accuracy would not increase in a way that would justify the computational effort.

In Figure 8 the right plot shows the error of the discharge coefficient relative to the ‘Extra-Fine’ discharge coefficient. For the ‘Fine’ grid the relative error is around 0.7%, which is below 1% and therefore in a tolerable range. As the relative error shows asymptotic behaviour, it is not expected that further refinement would result in considerably higher error for the finer mesh. Hence, the grid parameters of the ‘Fine’ mesh were chosen for the parameter study because it has sufficient accuracy within a reasonable computational expense.

5. Results and Discussion

5.1. Flow Behaviour for the Different Operating Points

Figure 9 shows the ensemble averaged velocity fields over 1000 single velocity fields of the measurements and the streamlines (in white) for different operating points. The nomenclature is chosen as in the following example: OP_4_39.6, where the first number (4) corresponds to the valve lift in millimetres and the second number (39.6) to the mass flow in kilograms per hour.

As to be expected, with increasing mass flow rate, the velocity magnitude increased at constant valve lift, while the dominant flow structures merely changed. All the velocity fields at 4 mm valve lift for different mass flows showed a straight jet exiting on the right side of the valve. On the left side, there were two flow phenomena observed. Firstly, the jet exited somewhat straight as well, but with a higher curvature compared to the right side. Secondly, the flow generated a recirculation zone right above the valve. The difference in flow behaviour between the left and the right jet and the asymmetrical flow is caused by the inclined surface at the valve exit, as shown in Figure 2.

With increasing valve lift at a constant mass flow rate, the maximum velocity decreased since the valve gap was increasing, though not only the velocity but also the dominant flow structures changed. At 1 mm and 2 mm, an almost symmetrical velocity field could be observed. At 3 mm and 4 mm, an asymmetry could be seen between the two jets, where the left jet showed a larger recirculation zone towards the middle of the valve compared to the right jet. Further increasing the valve lift resulted in a rising recirculation zone at the left jet, with a straight jet on the right and an alternating flow with recirculation towards the middle on the left. With 5 mm and 6 mm, this recirculating flow became more dominant while also getting flatter. At maximum valve lift (6 mm), the recirculation appeared almost perpendicular to the valve axis, which was supposed to cause an overall tumble motion within the flow field above the valve.

Using the operating points from Figure 9, an uncertainty study of the deviations resulting from the PIV method was carried out according to Wieneke [98] and Sciacchitano & Wieneke [99] using the DaVis software. The corresponding uncertainties are shown in Figure 10. For this purpose, the median of all relative percentage deviations (absolute deviation divided by the respective individual velocity fields) was determined for each operating point. The median was chosen because the calculation occasionally involved division by zero, which results in infinitely large values. This would result in unsuitable results when averaging. The results in Figure 10 show small deviations where the velocity magnitude is large and vice versa.

5.2. Validation of the LES

To validate the LES results with the results from the PIV measurements, the ensemble/time averaged velocity fields were used. Note that the colour coding for the velocities is adjusted in terms of minimum and maximum values to have the best resolution in terms of the comparison between the LES and the PIV experiments.

The variation in the mass flow rate from 19.8 kg/h to 59.5 kg/h with a constant valve lift of 4 mm is presented in Figure 11 for the simulation (LES) and experimental (PIV) results. These results show that the time-averaged velocity of the simulation agrees very well with the PIV results. The structure of the flow hardly changed depending on the mass flow; only the global velocity scaled with the mass flow rate. It is noticeable that the velocity fields of the simulation results are less smooth compared to those of the PIV results. The reason for this is the averaging over a shorter period of 40 ms.

Figure 12 shows a comparison between the LES and the PIV results for the valve lift variation with valve lifts ranging from 1 mm to 6 mm at a mass flow rate of 39.6 kg/h. Especially for valve lifts between 4 mm and 6 mm, a good agreement between LES results and PIV measurements can be seen. The simulation accurately captured the behaviour of the left jet as the valve lift increased. Even with a valve lift of 2 mm, the LES showed qualitatively the same behaviour as the PIV, with the left and right jets being slightly curved towards the middle, resulting in symmetrical flow behaviour. With valve lifts of 1 mm and 3 mm, the deviation was somewhat larger. With a valve lift of 1 mm, the experiment showed stronger fluctuations in the two jets (right and left) compared to the velocity field of the LES. The reason for better agreement in the larger valve lifts is assumed to be that the wall had a large influence on the flow with small valve lifts. Furthermore, with the fixed cell size around the valve, the number of cells in the valve gap was lower for smaller valve lifts. This was aggravated by the fact, that for the smallest valve gaps, the velocity was the highest and therefore the gradients were greater, which led to bigger numerical errors. Overall, a good agreement between the LES results and PIV measurements can be observed.

In addition to the more detailed validation of the flow structure around the valve, a global validation was also carried out by comparison of the measures and by the LES methodology’s calculated valve discharge coefficient according to Equation (20).

Figure 13 shows the variation in the mass for the operating points “OP_4_19.8”, “OP_4_39.6”, and “OP_4_59.5”. For the operating points “OP_1_39.6”, “OP_2_39.6”, “OP_3_39.6”, “OP_4_39.6”, “OP_5_39.6”, and “OP_6_39.6” the valve gap variation is shown with valve lifts ranging from 1 mm to 6 mm. An experimental uncertainty study according to “Guide to the Expression of Uncertainty in Measurement” (GUM) [100] was performed including the type A and type B uncertainties of the pressure, temperature, and mass flow rate quantities. The combined uncertainty of the discharge coefficient $C_d$ was then calculated and is shown with error bars in Figure 12. The overall physical behaviour is represented correctly. For the variation in the mass flow the discharge coefficient according to Equation (20) remained approximately constant because the discharge coefficient describes the pressure loss of a valve contour and is therefore driven by the shape and not the mass flow. For the valve lift variation the shape changed for each operating point, and therefore changing values could be observed. For the smallest valve gaps the pressure loss was the highest, resulting in the smallest discharge coefficient. With greater valve lift the pressure loss was reduced and the discharge coefficient increased to an accumulated value where the valve disc was not affecting the flow. This can be seen for the valve lifts of 5 mm and 6 mm, where the discharge coefficient remains roughly constant. Besides the plausibility check of the physical behaviour, the comparison of LES and experimental data shows good agreement. Only for the operating points “OP_4_39.6” and “OP_4_59.5” is there a bigger difference between simulation and experiment; however, this is still considered to be in a reasonable range.

5.3. Turbulent Kinetic Energy

The turbulent kinetic energy $k$ was computed from the velocity field for both the LES results and the PIV measurements to enable a statement about how much turbulent kinetic energy $k$ can be resolved in LES and to analyse the turbulence within the flow field. In Figure 14 $k$ at 4 mm valve lift and a mass flow rate of 39.6 kg/h of both PIV and LES is shown. There is an adequate agreement in mean turbulent kinetic energy between the simulation and the experimental results. There were two shear layers at the right- and left-hand side of the two jets. The turbulent kinetic energy was highest within these shear layers. This seems plausible according to the energy cascades since turbulence is produced by velocity gradients. It turns out that $k$ was larger in the right jet compared to the left one, which can be attributed to the intake port geometry being inclined to the right.

To quantify the fraction of turbulent motion resolved within the LES model, the modelled turbulent kinetic energy at 4 mm valve lift and a mass flow rate of 39.6 kg/h was analysed. Furthermore, the percentage amount of modelled turbulent kinetic energy compared to the total kinetic energy (modelled and resolved) was computed and is shown in Figure 15. The result shows that about 99.9% of turbulence was resolved within the time-dependent velocity field and there was only a small amount (about 0.1%) modelled within LES. It can be concluded that the LES is of good quality. The aim of LES is to resolve 80% of the turbulence, and a maximum of 20% should be modelled with the subgrid model.

5.4. Proper Orthogonal Decomposition (POD) of PIV Measurements and Turbulent Structures

The coefficients of Skewness $S_i$ and normalised Flatness $T_i$ as well as the Relevance index $R_i$ were used to find the cut-off modes between the dominant, coherent, and turbulent structures for quadruple decomposition. Those coefficients are shown in Figure 16 for all the operating points. The Flatness coefficient $T_i$ is normalised by subtracting the target value of 3 from it and dividing it by the maximum value of the first peak to be able to carry out quadruple decomposition for all nine operating points uniformly. Due to the large fluctuations, all curves were smoothed with a moving average, using a sliding mode window of 10, 200, and 50 for $S_i$, $T_i$, and $R_i$, respectively. According to this, the cut-off modes of the various operating points with different mass flow rates and the valve lifts were found as shown in

Table 3. Three cut-off modes of the POD quadruple decomposition for all the operating points measured using PIV.

OP	Mass Flow Rate in kg/h	Valve Lift in mm	First Cut-off Mode	Second Cut-off Mode	Third Cut-off Mode
OP_4_19.8	19.8	4	45	275	724
OP_4_39.6	39.6	4	41	237	672
OP_4_59.5	59.5	4	29	495	661
OP_1_39.6	39.6	1	26	543	1219
OP_2_39.6	39.6	2	26	264	697
OP_3_39.6	39.6	3	37	569	1163
OP_5_39.6	39.6	5	48	328	634
OP_6_39.6	39.6	6	40	307	672

.

After ranking and sorting the POD modes according to their energy, the energy per mode result is shown in Figure 17 on the left-hand side. Besides this, the cumulative curve of percent energy per mode is visualised, showing that, for example, the first 200 modes contained 80% of the turbulent kinetic energy. This agrees with the theory of the energy cascades and the behaviour of turbulent flow. Figure 18 shows an example of the POD modes, and the associated time coefficients are shown in Figure 17 on the right-hand side for 4 mm and 39.6 kg/h. The time coefficients show a decrease in the amplitude with higher POD modes due to the smaller amount of energy within the smaller structures at higher modes.

The exemplary POD modes in Figure 18 are a measure of the correlation at different points in the flow. This means that the fluctuating components of the velocity within the red region in one mode were the most strongly correlated. The flow structures within the POD modes were random due to the behaviour of turbulence. However, each of the modes visualised a certain proportion of spatial correlation and contributed its amount of turbulent kinetic energy to the flow. Summing up all the modes multiplied by their time coefficients leads to the original velocity fluctuations. Because of that, POD can be used to build a reduced order model of the velocity fluctuations. Since the small turbulent structures were of interest within the present study, a reduced order model was not considered here.

Figure 19 shows the turbulent kinetic energy $k$ of all modes of the operating points as well as the cut-off into dominant, coherent, turbulent, and remaining modes. It is shown that the division into the different structures was almost equal for all the operating points and was therefore universal. The largest part of the turbulent kinetic energy can be assigned to the dominant structures. The energy decreased over the dominant, coherent, turbulent, and remaining modes, which agrees well with theory. The mass flow variation in Figure 19 on the left-hand side shows that the turbulent kinetic energy increased quadratically with increasing mass flow rate. Furthermore, total $k$ decreased quadratically with increasing valve lift, which can be seen in Figure 19 on the right-hand side.

By summing up the POD modes over all the dominant, coherent, and turbulent modes separately by multiplying the modes with their corresponding time coefficients, the representation in Figure 20 is obtained. However, it should be mentioned that this is only a qualitative comparison and that it is difficult to compare the structures of different operating points due to false vectors in the individual velocity fields, the different image sections, and the different total number of modes of the various operating points. However, a qualitative analysis can be made regarding the behaviour of the turbulent flow. In Figure 20, the sum of dominant structures, which are large structures that contain a lot of energy, is shown on the left-hand side. The dominant structures were particularly present in the vicinity of the two jets, because the turbulent energy was generated by velocity gradients. It turns out that the dominant structures were more present in the right jet compared to the left one. The reason is assumed to be the intake port, which was inclined to the right-hand side. The coherent structures in the middle of Figure 20 are smaller structures compared to the dominant ones and have less energy. For reasons of momentum conservation, there are many more coherent structures compared to the dominant ones, because they result from the decay of large structures. The coherent structures were also more dominant in the vicinity of the jets than in the rest of the flow area, but they could be observed throughout the flow area. In general, dominant structures were formed by shear layers and, because of the vortex decay, the structures were getting smaller with higher modes were distributed more throughout the entire flow area. The turbulent structures are shown on the right-hand side in Figure 20. There were many small structures with little energy, which were homogeneously distributed in the entire flow area.

5.5. Reynolds Stress Transport at 39.6 kg/h and 4 mm

The results of the Reynolds stresses ${\tau }_{xx}$, ${\tau }_{yy}$ and ${\tau }_{xy}$ of the operating point at 4 mm and 39.6 kg/h are shown in Figure 21. The Reynolds stresses describe the momentum exchange related to the turbulent fluctuations of the velocity components. The two components of the Reynolds stress tensor ${\tau }_{xx}$ and ${\tau }_{yy}$ show that the velocity fluctuations in the x and y directions on the bottom and top of the jets, respectively, led to an increased turbulent momentum exchange. It is obvious that the two components together form the turbulent kinetic energy, which was maximal in the shear layers of the jet. The Reynolds stress component ${\tau }_{xy}$ is proportional to the covariance between the x and y velocity components. According to Figure 21, the fluctuations in velocity in the x and y directions correlated the most within the jets.

Figure 22 shows the different terms $T_{ij}$, $C_{ij}$, $D_{MT,ij}$, $D_{TT,ij}$, $P_{ij}$, $E_{ij}$ of the specific Reynolds stress transport equation from Equation (7) within the different rows from top to bottom. The columns in Figure 22 show the different components of ${\tau }_{xx}$, ${\tau }_{yy}$, and ${\tau }_{xy}$ from left to right, where the terms involving the third dimension $z$ are set to zero. The two terms $D_{PVF,ij}$ and ${\mathsf{\Pi }}_{ij}$ involving the pressure $p$ could not be calculated because there is no spatially resolved pressure data available from the PIV measurements. The transient terms $T_{ij}$ are presented in the first row in Figure 22 and are in all three equations (columns) with a maximum amount of $5 {\mathrm{m}}^2/{\mathrm{s}}^3$ very small compared to the other terms, which are in the order of magnitude of ${10}^3 {\mathrm{m}}^2/{\mathrm{s}}^3$ or larger. The reason for that is the large time steps of 0.25 s between the individual velocity fields, which are chosen to ensure statistical independence. Therefore, it is assumed that the transient terms $T_{ij}$ are greatly underestimated. Compared to $T_{ij}$, the diffusion due to molecular transport $D_{MT,ij}$ and the dissipation rate $E_{ij}$ in the third and sixth rows in Figure 22 are much larger. However, $D_{MT,ij}$ and $E_{ij}$ are also negligible compared to the other terms, because those scale linearly with the kinematic viscosity $\nu$, which is very small for air with $\nu =1.6\cdot {10}^{-5}\frac{{\mathrm{m}}^2}{\mathrm{s}}$. Basically, the convective flux term $C_{ij}$, the diffusion to turbulent transport $D_{TT,ij}$ and the production/generation term $P_{ij}$ play an important role within the specific Reynolds stress Equation (7). $C_{ij}$, $D_{TT,ij}$, and $P_{ij}$ were particularly large in the vicinity of the left and the right jets of the valve, and the Reynolds stresses were produced and transported by convection between the two jets, where the vortex was formed. However, the specific Reynolds stress transport equations cannot be approximated using the three terms $C_{ij}$, $D_{TT,ij}$, and $P_{ij}$ because calculating $C_{ij}-D_{TT,ij}-P_{ij}$ does not add up to value of zero. Apparently, some of the terms of $T_{ij}$, $D_{PVF,ij}$ and ${\mathsf{\Pi }}_{ij}$ are important to satisfy the equation.

The same results as in Figure 22 are presented in Figure 23 based on the CFD simulation results. Basically, there is a good qualitative agreement between the terms of the Reynolds stress transport equation from simulation and the experimental PIV results, and the order of magnitude also agrees well. However, the results in Figure 22 calculated from the PIV data look smoother compared to the simulation results because there were not enough timesteps used for averaging in the latter. Due to the sufficiently small timestep of $5\times {10}^{-5} \mathrm{s}$ in the LES, the transient term with a maximum value of ${10}^5 {\mathrm{m}}^2/{\mathrm{s}}^3$ is also important for the Reynolds stress transport equation. It was present in the entire flow area and was particularly dominant near the two jets. However, the LES data using smaller time steps are not statistically independent, which means that no POD can be performed with the LES data.

For Reynolds stress turbulence modelling, the three terms of $D_{PVF,ij}$, ${\mathsf{\Pi }}_{ij}$, and $D_{TT,ij}$ involving a product of three fluctuations are modelled. There is no information about $D_{PVF,ij}$ and ${\mathsf{\Pi }}_{ij}$, and, because of the unknown pressure field $p$ from the PIV data, modelling of the diffusion due to turbulent transport $D_{TT,ij}$ is analysed in Figure 24. To do so, the dissipation rate of the turbulent kinetic energy is needed. It is calculated from the two Equations (3) and (4); both results are shown in Figure 25. The dissipation rate $\epsilon$ is shown on the left-hand side and the dissipation rate using the Smagorinsky sub-grid model (SSM) ${\epsilon }_{SSM}$ with $\mathsf{\Delta }=0.3763 \mathrm{m}$ and $C_{\mathrm{SSM}}=0.21$ is shown on the right-hand side in Figure 25. Both $\epsilon$ and ${\epsilon }_{SSM}$ are qualitatively identical, but ${\epsilon }_{SSM}$ is larger by one order of magnitude compared to $\epsilon$. Because of not resolving the Kolmogorov scales, the dissipation rate of the turbulent kinetic energy is underpredicted. This is one of the uncertainties of modelling the diffusion due to turbulent transport $D_{TT,ij}$. In Figure 24, $D_{TT,ij}$ calculated from the PIV data is shown in the first row; the calculations of the models by “Daly & Harlow” (DH) [87], “Hanjalic & Launder” (HL) [88], and “Mellor & Herring” (MH) [89] $D_{TT,DH,ij}$, $D_{TT,HL,ij}$, and $D_{TT,MH,ij}$ are presented in the second to the fourth row. On the one hand, there is a good qualitative agreement within the flow pattern between the three models and with the result of $D_{TT,ij}$. This leads to the conclusion that all three models are equally suitable to predict the diffusion due to turbulent transport within the flow around an intake valve. On the other hand, $D_{TT,ij}$ is larger by a factor of about four compared to the models $D_{TT,DH,ij}$, $D_{TT,HL,ij}$, and $D_{TT,MH,ij}$. One reason for this might be an underprediction of $\epsilon$, which would lead to an overprediction of $D_{TT,DH,ij}$, $D_{TT,HL,ij}$, and $D_{TT,MH,ij}$. Another reason could be the two-dimensional flow calculations and the terms involving the third dimension $z$, which are omitted with two-dimensional data. Furthermore, the Smagorinsky constant of $C_{\mathrm{SSM}}=0.21$ is based on an assumption, and the exact value is not known. However, all three models can predict the diffusion due to turbulent transport $D_{TT,ij}$ very well in a qualitative manner.

6. Summary, Conclusions, and Outlook

A novel test bench using a single inlet valve of the Flex-OeCoS engine was designed to perform Particle Image Velocimetry (PIV) measurements in addition to the development of a Large Eddy Simulation (LES) model. Experiments and simulations were carried out for different valve lifts and mass flows and the simulation results were validated with the experimental results.

The results show that the mass flow has no influence on the flow structure, but that the velocity magnitude increases with an increasing mass flow rate. Compared to the mass flow variation, the valve lift affects the flow structure, starting from an almost symmetrical flow situation at 1 mm valve lift to an asymmetrical flow situation at 6 mm, where the left jet leads to strong oscillations and a recirculation zone above the valve. The corresponding LES results capture the flow structure very well and agree with the experimental results. On top of that, the turbulent kinetic energy shows that the LES is of high quality, and only about 0.1% of the turbulent kinetic energy is modelled, whereas 99.9% is resolved. Furthermore, the turbulent structures are investigated with POD quadruple decomposition, where the POD modes are split into the dominant, coherent, and turbulent structures. The dominant modes correspond to the large flow structures with a high amount of energy, which are particularly present within the two jets on the left and the right side of the valve. The small structures with a low amount of energy are homogenously distributed within the whole flow field. According to the amount of energy within the different structures, the division is almost equal (universal) for all the operating points. However, the turbulent kinetic energy is not equal, increasing quadratically with an increasing mass flow rate and with decreasing valve lift. Analysing the Reynolds stress transport equation, the most important terms are the rate of change $T_{ij}$, the convective flux $C_{ij}$, the turbulent diffusive flux $D_{TT,ij}$, and the production term $P_{ij}$. The molecular diffusive flux $D_{MT,ij}$ as well as the dissipation rate $E_{ij}$ can be neglected since their order of magnitude is about 100 times smaller. Moreover, the turbulent diffusive flux can be modelled, and three different models show a qualitative agreement between each other and with the experimental result.

In a future study, the existing test rig will be adapted for a dual valve study under steady state and transient conditions (moving valve) for similar boundary conditions to investigate the interaction between the flow around two valves.