# 4AIAA CFD High Lift Prediction Workshop results using metric-based anisotropic mesh adaptation

Frédéric Alauzet, Francesco Clerici, Adrien Loseille, Matthieu Maunoury,

Lucien Rochery, Cosimo Tarsia-Morisco, Lucille-Marie Tenkes, Julien

Vanharen

► To cite this version:

Frédéric Alauzet, Francesco Clerici, Adrien Loseille, Matthieu Maunoury, Lucien Rochery, et al. 4AIAA CFD High Lift Prediction Workshop results using metric-based anisotropic mesh adaptation. AIAA Fluid Dynamics Conference 2022, Jun 2022, Chicago, IL, USA, Unknown Region. 10.2514/6.2022-3521. hal-04401293

HAL Id: hal-04401293

https://inria.hal.science/hal-04401293v1

Submitted on 13 Dec 2024

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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
4th AIAA CFD High Lift Prediction Workshop results
using metric-based anisotropic mesh adaptation

Frédéric Alauzet*, Francesco Clerici†, Adrien Loseille‡,

Matthieu Maunoury§, Lucien Rochery¶, Cosimo Tarsia Morisco∫,

Lucille-Marie Tenkès;** and Julien Vanharen††

GAMMA Team, Inria Saclay Ile-de-France, Palaiseau, France

This article presents CFD results obtained with the metric-based mesh adaptation plat-
form for the 4th AIAA CFD High-Lift Prediction Workshop on the new CRM-HL geometry.
We emphasize the need for a strong solver to achieve iterative convergence for this complex
case and we propose a α-continuation mesh adaptation algorithm to improve the prediction.
Two test cases of the workshop have been carried out: the mesh convergence study at an
angle of attack of 7.05 degrees and the C$_{L,max}$ study.

I. Introduction

Mesh generation for accurate and reliable numerical simulations of Reynolds-Averaged Navier-Stokes (RANS) applications is a time consuming and tedious task. Traditional processes rely on the experience and intuition of a skilled engineer to predict the flow features and prescribe mesh requirements accordingly. Following meshing guidelines slows down the mesh generation process leading to a prohibitive cost in CPU time in the numerical simulation pipeline. This lack of automation is an impediment for many applications such as shape optimization. Moreover, as realistic applications consider complex geometries and flow patterns, the analysis of the various physical features, and thus, the prescription of an appropriate mesh, is almost impossible. The provided meshes for complex configurations generally rely on a priori refinement of large zones, leading to unnecessary large computational resource consumption. In this context, mesh adaptation strategies, that automatically analyze and adapt to the flow solutions can bring significant improvement in terms of CPU times, memory consumption and accuracy through error control. Indeed, the simulation will then be independent of the initial mesh, and more importantly, it is an automated process, which removes human intervention from the simulation pipeline.⁷,⁸

Metric-based anisotropic mesh adaptation has been very successful in predicting the lift for the High-
Lift Common Research Model (CRM-HL) provided at the 3rd AIAA CFD High-Lift Prediction Workshop
(HLPW3).<sup>2,3</sup> The obtained results are shown in Figures 1 and 2 for the 16 degrees angle of attack test case.

The convergence of the lift value throughout the whole mesh convergence study is shown in Figure 1 (left picture). It shows, in red, the evolution of the lift value for all the computations for all the complexities, and, in blue, the final retained values obtained for each complexity. We note that it is very advantageous to start and perform many iterations at small complexities because most of the convergence is achieved by the coarse adapted meshes, and these iterations are inexpensive compared to larger complexities, whereas only a minimum number of iterations is necessary for the finer adapted meshes. Moreover, this plot points out that we obtained mesh a converged solution for this 3D case.

<footnote>*Senior Researcher, Frederic.Alauzet@inria.fr</footnote>

<footnote>†PhD Student, Francesco.Clerici@inria.fr</footnote>

<footnote>‡Researcher, Adrien.Loseille@inria.fr</footnote>

<footnote>§Research Engineer, Matthieu.Maunoury@inria.fr</footnote>

<footnote>¶PhD Student, Lucien.Rochery@inria.fr</footnote>

<footnote>‖Research Engineer, cosimo.tarsia-morisco@inria.fr</footnote>

<footnote>**PhD Student, lucille-marie.tenkes@inria.fr</footnote>

<footnote>††Research Engineer, Julien.Vanahren@onera.fr</footnote>
Figure 1 (right) compares the results obtained with the metric-based anisotropic mesh adaptation to all the workshop results (green lines). Particularly, red lines show the results using the feature-based error estimate based on controlling the interpolation error of the Mach number in $L^4$ norm. Blue lines show the results using the viscous goal-oriented error estimate based on the control of the lift output functional. It points out the early capturing property of the mesh-adaptive solution platform (i.e., an accurate prediction is obtained for relatively coarse adapted meshes), and it emphasizes the obtention of a mesh-converged solution. Two main considerations can be done by observing this chart. The first one is the evident superiority of the goal-oriented approach compared to the feature-based one for this kind of application. The second one is the improvement of the results over the years thanks to the progress in numerical schemes (a less dissipative limiter and a better iterative convergence thanks to an improved Newton's method).

Figure 2 is interesting on several aspects. We compare the lift prediction obtained with the mesh-adaptive solution process (blue stars, each star being a different mesh size) to all the results obtained during the HLPW3 on the coarse (8M vertices), the medium (26M vertices), the fine (70M vertices) and the x-fine (206M vertices) workshop's meshes (red squares). We observe that a result similar to the fine grid is obtained with a 1.3M vertices adapted mesh, and the 2.73M and 5.40M vertices adapted meshes provide results similar to the x-fine mesh that have been run by only five participants. We are predicting a lift value of 2.382 which is above the workshop's value obtained on the x-fine mesh, i.e. 2.365. Looking at the overall solution trend, lift values tend to increase with the mesh size, so it is very likely that even the x-fine mesh solutions are not yet mesh-converged. In conclusion, we achieve the same accuracy as meshes manually generated following the best practice meshing guidelines, with 75 times less vertices for the lift prediction.

In this paper, we propose to analyze the results obtained for the 4th AIAA CFD High-Lift Prediction Workshop (HLPW4) on the new CRM-HL geometry. This new geometry is a lot more complicated than the one of the HLPW3 due to the presence of leading-edge brackets, trailing-edge support fairings, pylon, nacelle, and nacelle chine. In consequence, the physics is more complex, making this case a lot more challenging. Our objective is to solve such a complex problem using the mesh-adaptive solution platform with the following constraints:

- pure tetrahedra adapted meshes are considered (no pseudo structured boundary layer mesh with prisms are allowed)

- no *ad hoc* size prescription in the boundary layer region, i.e., the size prescription should only be provided by the error estimate

- using a minimum number of vertices to achieve the correct behavior. From the HLPW3 results, we chose to analyse the physical behavior predicted using adapted meshes composed of 10M vertices, which should

<figure><img src="image_2.png" /><figcaption>Figure 1. HLPW3 HLCRM 16° case. Left, convergence history of the total lift value $C_L$ for the viscous goal-oriented error estimate throughout the whole mesh convergence study. In red, the convergence of the total lift at each complexity and, in blue, the global convergence of the total lift by retaining the final lift value for each complexity. Right, convergence of the total lift for the mesh-adaptive solution process with the feature-based error estimate (red lines) and the viscous goal-oriented error estimate (blue lines) with respect to all workshop entries (green lines).</figcaption></figure>
be a reasonable mesh size for this configuration.

We have encountered several difficulties. Some are due to the above constraints, some have been solved, others are still ongoing. The first difficulty is that steady state for such a complex geometry and flow can be difficult to achieve. At the beginning, we observed oscillating aerodynamic coefficients (drag and lift), in other words the solution was not steady (vortex shedding was observed at several locations). This issue has been solved by developing a stronger solver in order to fulfill iterative convergence. The second difficulty is that the flow can be separated at several places on coarse adapted meshes. Some of these flow separations are physically correct but some of them are due to a lack of mesh refinement in these areas and are unphysical. We were able to reattach the flow in some of these regions but others are more complicated and stay separated during the process. To solve this issue, we have set-up a new α-continuation mesh adaptation algorithm where the angle of attack is increased during the mesh-adaptive solution process.

The paper is outlined as follow. Section II presents the *C*-continuation and the *α*-continuation mesh adaptation algorithms. Section III describes each component of the mesh-adaptive solution process. Section IV provides the numerical results. We end the paper with some conclusions and perspectives.

<figure><img src="image_3.png" /><figcaption>Figure 2. HLPW3 HLCRM 16° case. Comparison of the lift prediction results obtained with the mesh-adaptive solution platform using the viscous goal-oriented error estimate (blue stars) with respect to all the HLPW3 results (red squares) obtained on the coarse (top left), medium (top right), fine (bottom left) and x-fine (bottom right) meshes.</figcaption></figure>
## II. Anisotropic mesh adaptation algorithm with mesh-convergence analysis

Mesh adaptation is a non-linear problem where the couple formed by the mesh and the solution needs to be converged at the same time. The goal is to generate an optimal adapted mesh for the computed solution for a prescribed mesh size (i.e., a prescribed number of vertices or mesh complexity). Therefore, an iterative process is required which is achieved by means of a mesh adaptation loop starting from an initial mesh $H_0$, an initial solution $W_0$, an initial adjoint state $W_0^*$ if goal-oriented mesh adaptation is considered, and a given mesh complexity $\mathcal{C}$ (the continuous counterpart of the mesh size).

At each step of the mesh adaptation loop, a metric tensor $M_i$ is computed from the triple $(H_i, W_i, W_i^*)$ and the given mesh complexity $\mathcal{C}$, using the selected error estimate. The metric tensor field $M_i$ contains information on sizes and directions of the elements of the adapted mesh we seek. This information is then used by the remesher to generate a new adapted mesh $H_{i+1}$.¹⁵ Then $W_i$ is interpolated on $H_{i+1}$ to obtain $(W^0)_{i+1}$ which is then used as a restart solution for the next flow solution of the mesh adaptation loop.¹ In the case of goal-oriented mesh adaptation, the adjoint state $W_i^*$ can be also interpolated on the new mesh $H_{i+1}$ to obtain $(W^{*,0})_{i+1}$ which is used as a restart for the next adjoint solution. Restarts are important not to waste time in the adaptive process and to reuse the previous work done as much as possible. Finally, a new solution $W_i$ and a new adjoint $W_i^*$ fields are computed. Solution and adjoint fields are converged to machine zero at each iteration. This iterative process is depicted by the inner while loop in Algorithm 1.

The convergence criteria of step 1(f) is up to the expectations of the user, it specifies when the couple mesh-solution is considered as converged for the current complexity in the process. In this work, for aeronautics applications, we consider that the couple mesh-solution is converged at the given complexity if the lift coefficient, the pressure component of the drag and viscous component of the drag are not varying by a given percentage $\epsilon$ for three consecutive iterations. Usually, we choose $\epsilon$ between 0.001 (i.e., 0.1%) and 0.01 (i.e., 1%). It is important to check the convergence of all the components of the drag otherwise we may get early break because of compensation between the pressure and the viscous components on the total drag value.

In the following, two extensions of the mesh adaptation algorithm are presented. The first one is designed to perform convergence study and the second is tailored to high-lift application to compute $C_{L,max}$ at a given mesh size.

### A. Convergence study algorithm: $\mathcal{C}$-continuation method

In the context of a mesh convergence analysis this adaptation loop (step 1) has to be repeated for several increasing mesh complexities $\{C^j\}_{j=1..jmax}$ with $jmax$ the number of considered complexities, i.e., the total number of meshes for the mesh convergence study. An efficient strategy consists in converging the couple mesh-solution for a given complexity and reuse the final mesh, solution and adjoint state to initialize the next computations at an increased mesh complexity. Such a process enables a multiscale resolution of the flow by solving large scale features on coarse adapted meshes (at the smallest complexities) and the fine scale features of the flow on fine adapted meshes (at the largest complexities). This acts like a "multigrid effect" and enables faster convergence on fine adapted meshes. This process is represented by the outer while loop in Algorithm 1. We propose the following strategy to increase the complexity for the mesh convergence study. At each outer loop iteration, the complexity is increased by a factor $\alpha$. In this work, we have set $\alpha = 2$ to multiply the mesh size by a factor 2 when increasing the complexity. Moreover, we have found that it is very advantageous to start at small complexities and to perform many iterations as coarse meshes contribute dramatically in converging the solution, and these iterations are inexpensive in comparison to the largest complexities. The major large scale features of the flow such as shocks, shear layers and wakes are detected, advected and meshed on the coarse adapted meshes. They are then refined alongside the boundary layer as the discretization increases. And, as the turbulent boundary layer takes time to develop, it is very efficient to capture it on coarse adapted meshes.

We call this method $\mathcal{C}$-continuation because the parameters of the simulation stay the same, only the mesh complexity is varying. The $\mathcal{C}$-continuation mesh adaptation algorithm has been used for the mesh convergence study (Section IV A). This algorithm is also very interesting when a full polar is studied because each configuration can be run independently in parallel which is very attractive.

<footnote>a This is not a true multigrid because the coarse adapted meshes are not used in the flow solver to compute the solution.</footnote>
<caption>**Algorithm 1** General mesh adaptation algorithm with mesh-convergence analysis</caption>

Initial mesh $\mathcal{H}_0^0$, solution $W_0^0$, adjoint $W_0^{*,0}$, and complexity $\mathcal{C}^0$

//--- Outer loop to perform the convergence study

while $\mathcal{C}^j \le \mathcal{C}^{jmax}$ do

//--- Inner loop to converge the mesh adaptation at fixed complexity

1. **while** $i \le n_{\text{adap}}$ **do**

(a) Compute optimal metric for the considered error estimate and complexity $\Rightarrow \mathcal{M}_{i-1}^j$

(b) Generate new adapted mesh $\Rightarrow \mathcal{H}_i^j$

(c) Interpolate primal and adjoint states on the new mesh $\Rightarrow (W^0)_i^j$ and $(W^{*,0})_i^j$

(d) Compute primal state $\Rightarrow W_i^j$

(e) Compute adjoint state $\Rightarrow W_i^{*,j}$

(f) **if** (convergence check) **then**

$$i = n_{\text{adap}} + 1$$

**else**

$$i = i + 1$$

**fi**

**done**

2. $\mathcal{H}_0^{j+1} = \mathcal{H}_{n_{\text{adap}}}^j$ ; $W_0^{j+1} = W_{n_{\text{adap}}}^j$ ; $(W^*)_0^{j+1} = (W^*)_{n_{\text{adap}}}^j$ ; $\mathcal{C}^{j+1} = \alpha \cdot \mathcal{C}^j$ with $\alpha > 1$ \\
$j = j + 1$

**done**

However, for high-lift computations at high angle of attack, multiple solutions to the Reynolds-Averaged Navier-Stokes equation have been numerically observed.¹⁰ Both physical and numerical solutions exhibit an hysteresis cycle with regard to the angle of attack. It is generally considered that the correct sought solution is the one with the highest lift value. Most of the time, the flow is attached for the highest lift solution while it is separated for the other ones. This points out the high sensitivity of the flow for such applications and it explains some of the difficulties encountered with metric-based mesh adaptation for high-lift applications.

Algorithm 1 can be used as it is in most of the cases, for instance, it has been used successfully for HLPW3 geometry at 8 and 16 degrees angle of attack or to compute $C_{L,max}$ for the 2-D CRM-HL geometry. Unfortunately, we have encountered difficulties in using Algorithm 1 for the HLWP4 geometry at moderately high angle of attack (11 degrees or above). This is mainly due to geometric effects (such as the brackets) that are responsible of flow separations.

Consequently, many regions with separated flow were obtained on coarse adapted meshes. In 2D, the mesh-adaptive solution process was able to correct itself automatically by reattaching the flow when the mesh complexity increased. This was very positive because the algorithm was able to automatically recover the correct solution. But, in 3D, this occurred in some regions but other regions stay separated even when the mesh complexity increased. For instance, this was the case for the top of the nacelle and the wing tip at an angle of attack of 11 degrees where the flow stayed separated until the end of the adaptation loop, see section IV B. We think that this is due to the error estimate, which is unable to refine properly the mesh in these regions even when the mesh complexity increases. In consequence, a lower lift solution is obtained. In that context, it is not due to an hysteresis phenomena but to a defect of the mesh adaptation process. In such cases, this algorithm cannot be directly applied but a $\alpha$-continuation method can be considered to solve the problem.
## B. Tailored high-lift algorithm: α-continuation method

The idea of the α-continuation algorithm is, for a given complexity, to run all the studied angles of attack α. In that case, only α is varying. We start at a low angle of attack and then the angle of attack is increased during the process. This strategy proved to be more efficient to obtain the highest lift solution at very high angle of attack just before stall (even on coarse adapted meshes). The drawback is that each configuration cannot be run independently and make this process less efficient. In this case, the global algorithm is modified as described in Algorithm 2 where {$\alpha_k$}$_k$ are the prescribed angles of attack. In that case, each complexity can be run independently in parallel but not each angle of attack. The modified initialization using a previous angle of attack final mesh, solution and adjoint state is done at step 1(b).

Note that a combination of Algorithm 1 (for lower angles of attack) and Algorithm 2 (if higher angles of attack) can be designed for improved efficiency. For instance, Algorithm 1 can be applied to all angles of attack. If the flow is separated at a given complexity and angle of attack, then the solution at the same complexity and at the previous angle of attack can be used as a restart for this angle of attack like in Algorithm 2. If the flow re-attaches, then Algorithm 1 is executed for the next complexities using that solution as initialization. If the flow does not re-attach, then the α-continuation solution with the next complexity is considered.

The α-continuation mesh adaptation algorithm has been used for the $C_{L,max}$ study (Section IV B).

<caption>**Algorithm 2** General mesh adaptation algorithm with mesh convergence analysis and angle of attack study</caption>

//--- Outer loop to perform the convergence study

while $C^j \le C^{max}$ do

Initial mesh $\mathcal{H}_0^{j,\alpha_0}$, solution $W_0^{j,\alpha_0}$, adjoint $W_0^{*,j,\alpha_0}$, and angle of attack $\alpha_0$; set complexity $C^j$

//--- Second loop to compute all the angles of attack

1. while $\alpha_k \le \alpha_{max}$ do

//--- Inner loop to converge the mesh adaptation at fixed complexity

(a) while $i \le n_{adap}$ do

i. Compute optimal metric for the considered error estimate and complexity $\Rightarrow \mathcal{M}_{i-1}^{j,\alpha_k}$

ii. Generate new adapted mesh $\Rightarrow \mathcal{H}_i^{j,\alpha_k}$

iii. Interpolate primal and adjoint states on the new mesh $\Rightarrow (W^0)_i^{j,\alpha_k}$ and $(W^{*,0})_i^{j,\alpha_k}$

iv. Compute primal state $\Rightarrow W_i^{j,\alpha_k}$

v. Compute adjoint state $\Rightarrow W_i^{*,j,\alpha_k}$

vi. if (convergence check) then

$$i = n_{adap} + 1$$

else

$$i = i + 1$$

fi

done

(b) if (stall occurs) then

break

else

$$\mathcal{H}_0^{j,\alpha_{k+1}} = \mathcal{H}_{n_{adap}}^{j,\alpha_k} ; \quad W_0^{j,\alpha_{k+1}} = W_{n_{adap}}^{j,\alpha_k} ; \quad (W^*)_0^{j,\alpha_{k+1}} = (W^*)_{n_{adap}}^{j,\alpha_k} ; \quad \text{set } \alpha_{k+1} = \alpha_k + \Delta\alpha$$

fi

done

$$2. C^{j+1} = \beta \cdot C^j \text{ with } \beta > 1$$

done
### III. Mesh adaptation components

#### A. WOLF flow solver

WOLF is a vertex-centered (flow variables are stored at vertices of the mesh) mixed finite-volume - finite-element Navier-Stokes solver on unstructured meshes composed of triangles in 2D and tetrahedra in 3D.

The convective terms are solved by the finite-volume method on the dual mesh composed of median cells. It uses the HLLC approximate Riemann solver to compute the flux at the cell interface. Second order space accuracy is achieved through a piecewise linear interpolation based on the Monotonic Upwind Scheme for Conservation Law (MUSCL) procedure which uses a particular edge-based formulation with upwind elements. A specific low dissipation scheme is considered using combination of centered (edge gradient) and upwind gradients (element gradient). A dedicated slope limiter is employed to damp or eliminate spurious oscillations that may occur in the vicinity of discontinuities. This low dissipative limiter is fully detailed in.³ The viscous terms are solved by the $P^1$ Galerkin finite element method (FEM) which provides second order accuracy.

The implicit temporal discretization considers the backward Euler time-integration scheme. At each time step, the linear system of equations is approximately solved using a Symmetric Gauss-Seidel (SGS) implicit solver and local time stepping to accelerate the convergence toward the steady state. A Newton's method based on the SGS relaxation is very attractive because it uses an edge-based data structure which can be efficiently parallelized.

As for the turbulence model, the negative Spalart-Allmaras (SA-neg) is loosely-coupled to the mean-flow equations, where the mean-flow and turbulence model equations are relaxed in an alternating sequence. The flow solver WOLF is thoroughly detailed in³ with all the associated bibliography.

As stated in the introduction, it can be very difficult to achieve steady state when a complex geometry, such as the one of the HLPW4, is considered for high-lift applications. Indeed, unsteady solutions were obtained (vortex shedding was observed at several places) using the Newton's method based on the SGS relaxation leading to oscillating aerodynamic coefficients. This issue has been solved by developing a stronger solver to fulfill iterative convergence. In consequence, since paper,³ the main improvements have been made on the implicit time integration method.

More precisely, we have observed that force and moment values were oscillating periodically around a fixed values without converging toward that values, assuming that the RANS modeling provides a steady solution in that case and we are not considering an unsteady case. In other words, we were facing Lyapunov stability issues with the current implicit scheme on the HLPW4 geometry. This issue comes from the Symmetric Gauss-Seidel (SGS) implicit solver which is not strong enough to achieve machine-zero convergence in this case. We tried a GMRES implicit solver preconditioned with SGS relaxation but it did not help. We solved that issue by developing a non-linear SGS solver consisting in coupling SGS relaxation with fixed point Defect Correction iterations. This is illustrated in² where we were able to reproduce the same issue on the 2-D CRM-HL airfoil with slat, main, and flap elements on extremely fine adapted meshes where the numerical dissipation was very low. We shown that this issue is solved by using the non-linear SGS implicit solver.

In order to design a true Newton's method, we followed some of the lines of.⁴,⁶,¹⁷ In such process, the CFL evolves depending on the behavior of the Newton's method, and an under-relaxation can applied to the increment obtained after the implicit solve. We observe that zeroing some of the increments was an impediment to the converge to machine zero. Moreover, a key improvement is to check whether the current non-linear residual has been decreased by the implicit solve:

$$ \| \mathbf{R}(W^{n+1}) + \frac{|C|}{\Delta t} (W^{n+1} - W^n) \| \le \| \mathbf{R}(W^n) \|, $$

instead of checking the reduction of the residual: $\| \mathbf{R}(W^{n+1}) \| \le \| \mathbf{R}(W^n) \|$. This change enables the CFL to grow faster resulting in faster convergence to machine zero or, convergence to machine zero when the previous version stalls.

These improvements were another step forward in designing of a finite volume strong solver.

#### B. WOLF adjoint solver

As regards the adjoint state computation, needed for goal-oriented error estimates, the matrix of the linear system is simply the implicit matrix (without the mass matrix) transposed and the right hand-side of the
system is the chosen functional (for instance, drag, lift, ...) exactly differentiated. In particular, for viscous flows, $\mu$ and the stress tensor $\tau$ are exactly differentiated. In order to solve the adjoint system, we use a restarted GMRES preconditioned with SGS relaxation. Note that, it is important to converge the adjoint problem to machine zero to obtain an accurate adjoint state for mesh adaptation.

### C. FEFLO.A local adaptive remesher

FEFLO.A is a generic purpose adaptive mesh generator dealing with 2D, 3D and surface mesh generation. It belongs to the class of metric-based mesh generator which aims at generating a unit mesh with respect to a prescribed metric field $M$. A mesh is said to be unit when composed of almost unit-length edges and unit-volume element.

The adaptive remesher uses on a combination of generalized standard operators (insertion, collapse, swap of edges and faces). The generalized operators are based on recasting the standard operators in a cavity framework.¹³,¹⁵ Additional modifications on the cavity allow to either favor a modification, that would have been rejected with the standard operator, or to improve the final quality by combining automatically many standard operators at once. In addition, the CPU time is also improved and becomes independent of the current modification. The unit speed is around 20,000 points inserted or removed per second on Intel i7 architecture at 2.7 GHz. For robustness purpose, both the surface and the volume mesh are adapted simultaneously, and each local modification is checked to verify that a valid mesh is obtained. For the volume, the validity consists in checking that each newly created element has a strictly positive volume. For the surface, the validity is checked by ensuring that the deviation of the geometric approximation with respect to a reference surface mesh remains within a given tolerance.

The generation of a unit mesh is decomposed into two steps:

1. Generate a unit-mesh: the mesh modification operators are used in the goal to optimize the length of the edges in $M$.

2. Optimization: the mesh modification operators are used to improve the quality $Q_M$.

During surface remeshing, either a P3 background surface is used or a direct CAD queries is used. The CAD kernel is based on EGADS and OpenCascade.⁹

### D. RANS goal-oriented error estimate

A goal-oriented error estimate based on an a priori error analysis has initially been proposed for the inviscid Euler equations in.¹⁴ The main idea was to translate the error on the considered output functional into a weighted interpolation error estimate. Weights are given by derivatives of the adjoint state and interpolation errors are on the Euler fluxes. As we are left with weighted interpolation errors, we can use the continuous mesh framework¹¹,¹² to obtain an analytical expression of the optimal metric field. An extension of this goal-oriented error estimate has been proposed for the laminar Navier-Stokes equations in.⁵ The main advantage of these error estimates in comparison to other goal-oriented error estimates is that the anisotropy of the mesh appears naturally.

From the analysis of the behavior several error estimates for the Reynolds Averaged Navier-Stokes equations,¹⁶ we came up with the following new goal-oriented error estimate for RANS using integration by part and linearization:³

$$|J(W) - J(W_h)| \le \int_{\Omega} \left| -\frac{\partial J}{\partial W} + \sum_{i} \left(\frac{\partial \mathcal{F}_i^E}{\partial W}\right)^T \nabla_{x_i} W^* + \sum_{i,j} \left(\frac{\partial \mathcal{F}_i^V}{\partial \nabla_{x_j} W}\right)^T W_{x_i, x_j}^* \right| |W - \Pi_h W| d\Omega \quad (1)$$

where $J$ is the considered output functional, $W$ is the conservative variables vector, $W^*$ the associated adjoint state, $\mathcal{F}_i^E$ the convective fluxes, $\mathcal{F}_i^V$ the viscous fluxes, $\nabla_{x_i} W^*$ and $W_{x_i, x_j}^*$ the gradient and the hessian of the adjoint state. The error estimate is a weighted sum of $L^1$ interpolation error on the conservative variables where the weights depend on the gradient and the hessian of the adjoint state and on the convective and viscous fluxes. Therefore, we can directly apply the continuous mesh framework to obtain an analytical expression of the optimal metric field.³
IV. Numerical results

We present the computations on the high-lift version of the NASA CRM (CRM-HL) geometry used for the 4rd AIAA CFD High Lift Prediction Workshop (HLPW4) using the mesh-adaptive solution platform. The considered geometry is equipped with inboard/outboard trailing-edge (TE) flap deflection 40°/37° with a nominal 30°/30° inboard/outboard leading-edge (LE) slat setting, nacelle, pylon, nacelle chine, LE brackets, TE support fairings, but without landing gear, nor horizontal or vertical tail.

First, we present the mesh convergence study corresponding to Case 1b of the workshop. Second, we present the $C_{L,max}$ study corresponding to Case 2 of the workshop.

A. Case 1b: Mesh convergence study

We perform the mesh convergence study at an angle of attack of 7.05 degrees corresponding to case 1b of the workshop. The flow conditions are:

<table><thead><tr><td>Mach number</td><td>Angle of attack</td><td>Reynolds number based on MAC</td><td>Temperature (K)</td></tr></thead><tbody><tr><td>0.2</td><td>7.05</td><td>5.49e6</td><td>289.44</td></tr></tbody></table>

The case is considered in SI units. The geometry used is expressed in meters, so that the reference length, based on MAC, is 7.00532 *m* and the reference surface is 191.84477 *m*².

The convergence study is performed using the *C*-continuation mesh adaptation algorithm described in Section II A. For the mesh adaptation, the goal-oriented error estimate given by Relation (1) is considered with the lift coefficient as output functional: $J = C_L$. The convergence study is performed by considering five mesh sizes (number of vertices):

$$ \{1\,250\,000, 2\,500\,000, 5\,000\,000, 10\,000\,000, 20\,000\,000\}. $$

We compare two runs: the first one considering the Newton's method based on SGS relaxation and the second one considering the Newton's method based on SGS relaxation coupled with Defect Correction iterations.

The starting mesh is composed of 135 043 vertices, 685 648 tetrahedra and 82 858 triangles on the surface, see Figure 3. This is a very coarse inviscid mesh without any boundary layer mesh nor any specific refinement for viscous flows. In other words, it has been generated without any meshing guidelines, thus it is very easy and quick to generate. We chose to start from this coarse and clearly unresolved mesh to illustrate the non-dependency of the mesh-adaptive solution process to the initial data.

Figure 4 (left pictures) shows the convergence of the total lift and moment with respect to the mesh size (retaining the final value for each complexity) for both studies. In both cases, the lift value is approaching

<figure><img src="image_4.png" /><figcaption>Figure 3. HLPW4 CRM-HL case. Initial coarse and clearly unresolved mesh only composed of 135 043 vertices, 685 648 tetrahedra and 82 858 triangles on the surface.</figcaption></figure>
1.76, and the moment is approaching -0.35. The simulation using the non linear SGS method shows a better convergence behavior and provides more accurate predictions on the coarser adapted meshes. Figure 6 (left pictures) shows the convergence of the total drag, as well as its pressure and viscous components with respect to the mesh size for both studies. We note that the viscous drag component is similar in both cases pointing out a similar behavior in the boundary layer region, but the viscous drag is not yet mesh-converged, and finer adapted meshes are required. Conversely, a discrepancy is present on the pressure component of the drag. In fact, the solution obtained with the Newton's method based on SGS relaxation is approaching 0.170 while that resulting from the non-linear SGS relaxation is about 0.172 with a better convergence behavior. Nevertheless, finer adapted meshes are required to conclude about it. A similar trend was observed on the HLPW3 where the convergence of the total drag required finer meshes than the convergence of the total lift.³

The right pictures of Figures 4 and 6 depict the force and moment history concerning all of the runs performed with the mesh adaptation algorithm. Similarly, we observe a better convergence behavior with the non-linear SGS method, and more importantly there are less fluctuations in the force and moment prediction at each complexity thanks to the improved iterative convergence.

The skin friction distribution is given for the case using the non-linear SGS method in Figure 5. While the impact the trailing-edge support fairings are visible at all complexities, we note that 5M vertices are needed to see the incipient effect of the leading-edge brackets. However, adapted meshes composed of 20M vertices are required to see all the effects of the leading-edge brackets and their impact on the final predicted lift value. Similarly, it requires at least 10M vertices to compute accurately the flow around the nacelle, before

<figure><img src="image_5.png" /><figcaption>Figure 4. HLPW4 CRM-HL case 1b. Left, convergence of the total lift (top), and moment (bottom) by retaining the final value for each complexity. Right, convergence history of the total lift (top), and moment (bottom) throughout the whole mesh-convergence analysis. For all pictures, results with SGS are in blue and results with the non-linear SGS are in green.</figcaption></figure>
that it is fully separated. Finally, we observe comparable flow patterns for the skin friction between 10M and 20M vertices, for instance they are similar on the flaps and seem to be converged. The main differences are: a more attached flow on the nacelle and at the wing tip at 20M vertices.

The off-body vorticity is depicted in Figure 7 for several mesh sizes. It seems that most of the flow features are already visible at 2.5M vertices, while the main variation between all the complexities is the shape of the vortices. These considerations point out the early-capturing property of the mesh adaptation process, i.e., the main flow features are already captured on coarse adapted meshes.

Two views of the adapted mesh composed of 20M vertices are shown in Figure 8 with the associated Mach number field. We emphasize the high mesh resolution in the boundary layer, shear layer and wake regions, as well as the numerous tip vortices in the wake, due to the brackets, the fairings and the nacelle.

**Remark:** Similar results are obtained at angle of attack 2.78 degrees, the total lift and drag convergence are plotted in Figure 9.

<figure><img src="image_6.png" /><figcaption>Figure 5. HLPW4 CRM-HL case 1b. Evolution of the skin friction on the sequence of adapted meshes.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 6. HLPW4 CRM-HL case 1b. Left, convergence of the total drag (top), pressure drag (middle) and viscous drag (bottom) by retaining the final value for each complexity. Right, convergence history of the total drag (top), pressure drag (middle) and viscous drag (bottom) throughout the whole mesh-convergence analysis. For all pictures, results with SGS are in blue and results with the non-linear SGS are in green.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 7. HLPW4 CRM-HL case 1b. Evolution of the off body vorticity on the sequence of adapted meshes.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 8. HLPW4 CRM-HL case 1. Views of the adapted meshes composed of 20M vertices and its associated Mach number field.</figcaption></figure>

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

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

<figure><img src="image_12.png" /><figcaption>Figure 9. HLPW4 CRM-HL at angle of attack 2.78 degrees. Convergence history of the total lift (left), and total drag (right) throughout the whole mesh-convergence analysis. For all pictures, results with SGS are in blue and results with the non-linear SGS are in green.</figcaption></figure>
## B. Case 2: $C_{L,max}$ study

*$C_{L,max}$ study based on the C-continuation mesh adaptation algorithm*

We did not succeed in applying the C-continuation mesh adaptation algorithm to the $C_{L,max}$ study given the following constraints we set in the introduction:

- only pure tetrahedra adapted meshes are considered (no pseudo structured boundary layer mesh with prisms are allowed)

- no adhoc size prescription in the boundary layer region, i.e., the size prescription should only be provided by the error estimate.

Results at 2.78 and 7.05 degrees are consistent with the TFG Mesh Adaptation Group and some mesh-convergence is obtained, as shown in the previous section. However, at the angle of attack 11.29 degrees, we obtained a flow separation on the nacelle and at the tip of the wing leading to an under prediction of the lift and an over prediction of the moment, see Figure 10 (top left), while some results of the TFG Mesh Adaptation Group show attached flow in these regions. Besides, the flow did not re-attach on finer adapted meshes.

At angle of attack 17.05 and 19.57 degrees, flow patterns were consistent with the TFG Mesh Adaptation Group (similar patterns on the nacelle and at the wing tip) even if a lower lift prediction than the corrected experimental data is predicted. The lower lift and higher moment values at high angle of attack are due to the flow separation at the wing tip which results in the well-known "pizza slice" shape which seems to be unphysical, see Figure 10 (middle and bottom left). However, we were not able to achieve similar simulations at 20.55 and 21.47 degrees because of adjoint state convergence issues.

The polar associated with these results is plotted with the blue curve in Figure 13. The previous comments can be clearly drawn from the polar curves.

*$C_{L,max}$ study based on the α-continuation mesh adaptation algorithm*

Consequently, we chose to perform the $C_{L,max}$ study using the α-continuation mesh adaptation algorithm presented in Section II where the angle of attack is incremented by half a degree each time. For each angle of attack, we performed two runs (i.e., two mesh adaptations and two flow simulations). From HLPW3 results and the convergence study at 7.05 degrees, we chose to perform the $C_{L,max}$ study for adapted meshes composed of 10M vertices which should be a reasonable mesh size for this configuration. Two studies were done.

**FIRST STUDY.** For the first study, we started from the 10M vertices adapted mesh and solution obtained at an angle of attack of 14.29 degrees obtained with the C-continuation mesh adaptation algorithm. For this initial solution, the flow is already separated on the nacelle and at the wing tip with two "pizza slice" patterns. The skin friction of the starting solution is shown in Figure 11 (top left). Using this modified mesh adaptation algorithm, we were able to run up to 24.29 degrees and make the adjoint problem to converge.

We first remark that the vortices generated by the leading-edge brackets appear quickly with the α-continuation method while they are more difficult to obtain with the C-continuation approach.

The second remark is that, as we start from a lower lift and a higher drag/moment solution branch because the flow is separated above the nacelle and at the wing tip, the process stays on that lower lift and higher drag/moment solution branch. Indeed, between 14 and 18 degrees, the flow has almost the same features. At 20.29 degrees, we note that the flow separation increases at the wing tip and, then, the behavior stays the same at higher angle of attack. The lift drop is due to the large flow separation at the wing-body junction that can be observed at 22.29 degrees, see Figure 11 (bottom left). In fact, this flow separation occurs between 20.79 and 21.29 degrees as one can see in Figure 12. As between 20.79 and 21.29 degrees, the flow changes drastically in this region, we note the solution is a bit noisy at 21.29 degrees pointing out a lack of iterative convergence. At that angle of attack, we should have done more than two mesh adaptations and two flow simulations to obtain a more converged couple mesh-solution in the wing-body junction region. Then, at angle of attack above 21.29 degrees, we observe the growth of the wing-body junction flow separation.
<figure><img src="image_13.png" /><figcaption>Figure 10. HLPW4 CRM-HL at angle of attack 11.29 (top), 17.05 (middle), and 19.57 (bottom) degrees. Left, skin friction obtained for an adapted mesh of 10M vertices using the *C*-continuation mesh adaptation algorithm. Right, skin friction obtained for an adapted mesh of 10M vertices using the *α*-continuation mesh adaptation algorithm starting at angle of attack 14 degrees.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Figure 11. HLPW4 CRM-HL α-continuation. Skin friction obtained for adapted meshes of 10M vertices using the α-continuation mesh adaptation algorithm starting at angle of attack 14.29 degrees. From top to bottom and from left to right, angle of attack from 14.29 to 24.29 degrees, each 2 degrees.</figcaption></figure>
The polar curves associated with this study are represented by the red curves in Figure 13. A lower lift and higher drag/moment branch with respect to the experimental data is obtained and match pretty well with the results obtained with the *C*-continuation mesh adaptation algorithm (blue curves). Indeed, if we compare both solutions at 17.05 and 19.57 degrees, we note that nearly identical solutions are obtained, see Figure 10 (middle and bottom). Interestingly, the lift, drag and moment discontinuities occur at the same angle of attack as the experimental data pointing out the correct behavior of the mesh-adaptive solution process in predicting $C_{L,max}$. Thus, the global behavior of the solution is obtained even if quantitatively the force and moment values do not match the experiment. We thus have to figure out if it is a modeling and/or a mesh resolution and/or an iterative convergence problems.

**SECOND STUDY.** As it is not clear whether the wing tip and nacelle flow separations are due to the turbulence model or to a lack of mesh resolution or lack of iterative convergence leading to a lower lift branch, we performed a second study where we started from the 10M vertices adapted mesh and solution obtained at an angle of attack of 7.05 degrees obtained with the *C*-continuation mesh adaptation algorithm, knowing that this solution has no separation. At present, this study has been carried out up to 13 degrees and is still ongoing. It is interesting to compare the solutions obtained with the *C*-continuation and the *α*-continuation mesh adaptation algorithms at 11.29 degrees, see Figure 10 (top). Using the *α*-continuation mesh adaptation algorithm, we obtain a solution at 11.29 degrees where no separation occurs at the wing tip and on the nacelle (top pictures) leading to a higher lift and a lower drag/moment solution meaning that, at that angle of attack, the flow separation was due to the numerics and not to the turbulence model. The polar curves associated with this study are represented by the green curves in Figure 13. For these angle of attacks, we note that the drag value fits the experiment, and that the lift and moment values are closer to the experimental data than the ones at higher angles of attack. Thus, a higher lift and lower drag/moment branch has been obtained. However, we are not able to conclude for higher angles of attack until the end of the study.

<figure><img src="image_15.png" /><figcaption>**Figure 12.** HLPW4 CRM-HL lift drop. Skin friction obtained for adapted meshes of 10M vertices using the *α*-continuation mesh adaptation algorithm starting at angle of attack 14.29 degrees. Left, at an angle of attack of degrees 20.79 before the lift drop and, right, at an angle of attack of degrees 21.29 after the lift drop.</figcaption></figure>
<figure><img src="image_16.png" /></figure>

HLPW4 - CRM-HL - M = 0.2 - Rey = 5.49e6

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

HLPW4 - CRM-HL - M = 0.2 - Rey = 5.49e6

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

HLPW4 - CRM-HL - M = 0.2 - Rey = 5.49e6

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

HLPW4 - CRM-HL - M = 0.2 - Rey = 5.49e6

<figure><img src="image_20.png" /><figcaption>Figure 13. HLPW4 CRM-HL case 2. Polar curves obtained with the C-continuation mesh adaptation algorithm (blue curve) and the α-continuation mesh adaptation algorithm (green curve when starting at 7 degrees and red curve when starting at 14 degrees). For all cases, we plotted the polar obtained with adapted meshes composed of 10M vertices. Corrected experimental data are represented by the black curve.</figcaption></figure>
V. Conclusions

The studies made in this work lead to several conclusions.

Thanks to mesh adaptation, we are able to show mesh-convergence in the solution at low angles of
attack (2.78 and 7.05 degrees) despite the geometric complexity of the configuration. The convergence was
observed on the lift, moment and pressure component of the drag. The viscous component of the drag was
not mesh-converged at 20M vertices.

We have observed at medium angle of attack that the classical mesh adaptation algorithm (*C*-continuation) provides a lower lift solution which is certainly due to the mesh resolution or the iterative convergence as a higher lift solution was recovered with the *α*-continuation mesh adaptation algorithm. In our opinion, improvement of the error estimates should be done to increase the mesh resolution in these flow separation regions in order to recover the correct solution. This is mandatory to make the mesh-adaptive solution process more robust for that kind of applications.

However, it is more difficult to draw conclusions at higher angle of attack. Indeed, as other participants of the Mesh Adaptation Technology Focus Group obtained very similar flow patterns at the wing tip with adapted meshes composed of 50M vertices, these "pizza slices" patterns at the wing tip may be due to the turbulence modeling and being the true numerical solution. On the other side, as the physics becomes richer, we may need larger size meshes to recover a higher lift solution while having similar flow patterns.

It is very interesting to see that, with only 10M vertices, we were able to obtain the correct behavior of
the polar and that the main feature of the flow were represented in the solution (even if quantitatively we
got a lower lift prediction). Moreover, these solutions match with other participants of the Mesh Adaptation
Technology Focus Group pointing out the drastic reduction of the dependency of the solution to the mesh.
The α-continuation mesh adaptation algorithm worked perfectly with increments of half a degree and only
two runs at each angle of attack. Note that, when discontinuities occurs during the polar, more runs seem
to be needed to properly converge the couple mesh-solution.

A last point of interest, we have observed better (higher) lift prediction using the 2nd order convection
Spalart-Allmaras turbulence model on the 2D CRM-HL case when running the turbulence model verification
study. We may expect similar conclusion in 3D.

References

<footnote n="1">¹F. Alauzet. A parallel matrix-free conservative solution interpolation on unstructured tetrahedral meshes. *Comput. Methods Appl. Mech. Engrg.*, 299:116–142, 2016.</footnote>

<footnote n="2">²F. Alauzet, F. Clerici, A. Loseille, C. Tarsia-Morisco, and J. Vanharen. Some progress on CFD high lift prediction using metric-based anisotropic mesh adaptation. In *AIAA SCITECH 2022 Forum*, AIAA Paper 2022-0388, Jan 2022.</footnote>

<footnote n="3">³F. Alauzet and L. Frazza. Feature-based and goal-oriented anisotropic mesh adaptation for RANS applications in aero-nautic and aerospace. *J. Comp. Phys.*, 439(110340):77 pages, 2021.</footnote>

<footnote n="4">⁴S.R. Allmaras, J.E. Bussoletti, C.L. Hilmes, F.T. Johnson, R.G. Melvin, E.N. Tinoco, V. Venkatakrishnan, L.B. Wigton, and D.P. Young. Algorithm Issues and Challenges Associated with the Development of Robust CFD Codes. In *Variational Analysis and Aerospace Engineering*, Springer Optimization and Its Applications book series, chapter 1, pages 1–19. Springer, Dordrecht, Heidelberg, London, New York, 2019.</footnote>

<footnote n="5">⁵A. Belme, F. Alauzet, and A. Dervieux. An a priori anisotropic goal-oriented estimate for viscous compressible flow and application to mesh adaptation. *J. Comp. Phys.*, 376:1051–1088, 2019.</footnote>

<footnote n="6">⁶N.K. Burgess and R.S. Glasby. Advances in numerical methods for CREATE-AV analysis tools. In *52th AIAA Aerospace Sciences Meeting*, AIAA Paper 2014-0417, National Harbor, MD, USA, Jan 2014.</footnote>

<footnote n="7">⁷J. Dompierre, M.-G. Vallet, Y. Bourgault, M. Fortin, and W.G. Habashi. Anisotropic mesh adaptation: towards user-independent, mesh-independent and solver-independent CFD. Part III: Unstructured meshes. *Int. J. Numer. Meth. Fluids*, 39:675–702, 2002.</footnote>

<footnote n="8">⁸W.G. Habashi, J. Dompierre, Y. Bourgault, D. Ait-Ali-Yahia, M. Fortin, and M.-G. Vallet. Anisotropic mesh adaptation: towards user-independent, mesh-independent and solver-independent CFD. Part I: General principles. *Int. J. Numer. Meth. Fluids*, 32(6):725–744, 2000.</footnote>

<footnote n="9">⁹R. Haimes and M. Drela. On the construction of aircraft conceptual geometry for high-fidelity analysis and design. *50th AIAA Aerospace Sciences Meeting*, Jan 2012.</footnote>

<footnote n="10">¹⁰D.S. Kamenetskiy, J.E. Bussoletti, C.L. Hilmes, V. Venkatakrishnan, and L.B. Wigton. Numerical evidence of multiple solutions for the Reynolds-Averaged Navier-Stokes equations. *AIAA Journal*, 52(8):1686–1698, 2014.</footnote>

<footnote n="11">¹¹A. Loseille and F. Alauzet. Continuous mesh framework. Part I: well-posed continuous interpolation error. *SIAM J. Numer. Anal.*, 49(1):38–60, 2011.</footnote>

<footnote n="12">¹²A. Loseille and F. Alauzet. Continuous mesh framework. Part II: validations and applications. *SIAM J. Numer. Anal.*, 49(1):61–86, 2011.</footnote>

<footnote n="13">¹³A. Loseille, F. Alauzet, and V. Menier. Unique cavity-based operator and hierarchical domain partitioning for fast parallel generation of anisotropic meshes. *Comput. Aided Des.*, 85:53–67, 2017.</footnote>

<footnote n="14">¹⁴A. Loseille, A. Dervieux, and F. Alauzet. Fully anisotropic goal-oriented mesh adaptation for 3D steady Euler equations. *J. Comp. Phys.*, 229:2866–2897, 2010.</footnote>
<footnote n="15"><sup>15</sup>A. Loseille and R. Löhner. Cavity-based operators for mesh adaptation. In *51th AIAA Aerospace Sciences Meeting*, AIAA Paper 2013-0152, Dallas, TX, USA, Jan 2013.</footnote>

<footnote n="16"><sup>16</sup>T. Michal, D.S. Kamenetskiy, D. Marcum, F. Alauzet, L. Frazza, and A. Loseille. Comparing anisotropic error estimates for ONERA M6 wing RANS simulations. In *56th AIAA Aerospace Sciences Meeting*, AIAA Paper 2018-0920, Kissimmee, FL, USA, Jan 2018.</footnote>

<footnote n="17"><sup>17</sup>M.J. Pandya, B. Diskin, J.L. Thomas, and N.T. Frink. Improved convergence and robustness of USM3D solutions on mixed element grids. *AIAA Journal*, 54(9):2589-2610, 2016.</footnote>