# Steady-State Turbulent Flow Simulations of the NASA High Lift Common Research Model in an Immersed Boundary Context

Benjamin Constant, Antoine Michael Diego Jost, Stéphanie Péron, Christophe Benoit, Héloïse Beaugendre

► To cite this version:

Benjamin Constant, Antoine Michael Diego Jost, Stéphanie Péron, Christophe Benoit, Héloïse Beaugendre. Steady-State Turbulent Flow Simulations of the NASA High Lift Common Research Model in an Immersed Boundary Context. AIAA AVIATION Forum and ASCEND 2025, Jul 2025, Las Vegas, United States. 10.2514/6.2025-3791. hal-05190236

HAL Id: hal-05190236

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

Submitted on 28 Jul 2025

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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Steady-State Turbulent Flow Simulations of the NASA High Lift Common Research Model in an Immersed Boundary Context

Benjamin Constant*, Antoine Michael Diego Jost†, Stéphanie Péron‡, Christophe Benoit§

ONERA, Institut Polytechnique de Paris, F-92320, Châtillon, France

Héloïse Beaugendre¶

Institut de Mathématiques de Bordeaux (IMB), Université de Bordeaux, CNRS UMR 5251, Bordeaux INP, INRIA, F-33400, Talence, France

This paper examines several steady-state turbulent flow simulations of the NASA High Lift Common Research Model on Cartesian grids at increasing angles of attack. The Reynolds-Averaged Navier-Stokes equations are solved with the one-equation Spalart-Allmaras model, while the presence of a geometry is accounted for by an immersed boundary method based on a ghost cell direct forcing approach. A wall model is applied to address the inability of isotropic cells to accurately capture the near-wall flow physics at high Reynolds numbers. This work further validates the recent developments brought to the Cartesian methodology, which are based on a geometric regularization of the forcing point location around the immersed geometry. Additional developments are presented to improve the pre- and post-processing for this type of 3D applications. Comparisons are made with state-of-the-art solutions from the reference literature on body-fitted meshes, both structured and unstructured. A 2D multi-element airfoil is also studied to facilitate the analysis of the solutions obtained on the full aircraft in a high lift configuration. For higher angles of attack, more detailed investigations are proposed to better understand the physics closer to the stall phenomenon, especially near the wing tip.

## I. Introduction

In 2024, analyzing the complex flow around an aircraft during landing and takeoff maneuvers remains a major challenge for the aerodynamics industry, addressing important issues such as flight safety and emissions reduction. Landing and takeoff are critical phases during which an aircraft is particularly vulnerable [1]. A thorough assessment of the turbulent flow that develops around the geometry at high angles of attack is therefore essential to better understand the complex interactions between the aircraft and the air that can severely affect aerodynamic performance. Atmospheric disturbances and ground effects must also be studied in detail in the early stages of testing in order to design aircraft capable of reacting quickly and efficiently under extreme conditions [2, 3]. A better understanding of the aerodynamic characteristics of the aircraft during these phases of flight is also important for designing more efficient wings as well as identifying flight paths or takeoff angles that can help minimize fuel consumption at the beginning and end of each flight. Reducing fuel consumption also helps to reduce emissions of carbon dioxide (CO₂) and other greenhouse gases [4], which is important for an industry that aims to be carbon-neutral by 2050 [5-7]. In 2018, air transport was still responsible for 2.4% of global CO₂ emissions from fossil fuels [8].

The use of computational fluid dynamics (CFD) plays a strategic and essential role in achieving these objectives by enabling rapid exploration of advanced aerodynamic designs and evaluating their performance over a wide range of flow conditions [9–11]. CFD has become an integral part of both the private and public aerospace sectors [12], providing a cost-effective and efficient complement to experimental test campaigns [13]. Within the CFD community, international workshops such as the CFD High Lift Prediction Workshop are important platforms for discussing and evaluating the wide range of numerical tools being developed and used to simulate turbulent flows around modern aircraft [14]. In this way, it is possible to find the ensemble of modeling techniques best suited for a particular set of applications, rather than relying on a single overarching model.

<footnote>*Research engineer, DAAA, benjamin constant@onera.fr</footnote>

<footnote>†Research engineer, DAAA, antoine.jost@onera.fr</footnote>

<footnote>‡Research engineer, DAAA, stephanie.peron@onera.fr</footnote>

<footnote>§Research engineer, DAAA, christophe.benoit@onera.fr</footnote>

<footnote>¶Professor, CARDAMOM, heloise.beaugendre@math.u-bordeaux.fr</footnote>
Among the various simulation methods currently available, approaches based on the Reynolds-Averaged Navier-Stokes (RANS) equations have been widely used in the past due to their ability to simulate steady-state turbulent flows with moderate computational demands [15–17]. In fact, most CFD simulations performed for industrial applications still rely on these approaches, as they provide fast and robust results within a day or two. This makes them particularly useful in the design phase, where different wing, engine, or airframe configurations need to be tested quickly. However, RANS approaches, regardless of the model applied near walls, are inherently limited in capturing complex instabilities and separated flow phenomena [18]. Though unsteady RANS (URANS) simulations improve on RANS by accounting for time-dependent flow features, they still struggle to capture fine turbulence structures, as many of these are either smoothed out or approximated [19].

Direct Numerical Simulations (DNS) of full three-dimensional geometries remain inaccessible due to current computational limitations [20]. Thus, the development of more advanced methods, such as Large-Eddy Simulations (LES) [21, 22] and hybrid RANS-LES (HRLES) approaches [23, 24], is crucial to overcome the well-known limitations of RANS and URANS and to achieve a finer understanding of complex turbulent phenomena. LES approaches provide high-fidelity simulations, as close as possible to DNS simulations, by explicitly resolving medium-scale turbulent structures. However, they require finely resolved grids in all three spatial directions, locally approaching near-isotropic refinement [25]. This requirement, along with the high computational power and time needed to achieve simulation convergence for complete aircraft configurations, limits the use of these approaches in industry [20]. In contrast, hybrid RANS-LES models, such as Detached Eddy Simulations (DES) [26, 27] and Zonal DES (ZDES) [28, 29], combine the strengths of both approaches by applying RANS in near-wall regions and LES in the bulk of the flow, offering a promising compromise between accuracy and computational cost. Wall-Modeled LES (WMLES) [30, 31], though more computationally demanding, simplifies boundary layer modeling by approximating near-wall turbulence, enabling LES on complex geometries with reduced computational cost [20].

Finally, other methods such as the Lattice Boltzmann Method (LBM) are gaining increasing interest as alternatives to traditional Navier-Stokes-based methods [32, 33]. Although LBM has historically been used in fields such as automotive aerodynamics, its potential for simulating unsteady turbulent flows around aircraft is currently being explored due to its ability to efficiently handle fluid-structure interactions with Cartesian meshes and immersed boundaries [34–37]. A key advantage of LBM is its ability to accurately predict aeroacoustic behavior on coarser Cartesian meshes, significantly reducing computational requirements [38]. LBM solvers are now systematically represented at every HLPW event [39–41], and the predictive capabilities of LBM Very Large Eddy Simulations (LBM-VLES), which incorporate LES-like concepts within the LBM framework, continue to improve each year [42]. In the CFD community, LBM is recognized as a promising and efficient tool for tackling complex aerodynamic and aeroacoustic challenges, and has already integrated engineering services within international aerospace companies.

As turbulence modeling has evolved, different meshing methods have emerged to address specific model needs and to accommodate the increasing complexity of modern aircraft geometries. Some mesh typologies, such as body-fitted curvilinear meshes, are particularly well suited to represent the near-wall flow physics, but often require considerable manual intervention and expertise, especially for complex geometries. Conversely, unstructured and Cartesian meshes allow for rapid and automatic mesh generation, offering fast numerical setup ideal for early-phase simulations. Cartesian meshes, in particular, have regained popularity in the CFD community due to their efficiency and adaptability in fast CFD solvers, especially with the emergence of LBM. Immersed Boundary Methods (IBM) [43, 44] facilitate fluid-structure interactions on these nonconformal meshes by considering immersed obstacles within the volume mesh, either indirectly through external force fields [45] or directly through localized interpolations or reconstructions [46]. While isotropic Cartesian cells are advantageous for capturing wake-region turbulence, they face limitations in resolving steep gradients near boundaries without exponentially increasing the cell count. For this reason, IBM are often paired with wall models to efficiently handle near-wall turbulence [47]. Over the past decade, IBM have been successfully implemented on complete aircraft geometries in a wide range of flow regimes, including low- and high-speed conditions, and for a variety of CFD models: RANS [48, 49], hybrid RANS-LES [50], LBM [41], and WMLES [51].

In this context, Peron et al. [52, 53], followed by Constant et al. [49, 54], have introduced and developed an efficient Cartesian methodology for RANS applications, based on direct forcing IBM and fast structured solvers. Before extending this approach to more advanced models, such as hybrid RANS-LES or WMLES, this paper aims to evaluate the predictive capabilities of this methodology on complex aircraft geometries approaching the known limits of RANS applicability. The paper is organized as follows: first, the flow solver and IBM features are briefly reviewed. We then present recent developments that facilitate the numerical simulation of external flows over industrial configurations. This methodology section is followed by a numerical analysis of steady-state turbulent flows around 2D multi-element airfoils (2DMEA) and CRM-HL aircraft at increasing angles of attack, using the latest numerical tools at our disposal.
## II. Numerical Method

### A. Description of the flow solver

The steady-state mean-flow solution is obtained by solving the Reynolds-Averaged Navier-Stokes (RANS) equations, which are derived from the following compressible Navier-Stokes equations:

$$ \left\{ \begin{aligned} \frac{\partial \rho}{\partial t} + \frac{\partial \rho u_j}{\partial x_j} &= 0 \\ \frac{\partial \rho u_k}{\partial t} + \frac{\partial (\rho u_k u_j + \delta_{kj} p)}{\partial x_j} - \frac{\partial \tau_{kj}}{\partial x_j} &= 0 \quad (\text{with } k = 1, 2, 3), \\ \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 = (u_1, u_2, u_3)$ 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. Note that the Einstein notation is used on the indices $i$ and $j$ to simplify the previous notations. The system of equations 1 can finally be represented in the following matrix form:

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

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. The equations are then averaged to obtain the RANS equations and closed by the Spalart-Allmaras (SA) turbulence model [16]. The latter is a one-equation transport model that introduces a modified eddy viscosity $\tilde{\nu}$ such as:

$$ v_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 (3) $$

where $\nu$ and $v_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 [16].

The compressible RANS equations are solved using the Navier-Stokes solver FastS [55], which is based on a 2nd-order Finite-Volume Method (FVM) and optimized for block-structured meshes where the flow solution is known at the center of the cells. Using HPC techniques as well as hybrid MPI/OpenMP parallelization, this efficient solver is capable of updating up to $10 \times 10^6$ cells per second on a single Broadwell core [33, 56]. The additional presence of a dedicated Cartesian solver also helps to balance the important number of points required to accurately simulate the turbulent flow field around complex industrial-scale geometries on isotropic meshes.

In the following, the Roe scheme [57] is used to estimate the convective fluxes, using a 3rd-order MUSCL reconstruction [58] to construct the left and right states at each cell boundary. The viscous fluxes are simply approximated with a 2nd-order finite volume scheme based on the compact approach developed by Thomas & Walters [59]. Since the present work focuses only on steady-state RANS simulations, the time integration is handled by a simple and robust 1st-order accurate implicit Euler scheme together with a local time step based on a global CFL number. This approach relies on a Newton's method and uses the Jacobian approximations proposed by Yoon and Jameson [60] for the convective parts and Coakley [61] for the viscous parts.

<figure><img src="image_2.png" /><figcaption>Fig. 1 Stencils for the evaluation of convective (left) and viscous (right) fluxes at a given face $\Gamma_{k,2}$ of an elementary cell $\Omega_k$ with an outward-pointing normal $n_2$.</figcaption></figure>
## B. Description of the Immersed Boundary Method

The FastS solver can operate on both curvilinear meshes and Cartesian grids. The latter allows us to study geometrically complex configurations in a short period of time due to an automatic mesh generation where the mesh does not have to conform to the obstacle. In this case, the presence of the geometry is taken into account through an immersed boundary method (IBM) based on a ghost cell direct forcing approach [53]. The automatic Cartesian grid generation developed by Peron and Benoit [52] is also considered, with some adjustments discussed in the next section. This approach is based on an octree skeleton, which ultimately results in a set of structured uniform Cartesian grids with minimal overlap between adjacent blocks. Three parameters ultimately control the complete mesh generation from the surface geometry alone, namely the domain extent $d_{ext}$, the set of near-wall resolutions $h_{min}$, and the minimum number of cells per direction for each elementary Cartesian block $v_{min}$. The meshing difficulty is thus shifted from the generation of the volume mesh to the generation of the surface mesh, structured or not, from the CAD definition.

<figure><img src="image_3.png" /><figcaption>Fig. 2 Illustration of the direct forcing approach for a given computed point P and a set of IBM points (W, A, B) for RANS or hybrid RANS/LES simulations.</figcaption></figure>

Figure 2 summarizes the IBM developed by Peron et al. [53]. This approach consists in forcing the flow solution at some selected points close to the immersed geometry, called target points, using the interpolated flow solution obtained further away from the wall at different probes, called image points. One of the original features of this method remains the use of an irregular step front defined by the first computed points for the projection of these image points. In more recent work, Constant et al. [54] have shown that this technique can be used for accurate aerodynamic predictions while preserving a smooth behavior of the near-wall flow solution, especially on the skin coefficients. Given the numerical solver used in this work, which is based on a five-point stencil in each direction, at least two rows of IBM target points must be defined near the wall prior to simulation. Finally, the wall points, shown in gray in Figure 2, are used to recover the surface flow solution during or after the simulation using a third-order Moving Least-Square (MLS) interpolation technique.

To simulate turbulent flows at realistic Reynolds numbers on the order of $10^6$, this Cartesian method should be coupled with a wall model to limit the number of cells within the volume mesh. In this work, Musker's wall function [62] is applied in the modeling region to recover both the tangential velocity at the target points and the friction velocity $u_τ$ at the wall points. From the flow solution $\dot{W}_B^n = (\rho_B^n, u_B^n, T_B^n, \bar{v}_B^n)$ reconstructed by linear interpolation at a given image point B, the flow solution $W_A^n = (\rho_A^n, u_A^n, T_A^n, \bar{v}_A^n)$ at its corresponding target points A can be obtained by applying the following formulas. The tangential velocity $u_A$ is expressed by:

$$u_A = u_τ M(y_A^+), \quad \text{where} \quad y_A^+ = \frac{u_τ y_A}{\nu}, \qquad (4)$$
where M denotes the Musker’s algebraic wall function at equilibrium:

$$M(u_{\tau}) = 5.424 \underbrace{\arctan\left[\frac{2y^{+} - 8.15}{16.7}\right]}_{f_{1}^{M}(u_{\tau})} + \underbrace{\log\left[\frac{(y^{+} + 10.6)^{9.6}}{(y^{+2} - 8.15y^{+} + 86)^{2}}\right]}_{f_{2}^{M}(u_{\tau})=\log[g_{1}(u_{\tau})/g_{2}(u_{\tau})]} - 3.52. \quad (5)$$

The normal component of the velocity is then simply estimated by linear reconstruction:

$$v_A = \frac{y_A}{y_B} v_B. \quad (6)$$

For the pseudo-viscosity $\tilde{v}$, Peron et al. [53] showed that this quantity derived from the SA turbulence model can be evaluated by solving the following nonlinear equation using Ferrari’s method:

$$(\tilde{v}_A)^4 - \kappa u_{\tau} y_A D (\tilde{v}_A)^3 - \kappa u_{\tau} y_A D v^3 C_{v1}^3 = 0, \quad (7)$$

where $\kappa = 0.41$ is the Von Kármán constant and $D = [1 - \exp(-y_A^+/19)]$ is the Van Driest damping term. Finally, the temperature, density, and pressure information are obtained using a Crocco-Busemann relation [63], the ideal gas law, and a direct extrapolation of the pressure solution in the wall normal direction, respectively:

$$\left\{ \begin{array}{l} T_A = T_B + \frac{\sqrt{Pr}}{2C_P}(u_B^2 - u_A^2) \\ \rho_A = \frac{p_A}{RT_A} \\ p_A = p_B \end{array} \right. \quad (8)$$

where $Pr$ is the Prandtl number, $C_P$ is the specific heat coefficient at constant pressure, and $R$ is the ideal gas constant.

Given the complexity of the applications discussed in these pages, related to aerodynamic geometries in high lift configurations facing strong incidences, the recent developments brought to the IBM by Constant et al. [49] are considered and applied throughout the results section. This refers to an extension of the geometric approach introduced in 2021 [54], which further improves the near-wall flow prediction by adjusting the location of the target points to account for the real boundary layer thickness, as well as taking into account normal pressure gradients to accurately reconstruct the wall pressure in post-processing, thus reducing the error made in the prediction of the aerodynamic coefficients. We recall that both of these approaches operate with a modeling height defined around the immersed obstacles, which is related to the flow parameters by the following formula:

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

where $y_0$ is the height of the first cell on a flat plate for which $y^+ = 1$ and $f(Re) = 0.058Re^{-0.2}$.

### III. Simulation workflow for external aerodynamic studies on industrial level configurations

This section presents the latest developments and improvements brought to the IBM approach to simulate challenging geometries and flow behaviors. This work can be seen as complementary to the work done in 2021 and 2024 regarding the accuracy of the near-wall treatment for turbulent flow prediction within the boundary layer [49, 54]. In this sense, several specific pre- and post-processing operations are presented in the following subsections.

#### A. Test case setup with a symmetry plan

Simulating only half of a geometry is very common for complex 3D applications involving a large number of points inside the volume mesh, especially around wings or aircraft when a natural symmetry plane already exists. In 2021, a symmetry plane was considered to simulate transsonic flows past a ONERA M6 wing via a slip immersed boundary condition [54]. Although this approach was easy to implement, the IBM treatment would require a minimum resolution near the entire symmetry plan to properly apply the immersed boundary condition. This technique could also be responsible for numerical instabilities near the outer boundaries of the domain, where the immersed boundary condition is in contact with the far-field boundary condition handled directly by the solver, but also near the geometry, where the slip and wall conditions locally intersect.
These issues have been resolved by applying a conventional symmetry condition instead of a slip immersed boundary condition. Prior to mesh generation, the full geometry is first symmetrized along one of the three directions in space, classically the y-axis. A box is then automatically generated around the symmetrized part of the geometry, larger than the final domain, so that the symmetry condition crosses the adjacent boundaries. This symmetrized problem is finally used to generate the octree as well as to naturally close the aircraft configuration near the symmetry plan for the following blanking operation within the mesh. By not accounting for this boundary condition through immersed boundaries, the octree mesh can also naturally coarsen away from the wall, reducing the total number of cells. Finally, this second methodology completely eliminates the numerical oscillations that could occur near the wall and domain boundaries.

<figure><img src="image_4.png" /><figcaption>**Fig. 3 Set of boundary conditions for a given turbulent flow simulation around the half-plane Common Research Model in High-Lift configuration (CRM-HL) in a freestream environment.**</figcaption></figure>

Figure 3 provides an illustration of this type of domain configuration for the CRM-HL, displayed in full meters. The domain extent $d_{ext}$ is set to 100 times the length of the aircraft and has been applied in all three directions of space to create the initial octree. The symmetry plan is positioned at $y = 0$ to split the geometry along the fuselage. The Cartesian mesh is created on half of the initial octree skeleton, using a simple cell selection based on the y-coordinate. The final domain can be seen as a square cut in half, with a symmetry boundary condition on one face (in red) and non-reflecting far-field boundary conditions on the remaining faces (in blue). All wall boundary conditions are ultimately taken into account by immersed boundary conditions, using the modeling ingredients recalled in the previous section.

## B. Mesh optimization with refinement boxes

It is well known that the use of Cartesian grids tends to exponentially increase the number of cells in 3D, even more so in a block-structured framework. Although the Cartesian methodology presented above has been thoroughly optimized to efficiently handle large meshes, additional efforts to further reduce the number of cells are still important when the mesh resolution cannot be too coarse. This is particularly the case for complex flow topologies where the wall model is partially inadequate to represent the real flow physics near the wall. In this sense, recent developments have aimed at reducing the extent of the minimum near-wall spacing inside the mesh by using local mesh refinements around the geometry as well as refinement boxes to limit the changes in mesh resolution in some flow regions.

During the octree generation step, some near-wall grids may be predominantly positioned inside the obstacle. Without additional treatment, it is then possible to observe a rapid and unwanted change in resolution in the wall normal direction, which can negatively affect the application of our IBM. To overcome this problem, the generation step is usually followed by an extension of the lowest-level nodes from the geometry definition, which proved to be robust enough for the original IBM with two rows of target points. However, the recently introduced responsive IBM solution further extends the modeled region near the immersed geometries, which means that the extent must also be increased to ensure that all IBM points, as well as the first row of computed points, are ultimately located within the finest mesh resolution. The first and simplest idea is to increase the number of lowest-level nodes by two layers in the three directions of space from the geometry surface. This solution is illustrated in Figures 4-(b) and 5-(b), for different $v_{min}$ but constant $h_{min}$ around a circular obstacle represented with a solid line. For a small number of cells per elementary block (Figure 4), this approach is not able< to propagate the minimum resolution over the entire modeling domain bounded by a
dashed line. Conversely, for a larger $v_{min}$ value (Figure 4), the extent around the circle is now too important. To avoid these problems, a more adaptive solution is proposed, this time based on refinement boxes that encapsulate the entire modeling domain. These boxes can be automatically generated around the entire immersed geometry, or just some preselected zones, with an offset that uses a distance parameter based on the modeling height. As can be seen in Figures 4-(c) and 5-(c), the near-wall resolution can now be sufficiently propagated from the surface while keeping the number of new cells to a minimum.

Octree refinement via refinement boxes has been widely used in the Cassiopee/Fast framework for more than a decade. This operation is actually based on a blanking algorithm, the same algorithm that is used to blank the resulting Cartesian mesh with all IBC bodies during the IBM pre-processing [64]. Apart from the systematic use of these bodies around arbitrary geometries to sufficiently propagate the finest resolution away from the wall definition, the novelty of the current approach lies in applying the blanking operation at the octree nodes, rather than at the centers, to refine all cells within the refinement bodies. Further subdivision of any octree child is then toggled if and only if at least one node falls within a refinement body.

<figure><img src="image_5.png" /><figcaption>Fig. 4 Automatic octree generation around a circular obstacle with $h_{min} = R/100$ and $v_{min} = 7$.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 5 Automatic octree generation around a circular obstacle with $h_{min} = R/100$ and $v_{min} = 14$.</figcaption></figure>
### C. Automatic and robust projection regularization for IBM points around complex geometries

During the IBM pre-processing, the search for the image and wall points follows the methodology introduced by Peron et al. in 2020 [53]. Starting from the step front of the first computed points, the idea is to project each target point towards the geometry using of wall distance gradients to obtain the wall point, and then follow the resulting normal in the opposite direction to find the corresponding image point. However, this approach alone is not robust enough to avoid incorrect projections in some cases. In fact, after this first projection attempt to the geometry, some wall or image points may end up very far apart if their corresponding target point is located in a highly concave region, resulting in a poor reconstruction of the near-wall solution.

Figure 6 illustrates this problem, in the case of a geometry that resembles a backward facing step. While the majority of the target points are successfully projected onto the wall and the image front after this first attempt (Figure 6-(a)), the two marked target points have actually failed to find suitable candidates for their corresponding image points. In fact, their normal vectors are perfectly tangent to the wall, which significantly shifts the intersection with the front of the first calculated points (Figure 6-(b)). The solution here is to consider an orthogonal projection onto this front instead of the geometry, allowing the IBM triplets to be positioned in a restricted environment (Figure 6-(c)).

<figure><img src="image_7.png" /><figcaption>**Fig. 6** Illustration of the problems associated with concave regions when projecting target points onto the wall.</figcaption></figure>

In practice, a tolerance is considered for the maximum projection distance of the target points to the geometry and the image front. This tolerance is the same for both projections and can be determined either by the local near-wall resolution for fixed layers of target points ($d_1 = n_{layers} \cdot \sqrt{dim} \cdot h_{loc}$), or by the modeling height for adaptive methods ($d_2 = 2 \cdot H_{mod}$). The two projection techniques are then used consecutively to ensure that each target point is correctly projected in both directions. In the few cases where the two projections do not yield three aligned IBM points given the previously imposed tolerance, the minimum projection distances are used to locate both the wall and image points. This choice is based on the consideration that the error caused by three unaligned IBM points is actually less dramatic than applying the wall model at the target points with an incorrect wall distance, or using flow information from distant locations to feed the wall model.

### D. Efficient co- and post-processing in an immersed boundary context

The IBM post-processing has recently been optimized to allow efficient co-processing of aerodynamic loads during the simulation. This involved modifying the algorithm described in [49, 53] for extracting the surface flow solution onto the immersed geometries to actually operate in a shared memory context. Prior to simulation, the full surface projection is created from the case geometry by removing some of the zones based on their immersed boundary conditions or family names, before splitting and redistributing the initial triangular surface mesh among all MPI processes. It is important that each zone of the surface definition has access to the flow solution stored in its neighboring IBM wall points, in order to accurately reconstruct the flow solution via third-order MLS interpolation.

Since this new surface mesh distribution does not match that of the volume mesh, a new communication graph must be initialized to successfully perform the final projection of the flow solution onto the wall. Note that it would have been possible to avoid this extra step by splitting the surface geometries according to the pre-existing distribution in the volume mesh (Figure 7-(a)). However, this would have caused a significant imbalance in the projection step, as some
processes would have to handle larger portions of the surface geometry than others, ultimately slowing down the whole
process. Similarly, it is not entirely certain that the amount of communication between procs would have been reduced.
As can be seen in Figure 7-(a), many sections of the obstacle are actually crossed by multiple Cartesian zones located on
different processors, meaning that a global MPI send-receive communication would still have been necessary. The new
distribution is shown in Figure 7-(b), where each portion of the circular obstacle receives a copy of the IBM information
from the adjacent near-wall Cartesian blocks. Since the Cartesian blocks involved in this IBM step feature the finest grid
resolutions of the entire volume mesh, note that only the flow information located at the wall points is actually shared
between the processors to avoid the risk of memory overflow.

<figure><img src="image_8.png" /><figcaption>**Fig. 7** Cartesian block distributions over four MPI nodes around a circular obstacle. Distributions for simulation (a) and IBM post-processing (b) Each color represents a single MPI process. Only the Cartesian blocks that intersect the circle are shown.</figcaption></figure>

The bottleneck of the entire IBM post-processing is actually the search for interpolation coefficients to extract the surface flow solution. For this reason, we have added the ability to store the MLS interpolation coefficients prior to simulation to ensure fast and efficient co- and post-processing. Recall that in order to accurately interpolate the solution at a given vertex, at least 10 wall points must have been projected onto adjacent surfaces to perform the final MLS interpolation. Otherwise, the surface solution is simply extrapolated from the nearest wall point. Suitable candidates are usually found by using bounding boxes and orthogonal projections of each wall point onto the surface mesh, as well as a nearest neighbor algorithm using *k*-d tree structures to precondition the search. To avoid prolonging this important but tedious step, it is important that each process only holds adjacent regions of the immersed boundaries. Otherwise, the search for suitable candidates would take more time, as distance wall points would unnecessarily try to project onto a surface they didn't originally project onto. This problem, which seems to be inherent to *k*-d tree-based search, where the overall speed of the algorithm is actually limited by distant candidates, is also the reason why wall distance computation can often be an important part of the entire pre-processing step.

Compared to previous implementations, this new methodology also provides much more flexibility. During co- or post-processing, the flow solution can easily be reconstructed in parallel using the stored information at low computational cost. The capture of the aerodynamic coefficients is then divided into two steps. The first is the integration of the pressure and friction information in the body frame, which can be done separately per zone or per zone family. This can be particularly useful, for example, to understand which part of the geometry is responsible for most of the drag generation. The coefficients are then normalized using standard flow parameters, namely static pressure and reference surface or length information, and reintegrated into the wind frame.
In addition, the approach for extracting the wall pressure with first and second order formulations has been slightly modified and simplified. This additional post-processing procedure follows the methodology introduced in [49] and involves a second set of image points located further into the fluid domain. Since the IBM data is stored from the donor viewpoint, which is intended to simplify the interpolation step at the image points and reduce communication, the wall point distribution around the immersed geometry is not exactly the same for the first and second image points. Therefore, two MLS projections must be performed on the surface, one using the first image points for friction and velocity information and one using the second image points for pressure information. With the latest developments, the additional cost of this improved IBM post-processing actually comes from calculating the pressure gradients and updating the wall pressure solution using the following formulas:

$$ \left\{ \begin{aligned} p_{w,o1} &= p_{B_2} - \left. \frac{\partial p}{\partial \mathbf{n}} \right|_{B_2} \cdot \left\| \overrightarrow{WB_2} \right\| \\ p_{w,o2} &= p_{B_2} - \left. \frac{\partial p}{\partial \mathbf{n}} \right|_{B_2} \cdot \left\| \overrightarrow{WB_2} \right\| + \left. \frac{\partial^2 p}{\partial nn} \right|_{B_2} \cdot \frac{\left\| \overrightarrow{WB_2} \right\|^2}{2} \end{aligned} \right. \quad (10) $$

where $B_2$ refers to the second image point, while $\partial_n p$ and $\partial_{nn} p$ are the first and second derivatives of the pressure in the wall normal direction $\mathbf{n}$. Although this procedure is crucial for faithfully predicting the correct wall solution at the wall in our framework, especially near curved regions or in the presence of thin turbulent boundary layers [49], this additional extraction is not considered to follow the convergence of the flow solution through the evolution of the lift coefficient during the simulation. Instead, due to their incompressible cost, these higher-order extractions are only used during the IBM post-processing, while we simply observe $p_w = p_{w,o0} = p_{B_1}$ during co-processing. Figure 8 below illustrates the two steps of an improved IBM post-processing involving two image points.

<figure><img src="image_9.png" /><figcaption>Fig. 8 Interaction between two image points to recollect all flow information at the wall during post-processing.</figcaption></figure>

Finally, let us recall the definitions of the skin pressure and friction coefficients, hereafter denoted $C_p$ and $C_f$:

$$ 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}, \quad (11) $$

where $\cdot_w$ and $\cdot_\infty$ refer to the wall and freestream quantities, respectively. The friction coefficient is calculated using the friction velocity determined by the algebraic wall model. The wall temperature and density are calculated using the adiabatic Crocco-Busemann formula and the ideal gas law, assuming zero velocity at the wall. The integration of all aerodynamic forces and moments acting on the bodies is then performed by direct near-field integration. In fact, previous work has shown that the differences between near-field and far-field integration are actually marginal when considering higher-order wall pressure reconstructions, which is the reason why we will consider only the former to study the evolution of drag, lift, and pitching moments.
IV. Numerical Results

This section presents a detailed analysis of subsonic flow simulations around a complete aircraft geometry in a
high lift configuration, representative of takeoff and landing phases, using the numerical tools and recent developments
described in the previous sections. The numerical setup and flow parameters follow some of the test cases presented
at the fourth edition of the international AIAA CFD High Lift Prediction Workshop (HLPW) [65], which provides
a comprehensive overview of the current prediction capabilities of modern CFD solvers, based on realistic aircraft
geometries facing strong flow incidences, and establishes best practice guidelines for the future.

The first three editions focused on simplified wing geometries (three-element NASA Trapezoidal Wing [66]) and
several aircraft models from leading aerospace research centers (DLR-F11 [67], NASA-CRM [68], and JAXA Standard
Model [68]). For the fourth edition [69], the workshop introduced a revision of the well-known NASA High Lift
Common Research Model (CRM-HL), in its full version, which is the result of a close collaboration between Boeing and
NASA to create a reference open geometry that reflects current commercial aircraft [70]. In recent years, many variations
of the CRM-HL model have been tested both numerically and experimentally, mostly by modifying the deployment of
slats and flaps positioned around the wing, as well as their attachment systems. This application represents a major
challenge for most of numerical CFD solvers, especially when approaching the maximum lift coefficient, as the flow
becomes massively detached near the wing tip. Additional difficulties include inboard and outboard flow separations
and wake interactions from the nacelle pylon, which are also sources of discrepancies in all numerical solutions.

Clark et al. [71] and Duensing et al. [72], recalled that these applications are even more challenging for steady-state
RANS methods, since no turbulence model seems fully capable of representing the complex nature of the flow between
19° and 22°, i.e. near maximum lift ($C_{Lmax}$). In particular, the prediction of aerodynamic coefficients on the CRM-HL
is highly dependent on the fraction of separated flow measured on the wing suction side, which is notoriously difficult to
predict using RANS approaches, despite increasingly complex models. For this workshop, which focused on simpler
one-equation Spalart-Allmaras (SA) turbulence models, no closure formulation seemed to stand out clearly from the
others. Although some advanced SA corrections improved the overall flow predictions with respect to experimental
results in some specific regions, they also tended to underperform in others. As a result, most participants adopted the
standard SA or SA-neg turbulence models as their best practice choice for RANS computations, viewing them as the
best current compromise between robustness, accuracy, and computational efficiency.

Prior to the HLPW4, participants were organized into six Technology Focus Groups (TFG) to draw broad conclusions
about different families of numerical methods. These included a variety of meshing techniques (structured curvilinear
grids, unstructured grids, adapted unstructured grids, Cartesian grids, etc.), turbulent modeling strategies (RANS,
URANS, WMLES, HRLES, etc.), and discretization techniques (finite volume methods, high order finite element
methods, lattice Boltzmann methods, etc.). Comparative analyses have been performed separately [73–77] and
synthesized in a final meta-analysis [69]. This edition of the HLPW has confirmed the limitations of RANS approaches
when approaching $C_{Lmax}$, with clear examples of underconvergence and significant scatter in the solutions with respect
to flow topologies [73]. While the use of adapted grids improves the consistency of the solutions to a large extent, the
participants of the mesh adaptation TFG also showed deficiencies compared to experimental results at high angles
of attack, which inevitably affected the aerodynamic coefficients [69]. In addition, the specific influence of different
numerical parameters (CFL, spatial discretization, etc.), different mesh refinements, or different flow initializations is
still partly unknown and will require further testing to develop common best practice guides for this type of application.
On the other hand, scale resolving methods such as LBM, HRLES, and WMLES offered better consistency with respect
to the reference solutions and showed promising results that seemed to eliminate some of the shortcomings of RANS
methods [76, 77]. Future workshops are expected to confirm this conclusion with even more participants.

Although lattice Boltzmann methods also rely on Cartesian grids and immersed boundary methods coupled with
wall modeling techniques to represent high Reynolds number turbulent flows, the current study focuses solely on the
results provided by the RANS TFG. In fact, it was considered preferable to compare our own simulations with reference
steady-state RANS simulations using body-fitted wall-resolved meshes, both structured and unstructured, rather than
unsteady LBM-based simulations, in order to limit the differences to the treatment of near-wall flows. In addition, the
LBM contributions for this fourth edition were limited to two participants, with no clear consensus at high angles of
attack, while the fixed-grid RANS TFG received contributions from 18 participants, providing a broader and more
consistent reference dataset. In the following, we primarily limit our analysis to two angles of attack, 7.05°and 17.05°,
to evaluate the performance of our Cartesian methodology under both favorable and challenging flow conditions within
a RANS-SA framework. The upper limit of 17.05°is intended to remain within a range where reference body-fitted
wall-resolved RANS solutions from the literature have been shown to provide accurate and reliable results, i.e., away
from the stall phenomenon.
## A. Turbulent flow simulations past a 2D multi-element airfoil

To better understand the IBM results obtained on the full aircraft geometry, we first propose to briefly discuss the characteristics of the turbulent flow established around a 2D multi-element airfoil (2DMEA) at 16° of incidence. This two-dimensional case, documented on the NASA TMR website [78], is actually a cross section of the CRM-HL geometry taken about 16 meters from the wing root. It has been used to analyze several mesh effects in an unstructured context, mainly for compressible RANS approaches based on SA [79] and SA-RC [80] turbulence models. To make this study as representative of the real case as possible, we consider a Mach number of 0.2 and a Reynolds number of 5 × 10⁶, based on the chord length of the retracted profile. Note that this 2D test case was also part of the fourth HLPW, as the third test case, and was intended to perform a verification study on the standard Spalart-Allmaras turbulence model to highlight potential implementation differences that would cause greater discrepancies among the participants of the RANS TFG on the CRM-HL.

Since this test case has already been thoroughly studied by Constant *et al.* [49] to analyze the cross-benefits of a local modeling height and higher-order wall pressure reconstructions in post-processing, the present study will only analyze three IBM solutions with three different mesh resolutions and corresponding $y_{target}^{+}$, after modeling height adaptation and wall pressure correction via second-order formulation. In addition, we now observe two resolution changes on the main element to get as close as possible to the 3D test case, while demonstrating that the current Cartesian methodology is effortlessly able to manage 1:2 abutting grids without generating near-wall perturbations. The first near-wall resolution (IBM4) was chosen to position the first computed cells at $y_{target}^{+} = 400$, which is representative of the resolutions that can realistically be applied in 3D without drastically increasing the number of cells. The last two resolutions (IBM5 and IBM6) allow smaller $y_{target}^{+}$ of 200 and 100, respectively, and are intended to discuss the capabilities of the method under optimal conditions and to demonstrate mesh convergence. Unlike in [49], the number of cells has been fully optimized, using the latest developments and a smaller $v_{min}$ value. This explains why the mesh sizes are now 2.5 times smaller than the Cartesian meshes used in previous studies.

In this case, the reference solution comes from the unstructured FUN3D solver developed by the NASA LRC, which has similar numerical ingredients to FastS. In the following, we consider only the body-fitted solution obtained on the finest mesh, i.e. the hex-dominant L7 unstructured mesh, from the first family of meshes proposed on the NASA TMR website (with first-order upwinding for the advective terms) [78]. The main features of the four meshes are summarized in Table 1. The percentages of cell types are also given for our own Cartesian meshes, with little difference from one resolution to another. Each mesh contains approximately 20% of interpolated cells (IBM target cells and Chimera ghost cells), 20% of blanked cells inside the geometry, and 60% of computed cells. Figure 9 then proposes some views of the IBM4, IBM5, and IBM6 Cartesian meshes in the right column, along with equivalent views of their respective skeleton quadtrees in the left column. Figures 9-(a,c,e) show the refinement boxes automatically created around the 2DMEA to sufficiently propagate the minimum resolution from the geometry to the boundary layer. These boxes, located around the slat, flap, leading edge and trailing edge of the wing, are indicated by colored dashed lines. For ease of reading, the same color code is maintained in the following figures, with IBM4, IBM5, and IBM6 results shown in green, orange, and blue, respectively, while the reference solution from FUN3D is shown in gray. Local magnifications over the main wing near two resolution changes, from fine to coarse near the leading edge and then back to coarse to fine near the trailing edge, help visualize the IBM point distribution in these regions in figure 9-(a) for the IBM4 Cartesian mesh. As shown in [54], the use of a modeling height naturally regularizes the location of the image points without a clear jump from one resolution to another, which would have affected the near-wall solution.

<table><caption>Table 1 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>BF L7</th></tr></thead><tbody><tr><td>Grid topology</td><td>Cartesian</td><td>Cartesian</td><td>Cartesian</td><td>Unstructured (Hex)</td></tr><tr><td>Minimum spacing h<sub>min</sub>/c</td><td>4e-4</td><td>2e-4</td><td>1e-4</td><td>N/A</td></tr><tr><td>Corresponding y<sub>target</sub><sup>+</sup></td><td>400</td><td>200</td><td>100</td><td>0.25</td></tr><tr><td>Number of points</td><td>0.56M</td><td>1.05M</td><td>2.01M</td><td>8M</td></tr><tr><td>% of cell types (0 / 1 / 2)</td><td>17% / 64% / 19%</td><td>21% / 60% / 19%</td><td>23% / 59% / 18%</td><td>-</td></tr></tbody></table>
<figure><img src="image_10.png" /><figcaption>Fig. 9 2DMEA (Ma 0.2 - Re 5M - AoA 16°): IBM mesh views around the geometry.</figcaption></figure>

Satisfactorily, the IBM solution behavior observed in 2024 is recovered with these lighter Cartesian meshes. All IBM solutions are in very good agreement with the reference solution, whether at the wall (Figure 10) or across the boundary layer (Figure 11). Figure 10 shows the mesh convergence on the pressure and friction coefficients, from the coarsest to the finest Cartesian mesh. For each IBM solution, the transitions from one resolution to another are almost imperceptible. Only the transition from fine to coarse, located at $x/c = 0.1$, seems to have a small impact on the skin friction distribution, as it leads to a slight accentuation of the lack of relaxation towards the reference solution, which was already observed in [49]. On the other hand, the pressure profile remains absolutely unaffected by the change in resolution, despite the use of first and second order pressure gradient information in the final reconstruction. A closer look at the IBM6 solution shows that it is almost perfectly superimposed on the FUN3D solution for the pressure coefficients, missing the friction peak on the slat by only a few percent. For its part, the IBM5 solution accentuates the limitations observed with IBM6, with an additional underestimation of the friction peak on the main element. Overall, the IBM5 solution remains close to its finer counterpart. Finally, the IBM4 solution shows an underprediction of the friction near each leading edge, but still provides a satisfactory solution with respect to the reference. The use of a larger modeling height also tends to dampen the overall behavior of the skin pressure distribution, while still following the correct trend and successfully approaching all three pressure peaks. While the modeling height is locally adapted so that the real $y^+$ measured at the first computed cells does not exceed the initial $y_{target}^+$ value, we know that the boundary layer at the slat front is actually too thin to allow the use of wall models that far from the wall, thus overestimating the boundary layer thickness for the entire profile. Nevertheless, the accuracy remains quite good with the coarsest mesh, which keeps the error relatively small at a low computational cost.
<figure><img src="image_11.png" /><figcaption>(a) Skin pressure coefficient at $y^+_{target} = 400$</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>(b) Skin friction coefficient at $y^+_{target} = 400$</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(c) Skin pressure coefficient at $y^+_{target} = 200$</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>(d) Skin friction coefficient at $y^+_{target} = 200$</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>(e) Skin pressure coefficient at $y^+_{target} = 100$</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>(f) Skin friction coefficient at $y^+_{target} = 100$</figcaption></figure>

<caption>Fig. 10 2DMEA (Ma 0.2 - Re 5M - AoA 16°): skin pressure and friction distributions.</caption>
Away from the wall, the x-velocity profiles taken on the multi-element suction side at constant abscissa also shed some relevant light on the capabilities of each mesh refinement. In Figure 11, both IBM5 and IBM6 solutions appear to be in excellent agreement with the FUN3D solution, which also confirms previous observations. While the finest solution is able to further reduce the error made on the first profile, at $x/c = -0.03$, over the slat where the boundary layer is too thin to have a sufficiently developed logarithmic layer, the IBM5 mesh seems to be almost as good at representing the complex flow behavior of the slat wake interacting with the boundary layer on the last two sections. Analysis of the IBM4 solution confirms that the boundary layer at the top of the slat is overestimated, which naturally affects the solution further downstream. In fact, the differences with FUN3D on the first section seem to evolve linearly from one mesh refinement to the next, which means that the level of error remains under control. Nevertheless, the IBM4 solution still offers good qualitative results in the slat wake of each section, while approaching the reference satisfactorily when the boundary layer is finally thick enough to allow the imposition of the wall model within its validity range. This can be seen in the first 10% of the velocity profile for section 2 and in the first 20% for section 3. Under these flow conditions, it is clear that the IBM6 and IBM5 meshes remain our best practice choice, but the IBM4 is still a viable alternative when the finest resolutions cannot realistically be used.

Figure 12 supports the previous conclusions by quantifying the error in the aerodynamic coefficients. The finest IBM5 and IBM6 solutions show excellent agreement with the reference body-fitted wall-resolved solution for the lift, drag, and pitching moment coefficients, respectively. For the IBM4 mesh, the drag coefficient error remains significant due to its higher sensitivity, as drag is typically an order of magnitude smaller than lift. In this 2D test case, only the finest resolutions are capable of capturing drag with sufficient accuracy in the absence of compensating effects for the excess drag generated upstream, with the IBM6 Cartesian mesh yielding an error of only 5%. It remains unclear whether improved wall modeling alone could significantly reduce the error without additional mesh refinement, since most analytical wall laws assume a turbulent boundary layer at equilibrium, an assumption that is not satisfied at the leading edge of an airfoil where the boundary layer is still developing and locally infinitesimal.

Although the pitching moment error is smaller than the drag error, it remains relatively high for IBM4 as the dominant source of error lies on the slat, far from the center of integration. In contrast, the lift coefficient is well predicted across all mesh levels. Its error decreases slowly and linearly with mesh refinement, further supporting the use of the coarsest resolution to analyze the CRM-HL facing high angles of attack.

<figure><img src="image_17.png" /><figcaption>Fig. 11 2DMEA (Ma 0.2 - Re 5M - AoA 16°): x-velocity profiles adimensionalized by freestream velocity norm.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>Fig. 12 2DMEA (Ma 0.2 - Re 5M - AoA 16°): aerodynamic coefficients (lift, drag, and pitching moment).</figcaption></figure>

Finally, Figure 13 illustrates the mesh convergence criterion used for the subsonic flow simulations around the 2DMEA, based on the evolution of the lift coefficient monitored in co-processing. The simulation is performed in two steps. In the first step, a global modeling height is defined using Equation 9 to initialize the flow field and compute the boundary layer thickness. This value is then locally adapted in the second step to more accurately reflect the actual boundary layer evolution. The simulation is then resumed on the same mesh, starting from the previous steady-state solution, until a new convergence of the lift coefficient is reached.

Analysis of the convergence rate for these three Cartesian meshes provides valuable insight into the applicability of each near-wall resolution in 3D. On the finest mesh, convergence requires approximately twice as many iterations as on the coarsest mesh, due to both a finer near-wall resolution and a smaller modeled area around the geometry, as discussed in [81]. This convergence rate becomes another critical factor to weigh, alongside CPU time and memory footprint, when selecting a near-wall resolution for 3D applications. For all these reasons, we chose a $y_{target}^{+}$ of 500 for the $C_{Lmax}$ study of the CRM-HL as the best practical compromise between accuracy and computational efficiency. Although lower modeling heights would likely improve predictions in complex regions where the wall model performs poorly, the mesh resolution they require was not realistically applicable within our Cartesian methodology at the time of this study.

<figure><img src="image_19.png" /><figcaption>Fig. 13 2DMEA (Ma 0.2 - Re 5M - AoA 16°): iterative convergence for the lift coefficient. Circle markers represent the iteration for which the stop criterion based on the L2 norm was reached.</figcaption></figure>
## B. Turbulent flow simulations past a NASA CRM-HL

Following the analysis of the 2DMEA facing a 16° turbulent flow, this section investigates the subsonic turbulent flow around the NASA CRM-HL at increasing angles of attack. These simulations follow the HLPW4 guidelines for Case 2, which focuses on the $C_{Lmax}$ study of the aircraft, using nominal slat and flap deflections. The numerical setup uses a half span of 29.4 meters and includes all the elements that compose the CRM-HL, as complex as they may be. In fact, unlike the fifth edition, the fourth edition of the workshop didn't propose a technological effect study by gradually adding different components to isolate their impact on the global physics, such as the slats, the flaps or the nacelle. A schematic representation of the full configuration is provided in Figure 14. The inboard/outboard deflections are set to 30°/30° for the slats, and 40°/37° for the flaps. The reference static temperature and pressure are set to 289 K and 1.70093 bar, respectively, while the freestream velocity corresponds to a subsonic Mach number of 0.2. To achieve the target Reynolds number of $5.49 \times 10^6$, based on the mean aerodynamic chord (MAC = 7.005 m), the reference viscosity is also artificially adjusted to a non-physical value, although participants were also given the option to scale the geometry instead.

Two reference body-fitted wall-resolved solutions from the RANS TFG of the HLPW4 are used to benchmark the present IBM results. The first set of solutions originates from the LAVA CFD solver [72], developed at NASA Ames, and was obtained using multi-block structured grids coupled with overset techniques. The second set comes from the FUN3D CFD solver [82], developed at NASA Langley, and was computed on unstructured fixed grids provided directly by the workshop. Both datasets were presented at the 2022 AIAA AVIATION Forum and are publicly available on the HLPW website [65]. As previously mentioned, other Cartesian methods could have been used as additional references, but none were submitted in a RANS context. To ensure that the analysis remains as representative as possible of the current capabilities of both structured and unstructured solvers in a RANS-SA context, we selected each participant's best-practice numerical solutions as references: the C-level overset mesh for LAVA and the 2.2.D unstructured mesh for FUN3D, generated with Fidelity Pointwise using a combination of prism and tetrahedral elements. Comparing the IBM results with the two reference solutions, based on different solvers and mesh topologies, also provides valuable insight into flow prediction at high angles of attack across RANS approaches, where discrepancies among the RANS TFG become more pronounced [73]. Furthermore, both LAVA and FUN3D solutions rely on second-order finite volume methods using the classical Roe scheme, combined with either the standard SA or SA-neg turbulence models, which closely matches the numerical framework employed by FastS in the present study.

<figure><img src="image_20.png" /><figcaption>Fig. 14 Half-plane case description of the NASA CRM-HL for the HLPW4.</figcaption></figure>

Building on observations from simpler 2D cases, as well as preliminary results obtained on the CRM-HL in 2023 [81] and 2024 [49], the following discussion focuses solely on our best-practice IBM solutions, compared against both LAVA and FUN3D references. Within the F42 adapted methodology, a $y_{\tau target}^+$ of 500 was selected as the best practical compromise, corresponding to a minimum near-wall resolution of 3.5 mm. Since applying the finest near-wall resolution across the entire aircraft would unnecessarily inflate the final cell count, the minimum near-wall spacing $h_{min}$ is imposed only in the vicinity of the slats, slat brackets, flaps, and nacelle chine, where the modeling height is most
likely to be locally adapted (cf. Figure 15). Based on prior experience with this application, we also refine the leading and trailing edges of the wing, the wing tip, and the forward section of the nacelle to ensure proper IBM point projection and adequate resolution in zones where separation may occur. The rest of the wing, the majority of the nacelle, and all three flat fairings are meshed at twice this minimal spacing, while the fuselage is intentionally left underresolved. The surface mesh of the entire geometry has been generated in Fidelity Pointwise directly from the CAD definition, respecting both local curvature and the prescribed near-wall spacing.

<figure><img src="image_21.png" /><figcaption>Fig. 15 NASA CRM-HL (Ma 0.2 - Re 5.49M - AoA 7.05°/ 17.05°): near-wall resolutions.</figcaption></figure>

As discussed in Section III.A, the computational domain is defined around the half-plane using the symmetry plane location at $y = 0$, and extending the external boundaries in all three spatial directions by 100 times the half-span. Leveraging the modified octree generation technique with automatic creation of refinement boxes around the immersed geometry, as described in section III.B, the final number of cells slightly exceeds $2 \times 10^9$, as reported in Table 2. Although this number of cells appears large compared to the 3.C structured and 2.2.D unstructured meshes, it is important to emphasize that only half of these cells are actual interior, computed cells. It should also be noted that the resulting IBM4 mesh was automatically generated in parallel within 5 to 10 minutes, whereas the body-fitted meshes required up to several weeks of manual preparation. In addition, the isotropic nature of Cartesian cells inherently increases the number of cells in both spanwise and chordwise directions, dimensions where RANS solvers typically allow coarser resolution. To a lesser extent, these cell count differences between Cartesian and curvilinear meshes can also be observed in wall-modeled LES simulations, for which the close-wall cells approach isotropy in both cases [51].

The entire IBM pre-processing workflow, including full mesh generation, flow field initialization using freestream conditions, and interpolation data search for both IBM image points and inter-grid ghost cells, took approximately one hour to complete in parallel. The main challenge posed by this geometry lies in its numerous fine geometric details, with many sharp corners and concave regions. The most critical areas are near the flap and slat attachments, where the limited spacing offered by this CRM-HL configuration makes it difficult to properly project target points onto the wall or image fronts. As a result, some IBM target points ended up with mismatched projection normals, forming trios of unaligned IBM points. Although a smaller modeling height could reduce the number of these misprojected points, the issue would largely persist without additional mesh refinement. In fact, the small geometry gaps also caused issues for both unstructured and structured body-fit meshes. Fortunately, in our case, the occurrence of such unaligned IBM points remains highly localized, and the vast majority were successfully projected. For stability reasons, however, a minmod slope limiter was applied to each Cartesian zone that contained one or more unaligned IBM points.

Each simulation on the IBM4 Cartesian mesh was performed on 10 Intel Sapphire Rapids nodes on ONERA's SATOR computer cluster, distributed across 40 × 24 CPUS to operate in a hybrid MPI-openMP configuration. Convergence was typically achieved over four runs of 20,000 iterations, including the modeling height adaptation applied after the first two runs. As with the 2DMEA, each IBM solution is analyzed only after both the modeling height adaptation and the application of second-order wall pressure reconstruction. The benefits of these two steps have already been demonstrated in [49], allowing us to reduce the number of solutions plotted in each figure for improved clarity and conciseness. The flow solution is then examined in detail, both at the geometry surface using pressure and friction coefficients, and within the volume mesh using constant-x slices and constant-z velocity profiles.
<table><caption>Table 2 NASA CRM-HL (Ma 0.2 - Re 5.49M - AoA 7.05°/ 17.05°): mesh characteristics.</caption><thead><tr><th></th><th>IBM4</th><th>BF 3.C</th><th>BF 2.2.D</th></tr></thead><tbody><tr><td>Grid topology</td><td>Cartesian</td><td>Structured with Overset</td><td>Unstructured (Prism-Tet)</td></tr><tr><td>Minimum spacing h<sub>min</sub>/c</td><td>3.5e-3</td><td>N/A</td><td>N/A</td></tr><tr><td>Corresponding y<sup>+</sup><sub>target</sub></td><td>400</td><td>1</td><td>0.75</td></tr><tr><td>Number of points</td><td>2194M</td><td>223M</td><td>200M</td></tr><tr><td>% of cell types (0 / 1 / 2)</td><td>26% / 47% / 27%</td><td>-</td><td>-</td></tr></tbody></table>

The evolution of the aerodynamic coefficients with the angle of attack is presented in Figure 16. All three CFD simulations are compared with the experimental results acquired in the QinetiQ 5-meter pressurized low-speed wind tunnel at Farnborough (UK), which provides full polar data from -4° to 21.5° and served as the primary experimental reference for the HLPW4 [83]. The IBM results, shown in green, compare well with both experimental and numerical references for the lift and drag coefficients. Agreement is particularly good at lower angles of attack (2.78° and 7.05°), with lift coefficient errors kept below 3%. The discrepancy increases at AoA = 11.29°, and continues to grow slightly to 6.8% at AoA = 17.05°, as the flow topology predicted by the IBM solutions begins to deviate from the reference trends. That said, the error made on the lift coefficient remains constrained and consistent with the deviations observed on the 2DMEA at a similar $y^+_{target}$ value. Based on this observation, finer Cartesian meshes combined with smaller modeling heights should reliably reduce lift coefficient errors.

In terms of drag, the IBM solution remains close to the experimental data even at higher angles of attack, likely benefiting from compensating effects around the wing in view of the growing discrepancies in lift and especially in pitching moment. As observed with the previous multi-element profile, the largest sources of error lie far from the center of integration, primarily near the inboard and outboard separation regions, and around the flaps. This explains why the pitching moment coefficient shows the largest scatter among the fixed-grid RANS TFG participants, and why our own IBM solutions have difficulty converging toward the experimental reference. The two body-fitted wall-resolved solutions also diverge in their pitching moment predictions, although both remain close to the experimental data. While FUN3D tends to overpredict the moment, LAVA slightly underpredicts it, further highlighting the sensitivity of this quantity when the geometry faces strong flow incidences. The following detailed analysis of the results obtained at AoA = 7.05° and AoA = 17.05° further clarifies the differences between the three CFD simulations, both at the wall and within the flow region, that ultimately led to these discrepancies in the integral quantities.

<figure><img src="image_22.png" /><figcaption>Fig. 16 NASA CRM-HL (Ma 0.2 - Re 5.49M - AoA 7.05°/ 17.05°): aerodynamic coefficients (lift, drag, and pitching moment). Errors are shown for each CFD solution (LAVA, FUN3D, and FastS) with respect to the experimental results. Note that the experimental data have been corrected for wind tunnel effects.</figcaption></figure>
1. Results at AoA 7.05°

<figure><img src="image_23.png" /><figcaption>**Fig. 17 NASA CRM-HL (Ma 0.2 - Re 5.49M - AoA 7.05°): skin friction distributions over the wing surface.**</figcaption></figure>

At AoA = 7.05°, the IBM solution shows excellent agreement with the LAVA and FUN3D body-fitted results, which are nearly superimposed at this low incidence. Qualitatively, the friction coefficient distributions on the wing suction side (Figure 17) and the off-body vorticity contours in the flow field (Figure 20) both demonstrate very good prediction of the overall turbulent flow characteristics from the IBM solution compared to the body-fitted references, even in complex regions such as the wing tip and nacelle. The main differences occur in the vicinity of the flaps, where the IBM solution underpredicts the amount of detached flow compared to the body-fitted references. Another smaller deviation occurs in the immediate continuity of the nacelle, likely due to insufficient resolution to properly capture and propagate the turbulent wakes generated by the two concave corner flows emanating from the slats. Elsewhere, the major flow structures are accurately reproduced by the IBM simulation, as can be seen in each of the constant-x slices shown in Figure 20, which demonstrate excellent agreement across all three simulations. This is consistent with initial expectations, since unidirectional equilibrium wall models such as the one applied in these simulations have been specifically designed and validated for this type of aerodynamic application, however geometrically complex it may be, with predominantly attached boundary layers and negligible history effects. The validity range of the algebraic wall model can then extend up to the logarithmic boundary layer limit, which increases with the Reynolds number, allowing higher $y_{target}^+$ values and coarser mesh resolutions.

These initial observations are fully supported by the quantitative analysis of the pressure and friction coefficient distributions shown in sections A to H (excluding section C, located downstream of the nacelle), extracted along the wingspan and presented in Figures 18 and 19. The near-perfect superposition of the IBM pressure curves with those of LAVA and FUN3D further demonstrates the predictive capabilities of the approach at this angle of attack and justifies the good agreement observed in the aerodynamic coefficients in Figure 16. The friction coefficient distributions also exhibit strong overall agreement, despite a mild but systematic underestimation of the leading edge friction peaks on each element compared to LAVA. Additional differences appear in sections D, E, and F, where the IBM solution maintains mostly attached flow over the flaps, in contrast to the flow separations predicted by the body-fitted references at these locations. Further validation of the IBM solution can be observed deeper into the boundary layer through the x-velocity profiles shown in Figure 21, sampled from both the inboard (sections A*/B*) and outboard (sections C*/D*) regions of the wing suction side at constant-z abscissas. Across these locations, the IBM solution provides an accurate prediction of axial velocities that closely matches LAVA, FUN3D, or both, depending on the probe location.

These first results reinforce earlier conclusions regarding the ability of our IBM approach to deliver highly satisfactory approximations to reference wall-resolved, body-fitted RANS approaches for fully turbulent, bounded flows at modest incidence. In particular, the relatively coarse mesh resolution applied around the CRM-HL, along with the correspondingly large $y_{target}^+$ value of 500, proves sufficient for this configuration, offering an excellent balance between accuracy, stability, and computational efficiency. The next case, however, involves a less favorable scenario, in which the limitations previously observed on the 2DMEA become locally exacerbated.
<figure><img src="image_24.png" /><figcaption>Fig. 18 NASA CRM-HL (Ma 0.2 - Re 5.49M - AoA 7.05°): skin pressure distributions.</figcaption></figure>
<figure><img src="image_25.png" /><figcaption>Fig. 19 NASA CRM-HL (Ma 0.2 - Re 5.49M - AoA 7.05°): skin friction distributions.</figcaption></figure>
<figure><img src="image_26.png" /><figcaption>Fig. 20 NASA CRM-HL (Ma 0.2 - Re 5.49M - AoA 7.05°): off-body vorticity contours.</figcaption></figure>
<figure><img src="image_27.png" /><figcaption>Fig. 21 NASA CRM-HL (Ma 0.2 - Re 5.49M - AoA 7.05°): x-velocity profiles.</figcaption></figure>
## 2. Results at AoA 17.05°

<figure><img src="image_28.png" /><figcaption>**Fig. 22 NASA CRM-HL (Ma 0.2 - Re 5.49M - AoA 17.05°): skin friction distributions over the wing surface.**</figcaption></figure>

At AoA = 17.05°, the differences between the IBM and reference solutions are slightly more pronounced, as the flow naturally grows more unsteady and tends to separate across broader regions of the CRM-HL. Both qualitatively and quantitatively, the IBM results continue to provide strong flow prediction over most of the configuration, with some localized discrepancies compared to the body-fitted references. In Figure 27, all three solutions show excellent agreement across the flow field for sections I through IV, particularly on the wing section side, where all solvers consistently capture the major vortex structures, including all the swan-shaped vortices occurring in sections III and IV. Closer to the fuselage, the IBM solution begins to diverge from the body-fitted references, where the near-wall resolution appears too coarse to fully resolve the complex flow structures developing in these regions. The discrepancies extend into the immediate vicinity of the wing root, and may partly explain the inboard flow separation differences observed in Figure 22. Finally, the vorticity magnitudes plotted in section V indicate two distinct flow behaviors on the wing, whether the flow remains attached to the wall (LAVA) or separates significantly (FastS and FUN3D).

The skin friction distributions over the wing surface provide further insight into the discrepancies observed among the three simulations in the flow field. While the surface flow is largely consistent across the slats, flaps, and most of the wing, the intensity of the slat bracket wakes varies between simulations, ultimately leading to different flow behaviors near the wing tip. Consistent with experimental results, LAVA predicts a limited flow separation near the trailing edge, whereas FUN3D predicts a broader separated region forming a triangular pattern, likely induced by a specific slat bracket. A similar separation pattern is observed in the IBM results and among several other RANS TFG participants. Within the HLPW4 context, this unphysical flow separation has been known to occur sporadically, influenced by mesh resolution, flow initialization, or turbulence model choices. In the IBM simulation, misprojected target points may be another cause of this behavior. As previously mentioned, projection issues are more pronounced at the ends of the brackets in contact with the wing, where the local use of a slope limiter seems necessary to prevent the simulation from diverging. However, more validation tests are required to precisely identify the cause or combination of causes of these flow separations in our simulations.

Additionally, the CRM-HL geometries used during HLPW5 proposed new, simpler attachment systems that offer better spacing and seem more suitable for CFD simulations than those used during HLPW4, regardless of the meshing technique employed. To illustrate this aspect, sectional views of the same slat bracket are proposed for three different HLPW test cases in Figure 23, while the same is done for a flat fairing in Figure 24. Most of the geometric features are measured to accurately show the major differences between each configuration. Critical regions can be observed toward the wing, with a steep concave angle of 44°, and on the bottom part of the wing slat, where the brackets reattach. For the current geometry, it is evident that the spacing allocated to some brackets is insufficient to properly locate IBM target points and their corresponding image points, which is necessary for the proper application of the wall function. Conversely, the fifth edition appears to partially alleviate this problem, by providing more room in this area. These geometric changes to the wing attachments could improve the projection step in the pre-processing stage and enhance the outboard flow prediction in the IBM simulation. That said, these simplified geometries have also been associated with an increased incidence of artificial wedge separation patterns in RANS solutions during HLPW5, which further complicates the flow analysis in this region.
<figure><img src="image_29.png" /><figcaption>Fig. 23 NASA CRM-HL (Ma 0.2 - Re 5.49M - AoA 7.05°/ 17.05°): Comparison of slat bracket geometries between three CRM-HL configurations studied at HLPW4 and HLPW5. Cross sections of the 8th slat bracket from the wing tip. Lengths are measured in millimetres.</figcaption></figure>

<figure><img src="image_30.png" /><figcaption>Fig. 24 NASA CRM-HL (Ma 0.2 - Re 5.49M - AoA 7.05°/ 17.05°): Comparison of flap track fairing geometries between three CRM-HL configurations studied at HLPW4 and HLPW5. Cross sections of the 2nd flap track fairing from the wing tip. Lengths are measured in millimetres.</figcaption></figure>

On the other hand, inboard flow separation may be as critical as outboard flow separation at this angle of attack. Near the fuselage, the IBM simulation predicts a different corner flow wake than the body-fitted solutions. While the IBM wake follows an oscillating pattern along the wing root, the body-fitted solutions produce a more stable wake following a 10° angled path, which later interacts with the first two bracket-induced wakes. Given the earlier observations in the flow field, these inboard differences may result from insufficient resolution near the fuselage, leading to incorrect upstream flow predictions. Finally, all three simulations diverge around the nacelle, showing different levels of separated flow around its pylon. Further downstream, the IBM solution eventually follows the LAVA trend, although the chordwise change in resolution appears to diffuse the turbulent wake within the IBM Cartesian mesh.

Figures 25, 26, and 28 provide a more quantitative perspective on the observations discussed above. In Figure 25, sections A, B, and C, show excellent agreement across all solutions, including the accurate capture of all three leading edge pressure peaks. Further along the wing span, small deviations emerge in sections D and E, near the flaps, where the IBM solution continues to underpredict the flow separation in these regions. Additional discrepancies are generated on the wing suction side, where the IBM pressure trend seems more damped than the references, although in the same proportion to what was observed on the 2DMEA. Towards the wing tip, the IBM solution aligns more closely with FUN3D, as both deviate from LAVA by predicting significant flow separations near sections G and H. Overall, these skin pressure coefficient profiles support the previous observations and rationalize the combination of favorable and less favorable trends in the global aerodynamic coefficients.
<figure><img src="image_31.png" /><figcaption>Fig. 25 NASA CRM-HL (Ma 0.2 - Re 5.49M - AoA 17.05°): skin pressure distributions.</figcaption></figure>
<figure><img src="image_32.png" /><figcaption>Fig. 26 NASA CRM-HL (Ma 0.2 - Re 5.49M - AoA 17.05°): skin friction distributions.</figcaption></figure>
<figure><img src="image_33.png" /><figcaption>Fig. 27 NASA CRM-HL (Ma 0.2 - Re 5.49M - AoA 17.05°): off-body vorticity contours.</figcaption></figure>
<figure><img src="image_34.png" /><figcaption>Fig. 28 NASA CRM-HL (Ma 0.2 - Re 5.49M - AoA 17.05°): x-velocity profiles.</figcaption></figure>
In the absence of skin friction data from FUN3D on the workshop website, the IBM solution is compared only with LAVA in Figure 26. Excellent agreement is observed between the two for sections D, E, and F, where the flow remains mostly attached to the wall. The flow behavior near the slats is particularly well predicted by the IBM solution, and the overall trends in these mid-span sections are consistent with earlier predictions from the 2DMEA test case. In contrast, in sections A, B, G, and H, the discrepancies found in the inboard and outboard flow separations cause the IBM solution to diverge more significantly from the body-fitted reference on the main element.

Lastly, in Figure28, the IBM solution is compared with both LAVA and FUN3D results using velocity profiles extracted at constant abscissas on the wing suction side. As with the previous angle of attack, excellent agreement is observed between all CFD solutions for the A* and B* probes. While the IBM solution moderately overestimates the boundary layer thickness in section A1, this trend eventually corrects itself further downstream. Consistent with the earlier skin friction analysis, the IBM results follow more closely the LAVA solution at the B* probe loactions, which lie downstream of the nacelle. The opposite behavior is seen at the C* probes, where IBM solution now produces velocity profiles similar to FUN3D, at the center of the triangular separation region. However, at the tip of the wing, the IBM profiles fail to capture the same flow behavior predicted by both body-fitted, wall resolved solutions, despite maintaining a constant near-wall resolution to account for the complexity of the flow topology in this region.

V. Conclusion

This work is a direct continuation of the recent developments brought to the immersed boundary method (IBM) developed by Peron et al. [53] and Constant et al. [49, 54]. The main contribution of this work is a more in-depth and extensive validation of our recent methodology for RANS simulations on the complex NASA High Lift Common Research Model (CRM-HL), an application that still represents a significant challenge for most CFD solvers. By modifying the IBM pre- and post-processing to simplify case setup, reduce cell count, and allow easy tracking of aerodynamic coefficients, multiple simulations on large Cartesian meshes were enabled at low computational cost. Two angles of attack were studied to evaluate the performance of our IBM in both favorable and unfavorable situations, providing numerous insights that will undoubtedly inform future development. Following the numerical setup provided by the fourth edition of the international AIAA CFD High Lift Prediction Workshop (HLPW4), turbulent flow prediction at the wall and across the flow field was thoroughly investigated at different probe locations. Overall, good to excellent agreement was found between the IBM solutions and the two body-fitted wall-resolved NASA solutions in both structured and unstructured contexts. These positive results were obtained without any numerical perturbations in the vicinity of the aircraft or spurious oscillations at the wall, further validating the robustness and accuracy of our Cartesian methodology when coupled with a geometric regularization of the forcing point location around the geometry.

Despite these promising observations, the current IBM approach showed some limitations at higher angles of attack, particularly in predicting inboard and outboard flow separations. To stabilize the simulations under these conditions, a slope limiter was finally applied to any Cartesian block containing one or more misprojected target points so that the simulation didn't diverge due to the creation of unphysical and local density sinks or sources. As mentioned above, this primarily affected the wing attachments, namely the slat brackets and the flap fairings. It is not known whether not relying on the sporadic use of a limiter near these elements would have resulted in increased fidelity with both references, or different flow patterns in detached regions. If a complete redesign of the IBM point projection algorithm around the immersed geometries seems overly complicated at this time, simpler CAD definitions with larger geometry gaps and wider concave regions could greatly limit the this problem and reduce the use of slope limiter throughout the domain. Studying the CRM-HL using HLPW5 test cases with increasing geometric complexity could either confirm or refute these assumptions and guide future work.

In addition, while RANS approaches can have difficulty approaching the stall phenomenon, especially with the standard Splart-Allmaras model, it appears that the use of an equilibrium wall model tends to exacerbate these limitations at lower angles of attack. To compensate for this, a lower modeling height should theoretically be applied at 17.05°, in addition to a more resolved Cartesian mesh to better account for the complex reality of the flow. The study of the 2D multi-element airfoil actually helped to quantify the error in the aerodynamic coefficient predictions with respect to the near-wall resolution within our IBM approach at high angles of attack. Unfortunately, only the two finest meshes were able to reduce the prediction error below the accepted error threshold for this type of aerodynamic simulation. Since these resolutions cannot realistically be applied to the wall without significantly increasing the mesh size and computational time in 3D, other solutions will have to be envisioned for industrial-level cases to further improve the accuracy of our methodology under unfavorable flow conditions. These solutions could be based on improved wall modeling [84–86] or alternative IBM approaches allowing smaller modeling heights at coarser mesh resolutions [87, 88].
References

[1] Belcastro, C., and Foster, J., “Aircraft Loss-of-Control Accident Analysis,” *AIAA Guidance, Navigation, and Control Conference*, 2010. https://doi.org/10.2514/6.2010-8004.

[2] Abu Salem, K., Palaia, G., Chiarelli, M., and Bianchi, M., “A Simulation Framework for Aircraft Take-Off Considering Ground Effect Aerodynamics in Conceptual Design,” *Aerospace*, 2023. https://doi.org/10.3390/aerospace10050459.

[3] Liu, X., Abà, A., Capone, P., Manfriani, L., and Fu, Y., “Atmospheric Disturbance Modelling for a Piloted Flight Simulation Study of Airplane Safety Envelope over Complex Terrain,” *Aerospace*, 2022. https://doi.org/10.3390/aerospace9020103.

[4] Lee, D., Fahey, D., Skowron, A., Allen, M., Burkhardt, U., Chen, Q., Doherty, S., Freeman, S., Forster, P., Fuglestvedt, J., Gettelman, A., De León, R., Lim, L., Lund, M., Millar, R., Owen, B., Penner, J., Pitari, G., Prather, M., Sausen, R., and Wilcox, L., “The contribution of global aviation to anthropogenic climate forcing for 2000 to 2018,” *Atmospheric Environment*, 2021. https://doi.org/10.1016/j.atmosenv.2020.117834.

[5] Krein, A., and Williams, G., “Flightpath 2050: Europe’s Vision for Aeronautics,” *Innovation for Sustainable Aviation in a Global Environment*, 2012. https://doi.org/10.3233/978-1-61499-063-5-63.

[6] Yang, L., Hu, Y.-J., Wang, H., Li, C., Tang, B.-J., Wang, B., and Cui, H., “Uncertainty quantification of CO2 emissions from China’s civil aviation industry to 2050,” *Journal of Environmental Management*, 2023. https://doi.org/10.1016/j.jenvman.2023.117624.

[7] Kito, M., Takayabu, H., and Nansai, K., “Carbon-neutral pathways to 2050 for Japan’s aviation industry in the absence of a mass supply of sustainable aviation fuels,” *Journal of Industrial Ecology*, 2023. https://doi.org/10.1111/jiec.13443.

[8] Graver, B., Zhang, K., and Rutherford, D., “CO2 emissions from commercial aviation, 2018,” *International Council on Clean Transportation*, 2019. URL https://theicct.org/sites/default/files/publications/ICCT_CO2-commercl-aviation-2018_20190918.pdf.

[9] Cambier, Laurent, Heib, Sébastien, and Plot, Sylvie, “The Onera elsA CFD software: input from research and feedback from industry,” *Mechanics & Industry*, 2013. https://doi.org/10.1051/meca/2013056.

[10] Kroll, N., Langer, S., and Schwöppe, A., “The DLR Flow Solver TAU - Status and Recent Algorithmic Developments,” *52nd Aerospace Sciences Meeting*, 2014. https://doi.org/10.2514/6.2014-0080.

[11] Kiris, C., Barad, M., Housman, J., Sozer, E., Brehm, C., and Moini-Yekta, S., “The LAVA Computational Fluid Dynamics Solver,” *AIAA Scitech 2014 Forum*, 2014. https://doi.org/10.2514/6.2014-0070.

[12] Spalart, P. R., and Venkatakrishnan, V., “On the role and challenges of CFD in the aerospace industry,” *The Aeronautical Journal*, 2016. https://doi.org/10.1017/aer.2015.10.

[13] Slotnick, J., Khodadoust, A., Alonso, J., Darmofal, D., Gropp, W., Lurie, E., and Mavriplis, D., “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.

[14] Rumsey, C., “High-Lift Prediction Workshops: Retrospective, Lessons Learned, and Future Prospects,” *34th Congress of the International Council of the Aeronautical Sciences*, 2024. URL https://ntrs.nasa.gov/api/citations/20240006238/downloads/ICAS-rumsey-hilift_4.pdf.

[15] Chien, K.-Y., “Predictions of Channel and Boundary-Layer Flows with a Low-Reynolds-Number Turbulence Model,” *AIAA Journal*, 1982. https://doi.org/10.2514/3.51043.

[16] Spalart, P. R., and Allmaras, S. R., “A one-equation turbulence model for aerodynamic flows,” *30th Aerospace Sciences Meeting and Exhibit*, 1992. https://doi.org/10.2514/6.1992-439.

[17] Wilcox, D. C., “Formulation of the $k-\omega$ Turbulence Model Revisited,” *AIAA Journal*, 2008. https://doi.org/10.2514/1.36541.

[18] Kiris, C. C., Ghate, A. S., Duensing, J. C., Browne, O. M., Housman, J. A., Stich, G.-D., Kenway, G., Fernandes, L. S., and Machado, L. M., “High-Lift Common Research Model: RANS, HRLES, and WMLES perspectives for CLmax prediction using LAVA,” *AIAA Scitech 2022 Forum*, 2022. https://doi.org/10.2514/6.2022-1554.

[19] I. Celik, M. F., M. Klein, and Janicka, J., “Assessment measures for URANS/DES/LES: an overview with applications,” *Journal of Turbulence*, 2006. https://doi.org/10.1080/14685240600794379.
[20] Yang, X. I. A., and Griffin, K. P., “Grid-point and time-step requirements for direct numerical simulation and large-eddy simulation,” *Physics of Fluids*, 2021. https://doi.org/10.1063/5.0036515.

[21] Sagaut, P., “Large Eddy Simulation for Incompressible Flows. An Introduction,” *Measurement Science & Technology*, 2001. https://doi.org/10.1088/0957-0233/12/10/70.

[22] Garnier, E., Adams, N., and Sagaut, P., “Large Eddy Simulation for Compressible Flows,” *Springer*, 2009. https://doi.org/10.1007/978-90-481-2819-8.

[23] Labourasse, E., and Sagaut, P., “Reconstruction of Turbulent Fluctuations Using a Hybrid RANS/LES Approach,” *Journal of Computational Physics*, 2002. https://doi.org/10.1006/jcph.2002.7169.

[24] Spalart, P. R., “Chapter Four - Hybrid RANS-LES Methods,” *Advanced Approaches in Turbulence*, 2021. https://doi.org/10.1016/B978-0-12-820774-1.00010-0.

[25] Pope, S. B., “Turbulent Flows,” *Cambridge University Press*, 2000.

[26] Spalart, P., Deck, S., Shur, M., Squires, K. D., Strelets, M. K., and Travin, A., “A New Version of Detached-eddy Simulation, Resistant to Ambiguous Grid Densitie,” *Theoretical and Computational Fluid Dynamics*, 2006. https://doi.org/10.1007/s00162-006-0015-0.

[27] Spalart, P. R., “Detached-Eddy Simulation,” *Annual Review of Fluid Mechanics*, 2009. https://doi.org/10.1146/annurev.fluid.010908.165130.

[28] Deck, S., “Recent improvements in the Zonal Detached Eddy Simulation (ZDES) formulation,” *Theoretical and Computational Fluid Dynamics*, 2012. https://doi.org/10.1007/s00162-011-0240-z.

[29] Deck, S., and Renard, N., “Towards an enhanced protection of attached boundary layers in hybrid RANS/LES methods,” *Journal of Computational Physics*, 2020. https://doi.org/10.1016/j.jcp.2019.108970.

[30] Kawai, S., and Larsson, J., “Wall-modeling in large eddy simulation: Length scales, grid resolution, and accuracy,” *Physics of Fluids*, 2012. https://doi.org/10.1063/1.3678331.

[31] Larsson, J., Kawai, S., Bodart, J., and Bermejo-Moreno, I., “Large Eddy Simulation with Modeled Wall-Stress: Recent Progress and Future Directions,” *Mechanical Engineering Reviews*, 2016. https://doi.org/10.1299/mer.15-00418.

[32] Marié, S., Ricot, D., and Sagaut, P., “Comparison between lattice Boltzmann method and Navier-Stokes high order schemes for computational aeroacoustics,” *Journal of Computational Physics*, 2009. https://doi.org/10.1016/j.jcp.2008.10.021.

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

[34] Wilhelm, S., Jacob, J., and Sagaut, P., “An explicit power-law-based wall model for lattice Boltzmann method-Reynolds-averaged numerical simulations of the flow around airfoils,” *Physics of Fluids*, 2018. https://doi.org/10.1063/1.5031764.

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

[36] Degrigny, J., Cai, S. G., Boussuge, J. F., and Sagaut, P., “Improved wall model treatment for aerodynamic flows in LBM,” *Computers and Fluids*, 2021. https://doi.org/10.1016/j.compfluid.2021.105041.

[37] Husson, J., Terracol, M., Deck, S., and Le Garrec, T., “Critical assessment of wall model numerical implementation in LBM,” *Computers and Fluids*, 2023. https://doi.org/10.1016/j.compfluid.2023.105857.

[38] Suss, A., Mary, I., Le Garrec, T., and Marié, S., “A hybrid lattice Boltzmann - Navier-Stokes method for unsteady aerodynamic and aeroacoustic computations,” *Journal of Computational Physics*, 2023. https://doi.org/10.1016/j.jcp.2023.11209.

[39] Trapani, G., Brionnaud, R., and Holman, D., “XFlow Contribution to the Third High-Lift Prediction Workshop,” *2018 Applied Aerodynamics Conference*, 2018. https://doi.org/10.2514/6.2018-2847.

[40] Konig, B., Fares, E., Murayama, M., and Ito, Y., “PowerFLOW Simulations for the Third AIAA High-Lift Prediction Workshop,” *2018 AIAA Aerospace Sciences Meeting*, 2018. https://doi.org/10.2514/6.2018-1255.
[41] Konig, B., Duda, B. M., and Laskowski, G. M., "Lattice Boltzmann Simulations for the 4th AIAA High-Lift Prediction Workshop using PowerFLOW," *AIAA AVIATION 2022 Forum*, 2022. https://doi.org/10.2514/6.2022-3433.

[42] Kiris, C. C., Ghate, A. S., Browne, O. M. F., Slotnick, J., and Larsson, J., "HLPW-4: Wall-Modeled Large-Eddy Simulation and Lattice-Boltzmann Technology Focus Group Workshop Summary," *Journal of Aircraft*, 2023. https://doi.org/10.2514/1.C037193.

[43] Mittal, R., and Iaccarino, G., "Immersed Boundary Methods," *Annual Review of Fluid Mechanics*, 2005. https://doi.org/10.1146/annurev.fluid.37.061903.175743.

[44] Verzicco, R., "Immersed Boundary Methods: Historical Perspective and Future Outlook," *Annual Review of Fluid Mechanics*, 2023. https://doi.org/10.1146/annurev-fluid-120720-022129.

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

[46] Mohd-Yusof, J., "Combined immersed-boundary/B-spline methods for simulations of flow in complex geometries," *Annual Research Briefs, Center for Turbulence Research*, 1997.

[47] Capizzano, F., "Turbulent Wall Model for Immersed Boundary Methods," *AIAA Journal*, 2011. https://doi.org/10.2514/1.J050466.

[48] Tamaki, Y., and Imamura, T., "Turbulent Flow Simulations of the Common Research Model Using Immersed Boundary Method," *AIAA Journal*, 2018. https://doi.org/10.2514/1.J056654.

[49] Constant, B., Péron, S., Beaugendre, H., and Benoit, C., "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. https://doi.org/10.1016/j.jcp.2024.113441.

[50] Degrigny, J., "Towards the Computational Prediction of Low-Speed Buffet: Improved Wall Modeling for the Lattice Boltzmann Method," *PhD thesis (Université d'Aix-Marseille)*, 2022.

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

[52] Péron, S., and Benoit, C., "Automatic off-body overset adaptive Cartesian mesh method based on an octree approach," *Journal of Computational Physics*, 2013. https://doi.org/10.1016/j.jcp.2012.07.029.

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

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

[55] ONERA, "Fast - CFD Solver," https://onera.github.io/Fast/, Last accessed: June 2025.

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

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

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

[59] Thomas, J. L., and Walters, R. W., "Upwind relaxation algorithms for the Navier-Stokes equations," *AIAA Journal*, 1987. https://doi.org/10.2514/3.9659.

[60] Yoon, S., and Jameson, A., "A multigrid LU-SSOR scheme for approximate Newton iteration applied to the Euler equations," *NASA Contractor Report*, 1986.

[61] Coakley, T. J., "Implicit upwind methods for the compressible Navier-Stokes equations," *AIAA Journal*, 1985. https://doi.org/10.2514/3.8923.
[62] Musker, A., "Explicit expression for the smooth wall velocity distribution in a turbulent boundary layer," *AIAA Journal*, 1979. https://doi.org/10.2514/3.61193.

[63] White, F. M., "Viscous Fluid Flow," *McGraw-Hill*, 1974.

[64] Benoit, C., Péron, S., and Landier, S., "Cassiopee: A CFD pre- and post-processing tool," *Aerospace Science and Technology*, 2015. https://doi.org/10.1016/j.ast.2015.05.023.

[65] AIAA, "AIAA CFD High Lift Prediction Workshop (HLPW)," https://hiliftpw.larc.nasa.gov, Last accessed: June 2025.

[66] Rumsey, C. L., Slotnick, J. P., Long, M., Stuever, R. A., and Wayman, T. R., "Summary of the First AIAA CFD High-Lift Prediction Workshop," *Journal of Aircraft*, 2011. https://doi.org/10.2514/1.C031447.

[67] Rumsey, C. L., and Slotnick, J. P., "Overview and Summary of the Second AIAA High-Lift Prediction Workshop," *Journal of Aircraft*, 2015. https://doi.org/10.2514/6.2014-0747.

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

[69] Rumsey, C. L., Slotnick, J. P., and Woeber, C. D., "Fourth High-Lift Prediction/Third Geometry and Mesh Generation Workshops: Overview and Summary," *Journal of Aircraft*, 2023. https://doi.org/10.2514/1.C037168.

[70] Lacy, D. S., and Sclafani, A. J., "Development of the High Lift Common Research Model (HL-CRM): A Representative High Lift Configuration for Transonic Transports," *54th AIAA Aerospace Sciences Meeting*, 2016.

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

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

[73] Gooch, C. F. O., and Coder, J. G., "Lessons Learned by the Fixed-Grid RANS TFG for HLPW-4/GMGW-3," *AIAA AVIATION 2022 Forum*, 2022. https://doi.org/10.2514/6.2022-3211.

[74] Galbraith, M. C., and Karman, S. L., "HLPW-4/GMGW-3: High Order Discretization Technology Focus Group Workshop Summary," *AIAA AVIATION 2022 Forum*, 2022. https://doi.org/10.2514/6.2022-3292.

[75] Park, M. A., "HLPW-4/GMGW-3: Mesh Adaptation for RANS Technology Focus Group Workshop Summary," *AIAA AVIATION 2022 Forum*, 2022. https://doi.org/10.2514/6.2022-3210.

[76] Ashton, N., Batten, P., Cary, A. W., Holst, K. R., and Skaperdas, V., "HLPW-4/GMGW-3: Hybrid RANS/LES Technology Focus Group Workshop Summary," *AIAA AVIATION 2022 Forum*, 2022. https://doi.org/10.2514/6.2022-3293.

[77] Kiris, C. C., Ghate, A. S., Browne, O. M., Slotnick, J. P., and Larsson, J., "HLPW-4/GMGW-3: Wall-Modeled LES and Lattice-Boltzmann Technology Focus Group Workshop Summary," *AIAA AVIATION 2022 Forum*, 2022. https://doi.org/10.2514/6.2022-3294.

[78] NASA, "Langley Research Center Turbulence Modeling Resource," https://turbmodels.larc.nasa.gov/, Last accessed: June 2025.

[79] Allmaras, S., Johnson, F., and Spalart, P., "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.

[80] Shur, M. L., Strelets, M. K., Travin, A. K., and Spalart, P. R., "Turbulence Modeling in Rotating and Curved Channels: Assessing the Spalart-Shur Correction," *AIAA Journal*, 2000. https://doi.org/10.2514/2.1058.

[81] Constant, B., "An improved immersed boundary method for the simulation of compressible turbulent flows around complex geometries," *PhD thesis (Université de Bordeaux)*, 2023. URL https://theses.hal.science/tel-04146884.

[82] Ahmad, N. N., Wang, L., Anderson, W. K., Balakumar, P., Iyer, P. S., and Nielsen, E. J., "FUN3D Simulations for the 4th AIAA High-Lift Prediction Workshop," *AIAA AVIATION 2022 Forum*, 2022. https://doi.org/10.2514/6.2022-3436.
[83] Lacy, D. S., and Clark, A. M., “Definition of Initial Landing and Takeoff Reference Configurations for the High Lift Common Research Model (CRM-HL),” *AIAA AVIATION 2020 Forum*, 2020. https://doi.org/10.2514/6.2020-2771.

[84] Berger, M. J., and Aftosmis, M. J., “An ODE-based Wall Model for Turbulent Flow Simulations,” *AIAA Journal*, 2017. https://doi.org/10.2514/6.2017-0528.

[85] Romanelli, M., Beneddine, S., Beaugendre, H., Sipp, D., Mary, I., and Bergmann, M., “Deep learning-based wall models for aerodynamic simulations: A new approach inspired by classical wall laws,” *in review*, 2022.

[86] A. L.-D., and HJ, B., “Machine learning building-block-flow wall model for large-eddy simulation,” *Journal of Fluid Mechanic*, 2023. https://doi.org/10.1017/jfm.2023.331.

[87] Tamaki, Y., Harada, M., and Imamura, T., “Near-wall modification of Spalart–Allmaras turbulence model for immersed boundary method,” *AIAA Journal*, 2017. https://doi.org/10.2514/1.J055824.

[88] Terracol, M., and Manueco, L., “A dynamic linearized wall model for turbulent flow simulation: Towards grid convergence in wall-modeled simulation,” *Journal of Computational Physic*, 2025. https://doi.org/10.1016/j.jcp.2024.11355.