domain: spi phys

text:
Gérald Carrier

Guillaume Arnoult

Nicolo Fabbiane

Jean-Sébastien Schotté

Christophe David

Sébastien Defoort

Martin Delavenne

Emmanuel Bénard

Multidisciplinary analysis and design of strut-braced wing concept for medium range aircraft

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

II.Introduction

Reducing the environmental footprint of aviation and especially its impact on climate change has now become the major challenge of the aeronautical industry. In this context, the European Commission has set ambitious goals for aviation defined in the "Fligthpath 2050"[1], including 75% reduction of CO2 and 90% reduction of NOx emissions, as well as a 65% reduction of noise. Achieving such targets will likely be possible only through radical evolutions of the aircraft architecture, propulsion system and the introduction of disruptive technologies. For what concerns the aircraft architecture, the Blended Wing Body (BWB) and the Strut-Braced Wing (SBW) seem today to be the two most promising alternatives to the conventional "Twin Underwing Turbofan Tube-and-Wing" configuration that has become the predominant solution over the last decades. The BWB, and the Hybrid Wing Body, are currently investigated in several important research efforts in Europe, USA and China [START_REF] Collier | Status of Research and Development of the Hybrid Wing Body (HWB) Aircraft Concept[END_REF][3] [START_REF] Zhenli | Binqian ZHANG Assessment on Critical Technologies for Conceptual Design of Blended Wing Body Civil Aircraft[END_REF]. The present paper addresses investigations of the SBW configuration. The main rationale for the SBW configuration is to enable significant increase of the wing aspect ratio (AR) to improve overall aircraft performances through aerodynamic efficiency gains without the tremendous wing weight increase that a cantilever wing of similar AR would suffer from. This is achieved through a reduction of the bending moment in the inner wingbox permitted by a traction of the strut at its junction with the wing, as illustrated in Figure 1. This enables to downsize the wing box structural components in the inner wing that are primarily sized by the bending moment in the case of a cantilever wing. This concept can be used to increase the wing aspect ratio and take advantage of an aerodynamic efficiency improvement obtained by reducing the lift induced drag component, while alleviating the increase of the wing structural weight associated with higher span (involving higher bending moment in the central wing-box). Therefore adding a strut, or more generally a truss structure, opens an additional degree of freedom for the designer to better exploit high aspect ratio wings at overall aircraft level by achieving new and better trade-off between aerodynamic efficiency and wing weight, compared to conventional cantilevered-wing concept.

Such trade-offs for the SBW were already investigated at ONERA in the ALBATROS project [START_REF] Carrier | Investigation of a strut-braced wing configuration for future commercial transport[END_REF]. More recently, in the Clean-Sky 2 thematic topics project U-HARWARD [6] different aircraft-level solutions enabling wing aspect ratio increase (such as cantilever flexible wings, SBW and folding wingtips) are currently investigated to reduce the environmental impact and improve the overall aircraft performance. In this project, ONERA and ISAE-SUPAERO are exploring the potential benefits and eventual pitfalls of the HAR-SBW concept through a multi-level multidisciplinary design approach including overall aircraft conceptual design (OACD) and refined disciplinary design studies of the aerodynamic and structural key aspects of HAR-SBW concept by high fidelity analysis methods.

This paper describes the investigation of the HAR-SBW concept and tries to evaluate its potential benefits at overall aircraft level. Basically, our objective is to explore the impacts at aircraft level of the increase of the wing aspect ratio, with and without strut and identify what are the maximum benefits expectable with and without strut and for which values of AR (and wing design parameters) this occurs.

The conceptual design approach used to assess the overall aircraft impacts of HARW and SBW concepts is introduced in the first section. Then the high fidelity aerodynamic analysis and design approach developed and used for this project to evaluate the HAR-SBW aerodynamic benefits is described and illustrated with preliminary design results in the second section. The third section describes the development of a wing structural weight evaluation module for the SBW concept and illustrates its capabilities. Finally, the last section describes the first results obtained with this multilevel OAD process combining the different models describes in sections 1 to 3 for the HAR-SBW concept. The benefits at overall aircraft level of the SBW concept is evaluated by investigating different values of the wing aspect ratio ranging from 12 to 22 and comparing the resulting aircraft characteristics to those of a conventional tube-and-wing aircraft configuration designed for the same mission with the same tools and approach.

III.Conceptual design approach

A. Top Level Aircraft Requirements (TLAR)

The aircraft investigated are designed for the following mission characteristics, close to the Airbus A321-LR aircraft TLAR:

 a sizing range of 7400 km,  a sizing payload of 18 tons (~200 pax),  a maximum payload of 23 tons  a cruise Mach number of 0.78;

B. FAST-OAD tools for conceptual design

The ONERA-ISAE FAST-OAD aircraft conceptual design framework [START_REF] David | From FAST to FAST-OAD: An open source framework for rapid Overall Aircraft Design[END_REF] is used in this work to perform the overall aircraft design in order to maximize the SBW benefits in terms of overall aircraft performances. FAST-OAD is implementing a multidisciplinary sizing process of the aircraft (see Figure 2) and is designed so that models can be added or removed from the process, which allows to adapt the OAD problem to the studied configuration and introduce various fidelity models for the same discipline. First this framework has been applied to generate a conventional aircraft (tube and wing) that meets the chosen TLARs, with the intent of providing a reference configuration to be compared to the SBW designs that will be lately investigated.

The models initially in place in this process (before introducing specific models for the SBW configuration) lie in the family of the so-called "Level 0" (statistical/analytical) models, well adapted for rapid sizing of conventional configurations:

 Geometry: semi-empirical formulas for adapting general characteristics of aircraft to TLARs  Weight: semi-empirical formulas for estimating weights (detailed breakdown by components)  Handling Qualities: analytical formulas for sizing empennage and setting wing position along fuselage  Aerodynamics: analytical and empirical formulas (Oswald coefficient, flat plate friction drag)  Performances: estimation of fuel consumption using time-step integration of point-mass mission simulations  Propulsion: semi-empirical parametric model for turbofan engine (directly used by performance model).

The design process of FAST-OAD with this collection of models is validated against an A320-CEO-like configuration [8] (Table 1) and can easily be extended to the parameters of an A321-like aircraft. The next steps is to introduce in this OAD process the capability to perform a sizing loop on a SBW configuration. Low-fidelity aerodynamic and structural weight evaluation models, specific to SBW, are integrated into the FAST-OAD framework.

Weight/Structure

The wings and struts primary structure weights are computed using physical analytical models based on beam theory [8]. The wing structure is simplified to an equivalent spar plus skin model with the spar flanges supporting bending moment, the web supporting the shear and the skin supporting the torsion This methodology has been extended to the SBW considering they only support traction and compression forces, as it will be described in section V. This method requires aerodynamic forces and moment spanwise distributions that can be estimated from elliptical loading or retrieved from Vortex Lattice Method computations as explained below.

Also implemented in FAST-OAD, a beam-like finite element model of the struts and the wings can be substituted to the previous models to compute both primary structural weights and deformations. Sections properties (area and inertia) are computed considering web thickness, equivalent upper and lower skin thickness and stringers area. The weight simply derives from material properties and thickness necessary to sustain static manoeuvre loads (2.5g, -1g). This approach can be coupled to aerodynamic force computation to size the wings and the struts within an aerostructural process. This approach is relevant for high-aspect ratio wings with non-negligible deformations.

The weight of the other components (fuselage, tails, systems, …) is assessed through statistical formulations.

Aerodynamics

Originally, aerodynamic models implemented in FAST-OAD rely on analytical and empirical formulations. The friction drag comes from a flat plate approximation corrected by compressible, sweep, camber and thickness effects. The induced and wave drag are derived from empirical formulations of the Oswald coefficient and critical Mach number estimates. These formulations may be extended to the struts that are mainly accounted through their contribution to viscous drag and their interference with the wing.

Vortex Lattice Method computations can also be performed within FAST-OAD using AVL software [START_REF] Drela | [END_REF]. Particularly, loads distributions can be extracted to feed wings and struts weight computations and the aerostructural process.

Propulsion

This module, not specific to SBW configuration, is however tweaked with respect to the original turbofan model of the FAST-OAD distribution to mimic the performances of an EIS 2015 engine such as the LEAP-1A, because a future SBW aircraft will likely use that family of turbofans and has to be compared to a -NEO aircraft. Around 10% SFC savings compared to CEO can be expected from this new propulsion model.

IV.High-fidelity Aerodynamic Design and Analysis

A. Description of the parametric automated aerodynamic analysis and design framework

For the sake of performing the aerodynamic design of the SBW configuration and feed the aerodynamic model used in the OAD environment with high-fidelity aerodynamic data , the following Hi-Fi design framework has been set up (Figure 3). First, a fully parametric geometry model has been generated using the Engineering Sketch Pad (ESP) [START_REF] Haimes | The Engineering Sketch Pad: A Solid-Modeling, Feature-Based, Web-Enabled System for Building Parametric Geometry[END_REF] tool developed by MIT. This freely available and open source software allows the user to build its own geometry with different levels of fidelity. This geometry can be connected to Analysis Interface Modules (AIM) through the Computational Aircraft Prototype Synthesis (CAPS) [START_REF] Christopher | Parameterized, Multi-fidelity Aircraft Geometry and Analysis for MDAO Studies using CAPS[END_REF] software (for meshing, structural analysis,..). In the current framework, four AIM modules are used and assembled into an automated design analysis chain using plain Python scripts (Figure 3): Pointwise (meshing), SU2 (CFD calculations), Friction module.

The commercial mesh generation software Pointwise enables the engineer to create structured and unstructured meshes, covering a wide range of CFD calculation needs. It is used there in batch mode using the Glyph scripts capability that enables to automate the generation of unstructured CFD mesh around almost arbitrary aircraft geometries [START_REF] Karman | Automatic Unstructured Mesh Generation with Geometry Attribution[END_REF]. Then, the open-source CFD software SU 2 [START_REF] Palacios | Stanford University Unstructured (SU2): An Open-source Integrated Computational, Environment for Multi-Physics Simulation and Design[END_REF] written in C++ for the numerical solution of partial differential equations, provides the flow field around a geometry and thus, the aerodynamic coefficients (solving either the RANS or Euler equations). Finally, the ONERA Far-Field Drag post-processor FFD00 [START_REF] Destarac | Drag/Thrust Analysis of Jet-Propelled Transonic Transport Aircraft: Definition of Physical DragComponents[END_REF] is used to evaluate lift-induced, wave and viscous drag components and suppress the spurious irreversible drag inherent to numerical CFD solutions.

B. Aerodynamic analysis and design process set-up and demonstration on A321 aircraft configuration

As introduced in section III.A, the reference configuration considered in this work is a long-range (LR) version derived from a single-aisle (SA) short-to-medium range (SMR) aircraft similar to a A321-LR. Each geometrical parameter is written in the ESP/EGADS and describe in Figure 4 to Figure 6.

-A321-LR reference aircraft geometrical characteristics used for the reconstruction

An illustration of the reconstructed A321-LR is presented in Figure 7 and used in the process presented in the previous section to achieve CFD simulations and analysis. These computations are performed using SU² solver with target lift coefficients which allows drawing the polar curve of the aircraft. A far field analysis is also performed during the process to identify the drag components, as shown in Figure 9 where the influence of the Mach number on the viscous (triangle shape), wave (diamond shape) and induced (circle shape) components of the drag is illustrated.

C. Aerodynamic analysis and design of the SBW aircraft configuration

SBW aircraft geometry parameterization

The objective of this part of the study is to define a generic parameterization of the strut, which will be used in the future for optimization purposes. The geometry of the strut is thus defined by 9 top-level parameters, which are gathered in Table 2.

Parameter Units Description

Span_frac

Table 2 -SBW geometry parameters describing the strut

These last parameters are directly taken from an aerodynamics consideration since they are supposed to have an influence on the channel effect described in [START_REF] Carrier | Investigation of a strut-braced wing configuration for future commercial transport[END_REF] and induced by the junction between the wing and the strut, leading to an important contribution of the wave drag to the global drag of the aircraft.

Aerodynamic analysis of the reference Strut-Braced Wing configuration

The reference configuration for this part of the study is taken from a former project studied at ONERA, called ALBATROS and described in [START_REF] Carrier | Investigation of a strut-braced wing configuration for future commercial transport[END_REF]. This specific configuration, illustrated in Figure 11, is reconstructed using the ESP/EGADS tools in the same way as the A321-LR reference configuration for the cantilever wing. After these first computations, it appears that the conventional tail configuration is not well adapted for this kind of configurations. The engines could be located either under the wing (which induces high interaction with the strut) or at the rear of the fuselage. In both cases, the engine exhaust would impact directly the HTP, with the risk of reducing its efficiency. A T-Tail configuration is thus designed and used in the next steps of the study, as illustrated in Figure 13.

Figure 13 -Modelling of the T-tail geometry

Construction of a strut analytical modelling

Based on this last geometrical modelling, the influence of the strut on the drag components is studied. To perform this analysis, a strut-less aircraft configuration and a rear strut-braced configuration are considered. Both configurations are illustrated in Figure 14. The aim here is to isolate the contribution of the strut to define an analytical model of the additional drag generated by the strut that could be integrated in the FAST-OAD process. RANS meshes are automatically generated (Figure 14) with Pointwise and SU² computations are achieved for the target lift coefficients in [0.3; 0.4; 0.5; 0.6; 0.7; 0.8].

Figure 14 -Illustration of strut-less and rear strut-braced configurations

The rear strut-braced configuration is considered here as it allows reducing the interaction between the wing and the strut, thus lowering the impact of the channel effect and the wave drag generated by the configuration.

In order to generalize this modelling to any strut geometry that will be considered in the future, only the additional friction drag is considered here and multiplied by the ratio

𝑆𝑤𝑒𝑡 𝑤𝑖𝑛𝑔 𝑆𝑤𝑒𝑡 𝑠𝑡𝑟𝑢𝑡
. The evolution of this coefficient is shown in Figure 15 as a function of the lift coefficient at which the computations were achieved.

Figure 15 -Strut viscous correction factor

This is then directly usable in the FAST-OAD framework to rapidly evaluate the viscous drag induced by a strut and identify the most promising design of strut-braced wing configurations from the aero-structural point of view resulting in optimum overall aircraft performances.

The sweep values considered here are not optimum but the authors hope to be able to identify an optimum configuration in the future and examine its impact on the other disciplines involved in the design process of the strutbraced wing configuration.

Illustration of aerodynamic design optimization capability for the SBW configuration

The CFD-based aerodynamic analysis process described previously (Figure 3) can be used to perform automated aerodynamic design optimization by linking this parametric analysis process with a numerical optimization algorithm. The capability of conducting aerodynamic design based on shape optimization using the present set of tools is illustrated hereafter. For that a simplified wing design exercise on the SBW configuration is performed using gradient based optimization and Euler flow modelling.

The optimization problem is defined as a full SBW aircraft inviscid drag minimization in transonic cruise condition (Mach number 0.78) at constant lift coefficient (CL=0.5). A subset of the complete set of geometry parameters defined in the ESP model are used to control the airfoil geometry of the wing. Fourteen design parameters are used to control simultaneously the wing twist and the wing airfoil geometry using the implementation of the CST airfoil parameterization by Kulfan [START_REF] Kulfan | A Universal Parametric Geometry Representation Method -CST[END_REF]. The optimization was performed using the Modified Method of Feasible Direction algorithm of Vanderplaats [START_REF] Gn | A robust Feasible Directions algorithm for design synthesis[END_REF] and included a constraint on the minimum wing airfoil relative thickness of 10%, while the initial wing geometry is equipped with aifoil of 12% relative thickness.

The evolution of the 14 design variables and optimization functions are plotted in Figure 16. A comparison of the initial and optimized airfoil geometry with the corresponding pressure distributions on an outboard wing section at 75% of span is given in Figure 17. The optimization achieved a reduction of the inviscid drag of about 100 drag counts, thanks mostly to the reduction of the wave drag which was significant on the initial configuration due to the "channel flow" between the wing intrados and the strut in these transonic flow conditions as it is illustrated in Figure 18. The wing extrados pressure recovery is also significantly improved as illustrated in Figure 17 and Figure 18.

This demonstrates that the CFD-based parameterized aerodynamic model of the complete configuration can be successfully used to perform basic aerodynamic design tasks by automated optimization. This capability is particularly important for evaluating the aerodynamic efficiency in transonic flow condition of an aircraft concept such as the SBW and deriving aerodynamic data that can be used in an OAD process like in III.B

Illustration of engine integration design capability

The CFD-based process with parameterized geometry based on ESP [START_REF] Haimes | The Engineering Sketch Pad: A Solid-Modeling, Feature-Based, Web-Enabled System for Building Parametric Geometry[END_REF] was also used successfully to investigate engine integration aspect on the SBW configuration. A geometry model of a generic high by-pass ratio engine nacelle was developed with ESP and added to the glider configuration, hold by a basic (parallel flange) engine

V. Structural Design and Analysis of the SBW configuration A. Construction of a wing weight estimation model for the SBW configuration

As exposed above, the initial version of FAST-OAD relies on semi-empirical formulas for weight estimates and aerodynamics. In order to account for strut-braced wing configurations, new loads, weight estimates and aerodynamic models have been integrated within FAST-OAD framework.

Wing loads model

First, a loads evaluation module has been set up to compute both limit and ultimate load factor (𝑛 𝑧 ) considering a safety factor imposed by the user, typically 1.5. So far, only one load case is implemented that consists in a pull-up manoeuvre at MTOW. Then, specific sub-models are developed to compute the external and corresponding internal loads resulting from aerodynamic, fuel distribution, engines, structural weight itself and induced by the strut. Only the aerodynamic loads computation is mandatory, the consideration of other contributions is let to the user.

a. Geometric considerations

Prior to any loads computations, some geometric considerations have to be introduced. As illustrated in Figure 20, parameters are introduced to take into account the geometry of the wing box. Particularly, the local wingbox chord is deduced from the local aerodynamic chord considering the spars chord ratios and a projection factor, 𝑓 𝜙 𝑒 , to account for the ratio between the chord in aerodynamic axis and its projection onto elastic axis:

𝑐 𝑏𝑜𝑥(𝑦) = (𝑘 𝑟𝑒𝑎𝑟-𝑠𝑝𝑎𝑟 -𝑘 𝑓𝑟𝑜𝑛𝑡-𝑠𝑝𝑎𝑟 ) • 𝑐 𝑎𝑒𝑟𝑜 (𝑦) • 𝑓 𝜙 𝑒 • cos 𝜙 𝑒
The location of the elastic axis, is defined from the contribution of each spar to the inertia of the box. This contribution, x, is fixed by the user and the chord-wise location of the elastic axis is deduced from rear and front spar ratios as follows:

𝑘 𝑒𝑎 = 𝑥 • 𝑘 𝑓𝑟𝑜𝑛𝑡-𝑠𝑝𝑎𝑟 + (1 -𝑥) • 𝑘 𝑟𝑒𝑎𝑟-𝑠𝑝𝑎𝑟
The height of the wingbox is also approximated from airfoil characteristics considering the relative thickness and a coefficient that takes into account the non-rectangular shape of the box, 𝑘 𝑠𝑝𝑎𝑟 , whose typical value is around 0.93:

ℎ 𝑠𝑝𝑎𝑟 = 𝑘 𝑠𝑝𝑎𝑟 • 𝑐 𝑎𝑒𝑟𝑜 • 𝑒 𝑎𝑖𝑟𝑓𝑜𝑖𝑙

Figure 20 -Wingbox parametrisation for loads computation

This parametrization is valid whatever the wing type considered (cantilever or strut-braced). But, the consideration of the SBW configuration requires additional geometric parameters that have been introduced in section C.1 and are summarized for what concerns the present structural model in Figure 21. Particularly, the strut/wing junction relative location, 𝜂 𝑠𝑡𝑟𝑢𝑡 and the distance between wing and strut roots, ℎ 𝑠𝑡𝑟𝑢𝑡 , note also that the chord of the strut is assumed to be constant span-wise as well as the size of the box. Figure 21

b. Aerodynamic loads

In the current version of FAST-OAD, no assumption about lift distribution is made. Therefore, a new model capable of estimating the local lift has been developed, either from a simple elliptical loading for rapid evaluations or from a computation with a Vortex Lattice Method taking into account airfoil camber, wing planform and twist.

The aerodynamic shear is then computed at any load section considering the local lift distribution 𝑙(𝑦):

𝐹 𝑧,𝑎𝑒𝑟𝑜 (𝑦 𝑙𝑜𝑎𝑑 ) = ∫ 𝑙(𝑦)𝑑𝑦 𝑏/2 𝑦 𝑙𝑜𝑎𝑑
Then, the bending moment is deduced, taking into account the local sweep angle of the elastic axis with respect to aerodynamic reference axis 𝜙 𝑒 :

𝑀 𝑥,𝑎𝑒𝑟𝑜 (𝑦 𝑙𝑜𝑎𝑑 ) = ∫ 𝐹 𝑧 (𝑦) cos 𝜙 𝑒 𝑑𝑦 𝑏 2 𝑦 𝑙𝑜𝑎𝑑
The local torque is also considered through two contributions. The first one resulting from the natural pitching moment of the airfoil, the second due to the distance between the aerodynamic centre and the elastic axis. To simplify the computations, the aerodynamic centre is assumed to be located at 25% of the local chords independently of the spanwise location.

c. Fuel loads

The fuel distribution is also considered for internal loads computation. The wing is divided into three tanks (center, inner, outer) as illustrated in Figure 22. The 20% outermost part of the wing is considered empty. Besides, the filling of each tank can be controlled individually through filling coefficients (𝑘 𝑓𝑖𝑙𝑙,𝑐𝑡 , 𝑘 𝑓𝑖𝑙𝑙,𝑖𝑡 , 𝑘 𝑓𝑖𝑙𝑙,𝑜𝑡 ) that range from 0.0 to 1.0. The fuel density 𝜌 𝑓𝑢𝑒𝑙 is set to 0.78kg.m -3 . The local shear due to fuel is then straightforwardly obtained from:

𝐹 𝑧,𝑓𝑢𝑒𝑙 (𝑦 𝑙𝑜𝑎𝑑 ) = ∫ -𝑛 𝑧 • 𝑔 • 𝑤 𝑓𝑢𝑒𝑙 (𝑦)𝑑𝑦 𝑏/2 𝑦 𝑙𝑜𝑎𝑑
And the local bending moment is:

𝑀 𝑥,𝑓𝑢𝑒𝑙 (𝑦 𝑙𝑜𝑎𝑑 ) = ∫ 𝐹 𝑧,𝑓𝑢𝑒𝑙 (𝑦) cos 𝜙 𝑒 𝑑𝑦 𝑏/2 𝑦 𝑙𝑜𝑎𝑑
At this stage the torsion due to the fuel is not considered, the underlying hypothesis is that the fuel center of gravity locus is close to the elastic axis.

d. Engine loads

In the new developed model, the loads due to the engine are, to a certain extent, considered. Indeed, only the local shear and bending moment due to the weight are taken into account. While torques resulting from weight and thrust are not computed. The internal engine loads are then given by:

𝐹 𝑧,𝑒𝑛𝑔𝑖𝑛𝑒 (𝑦 𝑙𝑜𝑎𝑑 ) = { 0 𝑖𝑓 𝑦 𝑙𝑜𝑎𝑑𝑠 > 𝑦 𝑒𝑛𝑔𝑖𝑛𝑒 -𝑛 𝑧 • 𝑔 • 𝑊 𝑒𝑛𝑔𝑖𝑛𝑒 𝑒𝑙𝑠𝑒 𝑀 𝑥,𝑒𝑛𝑔𝑖𝑛𝑒 (𝑦 𝑙𝑜𝑎𝑑 ) = { 0 if 𝑦 𝑙𝑜𝑎𝑑𝑠 > 𝑦 𝑒𝑛𝑔𝑖𝑛𝑒 -𝑛 𝑧 • 𝑔 • 𝑊 𝑒𝑛𝑔𝑖𝑛𝑒 (𝑦 𝑒𝑛𝑔𝑖𝑛𝑒 -𝑦 𝑙𝑜𝑎𝑑 ) else

e. Structural weight loads

The structural weight distribution itself is the source of internal loads. This distribution can be either estimated as elliptic or be computed from the results of the structural sizing. The latter induces a loop between the loads computation and the structure sizing model that is solved in FAST-OAD.

Considering the weight distribution 𝑤 𝑠𝑡𝑟𝑢𝑐𝑡 (𝑦), the shear and bending moment are expressed straightforwardly as follows:

𝐹 𝑧,𝑠𝑡𝑟𝑢𝑐𝑡 (𝑦 𝑙𝑜𝑎𝑑 ) = ∫ -𝑛 𝑧 • 𝑔 • 𝑤 𝑠𝑡𝑟𝑢𝑐𝑡 (𝑦)𝑑𝑦 𝑏 2 𝑦 𝑙𝑜𝑎𝑑 𝑀 𝑥,𝑠𝑡𝑟𝑢𝑐𝑡 (𝑦 𝑙𝑜𝑎𝑑 ) = ∫ 𝐹 𝑧,𝑠𝑡𝑟𝑢𝑐𝑡 (𝑦) cos 𝜙 𝑒 𝑑𝑦 𝑏 2 𝑦 𝑙𝑜𝑎𝑑

f. Strut loads

In the case a strut-braced wing configuration is considered, the loads introduced by the strut within the wing must be accounted for. Because the problem is hyperstatic, the traction in strut in sizing manoeuvre conditions is assumed to represents a portion 𝑘 𝑙𝑖𝑓𝑡 of the total aircraft lift. It is then hypothesised that the strut could be preconstrained to reach this force in manoeuvre conditions.

The resulting internal loads are expressed as follows:

𝐹 𝑧,𝑠𝑡𝑟𝑢𝑡 (𝑦 𝑙𝑜𝑎𝑑 ) = { 0 if 𝑦 𝑙𝑜𝑎𝑑 > 𝜂 𝑠𝑡𝑟𝑢𝑡 • 𝑏 2 -𝑘 𝑙𝑖𝑓𝑡 𝐿 sin 𝛼 𝑠𝑡𝑟𝑢𝑡 else 𝑀 𝑥,𝑠𝑡𝑟𝑢𝑡 (𝑦 𝑙𝑜𝑎𝑑 ) = ∫ 𝐹 𝑧,𝑠𝑡𝑟𝑢𝑡 (𝑦) cos 𝜙 𝑒 𝑑𝑦 𝑏/2 𝑦 𝑙𝑜𝑎𝑑
The normal compressive force introduced in the wing must also be taken into account as it will influence spar flanges sizing as detailed in section 2 and may also be responsible for buckling of the inner wing.

𝐹 𝑦,𝑠𝑡𝑟𝑢𝑡 (𝑦

𝑙𝑜𝑎𝑑 ) = { 0 if 𝑦 𝑙𝑜𝑎𝑑 > 𝜂 𝑠𝑡𝑟𝑢𝑡 • 𝑏 2 -𝑘 𝑙𝑖𝑓𝑡 𝐿 cos 𝛼 𝑠𝑡𝑟𝑢𝑡 else
The torque due to the strut is not considered in the current implementation, meaning that the junction with the wing is considered to be located on the elastic axis.

g. Loads envelop

Once all contributions have been computed, the envelop used for the structural sizing can be issued for the considering load case: 𝐹 𝑧 (𝑦 𝑙𝑜𝑎𝑑 ) = 𝐹 𝑧,𝑎𝑒𝑟𝑜 (𝑦 𝑙𝑜𝑎𝑑 ) + 𝐹 𝑧,𝑓𝑢𝑒𝑙 (𝑦 𝑙𝑜𝑎𝑑 ) + 𝐹 𝑧,𝑒𝑛𝑔𝑖𝑛𝑒 (𝑦 𝑙𝑜𝑎𝑑 ) + 𝐹 𝑧,𝑠𝑡𝑟𝑢𝑐𝑡 (𝑦 𝑙𝑜𝑎𝑑 ) + 𝐹 𝑧, 𝑠𝑡𝑟𝑢𝑡 (𝑦 𝑙𝑜𝑎𝑑 ) 𝑀 𝑥 (𝑦 𝑙𝑜𝑎𝑑 ) = 𝑀 𝑥,𝑎𝑒𝑟𝑜 (𝑦 𝑙𝑜𝑎𝑑 ) + 𝑀 𝑥,𝑓𝑢𝑒𝑙 (𝑦 𝑙𝑜𝑎𝑑 ) + 𝑀 𝑥,𝑒𝑛𝑔𝑖𝑛𝑒 (𝑦 𝑙𝑜𝑎𝑑 ) + 𝑀 𝑥,𝑠𝑡𝑟𝑢𝑐𝑡 (𝑦 𝑙𝑜𝑎𝑑 ) + 𝑀 𝑥, 𝑠𝑡𝑟𝑢𝑡 (𝑦 𝑙𝑜𝑎𝑑 ) 𝑀 𝑦 (𝑦 𝑙𝑜𝑎𝑑 ) = 𝑀 𝑦,𝑎𝑒𝑟𝑜 (𝑦 𝑙𝑜𝑎𝑑 ) + 𝑀 𝑦,𝑒𝑛𝑔𝑖𝑛𝑒 (𝑦 𝑙𝑜𝑎𝑑 ) + 𝑀 𝑦, 𝑠𝑡𝑟𝑢𝑡 (𝑦 𝑙𝑜𝑎𝑑 )

A typical example of internal loads and envelops for a strut-braced wing configuration is provided in Figure 23:

Analytical wing structure sizing and weight estimate

Once the loads have been properly assessed, the wing primary structure can be sized. The wings and struts primary structure weights are computed using physical analytical models based on beam theory [8], while the weight of the other parts (fuselage, tails, systems, engines, …) are assessed through the semi-empirical formula already implemented in FAST-OAD.

For the new physics-based models developed here, the wing structure is simplified to an equivalent spar plus skin model with the spar flanges supporting bending moment, the web supporting the shear and the skin supporting the torsion (Figure 24). This methodology has been extended to the struts considering they only support traction. The models developed here are only valid for isotropic metallic materials. For those materials we consider the tensile yield stress 𝜎 𝑦𝑖𝑒𝑙𝑑,𝑡𝑒𝑛𝑠 , the compressive yield stress 𝜎 𝑦𝑖𝑒𝑙𝑑,𝑐𝑜𝑚𝑝 , the maximum shear stress 𝜏 𝑚𝑎𝑥 , the density 𝜌 and a minimum technological thickness for metallic sheets 𝑡 𝑚𝑖𝑛 , whose typical value is around 2mm.

a. Spar flange sizing

As already mentioned, spar flanges are supposed to support all the bending moment. As a consequence, in the case of a fully stressed structure and taking into account the minimum technological thickness, the surfaces of the upper and lower flanges are given by:

𝑆 𝑢𝑝𝑝𝑒𝑟-𝑓𝑙𝑎𝑛𝑔𝑒 (𝑦) = max (𝑡 𝑚𝑖𝑛 • 𝑐 𝑏𝑜𝑥 (𝑦), | -𝑀 𝑥 (𝑦) ℎ 𝑠𝑝𝑎𝑟 (𝑦) • 𝜎 𝑦𝑖𝑒𝑙𝑑,𝑐𝑜𝑚𝑝 |) 𝑆 𝑙𝑜𝑤𝑒𝑟-𝑓𝑙𝑎𝑛𝑔𝑒 (𝑦) = max (𝑡 𝑚𝑖𝑛 • 𝑐 𝑏𝑜𝑥 (𝑦), | 𝑀 𝑥 (𝑦) ℎ 𝑠𝑝𝑎𝑟 (𝑦) • 𝜎 𝑦𝑖𝑒𝑙𝑑,𝑡𝑒𝑛𝑠 |)
The introduction of a strut implies to consider also the normal compressive force (for positive load factors) acting on flanges:

𝑆 𝑢𝑝𝑝𝑒𝑟-𝑓𝑙𝑎𝑛𝑔𝑒 (𝑦) = max (𝑡 𝑚𝑖𝑛 • 𝑐 𝑏𝑜𝑥 (𝑦), | -𝑀 𝑥 (𝑦) ℎ 𝑠𝑝𝑎𝑟 (𝑦) • 𝜎 𝑦𝑖𝑒𝑙𝑑,𝑐𝑜𝑚𝑝 + 0.5 • 𝑇 𝑦 (𝑦) 𝜎 𝑦𝑖𝑒𝑙𝑑,𝑐𝑜𝑚𝑝 |) 𝑆 𝑙𝑜𝑤𝑒𝑟-𝑓𝑙𝑎𝑛𝑔𝑒 (𝑦) = max (𝑡 𝑚𝑖𝑛 • 𝑐 𝑏𝑜𝑥 (𝑦), | 𝑀 𝑥 (𝑦) ℎ 𝑠𝑝𝑎𝑟 (𝑦) • 𝜎 𝑦𝑖𝑒𝑙𝑑,𝑡𝑒𝑛𝑠 + 0.5 • 𝑇 𝑦 (𝑦) 𝜎 𝑦𝑖𝑒𝑙𝑑,𝑡𝑒𝑛𝑠 |)
The weight of the flanges for both swept wings can be easily derived:

𝑊 𝑓𝑙𝑎𝑛𝑔𝑒𝑠 = 2 • 𝜌 cos 𝜙 𝑒 • ∫ (𝑆 𝑢𝑝𝑝𝑒𝑟-𝑓𝑙𝑎𝑛𝑔𝑒 + 𝑆 𝑙𝑜𝑤𝑒𝑟-𝑓𝑙𝑎𝑛𝑔𝑒 )𝑑𝑦 𝑏 2 0 b.

Spar web sizing

The web is supposed to support only the shear. Because expected thicknesses are small enough, the shear flow can be assumed to be constant within the web and the minimal surface allowing a fully stressed structure is given by:

𝑆 𝑤𝑒𝑏 (𝑦) = max (𝑡 𝑚𝑖𝑛 • ℎ 𝑠𝑝𝑎𝑟 (𝑦), 𝐹 𝑧 (𝑦) 𝜏 𝑚𝑎𝑥 )
The weight of the web is then obtained straightforwardly:

𝑊 𝑤𝑒𝑏 = 2 • 𝜌 cos 𝜙 𝑒 • ∫ 𝑆 𝑤𝑒𝑏 𝑑𝑦 𝑏 2 0 c.

Wingbox skin sizing

The torsion shear flow within the wingbox skin is expressed as:

𝜙(𝑦) = - 𝑀 𝑦 (𝑦) 2𝑆 𝑏𝑜𝑥
Then the minimum allowable thickness for the skins satisfies: A correction factor 𝑘 𝑙 is introduced to take into account the non-rectangular shape of the box in the computation of perimeter. A typical value for this coefficient is 0.97.

𝑡 𝑠𝑘𝑖𝑛 (𝑦) ≥ |𝜙 (

d. Strut sizing

The strut is assumed to support only traction force (compression case is eluded at this stage of development). The chord of the strut being considered as constant, the surface to sustain loads is given by:

𝑆 𝑠𝑡𝑟𝑢𝑡 = max (𝑡 𝑚𝑖𝑛 • 2 • (𝑐 𝑏𝑜𝑥,𝑠𝑡𝑟𝑢𝑡 + ℎ 𝑏𝑜𝑥, 𝑠𝑡𝑟𝑢𝑡 ), 𝑘 𝑙𝑖𝑓𝑡 • 𝐿 𝜎 𝑦𝑖𝑒𝑙𝑑,𝑡𝑒𝑛𝑠 ) 20
And the weight for both side struts is directly derived from the length 𝑙 𝑠𝑡𝑟𝑢𝑡 of a strut:

𝑊 𝑠𝑡𝑟𝑢𝑡 = 2 • 𝜌 • 𝑆 𝑠𝑡𝑟𝑢𝑡 • 𝑙 𝑠𝑡𝑟𝑢𝑡

e. Total wing weight

Finally, the total wing weight is computed adding the contribution of skin, flanges and web and considering also ribs and secondary parts through empirical formulations. The ribs are supposed to be evenly spaced spanwise with a constant thickness fixed by the user. The secondary structure is computed with the following formula:

𝑊 𝑤𝑖𝑛𝑔,𝑠𝑒𝑐 = 0,3285 • 𝑘 𝑤𝑖𝑛𝑔 • 𝑀𝑇𝑂𝑊 0,35 • 𝑆 𝑐𝑎𝑛𝑡𝑖𝑙𝑒𝑣𝑒𝑟 • 𝑘 𝑚𝑣𝑜
With 𝑘 𝑤𝑖𝑛𝑔 and 𝑘 𝑚𝑣𝑜 that are respectively a correction coefficient depending on engine layout (its value is 1 for 4 engine aircraft, 1.05 for two engines and 1.1 for rear engines) and a "cultural" coefficient to take into account structural additional weights. 𝑆 𝑐𝑎𝑛𝑡𝑖𝑙𝑒𝑣𝑒𝑟 is the cantilevered surface of the wing (outside fuselage).

Verification of the analytical SBW structural weight model

As a verification of the response of the SBW structural weight estimation model, parametric studies of the response of the model to different parameters such as the aspect ratio of the wing, the spanwise location of the strut junction on the wing, strut traction load and wing box thickness have been conducted. These exercises were conducted for a total aerodynamic loading on the wing assumed to be constant and elliptically distributed along the span.

The impact of the spanwise location of the strut on the wing and of the traction force in the strut (expressed as a fraction of the total wing aerodynamic load, L) is illustrated in Figure 25 for a wing of aspect ratio fixed at 16. This clearly shows, for each strut traction value, an optimum strut attachment location and the global optimum of this parameter being around 65% of the span, from a pure structural weight perspective.

The investigation of the effect of the aspect ratio on the minimal structural wing+strut weight for the SBW is illustrated in Figure 26 (orange line in the bottom left plot). This figure also shows the weight of a cantilever wing structure, sized with the same approach and hypothesis, for different aspect ratios. The potential of the SBW to alleviate the structural wing weight increase with AR is clearly demonstrated by a much slower increase of the wing weight with aspect ratio for the SBW. In these last results, the wing weight was minimised by identifying the optimal value of strut junction location and strut traction. The bottom right plot of Figure 26 shows the evolution of the optimal values of these two parameters with aspect ratio.

B. High fidelity structural analyses

The high-fidelity structural model used to develop the design methodology is the FEM model of the ALBATROS project illustrated in Figure 27. The particularity of this model is to use a curved strut in order to mitigate buckling problems. Static computations of this strut braced wing configuration have been made using the FEM code MSC Nastran in linear and non linear approaches (see Figure 28 and Figure 29). Taking the aerodynamic forces into account in these computations will be the next step of this work.

VI. Results of the OAD process applied to investigate the potential benefits of HAR-SBW concept

A. Aircraft designs evaluations using FAST-OAD process

One of the main interests of the FAST-OAD process lies in the ability to quickly evaluate different aircraft designs. This opportunity is particularly given by the possibility to implement custom modules to evaluate the performances of unconventional configurations. To illustrate this feature, the reconstruction of a cantilever equivalent of the reference SBW configuration starting from the A321-LR configuration is detailed in the following.

Starting from the reconstruction of the A321-LR, the first modification, motivated both by the aerodynamics and the structural disciplines, consisted in moving the engines from their position under the wing to a rear position on the fuselage. In the same time, the geometry of the tail is modified from a conventional tail to a T-Tail geometry, preserving the effectiveness of the HTP. These modifications are shown in Figure 30.

Figure 30 -Geometry modifications modelled by FAST-OAD after rear positioning of the engines

Additional geometry modifications lead to the reproduction of the SBW reference configuration ALBATROS, illustrated in green in Figure 31.

Figure 31 -Reproduction of the ALBATROS SBW configuration with FAST-OAD

As illustrated in Figure 31, the wing geometry is changed between conf 1 (in blue) and conf 2 (in red) by modelling a simple trapezoidal wing with an aspect ratio 𝐴𝑅 = 16. The sweep angle of the configuration is also reduced to match the sweep angle of the ALBATROS configuration, which is supposed to extend the laminar portion of the wing. Finally, a generic strut is added to this last configuration in the configuration 2.1 (in green). The mass breakdown of these configurations is presented in Figure 32.

Figure 32 -Mass breakdown comparison for the reproduced ALBATROS configuration

The interest of the FAST-OAD process resides in the possibility to quickly evaluate and compare cantilever and SBW configurations. From this preliminary analysis, the strut would allow reducing the fuel needed by about 12% to achieve the Breguet-Leduc mission that was defined.

B. Fast evaluation of the aspect ratio influence on the SBW configuration using FAST-OAD

From the positive results of the previous paragraph, next step consists in starting to optimize the SBW configuration by sweeping over the AR parameter.

Several geometries of the SBW configuration are then considered, starting from the swept configuration with an AR = 16 introduced in part A of this section. The FAST-OAD analysis is performed for the aspect ratios 𝐴𝑅 ∈ [START_REF] Drela | [END_REF][START_REF] Park | Geometry Modeling for Unstructured Mesh Adaptation[END_REF][START_REF] Destarac | Drag/Thrust Analysis of Jet-Propelled Transonic Transport Aircraft: Definition of Physical DragComponents[END_REF][START_REF] Karman | Automatic Unstructured Mesh Generation with Geometry Attribution[END_REF][START_REF] Gn | A robust Feasible Directions algorithm for design synthesis[END_REF][START_REF] Ko | Transonic Aerodynamics of a Wing/Pylon/Strut Juncture[END_REF] in combination with the position of the junction between the strut and the wing written 𝜂 𝑠𝑡𝑟𝑢𝑡 = 𝑦 𝑗𝑢𝑛𝑐𝑡𝑖𝑜𝑛 𝑏 . For this analysis, the Breguet-Leduc equation is considered to determine the fuel weight needed by the aircraft to achieve the mission.

The fixed main geometrical parameters of the SBW configuration are the strut wetted surface and the aspect ratio while the variables of interest are the wing span and the wing planform area. The main quantities of interest determined in the end by the FAST-OAD process are the MTOW, block fuel and wing mass, which are expected to decreased if compared with a classical cantilever configuration. The evolutions of these quantities are illustrated in the following figures, where the values plotted for 𝜂 𝑠𝑡𝑟𝑢𝑡 = 0 corresponds to the cantilever configuration of the airplane: As can be seen from Figure 33, the presence of a strut presents a direct advantage compared to a cantilever configuration from a mass reduction perspective, and this for all the aspect ratio considered. From these images, it seems that an optimum position for the junction between the strut and the wing is appearing between 40% and 50% of the wingspan, depending on the aspect ratio. This particularity will be more precisely studied in the future.

VII.Conclusion and Perspectives

In the context of the EU funded Clean-Sky 2 project U-HARWARD [6], a multi-fidelity investigation of the SBW configuration has been initiated in 2020 by ONERA and ISAE-SUPAERO. A multi-disciplinary and multi-level design approach has been set up combining an overall aircraft conceptual design framework, FAST-OAD [START_REF] David | From FAST to FAST-OAD: An open source framework for rapid Overall Aircraft Design[END_REF] and highfidelity disciplinary analysis and design frameworks for aerodynamics and structures. The logic of this multi-level approach is to feed the conceptual design process with physics-based information from high-fidelity analyses and design tasks from the aerodynamics and structural disciplines.

The disciplinary models initially embedded in the FAST-OAD conceptual design framework were based on empirical and semi-analytical formulas that provide a consistent response to overall aircraft parameters for conventional tube-and-wing aircraft configurations. To investigate SBW configurations, these FAST-OAD legacy models have been complemented or substituted by an analytical structural sizing model for what concerns wing weight estimation and a strut drag penalty model derived from high-fidelity CFD simulations of SBW configurations.

The OAD process embedding these new models for the wing structure and aerodynamic specifically developed for SBW configuration has been applied to perform first investigations of the SBW concept. Although the first results of these OAD investigations include some conservative hypotheses and therefore do not yet take full benefits of the SBW concept (e.g. wing airfoil thickness reduction and natural laminar flow potential on the wing are not accounted for at this stage), the potential benefits of a HAR-SBW aircraft configuration aircraft are clearly identified, compared to a conventional tube-and cantilever-wing aircraft designed for the same mission corresponding to the one of a Airbus A321-LR. A parametric investigation of the effect of the wing aspect ratio showed an optimum benefit in term of block-fuel reduction for a value of the aspect ratio around 20. For such high aspect ratio values, thanks to the strut, the OEW is even reduced compared to a conventional tube-and-wing aircraft thanks to the snowball effects and therefore the MTOW is significantly reduced.

This SBW concept investigation will be consolidated and further elaborated in the future, with the exploration of additional wing design parameters such as the wing thickness distribution, or by accounting for the potential drag reduction benefits achievable by natural laminar flow on the low sweep, thin wing enabled by the SBW concept.

VIII. Acknowledgment

Figure 1 -

 1 Figure 1 -Reduction of inner wing flexural moment by the SBW concept.

Figure 2 -

 2 Figure 2-FAST-OAD process (XDSM)

Figure 3 -

 3 Figure 3 -Aerodynamic design framework based on CAPS/ESP/EGADS tools

Figure 4 -Figure 5 -Figure 6 -

 456 Figure 4 -Geometrical parameters associated to the fuselage

Figure 7 -

 7 Figure 7 -Reference aircraft geometry CFD simulations based on both Euler and RANS (Reynolds Averaged Navier-Stokes) modelling are achieved on this configuration for several Mach numbers. The analysis can be performed either in the near field, for instance as illustrated in Figure 8 with the surface pressure coefficients

Figure 8 -

 8 Figure 8 -Pressure coefficient distribution on the reference aircraft at M = 0.78, RANS (Spalart-Allmaras) modelling of the flow

Figure 9 -

 9 Figure 9 -Influence of the Mach number on the drag decomposition for the reference aircraft, CFD computations achieved with RANS modelling of the flow These results are considered as the starting point of the present study to conclude on the effect of the strut on the configuration performances.

  Z-direction of the junction between the strut and the fuselage Strut_emp_Y m Position in the spanwise-direction of the junction between the strut and the fuselage Strut_curvature_r adius m Radius joining the "vertical" and "horizontal" components of the strut Strut_Vert_sweep ° Sweep angle of the "vertical" part of the strut Strut_LE_sweep ° Sweep angle of the "horizontal" part of the strut Junction_TR -Ratio applied to the wing chord at the junction between the strut and the wing Strut_Vert_twist ° Twist angle of the "vertical" part of the strut Strut_twist ° Twist angle of the "horizontal" part of the strut

Figure 10 -

 10 Figure 10 -Illustration of the geometrical parameters associated to the strut

Figure 11 -

 11 Figure 11 -Reconstruction of the ONERA ALBATROS configuration using ESP/EGADS tools CFD computations are achieved at cruise conditions of the mission defined earlier using Euler flow modelling at a Mach number M = 0.78 and for several target lift coefficients. The Far Field Drag post-processing included in the process allows identifying the sources of the different drag components, as illustrated in Figure 12 for the wave drag.

Figure 12 -

 12 Figure 12 -Wave drag sources based on Euler simulations on the reference Strut-Braced Wing configuration

Figure 16 -Figure 17 -Figure 18 -

 161718 Figure 16 -History of evolution of airfoil geometry variables during the aerodynamic drag minimization by gradient-based optimisation

  length. This enabled straightforward parametric investigation of the impact of engine position and settings (toe and pitch angle) on the overall aircraft aerodynamic coefficients, as it is illustrated in Figure19. Each of the configurations corresponding to underwing and rear-fuselage mounted engine locations could be analysed by just changing the values defining the position, pitch and toe angles of the nacelle in the geometry model.

Figure 19 -

 19 Figure 19 -Illustration of different variants of engine integration investigated with the automated CFD-based analysis process

  also shows the selected convention of the internal forces orientation applied for loads computation.

Figure 21 -

 21 Figure 21-Strut-braced wing parametrization and internal forces orientation.

  𝑒𝑎 -0.25) cos 𝜙 𝑒 ∫ 𝑐 𝑎𝑒𝑟𝑜 𝑙(𝑦)𝑑𝑦 𝑏 2 𝑦 𝑙𝑜𝑎𝑑

Figure 22 -

 22 Figure 22-Wing tanks layout.Finally, the fuel distribution is given by:

Figure 23 -

 23 Figure 23-Typical internal loads and envelop for structural sizing for a Strut-braced wing configuration.

Figure 24 -

 24 Figure 24-Wing primary structure model for weight computation [8].

  𝑦)| 𝜏 𝑚𝑎𝑥 ≡ |-𝑀 𝑦 (𝑦) 2𝑆 𝑏𝑜𝑥 (𝑦) • 𝜏 𝑚𝑎𝑥 | With: 𝑆 𝑏𝑜𝑥 ≡ 𝑆 𝑟𝑖𝑏𝑠 = 𝑘 𝑆𝑟𝑖𝑏𝑠 × 𝑒 𝑎𝑖𝑟𝑓𝑜𝑖𝑙 × 𝑐 𝑎𝑒𝑟𝑜 2 × (𝑘 𝑟𝑒𝑎𝑟-𝑠𝑝𝑎𝑟 -𝑘 𝑓𝑟𝑜𝑛𝑡-𝑠𝑝𝑎𝑟 ) × cos(𝜙 𝑒 ) × 𝑓 𝜙 𝑒 A correction coefficient 𝑘 𝑆𝑟𝑖𝑏𝑠 that takes into account the non-rectangular shape of the box in the ribs surface (enclosed by skin) is introduced. Its typical value is around 0.93. Finally, the local surface of the wingbox skin can be simply derived by the multiplication of skin thickness and wingbox perimeter p. 𝑆 𝑠𝑘𝑖𝑛 (𝑦) = 𝑝(𝑦) • 𝑡 𝑠𝑘𝑖𝑛 (𝑦) ⇔ 𝑆 𝑠𝑘𝑖𝑛 (𝑦) = 2 • 𝑘 𝑙 • [𝑒 𝑎𝑖𝑟𝑓𝑜𝑖𝑙 (𝑦) + (𝑘 𝑟𝑒𝑎𝑟-𝑠𝑝𝑎𝑟 (𝑦) -𝑘 𝑓𝑟𝑜𝑛𝑡-𝑠𝑝𝑎𝑟 (𝑦)) 𝑓 𝜙 𝑒 cos 𝜙 𝑒 ] • 𝑐 𝑎𝑒𝑟𝑜 (𝑦) • |-𝑀 𝑦 (𝑦) 2𝑆 𝑏𝑜𝑥 (𝑦) • 𝜏 𝑚𝑎𝑥 |

Figure 25 -

 25 Figure 25-Investigation of the impact of spanwise strut/wing junction location on the primary wing (flange), strut and total wing+strut weights

Figure 27 -Figure 29 -

 2729 Figure 27-Finite Element Model of ALABATROS strut braced wing

  a) Influence of the position of the junction on the MTOW b) Influence of the position of the junction on the block fuel mass c) Influence of the position of the junction of the wing mass

Figure 33 -

 33 Figure 33 -Preliminary study of the influence of the Aspect Ratio of the SBW configuration

Table 1 -Weights from FAST-OAS sizing for CSR-01

 1

Table 1 . width_fus Heigth_fus c_wing_tip ϕ_wing_LE Span_wing c_wing_root c_HTP_root c_HTP_tip Span_HTP c_VTP_root c_VTP_tip ϕ_VTP_LE Span_VTP Parameter Units Value

 1

	Fuselage Length	m	44.51
	Fuselage Width	m	3.95
	Fuselage Height	m	4.05
	Wing Span	m	34.1
	Wing Area	m²	126
	Wing Aspect Ratio	-	9.23
	Wing Taper Ratio	-	0.25
	HTP Span	m	12.45
	HTP Area	m²	30.75
	HTP Aspect Ratio	-	5.04
	HTP Taper Ratio	-	0.37
	VTP Span	m	5.87
	VTP Area	m²	22.3
	VTP Aspect Ratio	-	1.55
	VTP Taper Ratio	-	0.26

Table 1

 1

This choice was found to be a good compromise between aircraft size (that could enable to derive conclusions for SMR and even regional mission) and scalability up to longer ranges, where the induced drag savings should provide the best positive effect.Additionally to above requirements, engines are fixed with a Maximum TakeOff thrust of 134kN.

Ass. Prof., ISAE-SUPAERO-DCAS Department.

The work presented in this paper has been funded by the European Union in H2020-EU.3.4.5.10 through the CleanSky2 project U-HARWARD, under grant Grant agreement ID: 886552

IX. References


Chawki Habchi
email: chawki.habchi@ifpen.fr

A GIBBS ENERGY RELAXATION (GERM) MODEL FOR CAVITATION SIMULATION

Keywords: cavitation, turbulence, two-fluid model, stiffened gas equation of state

In this paper, a comprehensive highly compressible and turbulent two-fluid multispecies model is presented. It involves an equation for the transport of the liquid volume fraction in addition to two different sets of partial differential equations for the gas and the liquid phase. The multicomponent gas phase is governed by an ideal gas equation of state (EOS) while the stiffened gas EOS is speci-fied to the single-component liquid phase. In this work, a Reynolds averaged Navier-Stokes (RANS) formulation is adopted. For the turbulence, a standard k -ε model is used for the gas phase while a turbulent viscositybased model is used for the liquid phase in order to improve the laminar-turbulent transition computation through the nozzle. In addition, the model equations include different relax-ation terms for mass, momentum, and energy exchanges at the liquid-gas interfaces. For the present cavitation modeling, an instantaneous relaxation procedure is used for the velocity, pressure, and temperature; while a slower procedure is adopted for the Gibbs free-energy relaxation model (GERM) at the interfaces. These models have been applied for the simulation of the cavitation inside a trans-parent single-hole nozzle. The obtained cavitation pocket has a similar shape as the experiments. Moreover, two different cavitation regimes have been identified. A gaseous cavitation regime appears in a region in which the static pressure is close to but above the liquid saturation pressure; a second cavitation regime may happen when the static pressure goes below the liquid saturation pressure. In the latter case, the liquid become superheated and leads to a vaporous cavitation regime. Also, the cooling of the fuel and the density variation due to the expansion of compressible liquid through the nozzle is among the more interesting findings of this paper.

INTRODUCTION

With increasingly stricter emission regulations and greater demand on fuel economy, the injector has become one of the most critical components of modern diesel engines. Since the nozzle geometry influences the spray characteristics and the mixture preparation in the engine, the design of injector nozzles has become an important factor for the improvement of the combustion performance and pollutant reduction. In addition, the physical processes involved in direct injection engines are complex. Consequently, it is important to understand the effects of the nozzle geometry (sac, orifice, etc.) on 1 NOMENCLATURE A 0 switch (0 or 1) for the turbulence µ, ρ viscosity, density models in Eqs. ( 2) and (3) Π equal to (P + A 0 2/3ρk) C k C k = 0.006 is a constant φ, ω, θ, ψ interfacial velocity, in Eq. ( 9) pressure, temperature and C s C s = 0.2 is a constant in Eq. ( 10) Gibbs relaxation to C µ C µ = 0.09 is a constant in Eq. ( 8) equilibrium characteristic C v liquid specific heat at constant frequency volume coefficient in Eq. ( 7) D diffusion coefficient in Eq. ( 5) Sub-and Superscripts e specific internal energy f fuel g specific free Gibbs energy i index of the Cartesian h specific enthalpy components k turbulent kinetic energy (TKE) I interfacial parameter l length scale l liquid P pressure L laminar condition P ∞,l , e ∞,l SG-EOS Eq. ( 7) fluid parameters m index of the gas species Q specific heat flux in Eq. (1) R perfect gas universal constant p, q q = 3 -p, index of the S ij , S strain tensor and strain intensity phases in Eqs. ( 1 the liquid injection, atomization, and combustion behavior, especially as the orifice size keeps getting smaller and the injection pressure higher. Moreover, the description of turbulent flow and its coupling with liquid atomization and dispersion is required, as well as the description of the mixing, the ignition, and the combustion itself. An accurate prediction of liquid injection is then critical in this context, and has to be modeled from the in-nozzle region up to the dispersed evaporating spray region, where combus-tion takes place. A number of studies have been conducted on the effect of the nozzle characteristics on the internal and external spray characteristics [START_REF] Andriotis | Influence of vortex flow and cavitation on near-nozzle diesel spray dispersion angle[END_REF][START_REF] Blessing | Analysis of flow and cavitation phenomena in diesel injection nozzles and its effects on spray and mixture formation[END_REF][START_REF] Soteriou | Direct injection diesel sprays and the effect of cavitation and hydraulic flip on atomization[END_REF]. Even if experiments of the complete process are a source of information and physical understanding, there is, at the present time, a growing willingness for the development of predictive numerical computations of such flows, which is particularly challenging for regions at the inner and near injector, as shown in [START_REF] Dumont | Cavitating flow in Diesel injectors and atomization: A bibliographical review[END_REF], [START_REF] Habchi | Modélisation Tridimensionnelle de l'injection et de la préparation du mélange gaz/carburant dans les moteurs combustion interne[END_REF], and [START_REF] Moreau | A numerical study of cavitation influence on diesel jet atomization[END_REF]. Currently, most cavitation models are based either on a single-fluid mixture approach using a homogeneous equilibrium model (HEM) [START_REF] Habchi | Multidimensional simulation of cavitating flows in diesel injectors by a homogeneous mixture modeling approach[END_REF][START_REF] Lee | A numerical investigation of transient flow and cavitation within minisac and valve-covered orifice diesel injector nozzles[END_REF][START_REF] Schmidt | A fully compressible, two-dimensional model of small, high-speed, cavitating nozzles[END_REF]Wang et al., 2011) or a simplified Rayleigh-Plesset (RP) bubble dynamics equation [START_REF] Giannadakis | Modelling of cavitation in diesel injector nozzles[END_REF]. In the HEM model, the cavitation is computed using a barotropic equation of state (EOS) derived from the Tait's correlation (Wang et al., 2011) or the Wood's speed of sound [START_REF] Wood | A Textbook of Sound[END_REF]. In the RP models, the cavitation is computed as the result of the pressure difference between the inner bubble and the surrounding liquid. Besides, it is usually assumed that cavitation should take place in regions in which the static pressure falls below the saturation pressure of the flowing liquid. Although this condition is not necessary for the growing and collapse of cavitation bubbles, both of the HEM and RP models assume that cavitation bubbles or pockets are full of vapor and do not consider gaseous cavitation which is observed experimentally [START_REF] Dunn | Experimental characterization of aviation-fuel cavitation[END_REF]. In addition, high velocity shears and stretching, pressure, or temperature variations that one may encounter in two-phase flows formed inside and in the close vicinity of injectors cannot be resolved by such HEM models, as only mass and momentum equations for the mixture are accounted for. Moreover, both HEM and RP models ignore the thermodynamic processes happening at the liquid-gas interfaces during phase changes. A more sophisticated mixture model has been recently developed by Schmidt and co-workers [START_REF] Shields | Cavitation as rapid flash boiling[END_REF]. The phase change is computed using a homogeneous relaxation model (HRM) which assumes a first-order rate equation for the evolution of the instantaneous vapor mass fraction toward its equilibrium value over a given empirical timescale. In a study performed by [START_REF] Battistoni | Comparison of mixture and multifluid models for in-nozzle cavitation prediction[END_REF], this mixture-HRM model has proved to be less accurate than a multifluid model using a RP bubble dynamic for phase change. In addition, Battistoni numerical investigations of gaseous cavitation have highlighted the importance of accounting for dissolved gases in the liquid. While these models have proved to be useful in many applications, they cannot be universally applied to all fuel injection cases, especially for highly compressible and non-isothermal conditions. More recently, a new family of Eulerian two-phase flow models has appeared. In these models, both liquid and gas are compressible and obey different set of transport equations for the mass, momentum, and energy in addition to an equation for the transport of the volume fraction. Presently, at least two different types of Eulerian models are under study in the literature. The most general model is the "seven equation model" of Saurel and co-workers (2008), from which a "five equation model" may be derived by assuming the pressures and velocities of the two phases to be equal. [START_REF] Wang | CFD simulation of fuel injection using a five equation Eulerian twophase flow model[END_REF] have applied the latter model to the case of fuel injection. They have shown that gaseous near-wall cavitation can significantly increase in the nozzle orifice due to mechanical expansion of dissolved gas, without the help of mass transfer. In this work, a new cavitation model based on the equilibrium at liquid-gas interface of the Gibbs free energy of liquid and gas is suggested in a general framework of highly compressible and turbulent seven equation model. This paper is organized as follows: First, the comprehensive highly compressible and turbulent two-fluid model is presented in Section 2, along with the GERM cavitation model. Next, the numerical implementation of these models in the IFP-C3D software [START_REF] Bohbot | IFP-C3D: An unstructured parallel solver for reactive compressible gas flow with spray[END_REF][START_REF] Velghe | A high efficiency parallel unstructured solver dedicated to internal combustion engine simulation[END_REF] is described in Section 3. The GERM cavitation model is then applied in Section 4 for the simulation of the cavitation inside a transparent single-hole nozzle for which experimental cavitation results are available [START_REF] Blaisot | NADIA bio: New advanced diagnosis for diesel injection analysis and bio fuels -A federative project around hydraulics of CR systems[END_REF]. The present study also attempts to improve the understanding of the different cavitation regimes by varying the temperature of the injected fuel from 300 to 450 K. The main results of this study are discussed in Section 5. Finally, the conclusion of this work is reported with a future work.

MATHEMATICAL MODELING

Flow regimes in internal and near-injector regions can be handled through different Eulerian/Eulerian approaches [START_REF] Herrmann | A parallel Eulerian interface tracking/Lagrangian point particle multi-scale coupling procedure[END_REF][START_REF] Lebas | Numerical simulation of primary break-up and atomization: DNS and modelling study[END_REF]. One may consider direct interface resolution using Navier-Stokes equations in each phase and jump conditions at the interfaces. The liquid jet atomization and the cavitation inception modeling are possible with sufficiently refined meshes and high-order numerical schemes in order to be able to capture the interface position accurately. The minimum liquid or gas structure size that can be accounted for is directly related to the grid resolution. Direct interface resolution approaches can be used for understanding studies and validation purposes, limited to simple geometries. Yet, they are not appropriate for complex configurations such as internal combustion engines and multihole injector simulations. In addition, since these approaches are usually derived using an incompressible liquid assumption, they are not able to properly compute the cavitation inside the nozzle, where the liquid velocity may have a value as high as 700 ms -1 and a Mach number greater than 0.5.

The Two-Fluid Model

In this paper, a comprehensive highly compressible and turbulent two-fluid multispecies model is presented. It is based on the previous work of [START_REF] Baer | A 2-phase mixture theory for the deflagration-to-detonation transition (DDT) in reactive antigranulocytes-materials[END_REF] model and [START_REF] Kapila | Two-phase modeling of deflagration-to-detonation transition in granular materials: Reduced equations[END_REF] as well as the more recent work of Saurel and co-workers (2008).

Their method is based on a compressible model involving seven partial differential equations, three for the gas phase and three for the liquid phase in addition to an equation for the liquid volume fraction transport. In this work, a multispecies model is also considered in order to be able to investigate non-condensable gas (such as N 2 ) in addition to condensable vapor. Besides, the equations are written in the Reynolds averaged Navier-Stokes (RANS) framework. The details of the averaging procedure may be found in [START_REF] Devassy | Atomization Modeling of Liquid Jets using an Eulerian-Eulerian Model and a Surface Density Approach[END_REF]. In addition, the specific internal energy equation is used instead of the total specific energy adopted in the work by [START_REF] Saurel | A multiphase Godunov method for compressible multifluid and multiphase flow[END_REF] and [START_REF] Saurel | Modelling phase transition in metastable liquids: Application to cavitating and flashing flows[END_REF]. Also, this two-fluid model is able to deal with liquid jet atomization as suggested by [START_REF] Devassy | Atomization modelling of liquid jets using a twosurface-density approach[END_REF]. Phase change at the interface is also considered in this paper.

For two-phase flows, the model consists of the following governing equations:

∂α p ρm,p ∂t + ∂α p ρm,p V p,i ∂x i = ∂ ∂x i ρp D p ∂ ∂x j ρm,p ρp + ψ (g q -gp ) δ mf (1) ∂ ∂t α p ρ p V p,i + ∂ ∂x i α p ρ p V p,i V p,j = - ∂α p Π p ∂x i + Π I ∂α p ∂x i + ∂α p τ L,T p,ij ∂x i + ϕ V q,i -V p,i (2)
∂ ∂t (α p ρ p e p ) + ∂ ∂x i α p ρ p V p,i e p = -α p P p ∂ V p,i ∂x i - ∂α p Q L,T p,i ∂x j + (1 -A 0 )α p τ L,T p,ij ∂ V p,i ∂x j + A 0 α p ρ p ε p + ϕ V I V q,i -V p,i + ωP I P q -P p + θ( T q -T p ) + ẽI ψ (g q -gp ) (3) ∂α p ∂t + V I i ∂α p ∂x i = ω P q -P p + θ κ ( T q -T p ) + ψ ρI (g q -gp ) (4)
Here p = {1, 2} and q = 3 -p; say p = 1 for the liquid phase (index l) and p = 2 for the gas phase (index g). The overbar (•) and tilde (•) designate the Reynolds and Favre average operation, respectively. Equation ( 1) is the continuity equation for species m.

The f th species index designates the fuel vapor. For the gas, a Fick law with a single turbulent diffusion coefficient D p is assumed, while the liquid is assumed single-component (D 1 = 0). In the system of Eqs. ( 1)-( 4) also intervenes interfacial quantities denoted with a superscript "I." The interfacial velocity, V I , and the interfacial pressure, P I , are assumed according to the [START_REF] Baer | A 2-phase mixture theory for the deflagration-to-detonation transition (DDT) in reactive antigranulocytes-materials[END_REF] model (V I = V l and P I = P g ). In Eq. ( 2), Π is equal to (P +A 0 2/3ρk), where A 0 = 1 for turbulent flows. τ L,T p,ij is the shear stress tensor in Eqs. ( 2) and (3). It can be written as τ

L,T p,ij = τ L p,ij + A 0 τ T p,ij
, where the superscripts L and T stand for laminar and turbulent flow regimes, respectively. Its modeling is given below in Section 2.3. In the energy equation ( 3), e p is the Favre averaged specific internal energy and Q L,T p,i the heat flux which is modeled as follows:

Q LT p,i = -λ p ∂ T p ∂x i -ρp D p m h m,p ∂ ∂x i ρm,p ρp (5)
The first term in the right-hand side (RHS) of Eq. ( 5) is due to the conduction (Fourier's law). For both phases, the turbulent contribution of the conductivity coefficient λ p is computed from a specified turbulent Prandtl number taken equal to 0.9. However, only the liquid has a significant laminar conductivity coefficient contribution for λ p . This contribution is computed using the Wilke and Lee correlation [START_REF] Reid | The Properties of Gases and Liquids[END_REF]. The second term in the RHS of Eq. ( 5) represents the total flux of the enthalpy transported by the gaseous species due to turbulent advection. This term is not considered in the liquid since this phase is assumed single component.

Equations of State

Equations of state (EOS) are used to close the system of Eqs. ( 1)-( 4). In order to circumvent the difficulty of models having negative squared speed of sound in the two-phase region, the present study uses two EOSs. Indeed, each fluid possesses a separate EOS.

In the present investigation, terms arising in the multispecies gas phase are closed by the perfect gas EOS where R is the universal constant and W g is the gas mixture molecular weight:

P g = ρ g R W g T g (6)
In addition, the stiffened gas EOS (SG-EOS) is used in the single-species liquid phase. Two different expressions of the SG-EOS are given as follows:

P l = ρ l (γ l -1) (e l -e ∞,l ) -γ l P ∞,l , T l = P l + P ∞,l C v,l ρ l (γ l -1) (7)
Table 1 summarizes the values of the SG-EOS parameters for the NF82BN18 fuel used below. This fuel is a mixture of 82% vol Normafluid and 18% vol Bromonaphthalene used in the optical investigation by [START_REF] Blaisot | NADIA bio: New advanced diagnosis for diesel injection analysis and bio fuels -A federative project around hydraulics of CR systems[END_REF]. This mixture allows concealing the fluid-PMMA interface, especially in the orifice where the cavitation occurs. Details concerning the SG-EOS and its parameters may be found in Le [START_REF] Métayer | Elaborating equations of state of a liquid and its vapor for two-phase flow models[END_REF].

Stress Tensors and Turbulence Modeling

The stress tensor τ L,T p,ij = τ L p,ij +A 0 τ T p,ij requires modeling in order to close the equations system (1)-(4). A standard Boussinesq approximation is used for the modeling of the Reynolds stress tensor (τ T p,ij = -ρ p V ′′ p,i V ′′ p,j ). As such, the total shear stress tensor can be written as

τ L,T p,ij = µ p ∂ V p,i ∂x j + ∂ V p,j ∂x i - 2 3 µ p ∂ V p,i ∂x i δ ij (8)
where µ p is the sum of the laminar and turbulent dynamic viscosities (µ p = µ L p +A 0 µ T p ). The laminar viscosity of gas phase, µ L g , is calculated according to Sutherland's formula, a function of the gas temperature. The laminar viscosity of liquid phase, µ L l , is given in a table as a function of the liquid temperature. The turbulent contribution, µ T g is given using a standard k g -ε g model [START_REF] Amsden | A computer program for chemically reactive flows with sprays[END_REF]) by

µ T g = ρg C µ k 2 g ε g (9)
with C µ = 0.09. For the liquid phase, a simple algebraic k l -l model based on a characteristic length, l taken of the order of the of grid cell, l = (Vol cell ) 1/3 has been used,

µ T l = ρ l (C s l) 2 S ( 10
)
where C s is a modeling constant taken equal to 0.2 and S = 2 Sij : Sij is the mean strain rate tensor intensity with Sij = 1/2{∂ V l,i /∂x j + ∂ V l,j /∂x i }. In addition, k l is modeled according to the relation k l = C k l 2 S2 where C k is taken equal to 0.006. Regarding turbulence modeling, the standard kε model and the Re-Normalization Group (RNG) kε model have also been implemented for the liquid phase in addition to the gas phase. Using these models in conjunction with the suggested two-phase "seven equation" model has proved to be very unstable in the liquid phase, especially during transient flow regimes. The reason for that behavior is that the kε model type is only valid for high Reynolds number conditions. The suggested turbulent viscosity based model is a first attempt to compute the needle opening flows and the nozzle flows including transition from laminar flow in the sac toward turbulent flow at the exit of the nozzle. In future work, further improvements in the turbulence modeling are required, maybe using the Large Eddy Similation (LES) approach.

GERM CAVITATION MODEL

The two-fluid model described in the previous section includes different terms for the relaxation of the velocity, the pressure, the temperature, and the Gibbs free energy at liquid-gas interfaces. The details of the relaxation methods are not given in this article because they are available elsewhere. However, the manner in which the relaxation terms are considered and numerically implemented is described in this section. In this work on the cavitation problem, the velocity and pressure in Eqs. ( 2)-( 4) are relaxed with an infinite rate by making φ and ω respectively, tend to infinity. The description of these relaxation procedures can be found in Saurel and Abrall (1999) and [START_REF] Saurel | A multiphase model for compressible flows with interfaces, shocks, detonation waves and cavitation[END_REF]. Similarly, an infinite relaxation of the temperature at the interfaces is applied in Eqs. ( 3) and ( 4) by tending the relaxation parameter θ to infinity (see [START_REF] Zein | Modeling phase transition for compressible two-phase flows applied to metastable liquids[END_REF]. Note that the pressure and temperature equilibrium are accomplished by computing the exact volume fraction variations [see the RHS of Eq. ( 4)]. These volume fraction changes are happening without any phase change since they are absent from the RHS of the mass conservation Eq. ( 1). Finally, the last terms in Eqs. ( 1) and ( 4) represent, respectively, the mass flow rate and volume changes due to phase change at the interfaces. These quantities have been computed using a Gibbs free-energy relaxation model (GERM). Unlike the model of [START_REF] Zein | Modeling phase transition for compressible two-phase flows applied to metastable liquids[END_REF], the GERM model adopts a finite relaxation rate for the Gibbs free energy at the interface. Indeed, in order to ensure the stability of the explicit scheme, the numerical procedure limits the vaporized mass and its corresponding heat transfer using a volume fraction Courant-Friedrichs-Lewy (CFL) number CFL α = (α eq,icellα icell )/α icell , where α eq,icell is the equilibrium gas volume fraction obtained by an infinite relaxation procedure for the Gibbs free energy. This means that the maximum variation of α icell is equal to CFL α α icell for each cell icell of the grid. A value of CFL α = 10 -3 has been used in this paper. A similar procedure may be found in [START_REF] Zein | Modeling phase transition for compressible two-phase flows applied to metastable liquids[END_REF] along with the expression for the parameters κ, ρI , and ẽI appearing in the RHS of Eqs. ( 3) and ( 4). As such, the GERM model may simulate different cavitation regimes. Indeed, the first two terms, ω P q -P p + (θ/κ)( T q -T p ), in Eq. ( 4) are dedicated for the modeling of the so-called "gaseous cavitation." This first regime of cavitation may appear in regions in which the static pressure is close to but above the liquid saturation pressure. In this case, not only the disequilibrium of pressure but also the disequilibrium of temperature between the liquid and the gas are responsible for the growth or collapse of the bubbles without the help of any phase change. In addition, a second cavitation regime may happen when the static pressure goes below the liquid saturation pressure. In this case, the liquid becomes superheated and leads to a "vaporous cavitation" regime. This vaporous regime is modeled by the last term, ṁp = ψ (g q -gp ) in Eq. ( 1) and its corresponding energy flux and volume fraction in Eqs. ( 3) and (4), respectively.

NUMERICAL METHODS

The two-fluid model along with the GERM cavitation model described in the previous sections has been implemented in the software IFP-C3D. The details of description of the numerical methods and how the equations are coupled and solved for the gas phase are detailed in [START_REF] Bohbot | IFP-C3D: An unstructured parallel solver for reactive compressible gas flow with spray[END_REF] and [START_REF] Velghe | A high efficiency parallel unstructured solver dedicated to internal combustion engine simulation[END_REF]. Since the liquid phase subsystem (1)-( 3) is similar to the gas one, the same numerical schemes have been used. Besides, the numerical schemes for the liquid volume fraction resolve Eq. ( 4) in the following modified form:

∂α p ∂t + ∂α p V I i ∂x i = α p ∂ V I i ∂x i + ω P q -P p + θ κ ( T q -T p ) + ψ ρI (g q -gp ) (11)
This modification allows using the same numerical scheme for the convection (i.e., rezoning) term as for the mass, momentum, and energy equations. Equation ( 11) is solved only for p = 1 (i.e., for the liquid volume fraction, α l = α 1 ). Next, the relaxation terms have been solved at the end of each time step, as described in the previous section. Finally, note that the transport equations in IFP-C3D are written in the arbitrary Lagrangian Eulerian (ALE) formalism in order to allow the description of moving components such as the injector needle.

COMPUTATIONAL CONFIGURATION

In this section, the GERM model is applied to the simulation of the cavitation inside a single-hole nozzle, for which experimental cavitation results are available [START_REF] Blaisot | NADIA bio: New advanced diagnosis for diesel injection analysis and bio fuels -A federative project around hydraulics of CR systems[END_REF]. The main features of the optical nozzle used in this work are given in Table 2 and Fig. 1. Besides, a simple linear interpolation of the vertices' coordinate between two grids is used for the needle movement. The grids correspond to the maximum and minimum needle lifts. A minimum needle lift of 10 µm is used as shown in Fig. 2. The three-dimensional mesh shown in Fig. 3 includes the nozzle and the combustion chamber. Only half of the 3D configuration has been meshed by 600,000 hexahedron cells. The plane of symmetry is shown in Fig. 3. A typical size of 60 mm has been chosen for the diameter and height of the combustion chamber. In order to make the specification of the inlet pressure boundary condition easier, the nozzle mesh has been extended to reach the inlet command (i.e., control) chamber of the injector. In addition, a relatively rough mesh is used in the combustion in order to limit the simulation cost. Note that the needle lift region is initially meshed by 10 hexahedron cells with a thickness lower than 1 µm. Furthermore, the cell size in the hole is taken relatively small (10 µm). They are even smaller at the orifice wall (3 to 4 µm) in order to resolve possible secondary flows in the cavitation inception regions. The initial and boundary conditions are also given in

FIG. 2:

Needle lift (ASI stands for "after start of injection").

Fig. 3. Initially, the whole injector is full of liquid while the combustion chamber is full of gas (nitrogen). Also, due to the two-fluid modeling, each phase contains a weak volume fraction α ε of the other phase. For the cavitation modeling, α ε of gas in the liquid is taken of the order of the dissolved gas value in typical fuels (α εg = 10 -2 ). However, a much smaller value of liquid volume fraction α εl in the gas is necessary and a value α εl = 10 -6 is adopted from previous works [START_REF] Saurel | A multiphase model for compressible flows with interfaces, shocks, detonation waves and cavitation[END_REF][START_REF] Zein | Modeling phase transition for compressible two-phase flows applied to metastable liquids[END_REF]. The injection pressure is 25 MPa. This value is also the initial liquid pressure in the injector above the needle seat region. However, the liquid filling the sac and the hole is set to 0.1 MPa, which is the same pressure as the gas phase in the combustion chamber. The temperature of the gas in the combustion chamber is set to 300 K, while the injected fuel has been varied from 300 to 450 K. Also, the computations are assumed to start with a zero velocity and very small values of turbulent kinetic energy (k g = 10 -12 m 2 /s 2 ) and dissipation rate (ε g = 10 -12 m 2 /s 3 ) in the computational domain. In addition, weak turbulent kinetic energy and dissipation rate values have been specified at the inlet. Finally, all the walls are considered with a slip condition.

RESULTS AND DISCUSSION

Two computations have been carried out using two different injected liquid temperatures: 300 and 450 K. For the low-temperature case, shadow images are available from the NADIA-BIO project [START_REF] Blaisot | NADIA bio: New advanced diagnosis for diesel injection analysis and bio fuels -A federative project around hydraulics of CR systems[END_REF]. Figure 4 tive agreement between the numerical and the experimental cavitation regions inside the optical nozzle. However, the computed cavitation pocket seems to be larger and closer to the orifice inlet than in the experimental image. Indeed, the cavitation pocket shape seems to be sensitive to the geometrical details at the orifice inlet (see the nozzle geometrical uncertainties reported in Table 2). In Fig. 5, the numerical results of the two computations are compared in terms of cavitation regimes. For the 450 K case, the cavitation has appeared faster than in the Identification of the gaseous and vaporous cavitation regimes. The presence of vapor inside the cavitation regions is depicted by its mass fraction iso-lines. 300 K case. In addition, there is no phase change for the 300 K case due to the extremely low saturation pressure (less than 10 Pa) of the fuel at 300 K (see Fig. 6). However, since the static pressure in the orifice has decreased to a value close to but above the liquid saturation pressure (Fig. 7), a gaseous cavitation regime has appeared during the injection period. Since the saturation pressure is significant at 450 K (greater than 10 kPa; see Fig. 6) and the pressure has decreased below 10 kPa in the orifice (see Fig. 7), a vaporous cavitation regime has appeared during the needle opening, until 330 s. In this case, the fuel vapormass fraction (shown in Fig. 5 by iso-lines) is up to 0.6 inside the cavitation pockets. Nevertheless, the results of both temperature cases are quite similar in terms of cavitation pocket size during the needle closing period. But, this is not the case during the injection opening transient, where the cavitation is expected to have a significant effect on the jet primary atomization. Also, the results are different in term of velocity magnitude as shown in Fig. 8. This is due to the decrease of liquid density with temperature increase. This topic is very important and will be among future studies. Finally, one may observe in Fig. 9 the variation of the liquid density as a function of the temperature in the sac and the orifice during injection for the 450 K case. Upon the needle opening, a shock wave develops due to the liquid pressure difference between the seat and the sac. The compression induced by the shock wave leads to an increase of liquid density and temperature, especially in the bottom of the sac. Next, a cooler fuel fills the sac progressively. Finally, the needle closing induces a rarefaction wave and leads to a decrease of liquid density and temperature in the hole, then in the sac. Although the temperature variation is relatively small (DT = ±10 K), the liquid density changes (±5%) could affect the injector flow rate significantly. This study shows that it is possible to reproduce the appearance and collapse of cavitation without using adjustable parameters. In addition, two regimes of cavitation have been identified and simulated by the GERM model. It has been shown that cavitation bubbles do not always contain vapor. Indeed, for the gaseous cavitation, it is not necessary that static pressure decreases below the liquid saturation pressure as assumed by most previous models. Also, the cooling of the fuel and the density variation due to the expansion of compressible liquid through the sac and the orifice is among the more interesting findings of this paper. Future works will attempts to improve the understanding of the different cavitation regimes and their effects on atomizing liquid jets. For instance, the influence of wall-heat transfer on the cavitation should be investigated along with wall-slip or no-slip conditions. In addition, more quantitative comparisons with accurate experimental data are necessary. Particularly, these comparisons will be based on the recent work performed by [START_REF] Manin | Boundary condition and fuel composition effects on injection processes of diesel sprays at the microscopic level[END_REF] and [START_REF] Kastengren | Engine combustion network (ECN): Measurements of nozzle geometry and hydraulic behavior[END_REF][START_REF] Kastengren | Timeresolved X-ray radiography of sprays from engine combustion network spray A diesel injectors[END_REF] in the framework of the Engine Combustion Network (ECN).

FIG. 3 :

 3 FIG. 3: 3D Mesh, initial, and boundary conditions.

FIG. 4 :

 4 FIG. 4: Qualitative comparison of the numerical and experimental cavitation regions. (a) Shadow image obtained in the framework of NADIA-BIO project. (b) A volume fraction iso-value equal to 0.5 is used for delimiting cavitation region.

  FIG. 5: (Color online) Liquid volume fraction and vapormass fraction distributions. Identification of the gaseous and vaporous cavitation regimes. The presence of vapor inside the cavitation regions is depicted by its mass fraction iso-lines.

FIG. 6 :

 6 FIG. 6: Pressure saturation (Psat) curve of the NF82NB18 fuel.

  FIG. 8: (Color online) Velocity field distributions. The cavitation region is depicted by an iso-surface of liquid volume fraction equal to 0.5. Cavitation region is located above this iso-surface. Velocity profiles in the orifice are also shown with the arrows.

  FIG. 9: (Color online) Temperature and liquid density fields in the sac and the hole.

TABLE 1 :

 1 SG-EOS parameters for the NF82BN18 fuel

	P ∞,l	γ l	C v,l	e ∞,l
	5.3 × 10 8 Pa	3.48	761 J/kg/K	-1004 kJ

TABLE 2 :

 2 Nozzle geometric features

	Units Average value Uncertainties

FIG. 1: Sac and orifice offset.

ACKNOWLEDGMENTS

The author would like to thank the French ANR agency for its support of this work in the framework of the NADIA-BIO project (contract no. 143886). The author is indebted to Professor E. Daniel from Aix-Marseille University for his valuable advice and to Mr. B. Reveillé (IFPEN) for providing the grid used in this article. Also, the author would like to thank Mr. R. Kumar (IFPEN) for his help with grid generation.