Implementation of an immersed boundary method in the Meso-NH v5.2 model: applications to an idealized urban environment

Abstract

This study describes the numerical implementation, verification and validation of an immersed boundary method (IBM) in the atmospheric solver Meso-NH for applications to urban flow modeling. The IBM represents the fluid–solid interface by means of a level-set function and models the obstacles as part of the resolved scales.

The IBM is implemented by means of a three-step procedure: first, an explicit-in-time forcing is developed based on a novel ghost-cell technique that uses multiple image points instead of the classical single mirror point. The second step consists of an implicit step projection whereby the right-hand side of the Poisson equation is modified by means of a cut-cell technique to satisfy the incompressibility constraint. The condition of non-permeability is achieved at the embedded fluid–solid interface by an iterative procedure applied on the modified Poisson equation. In the final step, the turbulent fluxes and the wall model used for large-eddy simulations (LESs) are corrected, and a wall model is proposed to ensure consistency of the subgrid scales with the IBM treatment.

In the second of part of the paper, the IBM is verified and validated for several analytical and benchmark test cases of flows around single bluff bodies with an increasing level of complexity. The analysis showed that the Meso-NH model (MNH) with IBM reproduces the expected physical features of the flow, which are also found in the atmosphere at much larger scales. Finally, the IBM is validated in the LES mode against the Mock Urban Setting Test (MUST) field experiment, which is characterized by strong roughness caused by the presence of a set of obstacles placed in the atmospheric boundary layer in nearly neutral stability conditions. The Meso-NH IBM–LES reproduces with reasonable accuracy both the mean flow and turbulent fluctuations observed in this idealized urban environment.

1 Introduction

Urbanization impacts the physical and dynamical structure of the atmospheric boundary layer, influencing both the local weather and the concentration and residence time of pollutants in the atmosphere, which in turn impact air quality. While the physical mechanisms driving these interactions and their connections to climate change are well understood (the urban heat island effect, anthropological effects), their precise quantification remains a major modeling challenge. Accurate predictions of these interactions require modeling and simulating the underlying fluid mechanics processes to resolve the complex terrain featured in large urban areas, including buildings of different sizes, street canyons and parks. For example, it is well known that pollution originates from traffic and industry in and around cities, but the actual dispersion mechanisms are driven by the local weather. Furthermore, fine-scale flow fluctuations influence nonlinear physicochemical processes. The present study addresses these issues by focusing on the numerical aspects of the problem.

With the progress in metrology, it is now possible to obtain reliable measurements of the atmospheric conditions over a city. For example, during the Joint Urban experiment (JU2003), scalar dispersion was measured experimentally over the streets of Oklahoma City (; ). Similarly, the CAPITOUL experiment (2004–2005), conducted in Toulouse, analyzed the turbulent boundary layer developed over the urban topography and evaluated the energy exchanges between the surface and atmosphere . More recently, the multiscale field study by provided meteorological observations and tracer concentration data in Salt Lake City. Other studies analyzed reduced-scale and/or idealized models to understand urban climate features as in the COSMO (Comprehensive Outdoor Scale Model Experiment for Urban Climate) and Kugahara projects . For example, and respectively used an array of cubic bodies and stone fields as small-scale models.

In order to use these experimental data in the future for model validation, the numerical models need to first be verified for academic test cases and simplified scenarios representative of atmospheric turbulent boundary layer flows. In particular, flow interaction with buildings or any generic obstacles plays a crucial role in urban flow modeling. The range of scales of objects acting as obstacles is huge in an urban setting, encompassing large buildings and small vegetation scales, and so is the range of the corresponding flow–obstacle interactions. Covering all possible cases is obviously impossible but we can rely on the invariance of certain flow characteristics at different scales. For example, the von Kármán streets are observed in the wake of a centimeter-scale cylinder as well as in the cloud layout behind an island. Following this, a wide selection of benchmark flows can be analyzed to verify and validate the numerical treatment of fluid–obstacle interaction with a view to atmospheric applications.

The physical application considered in this work is the atmospheric mesoscale reaction to perturbations induced by urban areas; the more the obstacles are considered part of the scales numerically resolved, the higher the accuracy of the results. To access this resolution, this study presents the development, implementation, verification and validation of an immersed boundary method (IBM) in the Meso-NH model (MNH) for applications to urban flow modeling

Meso-NH scientific documentation: http://mesonh.aero.obs-mip.fr/mesonh52/BooksAndGuides (last access: 13 May 2019).

. This choice was dictated by the fact that numerical solvers in MNH enforce conservation on structured grids and hence cannot handle body-fitted grids with steep topological gradients. The main idea behind IBM is the detection of an interface separating a fluid region (where conservation laws hold) from a solid region (corresponding to the obstacle volume) using different techniques (e.g., markers, level-set functions, local volume fraction, etc). Two main classes of IBM exist based on the continuous and discrete forcing approaches, respectively. The continuous forcing approach was developed by for biomechanics applications and consists of the addition of a continuous artificial force (acceleration indeed) in the momentum conservation equation that mimics the effect of the obstacles (heart linings) and drives the flow to relax to no-slip conditions at the wall of the obstacles. This approach and its variant developed by for a rigid interface can suffer from the lack of stiffness (fluid–solid interface is generally spread over few cells) and the time step restriction (spring and damper model with large natural frequency). Nevertheless, the continuous forcing approaches are very successful in many applications (penalization method as in , fictitious domain method, etc). In the second IBM class, the discrete approach, the boundary conditions are specified at the immersed interface. To simulate flows around nonmoving and rigid bodies, two subclasses of discrete approaches can be defined as in : direct or indirect approaches, depending on the forcing location . Many types of discrete forcing exist, e.g., direct forcing in the fluid region near the interface as in , the immersed interface method and the Cartesian grid method . Depending on how to resolve the partial differential equations, Cartesian grid methods are written for finite-volume discretizations (cut-cell technique, CCT) and for finite-difference discretizations (ghost-cell technique, GCT) as in . CCT reshapes the cell cut by the interface to preserve mass, momentum and energy. Using GCT, the local spatial reconstruction is done inside the solid region. Note that the latter technique has been successfully implemented in the Weather Research and Forecasting (WRF) model .

In this study, a discrete forcing approach is adopted wherein the fluid–solid interface is modeled by means of a level-set function . The motivation behind this choice is that we are primarily interested in modeling explicitly rigid and nonmoving bodies in a turbulent flow and with sufficiently fine resolution to avoid the large dissipation inherent in the presumed spread interface. The GCT does not introduce source terms in the conservation equations modeling the fluid region so that boundary conditions are imposed at the interface and/or in the solid region. The only corrections to the physical model in the fluid region come from subgrid turbulent parameterizations. The idea is that in future mesoscale applications, IBM will be used to resolve large obstacles (in the solid region), such as buildings, but also mountains, whereas a subgrid drag model will be used to handle unresolved obstacles such as vegetation .

The paper is organized as follows. Section  briefly describes the general features of MNH. Section  details the numerical implementation of the IBM. Inspired by the works of , , and and to close the turbulence problem, an immersed wall model is proposed in Sect. . Section  and  describe the validation of the method for academic flows, respectively potential and viscous flows. Finally, Section  and  discuss the results of turbulent flow simulations and comparisons with data from field experiments. Conclusions are drawn in Sect. . Additional tests and validations on potential and inviscid flows are respectively documented in Appendix A and B. The study of a viscous and thermodynamic case is given in the Supplement.

2 The Meso-NH code at a glance

MNH is an atmospheric non-hydrostatic research model. Its spatiotemporal resolution ranges from the large meso-alpha scale (hundred of kilometers and days) down to the micro-scale (meters and seconds). It is massively parallel on the nested and structured grids adapted on most international hosting computer platforms. Several parameterizations are available: radiation, turbulence, microphysics, moist convection with phase change, chemical reactions, electric scheme and externalized surface scheme. In the present study, only two subgrid parameterizations are approached: turbulence and surface schemes.

2.1 The conservation laws

The spatial discretization $\mathbf{x}$ is based on terrain-following coordinates . 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 $\mathrm{\Delta }z=\mathrm{\Delta }$ 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 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 $\mathit{\psi }$ 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. ). 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 }$ (${\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. ). The potential temperature is defined through the Exner function $\mathrm{\Pi }$ 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 ${\mathit{\rho }}_{\mathrm{r}}$ is the density of the reference state, $P$ is the local pressure, $P_{\mathrm{r}\mathrm{0}}$ the reference ground-level pressure, $R_{\mathrm{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 $$\begin{array}{}\text{1}& {\displaystyle \frac{\partial \left({\mathit{\rho }}_{\mathrm{r}}\stackrel{\mathrm{‾}}{\mathit{\theta }}\right)}{\partial t}}{|}_{\mathrm{csl}}=-\mathrm{\nabla }\left({\mathit{\rho }}_{\mathrm{r}}\stackrel{\mathrm{‾}}{\mathit{\theta }}\stackrel{\mathrm{‾}}{\mathit{u}}\right)-\mathrm{\nabla }\cdot \left({\mathit{\rho }}_{\mathrm{r}}\stackrel{\mathrm{‾}}{{\mathit{\theta }}^{\prime }{\mathbf{u}}^{\prime }}\right)+{\mathit{\rho }}_{\mathrm{r}}{\stackrel{\mathrm{‾}}{F}}_{\mathit{\theta }}^{\mathrm{\Pi }},\text{2}& {\displaystyle \frac{\partial \left({\mathit{\rho }}_{\mathrm{r}}\stackrel{\mathrm{‾}}{s}\right)}{\partial t}}{|}_{\mathrm{csl}}=-\mathrm{\nabla }\cdot \left({\mathit{\rho }}_{\mathrm{r}}\stackrel{\mathrm{‾}}{s}\stackrel{\mathrm{‾}}{\mathit{u}}\right)-\mathrm{\nabla }\cdot \left({\mathit{\rho }}_{\mathrm{r}}\stackrel{\mathrm{‾}}{s^{\prime }{\mathbf{u}}^{\prime }}\right),\end{array}$$ where ${\stackrel{\mathrm{‾}}{F}}^{\mathrm{\Pi }}$ corresponds to pressure effects. The transport of each prognostic scalar in Eqs. (), () and () is made by a piecewise parabolic method (PPM) with undershoot and overshoot limitation . The temporal algorithm of the advection term in these scalar transports is a forward-in-time scheme (noted FT). The momentum equations are

3$$\begin{array}{rl}{\displaystyle \frac{\partial \left({\mathit{\rho }}_{\mathrm{r}}\stackrel{\mathrm{‾}}{\mathit{u}}\right)}{\partial t}}{|}_{\mathrm{csl}}=& -\mathrm{\nabla }\cdot \left({\mathit{\rho }}_{\mathrm{r}}\stackrel{\mathrm{‾}}{\mathit{u}}\otimes \stackrel{\mathrm{‾}}{\mathit{u}}\right)+\mathrm{\nabla }\cdot \left({\mathit{\mu }}_{\mathrm{f}}\mathrm{\nabla }\stackrel{\mathrm{‾}}{\mathbf{u}}\right)\\ & -\mathrm{\nabla }\cdot \left({\mathit{\rho }}_{\mathrm{r}}\stackrel{\mathrm{‾}}{{\mathbf{u}}^{\prime }\otimes {\mathbf{u}}^{\prime }}\right)+{\mathit{\rho }}_{\mathrm{r}}{\stackrel{\mathrm{‾}}{\mathbf{F}}}_u^{\mathrm{\Pi }}+{\mathit{\rho }}_{\mathrm{r}}\mathbf{g}{\displaystyle \frac{\stackrel{\mathrm{‾}}{\mathit{\theta }}-{\mathit{\theta }}_{\mathrm{r}}}{{\mathit{\theta }}_{\mathrm{r}}}},\end{array}$$ where $\stackrel{\mathrm{‾}}{\mathit{u}}$ is the resolved wind, $\mathbf{g}$ the acceleration due to the gravity appearing in the buoyancy term, ${\mathit{\theta }}_{\mathrm{r}}$ is the potential temperature of the reference state, ${\mathit{\mu }}_{\mathrm{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. () can be done by second- or fourth-order centered schemes and third- or fifth-order weighted essentially non-oscillatory schemes . The temporal evolution of the resolved wind is achieved by a fourth-order explicit Runge–Kutta (ERK) algorithm . In the present study $\mathrm{\Delta }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 $\mathcal{O}({\mathit{\nu }}_{\mathrm{f}}/{\mathrm{\Delta }}^{\mathrm{2}})$ (${\mathit{\nu }}_{\mathrm{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 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.  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 ${\mathrm{\Psi }}^{\ast }$:

4$${\stackrel{\mathrm{‾}}{\mathit{u}}}^{n+\mathrm{1}}={\stackrel{\mathrm{‾}}{\mathit{u}}}^{\ast }-{\displaystyle \frac{\mathrm{\Delta }t}{{\mathit{\rho }}_{\mathrm{r}}}}\mathbf{\nabla }{\mathrm{\Psi }}^{\ast }.$$

${\mathrm{\Psi }}^{\ast }$ is obtained with the resolution of the pseudo-Poisson equation written as 5$$\mathrm{\nabla }\cdot \left({\mathit{\rho }}_{\mathrm{r}}^{-\mathrm{1}}\mathbf{\nabla }{\mathrm{\Psi }}^{\ast }\right)=\mathrm{\Delta }{t}^{-\mathrm{1}}\mathrm{\nabla }\cdot {\stackrel{\mathrm{‾}}{\mathit{u}}}^{\ast }.$$

The horizontal part of the operator to invert in the elliptic problem is treated in the Fourier space , and its vertical part leads to the classical tridiagonal matrix. The mathematical operator to invert $\mathrm{\nabla }\cdot \left(\mathbf{\nabla }\right)$ is exact for flat problems . 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 or a residual conjugate gradient 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. () are estimated. The LES closure is done by an eddy–diffusivity approach called 1.5TKE with a 1.5-order closure scheme . 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}}})$:

6$$\begin{array}{rl}{\displaystyle \frac{\partial \left({\mathit{\rho }}_{\mathrm{r}}e\right)}{\partial t}}{|}_{\mathrm{csl}}=& -\mathrm{\nabla }\left({\mathit{\rho }}_{\mathrm{r}}e\stackrel{\mathrm{‾}}{\mathit{u}}\right)-{\mathit{\rho }}_{\mathrm{r}}\mathbf{g}{\displaystyle \frac{\stackrel{\mathrm{‾}}{{\mathit{\theta }}^{\prime }{\mathbf{u}}^{\prime }}}{\stackrel{\mathrm{‾}}{\mathit{\theta }}}}-{\mathit{\rho }}_{\mathrm{r}}\stackrel{\mathrm{‾}}{{\mathbf{u}}^{\prime }\otimes {\mathbf{u}}^{\prime }}\cdot \mathrm{\nabla }\stackrel{\mathrm{‾}}{\mathit{u}}\\ & +\mathrm{\nabla }\cdot \left({\mathit{\rho }}_{\mathrm{r}}{K}_e{l}_m\sqrt{e}\mathrm{\nabla }e\right)-{\mathit{\rho }}_{\mathrm{r}}{K}_{\mathit{ϵ}}{e}^{\mathrm{3}/\mathrm{2}}/l_{\mathit{ϵ}},\end{array}$$ where $K_e$ and $K_{\mathit{ϵ}}$ are constants prescribed in the turbulence scheme, and $l_m$ and $l_{\mathit{ϵ}}$ are the length scales defining the turbulent viscosity. The dissipation term is directly estimated from $e$ and $l_{\mathit{ϵ}}$ (the left-hand term in Eq. ). 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$. $$\begin{array}{}\text{7}& {\displaystyle }\begin{array}{rl}\stackrel{\mathrm{‾}}{u_i^{\prime }{u}_j^{\prime }}=& +(\mathrm{2}/\mathrm{3}){\mathit{\delta }}_{ij}e-(\mathrm{1}/\mathrm{15})l_m\sqrt{e}\left({\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{u_i}}{\partial x_j}}+{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{u_j}}{\partial x_i}}\right.\\ & \left.-(\mathrm{2}/\mathrm{3}){\mathit{\delta }}_{ij}{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{u_m}}{\partial x_m}}\right)\end{array}\text{8}& {\displaystyle }\stackrel{\mathrm{‾}}{{\mathit{\theta }}^{\prime }{u}_i^{\prime }}=-(\mathrm{1}/\mathrm{6})l_m\sqrt{e}{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{\mathit{\theta }}}{\partial x_i}}{\mathrm{\Phi }}_i\text{9}& {\displaystyle }\stackrel{\mathrm{‾}}{{{\mathit{\theta }}^{\prime }}^{\mathrm{2}}}=+(\mathrm{5}/\mathrm{36})l_m^{\mathrm{2}}{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{\mathit{\theta }}}{\partial x_m}}{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{\mathit{\theta }}}{\partial x_m}}{\mathrm{\Phi }}_m\end{array}$$ The Einstein summation convention is applied and ${\mathrm{\Phi }}_{i,m}$ represents atmospheric stability functions . The ground condition can be modeled by the externalized surface scheme SURFEX . 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^{\ast }$ proportional to the turbulent wall shear and a roughness length $z_{\mathrm{0}}$, the vertical gradient of $\stackrel{\mathrm{‾}}{\mathit{u}}$ is recovered by specifying a logarithmic profile 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 term is neglected and functions do not appear in the previous formulation). SURFEX is employed in Sect. .

3 The IBM forcing in the Meso-NH code

The numerical domain is divided into two regions: where the equations of continuum mechanics hold and a solid region embedding the obstacle where they do not. After comparing the methods (Fig. a) based on a local volume fraction function and the LSF , it was decided to use the LSF as it was able to capture the interface between the fluid and solid regions more accurately. The $\mid \mathit{\varphi }\mid$ distance informs us about the minimal distance to the fluid–solid interface and the $\mathit{\varphi }$ sign about the region nature: sgn$\left(\mathit{\varphi }\right)>\mathrm{0}$ for the solid one; otherwise, sgn$\left(\mathit{\varphi }\right)<\mathrm{0}$. The vector $\mathit{n}$ normal to the interface and its local curvature $\mathit{\sigma }$ are defined as $\mathit{n}=\frac{\mathrm{\nabla }\mathit{\varphi }}{\left|\mathrm{\nabla }\mathit{\varphi }\right|}$ and $\mathit{\sigma }=-\mathrm{\nabla }\cdot \mathit{n}$. Figure a illustrates the continuous variation of LSF for an arbitrary bell-shaped interface. The LSF is estimated at the seven available point types per cell to limit the discretization errors (Fig. b): at the mass point $P$, where prognostic scalar variables are localized, at the three velocity nodes $U,V,W$ where each projection $\stackrel{\mathrm{‾}}{\mathit{u}}$ is characterized, and at the A, B, C vorticity nodes employed by turbulent variables. The points of the solid region act as external points of the computational grid (as do external points in a boundary-fitted method at the grid limit). An intensive study has been done to estimate the modeling of the vector normal to the interface and the local curvature using LSF. The forcing based on a ghost-cell technique (GCT; and CCT or cut-cell technique) is applied to the explicit-in-time schemes (and the pressure solver) and detailed in Sect.  (and Sect. ).

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

The ${\mathit{\psi }}^n$ value is estimated at the time $n\mathrm{\Delta }t$ ($\mathrm{\Delta }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. ) cannot be deduced from the conservation laws in the solid region. Expecting a correction due to IBM, wherein $\mathit{\varphi }\ge \mathrm{0}$, a general formulation of the tendencies is written as

10

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

[Figure omitted. See PDF]

In classical GCT 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 a shows the characterization of one ghost $G$ (of LSF value ${\mathit{\varphi }}_G$), its associated mirror $I$ (of LSF value ${\mathit{\varphi }}_I$) and the interface node $B$ ($GI=\mathrm{2}{\mathit{\varphi }}_G\mathit{n}$). Figure b 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 $\mathrm{\Delta }$ the space step), leading to a not well-posed condition.

The new GCT. To overcome this problem, we define image points (noted $I_{\mathrm{1}}$ and $I_{\mathrm{2}}$ in Fig. a; 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 a shows the images for one ghost. The new approach enforces a large enough value of the $\left|I_{l}B\right|$ distance. Figure b (c) illustrates classical (new) GCT for several ghosts. Figure b shows some mirror points associated with ghosts of the first layer in the vicinity of the interface. Figure c 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${}^{\mathrm{a}}$ and PLI${}^{\mathrm{b}}$ are tested to build the Lagrangian interpolation: $$\begin{array}{}\text{12}& {\displaystyle }\begin{array}{rl}{\stackrel{\mathrm{‾}}{\mathit{\psi }}}^{\mathrm{a}}\left(I\right)=& \left[\mathrm{2}{L}_G^{\mathrm{a}}\left(I\right)\stackrel{\mathrm{‾}}{\mathit{\psi }}\left(B\right)+L_{I_{\mathrm{1}}}^{\mathrm{a}}\left(I\right)\stackrel{\mathrm{‾}}{\mathit{\psi }}\left(I_{\mathrm{1}}\right)\right.\\ & \left.+L_{I_{\mathrm{2}}}^{\mathrm{a}}\left(I\right)\stackrel{\mathrm{‾}}{\mathit{\psi }}\left(I_{\mathrm{2}}\right)\right]{\left[\mathrm{1}+L_G^{\mathrm{a}}\left(I\right)\right]}^{-\mathrm{1}},\end{array}\text{13}& {\displaystyle }{\stackrel{\mathrm{‾}}{\mathit{\psi }}}^{\mathrm{b}}\left(I\right)=L_B^{\mathrm{b}}\left(I\right)\stackrel{\mathrm{‾}}{\mathit{\psi }}\left(B\right)+L_{I_{\mathrm{1}}}^{\mathrm{b}}\left(I\right)\stackrel{\mathrm{‾}}{\mathit{\psi }}\left(I_{\mathrm{1}}\right)+L_{I_{\mathrm{2}}}^{\mathrm{b}}\left(I\right)\stackrel{\mathrm{‾}}{\mathit{\psi }}\left(I_{\mathrm{2}}\right),\end{array}$$ where $L^{\mathrm{a}}\left(I\right)$ and $L^{\mathrm{b}}\left(I\right)$ are the Lagrangian polynomials, as follows.

$$\begin{array}{}\text{14}& {\displaystyle }\begin{array}{rl}& L_{I_{\mathrm{1}}}^{\mathrm{a}}\left(I\right)=\left({\displaystyle \frac{\mathrm{2}\mathrm{\Delta }-\mathit{\varphi }}{\mathrm{\Delta }}}\right)\left({\displaystyle \frac{\mathrm{2}\mathit{\varphi }}{\mathrm{\Delta }+\mathit{\varphi }}}\right)\\ & L_{I_{\mathrm{2}}}^{\mathrm{a}}\left(I\right)=\left({\displaystyle \frac{\mathit{\varphi }-\mathrm{\Delta }}{\mathrm{\Delta }}}\right)\left({\displaystyle \frac{\mathrm{2}\mathit{\varphi }}{\mathrm{2}\mathrm{\Delta }+\mathit{\varphi }}}\right)\\ & L_G^{\mathrm{a}}\left(I\right)=\left({\displaystyle \frac{\mathit{\varphi }-\mathrm{\Delta }}{\mathit{\varphi }+\mathrm{\Delta }}}\right)\left({\displaystyle \frac{\mathit{\varphi }-\mathrm{2}\mathrm{\Delta }}{\mathit{\varphi }+\mathrm{2}\mathrm{\Delta }}}\right)\end{array}\text{15}& {\displaystyle }\begin{array}{rl}& L_{I_{\mathrm{1}}}^{\mathrm{b}}\left(I\right)=\left({\displaystyle \frac{\mathrm{2}\mathrm{\Delta }-\mathit{\varphi }}{\mathrm{\Delta }}}\right)\left({\displaystyle \frac{\mathit{\varphi }}{\mathrm{\Delta }}}\right)\\ & L_{I_{\mathrm{2}}}^{\mathrm{b}}\left(I\right)=\left({\displaystyle \frac{\mathit{\varphi }-\mathrm{\Delta }}{\mathrm{\Delta }}}\right)\left({\displaystyle \frac{\mathit{\varphi }}{\mathrm{2}\mathrm{\Delta }}}\right)\\ & L_B^{\mathrm{b}}\left(I\right)=\left({\displaystyle \frac{\mathit{\varphi }}{\mathrm{\Delta }}}\right)\left({\displaystyle \frac{\mathit{\varphi }-\mathrm{2}\mathrm{\Delta }}{\mathrm{2}\mathrm{\Delta }}}\right)\end{array}\end{array}$$

Figure 3

Quadratic 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^{\mathrm{a}}$ ($L^{\mathrm{b}}$) polynomial results.

[Figure omitted. See PDF]

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. a) and $\mathit{\psi }={\mathit{\varphi }}^{\mathrm{1}/\mathrm{4}}$ (Fig. b) is studied. As illustrated, PLI${}^{\mathrm{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. $$\begin{array}{}\text{16}& {\displaystyle }\begin{array}{rl}& L_i^{\mathrm{LP}}\left(x_l\right)=\stackrel{N}{\prod _{p=\mathrm{1},p\ne i}}{\displaystyle \frac{x_l-x_p}{x_i-x_p}}\\ & L_j^{\mathrm{LP}}\left(y_l\right)=\stackrel{N}{\prod _{p=\mathrm{1},p\ne j}}{\displaystyle \frac{y_l-y_p}{y_j-y_p}}\\ & L_k^{\mathrm{LP}}\left(z_l\right)=\stackrel{N}{\prod _{p=\mathrm{1},p\ne k}}{\displaystyle \frac{z_l-z_p}{z_k-z_p}}\end{array}\text{17}& {\displaystyle }\begin{array}{rl}\stackrel{\mathrm{‾}}{\mathit{\psi }}(x_l,y_l,z_l)=& \stackrel{N}{\sum _{i=\mathrm{1}}}\stackrel{N}{\sum _{j=\mathrm{1}}}\stackrel{N}{\sum _{k=\mathrm{1}}}{L}_i^{\mathrm{LP}}\left(x_l\right)L_j^{\mathrm{LP}}\left(y_l\right)L_k^{\mathrm{LP}}\left(z_l\right)\\ & \cdot \stackrel{\mathrm{‾}}{\mathit{\psi }}(x_i,y_j,z_k)\end{array}\end{array}$$

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: 18$$\begin{array}{rl}& \stackrel{\mathrm{‾}}{\mathit{\psi }}\left(x_l\right)={\displaystyle \frac{\stackrel{N}{\sum _{i=\mathrm{1}}}{L}_i^{\mathrm{DW}}\left(x_l\right).\stackrel{\mathrm{‾}}{\mathit{\psi }}\left(x_i\right)}{\stackrel{N}{\sum _{i=\mathrm{1}}}{L}_i^{\mathrm{DW}}\left(x_l\right)}};\\ & \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}}},\end{array}$$ where $L_i^{\mathrm{DW}}\left(x\right)=\mid x_l-x_i{\mid }^{-\mathit{\alpha }}$ ($\mathit{\alpha }=\mathrm{1}$). This formulation diverges when $x_i\to 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{ϵ}$ ($\mathit{ϵ}$ is an arbitrary parameter depending on the mesh discretization, $\mathit{ϵ}\ll \mathrm{\Delta }$). 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 and tests on the $\mathit{\alpha }$ 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. b) 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$ ($\mathrm{\Delta }{n}_{\mathrm{1}}=\mathrm{\Delta }$ and $\mathrm{\Delta }{n}_{\mathrm{2}}=\mathrm{2}\mathrm{\Delta }$ in Fig. a), 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, ($\mathbf{t},\mathbf{c}$) represents two arbitrary tangent directions. The tangent direction $\mathbf{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 19$$\mathbf{c}\left(I_l\right)={\displaystyle \frac{\mathit{n}\otimes \stackrel{\mathrm{‾}}{\mathit{u}}\left(I_l\right)}{\mid \mid \mathit{n}\otimes \stackrel{\mathrm{‾}}{\mathit{u}}\left(I_l\right)\mid \mid }};\mathbf{t}\left(I_l\right)=\mathbf{c}\left(I_l\right)\otimes \mathit{n}.$$

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:

20$$\begin{array}{rl}& \mathbf{t}\left(B\right)=\mathbf{t}\left(I_{\mathrm{1}}\right)\text{ if no rotation};\\ & {\mathbf{e}}_{\mathbf{t}}\left(B\right)=\mathrm{2}{\mathbf{e}}_{\mathbf{t}}\left(I_{\mathrm{1}}\right)-{\mathbf{e}}_{\mathbf{t}}\left(I_{\mathrm{2}}\right)\text{ if linear evolution}.\end{array}$$

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. .

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 })$, $\mathit{\psi }\left(G\right)$ for a Dirichlet (Neumann) condition is written as 21$$\begin{array}{rl}& \mathit{\psi }\left(G\right)=\mathrm{2}\mathit{\psi }\left(B\right)-\mathit{\psi }\left(I\right)\text{ (Dirichlet)},\\ & \mathit{\psi }\left(G\right)=\mathrm{2}\mathit{\varphi }{\displaystyle \frac{d\mathit{\psi }}{dn}}{|}_B+\mathit{\psi }\left(I\right)\text{ (Neumann)}.\end{array}$$

Figure 4

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

[Figure omitted. See PDF]

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 a and b illustrate these interpolations by considering two analytical profiles (red lines): the quadratic interpolation $Q{I}_{\mathrm{1}}$ ($Q{I}_{\mathrm{2}}$) is based on the image values located at $\mathit{\varphi }=\mathrm{1}/\mathrm{2}\mathrm{\Delta }$ and $\mathit{\varphi }=\mathrm{\Delta }$ and plotted in green symbols (at $\mathit{\varphi }=\mathrm{\Delta }$ and $\mathit{\varphi }=\mathrm{2}\mathrm{\Delta }$ plotted in blue symbols). Depending on the analytical profile, Fig.  shows the influence of the image location choice. As expected, $Q{I}_{\mathrm{1}}$ ($Q{I}_{\mathrm{2}}$) appears to be less accurate than $Q{I}_{\mathrm{2}}$ ($Q{I}_{\mathrm{1}}$) 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)$). $Q{I}_C$ is the combination of $Q{I}_{\mathrm{1}}$ and $Q{I}_{\mathrm{2}}$ (purple line). $Q{I}_C$ preserves the advantage of each quadratic interpolation and when ${\mathit{\varphi }}_G<\mathrm{\Delta }$ (${\mathit{\varphi }}_G>\mathrm{\Delta }$), $Q{I}_{\mathrm{1}}$ ($Q{I}_{\mathrm{2}}$) 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 $Q{I}_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. (), 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. ) 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 ${\mathcal{Q}}_{\mathrm{cct}}$):

22$$\begin{array}{rl}\mathcal{V}\mathrm{\nabla }\cdot {\stackrel{\mathrm{‾}}{\mathit{u}}}^{\ast }& =\underset{{\mathcal{V}}_{\mathrm{f}}}{\int }\mathrm{\nabla }\cdot {\stackrel{\mathrm{‾}}{\mathit{u}}}^{\ast }\mathrm{d}\mathcal{V}+\underset{{\mathcal{V}}_{\mathrm{s}}}{\int }\mathrm{\nabla }{\stackrel{\mathrm{‾}}{\mathit{u}}}^{\ast }\mathrm{d}\mathcal{V}=\sum \pm {\stackrel{\mathrm{‾}}{u_i}}^{\ast }\cdot {\mathcal{S}}_i\\ & =\sum \pm \stackrel{\mathrm{̃}}{{\mathrm{\Delta }}^{\mathrm{2}}{\stackrel{\mathrm{‾}}{u_i}}^{\ast }},\end{array}$$ where $\mathcal{V}={\mathrm{\Delta }}^{\mathrm{3}}$ is the cell volume, ${\mathcal{V}}_{\mathrm{f}}$ (${\mathcal{V}}_{\mathrm{s}}$) the fluid (solid) part of $\mathcal{V}$ and ${\mathcal{S}}_i$ the cell surfaces for which $i$ is the index of each surface orientation [e, w, n, s, b, f] as illustrated in Fig. a.

Figure 5

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

[Figure omitted. See PDF]

According to the Green–Ostrogradski theorem, the ${\stackrel{\mathrm{‾}}{u_i}}^{\ast }{\mathcal{S}}_i$ calculation is the classical way of a CCT 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. a). Figure b 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 $\mathit{\varphi }={\mathit{\varphi }}_{\mathrm{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 ${\mathit{\delta }}_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) 23$$\begin{array}{rl}\stackrel{\mathrm{̃}}{{\mathrm{\Delta }}^{\mathrm{2}}{\stackrel{\mathrm{‾}}{u_i}}^{\ast }}=& {\displaystyle \frac{{\mathrm{\Delta }}^{\mathrm{2}}}{\mathrm{8}}}\left[\stackrel{8}{\sum _{p=\mathrm{1}}}\mathcal{H}(-{\mathit{\varphi }}_p)\mathcal{H}(-{\mathit{\varphi }}_i){\stackrel{\mathrm{‾}}{u_i}}^{\ast }\right.\\ & +\stackrel{8}{\sum _{p=\mathrm{1}}}\mathcal{H}(-{\mathit{\varphi }}_p{\mathit{\varphi }}_i)\cdot \left|\mathcal{H}\right(-{\mathit{\varphi }}_p)-\mathit{\pi }{\left({\displaystyle \frac{{\mathit{\delta }}_p}{\mathrm{\Delta }}}\right)}^{\mathrm{2}}|\\ & \left.\cdot \left(\mathcal{H}(-{\mathit{\varphi }}_p){\stackrel{\mathrm{‾}}{u}}_p^{\ast }+\mathcal{H}(-{\mathit{\varphi }}_i){\stackrel{\mathrm{‾}}{u_i}}^{\ast }\right)\right].\end{array}$$

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 ${\mathit{\varphi }}_p<\mathrm{0}$ and ${\mathit{\varphi }}_i<\mathrm{0}$ (Fig. a); a pure solid cell $\stackrel{\mathrm{̃}}{{\mathrm{\Delta }}^{\mathrm{2}}{\stackrel{\mathrm{‾}}{u_i}}^{\ast }}=\mathrm{0}$ when ${\mathit{\varphi }}_p>\mathrm{0}$ and ${\mathit{\varphi }}_i>\mathrm{0}$ (Fig. b); 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. c–d). Using Eq. (), Eq. () is solved and leads to the RHS computation of Eq. ().

Figure 6

(a–d) $\pm \stackrel{\mathrm{̃}}{{\mathrm{\Delta }}^{\mathrm{2}}{\stackrel{\mathrm{‾}}{u_i}}^{\ast }}$ calculations depending on the signs of ${\mathit{\varphi }}_i=\left(\mathit{\varphi }\right)$ and ${\mathit{\varphi }}_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. ).

[Figure omitted. See PDF]

Knowing ${\mathcal{Q}}_{\mathrm{cct}}$, the reflection now concerns the $\mathcal{P}$ matrix to invert. The classical interface condition on the potential ${\mathrm{\Psi }}^{\ast }$ 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 appears only on the border of a numerical domain. Using an IBM and without any impact of this interface condition on the $\mathcal{P}$ coefficients, the impermeability character of solid obstacles is not achieved. Due to the inversion of the horizontal part of $\mathcal{P}$ by a fast Fourier transform , the solution of calculating ${\mathcal{P}}_{\mathrm{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 ${\mathrm{\Psi }}^{\ast }$ condition in $\mathcal{P}$ 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. () 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. .

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_{\mathit{ϵ}}$ (hereafter noted $l_m$ and called the mixing length). The 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. ) is modified in the presence of immersed obstacles in the subgrid turbulent kinetic energy equation (Eq. ), mixing length computation and Reynolds stress diagnosis (Eqs. , and ).

The subgrid turbulent kinetic energy condition. The explicit-in-time resolution of Eq. () claims a GCT forcing and an interface condition on the STKE $e$. Commonly, the STKE profile is considered parabolically in the viscous sublayer and constant in the inertial and wake outer layers . 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 $({\mathit{\nu }}^{\mathrm{3}}{\mathit{ϵ}}^{-\mathrm{1}}{)}^{\frac{\mathrm{1}}{\mathrm{4}}}$. 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\approx \mathcal{O}\left(\mathrm{10}m\right)$; $\mathit{\text{Re}}\approx \mathcal{O}\left({\mathrm{10}}^{\mathrm{7}}\right)$. This thickness estimate is also proportional to $E\mathit{\nu }/u^{\ast }$ ($E\approx \mathrm{9.8}$ is commonly employed), whereby the friction velocity $u^{\ast }$ 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. ). 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^{\mathrm{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 })$.

24

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

[Figure omitted. See PDF]

Figure  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_{\mathrm{\infty }}$ and parallel streamlines (Fig. a). 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_{\mathrm{cyl}}=\mathrm{2}{R}_{\mathrm{cyl}}$ is achieved. Figure b shows the streamlines obtained with the MNH pressure solver modified to take into account the presence of immersed obstacles (Sect. ). Defining $\mathbf{x}$ as the direction parallel to the initial streamlines and $\mathbf{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, $(\mathit{\alpha };r)$ cylindrical coordinates). The numerical confinement is discussed hereafter, characterized by $L=L_{\mathrm{cyl}}/R_{\mathrm{cyl}}$, where $L_{\mathrm{cyl}}$ is the distance separating the lateral domain surfaces (Fig. a).

Figure 9

Potential 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=\mathrm{16}$); (b) evolution of the added mass coefficient $A(N;L)$ with the confinement $L=L_{\mathrm{cyl}}/R_{\mathrm{cyl}}$ ($N=\mathrm{16}$) and with the node number per radius cylinder $N$ ($L=\mathrm{16}$). The confinement is defined in Fig. a.

[Figure omitted. See PDF]

The Richardson (RICH) and the residual conjugate gradient iterative (RESI) methods are tested (Sect. ). Figure a plots the evolution of the dimensionless residue $R\left(k\right)$ (based on a characteristic divergence $U_{\mathrm{\infty }}/\mathrm{\Delta }$) with the iteration number $k$ and obtained with the confinement $L=\mathrm{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_{\mathrm{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\left(k\right)$. 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_{\mathrm{cyl}}$ location (Fig. a) 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_{\mathrm{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 ${\mathbf{a}}_{\mathrm{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 $\mathrm{\Delta }t\to \mathrm{0}$. $\mathcal{A}$ is the dimensionless coefficient and $m_{\mathrm{f}}$ the displaced fluid mass. ${\mathcal{A}}_{\mathrm{cyl}}$ theoretically equals $\mathrm{1}$ in the non-confined cylinder case . The red curve in Fig. b illustrates the effect of the confinement $L$ on ${\mathcal{A}}_{\mathrm{cyl}}$ for $N=\mathrm{16}$ resolution. Unsurprisingly, ${\mathcal{A}}_{\mathrm{cyl}}$ increases with the confinement . The weak dependence of ${\mathcal{A}}_{\mathrm{cyl}}$ with $L>\mathrm{16}$ allows us to consider the body to be isolated for $L\sim \mathrm{16}$. The green curve in Fig. b shows the impact of the space resolution for $L=\mathrm{16}$. The numerical added mass coefficient is in good agreement with the theoretical one, presenting a relative error of about $\mathrm{2}\mathit{\%}$ for $N>\mathrm{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  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 .

Figure 10

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

[Figure omitted. See PDF]

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 , who experimentally studied the nature of eddy structures.

Figure 12

Recirculating region at $\mathit{\text{Re}}=\mathrm{40}$ (MNH–IBM, $\mathrm{20}\mathrm{pts}/D_{\mathrm{cyl}}$): definition of the ${\mathit{\theta }}_{\mathrm{d}}$ (${}^{\circ }$) separation angle, $l_{\mathrm{r}}$ recirculating length and $(a;b)$ vortex core location. The distance is dimensionless by $D_{\mathrm{cyl}}$.

[Figure omitted. See PDF]

The $\mathrm{20}\mathrm{pts}/D_{\mathrm{cyl}}$ and $\mathrm{30}\mathrm{pts}/D_{\mathrm{cyl}}$ meshes present a good spatial convergence and weak differences at $\mathit{\text{Re}}=\mathrm{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 $\mathrm{10}\mathrm{pts}/D_{\mathrm{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. ). Table  compares the $\mathrm{20}\mathrm{pts}/D_{\mathrm{cyl}}$ results with a part of the results literature collected in .

Table 1

Description of the standing eddies in the wake of the solid cylinder ($\mathit{\text{Re}}=\mathrm{40}$): comparison of the separation angle ${\mathit{\theta }}_{\mathrm{d}}$ (${}^{\circ }$), recirculating length $l_{\mathrm{r}}$ (m) and vortex core location $(a;b)$ between the literature and MNH–IBM.

Authors	${\mathit{\theta }}_{\mathrm{d}}$ (${}^{\circ }$)	$l_{\mathrm{r}}/D_{\mathrm{cyl}}$	$b/D_{\mathrm{cyl}}$	$a/D_{\mathrm{cyl}}$
	53.8	2.13	0.76	0.59
	53.6	2.28	0.72	0.60
	53.7	2.30	0.73	0.60
	53.4	2.26	0.71	0.60
	53.6	2.24	0.71	0.59
	54.5	2.34	0.76	0.62
MNH–IBM $\mathrm{30}\mathrm{pts}/D_{\mathrm{cyl}}$	$\approx \mathrm{54}$	$\approx \mathrm{2.2}$	$\approx \mathrm{0.7}$	$\approx \mathrm{0.6}$

The focus is on the unsteady mode at $\mathit{\text{Re}}=\mathrm{140}$. The ratio between the characteristic time of inertial effects $D_{\mathrm{cyl}}/U_{\mathrm{\infty }}$ and the one related to the vortex shedding $\mathrm{1}/f$ defines the Strouhal number $\mathit{\text{St}}=\frac{fD}{U_{\mathrm{\infty }}}$. , and 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 . 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 $\mathrm{2}\mathit{\%}$ in regard to the 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 Turbulent flows around parallelepiped(s)

This section is devoted to turbulent flows approaching our perspective: the simulation of an atmospheric flow over a city. The turbulent flows around a cubic body vertically confined in a channel and over an urban-like roughness (set of obstacles) are here described. MNH–IBM is explicitly compared to experimental investigations in the two cases. Comparisons to other LESs from the literature will be mentioned.

Common hypothesis and methods. The fluid is considered as neutrally stratified. The Coriolis term is negligible due to the addressed space scales and timescales. The turbulent diffusion is modeled by the subgrid TKE1.5 scheme transported by PPM (Sect. ). All surfaces are considered non-permeable and the IBM wall model (Sect. ) is activated. An $(x,y,z)$ reference frame is defined ($z$, vertical direction) and the velocity vector is written as $\mathbf{u}\left(\mathbf{t}\right)=u\left(t\right)\mathbf{x}+v\left(t\right)\mathbf{y}+w\left(t\right)\mathbf{z}$. A time simulation is needed afterwards to establish the turbulence state (not shown here). The overline notation refers to the mean value in time in this section.

5.1 Flow over a surface-mounted cube

Using static pressure measurements, as well as laser-sheet and oil-film visualizations, and generated a large data set for the study of flows around a cubic body placed in a channel (Fig. a). RANS and LES have been used to explore in detail this physical case .

Figure 14

Vertical symmetry plane of the mean flow: (a–c) MNH–IBM time-averaged streamlines; (d) streamlines observed by ; (e) MNH–IBM instantaneous visualization of the Q criterion .

[Figure omitted. See PDF]

Results. Figure a, 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 , who used a far more precise grid near walls. The same figure obtained by the experimental investigation is given in Fig. d. The size of the front (rear) region is characterized by the recirculating length $x_{\mathrm{f}}/H$ ($x_{\mathrm{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 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 . This bifurcation point was also commented on in , even if their experimental uncertainty did not allow us to visualize it in Fig. d.

Figure e shows an MNH–IBM instantaneous flow field with the Q criterion , 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\left(t\right)$ 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 15

Mean vertical profiles of velocities (top), turbulent kinetic energy and Reynolds stress (bottom). The lines correspond to the MNH–IBM results. The symbols are the data except for TKE . 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}$.

[Figure omitted. See PDF]

Still in the vertical symmetry plane of the mean flow, Fig.  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 ($\mathrm{10}\mathrm{pts}/H$ in red; $\mathrm{20}\mathrm{pts}/H$ in green; $\mathrm{40}\mathrm{pts}/H$ in blue). Note that the $U_{\mathrm{b}}$ bulk velocity was set to unity and therefore presents variables in Fig.  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.  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 $S\mathrm{1}$ (Fig. b). An underestimation appears on the counterflow at $(x/H;z/H)\approx (\mathrm{2};\mathrm{0})$ and is frequently observed in the literature (Fig. c). Some discrepancies are found on two $\stackrel{\mathrm{‾}}{w}$ profiles (Fig. f–g). Note that and highlight the difficulty of recovering it. The turbulent kinetic energy and the Reynolds stress are correctly predicted at the cube roof (Fig. i–j). The vertical profiles of TKE and $\stackrel{\mathrm{‾}}{u^{\prime }{w}^{\prime }}$ downstream of the cube show an overestimate tendency (Fig. k, p, l). No experimental data are available on TKE at $x/H=\mathrm{4}$ (Fig. l), 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 .

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 $S\mathrm{1}$ 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 , 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 $S\mathrm{2}$ 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_{\mathrm{0}}$ roughness length fixed in the immersed wall model (Sect. ) 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_{\mathrm{tke}}$ is applied on the channel surfaces (non-slip condition), $x_{\mathrm{f}}$ highly increases and the vortex at $S\mathrm{2}$ appears. Otherwise, $K_{\mathrm{tke}}$ and $z_{\mathrm{0}}$ 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 , , 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_{\mathrm{tke}}\approx \mathrm{2}$ in Eq. ().

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 . 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 17

Mean 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 experimental measurements (MNH–IBM results).

[Figure omitted. See PDF]

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 $-\mathrm{41}$${}^{\circ }$ 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. a, 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${}^{-\mathrm{1}}$ and allows us to build the vertical profile of the mean incoming flow (Fig. a, blue line). This $u^{\ast }$ value is close to the experimental one $u_{\exp }^{\ast }=\mathrm{0.68}$ m s${}^{-\mathrm{1}}$ found by a sonic anemometer at the feet of the S tower. The same inlet condition is used in the numerical studies of , , and . 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 in Figure a (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\mathit{\lesssim }\mathrm{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 . A wind acceleration and an increase in the wind angle are observed by MNH–IBM for $z\mathit{\gtrsim }\mathrm{8}$ m as in the LES results of and 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 18

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

[Figure omitted. See PDF]

Figure a (b) illustrates the time-averaged horizontal wind field $\left|\stackrel{\mathrm{‾}}{\mathit{u}}\right|$ ($\left|\mathit{u}\right|$ at $t=\mathrm{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 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_{\mathrm{0}}$ 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 2

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

	$E_{\mathrm{k}}$ (m${}^{\mathrm{2}}$ s${}^{\mathrm{2}})$		TKE (m${}^{\mathrm{2}}$ s${}^{\mathrm{2}})$		$u^{\ast }$ (m s${}^{-\mathrm{1}})$
	MNH–IBM		MNH–IBM		MNH–IBM
$z=\mathrm{4}$ m	15.3	14.1	3.33	3.78	1.14	1.08
$z=\mathrm{8}$ m	36.7	29.7	1.70	3.28	0.68	0.83
$z=\mathrm{16}$ m	65.8	48.5	0.02	1.75	0.02	0.60

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 . All the variables are in good agreement with the experiment at $z=\mathrm{4}$ m (about 2 times the container height). The friction velocity at $z\left(\mathrm{T}\right)=\mathrm{4}$ m is at least 2 times larger than $u_{\exp }^{\ast }\left(S\right)=\mathrm{0.68}$ m s${}^{-\mathrm{1}}$ 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\left(\mathrm{T}\right)=\mathrm{8}$ m and $z\left(\mathrm{T}\right)=\mathrm{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=\mathrm{16}$ m is close to the upstream value.

Figure 19

Spectrum of the measured (blue line) and modeled (green line) $u$ wind component at the T tower at $z=(\mathrm{4};\mathrm{8};\mathrm{16})$ m. $f^{\ast }$ is dimensionless as a Strouhal number using $|\mathit{u}{|}_{\mathrm{inlet}}$($z=\mathrm{2.5}$ m), and $z=\mathrm{2.5}$ m is the container length.

[Figure omitted. See PDF]

The discrete Fourier transform of the $u$ temporal evolution (Fig. ) at $z\left(\mathrm{T}\right)=(\mathrm{4};\mathrm{8})$ m shows the coherence between the energetic cascade of the experimental investigations and that of the numerical ones. At $z\left(\mathrm{T}\right)=\mathrm{16}$ m, MNH–IBM underestimates the unsteady part of the solution for all wave numbers. The same behavior is observed on $v$ and $w$ (not shown here).

The MNH–IBM results are consistent with the experimental observations below $z\left(\mathrm{T}\right)\mathit{\lesssim }\mathrm{10}$ m. This modeling induces the turbulence to be mostly due to the container wakes upstream of the T tower and not to the incoming turbulence (not modeled in the simulation). The influence of the atmospheric turbulence grows with altitude and MNH–IBM diverges with the experiment. This divergence for $z\left(\mathrm{T}\right)\mathit{\gtrsim }\mathrm{10}$ m leads us to think that the thickness of the container has an influence about 4–5 times the height of the urban-like canopy.

6 Conclusions and perspectives

This study details the first implementation of an immersed boundary method (IBM) in the atmospheric Meso-NH (MNH) model currently based on mathematical formulations written for structured grids. The MNH–IBM aim is to explicitly model the fluid–solid interaction in the surface boundary layer developed over grounds presenting complex topographies such as cities or industrial sites.

A level-set function characterizes the geometric properties of the fluid–solid interface. Two original approaches of the ghost-cell and cut-cell techniques are implemented to correct the MNH numerical schemes. A newly proposed GCT recovers the fluid information in several image points by presenting a distance to the interface independent of the ghost's distance. The CCT consists of a new finite-volume approach of flux balance near the immersed interface. The GCT is applied to the numerical schemes based on explicit time integration and the CCT is employed in the implicit resolution of the Poisson equation, satisfying the incompressibility hypothesis. The adaptation and use of iterative procedures solves the pressure problem without any modification to the inverted matrix. The turbulence problem is closed at the fluid–solid interface by a pragmatic LES–RANS formulation based on the subgrid turbulent kinetic energy and the length of the smallest energetic eddies.

The pressure solver, adapted to the IBM and isolated from the rest of MNH, is used to model potential flows around several obstacles. Compared to analytical and theoretical solutions, the numerical results demonstrate the ability of the IBM adaptation to ensure that the momentum is preserved and the continuity equation is respected. Non-dissipative flows are simulated to test the IBM forcing of the wind advection scheme (the impact of interpolations collected in , classical and novel GCT comparison, and numerical diffusion near the interface). These tests validate the proposed GCT “three images and ghost points” using an inverse distance-weighting (trilinear) interpolation near (far from) the interface. The modeling of the advection term near the fluid–solid interface is ensured by a second-order centered scheme associated with an artificial viscosity. The artificial viscosity is calibrated after comparisons with a third-order WENO scheme. With these numerical choices, MNH–IBM demonstrates its ability to model wake instabilities past a circular cylinder placed in a viscous fluid. Then, large-eddy simulations of turbulent flows around bodies with sharp edges and corners are executed (i.e., a cube placed in a channel as in , and a near-neutral atmospheric application over an array of containers as in ). These two LESs validate the proposed immersed wall model, switching the characteristic space scale defining the turbulent Reynolds number between one obtained by either a viscous length scale (the cube case) or a roughness length (the container case).

Future work. This study constitutes a first robust step towards a better understanding of the interactions between “weather and cities” and better predictions of such interactions. The idealized character of the physical cases approached here offers some insights. One improvement would consist of a generalization of the IBM writing to terrain-following coordinates , allowing for the simulation of high-curvature bodies in the presence of non-flat ground. In the run-up of the resolution of multiscale problems, consistency with grid nesting would be pertinent as would coupling to a drag model . In the current paper, “simple” bodies are investigated; the modeling of real houses and buildings with arbitrary shapes in close proximity to each other is ongoing with work dedicated to a brief and intense pollutant episode due to a factory explosion in 2001 over Toulouse, assuming a dry neutral case with nonreactive gas dispersion. Such hypotheses involve a broad range of physical phenomena requiring numerical compliance with IBM, including chemical reactions, phase changes and radiative effects. Such compliance would allow for access to a large variety of atmospheric situations.

Code availability

The immersed boundary method has been implemented in the 5.2 version of the Meso-NH code. This reference version is under the CeCILL-C license agreement and freely available at http://mesonh.aero.obs-mip.fr/mesonh52 (last access: 13 May 2019).

The source files dedicated to IBM and the input files for the simulations in Sects.  and  can be downloaded from the CERFACS web page: https://cerfacs.fr/MNHIBM/Auguste-GMD-2019 . A Supplement to this article contains the tar ball (1Mo) of these files. The immersed boundary method will be integrated in a future Meso-NH version.

Appendix A Space convergence of the pressure solver

Array of vortices. A Poisson equation solution is investigated (Fig. a) by imposing in the RHS of Eq. () the divergence $\mathrm{\nabla }\cdot \stackrel{\mathrm{‾}}{{\mathit{u}}^{\ast }}=-\mathit{\pi }(l^{\mathrm{2}}+m^{\mathrm{2}})\cos \left(\mathit{\pi }lx\right)\sin \left(\mathit{\pi }my\right)$, where $l=m=\mathrm{cste}$. The error norms ($L_p=\sqrt[p]{\sum |P_{\mathrm{n}}-P_{\mathrm{t}}{|}^p}$, where $P_{\mathrm{n}}$ is the numerical pressure and $P_{\mathrm{t}}$ the theoretical one) are estimated in the presence (or not) of an immersed cylindrical body (Fig. ). The space second order of the pressure solver is recovered without IBM. The order decreases with IBM and stays consistent regarding the $L_{p=(\mathrm{\infty };\mathrm{1};\mathrm{2})}$ slopes. Note that an immersed square or sphere gives similar results.

Agnesi hill. The irrotational solution around two 2-D bell-shaped interfaces is investigated with IBM and the boundary-fitted method (BFM; terrain-following coordinates). The topography is characterized by a height $h_a$ and a shape $h\left(x\right)=\frac{h_a}{\mathrm{1}+{\left(\frac{k_a\cdot x}{h_a}\right)}^{\mathrm{2}}}$ ($k_a=\mathrm{4}$, bell 1; $k_a=\mathrm{8}$, bell 2). The bell slope is arbitrarily and respectively described here as gentle or steep. Figure  shows the pressure contours obtained with IBM (left) and BFM (right) for a gentle (top) and a steep (bottom) shape. The minimal pressure value is localized at the top of each bell and goes to zero far from this location. The reference BFM and IBM simulations with the fine resolution ($N=\mathrm{160}$ nodes per $h_a$, red) show a good agreement for each hill. The blue (green) corresponds to a coarser mesh employing $N/\mathrm{3}$ ($N/\mathrm{9}$) nodes per $h_a$. Weak differences appear between the $N$ and $N/\mathrm{3}$ meshes for both IBM and BFM, revealing a good space convergence (Fig. a–b). Numerical errors are visible with IBM near the interface, but the Venturi effect is well modeled. Differences become more significant with the $N/\mathrm{9}$ mesh, especially with the BFM BELL2 presenting the highest curvature value (Fig. d). IBM appears less accurate than BFM when the ground presents low curvature in regard to the space resolution. Otherwise, IBM seems more pertinent than BFM to model high interfaces such as sharp edges or corners. The minimum pressure that could reach $-\mathrm{\infty }$ is smoothed by IBM and allows the pressure solver not to diverge. Note that used a compressible WRF model and observed similar behaviors.

Figure A2

Potential flow function of the space resolution (color code) plotting the pressure contours around two bells (a, b: BELL1, gentle slope; c, d: BELL2, steep slope) and obtained with IBM (a, c) and the boundary-fitted method (b, d).

[Figure omitted. See PDF]

Appendix B Inviscid flow around a circular cylinder

For most atmospheric applications, the region size for which the fluid molecular viscosity ${\mathit{\nu }}_{\mathrm{f}}$ influences the dynamic is sufficiently small to be considered negligible (Sect. ). Solving the Euler equations, the impact of the numerical diffusion could be significant, especially near the fluid–solid interface. The adopted strategy with IBM is to model the advection term with a low-order scheme near the interface (Sect. ). The order decrease in IBM is not essential but allows us to limit the number of ghosts in the solid region, thereby limiting the communications during a parallel computation. Indeed, the chosen implementation implies that the associated images and ghost points have to be localized in the integration volume of each processor. The WEN3 third-order weighted essentially non-oscillatory and CEN2 second-order centered schemes are available in MNH. Far from the interface a CEN4 fourth-order centered scheme is employed. Table  summarizes the advection scheme nomenclature.

The vorticity equation for a 2-D inviscid flow reveals no production in time. The numerical vorticity production at the immersed surface of a cylindrical body is studied here by initializing the simulation with the potential solution. To fit the potential solution, a nontrivial condition is employed on the tangent velocity $\frac{{\partial }^{\mathrm{3}}{\stackrel{\mathrm{‾}}{u}}_{\mathrm{t}}}{\partial n^{\mathrm{3}}}=\mathrm{0}$. Expecting a numerical vorticity sufficiently controlled to avoid the flow separation, the effect of the artificial diffusion ${\mathit{\nu }}_{\mathrm{art}}\mathrm{\Delta }\stackrel{\mathrm{‾}}{\mathit{u}}$ injected with CEN2 is compared to the WEN3 intrinsically diffusive behavior. Furthermore, this study had estimated the 3-D interpolation impact (not detailed here).

Figure B1

Solving the Euler equations: (a) vorticity production ${\stackrel{\mathrm{‾}}{E}}_{\mathrm{s}}^{\ast }\left(t^{\ast }\right)$ influenced by the Courant number (CFL $=$ 0.8, line; CFL $=$ 0.2, symbol) and by the ${\mathit{\nu }}_{\mathrm{art}}$ artificial viscosity (red, green, blue, purple and cyan, respectively; ${\mathit{\nu }}_{\mathrm{art}}={\mathit{\nu }}_{\mathrm{art}}^{\mathrm{ref}}[\mathrm{1};\mathrm{2};\mathrm{4};\mathrm{16};\mathrm{256}{]}^{-\mathrm{1}}$) using an advection CEN2 second-order centered scheme near the interface; (b) velocity magnitude field obtained with a WEN3 third-order weighted essentially non-oscillatory scheme (top) and $\mathrm{CEN}\mathrm{2}+{\mathit{\nu }}_{\mathrm{art}}^{\mathrm{ref}}{\mathrm{256}}^{-\mathrm{1}}$ (bottom). CEN4 is the advection scheme used far from the interface. The mesh is the coarser one (MESH1).

[Figure omitted. See PDF]

Figure B2

Solving the Euler equations: (a) influence of the space resolution (red: MESH1; green: MESH2; blue: MESH3) on the vorticity production ${\stackrel{\mathrm{‾}}{E}}_{\mathrm{s}}^{\ast }\left(t^{\ast }\right)$ when the near-interface advection is modeled by WEN3 (line) and $\mathrm{CEN}\mathrm{2}+{\mathit{\nu }}_{\mathrm{art}}^{\mathrm{ref}}$ (symbol); (b) vorticity magnitude field obtained by WEN3 (top) and $\mathrm{CEN}\mathrm{2}+{\mathit{\nu }}_{\mathrm{art}}^{\mathrm{ref}}$ (bottom).

[Figure omitted. See PDF]

Figure a plots the evolution in time of the enstrophy ${\stackrel{\mathrm{‾}}{E}}_{\mathrm{s}}^{\ast }\left(t^{\ast }\right)=\frac{D_{\mathrm{cyl}}}{U_{\mathrm{\infty }}{\mathcal{V}}_{\mathrm{f}}}{\int }_{{\mathcal{V}}_{\mathrm{f}}}|\mathrm{\nabla }\times \stackrel{\mathrm{‾}}{\mathit{u}}|\mathrm{d}\mathcal{V}$ depending on the $\mathrm{\Delta }t$ time step and ${\mathit{\nu }}_{\mathrm{art}}$ using CEN2 near the interface (with $U_{\mathrm{\infty }}$ the velocity of the incoming flow and ${\mathcal{V}}_{\mathrm{f}}$ the integration volume in the fluid region). The enstrophy increases in time and reaches a mean value when the produced vorticity near the interface is evacuated from the numerical domain and in the body wake. Except for the simulations with low artificial diffusion (symbols and curves in cyan and purple), the vorticity production is weakly dependent on the physical time and CFL Courant number (symbols and curves in red, green and blue). It induces ${\mathit{\nu }}_{\mathrm{art}}$ proportional to $U_{\mathrm{\infty }}\mathrm{\Delta }x\approx \mathcal{O}\left(\frac{\mathrm{\Delta }{x}^{\mathrm{2}}\mathrm{CFL}}{\mathrm{\Delta }t}\right)$. A reference value of the artificial viscosity is also defined as ${\mathit{\nu }}_{\mathrm{art}}^{\mathrm{ref}}=\frac{\mathrm{\Delta }{x}^{\mathrm{2}}}{\mathrm{\Delta }t}\mathrm{CFL}$.

Figure b illustrates the vorticity field in the vicinity of the interface between the intrinsically diffusive WEN3 and $\mathrm{CEN}\mathrm{2}+{\mathit{\nu }}_{\mathrm{art}}$ with ${\mathit{\nu }}_{\mathrm{art}}={\mathit{\nu }}_{\mathrm{art}}^{\mathrm{ref}}{\mathrm{256}}^{-\mathrm{1}}$. The streamlines are maintained without detachment near the interface with WEN3. Otherwise, the CEN2 solution with low artificial diffusion presents numerical instabilities and vortex shedding.

Table B1

Summary of the mean wind advection scheme used (WENO: weighted essentially non-oscillatory).

	Advection	Space	Intervention location	Artificial
	scheme	order		diffusion
CEN4	centered	4	far from the interface	no
CEN2	centered	2	near the interface	yes
WEN3	WENO	3	near the interface	no

Figure a plots the enstrophy evolution for three meshes (color code) for WEN3 (lines) and CEN2+${\mathit{\nu }}_{\mathrm{art}}$ with ${\mathit{\nu }}_{\mathrm{art}}={\mathit{\nu }}_{\mathrm{art}}^{\mathrm{ref}}$ (symbols). MESH1 (10 nodes per $D_{\mathrm{cyl}}$), MESH2 (20 nodes per $D_{\mathrm{cyl}}$) and MESH3 (40 nodes per $D_{\mathrm{cyl}}$) are respectively the coarse, intermediate and fine mesh. The border of the numerical domain is always distant from the cylinder of more than $\mathrm{10}{R}_{\mathrm{cyl}}$. The CEN2+${\mathit{\nu }}_{\mathrm{art}}^{\mathrm{ref}}$ vorticity production appears fairly close to the WEN3 one for the three space resolutions. Figure b corroborates the last comment, presenting the vorticity contours dimensionless by $\frac{D}{U_{\mathrm{\infty }}}$. A suitable ${\mathit{\nu }}_{\mathrm{art}}$ combined with the CEN2 choice is also in the range of the too-diffusive WEN3 results and the growth of numerical instabilities. CEN2+${\mathit{\nu }}_{\mathrm{art}}^{\mathrm{ref}}{\mathrm{4}}^{-\mathrm{1}}$ is retained as the advection scheme of the mean wind near an immersed interface.

The supplement related to this article is available online at: https://doi.org/10.5194/gmd-12-2607-2019-supplement.

Author contributions

All authors contributed to the development of the source code. The execution and the exploitation of the presented simulations were conducted by FA and GR. All authors contributed to the writing and editing of the paper.

Competing interests

The authors declare that they have no conflict of interest.

Acknowledgements

We thank the following for stimulating and fruitful discussions: Jacques Magnaudet (Institut de Mécanique des Fluides, IMFT, Toulouse), Yannick Hallez (Laboratoire de Génie Chimique de Toulouse, LGC, Toulouse), Jean-Lou Pierson (Institut Francais du Pétrole Energies Nouvelles, IFPEN, Lyon), Jean-Luc Redelsperger (Laboratoire de Physique des Océans, IFREMER, Brest), Jean-Pierre Pinty and Juan Escobar (Laboratoire d'Aérologie, LA, Toulouse), Odile Thouron, Thibaut Lunet, Luc Giraud, Isabelle d'Ast and Gerard Dejean (Centre Européen de Recherche et Formation Avancée en Calcul Scientifique, CERFACS, Toulouse). The simulations were performed on the Neptune/Nemo-CERFACS supercomputers on the Occigen-CINES Bull cluster (c2016017724 and a0010110079 GENCI projects).

Financial support

Part of this work was supported by Région Midi-Pyrénées funding.

Review statement

This paper was edited by Ignacio Pisso and reviewed by two anonymous referees.