A Numerical Study of the Minimum-Dissipation Model for Large-Eddy Simulation in OpenFOAM

Abstract

The minimum-dissipation model (QR) has been utilized in studying turbulent channel flows at Reynolds numbers $\left({\mathit{Re}}_{\tau }\right)$ of up to 2000, as well as in investigating the flow past a circular cylinder at a Reynolds number (Re) of 3900, and flow over periodic hills at Re = 10,595. In our investigations, we have employed both symmetry-preserving discretizations and standard second-order accurate discretization methods within the OpenFOAM framework. The outcomes are compared with Direct Numerical Simulation (DNS) and experiment, indicating a favorable alignment between the QR results and the reference data. The findings suggest that the static QR model achieves comparable performance to dynamic models while cutting computational costs by a factor of three. The model coefficient C = 0.024 produces the most precise predictions, and as the mesh resolution increases, the influence of the subgrid model decreases, dropping to less than 20% of the molecular viscosity at the finest mesh. Furthermore, the QR model can predict the mean and root-mean-square velocity accurately up to ${\mathit{Re}}_{\tau }=2000$ without a wall damping function. The characteristics of turbulence strongly rely on spatial discretization methods. Various comparisons demonstrate the QR model conjugated with symmetry-preserving discretization performs better than the standard OpenFOAM discretization. Within the realm of OpenFOAM discretization, central difference schemes outperform other approaches.

1. Introduction

Turbulent flows are prevalent in numerous engineering applications, yet simulating them using Direct Numerical Simulation (DNS) can be prohibitively expensive, especially for high-Reynolds-number flows. The Unsteady Reynolds-Averaged Navier-Stokes (RANS) method, while reducing computational costs, introduces significant simplifications, leading to decreased accuracy due to its treatment of large-scale unsteadiness. To tackle these challenges, large-eddy simulation (LES) has been developed. LES directly resolves the larger-scale, unsteady turbulent motions while modeling the effects of smaller-scale motions. This is achieved by representing the unresolved scale of motion using a subgrid model.

Among the subgrid models available for LES, the Smagorinsky model stands out as the most widely used one [1]. While the Smagorinsky model demonstrates satisfactory performance in decaying homogeneous turbulence simulations [2], it tends to dissipate eddies inappropriately for laminar and transitional flows. The performance of the Smagorinsky model can be enhanced by dynamically computing the model coefficient, but this approach comes with increased computational costs. Another alternative is the WALE model [3], which addresses behavior near walls by utilizing the square of the velocity gradient tensor. Additionally, there is the Vreman model, which remains insensitive to pure shear but may lead to eddy dissipation in cases of back-scatter and solid body rotation [4].

Minimum-dissipation models offer a straightforward alternative to Smagorinsky-type approaches for parameterizing subfilter turbulent fluxes in large-eddy simulation (LES). The QR model, proposed by Verstappen [5, 6], is the first minimum-dissipation model introduced in this context. The QR model possesses several desirable characteristics: it is more computationally efficient than dynamic models, effectively deactivates in laminar and transitional flows, and maintains consistency with the exact subfilter stress tensor on isotropic grids. Subsequently, the anisotropic minimum-dissipation model (AMD) was developed by Rozema et al. specifically for flows on anisotropic grids [7]. Abkar and Moin utilized the AMD model to examine high-Reynolds-number rough-wall boundary layer flow [8]. Zahiri et al. further extended the application of the AMD model by implementing it into OpenFOAM and conducting tests on both single-phase and multi-phase flows [9].

However, there have been limited investigations into the QR model, particularly within open-source software. In this study, we introduce the QR model into OpenFOAM and conduct simulations in scenarios involving high Reynolds numbers and complex geometries, marking the first instance of such exploration.

The discretization error discretizing the partial differential equations into algebraic equations. Verstappen and Veldman [10] put forward a proposition to precisely maintain the symmetry properties of the underlying differential operators. This involves representing the convective operator with a skew-symmetric matrix and the diffusive operator with a symmetric, positive-definite matrix. Trias et al. [11] extended this concept to unstructured collocated grids and introduced an approach based on fully conservative regularization of the convective term to address the issue of checkerboard spurious modes. Expanding upon these concepts, Komen et al. [12] devised a symmetry-preserving second-order time-accurate PISO-based method for solving incompressible flow equations on collocated grids. Hopman et al. made the code available to the public [13].

This chapter starts by introducing the QR model, the symmetry-preserving discretization, and OpenFOAM standard discretization methods. Then, simulations are conducted for channel flow, flow over periodic hills, and flow over a cylinder. The model coefficient is optimized and an investigation of the numerical discretization schemes is carried out. The obtained results are compared with the references [14, 15, 16, 17, 18, 19, 20, 21]. The chapter concludes with a summary.

2. Numerical schemes

2.1 Minimum-dissipation model

The dynamics of large eddies in incompressible fluid flow are governed by the following momentum and continuity equations

$${\partial }_t\upsilon +\left(\upsilon \cdot \nabla \right)\upsilon +\nabla p-2\upsilon \nabla \cdot S\left(\upsilon \right)=-\nabla \cdot \tau \left(v\right),\nabla \cdot \upsilon =0,$$E1

In this context, v represents the spatially filtered velocity, $\nu$ denotes the kinematic molecular viscosity, and p represents the pressure. $S\left(\upsilon \right)=\left(\nabla \upsilon +\nabla {\upsilon }^T\right)/2$ refers to the symmetric portion of the velocity gradient. The subgrid tensor $\tau \left(\upsilon \right)$ can be expressed as:

$${\tau }_{\mathit{ij}}=\frac{1}{3}{\tau }_{\mathit{kk}}{\delta }_{\mathit{ij}}+\left({\tau }_{\mathit{ij}}-\frac{1}{3}{\tau }_{\mathit{kk}}{\delta }_{\mathit{ij}}\right)=\frac{2}{3}{k}_{\mathit{sgs}}{\delta }_{\mathit{ij}}+\left({\tau }_{\mathit{ij}}-\frac{1}{3}{\tau }_{\mathit{kk}}{\delta }_{\mathit{ij}}\right).$$

where $k_{\mathit{sgs}}=\frac{1}{2}{\tau }_{\mathit{kk}}$ is the subgrid kinetic energy. The subgrid stress tensor ${\tau }_{\mathit{ij}}$ is split into an isotropic part $\frac{1}{3}{\tau }_{\mathit{kk}}{\delta }_{\mathit{ij}}$ and anisotropic part ${\tau }_{\mathit{ij}}-\frac{1}{3}{\tau }_{\mathit{kk}}{\delta }_{\mathit{ij}}$. The isotropic part can be included into the pressure. The eddy-viscosity model describes the anisotropic part of the subgrid tensor as

$${\tau }_{\mathit{ij}}-\frac{1}{3}{\tau }_{\mathit{kk}}{\delta }_{\mathit{ij}}=-2{\nu }_{e}S\left(\upsilon \right)$$E2

Note that the trace of $S\left(\upsilon \right)$ is zero due to the divergence-free condition, $\nabla \cdot \upsilon =0$. The coefficient ${\nu }_e$ is termed the eddy viscosity of the model. This subgrid model is time irreversible (for ${\nu }_e>0$), meaning it contributes to dissipation forward in time. In the classical Smagorinsky eddy-viscosity model [1], the eddy viscosity is set equal to:

$${\nu }_e=C_S^2{\mathrm{\Delta }}^2\sqrt{4q}.$$E3

Here $q\left(\upsilon \right)=\frac{1}{2}\mathit{tr}\left(S^2\left(\upsilon \right)\right)$ represents the second invariant of the strain-rate tensor $S\left(\upsilon \right)$. It’s worth noting that $\mid S\left(\upsilon \right)\mid =\sqrt{2\mathit{tr}\left(S{\left(\upsilon \right)}^2\right)}=\sqrt{4q}$. Various Smagorinsky coefficient values $C_s$ have been suggested for different scenarios, spanning from $C_s=0.1‐0.17$ for channel flow and temporal mixing layers, to $C_s=0.20-0.22$ for decaying homogeneous isotropic turbulence [2].

The first minimum-dissipation model, introduced by Verstappen [5], is based on the invariants of the strain-rate tensor, and is designed to deactivate in laminar flow and flows exhibiting negative eddy dissipation. This model operates under the assumption that the eddy-viscosity model should prevent the residual field, denoted as ${\upsilon }^{\prime }=\upsilon -\overline{\upsilon }$ where the bar represents the LES filter, from being dynamically significant. This criterion is formalized by constraining the subgrid kinetic energy using Poincaé’s inequality. According to Poincaré’s inequality, there exists a constant $C_{\mathrm{\Delta }}$, contingent solely on ${\mathrm{\Omega }}_{\mathrm{\Delta }}$, such that for every function $\upsilon$ in the Sobolev space $W^{1,2}\left({\mathrm{\Omega }}_{\mathrm{\Delta }}\right)$

$${\int }_{{\mathrm{\Omega }}_{\mathrm{\Delta }}}\parallel \upsilon -\overline{\upsilon }{\parallel }^2\mathit{dx}\le C_{\mathrm{\Delta }}{\int }_{{\mathrm{\Omega }}_{\mathrm{\Delta }}}\parallel \nabla \upsilon {\parallel }^2\mathit{dx}$$E4

where the residual field ${\upsilon }^{\prime }=\upsilon -\overline{\upsilon }$ contains the eddies of size smaller than the length of the filter $\mathrm{\Delta }$, and $\parallel \cdot \parallel =\sqrt{⟨,⟩}$ is the standard norm of the inner product $⟨,⟩$ on the space of real-valued $L^2\left({\mathrm{\Omega }}_{\mathrm{\Delta }}\right)$ functions, where ${\mathrm{\Omega }}_{\mathrm{\Delta }}$ is a domain with diameter ${\mathrm{\Delta }}_j$ here it will be taken equal to a grid cell. The Poincaré constant $C_{\mathrm{\Delta }}$, which remains unaffected by $\upsilon$, is equal to the inverse of the smallest non-zero eigenvalue of the dissipative operator $-\mathrm{\Delta }=-\nabla \cdot \nabla ={\nabla }^T\nabla$ on a grid cell ${\mathrm{\Omega }}_{\mathrm{\Delta }}$ [22]. $-\nabla ={\nabla }^T$ is utilized for the $L^2\left({\mathrm{\Omega }}_{\mathrm{\Delta }}\right)$ inner product and periodic domain ${\mathrm{\Omega }}_{\mathrm{\Delta }}$. For convex domains, the Poincaré constant is denoted as $C_{\mathrm{\Delta }}={\left(\mathrm{\Delta }/\pi \right)}^2$ [23]. Poincaré’s inequality dictates that the kinetic energy of the residual field ${\upsilon }^{\prime }$ is bounded by a constant multiplied by the velocity gradient energy

$${\int }_{{\mathrm{\Omega }}_{\mathrm{\Delta }}}\frac{1}{2}\parallel {\upsilon }^{\prime }{\parallel }^2\mathit{dx}\le C_{\mathrm{\Delta }}{\int }_{{\mathrm{\Omega }}_{\mathrm{\Delta }}}\frac{1}{2}\parallel \nabla \upsilon {\parallel }^2\mathit{dx}$$E5

The equation governing the evolution of velocity gradient energy can be obtained by computing the inner product of $L^2$ with ${\nabla }^2\upsilon$. Integration by part yields

$$\frac{d}{\mathit{dt}}{\int }_{{\mathrm{\Omega }}_{\mathrm{\Delta }}}\frac{1}{2}\parallel \nabla \upsilon {\parallel }^2\mathit{dx}=-\nu {\int }_{{\mathrm{\Omega }}_{\mathrm{\Delta }}}\parallel {\nabla }^2\upsilon {\parallel }^2\mathit{dx}+{\int }_{{\mathrm{\Omega }}_{\mathrm{\Delta }}}\left(\upsilon \cdot \nabla \right)\upsilon \cdot \mathrm{\Delta }\mathit{\upsilon dx}-{\nu }_e{\int }_{{\mathrm{\Omega }}_{\mathrm{\Delta }}}\parallel {\nabla }^2\upsilon {\parallel }^2\mathit{dx}.$$E6

The boundary terms arising from integration by parts vanish since ${\mathrm{\Omega }}_{\mathrm{\Delta }}$ is a periodic box. The second term on the right-hand side of Eq. (6) signifies the generation of velocity gradient energy by the convective term. It can be represented as $r\left(\upsilon \right)=-\mathit{tr}\left(S^3\left(\upsilon \right)\right)/3=-\mathit{det}\left(S\left(\upsilon \right)\right)$, which corresponds to the third invariant of the strain-rate tensor $S\left(\upsilon \right)$ (for further details, see Ref. [5, 6, 7]). We make the assumption that the molecular and eddy viscosity remains constant over a grid cell. The third term in Eq. (6) accounts for the dissipation caused by eddy viscosity and can be expressed as $q\left(\omega \right)=\mathit{tr}\left(S^2\left(\omega \right)\right)/2$, where $\omega$ denotes the vorticity, $\omega =\nabla \times \upsilon$.

Incorporating $r\left(\upsilon \right)$ and $q\left(\omega \right)$ into the evolution equation of velocity gradient energy (Eq. (6))

$$\frac{d}{\mathit{dt}}{\int }_{{\mathrm{\Omega }}_{\mathrm{\Delta }}}\frac{1}{2}\parallel \nabla \upsilon {\parallel }^2\mathit{dx}=-\nu {\int }_{{\mathrm{\Omega }}_{\mathrm{\Delta }}}\parallel {\nabla }^2\upsilon {\parallel }^2\mathit{dx}+4{\int }_{{\mathrm{\Omega }}_{\mathrm{\Delta }}}r\left(\upsilon \right)\mathit{dx}-4{\nu }_e{\int }_{{\mathrm{\Omega }}_{\mathrm{\Delta }}}q\left(\omega \right)\mathit{dx}.$$E7

Let us consider that the eddy viscosity is chosen in such a way that the last two terms on the right-hand side of Eq. (7) counterbalance each other

$${\int }_{{\mathrm{\Omega }}_{\mathrm{\Delta }}}r\left(\upsilon \right)\mathit{dx}={\nu }_e{\int }_{{\mathrm{\Omega }}_{\mathrm{\Delta }}}q\left(\omega \right)\mathit{dx}$$E8

Then, we obtain

$$\frac{d}{\mathit{dt}}{\int }_{{\mathrm{\Omega }}_{\mathrm{\Delta }}}\frac{1}{2}\parallel \nabla \upsilon {\parallel }^2\mathit{dx}=-\nu {\int }_{{\mathrm{\Omega }}_{\mathrm{\Delta }}}\parallel {\nabla }^2\upsilon {\parallel }^2\mathit{dx}.$$E9

Utilizing Poincaré’s inequality and Gronwall’s lemma on the right-hand side of the aforementioned equation yields:

$${\int }_{{\mathrm{\Omega }}_{\mathrm{\Delta }}}\parallel {\upsilon }^{\prime }{\parallel }^2\left(x,t\right)\mathit{dx}\le C_{\mathrm{\Delta }}{\int }_{{\mathrm{\Omega }}_{\mathrm{\Delta }}}\parallel \nabla \upsilon {\parallel }^2\left(x,t\right)\mathit{dx}\le C_{\mathrm{\Delta }}{e}^{\frac{-2\mathit{\nu t}}{C_{\mathrm{\Delta }}}}{\int }_{{\mathrm{\Omega }}_{\mathrm{\Delta }}}\parallel \nabla {\upsilon }^{\prime }{\parallel }^2\left(x,0\right)\mathit{dx}.$$

The energy of the subgrid-scale decays at least as fast as $C_{\mathrm{\Delta }}{e}^{\frac{-2\mathit{\nu t}}{C_{\mathrm{\Delta }}}}$, for any filter length $\mathrm{\Delta }$. Therefore, we are able to keep the subfilter component ${\upsilon }^{\prime }$ under control through the utilization of Eq. (8). The minimum amount of eddy dissipation must adhere to the dissipation condition (9).

The expression for $q\left(\omega \right)$ in Eq. (8) can be equivalently written in terms of $q\left(v\right)$:

$${\int }_{{\mathrm{\Omega }}_{\mathrm{\Delta }}}q\left(\omega \right)\mathit{dx}=\frac{1}{4}{\int }_{{\mathrm{\Omega }}_{\mathrm{\Delta }}}{\left|\nabla \omega \right|}^2\mathit{dx}=\frac{1}{4}{\int }_{{\mathrm{\Omega }}_{\mathrm{\Delta }}}\omega \cdot \left(-\mathrm{\Delta }\right)\mathit{\omega dx}=\frac{{\int }_{{\mathrm{\Omega }}_{\mathrm{\Delta }}}\omega \cdot \left(-\mathrm{\Delta }\right)\mathit{\omega dx}}{{\int }_{{\mathrm{\Omega }}_{\mathrm{\Delta }}}\omega \cdot \mathit{\omega dx}}\cdot {\int }_{{\mathrm{\Omega }}_{\mathrm{\Delta }}}q\left(v\right)\mathit{dx}.$$E10

Thus, Eq. (8) becomes

$${\nu }_e=\frac{{\int }_{{\mathrm{\Omega }}_{\mathrm{\Delta }}}\omega \cdot \mathit{\omega dx}}{{\int }_{{\mathrm{\Omega }}_{\mathrm{\Delta }}}\omega \cdot -\mathrm{\Delta }\mathit{\omega dx}}\cdot \frac{{\int }_{{\mathrm{\Omega }}_{\mathrm{\Delta }}}r\left(\upsilon \right)\mathit{dx}}{{\int }_{{\mathrm{\Omega }}_{\mathrm{\Delta }}}q\left(\upsilon \right)\mathit{dx}}$$E11

The first fraction on the right-hand side above is at most $C_{\mathrm{\Delta }}$, which corresponds to the reciprocal of the smallest eigenvalue of $-\mathrm{\Delta }$ on ${\mathrm{\Omega }}_{\mathrm{\Delta }}$.

Therefore, we establish ${\nu }_e{\int }_{{\mathrm{\Omega }}_{\mathrm{\Delta }}}q\left(\upsilon \right)\mathit{dx}=C_{\mathrm{\Delta }}{\int }_{{\mathrm{\Omega }}_{\mathrm{\Delta }}}r\left(\upsilon \right)\mathit{dx}$. This equation guarantees that the subgrid scales remain dynamically insignificant, implying that their energy is bounded by Eq. (9), where the subgrid scales energy ${\int }_{{\mathrm{\Omega }}_{\mathrm{\Delta }}}\parallel {\upsilon }^{\prime }{\parallel }^2\left(x,t\right)\mathit{dx}$ decays at least as rapidly as $C_{\mathrm{\Delta }}{e}^{\frac{-2\mathit{\nu t}}{C_{\mathrm{\Delta }}}}$, for any filter length $\mathrm{\Delta }$. Consequently, the minimum eddy viscosity required to counteract the nonlinear production is determined by

$${\nu }_e=C_{\mathrm{\Delta }}\frac{\overline{\mid r\left(\upsilon \right)\mid }}{\overline{q\left(\upsilon \right)}}.$$E12

To ensure non-negativity of the eddy viscosity, we employ the absolute value of $r\left(\upsilon \right)$. The quantities $\overline{q\left(\upsilon \right)}$ and $\overline{r\left(\upsilon \right)}$ represent the grid cell averages of the second and third invariant of the strain-rate tensor, respectively. In practical applications, grid cell averaging is typically approximated using midpoint integration. This characterizes the QR model

$$\tau -\frac{1}{3}\mathit{tr}\left(\tau \right)I=-2{\nu }_{e}S\left(\upsilon \right)=-2C_{\mathrm{\Delta }}\frac{\mid r\left(\upsilon \right)\mid }{q\left(\upsilon \right)}S\left(\upsilon \right)$$E13

2.2 Numerical schemes in OpenFOAM

2.2.1 Spatial and temporal discretization

The discretization schemes typically maintain a second-order accuracy level. OpenFOAM employs a collocated arrangement where all variable values are computed and stored at the center $x_P$ of the control volume $V_P$. These values are depicted by a piecewise constant profile, commonly known as the mean value

$$\varphi \left(x_P\right)={\varphi }_P\approx \frac{1}{V_P}{\int }_{V_P}\varphi \left(x\right)\mathit{dV},$$E14

where $\varphi$ denotes a discretized quantity. By employing the Divergence theorem (or Gauss theorem), the volume integrals present in the governing equations are transformed into surface integrals. Consequently, the problem shifts toward interpolating cell-centered values to face-centered values. The face values involved in the convective and diffusive fluxes must be determined through some form of interpolation from the centroid of the control volumes to their respective faces. The Gauss linear interpolation scheme in OpenFOAM is given by ${\varphi }_f={\mathit{\alpha \varphi }}_P+\left(1-\alpha \right){\varphi }_N$, where $\alpha =\frac{f-N}{P-N}$, which yields a central difference scheme on a uniform mesh. The solver we use is pimpleFOAM if no other specification is provided.

If no specific instruction is provided, the continuous temporal derivative $\partial /\partial t$ is discretized using an implicit backward scheme, as indicated by Eq. (15)

$$\frac{3{\varphi }^{n+1}-4{\varphi }^n+{\varphi }^{n-1}}{2\mathrm{\Delta }t}\approx \frac{\partial }{\partial t}\left(\varphi \right),$$E15

where the n + 1 represents the value at the subsequent time level $t+\mathrm{\Delta }t$, n denotes the value at the current time level, and n – 1 signifies the value at the preceding time level $t-\mathrm{\Delta }t$.

2.3 Symmetry-preserving discretization scheme

2.3.1 Navier-stokes equation

The governing equation for incompressible Navier-Stokes (NS) is expressed as

$${\partial }_t\mathbf{u}+\left(\mathbf{u}\cdot \nabla \right)\mathbf{u}-\frac{1}{\mathit{Re}}\nabla \cdot \nabla \mathbf{u}+\nabla p=0,\nabla \cdot \mathbf{u}=0.$$E16

Here, u stands for the velocity (no filter applied), p is the pressure, and Re represents the Reynolds number.

When there are no external sources, for instance, body or boundary forces, the change rate of the total energy remains unaffected by both pressure disparities and convective transport. Instead, it is solely dictated by dissipation. This fundamental physical characteristic can be easily inferred from the symmetric properties of the differential operators in the NS Eq. (16).

The total energy of the flow (u, u) is defined using the conventional scalar product. Its temporal evolution can be derived by differentiating (u, u) with respect to time and then reformulating ${\partial }_t\mathbf{u}$ with the assistance of Eq. (16). This process yields:

$$\frac{d}{\mathit{dt}}\left(\mathbf{u},\mathbf{u}\right)=-\left(\left(\mathbf{u}\cdot \nabla \right)\mathbf{u},\mathbf{u}\right)-\left(\mathbf{u},\left(\mathbf{u}\cdot \nabla \right)\mathbf{u}\right)+\frac{1}{\mathit{Re}}\left(\left(\nabla \cdot \nabla \mathbf{u},\mathbf{u}\right)+\left(\mathbf{u},\nabla \cdot \nabla \mathbf{u}\right)\right)-\left(\nabla p,\mathbf{u}\right)-\left(\mathbf{u},\nabla p\right).$$

By performing integration by parts on the linear and trilinear forms on the right-hand side, disregarding any boundary contributions, we

$${\left(\mathbf{u},\nabla \right)}^{\ast }=-\left(\mathbf{u},\nabla \right)\text{and}{\nabla }^{\ast }=-\nabla$$E17

As a result of these (skew-) symmetries, the convective- and pressure-dependent terms cancel, leading to a simplification in the change rate of the total energy to

$$\frac{d}{\mathit{dt}}\left(\mathbf{u},\mathbf{u}\right)=-\frac{2}{\mathit{Re}}\left(\nabla \mathbf{u},\nabla \mathbf{u}\right)\le 0$$E18

In the discrete context, the evolution of energy also adheres to Eq. (18), where u is substituted with the discrete velocity, and $\nabla$ is replaced by its discrete counterpart. This holds true provided that the discretization of the differential operator upholds the (skew-) symmetries outlined in Eq. (17). Under this condition, the energy of any discrete solution remains conserved in the absence of viscosity, and diminishes over time in the presence of dissipation. Essentially, a symmetry-preserving spatial discretization of the NS equation exhibits unconditional stability and conservativeness.

2.3.2 Symmetry-preserving discretization

Consider the discretization of a first-order derivative in one spatial dimension. The Lagrangian interpolation (minimizing local truncation error) violates the skew-symmetry of the convective operator on nonuniform grids. Consequently, quantities conserved in the continuous formulation, such as kinetic energy, are not preserved at the discrete level. The energy undergoes systematic damping; for instance, in upwind methods, where the convective term introduces artificial dissipation that dampens the kinetic energy. Alternatively, explicit damping may be required to maintain stability. However, artificial dissipation disrupts the delicate balance between convective transport and physical dissipation, particularly at the smallest scales of motion, thereby distorting the essence of turbulence. Hence, symmetry-preserving discretization is employed.

Consider the continuing equation and linear momentum in one spatial dimension

$${\partial }_{t}u+\overline{u}{\partial }_{x}u-{\partial }_{\mathit{xx}}u/\mathit{Re}+{\partial }_{x}p=0,{\partial }_{x}u=0.$$E19

For simplicity, the convective transport velocity $\overline{u}$ is assumed to be constant. In matrix-vector notation, the spatial discretization of Eq. (19) can be expressed as:

$${\mathrm{\Omega }}_0\frac{d{\mathbf{u}}_h}{\mathit{dt}}+{\mathbf{C}}_{\mathbf{0}}\left(\overline{u}\right){\mathbf{u}}_h+{\mathbf{D}}_0{\mathbf{u}}_h+{\mathrm{\Omega }}_0{\mathbf{G}}_0{\mathbf{p}}_h=0,{\mathbf{M}}_0{\mathbf{u}}_h=0$$E20

where the diagonal matrix ${\mathrm{\Omega }}_0$ is built of the spacing of mesh: ${\left({\mathrm{\Omega }}_0\right)}_{i,i}=\frac{1}{2}\left(x_{i+1}-x_{i-1}\right)$, the discrete velocities $u_i\left(t\right)\approx u\left(x_i,t\right)$ constitute the vector ${\mathbf{u}}_{\mathbf{h}}$; the tridiagonal matrices ${\mathbf{C}}_{\mathbf{0}}\left(\overline{u}\right)$ and ${\mathbf{D}}_{\mathbf{0}}$ represent the convective and diffusive operate, respectively.

The continuity equation ${\partial }_{x}u=0$ is discretized as $\frac{u_{i+1}-u_{i-1}}{x_{i+1}-x_{i-1}}=0$. This yields ${\mathbf{M}}_0{\mathbf{u}}_h=0$.

The discretization of the pressure gradient is facilitated by utilizing the symmetry relation outlined in Eq. (17). According to Eq. (17), the continuous gradient operator equals the negative transpose of the divergence, implying that any velocity and pressure fields adhere to the relation $\left(\nabla p,\mathbf{u}\right)=-\left(p,\nabla \cdot \mathbf{u}\right)$. This relation also applies to the discrete pressure ${\mathbf{p}}_h$ and the discrete pressure gradient ${\mathbf{G}}_0{\mathbf{p}}_h$, expressed as

$${\left({\mathbf{G}}_0{\mathbf{p}}_h\right)}^{\ast }{\mathrm{\Omega }}_0{\mathbf{u}}_h={\mathbf{p}}_h^{\ast }{\mathbf{G}}_0^{\ast }{\mathrm{\Omega }}_0{\mathbf{u}}_h=-{\mathbf{p}}_h^{\ast }{\mathbf{M}}_0{\mathbf{u}}_h,$$E21

if the gradient operator is approximated by

$${\mathbf{G}}_0=-{\mathrm{\Omega }}_0^{-1}{\mathbf{M}}_0^{\ast }.$$E22

Note that the gradient matrix, describing the integration of the pressure over the control volumes, is represented by $-{\mathbf{M}}_0^{\ast }$. Since the discrete pressure gradient also inherits boundary conditions from the discrete divergence, there is no need to explicitly specify boundary conditions for the pressure.

In the absence of diffusion, that is for ${\mathbf{D}}_{\mathbf{0}}=0$, the energy ${‖{\mathbf{u}}_{\mathbf{h}}‖}^2={\mathbf{u}}_{\mathbf{h}}^{\ast }{\mathrm{\Omega }}_0{\mathbf{u}}_{\mathbf{h}}$ of any solution ${\mathbf{u}}_{\mathbf{h}}$ of the dynamic system of Eq. (20) is conserved if and only if the right-hand side of

$$\frac{d}{\mathit{dt}}{‖{\mathbf{u}}_{\mathbf{h}}‖}^2=-{\mathbf{u}}_{\mathbf{h}}^{\ast }\left({\mathbf{C}}_{\mathbf{0}}\left(\overline{u}\right)+{\mathbf{C}}_{\mathbf{0}}^{\ast }\left(\overline{u}\right)\right){\mathbf{u}}_{\mathbf{h}}+{\mathbf{u}}_h^{\ast }\left({\mathbf{M}}_0^{\ast }{\mathbf{p}}_h\right)+{\left({\mathbf{M}}_0^{\ast }{\mathbf{p}}_h\right)}^{\ast }{\mathbf{u}}_h$$

is zero. This property remains valid (for any ${\mathbf{u}}_h$) if and only if the coefficient matrix ${\mathbf{C}}_0\left(\overline{u}\right)$ exhibits skew-symmetry,

$${\mathbf{C}}_0\left(\overline{u}\right)+{\mathbf{C}}_0^{\ast }\left(\overline{u}\right)=0,$$E23

In other words, the discrete operator must retain the skew-symmetry observed in the continuous convective derivative $\left(\mathbf{u}\cdot \nabla \right)$.

To meet the skew-symmetry condition stated in Eq. (23), the interpolation weights of the neighboring discrete variables should be set to $1/2$. Consequently, symmetry-preserving discretization yields

$$\overline{u}{\partial }_{x}u\left(x_i\right)\approx \overline{u}\frac{u_{i+1}-u_{i-1}}{x_{i+1}-x_{i-1}}={\left({\mathrm{\Omega }}_0^{-1}{\mathbf{C}}_{\mathbf{0}}\left(\overline{u}\right){\mathbf{u}}_{\mathbf{h}}\right)}_i$$E24

The elements of the tridiagonal matrix ${\mathbf{C}}_{\mathbf{0}}\left(\overline{u}\right)$ are defined as follows: ${\mathbf{C}}_{\mathbf{0}}{\left(\overline{u}\right)}_{i,i-1}=-\frac{1}{2}\overline{u}$, ${\mathbf{C}}_{\mathbf{0}}{\left(\overline{u}\right)}_{i,i}=0$_, and ${\mathbf{C}}_{\mathbf{0}}{\left(\overline{u}\right)}_{i,i+1}=\frac{1}{2}\overline{u}$. It has been rigorously proven that Eq. (24) yields second-order accurate solutions on both uniform and nonuniform grids [24].

Diffusion is discretized in the same vein. The resulting coefficient matrix ${\mathbf{D}}_0$ is positive-definite, like the underlying differential operator $-{\partial }_{\mathit{xx}}$

$${\mathbf{D}}_0=\frac{1}{\mathit{Re}}{\mathrm{\Delta }}_0^{\ast }{\mathrm{\Lambda }}_0^{-1}{\mathrm{\Delta }}_0,$$

where the difference matrix ${\mathrm{\Delta }}_0$ is defined as ${\left({\mathrm{\Delta }}_0{\mathbf{u}}_h\right)}_i=u_i-u_{i-1}$, and the non-zero entries of the diagonal matrix ${\mathrm{\Lambda }}_0$ are given by ${\left({\mathrm{\Lambda }}_0\right)}_{i,i}=\delta x_i$. Consequently, the symmetric component of ${\mathbf{C}}_{\mathbf{0}}\left(\overline{u}\right)+{\mathbf{D}}_{\mathbf{0}}$ is solely influenced by diffusion, and hence is positive-definite. The energy of any solution ${\overline{\mathbf{u}}}_h$ of the semi-discrete system Eq. (20) evolves in a manner analogous to the continuous case; compare Eq. (18) to

$$\frac{d}{\mathit{dt}}\left({\mathbf{u}}_h^{\ast }{\mathrm{\Omega }}_0{\mathbf{u}}_h\right)\stackrel{\left(20\right)+\left(23\right)}{=}-{\mathbf{u}}_h^{\ast }\left({\mathbf{D}}_0+{\mathbf{D}}_0^{\ast }\right){\mathbf{u}}_h\le 0.$$

The right-hand side equals zero if and only if ${\mathbf{u}}_h$ lies within the null space of ${\mathbf{D}}_0+{\mathbf{D}}_0^{\ast }$. Thus, in summary, as the energy ${‖{\mathbf{u}}_h‖}^2=\left({\mathbf{u}}_h^{\ast }{\mathrm{\Omega }}_0{\mathbf{u}}_h\right)$ remains constant over time, a stable solution can be achieved on any grid. Similarly, higher-order symmetry-preserving discretization can be derived in a comparable manner (refer to [10] for details). Taking all ingredients together yields the symmetry-preserving discretization of the NS equation.

2.3.3 Symmetry-preserving discretization of Navier-Stokes equation

the incompressible Navier-Stokes equations in their semi-discrete form are expressed as

$$\mathrm{\Omega }\frac{d{\mathbf{u}}_h}{\mathit{dt}}+\mathbf{C}\left({\mathbf{u}}_h\right){\mathbf{u}}_h+{\mathbf{Du}}_h-{\mathbf{M}}^{\ast }{\mathbf{p}}_h=0,{\mathbf{Mu}}_h=0$$E25

Global conservation laws involve integrals across the flow domain, which are transformed into scalar products upon discretization of the flow. For example, the change of the total mass of the flow is discretized as the scalar product of the constant vector 1 (with dimensions equal to the number of grid cells) and the discrete mass flux ${\mathbf{Mu}}_h$. Due to this scalar product equals zero (${\mathbf{Mu}}_h=0$) the total mass is conserved.

To determine the total momentum, the scalar product of velocity vector ${\mathbf{u}}_{\mathbf{h}}$ with the vector $\mathrm{\Omega }\mathbf{1}$ is taken, where the constant vector now encompasses as many entries as there are control volumes for the discrete velocity components $u_{i,j}$ and $v_{i,j}$). The evolution of total momentum straightforwardly follows from Eq. (25)

$$\frac{d}{\mathit{dt}}\left({\mathbf{1}}^{\ast }\mathrm{\Omega }{\mathbf{u}}_h\right)=-{\mathbf{1}}^{\ast }\left(\mathbf{C}\left({\mathbf{u}}_h\right)+\mathbf{D}\right){\mathbf{u}}_h+{\mathbf{1}}^{\ast }{\mathbf{M}}^{\ast }{\mathbf{p}}_h=0.$$E26

Thus, momentum conservation is ensured if ${\left(\mathbf{C}\left({\mathbf{u}}_h\right)+\mathbf{D}\right)}^{\ast }1=0$, and the law of mass conservation is consistently discretized as $\mathbf{M}\mathbf{1}\mathbf{=}\mathbf{0}$. The former condition can be divided into two conditions: one for convective discretization, ${\mathbf{C}}^{\ast }\left({\mathbf{u}}_h\right)\mathbf{1}=0$, and one for diffusive discretization, ${\mathbf{D}}^{\ast }\mathbf{1}\mathbf{=}\mathbf{0}$. Moreover, we can omit the ${}^{\ast }$'s, yielding $\mathbf{C}\left({\mathbf{u}}_h\right)\mathbf{1}=0$ and $\mathbf{D}\mathbf{1}=\mathbf{0}$, since the convective matrix $\mathbf{C}\left({\mathbf{u}}_h\right)$ is skew-symmetry and the diffusive matrix D is symmetric. Thus, it suffices to ensure that the constant vector lies in the null space of the approximate, convective, and diffusive operators.

The discretization is configured to ensure that the evolution of the (kinetic) energy ${\mathbf{u}}_h^{\ast }\mathrm{\Omega }{\mathbf{u}}_h$ of any solution to Eq. (25) adheres to

$$\frac{d}{\mathit{dt}}\left({\mathbf{u}}_h^{\ast }\mathrm{\Omega }{\mathbf{u}}_h\right)=-{\mathbf{u}}_h^{\ast }\left(\mathbf{C}\left({\mathbf{u}}_h\right)\right)+{\mathbf{C}}^{\ast }\left({\mathbf{u}}_h\right)){\mathbf{u}}_h-{\mathbf{u}}_h^{\ast }\left(\mathbf{D}+{\mathbf{D}}^{\ast }\right){\mathbf{u}}_h+{\mathbf{u}}_h^{\ast }{\left({\mathbf{M}}^{\ast }\mathbf{p}\right)}_h+{\left({\mathbf{M}}^{\ast }{\mathbf{p}}_h\right)}^{\ast }{\mathbf{u}}_h,$$

where the right-hand side is negative for all ${\mathbf{u}}_h$'s, except those within the null space of $\mathbf{D}+{\mathbf{D}}^{\ast }$. The convective component cancels out due to $\mathbf{C}{\left(\mathbf{u}\right)}_h$‘s skew-symmetry, while the pressure terms nullify on the staggered grids (thus, cannot destabilize the spatial discretization) because the discrete pressure gradient correlates with the transpose of M, as seen in Eq. 22.

In summary, energy is conserved for inviscid flow, while for viscous flow, the energy ${‖{\mathbf{u}}_{\mathbf{h}}‖}^2={\mathbf{u}}_{\mathbf{h}}^{\ast }\mathrm{\Omega }{\mathbf{u}}_{\mathbf{h}}$ does not escalate in time. This indicates that the symmetry-preserving discretization Eq. 25 is stable and upholds conservation laws for mass, momentum, and energy.

2.3.4 Addressing the pressure-velocity coupling on the collocated grid

n a collocated grid setup the actual velocity ${\mathbf{u}}_{\mathbf{c}}$ is stored in the cell center. The coupling between velocity and pressure introduces an additional error term proportional to the third-order derivative of pressure ${\tilde{\mathbf{p}}}_c^{\prime }$ into the momentum equation. This issue is commonly referred to as the checkerboard problem [25]. Trias et al. [11] proposed a method to eliminate the checkerboard spurious mode without introducing any nonphysical dissipation. The concept behind this approach involves utilizing a linear shift operator to transform a cell-centered velocity ${\mathbf{u}}_{\mathbf{c}}$ into a staggered one ${\mathbf{u}}_{\mathbf{s}}$ and employing a fully conservative regularization of the convective term to restrain the generation of small motions.

The linear shift operator serves to establish a relationship between the staggered velocity field and the cell-centered counterparts and vice versa. In this context, the subscript c denotes variables centered on the collated mesh, while s indicates variables staggered on the faces. The cell-to-face linear shift operator, denoted as ${\mathrm{\Gamma }}_{c\to s}\in {\mathbb{R}}^{m\times 3n}$, converts a cell-centered velocity field into a staggered one

$${\mathbf{u}}_s={\mathrm{\Gamma }}_{c\to s}{\mathbf{u}}_c.$$E27

Conversely, the cell-centered fields are linked to the staggered ones through the linear shift transformation ${\mathrm{\Gamma }}_{s\to c}\in {\mathbb{R}}^{3n\times m}$,

$${\mathbf{u}}_c={\mathrm{\Gamma }}_{s\to c}{\mathbf{u}}_s.$$E28

Note the general relation ${\mathrm{\Gamma }}_{s\to c}{\mathrm{\Gamma }}_{c\to s}=I$ only holds approximately, meaning ${\mathbf{u}}_c\approx {\mathrm{\Gamma }}_{s\to c}{\mathrm{\Gamma }}_{c\to s}{\mathbf{u}}_c$.

The face-to-cell shift operator ${\mathrm{\Gamma }}_{s\to c}$ is constrained by Eq. (22) to the contribution of the pressure gradient term to the global kinetic energy. The expression for is as follows

$${\mathrm{\Gamma }}_{s\to c}={\left(I_3\otimes {\mathrm{\Omega }}_c\right)}^{-1}{\mathrm{\Gamma }}_{c\to s}^{\ast }{\mathrm{\Omega }}_s,$$E29

where $I_3\in {\mathbb{R}}^{3\times 3}$ is the identity matrix.

The linear shift operator is given by

$${\mathrm{\Gamma }}_{c\to s}=N_s\mathrm{\Pi },$$E30

where matrices $\mathrm{\Pi }\in {\mathbb{R}}^{3m\times 3n}$ and $N_s\in {\mathbb{R}}^{m\times 3m}$ are defined as $\mathrm{\Pi }=I_3\otimes {\mathrm{\Pi }}_{c\to s}$ and $N_s=\left(N_{s,1},N_{s,2},N_{s,3}\right)$, where ${\mathrm{\Pi }}_{c\to s}\in {\mathbb{R}}^{m\times n}$ is the operator that interpolates a cell-centered scalar field to the faces, and $N_{s,i}\in {\mathbb{R}}^{m\times m}$ are diagonal matrices containing the $x_i$-spatial components of the face-normal vectors.

A classical fractional step projection method [26, 27, 28] is employed to address the velocity-pressure coupling. The staggered velocity field, denoted as ${\mathbf{u}}_s$, a velocity ${\mathbf{u}}_s^p$ can be uniquely decomposed into a curl-free vector represented as the gradient of a scalar field, $\mathbf{G}{\tilde{\mathbf{p}}}_c^{\prime }$ and a solenoidal vector, ${\mathbf{u}}_s^{n+1}$. This decomposition is formulated as

$${\mathbf{u}}_s^p={\mathbf{u}}_s^{n+1}+\mathbf{G}{\tilde{\mathbf{p}}}_c^{\prime },$$E31

taking the divergence of Eq. (31) yields a discrete Poisson equation for ${\tilde{\mathbf{p}}}_c^{\prime },$

$${\mathbf{Mu}}_s^p={\mathbf{Mu}}_s^{n+1}+\mathbf{MG}{\tilde{\mathbf{p}}}_c^{\prime }\Rightarrow {\mathbf{Mu}}_s^p=\mathbf{MG}{\tilde{\mathbf{p}}}_c^{\prime }.$$E32

The cell-center predicted velocity field, ${\mathbf{u}}_c^p$, is computed with the projection method and then corrected to obtain the velocity at the next time step, ${\mathbf{u}}_{\mathbf{c}}^{\mathbf{n}\mathbf{+}\mathbf{1}}$. Assuming ${\mathbf{u}}_c^p\approx {\mathrm{\Gamma }}_{s\to c}{\mathrm{\Gamma }}_{c\to s}{\mathbf{u}}_c^p$ is satisfied, the overall correction algorithm can be explicitly written by combining the expression of Eq. (31) and the linear shift operator, i.e.,

$${\mathbf{u}}_c^{n+1}={\mathbf{u}}_c^p+{\mathrm{\Gamma }}_{s\to c}{\mathrm{\Omega }}_s^{-1}{\mathbf{M}}^{\ast }{\mathbf{L}}^{-1}\mathbf{M}{\mathrm{\Gamma }}_{c\to s}{\mathbf{u}}_c^p$$E33

$$={\mathbf{u}}_c^p+{\left(I_3\otimes {\mathrm{\Omega }}_c\right)}^{-1}{\mathrm{\Gamma }}_{c\to s}{\mathbf{M}}^{\ast }{\mathbf{L}}^{-1}\mathbf{M}{\mathrm{\Gamma }}_{c\to s}{\mathbf{u}}_c^p,$$E34

where the ${\mathrm{\Omega }}_s\in {\mathbb{R}}^{m\times m}$ and ${\mathrm{\Omega }}_c\in {\mathbb{R}}^{n\times n}$ are diagonal matrices with staggered and cell-centered control volumes, respectively. The discrete Laplacian operator $\mathbf{L}\in {\mathbb{R}}^{n\times n}$, defined as$L\equiv -\mathbf{M}{\mathrm{\Omega }}^{\mathbf{-}\mathbf{1}}{\mathbf{M}}^{\ast }$, is symmetric and negative-definite. The checkerboard problem stems from the unrealistic component of the cell-centered velocity field that cannot be eliminated by the pseudo-projection matrix,

$${\mathbf{u}}_c^{\ominus }={\mathbf{u}}_c^{\ominus }+\left({\left(I_3\otimes {\mathrm{\Omega }}_c\right)}^{-1}{\mathrm{\Gamma }}_{c\to s}{\mathbf{M}}^{\mathbf{\ast }}{L}^{-1}\mathbf{M}{\mathrm{\Gamma }}_{c\to s}\right){\mathbf{u}}_c^{\ominus },$$E35

where ${\mathbf{u}}_c^{\ominus }$ denotes the spurious modes. These’unrealistic’ velocity modes persist unless explicitly removed. Trias et al. [27] suggested regularizing (smoothing approximation) the convective term to address the origin of these ‘unrealistic’ velocity modes while maintaining the numerical solution devoid of nonphysical oscillations.

2.3.5 Constructing the discrete operators on unstructured collocated grids

Skew-symmetry convective operator $\mathbf{C}\left({\mathbf{u}}_{\mathbf{s}}\right)$ is obtained in two stages [10]. Firstly, attention is given to the off-diagonal elements. The matrix $\mathbf{C}\left({\mathbf{u}}_{\mathbf{s}}\right)$-dia $\left(\mathbf{C}\left({\mathbf{u}}_{\mathbf{s}}\right)\right)$ is skew-symmetry when the interpolation weights of the adjacent discrete variables are set equal to $1/2$ hence the discrete normal velocity ${\left[{\mathbf{u}}_s\right]}_f\approx {\mathbf{u}}_f\cdot {\mathbf{n}}_f$, situated at the centroid of the cell faces f, is expressed as

$${\left[{\mathbf{u}}_s\right]}_f={\left[{\mathbf{\Gamma }}_{c\to s}{\mathbf{u}}_c\right]}_f=\frac{1}{2}\left({\left[{\mathbf{u}}_c\right]}_{c1}+{\left[{\mathbf{u}}_c\right]}_{c2}\right)\cdot {\mathbf{n}}_f,$$E36

where $c1$ and $c2$ represent the cells neighboring the face f. The entries of the matrix ${\left[C\left({\mathbf{u}}_s\right)\right]}_k$ correspond to half of the flux through face f between the adjacent cells i and j, formulated as

$${\left[C\left({\mathbf{u}}_s\right)\right]}_k=\frac{1}{2}{\left[{\mathbf{u}}_s\right]}_f{A}_f,$$E37

where $A_f$ denotes the face area of f. In addition, for skew-symmetry, the diagonal elements should be zero,

$${\left[C\left({\mathbf{u}}_s\right)\right]}_{k,k}=\frac{1}{2}\sum _{f\in F_f\left(i\right)}{\left[{\mathbf{u}}_s\right]}_f{A}_f=0,$$E38

where the $F_f\left(i\right)$ represents the set of faces neighboring face f. This criterion is met due to the vanishing discrete divergence of ${\mathbf{u}}_s$. Thus, the convective operator unstructured collocated grids

$${\left[C\left({\mathbf{u}}_s\right){\mathbf{u}}_c\right]}_k=\sum _{f\in F_f\left(k\right)}\frac{1}{2}\left({\left[{\mathbf{u}}_c\right]}_{c1}+{\left[{\mathbf{u}}_c\right]}_{c2}\right){\left[{\mathbf{u}}_s\right]}_f{A}_f,$$E39

By integrating the continuity equation given in Eq. (16) over an arbitrary cell k with volume ${\left({\mathrm{\Omega }}_c\right)}_{\mathit{kk}}$, we obtain

$${\int }_{{\left({\mathrm{\Omega }}_c\right)}_{k,k}}\nabla \cdot \mathbf{u}\mathit{dV}={\int }_{\partial {\left({\mathrm{\Omega }}_c\right)}_{k,k}}\mathbf{u}\cdot \mathbf{n}\mathit{dS}=\sum _{f\in F_f\left(k\right)}{\int }_{S_f}\mathbf{u}\cdot \mathbf{n}\mathit{dS}.$$E40

A discretization of Eq. (40) that is second-order accurate is

$$\sum _{f\in F_f\left(k\right)}{\int }_{S_f}\mathbf{u}\cdot \mathbf{n}\mathit{dS}\approx \sum _{f\in F_f\left(k\right)}{\left[{\mathbf{u}}_s\right]}_f{A}_f.$$E41

Hence, the divergence operator is as

$${\left[{\mathbf{Mu}}_s\right]}_k\equiv \sum _{f\in F_f\left(k\right)}{\left[{\mathbf{u}}_s\right]}_f{A}_f=0.$$E42

According to Eq. (22) the integrated pressure gradient operator, ${\mathrm{\Omega }}_s\mathbf{G}$, equals the negative of the transpose of the divergence operator $-{\mathbf{M}}^{\mathbf{\ast }}$. Hence, the discretization of the pressure gradient at the face f follows from Eq. (42)

$${\left[\mathrm{\Omega }{\mathbf{Gp}}_c\right]}_f=\left(p_{c1}-p_{c2}\right)A_f,$$E43

The discrete gradient inherits boundary conditions from the discrete divergence, thus negating the need for specifying boundary conditions for the pressure. Finally, we compute the pressure from a Poisson equation that from the incompressibility constraint.

The Laplacian operator is represented by the matrix,

$$\mathbf{L}=-\mathbf{M}{\mathrm{\Omega }}_s^{-1}{\mathbf{M}}^{\mathbf{\ast }}.$$E44

negative-definite and symmetric matrix resembles the continuous Laplacian operator $\mathrm{\Delta }\equiv \nabla \cdot \nabla$.

The approach used to discretize the Laplacian operator in the pressure Poisson equation is also employed to discretize the diffusive term in NS equations. The diffusive operator is seen as the product of two first-order differential operators, a divergence and a gradient. The divergence operator is discretized and the gradient operator becomes the transpose of the discrete divergence

$${\mathbf{D}}_c=-\frac{1}{\mathit{Re}}\mathbf{MG}=\frac{1}{\mathit{Re}}\left(\mathbf{M}{\mathrm{\Omega }}_s^{-1}\right){\mathbf{M}}^{\mathbf{\ast }}.$$E45

This procedure results in a symmetric, positive-definite, approximation of the diffusive operator $-\nabla \cdot \nabla$. The collocated diffusive operator for a cell-centered variable ${\varphi }_c$ is expressed as,

$$\left[{\mathbf{D}}_c{\varphi }_c\right]=\frac{1}{\mathit{Re}}\sum _{f\in F_f\left(k\right)}\frac{\left({\varphi }_{c1}-{\varphi }_{c2}\right)A_f}{\delta n_f},$$E46

where the length $\delta n_f$ approximates the distance between the centroid of cells $c1$ and $c2$, calculated as $\delta n_f=\mid {\mathbf{n}}_f\cdot \overleftarrow{c1c2}\mid$. Afterwards, the volume of the cell associated with the face-normal velocity at face f is denoted as ${\left({\mathrm{\Omega }}_s\right)}_f=\delta n_f{A}_f$, where the $F_f\left(k\right)$ is the set of faces boarding the face f.

2.3.6 Solver

The RKSymFoam solver’s primary algorithm comprises three hierarchical iterative levels in the case of implicit time discretization [12, 13]:

1. The outer loop iterates over each RungeKutta (RK) stage i, ranging from 1 to s (see Section 2.3.8 for details).

2. A loop manages to update the nonlinear convective term.

3. An inner PISO loop handles the pressure-velocity coupling.

For explicit temporal discretization, a single projection step suffices for the pressure-velocity coupling, and there’s no need to update the convective term.

2.3.7 Temporal discretization

To maintain the favorable conservation and stability properties in discrete time, the time discretization in Eq. (25) is a skew-symmetric operator, implicit temporal integration. In high-Reynolds-number flow simulations, the implicit method may incur higher computational costs compared to the explicit counterpart. Yet in an explicit temporal discretization scheme for convection-diffusion equations, the time step $\mathrm{\Delta }t$ is typically constrained by convective stability conditions such as $U\mathrm{\Delta }t<\mathrm{\Delta }y$ (where $\mathit{U}$ represents the maximum velocity magnitude and $\mathrm{\Delta }y$ denotes the spatial grid size), diffusive stability conditions akin to $2\mathrm{\Delta }t<\mathit{Re}\mathrm{\Delta }{y}^2$. The implicit method does not impose such stringent restrictions, allowing for a time step that can be larger. To strike a balance between accuracy and stability, ensuring that the simulation remains globally second-order accurate, the Crank-Nicolson and diagonal-implicit Runge–Kutta methods are chosen.

2.3.8 Butcher tables specify Runge-Kutta methods

After applying the spatial discretization the momentum equations Eq. (16) become

$$\frac{d{\mathbf{u}}_c}{\mathit{dt}}={\mathbf{F}}_c\left(t,{\mathbf{u}}_c\right)-{\left(\nabla p\right)}_c\text{with}{\mathbf{u}}_c\left(t^n\right)={\mathbf{u}}_c^n.$$E47

Here, the subscript c denotes the cell-centered values, while the pressure gradient term remains to be discretized. Runge-Kutta (RK) schemes, specified by a Butcher table [29], are used to discretize Eq. (47) in time. Using $a_{\mathit{ij}}$ to denote the stage weights for stage i, and $c_i$ as the quadrature nodes for the schemes with $c_i={\sum }_{j=1}^s{a}_{\mathit{ij}}$, where $i=1,\dots ,s$, and $t^i=t^n+c_i\mathrm{\Delta }t$, s denotes the number of stages, and $b_j$ denotes the primary weights of the utilized RK scheme, satisfying ${\sum }_{j=1}^s{b}_j=1$. Explicit schemes are characterized by having only non-zero entries in the lower triangle of the $a_{i,j}$ portion of the table. Entries on or above the diagonal lead to implicit methods as they involve ${\mathbf{u}}_c^{n+1}$ on the right-hand side of Eq. (47).

The intermediate solution ${\tilde{\mathbf{u}}}_c^i$ for stage i at time $t^i$ is expressed as

$${\tilde{\mathbf{u}}}_c^i={\mathbf{u}}_c^n+\mathrm{\Delta }t\left(\sum _{j=1}^i{a}_{\mathit{ij}}{\mathbf{F}}_c\left(t^j,{\tilde{\mathbf{u}}}_c^j\right)-c_i{\left(\nabla \tilde{p}\right)}_c^i\right),$$E48

and the final solution ${\mathbf{u}}_c^{n+1}$ at time $t^{n+1}=t^n+\mathrm{\Delta }t$ by

$${\mathbf{u}}_c^{n+1}={\mathbf{u}}_c^n+\mathrm{\Delta }t\left(\sum _{j=1}^i{b}_j{\mathbf{F}}_c\left(t^j,{\tilde{\mathbf{u}}}_c^j\right)-{\left(\nabla p\right)}_c^{n+1}\right).$$E49

Here we concentrate on explaining the forward-Euler method as an example, as higher-order explicit RK methods essentially involve a series of forward-Euler stages. We compute an intermediate velocity field ${\mathbf{u}}_c^p$ through the following predictor step

$$\frac{{\mathbf{u}}_c^p-{\mathbf{u}}_c^n}{\mathrm{\Delta }{t}^n}=\mathbf{H}\left({\mathbf{u}}_s^n\right){\mathbf{u}}_c^n-{\mathbf{G}}_c{\mathbf{p}}_c^p.$$E50

Here $\mathbf{H}\left({\mathbf{u}}_s^n\right){\mathbf{u}}_c^n\equiv -{\mathrm{\Omega }}_c^{-1}\left(\mathbf{C}\left({\mathbf{u}}_s^n\right){\mathbf{u}}_c^n+{\mathbf{Du}}_c^n\right)$, where $\mathbf{H}\left({\mathbf{u}}_s^n\right)\in {\mathbb{R}}^{3n\times 3n}$, and $\mathrm{\Delta }{t}^n$ denotes the time difference between time levels n and $n+1$. As the intermediate velocity is not divergence-free, we correct the final values for the time step $n+1$ by adding the following corrections to the intermediate values: ${\mathbf{u}}_c^{n+1}={\mathbf{u}}_c^p+{\mathbf{u}}_c^{\prime }$ and ${\mathbf{p}}_c^{n+1}={\mathbf{p}}_c^p+{\tilde{\mathbf{p}}}_c^{\prime }$, the velocity and pressure correction

$${\mathbf{u}}_c^{\prime }=-\mathrm{\Delta }t{\mathbf{G}}_c{\tilde{\mathbf{p}}}_c^{\prime }.$$E51

After obtaining the pressure correction ${\tilde{\mathbf{p}}}^{\prime }$ from the Poisson equation Eq. (32), we can calculate the velocity correction ${\mathbf{u}}^{\prime }$ using Eq. (51). Subsequently, the new velocity and pressure fields for the next time step can be computed using the velocity and pressure corrections ${\mathbf{u}}_c^{\prime }$ and ${\tilde{\mathbf{p}}}_c^{\prime }$.

The RKSymFoam solver employs the PISO approach for implicit time integration. Let us illustrate this with the backward Euler method. The PISO algorithm comprises a predictor step followed by $n_{\mathit{\text{corr}}}$ corrector steps (or inner iterations), as outlined below [30]:

1. Predictor step: The initial intermediate velocity ${\mathbf{u}}_c^1$ is calculated using the predictor equation

   $$\frac{{\mathbf{u}}_c^1-{\mathbf{u}}_c^n}{\mathrm{\Delta }{t}^n}=\mathbf{H}\left({\mathbf{u}}_s^p\right){\mathbf{u}}_c^p-{\mathbf{G}}_c{\mathbf{p}}_c^p,$$E52

   where $\mathbf{H}\left({\mathbf{u}}_s^p\right){\mathbf{u}}_c^p=\mathbf{H}\left({\mathbf{u}}_s^n\right){\mathbf{u}}_c^n$ is utilized to refrain from an implicit treatment of the convective and diffusive terms. The resulting initial intermediate velocity may not be divergence-free, necessitating subsequent corrector steps.

2. Corrector steps: In each correct step, a new pressure field ${\mathbf{p}}_c^k$ and a corresponding revised velocity ${\mathbf{u}}_c^{k+1}$ (which is divergence-free $\mathbf{M}{\mathrm{\Gamma }}_{c\to s}{\mathbf{u}}_c^{k+1}=0_{c,n}$) are. To enhance the stability of the momentum equation during the corrector step, the operator is partitioned into a diagonal component $A^d\left({\mathbf{u}}_s^n\right)$ acting on ${\mathbf{u}}_c^{k+1}$ and an off-diagonal component ${\mathbf{H}}^{\mathit{od}}\left({\mathbf{u}}_s^n\right)$ acting on ${\mathbf{u}}_c^k$, resulting in

$${\mathbf{u}}_c^{k+1}={\mathbf{B}}^{-1}\left(\mathrm{\Delta }t{\mathbf{H}}^{\mathit{od}}\left({\mathbf{u}}_s^n\right){\mathbf{u}}_c^k+{\mathbf{u}}_c^n\right)-\mathrm{\Delta }t{\mathbf{B}}^{-1}{\mathbf{G}}_c{\mathbf{p}}_c^k,\text{with}k=1,2,\dots ,n_{\mathit{\text{corr}}},$$E53

where $\mathbf{B}=\mathbf{I}-\mathrm{\Delta }{t}^n{A}^d\left({\mathbf{u}}_s^n\right)\in {\mathbb{R}}^{3n\times 3n}$. A Poisson equation for pressure ${\mathbf{p}}_c^k$ by taking the divergence of Eq. (53) and employing $\mathbf{M}{\mathrm{\Gamma }}_{c\to s}{\mathbf{u}}_c^{k+1}=0_{c,n}$

$${\mathbf{MB}}_S^{-1}{\mathrm{\Omega }}_s^{-1}{\mathbf{M}}^{\mathbf{\ast }}{\mathbf{p}}_c^k=\frac{1}{\mathrm{\Delta }t}\mathbf{M}{\mathrm{\Gamma }}_{c\to s}{\mathbf{B}}^{-1}\left(\mathrm{\Delta }t{\mathbf{H}}^{\mathit{od}}\left({\mathbf{u}}_s^n\right){\mathbf{u}}_c^k+{\mathbf{u}}_c^n\right),k=1,2,\dots ,n_{\mathit{\text{corr}}},$$E54

where ${\mathbf{B}}_s^{-1}\approx {\mathrm{\Gamma }}_{c\to s}{\mathbf{B}}^{-1}{\mathrm{\Gamma }}_{s\to c}$, is a square diagonal matrix. The implicit scheme is$\mathbf{L}=-{\mathbf{MB}}_s^{-1}{\mathrm{\Omega }}_s^{-1}{\mathbf{M}}^{\mathbf{\ast }}$. The pressure correction ${{\tilde{\mathbf{p}}}^{\prime }}_c^k$ is obtained from the following Poisson equation

$${\mathbf{MB}}_S^{-1}{\mathrm{\Omega }}_s^{-1}{\mathbf{M}}^{\ast }{{\mathbf{p}}^{\prime }}_c^k=\frac{1}{\mathrm{\Delta }t}\mathbf{M}{\mathrm{\Gamma }}_{c\to s}{\mathbf{B}}^{-1}\left(\mathrm{\Delta }t{\mathbf{H}}^{\mathit{od}}\left({\mathbf{u}}_s^n\right){\mathbf{u}}_c^k+{\mathbf{u}}_c^n-{\mathbf{Gp}}_c^n\right),k=1,..,n_{\mathit{\text{corr}}},$$E55

After obtaining the pressure field ${\mathbf{p}}_c^k$ from the Poisson equation Eq. (54), the corresponding updated velocity can be computed using Eq. (53). This iterative process continues for $n_{\mathit{\text{corr}}}$ iterations until an inner iteration criterion is met. The inner PISO iteration procedure concludes by updating the velocity, which is utilized to solve the non-linearity in the convective the operator H, with the new velocity ${\mathbf{u}}_c^{n_{\mathit{\text{corr}}}}$ for the subsequent outer iteration.

3. Verification and validation

The test cases encompass channel flow across a range of friction Reynolds numbers $\left({\mathit{Re}}_{\tau }\right)$, including 180, 550, 1000, and 2000, along with flow over periodic hills and flow past a circular cylinder. Both static and dynamic minimum-dissipation models are subjected to numerical examinations in channel flow scenarios. The investigation delves into the sensitivity of the QR model concerning spatial discretization schemes, mesh resolution, and Reynolds number variations, comparing standard OpenFOAM discretization against symmetry-preserving discretization methods. This section focuses on the main findings; for comprehensive details, see [31].

3.1 Optimal QR model coefficient

The simulations were performed using the QR model with different model constants for a channel flow at ${\mathit{Re}}_{\tau }=180$ and periodic hills, aiming to identify the optimal model constant. Figure 1 depicts the normalized mean streamwise velocity against the wall. It demonstrates that the model coefficient C = 0.024 aligns precisely with DNS data. Additionally, utilizing five different error measurements (mean square, absolute, maximum absolute, and slope errors), the constant C = 0.024 was identified as yielding the smallest error. This model coefficient yields accurate results under mesh convergence. Furthermore, the various comparisons carried out for flows over periodic hills demonstrate that the model coefficient of C = 0.024 is again the best one. In literature, the Smagorinsky model coefficient is typically considered optimal within the range of $C_s=0.1$ to $C_s=0.2$. Consequently, the ideal coefficient for the QR model falls within the range of $C_s^2$. Note that the model coefficient utilized in the QR model is the Smagorinsky.

Figure 1.

Normalized mean streamwise velocity versus wall-normal distance in channel flow simulation at ${\mathit{Re}}_{\tau }=\mathit{180}$_.

3.2 Static QR model compared to dynamic models in channel flow at ${\mathbf{Re}}_{\mathbf{\tau }}\mathbf{=}\mathbf{180}$

The static QR model is compared with four dynamic models. The dynamic procedure is based on the Lilly-Germano identity [2]. All subgrid models depicted in Figure 2 closely approximate the DNS results in the near-wall region, particularly the dynamic QR model, which agrees with the DNS data. A detailed examination of the figure highlights that in the logarithmic wall range $y^{+}>30$, the static QR, dynamic minimum-dissipation (DMD), and dynamic QR (DQR) models closely match the DNS results. Conversely, the hybrid dynamic (HDM) and dynamic Smagorinsky (DSM) models nearly overlap and deviate from the DNS data in the logarithmic wall region.

Figure 2.

Average streamwise velocity $u^{+}$ obtained from static QR and dynamic models compared to DNS data.

, both the static QR and DQR models exhibit similar performance, with negligible differences in errors. On the other hand, the DSM and HDM demonstrate comparatively lower accuracy, while the DMD model exhibits the largest error.

In summary, this study underscores the reliability of the static QR model. A suitably chosen model coefficient can yield results akin to a dynamic model while cutting computational costs by a factor of three.

3.3 Sensitivity to expansion ratio in channel flow at ${\mathbf{Re}}_{\mathbf{\tau }}\mathbf{=}\mathbf{180}$

This investigation focuses on the expansion ratio of the mesh in the wall-normal direction. The expansion ratio is defined as follows

$$\text{Expansion\hspace{0.17em}\hspace{0.17em}ratio}=\frac{\mathrm{\Delta }{y}_c^{+}}{\mathrm{\Delta }{y}_w^{+}},$$E56

where $\mathrm{\Delta }{y}_c^{+}$ represents the size of the cell at the channel center, while $\mathrm{\Delta }{y}_w^{+}$ denotes the size of the first cell adjacent to the wall. Figure 3 illustrates the normalized mean streamwise and rms velocity against wall distance in wall unit at ${\mathit{Re}}_{\tau }=180$, employing a grid resolution of $40\times 50\times 30$. The normalized uniform grid spacing in the streamwise direction is $\mathrm{\Delta }{x}^{+}=18$ and in spanwise direction is $\mathrm{\Delta }{z}^{+}=12$. The expansion ratios in the wall-normal direction are (approximately) 5, 10, 15, 20, and 25. The first normalized wall-normal grid height is $\mathrm{\Delta }{y}_w^{+}=1.728$. The time step for simulation is $\mathrm{\Delta }{t}^{+}=36$. As the expansion ratio increases, the near-wall mesh, while the mesh in the channel center.

Figure 3.

Normalized mean streamwise and rms velocity using QR model with C = 0.023 at ${\mathit{Re}}_{\tau }=\mathit{180}$ expansion ratios: 5,10,15, 20, and 25.

For the mean streamwise velocity $u^{+}$, reducing the expansion ratio from 15 to 5 leads to increasingly pronounced overestimation of the mean velocity beyond $y^{+}>10$. However, this overestimation exacerbates further with expansion ratios of 20 and 25.

For the velocity fluctuation in the streamwise direction $u^{\prime }{u}^{\prime }$, lowering the expansion ratio from 25 to 5 results in notably exaggerated peaks in the near-wall region and relatively diminished overestimation in the central channel. For the velocity fluctuation in the wall-normal ${\upsilon }^{\prime }{\upsilon }^{\prime }$ and spanwise $w^{\prime }{w}^{\prime }$ directions (not shown here), reducing the expansion ratio from 20 to 5 causes overestimation of peak values near the wall and comparatively smaller underestimations in the central channel. The highest expansion ratio, 25, underestimates the near-wall peak but provides more accurate values in the channel center.

3.4 Grid Independence study in channel flow at ${\mathbf{Re}}_{\mathbf{\tau }}\mathbf{=}\mathbf{550}$

This section of the survey investigates how different grid configurations affect the mean and root-mean-square velocities at a Reynolds number of ${\mathit{Re}}_{\tau }=550$, utilizing a QR model with C = 0.023. Four grid setups are considered: ${128}^3$, ${64}^3$, ${32}^3$, and ${16}^3$. Figure 4 illustrates the normalized mean streamwise velocity as a function of wall distance in wall units, with the channel dimensions being $2\pi \times 2\times \pi$. The expansion in the wall-normal direction is approximately 10, and the simulation employs a time step of $\mathrm{\Delta }{t}^{+}=55$.

Figure 4.

Comparison of streamwise velocities using QR Model (C = 0.023) at ${\mathit{Re}}_{\tau }=\mathit{550}$ on various mesh resolutions (${128}^3$, ${64}^3$, ${32}^3$, ${16}^3$) with DNS results at ${\mathit{Re}}_{\tau }=\mathit{500}$.

Initially, the analysis focuses on the results of the simulation conducted using the ${16}^3$ grid points (depicted by the blue line on the left side of Figure 4). The first mesh point is situated around $y^{+}=7$, indicating a significant distance from the wall. Despite the coarse mesh, the obtained results surprisingly exhibit reasonable agreement with the overall trend observed in the other simulations.

A noteworthy observation is the discernible trend that emerges with increasing mesh resolution: the mean velocity profile progressively aligns more closely with DNS data, particularly evident when utilizing ${128}^3$ mesh points. As resolution increases, the proximity of the initial mesh points to the wall decreases, resulting in a greater concentration of mesh in the near-wall region. Consequently, the ${128}^3$ case demonstrates complete agreement with DNS data in the near-wall region.

3.5 Symmetry-preserving compared to standard OpenFOAM discretization

In this section, we illustrate the effect of the numerical schemes by comparing the symmetry-preserving discretization with the Pimple algorithm conjugated with Gauss linear schemes in three benchmarks, that is, periodic channel flow, periodic hill flow, and flow over a cylinder.

The channel flow outcomes depicted in Figure 5, the results obtained symmetry-preserving techniques closely resemble those of the DNS reference. Notably, in the region where $y^{+}\ge 10$, the symmetry-preserving scheme with the QR model precisely aligns with the DNS results.

Figure 5.

The mean streamwise velocity obtained from the QR model with standard OpenFOAM and symmetry-preserving (SP) schemes for fully developed turbulent channel flow at ${\mathit{Re}}_{\tau }=\mathit{1000}$.

The symmetry-preserving method enhances velocity fluctuation predictions in the region where $y^{+}\le 220$. In other areas of the computational domain, both methods yield comparable results. For the Reynolds shear stress $u^{\prime }{\upsilon }^{\prime }$, the two methods give the same accurate prediction.

Additionally, it has been observed that the contribution of the subgrid-scale model to the diffusive flux in both methods is lower than 20% of the contribution of molecular viscosity [31].

3.6 Application to periodic hill flow

The phenomenon of flow separation from curved surfaces followed by reattachment is frequently encountered in engineering applications. To evaluate the effectiveness of the proposed minimum-dissipation model in OpenFOAM for simulating separated flows, three-dimensional flow simulations over periodic hills at Re = 10,595 as, illustrated in Figure 6, have been conducted.

Figure 6.

The side view of the periodic hill flow.

The mean velocity and Reynolds stress calculated by the QR model, without incorporating a wall damping function, are compared with experimental data [32] and three numerical experiments [33]. It is important to note that all large-eddy simulations are based on second-order accurate numerical techniques both in space and time. Additionally, the QR model grid comprises six times fewer grid points compared to the grid used in the reference CFD data. The QR simulation mesh comprises 2.56 million grid points and necessitates 6 hours for simulation on two nodes, each equipped with 128 cores. This mesh size proves adequate for modeling this case, as further refinement does not yield improved results.

3.6.1 Separation and reattachment lengths

The separation and reattachment points are obtained at Re = 11,230, that is, the bulk velocity $u_b=1.06$, using Gauss linear spatial discretization and backward temporal discretization. The separation point is accurately pinpointed by numerically solving boundary layer equations under pressure-adverse conditions. The QR model predicts the separation point, where the wall shear stress reaches zero, to be approximately $x/H\approx 0.175$, slightly smaller than the reference value of $x/H\approx 0.19$. This deviation is expected, as the separation point tends to shift upstream with increasing Reynolds numbers. The precise location of the separation point significantly influences the point of reattachment. Recirculation initiates around $x/H\approx 0.27$ and extends until approximately $x/H\approx 4.71$, with the main recirculation bubble spanning a length of roughly 4.48. The reattachment position, where the dividing streamline reconnects with the wall, is determined to be at $x/H=4.71$, closely aligning with the reference value of 4.69 [32].

3.6.2 Sensitivity to the numerical schemes of standard OpenFOAM

A significant discrepancy has been identified in the prediction of the streamwise velocity component $⟨u⟩/u_b$ within the range $1.5<y/H<2.8$, despite employing the same Reynolds number as the reference case (Re = 10, 595).

In an effort to understand the cause of this discrepancy, various interpolation numerical schemes for solving the divergence term and pressure gradient have been evaluated. These schemes include central difference, filtered and limited central difference, linear upwind blended with central difference, least squares, and detached eddy discretization schemes. However, the outcomes of these different numerical schemes fail to accurately predict the mean velocity, as illustrated in Figure 7. Consequently, it has been determined that increasing the Reynolds number to better align with the reference data is necessary. Furthermore, despite minor differences observed between the simulations, the central difference schemes are ultimately selected due to their satisfactory performance in previous validation cases. These discrepancies, while present, are deemed insignificant in distinguishing one technique from another.

Figure 7.

Streamwise velocity at five different positions with varying finite volume discretization schemes.

3.6.3 Sensitivity to the Reynolds number

In order to address the underprediction of the streamwise velocity component $⟨u⟩/u_b$ within the range of $1.5<y/H<2.8$ while maintaining the Reynolds number equal to the reference value of Re = 10,595, simulations at higher Reynolds numbers were explored. The most notable revelation from the analysis, as depicted in Figure 8, is that the mean velocity and Reynolds stress obtained with a 6.5% higher Reynolds number align with the B-spline reference. Furthermore, a 7% increase in Reynolds number corresponds to agreement with experimental data (not illustrated here). These findings confirm that augmenting the Reynolds number by 3−7% successfully aligns the simulations with the reference data. Note that despite the experimental data and three LES simulations used for comparison being conducted at the same Reynolds number of Re = 10,595, they exhibit discrepancies among themselves.

Figure 8.

Streamwise velocity (upper) and Reynolds stress (lower) at ten distinct locations, with a 6.5% increment in Reynolds number compared to the B-spline data. 11,284 refers to the QR results and LES-1 indicates the reference B-spline results.

3.7 Application to flow over a cylinder

3.7.1 Comparison between standard OpenFOAM spatial discretization

To examine the impact of different numerical schemes employed in standard OpenFOAM spatial discretization on turbulence, simulations are conducted for both periodic hills and flow over a cylinder.

Figure 9 displays the mean streamwise velocity u resulting from three distinct numerical discretizations applied to flow over a cylinder. Notably, discrepancies in mean streamwise velocity u are observed in the vicinity of the wake $\left(x/D<2.02\right)$. Specifically, the upwind-blended scheme predicts a U-shaped mean velocity profile at $x/D=1.06$, transitioning to a markedly lower V-shaped profile at $x/D=1.54$ and 2.02, compared to the other two schemes.

Figure 9.

The mean streamwise velocity at six locations (left) and the Reynolds stress in the cross-section at $x/D=\mathit{1.54}$ (right) are depicted for flow over a cylinder at ${\mathit{Re}}_D=\mathit{3900}$. The QR model is utilized with filtered central difference (Flinear), pure central difference (Linear), and upwind-biased (LUST) schemes.

Analysis from the central difference simulation reveals that the onset of turbulence in the separating shear layers occurs nearer to the cylinder, leading to the formation of a shorter vortex region and the development of a V-shaped velocity profile. Consequently, the shear layers exhibit reduced length, and the recirculation region diminishes in size [31]. In the downstream region, variations in numerical schemes appear to have minimal impact on mean velocity. However, simulations of periodic hills demonstrate that the selection of numerical schemes substantially influences both mean and root-mean-square variables.

4. Conclusions

A comprehensive comparison between the minimum-dissipation model (QR) and experimental/DNS/LES results across channel flow, flow past a circular cylinder, and flow over periodic hills reveals generally favorable agreement, particularly on coarse meshes.

The findings suggest that the static QR model is comparable to dynamic models while reducing computational costs. Notably, a model coefficient of C = 0.024 produces the most accurate predictions, with the contribution of the subgrid model diminishing to less than 20% of the molecular viscosity as mesh resolution increases to its finest level.

The characteristics of turbulence exhibit dependence on spatial discretization methods. Various comparisons highlight the superiority of the QR model when paired with symmetry-preserving discretization over standard OpenFOAM discretization techniques. Among standard OpenFOAM discretization methods, central difference schemes emerge as the most effective. Finally, for the periodic hill flow, augmenting the Reynolds number by 3−7% successfully aligns the simulations with the reference data.

Acknowledgments

This work was supported by the Chinese Scholarship Council and the University of Groningen under Grant (CSC No. 202006870021). The computing time was granted by the Dutch Research Council (NWO) under project 2022.009.

The first version was published as a preprint titled ‘Minimum-dissipation model for large-eddy simulation in OpenFOAM – A study on channel flow, periodic hills, and flow over a cylinder’ on arXiv.org in 2023.