Implementation of Immersed Boundaries through Volume Penalization in the Industrial Aeronautical Solver CODA
Jonatan Nunez David Huergo Suyash Shrestha Esteban Ferrer Eusebio Valero Diego Lodares Juan Guerra Juan Florenciano





Performant near-term quantum combinatorial optimization
    Phillip C. Lotshaw
    April 29, 2024
=======================================================







  In this work, we present the implementation and validation of the immersed boundary volume penalization methods in the CFD solver CODA. Our goal is the modelling and simulation of turbulent fluid flowsin complex 3D aerodynamic configuration through the numerical solution of the  Reynolds-Averaged Navier–Stokes equations using the Spalart–Allmaras turbulent model. To do that, an efficient immerse boundary method, together with an efficient preprocessing tool for the construction of unstructured hexahedral meshes with adaptive mesh refinement around immersed geometries has been developed.
  We report several numerical computations, including subsonic flow past a NACA0012 airfoil, subsonic flow past an MDA30P30N multi-element airfoil and subsonic flow around the NASA high-lift CRM aircraft. These simulations have been performed in the CODA solver with a second-order finite volume scheme as spatial discretization and an implicit backward-Euler scheme based on the matrix-free GMRES block-Jacobi iterative method. The reported numerical simulations are in very good agreement with their corresponding experimental data.







§ INTRODUCTION


In the ever-evolving landscape of computational fluid dynamics (CFD) and industrial simulations, the accurate representation of complex geometries and their interaction with fluid flow remains a critical challenge. The ability to model and analyze intricate structures in a computationally efficient manner is paramount, as it directly influences the design, performance, and optimization of various industrial systems. This challenge is particularly relevant in the context of the industrial solver CODA <cit.>, where the interaction between fluid flow and complex geometries plays a central role.CODA is a parallel software framework for multidisciplinary analysis and optimization of aircraft and helicopters based on advanced and accurate numerical methods <cit.>. The solver is being developed as part of a collaboration between the French National Aerospace Research Center (ONERA), the German Aerospace Center (DLR) and Airbus. CODA provides a robust, scalable and computational efficient integrated design process for aerodynamics and structural analysis. To perform an efficient analysis and optimization of aircraft on state-of-the-art HPC systems, the Navier–Stokes and the Reynolds-averaged Navier–Stokes (RANS) equations are solved for high Reynolds-number flow on unstructured grids with second-order finite-volume and higher-order discontinuous-Galerkin discretizations.The implementation of immersed boundaries through volume penalization has emerged as a promising technique to address the difficulty to generate body fitted meshes for complex geometries and simulate moving geometries. Immersed boundaries allow for the inclusion of detailed and intricate geometries within the computational domain, even in cases where traditional structured grids might be impractical or prohibitively expensive to generate. In the realm of CFD simulations, the term "immersed boundary" refers to the representation of complex structures, such as solid objects or porous media, as part of the fluid computational domain, without the need for a grid that conforms to their intricate shapes. This not only enhances the capability to model real-world industrial scenarios, but also significantly improves computational efficiency.
 The primary objective of this research is to provide an in-depth exploration of the methodology, mathematical formulations, and computational considerations involved in the implementation of immersed boundaries through volume penalization in the solver CODA. By doing so, this study aims to facilitate a deeper understanding of this technique and its potential applications in aeronautics. It is our hope that the insights and findings presented in this paper will contribute to the advancement of CFD simulations, offering engineers and researchers a valuable tool for the accurate modeling of complex industrial systems and the optimization of their performance.This paper is organized as follows: In <ref> we summarize the main numerical schemes implemented in CODA for simulating turbulent fluid flow with immersed boundary methods. In <ref> we discuss the preprocessing and postprocessing tools developed for generating refined meshes around immersed bodies and the analysis of simulation data. Next, in <ref> we present several numerical computations done with CODA with immersed boundary volume penalization methods: subsonic flow past a NACA0012 airfoil, subsonic flow past an MDA30P30N multi-element airfoil and subsonic flow around the NASA high-lift CRM aircraft. Finally, in <ref>, a summary of this work is presented.






§ NUMERICAL METHODS


CODA is a CFD solver, which  integrates an extensive set of algorithms for solving partial differential equations in a multi-physics context. In this work, we focus on the the implementation of an immersed boundary methodology for the modelling of turbulent fluid flows by the well-known Reynolds-Averaged Navier–Stokes equations. In our formulation, we have selected the Spalart–Allmaras turbulent model that we briefly reproduce in <ref> for completeness of the paper. Regarding the spatial and time discretization employed in CODA, we summarize the key ingredients in the next section.



 §.§ Spatial and Time Discretization


The spatial discretization schemes used in CODA to solve the Navier–Stokes equations and the Reynolds-Averaged Navier–Stokes equations are the finite volume method <cit.>, the modal discontinuous Galerkin method <cit.>, and the discontinuous Galerkin spectral element method <cit.>. The partial differential equations of interest can be written as a hyperbolic-parabolic system of conservation laws in differential form (see <ref> for more details)

    ut+∇·f(u,∇u)=s(u),

where u=u(x,t) is the vector of conserved quantities, and f=f(u,∇u) is the tensor of physical fluxes, which has two terms, the advective flux and the diffusive term, namely

    f(u,∇u)=
        f^A(u)-f^D(u,∇u),

and s=s(u) represents the source terms. The conservation laws can also be written in integral form as follows

    t∫_ΩuΩ +
      ∮_∂Ω(f^A(u)-f^D(u,∇u))·nσ = ∫_ΩsΩ.

Here, Ω is the spatial domain and n is the normal vector on the surface ∂Ω enclosing the spatial domain. In this work we are interested only in the finite volume schemes. A short description of these methods is outlined; for the discontinuous Galerkin schemes, see <cit.>. In the finite volume framework, the computational domain is divided into a set of non-overlapping polyhedral control volumes or cells, that is Ω=⋃_i^NΩ_i, and the integral form of the equations is discretized for each control volume Ω_i.
Defining the mean value of the function u in the cell Ω_i by

    u_i:=1|Ω_i|∫_Ω_iu(x,t)Ω,

we get the so-called semi-discrete formulation

    u_it:=-1|Ω_i|{∑_l[f^A·n̂]^*σ_l-∑_l[f^D·n̂]^*σ_l
        - s|Ω_i|
        }.

The surface integral has been approximated by a sum of the fluxes through the faces ∂Ω_l of the cell Ω_i. The σ_l stands for the area of the face ∂Ω_l, and |Ω_i| is the volume of the cell Ω_i.The finite volume discretization implemented in CODA is second-order accurate. The advective fluxes are approximated with the Roe schemes with an entropy fix. In this scheme, the left and right states at the cell interfaces are piecewise linearly reconstructed, and therefore gradient-based. More details on how are computed the element and face gradients can be found in <cit.>.
Regarding the time discretization, in CODA are available explicit, implicit and implicit-explicit integrators. In this work we have used an implicit backward-Euler scheme based on preconditioned matrix-free GMRES. 
Preconditioners include incomplete LU and line-inversion <cit.>.



 §.§ Immersed Boundary Methods


Immersed Boundary Methods are numerical techniques for simulating fluid-structure interaction problems involving complex geometries, providing an alternative to numerical methods based on body-fitted meshes. Immersed Boundary Methods were first proposed by <cit.> to simulate blood flow in the human heart. These methods provide a flexible framework for modeling fluid flow over complicated immersed bodies on grids that do not conform to the surface of the body by treating these as immersed boundaries within the computational domain. In this way, the Immersed Boundary Methods overcome the challenges of generating high-quality meshes around complex geometries, which can be time-consuming and computationally demanding, while accurately resolving complex flows using simple grids. The geometric components are discretized using a separate representation, such as a surface or volume mesh, which is then immersed within the computational grid. Immersed boundary conditions are applied to take into account the geometric components lying in the flow.Because the meshes used in the immersed boundary method are, in principle, non-conforming with the geometry of the obstacle, a very high resolution of this background mesh is required to obtain reliable numerical solutions. The employment of very fine meshes in the whole computational domain makes unpractical the use of the immersed boundary methods for industrial applications, like the study of the aerodynamics of aircrafts. In fact, very fine meshes are only necessary in specific places in the computational domain, like around the immersed geometry and the wake regions. Fine meshes far away from the immersed bodies are unnecessary. For these reasons, mesh refinement algorithms are used to refine the mesh only in the regions where it is actually essential for accurate computations.



  §.§.§ Immersed Boundary Volume Penalization Method


In this family of immersed boundary volume penalization methods, the flow equations are solved on the whole computational domain, while a source term is introduced to represent the body as a porous medium with very low permeability
<cit.>. This method is based on the idea of penalizing the integration points which are located inside the immersed body using source terms, instead of using a boundary condition over the surface of the body, as it is the case when using conforming meshes. In this way, a simple Cartesian mesh can be used to solve the Navier–Stokes equations over arbitrary geometries, reducing the time required for the mesh generation.
In the IBVP method, a mask function is required to discriminate between those points that are inside or outside the immersed body, as it is represented in <ref>. Then, source terms are applied to modify the behavior of the flow field inside the immersed body. Given the governing equations for a compressible viscous fluid:

    ut+∇·f(u,∇u)=s(u)

where s(u) is the IBM source term. This source term is written as

    s(u;χ,η)=
      χ/η[                0;         ρ(v-v_s); 12ρ(v·v-v_s·v_s) ],

where v_s is the velocity of the moving geometry (note that for static object v_s=0 andχ represents the mask function and distinguishes between the fluid region Ω_f and the solid region Ω_s:

    χ(x,t) =
      
        1,    if x∈Ω_s,
    
        0,    otherwise.

and 0<η≪1 is the penalization parameter. If a RANS Spalart–Allmaras model is used, the eddy viscosity, ν̃, is set to zero inside the body <cit.>:

    s_SA(u;χ,η)=-χ/ην̃.

Furthermore, a wall model has been developed to improve the accuracy of RANS simulations. This model selects a set of points near the immersed body where the wall function will be applied, which are called wall points. For each wall point, the normal vector to the body's surface is computed and a new point is defined, called image point. The image points will provide the required information to the wall function to update the state vector at the wall points. However, image points (in general) are not integration points of the computational mesh. Therefore, the state vector at the image points has to be interpolated from a set of interpolation points. These interpolation points, which are computed through a kd-tree algorithm, are the N-nearest integration points to each image point. An overview of this approach is represented schematically in <ref>. The distance which is used to discriminate between the regions where the wall points and the image points live, is defined as the length of the smallest edge of the smallest cell of the mesh.
To interpolate the state vector at each image point, a Radial Basis Function (RBF) interpolator with a Gaussian kernel is used. Then, a wall function is used to compute the friction velocity at each image point <cit.>:

    v_||^+(y^+)=1/κlog(1+κy^+) + (C-1/κlog(κ))(1-exp(-y^+/11)-y^+/11exp(-y^+/3)),

and

    v_IP,t = v_τv_||^+(y^+_IP),  with  
        y^+_IP = v_τd_IP/ν_IP,

where the subscript IP refers to each image point and the subscript t refers to the tangential component of the velocity vector.Finally, the resulting friction velocity v_τ is considered to be constant along the normal direction and, hence, it can be used at the wall points (subscript WP). Taking into account this information, the velocity, the eddy viscosity and the non-dimensional wall distance can be obtained at the wall points:

    y^+_WP   = v_τd_WP/ν_WP, 
    
      v_WP,t   = v_τv^+(y^+_WP), 
    
      v_WP,n   = d_WP/d_IPv_IP,n, 
    ν̃_WP   = κv_τd_WP,

which can be solved including an additional condition of parallelism:

    v_WP,t∥ v_IP,t

The subscript n refers to the normal direction of the velocity vector, κ=0.41 is the von Kármán constant and d is the distance from a point to the immersed body. The new velocity and turbulent viscosity at the wall points are considered by including additional source terms for the momentum, energy and turbulent viscosity equations. Finally, this wall model as implemented in CODA along with the IBVP scheme has proved to improve the results of the standard RANS simulation if the y^+ at the wall points has a value (approximately) of 60 or less.



  §.§.§ IBVP Workflow in CODA


The implemented IBVP workflow in CODA consists in the following steps:

  
  Mesh preprocessing In the mesh preprocessing stage, an octree-based Cartesian mesh with refinement around the immersed body is created. This task is done by an external library. The details of how this mesh is generated are given in <ref>.
  
  IBVP preprocessing The immersed bodies are imported by the function immersedBody.ImmersedBodyGroup from the FSImmersedBody plugin. This plugin contains utilities to place solid bodies inside a CFD mesh, mark points inside them, calculate distances, and generate integration points and weights over them.
  
  CFD simulation CODA loads the octree mesh and calls the preprocessing functions for generating the IBVP mask, and export it with the function ExportIBVPMaskNSToFSMeshOfExportedIntegrationPoints if the Navier–Stokes equations are used, and with the function ExportIBVPMaskRANSToFSMeshOfExportedIntegrationPoints if the RANS equations are employed. Both functions are part of
  immersedBody.ImmersedBodyGroup and then performs the CFD simulation.
  
  IBVP postprocessing The function Postprocessing.ComputeLiftAndDrag computes the lift and drag coefficients. More details are discussed in <ref>.







§ PREPROCESSING AND POSTPROCESSING




 §.§ Preprocessing


Immersed boundary methods use a background mesh for solving a fluid-structure interaction problem. In 3D applications, this background mesh is typically a simple Cartesian mesh made of hexahedral elements. However, it can be also unstructured and made of hybrid elements (hexahedra, prisms, tetrahedra, etc.), or even a body-fitted mesh. In this section we focus on the construction of a Cartesian background mesh built out of hexahedra. This kind of mesh is assembled in an unstructured fashion by the preprocessing tool.



  §.§.§ Background Mesh Generation


The preprocessing tool has two ways for generating the initial background mesh: the first way consists in creating a new Cartesian mesh from scratch, and the second way consists in importing a mesh generated by the software GMSH (strictly speaking, the preprocessing tool is not generating a mesh but importing an already generated one by an external tool). The preprocessing tool is only capable of handling meshes made of hexahedral elements.
Next, we will describe the tasks the preprocessing tool does for constructing an unstructured and adaptively refined Cartesian mesh.



  §.§.§ Loading Parameters for Background Mesh Generation


The first task to be performed is to load from file the parameters of the mesh to be constructed. Depending on the type of mesh to be generated (built-in or imported mesh), different parameters are read-in.For the built-in mesh type, the preprocessing tool requires some parameters that specify the type of built-in box (standard box, curved box, box deformation function, type of target geometry, mesh stretching factors in each spatial direction, etc.). The number of elements in each direction of the Cartesian box, the degree of the mesh (for high-order curved meshes), and the kind of curved boundary are also required.For the imported mesh type, the filename of the GMSH file is required. Additional parameters are necessary for setting up the name and location of boundary conditions and the output format of the mesh file (formats supported: HDF5, GMSH, and TECPLOT).



  §.§.§ Creation of the Background Mesh


Once the mesh parameters are loaded, the preprocessing tool proceeds to create the background mesh, either a built-in or an imported mesh. See Figure 5 for a built-in and an imported mesh generated by the preprocessing tool. The typical workflow of the preprocessing tool is as follows:

  
  * Creation of interpolation matrices between different polynomial bases.
  
  * Creation of a reference box.
  
  * Creation of a curved box from the reference box.
  
  * Application of the stretching functions on each spatial direction, if required.
  
  * Application of mesh deformation functions, if required.
  
  * Application of mesh rotation matrices, if required.
  
  * Storing the final mesh in arrays.
  
  * Numbering of nodes, elements, elements faces and boundary faces.
  
  * Creation of the main data structure of the mesh: the list of elements with their corresponding information (nodes ID, elements ID, local faces definition, boundary conditions). This elements list is a dynamical data structure, namely a linked list. Taking into account that the number of elements increases dynamically according to the desired target mesh size around the immersed geometry, this data structure makes possible a straightforward application of the mesh refinement algorithm.




  §.§.§ Mesh Refinement around Immersed Geometries


Once the background mesh is generated and stored in the element list data structure, the preprocessing tool will refine this mesh provided the geometry of the body in STL format.



  §.§.§ Loading Parameters for the Mesh Refinement


From the same parameter file used in the generation of the background mesh, the preprocessing tool will load all the information necessary for the mesh refinement task: the maximum level of refinement around the immersed geometry, the number and location of special regions that need to be refined, and their corresponding levels of refinement, the type of refinement (either isotropic or anisotropic), and the filenames of the STL files used to describe the geometry of the immersed bodies. After reading this information, the preprocessing tool will perform the mesh refinement around the immersed geometry.



  §.§.§ Importing Geometries from STL Files


The immersed geometries are described by STL files. These files contain an unstructured triangulated surface (triangles unit vectors and the coordinates of their corresponding vertices), as they are generated by a CAD software. The preprocessing tool is responsible of loading the tessellated surface of the immersed body for its further employment in the refinement stage.



  §.§.§ Mesh Refinement around Immersed Geometries


The preprocessing tool has the capability of refining the background hexahedral mesh in special regions within the computational domain (for instance, wake regions) and also around the triangulated surface of the immersed geometry. The preprocessing tool will follow the next steps for refining the background mesh:

  
  * Flagging the elements that require refinement. For the refinement of special regions in the computational domain, the flagging procedure marks all hexahedral cells inside the region coordinates provided in the parameter file. For the refinement of the hexahedral elements intersecting the triangulated surface of the immersed body, an efficient algorithm for the triangle-box overlap is implemented. Before the mesh refinement takes place, this triangle-box overlap is tested for all triangles of the surface triangulation, and if such overlap takes place, the involved elements store the overlapped triangles. During the mesh refinement loop (refinement level l>0), this triangle-box overlap is tested only for those hexahedra storing overlapped triangles.
  
  * Flagging the elements for balancing. The FSDM, and therefore CODA can handle hanging nodes for different face element types (triangular and quadrilateral faces). This hanging nodes capability is restricted in FSDM to elements with a 2:1 refinement level ratio. This means that for two neighboring elements with different refinement levels, only a difference of level of at most 1 between the neighboring cells is allowed. The preprocessing tool is aware of this restriction and therefore an algorithm for preserving the 2:1 balancing has been implemented.
  
  * A feature implemented in the preprocessing tool is the refinement of elements around the immersed body up to certain number of neighboring cells far away from the truly overlapped elements with the geometry. This allow us to have a thick layer of small elements around the immersed body, which is useful in the simulation of flow problems with the RANS equations and the employment of wall functions.



  
  * Refinement of hexahedral elements. All flagged cells are splitted into 4 or 8 new children cells (for anisotropic or isotropic mesh refinement, respectively). Hanging nodes can appear and the relation between the face nodes of neighboring elements with different level of refinement has to be computed and exported in such a way that FSDM can handle them. For the hexahedral meshes used within the preprocessing tool and when hanging nodes are present, the Quad2Quad and Quad4Quad interfaces are computed (see <ref>).
  
  * Storing the final mesh in arrays. In these arrays are stored the nodes, elements, conforming faces and non-conforming faces.
  
  * Exporting the mesh arrays to HDF5 format, following the FSDM guidelines for the HDF5 format.




 §.§ Postprocessing


The postprocessing for the IBVP methodology is based on the computation of the lift and drag of an immersed body, known the geometry (as an STL file) and the converged flow field. Two different approaches are considered: integration over the STL and a Voronoi integration.



  §.§.§ Integration over the Original Immersed Body Tessellation


This approach takes advantage of the information within the STL file to perform the integration. First, for each triangle of the tessellation, the fluid variables are interpolated at the barycenter, at the vertices and at the middle point of each edge of the triangle. The interpolation can be performed using different approaches: Radial Basis Functions (RBF), inverse distance or linear interpolation. In general, a simple linear interpolation has shown a good behavior and robustness.Then the viscous stress tensor (including the pressure) τ̃=τ-pI is computed on those points. Finally, the integration over each triangle is performed as follows:

    F_t = S/60( 27τ̃_bary + 3∑_i=1^3τ̃_vert,i + 8∑_i=1^3τ̃_mid,i).

In this case, the integrated force on the triangle t, F_t is computed as a weighted average of the stress tensor (including the pressure), τ̃, projected over the normal vector of the local triangle. The value of τ̃_bary, τ̃_vert,i, and τ̃_mid,i correspond to the projected force at the barycenter, at the vertex i and at the middle point of the side i respectively.



  §.§.§ Integration over a Modified Immersed Body Tessellation


This approach provides an alternative to compute the lift and drag when the STL is too coarse (and the integration over the STL shows a bad quality) or too fine (and the computational cost has to be reduced). In these cases, a new set of points over the surface of the immersed body are defined to perform the integration.First, a region around the body is defined and every integration point of the mesh within that region is projected over the body's surface. Then, an algorithm is used to compute the weights for each point. Once the points over the surface (and their weights) are known, the fluid variables are interpolated on those points.
Finally, the stress tensor (including the pressure) is computed following the same equations shown above and the projected force is integrated as follows:

    F = ∑_i=0^nPointsτ̃_iω_i.

In this case, the global force is computed as the dot product between the projected force at each point, τ̃, and the weights, ω.






§ NUMERICAL COMPUTATIONS


To prove the validity and performance of this approach,  we  present in this section several numerical simulations of different fluid flow problems of increasing dificulty. These tests are the subsonic flow around a NACA0012 airfoil at different angles of attack, the subsonic flow around the MDA30P30N multi-element airfoil, and the more challenging NASA High-Lift Common Research Model (CRM-HL) aircraft (Figure <ref>-top-left, top-right and bottom respectively). Experimental results of most of those problems are reported in the literature and thus will serve as validation of our methodology.  For all simulations, the RANS equations were used for modeling the fluid flow, and the numerical methods employed by CODA were, for the spatial discretization, the second order finite volume method coupled with the immersed boundary volume penalization. A linearized implicit Euler scheme is employed for time discretization. All simulations were executed using 600 cores in the Magerit Supercomputer at Supercomputing and Visualization Center of Madrid (CeSViMa) and 3840 cores in the Marenostrum4 Supercomputer at the Barcelona Supercomputing Center.




 §.§ Flow around a NACA0012 airfoil


The computational domain for the fluid flow around a NACA0012 airfoil is the box with dimensions [-20,+20]×[-20,+20]×[0,+1]. The span of the domain is sufficient to avoid domain confinement effects and wave reflections from the domain boundaries, which could lead to a significant error when computing the lift and the drag force around the airfoil. The extent of the computational domain corresponds to 40c, where c is the chord length and is equal to 1. The domain is discretized initially with n_xn_yn_z=1600 hexahedral elements, where n_x=40, n_y=40, and n_z=1 are the number of elements in x, y, and z direction, respectively. This initial mesh is further refined around the immersed geometry and also in the wake region. The grid spacing obeys h=2^-l, where l is the level of refinement. Sufficient grid resolution around the airfoil is crucial to obtain accurate results for the lift and drag coefficients. The mesh around the immersed geometry is refined up to level of refinement l=12, we get a mesh with around 2.4E6 elements. The wake region is resolved with a refined subregion with refinement level l=8. This region is approximately 20 airfoil chords downstream of the lead edge of the airfoil. The characteristics of the meshes used for the mesh convergence study are shown in <ref>; ranging from a very coarse mesh up to a very fine one. The refinement level l_max listed in the table is the highest level of refinement around the immersed body.

The immersed geometry is located in the center of the computational domain, with its cross-section on plane z, and this corresponds to the NACA0012 profile with chord length c=1. The triangulation of the surface geometry is made of 18524 triangles. The geometry is shown in <ref>, and in <ref> is depicted the mesh around the immersed geometry.



  §.§.§ Subsonic flow at M_∞=0.15 and α=0 and α=10


The first test corresponds to the subsonic flow around an NACA0012 airfoil. This problem has become a classic test case for RANS solvers due to the simple geometry and large amount of available numerical and experimental data <cit.>. It is used primarily for the analysis of turbulence models, testing their convergence properties and their effect on the accuracy. We perform the simulations for angles of attack α=0 and α=10. The flow satisfies the following conditions: a gas with adiabatic index γ=1.4 and Reynolds number Re=6E6 flows with a freestream Mach number M_∞=0.15. The non-dimensional density is set to ρ=1, and the non-dimensional pressure to p=1. On the boundaries of the computational domain were set the following boundary conditions: left face is set to inflow, the right face to outflow, the front and back faces were set to periodic, and the top and bottom to far field. The inflow pressure assumes the value of the stagnation pressure and the outflow pressure to the static pressure.

In the <ref> are depicted the contour plots of the velocity magnitude for the flow past the NACA0012 profile at angle of attack α=0. The mesh employed in the simulation is the very fine mesh with element's size around the body h_body=2.44e-4, element's size in the wake region h_wake=3.91e-3. The Spallart–Allmaras model has been activated. We observe the flow field is symmetric as the fluid past the airfoil. As the fluid flows over the blunt trailing edge, a minor amount of flow separation occurs, but it reattaches shortly after it passes the trailing edge and recovers the mean velocity in the distant far field.In the <ref> are depicted the contour plots of the velocity magnitude for the flow past the NACA0012 profile at angle of attack α=10. The mesh employed is the same as for angle of attack α=0. The Spallart–Allmaras model has been activated. We observe that the stagnation point is shifted downward (relative to the airfoil), reflective of the shifted angle of attack. The fluid velocity above the airfoil is clearly different than that below the airfoil, showing a low-velocity region toward the upper side of the trailing edge, which means flow separation.

Now we analyze the pressure coefficient plots. The <ref> and <ref> show pressure coefficient C_p over the airfoil surface for all considered meshes and also for the cases with wall model deactivated and activated. Experimental data are also plotted and they were obtained from <cit.>. For an angle of attack α=0, and for all considered meshes, the pressure coefficient curve is very close to the experimental data, but some unphysical oscillations are clearly perceptible in all simulations, both with the wall model switched off and switched off. These oscillations tend to diminish as the mesh becomes finer. For an angle of attack α=10, the pressure coefficient curves for the very coarse, coarse and medium meshes show a large discrepancy regarding the experimental data, especially on the leading edge upper surface. Only the fine and very fine meshes present a good agreement with the experiment. The activation of the wall model helps to reduce this disagreement between the numerical results and experimental data for the coarse meshes.Values of the lift and drag are shown in <ref> for angle of attack α=0 and in <ref> for angle of attack α=10. For the simulation with angle of attack α=10, the lift coefficient C_L increases monotonically as the grid is refined, getting closer to the experimental value, while the drag coefficient C_D appears to decrease monotonically and also it is approaching its numerical value to the experimental one. This behavior occurs for the simulations with the wall model deactivated and activated. For the simulation with angle of attack α=0, we do not observe monotonicity in the evolution of the lift coefficient C_L and the drag coefficient C_D as the mesh is refined.





 §.§ Subsonic flow around an MDA30P30N multi-element airfoil


The third test in the benchmarking is a three-element configuration, specifically the McDonnell Douglas 30P-30N landing configuration (MDA30P30N). Many experimental and computational studies have been performed for the flow past multi-element airfoils in the last decades <cit.>. Accurate prediction of the flow over multi-element airfoils during high-lift operations can improve the performance and the safety factor of aircrafts. The flow around multi-element airfoils is complex and it is well known to be dominated by different flow mechanisms at different operating conditions, making rather difficult to accurately predict high-lift flow fields. The geometry of the configuration is depicted in <ref>. The leading-edge slat and the trailing-edge flap have a deflection angle of 30. This airfoil has been tested extensively at NASA Langley Low Turbulence Pressure Tunnel at various Reynolds and Mach numbers and also has been numerically simulated using a wide range of numerical techniques for solving the Navier–Stokes equations along with different turbulence models.The computational domain is similar to the employed in the NACA0012 case, that is, a box with dimensions [-20,+20]×[-20,+20]×[0,+1]. The domain is discretized initially with n_xn_yn_z=1600 hexahedral elements, where n_x=40, n_y=40, and n_z=1. The mesh is then refined around the immersed geometry and in the wake region. The mesh around the immersed geometry is refined up to level of refinement l=12, and the wake region is resolved with a refined subregion with refinement level l=8. This region is approximately 20 airfoil chords downstream of the lead edge of the airfoil. In <ref> are shown the characteristics of the meshes used for the mesh convergence study, ranging from a medium mesh up to a very fine mesh. The refinement level l_max listed in the table is the highest level of refinement around the immersed body.

The immersed geometry is located in the center of the computational domain, with its cross-section on plane z, and this corresponds to the MDA30P30N profile with chord length c=1.2. The triangulation of the surface geometry is made of 66320 triangles. In <ref> is depicted the mesh around the immersed geometry. The inner elements are not shown.
Regarding the flow conditions, the gas has adiabatic index γ=1.4 and it flows with Reynolds number Re=9E6 and freestream Mach number M=0.2. Simulations were run with the airfoil at an angle of attack α=8. The non-dimensional density is set to ρ=1, and the non-dimensional pressure p=1. On the boundaries of the computational domain were set the following boundary conditions: left face is set to inflow, the right face to outflow, the front and back faces were set to periodic, and the top and bottom to far field. The inflow pressure assumes the value of the stagnation pressure and the outflow pressure to the static pressure.

Contour plots of velocity magnitude and streamtraces from the numerical solution on the very fine mesh are shown in <ref> and <ref>. A small separation region is present on the suction side of the upstream side of the trailing edge of the flap.

In the <ref> is presented the grid convergence in terms of the pressure coefficient C_p over the multi-element airfoil surface, with wall model deactivated and activated. Experimental data were taken from <cit.>. In the main part of the airfoil, the pressure coefficient curve is in good agreement with the experimental data, except at its trailing edge. We can observe a remarkable discrepancy with respect to the experimental data on the upper surface of the slat and the flap.



 §.§ Subsonic flow around NASA high-lift CRM


The last test we consider in this work is the NASA CRM high-lift configuration (CRM-HL), selected from the Fourth AIAA CFD High Lift Prediction Workshop <cit.>. The NASA Common Research Model has been used mainly in the Drag Prediction Workshop (DPW) and the High Lift Prediction Workshop (HLPW) assess the accuracy of numerical methods in the prediction of aircraft forces and moments. The database from these workshops allow to assess the robustness and effectiveness of state-of-the-art numerical programs and turbulence modeling techniques using Navier–Stokes solvers.The CRM configuration was initially designed by Boeing and further manufactured and tested by NASA <cit.>.
The NASA CRM consists of a contemporary supercritical transonic wing and a fuselage that is representative of a wide-body commercial transport aircraft. Several experiments of the CRM have been carried out in several research facilities <cit.>.The computational domain used for this simulation is a box with dimensions [-20,+20]×[-20,+20]×[0,+20]. The domain is discretized initially with n_xn_yn_z=32000 hexahedral elements, where n_x=40, n_y=40, and n_z=20. The mesh is then refined around the immersed geometry and in the wake region. The mesh around the immersed geometry is refined up to level of refinement l=10, and the wake region is resolved with a refined subregion with refinement level l=8. In <ref> are shown the characteristics of the meshes used for the mesh convergence study, ranging from a very coarse mesh up to a medium mesh. The refinement level l_max listed in the table is the highest level of refinement around the immersed geometry. The geometry of the configuration is depicted in <ref>.

The immersed geometry is located in the center of the computational domain. The triangulation of the surface geometry is made of 74331 triangles. In <ref> is depicted the mesh around the immersed geometry.
The Reynolds number based on the mean aerodynamic chord (MAC=7.005) is Re=5.6 million, and the angle-of-attack is α=8 with an incoming Mach number of M_∞=0.2. On the boundaries of the computational domain were set the following boundary conditions: left face is set to inflow, the right face to outflow, the front and back faces were set to periodic, and the top and bottom to far field. The inflow pressure assumes the value of the stagnation pressure and the outflow pressure to the static pressure.






§ CONCLUSIONS


The immersed boundary volume penalization method for the Navier–Stokes and the RANS equations has been implemented in the CFD solver CODA. The immersed boundary method is compatible with the finite volume, modal discontinuous Galerkin and discontinuous Galerkin spectral element methods and with the explicit and implicit temporal schemes available in CODA.
An efficient preprocessing tool for the construction of unstructured hexahedral meshes with adaptive mesh refinement around immersed geometries has been developed. The octree meshes generated with the tool are imported in CODA and used for simulating complex aerodynamics flow problems. Immersed boundary methods require a very simple meshing process: only the box dimension and number of elements of the Cartesian mesh and the geometry of the body are required as an input. The automatic mesh generation process is able to automatically refine the mesh near the body. Besides, the immersed boundary technology can be used straightforward with an external predefined mesh, without any additional modification.Several numerical computations were presented and discussed: subsonic flow past a NACA0012 airfoil, subsonic flow past an MDA30P30N multi-element airfoil and subsonic flow around the NASA high-lift CRM aircraft. The computations are in very good agreement with their corresponding experimental data for the pressure coefficient curve, and the drag and lift coefficients when a wall model for the RANS equations is used.Even though the immersed boundary methodology works within the CFD solver CODA with very good results, a critical disadvantage of these techniques with respect to a body-fitted approach is that the mesh refinement around the immersed geometry leads to a very high number of elements when a high resolution of the boundary layer is desired. Anisotropic mesh refinement reduce the number of elements in comparison with isotropic mesh refinement, but the mesh size is still prohibitive.As a consequence of these considerations, we will explore also alternative approaches to automatic mesh generation for IBM. In particular, future work will focus on hybrid IBM approaches, in which only geometrical details (control surfaces, ice shapes, etc) will be immersed in a body-fitted, wall-resolved mesh of a clean aerodynamic surface. This approach would benefit from IBM advantages (several configurations of the geometrical details can be tested on the same body-fitted background mesh) without handling meshes of very big sizes.






§ ACKNOWLEDGMENTS


Jonatan Núñez, David Huergo, Esteban Ferrer and Eusebio Valero acknowledge the funding received by the Grant NextSimAEI10.13039501100011033 and H2020, GA-956104. This project has received funding from the Clean Aviation Joint Undertaking under the European Union’s Horizon Europe research and innovation programme under Grant Agreement HERA (Hybrid-Electric Regional Architecture) no. 101102007. Views and opinions expressed are, however, those of the author(s) only and do not necessarily reflect those of the European Union or CAJU. Neither the European Union nor the granting authority can be held responsible for them.
Esteban Ferrer and Eusebio Valero acknowledge the funding received by the Grant DeepCFD (Project No. PID2022-137899OB-I00) funded by MCINAEI10.13039501100011033 and by ERDF A way of making Europe.
Esteban Ferrer would like to thank the support of the Comunidad de Madrid and Universidad Politécnica de Madrid for the Young Investigators award: APOYO-JOVENES-21-53NYUB-19-RRX1A0.
Finally, all authors gratefully acknowledge Universidad Politécnica de Madrid (<www.upm.es>) for providing computing resources on Magerit Supercomputer. The authors also thankfully acknowledge the computer resources at MareNostrum and the technical support provided by Barcelona Supercomputing Center (RES-IM-2022-3-0023).
















§ GOVERNING EQUATIONS




 §.§ Navier–Stokes Equations


The Navier–Stokes equations can be written as a hyperbolic-parabolic system of conservation laws in the following way

    ut+∇·f(u,∇u)=0.

The vector of conserved quantities is defined by u(x,t)=(ρ,ρv,ρE), where ρ is the mass density, v=(v_x,v_y,v_z) is the velocity vector, and E is the total energy. The physical flux is defined by

    f(u,∇u)=
        f^A(u)-f^D(u,∇u),

with the convective and viscous fluxes given, respectively, by

    f^A(u)=
        [      ρ v; ρ v⊗v+pI;  v(ρE+p) ],  f^D(u,∇u)=
        [     0;    -τ; τ·v-q ].

The viscous stresses are described by the stress tensor τ, defined by

    τ=2μS^D=μ(∇v+(∇v)^⊺-2/3(∇·v)I),

where S^D is the deviatoric component of the strain-rate tensor

    S=1/2(∇v+(∇v)^⊺)

and 

    q=-k∇T,   k=μc_p/Pr.

The equation of state

    p=(γ-1)(ρE-1/2ρv^2),

where p is the static pressure, γ is the ideal gas index, μ is the dynamic viscosity, T is the temperature, k is the thermal conductivity, and μ denotes the dynamic viscosity coefficient. The kinematic viscosity coefficient is defined by the formula

    ν = μ/ρ.




 §.§ Reynolds-Averaged Navier–Stokes Equations


The Navier–Stokes equations govern the motion of fluids in both laminar and turbulent regimes. Because the large range of spatial and temporal scales present in turbulent flows, solving directly the Navier–Stokes equations is prohibitively expensive. Therefore the Reynolds-averaged Navier–Stokes (RANS) equations are solved instead to model steady turbulent mean flows. The RANS equations couple the mean flow equations with the one-equation turbulence model of Spallart–Allmaras. These equations can be written also as a hyperbolic-parabolic system of conservation laws with source term in the following way

    ut+∇·f(u,∇u)=s(u).

Here, u is the vector of time-averaged conservative variables over a given time interval, and has the following components: u(x,t)=(ρ,ρv,ρE,ρν̃), where ρ is the time-averaged mass density, v=(v_x,v_y,v_z) is the time-averaged velocity vector, E is the time-averaged total energy, and ρν̃ is a new conservative variable relating the time-averaged mass density and the eddy viscosity ν̃. The advective and diffusive components of the physical flux associated to the quantity ρν̃ are given by

    f^A[ρν̃]=ρν̃v,  f^D[ρν̃]=1/σ(μ+f_n1ρν̃)∇ν̃

In the diffusive fluxes, the turbulent stress tensor τ_t and the turbulent
heat fluxes q_t are added respectively to τ and q appearing in the mean flow equations:

    τ_t=2μ_tS^D,  q_t=-μ_t/Pr_tc_p∇T

where Pr_t=0.9 is the turbulent Prandtl number and μ_t is the turbulent dynamic viscosity

    μ_t=
        ρν̃f_v1(χ)    for ν̃≥0,
    
          0    for ν̃<0,
    ,  
      f_v1(χ)=χ^3/χ^3+c_v1^3,  χ=ρν̃/μ.

The source terms only have nonzero component in the equation for the turbulent variable ρν̃:

    S[ρν̃]=-ρ(P-D)-c_b2/σρ∇ν̃+1/σ(ν+f_n1ν̃)∇ρ·∇ν̃,

where the production and destruction terms, P and D, are defined by:

    P   =
        
          c_b1(1-f_t2)ω̃ν̃   for ν̃≥0,
    
          c_b1(1-c_t3)ων̃   for ν̃<0,
    
    
      D   =
        (c_w1f_w-c_b1/κ^2f_t2)(ν̃/d)^2   for ν̃≥0,
    
          -c_w1(ν̃/d)^2   for ν̃<0,

and

    f_n1=
        
          1    for ν̃≥0,
    c_n1+χ^3/c_n1-χ^3   for ν̃<0,
    ,  
      f_t2=c_t3exp(-c_t4χ^2),   f_w=g(1+c_w3^6/g^6+c_w3^6)^1/6

with

    g=r+c_w2(r^6-r),   r=min(r_lim,ν̃/ωκ^2d^2),

and d is the distance to the nearest wall and ω the vorticity magnitude. The modified vorticity magnitude ω̃ is given by

    ω̃=
        ω+ω̅   for ω̅>-c_v2ω,
    ω+ω( c_v2^2+c_v3ω̅)/(c_v3-2c_v2)ω-ω̅   for ω̅<-c_v2ω,

where ω̅ and f_v2 are given by

    ω̅=ν̃/κ^2d^2f_v2,   f_v2=1-χ/1-f_v1

For the sake of completeness, we give the values of the constants in the above expressions:
c_v1=7.1,
σ=2/3,
c_b1=0.1355,
c_b2=0.622,
κ=0.41,
c_w2=0.3,
c_w3=2,
r_lim=10,
c_t3=1.2,
c_t4=0.5,
c_v2=0.7,
c_v3=0.9,
c_n1=16.




