Inês Da Costa Cardoso

Gaspard Berthelin

Sylvain Dubreuil

Michel Salaün

Nathalie Bartoli

Christian Gogu

I Cardoso
email: ines.cardoso@onera.fr

Model Order Reduction and Bayesian Optimization for MDO problems

Keywords: Multidisciplinary analysis and optimization, model order reduction, surrogate models, Gaussian processes, overall aircraft design

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Multidisciplinary Analysis (MDA) and Optimization (MDO) are useful tools when addressing problems which involve several disciplines or subsystems. For instance, in aircraft design, the aerodynamic and structural disciplines must be taken into account in order to fully characterize the aeroelastic behaviour of a wing [START_REF] Grossman | Integrated aerodynamic/structural design of a sailplane wing[END_REF][START_REF] Jasa | Open-source coupled aerostructural optimization using Python[END_REF]. With the aid of disciplinary solvers, it is possible to model the physical phenomena, as well as the coupling between them. These solvers then interact with one another via the coupling variables, resulting in a non-linear system of equations. To solve this coupled system, each disciplinary solver must be evaluated iteratively by a global solver until a given convergence criterion is reached. The role of MDA is thus to find the equilibrium of the system at a given design space point.

Multidisciplinary Optimization, on the other hand, may be used to iteratively change the design parameters towards a given goal. Different architectures are available for multidisciplinary optimization. The Individual Disciplinary Feasible (IDF) and Multidisciplinary Feasible (MDF) approaches attempt to find the optimal point via a single optimization problem [START_REF] Martins | Multidisciplinary Design Optimization: A Survey of Architectures[END_REF][START_REF] Cramer | Problem Formulation for Multidisciplinary Optimization[END_REF]. In the IDF approach, the optimization algorithm handles both the design and coupling variables, by uncoupling the disciplinary solvers and instead imposing consistency constraints. Despite the reduction in the computational cost achieved by avoiding the resolution of the coupled MDA, this formulation can become difficult to set up if the number of coupling variables is significant. Contrarily, in the MDF approach, the coupled MDA is solved via fixed point iteration, usually by means of a non-linear Gauss-Seidel solver. This formulation is easier to set up and returns a physically feasible design at every optimization iteration. Furthermore, the high cost of solving the coupled MDA may be improved if gradient information is considered, for instance, by employing Newton-based solvers.

The computational burden of solving the multidisciplinary analysis might be further reduced by constructing surrogate models of the objective and constraint functions [START_REF] Jones | Efficient Global Optimization of Expensive Black-Box Functions[END_REF][START_REF] Bartoli | An adaptive optimization strategy based on mixture of experts for wing aerodynamic design optimization[END_REF], thus limiting the number of whole MDA resolutions. However, without an adequate enrichment criterion, it might become necessary to have accurate surrogate models throughout the entire design space, even in areas where the local or global optima are unlikely to be. Similarly, to improve the computational time of the multidisciplinary optimization, gradient-based algorithms may be employed. Nonetheless, these are designed to converge to local minima and may require attempting different starting points in order to find a global optimum.

An alternative strategy for solving the parametrized MDA, denominated Disciplinary Proper Orthogonal Decomposition and Interpolation (DPOD+I), is proposed in [START_REF] Berthelin | Disciplinary Proper Orthogonal Decomposition and Interpolation for the resolution of parametrized Multidisciplinary Analysis[END_REF], where the issue of using high dimensional coupling variables is addressed by uncoupling the disciplinary solvers and replacing each one by a reduced order model [START_REF] Chatterjee | An introduction to the proper orthogonal decomposition[END_REF][START_REF] Zimmermann | Gradient-enhanced surrogate modeling based on proper orthogonal decomposition[END_REF]. These reduced order models are then trained using Gaussian processes, which replace the disciplinary solvers in the multidisciplinary analysis. The surrogate disciplines are enriched iteratively throughout the optimization process by means of sensitivity analysis, and according to the methodology proposed in [START_REF] Dubreuil | Towards an efficient global multidisciplinary design optimization algorithm[END_REF]. As other surrogate-based approaches, this methodology avoids the resolution of the MDA via the potentially costly solvers, but requires enrichment of the surrogates models only when the uncertainty is high, thus resulting in a lower number of expensive solver calls. To address global optimization, the authors of [START_REF] Dubreuil | Towards an efficient global multidisciplinary design optimization algorithm[END_REF] propose an iterative procedure denominated Efficient Global Multidisciplinary Design Optimization (EGMDO). This approach enriches the disciplinary models solely in the areas of the design space where the global optimum is more likely to be. To do so, an Expected Improvement criterion, such as the one introduced in [START_REF] Jones | Efficient Global Optimization of Expensive Black-Box Functions[END_REF], is proposed, based on the uncertainty induced by the disciplinary surrogates on the objective function. By iteratively enriching the disciplinary models in the relevant design space areas, this uncertainty is reduced and the global optimum is found.

The present article proposes to combine the replacement of the disciplinary solvers by their respective surrogates obtained by means of DPOD+I, with the use of EGMDO procedure as optimization algorithm. The proposed framework is then used to solve an aeroelastic optimization problem associated with an aircraft wing, by coupling a potential flow model with a linear elastic structural model. An inverse problem is chosen as a test case, where the goal is to retrieve the set of design parameters which results in the desired reference solution.

The present work is organized as follows. Section 2 presents the chosen test case as well as its resolution via the computationally expensive disciplinary solvers. This will serve as a reference solution against which the remaining implementations will be compared. Section 3 will address the replacement of the disciplinary solvers by their DPOD+I surrogates, while Section 4 will combine the DPOD+I with the EGMDO optimizer. Finally, some conclusions and perspectives are given.

Static aeroelastic optimization of an aircraft wing

We consider the aeroelastic analysis and optimization problem associated with an aircraft wing, whose coupled behaviour is described by the aerodynamic and structural disciplines. Together, these disciplines define the following non-linear system of equations:

   u s = M s (x, f a ) f a = M a (x, u s ) (1)
Here, u s is the structural displacement of the wing, f a is the vector of aerodynamic forces acting upon the wing, M s and M a are, respectively, the chosen structural and aerodynamics models and x is the set of the design variables. Together, u s and f a define the set of coupling variables considered in our problem. As will be seen in the following section, they are defined as high-dimensional vectors.

Disciplinary solvers

To solve the MDA it is necessary to model the participating disciplines. In this work, the disciplinary solvers are assumed to be black-box solvers, meaning that no intrusive methods can be used. Additionally, we assume that the derivatives of their outputs with respect to the design variables are not available.

The aerodynamics solver is then chosen as the linear, potential flow model, solved by Vortex Lattice Method (VLM) [START_REF] Katz | Low-Speed Aerodynamics[END_REF]. The VLM models the surface of the wing as a thin plate, neglecting both the thickness and the viscosity of the wing. If the wing geometry is discretized into a structured mesh, it is possible to compute the aerodynamic loads by modelling its surface using horseshoe vortices. Here, the considered planform geometry of the wing is that of the Common Research Model (configuration uCRM-9) [START_REF] Brooks | Undeflected Common Research Model (uCRM): An Aerostructural Model for the Study of High Aspect Ratio Transport Aircraft Wings[END_REF] and the aerodynamic mesh discrezation is such that the load vector f a is of dimension d a = 2100. A representation of the aerodynamic mesh is provided in Figure 1(a).

The structural discipline is modelled via linear elastic theory, and solved by the finite element solver Code Aster [START_REF] Edf | Finite element code_aster , Analysis of Structures and Thermomechanics for Studies and Research[END_REF], on which a static linear analysis is conducted to determine the wing's displacements. The finite element model is composed of the upper and lower skins, 3 spars and evenly spaced ribs, all modelled by thin plate elements. In addition, each of the skins is reinforced by 9 stiffners, modelled by beam elements. In total, the model consists of d s = 43416 degrees of freedom, corresponding to the dimension of u s . The structural mesh is presented in Figure 1(b), where the upper skin was removed to facilitate visualization. In addition to modelling the disciplines, it is necessary to transfer the loads and displacements between the non-coincident structural and aerodynamic meshes. This is made by radial basis interpolation [START_REF] Rendall | Unified fluid-structure interpolation and mesh motion using radial basis functions[END_REF], through which an interpolation matrix H ∈ R d a ×d s is built. This matrix allows to obtain the displacements of the aerodynamic nodes as u a = Hu s and the structural loading as f s = H T f a .

Problem statement

To facilitate comparison between the different implementations, an inverse problem may be defined, for which 4 parameters are chosen as design variables. Two of these parameters, namely the angle of attack α and the freestream velocity V ∞ , influence only the aerodynamic analysis. Contrarily, the remaining parameters impact solely the structural model. These are the skin thicknesses t sk , assumed to be the same for both lower and upper skins, and the spar thicknesses t sp , considered to be the same for all three spars. The thickness of the wing ribs and stiffners is kept constant throughout the resolution of the problem. Table 1 presents the range of variation of each of the considered design variables. It is worth noting, however, that for a more practical implementation, these variables have been scaled to take values in the range [0, 1].

Variable α[ • ] V ∞ [m/s] t sk [m] t sp [m] [1, 9] [220, 250] [0.003, 0.01] [0.01, 0.1]
Table 1: Unscaled design variable range of variation As previously stated, an inverse problem was chosen. This means that, throughout the optimization process, the goal is to find the set of design parameters which corresponds to a given reference solution obtained at x = x ref .

The objective function is then described by the sum of the relative error committed by each of the coupling variables at any evaluated design point x, and the optimization problem may be written as,

x * = arg min x∈X f obj (x) with f obj (x) = ∥f a (x ref ) -f a (x)∥ 2 ∥f a (x ref )∥ 2 + ∥u s (x ref ) -u s (x)∥ 2 ∥u s (x ref )∥ 2 (2)
where X ⊂ R 4 is the subspace defined by the design variable bounds. For the remainder of this work, x ref will be chosen as the corner of the design space for which the maximum wing tip displacement is achieved, namely, for the scaled design variables, x ref = {α, V ∞ , t sk , t sp }, which corresponds to the set {1, 1, 0, 0} in scaled design variables. The wing vertical displacement field obtained at the chosen reference design point is presented in Figure 2. As is shown, the structural displacement is maximum at the wing tip, where a vertical displacement of 2.05 meters is observed. Contrarily, and due to the imposed boundary conditions, no displacement is observed at the wing root. Although not shown here, the minimum displacement of the wing is obtained at x = {0, 0, 1, 1} for scaled variables, for which a vertical displacement of 0.067 meters at the wing tip is observed.

Finally, it is worth noting that, in the present work, no constraint was considered and the optimization problem is fully defined by Eq. (2).

Optimization framework

The previously defined aeroelastic problem was implemented within NASA OpenMDAO framework [START_REF] Heath | OpenMDAO: Framework for Flexible Multidisciplinary Design, Analysis and Optimization Methods[END_REF][START_REF] Gray | OpenMDAO: An open-source framework for multidisciplinary design, analysis, and optimization[END_REF], an open-source code for multidisciplinary analysis and optimization, which facilitates gradient computation. OpenMDAO allows the resolution of the non-linear MDA through the non-linear block Gauss-Seidel or Newton solvers, and makes use of SciPy optimizers to solve the MDO. OpenMDAO presents the benefit of a modular implementation, which allows for an easy reformulation of the problem. Furthermore, thanks to ONERA WhatsOpt environment [START_REF] Lafage | WhatsOpt: a web application for multidisciplinary design analysis and optimization[END_REF], which generates the OpenMDAO skeleton code, the implementation of the disciplines is facilitated. Figure 3 presents the implemented optimization framework for the reference case where the real solvers are used to model the disciplines.

Some additional remarks should be made concerning the presented implementation. Firstly, to take into account the high-dimensional coupling variables, only the MDF approach could be used, hence the MDA is solved at each optimization iteration. Secondly, the resolution of the MDA is performed by the nonlinear block Gauss-Seidel solver, for which the use of Aitken acceleration [START_REF] Irons | A version of the Aitken accelerator for computer iteration[END_REF] showed to significantly reduce the number of necessary Gauss-Seidel iterations. Lastly, at the MDO level, the solver used was the gradient-based Sequential Least SQuares Programming (SLSQP) solver, which makes use of OpenMDAO capability to calculate the derivatives of the objective function with respect to the design variables via finite difference approximation. Regarding the convergence criteria, a tolerance of 10 -8 was set for the nonlinear solver, while a tolerance of 10 -6 was set for the optimizer.

Concerning the initial guess, the center of design space was chosen, i.e., x 0 = {0.5, 0.5, 0.5, 0.5}. This point was kept for the DPOD+I implementation to allow for the comparison of the computational cost.

Reference problem solution

For the reference problem, the optimizer was able to converge to the reference solution in 17 iterations. Table 2 presents the obtained x * , the corresponding objective function value, and the number of disciplinary solver calls performed throughout the optimization. Note that while the chosen disciplinary solvers are not particularly expensive, if higher-fidelity, such as CFD, or non-linear solvers were used, their computational cost would be significantly greater.

α * V * ∞ t * sk t * sp f (x * ) n a n s
1.0 1.0 0.0 0.0 0.0 286 286 Table 2: Optimization results for the reference problem, where n a is the number of aerodynamic solver calls and n s is the number of structural solver calls.

As is shown, to converge exactly to the reference solution, a total of 286 solver calls to each disciplinary solver were needed. Of the 286 calls, 33.6% were made at the gradient computation step. This is an expected result, as the derivatives were calculated via finite difference approximation, which requires additional objective function evaluations.

DPOD+I implementation

According to the methodology proposed in [START_REF] Berthelin | Disciplinary Proper Orthogonal Decomposition and Interpolation for the resolution of parametrized Multidisciplinary Analysis[END_REF], a surrogate model is built for each disciplinary solver. However, due to the high dimensionality of the coupling variables, a model order reduction must be performed before the surrogate interpolation. This reduction is carried out by Disciplinary Proper Orthogonal Decomposition (DPOD) and the surrogate interpolation is made by Gaussian Processes to approximate the POD coefficients. This is considered a pre-processing step which takes place before the beginning of the optimization (see [START_REF] Berthelin | Disciplinary Proper Orthogonal Decomposition and Interpolation for the resolution of parametrized Multidisciplinary Analysis[END_REF] for details on the initial training phase).

It remains, nonetheless, possible to enrich the disciplinary surrogate models throughout the optimization. The decision of whether or not to enrich is made via sensitivity analysis, where, in the case where enrichment is deemed necessary, the less accurate discipline is identified and enriched. The adapted optimization framework is presented in Figure 4, where CoV is the coefficient of variation of the objective function and ϵ is a user defined threshold, used to decide whether or not to enrich the disciplinary surrogates. For the current implementation, a threshold of 10 -3 was set, while the convergence criteria for the MDA and MDO were the ones defined for the reference problem. The gradient computation is made using the Gaussian process derivatives.

Lastly, it is worth noting that this implementation introduces randomness in the obtained solution. This is due to the DPOD+I training stage, where the DPOD coefficients are replaced by Gaussian processes, resulting in a corresponding random MDA (once more, readers are referred to [START_REF] Berthelin | Disciplinary Proper Orthogonal Decomposition and Interpolation for the resolution of parametrized Multidisciplinary Analysis[END_REF] for details). For this reason, several runs must be performed so as to estimate the value of each parameter in the solution.

The number of performed runs was set to 10.

DPOD+I optimization results

As previously stated, the construction of the disciplinary surrogates is a pre-processing step and therefore is not a part of the represented framework. Nonetheless, and to allow for a fair comparison, the associated computational cost should be taken into account. Concerning the size of the obtained POD basis, the average of the 10 performed runs resulted in a structural basis composed of 6 coefficients, while the aerodynamic basis was made of 5 coefficients.

The obtained results for the DPOD+I optimization framework are presented in Table 3, for the 10 performed runs. As is shown, the replacement of the disciplinary solvers by their DPOD+I surrogates induces a significant reduction in the overall computational cost of the optimization. Namely, for the aerodynamic discipline, an average of 60.7 solver calls were needed, of which around 52 were made in the training stage. Similarly, for the structural discipline, an average of 51.5 solvers calls were made, of which 42 corresponded to the training stage. In comparison to the reference problem this represents a reduction by a factor of 5 on the number of disciplinary solver calls. We further remark that, since the training phase is independent of the optimization, for a problem with multiple local minima where multi-start strategies are required for gradient-based solvers, the same basis may be kept regardless of the starting point. Although not explored in the present work, this is a clear advantage of the proposed DPOD+I implementation.

Finally we note that, due to the performed approximations, a small error is committed at the optimum point, both at the level of the design variables, at an order of magnitude of 10 -4 , and of the objective function, at an order of magnitude of 10 -2 , as is shown in table 3. Figure 5 compares the number of disciplinary solver calls, as a function of the optimizer iteration, for one of the 10 performed runs of DPOD+I implementation against that of the reference problem. As is shown, the DPOD+I approach starts at a higher number of disciplinary solver calls, due to the training stage. Nonetheless, throughout the optimization not only fewer calls are made, but also a lower number of optimizer iterations is necessary, resulting in the overall lower computational cost.

EGMDO implementation

In addition to replacing the disciplinary solvers by their DPOD+I surrogates, we now consider the Efficient Global Multidisciplinary Design Optimization (EGMDO) algorithm proposed in [START_REF] Dubreuil | Towards an efficient global multidisciplinary design optimization algorithm[END_REF]. Here, we depart from an initial Design of Experiments (DoE) which is used to obtain an approximation of the objective function, itself a random function. Contrarily to what is proposed in the classical Efficient Global Optimization (EGO) algorithm [START_REF] Jones | Efficient Global Optimization of Expensive Black-Box Functions[END_REF], however, in EGMDO the objective function cannot be modelled as a Gaussian Process due to the non-linearity of the MDA. Instead, its value at the DoE points is obtained by means of a Karhunen Loève expansion, and then extended to the rest of the domain via Gaussian process interpolation.

As in the previous implementation, the disciplinary surrogate models may be enriched at a given x i point, if the coefficient of variation of the objective function is greater than a given threshold. However, in the EGMDO approach, an additional condition for enrichment is imposed -the point must have some likelihood of being the minimum, namely, P min (x i ) ≥ 1 n DoE . This avoids spending computational effort in areas of the design space where the global optimum is unlikely to be.

New points may also be added to the DoE. These points are determined by a modified Expected Improvement (EI) criterion, which, due to the non-Gaussian nature of the objective function, has no analytical expression and must be estimated via Monte Carlo simulation (see [START_REF] Dubreuil | Towards an efficient global multidisciplinary design optimization algorithm[END_REF] for details on EI calculation). Regarding the stopping criterion, while in EGMDO the optimization is stopped when the change of maximum value of EI between two consecutive iterations is less than a given threshold, here we impose a maximum of 200 iterations. The computational framework of the EGMDO implementation is presented in Figure 6, where an initial DoE x DoE of 20 points was considered.

Lastly, and before proceeding to the results, it is recalled that, since the combined DPOD+I and EGMDO implementation also introduces randomness in the solution, 10 runs were also performed for the present optimization.

EGMDO optimization results

As for the previous case, the training of the DPOD+I surrogates is considered as a pre-processing step and is not featured in the optimization framework. However, its computational cost is taken into account. The obtained results for the EGMDO implementation are presented in Table 4. As is shown, an average of 61 solver calls were performed for the aerodynamics discipline, of which 59 were made during the DPOD+I surrogate training. Similarly, an average of 51 solver calls were made for the structural discipline, of which 49 corresponded to the pre-processing stage. This signifies that an average of only 2 calls were made to each disciplinary solver during the actual MDO resolution, implying that very few disciplinary solver enrichments were necessary to find the optimal point. In fact, comparing to the DPOD+I implementation, a reduction by a factor of 5 is observed in the number of disciplinary solver enrichments. A possible explanation for why very few enrichments are made is that, at the level of the DPOD+I surrogate training, a too strict convergence criterion was chosen, thus reducing the need for enrichment. For the present problem, it might be possible to further reduce the number of disciplinary solver calls if a less strict criterion was chosen. However, more complex problems will likely benefit from the well trained DPOD+I surrogates.

α * V * ∞ t *
Finally, and as was the case for the DPOD+I implementation, due to the performed approximations, a small error, of the order of 10 -2 , is committed both for the design variables and the objective function, as is evidenced in Table 4.

Conclusion and perspectives

In this work, different frameworks for multidisciplinary optimization were implemented and compared for an aeroelastic problem with two participating disciplines. The objective function was chosen as an inverse problem, where the goal was to find a known reference solution. In addition, no constraint was considered.

The initial implementation, and the most costly, consisted on modelling the disciplines by possibly expensive black-box solvers and solving the optimization by using the gradient-based SLSQP algorithm, in an MDF approach. This implementation, in addition to making several disciplinary solver calls to solve the MDA, requires gradient computation via finite difference approximation, which is itself an expensive computation.

A second implementation replaced the disciplinary solvers by their DPOD+I surrogates, according to the methodology proposed in [START_REF] Berthelin | Disciplinary Proper Orthogonal Decomposition and Interpolation for the resolution of parametrized Multidisciplinary Analysis[END_REF]. This implementation resulted in significantly less disciplinary solver calls, as, throughout the optimization, the computationally expensive solvers were only called when the disciplinary Gaussian processes needed to be enriched. In addition, this approach eliminated the need to perform gradient computation via finite difference approximation, as the Gaussian process derivatives were used. Nonetheless, the training of the DPOD+I surrogates had an associated computational cost which could not be disregarded. As was remarked, however, in a multi-start context, the training phase would be run a single time, as it is independent of the optimization problem.

A last implementation replaced the optimizer by the EGMDO algorithm proposed in [START_REF] Dubreuil | Towards an efficient global multidisciplinary design optimization algorithm[END_REF], in an effort to address global optimization, while keeping the DPOD+I surrogates to approximate the disciplines. This approach modelled the objective function as a non-Gaussian random field, approximated by a combination of Karhunen Loève expansion and Gaussian process interpolation, and used an adapted Expected Improvement criterion to enrich the disciplinary solvers solely where the global minimum was more likely to be. It was shown that, through this approach, very few enrichments were needed in order to find the optimum point, possibly due to the easy to find global optimum.

Having demonstrated the capabilities of the proposed framework, some perspectives are envisioned. Namely, the application of this strategy to more complex problems, such as the ones involving one or multiple non-linear disciplinary solvers, as these could prove to be a challenge for the POD approximation. Similarly, the inclusion of a greater number of design variables could raise a problem at the level of the disciplinary GP interpolation, currently done by the universal kriging model. Fortunately, some work has already been done in this field, from which we may benefit. For instance, in [START_REF] Bouhlel | Improving kriging surrogates of highdimensional design models by Partial Least Squares dimension reduction[END_REF] Partial Least Squares is used for dimension reduction. Finally, it would be interesting to apply the proposed strategy to a problem with multiple local minima, which usually accompanies the aforementioned increase of complexity and dimensionality.

Figure 1 :

 1 Figure 1: Wing meshes used in the disciplinary solvers.

Figure 2 :

 2 Figure 2: Wing vertical displacement field at reference solution.

Figure 3 :

 3 Figure 3: Reference problem optimization framework.

Figure 4 :

 4 Figure 4: DPOD+I optimization framework.

Table 3 :

 3 Optimization results for the DPOD+I implementation. Average and coefficient of variation values obtained from a total of 10 runs. (a) Aerodynamic solver calls (b) Structural solver calls

Figure 5 :

 5 Figure 5: History of disciplinary solver calls for the reference and DPOD+I implementations.

Figure 6 :

 6 Figure 6: EGMDO optimization framework.

Table 4 :

 4 Optimization results for the EGMDO implementation. Average and coefficient of variation values obtained from a total of 10 runs.

			sk	t * sp	f (x * )	n a	n s
	E	0.999 0.987 0.006 3 × 10 -7 0.044	61	51
	CoV 0.003 0.02	0.9	3.16	0.44 0.1117 0.1225

Acknowledgments

This work is a part of the activities of the ONERA-ISAE-ENAC joint research group.


Agathe De La Hougue

Guillaume Androuin

Laetitia Pernod

Matthieu Sacher
email: matthieu.sacher@ensta-bretagne.frd.envsn

Benoît Augier
email: benoit.augier@ifremer.fr

Patrick Bot
email: patrick.bot@ecole-navale.fr

Paul Iachkine
email: paul.iachkine@envsn.sports.gouv.fr

A De

L A Hougue Abc

Influence of shaft stiffness on the performance of a plane-like hydrofoil

Keywords: Hydrofoils, Velocity Prediction Program, Free Surface Interaction, FSI

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

Introduction

Hydrofoils are used to partially or fully support a ship's weight, thus reducing the hydrodynamic resistance as a significant fraction of the ship's hull is lifted out of the water. The use of lifting hydrofoils as a means of drag reduction consequently caused a jump in yacht performance and increased velocities over the recent years. However, as flying crafts are relatively new, these series have not yet reached the same level of expertise as older yacht racing series. In particular, designing hydrofoils for sailing races such as the America's Cup or the Olympics remains a complex challenge since the hydrofoil has to achieve boat flying conditions as much as possible, while i) minimizing drag over a large range of boat speeds and while ii) keeping the craft as stable and manoeuvrable as possible. For small foiling crafts, T-foils or plane-like hydrofoils are now widely used, as they offer a good trade-off between the generated lift, the total drag and the equilibrium of the submersed system. These foils are characterized by a horizontal lifting surface connected to the craft by a vertical shaft. The performance can be studied thanks to Velocity Prediction Programs (VPP) where the equilibrium of the foiling craft is solved for different sailing conditions. VPPs have been developed for different sailing vessels, from classical highperformance boats [START_REF] Philpott | Yacht velocity prediction using mathematical programming[END_REF]] to foiling catamarans [START_REF] Bagué | Dynamic stability analysis of a hydrofoiling sailing boat using CFD[END_REF][START_REF] Graf | VPP-driven sail and foil trim optimization for the Olympic Nacra 17 foiling catamaran[END_REF], including dinghies such as lasers [START_REF] Day | Performance prediction for sailing dinghies[END_REF]]., with an increased complexity in the physical phenomena considered and the numerical methods used. For instance, [START_REF] Day | Performance prediction for sailing dinghies[END_REF]] takes into account the sail power and the position of the crew, and [START_REF] Bagué | Dynamic stability analysis of a hydrofoiling sailing boat using CFD[END_REF][START_REF] Graf | VPP-driven sail and foil trim optimization for the Olympic Nacra 17 foiling catamaran[END_REF]] use numerical methods from lifting line to full CFD. The influence of the structure deformation may also be considered thanks to the introduction of FSI effects in the VPP, such as in the case of foiling 60ft Imoca [START_REF] Horel | Application of system-based modelling and simplified-fsi to a foiling open 60 monohull[END_REF] or L foils on catamarans [START_REF] Graf | Prediction and Optimization of Aerodynamic and Hydrodynamic Forces and Boat Speed of Foiling Catamarans with a Wing Sail and a Jib[END_REF]].

An important parameter in the dynamics of foiling crafts that needs to be considered in Velocity Prediction Programs (VPP) is the height of flight (Hagemeister and Flay 2019), as it plays a critical role on the side force and the drag through the wetted area. Indeed, the ride height is often one of the key factors in the foiling equilibrium and stability. The issue of the optimum height of flight for the sail craft "moth" equipped with T-Foils is studied both experimentally [Kebbell et al., Binns et al.] and numerically with the evaluation of a lifting line method [START_REF] Day | Performance prediction for sailing dinghies[END_REF]].

This work focuses on the behaviour of windfoils (foiling winfsurfs) which are equipped with a planelike foil package composed of a shaft, a main wing and a smaller rear wing called stabilizer, connected by a fuselage (Figure 1). When foiling, the craft flies around 0.7m above the water surface with a speed from 15 up to 30 knots (7.5 up to 15 m/s) and with great fluctuations in attitudes (for example during a "stable" flight, the heel angle may vary by up to 10°). Important issues on the windfoil are thus the effect of the free surface and the deformation of the shaft, as this slender structure is heavily loaded. The interaction with the free surface is indeed an important characteristic of the windfoil with a surface piercing shaft and horizontal lifting surfaces being at low depth or sometime surface piercing. The performance of hydrofoils close to a free surface has been presented in the literature both experimentally and numerically for a single extruded NACA foil [START_REF] Ni | Performance of a Hydrofoil Operating Close to a Free Surface over a Range of Angles of Attack[END_REF][START_REF] Daskovsky | The hydrofoil in surface proximity[END_REF]] or in tandem configuration [START_REF] Bal | High-speed submerged and surface-piercing cavitating hydrofoils, including tandem case[END_REF]] with a diminution of the lift and increase of the drag when getting closer to the surface. no change of efficiency when the heel angle is increased from 0° to 30°. They gave results on the variations of lift and drag with the depth-to-chord ratio. [START_REF] Day | Experimental measurement and simplified prediction of T-foil performance for monohull dinghies[END_REF] performed tank testing on monohull dinghy T-foils and showed that a lifting line approach was correct for low angles of attack. They tested two different submergence heights and found a reduced drag close to the surface. [START_REF] Daskovsky | The hydrofoil in surface proximity[END_REF] found a decrease in lift when the foil gets closer to the free surface for moderate angles of attack. [START_REF] Ni | Performance of a Hydrofoil Operating Close to a Free Surface over a Range of Angles of Attack[END_REF] showed the same trend with a sharp decrease for submergence-to-depth ratio lower than one. The fluid structure interaction (FSI) of foils has been studied by numerous authors. For example, Temtching Temou et al [2021] investigated the hydro-elastic response of a composite hydrofoil with models of different levels of fidelity and highlighted the importance of the bend-twist coupling (see also [START_REF] Pernod | Experimental and Numerical Investigation of the Fluid-Strucure Interaction on a Flexible Composite Hydrofoil under Viscous Flows[END_REF]). However, FSI is rarely considered by the commonly used VPPs. As a notable exception, Horel and Durand (2019) studied a foiling Open60 monohull and found that the foil flexibility had little influence on the lift and drag, but significantly affects the side force and changes how the loading is distributed among the different appendages.

Fs Fs

This study focuses on the performance prediction of a foiling windsurf, namely a windfoil. The sailing craft is composed of a 9m² sail, a 2.2m long board and a plane-like hydrofoil. The hydrofoil, presented on Figure 1, is composed of a plane based on a main wing (0.095m²) and a stabilizer (0.025m²) connected by a fuselage. The plane is connected to the board by the help of a vertical shaft, itself a 1m long freesurface piercing composite structure with a symmetric section.

In the present work, it is considered an upwind sailing point of a rigid body formed by the sail, the board and the athlete (see Figure 2), while the shaft of the foil is the only structure allowed to deform. As shown in Figure 2, the craft is sailing at an apparent wind angle AWA which is obtained from a vector construction of the true wind and the speed V of the craft. The craft is sailing with a heel angle ϴg to windward (opposite to common non-foiling yachts), a small trim angle ϴa which mainly affects the angle of attack of the foil (main wing and stabilizer), and finally a small drift angle ϴd is also considered, which is commonly called the leeway angle, that mainly affects the angle of attack of the submerged portion of the shaft. In order to study the influence of different tuning parameters on the performances of the windfoil, a 5 degree-of-freedom (DOF) VPP has been implemented to solve: 1. Speed V 2. Leeway angle ϴd 3. Heel angle ϴg 4. Height of flight Tz 5. Trim angle ϴa The yaw balance is the only unsolved degree of freedom. To ease the VPP convergence, some parameters remain also fixed. Indeed, the opening angle of the sail (around Z-axis when ϴg=ϴa=0°) is set to 13.5° while its trim angle (around Y-axis when ϴg=ϴa=0°) is set to -10°. It should be noted that since a rigid body (sail, board and athlete) is considered, the craft angles ϴd, ϴg and ϴa also affect the forces and moments of the sail. The sailor position from the reference point in Figure 2 is also set to -0.4m, 1.4m and 1m for respectively X, Y and Z axis when ϴg=ϴa=0°. The considered nonlinear system of forces and moments equations to be solved is given by,

∑ 𝐹 ⃗ = 𝐹 𝑓 ⃗⃗⃗⃗ + 𝐹 𝑠 ⃗⃗⃗⃗ + 𝑃 ⃗⃗ = 0 ⃗⃗ (1) ∑ 𝑀 ⃗⃗⃗ = 𝑀 𝑓 ⃗⃗⃗⃗⃗⃗ + 𝑀 𝑠 ⃗⃗⃗⃗⃗⃗ + 𝑀 𝑃 ⃗⃗⃗⃗⃗⃗ = 0 ⃗⃗ (2)
where, 𝐹 𝑠 ⃗⃗⃗⃗ is the aerodynamic force from the sail and 𝐹 𝑓 ⃗⃗⃗⃗ is the hydrodynamic force of the foil. The weight of the athlete and the craft 𝑃 ⃗⃗ is also considered. In equation ( 2), moments are computed at the reference point (i.e. bottom of the sail mast). The aerodynamic loads of the board and the athlete are not considered. The system of equations ( 1) and ( 2) is solved using the Newton-Raphson method and the convergence of the VPP is reached when the residual of both forces and moments are below a convergence threshold fixed to 1N and 1N.m.

Numerical force models

In the present section, numerical models to compute aerodynamic and hydrodynamic loads that are used in the VPP are presented.

Aerodynamic loads

In order to model the aerodynamics of the sail, sail shape measurements have been carried out by the ESPCI -Paris. More precisely, the sail has been hoisted on its mast without wind and each sail batten has been measured with 10 points to then reconstruct the whole sail surface. In the present work, this sail shape is used in the VPP and FSI on the sail and mast are neglected. Aerodynamic loads on the sail are determined with a Vortex Lattice Method (VLM). The following hypotheses are considered: inviscid fluid, incompressible and irrotational flow, profile thickness is neglected and turbulence is not modelled. The surface of the sail is discretized into several quadrilateral panels, each one composed of an unknown horseshoe vortex and a control point to prescribe that the normal velocity across the panel is zero. The numerical solution of the lifting surface problem computes each vortex strength then to determine the forces and moments generated by the fluid on the studied thin body. The VortexLattice Julia-based code of Taylor McDonnell and Andrew Ning [START_REF] Day | Performance prediction for sailing dinghies[END_REF] is presently used to compute the sail aerodynamic loads. VortexLattice has been designed and validated [START_REF] Day | Performance prediction for sailing dinghies[END_REF] to obtain fast and consistent results for stationary and unsteady flows.

A VLM mesh convergence analysis has been carried out and has led to discretize the sail with 1250 elements: 50 along the span and 25 along the chord, as shown on Figure 3. The discretization along the span is refined near the top and bottom, using a cosine distribution, to better model the circulation drop.

Hydrodynamic loads

The geometry of the foil (i.e. profile, chord and span) has been built from an average of commercially available windfoil geometries. Hydrodynamic loads on the three lifting surfaces of the foil (shaft, main wing and stabilizer) are computed with the Athena Vortex Lattice (AVL) code [Drela and Youngren 2011], also based on the VLM. A viscous correction is added to account for the viscous drag, which is estimated with the XFOIL code [START_REF] Drela | Xfoil: An analysis and design system for low Reynolds number airfoils[END_REF], based on a 2D panel method where a coupled viscousinviscid interaction approach is employed. The free surface effects can be modelled in a simplified way but at a low computation cost by using a symmetry or an anti-symmetry condition on the free surface plane. It is generally admitted that the symmetry condition is more adapted to low-Froude number situations [Temtching Temou et al. 2021] while the anti-symmetry condition is more suitable for high Froude numbers [Faltinsen 2005, Dutton 2017] (see also [START_REF] Perali | Comparaison de conditions de surface libre linéarisées pour l'étude d'un Hydrofoil submerge à l'aide d'une approche potentielle[END_REF]). In the present work, the foil chord-based Froude number is of order 10 which leads to the anti-symmetry condition as free surface model. This anti-symmetric plane is set at the position corresponding to the flying height and the part of the shaft above the free surface is not considered.

Figure 4 shows the influence of the free surface on the hydrodynamic forces of the foil. Total forces of the foil (main, shaft and stabilizer) are reported and given along X (drag), Y (lateral force) and Z (upward force) components. Forces are scaled by hydrodynamic forces at infinite submergence and are plotted as functions of h/c, where h is the shaft submergence and c is the shaft chord length. Specifically, when h=1m, the shaft is fully submerged whereas when h=0m, the shaft is completely out of the water. Results observed in Figure 4 The free surface effects, as functions of the flying height, on hydrodynamic forces seem to be correctly predicted and this model can thus be used in the VPP to account for effects of the craft flying height on hydrodynamic forces of the foil. For instance, the VPP should be able to model the compromise between flying higher to reduce drag, and keeping enough submersion to achieve the necessary vertical lift and side force.

Fluid Structure Interaction on the shaft

While sailing, the hydrofoil is subjected to a lateral force, which is mainly distributed on the immersed part of the shaft. Thus, the shaft is deformed which, in the end, might influence the lift created and the way of sailing. In this section, a model to determine the deformation of the shaft while sailing is proposed. This model is based on bending tests and the beam theory.

Punctual load

First, a shaft is tested as a cantilevered beam subjected to a point load on its free end. This type of load is easier to apply than a distributed load and is suitable to determine the mechanical characteristics. As the analytical expression is known from beam theory [START_REF] Fanchon | Guide de la mécanique[END_REF]], the identification of stiffness properties is possible (Eq 3).

𝜐 𝑝𝑜𝑛𝑐𝑡𝑢𝑎𝑙 (𝑥) =

-𝐹𝑥² 6 𝐸𝐼 (3𝐿 -𝑥) (3)
where vponctual (x) is the deflection of the beam at the section x, F is the force, L is the total length of the beam, E is the Young's modulus and I is the quadratic moment. This form assumes that E and I are constant along the beam, which in our case is not rigorously exact. However, if the deflection calculated with this model is fitted to the experimental values measured, the average EI product can be identified as can be seen on Figure 5. It appears that the deflection shape is rather well represented with a maximum deviation of 5,11% at 500mm. In this study, the nonlinear effects are neglected and the EI calculated with the experiment is supposed to stay the same whatever the value of the force applied. For this study, a punctual load of 20kg is applied on the shaft at x=900mm. The deflection at this point is measured experimentally and the EI corresponding to the shaft is calculated thanks to equation (3). Twenty shafts were tested and a range of values for EI are found.

Distributed load

While sailing, the load is not punctual but distributed on the immersed part. Here, the load is considered uniform between the water surface and the end of the shaft. Using beam theory, the deflection shape can be calculated (Eq. 4). [ -(𝐿 -𝐴) 4 + 4𝑥(𝐴 3 -𝐿 3 ) + 𝐿 4

-𝐴 4 ] (4)
where :

-A is the first section immersed -L is the total length of the shaft -q is the distributed load, also defined using the first model as q = F/(L-A) -E is the Young's modulus -I is the quadratic moment of the section. The deflection of the whole shaft is then computed using equation ( 4) with the EI previously identified from the tests and considering a distributed load q=F/(L-A) with F the lateral force on the shaft and L-A=Tz, the immerged part of the shaft depending on the height of flight. 4 Results

Deflection of the shaft for different stiffness

The VPP is run to find the operating point of sail at equilibrium for different values of the shaft stiffness.

The stiffness range investigated is larger than what is measured on actual foils, in order to highlight the effects of deformation and test if optimum values may exist even out of the existing range. VPP runs are done for the following cases:

-Perfectly rigid shaft, called "No deformation" -The highest EI value measured, called "stiff" -The lowest EI value measured, called "flexible" -Intermediary EI value equal to 85% of maximum, called 85% -More flexible shafts with EI equal to 75% and 50% of maximum, called 75% and 50%. Figure 7 represents the deflection of the shaft when the VPP has reached the equilibrium in each of these cases. The deflection in mm is shown as a function of the z coordinate along the shaft (1000 mm is the top of the shaft where the board is and 0 mm is the bottom, where the fuselage is). As expected, the more flexible the shaft, the higher the deflection at the end. For a shaft with EI 50% lower than the stiffest shaft tested, the deflection at the bottom of the shaft is almost twice as high. It should be noted that this most flexible case is close to the limit of validity of the linear deformation theory used here. In the following graphs presented in Figure 8, the converged values of the different parameters of the VPP are plotted for the different EI simulated. These values are compared to results obtained for the shaft with no deformation. It can be observed that a difference exists between a non-elastic shaft and the others where the balance point is not the same. In addition, the trend of the evolution of the balance values is consistent with the deformation of the shaft. leeway angle |d| and the heel angle |g| grow as the flexibility increases whereas the height of flight (L-Tz) and the craft velocity decrease. However, in the range of stiffness measured on actual foils, the variation in balance point between the most flexible and the stiffest foils is very low. Figure 9 illustrates the equilibrium point reached for a perfectly rigid shaft and the most flexible investigated case -EI=50%with a shaft half as stiff as the stiffest measured.

In this study, where only a static equilibrium is computed and the stability is not considered, the effect of the stiffness is thus not that important. However, as illustrated in Fig. 9, the stiffness of the shaft does affect the projections of forces and the contribution of the main wing and stab to the lateral force along the Y axis.

Distribution of the lateral force

if the point of equilibrium and the performance are barely affected by the deformation of the foil, it is noticeable that the lateral hydrodynamic force (anti-drift) which is necessarily produced to balance the aerodynamic side force is distributed differently among the different parts of the foil, when the shaft is deformed. Figure 10 shows the distribution of lateral force for the different considered stiffness of the shaft. On this graph (Fig. 10), the total anti-drift force (Py_tot) is plotted as well as the contribution from each part: the main wing, the stabilizer (negative contribution) and the shaft, for different values of shaft stiffness. What can be observed is that as the shaft is more flexible, the lateral force from the main wing is decreasing as the wing is more horizontal, and the contribution of the shaft to the side force increases. This is also correlated to the increase of the leeway angle which results in a higher angle of attack for the shaft. For a perfectly rigid shaft, most of the anti-drift force is given by the main wing thanks to the windward heel, which allows the craft to foil higher, meaning low submerged area of the shaft. For the most flexible shaft studied (50%), the contributions of the shaft and the main wing to the side force are similar, compensated by a lower height of flight.

Conclusions

This paper presents the development of a novel five degree-of-freedom (DoF) VPP for foiling crafts, which includes the FSI effects on the shaft of the vessel. Hydrodynamic forces are computed using a vortex lattice model with a viscous correction, solved by the panel method implemented in the software Athena Vortex Lattice (AVL). The structural behavior of the shaft is introduced in the VPP in the form of an analytical model derived from experimental tests performed on twenty shafts. Shafts of different stiffness are tested, and their respective deflections, as well as the equilibrium state computed by the VPP, are compared. Results show that, even if the deflection varies with the stiffness of the shaft, for a static analysis the equilibrium state and the performance are barely affected. However, the use of a flexible shaft leads to a different distribution of the lateral hydrodynamic force (anti-drift) on the different parts of the foil. The stiffer the shaft, the more the anti-drift force is created by the main wing allowing the craft to fly higher. Further work is underway in two directions. First, in the model for the deformation of the foil, only the lateral force on the shaft was considered and its effect to deflect the shaft. In reality, the horizontal foil can also be deformed and the shaft is submitted to an horizontal compressive force. In addition, during this study, the sailor was kept in a fixed position. It would be interesting to see if, by changing the position of the center of gravity of the sailor, a more elastic shaft could be interesting in terms of performance compared to a stiffer shaft.Last but not least, only the upwind sailing condition is simulated. The downwind situation, with a typical speed around 15m/s (30kts) and heel angle around 5° could also be investigated.

Figure 1 :

 1 Figure 1: Olympic windfoil plane-like hydrofoil and structure displacement under loading (artist's view)

  This work investigates the effect of the deflection of the shaft on the performance of a generic planelike hydrofoil typical of a windfoil using of 5 Degrees of Freedom (DOF) Velocity Prediction Program (VPP: resolving the balance between aero and hydrodynamic forces to find the operating point). The hydrofoils are heavily loaded and thus deformed depending on the stiffness of the shaft and the height of flight. A simplified structural model is proposed thanks to experimental test on shaft to represent the Fluid Structure behaviour of the foil subjected to hydrodynamic loading. Results obtained with a rigid shaft are compared with results obtained with flexible shafts of different stiffness values.

Figure 2 :

 2 Figure 2: Forces on the windfoil. AWA is the Apparent Wind Angle, P is the total weight (sailor+ equipment), Fs is the aerodynamic force from the sail, Px, Py, Pz are the three components of the total hydrodynamic force Ff from the foil package (shaft+main wing+stabilizer).

Figure 3 :

 3 Figure 3: Points measured on the real sail (left) and meshed sail (right)

  Figure4shows the influence of the free surface on the hydrodynamic forces of the foil. Total forces of the foil (main, shaft and stabilizer) are reported and given along X (drag), Y (lateral force) and Z (upward force) components. Forces are scaled by hydrodynamic forces at infinite submergence and are plotted as functions of h/c, where h is the shaft submergence and c is the shaft chord length. Specifically, when h=1m, the shaft is fully submerged whereas when h=0m, the shaft is completely out of the water. Results observed in Figure4are in line with what is expected [Binns et al. 2008, Pernod et al. 2022]. The drag (X) and the lateral (Y) forces decrease as the foil is getting closer to the free surface which can be explained by the reduction of the wetted surface of the shaft. The lift force (Z) firstly increases as the foil gets closer to the free surface because of the acceleration of the flow between the horizontal wing and the free surface. When the foil gets too close to the free surface (h/c below ≈ 1), the lift brutally drops.

Figure 4 :Figure 5 :

 45 Figure 4: h/c versus hydrodynamic forces on the foil scaled by forces at infinite submergence

Figure 1 :

 1 Figure 1: Application of a distributed load between section A and L

Figure 7 :

 7 Figure 7: Shaft deflection at equilibrium for the different considered stiffness of the shaft for the upwind sailing condition

Figure 8 :Figure 9 :

 89 Figure 8: Effect of shaft stiffness on the angles, height of flight and speed of the windfoil

Figure 10 :

 10 Figure 10: Effect of shaft stiffness on the side force Py distribution between the different lifting surfaces of the hydrofoil

Acknowledgments

This work was supported by Region Bretagne and the Agence Nationale de La Recherche (ANR) through grant n°ANR-19-STHP-0002. The authors are grateful to the technical staff at IFREMER and ENVSN for their help in the experiments, ESPCI for the sail shape and the FFV for the contribution in the shaft characterization study.


Résumé :

Les écoulements multi-phasiques naturels ou industriels (glissements de terrain, mélanges de poudres, ...) ont des propriétés rhéologiques non triviales, qu'il est utile de mieux comprendre à partir d'expériences modèles. Les matériaux granulaires humides, à mi-chemin entre les matériaux granulaires secs et les suspensions denses concentrées, sont un exemple de matériaux multi-phasiques, dont la rhéologie reste assez peu documentée. Nous réalisons pour la première fois à notre connaissance, des expériences d'écoulement de matériaux granulaires humides, sur plan incliné : celles-ci révèlent un régime d'écoulement stationnaire uniforme pour une certaine gamme de paramètres expérimentaux (angle de la pente, hauteur de l'ouverture du silo). Nous mettons en évidence la pertinence de la contrainte de cohésion τ c décrivant la limite de stabilité statique via le critère de Mohr-Coulomb avec cohésion τ = τ c + µ 0 P , avec τ et P , les contraintes maximales de cisaillement et de pression et µ 0 le coefficient de frottement interne statique, en écoulement rapide. Ce critère semble ainsi adapté pour décrire la rhéologie des matériaux granulaires humides en écoulement rapide, à condition de caractériser le frottement interne µ = µ 0 en situation d'écoulement rapide et de comprendre ses variations avec les paramètres expérimentaux.

Abstract :

Multi-phase flows encountered in nature or industry (landslides, mudflows, powder mixtures, ...) exhibit non trivial rheological properties, that can be understood better thanks to model materials and appropriate rheometers. Unsaturated wet granular materials, between dry granular and immersed ones (dense granular suspensions), are one example of multi-phase materials, which rheology was not studied much. We realize for the first time to our knowledge, experiments of wet granular flows over a rough inclined plane : our results show steady uniform flows for a wide range of parameters (the inclination angle and the mass flow-rate). We demonstrate the relevance of the cohesion stress τ c well defined in the cohesive Mohr-Coulomd yield criterion τ = τ c + µ 0 P , with τ and P , the maximal values of the shear and pressure stresses respectively, and µ 0 the static internal friction coefficient, for rapid flows. Thus this yield threshold should be extended to the rheology of rapid wet granular flows, by allowing the internal friction µ = µ 0 to differ from its static value µ 0 , after having understood its variations with the experimental parameters.

Mots clefs : Wet granular flow ; Inclined plane ; Rheology ; Cohesion ; Friction 1 Introduction

Wet (unsaturated) granular materials are triphasic systems composed of an assembly of grains which pore space is partially filled by a liquid and air. These model granular systems may be useful for the understanding of multi-phasic materials encountered in natural (landslides, debris flows, rock or snow avalanches) and industrial fields (civil engineering, powder-process such as wet granulation, fertilizer production). They exhibit striking behaviours and unique rheologies, however the understanding of their flow and even their description are still fragmental, due to the many parameters playing important roles and the difficulty of experimental instrumentations. Thus, a unified approach, successfully comparing theoretical and experimental results in different configurations, is still lacking, appealing for much more experimental studies.

It is commonly known that the addition of liquid among grains, unless saturation of the pores, gives cohesion properties to the wet granular material due to surface tension, when liquid form bridges between grains, that attract them, possibly merging in liquid/grain/air clusters. A traditional description at the macroscopic scale is possible by the cohesive Mohr-Coulomb yield criterion, with τ c , the cohesion stress and µ 0 , the static internal friction coefficient :

τ = τ c + µ 0 P, (1)
where τ and P are the shear stress and the normal confining stress (or pressure) respectively. It successfully describes plastic flows of wet grains, in the quasi-static limit of slow flows, with the same static internal friction coefficient µ 0 for wet and dry grains [START_REF] Pierrat | Effect of moisture on the yield locus of granular materials : theory of shift[END_REF][START_REF] Badetti | Rheology and microstructure of unsaturated granular materials : Experiments and simulations[END_REF]. This means that in the quasi-static limit, the internal friction is not modified by the presence of liquid, suggesting the absence of modification of the inter-grain friction [START_REF] Tapia | Influence of surface roughness on the rheology of immersed and dry frictional spheres[END_REF], a question not settled yet [START_REF] Rognon | Rheophysics of cohesive granular materials Europhys[END_REF]. What about the case of non quasi-static but rapid flows [START_REF] Berger | Scaling behaviour of cohesive granular flows[END_REF][START_REF] Mandal | Insights into the rheology of cohesive granular media[END_REF][START_REF] Vo | Additive rheology of complex granular flows[END_REF] ? In a non-uniform shear stress configuration, is this criterion useful at the yield plane only or above the yield threshold in dynamic conditions too ? Can we extend the yield criterion (1) to rapid flows, even above the yield threshold, by introducing the internal friction µ different from its static value µ 0 :

τ = τ c + µ P. (2)
The influence of a liquid wetting the grains on the wet granular material rheology sensitively depends on the liquid content, its spatial distribution around grains and within pores, and the morphology of the liquid phase [START_REF] Mitarai | Wet granular materials[END_REF][START_REF] Scheel | Liquid distribution and cohesion in wet granular assemblies beyond the capillary bridge regime[END_REF]. Due to the high solid fraction of grains and close distances between bonded grains, direct inter-grain contacts play an important role in their rheological behaviour, depending on the shape, surface and grain properties (roughness, friction, elasticity, ...). To those micromechanical features, one should add the physical properties of the liquid (viscosity, surface tension, contact angle, ...) and the geometrical microstructures of the granular assembly [START_REF] Mitarai | Wet granular materials[END_REF][START_REF] Scheel | Liquid distribution and cohesion in wet granular assemblies beyond the capillary bridge regime[END_REF]. Indeed, in addition to the cohesion induced by the presence of liquid interfaces, possible dissipative forces (viscous, capillary, inertial, ...) may occur in the liquid phase. Consequently, the role of some liquid among grains is not obvious at all, thus wet (unsaturated) grains are not as simple as model cohesive materials, with the addition of an inter-grain force interaction [START_REF] Rognon | Rheophysics of cohesive granular materials Europhys[END_REF][START_REF] Rognon | Dense flows of cohesive granular materials[END_REF][START_REF] Vo | Additive rheology of complex granular flows[END_REF].

All this suggests us to realize experiments of wet granular flows over a rough inclined plane. First, the experimental methods and the configuration of the inclined plane are presented in sections 2 and 3 respectively. Second, we demonstrate in section 4 that our measurements of the maximal (shear and normal) stresses in wet granular layers flowing in the steady and uniform regime are compatible with the emergence of a cohesion stress, that is the same than the one appearing in the cohesive Mohr-Coulomb yield criterion (1), suggesting to extend this latter as a rheology (2) for rapid flows. Finally, we conclude in section 5.

Experimental methods

Our granular materials consist of glass beads (density ρ = 2500kg.m -3 , diameter d 250µm ± 50µm) either dry or wet by silicon oil (density ρ l = 950kg.m -3 , viscosity η = 20mPa.s, surface tension Γ = 20mN.m -1 ). Cohesive forces exerted by liquid bridges on bounded grains become important compared to their weight if the Bond number Bo = ρgd 2 /(πΓ) 1 : here, Bo 10 -2 , meaning that grain diameter d is much smaller than the capillary length πΓ/(ρg). The liquid content = m l /m s , defined as the mass ratio of liquid over grains, is either 0 or = 0.5%, so that we focus on the dry state and on the wet pendular state [START_REF] Mitarai | Wet granular materials[END_REF], in which liquid bridges are of small volumes and involve mostly two grains [START_REF] Scheel | Liquid distribution and cohesion in wet granular assemblies beyond the capillary bridge regime[END_REF]. The solid fraction in steady rapid flows is supposed to be φ ≈ 0.5 for wet grains [START_REF] Raux | Cohesion and agglomeration of wet powders[END_REF], while φ ≈ 0.6 for dry grains.

A mass of about 20kg of (wet or dry) grains is stored into the hopper located at the top of the inclined plane, before opening the hopper frontal gate of the height f to release them (figure 1a). The 1.6m-long and 34cm-wide plane is made of a rough bottom plate (sand-paper of roughness 350µm) and two smooth lateral Plexiglass walls, to ensure a no-slip boundary condition at the bottom and negligible sidewall effects [START_REF] Jop | Crucial role of sidewalls in granular surface flows : consequences for the rheology[END_REF] : our experimental parameters lead to shear flows deep inside the granular layer up to the bottom plate, and not only at the surface like in a heap flow. The plane can be inclined at the imposed slope angle θ. The hopper is elevated (∼ 10cm) with respect to the inclined plane, so that when opening its gate, the beads first fall down continuously on the rough plane before flowing over it. This thin and dilute 'rain' allows a good reproducible initial condition and a constant input mass flow rate Q m [START_REF] Freyssingeas | Flowers in flour : avalanches in cohesive granular matter[END_REF].

The mass flow rate Q m is measured from the linear variations with time of the weighted mass m of (wet or dry) grains flowing out of the inclined plane. A laser sheet is projected on the plane with a small angle with respect to the plane normal to the camera axis, so that it is shifted by the presence of a granular layer proportionally to its thickness h(x, t) at the position x and the time t [START_REF] Pouliquen | Scaling laws in granular flows down rough inclined planes[END_REF]. The camera positioned at the normal of the plane, record the central region of the inclined plane at a 20Hz frequency, which images are processed using ImageJ and Matlab.

The inclined plane as a rheometer

It is well established that steady and uniform flows are observed over a rough wide inclined plane for saturated (either dry or immersed) granular materials [START_REF] Pouliquen | Scaling laws in granular flows down rough inclined planes[END_REF][START_REF] Gdr Midi | On dense granular flows[END_REF][START_REF] Cassar | Submarine granular flows down inclined planes[END_REF]. With our experimental set-up and materials, we obviously observe such steady uniform flows with dry grains, from which we will measure in the following (figure 2) the static internal friction µ 0 appearing in equation [START_REF] Pierrat | Effect of moisture on the yield locus of granular materials : theory of shift[END_REF]. The mass flow rates Q m change between 1kg.s -1 and 10kg.s -1 .

For wet (unsaturated) grains, for the first time, to our knowledge, we observe experimentally a steady state and uniform regime too with a finite and (quasi-)constant thickness layer of wet grains propagating along the rough inclined plane, which angle of the free surface is equal to the slope angle θ, for high

a) b) F 1 -a)
The experimental setup consists of a wide rough plane that can be inclined at the imposed slope angle θ. Wet or dry grains, stored in a hopper at the top of the plane, are released through a gate, which aperture thickness f can be adjusted in order to control the mass flow rate Q m , measured by a scale located underneath. The granular layer free surface is recorded by a camera, allowing for the measurement of the height profile h(x, t) thanks to a LASER diode. b) Schematic representation of the vertical section of a wet granular flow : the longitudinal direction is referred to as x and the direction normal to the plane as z.

enough mass flow rates Q m and slope angles θ. The mass flow rates Q m change between 1kg.s -1 and 10kg.s -1 .

We refer to the longitudinal direction as x and the direction normal to the plane as z (figure 1b). The momentum equations for a steady free surface flow down an inclined plane at an angle θ give :

P = φρg(h -z) cos θ, (3)
τ = φρg(h -z) sin θ, (4)
with τ and P , the shear and normal stresses. Thus, in case of a steady uniform flow, the inclined plane geometry, whatever the rheology is, allows to control the apparent friction throughout the whole layer thickness :

µ app = τ P = tan θ. (5)
It is interesting to note that with the extended cohesive Mohr-Coulomb criterion in equation ( 2), the apparent friction coefficient µ app in (5) and the internal friction µ = dτ /dP in (2) differ :

µ app = τ c P + µ. (6)
This may imply some confusions, such as this one : the same internal friction µ for dry and wet grains leads nevertheless to a larger apparent friction as soon as a non-zero cohesion stress τ c appears, leading to larger slope angles for example [START_REF] Bocquet | Moisture-induced ageing in granular media and the kinetics of capillary condensation[END_REF][START_REF] Strauch | Wet granular matter : a truly complex fluid[END_REF].

We recall that it is well known that the friction of dry granular flows (the internal friction being equal to the apparent one) is not exactly constant, and its variations are well described by a function of the inertial number [START_REF] Gdr Midi | On dense granular flows[END_REF] : I = γd/ P/ρ, with γ the shear rate. The common interpretation of the I-dependence F 2 -Parameters space (θ, h) of the experiments reported here, with θ the inclined plane slope angle and h the constant height of the wet (circles) and dry (squares) granular flows. The colors encode for the values of θ. The continuous curve h c (θ) is computed from equation ( 7) with τ c = 70Pa and µ 0 = 0.45, for the height threshold for a cohesive steady uniform free-surface gravitational flow. The dashed line corresponds to the angle threshold θ 0 = arctan(µ 0 ). lies in the comparison of two time-scales at the grain level : the typical time of deformation 1/ γ and the free-fall relaxation time d ρ/P , when the fluid does not apply any drag on the grain. Here the inertial number typically changes between 10 -2 and 1, corresponding to the inertial flow regime, well above the quasi-static regime, for which the friction tends to its quasi-static value.

When some cohesion is introduced in a granular material, the velocity profile of a gravity-induced uniform steady flow is modified [START_REF] Brewster | Plug flow and the breakdown of Bagnold scaling in cohesive granular flows[END_REF][START_REF] Rognon | Dense flows of cohesive granular materials[END_REF][START_REF] Mandal | Insights into the rheology of cohesive granular media[END_REF] compared to the Bagnold profile (with a 3/2 power zdependence). A plug region with no shear appears at the top of the sheared layer, due to the existence of a minimal value for the shear stress to induce some shear flow, as it is for a yield-stress fluid (figure 1b). In general, a yield-stress τ Y (so that the shear rate γ = 0 for a shear stress τ ≥ τ Y ) leads to the emergence of a length-scale h Y , as a height threshold for shear in a steady uniform gravitational flow of density ρ : h Y = τ Y /(ρg sin θ). In a similar way, cohesion in a granular material induces a yield-stress τ Y = τ c + µ 0 P , leading to a height threshold h c for shear in a steady uniform gravitational flow of density φρ :

h c = τ c φρg cos θ 1 (tan θ -µ 0 ) , (7)
with µ 0 , the static internal friction of the granular material, that was shown to be the same with or without liquid (saturated or unsaturated) [START_REF] Pierrat | Effect of moisture on the yield locus of granular materials : theory of shift[END_REF][START_REF] Badetti | Rheology and microstructure of unsaturated granular materials : Experiments and simulations[END_REF]. This means that a steady uniform shear flow should occur only for cohesive granular layers thicker that the height threshold h c ; such a cohesive granular layer (of thickness h ≥ h c ) will experience some shear at the bottom on the height h -h c and a plug velocity profile of height h c above it.

If we wanted to do further predictions on cohesive or wet granular flows, we would need more information about their rheology and thus the measurement of the shear rates within our flows, that is not obvious as soon as shear is localized, which will be devoted to a future work. Instead, here, we will investigate if our experimental measurements are compatible with a cohesive Mohr-Coulomb-type rheology, with some cohesion stress relevant both for the yield criterion and for rapid flows.

The cohesion stress

The cohesion stress at the yield threshold

Figure 2 shows the thickness h of steady-uniform flows for the different slope angles θ explored for both the dry (squares) and the wet (circles) granular samples ( = 0.5%). Steady-uniform flows are achievable for large enough values of h (or f ) and θ. One can see that wet granular flows are much thicker than dry ones. Indeed minima of h (and θ) for wet granular flows are much larger than the ones for dry granular flows.

As we did not observe any steady uniform dry granular flows for slopes smaller than 24.6 ± 0.5 • , we identify the angle threshold θ 0 = 24 • and get for the static internal friction µ 0 = 0.45, represented by the dashed line in figure 2. The relation [START_REF] Vo | Additive rheology of complex granular flows[END_REF] for the height threshold is drawn as a continuous line for τ c = 70Pa and µ 0 = 0.45 : h c (θ) is below the values (θ, h) of the steady uniform wet granular flows reported here, seeming to describe well the yield threshold. As the minimal flowing height h c (θ) increases and diverges when the slope approaches the angle threshold θ 0 , getting steady uniform wet flows at smaller θ would need much larger thicknesses h, gate heights f and mass flow rate Q m , not achieved in our set-up. That is why minima of θ for our wet granular flows seem larger than the ones for our dry granular flows. Now we will check whether or not this cohesion stress τ c = 70Pa describing the yield threshold for our steady uniform wet granular flows ( = 0.5%) over an inclined plane is useful to describe our rapid flows.

The cohesion stress relevant for rapid wet granular flows

In our configuration of a gravity-induced uniform steady flow, under the assumption of a no-slip boundary condition at the bottom, the material is sheared from the bottom (z = 0) to the interface with the plug (z = h -h c ), we choose to plot the maximal values of the pressure and the shear stress (max(P ), max(τ )) as coordinates ; each point (both for dry and wet flows) is obviously computed from equations (3) and ( 4) for z = 0 : (P (0), τ (0)). Thus, each set of data realized at a constant slope angle θ is aligned according a straight line of slope tan θ passing through the origin, as shown in figure 3a for θ = 35 • , by wet (circles) and dry (squares) flows. As expected from figure 2 showing steady uniform dry flowing layers thinner than wet ones, (P (0), τ (0)) data of the dry glass beads have smaller values than the wet ones : dry data are closer to the origin than wet ones, these latter must have the shear stress larger than the yield stress (see Section 4.1).

However, when we plot all of our experimental wet data, for slope angles θ ∈ [29-37] • , as shown by the circles in figure 3b and4, surprisingly, they all lie along a single line, demonstrating that the wet material can be straightforwardly described by an affine law, indicating a finite cohesion stress τ c and an approximately constant internal friction µ 0.53. In this figure 3b, all of our dry granular flows, for slope angles θ ∈ [24-39] • , are also plotted (squares), but they have not the same friction (the internal friction being equal to the apparent one) µ = µ app ∈ [0.45-0.81], as shown by the non parallel lines delimited by the colored lines and squares : this is a consequence of the dependence µ(I) with the inertial number. That's why we restrict data from our dry granular flows to a small range of slopes θ ∈ [26-30] • in order to get approximately the same internal friction µ 0.53 in figure 4.

Finally, our major result is visible in figure 4, where are plotted the maximal values of shear τ (0) and pressure P (0) stresses for our rapid granular flows, either dry or wet, but having approximately the same internal friction, so that the two sets of data lie along two parallel lines. Whereas dry data pass through the origin indicating no cohesion (τ c = 0), all of our wet data while steadily and uniformly flowing over the inclined plane lie along a single line, whose shift with respect to the origin demonstrate the emergence of a cohesion stress τ c , even above the yield threshold. The measurement of this cohesion stress τ c 70Pa is compatible with the one characterizing the yield threshold (see Section 4.1). Besides, this shift between dry and wet data (at the same internal friction) indicate that the apparent friction coefficient of the wet beads is higher than the dry ones, only because of a finite cohesion stress.

Conclusion

For the first time to our knowledge, we realized experiments of wet (unsaturated) granular flows over a wide rough inclined plane and we observe experimentally a steady state and uniform regime with a finite and (quasi-)constant thickness propagating at a constant velocity for high enough mass flow rates Q m and slope angles θ. As expected from other configurations [START_REF] Pierrat | Effect of moisture on the yield locus of granular materials : theory of shift[END_REF][START_REF] Badetti | Rheology and microstructure of unsaturated granular materials : Experiments and simulations[END_REF][START_REF] Abramian | The slumping of a cohesive granular column : Continuum and discrete modeling[END_REF], the minimal thickness h c for a wet granular layer to be sheared uniformly and steadily is well described thanks to the cohesive Mohr-Coulomb yield criterion, in equation ( 1), with τ c , the cohesion stress and µ 0 , the static internal friction coefficient, being the same for wet and dry grains. Surprisingly, extended equation (1) by allowing some variations of the internal friction to equation (2) allows to describe a whole set of data of wet granular flows (at a constant liquid content and for a given liquid), even if far above the yield threshold, and without needing to be at the yield plane [START_REF] Mandal | Insights into the rheology of cohesive granular media[END_REF]. This demonstrates the relevance of the cohesion stress τ c in equation ( 2) for rapid wet granular flows. Furthermore, the same value of the cohesion stress τ c holds for these two regimes.

In the future, we need to characterize the internal friction of wet granular flows and to understand its variations, especially when changing the liquid content and the properties of the liquid. To this aim, we will have to be able to measure the shear rate characteristic of such localized flows.

  Wet granular flows over a rough inclined plane S. DEBOEUF a , C. NOÛS b , A. FALL c a. Sorbonne Université, CNRS, UMR 7190, Institut Jean Le Rond d'Alembert, 75005 Paris + sdeboeuf@dalembert.upmc.fr b. Laboratoire Cogitamus c. Labo. Navier, UMR 8205, CNRS, Ecole des Ponts ParisTech, Univ. Gustav Eiffel, 77420 Champs sur Marne

F 3 -

 3 Shear and normal stresses τ (0) and P (0) computed at the the bottom (z = 0) where they are maximal : a) for our dry and wet flows at θ = 35 • ; b) for all of our wet flows at θ ∈ [29-37] • and for all of our dry flows at θ ∈ [24-39] • , having neither the same apparent friction µ app nor the same internal friction µ. The colors encode for the inclined plane slope angle θ as in figure 2.

F 4 -

 4 Shear and normal stresses τ (0) and P (0) computed at the the bottom (z = 0) where they are maximal, for all of our wet flows and for dry flows at θ ∈ [26-30] • , equivalent to µ ∈ [0.5-0.6], such that they have approximately both the same internal friction µ. The colors encode for the inclined plane slope angle θ as in figure 2.


Jean-Pierre Dalmont

Rémi Mattéoli

Sylvain Maugeais

Soizic Terrien

Christophe Vergez

THE BICHROMATOR AS THE QUINTESSENCE OF A SAXOPHONE AND THE QUESTION OF HARMONICITY

Keywords:

Reed wind instruments are divided into two families: those with a cylindrical bore (clarinet) and those with a conical bore (saxophone, oboe, bassoon). Apart from the timbre, the main features of the behavior of a clarinet can be obtained with a resonator modeled with a single mode (monochromator), while that of a saxophone can be obtained with a resonator modeled with two modes (bichromator), whose proper frequencies ratio is close to 2. This system has been largely studied by J. Gilbert et al. In particular, it has been demonstrated that when the inharmonicity is small, the oscillation on the fundamental frequency is obtained by an inverse bifurcation. Moreover, perfect harmonicity represents an optimum in terms of threshold pressure. Now a question remains: to what extent does the inharmonicity impact playability? This question is investigated in the present paper using both experiments and advanced numerical modeling.

INTRODUCTION

In the 1980s, Jean Kergomard, with Pierre Boulez and Ircam, worked on a quarter-tone wind instrument project.

---------* Corresponding author: jean-pierre-dalmont@univlemans.fr Copyright: ©2023 Jean-Pierre Dalmont et al. This is an openaccess article distributed under the terms of the Creative Commons Attribution 3.0 Unported License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

The principle was to branch on the instrument, near the mouthpiece, a small narrow tube whose effect was to shift the range by a quarter tone, thus allowing easy access to micro intervals without using complicated fingerings. This system worked rather well for the flute and the clarinet and had resulted in a patent. Unfortunately, the device did not experience the expected development. Nevertheless, Jean Kergomard had commissioned us (i.e. Joël Gilbert and Jean-Pierre Dalmont) to design a similar device for the saxophone and the bassoon. The problem turned out to be much more difficult than for cylindrical instruments because it quickly became clear that shifting the first impedance peak would not be enough because that would increase the inharmonicity to a point that makes the instrument unplayable. This produced, in certain situations, pseudo-periodic sounds combining the first two eigenfrequencies. This work allowed us to wonder about the functioning of reed instruments with conical bores. Examination of the impedance curves shows that conical instruments have resonance frequencies that are almost harmonic, following the series f 1 , 2f 1 , 3f 1 etc. with a first impedance peak of amplitude often lower than the second. On the other hand, the clarinet works on a series f 1 , 3f 1 , etc. and the first peak has a larger amplitude than all the others. We have come to the conclusion that, just as many aspects of clarinet operation can be explained with a single-mode resonator model (the monochromator), many aspects of saxophone operation can be explained with a two-mode resonator model: the so called bichromator. Investigating the dynamics of this model ( [START_REF] Grand | Oscillation threshold of woodwind instruments[END_REF], then [START_REF] Dalmont | Reed instruments, from small to large amplitude periodic oscillations and the Helmholtz motion analogy[END_REF] and [START_REF] Gilbert | Minimal blowing pressure allowing periodic oscillations in a simplified reed musical instrument model: Bouasse-Benade prescription assessed through numerical continuation[END_REF]) has made it possible to understand part of the specificities of the saxophone and more generally of reed instruments with conical bore.

BICHROMATOR WITH PERFECTLY HARMONIC RESONANCES

Analysis of the impedance curves of reed instruments with conical resonators shows that the amplitude of the first impedance peak is higher than all the others for the notes at the top of the first register (beyond the fingering of A on a saxophone or an oboe) while it is lower than the second for the fingerings at the bottom of the first register. It was therefore relevant to focus on the importance of the relative amplitude of these two peaks on the oscillations.

In [START_REF] Grand | Oscillation threshold of woodwind instruments[END_REF], appears a diagram of first importance (Fig. 1). This diagram is a plot of the following equation which gives the condition for a direct bifurcation [START_REF] Grand | Oscillation threshold of woodwind instruments[END_REF]:

< ² ( ) , ( 1
)
where and are the second and third coefficient of the Taylor expansion of the non-linear characteristics ( ) with the volume velocity and the pressure in the mouthpiece. = and = are the admittances of the two modes. This shows that when >> and < 0 (which is always the case for the saxophone), the Hopf bifurcation responsible for the emergence of oscillations is direct (top corner left on Fig. 1) as in the case of the clarinet. This means that when the pressure is gradually increased, the reed destabilizes for a given pressure called "direct pressure threshold" and begins to oscillate first with a low and then growing amplitude (see Fig. 2a). The diagram in Fig. 1 shows that when is only slightly larger than , the bifurcation becomes inverse, i.e. when the pressure is gradually increased, the reed destabilizes at the "direct pressure threshold" but begins to oscillate with a large amplitude (blue arrow on Fig. 2b). It is then possible to reduce the pressure until the sound dies out at a pressure lower than the direct threshold pressure, known as the inverse pressure threshold or subcritical threshold (red arrow on Fig. 2b), The analysis in [START_REF] Grand | Oscillation threshold of woodwind instruments[END_REF] shows that if > the bifurcation is always direct because < 0 (Eqn. ( 1)). This result is surprising to say the least: indeed, this situation corresponds to the notes of the bottom of the first register and there is no reason for the bifurcation to become direct in a sudden manner when = . The question was investigated in [START_REF] Dalmont | Reed instruments, from small to large amplitude periodic oscillations and the Helmholtz motion analogy[END_REF], where it has been shown that the direct bifurcation actually corresponds to the branch of the inverted Helmholtz motion which is only observed at high level, is the admittance difference between the second and first modes (from [START_REF] Dalmont | Reed instruments, from small to large amplitude periodic oscillations and the Helmholtz motion analogy[END_REF]). HM stands for Helmholtz Motion. beyond the extinction threshold of the Helmholtz motion [START_REF] Ollivier | Idealized models of reed woodwinds. Part I: analogy with the bowed string[END_REF][START_REF] Ollivier | Idealized models of reed woodwinds. Part II: on the stability of 'two step' oscillations[END_REF][START_REF] Colinot | Numerical Simulation of Woodwind Dynamics: Investigating Nonlinear Sound Production Behavior in Saxophone-like Instruments[END_REF]. The branch of the Helmholtz motion is connected to the branch of the second regime (at the octave) whose direct threshold is

. In this situation is lower than the direct threshold of the fundamental (Fig. 3): indeed the higher the impedance, the lower the threshold. The Helmholtz motion thus emerges from an inverse bifurcation with a brief passage over the octave. These results have an obvious practical consequence: it is difficult or even impossible to play pianissimo in the low register with a saxophone.

Figure 3. Bifurcation diagram for

> . RMS pressure with respect to the mouthpressure . All quantities are normalized by the minimum closing pressure . +++ unstable solutions (from [START_REF] Dalmont | Reed instruments, from small to large amplitude periodic oscillations and the Helmholtz motion analogy[END_REF]).

SLIGHTLY INHARMONIC BICHROMATOR

In practice, the harmonicity of the resonance frequencies is never perfectly satisfied, especially at the top of the first register, because the length of the truncation is no longer small compared to the wavelength. The truncated cone with equivalent volume is then no longer a good approximation of the stepped cone. The problem of the slightly inharmonic bichromator has been treated in [START_REF] Gilbert | Minimal blowing pressure allowing periodic oscillations in a simplified reed musical instrument model: Bouasse-Benade prescription assessed through numerical continuation[END_REF] where it is also observed that the lowest oscillation thresholds are obtained for perfectly harmonic resonances thus confirming the prescription of Bouasse-Benade [START_REF] Bouasse | Instruments à vent tomes I et II[END_REF][START_REF] Benade | Sound production in wind instruments[END_REF]. Moreover, it appears that the almost non-existent octave stability range in the harmonic case takes on a significant proportion. It is observed that there is a pressure range in which the two regimes are stable (Fig. 4). In practice, this means a higher risk of inadvertent emission of the octave. On the other hand, this phenomenon can be advantageously used to make it possible to play the note and its octave with the same fingering as is the case for the baroque oboe for example. Indeed, on this instrument, there is no register key and the octave is obtained by a subtle modification of the embouchure.

QUANTITATIVE ANALYSIS OF THE INFLUENCE OF THE INHARMONICITY AND THE AMPLITUDE RATIO OF THE TWO MODES

Beyond the qualitative observations in [START_REF] Gilbert | Minimal blowing pressure allowing periodic oscillations in a simplified reed musical instrument model: Bouasse-Benade prescription assessed through numerical continuation[END_REF], we would like to be able to easily deduce the possible playing frequency from the impedance curve. In the present work, with the methods used in [START_REF] Gilbert | Minimal blowing pressure allowing periodic oscillations in a simplified reed musical instrument model: Bouasse-Benade prescription assessed through numerical continuation[END_REF], we aim at evaluating the playing frequency as a function of both the inharmonicity and the amplitude ratio of the two modes. The results will be compared with those obtained with the "sum function" of Wogram and Worman [START_REF] Wogram | Ein Beitrag zur Ermittlung der Stimmung von Blechblasinstrumenten[END_REF][START_REF] Worman | Self-sustained non-linear oscillations of medium amplitude in clarinet-like systems[END_REF]. We can also compare with Grothe's proposal who weights the reduced frequency / of each mode by their amplitudes. We also investigate to what extent studying the dependence of stability ranges of the sound regimes on the harmonicity and relative amplitude of modes can help predict playing difficulties.

EXPERIMENTAL VALIDATION

In [START_REF] Mattéoli | Étude de la facilité d'émission des cuivres et des instruments à anche au moyen des diagrammes de bifurcation[END_REF], a first experimental validation was attempted with an artificial mouth and a saxophone body extended with cylindrical tubes of different lengths. It turns out that the hypothesis of a linear resonator is not tenable with a saxophone because the non-linear losses at the side holes [START_REF] Atig | Saturation mechanism in clarinet-like instruments, the effect of the localized non-linear losses[END_REF][START_REF] Dalmont | Oscillation and extinction thresholds of the clarinet: Comparison of analytical results and experiments[END_REF] are in practice never negligible. Indeed the comparison between theory and experiment on Fig. 5 emphasize the fact that the dynamics is more reduced in practice than in the theory. This work will be resumed with resonators made up of two cylinders, the last cylinder being extended by a network of large lateral holes. Whatever the results of these experiments, the existence of non-linear losses in the instruments with lateral holes must lead us to be cautious about the results of the numerical calculations and in particular for instruments with small lateral holes. Light/dark green: experiment with/without reed dynamics. Blue/orange marks: experiments with increasing/decreasing mouthpressure (from [START_REF] Mattéoli | Étude de la facilité d'émission des cuivres et des instruments à anche au moyen des diagrammes de bifurcation[END_REF]).

CONCLUSION

The approach proposed by Joël Gilbert, consisting in approaching the saxophone by a two-mode model, has proven to be extremely fruitful for the understanding of its operation. The model has not been fully exploited yet, and with current digital tools can be used to unveil behavior laws to extrapolate playing frequencies and anticipate some emission difficulties.

Figure 1 .

 1 Figure 1. Diagram showing the regions where the bifurcation is direct as well as those regions where it is inverse (Eqn. (1)). Abscissa is for coefficient C, cubic term of the Taylor expansion of the non-linear characteristics ( ).is the admittance difference between the second and first modes (from[START_REF] Dalmont | Reed instruments, from small to large amplitude periodic oscillations and the Helmholtz motion analogy[END_REF]). HM stands for Helmholtz Motion.

Figure 2 .

 2 Figure 2. Diagram showing the amplitude | | of the first harmonic as a function of the mouth pressure for two situations: a) direct bifurcation b) inverse bifurcation (dotted line: unstable solution) (from [1]).

Figure 4 .

 4 Figure 4. Bifurcation diagram for = and ℎ = -1 = 0.03. Thick/thin lines: stable/unstable oscillations. Blue: fundamental Helmholtz motion. Red: octave. Green: inverted fundamental Helmholtz motion (from [3]).

Figure 5 .

 5 Figure 5. Experimental and theoretical bifurcation diagram for fingering C of an alto saxophone.Light/dark green: experiment with/without reed dynamics. Blue/orange marks: experiments with increasing/decreasing mouthpressure (from[START_REF] Mattéoli | Étude de la facilité d'émission des cuivres et des instruments à anche au moyen des diagrammes de bifurcation[END_REF]).