text:
Vincent Liguori
email: vincent.liguori@onera.fr

Frederic Moens
email: frederic.moens@onera.fr

Jacques Peter
email: jacques.peter@onera.fr

Eric Laurendeau
email: eric.laurendeau@polymtl.ca

AERO2023-81-LIGUORI

The implementation of the Non-Linear Vortex Lattice Method in the CHAMPS framework is presented. Developments are made with the aim of improving the prediction of pressure distributions in the region near the root of the wing with physics-based modifications. The viscous coupling algorithm is modified to better take into account compressibility effects on swept wings. A methodology to better capture the spanwise variation of the flow caused by the wing-body geometry while computing the database used for the viscous coupling is tested. Results for the isolated modifications are shown. The results of the complete methodology applied on classical test cases of the DLR F4 wing-body and the High-Speed Common Research Model wing-body are shown. Improvements in shock prediction are observed however limitations associated to drag overprediction are identified.

INTRODUCTION

Potential methods based on a distribution of discrete perturbations to represent the geometry of a body are useful tools in the context of early aircraft design phases because of their low computational cost compared with other higher fidelity methods such as Reynolds-Averaged Navier-Stokes (RANS) equations solver [START_REF] Taylor | Low-fidelity method for rapid aerostructural optimisation and design-space exploration of planar wings[END_REF]. They can in particular model adequately the spanwise distribution of the induced angle of attack on a wing caused by a lift distribution. These methods in their original forms are limited by their underlying assumptions of inviscid, incompressible and irrotational flow. These limitations can be partially overcome by adding corrections to the original potential equations to take into account the phenomena that would otherwise be neglected. Many such approaches have been tested in the literature. Notably, to consider some compressibility effects up to the transonic regime the Prandtl-Glauert transformation can be used. To consider additional phenomena such as viscosity and shocks, coupling the potential method with higher fidelity data is a popular method that has been applied by numerous authors with various methodologies, as compiled in the literature [START_REF] ¸ugar Gabor | Fast and accurate quasi-3D aerodynamic methods for aircraft conceptual design studies[END_REF].

The work presented in this paper uses a methodology of this type, based on coupling a Vortex Lattice Method (VLM) with two-dimensional viscous data computed from a RANS solver to correct the local lift and drag coefficients using strip-theory. This method is referred to as Non-Linear Vortex Lattice Method (NL-VLM) and is implemented in the CHAMPS framework developed at École Polytechnique de Montréal [START_REF] Parenteau | Development of Parallel CFD Applications with the Chapel Programming Language[END_REF] [START_REF] Liguori | Validation of a nonlinear Vortex Lattice Method for multiple wing sweep angle configurations[END_REF].

While the method has demonstrated a good agreement with higher fidelity solvers and experimental data in lowspeed and transonic conditions, the prediction of shocks in transonic conditions near the fuselage in the root region of the wing has been identified as a limitation of the method. In one previous work [START_REF] Gallay | Preliminary-Design Aerodynamic Model for Complex Configurations Using Lifting-Line Coupling Algorithm[END_REF], the method was seen predicting shocks in this region when a gradual pressure decrease was observed experimentally. In another previous work [START_REF] Kontogiannis | Viscous-Inviscid Analysis of Transonic Swept Wings using 2.5D RANS and Parametric Shapes[END_REF], the method predicted the appearance of large shocks that caused stalling of the wing sections near the fuselage in transonic conditions, which did not match with numerical results computed using a more computationally intensive 3D RANS method.

Based on this observation, modifications to the implementation of NL-VLM were tested in this paper with the goal of improving the prediction of shocks on wing-body configurations in the wing region near the fuselage in transonic flow conditions.

The numerical implementation of the NL-VLM algorithm without the modifications is first summarized. A more detailed description of the modifications to the methodology compared with previous implementations is then presented. Numerical results obtained on the wingbody configurations of the DLR-F4 and Common Research Model geometries in transonic conditions are finally shown. The effects of the modifications on these test cases are then discussed.

INITIAL NON-LINEAR VORTEX LATTICE METHOD IMPLEMENTA-TION

The NL-VLM implementation is based on a potential three-dimensional method, a two-dimensional higher fidelity method to create a viscous database, a coupling algorithm to correct the potential method with the viscous database and a methodology to compute the total forces over the body from the coupled system. The initial implementation of these main steps is presented in this section.

The NL-VLM in the CHAMPS framework is already described in the literature [START_REF] Liguori | Validation of a nonlinear Vortex Lattice Method for multiple wing sweep angle configurations[END_REF], as such unchanged parts of the methodology are briefly summarized to the extent required for the subsequent description of the modifications applied in the context of this work in Sec. 3.

Three dimensional potential method

To model the three-dimensional geometry, a potential method with hybrid elements is used. The geometry is considered to be an assembly of two categories of bodies. The first category is made of lifting surfaces, which are considered to be thin bodies. In the case of this work, this category is used for wing surfaces. These bodies are modeled using a classical vortex-lattice method [START_REF] Katz | Low-Speed Aerodynamics[END_REF]. A structured grid is placed over the planform of the wing, simplified as being flat between the leading edge and trailing edge. Vortex ring elements are then placed at the quarter length of the structured grid elements.

The second type of bodies are non-lifting bodies, which are considered to be bluff bodies. In the case of this work, this category is used for fuselages. These bodies are modeled using an unstructured constant source panels method, the panels being placed over the surface of the bodies.

For both types of elements, a Neumann boundary condition of velocity tangent to the element is imposed. This condition is verified at a control point placed at the center of the element in the case of the source panels and at 3/4 of the length of the elements for the vortex rings, which allows the vortex rings to create lift with a slope corresponding to thin airfoil theory [START_REF] Deyoung | Historical evolution of vortex-lattice methods[END_REF] [START_REF] Pistolesi | Considerations on the Mutual Interference of Aerofoil Systems[END_REF]. The complete system considering both vortex ring and source panel elements is solved to compute the intensity of the perturbations. To improve the accuracy of the model in transonic conditions, the Prandtl-Glauert transformation is applied to both types of elements [START_REF] Drela | Flight Vehicle Aerodynamics[END_REF].

Viscous database

The viscous database consists of tables of lift, drag and moment coefficients of an airfoil for a range of angles of attack. These coefficients are obtained with the CHAMPS RANS equations solver, which is a solver using a cell-centered finite-volume approach. A more detailed description of the solver and its capabilities is outside the scope of this work but is available in the literature [START_REF] Parenteau | Development of Parallel CFD Applications with the Chapel Programming Language[END_REF]. These tables of coefficients are computed for airfoils representing chosen sections along the span of the wing. The database is then interpolated along the wing span so that a corresponding local viscous database is allocated to every wing strip.

Since wing sweep has an important effect on flow fields around airfoils, it has been shown that simply using a 2D RANS simulation is insufficient to obtain solutions that match experimental results both in force coefficients and pressure distribution [START_REF] Gallay | Preliminary-Design Aerodynamic Model for Complex Configurations Using Lifting-Line Coupling Algorithm[END_REF]. Wing sweep effects are considered in the viscous database calculation through a 2.5D approach to consider an infinitely swept wing of constant geometry instead of a straight wing. This is done using a "sheared-mesh" approach in which a 2D mesh is extruded in a third dimension over 1 cell. The extrusion direction is aligned with the quarter-chord sweep angle of the wing. A periodic boundary condition is then applied on the side faces of the mesh created during the extrusion. An alternative is to directly implement the 2.5D hypothesis in the RANS equations, removing the need for halocells [START_REF] Ghasemi | A Two-Dimensional/Infinite Swept Wing Navier-Stokes Solver[END_REF]. The former method can be viewed as a discrete form, the latter as a continuous form. In this paper, the discrete form is used.

Viscous Coupling

The method uses a coupling between the inviscid method and a viscous database to include non-linear effects. This coupling is done only on lifting elements by forcing the chordwise sections of the wing in the potential flow to match the lift of their corresponding viscous database at their effective angle of attack. To enforce this, the algorithm used is inspired by Van Dam's α-coupling method [START_REF] Van Dam | On the aerodynamic design of multi-element high-lift systems[END_REF], which consists in rotating the normal vector of the elements by a viscous coupling correction angle to influence the vortex intensity of the elements. The coupling methodology, similar to the loosely coupled algo-rithm used by Gallay and Laurendeau [START_REF] Gallay | Preliminary-Design Aerodynamic Model for Complex Configurations Using Lifting-Line Coupling Algorithm[END_REF], is summarized as 1. Compute the sectional lift coefficient of the local strip in the inviscid potential flow using the Kutta-Joukowski theorem [START_REF] Katz | Low-Speed Aerodynamics[END_REF] 2. Compute the local effective angle of attack α e from the theoretical lift slope of the wing section.

3. Obtain through interpolation the lift coefficient at the computed effective angle of attack from the viscous database.

4. Update a viscous coupling correction angle based on the difference between the two lift coefficients.

5. Apply a rotation to the normal of the elements in the section corresponding to the viscous coupling correction angle.

These steps are repeated iteratively for each chordwise strip of elements until the difference between the viscous and inviscid lift coefficients is sufficiently small, at which point the update to the viscous correction angle will tend toward zero and convergence is considered achieved.

Force Coefficients Computation

The lift and moment coefficients are computed from the potential solution using the Kutta-Joukowski theorem. The drag coefficient is computed in two parts. The induced drag is computed using a Trefftz plane [START_REF] Katz | Low-Speed Aerodynamics[END_REF] and the viscous and pressure components of the drag are computed by integrating the corresponding drag coefficients at the effective angle of attack of the viscous database for each wing section. Both contributions are summed to obtain the total drag. A more detailed description of the methodology is available in a previous presentation of the method [START_REF] Liguori | Validation of a nonlinear Vortex Lattice Method for multiple wing sweep angle configurations[END_REF].

MODIFICATIONS TO THE NL-VLM IMPLEMENTATION

The three modifications applied to the NL-VLM are the use of a new equation to compute the local effective angle of attack, taking into account the taper of the wing in the viscous database computation and adding a new boundary condition in the computation of the viscous database to model the perturbation effect of the fuselage.

Calculation of the local effective angle of attack

The coupling algorithm presented in section 2.3 computes the local effective angle of attack by dividing the inviscid lift coefficient of the section by a theoretical local lift slope value.

The collocation points in the VLM method are placed to ensure a local slope of 2π in the case of an infinite unswept wing with no other transformation being applied in accordance with the 1/4-3/4 rule of Pistolesi [START_REF] Pistolesi | Considerations on the Mutual Interference of Aerofoil Systems[END_REF]. The Prandtl-Glauert correction and the presence of sweep however have an influence on the lift slope of the elements.

In a previous work [START_REF] Parenteau | VLM Coupled with 2.5D RANS Sectional Data for High-Lift Design[END_REF] these effects are considered with two distinct multiplying factors, so that

C Lα,inv = 2πcos(Φ) (β ) (1)
β = 1 -M 2 (2)
where Φ is the sweep angle of the quarter-chord of the wing, β is the Prandtl-Glauert transformation factor and M is the Mach number of the freestream. In another previous work the lift slope used is taken for each section directly from the linear part of the lift slope in the viscous database [START_REF] Kontogiannis | Adjoint State of Nonlinear Vortex-Lattice Method for Aerodynamic Design and Control[END_REF]. In this work, the lift slope used is the theoretical lift slope of an infinite swept wing in a potential field where the Prandtl-Glauert transformation has been applied, derived in the literature [START_REF] Drela | Flight Vehicle Aerodynamics[END_REF] as

C Lα,inv = 2πcos(Φ) β 2 cos 2 (Φ) + sin 2 (Φ) . (3)
In the two limiting cases of an unswept wing and a Mach number of zero, Eq. 1 and Eq. 3 are equivalent. They however differ when both effects are present simultaneously. To validate this approach, the VLM code without viscous coupling is used on the case of a swept wing with an aspect ratio of 1000, which is close to an infinite wing. The sweep angle used is 26.57 degrees and the Mach number used for the Prandlt-Glauert transformation is 0.85.

It can be seen in Fig. 1 that the new theoretical lift slope is closer to the numerical results than the previous method. Using the approach of estimating the effective angle of attack by assuming that the local lift slope corresponds to the infinite swept wing case with the local section geometry, the new method is considered more accurate than the previous formula to compute the effective angle of attack. In contrast with the method of using the linear section of the viscous slope shown in the literature, this method is only based on inviscid flow and as such is considered more representative of the numerical method. It also does not require the delicate step of identifying the linear part of the viscous slope for every wing section.

Consideration of tapered geometry in viscous database

As mentioned in Sec. 2.2, a 2.5D sheared-cell approach is used in the initial method. This deformation of the mesh only considered the quarter-chord sweep of the wing section. However, the difference in sweep between the leading edge and the trailing edge in a tapered wing could lead to differences in flow solutions, particularly near the stall angle.

The region of the wing near the fuselage is usually more tapered than the outboard sections of the wing for transonic transport aicraft, this difference being visually identifiable at a kink in the wing trailing edge. For this reason, a new mesh deformation method was tested to model the taper effect during the viscous database computation for sections in this region. This new method consists in varying the extrusion direction of the 2D mesh based on the location relative to the airfoil.

The mesh points located upstream of the leading edge of the airfoil are moved in the direction of the leading edge sweep. The extrusion direction of the mesh points located between the leading edge and trailing edge varies linearly between the leading edge sweep and trailing edge sweep so that the taper ratio is applied to the extruded mesh. Finally, the points located downstream of the trailing edge are moved in the direction of the trailing edge sweep. To keep the airfoil section consistent, a vertical scaling corresponding to the taper ratio is also applied to the mesh points located between the lowest and highest points of the airfoil. A schematic representation of this extrusion method from a top-view can be seen in Fig. 2. The airfoil is shown with continuous lines and the mesh extrusion direction is shown in dashed lines. Using this methodology, the taper effects are considered di-rectly by the physics of the solver similarly to recent 2.5D approaches [START_REF] Ghasemi | A Two-Dimensional/Infinite Swept Wing Navier-Stokes Solver[END_REF][3], in contrast with other approaches for taper effects in the literature based on subsequent correction of non-tapered results [START_REF] Xu | An improved 2.75D method relating pressure distributions of 2D airfoils and 3D wings[END_REF]. This varying deformation has the effect of creating a difference between the geometries at the end of the extrusion. In the limiting case of a large enough extrusion, mesh cells in the airfoil region would eventually cross and create negative volumes. Even before this limiting case, a large variation of chord could have a non-negligible impact on the average Reynolds number of the flow over the airfoil. To avoid these effects, a parametric study is done on the extrusion length of the airfoil and its impact on the obtained aerodynamic coefficients. The results, shown in Fig. 3, indicate that the results tend to converge when a sufficiently small extrusion length is used. Based on these results the extrusion length used in this work for a tapered geometry is fixed at five percent of the chord.

To validate the hypothesis that this modification allowed the solver to capture physical phenomena caused by the tapered geometry, a comparison is made between the surface solution obtained with the tapered extrusion and the original 2.5D method. The test case corresponds to the airfoil at the wing root of the High-Speed Common Research Model (HS-CRM) geometry [START_REF] Vassberg | Development of a Common Research Model for Applied CFD Validation Studies[END_REF], used notably in the fifth AIAA drag prediction workshop [START_REF] Levy | Summary of Data from the Fifth Computational Fluid Dynamics Drag Prediction Workshop[END_REF]. The flow conditions are those of the first test case of the workshop, namely a Mach number of 0.85 and a global Reynolds number of 5 000 000 which corresponds to a Reynolds number of 8 461 600 relative to the local chord since the global Reynolds number uses the mean aerodynamic chord as a reference length.

A top-view of the wall shear stress lines in Fig. 4 and Fig. 5 before stall for each configuration reveals that there is higher convection by a crossflow in the 2.5D case compared with the 2.5D tapered case near the trailing edge

φ LE =37.1 • , φ T E =9.8 • )
of the airfoil. This behavior is expected and shows that the method is able to capture the difference caused by the geometry. A comparison of the lift coefficients in Fig. 6 shows the impact of this difference on the stall behavior. Notably, a steeper loss of lift is observed in the 2.5D tapered case after the onset of stall. It is also noted during testing that this mesh modification significantly hampers the convergence of the solver. In some instances, stabilization techniques such as Selective Frequency Damping [START_REF] Liguori | Implementation of an efficient Selective Frequency Damping method in a RANS solver[END_REF] are used to obtain converged solutions.

Consideration of fuselage effect on crossflow in viscous database

The last modification aims to model the fact that, near the root of the wing, the presence of the fuselage reduces the crossflow over the wing. In the initial implementation, this variation in crossflow is not considered in the computation of the viscous database since the 2.5D approach assumes an infinite swept wing and as such no spanwise variation of crossflow.

To model this variation of crossflow, the boundary conditions applied on the extruded 2.5D mesh are modified. Instead of a periodic boundary condition on the faces of the mesh in the extrusion direction, a symmetry condition is applied on the face at the inboard side of the extruded airfoil. A zero-order extrapolation boundary condition is applied to the outboard side face so that the values of each conservative variable are copied in their corresponding halo cell. This approach creates a spanwise variation in the crossflow between the two extrusion faces of the mesh.

This approach is tested on a near-root airfoil section of the DLR-F4 wing body. This geometry was used in the first drag prediction workshop [START_REF] Levy | Summary of data from the first AIAA CFD Drag Prediction Workshop[END_REF] and has also already been tested in a previous work with a NL-VLM approach using a RANS 2.5D database [START_REF] Gallay | Preliminary-Design Aerodynamic Model for Complex Configurations Using Lifting-Line Coupling Algorithm[END_REF]. The flow conditions correspond to the first test case of the workshop, which are a lift coefficient of 0.5, a Reynolds number of 3 000 000 using the mean aerodynamic chord as the reference length and a Mach number of 0.75. In this previous work the shock prediction in the section with pressure measurement nearest to the fuselage did not match the experimental results.

The test case is first solved with the initial NL-VLM implementation to obtain a spanwise lift distribution. Four 2.5D polars are used to create the original viscous database, using the root section, the kink section with two different sweeps corresponding to the wing sweep before and after the kink and the wing tip section. The airfoil section at 18.5% of the span away from the symmetry plane is then extracted and the local lift coefficient obtained through the initial NL-VLM solver is used as a reference for the section. This section is selected because it is the nearest section to the fuselage for which pressure measurements were taken in experimental testing.

The section is then tested with the CHAMPS RANS solver in its local flow conditions using the original 2.5D approach and the modified boundary condition approach, each time matching the same section lift coefficient. In both cases, an extrusion distance of 1 chord is used. The goal of this test is to investigate if the modification of the boundary condition improves the agreement between the pressure coefficients of the RANS solver and the experimental result. The pressure coefficients distribution obtained in Fig. 7 shows that the modified boundary conditions approach predicts a distribution closer to the experimental results. In particular, the shock over the upper surface predicted in the original 2.5D approach is not present with the new approach.

Since the new boundary conditions create a variation in crossflow between the two extrusion faces of the mesh, it is expected that varying the extrusion length would have an impact on the result. This is confirmed in Fig. 8 and Fig. 9, showing a parametric study of the extrusion length when using the one-sided symmetry approach. This parametric study is on the same HS-CRM root airfoil test case shown in Sec. 3.2. The same results with the original 2.5D approach are also shown for comparison. The aero-Figure 7: Pressure coefficients comparison of the 2.5D and 2.5D with symmetry condition approach at equivalent lift with experimental results (Ma=0.75, Re=3.92e6, φ c/4 =21.0 • ) dynamic coefficients tend toward the 2.5D results when the extrusion length becomes large. This is expected and reflects the fact that since the extrusion is done over only one cell in which the solver assumes a linear variation and that one face imposes a null crossflow, increasing the distance between the two faces reduces the spanwise component of the crossflow gradient since the variation is performed over a larger distance. In the limiting case of an infinite-length extrusion, the results would converge toward the original 2.5D case. When applying the one-sided symmetry methodology, two behaviors are observed. The first one is an increase in C Lmax and of the stall angle of attack. The second one is a global increase in drag. The extrusion length impact is problematic for two rea- sons. The first one is that it was observed in the tapered geometry database investigation in Sec. 3.2 that the extrusion length has an impact on the results and that a small extrusion length relative to the chord is desirable. The second reason is practicality, since modifying a mesh is a more time-demanding procedure compared to varying a numerical coefficient, particularly if multiple values must be tested for a parametric study. To avoid these difficulties, a new boundary condition is implemented in the RANS solver. For this boundary condition, the density and pressure of the interior cell are copied in the halo cell, while the velocity is set so that a predetermined percentage of the normal velocity through the face of the boundary is canceled. The implementation is summarized as:

ρ halo = ρ in p halo = p in v x,halo = v x,in -2k(⃗ v in •⃗ n) • n x v y,halo = v y,in -2k(⃗ v in •⃗ n) • n y v z,halo = v z,in -2k(⃗ v in •⃗ n) • n z (4)
where ρ is the density, p is the pressure, ⃗ v is the velocity at the cell center and ⃗ n is the unit normal vector for the boundary face. The halo cell and the interior cell are identified with their respective indices and the x,y and z indices identify the directional components of vectors. The parameter k is the normal forcing factor and is set by the user. In the limiting case where k = 1.0, the boundary condition is equivalent to a symmetry condition. In the limiting case where k = 0.0, the boundary becomes a zero-order extrapolation like the periodicity condition used in the original 2.5D approach.

A parametric study of the variation of this normal forc-ing factor in Fig. 10 with a constant extrusion length of 5 percent of the airfoil chord shows that as expected, the behavior of the aerodynamic coefficients follows the same trend when increasing the normal forcing factor as when reducing the extrusion length with the root symmetry condition in Fig. 9. Based on these results the new normal forcing boundary condition with an extrusion length of 5 percent of the chord is used to generate the viscous databases in the following results.

NUMERICAL RESULTS ON WING-BODY CONFIGURATIONS

To analyze the effect of the new modifications on the complete NL-VLM solver, the NL-VLM with the modifications presented in Sec. 3 is used on wing-body configurations. The same test cases are also run with the NL-VLM implementation without the modifications to the viscous database for comparison. It should be noted that the modification to the inviscid slope used in the effective angle computation and presented in Sec. 3.1 is used in both cases. This allows the analysis to focus on the effect of the viscous database modification, the improvement in effective angle computation by the slope modification already having been showcased.

The first considered test case is case 1 of the first drag prediction workshop, which consists in enforcing a lift coefficient of 0.5 for the DLR-F4 wing-body geometry [START_REF] Levy | Summary of data from the first AIAA CFD Drag Prediction Workshop[END_REF]. The flow conditions are a Reynolds number of 3 000 000 and a Mach number of 0.75, which are transonic conditions. For the NL-VLM computation, the viscous database is assembled from 4 RANS polars.

One is computed at the root of the wing, two are computed on each side of the kink of the wing to account for the sweep variation and one is computed for the tip airfoil section. For the new implementation, a normal forcing factor of 0.05 is used for the inboard normal forcing boundary condition in the root section database computation, based on a parametric study to identify an adequate value to avoid an increase of drag that would be too important. The wing surface solution is generated from the NL-VLM results through an interpolation of the surface solution of the RANS results used to assemble the database. The pressure coefficients of these reconstructed solutions are compared to the experimental pressure measurements from the ONERA wind tunnel for the first drag prediction workshop in Fig. 11 and Fig. 12. The sections considered are the two sections closest to the fuselage available. The results show that the new implementation avoids overpredicting a shock near the wing root. This overprediction of the shock by the unmodified method is consistent with results from a similar previous work [START_REF] Gallay | Preliminary-Design Aerodynamic Model for Complex Configurations Using Lifting-Line Coupling Algorithm[END_REF]. Some discrepancies with experimental measurements are still present, notably in the second section measured along the span.

To analyze the effect of the modifications on the aerodynamic coefficients predicted by the NL-VLM, an angle of attack sweep is done on the DLR-F4 geometry in the same flow conditions and compared with experimental results from ONERA. For the drag computation, a constant value of 0.009 corresponding to wind tunnel measurements is added to the NL-VLM C D result to model the fuselage drag since no viscous coupling is applied on the body. The results, in Fig. 13 and Fig. 14, show that while the modifications allow the method to predict a Figure 12: Comparison of the reconstructed solution of the initial NL-VLM and the modified implementation with experimental measurement as reference for the DLR-F4, section at 23.8% of the span higher stall angle and maximum C L which is closer to the experimental results, the drag seems to be overpredicted compared with experimental data. The drag predicted by the initial method is closer to the experimental results for lower angles of attack. As it was observed that the normal forcing modification to the viscous database caused an increase in the drag coefficient, these results seem to indicate that the overprediction of drag is a downside of this modification. A second test case is considered consisting of the high-speed Common Research Model wing-body geometry [START_REF] Vassberg | Development of a Common Research Model for Applied CFD Validation Studies[END_REF]. The first case of the fifth drag prediction workshop is used for the flow conditions, which are a lift coefficient enforced at 0.5, a Reynolds number of 5 000 000 Figure 14: Drag coefficient against angle of attack of the original and modified NL-VLM implementation and experimental data for the DLR-F4 geometry in transonic conditions and a Mach number of 0.85. The CRM wing shows significant airfoil variation along its span, and for this reason, in this case, 12 sections are considered to assemble the viscous database, distributed roughly at every 10% of the span from the root to the tip with the addition of the wing kink for which a section on either side is considered to take the sweep and taper change into account. Due to the instability in the nonlinear RANS solver caused by the tapered geometry modification to the database mentioned in Sec. 3.2, the tapered geometry modification is only applied to the first viscous section. The normal forcing factor used for the RANS boundary condition varies from 0.1 for the first section to 0.04 to the second and 0.01 for the third. The other sections are considered sufficiently far from the fuselage for this modification to be unnecessary.

The pressure coefficients over the reconstructed solution are again compared with experimental results from the National Transonic Facility for the same case for two sections near the fuselage in Fig. 15 and Fig. 16. It is again observed that while discrepancies are present between the NL-VLM and experimental results, the modified method is closer to the experimental pressure coefficients. In particular, the strong shocks predicted by the initial method near the wing root are not present in the modified method. The presence of these shocks in the original method is consistent with results from previous work on the same case, where it was reported that this wrong prediction near the wing root caused an early stalling of the inboard wing section [START_REF] Kontogiannis | Viscous-Inviscid Analysis of Transonic Swept Wings using 2.5D RANS and Parametric Shapes[END_REF].

CONCLUSION

The aim of this paper was the exploration of methods to improve the prediction of shocks by the NL-VLM for transonic conditions in the region of the wing near the fuselage. To achieve this, three modifications to the original implementation were explored and integrated. The results obtained indicate that the modifications allow a more reliable prediction of flow physics. However, some discrepancies with experimental results remain and while the normal forcing correction seems to improve the shock prediction for sections near the fuselage, caution should be taken about the drag increase it causes.

Figure 1 :

 1 Figure 1: Comparison of the initial and modified lift slope with the numerical results of the VLM on an near-infinite swept wing (Ma=0.85, φ =26.57 • )

Figure 2 :

 2 Figure 2: Schematic top view of the extruded mesh with 2D, 2.5D and 2.5D tapered geometries

Figure 3 :

 3 Figure 3: Lift coefficient for various angle of attacks using tapered extrusion method relative to the inverse of the extrusion length (Ma=0.85, Re=8.46e6, φ LE =37.1 • , φ T E =9.8 • )

Figure 4 :Figure 5 :Figure 6 :

 456 Figure 4: Wall shear stress lines and y-direction shear stress contour for the 2.5D case at an angle of attack of 1.0 (Ma=0.85, Re=8.46e6, φ LE =37.1 • , φ T E =9.8 • )

Figure 8 :

 8 Figure 8: Parametric study of extrusion length when using a symmetry condition at root for lift coefficient against angle of attack and comparison with 2.5D approach (Ma=0.85, Re=8.46e6, φ c/4 =31.4 • )

Figure 9 :

 9 Figure 9: Parametric study of extrusion length when using a symmetry condition at root for lift coefficient against drag coefficient and comparison with 2.5D approach (Ma=0.85, Re=8.46e6, φ c/4 =31.4 • )

Figure 10 :

 10 Figure 10: Parametric study of the k normal forcing factor when using the modified boundary condition for lift coefficient against drag coefficient and comparison with 2.5D approach (Ma=0.85, Re=8.46e6, φ c/4 =31.4 • )

Figure 11 :

 11 Figure 11: Comparison of the reconstructed solution of the initial NL-VLM and the modified implementation with experimental measurement as reference for the DLR-F4, section at 18.5% of the span

Figure 13 :

 13 Figure 13: Lift coefficient against angle of attack of the original and modified NL-VLM implementation and experimental data for the DLR-F4 geometry in transonic conditions

Figure 15 :Figure 16 :

 1516 Figure 15: Comparison of the reconstructed solution of the initial NL-VLM and the modified implementation with experimental measurement as reference for the HS-CRM, section at 13.1% of the span

ACKNOWLEDGEMENTS

This work benefits from FRQNT and NSERC doctoral scholarships as well as from the Canada Research Chair program.


Edward Yeboah
email: yeboah5@hotmail.com

Smallholder agriculture: Factors (causes) generating a change

Keywords: Smallholder agriculture, Factors generating a change, Integrated Soil Fertility Management (ISFM), Monitoring and Evaluation; Semi-deciduous forest zone of Ghana

Context

Food security in Sub-Saharan Africa (SSA) is adversely affected by climate change related to decreases in rainfall, nutrient and carbon mining, varying rainfall amount, low use of inorganic fertilizers and limited amount and varying quality of organic resources. This study concerns the smallholder agriculture for which food insecurity, poverty, land degradation, and climate change were the main factors degrading livelihood. Integrated Soil Fertility Management (ISFM) which is the combined application of organic and inorganic fertilizers for sustainable crop production has been proposed to address the inherently low fertility of tropical soils particularly in SSA.

Prospects

Lessons learnt for the future (future evolution scenario) are a change of scale, dissemination, study topic, etc…with aim toup-scaling of technologies, mintaining the study as Long term experiment site and linking smallholder farmers to market

Objectives

The main objective is to evaluate organic resource quality and quantity on sustainable maize crop yields in the semi-deciduous forest zone of Ghana where the soil fertility is decreasing.

Description of the Project

The project started on September 2002 to date. In a farmer participatory on farm trial in the semi-deciduous forest zone of Ghana, six organic resources commonly available to smallholder farmers and ranging from high, medium and low qualities with and without inorganic fertilizers were evaluated. High yielding maize variety "mamaba" was chosen for test crop. The experiment is fully replicated and repeated on another site of contrasting soil type.

The key stakeholders and beneficiaries are smallholder farmers, International Scientific community, International Donor Organizations, Ministry of Food and Agriculture, Extension Agents, Non Governmental Organizations (NGO'S), Universities and Agricultural Colleges as well as District Assemblies.

Results and Impacts

Quantitative gains (production, productivity, income, etc.) Technological options to produce continuous maize grain yield between 3 to 6 t ha -1 for two seasons in year depending on farmer's resource endowment Qualitative gains (living conditions, well-being, etc.) Ability to pay children school fees and registration for National Health Insurance Scheme (NHIS) Sustainability (economic, social, environmental) Carbon sequestration with innovations, food security and ability to support extended family members

Dynamic analysis of factors determining change

Use a time scale (e.g. below) to identify:

-The stage of development of the process described The criteria used by farmers were: availability of organic inputs; labour involved; transportation; bulkiness; yield and soil improvement


domain: spi phys

text:
Hugo Dornier
email: hugo.dornier@onera.fr

Julien Marty
email: julien.marty@onera.fr

Sébastien Bourasseau
email: sebastien.bourasseau@onera.fr

Itham Salah
email: itham.salaheldin@onera.fr

AERO2023-35-DORNIER ROBUST AND EFFICIENT CFD SIMULATIONS OF THE ARL-SL19 SUPERSONIC CASCADE THROUGH ADAPTIVE MESH REFINEMENT

This paper presents an efficient adaptive mesh refinement method able to generate automatically high quality grids that allow accurate CFD simulations for supersonic turbomachinery flow configurations. The chosen test case is the well documented ARL-SL19 [10] supersonic cascade. First, a numerical error analysis and a validation against reference results are performed to quantify the reliability of the numerical method. The next part consists in studying two mesh adaptation approaches based on a shock sensor and another one based on wake and recirculation. The results are compared to those obtained with the one relying on the more global Mach number sensor. After defining a relevant adaptation strategy, a comparison is made of the different approaches performance. Finally, the shock and the turbulent features sensors are used simultaneously to find the best combination ratios.

INTRODUCTION

Supersonic flows in axial compressors and turbines produce complex shock structures between blades. On the one hand, these phenomena can induce strong pressure losses and flow separations. But on the other hand, a supersonic compressor can use faster rotors resulting in higher pressure ratios for each stage than its subsonic counterpart [START_REF] Kantrowitz | The supersonic axial-flow compressor[END_REF]. Consequently, for a given target compression ratio, less stages are necessary which makes the whole system smaller and generally lighter. The industry has been driven by this advantages to design mostly transonic and supersonic compressors. Nevertheless, developing blade geometries for supersonic compressors and turbines remains a challenge because small changes in the upstream flow field can have high impacts on the shock structures. Thus accurately resolving such flows where shocks can interact with each others or with viscous layers is mandatory. A major difficulty is that for complex geometries the shock patterns cannot be known a priori and generating a mesh density able to capture any flow topology possible would be computationally too expensive. According to the literature [START_REF] Alauzet | A decade of progress on anisotropic mesh adaptation for computational fluid dynamics[END_REF][1] [START_REF] Vivarelli | Efficient adjoint-based mesh adaptation applied toturbo-machinery flows[END_REF] [START_REF] Benard | Mesh adaptation for large-eddy simulations in complex geometries[END_REF], mesh adaptation appears to be one of the best strategies to reach both accurate results and computationally acceptable meshes. Hence the effort in the initial mesh generation process can be significantly reduced. This paper details an adaptive mesh refinement method to automatically generate grids in order to improve CFD calculation accuracy for the aforementioned supersonic turbomachinery flows. The adaptation uses existing feature based sensors and criteria gathered from the literature along with new ones conceived during this study. Besides, the mesh modification includes a remeshing phase based on an isotropic insertion of nodes by Cassioppée's [START_REF] Benoit | Cassiopee: A CFD pre-and post-processing tool[END_REF] tool : adaptCells. The performance of the adaptation process has been assessed on the well documented ARL-SL19 [START_REF] Fottner | Test cases for computation of internal flows in aero engine components[END_REF] supersonic cascade. This test case is representative of complex supersonic turbomachinery flow configurations. It features phenomena such as oblique and strong shocks, shock/shock and shock/boundary-layer interactions as well as induced separations and wakes.

In the first part the down-selected methods and test cases are presented. The second part focuses on assessing the reliability of the numerical methods with a numerical error study and a validation against experimental and CFD references. Then the impact of the main parameters of the adaptation process is investigated with a prescribed shock sensor. The resulting know-how is the used to compare the efficiency regarding the other sensors as well as a combination of the shock and turbulence based ones. Finally the main outcomes of our study are summarized and some perspective given for future work.

METHOD

2.1 ARL-SL19 configuration and numerical setup

Test case

The test case for this study is the ARL-SL19 supersonic cascade. This case has been conceived as a reference case for supersonic compressors cascades and has been tested in three different supersonic wind tunnels. The first test campaign was made by Detroit Diesel Allison (DDA) [START_REF] Fleeter | Experimental Investigation of a Supersonic Compressor Cascade[END_REF] in 1974 while the test campaigns at ONERA [START_REF] Fourmaux | Test results and theoretical investigation on the ARL 19 supersonic blade cascade[END_REF] and DFVLR [START_REF] Tweedt | Experimental investigation of the performance of a supersonic compressor cascade[END_REF] were made in 1984. Each of the test has been performed for slightly different flow conditions and a wide variety of data are available from them for validation:

• Inlet and outlet quantities of interest: Mach number, static pressure and incidence, allowing to compute the different turbomachinery quantities such as pressure ratio and deviation (see table 2).

• On blade pressure measurements at different locations on both suction and pressure sides. These data allow to compute and compare the isentropic Mach number at blade wall.

• Wake traverse measurements of total pressure, static pressure and incidence, downstream of the cascade.

Numerical setup

Domain For the numerical CFD study, the ARL-SL19 cascade has been considered to be composed of an infinite number of blades each of infinite length, this allows to consider by translation periodic boundary condition both span-wise and pitch-wise. Therefore, only the domain described in the following figure needs to be meshed for the numerical study. The computation is made 3D but the flow is assumed to be 2D, as in the numerical simulations which will be used for the validation process. During the mesh adaptation process no subdivisions is made in span direction, thus the cells counts will be reported by layer.

Blade Geometry

The blade geometry has been reconstructed using the list of suction and pressure side points coordinates given in [START_REF] Fottner | Test cases for computation of internal flows in aero engine components[END_REF] and closed at the leading and trailing edges by inserting tangent ellipses. Boundary conditions The flow is axially subsonic everywhere in the cascade at normal operating conditions therefore the boundary conditions used at inlet is a subsonic inflow condition, for this condition the tangential (along Y-axis) velocity stagnation pressure and stagnation enthalpy need to be set. For outlet a subsonic outflow is applied, it requires providing a static pressure constraint. Along the blade, viscous wall conditions are used and periodic conditions are set to the remaining boundaries. Finally, the main parameters to set the cascade operating condition are:

• The inlet tangential velocity (derived from the inlet Mach number and incidence) V 1y

• The outlet pressure P 2

Solver

The simulation were performed with RANS modelling using the elsA [START_REF] Cambier | The ONERA elsA CFD software : Input from research and feedback from industry[END_REF] finite volume and cell centered solver (ONERA-SAFRAN property). The turbulence model used is the negative Spalart-Allmaras [START_REF] Allmaras | Modifications and Clarifications for the Implementation of the Spalart-Allmaras Turbulence Model[END_REF] model, this choice has been driven by the good results reported by Piovesan [START_REF] Piovesan | Accurate 2-D modelling of transonic compressor cascade aerodynamics[END_REF] and Zhang [START_REF] Zhang | Investigation of hysteresis loop in thestarting process of supersonic cascade[END_REF] using this model. For the discretisation of the equations the Roe scheme [START_REF] Roe | Approximate Riemann solvers, parameter vectors, and difference schemes[END_REF] is used, since it creates minimal numerical dissipation it is consequently well suited for configurations with shocks.

The mesh adaptation strategy used in this paper works only on unstructured grids, therefore all the calculations are performed using elsA's unstructured capabilities.

Operating conditions

The test points from ONERA's [START_REF] Fourmaux | Resultats d'Essais à S5 CH d'une grille d'Aubes ARL[END_REF] campaign have been used in this study, since the authors had access to extended test data and to the actual blades used for the tests. Therefore, the standard point chosen is at M 1 = 1.66 and P 2 /P 1 = 1.95 (see figure 1). It will be used for the main part of this study. The wind tunnel measurement's uncertainty has been reported along with the data when it was available.

The flow field in the cascade for these operating conditions is represented on the figure 2.

Figure 2: Mach number contours trough the ARL-SL19 cascade, with features highlighted, on a fine mesh (2M elements/layer) An important network of shock waves can be observed in the ARL-SL19 cascade (dashed lines on figure 2). It features a detached shock at the leading edge of the profile spreading upstream for its upper part and reflecting on the other blades for its lower part. A strong shock is also visible on the pressure side interacting with: the leading edge shock reflections, the trailing edge shock and the pressure side boundary layer. Because of this shock's interaction with the boundary layer a recirculation area is formed on the pressure side (dotted lines on figure 2). Finally, a wake is created downstream of the blade and interacts with the shocks (dashed-dotted lines on figure 2).

To compare the results of this study with other CFD studies from the literature, a second test point at M 1 = 1.58 and P 2 /P 1 = 2.12 as been selected. For this test point the results are reported in several papers [START_REF] Piovesan | Accurate 2-D modelling of transonic compressor cascade aerodynamics[END_REF][21] [START_REF] Liu | A new method for rapid shock loss evaluation and reduction for the optimization design of a supersonic compressor cascade[END_REF].

Feature sensors

It has been seen in the previous section that a wide variety of features and interactions are present within the ARL-SL19 flow field. As the goal here is to use these features for the mesh adaptation process, a first step would be to identify sensor that can efficiently isolate one or several features from the flow field. After a bibliographic research and preliminary studies, three types of sensors have been selected: a sensors to detect shock waves, a sensor to detect recirculation area and wakes and a more generalist sensor that can detect all of these features at once. This selection of sensors aims to include most of the ARL-SL19 specific phenomena while providing a representative sample of the different sensors in the literature.

Shock sensor To detect and isolate shock waves and their interactions, the Lovely & Haimes [START_REF] Lovely | Shock detection from computational fluid dynamics results[END_REF] sensor has been used. This sensor is based on the pressure gradient and the Mach number vector (computed from the Mach number and the velocity vector). The idea here is not to detect the zones with strong pressure gradient, as it leads to the detection of unwanted features such as strong expansions. Instead, the goal is to compute the component of the Mach number vector tangent to the pressure gradient, which at a shock, is always normal to the shock. Therefore the computed Mach number is the Mach number normal to the shock which is known to always be greater than one upstream of the shock and lower than one downstream. Consequently any value of the Lovely & Haimes sensor close to one means that a shock wave is detected. In the original version of the sensor, the authors propose to scale the pressure gradient with its norm, however this method has shown to create singular values where the pressure variations are almost null. To avoid this problem the scal-ing have been replaced with the following:

θ lh = -→ ∇p ρ 1 ∆L - → V 2 . - → M (1)
This is the final form of the Lovely & Haimes sensor that will be use in this study. While it does not preserve the full physical meaning of the original sensor it is still very efficient as it detects the zone where the pressure gradient is strong and where the Mach vector is mainly aligned with it.

Wakes and recirculation As RANS models are used in this study, the turbulence in the flow field is modelled by the equations of the chosen model, here the negative Spalart-Allmaras [START_REF] Allmaras | Modifications and Clarifications for the Implementation of the Spalart-Allmaras Turbulence Model[END_REF] model. Given that production of turbulence is especially high in zones with high vorticity and shear stress. It's possible to conceive a sensor, from the turbulence equation, dedicated to the second group of features present in the cascade: wakes and recirculation. To this end, the term of production of the pseudo viscosity equation has been used in this study, it can be written as follows:

θ pk = Cb1ρ ν Ω Ω = - → ω + ν (κ.dw) 2 f v2 (2)
With Cb1 = 0.1355, κ = 0.41, Cv1 = 7.1 and f v2 a damping function.

This sensor can be computed for several other turbulence model such as Menter's [START_REF] Menter | Zonal Two Equation k-w Turbulence Models For Aerodynamic Flows[END_REF] using their specific turbulent production term. Or with the generalist approach featured in Driver and Seegmiler's [START_REF] Driver | Features of a reattaching turbulent shear layer in divergent channel flow[END_REF] work.

General sensor A generalist sensor for the ARL-SL19 should be able to detect at the same time the shocks, the recirculation bubble and the wake present in the cascade. As all of these features show significant velocity variations, a sensor based on velocity gradient should be relevant. Hence, the Mach number gradient has historically been used for this task. In this study, its norm will be considered to construct the following sensor:

θ ∇M = ∇M (3)
It is important to note that this sensor tends to overdetect the strong gradients in boundary layers whereas these zones are known to be correctly refined in the normal to the wall direction for the initial mesh. Moreover, uneven refinement and strong size transition close to the wall can cause convergence issues and harm the quality of the solution. Therefore, it has been chosen to "protect" the boundary layers via a weighting of the sensor value by the smallest cell edge length. This way the thin cells in the boundary region wont be as much represented by the sensor. The authors are aware that this solution is not optimal because with the isotropic refinement used in this study it also prevents the refinements in the tangent to the wall direction whereas they aren't problematic and can highly improve the walls curvature definition.

Adaptation process

The iterative adaptation process used in this study is summarized by the following figure : The mesh adaptation is based on an isotropic insertion of nodes with Cassiopee's [START_REF] Benoit | Cassiopee: A CFD pre-and post-processing tool[END_REF] tools adaptCells. The tools is able to subdivide any given basic cell (parent) into a given number of equivalent new cells (children), see figure 4. As the grids resulting from this subdivision, when conformized, wont be fully hexahedral anymore, adaptCells uses a hierarchical representation of the grids to perform successive mesh adaptations. The tool performs smoothing when necessary to limit the strong size transition ratio induced by the adaptation. It is also able to perform directional refinement, which as been used here to prevent the insertion of new cells span wise. Finally, adaptCells is able to enforce coincidence at periodic interface by splitting the required opposites faces.

To transform the sensors values into integer subdivisions count for adaptCells, a quantile criteria has been used. It work by taking a given fraction (or quantile) of the For example choosing the 75 th percentile (or 3 rd quartile) will flag the 25% of cells having the highest values of sensor. As the parameter of the criteria can have a huge impact on the adaptation process, it will be studied in detail in section 4.

After the subdivision the resulting mesh is re-projected onto the high resolution geometry of the ARL-SL19 using Ael's Quantum mesh deformation tool. This step is mandatory to ensure that the resolution of wall curvatures is improved when the available number of points along them increases from the adaptation.

The final step of the adaptation is a prolongation of the solution of the parent cells to their children in order to create a good initialisation state for the next step.

Furthermore, at each step a convergence criterion is computed to determine if the adaptation loop should continue or not.

For this study, as the goal is to better understand how the adaptation process behaves no convergence criterion have been used, instead, the adaptation are all run until they either reach a prescribed number of adaptation iterations or a maximum number of elements per layer.

NUMERICAL UNCERTAINTY AND EXPERIMENTAL VALIDATION

The first step of this study has been to validate the numerical method used for the simulation of the ARL-SL19 cascade and to evaluate the solution sensitivity to the mesh.

Validation

For the validation part of the study, the results using a fine grid have been compared at the reference M 1 = 1.58, P 2 /P 1 = 2.12 test point with the data from Tweedt's testing [START_REF] Tweedt | Experimental investigation of the performance of a supersonic compressor cascade[END_REF] and with the numerical study of Piovesan [START_REF] Piovesan | Accurate 2-D modelling of transonic compressor cascade aerodynamics[END_REF] on the ARL-SL19 cascade. The first validation aims at showing that the numerical set-up can efficiently represents the flow field at the desired test point. While the second is done to ensure that the CFD results are coherent with the ones of Piovesan which have also been made with a Spalart-Allmaras model and a Roe scheme but with a purely 2D configuration. 12. From this study and from the CFD study of Piovesan [START_REF] Piovesan | Accurate 2-D modelling of transonic compressor cascade aerodynamics[END_REF] and the experimentation of Tweedt [START_REF] Tweedt | Experimental investigation of the performance of a supersonic compressor cascade[END_REF] On figure 5 a general good agreement between the results of this study and Piovesan's is found. Both feature the same small separation at 80% of chord on the suction side induced by the leading edge shock reflection and the separation on the pressure side due to the strong shock is also consistent.

Nevertheless, both CFD studies struggle to accurately represent the plateau of isentropic Mach number before the strong pressure side shock measured during Tweedt's experiment.

According to the results shown is this section, the numerical method seems to be accurate enough to give a good representation of the ARL-SL19 flow field at the desired test points. For the following part it has been decided, using the mesh sensitivity analysis, to start all the adaptation computations from the first grid and solution reported on figure 6 (5.10 3 elements per layer).

Numerical error and uncertainty

The procedure used in this paper for numerical error and uncertainty approximation is from Ec ¸a [START_REF] Ec | A procedure for the estimation of the numerical uncertainty of CFD calculations based on grid refinement studies[END_REF]. It relies on the following error formulation :

ε φ ≈ φ i -φ 0 = αh p i (4)
A fitting is performed on at least 4 grids to find the unknown parameters : φ 0 the approximated exact solution, p the observed order of grid convergence and α a constant. From these values, the numerical uncertainty can be computed (see equations ( 20) and ( 21) in Ec ¸a [START_REF] Ec | A procedure for the estimation of the numerical uncertainty of CFD calculations based on grid refinement studies[END_REF]). This method is especially interesting because it allows to perform the analysis with more than the minimum of four grids, therefore the estimation's quality can be improved from additional grids. Ec ¸a also gives a clear methodology to apply when the order of convergence of the method underlines poor convergence. To use this method two assumptions are made. First, the different grids must be in the "asymptotic range" then the grids must be geometrically similar in order to be representable by a single size parameter.

Two methods have been used and compared for this study, first adaptCells has been used to create a sequence of self similar fully refined grid. Then the authors manually created grids that match the cell count of the different grids of the previous sequence without isotropically refining the mesh. The main difference are that first cell height close to wall was not changed and that the refinement in the upstream and far downstream regions of the computational domain were limited.

For both of these methods the compression ratio has been computed at each step, along with its numerical error and uncertainty. Five steps have been made starting with a 5.10 3 elements by layer mesh and finishing with close to 2.10 6 elements by layer meshes as can be seen on figure 6. On this figure a good asymptotic convergence toward a compression ratio of 1.957 can be observed when the mesh size increases, with the self similar grids. The outlet pressure P 2 derives from the ratio P 2 /P 1 and the inlet pressure P 1 computed from the stagnation pressure and the inlet Mach number M 1 . In CFD, M 1 is not specified by user but results from the computation. Thus, any discrepancy in M 1 prediction results in a deviation on P 2 and leads to the small difference encountered. With the manual grids the static pressure ratio at the end is close to the one obtained before but the convergence is degraded by some overshooting made with the second and third grid. Therefore, the numerical error is higher resulting in noticeably more important numerical uncertainty even with the finest grids. This could be induced by lower mesh quality at important locations of the domain, or by the density of the mesh not being sufficient in those same key areas.

For computational time issues and considering the error level in static pressure ratio of the fourth isotropically refined grid, its cell count (4.10 5 elements per layer) is

SENSOR BASED MESH ADAPTA-TION INVESTIGATION

The different sensors and parameters used in the adaptive mesh refinement procedure can have a huge impact on the flow solution. Therefore, to evaluate and improve the robustness of the method several parametric studies were performed. First, the impact of the quantile criterion on the adaptation using the Lovely& Haimes shock sensor was examined. Then, the performance of the different sensors were compared with a fixed quantile criterion. Finally, the efficiency of sensor combination was assessed by combining the Lovely& Haimes shock sensor with the turbulent production sensor in different proportions.

Influence of quantile on the adaptation

The choice of an adequate quantile criterion for adaptation is not an easy task, choosing a to low parameter could make the adaptation miss certain features. Whereas choosing a to high quantile parameter could make the adaptation inefficient. On the figure 7 the static pressure ratio evolution during the adaptation process has been plotted for eight different quantile criteria.

A strong link between the quantity of refinement done at each step and the intensity of the induced variation on the static pressure ratio is found. Except for the 95 th and 70 th criterion all the studied criteria induce a strong Figure 7: Error relative to reference on the static pressure ratio trough the ARL-SL19 cascade for several adaptation iterations. Using the Lovely & Haimes [START_REF] Lovely | Shock detection from computational fluid dynamics results[END_REF] shock sensor and a quantile type refinement criterion ranging from the 95 th to the 60 th augmentation of the static pressure ratio during the first adaptation iterations. During these steps the criterion converges significantly faster toward a final plateau than the reference isotropic refinement. During the adaptation process elements are added in regions that requires more refinement than the simple isotropic method can provide. The adaptation made with the 70 th quantile criterion does not follow the general trend observed earlier, indeed, it converges toward a final value between the 85 th and 80 th criteria's final value. Moreover, its initial convergence rate is lower than the others and it does not top out either. The figure 8 gives some more information to explain this difference. On the figure, one can see that during the adaptation the zones refined with each criterion are interlocked, the biggest quantile contain the smallest ones. Also, the zone refined during the later iterations of an adaptation process are generally contained inside zones that have already been refined. This behaviour is due to the nature of shocks, because they are infinitely thin discontinuities, the smaller the mesh gets through the shock, the stronger the gradients get. Therefore, if the shock sensor is not limited the refined zone is getting more and more focus on the exact locations of the shocks. Yet, the 70 th criterion does not follow this pattern. At both iterations 3 and 6 it can be observed that some of the zones refined with it (in yellow) are different than with the other criteria. The discrepancies in refinement area seems to appear on the suction side near the shock boundary layer interaction at 80% of the chord. Before this point the reflections of the leading edge shock are similar for all the criteria, but after, the resulting shock have a slight offset with the 70 th criterion. This difference directly impact the strong shock and all Figure 8: Refined zone evolution on the ARL-SL19 cascade for the quantile criterion ranging from the 95 th to the 60 th and for different adaptation steps. The refined zone are coloured by quantile value and nested from smaller quantile to bigger ones. the downstream shock structures. An observation of the refinement made in this zone shown on figure 9 underlines that with the criteria lower than the 70 th the boundary layer is not refined while with the ones greater than the 70 th the boundary layer is refined down to the wall uniformly. In contrast, with the 70 th criterion the refinement seems to stop inside the boundary layer, in the ten first layers of cells. The induced size transition between un-refined and refined regions in this area could explain the difference observed. Knowing that the shock and the boundary layer are highly coupled and that the mesh refinement near the wall has a strong impact on the boundary layer, it can induce a modification of the shock position. As a result, the protection of already well resolved boundary layers seems to be essential to the robustness of the adaptation process when using isotropic refinement as is done here.

Finally, there is a close proximity between the results with a 65 th and 60 th quantile criterion, thus increasing the quantity of refined cells more than that will be inefficient as it wont change the final solution significantly. It would be interesting for a future study to understand how the adaptation process can lead to value greater than the supposedly converged reference.

Sensors performance comparison

To better understand how the specificities of each sensor influence their behaviour and performance, the three sensors have been compared using, for each, the best performing quantile criterion observed in the study. The resulting static pressure ratio evolution have been reported on figure 10. The static pressure ratio evolution during the adaptation, observed on figure 10 is very different from one sensor to an other. As seen in the previous section, the shock sensor allows a faster convergence of the static pressure ratio toward a final plateau than the isotropic solution. Meanwhile, using the Mach number gradient sensor does not result in such a strong asymptotic convergence on the static pressure ratio. With this sensor the quantity of interest first dips lower than its initial value for two iteration before slowly increasing toward the reference value. The static pressure ratio peaks three iterations before the end of the process and decreases for the last iterations. This behaviour on the first and last iterations is due to the shock system moving toward the leading edge of the blade and quitting the already refined zone. As for the first iterations, the decline observed at the end is likely to be followed by a new increase sequence when the adaptation "catches up" with the new position of the shocks. This phenomenons is due to the fact that the Mach number sensor is a global sensor, therefore is refines simultaneously several features. On the ARL-SL19 most of the feature are linked meaning that refining one of them will affect the others. Here, refining the recirculation area and wake will affect the position of the shocks, which themselves will affect the recirculation and wake. As a consequence, the adaptation process will oscillate around its final state before converging. In this case the approach chosen here could be improved. As no de-refinement is performed, the zones refined where the features are not located any-more remain, increasing the size of the successive meshes.

The turbulent production sensor does not seem to have any positive effect, when applied alone, on the static pressure ratio. Because, with this sensor, a decrease is observed at first followed by a re-increase and a plateau near the initial value. That is expected because the physics of the ARL-SL19 is dominated by the shocks, as a result not refining them limits the quality improvement that can be made. Nevertheless, the results obtained with the shock sensor show that refining only the shock isn't sufficient either to reach the reference solution.

Sensors combination

In the previous section it was concluded that neither the shock sensor nor the turbulence sensor alone is sufficient to reach the desired solution, therefore the combination of these two sensor was investigated in this part. Here, the zone flagged for refinement by each sensor is combined (by union) at each step. The goal of this part of the study is to find the best quantile criterion to apply to each sensor before the combination to get the most accurate results. To do so a parametric study has been made by combining both sensor with quantile criteria ranging from the 95 th to the 60 th (only the combinations creating a criterion greater or equal to the 60 th have been studied). The error relative to the reference value at the end of each adaptation have been reported on figure 11. It can be observed that the difference with the reference for the initial grid is close Figure 11: Relative difference between the static pressure ratio at the end of the simulation of combined Lovely & Haimes and turbulent production sensors and the reference simulation (on figure 6). For combinations of quantile parameter ranging between 95 th and 60 th for both sensors to 1.6%, this means that every adaptation resulting in a greater error value has ultimately worsen the simulation of the flow field, either because it converged toward a worse value or because it has not done enough iteration to fully converge. In contrast the adaptation featuring an error value lower to 1.6% have improved the representation of the flow filed. However, it doesn't mean that these adaptations were more efficient than the isotropic refinement. On the map the worse adaptations are mainly found for small quantile of the shock sensor but for all turbulence sensor's quantile value. This illustrate the strong important of shock resolution for the ARL-SL19, if none or to few refinement is done along the shocks the solution can't converge toward the value. Yet, the refinement of the recirculation and wake area isn't useless. The results show that the adaptations with only the shock sensor (horizontal axis) are not the ones with the lowest error. Indeed, coupling a high criterion for the shocks with a moderate criterion for the turbulence seems to give the best results. The adaptation allowing the reach the lowest error of the experience set, is the one made with a 75 th criterion for the shocks and a 85 th criterion for the turbulence. Generally, it seems that choosing a combination that allows for high number of total flagged cells with at least an 80 th shock criterion is the most efficient way of combining these two sensors. Strong refinements have some advantages for the robustness of the adaptation process. Because the mesh is usually refined further away from the features their movement is less likely to cause convergence rate issue. Besides, refining larger zones around the features has the side effect to allow the refinement of other nearby features that wouldn't have been refined otherwise. Moreover, refining faster the mesh ultimately means that less iteration are to be performed to reach a cell count target or a convergence criterion, meaning that the adaptation is likely to be less expensive. However there is also some drawbacks to this approach, as seen earlier refining far from the features can induce uneven refinement in the boundary layers and strong flow destabilisations.

Finally, one should keep in mind that the error measure provided in the study is not perfect, because it does not give any insight on convergence, a low value of error only means that the adaptation ended close to the reference value, but not that it converged toward this value. Also, because of the union made for the combination of criteria the map on figure 11 is not accurate: combining a 75 th criterion and a 85 th for example result in a refined zone smaller or equal to the zone refined by a 60 th criterion.

CONCLUSIONS AND PERSPEC-TIVES

A fully automated mesh adaptation process has been presented and applied to the ARL-SL19 supersonic cascade.

Once the numerical settings validated, a thorough convergence and error analysis have been carried out, providing insights on the most efficient and robust mesh adaptation strategy. First of all, a smooth geometric evolution of the cell size normal to the wall in the boundary layer seems to be play an important role. If this condition is not ensured in the key areas of the flow field such as shock/boundary-layer interaction and flow separations,the analysis of the quantities of interest have undergone significant alterations. Particular attention must be paid regarding the adaptation mechanism within the boundary layer. Then, when dealing with interacting features such as shock/boundary-layer interaction the efficiency of the adaptation process can be altered, reducing the order of convergence of the method and increasing the number of iterations needed. Being able to de-refine the mesh would be an interesting solution to keep the adaptation cost acceptable. Finally, combining sensors seemed to be a good practice to catch simultaneously several features. It has been observed that there can be an optimal combination of the sensors which required a lot of computational effort to be found. Further investigations would be necessary for a better understanding of this combination sensitivity with regard to the operating conditions. The methodologies and knowledge developed during this study is essential for future developments of adaptive mesh refinement, with a possible extension to the domain of parametric uncertainty. To do so the adaptation should prove to be robust enough to take into account the uncertainties when applied to different operating flow conditions.

Figure 1 :

 1 Figure 1: Schematic representation of the ARL-SL19 CFD configuration

Figure 3 :

 3 Figure 3: Main steps of the adaptation process

Figure 4 :

 4 Figure 4: Main steps of the adaptation process

Figure 5 :

 5 Figure 5: Isentropic Mach Number along the ARL-SL19 blade for the test point at M 1 = 1.58, P 2 /P 1 = 2.12. From this study and from the CFD study of Piovesan[START_REF] Piovesan | Accurate 2-D modelling of transonic compressor cascade aerodynamics[END_REF] and the experimentation of Tweedt[START_REF] Tweedt | Experimental investigation of the performance of a supersonic compressor cascade[END_REF]

Figure 6 :

 6 Figure 6: Compression ratio and numerical uncertainty versus grids size for successive self similar grids and their equivalent in size manually made grids.The convergence order of the self similar grids is 1.67

Figure 9 :

 9 Figure 9: Refinement differences in the inner boundary layer, on the suction side at 80% of chord

Figure 10 :

 10 Figure 10: Error relative to reference on the static pressure ratio trough the ARL-SL19 cascade for several adaptation iterations. Using the Lovely & haimes [15] shock sensor, the Mach number gradient sensor and the turbulent production sensor. Respectively with quantile 75 th , 85 th and 75 th refinement criteria.

Table 2 :

 2 Experimental data from Fourmaux's testing[START_REF] Fottner | Test cases for computation of internal flows in aero engine components[END_REF][12] * Axial Velocity Density Ratio * * Total pressure loss

			The geometric configuration and ex-
	perimental quantities associated with this test point are
	reported in the following tables:
		Parameter	Value
		c	85.0 mm
		t	55.6 mm
		β s	146.93 •
	Table 1: Experimental parameters from Fourmaux's test-
	ing [10][12]
	Inlet		Outlet
	M 1	1.66 ± 0.02 M 2	1.098 ± 0.02
	β 1	57.5 • ± 0.2 β 2	60.5 • ± 0.5
	AVDR * 1.057	ω * *	0.1363
			P 2 /P 1 1.95

ACKNOWLEDGEMENT

This work was partially supported by the SONICE project, granted by the French Directorate General for Civil Aviation (DGAC).