2.1 Wall-boundary conditions

As our study focuses on a geometry-driven flow, we shall briefly review the topic of wall-boundary conditions. In order to complete the advection step of the LB algorithm, the set of populations $\bar{D}$ advected from the solid into the fluid region are unknown and need to be specified. To that end, a variety of realizations can be found in the literature ranging from simple bounce-back boundary conditions using a crude staircase approximation to more sophisticated formulations trying to incorporate the curvature of the object. It is clear that for under-resolved simulations, the staircase approximation is not the best option as the geometry is not accurately represented. For curved boundary conditions, a common approach is to employ interpolation or extrapolation schemes onto the populations (see, e.g., Guo, Zheng & Shi Reference Guo, Zheng and Shi2002; Lallemand & Luo Reference Lallemand and Luo2003; Wang et al. Reference Wang, Shu, Teo and Wu2015). However, their usage is limited to low Reynolds number flows as spurious shocks are triggered at the boundary for turbulence simulations (Lammers et al. Reference Lammers, Beronov, Volkert, Brenner and Durst2006; Spasov, Rempfer & Mokhasi Reference Spasov, Rempfer and Mokhasi2009).

An alternative, recently proposed by Dorschner et al. (Reference Dorschner, Chikatamarla, Bösch and Karlin2015), is to impose the boundary conditions not on the highly fluctuating populations directly but rather on the slow varying moments, namely density, momentum and pressure. Here, the missing populations are approximated by an analogue of Grad’s distribution function, which results in a parametrization of the distribution in terms of relevant moments, including not only locally conserved quantities but also other pertinent moments. A derivation using maximum entropy or quasi-equilibrium considerations may be found in the work of Gorban & Karlin (Reference Gorban and Karlin2004). In the athermal case, it is sufficient to include the pressure tensor $\boldsymbol{{\it\Pi}}$ besides the conserved quantities, yielding

(2.13) $$\begin{eqnarray}\displaystyle f_{i}^{\ast }({\it\rho},\boldsymbol{u},\boldsymbol{{\it\Pi}})=W_{i}\left[{\it\rho}+\frac{{\it\rho}}{c_{s}^{2}}\boldsymbol{c}_{i}\boldsymbol{\cdot }\boldsymbol{u}+\frac{1}{2c_{s}^{4}}\left(\boldsymbol{{\it\Pi}}-{\it\rho}c_{s}^{2}\unicode[STIX]{x1D644}\right)\boldsymbol{ : }(\boldsymbol{c}_{i}\otimes \boldsymbol{c}_{i}-c_{s}^{2}\unicode[STIX]{x1D644})\right], & & \displaystyle\end{eqnarray}$$

where the pressure tensor $\boldsymbol{{\it\Pi}}$ is approximated by

(2.14) $$\begin{eqnarray}\displaystyle \boldsymbol{{\it\Pi}}=\boldsymbol{{\it\Pi}}^{eq}+\boldsymbol{{\it\Pi}}^{neq}, & & \displaystyle\end{eqnarray}$$

with

(2.15) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{{\it\Pi}}^{eq}={\it\rho}c_{s}^{2}\unicode[STIX]{x1D644}+{\it\rho}\boldsymbol{u}\otimes \boldsymbol{u}, & \displaystyle\end{eqnarray}$$

(2.16) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{{\it\Pi}}^{neq}=-\frac{{\it\rho}c_{s}^{2}}{2{\it\beta}}(\boldsymbol{{\rm\nabla}}\boldsymbol{u}+\boldsymbol{{\rm\nabla}}\boldsymbol{u}^{\dagger }). & \displaystyle\end{eqnarray}$$

The macroscopic variables to be used in $f_{i}^{\ast }$ (2.13) are specified by introducing the concept of target values. The target velocity is obtained using an interpolation scheme including the wall velocity and the velocity at the neighbouring fluid nodes. The target density is evaluated based on the bounce-back density. An additional term needs to be included into the target density in the case of a moving wall to account for the mass swept by the obstacle moving with the velocity $\boldsymbol{u}_{w}$ . This may be derived by introducing a forcing term necessary for the displacement and amounts to

(2.17) $$\begin{eqnarray}\displaystyle {\it\rho}_{dyn}=\mathop{\sum }_{i\in \bar{D}}6W_{i}{\it\rho}\boldsymbol{c}_{i}\boldsymbol{\cdot }\boldsymbol{u}_{w,i}. & & \displaystyle\end{eqnarray}$$

Finally, the pressure tensor is determined using a finite difference scheme for $\boldsymbol{{\rm\nabla}}\boldsymbol{u}$ and the velocity values from the previous time step. Note that for a moving object lattice sites are uncovered as the object passes by and require a reinitialization of the corresponding populations. In this case, we use the equilibrium distribution with the wall-velocity and a local density average. It needs to be mentioned that the wall-velocity may not be too large compared to the speed of sound $c_{s}$ . As shown in Dorschner et al. (Reference Dorschner, Chikatamarla, Bösch and Karlin2015), this yields a boundary condition compatible with the entropic considerations and thus the KBC model, where the stability issues of previously proposed schemes are resolved and second-order accuracy is retained. An extended discussion on the boundary conditions including all technical details is given in appendix B.

3.1 Turbulence in a periodic box

As the main flow in the bulk of the engine-like assembly is driven by vortex–vortex interaction and small-scale turbulence, we consider the Kida vortex flow as a classical benchmark simulation (see figure 1 a). This well-studied flow evolves from a deterministic and symmetric initial condition to a state that resembles a fully developed turbulent flow featuring a corresponding energy cascade and has been analysed extensively using DNS (Kida Reference Kida1985; Kida & Murakami Reference Kida and Murakami1987; Keating et al. Reference Keating, Vahala, Yepez, Soe and Vahala2007; Chikatamarla et al. Reference Chikatamarla, Frouzakis, Karlin, Tomboulides and Boulouchos2010). The initial conditions are given by

(3.1) $$\begin{eqnarray}\displaystyle \begin{array}{@{}rcl@{}}u_{x}(x,y,z)\; & =\; & U_{0}\sin x(\cos 3y\cos z-\cos y\cos 3z),\\ u_{y}(x,y,z)\; & =\; & U_{0}\sin y(\cos 3z\cos x-\cos z\cos 3x),\\ u_{z}(x,y,z)\; & =\; & U_{0}\sin z(\cos 3x\cos y-\cos x\cos 3y),\end{array} & & \displaystyle\end{eqnarray}$$

where $x,y,z\in [0,2{\rm\pi}]$ and periodic boundary conditions are imposed in all directions. The Reynolds number is defined as $Re=U_{0}N/{\it\nu}$ where $N$ is the domain size. While the kinetic energy is decaying, the evolution of enstrophy shows a steep increase in the early stage of the simulation and reaches a maximum value before it starts to decay. Just after the peak of enstrophy the flow reaches the most turbulent state producing large gradients and small-scale structures. While large gradients on the one hand may cause numerical instabilities, it is of paramount interest not to over-damp the dynamics on the other hand, which will lead to a corrupted and non-physical result lacking the small scales. Without employing explicit turbulence models one is usually restricted to increasing the resolution such that the smallest eddies are resolved. This is typically satisfied when the grid spacing is smaller than the Kolmogorov scale ${\it\eta}=({\it\nu}^{3}/{\it\epsilon})^{1/4}$ with kinematic viscosity ${\it\nu}$ and rate of energy dissipation ${\it\epsilon}$ .

In order to study the accuracy of the KBC model, a detailed investigation was recently conducted by Bösch et al. (Reference Bösch, Chikatamarla and Karlin2015). The Reynolds number here is $Re=6000$ , which is slightly higher than what can be expected for the valve/piston assembly considering the cylinder diameter and the maximum piston velocity as characteristic scales. A sufficiently resolved reference simulation ( ${\it\eta}\approx 1.2~{\rm\Delta}x$ ) is conducted with a box length of $N=600$ using the LBGK collision model. Further simulations with $N=\{100,200,400\}$ and the same Reynolds number using the KBC model are the carried out and compared to both the reference solution and theoretical limits (see Bösch et al. Reference Bösch, Chikatamarla and Karlin2015).

Figure 2.

Statistics for the Kida vortex flow at $Re=6000$ and resolutions of $N=\{100,200,400\}$ and $N=600$ for the simulations with KBC and the reference simulation with LBGK, respectively. The theoretical Kolmogorov scaling is indicated by the dotted line.

An important global characteristic is the evolution of the turbulent kinetic energy $k$ as shown in figure 2(a). For all resolutions in this study the energy decay seems to be captured well despite the rather severe under-resolution in the coarsest simulation ( ${\it\eta}\approx 0.2~{\rm\Delta}x$ ). However, a more meaningful insight is given by the energy distribution across the scales of the flow as shown by means of the normalized energy spectrum along with the theoretical Kolmogorov scaling with a slope of $-5/3$ in the inertial subrange in figure 2(b). It is apparent that the energy scales with marginal difference for all resolutions with a sharp cut-off at its smallest scale as expected for a well-behaved subgrid model.

A thorough convergence study of various statistical quantities sampled at time points around the peak of enstrophy demonstrates second-order accuracy (see Bösch et al. Reference Bösch, Chikatamarla and Karlin2015) as is expected for a lattice Boltzmann method.

While convergence towards the reference solution is established, it is of interest to quantify the recovery of the Navier–Stokes equations at small scales. To that end, let us remind the reader that the incompressible Navier–Stokes equation implies the following balance equations for the averaged momentum, vorticity, energy and enstrophy which yield for statistically homogeneous flows (Lamb Reference Lamb1932; Batchelor Reference Batchelor2000),

(3.2) $$\begin{eqnarray}\displaystyle & \displaystyle \partial _{t}\langle \boldsymbol{u}\rangle =0, & \displaystyle\end{eqnarray}$$

(3.3) $$\begin{eqnarray}\displaystyle & \displaystyle \partial _{t}\langle {\bf\omega}\rangle =0, & \displaystyle\end{eqnarray}$$

(3.4) $$\begin{eqnarray}\displaystyle & \displaystyle \partial _{t}k=-2{\it\nu}{\it\Omega}, & \displaystyle\end{eqnarray}$$

(3.5) $$\begin{eqnarray}\displaystyle & \displaystyle \partial _{t}{\it\Omega}=\langle {\bf\omega}\boldsymbol{\cdot }\unicode[STIX]{x1D668}\boldsymbol{\cdot }{\bf\omega}\rangle -2{\it\nu}P, & \displaystyle\end{eqnarray}$$

where

(3.6) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D668}={\textstyle \frac{1}{2}}(\boldsymbol{{\rm\nabla}}\boldsymbol{u}+\boldsymbol{{\rm\nabla}}\boldsymbol{u}^{\dagger }) & & \displaystyle\end{eqnarray}$$

is the rate-of-strain tensor and $P$ is the palinstrophy,

(3.7) $$\begin{eqnarray}\displaystyle P={\textstyle \frac{1}{2}}\langle \boldsymbol{{\rm\nabla}}{\bf\omega}\boldsymbol{ : }\boldsymbol{{\rm\nabla}}{\bf\omega}\rangle . & & \displaystyle\end{eqnarray}$$

Table 1.

Effective viscosity ratios at non-dimensional time $t/(N/U_{0})=0.75$ for simulations with different resolutions $N$ .

While the global conservation of average momentum (3.2) and vorticity (3.3) are satisfied up to machine precision for all times and all resolutions considered in Bösch et al. (Reference Bösch, Chikatamarla and Karlin2015), the balance of various terms in the energy (3.4) and enstrophy (3.5) equations is directly probing the recovery of the Navier–Stokes equation at small scales by the KBC model. To that end, we recast the balance equations (3.4) and  (3.5) in terms of the effective viscosity,

(3.8) $$\begin{eqnarray}\displaystyle & \displaystyle {\it\nu}_{eff,k}=-\frac{\partial _{t}k}{2{\it\Omega}}, & \displaystyle\end{eqnarray}$$

(3.9) $$\begin{eqnarray}\displaystyle & \displaystyle {\it\nu}_{eff,{\it\Omega}}=\frac{\langle {\bf\omega}\boldsymbol{\cdot }\unicode[STIX]{x1D668}\boldsymbol{\cdot }{\bf\omega}\rangle -\partial _{t}{\it\Omega}}{2P}. & \displaystyle\end{eqnarray}$$

In the simulation, the Navier–Stokes equation will be verified at small scales if the ratio ${\it\nu}_{eff}/{\it\nu}\approx 1$ . Thus, the evaluation of effective viscosities as in (3.8) and (3.9) is an important check of the accuracy and is listed in table 1. By increasing the resolution, the values approach ${\it\nu}_{eff,k}/{\it\nu}\approx 1$ . It is apparent that even for the coarsest run the additional dissipation is rather small, which is consistent with the evolution of turbulence kinetic energy $k$ shown in figure 2(a). The second effective viscosity ${\it\nu}_{eff,{\it\Omega}}$ is somewhat larger for simulations for the coarse grids, which is consistent with the under-prediction of the peak in enstrophy. For larger resolutions, however, the values are close to the nominal viscosity. Thus, we conclude that the KBC scheme recovers the Navier–Stokes equations well (in the absence of boundaries) while introducing only small additional dissipation on coarse grids. Further simulations of decaying turbulence are presented in Bösch et al. (Reference Bösch, Chikatamarla and Karlin2015). Note that due to the temporally and spatially varying parameter ${\it\gamma}$ , a fluctuating bulk viscosity is obtained, which in turn was found to reduce artificial compressibility effects in comparison to LBGK and other KBC models (see Bösch et al. Reference Bösch, Chikatamarla and Karlin2015).

In general, the family of KBC models has shown outstanding stability allowing for the operating range to be extended by orders of magnitude in terms of the Reynolds number compared to standard MRT- or LBGK-type models. Further, it has been shown to recover the well-established LBGK model for fully resolved simulations (Bösch et al. Reference Bösch, Chikatamarla and Karlin2015). In the next section, we aim to go beyond the periodic set-up in order to test for the next conceptual building block identified above.

3.2 Turbulent flow in a pipe

The chamber of the engine-like geometry is rotationally symmetric and it may be expected that the cylinder walls effect the dynamics of the flow to a large extent. Therefore, the turbulent flow through a round pipe is chosen as a validation of the second building block (see figure 1 b). This problem has been studied extensively in the literature experimentally, analytically and numerically. While for the flat channel there is consensus about the scaling of the mean velocity profile, it is less clear for the turbulent flow through a pipe and is being discussed in the literature (see, e.g., Barenblatt, Chorin & Prostokishin Reference Barenblatt, Chorin and Prostokishin1997; Zagarola, Perry & Smits Reference Zagarola, Perry and Smits1997; Zagarola & Smits Reference Zagarola and Smits1998; Wosnik, Castillo & George Reference Wosnik, Castillo and George2000; Perry, Hafez & Chong Reference Perry, Hafez and Chong2001; McKeon et al. Reference McKeon, Swanson, Zagarola, Donnelly and Smits2004; Morrison et al. Reference Morrison, McKeon, Jiang and Smits2004; Wei et al. Reference Wei, Fife, Klewicki and McMurtry2005; Wu & Moin Reference Wu and Moin2008). Nevertheless, there exist reliable DNS and experimental data. Here we choose a Reynolds number $Re_{D_{p}}=5300$ based on the pipe diameter $D_{p}$ and the mean bulk velocity $\overline{u_{bulk}}$ , the same as in the DNS of Wu & Moin (Reference Wu and Moin2008). This number is well within the range of what is expected in the chamber or the valve/piston assembly considered below.

As the problem is axially symmetric, it is conveniently formulated in cylindrical coordinates and one typically uses a corresponding computational mesh. The classical LB method, however, is restricted to a rectilinear Cartesian mesh (which is also employed for the engine-like geometry below). Thus, this benchmark problem is probing the performance of the boundary condition for curved walls to its full extent as the flow is wall bounded.

Three simulations are conducted at diameter $D_{p}=\{49,99,199\}$ lattice units (runs A, B and C). The domain length in the stream-wise direction is $L=16\,R$ , where $R=D_{p}/2$ is the pipe radius. The flow is initialized with a random velocity field and evolved for $200\,T$ , where the turnover time is given by $T=R/\overline{u_{bulk}}$ . After this initial transient, statistics are collected for another $200\,T$ , yielding a total run-time of $400\,T$ . The pressure gradient was adjusted during the simulation to reach the desired Reynolds number, which was realized through a body force. The corresponding Kàrmàn number is $R^{+}=u_{{\it\tau}}R/{\it\nu}=180$ with the wall friction velocity $u_{{\it\tau}}$ and the kinematic viscosity ${\it\nu}$ .

The distance from the pipe wall is given by $R-r$ with $r=\sqrt{x^{2}+y^{2}}$ , where $z$ denotes the spatial coordinate in the stream-wise direction. The non-dimensional wall units employed hereafter are defined as $x^{+}=xu_{{\it\tau}}/{\it\nu}$ and $u^{+}=u/u_{{\it\tau}}$ for space and velocity, respectively. Thus, the non-dimensional distance to the wall is $(R-r)^{+}=R^{+}-r^{+}$ . A natural measure for spatial resolution is the non-dimensional and uniform grid spacing ${\rm\Delta}x^{+}$ here, while for the DNS of Wu & Moin (Reference Wu and Moin2008) radial ( ${\rm\Delta}r^{+}$ ), azimuthal ( ${\rm\Delta}(r{\it\theta})^{+}$ ) and stream-wise ( ${\rm\Delta}z^{+}$ ) directions vary non-uniformly. The finest resolution for the DNS is typically found at the wall in the wall-normal direction. In Wu & Moin (Reference Wu and Moin2008) this amounts to ${\rm\Delta}r^{+}|_{r=R}=0.167$ (with a maximum ${\rm\Delta}r^{+}|_{r=0.409R}=1.647$ ), while in our simulations ${\rm\Delta}x^{+}=7.3$ (run A), ${\rm\Delta}x^{+}=3.6$ (run B) and ${\rm\Delta}x^{+}=1.8$ (run C).

Figure 3 shows a comparison of the mean stream-wise velocity component of runs A–C to the reference DNS results. Despite severe under-resolution, excellent agreement can be observed for runs B and C. Run A obviously employs a mesh which is too coarse to capture the scaling of the mean velocity correctly.

The next order statistical moments are shown in figure 4(a–d). The root-mean-square (r.m.s.) fluctuations of the axial, radial and azimuthal velocity component show the same trend as seen in figure 3; the coarsest simulation does not reproduce the expected values while runs B and C are very close to the DNS results. Moreover, the cross-correlations of the axial and radial fluctuations, figure 4(d), show excellent agreement for simulations B and C as well.

Figure 3.

Mean velocity component in flow direction for the turbulent pipe flow.

Figure 4.

Root mean square velocity profiles for the turbulent pipe flow. The legend is identical to figure 3.

3.3 Sedimenting sphere

The flow in the valve/piston assembly is induced through the moving piston and it is expected that the treatment of the moving boundary has a significant influence on the flow field inside the chamber. Stability and accuracy of the wall-boundary treatment as outlined in § 2.1 were investigated for various flow set-ups including moving objects in two dimensions in Dorschner et al. (Reference Dorschner, Chikatamarla, Bösch and Karlin2015). Here, we present a three-dimensional validation, where a sedimenting rigid sphere is considered. The sphere with density ${\it\rho}_{s}$ is resolved with a diameter of $D_{s}=30$ and the particle Reynolds number was taken to be $Re=100$ . When released from rest in a quiescent fluid with density ${\it\rho}_{f}$ , the sphere accelerates to its asymptotic settling velocity $u_{t}$ , where the gravitational force $F_{g}={\rm\pi}D_{s}^{3}{\it\rho}_{s}g/6$ is balanced by the buoyancy force $F_{b}={\rm\pi}D_{s}^{3}{\it\rho}_{f}g/6$ and the drag force $F_{d}=({\it\rho}_{f}u^{2}{\rm\pi}D_{s}^{2}C_{d})/8$ . The force balance immediately yields

(3.10) $$\begin{eqnarray}\displaystyle \frac{{\it\rho}_{s}}{{\it\rho}_{f}}=1+\frac{3u_{t}^{2}C_{d}}{4D_{s}g}, & & \displaystyle\end{eqnarray}$$

where $C_{d}$ denotes the drag coefficient. We set the nominal settling velocity to $u_{t,n}=0.01$ in lattice units and choose the particle density in accordance with (3.10). The drag coefficient was measured in a separate simulation for a stationary sphere and a mean flow velocity, which resulted in $C_{d}=1.1$ in agreement with literature values (see, e.g., Roos & Willmarth Reference Roos and Willmarth1971; Johnson & Patel Reference Johnson and Patel1999; Kim, Kim & Choi Reference Kim, Kim and Choi2001). In the dynamic case, we solve Newton’s equations for the particle motion

(3.11) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}\boldsymbol{x}_{s}}{\text{d}t}=\boldsymbol{u}_{s}, & \displaystyle\end{eqnarray}$$

(3.12) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}\boldsymbol{u}_{s}}{\text{d}t}=\frac{\boldsymbol{F}}{m_{s}}+\left(1-\frac{{\it\rho}_{f}}{{\it\rho}_{s}}\right)\boldsymbol{g}, & \displaystyle\end{eqnarray}$$

where the force $\boldsymbol{F}$ acting on the particle was evaluated using the momentum exchange method (see, e.g., Mei et al. Reference Mei, Yu, Shyy and Luo2002). We measure the terminal settling velocity to be $u_{t}=0.0100033$ , which amounts to less than $0.033\,\%$ error. Its equivalence to the nominal settling velocity establishes Galilean invariance between the static and the dynamic case and therefore validates the moving boundary condition.

4.1 Numerical set-up

The numerical set-up is identical to the experimental work of Morse et al. (Reference Morse, Whitelaw and Yianneskis1979) and the DNS of Schmitt et al. (Reference Schmitt, Frouzakis, Tomboulides, Wright and Boulouchos2014a ) for which the schematic and all geometrical specifications are shown in figure 5. The axis-symmetric assembly consists of a cylinder with diameter $D_{c}=75$  mm, a static centred valve and a flat piston for which a sinusoidal motion corresponding to a speed of $200$  rpm, a stroke of $S=60$  mm and a clearance at top dead centre (TDC) of $30$  mm is imposed. The maximum piston speed and the viscosity were chosen to obtain a Reynolds number of $Re=v_{p,max}D_{c}/{\it\nu}=3070$ . The valve has the radius $r_{v}=16.8$  mm and an angle of $30^{\circ }$ , which results in a uniform valve gap of approximately $4$  mm. In the experiment, the system is open and was mimicked in our simulation by a large reservoir on top of this assembly (not shown in the schematic) with a volume roughly three times as large as displaced by the cylinder. Due to the open system, compressibility effects are negligible and the flow may be regarded as incompressible without any phase change. Initially, a quiescent flow field is imposed where the piston is located at TDC. Further, the inlet condition at the top of the reservoir is evaluated instantaneously using the current piston velocity and the ratio of the reservoir and piston area. For all other boundary conditions, the no-slip boundary condition as outlined in § 2.1 is imposed.

Figure 5.

Schematic of the valve/piston assembly (all measures are given in mm).

In total, eight cycles were simulated (see the supplementary movie available athttp://dx.doi.org/10.1017/jfm.2016.448) and the first two were neglected in the accumulation of statistics to avoid initialization effects. Statistical quantities were computed using azimuthal and ensemble (across cycles) averaging. To study the nature of the subgrid model, simulations were carried out for a cylinder resolution of $D_{c,lb}=\{100,150,300\}$ points, respectively. This yields a total of approximately $N=\{8.5,28,228\}\text{Mio}$ nodes for each simulation and a cubic domain, respectively. The resolutions were chosen based on our preliminary studies in § 3 for which the flow regimes in terms of Reynolds number are similar to the one exhibited here.

Using a CRAY XC40 system the computational resources for $D_{c,lb}=300$ amount to $8.6\times 10^{3}\text{CPUh}$ per cycle, which is roughly $20$ times less than the corresponding requirement for the DNS in Schmitt et al. (Reference Schmitt, Frouzakis, Tomboulides, Wright and Boulouchos2014a ). The computational cost scales as $p^{D+1}$ for a refinement ratio $p$ and the dimensionality $D$ . Note that since a cubic domain and a regular grid without any local mesh refinement is used, significant optimization options remain. This is focus of future work as the scope of the study herein is to validate the implicit subgrid model of KBC for complex flows beyond classical benchmarks.

Figure 6.

Comparison of streamlines of the averaged velocity field. From (a) to (c): present results, DNS and experimental data are shown for different crank angles, respectively.

4.2 Velocity field

Resolved simulation and the DNS limit

In this section, we report the radial velocity profiles obtained for the intake stroke at crank angles (CA) of $36^{\circ }$ , $90^{\circ }$ and $144^{\circ }$ on different axial planes and make a comparison with the DNS and the experimental data. In the first instance, the profiles for the highest resolution are compared to establish the correct convergence of KBC models. Later in this section, the influence of under-resolution and its convergence is investigated. The main features of the mean flow field may be identified by the averaged streamlines shown in figure 6. In good agreement with both experimental and DNS data, the first phase of the intake stroke consists of fluid drawn into the chamber by the accelerated motion of the piston, forming a hollow cone jet. As a consequence, small vortex rings on both sides of the jet are generated, where the inner ring develops into the main feature as the piston moves further down and deflects the jet inwards before it is reflected by the piston. At this stage, one may observe jet breakup processes and the transition to turbulence caused by the interaction with the flow field. For instance, the Kelvin–Helmholtz instabilities in the shear layer are captured in the velocity snapshot of figure 7 for the first cycle. Subsequently, the jet is deflected towards the cylinder wall and a third vortex becomes apparent in the averaged flow field (see figure 6 c).

Note that at $90^{\circ }$ CA the flow undergoes a transition phase as the jet has reached its highest momentum and a fully turbulent flow field, dominated by the large vortical structure, is developed by interaction of the jet with the flow field.

Figure 7.

Instantaneous velocity magnitude snapshot at $77^{\circ }$ CA.

Quantitatively we compare the corresponding mean and r.m.s. velocity profiles of the axial velocity component $v_{z}$ in radial direction on different axial planes with respect to the cylinder head in figures 8–10). Consistent with the previous observation, the flow at $36^{\circ }$ CA is mainly dominated by the incoming jet, which results in a peak velocity at $r/R_{c}\approx 0.6$ . The r.m.s. velocity profiles indicate that the turbulence generation may be localized to jet shear layer on both sides. It is clear that the results predicted by the presented KBC model compare well to both DNS and experiment in terms of mean jet velocities and location. The r.m.s. velocities near the cylinder head at $z=10$ mm are also on top of the DNS data. It is apparent that more turbulence is generated compared to DNS and experiment close to the cylinder axis for the $z=20$  mm plane, resulting in a broadening of the r.m.s. velocity profile towards the centre. However, this does not seem to affect the good agreement for different crank angles as shown below.

At $90^{\circ }$ CA, the deflection of the jet at the cylinder wall leads to an increase of the jet radius, which is manifested by a shift of the average velocity profiles towards the cylinder wall. It should be noted that, as observed in the LES study of Liu & Haworth (Reference Liu and Haworth2010), the flow transition at this stage of highest momentum is difficult to capture for various turbulence modelling approaches. In their work, better agreement was achieved by adapting the model parameters but led to more discrepancies for $36^{\circ }$ CA and $144^{\circ }$ CA. Despite that, figure 9 clearly shows excellent agreement with the KBC model for both mean and r.m.s. velocity profiles compared with the reference data. A small shift in radial direction compared to the DNS is visible for the upper axial planes and the mean profiles but matches the experimental data. In planes beyond $z=20$  mm, our simulation is almost indistinguishable from the DNS result. Analogously, this holds for the r.m.s. velocity profiles.

Figure 8.

Comparison of the ensemble and azimuthally averaged axial mean (a) and r.m.s. (b) velocities at $36^{\circ }$ CA for a resolution of $D_{c,lb}=300$ .

Figure 9.

Comparison of the ensemble and azimuthally averaged axial mean (a) and r.m.s. (b) velocities at $90^{\circ }$ CA for a resolution of $D_{c,lb}=300$ .

Figure 10.

Comparison of the ensemble and azimuthally averaged axial mean (a) and r.m.s. (b) velocities at $144^{\circ }$ CA for a resolution of $D_{c,lb}=300$ .

At later crank angles, the piston is decelerating and we compare the corresponding profiles for $144^{\circ }$ CA in figure 10. Similarly, the discrepancies are marginal for all data.

Next, let us consider the turbulence statistics of the flow field by means of the Reynolds stress components and the turbulent kinetic energy $k$ as shown in figures 11–12. As the comparison is very similar during all phases, we choose the most interesting transition phase at $90^{\circ }$ CA as representative. It can clearly be seen that the turbulent kinetic energy follows the trend of the reference data. Its shift is analogous to the mean and r.m.s. velocities as shown above and minor discrepancies in terms of magnitude are visible. This, however, becomes more pronounced for lower values of the kinetic energy, suggesting that the contribution of the smallest scales is not fully accounted for. Similarly, the radial component $\overline{v_{r}^{\prime }v_{r}^{\prime }}$ shows a good agreement for the region near the jet entry, but dissipates slightly more. Analogous behaviour may be found for the Reynolds stress component $\overline{v_{{\it\varphi}}^{\prime }v_{{\it\varphi}}^{\prime }}$ and $\overline{v_{z}^{\prime }v_{r}^{\prime }}$ .

Under-resolved simulations

The above comparison demonstrated an overall excellent agreement between DNS, experimental data and our simulation using the entropic multi-relaxation time lattice Boltzmann model and the highest resolution of $D_{c,lb}=300$ . This provides evidence for the convergence of our simulations to DNS. Having established the correct limit, we can now study the effect of under-resolution.

Figure 11.

Comparison of the ensemble and azimuthally averaged turbulent kinetic energy $k$ and the Reynolds stress component $\overline{v_{r}^{\prime }v_{r}^{\prime }}$ at $90^{\circ }$ CA for a resolution of $D_{c,lb}=300$ .

Figure 12.

Comparison of the ensemble and azimuthally averaged Reynolds stress components $\overline{v_{{\it\varphi}}^{\prime }v_{{\it\varphi}}^{\prime }}$ and $\overline{v_{z}^{\prime }v_{r}^{\prime }}$ at $90^{\circ }$ CA for a resolution of $D_{c,lb}=300$ .

The discussion is similar for all crank angles and we again choose the most critical phase at $90^{\circ }$ CA as representative. In figure 13, the mean and r.m.s. velocity profiles are shown for cylinder resolutions of $D_{c,lb}=\{100,150,300\}$ points. It may be observed that the differences in terms of the average axial velocity profiles between $D_{c,lb}=150$ and $D_{c,lb}=300$ are marginal. The r.m.s. velocity profiles are very similar as well and only a small over-prediction at the lower resolution may be observed. At first sight, this is not intuitive as one would expect the turbulence production to be smoothed out, which would lead to lower values of the r.m.s. velocities. However, in this case, it may be traced back to very subtle inaccuracies of jet breakup processes, which affect the formation of the large-scale structures of the flow in a nonlinear manner and in turn affect the turbulence production leading to slightly higher r.m.s. values. On the other hand, for the case of $D_{c,lb}=100$ , the effect of under-resolution becomes noticeable. The incoming jet is hardly broken up by the surrounding flow field but rather by the confining cylinder walls, which leads to a higher penetration depth and jet radius on average. This may also be seen when considering the corresponding r.m.s. profiles, where the consistent under-prediction for $z=10{-}30$  mm leads to less turbulence production and therefore a delayed jet breakup with an increased diameter. Note that this is well within expectations as the valve gap is resolved by only three fluid points for this case. A comparison to our previous simulations of the pipe would suggest that the resolution is sufficient even at $D_{c,lb}=100$ . However, it needs to be kept in mind that due to the impinging jet the Kàrmàn number varies in the range of $R^{+}=u_{{\it\tau}}R/{\it\nu}\approx 90{-}780$ (measured for $D_{c,lb}=300$ ) and that for the local maximum of $R^{+}\approx 780$ the effect of the small-scale structures cannot be captured accurately at this resolution. On the other hand, the Kida vortex simulation reveals that the small structures in the bulk of the chamber are not correctly represented for a resolution as coarse as $D_{c,lb}=100$ , which further suggests that their contribution to the jet breakup is not fully accounted for. Hence, as expected, the complex flow in the cylinder as a combination of distinct flow features requires a higher resolution than the individual building blocks. A slight increase in resolution to $D_{c,lb}=150$ is necessary. Note that as indicated by the precursor simulations in §§ 3.1 and 3.2, this is still an under-resolved simulation in which the built-in subgrid model of KBC captures almost all flow features accurately.

Figure 13.

Comparison of the ensemble and azimuthally averaged axial mean (a) and r.m.s. (b) velocities at $90^{\circ }$ CA for different resolutions.

4.3 Cyclic variability

Capturing cyclic variability in internal combustion engine simulations is of crucial importance for its design as it may significantly influence its efficiency and pollutant formation among various other effects. Such variability has also been observed both numerically and experimentally in the valve/piston assembly (Haworth Reference Haworth1999; Rutland Reference Rutland2011; Schmitt et al. Reference Schmitt, Frouzakis, Tomboulides, Wright and Boulouchos2014a ). Its origin was discussed in the DNS work of Schmitt et al. (Reference Schmitt, Frouzakis, Wright, Tomboulides and Boulouchos2014b ), where it was mainly attributed to the remaining radial velocity at TDC from the previous cycle, the consequential jet location as well as the vortex ring orientation at bottom dead centre (BDC).

The aim of this section is to investigate the predictive capabilities of KBC models concerning the cyclic variability. An indicator of cyclic variability is given by the large vortex ring formed at BDC, which is visualized for all simulated cycles in figure 14 by means of the pressure iso-surfaces. Variation in shape, size and orientation is observed. Particularly pronounced is the distortion of the vortex ring for cycles five and six.

Quantitatively, we compare the correlation plot of the normalized, average radial velocity $\overline{v_{r,T}}/v_{p,mean}$ at TDC versus the normalized, average jet radius $\overline{r_{jet}}/r_{c}$ at $45^{\circ }$ CA in figure 15. Both quantities are evaluated as in Schmitt et al. (Reference Schmitt, Frouzakis, Wright, Tomboulides and Boulouchos2014b ), where the radial velocity is averaged at TDC azimuthally and radially in an axial window of $11.5$ mm starting from the cylinder head. The mean jet radius is defined as the average of the jet radii contained in the iso-surface of $u_{z}=-2.5~\text{m s}^{-1}$ . Our results show that the phenomenon of a cyclic varying jet radius is observed for all resolutions and confirm the result of Schmitt et al. (Reference Schmitt, Frouzakis, Wright, Tomboulides and Boulouchos2014b ) that it may be correlated to the residual turbulence at TDC from the previous cycle, where an increased averaged radial velocity at TDC results in an increased average jet radius and vice versa. It is notable that even the lowest resolution, for which the average axial velocity profiles have shown effects of under-resolution, exhibits pronounced cycle-to-cycle variation with the qualitatively correct cause and effect relation. This is indicated in figure 15, where the linear curve fit has an almost identical slope for all resolutions and therefore establishes the correlation of the normalized, average radial velocity $\overline{v_{r,T}}/v_{p,mean}$ at TDC versus the normalized, average jet radius $\overline{r_{jet}}/r_{c}$ . In addition, figure 15 reflects the dominant contributions to the subsequent jet breakup and is consistent with our previous observations in figure 13. For the lowest resolution, one can observe the most severe increase of the mean jet radius, which again may be attributed to the fact that the incoming jet is less influenced by the unresolved small-scale flow structures but rather is broken up by the interaction with the chamber walls leading to a broader jet with higher penetration depth. This was also observed for mean velocity profiles in figure 13. As the resolution is increased, the contributions of the small-scale structures become more pronounced and their interaction with the jet do not allow for such large jet radii as in the under-resolved case. On average this yields a lower penetration depth and smaller jet radii, consistent with both DNS and experimental data. As expected, the highest resolution yields the best match in terms of the average jet radius. More simulated cycles will presumably increase the spread on both the jet radius variability as well as the mean radial velocity.

Figure 14.

Cyclic variation of the vortex ring at BDC visualized by pressure iso-surfaces for $D_{c,lb}=300$ .

Figure 15.

Cyclic variation for all resolutions of $D_{c,lb}=\{100,150,300\}$ and the DNS results, quantified by mean jet radius $\overline{r_{jet}}/r_{c}$ at $45^{\circ }$ CA and correlated to the mean radial velocity $\overline{v_{r,T}}/v_{p,mean}$ at TDC. Lines represent a linear fit through the corresponding data points.