# An improved immersed boundary method for turbulent flow simulations on Cartesian grids: extension of a global geometric approach for thin boundary layers and strong flow incidence

Benjamin Constant, Stéphanie Péron, Heloise Beaugendre, Christophe Benoit

<figure><img src="image_1.png" /></figure>

## To cite this version:

Benjamin Constant, Stéphanie Péron, Heloise Beaugendre, Christophe Benoit. An improved immersed boundary method for turbulent flow simulations on Cartesian grids: extension of a global geometric approach for thin boundary layers and strong flow incidence. Journal of Computational Physics, 2024, 519, pp.113441. 10.1016/j.jcp.2024.113441. hal-04731107

HAL Id: hal-04731107

https://hal.science/hal-04731107v1

Submitted on 10 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
# An improved immersed boundary method for turbulent flow simulations on Cartesian grids: extension of a global geometric approach for thin boundary layers and strong flow incidence

Benjamin Constant<fnref n="1" />, Stéphanie Péron<fnref n="1" />, Héloïse Beaugendre<fnref n="2" />, Christophe Benoit<fnref n="1" />

## Abstract

This paper presents further improvements to the immersed boundary method introduced in [1]. The proposed developments take place during the pre- and post-processing stages of the simulation workflow and aim to further increase the accuracy of the approach for steady-state simulations of high Reynolds number turbulent flows around aerodynamic geometries facing strong incidence. To this end, the location of the forcing points is further optimized prior to simulation, with an adaptive and local modeling height that accounts for the evolution of the turbulent boundary layer thickness, especially at the leading edge. In addition, the direct extrapolation of the pressure solution at the wall is replaced by first- and second-order reconstructions using the normal pressure gradients interpolated at a new set of image points. This second approach is used only in post-processing, after the simulation, and prevents the degradation of the wall pressure in the presence of strong curvatures or thin boundary layers. These developments have been validated by simulating subsonic turbulent flows around the NACA0012 profile, 2D multi-element airfoil (2DMEA), and HL-CRM half-plane at significant angles of attack. Smooth and accurate skin pressure and friction coefficients are observed, in excellent agreement with body fitted wall-resolved solutions, even for coarser Cartesian meshes. Better drag and lift coefficients calculated by direct near-field integrations are obtained, along with accurate predictions of the near-wall flow physics throughout the turbulent boundary layer.

**Keywords:** Immersed Boundary Method, Wall Modeling, Cartesian Mesh, RANS Simulations

<footnote n="1">¹ONERA, Université Paris Saclay, F-92322, Châtillon, France</footnote>

<footnote n="2">²Univ. Bordeaux, INRIA, CNRS, Bordeaux INP, IMB, UMR 5251, F-33400, Talence, France</footnote>
# 1. Introduction

In the aerospace industry, the use of Computational Fluid Dynamics (CFD), along with other analytical and experimental methods, has increased significantly to simulate realistic and complex configurations with more and more fine geometric details. While the growth in computational power over the past decades is undeniable, and even more so with the upcoming expansion of exascale computing, the mesh generation step still represents a major bottleneck in the CFD workflow, as it typically requires manual interaction and expertise [2]. For this reason, the popularity of Immersed Boundary Methods (IBM) has grown exponentially in recent years due to their ability to greatly simplify mesh generation through the use of simple Cartesian grids that do not have to conform to the obstacles [3]. Following the initial work of Peskin [4], many variations have been proposed in the literature for a wide range of flows and applications [5]. Currently, immersed boundaries combined with Cartesian grids and local grid refinement seem to be particularly well-suited for high levels of automation and computational efficiency [6]. However, in order to limit the number of points within the boundary layer in the resulting Cartesian mesh, wall models are often required to solve turbulent flows at high Reynolds numbers [7, 8].

As a matter of fact, the coupling between an immersed boundary treatment and a wall model is not an easy task, especially since the latter tends to generate spurious oscillations in the vicinity of the immersed obstacle. Over the last five years, many authors have offered different solutions to this open problem, both for RANS (Reynolds-Averaged Navier-Stokes) and WMLES (Wall-Modeled Large Eddy Simulation) contexts. Following the initial work of Capizzano [7], Tamaki et al. have worked on a linearization of the velocity, temperature and turbulent viscosity profiles below the first computed points, thus acting directly on the wall model itself to circumvent the strong non-linearities that often occur within the turbulent boundary layer [9, 10]. For their part, Cai et al. have introduced an innovative way to compute the tangential velocity while remaining consistent with the near-wall flow physics, by interpolating this value with different probes found in the wall normal direction [11, 12]. Others have worked on geometric solutions, i.e. directly on the positioning of the IBM points around the immersed geometry to prevent any of them from lying too low inside the turbulent boundary layer. This correction has been proposed systematically [13], according to an absolute distance criterion [14] or according to $y^+$ values measured at the IBM points [15]. In 2021, Constant et al. also offered a geometric solution to this well-known problem, using a global modeling height around the immersed geometry, which allowed them to achieve
proper mesh convergence as well as a complete removal of all near-wall oscillations [1].

Another major shortcoming of the immersed boundary methods is the impossibility to directly analyze the flow fields at the wall boundaries without performing an additional reconstruction step. In this sense, Capizzano [16], followed by Manueco et al. [17], have introduced new and efficient ways to extract the wall information in an immersed boundary context, using watertight surfaces automatically defined near the wall during the pre-processing step to integrate the global loads at a low computational cost. Conversely, Tamaki and Imamura [9] rather considered a far-field integration domain defined by the front of the first computational cell faces, to simplify the computation of the total aerodynamic forces exerted on the obstacle. It is therefore undeniable that a great deal of progress has been made from a geometric point of view, which has considerably facilitated the on-the-fly calculation of aerodynamic loads. However, we believe that certain aspects of near-wall flow modeling can still be improved. In particular, in the aforementioned studies, as well as in most of the recent IBM literature devoted to the simulation of high Reynolds number turbulent flows, the pressure solution is generally assumed to be constant in the wall normal direction and throughout the entire turbulent boundary layer, which can ultimately lead to significant errors and even complicate the direct and accurate integration of the aerodynamic forces on the wall.

To address this shortcoming, we propose to study higher-order extractions of the wall pressure during the post-processing stage using normal pressure gradients interpolated further away from the geometry. The new methodology introduced in [1] is adopted and further improved to account for the real evolution of the turbulent boundary layer around profile geometries, with a more adaptive modeling height. The final objective of this paper is to demonstrate that the immersed boundary method we propose can be a viable alternative to conventional structured approaches, particularly for aerodynamic geometries facing strong flow incidence, while maintaining a realistic mesh size.

The paper is organized as follows: the governing equations and the flow solver are introduced in section 2. The entire IBM workflow is then summarized in section 3, from pre-processing to post-processing, to facilitate both reading and examination of subsequent developments. In section 4, a local modeling height is introduced that allows the user to reduce the portion of the modeled region near the wall in the presence of thin boundary layers, without compromising the stability and accuracy of the IBM methodology. In section 5, higher-order corrections to the wall pressure are adopted, with particular emphasis on the introduction of a new set of image points into our framework. Simulations of subsonic turbulent flows under strong flow incidences around the
NACA0012 profile, 2D multi-element airfoil, and HL-CRM half-plane are studied in section 6 to evaluate the cross-benefits of the two improvements brought to our methodology, resulting in more accurate drag and lift coefficients calculated by near-field integrations.

## 2. Description of the flow solver

### 2.1. Governing equations

The steady-state mean-flow solution is obtained by solving the compressible Navier-Stokes equations:

$$ \left\{ \begin{aligned} & \frac{\partial \rho}{\partial t} + \frac{\partial \rho u_k}{\partial x_k} = 0 \\ & \frac{\partial \rho u_1}{\partial t} + \frac{\partial (\rho u_1 u_j + \delta_{1j} p)}{\partial x_j} - \frac{\partial \tau_{1j}}{\partial x_j} = 0 \\ & \frac{\partial \rho u_2}{\partial t} + \frac{\partial (\rho u_2 u_j + \delta_{2j} p)}{\partial x_j} - \frac{\partial \tau_{2j}}{\partial x_j} = 0 \\ & \frac{\partial \rho u_3}{\partial t} + \frac{\partial (\rho u_3 u_j + \delta_{3j} p)}{\partial x_j} - \frac{\partial \tau_{3j}}{\partial x_j} = 0 \\ & \frac{\partial \rho E}{\partial t} + \frac{\partial (\rho E + p) u_j}{\partial x_j} + \frac{\partial q_j}{\partial x_j} - \frac{\partial \tau_{ij} u_i}{\partial x_j} = 0 \end{aligned} \right. \quad (1) $$

where $\rho$ refers to the fluid density, $u = (u1, u2, u3)$ is the velocity vector, $p$ is the pressure, $E$ is the total energy, $q$ is the heat flux vector, and $\bar{\tau}$ is the viscous stress tensor.

The equations are then averaged to obtain the Reynolds-Averaged Navier-Stokes equations (RANS) and closed by the Spalart-Allmaras (SA) turbulence model [18]. The latter is a one-equation transport model that introduces a modified eddy viscosity $\tilde{\nu}$ such as:

$$ \nu_t = \tilde{\nu} f_{v1}, \quad \text{with} \quad f_{v1} = \frac{\chi^3}{\chi^3 + c_{v1}^3}, \quad \chi = \frac{\tilde{\nu}}{\nu}, \quad \text{and} \quad c_{v1} = 7.1, \quad (2) $$

where $\nu$ and $\nu_t$ are the kinematic viscosity of the fluid and the eddy viscosity respectively. The damping function $f_{v1}$ goes to 0 near-wall and 1 in the logarithmic layer, in such a way that the pseudo viscosity $\tilde{\nu}$ matches the eddy viscosity in the latter region. The SA turbulence model adopted here is the standard model, as originally introduced in 1992 [18].
## 2.2. Numerical methods

The system of equations 1 can be represented in the following matrix form:

$$ \frac{\partial Q}{\partial t} + \nabla \cdot F_c(Q) - \nabla \cdot F_v(Q) = 0, \qquad (3) $$

where $F_c$ and $F_v$ describe the convective and viscous fluxes. $\nabla$ is the usual divergence operator and $Q$ is the vector of conservative variables. These quantities can then be written as:

$$ F_{c_j}(Q) = \begin{pmatrix} \rho u_j \\ \rho u_1 u_j + \delta_{1j} p \\ \rho u_2 u_j + \delta_{2j} p \\ \rho u_3 u_j + \delta_{3j} p \\ (\rho E + p) u_j \end{pmatrix}, \quad F_{v_j}(Q) = \begin{pmatrix} 0 \\ \tau_{1j} \\ \tau_{2j} \\ \tau_{3j} \\ \tau_{ij} u_i - q_j \end{pmatrix}, \quad \text{with } Q = \begin{pmatrix} \rho \\ \rho u_1 \\ \rho u_2 \\ \rho u_3 \\ \rho E \end{pmatrix} \qquad (4) $$

In this work, the system 3 is solved using the structured Navier-Stokes solver from the Fast software suite [19], hereinafter referred to as FastS and developed at ONERA since 2015. This solver is based on a 2nd-order accurate Finite-Volume Method (FVM) for block-structured meshes, for which the flow solution is known at the center of the cells. As stated in [20] and recalled in [21], one major feature of this solver lies in its ability to efficiently handle a broad range of unsteady simulations (URANS, ZDES, DDES or LES), by allowing the update of up to 10 million cells per second on a single Intel Broadwell core. These performances are made possible by the use of a hybrid MPI/OpenMP parallelization, coupled with various advanced HPC techniques (cache blocking, vectorization, etc.) thoroughly described by Alferez et al. [20]. FastS also features a dedicated Cartesian solver, 2.5 times faster than its curvilinear counterpart in terms of memory and CPU cost, as it simplifies the numerical fluxes and eliminates the need to store most of the mesh-related information. Ultimately, its use counterbalances the presence of a large number of points in our multi-block Cartesian meshes, which can be of the order of $10^9$ for three-dimensional industrial applications.

For RANS computations, two numerical schemes are considered to estimate the convective fluxes, depending on the flow regime. In practice, a reduced AUSM (Advection Upstream Split-
ting Method) scheme is applied when the Mach number is lower than 0.2. This first scheme is a modification of the AUSM+(P) scheme developed by Edwards & Liou [22], deprived of its shock capture capabilities and solely adapted to quasi-incompressible subsonic flows [23]. For higher Mach numbers, the AUSM scheme is then replaced by the classical Roe scheme [24], which is more dissipative. In the following, a 3rd-order accurate MUSCL reconstruction [25] is also taken into account to accurately recover the left and right states at each interface.

The viscous fluxes are discretized using a 2nd-order finite volume scheme. To this end, the approach developed by Thomas & Walters [26] is followed, which allows a direct evaluation of all gradient information at the interfaces in a robust and accurate way, resulting in a more compact stencil in spite of a higher cost. The use of this approach requires the knowledge of the velocity and temperature gradients at the interfaces of each cell, in the three directions of space. Figure 1 illustrates the stencil used to evaluate the convective, viscous, and total flux balances for all the faces composing an elementary volume $Ω_k$. For 3D applications, it can be seen that the entire stencil consists of 25 cells.

<figure><img src="image_3.png" /><figcaption>Figure 1: Sum of stencils used to evaluate both convective and viscous flux balances for a given elementary volume $Ω_k$ (here drawn in red).</figcaption></figure>

FastS supports different time integration methods, implicit [27] or explicit [28], to suit different types of numerical simulations, unsteady or not. Since the present work focuses only on steady-state RANS simulations, the temporal integration is performed using a simple and robust 1st order accurate implicit Euler scheme together with a local time step based on a CFL number imposed beforehand.
For the simulation of a turbulent flow in a free environment (cf. Figure 2), a far-field boundary condition is classically considered at the outer edges of the computational domain. This non-reflecting boundary condition is computed in FastS by a characteristic reasoning according to Thomson’s simplified approach [29]. Finally, all boundary conditions involving a solid wall are systematically handled by an immersed boundary condition (IBC).

<figure><img src="image_4.png" /><figcaption>Figure 2: Set of boundary conditions for a given turbulent flow simulation around a profile geometry in a freestream environment. The far-field boundary conditions are represented in red ($\Gamma_{FF}$), while the immersed boundary condition applied to the obstacle is displayed in black ($\Gamma_W$). Finally, the exterior surfaces used for a far-field integration of the aerodynamic forces are shown in green ($\Gamma_F$).</figcaption></figure>

To conclude this section, it is important to note that FastS actually works with the primitive variables rather than the conservative ones. This choice is driven by HPC considerations and aims at maximizing the number of multiplications over divisions, which are more computationally expensive. Thus, from now on $W = (\rho, u_1, u_2, u_3, T, \tilde{\nu})$ will denote the flow solution vector, to which is added the pseudo-viscosity derived from the Spalart-Allmaras turbulence model.

### 3. Description of the immersed boundary method

As recalled in the introduction, immersed boundary methods allow the use of automatic and fast mesh generation techniques such as the one introduced by Peron and Benoit [30]. Their meshing strategy is adopted here, which consists in generating of a set of structured uniform Cartesian grids with increasing levels of refinement. The different levels are entirely managed by an octree skeleton, which guarantees the correct mesh resolution near the walls while limiting the final number of cells inside the mesh. The connectivity between adjacent structured blocks is handled with minimal overlap between said blocks and a second-order Chimera method that performs inter-grid interpolations during the simulation [31, 32]. The resulting Cartesian mesh is thus structured by
blocks, allowing the use of an efficient structured solver such as the one provided by FastS.

The immersed boundary method studied and extended in this paper is based on a ghost-cell direct forcing approach. It consists of imposing the flow solution **W** at some preselected points around the obstacle, called target or forcing points. The solution at these points is usually obtained with the help of additional image points, whether fictitious or not, and positioned further away in the fluid region. In the literature, the latter can also be referred to as reference points [8].

To be fully consistent with the FastS Cartesian solver, which is based on a five-point stencil in each direction, at least two layers of target points must be defined near the obstacle so that an immersed boundary condition can be properly imposed. This is highlighted in Figure 3-(a), which shows the stencil involved in the evaluation of all numerical fluxes for a given calculated point *P*. It is also important to note that all target points are positioned outside the obstacle. This choice is motivated by the use of a wall model in combination with IBM, for RANS or hybrid RANS/LES simulations.

<figure><img src="image_5.png" /><figcaption>Figure 3: Illustration of the immersed boundary method implemented in Cassiopee, for RANS or hybrid RANS/LES simulations.</figcaption></figure>

Figure 3-(b) illustrates the reconstruction of the flow solution at a target point *A* (●) and a wall point *W* (○) from their corresponding image point *B* (○). Since the image points are not discretization points in the current method, the solution at these points is simply recovered by second-order linear interpolations from the surrounding computed donor points. The wall points, which are located directly on the immersed obstacles, are simply used to store all the information at the wall during the simulation for post-processing purposes. Like the image points, the wall points
are not actually mesh discretization points, and an additional operation is required to visualize the flow solution on the geometry surface.

Ultimately, the entire workflow of the immersed boundary method implemented within the Cassiopee/Fast framework can be divided into three main steps, summarized in the following Figure 4. It can be seen that the immersed boundary treatment is completely handled by Cassiopee [33, 34], from pre-processing to post-processing. During the simulation, the solver actually communicates with Cassiopee's modules to update both IBM target points and Chimera points, which are located inside the mesh, either near the immersed boundary (IBM target points) or in the overlapping regions between multiple Cartesian grids (Chimera points), and act as ghost cells from the solver's point of view [35]. Thus, although the current immersed boundary method is optimized for Finite Volume solvers such as FastS, this approach can theoretically be applied to many other situations, in both structured and unstructured contexts.

<figure><img src="image_6.png" /><figcaption>Figure 4: IBM simulation workflow within the Cassiopee/Fast environment. Cassiopee (---), Fast (----).</figcaption></figure>

The following three subsections summarize the specific IBM treatment that has been developed in Cassiopee by Peron *et al.* [36]. This summary also includes more recent developments by Constant *et al.* [1] and is important for a better understanding of sections 4 and 5, since some of the information recalled below has not been thoroughly introduced in previous articles.
### 3.1. IBM pre-processing

The pre-processing step begins with the automatic generation of an octree-based Cartesian mesh, where the resulting mesh is a collection of overset isotropic grids with minimal overlap. The detailed sub-steps are shown in Figure 5, from the generation of the quadtree skeleton around a spherical obstacle (Figure 5-(a)), through the resulting Cartesian grids (Figure 5-(b)), to the merging of elementary grids to increase the amount of intra-block communication in a shared memory context and thus take full advantage of the hybrid MPI/OpenMP parallelization of FastS (Figure 5-(c)). A hole-cutting algorithm is applied to blank interior cells, which also removes fully blanked blocks for memory reasons (Figure 5-(d)). A cell nature field is then used to store the indices of interior blanked cells (*cellN*=0), calculated cells (*cellN*=1), and interpolated cells (*cellN*=2).

<figure><img src="image_7.png" /><figcaption>Figure 5: Illustration of an automatic mesh generation around a circle. The parameter $d_{ff}$ defines the global domain extent, from the geometry bounding-box. $h_{min}$ is the minimum near-wall spacing. $v_{min}$ is the minimum number of cells per direction that compose each Cartesian grids of the multi-block mesh.</figcaption></figure>
The meshing operation is immediately followed by the IBM pre-processing stage. Using the original algorithm described in [36], the first two layers of cells near the obstacle are first marked as target points (Figure 6-(a)). In order to find suitable candidates for the image point that are both close enough to remain in the boundary layer and far enough to interpolate the flow solution only from interior computed points, a watertight surface is created from the stair-step fringe of the first computed cell centers (Figure 6-(e)). The image points are simply obtained by projecting the target points onto this surface front, following the direction of the wall normal, while the corresponding wall points are obtained by the same projection in the opposite direction (Figure 6-(f,g)).

In a previous paper [1], we have also proposed an improved algorithm for the location of the target and image points, based on geometrical considerations and absolutely independent of the chosen turbulence and wall models. This approach aimed at removing the spurious oscillations that could occur on the skin pressure and friction coefficients with the original IBM treatment developed by Peron et al. [36], as well as many others [37, 11]. To that end, a value of $y_{target}^{+}$ is provided for the first computed cells, ideally set between 50 and 500 to match the validity range of the algebraic wall function and the chosen minimum near-wall spacing. A flat plate approximation is then used to evaluate the distance of the first computed cells to the wall, from the Reynolds number $Re$, the characteristic length $L$ of the geometry, and the flow conditions. The resulting modeling height is denoted as $H_{mod}$ and is calculated from the following formula:

$$H_{mod} = y_0 \times y_{target}^{+} = \sqrt{2} \frac{y_{target}^{+} L}{Re \sqrt{f(Re)}}, \quad (5)$$

where $y_0$ is the height of the first cell on a flat plate for which $y^+ = 0$ and $f(Re) = 0.058 Re^{-0.2}$. This new methodology primarily ensures that the mesh is fine enough at the first computed cells to correctly capture the strong variations of the tangential velocity in the wall normal direction. By extending the modeling area around the obstacle, it also encompasses neighboring target points in the same region of the turbulent boundary layer, which greatly contributes to the reduction of spurious oscillations in the vicinity of the geometry.
<figure><img src="image_8.png" /><figcaption>Figure 6: Description of the IBM pre-processing workflow. Original algorithm (a,e,f,g) and additional steps (b,c,d) introduced by the improved methodology F42 (cf. [1]).</figcaption></figure>
In practice, this improved algorithm modifies the IBM pre-processing only after the first two layers of cells at the edge of the immersed obstacle have been marked as target points (Figure 6-(a)). After calculating $H_{\text{mod}}$ from Equation 5 (Figure 6-(b)), all points located below this threshold are now marked as target points (Figure 6-(c)). Since only the last two layers of target points are taken into account by the flow solver, a blanking operation can also be performed to reduce both the memory consumption and computational time (Figure 6-(d)). With this last step added to our framework, little to no difference was found between the original and the improved methods in terms of CPU performance, as stated in [1]. The rest of the IBM pre-processing, which is identical to the original IBM treatment, can then proceed (Figure 6-(e,f,g)).

## 3.2. Wall modeling in an immersed boundary context

<figure><img src="image_9.png" /><figcaption>Figure 7: Left: Asymptotic behavior of a turbulent boundary layer developing on a flat plate at equilibrium. Comparison with the analytical wall laws of Musker [38] and Spalding [39]. Right: Direct forcing approach between a target point A and its corresponding image point B.</figcaption></figure>

In the final step of the IBM reconstruction, the flow solution at the target points must be updated using the information available at their corresponding image points. This operation is illustrated in Figure 7-(right). Let us denote $\hat{W}_B = (\rho_B, u_B, T_B, \tilde{\nu}_B)$, the flow solution interpolated at a given image point B, and $W_A = (\rho_A, u_A, T_A, \tilde{\nu}_A)$, the flow solution reconstructed at a given target point A. To simplify the following notations, let us also define $u = u \cdot t$ and $v = u \cdot n$ as the normal and tangential components of a velocity vector $u$.

The current immersed boundary method discussed in this paper is mainly dedicated to aero-dynamic simulations at Reynolds numbers of the order of a million. Therefore, the close-wall flow physics cannot be solved directly with Cartesian grids, as this would lead to prohibitive cell counts.
To overcome this problem, the immersed boundary methods usually operate with wall functions to reconstruct the velocity at the target points located in the turbulent boundary layer, either based on algebraic models [8, 36, 37] or TBLE (Thin Boundary Layer Equation) based models [40, 41, 42].

In this framework, the algebraic wall function proposed by Musker [38] is adopted to recover both the friction velocity $u_τ$, which is assumed to be constant in the wall normal direction, and the tangential velocity at the target points. It is worth noting that this wall function, taken in its equilibrium version, is designed to fit the viscous and buffer sublayers in addition to the logarithmic sublayer. This aspect is demonstrated in Figure 7-(left).

Given the tangential velocity vector at the image point B, the friction velocity is obtained through a Newton-Raphson method on Musker's algebraic wall function<fnref n="3" />. The tangential velocity $u$ at the target point A then reads:

$$u_A = u_τ M(y_A^+), \quad \text{where} \quad y_A^+ = \frac{u_τ y_A}{ν}, \qquad (6)$$

while the normal velocity $v$ at the target point A is simply estimated by linear reconstruction, assuming its value to be zero at the wall:

$$v_A = \frac{y_A}{y_B} v_B. \qquad (7)$$

Then, since the definition $ν_t = κu_τy_D$ is valid in the log layer and down to the wall, the pseudo-viscosity $\tilde{ν}$ defined in Equation 2 can be evaluated at the target point A by solving the following nonlinear equation:

$$(\tilde{ν}_A)^4 - κu_τy_A D(\tilde{ν}_A)^3 - κu_τy_A Dν^3 C_{v1}^3 = 0, \qquad (8)$$

where $κ = 0.41$ is the Von Kármán constant and $D = [1 - \exp(-y_A^+/19)]$ is the Van Driest damping

<footnote n="3">³A strict convergence criterion is applied to both the function itself and its derivative. In practice, this iterative algorithm converges in a few steps.</footnote>
term. The resolution of Equation 8 by the Ferrari’s method is introduced and thoroughly described in [36]. To account for low compressibility effects within our wall model, a Crocco-Busemann relation [43] is used to reconstruct the temperature:

$$T_A = T_B + \frac{\sqrt[3]{Pr}}{2C_P}(u_B^2 - u_A^2), \qquad (9)$$

where Pr is the Prandtl number and $C_P$ is the specific heat coefficient at constant pressure. It is important to note, however, that Equation 9 is only valid for adiabatic walls. Knowing these quantities, the density is finally determined from the ideal gas law:

$$\rho_A = \frac{p_A}{RT_A}, \qquad (10)$$

where R is the ideal gas constant and by considering the pressure p constant along the wall normal direction within the turbulent boundary layer:

$$p_A = p_B. \qquad (11)$$

This last hypothesis will be discussed in the next sections. As recalled before, the IBM reconstruction is immediately followed by the update of all Chimera points, either for 1-to-1 or 1-to-2 abutting grids. The complete workflow of an IBM simulation using Cassiopee/Fast is detailed in Figure 8, which illustrates the case of two overlapping grids near the immersed boundary, with different refinement levels. In this situation, the fine and the coarse grids are extended by two and three cells respectively, to ensure that the donor points are always interior computed points and avoid any kind of extrapolation.

This figure is also useful for understanding how the two grids interact at each substep, especially in a situation as complex as this one. It is clear from studying this figure that the Chimera points cannot be donor points for IBM points, be they image or target points. Therefore, the donor points for some image points can be positioned in adjacent grids, while the IBM treatment must be performed systematically before the Chimera transfers.
<figure><img src="image_10.png" /><figcaption>Figure 8: Workflow of an IBM simulation performed in Cassiopee/Fast for 1-to-2 abutting grids close to a wall</figcaption></figure>
### 3.3. IBM post-processing

<figure><img src="image_11.png" /><figcaption>Figure 9: Description of the IBM post-processing workflow.</figcaption></figure>

As detailed in [36] and recalled in [1], the reconstruction of the surface solution is performed with the help of a third-order accurate Moving Least Squares (MLS [44, 45]) interpolation on the triangular mesh defining the obstacle, from the flow quantities computed and stored at the wall points at each iteration of the flow solver. This operation is illustrated in Figure 9.

During the post-processing step, each wall point is orthogonally projected onto the set of the $nT$ triangles $(T_m)_{m \in [1, nT]}$ that form the triangulation $T$ of the geometry. For a vertex $V_n$ in $T$, a minimum set of 10 wall points is considered to perform the MLS reconstruction, chosen from the adjacent $T_m$ triangles. To avoid the introduction of interpolation errors, the final triangulation of the immersed obstacle must also maintain a certain degree of consistency with the discretization of the Cartesian grids positioned in the vicinity of the obstacle.

After this important step of the immersed boundary post-processing, it is immediately possible to calculate the skin pressure and the coefficient of friction at the geometry surface:

$$C_p = \frac{p_w - p_\infty}{\frac{1}{2} \rho_\infty u_\infty^2} \quad \text{and} \quad C_f = \frac{\tau_w}{\frac{1}{2} \rho_\infty u_\infty^2} = \frac{\rho_w u_\tau^2}{\frac{1}{2} \rho_\infty u_\infty^2}, \qquad (12)$$

where $p_w = p_B$, $T_w = T_B + \frac{\sqrt{Pr}}{2C_p}u_B^2$, and $\rho_w = \frac{p_w}{RT_w}$. The drag and lift coefficients are then defined as follows:
$$C_D = \frac{F_D}{\frac{1}{2} \rho_\infty U_\infty^2 A_{ref}} \quad \text{and} \quad C_L = \frac{F_L}{\frac{1}{2} \rho_\infty U_\infty^2 A_{ref}}, \qquad (13)$$

where $F_D$ and $F_L$ are obtained by integrating all the aerodynamic forces acting on the body, parallel and perpendicular to the flow directions, respectively.

According to the work of Destarac [46], this integration can be done in two different ways, either by a near-field integration, which consists in integrating the forces directly on the body surface $\Gamma_W$, or by a far-field integration, which involves selecting a control volume surrounding the body to integrate the aerodynamic forces on its outer surfaces $\Gamma_F$ (cf. Figure 2). The next paragraphs describe our methodology for both approaches.

### 3.3.1. Near-field integration of the aerodynamic forces

The near-field approach for the integration of the aerodynamic forces on a body surface $\Gamma_W$ reads:

$$F_{i,W} = \int_{\Gamma_W} \left( (p - p_{\infty})\delta_{ij} - \tau_{ij} \right) n_j \, dS. \qquad (14)$$

In an immersed boundary context, this approach is particularly direct, since the wall solution is already computed at each iteration of the flow solver and stored at the wall points. Only the stress tensor needs to be recovered in post-processing after interpolating the flow quantities at the wall boundaries using the methodology described in the previous section.

Let us denote $n$, the unit normal vector to the wall, and $u$, the velocity vector defined at the target points. The unit tangent vector to the surface then becomes:

$$t = \frac{u - (u \cdot n)n}{\|u - (u \cdot n)n\|}. \qquad (15)$$

Knowing the wall shear stress $\tau_w$ from the wall density $\rho_w$ and the friction velocity $u_\tau$, and assuming that the shear force is equal to $\tau_w t$, the stress tensor $\bar{\tau}$ can be simplified as follows:
295

The pressure and friction forces are then integrated separately on the surface $\Gamma_W$ to quantify the effect of each component on the resulting drag and lift coefficients.

300

Thus, the key to obtaining the correct drag and lift coefficients by near-field integration is to have smooth and accurate skin distributions around the obstacle, which depends heavily on the near-wall treatment as well as the mesh refinement. Another disadvantage of this approach may be the inability to decompose the aerodynamic forces other than mechanically, i.e. with a pressure part and a viscous part. On the other hand, far-field integration allows a physical decomposition of the total drag into different components: spurious drag, wave drag, and lift-induced drag, which can be very useful for improving the performance of an aircraft [47]. For this reason, a considerable amount of research has been carried out over the last twenty years to improve our understanding of the underlying phenomena of the aerodynamic forces acting on a body under different flow conditions, as well as to provide an accurate decomposition of these forces by far-field integration [48, 49].

305

### 3.3.2. Far-field integration of the aerodynamic forces

The far-field approach for the integration of the aerodynamic forces on an exterior surface $\Gamma_F$ reads [50]:

310

First, we need to define a control around the obstacle or the set of obstacles that defines the immersed boundaries. In our framework, we start by creating a square box of length $\delta$, whose center is chosen as the barycenter of the obstacle, as shown in Figure 2. To minimize the effect of the choice of $\delta$ in the final integration, we impose $\delta = 2c$, where c is the chord length. In fact, for all the bidimensional cases tested so far, it has been found that the final integration does not

$$ \bar{\bar{\tau}} = \begin{pmatrix} \tau_{xx} & \tau_{xy} & \tau_{xz} \\ \tau_{yx} & \tau_{yy} & \tau_{yz} \\ \tau_{zx} & \tau_{zy} & \tau_{zz} \end{pmatrix} = \begin{pmatrix} 2\tau_w t_x n_x & \tau_w (t_x n_y + t_y n_x) & \tau_w (t_x n_z + t_z n_x) \\ \tau_w (t_x n_y + t_y n_x) & 2\tau_w t_y n_y & \tau_w (t_y n_z + t_z n_y) \\ \tau_w (t_x n_z + t_z n_x) & \tau_w (t_y n_z + t_z n_y) & 2\tau_w t_z n_z \end{pmatrix} \quad (16) $$
strongly depend on the choice of this distance parameter once the condition $\delta > c$ is satisfied.

Once the box is created, the steady-state flow solution can be interpolated directly from the volume mesh to the exterior surfaces of the box for the final computation. Provided that the far-field integration is performed far enough from the wall, especially downstream, the viscous stress can actually be ignored, as stated in [47]. This last hypothesis greatly simplifies the computation, since the reconstruction of the stress tensor at such a distance from the wall can be the source of many numerical errors, due to poor gradient evaluation in these distant regions. Equation 17 finally becomes:

$$F_{i,F} = - \int_{\Gamma_F} \left( \rho(u_i - u_{\infty,i})u_j + (p - p_{\infty})\delta_{ij} \right) n_j dS. \quad (18)$$

Ideally, the far-field integration should be performed directly on the original mesh, following a surface of the same level of refinement, to avoid numerical noise as much as possible. Successful attempts in this direction have been observed in the recent literature [9]. However, it can be quite difficult to create such an optimal front sufficiently far away from the wall to avoid as much as possible the influence of the immersed boundary treatment in the final integration, mainly because this type of closed watertight surface cannot be automatically retrieved from the positioning of the IBM points alone, and can easily cross grids of different refinement levels. In addition, the minimal overlap between adjacent Cartesian blocks is another difficulty that would need to be addressed in future work to avoid the tedious interpolation step in the current implementation.

# 4. An improved pre-processing for the simulation of high Reynolds number turbulent flows in an immersed boundary context

## 4.1. Motivation

While the improved F42 methodology proposed in 2021 [1] has been successful in removing near-wall oscillations for turbulent flow simulations on Cartesian meshes, a limitation has eventually emerged regarding the use of a single modeling height over the entire geometry. Recent work has shown that the modeling height cannot remain invariant for all flow configurations due to the large fluctuations observed in the thickness of the turbulent boundary layer. This is particularly true for airfoils facing strong flow incidences, where this threshold height should not be imposed
identically at the leading edge and further downstream, as the turbulent boundary layer evolves along the chord. In this section, an extension of the F42 approach is proposed to address the specific constraint imposed by thin turbulent boundary layers, without compromising the universality of this geometric methodology.

The standard flat plate approximation for the skin friction coefficient can sometimes lead to strong over-predictions of the $y^+$ distribution at the target and image points. As shown in Figure 10, the initial estimation of the modeling height is usually satisfactory around airfoil at zero angle of attack (Figure 10-(a)). For higher angles of attack (Figure 10-(b,c)), the approximation is still valid on the pressure side, where the flow behavior is close to that of a flow past a flat plate. On the suction side, however, the $y^+$ values captured at the image points differ significantly from the initially imposed $y_{target}^+$ value, especially at the leading edge. For example, the relative error between the $y_{target}^+ = 300$ value and the actual maximum measured at the leading edge for the corresponding solution is close to 300% at AoA=15°.

<figure><img src="image_12.png" /><figcaption>Figure 10: NACA0012 airfoil (Ma 0.15 - Re 6M): $y^+$ distributions at the image points for increasing angles of attack. Impact of the initial choice of $y_{target}^+ \in [100, 300]$.</figcaption></figure>

The discrepancies observed near the leading edge can be explained by the fact that the turbulent boundary layer is not fully developed in these regions. As a result, not locally adjusting the modeling height to actually reflect the physical reality of the near-wall physics carries the risk of applying the wall model outside its validity domain, potentially leading to significant numerical errors, even though the surface solution remains free of oscillations [1].
## 4.2. Modeling height adaptation

To solve this problem, the present work proposes a transition from a global to a local modeling height, while maintaining a good regularity with respect to the positioning of the front of the first computed points. After the convergence of a previous RANS simulation performed with the F42 approach and a global modeling height, an additional step is performed within the IBM workflow described in Figure 6, during which the $y^+$ values stored at the target point are extracted towards the Cartesian mesh. The collected information is then used to locally adjust the position of the target points in the vicinity of the immersed obstacle based on a comparison with the initial $y_{target}^+$ value. After this local update of the cell nature field, a small part of the IBM preprocessing must be repeated a second time to be able to run another simulation. This additional pre-processing step mainly involves finding all the interpolation data for the Chimera transfers and the IBM reconstructions prior to the simulation, but excludes any modification of the mesh as well as another costly computation of the turbulent wall distance. Since this entire operation is performed directly on the solution tree, the simulation can naturally restart from the previous converged solution without requiring mesh-to-mesh interpolations of all flow fields. The advantages of this warm start will be discussed later.

To illustrate this solution, a diagram of the new adaptation process is proposed in Figure 11, for a simplified one-dimensional mesh. Starting from the initial positioning of the target points according to a global modeling height $H_{mod}$ calculated a priori, a new distribution is finally obtained that takes into account the real variation of the turbulent boundary layer thickness. In practice, some target points are simply converted to calculated points inside the mesh, based on the $y^+$ value found one layer below, in the wall normal direction. The purpose of this safety is to limit the final adaptation in some regions of the turbulent flow where the $y^+$ distribution exhibits strong variations, so that this new process does not counteract the initial contribution of the F42 approach. It is also important to mention that the nature of the first two rows of target points remains unchanged, even though these points may locally satisfy the adaptation condition. This is a requirement of the flow solver, as recalled in the previous sections.

Let us then define:

$$y_{n-1}^{+} = y_{n}^{+} \cdot \frac{y_{n-1}}{y_{n}} \approx y_{n}^{+} \cdot \left(1 - \frac{h}{y_{n}}\right), \qquad (19)$$
where *h* is the minimum cell size within the mesh. The proposed approximation of the ratio $y_{n-1}/y_n$ is based on the favorable scenario where the Cartesian mesh is perfectly aligned with the obstacle.

<figure><img src="image_13.png" /><figcaption>Figure 11: Algorithm for adapting the position of the IBM target point using $y^+$ information.</figcaption></figure>

For many cases, such as those presented in [1], the initial flat plate approximation provides a satisfactory estimation of the flow behavior that does not tend to over-estimate the boundary layer thickness. However, for an equal number of cases, if not more, we believe that the sensitivity of the wall model in the vicinity of thin boundary layers justifies this local adaptation of the modeling height. Recent developments have also succeeded in limiting the cost of this approach, which does not require another complete IBM pre-processing workflow and allows the user to easily restart a simulation from the previous steady-state solution. The cost/benefit ratio of this adaptation is discussed in detail in the validation section.

Finally, unlike other solutions in the literature, such as Xu *et al.* [15], which propose to dynamically modify the location of the forcing points around an immersed boundary, the methodology introduced in this section is purely static and actually converges in a single pass.
# 5. An improved post-processing for the simulation of high Reynolds number turbulent flows in an immersed boundary context

## 5.1. Motivation

In previous work [51, 36, 1], the pressure solution has been systematically extrapolated from the image points, as the pressure is classically assumed to be constant in the wall normal direction throughout the referenced literature devoted to the simulation of weakly compressible turbulent flows. This hypothesis is also one of the main assumptions underlying the definition of the classical logarithmic law on which Musker's analytical law is based, and the reason why the pressure solution at the target points should remain directly extrapolated from the image points, i.e. to be absolutely consistent with the reconstruction of the velocity profile near the wall. However, the situation changes when it comes to recovering this solution at the wall points.

In this paper, we propose not to include this modeling hypothesis in the IBM post-processing stage and instead consider different higher-order extractions of the wall pressure using the pressure gradient information computed at the image points. As will be shown in the following sections, the zero-pressure-gradient assumption can indeed be quite erroneous in the vicinity of curved geometries and where the turbulent boundary layer is not yet fully developed, especially on the suction side of airfoils when facing strong flow incidence. In these regions, the pressure within the turbulent boundary layer varies significantly in the wall normal direction regardless of the distance to the geometry, which can ultimately lead to an incorrect prediction of the suction peak location. This error typically increases with the distance of the IBM points from the geometry, as determined by the near-wall mesh spacing alone or by the choice of a $y_{target}^{+}$ for the first computational cells, according to the various IBM methods at our disposal.

After numerous observations and tests, we have come to the conclusion that imposing the same pressure at the wall and at the image points may prevent an accurate evaluation of the pressure drag with a near-field integration, even if the deviations on the skin pressure distribution appear to be small. This would be especially the case in subsonic regimes where the pressure peak is intense and strongly confined to the leading edge as soon as the flow incidence increases. More than the correct maximum value itself, the correct location of this point seems to be a key element for a successful pressure drag assessment.

As an example, Figure 12 highlights the errors obtained on the maximum pressure value and the stagnation point location for increasing $y_{target}^{+}$, resulting in large errors in terms of pressure drag
prediction. The test-case is the simulation of a subsonic turbulent flow with a ten-degree angle of attack around a NACA0012 profile, where body-fitted wall-modeled solutions are compared with a reference body-fitted wall-resolved solution.

<figure><img src="image_14.png" /><figcaption>Figure 12: NACA0012 airfoil (Ma 0.15 - Re 6M - AoA 10°): evolution of the relative error on the total pressure drag and the pressure peak (location and intensity) for increasing modeling heights.</figcaption></figure>

Conversely, further investigation of this test-case has shown that a far-field integration would provide a better estimation of the pressure drag. While these observations are consistent with those of Cai *et al.* [11], our conclusions diverge. In fact, by closely examining the evolution of the pressure solution across the turbulent boundary layer, and in the absence of spurious oscillations near the wall, the pressure clearly appears to be degraded during the post-processing step, due to the current reconstruction of this flow variable. It has also become clearer that the current skin pressure distribution does not fully reflect the quality of the solution obtained further into the fluid domain, which shows a better agreement with the reference solution, ultimately leading to a better prediction of the integrated forces further away from the wall.

Moreover, based on the results published in [1] for the NACA0012 test-case on increasingly finer Cartesian meshes while maintaining the same modeling height ($y_{target}^{+} = 100$), it is actually possible to demonstrate how little impact a better reconstruction of the pressure solution at the target points would have during the simulation, other than the risk of jeopardizing the stability of the methodology as a whole, as well as increasing the computational cost of each iteration. Thus, on the left side of Figure 13, the evolution of the relative error on the pressure solution shows no significant difference from one solution to another, although the maximum distance ($d_{max}^{PI}$) at which
the zero-pressure-gradient assumption is applied is divided by four between the coarsest and the finest meshes, as illustrated on the right side of the same figure. In addition, the pressure profiles of these three simulations overlap in the vicinity of the suction peak, providing further evidence for a modification of the post-processing step alone.

<figure><img src="image_15.png" /><figcaption>Figure 13: NACA0012 airfoil (Ma 0.15 - Re 6M - AoA 10°): evolution of the relative error on the total pressure drag and the pressure peak (location and intensity) for increasing near-wall resolutions. Each IBM solution is obtained on a Cartesian mesh, using the F42 approach as well as a constant modeling height based on $y_{target}^{+} = 100$.</figcaption></figure>

## 5.2. Higher-order wall pressure reconstructions

To improve our IBM methodology for the post-processing stage, a second layer of image points has recently been added to our framework to extrapolate the wall pressure with first and second order formulations. The idea is to recover both the pressure and pressure gradient information by interpolating the flow solution further away from the wall and the modeled region to mitigate the influence of the target points on the computation of a corrected wall pressure. In fact, since the pressure at the target points is extrapolated directly from their corresponding image points, the pressure gradient distributions at the first cells are actually inaccurate and even oscillate when strong pressure gradients occur, which would eventually have a negative impact on the skin pressure distribution.

Instead, in order to find suitable donor points for the interpolation of accurate pressure information, we decided to compute a minimum normal distance to the wall to safely position these new image points. In line with the F42 approach, this distance was finally chosen as the modeling height $H_{mod}$, to which is added a smaller distance $\epsilon$, empirically defined as follows:
$$ \epsilon = 2 \cdot \sqrt{\dim \cdot h}, \qquad (20) $$

where *dim* is the dimension of the test-case and *h* is the minimum near-wall spacing. From there, two different approaches have been considered for the final location of the second layer of image points (●), summarized in the following figure:

470

<figure><img src="image_16.png" /><figcaption>Figure 14: Location of the target and wall points around a curved geometry, with their corresponding first and second image points. Iso-distance (left) and stair-step (right) fronts for the new set of image points.</figcaption></figure>
In order to observe the most uniform pressure gradient distribution possible, this second front of image points can be created at a given iso-distance from the wall, following the orange dashed lines plotted in Figure 14-(a). The main challenge with this first choice is the determination of a global distance, which can be particularly difficult to determine when a change in resolution occurs around the geometry. Under these conditions, the mesh parameter used to calculate the maximum projection distance (cf. Equation 20) tends to become a major issue. On the one hand, assuming a constant $h$ runs the risk of bringing the front too close or too far away from the wall, as it transitions from one mesh resolution to another. On the other hand, considering a local $h$ per zone would require a regularization of the overlapping regions between Cartesian grids of different resolutions, which would prove particularly difficult to implement, especially around complex geometries.

Conversely, it is also possible to consider a stair-step front (Figure 14-(b)), based on the same algorithm developed by Person et al. [36] for the projection of the first image points onto the fringe of the first computed cells, here recalled by blue dashed lines. This second choice cleverly avoids all the projection difficulties mentioned above, making it both more robust and more efficient, even with the use of a local mesh parameter $h$ in the vicinity of a mesh resolution switch near the wall. Furthermore, in order to highlight the minimal contribution of a perfectly smooth positioning of the second image points, against the gains in terms of stability and pre-processing time offered by the second approach and its irregular stair-step front, the two approaches are finally compared in the validation section. After this brief introduction of the new set of image points for this modified IBM post-processing, the different types of wall pressure reconstruction for increasing orders will be presented in the following paragraphs.

### 5.2.1. 0th-order formulation: direct extrapolation

As mentioned above, the default reconstruction of the wall pressure corresponds to a zeroth-order extrapolation of the solution interpolated at the first image points (Figure 15-(a)). Without considering any pressure gradient information, the number of interior calculated points involved in this reconstruction is inherently equal to the number of donor points (●) used to interpolate the solution at the image point, with up to four donors in 2D and twice as many in 3D. This first formulation then reads, for a given wall point $W$ and its corresponding first image point $B_1$:

$$p_{w,o0} = p_{B_1} \qquad (21)$$
### 5.2.2. 1st-order formulation: linear reconstruction

In order to propose a first-order extrapolation of the wall pressure, with the aim of taking into account the local variations of this solution in the wall normal direction, we then propose to consider the pressure and pressure gradient information interpolated at the second image points (Figure 15- (b)). In practice, this is done by using the same gradient calculation method employed by the solver, this time applied to the pressure solution at the closest calculated donor points surrounding each second image point. For this second formulation, which includes the first pressure derivatives, it can be seen that the number of donor points now exceeds the initial four calculated points. The first-order reconstruction of the wall pressure between a given wall point W and its corresponding second image point $B_2$ is performed according to the following expression:

$$p_{w,o1} = p_{B_2} - \frac{\partial p}{\partial \mathbf{n}}|_{B_2} \cdot \overrightarrow{||WB_2||}, \quad (22)$$

where $\partial_n p$ is the first directional derivative of the pressure in the wall normal direction, such as:

$$ \begin{aligned} \partial_n p &= \frac{\partial p}{\partial \mathbf{n}} = \mathbf{n}^T \cdot \nabla p \\ &= n_x \frac{\partial p}{\partial x} + n_y \frac{\partial p}{\partial y} + n_z \frac{\partial p}{\partial z}. \end{aligned} \quad (23) $$

Alternatively, it would be possible to estimate the value of $\partial_n p$ using a Finite Difference approach between the two image points, thus reducing the computational cost of a direct evaluation. However, for the purposes of this study, this solution cannot be considered entirely satisfactory. In fact, as recalled earlier, the first image points are projected onto the front of the first computed points and appear too close to the modeled area for their interpolated pressure solutions to be unaffected by this near-wall region. The distance between the first and second image points may also appear relatively large, especially for three-dimensional applications or after adjusting of the modeling front as described in the previous section. In addition, it is important to note that this deviation tends to increase linearly with the near-wall resolution, which may eventually introduce even more numerical errors. For all these reasons, we have chosen to calculate the pressure gradient information only at the second image points, without taking into account the first image points.
### 5.2.3. 2nd-order formulation: quadratic reconstruction

In some cases, it may also be necessary to include second-order pressure gradient information to capture slight curvature effects. For example, if the near-wall resolution is too coarse and the second image point is too far away from the geometry (cf. Equation 20), a first-order reconstruction may well prove insufficient and fail to produce the most accurate wall pressure solution. To overcome this problem, the second-order reconstruction of the wall pressure between a given wall point W and its corresponding second image points $B_2$ can be obtained by applying the following formula:

$$p_{w,o2} = p_{B_2} - \left. \frac{\partial p}{\partial \mathbf{n}} \right|_{B_2} \cdot \left| \overline{\overline{W B_2}} \right| + \left. \frac{\partial^2 p}{\partial \mathbf{n} \mathbf{n}} \right|_{B_2} \cdot \frac{\left| \overline{\overline{W B_2}} \right|^2}{2}, \quad (24)$$

where $\partial_{nn}p$ is the second directional derivative of the pressure in the wall normal direction:

$$\begin{align}
\partial_{nn} p = \frac{\partial^2 p}{\partial \mathbf{n} \mathbf{n}} &= \mathbf{n}^T \cdot \nabla^2 p \cdot \mathbf{n} \nonumber \\
&= n_x (n_x \partial_{xx} p + n_y \partial_{xy} p + n_z \partial_{xz} p) \nonumber \\
&\quad + n_y (n_x \partial_{yx} p + n_y \partial_{yy} p + n_z \partial_{yz} p) \nonumber \\
&\quad + n_z (n_x \partial_{zx} p + n_y \partial_{zy} p + n_z \partial_{zz} p). \tag{25}
\end{align}$$

The use of a Cartesian mesh eliminates the need to consider any particular metric when calculating the second derivatives $\partial_{ij}p$ involved in Equation 24. Instead, a new gradient calculation is simply performed, this time on the first derivatives $\partial_i p$. This quadratic reconstruction also takes into account a larger stencil of interior computed points for interpolating the second pressure derivatives at the new set of image points. For some configurations, such as the one depicted in Figure 15-(c), this stencil would also include some of the first computed points positioned near the modeled area, which could be a problem. However, in the context of this study, the curvature effects are actually highly localized at the leading edge of the geometries and appear to be negligible elsewhere. As explained in the previous section, it is also these regions that are most susceptible to modeling height adaptation. Therefore, it seems possible to safely apply this higher-order wall pressure reconstruction without using any of the pressure information obtained at the first calculated points, while maintaining the current location of the second image points.
<figure><img src="image_17.png" /><figcaption>Figure 15: Stencils of computed points whose pressure solution is taken into account for updating the pressure and pressure gradient information at the computed donor points, for increasing formulation orders.</figcaption></figure>

## 6. Numerical results

This section reviews all validation cases used to demonstrate that the modified IBM pre- and post-processing provides a better capture of skin pressure and friction distributions, as well as a better prediction of drag and lift coefficients with a near-field integration of all aerodynamic forces on an immersed obstacle. For this purpose, steady-state RANS simulations of subsonic turbulent flows around two airfoils and an aircraft were performed on Cartesian grids.

IBM results are systematically compared with body-fitted approaches, where the turbulent boundary layer is resolved directly by the mesh. Wherever possible, the same numerical parameters are used, including the turbulent model, the time integration method, and the numerical scheme used to estimate the convective fluxes. For the sake of clarity, a common standard is also adopted for the following figures and tables:

<figure><img src="image_18.png" /></figure>

In the following, each $y_{target}^{+}$ is given as a function of the near-wall resolution of the different Cartesian meshes considered, so that the mesh spacing measured at the first calculated points always appears perfectly in line with the required spacing found on a conformal mesh, at the same
distance from the wall. These values are chosen between 100 and 500, which ensures that the wall model is theoretically applied within its range of validity, given the Reynolds numbers that characterize the turbulent flows studied in this section.

In contrast to previous studies, more emphasis is now placed on capturing the wall pressure at the leading edge. The convergence rate of the IBM approach is thoroughly investigated, using the three reconstructions proposed in the last section, as well as increasing mesh resolutions in the vicinity of each wall. The behavior of the pressure solution throughout the entire boundary layer is also analyzed in great detail, using different profiles extracted from the upper surface of each bidimensional test-case. Finally, in order to study the most critical cases and to avoid any redundancy with the results already discussed in [1], only large angles of attack are considered.

## 6.1. Simulation of subsonic turbulent flows around a NACA0012 profile

The first test-case is the simulation of a subsonic turbulent flow around a NACA0012 airfoil with a critical flow incidence of fifteen degrees. A freestream Mach number of 0.15 is considered, as well as a Reynolds number of 6 million, based on a chord length of $c = 1$. This bidimensional test-case has already been studied in [1], under the same flow conditions and with the same simulation parameters, summarized in Table 1 below:

<table><caption>Table 1: NACA0012 airfoil (Ma 0.15 - Re 6M - AoA 15°): simulation parameters.</caption><thead><tr><td>Turbulence Model</td><td>Wall Model</td><td>Numerical Scheme</td><td>Time Integration</td><td>CFL</td></tr></thead><tbody><tr><td>RANS-SA</td><td>Musker</td><td>AUSM+(P)</td><td>BDF1</td><td>5</td></tr></tbody></table>

In 2021, good agreement was found between the improved IBM method and the reference body-fitted wall-resolved solution, whether on Cartesian or body-fitted meshes with Musker's wall function applied at the same threshold distance. Better results were obtained when the wall model was strictly applied in the logarithmic sub-layer, by imposing $y_{target}^{+} \in [100, 500]$, although this could lead to significant mesh refinement around the geometry.

This test case will now be used to investigate the impact of the immersed boundary method on the evaluation of the aerodynamic forces, with a particular emphasis on the drag coefficient. Since we have already discussed the negligible effect of the Cartesian mesh itself on the resulting skin pressure and friction distributions when imposing a proper near-wall spacing according to the positioning of the target and image points, the following parametric study will only take place on a
fine Cartesian mesh, using an ideal $y_{target}^{+} = 100$. The specific contribution of the modeling height adaptation on the one hand and the reconstruction order of the wall pressure on the other hand will be analyzed both separately and crosswise.

Information on the mesh properties can be found in Table 2, while views of the body-fitted mesh and the Cartesian mesh are presented in Figure 16. The IBM6 Cartesian mesh was automatically generated with Cassiopee, based on a near-wall spacing of $h_{min} = 1e-4$ and a domain extent of $d_{ff} = 25c$. Special care was taken to extend the width of the first Cartesian grid around the geometry to position both layers of image points and their interpolation donors within the finest mesh resolution. To this end, the minimum number of cells per direction composing each elementary Cartesian grid was set to $v_{min} = 25$. A local refinement box was also placed at the trailing edge to successfully propagate the minimum near-wall resolution around it.

<table><caption>Table 2: NACA0012 airfoil (Ma 0.15 - Re 6M - AoA 15°): mesh characteristics.</caption><thead><tr><th></th><th>IBM6</th><th>BF</th></tr></thead><tbody><tr><td>Minimum spacing $h_{min}/c$</td><td>1e-4</td><td>1e-6</td></tr><tr><td>Corresponding $y_{target}^{+}$</td><td>100</td><td>1</td></tr><tr><td>Number of cells</td><td>3.8M</td><td>8M</td></tr><tr><td>Proportion of cell types (blanked / computed / interpolated)</td><td>23% / 65% / 12%</td><td>-</td></tr></tbody></table>

<figure><img src="image_19.png" /><figcaption>Figure 16: Mesh views around the NACA0012 airfoil.</figcaption></figure>

The final location of the IBM points around the profile is examined in Figure 17, before and after the modeling height adaptation. As expected, the greatest effect of the approach proposed in section 4 is seen at the leading edge, where the first computed points now appear much closer to the
595

wall. Further downstream, the distribution is virtually unaffected by this adaptive process. This figure also confirms that the second image points are not actually shifted towards the geometry like the rest of the IBM points, to ensure an accurate capture of the first- and second-order pressure gradients, sufficiently far away from the modeled area.

600

The present case seems ideal to study the contribution of a local adaptation of the modeling height after a first simulation with a global threshold around the obstacle. While such a treatment seems mostly unnecessary in the absence of flow incidence, as demonstrated in [1], an adaptation of the modeling height greatly reduces the over-estimation of the $y^+$ values obtained at the first calculated points for a turbulent flow facing an incidence of fifteen degrees, thus limiting the error due to the wall model itself when the latter is applied outside its domain of validity.

<figure><img src="image_20.png" /><figcaption>Figure 17: NACA0012 airfoil (Ma 0.15 - Re 6M - AoA 15°): IBM points location for the IBM6 Cartesian mesh, before and after the modeling height adaptation.</figcaption></figure>

605

Figure 18 illustrates this aspect by showing the $y^+$ distribution around the airfoil before (F42: —) and after (F42-A: —) the adaptation. The strongest values on the suction side are almost divided by two and clearly justify this additional operation. As mentioned above, the modeling height distribution actually converges in a single step, which eliminates the need to dynamically modify the location of the IBM points during the simulation and reduces additional cost. Unfortunately, the maximum values monitored at the leading edge that exceed the initial $y^+_{target} = 100$ cannot be further reduced due to the safety previously imposed on this operation, which prevents any computed point from lying too low in the turbulent boundary layer. Only a finer Cartesian mesh would allow us to bring the first computed points so close to the geometry.

610
<figure><img src="image_21.png" /><figcaption>Figure 18: NACA0012 airfoil (Ma 0.15 - Re 6M - AoA 15°): y⁺ distributions at the first image points.</figcaption></figure>

Figures 19 to 21 examine the flow solution behavior across the turbulent boundary layer, for different sections of the profile, taken at 0.1%, 1%, 10%, 50% and 90% of the chord length. The focus here is solely on the upper side of the geometry, where the solution is most likely to deviate from the reference body-fitted solution. It should also be noted that each profile is plotted from the position of the first calculated point.

For the first 10% of the chord length, it is clear from the dimensionless velocity profiles (Figure 19) that the equilibrium wall function is unable to adequately represent the flow behavior, as it quickly leaves the inner region of the turbulent boundary layer. Conversely, beyond the first 10% of the chord, the boundary layer is now sufficiently developed, with a well-defined logarithmic region, for which the wall model is once again perfectly suited.

In this case, the use of a local modeling height seems essential to reduce the risk of over-estimating the boundary layer thickness at the leading edge. Otherwise, the non-adapted solutions are more prone to deviate from the reference solution and even exhibit a different flow behavior, as observed on the turbulent viscosity profiles near the trailing edge (Figure 21), where the error made at the leading edge (see the zooms in Sections 1 and 2) is carried further downstream and eventually added up. In addition, the error made at the wall boundary also appears to affect the flow prediction further away from the geometry, as shown by the velocity profiles (Figure 20), making the adaptation operation even more critical for this test-case.
<figure><img src="image_22.png" /><figcaption>Figure 19: NACA0012 airfoil (Ma 0.15 - Re 6M - AoA 15°): velocity profiles (logarithmic scale).</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>Figure 20: NACA0012 airfoil (Ma 0.15 - Re 6M - AoA 15°): velocity profiles (linear scale).</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>Figure 21: NACA0012 airfoil (Ma 0.15 - Re 6M - AoA 15°): turbulent viscosity profiles.</figcaption></figure>
Figure 22 gathers the pressure gradient distributions obtained in the wall normal and tangent directions. These results highlight magnitudes that should not be overlooked when modeling near-wall turbulence, especially over the first 10 % of the chord length when the boundary layer is the thinnest. For the pressure gradients measured in the wall normal direction, a comparison is then made between the values obtained at the first image points (___) and the same values obtained at the second image points, either positioned at iso-distance from the wall (---) or following a stepped front (----).

As expected, the gradients recorded at the first image points show significant oscillations, particularly at the leading edge where the values are actually the strongest. These distributions are obviously too noisy to be used for higher-order wall pressure reconstructions. In contrast, the pressure gradients obtained at the second image points confirm the benefits of their introduction into our methodology. Without any spatial smoothing, the distributions no longer exhibit spurious oscillations and are therefore perfectly suitable for first- and second-order extrapolations of the pressure solution. Furthermore, it appears that the positioning of the second set of image points around the geometry does not have a significant impact on the final solution. Although the irregular stepped front is still partially visible on the pressure gradient distribution, the zoom required to perceive it does not justify the application of a more regular front, which is more difficult to implement. For this reason, only the positioning of the second image points described in Figure 14-(b) will be considered in the following.

In Figure 22-(a), the tangential pressure gradient distributions around the airfoil also show that a strong adverse pressure gradient ($\partial_t p > 0$) develops on the suction side of the wing at large angles of attack. Similarly, favorable pressure gradients ($\partial_t p < 0$) occur near the leading edge. Since these pressure gradients are located in areas already identified as sources of discrepancies between the reference and the IBM solutions, this certainly implies that our wall model is missing some non-negligible information in these zones when reconstructing the velocity field in the turbulent boundary layer. Therefore, a wall model that takes into account the pressure gradients, both favorable and adverse, would certainly improve the solution in these regions.
<figure><img src="image_25.png" /><figcaption>Figure 22: NACA0012 airfoil (Ma 0.15 - Re 6M - AoA 15°): pressure gradient distributions.</figcaption></figure>

Looking at the surface solution, Figure 23-(b) shows an excellent agreement between the IBM solutions and the reference for the skin friction distribution. While the non-adapted solution already successfully approaches the body-fitted solution, as already demonstrated in [1], the contribution of the modeling height adaptation is clearly visible at the leading edge. Thanks to this additional operation, the friction peak is now much better represented in the first 5% of the chord length, a region where the wall model is globally deficient. To a lesser extent, the friction solution is also better predicted over the entire profile, especially towards the trailing edge, as shown in the second enlargement. Although it seems very difficult to further improve the friction prediction at the leading edge on the current Cartesian mesh, without making major modifications to the wall model, the pressure distributions shown in Figure 23-(a) demonstrate that it is nevertheless possible to fully exploit all the information derived from the converged computation. In an effort to improve the prediction of this quantity at the wall, three types of wall pressure reconstructions are then considered, of order 0 (—), 1 (---), and 2 (···), respectively.

The direct extrapolation of the pressure solution through the turbulent boundary layer already shows a relatively noiseless pressure peak at the leading edge for the adapted IBM solution (in blue), which also appears quite close to the reference solution (in black). On the contrary, this type of reconstruction clearly leads to a significant under-estimation for the non-adapted IBM solution (in purple), which can be easily explained by the variation in the location of the first image points before and after the modeling height adaptation (cf. Figure 17). To overcome this problem, the linear reconstruction of the wall pressure greatly improves the prediction of the pressure peak on the suction side, both in terms of position and intensity. More importantly, far from destabilizing
680

the solution, the use of pressure and pressure gradient information further down the mesh actually smooths the pressure coefficient distribution, eliminating the small amount of numerical noise observed with the previous reconstruction. For this fine Cartesian mesh, order convergence is achieved fairly quickly, as can be seen by superimposing the pressure solutions obtained by the first- and second-order reconstructions. Thus, the use of local curvature information seems unnecessary for the current mesh resolution and the chosen $y_{target}^{+}$ parameter.

685

<figure><img src="image_26.png" /><figcaption>Figure 23: NACA0012 airfoil (Ma 0.15 - Re 6M - AoA 15°): skin pressure and friction distributions.</figcaption></figure>

690

To better understand the improvements made to the skin pressure distributions, Figure 24 shows the behavior of the pressure solution at different sections of the suction side. This figure proves that the zero-pressure-gradient assumption does not quite hold near the leading edge at $x/c = 0.001$, where the pressure clearly varies in the wall normal direction, close to the geometry. As stated earlier, although the pressure profiles obtained with the wall-modeled solutions agree well with the reference solution above the modeling height and even almost perfectly with the IBM F42-A solution, the wall pressure is actually degraded during the post-processing with a standard extrapolation of this solution from the image points. The close-up view near the leading edge illustrates this by showing how the final wall pressures, indicated by square-shaped markers (■), are ultimately positioned further away from the baseline reference as the modeling height increases, from one IBM solution to another. Moving along the chord, the original assumption seems to hold again, even in the case of strong flow incidences, which is consistent with the previous observations

695
made on the velocity profiles. For the next 90% of the chord length, the pressure gradient in the normal direction is actually close to zero, as shown in Figure 22-(a), allowing us not to degrade the first prediction of the wall pressure solution, but to improve it wherever possible.

This problem is addressed by using Equations 22 and 24 to extract the wall pressure using both linear (▲) and quadratic (●) reconstructions. Near the image point location, the pressure solution tends to follow a linear evolution in the wall normal direction, which explains the rapid order convergence observed on the skin pressure distribution for the IBM6 Cartesian mesh. The corrected pressures now cluster around the expected values, and even more so after the modeling height adaptation. At a greater distance from the wall, the profile as a whole begins to exhibit significant curvature effects that cannot be accounted for without at least a second-order formulation. Therefore, the need to consider an additional reconstruction order for the wall pressure would become stronger with coarser Cartesian meshes coupled with larger $y_{target}^{+}$ values, as the image points will end up further away from the geometry. This aspect will be discussed in detail in the next applications.

<figure><img src="image_27.png" /><figcaption>Figure 24: NACA0012 airfoil (Ma 0.15 - Re 6M - AoA 15°): pressure profiles.</figcaption></figure>

The integration of the pressure and friction forces exerted on the body is finally studied in Figure 25, which collects all the drag and lift coefficients obtained with the IBM solutions using the near-field and far-field approaches described in section 3, with respect to the baseline body-fitted wall-resolved solution. In the absence of any wall pressure correction during the IBM post-processing, the drag coefficients are greatly overestimated with a near-field integration, especially
when the flow faces strong flow incidences. Further investigation shows that the pressure drag is mainly responsible for these deviations, while frictional drag remains close to the baseline solution. It appears that the discrepancies observed on both the skin pressure distributions (Figure 23-(a)) and the pressure profiles (Figure 24) eventually lead to even larger errors after integrating the pressure solution on the immersed geometry. In fact, since the suction peak is highly localized near the leading edge at higher angles of attack, most of the pressure drag is also generated in this zone and even small differences in its location ultimately lead to very large errors in the drag coefficient. As before, the induced error also increases linearly with the distance at which the wall model is applied. The relative error is then close to 100% for a flow incidence of fifteen degrees before the modeling height adaptation, but only half as much after.

As expected, the encouraging results observed on the pressure profiles due to the post-processing corrections systematically lead to a better capture of the drag coefficients with the near-field approach, as shown in Figures 25-(c) and 25-(d). There is also little to no difference between the near-field and far-field approaches when the wall pressure is correctly estimated using first- or second-order reconstructions, depending on the Cartesian mesh resolution or flow incidence. In the end, both approaches provide quasi-identical results, in good agreement with the reference solution, and the final drag coefficient is no longer dependent on the integration mean, as it should be.

Figure 25 also confirms the previous analysis regarding the mutual interest of better IBM pre- and post-processing. While the latter unlocks the full potential of a given IBM simulation in terms of surface solution reconstruction, the former is crucial for reducing the error due to the wall model itself and improving the flow solution as a whole. The combined use of the two developments presented in this paper appears to be the key to observing the best drag and lift coefficients within our current IBM methodology. This is highlighted in Table 3, which provides a detailed view of the drag and lift coefficients at a fifteen-degree flow incidence, with and without modeling adaptation, with and without higher-order pressure reconstructions. The worst and best case scenarios are circled in red and green, respectively. The aerodynamic coefficients are consistently better approximated after the adaptation process, regardless of the wall pressure reconstruction, which is perfectly in line with previous observations. For the adapted IBM solution, the relative error of the drag coefficients even drops from 47% to 17% with respect to the reference, as the extrapolation order increases.
<figure><img src="image_28.png" /><figcaption>(a) Lift - Before the modeling height adaptation</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>(b) Lift - After the modeling height adaptation</figcaption></figure>

<figure><img src="image_30.png" /><figcaption>(c) Drag - Before the modeling height adaptation</figcaption></figure>

<figure><img src="image_31.png" /><figcaption>(d) Drag - After the modeling height adaptation</figcaption></figure>

<caption>Figure 25: NACA0012 airfoil (Ma 0.15 - Re 6M): lift and drag coefficients for increasing angles of attack.</caption>

<table><caption>Table 3: NACA0012 airfoil (Ma 0.15 - Re 6M - AoA 15°): detailed drag and lift coefficients.</caption><thead><tr><th rowspan="2"></th><th colspan="3">IBM6 F42</th><th colspan="3">IBM6 F42-A</th><th rowspan="2">BF</th></tr><tr><th>0th-order</th><th>1st-order</th><th>2nd-order</th><th>0th-order</th><th>1st-order</th><th>2nd-order</th></tr></thead><tbody><tr><th>C<sub>D,v</sub></th><td>5.02e-3</td><td>5.02e-3</td><td>5.02e-3</td><td>5.46e-3</td><td>5.46e-3</td><td>5.46e-3</td><td>5.29e-3</td></tr><tr><th>C<sub>D,p</sub></th><td>35.42e-3</td><td>26.10e-3</td><td>26.89e-3</td><td>25.11e-3</td><td>19.31e-3</td><td>19.05e-3</td><td>15.56e-3</td></tr><tr><th>C<sub>D</sub></th><td>40.44e-3</td><td>31.12e-3</td><td>31.91e-3</td><td>30.57e-3</td><td>24.77e-3</td><td>24.51e-3</td><td>20.85e-3</td></tr><tr><th>C<sub>L</sub></th><td>1.416</td><td>1.421</td><td>1.421</td><td>1.471</td><td>1.476</td><td>1.476</td><td>1.516</td></tr></tbody></table>

err 94%/5.6%

err 47%/2.9%

err 17%/2.6%
Despite these significant improvements, some differences in drag coefficient prediction remain. This is particularly evident at fifteen degrees, near the stall angle, where a relative error of almost 20% is observed for the optimal IBM solution. It should be noted that this problem has already been raised in the literature, notably by Berger and Aftosmis [41]. In their study, the authors observed deviations comparable to ours with respect to this integral quantity, regardless of the excellent behavior of their solutions, using a modified cut-cell approach coupled with advanced wall modeling techniques. In order to explain these apparently contradictory results, they were able to cite the work of Diskin [52], which highlighted the sensitivity of the present test-case. For our part, it would then be necessary to apply the wall model nowhere but in the viscous sublayer, by imposing very low values of $y_{target}^{+}$ at the leading edge, which would go against the philosophy behind the use of wall laws.

To conclude this first study, Figures 25-(a) and 25-(b) reveal that the lift coefficient does not seem to vary with the integration method or the way the wall pressure is extracted, unlike the drag coefficient. As a matter of fact, in both cases, the wall-modeled solutions agree well with the body-fitted reference for this integral quantity. Once again, the prediction at fifteen degrees is the most critical, for which the non-adapted solution carries a relative error of almost 10%. The adapted solution is the only one able to keep the relative error below 3% for the lift coefficient at each flow incidence considered.

## 6.2. Simulation of subsonic turbulent flows around a 2D multi-element airfoil

The second test-case is the simulation of a subsonic turbulent flow around the 2D multi-element airfoil (2DMEA) facing a strong sixteen-degree flow incidence. This bidimensional case, extensively documented on the NASA TMR website [53], has been used to analyze the impact of the mesh in an unstructured context, mainly for compressible RANS approaches based on the SA [54] and SA-RC [55] turbulence models. This profile geometry is more representative of a real industrial configuration than the NACA0012 with three separate elements, and provides the perfect opportunity to further validate our recent developments with more challenging flow behavior before tackling a tridimensional test-case at the end of this section. The 2DMA profile can be seen as a cross-section of the High-Lift Common Research Model (HL-CRM [56]), taken about 16 meters from the wing root. A Mach number of 0.2 is considered, with a Reynolds number of 5 million, based on the chord length of the retracted airfoil. The main simulation parameters are summarized in Table 4:
<table><caption>Table 4: 2DMEA (Ma 0.2 - Re 5M - AoA 16°): simulation parameters.</caption><thead><tr><td>Turbulence Model</td><td>Wall Model</td><td>Numerical Scheme</td><td>Time Integration</td><td>CFL</td></tr></thead><tbody><tr><td>RANS-SA</td><td>Musker</td><td>Roe</td><td>BDF1</td><td>5</td></tr></tbody></table>

This time, the reference solution is provided by FUN3D, a Finite-Volume unstructured solver developed at NASA [57]. In the following, only the finest body-fitted mesh from the first family of hex-dominant unstructured meshes proposed on the NASA TMR website (with first-order turbulent advection only) will be considered for comparison with the best available reference solution. For this test-case, the previous AUSM+(P) scheme is then replaced by the Roe scheme introduced in section 2 to get as close as possible to the numerical parameters of the reference, which is also based on a flux-difference splitting scheme.

In contrast to the validations conducted on the NACA0012 profile, this subsection examines four quadtree Cartesian meshes whose information is collected in Table 5. Each of them is separated by one level of refinement in order to study the evolution of the IBM solutions as a function of the near-wall resolution. The idea is to evaluate the immersed boundary method under optimal conditions, with a local modeling height and for increasing $y_{target}^{+}$ values, which are perfectly suited for each Cartesian mesh. The previous parametric study is thus replaced by a mesh convergence test on a more complex geometry. However, it is important to note that of these four Cartesian meshes, only the IBM4 resolution can be realistically applied to an industrial-level test-case. While the other meshes are included to illustrate what we could theoretically achieve with unlimited computing resources, reducing the error introduced by our immersed boundary method on the coarsest mesh is actually the primary objective of this second test-case.

To better understand the geometry at hand, several views of the L7 body-fitted and IBM6 Cartesian meshes are shown in Figure 26. Both have approximately the same number of cells. For this particular airfoil, the near-wall resolution is severely limited by the small spacing between the three elements that make up the 2DMEA. In fact, a large near-wall spacing would, at some point, prevent the correct positioning of the IBM points around the geometry, a fortiori when considering the IBM F42 approach and its additional rows of target points. For this specific reason, but also due to the inherent limitation of the wall model currently implemented in our framework, coarser Cartesian meshes than IBM4 were not considered.
<table><caption>Table 5: 2DMEA (Ma 0.2 - Re 5M - AoA 16°): mesh characteristics.</caption><thead><tr><th></th><th>IBM4</th><th>IBM5</th><th>IBM6</th><th>IBM7</th><th>BF</th></tr></thead><tbody><tr><td>Minimum spacing h<sub>min</sub>/c</td><td>4e-4</td><td>2e-4</td><td>1e-4</td><td>5e-5</td><td>1e-6</td></tr><tr><td>Corresponding y<sup>+</sup><sub>target</sub></td><td>400</td><td>200</td><td>100</td><td>100</td><td>1</td></tr><tr><td>Number of cells</td><td>1.4M</td><td>2.7M</td><td>5.2M</td><td>10.4M</td><td>5.9M</td></tr></tbody></table>

<figure><img src="image_32.png" /><figcaption>Figure 26: Mesh views around the 2DMEA.</figcaption></figure>

Figure 27 below shows the location of all IBM points (wall points, target points, first and second image points) around the multi-element airfoil for the IBM6 Cartesian mesh. In the same figure, where the adaptation process has been successfully completed, it can be seen in each local zoom that the positioning of the IBM points has been particularly affected at the leading edge of each element and at the top surface of the geometry, while the bottom surface and the trailing edges remain unchanged after this operation. At first glance, it would seem that the methodology introduced in section 4 is naturally applicable to more complex geometries.

<figure><img src="image_33.png" /><figcaption>Figure 27: 2DMEA (Ma 0.2 - Re 5M - AoA 16°): IBM points location for the IBM6 Cartesian mesh, after the modeling height adaptation.</figcaption></figure>
Figure 28 examines the mesh convergence speed of each Cartesian mesh, from coarsest to finest, before and after the modeling height adaptation described in Section 3. The values of the lift coefficients computed on-the-fly during the CFD simulations are plotted as a function of the number of iterations. To extend this analysis, the convergence information of the initial RANS simulation is plotted ahead, using light solid lines. This convention (F42: — & F42-A: —) is followed throughout the subsection. For this first convergence with a global modeling height, it can be seen that the IBM7 mesh requires four times as many iterations as the IBM4 mesh, which could be a problem in a more complex configuration. In fact, the convergence rate for this initialization is quasi-proportional to the near-wall resolution. After the adaptation process, which is explicitly localized on each plot, this figure shows how important a warm start is for the method to reduce its inherent cost. Although the lift coefficients are significantly improved afterwards, a cold start would simply double the whole process and partially compromise the competitiveness of the approach we defend. On the other hand, restarting from the previous steady-state solution requires only an additional number of iterations of the order of a third of the initial simulation. In the end, both the adapted solutions obtained through warm and cold starts converge towards the same improved CL values, which also demonstrates the stability of our approach for this test-case.

<figure><img src="image_34.png" /><figcaption>Figure 28: 2DMEA (Ma 0.2 - Re 5M - AoA 16°): mesh convergence information.</figcaption></figure>
The importance of the adaptation process is even more apparent in this application, as the
differences are now directly visible in the flow field, even at a significant distance from the airfoil.
To illustrate this, Figure 29 describes the evolution of the eddy viscosity contours around the
2DMEA, for which the full set of IBM solutions provides very good qualitative results compared
to the reference body-fitted result after the modeling height adaptation, as can be seen on the
right side of this figure. On the other hand, the unadapted IBM solutions clearly overestimate the
thickness of the boundary layers created by the slat and the main element. When the turbulent
wakes interact behind the airfoil, near the slat trailing edge, it is actually possible to measure
the discrepancies in terms of wake size with the information taken from the reference. Before the
adaptation, the error decreases proportionally to the modeling height, from almost 80% with the
coarsest Cartesian mesh to a little more than 20% with the finest. In other words, in this case, even
the best IBM solution does not come close to the fitted solution. After the adaptation, the error is
now fully contained, even for the coarsest IBM solution, which produces an error close to the finest
IBM solution taken from the left side of the figure. This can also be observed inside the slat cavity,
where the turbulent intensity at the center of the fluid recirculation seems to be over-estimated for
all IBM solutions using a global modeling height. Finally, this analysis will be later confirmed using
unidirectional probes in each region of interest of the fluid domain, i.e. on the suction side of each
element.

For a more quantitative analysis, Figure 30 presents a detailed view of the integrated aero-
dynamic forces acting on the 2DMEA for increasing orders of wall pressure reconstruction. The
maximum relative errors on the lift and drag coefficients (expressed as percentages) are given for
each dataset and compared to the FUN3D solution using the coarsest IBM solution, before and
after the modeling adaptation. As expected, the following results are consistent with those obtained
for the NACA0012 profile. In particular, it appears that the pressure contributions dominate in
the estimation of the total drag, since the viscous contributions are much more negligible and, in
fact, very close to the reference solution. In light of these observations, the error obtained for the
drag coefficient is once again induced and driven by the error made on the pressure solution. More
importantly, these results also highlight the cross-benefits of the two improvements brought to our
methodology, since both the modeling height adaptation and the higher-order extractions of the
wall pressure are ultimately necessary to reduce the error made on the aerodynamic coefficients to
acceptable levels, regardless of the mesh resolution.
<figure><img src="image_35.png" /><figcaption>Figure 29: 2DMEA (Ma 0.2 - Re 5M - AoA 16°): eddy viscosity contours.</figcaption></figure>
Unsurprisingly, the solutions obtained for $y_{target}^{+} = 100$ (IBM6 and IBM7 Cartesian meshes) are the closest to the body-fitted reference solution, despite a relative drag error of the order of 20% at best, similar to that observed on the NACA0012 profile at a fifteen-degree flow incidence and for comparable flow parameters. This error is mainly due to the high sensitivity of the solution at the leading edge, which is not necessarily representative of the difficulties encountered in three dimensions for this type of high-lift configuration.

In this case, it is also interesting to study the effect of the wall pressure reconstruction for larger modeling heights. For the finest meshes, a linear reconstruction of the pressure solution from the second image points is sufficient to satisfy the convergence conditions for the drag coefficient. Conversely, for coarser meshes (IBM4 and IBM5 Cartesian meshes), a quadratic reconstruction of the wall pressure may actually be useful. Indeed, for $y_{target}^{+} = 200$ and $y_{target}^{+} = 400$, taking into account local curvature effects further away from the immersed geometry can still halve the relative error. This second-order wall pressure reconstruction even allows the adapted solution obtained at $y_{target}^{+} = 200$ to practically converge to the adapted solutions obtained at $y_{target}^{+} = 100$.

Again, the lift coefficient (Figure 30-(a)) seems to be less sensitive to the order of the wall pressure reconstruction, as the relative error hardly varies between the datasets related to each Cartesian mesh. For this coefficient, the use of a local modeling height seems to have a much larger impact, inducing the same shift towards the reference for all IBM solutions. Quite satisfactorily, relative errors of less than 5% are observed on the lift coefficient for the IBM5, IBM6, and IBM7 meshes. For the IBM4 mesh, however, it seems difficult to imagine a larger reduction of the current minimum error. Upon further investigation, this appears to be due to the problems associated with the leading edges of profile geometries, i.e. the inability of our wall model to perfectly reproduce the behavior of a turbulent boundary layer without a sufficiently established logarithmic region.

Finally, the following Table 6 provides a more detailed view of the drag and lift coefficients, only for a second-order reconstruction of the wall pressure and after the adaptation of the modeling height. The convergence rate of our approach can be calculated using the four Cartesian meshes considered in this section, with respect to the reference body-fitted solution. For the drag coefficients, the relative error decreases rapidly from 42% (IBM4) to 27% (IBM5) and ends at 23% (IBM6) and 22% (IBM7). For the lift coefficients, the convergence is even faster, with the error being divided by two for each new Cartesian mesh.
<figure><img src="image_36.png" /><figcaption>Figure 30: 2DMEA (Ma 0.2 - Re 5M - AoA 16°): lift and drag coefficients.</figcaption></figure>

<table><caption>Table 6: 2DMEA (Ma 0.2 - Re 5M - AoA 16°): detailed drag and lift coefficients.</caption><thead><tr><th></th><th>IBM4 F42-A</th><th>IBM5 F42-A</th><th>IBM6 F42-A</th><th>IBM7 F42-A</th><th>BF</th></tr><tr><th></th><th>2nd-order</th><th>2nd-order</th><th>2nd-order</th><th>2nd-order</th><th>-</th></tr></thead><tbody><tr><th>C<sub>D,v</sub></th><td>0.92e-2</td><td>0.96e-2</td><td>0.96e-2</td><td>0.94e-2</td><td>0.93e-2</td></tr><tr><th>C<sub>D,p</sub></th><td>7.93e-2</td><td>6.96e-2</td><td>6.70e-2</td><td>6.63e-2</td><td>5.29e-2</td></tr><tr><th>C<sub>D</sub></th><td>8.85e-2</td><td>7.92e-2</td><td>7.66e-2</td><td>7.57e-2</td><td>6.22e-2</td></tr><tr><th>C<sub>L</sub></th><td>3.472</td><td>3.620</td><td>3.721</td><td>3.734</td><td>3.779</td></tr></tbody></table>

To complete this analysis, a detailed study of the pressure and friction coefficient distributions is presented in Figures 31 and 32 to ensure that the aerodynamic coefficient capture does not benefit from compensatory effects. These results also confirm the extension of our new IBM methodology to multi-element airfoils, given the absence of spurious oscillations on all IBM solutions, as well as the good (IBM4, IBM5) or even excellent (IBM6, IBM7) agreement observed with the reference. As in the previous test-case, local enlargements of the skin pressure distribution indicate the regions that are particularly affected by the wall pressure correction, which happen to be the leading edges of each element. These distributions finally shed some relevant light on the integral results discussed above, especially for the solution obtained on the IBM4 mesh, where the curvature effects are much more clearly taken into account with the second-order reconstruction. For the latter, however, a slight under-estimation of the total pressure distribution on the upper surface of the wing is observed, which could explain the discrepancies in terms of drag coefficients, despite an optimal reconstruction of the solution at the wall.
<figure><img src="image_37.png" /><figcaption>Figure 31: 2DMEA (Ma 0.2 - Re 5M - AoA 16°): skin pressure distributions.</figcaption></figure>

Figure 32 shows a good capture of the skin friction information for all IBM solutions, especially for the fine IBM6 and IBM7 meshes, which even manage to capture the correct friction maxima near the wing and flap leading edges, failing only to perfectly represent the flow around the slat where the boundary layer is the thinnest. The IBM5 mesh shows similar behavior to the more resolved results, despite a slight under-estimation of the friction maxima, as well as a more damped behavior on the suction side of the 2DMEA. The same trend is observed for the IBM4 mesh, with some additional difficulties near the trailing edges, mainly due to the initial over-estimation of the boundary layer thickness. However, the accuracy obtained on the coarsest Cartesian mesh is still largely positive and validates the current IBM methodology for higher $y^+$ values without losing too much information, which would prove useful for 3D applications.
<figure><img src="image_38.png" /><figcaption>Figure 32: 2DMEA (Ma 0.2 - Re 5M - AoA 16°): skin friction distributions.</figcaption></figure>

To conclude this study of the 2D multi-element airfoil, Figures 33 and 34 describe the evolution of the flow solution within the turbulent boundary layer, for different sections close to the airfoil, all chosen at constant abscissa. Satisfactorily, the IBM solutions obtained on the finest meshes correspond almost perfectly to the reference body-fitted solution for both velocity and turbulent viscosity solutions. Only small deviations are observed in the last two sections, where the slat wake interacts with the near-wall flow, either on the wing or on the flap top surface. For the coarsest meshes, these small deviations eventually lead to significant gaps with the reference solution as we move away from the wall, but the solutions still remain close to FUN3D. On the contrary, the IBM solutions obtained before the adaptation process deviate much more from the baseline. In fact, the effect of this adaptation is so important that the adapted IBM4 solution almost matches the unadapted IBM6 solution, which is consistent with the observations made earlier.

Naturally, the solution taken at $y_{target}^{+} = 400$ suffers the most from this initial over-estimation of the turbulent boundary layer around the slat. This is particularly evident for the turbulent viscosity solution, where a sum of upstream errors eventually causes this IBM solution to deviate from the reference, well beyond the modeled zone near the wall. In the end, this simply raises the question of modeling a thin boundary layer without further refinement of the Cartesian mesh.
<figure><img src="image_39.png" /><figcaption>Figure 33: 2DMEA (Ma 0.2 - Re 5M - AoA 16°): x-velocity (linear scale) profiles.</figcaption></figure>
<figure><img src="image_40.png" /><figcaption>Figure 34: 2DMEA (Ma 0.2 - Re 5M - AoA 16°): turbulent viscosity profiles.</figcaption></figure>
### 6.3. Simulation of subsonic turbulent flows around a NASA HL-CRM

The third and final application considered in this article is the simulation of subsonic turbulent flows around a full plane in high-lift configuration. The setup of this last application follows one of the cases presented at the fourth AIAA CFD High-Lift Prediction Workshop (HLPW [58]), which focused specifically on the High-Lift Common Research Model (HL-CRM) during the 2022 edition. In the following, only nominal inboard/outboard deflections are considered, with 30°/30° for the slats and 40°/37° for the flaps.

The various numerical simulations proposed in the fourth edition of this international workshop do represent a major challenge for all steady-state RANS methods, as recalled by Clark *et al.* [59] as well as Duensing *et al.* [60]. When approaching the maximum lift coefficient ($C_{Lmax}$), i.e. just before the stall phenomenon, it is generally accepted that no turbulence model is fully capable of representing the complex nature of the flow. In particular, the prediction of the aerodynamic coefficients is highly dependent on the fraction of the flow separated near the wall boundaries, which is difficult to predict with RANS approaches, and even more so with a wall model that attempts to approximate most of the boundary layer behavior. For this reason, this study focuses on the HL-CRM under a flow incidence of 17.05 degrees, which is still very challenging for immersed boundary methods coupled with wall modeling techniques such as ours, but avoids the additional difficulty associated with massively detached flow near the critical $C_{Lmax}$ value. Moreover, this choice of angle of attack is perfectly consistent with the previous test-cases.

The numerical setup, now operating in a fully dimensionalized context, considers a half wingspan of 29.4 meters, as recommended by the committee in charge of the workshop. The reference static temperature is set to 289 kelvin, the reference static pressure is set to 1.70093 bar, and the freestream velocity is imposed from a subsonic Mach number of 0.2. To achieve the required Reynolds number of 5.49 million, this time based on the mean aerodynamic chord (MAC = 7.005m), the reference viscosity is also artificially modified to a non-physical value. In the end, these two dimensionless numbers are therefore equivalent to those used to study the 2DMEA airfoil. The simulation parameters are summarized in the table below:

<table><caption>Table 7: NASA CRM-HL (Ma 0.2 - Re 5.49M - AoA 17.05°): simulation parameters.</caption><thead><tr><td>Turbulence Model</td><td>Wall Model</td><td>Numerical Scheme</td><td>Time Integration</td><td>CFL</td></tr></thead><tbody><tr><td>RANS-SA</td><td>Musker</td><td>Roe</td><td>BDF1</td><td>5</td></tr></tbody></table>
The body-fitted structured solutions provided by the LAVA CFD solver [61] were selected to compare the latest IBM solutions to a baseline reference. All LAVA results are currently freely available and compiled on the AIAA HLPW-4 seminar website, which greatly facilitated the following comparisons. The Launch Ascent and Vehicle Aerodynamics (LAVA) solver, developed at NASA Ames Research Center, is based on a cell-centered second-order Finite Volume approach that is sufficiently close to the numerical core of FastS. The choice was made to compare with the finest body-fitted solutions available, and therefore only the solution obtained on the finest *D* structured grid is considered in the following. It is also worth mentioning that this grid was made available by the workshop to all participants belonging to the study group dedicated to RANS approaches.

<figure><img src="image_41.png" /><figcaption>Figure 35: NASA CRM-HL (Ma 0.2 - Re 5.49M - AoA 17.05°): half-plane case description.</figcaption></figure>

Given the flow parameters of this application, a $y_{target}^{+}$ of 500 was considered as the best compromise between mesh size and accuracy of the flow solution, given the previous in-depth analysis performed on the 2DMEA. According to the F42 methodology, this initial choice imposes a fine near-wall resolution of 3.5mm. Unlike previous two-dimensional test-cases, this minimum resolution cannot be applied around the entire HL-CRM geometry without significantly increasing the final cell count. Therefore, no less than three resolution levels were taken into account during the mesh generation, as shown in Figure 35. The minimum near-wall spacing $h_{min}$ is only imposed in the vicinity of the slats, the flaps, and the nacelle pylon, where the modeling height is most likely to be adapted. The rest of the wing and the nacelle have twice the near-wall resolution, while the fuselage is deliberately left under-meshed. This data set resulted in an octree Cartesian mesh of about 2.4 $10^9$ points, including all blanked and interpolated cells. While this mesh is four times larger than the reference mesh (see Table 8), it is important to emphasize the speed at which it was generated.
Taking into account the specific treatment related to our immersed boundary method, as well as the expensive calculation of the turbulent distance flow field, the entire pre-processing stage was completed in less than seven hours in a massively parallel environment, from the geometry surface mesh alone. Conversely, it can take several weeks to generate a curvilinear grid as resolved as the *D* grid, and up to several months to calibrate it for unsteady WMLES applications [62].

<table><caption>Table 8: NASA CRM-HL (Ma 0.2 - Re 5.49M - AoA 17.05°): mesh characteristics.</caption><thead><tr><th>-</th><th>IBM4</th><th>BF</th></tr></thead><tbody><tr><td>Minimum spacing h<sub>min</sub></td><td>3.5e-3</td><td>N/A (y<sub>target</sub><sup>+,D</sup> = 0.75)</td></tr><tr><td>Corresponding y<sub>target</sub><sup>+</sup></td><td>500</td><td>1</td></tr><tr><td>Number of cells</td><td>2400M</td><td>550M</td></tr><tr><td>Cell natures (blanked / computed / interpolated)</td><td>21% / 54% / 25%</td><td>-</td></tr></tbody></table>

Figure 36 illustrates the number of iterations required to achieve convergence with a criterion based on the lift coefficient. Using the IBM approach proposed in this article, this simulation used 20x28 Intel Broadwell CPUs on ONERA's SATOR computer cluster for five runs of 15,000 iterations, including the modeling height adaptation that occurred after the first three runs.

<figure><img src="image_42.png" /><figcaption>Figure 36: NASA CRM-HL (Ma 0.2 - Re 5.49M - AoA 17.05°): mesh convergence information.</figcaption></figure>

Again, this figure highlights the clear benefit of using a local modeling height, as the lift coefficients shift towards the reference value. Table 9 in turn confirms the importance of using first- and second-order gradient information for the wall pressure reconstruction, as it divides the error made on the drag coefficient by six and reduces the error on the lift coefficient to 8%.
<figure><img src="image_43.png" /><figcaption>Table 9: NASA CRM-HL (Ma 0.2 - Re 5.49M - AoA 17.05°): detailed drag and lift coefficients.</figcaption></figure>

The following Figures 37 and 38 show three slices of the geometry, following the HLPW post-processing guidelines for authors. Using the quadratic extrapolation of the pressure solution at the surface, the skin pressure distribution obtained with the improved immersed boundary method actually converges towards the reference body-fitted for all three slices considered. In fact, these observations are very close to those already made on the previous 2DMEA airfoil using the coarsest Cartesian mesh. Both the slat and the main element of the wing are well predicted by our methodology, even when considering the pressure peaks at the leading edges. In the vicinity of the flaps, the solution behavior is slightly more damped, which may be caused by the inherently unsteady behavior of the flow in these regions as we approach the wing tip, where the flow begins to separate from the aircraft. The same analysis can be done for the skin friction coefficient. Although this flow quantity is more difficult to predict, the IBM solution manages to come close to the baseline solution. Furthermore, both skin solutions are completely free of any oscillations, confirming that the latest developments in our geometric approach do not destabilize its overall stability, but only improve its accuracy.

For a more qualitative analysis of the turbulent flow developing around the HL-CRM, different sections of the domain are displayed in Figure 39 for vorticity magnitude values taken at constant x-abscissa. Even in the absence of local refinements within the Cartesian grid further away from the plane, the turbulent flow structures are well predicted with the improved IBM approach, confirming the ability of the current methodology to accurately predict challenging flow behavior around complex geometries, even for strong flow incidence. However, approaching the stall phenomenon at higher angles would certainly require a finer resolution near the wall to allow the use of an even smaller modeling height when the flow begins to separate, which cannot be properly predicted by classical wall functions.
<figure><img src="image_44.png" /><figcaption>Figure 37: NASA CRM-HL (Ma 0.2 - Re 5.49M - AoA 17.05°): skin pressure distributions.</figcaption></figure>

<figure><img src="image_45.png" /><figcaption>Figure 38: NASA CRM-HL (Ma 0.2 - Re 5.49M - AoA 17.05°): skin friction distributions.</figcaption></figure>
<figure><img src="image_46.png" /></figure>

<figure><img src="image_47.png" /></figure>

<figure><img src="image_48.png" /></figure>

<figure><img src="image_49.png" /></figure>

<figure><img src="image_50.png" /></figure>

<figure><img src="image_51.png" /><figcaption>Figure 39: NASA CRM-HL (Ma 0.2 - Re 5.49M - AoA 17.05°): off-body vorticity contours.</figcaption></figure>
1010

This last test-case demonstrates the natural extension of recent developments to tridimensional industrial applications. By achieving comparable errors to the IBM4 quadtree Cartesian mesh on the 2D multi-element airfoil in terms of skin pressure and friction distributions, as well as lift coefficients, this analysis also implies that further refinement of the mesh around the HL-CRM would certainly improve the IBM solution in the same way as the previous bidimensional test-cases, since this new IBM methodology has been shown to converge towards the reference solutions. Finally, this last application will be further analyzed in a separate article, covering more angles of attack as well as more details of the flow behavior around the geometry. Lower $y_{target}^{+}$ could also be considered based on available computing resources.

1015

## 7. Conclusions

1020

This article has introduced two major improvements to our immersed boundary method. While the work done in [1] had significantly improved the robustness of the approach, these new developments have improved its overall accuracy without compromising its stability or significantly increasing its computational cost. First, an adaptation of the modeling height was proposed using the initial threshold distance provided by the F42 methodology. This new adaptive procedure is based on physical criteria and aims to reduce the proportion of modeled area near the wall wherever possible and necessary. The idea is to better account for the real evolution of the turbulent boundary layer and to prevent the wall model from being applied outside its validity domain as much as possible when the boundary layer is too thin and features an under-undeveloped logarithmic region. In a second part, we exposed the inherent limitations of the common hypothesis regarding the direct extrapolation of the pressure solution from the image point down to the wall. We showed that this assumption is indeed responsible for a non-negligible lack of accuracy in the final pressure extraction, especially near curved regions. It was also observed that these inaccuracies, mainly introduced in the post-processing stage, were actually responsible for larger discrepancies observed in the integral quantities, as they prevented a near-field integration from capturing the correct pressure drag and lift. To overcome this, we developed higher-order linear and quadratic reconstructions of the wall pressure at the end of the simulation using normal pressure gradient information acquired at a new set of image points.

1025

1030

1035

These two developments have been thoroughly validated in section 6. The results obtained on a NACA0012 profile demonstrate an excellent agreement between the IBM solutions and the reference
body-fitted wall-resolved solution for all measured flow quantities. In particular, we have shown that the combined use of the improved pre- and post-processing allows a better capture of the skin pressure distribution near the leading edge, where the local curvature is at its maximum and the boundary layer is at its thinnest. We also demonstrated that the final drag and lift coefficients no longer depend on the integration mean (i.e. near-field or far-field), regardless of the selected angle of attack. The 2DMEA test-case confirmed these promising results, this time on a more complex geometry and for different Cartesian mesh resolutions. While the finest meshes showed excellent agreement with the reference solution, the coarser meshes also validated the contribution of recent developments, with more accurate IBM solutions at lower computational cost. It is also interesting to note that good quantitative results, integrated or not, can be obtained even for large $y_{target}^{+}$ values, which raises the general question of the benefit-cost ratio in an immersed boundary context. These observations were finally confirmed with the study of the HL-CRM in high-lift configuration, for which the current methodology delivered convincing results despite the complexity of this application. Future work will now consist in transposing the current approach to hybrid RANS-LES simulations, as well as in extending the current wall model based on an equilibrium formulation, which inherently limits the current $y_{target}^{+}$ range that can be reasonably imposed.

## Acknowledgement

The authors would like to thank Ivan Mary and Antoine Michael Diego Jost for their many helpful comments and suggestions.

## References

[1] B. Constant, S. Péron, H. Beaugendre, C. Benoit, An improved immersed boundary method for turbulent flow simulations on Cartesian grids, Journal of Computational Physics (2021). doi:10.1016/j.jcp.2021.110240.

[2] J. Slotnick, A. Khodadoust, J. Alonso, D. Darmofal, W. Gropp, E. Lurie, D. Mavriplis, CFD Vision 2030 Study: A Path to Revolutionary Computational Aerosciences, NASA Langley Research Center, CR-2014-218178 (2014). URL https://api.semanticscholar.org/CorpusID:215813979
[3] R. Verzicco, Immersed Boundary Methods: Historical Perspective and Future Outlook, Annual Review of Fluid Mechanics (2023). doi:10.1146/annurev-fluid-120720-022129.

1070

[4] C. S. Peskin, Flow patterns around heart valves: A numerical method, Journal of Computational Physics (1972). doi:10.1016/0021-9991(72)90065-4.

[5] R. Mittal, G. Iaccarino, Immersed Boundary Methods, Annual Review of Fluid Mechanics (2005). doi:10.1146/annurev.fluid.37.061903.175743.

1075

[6] C. Brehm, M. F. Barad, C. C. Kiris, Open rotor computational aeroacoustic analysis with an immersed boundary method, 54th AIAA Aerospace Sciences Meeting (2016). doi:10.2514/6.2016-0815.

[7] F. Capizzano, Turbulent Wall Model for Immersed Boundary Methods, AIAA Journal (2011). doi:10.2514/1.J050466.

1080

[8] S. Wilhelm, J. Jacob, P. Sagaut, An explicit power-law-based wall model for lattice Boltzmann method-Reynolds-averaged numerical simulations of the flow around airfoils, Physics of Fluids (2018). doi:10.1063/1.5031764.

[9] Y. Tamaki, T. Imamura, Turbulent Flow Simulations of the Common Research Model Using Immersed Boundary Method, AIAA Journal (2018). doi:10.2514/1.J056654.

[10] Y. Tamaki, S. Kawai, Wall-Resolved Large-Eddy Simulation of Near-Stall Airfoil Flow at $Re_c = 10^7$, AIAA Journal (2022). doi:10.2514/6.2022-1060.

1085

[11] S.-G. Cai, J. Degrigny, J.-F. Boussuge, P. Sagaut, Coupling of turbulence wall models and immersed boundaries on Cartesian grids, Journal of Computational Physics (2021). doi: 10.1016/j.jcp.2020.109995.

[12] S.-G. Cai, J. Jacob, P. Sagaut, Immersed boundary based near-wall modeling for large eddy simulation of turbulent wall-bounded flow, Computers and Fluids (2023). doi:10.1016/j.compfluid.2023.105893.

1090

[13] N. S. Dhamankar, G. A. Blaisdell, A. S. Lyrintzis, Implementation of a Wall-Modeled Sharp Immersed Boundary Method in a High-Order Large Eddy Simulation Tool for Jet Aeroacoustics, 54th AIAA Aerospace Sciences Meeting (2016). doi:10.2514/6.2016-0257.
1095

[14] J. Degrigny, S. G. Cai, J. F. Boussuge, P. Sagaut, Improved wall model treatment for aerodynamic flows in LBM, Computers and Fluids (2021). doi:10.1016/j.compfluid.2021.105041.

[15] Y. Xu, X. Liu, An immersed boundary method with $y^+$-adaptive wall function for smooth wall shear, International Journal for Numerical Methods in Fluids (2021). doi:10.1002/fld.4960.

1100

[16] F. Capizzano, Automatic generation of locally refined Cartesian meshes: Data management and algorithms, International Journal for Numerical Methods in Engineering (2018). doi: 10.1002/nme.5636.

[17] L. Manueco, P.-E. Weiss, S. Deck, On the estimation of unsteady aerodynamic forces and wall spectral content with Immersed Boundary Conditions, Computers and Fluids (2020). doi:10.1016/j.compfluid.2020.104471.

1105

[18] P. R. Spalart, S. R. Allmaras, A one-equation turbulence model for aerodynamic flows, 30th Aerospace Sciences Meeting and Exhibit (1992). doi:10.2514/6.1992-439.

[19] Fast - CFD Solver (ONERA), https://fast.onera.fr.

1110

[20] N. Alferez, I. Mary, J.-M. L. Gouez, A. Farjallah, C. andreolli, Intel Xeon and Xeon Phi optimizations of an industry-oriented Computational Fluid Dynamics solver, Intel HPC Developer Conference (2017).

[21] A. Suss, I. Mary, T. Le Garrec, S. Marié, Comprehensive comparison between the lattice Boltzmann and Navier-Stokes methods for aerodynamic and aeroacoustic applications, Computers and Fluids (2023). doi:10.1016/j.compfluid.2023.105881.

1115

[22] J. R. Edwards, M.-S. Liou, Low-Diffusion Flux-Splitting Methods for Flows at All Speeds, AIAA Journal (1998). doi:10.2514/2.587.

[23] I. Mary, P. Sagaut, Large Eddy Simulation of Flow Around an Airfoil Near Stall, AIAA Journal (2002). doi:10.2514/2.1763.

[24] P. Roe, Approximate Riemann solvers, parameter vectors, and difference schemes, Journal of Computational Physics (1981). doi:10.1016/0021-9991(81)90128-5.
1120

[25] B. van Leer, Towards the ultimate conservative difference scheme. V. A second-order sequel to Godunov's method, Journal of Computational Physics (1979). doi:10.1016/0021-9991(79) 90145-1.

[26] J. L. Thomas, R. W. Walters, Upwind relaxation algorithms for the navier-stokes equations, AIAA Journal (1987). doi:10.2514/3.9659.

1125

[27] F. Daude, I. Mary, P. Comte, Self-Adaptive Newton-based iteration strategy for the LES of turbulent multi-scale flows, Computers and Fluids (2014). doi:10.1016/j.compfluid.2014. 04.028.

[28] G. Jeanmasson, I. Mary, L. Mieussens, On some explicit local time stepping finite volume schemes for CFD, Journal of Computational Physics (2019). doi:10.1016/j.jcp.2019.07. 017.

1130

[29] K. W. Thompson, Time dependent boundary conditions for hyperbolic systems, Journal of Computational Physics (1987). doi:10.1016/0021-9991(87)90041-6.

[30] S. Péron, C. Benoit, Automatic off-body overset adaptive Cartesian mesh method based on an octree approach, Journal of Computational Physics (2013). doi:10.1016/j.jcp.2012.07. 029.

1135

[31] G. Jeanfaivre, C. Benoit, M.-C. L. Pape, Improvement of the Robustness of the Chimera Method, 32nd AIAA Fluid Dynamics Conference and Exhibit (2002). doi:10.2514/6. 2002-3290.

[32] C. Benoit, G. Jeanfaivre, Three-Dimensional Inviscid Isolated Rotor Calculations Using Chimera and Automatic Cartesian Partitioning Methods, Journal of the American Helicopter Society (2003). doi:10.4050/JAHS.48.128.

1140

[33] Cassiopee - CFD pre-, -co and post-processing tool (ONERA), https://cassiopee.onera.fr.

[34] C. Benoit, S. Péron, S. Landier, Cassiopee: A CFD pre- and post-processing tool, Aerospace Science and Technology (2015). doi:10.1016/j.ast.2015.05.023.

1145

[35] S. Péron, C. Benoit, V. Gleize, I. Mary, M. Terracol, A mixed overset grid/immersed boundary approach for CFD simulations of complex geometries, 54th AIAA Aerospace Sciences Meeting (2016). doi:10.2514/6.2016-2055.
1150

[36] S. Péron, C. Benoit, T. Renaud, I. Mary, An immersed boundary method on Cartesian adaptive grids for the simulation of compressible flows around arbitrary geometries, Engineering with Computers (2020). doi:10.1007/s00366-020-00950-y.

[37] Y. Tamaki, M. Harada, T. Imamura, Near-wall modification of Spalart-Allmaras turbulence model for immersed boundary method, AIAA Journal (2017). doi:10.2514/1.J055824.

[38] A. Musker, Explicit expression for the smooth wall velocity distribution in a turbulent boundary layer, AIAA Journal (1979). doi:10.2514/3.61193.

1155

[39] D. B. Spalding, A Single Formula for the “Law of the Wall”, Journal of Applied Mechanics (1961). doi:10.1115/1.3641728.

[40] F. Capizzano, Coupling a wall diffusion model with an immersed boundary technique, AIAA Journal (2016). doi:10.2514/1.J054197.

1160

[41] M. J. Berger, M. J. Aftosmis, An ODE-based Wall Model for Turbulent Flow Simulations, AIAA Journal (2017). doi:10.2514/6.2017-0528.

[42] L. Manueco, P.-E. Weiss, S. Deck, On the coupling of wall-model immersed boundary conditions and curvilinear body-fitted grids for the simulation of complex geometries, Computers and Fluids (2021). doi:10.1016/j.compfluid.2021.104996.

[43] F. M. White, Viscous Fluid Flow, McGraw-Hill (1974).

1165

[44] P. Lancaster, K. Salkauskas, Surfaces generated by moving least squares methods, Mathematics of computation (1981). doi:10.2307/2007507.

[45] Z.-Q. Cheng, Y. Wang, B. Li, K. Xu, G. Dang, S. Jin, A Survey of Methods for Moving Least Squares Surfaces, IEEE/EG Symposium on Volume and Point-Based Graphics (2008). doi:10.2312/VG/VG-PBG08/009-023.

1170

[46] D. Destarac, CFD-based aircraft drag prediction and reduction, von Kármán Institute for Fluid Dynamics, Lecture Series. Hampton, VA (2003).

[47] C. van Dam, Recent experience with different methods of drag prediction, Progress in Aerospace Sciences (1999). doi:10.1016/S0376-0421(99)00009-3.
1175

[48] B. Mele, R. Tognaccini, Aerodynamic force by Lamb vector integrals in compressible flow, Physics of Fluids (2014). doi:10.1063/1.4875015.

1180

[49] L. L. Kang, L. Russo, R. Tognaccini, J. Z. Wu, W. D. Su, Aerodynamic Force Breakdown in Reversible and Irreversible Components by Vortex Force Theory, AIAA Journal (2019). doi:10.2514/1.J058379.

[50] R. M. Cummings, M. B. Giles, G. N. Shrinivas, Analysis of the elements of drag in three-dimensional viscous and inviscid flows, 14th Applied Aerodynamics Conference (1996). doi:10.2514/6.1996-2482.

1185

[51] T. Renaud, C. Benoit, S. Péron, I. Mary, N. Alferez, Validation of an immersed boundary method for compressible flows, AIAA Scitech 2019 Forum (2019). doi:10.2514/6.2019-2179.

[52] B. Diskin, J. L. Thomas, C. L. Rumsey, A. Schwöppe, Grid convergence for turbulent flows, 53rd AIAA Aerospace Sciences Meeting (2015). doi:10.2514/6.2015-1746.

1190

[53] Langley Research Center Turbulence Modeling Resource, https://turbmodels.larc.nasa.gov/.

[54] S. Allmaras, F. Johnson, P. Spalart, Modifications and clarifications for the implementation of the Spalart-Allmaras turbulence model, Seventh International Conference on Computational Fluid Dynamics (ICCFD7) (2012). URL https://www.iccfd.org/iccfd7/assets/pdf/papers/ICCFD7-1902_paper.pdf

[55] M. L. Shur, M. K. Strelets, A. K. Travin, P. R. Spalart, Turbulence modeling in rotating and curved channels: Assessing the spalart-shur correction, AIAA Journal (2000). doi:10.2514/2.1058.

1195

[56] C. L. Rumsey, J. P. Slotnick, A. J. Sclafani, Overview and Summary of the Third AIAA High Lift Prediction Workshop, Journal of Aircraft (2019). doi:10.2514/1.C034940.

[57] FUN3D - CFD Solver (NASA), https://fun3d.larc.nasa.gov.

[58] AIAA CFD High Lift Prediction Workshop (HLPW), https://hiliftpw.larc.nasa.gov.
1200

[59] A. M. Clark, J. P. Slotnick, N. J. Taylor, C. L. Rumsey, Requirements and Challenges for CFD Validation within the High-Lift Common Research Model Ecosystem, AIAA Aviation Forum 2020 (2020). doi:10.2514/6.2020-2772.

[60] J. C. Duensing, J. A. Housman, L. Fernandes, L. Machado, C. C. Kiris, A Reynolds-Averaged Navier-Stokes Perspective for the High Lift-Common Research Model Using the LAVA Framework, AIAA Aviation Forum 2022 (2022). doi:10.2514/6.2022-3742.

1205

[61] C. Kiris, M. Barad, J. Housman, E. Sozer, C. Brehm, S. Moini-Yekta, The LAVA Computational Fluid Dynamics Solver, AIAA Scitech 2014 Forum (2014). doi:10.2514/6.2014-0070.

[62] A. S. Ghate, G.-D. Stich, G. K. Kenway, J. A. Housman, C. C. Kiris, A Wall-Modeled LES Perspective for the High Lift Common Research Model Using LAVA, AIAA Aviation Forum 2022 (2022). doi:10.2514/6.2022-3434.