A Finite Volume approximation of the Navier-Stokes equations with nonlinear filtering stabilization

TL;DR

This paper introduces a finite volume approach combined with a nonlinear filtering stabilization for simulating incompressible fluid flows at moderate Reynolds numbers, validated against benchmark and experimental data.

Contribution

It proposes integrating the Evolve-Filter-Relax algorithm with a finite volume method, providing practical guidance for parameter tuning and mesh feature effects.

Findings

Finite volume method yields different results than finite element method at similar mesh levels.

Parameter tuning can improve model accuracy and stability.

Mesh features significantly influence simulation accuracy.

Abstract

We consider a Leray model with a nonlinear differential low-pass filter for the simulation of incompressible fluid flow at moderately large Reynolds number (in the range of a few thousands) with under-refined meshes. For the implementation of the model, we adopt the three-step algorithm Evolve-Filter-Relax (EFR). The Leray model has been extensively applied within a Finite Element (FE) framework. Here, we propose to combine the EFR algorithm with a computationally efficient Finite Volume (FV) method. Our approach is validated against numerical data available in the literature for the 2D flow past a cylinder and against experimental measurements for the 3D fluid flow in an idealized medical device, as recommended by the U.S. Food and Drug Administration. We will show that for similar levels of mesh refinement FV and FE methods provide significantly different results. Through our numerical experiments, we are able to provide practical directions to tune the parameters involved in the model. Furthermore, we are able to investigate the impact of mesh features (element type, non-orthogonality, local refinement, and element aspect ratio) and the discretization method for the convective term on the agreement between numerical solutions and experimental data.

Figures

Tables

Table 1. Table 1: Name, minimum diameter h m ​ i ​ n subscript ℎ 𝑚 𝑖 𝑛 h_{min} , maximum diameter h m ​ a ​ x subscript ℎ 𝑚 𝑎 𝑥 h_{max} , and number of cells for all the prismatic ( P 𝑃 P ) or hexahedral ( H 𝐻 H ) and meshes used for the 2D flow past a cylinder.

mesh name	$h_{min}$	$h_{max}$	No. of cells	mesh name	$h_{min}$	$h_{max}$	No. of cells
$16k_P$	7.4e-3	1.5e-2	1.6e4	$16k_H$	4.2e-3	1.1e-2	1.59e4
$25k_P$	5.9e-3	1.3e-2	2.5e4	$25k_H$	3.5e-3	9.6e-3	2.45e4
$63k_P$	3.7e-3	8.1e-3	6.33e4	$63k_H$	1.9e-3	5.8e-3	6.28e4
$120k_P$	2.7e-3	6.2e-3	1.21e5	$120k_H$	1.4e-3	4.2e-3	1.22e5
$200k_P$	2e-3	4.8e-3	1.98e5	$200k_H$	1e-3	3.4e-3	1.98e5

Table 2. Table 2: Maximum lift and drag coefficients given by NSE and times at which the maxima occur on all the meshes under consideration. The bottom row reports the results from [ 26 ] .

Mesh	$t\left(c_{l,max}\right)$	$c_{l,max}$	$t\left(c_{d,max}\right)$	$c_{d,max}$	Mesh	$t\left(c_{l,max}\right)$	$c_{l,max}$	$t\left(c_{d,max}\right)$	$c_{d,max}$
$16k_H$	5.79	0.574	3.947	3.074	$16k_P$	5.769	0.505	3.928	2.906
$25k_H$	5.743	0.533	3.929	3.057	$25k_P$	5.80	0.558	3.946	3.066
$63k_H$	5.739	0.515	3.945	3.028	$63k_P$	5.737	0.561	3.943	3.043
$120k_H$	5.724	0.507	3.936	3.012	$120k_P$	5.737	0.484	3.935	3.012
$200k_H$	5.716	0.495	3.934	2.986	$200k_P$	5.726	0.478	3.939	2.994
$Q_2$-$P_1^d$	5.694	0.478	3.936	2.951	$P_2$-$P_1$	5.694	0.478	3.936	2.951

Table 3. Table 3: Maximum lift and drag coefficients and times at which the maxima occur for the EF L, EF NL and EFR algorithms with α = h m ​ i ​ n , η 𝛼 subscript ℎ 𝑚 𝑖 𝑛 𝜂 \alpha=h_{min},\eta on all the meshes under consideration.

Mesh name	Algorithm	$t\left(c_{l,max}\right)$	$c_{l,max}$	$t\left(c_{d,max}\right)$	$c_{d,max}$	$E_{c_l}$	$E_{c_d}$
$16k_H$	EF L, $\alpha$ = $h_{min}$	4.024	0.258	3.735	0.984	-47.87%	-67.04%
$16k_H$	EF NL, $\alpha$ = $h_{min}$	2.303	0.022	3.943	3.358	-95.55%	12.45%
$16k_H$	EF L, $\alpha$ = $\eta$	4	0.209	3.784	1.144	-57.77%	-61.68%
$16k_H$	EF NL, $\alpha$ = $\eta$	3.511	0.041	3.942	3.342	-91.71%	11.92%
$16k_H$	EFR, $\alpha$ = $h_{min}$	5.777	0.544	3.94	3.074	9.89%	2.94%
$16k_H$	EFR, $\alpha$ = $\eta$	5.777	0.545	3.943	3.074	10.1%	2.94%
$16k_P$	EF L, $\alpha$ = $h_{min}$	3.989	0.24	3.581	0.734	-49.79%	-75.48%
$16k_P$	EF NL, $\alpha$ = $h_{min}$	3.895	0.049	3.864	2.7	-89.74%	-9.82%
$16k_P$	EF L, $\alpha$ = $\eta$	4	0.173	3.762	1.026	-63.8%	-65.73%
$16k_P$	EF NL, $\alpha$ = $\eta$	0.905	0.0025	3.918	2.914	-99.47%	-2.67%
$16k_P$	EFR, $\alpha$ = $h_{min}$	5.76	0.474	3.931	2.905	-0.83%	-2.97%
$16k_P$	EFR, $\alpha$ = $\eta$	5.761	0.481	3.926	2.906	0.62%	-2.94%
$25k_H$	EF L, $\alpha$ = $h_{min}$	4.008	0.168	3.738	0.993	-66.06%	-66.74%
$25k_H$	EF NL, $\alpha$ = $h_{min}$	3.918	0.038	3.942	3.683	-92.32%	23.34%
$25k_H$	EF L, $\alpha$ = $\eta$	4.012	0.157	3.763	1.039	-68.28%	-65.2%
$25k_H$	EF NL, $\alpha$ = $\eta$	4.249	0.036	3.974	3.641	-92.72%	21.93%
$25k_H$	EFR, $\alpha$ = $h_{min}$	5.734	0.509	3.936	3.058	2.82%	2.41%
$25k_H$	EFR, $\alpha$ = $\eta$	5.734	0.509	3.935	3.058	2.82%	2.41%
$25k_P$	EF L, $\alpha$ = $h_{min}$	4.003	0.234	3.608	0.756	-51.04%	-74.75%
$25k_P$	EF NL, $\alpha$ = $h_{min}$	3	0.021	3.918	3.526	-95.6%	17.76%
$25k_P$	EF L, $\alpha$ = $\eta$	4.011	0.205	3.732	0.954	-57.11%	-68.13%
$25k_P$	EF NL, $\alpha$ = $\eta$	6.078	0.01	3.936	3.546	-97.9%	18.43%
$25k_P$	EFR, $\alpha$ = $h_{min}$	5.792	0.534	3.942	3.067	11.71%	2.43%
$25k_P$	EFR, $\alpha$ = $\eta$	5.792	0.538	3.931	3.067	12.55%	2.43%
$63k_H$	EF L, $\alpha$ = $h_{min}$	4.03	0.113	3.803	1.166	-77.17%	-60.95%
$63k_H$	EF NL, $\alpha$ = $h_{min}$	5.994	0.246	3.961	3.573	-50.3%	19.65%
$63k_H$	EF L, $\alpha$ = $\eta$	4.031	0.186	3.701	0.89	-62.42%	-70.19%
$63k_H$	EF NL, $\alpha$ = $\eta$	3.571	0.05	3.944	3.888	-89.89%	30.2%
$63k_H$	EFR, $\alpha$ = $h_{min}$	5.733	0.502	3.942	3.028	1.41%	1.4%
$63k_H$	EFR, $\alpha$ = $\eta$	5.734	0.502	3.941	3.029	1.41%	1.44%
$63k_P$	EF L, $\alpha$ = $h_{min}$	4.006	0.222	3.656	0.806	-53.55%	-73.08%
$63k_P$	EF NL, $\alpha$ = $h_{min}$	3.947	0.042	3.95	4.118	-91.21%	37.54%
$63k_P$	EF L, $\alpha$ = $\eta$	4	0.217	3.688	0.853	-54.6%	-71.5%
$63k_P$	EF NL, $\alpha$ = $\eta$	3.98	0.035	3.95	4.062	-92.67%	35.67%
$63k_P$	EFR, $\alpha$ = $h_{min}$	5.735	0.546	3.945	3.044	14.22%	1.67%
$63k_P$	EFR, $\alpha$ = $\eta$	5.734	0.547	3.935	3.044	14.43%	1.67%
$120k_H$	EF L, $\alpha$ = $h_{min}$	4.016	0.112	3.813	1.2	-77.37%	-59.81%
$120k_H$	EF NL, $\alpha$ = $h_{min}$	5.891	0.391	4	3.37	-21.01%	12.86%
$120k_H$	EF L, $\alpha$ = $\eta$	3.994	0.218	3.622	0.79	-55.96%	-73.54%
$120k_H$	EF NL, $\alpha$ = $\eta$	2.833	0.039	3.948	4.105	-92.12%	37.47%
$120k_H$	EFR, $\alpha$ = $h_{min}$	5.72	0.498	3.936	3.012	0.6%	0.87%
$120k_H$	EFR, $\alpha$ = $\eta$	5.72	0.498	3.943	3.013	0.6%	0.9%
$120k_P$	EF L, $\alpha$ = $h_{min}$	4.038	0.215	3.674	0.832	-55.02%	-72.21%
$120k_P$	EF NL, $\alpha$ = $h_{min}$	6.883	0.008	3.966	4.04	-98.32%	34.93%
$120k_P$	EF L, $\alpha$ = $\eta$	4.04	0.224	3.658	0.791	-53.13%	-73.58%
$120k_P$	EF NL, $\alpha$ = $\eta$	0.926	0.0014	3.977	4.189	-99.7%	39.91%
$120k_P$	EFR, $\alpha$ = $h_{min}$	5.734	0.474	3.94	3.013	-0.83%	0.63%
$120k_P$	EFR, $\alpha$ = $\eta$	5.734	0.474	3.944	3.014	-0.83%	0.66%

Table 4. Table 4: Throat Reynolds number R ​ e t 𝑅 subscript 𝑒 𝑡 Re_{t} , inlet Reynolds number R ​ e i 𝑅 subscript 𝑒 𝑖 Re_{i} , flow rate Q 𝑄 Q , Kolmogorov length scale η 𝜂 \eta and number of cells for a DNS for all the flow regimes under consideration.

$R{e}_t$	$R{e}_i$	$Q$ (m3/s)	$\eta$ (m)	No. of cells for DNS
2000	667	2.0825e-5	4.0124e-5	3.6285e8
3500	1167	3.6444e-5	2.6371e-5	1.2781e9
5000	1667	5.2062e-5	2.0182e-5	2.8513e9
6500	2167	6.7681e-5	1.6577e-5	5.1455e9

Table 5. Table 5: Mesh name, minimum diameter h m ​ i ​ n subscript ℎ 𝑚 𝑖 𝑛 h_{min} , maximum diameter h m ​ a ​ x subscript ℎ 𝑚 𝑎 𝑥 h_{max} , number of cells, and average number of cells along a throat diameter D t / h subscript 𝐷 𝑡 ℎ D_{t}/h for all the meshes under consideration.

mesh name	$h_{min}$	$h_{max}$	No. of cells	$D_t/h$
$100k$	2.6784e-4	1.1e-3	1.03e5	11
$320k$	1.8186e-4	7.7958e-4	3.21e5	16
$770k$	1.3611e-4	5.9389e-4	7.72e5	21
$1600k$	1.0613e-4	4.6966e-4	1.59e6	27
$2700k$	8.8289e-5	3.947e-4	2.66e6	32
$6400k$	6.529e-5	2.957e-4	6.38e6	43

Table 6. Table 6: Mesh name, minimum diameter h m ​ i ​ n subscript ℎ 𝑚 𝑖 𝑛 h_{min} , maximum diameter h m ​ a ​ x subscript ℎ 𝑚 𝑎 𝑥 h_{max} , and average number of cells along a throat diameter D t / h subscript 𝐷 𝑡 ℎ D_{t}/h for the tetrahedral meshes used in [ 8 ] (left) and [ 54 ] (right).

mesh name	$h_{min}$	$h_{max}$	$D_t/h$	mesh name	$h_{min}$	$h_{max}$	$D_t/h$
$140k$	3.39e-4	3.09e-3	10	$5000k$	3.2e-4	3.6e-4	13
$330k$	2.23e-4	1.93e-3	11	$15000k$	1.9e-4	2.3-4	21
$900k$	1.09e-4	1.87e-3	19	$50000k$	1.2e-4	1.7e-4	34
$1200k$	1.08e-4	1.63e-3	19
$1900k$	1.06e-4	1.49e-3	20
$3000k$	1.17e-4	9.64e-4	21

Equations

$$\displaystyle\rho\,\partial_{t}{\mbox{\boldmath{$u$}}}+\rho\,\nabla\cdot\left({\mbox{\boldmath{$u$}}}\otimes{\mbox{\boldmath{$u$}}}\right)-\nabla\cdot\boldsymbol{\sigma}$$

$$\displaystyle\nabla\cdot{\mbox{\boldmath{$u$}}}$$

$$u$$

$$\displaystyle\boldsymbol{\sigma}\cdot{\mbox{\boldmath{$n$}}}$$

$$\boldsymbol{\sigma}({\mbox{\boldmath{$u$}}},p)=-p\mathbf{I}+\mu(\nabla{\mbox{\boldmath{$u$}}}+\nabla{\mbox{\boldmath{$u$}}}^{T}),$$

$$\displaystyle\rho\,\partial_{t}{\mbox{\boldmath{$u$}}}+\rho\,\nabla\cdot\left({\mbox{\boldmath{$u$}}}\otimes{\mbox{\boldmath{$u$}}}\right)+\nabla p-2\mu\,\Delta{\mbox{\boldmath{$u$}}}={\mbox{\boldmath{$f$}}}\mbox{ in }\Omega\times(t_{0},T).$$

$$Re=\frac{UL}{\nu },$$

$$\eta =R{e}^{-3/4}L.$$

$$\displaystyle\rho\,\partial_{t}{\mbox{\boldmath{$u$}}}+\rho\,\nabla\cdot\left(\overline{{\mbox{\boldmath{$u$}}}}\otimes{\mbox{\boldmath{$u$}}}\right)-2\mu\Delta{\mbox{\boldmath{$u$}}}+\nabla p={\mbox{\boldmath{$f$}}}$$

$$\displaystyle\nabla\cdot{\mbox{\boldmath{$u$}}}=0$$

$$\displaystyle-2\alpha^{2}\nabla\cdot\left(a({\mbox{\boldmath{$u$}}})\nabla\overline{{\mbox{\boldmath{$u$}}}}\right)+\overline{{\mbox{\boldmath{$u$}}}}+\nabla\lambda={\mbox{\boldmath{$u$}}}$$

$$\displaystyle\nabla\cdot\overline{{\mbox{\boldmath{$u$}}}}=0$$

$$\displaystyle a({\mbox{\boldmath{$u$}}})\simeq 0$$

$$\displaystyle a({\mbox{\boldmath{$u$}}})\simeq 1$$

$$\displaystyle\overline{{\mbox{\boldmath{$u$}}}}$$

$$\displaystyle(2\alpha^{2}a({\mbox{\boldmath{$u$}}})\nabla\overline{{\mbox{\boldmath{$u$}}}}-\lambda\mathbf{I}){\mbox{\boldmath{$n$}}}$$

$$\displaystyle\rho\,\frac{3}{2\Delta t}\,{\mbox{\boldmath{$u$}}}^{n+1}+\rho\,\nabla\cdot\left(\overline{{\mbox{\boldmath{$u$}}}}^{n+1}\otimes{\mbox{\boldmath{$u$}}}^{n+1}\right)-2\mu\Delta{\mbox{\boldmath{$u$}}}^{n+1}+\nabla p^{n+1}$$

$$\displaystyle\nabla\cdot{\mbox{\boldmath{$u$}}}^{n+1}$$

$$\displaystyle-\alpha^{2}\nabla\cdot\left(a({\mbox{\boldmath{$u$}}}^{n+1})\nabla\overline{{\mbox{\boldmath{$u$}}}}^{n+1}\right)+\overline{{\mbox{\boldmath{$u$}}}}^{n+1}+\nabla\lambda^{n+1}$$

$$\displaystyle\nabla\cdot\overline{{\mbox{\boldmath{$u$}}}}^{n+1}$$

$$\displaystyle\rho\,\frac{3}{2\Delta t}\,{\mbox{\boldmath{$v$}}}^{n+1}+\rho\,\nabla\cdot\left({\mbox{\boldmath{$u$}}}^{*}\otimes{\mbox{\boldmath{$v$}}}^{n+1}\right)-2\mu\Delta{\mbox{\boldmath{$v$}}}^{n+1}+\nabla q^{n+1}$$

$$\displaystyle\nabla\cdot{\mbox{\boldmath{$v$}}}^{n+1}$$

$$\displaystyle-\alpha^{2}\nabla\cdot\left(a({\mbox{\boldmath{$v$}}}^{n+1})\nabla\overline{{\mbox{\boldmath{$v$}}}}^{n+1}\right)+\overline{{\mbox{\boldmath{$v$}}}}^{n+1}+\nabla\lambda^{n+1}$$

$$\displaystyle\nabla\cdot\overline{{\mbox{\boldmath{$v$}}}}^{n+1}$$

$$\displaystyle{\mbox{\boldmath{$u$}}}^{n+1}$$

$$p^{n+1}$$

$$\displaystyle\frac{\rho}{\Delta t}\overline{{\mbox{\boldmath{$v$}}}}^{n+1}-\nabla\cdot\left(\overline{\mu}\nabla\overline{{\mbox{\boldmath{$v$}}}}^{n+1}\right)+\nabla\overline{q}^{n+1}$$

$$a({\mbox{\boldmath{$v$}}})=a_{D}({\mbox{\boldmath{$v$}}})=\left|{\mbox{\boldmath{$v$}}}-D(F({\mbox{\boldmath{$v$}}}))\right|^{2},$$

$$D_N=\sum _{n=0}^N(I-F{)}^n.$$

$$\displaystyle a_{D_{0}}({\mbox{\boldmath{$v$}}})=\left|{\mbox{\boldmath{$v$}}}-F({\mbox{\boldmath{$v$}}})\right|.$$

$$F=F_H={\left(I-{\alpha }^2\Delta \right)}^{-1}.$$

$$\displaystyle{\mbox{\boldmath{$v$}}}-D_{N}(F_{H}({\mbox{\boldmath{$v$}}}))=(-1)^{N+1}\alpha^{2N+2}\Delta^{N+1}F_{H}^{N+1}{\mbox{\boldmath{$v$}}}.$$

$$\tilde{{\mbox{\boldmath{$v$}}}}^{n+1}-\alpha^{2}\Delta\tilde{{\mbox{\boldmath{$v$}}}}^{n+1}={\mbox{\boldmath{$v$}}}^{n+1}.$$

$$\displaystyle\rho\,\frac{3}{2\Delta t}\,\int_{\Omega_{i}}{\mbox{\boldmath{$v$}}}^{n+1}d\Omega$$

$$\displaystyle+\int_{\Omega_{i}}\nabla q^{n+1}d\Omega=\int_{\Omega_{i}}{\mbox{\boldmath{$b$}}}^{n+1}d\Omega.$$

$$\displaystyle\rho\,\frac{3}{2\Delta t}\,\int_{\Omega_{i}}{\mbox{\boldmath{$v$}}}^{n+1}d\Omega$$

$$\displaystyle+\int_{\partial\Omega_{i}}q^{n+1}d\textbf{A}=\int_{\Omega_{i}}{\mbox{\boldmath{$b$}}}^{n+1}d\Omega.$$

$${\int }_{\partial {\Omega }_i}{q}^{n+1}d\text{A}\approx \sum _j{q}_j^{n+1}{\text{A}}_j,$$

Full text

1 Introduction

The Direct Numerical Simulation (DNS) of the Navier-Stokes Equations (NSE) computes the evolution of all the significant flow structures by resolving them with a properly refined mesh. Unfortunately, when the convection dominates the dynamics, i.e. the Reynolds number (the dimensionless number that weights the importance of inertial forces versus viscous ones) is “large”, this requires very fine meshes, making DNS computationally unaffordable for practical purposes. Therefore, for many realistic problems the simulation of turbulent flows is performed by using different models.

The NSE can be properly averaged (in different ways, quite often in time), leading to the so-called Reynolds-Averaged Navier-Stokes equations (RANS), or filtered (usually in space), leading to Large Eddy Simulation (LES) techniques (see, e.g., [18]). In this work, we focus on the latter approach. We consider a variant of the so-called Leray model [31], where the small-scale effects are described by a set of equations to be added to the discrete NSE formulated on the under-refined mesh. The extra-problem can be devised in different ways, for instance by a functional splitting of the solved and unresolved scales [6] or by resorting to the concept of “suitability” of weak solutions [20]. We consider a variation of the Leray model for which the extra problem acts as a differential low-pass filter [13]. For its actual implementation, we use the Evolve-Filter-Relax (EFR) algorithm proposed in [30]. For a formal rewriting of the EFR algorithm as a LES technique we refer to [36], while we refer to [8] for a reformulation in an operator-splitting framework. One of the advantages of this approach is that it is easily implemented in a legacy Navier-Stokes solver, since the filtering step requires the solution of a Stokes-like problem.

The Leray model has been extensively applied within a Finite Element (FE) framework, although obviously other space approximations are possible. Here, we focus on Finite Volume (FV) methods, which have been widely used in the LES context. So far, the application of the Leray model in a FV framework has been unexplored. In this manuscript, we intend to fill this gap by proposing a computationally efficient FV method for the EFR algorithm. We target applications involving incompressible fluid flow in both two and three dimensions, with the 3D applications featuring moderately large Reynolds numbers (in the range of a few thousands). We will show that that for similar levels of mesh refinement FV and FE methods provide significantly different results. Thus, changing the space discretization method for the EFR algorithm is not a trivial exercise.

We choose two benchmarks to showcase the features of our approach. The first is a classical academic benchmark: 2D flow past a cylinder [26, 47] at Reynolds number 100. Although this flow field is not turbulent, it turns out to be very challenging for most numerical models and methods, especially on coarser computational grids. This first test serves the purpose of getting a preliminary understanding of how the EFR algorithm works in combination with a FV method, in terms of element type, level of mesh refinement, and choice of critical parameter values. The second benchmark has been issued by the U.S. Food and Drug Administration (FDA) and it involves 3D flow at Reynolds numbers up to 6500 through a nozzle. Three independent laboratories were requested by FDA to perform flow visualization experiments on fabricated nozzle models for different flow rates [21]. This resulted in benchmark data available online to the scientific community for the validation of Computational Fluid Dynamics (CFD) simulations [1]. Available experimental measurements enable us to check the effectiveness of our FV-based EFR algorithm in simulating average macroscopic quantities. Through our numerical experiments, we are able to provide practical directions to tune the parameters involved in the model. Moreover, we are able to investigate the impact of mesh features (element type, non-orthogonality, local refinement, and element aspect ratio) and the discretization method for the convective term on the agreement between numerical solutions and experimental data.

All the computational results presented in this article have been performed with OpenFOAM® [52], an open source finite volume C++ library widely used by commercial and academic organizations. See, e.g., [19, 32] for numerical results obtained with LES techniques implemented in OpenFOAM®. An important outcome of this work is that the code created for it is incorporated in an open-source library111https://mathlab.sissa.it/cse-software and therefore is readily shared with the community.

This work is outlined as follows. In Sec. 2, we introduce the continuous Leray model as well as the numerical approximation proposed in [30]. In Sec. 3, we detail our strategy for space discretization, which combines an operator splitting and a Finite Volume method. Additionally, we show how to tune certain model parameters using physical and discretization quantities. The numerical results for the flow around the cylinder are reported in Sec. 4.1, while the comparison between numerical results and experimental data provided by the FDA is presented in Sec. 4.2. Conclusions are drawn in Sec. 5.

2 Problem definition

2.1 The Navier-Stokes equations

We consider the motion of an incompressible viscous fluid in a time-independent domain $\Omega$ over a time interval of interest $(t_{0},T)$. The flow is described by the incompressible Navier-Stokes equations

[TABLE]

endowed with the boundary conditions

[TABLE]

and the initial data {\mbox{\boldmath{u}}}={\mbox{\boldmath{u}}}_{0} in $\Omega\times\{t_{0}\}$. Here $\overline{\partial\Omega_{D}}\cup\overline{\partial\Omega_{N}}=\overline{\partial\Omega}$ and $\partial\Omega_{D}\cap\partial\Omega_{N}=\emptyset$. In addition $\rho$ is the fluid density, $u$ is the fluid velocity, $\partial_{t}$ denotes the time derivative, $\boldsymbol{\sigma}$ is the Cauchy stress tensor, $f$ accounts for possible body forces (such as, e.g., gravity), {\mbox{\boldmath{u}}}_{D},{\mbox{\boldmath{g}}} and {\mbox{\boldmath{u}}}_{0} are given. Equation (1) represents the conservation of the linear momentum, while eq. (2) represents the conservation of the mass. For Newtonian fluids $\boldsymbol{\sigma}$ can be written as

[TABLE]

where $p$ is the pressure and $\mu$ is the constant dynamic viscosity. Notice that by plugging (5) into eq. (1), eq. (1) can be rewritten as

[TABLE]

In order to characterize the flow regime under consideration, we define the Reynolds number as

[TABLE]

where $\nu=\mu/\rho$ is the kinematic viscosity of the fluid, and $U$ and $L$ are characteristic macroscopic velocity and length, respectively. For an internal flow in a cylindrical pipe, $U$ is the mean sectional velocity and $L$ is the diameter. For moderately large Reynolds numbers the effects of flow disturbances cannot be neglected, and yet Reynolds-averaged Navier-Stokes (RANS) models [39] are generally too crude.

2.2 Leray model

In the framework of the Kolmogorov 1941 theory [28, 27], the turbulent kinetic energy, which is the kinetic energy associated with eddies in the turbulent flow, is injected in the system at the large scales (low wave numbers). Since the large scale eddies are unstable, they break down, transferring the energy to smaller eddies. Finally, the turbulent kinetic energy is dissipated by the viscous forces at the small scales (high wave numbers). This process is usually referred to as energy cascade. The scale at which the viscous forces dissipate energy is referred to as Kolmogorov scale. For a flow in developed turbulent regime and at statistical equilibrium, the Kolmogorov scale can be expressed as

[TABLE]

In order to correctly capture the dissipated energy, DNS needs a mesh with spacing $h\sim\eta$. As the Reynolds number increases, DNS leads to a huge number of unknowns and prohibitive computational costs. On the other hand, when the mesh size $h$ fails to resolve the Kolmogorov scale, the under-diffusion in the simulation leads to nonphysical computed velocities. In some cases, this is detectable simply looking at the velocity field, which features nonphysical oscillations eventually leading to a simulation break down. However, in some cases the velocity field does not display oscillations, yet it does not correspond to the physical solution. A possible remedy to this issue is to introduce a model which filters the nonphysical oscillations in the velocity field and conveys the energy lost to resolved scales.

The so called Leray model couples the Navier-Stokes equations (6),(2) with a differential filter. The resulting system reads

[TABLE]

Here, \overline{{\mbox{\boldmath{u}}}} is the filtered velocity, $\alpha$ can be interpreted as the filtering radius (that is, the radius of the neighborhood where the filter extracts information from the unresolved scales), the variable $\lambda$ is a Lagrange multiplier to enforce the incompressibility constraint for \overline{{\mbox{\boldmath{u}}}} and $a(\cdot)$ is a scalar function such that:

[TABLE]

This function, called indicator function, is crucial for the success of the Leray model. In Sec. 2.4, we will discuss our choice of a({\mbox{\boldmath{u}}}). Here, we mention that the choice a({\mbox{\boldmath{u}}})\equiv 1 corresponds to the classic Leray-$\alpha$ model [31]. This model has the advantage of making the operator in the filter equations linear and constant in time, but its effectivity is rather limited, since it introduces the same amount of regularization everywhere in the domain, hence causing overdiffusion.

Equations (11)-(12) require suitable boundary conditions. These are chosen to be

[TABLE]

Even though (9)-(10) are linear in ({\mbox{\boldmath{u}}},p) and the filter problem is linear in (\overline{{\mbox{\boldmath{u}}}},\lambda), the coupling is non-linear, due to the term \nabla\cdot\left(\overline{{\mbox{\boldmath{u}}}}\otimes{\mbox{\boldmath{u}}}\right) in eq. (9) and the term a({\mbox{\boldmath{u}}})\nabla\overline{{\mbox{\boldmath{u}}}} in eq. (11) (when $a(\cdot)$ is not constant).

2.3 Time discrete problem

To discretize in time problem (9)-(12), let $\Delta t\in\mathbb{R}$, $t^{n}=t_{0}+n\Delta t$, with $n=0,...,N_{T}$ and $T=t_{0}+N_{T}\Delta t$. Moreover, we denote by $y^{n}$ the approximation of a generic quantity $y$ at the time $t^{n}$.

For the time discretization of system (9)-(12), we adopt a Backward Differentiation Formula of order 2 (BDF2), see e.g. [40]. The Leray system discretized in time reads: given {\mbox{\boldmath{u}}}^{0}, for $n\geq 0$ find the solution ({\mbox{\boldmath{u}}}^{n+1},p^{n+1},\overline{{\mbox{\boldmath{u}}}}^{n+1},\lambda^{n+1}) of system:

[TABLE]

where {\mbox{\boldmath{b}}}^{n+1}={\mbox{\boldmath{f}}}^{n+1}+(4{\mbox{\boldmath{u}}}^{n}-{\mbox{\boldmath{u}}}^{n-1})/(2\Delta t). Obviously, other discretization schemes are possible. However, for clarity of exposition we will restrict the description of the approach to the case of BDF2.

A monolithic approach for problem (15)-(18) would lead to high computational costs, making the advantage compared to DNS questionable. To decouple the Navier-Stokes system (15)-(16) from the filter system (17)-(18) at the time $t^{n+1}$, we have two options:

Filter-then-solve: Solve the filter equations (17)-(18) first, with {\mbox{\boldmath{u}}}^{n+1} replaced by a suitable extrapolation {\mbox{\boldmath{u}}}^{*} and a({\mbox{\boldmath{u}}}^{n+1}) replaced by a({\mbox{\boldmath{u}}}^{*}), and then solve equations (15)-(16) with advection field given by the filtered velocity previously computed. 2. 2.

Solve-then-filter: Solve equations (15)-(16) first, replacing the advection field \overline{{\mbox{\boldmath{u}}}}^{n+1} with a suitable extrapolation {\mbox{\boldmath{u}}}^{*}, and then solve the filter problem (17)-(18).

In this work, we will focus on approach 2. To keep the computational costs low, we adopt a semi-implicit approach, i.e. we perform only one iteration per time step. In particular, we consider a modified version of approach 2 called Evolve-Filter-Relax (EFR), which was first proposed in [30]. This algorithm reads as follows: given the velocities {\mbox{\boldmath{u}}}^{n-1} and {\mbox{\boldmath{u}}}^{n}, at $t^{n+1}$:

1. i)

evolve: find intermediate velocity and pressure ({\mbox{\boldmath{v}}}^{n+1},q^{n+1}) such that

[TABLE]

where {\mbox{\boldmath{u}}}^{*}=2{\mbox{\boldmath{u}}}^{n}-{\mbox{\boldmath{u}}}^{n-1}. 2. ii)

filter: find (\overline{{\mbox{\boldmath{v}}}}^{n+1},\lambda^{n+1}) such that

[TABLE] 3. iii)

relax: set

[TABLE]

where $\chi\in(0,1]$ is a relaxation parameter.

Remark 2.1**.**

Filter problem (21)-(22) can be considered a generalized Stokes problem. In fact, by dividing eq. (21) by $\Delta t$ and rearranging the terms we obtain:

[TABLE]

*where $\overline{q}^{n+1}=\rho\lambda^{n+1}/\Delta t$. Problem (25),(22) can be seen as a time dependent Stokes problem with a non-constant viscosity $\overline{\mu}$, discretized by the Backward Euler (or BDF1) scheme. A solver for problem (25),(22) can then be obtained by adapting a standard linearized Navier-Stokes solver. Notice that for $\alpha^{2}/\Delta t\simeq h^{2}/\Delta t\to 0$ we have that $\overline{\mu}\rightarrow 0$ and \overline{{\mbox{\boldmath{v}}}}_{h}^{n+1}\to{\mbox{\boldmath{v}}}_{h}^{n+1}. *

Remark 2.2**.**

In the EFR algorithm proposed in [30] there is no relaxation for the pressure, i.e. the end-of-step pressure is set equal to the pressure of the Evolve step. In [8], two relaxations for the pressure were considered: $p^{n+1}=q^{n+1}+\frac{3}{2}\chi\overline{q}^{n+1}$ or $p^{n+1}=(1-\chi)q^{n+1}+\chi\overline{q}^{n+1}$. Notice while $\overline{q}^{n+1}$ has the same dimensional units as $q^{n+1}$, $\lambda^{n+1}$ does not.

In the rest of the paper, we will call EFR algorithm (19), (20), (25), (22)-(24). We will call the algorithm simply Evolve-Filter (EF) when $\chi=1$, since there is no actual relaxation step. In the next subsection, we will consider an indicator function that leads to a nonlinear alternative to the Leray-$\alpha$ model.

2.4 The indicator function

Different choices of $a(\cdot)$ have been proposed and compared in [11, 30, 22, 49, 12]. Here, we focus on a class of deconvolution-based indicator functions:

[TABLE]

where $F$ is a linear filter (an invertible, self-adjoint, compact operator from a Hilbert space to itself) and $D$ is a bounded regularized approximation of $F^{-1}$. A popular choice for $D$ is the Van Cittert deconvolution operator $D_{N}$, defined as

[TABLE]

The evaluation of $a_{D}$ with $D=D_{N}$ (deconvolution of order $N$) requires then to apply the filter $F$ a total of $N+1$ times. Since $F^{-1}$ is not bounded, in practice $N$ is chosen to be small, as the result of a trade-off between accuracy (for a regular solution) and filtering (for a non-regular one). In this paper we consider $N=0$, corresponding to $D_{0}=I$. For this choice of $N$, the indicator function (26) becomes

[TABLE]

We select $F$ to be the linear Helmholtz filter operator $F_{H}$ defined by

[TABLE]

It is possible to prove [16] that

[TABLE]

Therefore, a_{D_{N}}({\mbox{\boldmath{v}}}) is close to zero in the regions of the domain where $v$ is smooth. We remark that finding F_{H}({\mbox{\boldmath{v}}}^{n+1})=\tilde{{\mbox{\boldmath{v}}}}^{n+1} is equivalent to finding \tilde{{\mbox{\boldmath{v}}}}^{n+1} such that:

[TABLE]

We will refer to eq. (9)-(12) with indicator function given by (27) as Leray-$\alpha$-NL model.

3 Space discrete problem: the Finite Volume approximation

In this section we discuss the space discretization of problems (19)-(20) and (25),(22) for both the classic Leray-$\alpha$ and the Leray-$\alpha$-NL model.

We adopt the Finite Volume (FV) approximation that is derived directly from the integral form of the governing equations. We have chosen to implement the EFR algorithm within the finite volume C++ library OpenFOAM® [52]. We partition the computational domain $\Omega$ into cells or control volumes $\Omega_{i}$, with $i=1,\dots,N_{c}$, where $N_{c}$ is the total number of cells in the mesh. Let Aj be the surface vector of each face of the control volume, with $j=1,\dots,M$. The value of $M$ depends on the dimension of the domain and the type of mesh (hexahedral vs prismatic).

3.1 Numerical discretization for the evolve step of the problem

The integral form of eq. (19) for each volume $\Omega_{i}$ is given by:

[TABLE]

By applying the Gauss-divergence theorem, eq. (3.1) becomes:

[TABLE]

Each term in eq. (3.1) is approximated as follows:

Gradient term:

[TABLE]

where $q_{j}$ is the value of the pressure relative to centroid of the $j^{\text{th}}$ face. In OpenFOAM® solvers, the face center pressure values $q_{j}$ are typically obtained from the cell center values by means of a linear interpolation scheme. Such scheme is rigorously second-order accurate only on structured meshes [46].

• -

Convective term:

[TABLE]

where {\mbox{\boldmath{u}}}^{*}_{j} and {\mbox{\boldmath{v}}}^{n+1}_{j} are respectively the extrapolated convective velocity and the fluid velocity relative to the centroid of each control volume face. In (32), $\varphi^{*}_{j}$ is the convective flux associated to {\mbox{\boldmath{u}}}^{*} through face $j$ of the control volume. In OpenFOAM® solvers, the convective flux at the cell faces is typically a linear interpolation of the values from the adjacent cells. Also {\mbox{\boldmath{v}}}^{n+1} needs to be approximated at cell face $j$ in order to get the face value {\mbox{\boldmath{v}}}^{n+1}_{j}. Different interpolation methods can be applied: central, upwind, second order upwind and blended differencing schemes [25].

• -

Diffusion term:

[TABLE]

where (\nabla{\mbox{\boldmath{v}}}^{n+1})_{j} is the gradient of {\mbox{\boldmath{v}}}^{n+1} at face $j$. We are going to briefly explain how (\nabla{\mbox{\boldmath{v}}}^{n+1})_{j} is approximated with second order accuracy on a structured, orthogonal mesh. Let $P$ and $Q$ be two neighboring control volumes. The term (\nabla{\mbox{\boldmath{v}}}^{n+1})_{j} is evaluated by subtracting the value of velocity at the cell centroid on the $P$-side of the face, denoted with {\mbox{\boldmath{v}}}^{n+1}_{P}, from the value of velocity at the centroid on the $Q$-side, denoted with {\mbox{\boldmath{v}}}^{n+1}_{Q}, and dividing by the magnitude of the distance vector $\textbf{d}_{j}$ connecting the two cell centroids:

[TABLE]

For non-structured, non-orthogonal meshes (see Fig. 1), an explicit non-orthogonal correction has to be added to the orthogonal component in order to preserve second order accuracy. See [25] for details.

Let us denote with {\mbox{\boldmath{v}}}^{n+1}_{i} and {\mbox{\boldmath{b}}}^{n+1}_{i} the average velocity and source term in control volume $\Omega_{i}$, respectively. Moreover, we denote with {\mbox{\boldmath{v}}}^{n+1}_{i,j} and $q^{n+1}_{i,j}$ the velocity and pressure associated to the centroid of face $j$ normalized by the volume of $\Omega_{i}$. Then the discretized form of eq. (3.1), divided by the control volume $\Omega_{i}$, can be written as:

[TABLE]

Following [42], we now write eq. (33) in semi-discretized form, i.e. with the pressure term in continuous form while all the other terms are in discrete form

[TABLE]

with

[TABLE]

Next, we take the divergence of the eq. (34) and make use of eq. (20) to obtain

[TABLE]

By integrating eq. (35) over the control volume $\Omega_{i}$ and applying the Gauss-divergence theorem, we get

[TABLE]

So, the space discretized eq. (35), divided by the control volume $\Omega_{i}$, can be expressed as:

[TABLE]

In eq. (36), $(\nabla q^{n+1})_{j}$ is the gradient of $q^{n+1}$ at faces $j$ and it is approximated in the same way as (\nabla{\mbox{\boldmath{v}}}^{n+1})_{j}.

Finally, the fully discretized form of problem (19)-(20) is given by system (33),(36). In this work, we choose a partitioned approach to deal with the pressure-velocity coupling. The partitioned algorithms available in OpenFOAM® are SIMPLE [38] for steady-state problems, and PISO [23] and PIMPLE [33] for transient problems. For the results reported in Sec. 4, we used the PISO algorithm. The splitting of operations in the solution of the discretised momentum and pressure equations gives rise to a formal order of accuracy of the order of powers of $\Delta t$ depending on the number of operation-splittings used (see [23] for more details). We remark that in the OpenFOAM® implementation of the PISO solver the mass flux is modified through an additional term than can cause artificial dissipation and is not a part of the original PISO algorithm in [23]. See [50] for more details.

3.2 Numerical discretization for the filter step of the problem

In this subsection, we present the discretization of the filter problem (25),(22). The integral form of the eq. (25) for each volume $\Omega_{i}$ is given by:

[TABLE]

By making use of the Gauss-divergence theorem, eq. (37) becomes:

[TABLE]

Like in Sec. 3.1, we proceed with providing the approximation of each term in eq. (38):

Gradient term:

[TABLE]

where $\overline{q}_{j}^{n+1}$ is the value of the auxiliary pressure associated to the centroid of the $j^{\text{th}}$ face.

• -

Diffusion term:

[TABLE]

where (\nabla\overline{{\mbox{\boldmath{v}}}}^{n+1})_{j} is the gradient of \overline{{\mbox{\boldmath{v}}}}^{n+1} at face $j$. It is approximated in the same way as (\nabla{\mbox{\boldmath{v}}}^{n+1})_{j}; see Sec. 3.1.

Upon division by the volume of $\Omega_{i}$, the discretized form of the eq. (38) can be written as

[TABLE]

In eq. (41), we denoted with \overline{{\mbox{\boldmath{v}}}}^{n+1}_{i} the average filtered velocity in control volume $\Omega_{i}$, while $\overline{q}^{n+1}_{i,j}$ is the auxiliary pressure at the centroid of face $j$ normalized by the volume of $\Omega_{i}$.

Next, we rewrite eq. (41) in semi-discretized form, i.e. with the pressure term in continuous form while all the other terms are in discrete form, take its divergence and use eq. (22) to get:

[TABLE]

By integrating eq. (42) over the control volume $\Omega_{i}$ and by applying the Gauss-divergence theorem, we obtain

[TABLE]

Finally, we divide by the volume of $\Omega_{i}$ to get

[TABLE]

The fully discrete problem associated to the filter is given by (41), (44). Also for this problem, we choose a partitioned algorithm: a slightly modified version of the SIMPLE algorithm, called SIMPLEC algorithm [48], that features improved and accelerated convergence towards a steady state solution. We found that the SIMPLEC algorithm is a necessary choice for the filter problem in the Leray-$\alpha$ model, because the standard SIMPLE algorithm does not converge. On the other hand, the SIMPLE algorithm does converge for the filter problem in the Leray-$\alpha$-NL model.

For the Leray-$\alpha$-NL model, we also need to approximate the solution to problem (28). After using Gauss-divergence theorem, the integral form of eq. (28) reads:

[TABLE]

We approximate the diffusion term as:

[TABLE]

where (\nabla\tilde{{\mbox{\boldmath{v}}}}^{n+1})_{j} is treated in the same way as (\nabla{\mbox{\boldmath{v}}}^{n+1})_{j} and (\nabla\overline{{\mbox{\boldmath{v}}}}^{n+1})_{j}. Once we divide by the volume of $\Omega_{i}$, the space discretized form of eq. (45) is

[TABLE]

where \tilde{{\mbox{\boldmath{v}}}}^{n+1}_{i} is the average value of \tilde{{\mbox{\boldmath{v}}}}^{n+1} in control volume $\Omega_{i}$.

3.3 Setting of the relaxation parameter $\chi$

In [8], a heuristic formula for estimating the value of $\chi$ is provided. The proposed idea is to set $\chi$ such that the total viscous stress on an under-refined mesh of size $h$ provides the same amount of dissipation of the viscous stress on a fully refined mesh of size $\eta$. This leads to:

[TABLE]

We propose here an alternative approach. We set $\chi$ such that the viscous contribution on a fully refined mesh of size $\eta$ is equivalent to the convex linear combination of the viscous contribution on an under-resolved mesh of size $h$ and the viscous contribution due to the artificial viscosity introduced by the Leray model:

[TABLE]

Eq. (48) can be interpreted as a relaxation step for the velocity gradient. With the approximation $\nabla_{\xi}\approx\xi^{-1}$, we obtain

[TABLE]

In practice, we may set $||a||_{\infty}=1$, which is a good approximation on coarse meshes, and obtain:

[TABLE]

We observe that both estimates of $\chi$ defined in eq. (49) are clearly positive. In fact, they would be negative when $h<\eta$, which means that the mesh is refined enough for a DNS and the filtering step is not needed.

4 Numerical results

In this section, we present several numerical results for the NSE (i.e. no turbulence model), Leray-$\alpha$ and Leray-$\alpha$ NL models. Two test cases are investigated: 2D flow past a cylinder and a 3D benchmark from the FDA.

The number of PISO loops and non-orthogonal correctors has been fixed to 2 for all the simulations. The following solvers have provided a good compromise between stability, accuracy, and numerical cost. The linear algebraic system associated with eq. (33) is solved using an iterative solver with symmetric Gauss-Seidel smoother. For eq. (41) and (46), we use the Diagonal Incomplete Cholesky Preconditioned Conjugate Gradient. Finally, for Poisson problems (36) and (44) we use Geometric Agglomerated Algebraic Multigrid Solver GAMG with the Gauss-Seidel smoother. The required accuracy is 1e-6 at each time step.

4.1 2D benchmark: channel flow around a cylinder

The first test we consider is a well-known benchmark [26, 47]. The computational domain is a 2.2 $\times$ 0.41 rectangular channel with a cylinder of radius 0.05 centered at (0.2, 0.2), when taking the bottom left corner of the channel as the origin of the axes. Fig. 2 shows part of the computational domain. We impose a no slip boundary condition on the upper and lower wall and on the cylinder. At the inflow and the outflow we prescribe the following velocity profile:

[TABLE]

The partitioned algorithms we use (see Sec. 3.1 and 3.2) require a boundary condition for the pressure too. We choose {\partial q}/{\partial{\mbox{\boldmath{n}}}}={\partial\overline{q}}/{\partial{\mbox{\boldmath{n}}}}=0 on each boundary where $n$ is the outward normal. We set density $\rho=1$ and viscosity $\mu=10^{-3}$. There is no external force, so {\mbox{\boldmath{f}}}=\mathbf{0}. We start the simulations from fluid at rest. Note that the Reynolds number is time dependent, with $0\leq Re\leq 100$ [47]. The quantities of interest for this benchmark are the drag and lift coefficients:

[TABLE]

where $U_{r}=1$ is the maximum velocity at the inlet/outlet, $L_{r}=0.1$ is the diameter of the cylinder, $S$ is the surface of the cylinder, and $\boldsymbol{t}$ and $\boldsymbol{n}$ are the tangential and normal unit vectors, respectively.

We consider several meshes with prismatic and hexahedral elements generated using Gmsh [3]. Table 1 reports mesh name, minimum and maximum diameter, and number of cells for each mesh. The name of each grid refers to the number of cells and the subscript denotes the kind of element, prismatic ($P$) or hexahedral ($H$). See Fig. 2 (a) shows part mesh $16k_{H}$. The quality of all the meshes under consideration is high. Hexahedral and prismatic meshes have very low values of maximum non-orthogonality (33∘ to 39∘), average non-orthogonality (4∘ to 7∘), skewness (0.5 to 1.2), and maximum aspect ratio (up to 2). Time step used for the simulations was governed by the maximum Courant-Friedrichs-Lewy number ($CFL_{max}$), set to $CFL_{max}=0.2$. For all the simulations, we use the PISO algorithm with flux correction term [50] and second-order accurate Linear Upwind Differencing (LUD) [51] scheme for the convective term.

Remark 4.1**.**

As discussed in Sec. 2.2, DNS needs a mesh with spacing $h\sim\eta$. For the current test, the Reynolds number is time dependent. One could either calculate the Kolmogorov scale (8) based on the maximum Reynolds number, which gives $\eta\approx 3.2e-3$, or introduce a mean Kolmogorov scale $\bar{\eta}$:

[TABLE]

*where $T=8$ is the final computational time. Notice that $\eta\approx\bar{\eta}$, thus a DNS would roughly need a number of cells of the order of 100k. *

Fig. 3 compares the NSE solution obtained with meshes $200k_{P}$ and $200k_{H}$ with the reference curves reported in [26]. We see that the solutions are in good agreement, so a DNS is possible both with meshes $200k_{P}$ and $200k_{H}$. Notice that the level of refinement of these meshes is in agreement with that one predicted in Remark 4.1. Hereinafter, we will refer to the solutions computed with a DNS on meshes $200k_{P}$ and $200k_{H}$ as the true solution and we will refer to meshes $200k_{P}$ and $200k_{H}$ as DNS meshes.

Fig. 4 displays the lift and drag coefficients computed from the NSE solution on all the meshes in Table 1 that are coarser than the DNS meshes. We observe that the lift coefficients computed from the NSE solutions on the coarser meshes show a large difference with respect to the true $c_{l}$, in terms of both amplitude and phase. To quantitate this difference, Table 2 reports maximum lift and drag coefficients obtained from the NSE solutions and times at which the maxima occur on all the meshes, together with the corresponding values from [26]. We observe a monotonic convergence for the maximum values of $c_{l}$ and $c_{d}$ computed on hexahedral meshes, while for prismatic meshes the convergence is not monotonic. In Table 2, our results computed with hexahedral (resp., prismatic) meshes are compared with the results in [26] computed with $Q_{2}$-$P_{1}^{d}$ (resp., $P_{2}$-$P_{1}$) finite elements and the Crank-Nicolson scheme for time discretization. Notice that the two sets of value from [26] reported in Table 2 coincide up to the reported number of digits.

For all simulations, we evaluated also the error for the maximum values of aerodynamic coefficients compared to the true values:

[TABLE]

Fig. 5 confirms what already noted from Table 2: for hexahedral meshes the errors (53) decrease monotonically as the mesh is refined, while that is not the case for prismatic meshes. In particular, we notice that the $c_{l,max}$ computed with mesh $16k_{P}$ is surprisingly good, while it worsens for meshes $25k_{P}$ and $63k_{P}$. It is only with mesh $120k_{P}$ that we obtain a better $c_{l,max}$ than the one computed with mesh $16k_{P}$.

Now that we have identified the right level of refinement for a DNS (i.e.,  roughly $200k$ cells), we test the EFR algorithm on all meshes in Table 1 that have a lower level of refinement. We aim at comparing the results given by the Leray-$\alpha$ and Leray-$\alpha$ NL models, and understanding the role of $\chi$ and $\alpha$ for both models. We have chosen to set $\chi=0$, i.e. no relaxation, and $\chi=\Delta t$, as suggested in [30] to keep the numerical dissipation low. When $\chi$ = 0, we will refer to the algorithm as Evolve-Filter (EF), and in particular EF L (resp., EF NL) if the Leray-$\alpha$ (resp., Leray-$\alpha$ NL) model is used. As for $\alpha$, two characteristic values are considered: the length $h_{min}$ of the shortest edge in the mesh (see [10] and Remark 4.3 in [8]) and the Kolmogorov scale $\eta$ based on the maximum Reynolds number (see Remark 4.1). Notice that the following “coarse” meshes have $h_{min}<\eta$: $63k_{H}$, $120k_{P}$, $120k_{H}$. Moreover, meshes $25k_{H}$ and $63k_{P}$ have $h_{min}$ which is comparable to $\eta$. Thus, the choice $\alpha=\eta$ is reasonable. For all the different cases under consideration we will report $c_{l}$, which is the more critical coefficient, while for $c_{d}$ we will only report the maximum value in Table 3.

Fig. 6 and 7 show the evolution of $c_{l}$ over time computed by the EF L and EF NL algorithms (i.e, $\chi=0$) on all the meshes coarser than the DNS meshes with hexahedral and prismatic elements, respectively. Fig. 6 and 7 report also the NSE results, i.e. the same results shown in Fig. 4, which were obtained with no filtering step. From Fig. 6 and 7 we see that the $c_{l}$ computed on a given mesh with both EF L and EF NL algorithms is very far from the true $c_{l}$, and worse than then $c_{l}$ given by NSE regardless of the value of $\alpha$. In fact, as expected, the EF algorithm reduces the peaks due to excessive artificial diffusion. Moreover, the linear filter completely dampens the high frequency modes for either choice of $\alpha$ and for both kind of meshes. Fig. 6 shows that also when using the EF NL algorithm on the hexahedral meshes the computed results obtained get closer and closer to the true solution as the mesh gets refined and for $\alpha=\text{min}\{h_{min},\eta\}$. This improvement of the computed $c_{l}$ is lost on the prismatic meshes. Focusing on the hexahedral meshes, we see that the choice of $\alpha$ makes a visible difference in the computed $c_{l}$ on meshes finer that $25k$. In conclusion, we learned that the NSE algorithm is to be preferred to the EF algorithms and hexahedral meshes are superior to prismatic meshes.

Since the EF L algorithm systematically introduces excessive numerical diffusion, we test the only the nonlinear version of the EFR algorithm, denoted hereafter simply by EFR. We set $\chi=\Delta t$. Fig. 8 shows the evolution of $c_{l}$ over time computed by the EFR algorithm on all the meshes coarser than the DNS meshes. Fig. 8 reports also the NSE results. We notice that the EFR algorithm gives a lift coefficient very close to the true solution on all the meshes and for all values of $\alpha$. We also observe that the filter helps reduce the phase advancement and the choice of $\alpha$ makes little-to-no visible difference. Moreover, it is confirmed that the results obtained with hexahedral meshes converge monotonically to the true solution with increase mesh refinement.

Finally, in Table 3 we report a quantitative comparison of all the simulations in this section in terms of maximum lift and drag coefficients and times at which the maxima occur. Notice that for all the meshes the EFR algorithm with either $\alpha=h_{min}$ or $\alpha=\eta$ is the choice that minimizes $|E_{c_{l}}|+|E_{c_{d}}|$.

4.2 3D benchmark: FDA

The second test we consider is a benchmark from the U.S. Food and Drug Administration (FDA), issued within the “Critical Path Initiative” program [2]. The objective of this benchmark is to simulate the flow of an incompressible and Newtonian fluid in nozzle geometry for different flow regimes, from laminar to fully turbulent. Despite its relative simplicity of the geometry, the nozzle contains all the features commonly encountered in medical devices, i.e. flow contraction and expansion, recirculation zones etc. See Fig. 9. Three independent laboratories were requested by FDA to perform flow visualization experiments on fabricated nozzle models for different flow rates up to Reynolds number 6500 [21]. This resulted in benchmark data available online to the scientific community for the validation of CFD simulations [1]. Available experimental measurements enable us to check the effectiveness of the EFR algorithm in simulating average macroscopic quantities. The results of the published inter-laboratory experiments refer to values of the Reynolds numbers in the throat (defined as in (7)) of ${Re_{t}=500,2000,3500,5000,6500}$, $Re_{t}=2000$ being the critical Reynolds number for transitional flow in a straight pipe [41]. To test our methodology, we focus on Reynolds numbers $Re_{t}=2000,3500,5000,6500$.

In Table 4, we report the throat Reynolds number $Re_{t}$, the corresponding inlet Reynolds number $Re_{i}$, the flow rate $Q$, the Kolmogorov length scale $\eta$ and the number of cells required by DNS for the considered flow regimes. The value of $\eta$ was found by plugging into (8) the value of $Re_{t}$ and the diameter of the expansion channel $D_{i}$ as characteristic length. As for the 2D benchmark, the number of cells for the DNS was calculated by assuming a mean grid width equal to $\eta$ as spatial resolution. We notice that the number of cells reported in Table 4 assumes a uniform mesh that has the highest level of refinement possible, since $\eta$ is calculated using the highest Reynolds number in the domain. A mesh with a local refinement dictated by the local Reynolds number would have a more manageable number of cells.

The fluid used in the experiments has $\rho=1056$ Kg/m3 and $\mu=0.0035$ Pa/s. Following [8], the length of the inlet chamber $L_{i}$ was set to four times its diameter, while the length of the expansion channel ($L_{o}$ in Fig. 9) was set to 12 times its diameter. On the lateral surface of the computational domain we prescribe a no-slip boundary condition. At the inflow we prescribe a Poiseuille velocity profile:

[TABLE]

where $r$, $\theta$ and $z$ are the radial, polar and axial locations respectively, $v_{mean}$ is the mean inlet velocity magnitude to obtain the desired flow rate, and $R_{i}=D_{i}/2$. At the outlet section, we prescribe an advective outflow condition:

[TABLE]

where $a$ is the advection velocity magnitude computed so that the total mass is conserved. Boundary condition (55) preserves the numerical stability in case the jet starts to break down too close to the outlet section [54]. The actual experimental set up of the FDA benchmark is a closed loop [21] and it is not reflected by our condition. However, this is expected to introduce a minimal error localized only in the close neighborhood of the outlet section. The results of the flow analysis are not affected, provided that the expansion channel is long enough. We recall that the partitioned algorithms we use require a boundary condition for the pressure as well: we impose {\partial q}/{\partial{\mbox{\boldmath{n}}}}={\partial\overline{q}}/{\partial{\mbox{\boldmath{n}}}}=0 at the inlet and the wall and $p=\overline{q}=0$ at the outlet.

We start all the simulations from fluid at rest and use a smooth transition of the inlet velocity profile to regime flow conditions. We let the simulations run till $T=0.8$ s and then we compare with experimental data the solution obtained by averaging the snapshots collected between $t=0.4$ s and $t=0.8$ s [37]. The comparison between computational and experimental data is made in terms of normalized axial component of the velocity and normalized pressure difference along the centerline. The axial component of the velocity $u_{z}$ is normalized with respect to the average axial velocity at the inlet $\bar{u}_{i}$:

[TABLE]

The pressure difference is normalized with respect to the dynamic pressure in the throat:

[TABLE]

where $p_{z}$ denotes the wall pressure along the $z$-axis and $p_{z=0}$ is the wall pressure at $z=0$, i.e. at the sudden expansion.

In the following, we compare our computational results with the experimental data from [21, 1]. Moreover, we compare our findings with [8] and [54]. In [8], the same LES method is used but it is combined with a Finite Element method ($\mathbb{P}_{2}$-$\mathbb{P}_{1}$ elements) implemented in LifeV [4]. In [54] the authors use a different LES method, called $\sigma$-model [35], and combine it with a central fourth-order Finite Volume scheme implemented in YALES2BIO [5]. In [8] BDF2 is used for the time discretization, while an explicit Runge-Kutta scheme of the fourth order (RK4) and an advanced fourth-order accurate time scheme (TFV4A) are used in [54].

Based on the results in Sec. 4.1, we consider only meshes with hexahedral elements. For the mesh generation we used the mesh generation utility, blockMesh within OpenFOAM. Table 5 reports the mesh name, minimum and maximum diameter, number of cells, and average number of cells along a throat diameter for all the meshes under consideration. Just like in Sec. 4.1, the name of each grid refers to the number of cells (Figure  9). We generated the meshes to have a similar refinement to the tetrahedral meshes used in [8], while the tetrahedral meshes considered in [54] have a much larger number of cells. See Table 6. Also for this benchmark, we made sure that all the meshes feature high quality: low values of maximum non-orthogonality (20∘ to 25∘), average non-orthogonality (around 4∘), skewness (around 1), and maximum aspect ratio (up to 7). Meshes were selectively refined in the convergent, throat, and sudden expansion.

For the sake of completeness, we performed additional tests with the tetrahedral grids used in [8] (see Table 6 left) that worked well within a finite element framework. We found that the high non-orthogonality of those meshes (up to 55∘) significantly affects the computational time and the accuracy of the computed solutions. Thus, we decided to disregard those results and focus on low non-orthogonality hexahedral meshes that generally represent the best choice in the finite volume context [25]. However, we observe that tetrahedral grids with very good features in terms of skewness and aspect ratio were adopted for finite volume simulations in [54]. Moreover, the authors state that additional tests on hexahedral meshes (not reported in the manuscript) did not show substantial differences, confirming that the quality of their tetrahedral meshes is high.

We start to investigate the cases of $Re_{t}=3500$ and above. We present the results for the transitional case ($Re_{t}=2000$) last, since it is known being a tough test both from the experimental and the numerical point of view.

4.2.1 Case $Re_{t}=3500$

In [37, 8], it was shown that tetrahedral mesh $1200k$ (see Table 6 (left)) has a level of refinement that gives numerical results obtained with a finite element method in excellent agreement with the experimental data, despite the average mesh size being roughly 20 times larger than the Kolmogorov length scale at $Re_{t}=3500$. Thus, we are going to consider meshes with a similar level of refinement or coarser meshes.

It is known that the numerical results for the FDA benchmark test are particularly sensitive to small perturbations in the mesh, geometry, discretization parameters, and numerical schemes [45, 37, 8, 54, 7, 24, 15, 9, 14, 17, 34, 53]. We start with investigating the effect of interpolation scheme for the convective term. In particular, we consider two second-order accurate schemes: Linear Upwind Differencing (LUD) [51], which was used in Sec. 4.1, and Central Differencing (CD) [29]. Fig. 10 reports the normalized axial velocity (56) and the normalized pressure drop (57) along the $z$ axis from the NSE solution on meshes $100k$ and $320k$ (see Table 5) obtained by using the LUD and the CD convective schemes. We set $CFL=0.6$ as in [54]. The jet given by the NSE algorithm with the LUD convective scheme is much longer than the experimental jet for both meshes, while the NSE algorithm with the CD convective scheme provides a smaller jet length. It is known that upwind schemes introduce large diffusive errors and have a significant damping effect on the energy spectrum of the flow [43]. Therefore, central difference schemes are to be preferred for LES [43]. From now on, we will use the CD scheme for both NSE and EFR algorithms. In addition, we will use the original version of PISO since the flux correction term has been shown to introduce artificial dissipation [50]. From Fig. 10 (left) we observe that the pressure computed with either scheme and mesh is in excellent agreement with the experimental data for $z<0.3$. For $z>0.3$ the difference in the pressures computed with the two schemes for the convective term becomes large and experimental data are not available.

Next, we vary the $CFL$ number. Fig. 11 shows the normalized axial velocity (56) and the normalized pressure drop (57) along the $z$ axis from the NSE solution on meshes $100k$ and $320k$ (see Table 5) obtained with $CFL$ number 0.2 and 0.6. We see that a lower $CFL$ number produces a slightly shorter jet. In [54] the same low sensitivity to the $CFL$ number was observed when the RK4 scheme is used for time differentiation, while the numerical results obtained with the TFV4A scheme are extremely sensitive to the choice of $CFL$ number. From now on, we will set the $CFL$ number to 0.6.

We now consider all the meshes in Table 5. The associated results from the NSE solution are displayed in Fig. 12. We see that the axial velocities computed with meshes $1600k$ and $2700k$ match quite well the measurements all along the portion of the $z$-axis under consideration ($-0.088<z<0.08$). These results are in agreement with those in [8], where it is shown that a DNS is possible with the tetrahedral meshes $1200k$ and $1900k$ (see Table 6 (left)). Such meshes have number of cells, average number of cells along a throat diameter, and values of $h_{min}$ comparable to meshes $1600k$ and $2700k$ used here. However, while in [8] it is observed that instabilities in the NSE algorithm lead to a simulation breakdown for the meshes coarser than mesh $1200k$, no instability arises with coarser meshes when a finite volume method is adopted. From Fig. 12, we see that the results obtained with meshes $1600k$ and $2700k$ are in good agreement with the experimental data. However, we do not observe a clear convergence for the average velocity, i.e. the curves do not get closer to each other as the mesh is refined. This suggests that at $Re_{t}$ = 3500 the flow is still very sensitive to perturbations, since it is close to the transitional regime. See also Sec. 4.2.4. Thus, much finer meshes are needed for a DNS, as reported in [54]. As for the pressure, we observe excellent agreement with the experimental data for all the meshes but the coarsest, which underestimates the pressure drop in the entrance ragion. In particular, Fig. 12 (right) shows great agreement in the conical convergent, which could never be achieved with the finite element method used in [8].

Remark 4.2**.**

Additional tests (not reported for brevity) showed that small perturbations in the mesh strongly affect the results. For example, meshes very similar to $1600k$ and $2700k$ in terms of maximum non-orthogonality, average number of cells along the throat diameter, and mesh size but with slightly different refinement near the sudden expansion produced very long jets. Thus, we stress the importance of choosing the correct features in mesh generation, as already pointed out in [45, 37, 8, 54, 7, 24, 15, 9, 14, 17, 34, 53].

We have noticed that the choice $\chi=\chi_{2}$ introduces slightly more diffusion that $\chi=\chi_{1}$ for all the coarser meshes under consideration, thereby providing results in better agreement with the experimental data. Thus, we set $\chi=\chi_{2}$ whenever we use the EFR algorithm. In Fig. 13, we report the comparison between EFR solutions and experimental data in terms of normalized axial velocity (56) and pressure difference (57) for all the meshes in Table 5 coarser than mesh $1600k$. We see that the results computed with mesh $770k$ are in very good agreement with the measurements all along the $z$-axis. With meshes $100k$ and $320k$, the EFR algorithm provides a jet that starts to break down closer to the sudden expansion. Nevertheless, for all three meshes we notice that the EFR algorithm significantly improves the agreement with the experimental data with respect to NSE: compare Fig. 12 with 13.

Despite the fact that a Finite Volume approximation is used also in [54], a direct comparison between our results and those in [54] is complicated because we use difference LES approaches. In [54], the authors chose to use meshes tetrahedral elements and avoid local refinement. As a result, their meshes feature a higher level of refinement than ours. Moreover, the results in [54] are very sensitive to the time discretization scheme, $CFL$ condition and grid resolution. However, it is shown that introducing fluctuations in the upstream boundary condition drastically reduces this sensitivity and makes the prediction of the jet break-down point very robust.

4.2.2 Case $Re_{t}=5000$

At $Re_{t}=5000$ all the experiments in [21] observed turbulence downstream of the sudden expansion with a reproducible jet breakdown point. Just like for the $Re_{t}=3500$ case, we will first check the results given by the NSE algorithm and then consider the EFR algorithm.

Fig. 14 shows the normalized axial velocity (56) and the normalized pressure drop (57) along the $z$ axis computed from the NSE solution on all the meshes in Table 5 coarser than mesh $2700k$. First, we observe that that the axial velocities computed with meshes $1600k$ and $2700k$ match quite well the measurements all along the portion of the $z$-axis under consideration ($-0.088<z<0.08$), just like the $Re_{t}$ = 3500 case. Unlike the $Re_{t}$ = 3500 case, we do observe convergence as the mesh is refined. We remark that finite element method in [8] did not provide results in good agreement with the experimental data with a tetrahedral mesh $3000k$ (see Table 6). We believe this is due to the fact that mesh $3000k$ in [8] has a lower average number of cells throat diameter and a greater $h_{max}$ than meshes $1600k$ and $2700k$. See Tables 5 and 6. As for the pressure, from Fig. 14 (right) we see excellent agreement with the experimental data for all the meshes but the coarsest, which slightly underestimates the pressure drop in the entrance ragion.

Next, we consider the EFR algorithm. We report the comparison between computed and measured normalized axial velocity (56) and pressure drop (57) in Fig. 15. We see that the results computed with meshes $320k$ and $770k$ are in very good agreement with the measurements all along the $z$-axis. It is remarkable that the results given by the EFR algorithm on a very coarse mesh like $320k$ match so well the experimental data. As we will see in the next subsection, this is the case also for $Re_{t}=6500$. A possible reason why the results given by the EFR algorithm with mesh $320k$ were not as good for $Re_{t}=3500$ might be that $Re_{t}=3500$ is still close to the transitional regime. From Fig. 15 (left) we see that the EFR algorithm with mesh $100k$ provides a jet that starts to break down closer to the sudden expansion, just like the $Re_{t}=3500$ case.

From the tests performed for the $Re_{t}=3500$ and $Re_{t}=5000$ cases, we learned that the computed jet breakdown location is very sensitive to grid features (kind of element, non-orthogonality, local refinement, aspect ratio) and interpolation scheme for the convective term. Based on our experience, a robust FV-based CFD framework requires: (i) non-orthogonality and low skewness hexahedral meshes having cells with a contained element aspect ratio, and (ii) the CD scheme for the interpolation of the convective term. With this framework, we obtained favorable comparisons with the experimental data without introducing inlet perturbations as in [54].

4.2.3 Case $Re_{t}=6500$

The third flow regime we focus on features the highest Reynolds number considered by FDA: $Re_{t}=6500$. Also in this case, turbulence downstream of the sudden expansion was observed in all the experiments reported in [21] with a reproducible jet breakdown point.

The $Re_{t}=6500$ case is thoroughly investigated in [24]. The author uses the Smagorinsky subgrid-scale model [44] and combines it with a FV method implemented in the commercial CFD code ANSYS Fluent. Great agreement with the experimental data is obtained with a very fine hexahedral block-structured mesh that has roughly 9 millions elements ($h_{min}$ =2.5e-5) and excellent features in terms of non-orthogonality and skewness. We notice that the meshes used in [24] are very similar to those ones we use in this work.

Given the good EFR results shown in Sec. 4.2.1 and 4.2.2, for $Re_{t}=6500$ we limit the investigation to the EFR algorithm with $\chi=\chi_{2}$ and all the meshes in Table 5 coarser than mesh $1600k$. We report the comparison between computed and measured normalized axial velocity (56) and pressure drop (57) in Fig. 16. From Fig. 16 (left) we see that the axial velocity computed with meshes $330k$ and $770k$ is in very good agreement with the measurements all along the $z$-axis. As for mesh $100k$, we observe the same trend as for the $Re_{t}=3500,5000$ cases: the computed jet starts to break down closer to the sudden expansion, while still in good agreement with the experimental data in the rest of the domain. As for the pressure, from Fig.16 (right) we observe a significant overestimation in the conical convergent respect to the experimental data for all the meshes. A relaxation of the pressure (see Remark 2.2) might lead to a better agreement.

4.2.4 Case $Re_{t}=2000$

The transitional case ($Re_{t}=2000$) proved to be challenging from both the experimental and computational sides, thus we present it last. From the experimental side, the interlaboratory velocity data in [21] showed significant differences downstream of the sudden expansion, in particular concerning the jet breakdown point. This was attributed mainly to a 10% higher flow rate (and consequently higher $Re_{t}$), which caused premature jet breakdown in two experiments out of five [21]. However, minor differences in the fabricated geometrical models and inlet perturbation levels played a role also. From the numerical point of view, we found the results to be very sensitive to the mesh refinement level.

For this test, we look first at the results obtained with the NSE algorithm. Fig. 17 reports the normalized axial velocity (56) and the normalized pressure drop (57) along the $z$ for mesh $2700k$ and coarser (see Table 5). From Fig. 17 (left), we notice that the axial velocity computed with mesh $1600k$ (resp., $2700k$) match quite well the measurements giving a shorter (resp., longer) jet. This agrees with what reported in [37], where it is shown that the NSE results obtained with a mesh with 2.5 millions tetrahedra compare well with the experimental data giving a longer jet. However, the numerical axial velocity computed with meshes $1600k$ and $2700k$ are much further apart for $Re_{t}=2000$ than for $Re_{t}=3500$. Compare Fig. 17 (left) with Fig. 12 (left). Thus, we cannot argue that either mesh is refined enough for this test case. From Fig. 17 (right), we see that the simulated pressure drop is in very good agreement with the experimental data all along the $z$-axis and regardless of the mesh.

Next, we consider the EFR algorithm. We report the comparison between computed and measured normalized axial velocity (56) and pressure drop (57) in Fig. 18. We see that the axial velocity computed with mesh $770k$ is in very good agreement with the measurements giving a shorter jet. To reach such an agreement with the NSE algorithm we had to use mesh $1600k$. However, just like for the NSE algorithm, the numerical axial velocity computed with the different meshes are further apart for $Re_{t}=2000$ than for all the Reynolds numbers considered in the previous sections. Compare Fig. 18 (left) with Fig. 13 (left), 15 (left), and 16 (left). These results confirm the difficulty in the numerical simulation of the transitional regime.

5 Conclusions and Perspectives

We showed the effectiveness of a FV-based EFR algorithm in simulating flow problems at moderately large Reynolds numbers. Nonlinear filtering stabilizes marginally resolved scales without over-diffusing, thereby allowing to use less degrees of freedom than required by a DNS. To select the regions of the domain where filtering is needed, we employ a nonlinear differential low-pass filter. The interest in the Finite Volume approximation is due to the fact that it has been widely used in the LES context. However, the application of the Leray model in a FV framework has been unexplored.

In order to showcase the features of our approach, we presented a computational study related to two benchmarks: 2D flow past a cylinder and a 3D benchmark from the FDA. With reference to the FDA benchmark, we performed a complete characterization of the flow at Reynolds numbers starting from 2000 up to 6500. We proposed a new formula based on physical and numerical arguments to tune a value of $\chi$ that leads to very good agreement with the experimental measurements. Several meshes were considered to understand how under-refined the mesh can be while still capturing the physical average quantities. Furthermore, we investigated the impact of the mesh features as well as the discretization of the convective term on the results obtained.

As a follow-up of the present work, we are going to develop a Leray Reduced Order Model (ROM) within a Finite Volume framework. We are also interested in coupling the Leray model with an elasticity model to simulate fluid-structure interaction problems which are ubiquitous in science and engineering.

6 Acknowledgements

We acknowledge the support provided by the European Research Council Executive Agency by the Consolidator Grant project AROMA-CFD “Advanced Reduced Order Methods with Applications in Computational Fluid Dynamics” - GA 681447, H2020-ERC CoG 2015 AROMA-CFD and INdAM-GNCS projects. This work was also partially supported by US National Science Foundation through grant DMS-1620384.