2.1 The conservation laws

The spatial discretization x is based on terrain-following coordinates (Gal-Chen and Somerville, 1975). The staggered mesh is regular $\mathrm{\Delta }x=\mathrm{\Delta }y=\mathrm{\Delta }$ in the horizontal directions and a transformation of the vertical one is available in order to fit a non-plane surface. The release of the vertical space step is available wherein a fine resolution is unnecessary. In the current study, only flat problems are considered with a Δz=Δ restriction for altitudes in the presence of immersed obstacles.

The core of the MNH dynamic in its dry version is based on the resolution of the Euler and thermodynamic equations (energy preserving). The anelastic approximation (Lipps and Hemler, 1982; Durran, 1989) is assumed; the reference state is stratified, and the density deviation to the hydrostatic case in the buoyancy term is considered. The system can be simplified into the Boussinesq approximation when considering a uniform reference state. The tendencies of each prognostic variable ψ satisfying the usual conservation laws in MNH are expressed as $\frac{\partial \mathit{\psi }}{\partial t}{|}_{\mathrm{csl}}$, where the subscript csl is used to distinguish these tendencies from those arising from the application of the IBM procedure (Sect. 3.1). The prognostic variables are the resolved momentum, the potential temperature and if necessary an arbitrary passive scalar. The prognostic variable is decomposed into a resolved component, $\stackrel{\mathrm{‾}}{\mathit{\psi }}$, and an unresolved component, ψ^′ (${\mathit{\psi }}^{\prime }=\mathrm{0}$ in a direct numerical simulation, DNS). An additional prognostic equation on the subgrid turbulent kinetic energy (STKE) is solved for a large-eddy simulation (LES; Sect. 2.3). The potential temperature is defined through the Exner function Π and the absolute temperature $T=P/\left({\mathit{\rho }}_{\mathrm{r}}{R}_{\mathrm{d}}\right)=\mathit{\theta }(P/P_{\mathrm{r}\mathrm{0}}{)}^{R_{\mathrm{d}}/C_p}=\mathit{\theta }\mathrm{\Pi }$, where ρ_r is the density of the reference state, P is the local pressure, P_r0 the reference ground-level pressure, R_d the gas constant and C_p the specific heat capacity at a constant pressure for dry air. The thermodynamic equations and an additional passive scalar equation are

where ${\stackrel{\mathrm{‾}}{F}}^{\mathrm{\Pi }}$ corresponds to pressure effects. The transport of each prognostic scalar in Eqs. (1), (2) and (6) is made by a piecewise parabolic method (PPM) with undershoot and overshoot limitation (Colella and Woodward, 1984; Lin and Rood, 1996). The temporal algorithm of the advection term in these scalar transports is a forward-in-time scheme (noted FT). The momentum equations are

where $\stackrel{\mathrm{‾}}{\mathit{u}}$ is the resolved wind, g the acceleration due to the gravity appearing in the buoyancy term, θ_r is the potential temperature of the reference state, μ_f the dynamic viscosity and $\mathrm{\nabla }\cdot \left({\mathit{\rho }}_{\mathrm{r}}\stackrel{\mathrm{‾}}{{\mathbf{u}}^{\prime }\otimes {\mathbf{u}}^{\prime }}\right)$ the Reynolds stresses. The spatial discretization of $\mathrm{\nabla }\cdot ({\mathit{\rho }}_{\mathrm{r}}\stackrel{\mathrm{‾}}{\mathit{u}}\otimes \stackrel{\mathrm{‾}}{\mathit{u}})$ in Eq. (3) can be done by second- or fourth-order centered schemes and third- or fifth-order weighted essentially non-oscillatory schemes (Jiang and Shu, 1996). The temporal evolution of the resolved wind is achieved by a fourth-order explicit Runge–Kutta (ERK) algorithm (Shu and Osher, 1989; Lunet et al., 2017). In the present study Δt is fixed to respect the Courant number $\frac{\mid {\stackrel{\mathrm{‾}}{\mathit{u}}}^n\mid \mathrm{\Delta }t}{\mathrm{\Delta }}<\mathrm{1}$ (n, the time step index) and no additional time splitting is implied. The temporal viscous stability condition 𝒪(ν_f∕Δ²) (ν_f the kinematic viscosity) imposes an additional restriction when the viscous term is explicitly resolved in time.

The bottom, lateral wall and top surfaces take a free-slip, impermeable and adiabatic behavior without the call of an externalized surface scheme. The open boundary condition is a Sommerfeld equation defined as wave radiation (Carpenter, 1982) to enforce the large scales and allow for the reflection wave damping.

2.2 The incompressibility condition

The wind of the resolved scales has to satisfy the continuity equation $\mathrm{\nabla }\cdot \left({\mathit{\rho }}_{\mathrm{r}}{\stackrel{\mathrm{‾}}{\mathit{u}}}^{n+\mathrm{1}}\right)=\mathrm{0}$. The method consists of the projection of the predicted velocity field ${\stackrel{\mathrm{‾}}{\mathit{u}}}^{\ast }$ (solution of Eq. 3 without the pressure term) into the null divergence subspace. This projection estimates the irrotational correction to apply to ${\stackrel{\mathrm{‾}}{\mathit{u}}}^{\ast }$ through a potential scalar Ψ^∗:

Ψ^∗ is obtained with the resolution of the pseudo-Poisson equation written as

The horizontal part of the operator to invert in the elliptic problem is treated in the Fourier space (Schumann and Sweet, 1988), and its vertical part leads to the classical tridiagonal matrix. The mathematical operator to invert ∇⋅(∇) is exact for flat problems (Bernadet, 1995). When the mesh is built with terrain-following coordinates over a flat surface, the solution of the pressure problem becomes inaccurate. In this orography presence case, an iterative procedure is employed such as a Richardson, a conjugate gradient (Young and Jea, 1980) or a residual conjugate gradient (Skamarock et al., 1997) algorithm.

2.3 The turbulent subgrid scales

To execute LES, the Reynolds stress $\mathrm{\nabla }\cdot \left({\mathit{\rho }}_{\mathrm{r}}\stackrel{\mathrm{‾}}{{\mathbf{u}}^{\prime }\otimes {\mathbf{u}}^{\prime }}\right)$ appearing in Eq. (3) are estimated. The LES closure is done by an eddy–diffusivity approach called 1.5TKE with a 1.5-order closure scheme (Cuxart et al., 2000). The isotropic part of the subgrid turbulence is given by the prognosis of the subgrid turbulent kinetic energy $e=\frac{\mathrm{1}}{\mathrm{2}}(\stackrel{\mathrm{‾}}{{u^{\prime }}^{\mathrm{2}}}+\stackrel{\mathrm{‾}}{{v^{\prime }}^{\mathrm{2}}}+\stackrel{\mathrm{‾}}{{w^{\prime }}^{\mathrm{2}}})$:

where K_e and K_ϵ are constants prescribed in the turbulence scheme, and l_m and l_ϵ are the length scales defining the turbulent viscosity. The dissipation term is directly estimated from e and l_ϵ (the left-hand term in Eq. 6). The anisotropic part of the subgrid turbulence is diagnosed from the $\stackrel{\mathrm{‾}}{\mathit{\psi }}$ gradient and e. The diagnosis of the anisotropic part of the subgrid turbulence is obtained using $\stackrel{\mathrm{‾}}{\mathit{\psi }}$ and e.

The Einstein summation convention is applied and Φ_i,m represents atmospheric stability functions (Redelsperger and Sommeria, 1981). The ground condition can be modeled by the externalized surface scheme SURFEX (Masson et al., 2013). In the dry version of MNH and with the hypothesis of zero thermal flux at the ground and buildings, only the turbulent friction is used. To compute the nonzero values of $\stackrel{\mathrm{‾}}{u_i^{\prime }{u}_j^{\prime }}$ at the ground, the SURFEX call employed in this paper consists of the simple activation of a dynamic wall model related to the Prandtl theory (eddy viscosity concept). The form of the surface turbulent fluxes is $\stackrel{\mathrm{‾}}{u_i^{\prime }{u}_j^{\prime }}{|}_{\mathrm{surf}}\approx -l_m^{\mathrm{2}}\mid \frac{d\stackrel{\mathrm{‾}}{u_i}}{d{x}_j}\mid \frac{d\stackrel{\mathrm{‾}}{u_i}}{d{x}_j}$. Defining a friction velocity u^∗ proportional to the turbulent wall shear and a roughness length z₀, the vertical gradient of $\stackrel{\mathrm{‾}}{\mathit{u}}$ is recovered by specifying a logarithmic profile (von Kármán, 1930) as $\stackrel{\mathrm{‾}}{u}\left(z\right)=\frac{u^{\ast }}{k}\ln (\mathrm{1}+\frac{z}{z_{\mathrm{0}}})$ (note that the atmospheric stability conditions are neutral or near neutral in this paper; therefore, the additional Monin and Obukhov (1954) term is neglected and Businger et al. (1971) functions do not appear in the previous formulation). SURFEX is employed in Sect. 5.2.

3.1 Ghost-cell technique and explicit-in-time schemes

The ψⁿ value is estimated at the time nΔt (Δt, the time step). The tendencies of the prognostic variables $\stackrel{\mathrm{‾}}{\mathit{\psi }}=[\stackrel{\mathrm{‾}}{\mathit{u}},\stackrel{\mathrm{‾}}{\mathit{\theta }},\stackrel{\mathrm{‾}}{s},(e\left)\right]$ (Sect. 2.1) cannot be deduced from the conservation laws in the solid region. Expecting a correction due to IBM, wherein ϕ≥0, a general formulation of the tendencies is written as

The right-hand-side (RHS) first term of Eq. (10) is given by the conservation laws (Sect. 2.1). The $\frac{\partial }{\partial t}{|}_{\mathrm{ibm}}$ tendency is the correction due to the GCT in the solid region and near the immersed interface, satisfying the $\stackrel{\mathrm{‾}}{\mathit{\psi }}$ desired boundary conditions at ϕ=0:

Note that $\frac{\partial \stackrel{\mathrm{‾}}{\mathit{u}}}{\partial t}{|}_{\mathrm{ibm}}$ is taken into account in the ERK temporal algorithm. The freeze of the immersed wind conditions in the ERK algorithm has also been implemented; it has shown more unstable behavior for a large Courant number.

The forced points are called ghost points and are renamed ghosts. To estimate the variable $\stackrel{\mathrm{‾}}{\mathit{\psi }}$ and for each ghost, the physical information is extracted near the interface and from the fluid region. The extension (grid stencil) of the forcing zone depends on the spatial accuracy of the numerical scheme. For example, Fig. 2b–c show the case of a two-layer stencil in a two-dimensional grid. The characterization of the layer is done by a conditional loop applied direction by direction on the LSF. For a 2-D case, the sign of $\mathit{\varphi }(i,j)\cdot \mathit{\varphi }(i,[j-k_l:j+k_l\left]\right)$ and $\mathit{\varphi }(i,j)\cdot \mathit{\varphi }\left(\right[i-k_l:i+k_l],j)$ is estimated. The integer value k_l determines the cells truncated by the interface: k_l=1 (k_l=2) defines the first (second) layer. The calculation of these ghost layers has a computational overhead due to data exchange among processors in parallel simulations. The stencil of the numerical scheme modeling the interface defines the k_l value. In order to limit this overhead, a low-order version of a centered explicit-in-time scheme (Sect. 2.1) is employed when $\mathit{\varphi }>-\mathrm{\Delta }$. The CPU cost of the “hybrid” advection scheme is largely compensated for by the decrease in the ghost number and parallel exchanges. Appendix B reports a comparison analysis between third-order weighted essentially non-oscillatory (WENO) and second-order centered schemes used in the vicinity of the interface; the studied case is the inviscid flow around a circular cylinder.

Figure 2(a) Node definitions acting in the ghost-cell technique: the ghost (G), interface (B), normal vector (n), mirror (I) and images (I₁,I₂). (b, c) Illustration of classical (new) GCT using the mirror (images). Triangles correspond to one of the node types (see Fig. 1b).

Download

In classical GCT (Tseng and Ferziger, 2003) the fluid information is obtained at a mirror point (noted I, renamed mirror) found in the normal direction to the interface in such a way that the interface node B is equidistant to I and G. Figure 2a shows the characterization of one ghost G (of LSF value ϕ_G), its associated mirror I (of LSF value ϕ_I) and the interface node B (GI=2ϕ_Gn). Figure 2b illustrates several ghosts and mirrors. The $\left|IB\right|$ distance depends on the forcing stencil, and a problematic case regularly met in the mirror interpolation is the vicinity of ghosts with the interface (${\mathit{\varphi }}_G=-{\mathit{\varphi }}_I\ll \mathrm{\Delta }$, with Δ the space step), leading to a not well-posed condition.

The new GCT. To overcome this problem, we define image points (noted I₁ and I₂ in Fig. 2a; renamed images) having a distance to the interface that depends only on the grid spacing: $G{I}_l=l\mathrm{\Delta }+{\mathit{\varphi }}_G\mathit{n}$ with $l=(\mathrm{1};\mathrm{2})$. Figure 2a shows the images for one ghost. The new approach enforces a large enough value of the $\left|I_{l}B\right|$ distance. Figure 2b (c) illustrates classical (new) GCT for several ghosts. Figure 2b shows some mirror points associated with ghosts of the first layer in the vicinity of the interface. Figure 2c shows that the new approach ensures the image points are always located in the fluid region, irrespective of the ghost location. The definition of several images per ghost allows us to build a profile of the $\stackrel{\mathrm{‾}}{\mathit{\psi }}$ fluid information normal to the interface. $\stackrel{\mathrm{‾}}{\mathit{\psi }}\left(I\right)$ is therefore recovered by a quadratic reconstruction f using the $(B,I_{\mathrm{1}},I_{\mathrm{2}})$ points. Two distinct calculations of $f\left(\stackrel{\mathrm{‾}}{\mathit{\psi }}\right(B),\stackrel{\mathrm{‾}}{\mathit{\psi }}(I_{\mathrm{1}}),\stackrel{\mathrm{‾}}{\mathit{\psi }}(I_{\mathrm{2}}\left)\right)$ noted PLI^a and PLI^b are tested to build the Lagrangian interpolation:

where L^a(I) and L^b(I) are the Lagrangian polynomials, as follows.

Figure 3Quadratic interpolations of two analytical profiles $\stackrel{\mathrm{‾}}{\mathit{\psi }}={\mathit{\varphi }}^n$ (red lines) using two image points at $\mathit{\varphi }=-[\mathrm{1};\mathrm{2}]\mathrm{\Delta }$ and the interface node. Green (blue) corresponds to L^a (L^b) polynomial results.

Download

The accuracy of an interpolation depends on the $\stackrel{\mathrm{‾}}{\mathit{\psi }}$ profile. For example, a logarithmic evolution of the tangent velocity is expected in LES. Otherwise, when the viscous layer is modeled, a linear evolution is expected. To compare the ability of each quadratic interpolation to approach a wide variety of profiles, the recovery of power laws such as $\mathit{\psi }={\mathit{\varphi }}^{\mathrm{3}/\mathrm{2}}$ (Fig. 3a) and $\mathit{\psi }={\mathit{\varphi }}^{\mathrm{1}/\mathrm{4}}$ (Fig. 3b) is studied. As illustrated, PLI^a fits the two analytical solutions better and is therefore adopted. The classical and new GCTs have been compared thoroughly, and part of this analysis is deferred to the Appendix. The interpolated field of the potential flow around a single cylinder or a sphere was compared to the theoretical solution; the sensitivity of the inviscid flow around the same bodies (Appendix B) to the type of GCTs has been studied. The new GCT has given the best results, especially in the symmetry preservation in the inviscid flow cases. Note that these results are also dependent on the 3-D interpolation choice detailed in the following paragraph. The proposed GCT is employed in the rest of this study. The GCT implementation is divided into four main steps: the fluid information recovery, the interface basis change, the interface condition and the ghost value.

The fluid information recovery. ${\stackrel{\mathrm{‾}}{\mathit{\psi }}}_{I_l}$, for the images contained in a pure fluid cell (all corner nodes are in the fluid region), is recovered by a trilinear interpolation based on Lagrangian polynomials (LP), as follows.

For truncated cells (at least one corner node is in the solid region), ${\stackrel{\mathrm{‾}}{\mathit{\psi }}}_{I_l}$ is recovered using an inverse distance-weighting (DW) interpolation:

where $L_i^{\mathrm{DW}}\left(x\right)=\mid x_l-x_i{\mid }^{-\mathit{\alpha }}$ (α=1). This formulation diverges when x_i→x_l and it is commonly adopted to impose $\stackrel{\mathrm{‾}}{\mathit{\psi }}\left(x_l\right)=\stackrel{\mathrm{‾}}{\mathit{\psi }}\left(x_i\right)$ when $\exists (x_i-x_l)\le \mathit{ϵ}$ (ϵ is an arbitrary parameter depending on the mesh discretization, ϵ≪Δ). The 3-D extension is direct with $\mid {\mathbf{x}}_{\mathbf{l}}-{\mathbf{x}}_{\mathbf{i}}\mid =\sqrt{(x_l-x_i{)}^{\mathrm{2}}+(y_l-y_i{)}^{\mathrm{2}}+(z_l-z_i{)}^{\mathrm{2}}}$. The use of these interpolations was decided after comparisons with barycentric Lagrangian and modified distance-weighting interpolations (Franke, 1982) and tests on the α coefficient. As the boundary condition is expressed in the interface frame and the grid is staggered, the non-collocation of the $\stackrel{\mathrm{‾}}{\mathit{u}}$ components implies the interpolation of three different classes of cells (with $U,V,W$ corners; Fig. 1b) for each $U,V,W$ ghost and building the change of frame matrix for which the proposed GCT presents an interest during the characterization of the direction tangent to the interface.

The interface basis change. Velocity vector $\stackrel{\mathrm{‾}}{\mathit{u}}$, known in the Cartesian mesh basis at the images I_l (Δn₁=Δ and Δn₂=2Δ in Fig. 2a), is projected in the basis of the interface $\left(\mathit{n}\right(B),\mathbf{t}(B),\mathbf{c}(B\left)\right)$ in which the boundary conditions on each vector component are imposed. Computing the LSF gradient, the normal direction is defined. Otherwise, (t,c) represents two arbitrary tangent directions. The tangent direction t is considered the predominant tangent direction of the flow along the fluid–solid interface depending on the image values and defining the velocity vector as $\stackrel{\mathrm{‾}}{\mathit{u}}\left(I_l\right)={\stackrel{\mathrm{‾}}{u}}_n\left(I_l\right)\mathit{n}+{\stackrel{\mathrm{‾}}{u}}_t\left(I_l\right)\mathbf{t}\left(I_l\right)$. The cotangent direction is defined as

The $(\mathit{n},\mathbf{t},\mathbf{c})$ basis at the interface is defined by considering (or not) the rotation of the tangent velocity with the distance to the interface:

Finally, the third component is $\mathbf{c}\left(B\right)=\mathit{n}\otimes {\mathbf{e}}_{\mathbf{t}}\left(B\right)$ and (inverse) projection is known.

The interface condition. Let ${\stackrel{\mathrm{‾}}{\mathit{\psi }}}_B$ and $\mathrm{\Delta }\frac{\partial \stackrel{\mathrm{‾}}{\mathit{\psi }}}{\partial n}{|}_B$ be the Dirichlet and Neumann conditions on $\stackrel{\mathrm{‾}}{\mathit{\psi }}$. The general formulation of the boundary condition $\stackrel{\mathrm{‾}}{\mathit{\psi }}(\mathit{\varphi }=\mathrm{0})$ is written as a Robin condition: $\stackrel{\mathrm{‾}}{\mathit{\psi }}(\mathit{\varphi }=\mathrm{0})=k_{\mathrm{r}}{\stackrel{\mathrm{‾}}{\mathit{\psi }}}_B+(\mathrm{1}-k_{\mathrm{r}}).\left(\stackrel{\mathrm{‾}}{\mathit{\psi }}\right(\mathit{\varphi }=-l\mathrm{\Delta })-l\mathrm{\Delta }\frac{\partial \stackrel{\mathrm{‾}}{\mathit{\psi }}}{\partial n}{|}_B$). The switch between the Dirichlet condition and the Neumann condition is done through the coefficient $k_{\mathrm{r}}\in [\mathrm{0}:\mathrm{1}]$. To give some examples of a Dirichlet condition, $(k_{\mathrm{r}};{\stackrel{\mathrm{‾}}{\mathit{\psi }}}_B)=(\mathrm{1};\mathrm{0})$ is imposed on the $\stackrel{\mathrm{‾}}{\mathit{u}}\cdot \mathit{n}$ velocity component normal to the interface arising from the impermeability hypothesis and on the $\stackrel{\mathrm{‾}}{\mathit{u}}\cdot \mathbf{t}$ component tangent to the interface for a no-slip hypothesis. To give some examples of the Neumann condition imposed by $(k_{\mathrm{r}};\frac{\partial \stackrel{\mathrm{‾}}{\mathit{\psi }}}{\partial n}{|}_B)=(\mathrm{0};\mathrm{0})$: a no-flux condition on the potential temperature (as well on a passive scalar, subgrid kinetic energy) and a free-slip case applied to $\stackrel{\mathrm{‾}}{\mathit{u}}\cdot \mathbf{t}$. Note the $\frac{l\mathit{\varphi }}{\mathrm{2}}$ approximation in the location of the derivative term and the Neumann condition depending on the chosen image (in practice the selected image I_l is the closest one to the interface).

An interface condition depending on the characteristics of the surrounding fluid such as $\stackrel{\mathrm{‾}}{\mathit{\psi }}(\mathit{\varphi }=\mathrm{0})=F({\stackrel{\mathrm{‾}}{\mathit{\psi }}}_{I_l};\frac{\partial \stackrel{\mathrm{‾}}{\mathit{\psi }}}{\partial n}{|}_{I_l})$ is a wall model. Using two (three) images, simple wall models such as the constant (linear) extrapolation of the $\stackrel{\mathrm{‾}}{\mathit{\psi }}$ gradient is reached by the $\frac{{\partial }^{\mathrm{2}}\stackrel{\mathrm{‾}}{\mathit{\psi }}}{\partial n^{\mathrm{2}}}{|}_{I_l}=\mathrm{0}$ ($\frac{{\partial }^{\mathrm{3}}\stackrel{\mathrm{‾}}{\mathit{\psi }}}{\partial n^{\mathrm{3}}}{|}_{I_l}=\mathrm{0}$) computation. The consistency between the tangent component to the interface of the resolved wind and the subgrid turbulence is the subject of Sect. 2.3.

The ghost value. Knowing $\stackrel{\mathrm{‾}}{\mathit{\psi }}(\mathit{\varphi }=\mathrm{0})$ and ${\stackrel{\mathrm{‾}}{\mathit{\psi }}}_{I_l}=\stackrel{\mathrm{‾}}{\mathit{\psi }}(\mathit{\varphi }=-l\mathrm{\Delta })$, ψ(G) for a Dirichlet (Neumann) condition is written as

Figure 4Profile normal to the interface of two points of fluid information $\stackrel{\mathrm{‾}}{\mathit{\psi }}$: analytical solution (red lines), quadratic interpolation QI₁ using ${\stackrel{\mathrm{‾}}{\mathit{\psi }}}_{-\mathit{\varphi }=[\mathrm{1}/\mathrm{2};\mathrm{1}]\mathrm{\Delta }}$ (green symbols), QI₂ using ${\stackrel{\mathrm{‾}}{\mathit{\psi }}}_{-\mathit{\varphi }=[\mathrm{1};\mathrm{2}]\mathrm{\Delta }}$ (blue symbols), QI_C as a combination of QI₁ and QI₂ (purple lines).

Download

In practice, three I_l images are defined for which the locations are ${\mathit{\varphi }}_{I_l}=-l\mathrm{\Delta }$ with $l=[\mathrm{1}/\mathrm{2};\mathrm{1};\mathrm{2}]$. The choice of the image distance to the interface affects the results. To best approach the expected solution, two quadratic interpolations depending on the images used and one combination of these quadratic interpolations are tested. Figure 4a and b illustrate these interpolations by considering two analytical profiles (red lines): the quadratic interpolation QI₁ (QI₂) is based on the image values located at $\mathit{\varphi }=\mathrm{1}/\mathrm{2}\mathrm{\Delta }$ and ϕ=Δ and plotted in green symbols (at ϕ=Δ and ϕ=2Δ plotted in blue symbols). Depending on the analytical profile, Fig. 4 shows the influence of the image location choice. As expected, QI₁ (QI₂) appears to be less accurate than QI₂ (QI₁) for $\stackrel{\mathrm{‾}}{\mathit{\psi }}(\mathit{\varphi }\in [-\mathrm{2}\mathrm{\Delta }:-\mathrm{\Delta }\left]\right)$ (for $\stackrel{\mathrm{‾}}{\mathit{\psi }}(\mathit{\varphi }\in [-\mathrm{\Delta }:\mathrm{0}\left]\right)$). QI_C is the combination of QI₁ and QI₂ (purple line). QI_C preserves the advantage of each quadratic interpolation and when ϕ_G<Δ (ϕ_G>Δ), QI₁ (QI₂) is used in the rest of the study. Knowing ${\stackrel{\mathrm{‾}}{\mathit{\psi }}}^{n+\mathrm{1}}\left({\mathit{\varphi }}^{-}\right)$ at the end of the MNH temporal loop with QI_C, the ${\stackrel{\mathrm{‾}}{\mathit{\psi }}}^{n+\mathrm{1}}\left({\mathit{\varphi }}^{+}\right)$ profile is extrapolated from the fluid region to the solid region by applying an anti-symmetry ${\stackrel{\mathrm{‾}}{\mathit{\psi }}}^{n+\mathrm{1}}\left({\mathit{\varphi }}^{+}\right)=\mathrm{2}{\stackrel{\mathrm{‾}}{\mathit{\psi }}}^{n+\mathrm{1}}\left(\mathrm{0}\right)-{\stackrel{\mathrm{‾}}{\mathit{\psi }}}^{n+\mathrm{1}}\left({\mathit{\varphi }}^{-}\right)$. The ghost value is estimated and the $\stackrel{\mathrm{‾}}{\mathit{\psi }}$ gradient at the interface is also recovered.

3.2 Cut-cell technique and pressure solver

First looking at the RHS of Eq. (5), the $\frac{\partial \left({\mathit{\rho }}_{\mathrm{r}}{\stackrel{\mathrm{‾}}{\mathit{u}}}^{\ast }\right)}{\partial t}{|}_{\mathrm{csl}}$ coming from the resolution of the explicit-in-time schemes near the interface and in the solid regions badly affects the $\mathrm{\nabla }\cdot {\stackrel{\mathrm{‾}}{\mathit{u}}}^{\ast }$ computation (note that the GCT operates after the step projection). Therefore, the fictive wind of the solid region can spread errors in the fluid region during the pressure resolution. To avoid it, a correction of the pressure solver is proposed.

The elliptic problem (Eq. 5) is rewritten as a resolution of the linear system $\mathcal{P}\cdot {\mathrm{\Psi }}^{\ast }=\mathcal{Q}$. In the standard MNH version, $\mathrm{\nabla }\cdot {\stackrel{\mathrm{‾}}{\mathit{u}}}^{\ast }=\mathcal{Q}$ is estimated using a finite-difference approach. To uncouple the solid region from the fluid region our revisited version enforces a null divergence for pure solid cells and estimates the balance of momentum fluxes by a finite-volume approach for truncated cells (noted 𝒬_cct):

where 𝒱=Δ³ is the cell volume, 𝒱_f (𝒱_s) the fluid (solid) part of 𝒱 and 𝒮_i the cell surfaces for which i is the index of each surface orientation [e, w, n, s, b, f] as illustrated in Fig. 5a.

Figure 5(a) Momentum flux balance for an arbitrary truncated cell of volume 𝒱, where the ${\stackrel{\mathrm{‾}}{u_i}}^{\ast }$ velocities (U_i in the figure, $i=[e,w,n,s,b,f]$) are supported by the 𝒮_i surfaces in grey; the transparent volume is a part of the solid body. (b) Segmentation of the 𝒮_w arbitrary surface (red border) in eight $P_p^{\ast }$ pieces of cake (the border of $P_{\mathrm{3}}^{\ast }$ is indicated in green).

Download

According to the Green–Ostrogradski theorem, the ${\stackrel{\mathrm{‾}}{u_i}}^{\ast }{\mathcal{S}}_i$ calculation is the classical way of a CCT (Yang et al., 1997; Kim et al., 2001) to estimate the velocity divergence. A similar approach is performed here by rebuilding the flux $\stackrel{\mathrm{̃}}{{\mathrm{\Delta }}^{\mathrm{2}}{\stackrel{\mathrm{‾}}{\mathit{u}}}^{\ast }}$ for truncated and solid cells.

The $\pm \stackrel{\mathrm{̃}}{{\mathrm{\Delta }}^{\mathrm{2}}{\stackrel{\mathrm{‾}}{u_i}}^{\ast }}$ calculation consists of a weighting of the outflux and influx function of the fluid and cell surface ratios (Fig. 5a). Figure 5b gives an example of the west surface (i=w, red border) in which $\stackrel{\mathrm{̃}}{{\mathrm{\Delta }}^{\mathrm{2}}{\stackrel{\mathrm{‾}}{u_{\mathrm{w}}}}^{\ast }}(j,k)$ is calculated using the LSF value ϕ=ϕ_w and the ones of the eight adjacent nodes ${\mathit{\varphi }}_p(j\pm \mathrm{1},k\pm \mathrm{1})$. A disk of radius $\sqrt[-\mathrm{1}]{\mathit{\pi }}\mathrm{\Delta }$ is split into eight “piece-of-cake” segments $P_p^{\ast }$ ($p=[\mathrm{1}:\mathrm{8}]$). An LSF linear interpolation detects (or not) the interface location. In the presence of an interface, its distance from the studied node is $\mathrm{0}<{\mathit{\delta }}_p<\sqrt[-\mathrm{1}]{\mathit{\pi }}\mathrm{\Delta }$. Knowing δ_p, the momentum flux balance is formulated for a nonmoving body as (p is the index of the piece of cake and i the index of the cell surface)

The four encountered cases correspond to a pure fluid cell $\stackrel{\mathrm{̃}}{{\mathrm{\Delta }}^{\mathrm{2}}{\stackrel{\mathrm{‾}}{u_i}}^{\ast }}=\frac{{\mathrm{\Delta }}^{\mathrm{2}}}{\mathrm{8}}\stackrel{8}{\sum _{p=\mathrm{1}}}{\stackrel{\mathrm{‾}}{u_i}}^{\ast }$ when ϕ_p<0 and ϕ_i<0 (Fig. 6a); a pure solid cell $\stackrel{\mathrm{̃}}{{\mathrm{\Delta }}^{\mathrm{2}}{\stackrel{\mathrm{‾}}{u_i}}^{\ast }}=\mathrm{0}$ when ϕ_p>0 and ϕ_i>0 (Fig. 6b); and two types of truncated cells depending on the fluid–solid nature of the main node for which ${\mathit{\varphi }}_p.{\mathit{\varphi }}_i<\mathrm{0}$ (Fig. 6c–d). Using Eq. (23), Eq. (22) is solved and leads to the RHS computation of Eq. (5).

Figure 6(a–d) $\pm \stackrel{\mathrm{̃}}{{\mathrm{\Delta }}^{\mathrm{2}}{\stackrel{\mathrm{‾}}{u_i}}^{\ast }}$ calculations depending on the signs of ϕ_i=(ϕ) and ϕ_p on an arbitrary piece of cake. The white (grey) region corresponds to the solid (fluid) one of $P_p^{\ast }$ (same color code as in Fig. 5).

Download

Knowing 𝒬_cct, the reflection now concerns the 𝒫 matrix to invert. The classical interface condition on the potential Ψ^∗ is a homogeneous Neumann condition $\frac{\partial {\mathrm{\Psi }}^{\ast }}{\partial \mathit{\varphi }}=\mathrm{0}$. Using a boundary-fitted method (BFM), the interface condition of the moving or nonmoving body (Auguste, 2010) appears only on the border of a numerical domain. Using an IBM and without any impact of this interface condition on the 𝒫 coefficients, the impermeability character of solid obstacles is not achieved. Due to the inversion of the horizontal part of 𝒫 by a fast Fourier transform (Schumann and Sweet, 1988), the solution of calculating 𝒫_cct appears to be problematic. The adopted solution consists of an iterative procedure as used in MNH for non-flat problems. The non-respect of the Ψ^∗ condition in 𝒫 leads to a not well-posed system, and the iterative procedure spreads to the entire fluid domain the enforcement of the null divergence imposed on solid cells. The resolution of the pseudo-Poisson Eq. (5) leads to ${\mathrm{\Psi }}^{\ast }\to {\mathrm{\Psi }}^{*M}=\stackrel{M}{\sum _{m=\mathrm{1}}}{\mathcal{P}}^{-\mathrm{1}}.{\mathcal{Q}}_{\mathrm{cct}}^m$, where M is the number of iterations. This number is limited by a convergence criterion (compromise between incompressibility satisfaction and CPU cost). Many iterative procedures are available in MNH originally developed for non-Cartesian grids. Richardson and preconditioned conjugate residual algorithms have been adapted here to the obstacle immersion. The newly modified pressure solver is tested and validated in Sect. 4.1.

3.3 Consistency with the turbulence scheme

It is known that $\frac{l_m}{l_{\mathit{ϵ}}}\to \mathrm{1}$ is a reasonable approximation in nonhomogeneous, non-isotropic turbulence such as in the near-wall region. This approximation is indeed retained in the present IBM implementation, which assumes l_m=l_ϵ (hereafter noted l_m and called the mixing length). The Redelsperger et al. (2001) corrections near the ground are to match the similarity laws and the free-stream model constants are not activated. l_m is equal to the numerical cutoff space scale sufficiently far from the ground, leading to a $\mathrm{\Delta }\sqrt{e}$ turbulent viscosity. Near the ground and following the Prandtl idea consisting of the assumption of the linear variation of l_m in the near-wall region, the upper limit of the mixing length is $min(kz,\mathrm{\Delta })$ (k is the von Kármán constant and z is the altitude).

The turbulent characteristics are highly affected by surface interaction. As a consequence and for LES, the subgrid turbulence scheme (Sect. 2.3) is modified in the presence of immersed obstacles in the subgrid turbulent kinetic energy equation (Eq. 6), mixing length computation and Reynolds stress diagnosis (Eqs. 7, 8 and 9).

The subgrid turbulent kinetic energy condition. The explicit-in-time resolution of Eq. (6) claims a GCT forcing and an interface condition on the STKE e. Commonly, the STKE profile is considered parabolically in the viscous sublayer (Craft et al., 2002; Bredberg, 2000) and constant in the inertial and wake outer layers (Kalitzin et al., 2005; Capizzano, 2011). Due to the high turbulent Reynolds number ${\mathit{\text{Re}}}_{\mathrm{t}}\approx \mathcal{O}({\mathrm{10}}^{\mathrm{4}}-{\mathrm{10}}^{\mathrm{5}})$ encountered, a homogeneous Neumann condition is applied at the immersed interface. The equilibrium between the production and dissipation of STKE could be discussed and controverted; this choice acts as a first stage in IBM development.

The near-wall correction of the mixing length. The von Kármán limitation due to immersed walls acts through the LSF, and the upper limit on the mixing length l_m near the interface becomes $min(kz,-\mathit{\varphi },\mathrm{\Delta })$ with a banning of negative values in the solid region. Whatever the production of STKE and the turbulent shear, the lower limit $l_m(-\mathit{\varphi }\le \mathrm{0})$ induces a null value of the diagnosed surface fluxes. In addition, a singularity appears in the dissipative term ${\mathit{\rho }}_{\mathrm{r}}{K}_{\mathit{ϵ}}{e}^{\mathrm{3}/\mathrm{2}}{l}_m^{-\mathrm{1}}$. Through pragmatic reasoning, the singularity due to $l_m^{-\mathrm{1}}(\mathit{\varphi }\to {\mathrm{0}}^{-})\to \mathrm{\infty }$ amounts to the modeled length scales being smaller than the Kolmogorov scale (ν³ϵ⁻¹)^¼. Considering the Kolmogorov scale as modeled, the turbulence should vanish, which is in contradiction to the dissipative term. In order to overcome this ill-posed problem, a l_m lower limit has to be specified. In the study of atmospheric flows around buildings, a characteristic thickness of the viscous layer $H/\sqrt{\mathit{\text{Re}}}$ can be defined around an H bluff body for a Reynolds number based on the obstacle scale: H≈𝒪(10m); Re≈𝒪(10⁷). This thickness estimate is also proportional to $E\mathit{\nu }/u^{\ast }$ (E≈9.8 is commonly employed), whereby the friction velocity u^∗ is about 1 cm per second. Following these estimates, the length scale due to the viscous effects $z_{\mathit{\nu }}^{\mathrm{ib}}$ belongs to the millimeter domain in the expected atmospheric cases. Looking after a building surface and its large heterogeneity (door, windows, surface characteristics), its roughness length $z_{\mathrm{0}}^{\mathrm{ib}}$ is at least in the decimeter domain and $z_{\mathrm{0}}^{\mathrm{ib}}>z_{\mathit{\nu }}^{\mathrm{ib}}$ (Illustration in Fig. 7). For low Re and smooth surfaces, $z_{\mathit{\nu }}^{\mathrm{ib}}>z_{\mathrm{0}}^{\mathrm{ib}}$ could be encountered. Therefore, we assume $z^{\mathrm{ib}}=max(z_{\mathrm{0}}^{\mathrm{ib}},z_{\mathit{\nu }}^{\mathrm{ib}})$ and that z^ib is related to the size of smallest unresolved eddies near walls (i.e., dissipative scale). The mixing length near the wall is $z^{\mathrm{ib}}<l_m<min(\mathit{\text{k}}z,-\mathit{\varphi },\mathrm{\Delta })$.

Figure 7Illustration of the unresolved physical processes near a nonidealized solid wall (black line) in an atmospheric context: the length scale based on the viscous effects (grey line) is drastically smaller than the roughness length. The roughness length approaches the scale of smallest eddies and governs the log-law profile.

Download

The turbulent fluxes correction. The $\stackrel{\mathrm{‾}}{\mathit{\psi }}$ gradient and the turbulent diffusion $\mathcal{O}\left(z^{\mathrm{ib}}\sqrt{e}\right)$ prescribe the turbulent fluxes at the immersed interfaces (Eqs. 7, 8 and 9). As a first step in the MNH–IBM implementation, a no-flux condition on the mean potential temperature is imposed, leading to a zero value of the sensible heat flux. Writing the mean velocity field at B as $\stackrel{\mathrm{‾}}{\mathit{u}}=\stackrel{\mathrm{‾}}{u_t}\mathbf{t}$, $\stackrel{\mathrm{‾}}{u_t}\left(B\right)$ is needed to recover a gradient consistent with the turbulent shear. Considering the Prandtl (1925) or von Kármán (1930) theories, the logarithmic profile is assumed in the vicinity of the wall according to ${\stackrel{\mathrm{‾}}{u}}_t\left(z\right)=\frac{u^{\ast }}{k}\ln \left(\mathrm{1}+\frac{z}{z^{\mathrm{ib}}}\right)$. Considering Δ to be the limit of the resolved scales, most of the turbulent kinetic energy $\frac{\mathrm{1}}{\mathrm{2}}(\stackrel{\mathrm{‾}}{{u^{\prime }}^{\mathrm{2}}}+\stackrel{\mathrm{‾}}{{v^{\prime }}^{\mathrm{2}}}+\stackrel{\mathrm{‾}}{{w^{\prime }}^{\mathrm{2}}})$ is contained in the subgrid when $-\mathit{\varphi }<\mathrm{\Delta }$ and as $K_{\mathrm{tke}}\sqrt{e}$ with a constant K_tke≳1. This assumption is reinforced by the homogeneous Neumann condition applied on e. This approach derives from RANS (Reynolds-averaged Navier–Stokes) approaches, and the velocity friction is formulated as $u^{\ast }=K_{\mathrm{tke}}\sqrt[\mathrm{4}]{C_{\mathit{\mu }}}\sqrt{e}$, where C_μ is a constant evolving between 0.03 (atmospheric applications) and 0.09 (fluid mechanics applications). Adding a damping function for the viscous cases (low turbulent Reynolds number, Re_t<20), the tangent wall velocity at the interface is written as

Finally, the pragmatic limitation ${\stackrel{\mathrm{‾}}{u}}_t\left(B\right)\le \stackrel{\mathrm{‾}}{u_t}(\mathit{\varphi }=-\mathrm{\Delta }/\mathrm{2})$ operates if the STKE value is too high. The proposed dynamic wall model evolves between the no-slip and free-slip conditions. If the subgrid turbulence is weak or if the physical problem is fully resolved, the viscous layer is well modeled and $\stackrel{\mathrm{‾}}{u_t}\left(B\right)\to \mathrm{0}$. Otherwise, for an intense subgrid turbulence or a fully unresolved problem, the shear due to the wall presence is not perceived and $\frac{\mathrm{d}\stackrel{\mathrm{‾}}{u_t}}{\mathrm{d}n}{|}_B\to \mathrm{0}$. The wall model establishes an equilibrium between the production of STKE and the mean parietal friction. Note that the use of a log-law model near a singularity such as sharp edges and corners could be called into question. Nevertheless, Section 5.1 and 5.2 show LES results employing this proposition. After numerical investigations done during the single-cube study, $K_{\mathrm{tke}}\sqrt[\mathrm{4}]{C_{\mathit{\mu }}}\approx \mathrm{1}$ appears to be a suitable choice.

4.1 Potential flow

Isolated from the rest of the code, the resolution of the pseudo-Poisson Eq. (5) leads to potential solutions (Sect. 2.2). Theoretical ones are available for flow developed around a nondeformable obstacle such as an infinite cylinder or a sphere (Milne-Thomson, 1968; Batchelor, 2000). The two bodies are investigated here. The flow around the infinite cylinder is predominantly presented.

Figure 8Potential flow around a cylinder: (a) initial state around the body of diameter D_cyl=2R_cyl; (b) streamlines obtained after the Poisson equation resolution. The confinement is defined as $L=L_{\mathrm{cyl}}/R_{\mathrm{cyl}}$).

Download

Figure 8 illustrates the cylinder case. The fluid density is considered constant in time and in space. The flow is initially imposed as spatially homogeneous with a constant module of velocity U_∞ and parallel streamlines (Fig. 8a). This initialization does not respect the conservation of the momentum flux, and the irrotational correction of the projection method goes to recover this conservation. At the same time, the impermeability of the cylinder of diameter D_cyl=2R_cyl is achieved. Figure 8b shows the streamlines obtained with the MNH pressure solver modified to take into account the presence of immersed obstacles (Sect. 3.2). Defining x as the direction parallel to the initial streamlines and y as the perpendicular one, the expected solution is $\mathit{u}\cdot U_{\mathrm{\infty }}^{-\mathrm{1}}=\cos \mathit{\alpha }(\mathrm{1}-\frac{R_{\mathrm{cyl}}^{\mathrm{2}}}{r^{\mathrm{2}}})\mathbf{x}-\sin \mathit{\alpha }(\mathrm{1}+\frac{R_{\mathrm{cyl}}^{\mathrm{2}}}{r^{\mathrm{2}}})\mathbf{y}$ (single and non-confined body, (α;r) cylindrical coordinates). The numerical confinement is discussed hereafter, characterized by $L=L_{\mathrm{cyl}}/R_{\mathrm{cyl}}$, where L_cyl is the distance separating the lateral domain surfaces (Fig. 8a).

Figure 9Potential flow around a cylinder: (a) velocity convergence of two iterative methods (residual conjugate gradient, Richardson) for different spatial resolutions $N=[\mathrm{4}:\mathrm{32}]$(L=16); (b) evolution of the added mass coefficient A(N;L) with the confinement $L=L_{\mathrm{cyl}}/R_{\mathrm{cyl}}$ (N=16) and with the node number per radius cylinder N (L=16). The confinement is defined in Fig. 9a.

Download

The Richardson (RICH) and the residual conjugate gradient iterative (RESI) methods are tested (Sect. 3.2). Figure 9a plots the evolution of the dimensionless residue R(k) (based on a characteristic divergence U_∞∕Δ) with the iteration number k and obtained with the confinement L=16. The two algorithms converge with a weak dependence on the spatial discretizations ($N=\left[\mathrm{4}\right(\mathrm{red});\mathrm{8}(\mathrm{green});\mathrm{16}(\mathrm{blue});\mathrm{32}(\mathrm{purple}\left)\right]$ nodes per R_cyl). The slope coefficient $\frac{\mathrm{d}R\left(k\right)}{\mathrm{d}k}\left(\mathrm{RESI}\right)\mathit{\lesssim }\mathrm{3}\frac{\mathrm{d}R\left(k\right)}{\mathrm{d}k}\left(\mathrm{RICH}\right)$, so RESI demonstrates the highest velocity convergence. Even if RICH is about 20 % faster per iteration than RESI, the global CPU cost of the last one is lowest for the same solver residue R(k). For this reason and due to a higher a priori radius convergence, RESI is adopted. Note that the momentum flux computed after the solver convergence at the x_cyl location (Fig. 8a) shows good mass preservation with a relative error of $\left[\mathrm{0.48}\mathit{\%}\right(N=\mathrm{4});\mathrm{0.20}\mathit{\%}(N=\mathrm{8});\mathrm{0.18}\mathit{\%}(N=\mathrm{16});\mathrm{0.14}\mathit{\%}(N=\mathrm{32}\left)\right]$ in regard to the incoming flux localized by its x_inlet longitudinal coordinate. Similar results were obtained with a spherical body (not shown here).

With a change of Galilean reference frame, this study corresponds to a uniform body acceleration a_b in a fluid initially at rest. However, a possible viscous term, the hydrodynamic force exerted on the body, is reduced to the added mass effect $\mathcal{A}{m}_{\mathrm{f}}{\mathbf{a}}_{\mathrm{b}}={\int }_{{\mathcal{V}}_{\mathrm{s}}}\frac{\partial {\mathit{\rho }}_{\mathrm{f}}\mathit{u}}{\partial t}\mathrm{d}\mathcal{V}$ for Δt→0. 𝒜 is the dimensionless coefficient and m_f the displaced fluid mass. 𝒜_cyl theoretically equals 1 in the non-confined cylinder case (Lamb, 1932). The red curve in Fig. 9b illustrates the effect of the confinement L on 𝒜_cyl for N=16 resolution. Unsurprisingly, 𝒜_cyl increases with the confinement (Brennen, 1982). The weak dependence of 𝒜_cyl with L>16 allows us to consider the body to be isolated for L∼16. The green curve in Fig. 9b shows the impact of the space resolution for L=16. The numerical added mass coefficient is in good agreement with the theoretical one, presenting a relative error of about 2 % for N>16. It induces a respect for the impermeability hypothesis at the immersed interface. A similar study for a spherical body gives ${\mathcal{A}}_{\mathrm{sph}}=\frac{\mathrm{1}}{\mathrm{2}}+\mathrm{0.4}\mathit{\%}$. Figure 10 illustrates the contours of the kinetic energy around the sphere in an arbitrary symmetry plane. The green contours (numerical solutions) fit well with the red contours (theoretical solutions). A convergence study of the pressure solver is discussed in Appendix A.

Figure 10Potential solution around a sphere: (a) kinetic energy in an arbitrary symmetry plane; (b) zoom.

Download

4.2 Viscous flow

A pure dynamic and well-documented case that naturally follows previous ones is studied here. This physical case is the wake past a circular cylinder (non-stratified flow) at two moderate Reynolds numbers $\mathit{\text{Re}}=(\mathrm{40};\mathrm{140})$. One of the forerunners is Taneda (1956), who experimentally studied the nature of eddy structures.

Figure 11Eddy structure in a viscous fluid: steady (left, Re≈40) and unsteady (right, Re≈140) solutions obtained by the current numerical investigation (a and b, MNH–IBM and 30pts∕D_cyl) and by the Taneda (1956) experiments (c and d). The visualization is due to the presence of a passive tracer injected on the body surface and transported by the flow.

Taneda (1956) found a regular Hopf bifurcation at a critical Reynolds number ${\mathit{\text{Re}}}_{\mathrm{c}}=\frac{U_{\mathrm{\infty }}{D}_{\mathrm{cyl}}}{{\mathit{\nu }}_{\mathrm{f}}}\approx \mathrm{45}$. Below Re_c and above Re>5, a boundary layer separation leads to a steady recirculating region in the near wake (Fig. 11c). Above Re_c≈45, an unsteady mode breaks the planar symmetry and the body wake presents an alternate vortex shedding (Fig. 11-d). The standing eddy (the von Kármán street) obtained by MNH–IBM at Re=40 (140) is visualized in Fig. 11a (b) by the injection of a passive tracer on the body surface.

The standing eddies at Re<Re_c are commonly described with a θ_d detachment angle, l_r recirculating length and (a;b) location of the vortex core (Fig. 12). The limit of the numerical domain is 10D_cyl upstream of the obstacle for the inlet condition (U_∞, the uniform incoming velocity) and lateral condition (slip condition) and 15D_cyl for the outlet condition, allowing for the vorticity evacuation. As Cai et al. (2017) mention, this domain can induce a low numerical confinement effect. Three regular Cartesian meshes are built with 10, 20 and 30 nodes per D_cyl.

Figure 12Recirculating region at Re=40 (MNH–IBM, 20pts∕D_cyl): definition of the θ_d (^∘) separation angle, l_r recirculating length and (a;b) vortex core location. The distance is dimensionless by D_cyl.

Download

The 20pts∕D_cyl and 30pts∕D_cyl meshes present a good spatial convergence and weak differences at Re=40. ${\mathit{\theta }}_{\mathrm{d}}\sim \mathrm{53}{}^{\circ }\pm \mathrm{2}{}^{\circ }$ and the recirculation length $l_{\mathrm{r}}/D_{\mathrm{cyl}}\sim \mathrm{2.2}\pm \mathrm{0.05}$. The 10pts∕D_cyl mesh shows more discrepancies, which are attributable to the nonability of the coarsest resolution to capture the viscous boundary layer for which the thickness evolves in $D_{\mathrm{cyl}}\sqrt[-\mathrm{1}]{\mathit{\text{Re}}}$. Note that the impact of the low-order (centered or WENO) modeling of the advection at the immersed interface is weak for this viscous case (Sect. 3.1). Table 1 compares the 20pts∕D_cyl results with a part of the results literature collected in Gautier et al. (2013).

Table 1Description of the standing eddies in the wake of the solid cylinder (Re=40): comparison of the separation angle θ_d (^∘), recirculating length l_r (m) and vortex core location (a;b) between the literature and MNH–IBM.

Download Print Version | Download XLSX

The focus is on the unsteady mode at Re=140. The ratio between the characteristic time of inertial effects D_cyl∕U_∞ and the one related to the vortex shedding 1∕f defines the Strouhal number $\mathit{\text{St}}=\frac{fD}{U_{\mathrm{\infty }}}$. Brazza et al. (1986), Park et al. (1998) and Stålberg et al. (2006) confirm the equation $\mathit{\text{St}}\left(\mathit{\text{Re}}\right)=-\mathrm{3.3265}/\mathit{\text{Re}}+\mathrm{0.1816}+\mathrm{1.6}\cdot {\mathrm{10}}^{-\mathrm{4}}/\mathit{\text{Re}}$ proposed by Williamson (1989). MNH–IBM obtains $\mathit{\text{St}}(\mathit{\text{Re}}=\mathrm{140})\in [\mathrm{0.177}:\mathrm{0.179}]$ and an absolute maximum relative error lower than 2 % in regard to the Williamson (1989) formulation with the two finer resolutions. Our results are in good agreement with those presented in the abovementioned and more extensive studies. Details on DNS validation in a viscous buoyancy-driven flow are also presented in the Supplement.

5.1 Flow over a surface-mounted cube

Using static pressure measurements, as well as laser-sheet and oil-film visualizations, Martinuzzi and Tropea (1993) and Hussein and Martinuzzi (1996) generated a large data set for the study of flows around a cubic body placed in a channel (Fig. 13a). RANS and LES have been used to explore in detail this physical case (Breuer et al., 1996; Shah and Ferziger, 1997; Rodi et al., 1997; Frank, 1999; Krajnovic and Davidson, 2002; Farhadi and Rahnama, 2006).

Figure 13(a) Top visualization of the flow around a cube (Hussein and Martinuzzi, 1996); (b) schematic representation of the time-averaged vortex structure around a cube (Martinuzzi and Tropea, 1993) and index of the recirculating regions.

Physical details. A cube (H side) is placed in a channel of 2H height. The channel is sufficiently large in the spanwise direction to consider the cube to be single in that direction. Turbulent flow is generated in the channel upstream of the cube with a mean bulk velocity U_b. Defining the dimensionless wall coordinate $z^{+}=u^{\ast }\cdot z/{\mathit{\nu }}_{\mathrm{f}}$, the stream-wise upstream velocity corresponds to a log law for smooth walls $\stackrel{\mathrm{‾}}{u}\left(z^{+}\right)\cdot u^{\ast -\mathrm{1}}=\mathrm{5.54}+\frac{\mathrm{1}}{\mathit{\kappa }}\log \left(z^{+}\right)$ as described in Hussein and Martinuzzi (1996). The Reynolds number, as defined by the mean bulk velocity, the cube height and the molecular diffusion, is Re≈40 000.

The mean flow around the cube presents a set of five recirculating regions (Fig. 13b). Each cube surface is associated with one of these regions: the A–B vortex separations in front of the cube, which spread laterally in a horseshoe D, two vortices near side walls E, one F on the roof and main arch vortex G downstream.

Numerical details. The top and bottom surfaces of the cubic body are modeled by the IBM. A small value of the roughness length $z_{\mathrm{0}}/H\approx {\mathrm{10}}^{-\mathrm{6}}$ is imposed (low value to model a smooth interface, viscous-scale intervention in the z⁺ calculation). The stream-wise (spanwise) direction is x (y). The size of the grid is set as $(x,y,z)=(-\mathrm{24}H:\mathrm{8}H,-\mathrm{4}H:\mathrm{4}H,\mathrm{0}H:\mathrm{2}H)$ with a location of the cube center at $(\frac{H}{\mathrm{2}};\frac{H}{\mathrm{2}};\frac{H}{\mathrm{2}})$. Three regular Cartesian meshes are employed with a respective space step $H/\mathrm{\Delta }=[\mathrm{10};\mathrm{20};\mathrm{40}]$. $x/H\in [-\mathrm{24}:-\mathrm{4}]$ is a region employed to model the fully turbulent character of the incoming flow. The incoming turbulent state is obtained by the IBM and a pseudo-recycling method inspired from the works of Lund (1998), Mayor et al. (2002), and Yang and Meneveau (2016) (not detailed here). The vertical profiles of the stream-wise velocity and turbulent intensity $\sqrt{\stackrel{\mathrm{‾}}{u^{\prime \mathrm{2}}}}/U_{\mathrm{b}}^{\mathrm{2}}\approx {\mathrm{2.10}}^{-\mathrm{2}}$ (Hussein and Martinuzzi, 1996) are recovered at $x/H\sim -\mathrm{4}$. We note that the turbulence generation deserves more attention, but we prefer to concentrate only on the cube wake.

Figure 14Vertical symmetry plane of the mean flow: (a–c) MNH–IBM time-averaged streamlines; (d) streamlines observed by Martinuzzi and Tropea (1993); (e) MNH–IBM instantaneous visualization of the Q criterion (Hunt et al., 1988).

Results. Figure 14a, b and c show the time-averaged streamlines in the vertical symmetry plane of the cube obtained by MNH–IBM for the three space resolutions. The streamlines of the coarse, medium and fine resolution are respectively in red, green and blue. The discretization order of the fine resolution is close to that of most literature LESs except Shah and Ferziger (1997), who used a far more precise grid near walls. The same figure obtained by the experimental investigation is given in Fig. 14d. The size of the front (rear) region is characterized by the recirculating length x_f∕H (x_r∕H). The experiment gives $x_{\mathrm{f}}/H\approx [\mathrm{1.04}:\mathrm{1.05}]$ and $x_{\mathrm{r}}/H\approx [\mathrm{1.64}:\mathrm{1.67}]$. The LES reference results give the ranges $x_{\mathrm{f}}/H\in [\mathrm{0.81}:\mathrm{1.28}]$ and $x_{\mathrm{r}}/H\in [\mathrm{1.38}:\mathrm{2.25}]$. For the two finest resolutions (green and blue), the overall prediction of MNH–IBM recovers a consistent mean topological structure. MNH–IBM obtains for the two finest resolutions $x_{\mathrm{f}}/H\in [\mathrm{0.99}:\mathrm{1.21}]$ and $x_{\mathrm{r}}/H\in [\mathrm{1.48}:\mathrm{1.55}]$. MNH–IBM does not capture, as with most LESs, the two dividing lines A∕B mentioned by Martinuzzi and Tropea (1993) but only captures a flattened vortex. However, the bifurcation point near the rear edge and ground is not detected by MNH–IBM, while this point was modeled in Rodi et al. (1997). This bifurcation point was also commented on in Martinuzzi and Tropea (1993), even if their experimental uncertainty did not allow us to visualize it in Fig. 14d.

Figure 14e shows an MNH–IBM instantaneous flow field with the Q criterion (Hunt et al., 1988), and as the LES reference results mention, it presents a highly intermittent character clearly visible with the quasi-disappearance of the D horseshoe and G arch. A frequency f of vortex shedding dominates the highly intermittent activity in the body wake, leading to the experimental Strouhal number $\mathit{\text{St}}=\frac{f\cdot H}{U_{\mathrm{b}}}\approx \mathrm{0.145}$. An MNH–IBM energy spectrum (discrete Fourier transform of w(t) in the body wake) finds a peak $\mathit{\text{St}}\in [\mathrm{0.10}:\mathrm{0.12}]$ for all the studied resolutions and a $\sim -\mathrm{5}/\mathrm{3}$ energetic cascade slope for larger wave numbers (not shown). The St values obtained by MNH–IBM are slightly lower than the experimental one but stay consistent with the $\mathit{\text{St}}\in [\mathrm{0.10}:\mathrm{0.15}]$ range of other LESs.

Figure 15Mean vertical profiles of velocities (top), turbulent kinetic energy and Reynolds stress (bottom). The lines correspond to the MNH–IBM results. The symbols are the Martinuzzi and Tropea (1993) data except for TKE (Hussein and Martinuzzi, 1996). The profiles are given at four longitudinal locations: (a, e, i, m) $x/H=\mathrm{1}/\mathrm{2}$; (b, f, j, n) $x/H=\mathrm{1}$; (c, g, k, o) $x/H=\mathrm{2}$; (d, h, l, p) $x/H=\mathrm{4}$.

Download

Still in the vertical symmetry plane of the mean flow, Fig. 15 plots at four longitudinal positions $x/H=(\mathrm{1}/\mathrm{2};\mathrm{1};\mathrm{2};\mathrm{4})$ the vertical profiles of time-averaged quantities related to the steady ($\stackrel{\mathrm{‾}}{u}$, $\stackrel{\mathrm{‾}}{w}$) and unsteady (TKE, $\stackrel{\mathrm{‾}}{u^{\prime }{w}^{\prime }}$) parts of the solution. The color code corresponds to the spatial resolution (10pts∕H in red; 20pts∕H in green; 40pts∕H in blue). Note that the U_b bulk velocity was set to unity and therefore presents variables in Fig. 15 that are dimensionless. In most of the sub-figures, the higher the space resolution, the lower the gap with the experiment. The top of Fig. 15 leads to a similar conclusion as the literature LESs: the stream-wise velocity is well recovered. The counterflow at the rear and roof of the cube near $(x/H;z/H)\approx (\mathrm{1};\mathrm{1})$ informs us about the existence of the bifurcation point S1 (Fig. 15b). An underestimation appears on the counterflow at $(x/H;z/H)\approx (\mathrm{2};\mathrm{0})$ and is frequently observed in the literature (Fig. 15c). Some discrepancies are found on two $\stackrel{\mathrm{‾}}{w}$ profiles (Fig. 15f–g). Note that Shah and Ferziger (1997) and Krajnovic and Davidson (2002) highlight the difficulty of recovering it. The turbulent kinetic energy and the Reynolds stress are correctly predicted at the cube roof (Fig. 15i–j). The vertical profiles of TKE and $\stackrel{\mathrm{‾}}{u^{\prime }{w}^{\prime }}$ downstream of the cube show an overestimate tendency (Fig. 15k, p, l). No experimental data are available on TKE at $x/H=\mathrm{4}$ (Fig. 15l), but using the $\stackrel{\mathrm{‾}}{u^{\prime \mathrm{2}}}(x/H=\mathrm{4})$ experimental value (not shown here) and the $\stackrel{\mathrm{‾}}{u^{\prime \mathrm{2}}}/\mathrm{TKE}$ ratio obtained by MNH–IBM, we suspect that TKE$(x/H=\mathrm{4})$ is still overestimated. This turbulence diagnosis is relatively similar to that of Farhadi and Rahnama (2006).

Sensitivity study. Some tests have shown a sensitivity of the solution with both the incoming turbulence and constants in the immersed wall model. Indeed, the existence of the small vortex near the saddle node S1 (Depardon et al., 2006) turns out to be strongly dependent on the inlet turbulence (the lower the turbulent intensity, the bigger the size of this vortex); the reattachment (or not) of the roof vortex F with the main arch G is also observed. This sensitivity follows the observations of Castro and Robins (1977), who studied the cube placed in a uniform or turbulent incoming flow. In the same way the existence of the vortex near the saddle node S2 (Frank, 1999) is dependent on the ground boundary condition. The $u^{\ast }\sqrt[-\mathrm{1}]{e}=K_{\mathrm{tke}}{\sqrt[\mathrm{1}/\mathrm{4}]{C}}_{\mathit{\mu }}$ ratio and the z₀ roughness length fixed in the immersed wall model (Sect. 3.3) impact the nature of the incoming turbulent state for which the surface shear plays an important role. To give a significant example and if a null value of K_tke is applied on the channel surfaces (non-slip condition), x_f highly increases and the vortex at S2 appears. Otherwise, K_tke and z₀ do not crucially affect the nature of the dynamic in the vicinity of the cube surfaces; the pressure gradient between front and back faces governs.

To conclude this section and despite the uncertainties of the inlet condition, MNH–IBM is in good agreement with the experiments of Martinuzzi and Tropea (1993), Hussein and Martinuzzi (1996), and other LESs using a resolution of about 10 points per cube length. Even if the coarsest resolution loses a part of the expected physics, it maintains a suitable modeling of the largest structures of the flow. A parametric study on inlet turbulence generation has led us to fix K_tke≈2 in Eq. (24).

5.2 The Mock Urban Setting Test (MUST) experiment

The MUST is an experimental campaign organized during early autumn 2001 in Utah's West Desert (Biltoft, 2001; Biltoft et al., 2002). Its objective was to quantify the dispersion of a passive tracer (propylene) in a dry atmospheric context over a topography reproducing a near-urban canopy. The main interest lies in the similitude between this experiment and a pollution episode due to toxic gas propagation over a city with high population densities. It provides extensive measurements of meteorological variables and scalar dispersion information.

Figure 16(a) Photograph of the MUST containers in Utah's West Desert; (b) schematic representation of the container layout. The locations of the concentration (detectors) and wind sensors (S and T towers) are indicated, as are the position (red cross) of the pollutant release and the direction of the incoming flow (green arrow) for case 2681829. Source: Biltoft (2001) and Yee and Biltoft (2004).

Physical details. The near-regular array is composed of 120 containers. Figure 16 gives a photograph and a picture of the topography. The containers are equivalent in volume and shape. Their spatial dimensions are $(L_x,L_y,L_z)=\sim (\mathrm{2.4},\mathrm{12.2},\mathrm{2.5})$ m and the horizontal distance between containers is 𝒪(10) m. Following Table II in Yee and Biltoft (2004), the 2681829 case (25 September 2001 at 18:29 UTC) is selected. The Monin–Obukhov length is L_mo≈28 000 m and the stability condition is supposed neutral; the buoyancy effects and the sensible heat flux are negligible in regard to the inertia effects and the turbulent shear. The incoming flow shows a mean horizontal angle with the container layout (green arrow, Fig. 16b). The MNH–IBM results are compared to the experimental measurements reachable at several altitudes (4, 8, 16 m) at the south tower and at the main T tower placed in the array center. The locations of the towers are indicated in Fig. 16b. A roughness length z₀=0.045 m is given by the experiment and related to the surrounding desert vegetation.

Numerical details. The externalized scheme SURFEX (Masson et al., 2013) models the ground friction. LSF is generated to represent the topography and IBM is used to model the containers. The smallest characteristic container length is discretized by 𝒪(10) cells. The mesh is a Cartesian and regular grid for z<10 m with a space step Δ=0.2 m. Above 10 m, the vertical space step is released in a geometric progression of 1.08 ratio. The altitude of the numerical domain top is about 8 times the container height.

Figure 17Mean vertical profiles in the MUST experiment: (a) $\left|\stackrel{\mathrm{‾}}{\mathit{u}}\right|$ horizontal velocity; (b) $\mathrm{atan}(\stackrel{\mathrm{‾}}{v}/\stackrel{\mathrm{‾}}{u})$ wind direction at the S (blue) and T (black) towers. The symbols (lines) are the Yee and Biltoft (2004) experimental measurements (MNH–IBM results).

Download

The distance between the horizontal limit of the computational domain and the array is about 20 times the container height. The large-scale flow is forced by an open boundary condition. A mean horizontal angle of −41^∘ with the x direction is fixed for all altitudes; note that the low angle deviation and the turbulence observed upstream of the containers in the experiment are not numerically considered. Following the experimental data given at the S tower (Fig. 17a, blue symbols) and assuming a log law $\left|\stackrel{\mathrm{‾}}{\mathit{u}}\right|=\frac{u^{\ast }}{\mathit{\kappa }}\log (\mathrm{1}.+z/z_{\mathrm{0}})$, a least-squares regression estimates a friction velocity $u^{\ast }=\mathrm{0.71}$ m s⁻¹ and allows us to build the vertical profile of the mean incoming flow (Fig. 17a, blue line). This u^∗ value is close to the experimental one $u_{\exp }^{\ast }=\mathrm{0.68}$ m s⁻¹ found by a sonic anemometer at the feet of the S tower. The same inlet condition is used in the numerical studies of Hanna et al. (2004), Milliez and Carissimo (2007), and Donnelly et al. (2009). Note that an additional term $\mathcal{O}\left(L_{\mathrm{mo}}^{-\mathrm{1}}\right)$ appears in their formulation but is negligible in the 2681829 case.

Results. The black line (MNH–IBM) and symbols (Yee and Biltoft, 2004) in Figure 17a (b) show the impact of the near-urban canopy on the vertical profile of $\left|\stackrel{\mathrm{‾}}{\mathit{u}}\right|$ (wind angle) in the core of the array (T tower). Not surprisingly, the canopy induces a global slowdown of $\left|\stackrel{\mathrm{‾}}{\mathit{u}}\right|$ near the ground and up to z≲8 m. A decrease in the mean horizontal wind angle is found at the same altitudes. This deviation is related to the container orientation, which tends toward the flow being aligned with the y direction. The same pattern is discussed in Milliez and Carissimo (2007). A wind acceleration and an increase in the wind angle are observed by MNH–IBM for z≳8 m as in the LES results of König (2014) and Dejoan et al. (2010) on similar MUST cases. These few degrees of deviation are observed by the experiment but not the acceleration. A part of this acceleration may be explained by a Venturi effect and a too-closed top limit of the computational domain (numerical confinement, under investigation).

Figure 18Wind at the horizontal cut at z=1.6 m: (a) time-averaged on 200 s; (b) instantaneous at t=200 s. The black squares indicate the location of the T and S towers (MNH–IBM results).

Download

Figure 18a (b) illustrates the time-averaged horizontal wind field $\left|\stackrel{\mathrm{‾}}{\mathit{u}}\right|$ ($\left|\mathit{u}\right|$ at t=200 s) at the 1.6 m altitude. These figures highlight the fact that the incoming flow is not turbulent in the simulation. This assumed gap constitutes a perspective (Camelli et al., 2005) not directly linked to IBM. It also allows us to introduce here some comments on the turbulence state. The atmospheric turbulence is dependent on the roughness length (z₀ considered constant due to the homogeneous and flat desert over a few miles upstream). The first container rows are the scene of the boundary layer transition and act as a region of strong roughness change. The turbulence observed in the urban-like canopy has two origins: the incoming turbulence and that induced by the container presence. The contribution of both turbulence types varies as a function of the altitude and distance to the first row.

Table 2Kinetic energy, turbulent kinetic energy and friction velocity obtained by the experimental (Yee and Biltoft, 2004) and numerical (MNH–IBM) investigations at three altitudes of the tower T.

Download Print Version | Download XLSX

The mean kinetic energy $E_{\mathrm{k}}=\frac{\mathrm{1}}{\mathrm{2}}({\stackrel{\mathrm{‾}}{u}}^{\mathrm{2}}+{\stackrel{\mathrm{‾}}{v}}^{\mathrm{2}}+{\stackrel{\mathrm{‾}}{w}}^{\mathrm{2}})$, friction velocity $u^{\ast }=\sqrt[\mathrm{4}]{{\stackrel{\mathrm{‾}}{u^{\prime }{w}^{\prime }}}^{\mathrm{2}}+{\stackrel{\mathrm{‾}}{v^{\prime }{w}^{\prime }}}^{\mathrm{2}}}$ and turbulent kinetic energy $\mathrm{TKE}=\frac{\mathrm{1}}{\mathrm{2}}(\stackrel{\mathrm{‾}}{{u^{\prime }}^{\mathrm{2}}}+\stackrel{\mathrm{‾}}{{v^{\prime }}^{\mathrm{2}}}+\stackrel{\mathrm{‾}}{{w^{\prime }}^{\mathrm{2}}})$ are estimated at the T tower and indicated in Table 2. All the variables are in good agreement with the experiment at z=4 m (about 2 times the container height). The friction velocity at z(T)=4 m is at least 2 times larger than $u_{\exp }^{\ast }\left(S\right)=\mathrm{0.68}$ m s⁻¹ observed at the S tower feet. That increase is the signature of the turbulence developed by the urban-like canopy. Looking at the results at z(T)=8 m and z(T)=16 m, the higher the altitude, the more discrepancies appear between the experimental and numerical results. The experimental measurement of the friction velocity at z=16 m is close to the upstream value.