Zenjebil Jouini

Zarel Valdez-Nava

David Malec

Failure Analysis of Ceramic Substrates Used in High Power IGBT Modules

Keywords: High Voltage Power Module, Ceramic Substrate, Dielectric Strength, Electro-Mechanical Breakdown

High voltage power modules are used in numerous applications to build high power converters. Technically, these modules are made of different materials and among them, dielectric materials are organic and inorganic. Organic insulators (gels) are used to avoid corona discharges in the vicinity of connecting wires and high voltage dies (diodes and transistors) and to protect them from moisture and contaminants. Inorganic insulators (ceramic substrates) are used to insulate the high voltage which dies from the grounded elements and to transfer heat to the heat sink. Despite being used since the late 90s, there is a lack of fundamental knowledge about the electrical properties of these substrates. Consequently, manufacturers tend to assure the reliability by over sizing them. As there are no clear rules for how to do that, failures occur, leading to the converter shutdown. The aim of this study is to bring new information about the understanding of the dielectric strength of ceramic materials used in these modules. We have focused our work on the correlation between the mechanical and the dielectric properties of ceramics by using relevant experiments. We provide new information about the impact of existing cracks on the ceramic dielectric failure, according to the electromechanical breakdown model. Our conclusions bring crucial information about the precautions to be taken during manufacturing and implementation of these substrates in power modules to reduce the likelihood of the particular causes of failure.

Introduction 1.Modules and Ceramics

Currently used high voltage power electronics systems are based on an elementary brick called "power module", where the active semiconductors are packaged. Inside this assembly, one finds different materials, whose functions are to assure the interconnection, the insulation and the heat exchange between the semiconductor and its environment. A simplified schema of a power module is shown in Figure 1.

The interconnection of the semiconductors and their insulation from the heat sink (connected to the ground) is mainly achieved by a single element called "substrate". Furthermore, one of its main functions is to warrant the electrical insulation between the semiconductors and their environment. Moreover, it must also ensure the evacuation of the heat generated by the semiconductors. Theses substrates are composed of metallic conductors and an insulating material. When the thermal requirements are high, ceramics are the first choice due to their high thermal conductivity (>20 W/m•K).

Different ceramic materials are used for substrates elaboration (e.g. BeO, Si 3 N 4 , AlN) [START_REF] Chasserio | Ceramic Substrates for High-Temperature Electronic Integration[END_REF]. Different technologies allow joining the ceramic to a conducting metal, such as DBC (Direct Copper Bonding) or AMB (Active Metal Brazing). Reliability of the substrate is often related to the mechanical properties of the ceramic material coupled with the metal being used. With the development of the new wide bandgap semiconductors, such as SiC or GaN, higher voltages along with higher current densities could potentially increase the insulating requirements imposed on the ceramic materials inside the substrates [START_REF] Diaham | Dielectrics for High Temperature SiC Device Insulation: Review of New Polymeric and Ceramic Materials[END_REF].

Different causes may lead to the module failure:

-dies failures [3] [4],
bonding wires cut [START_REF] Ciappa | Selected Failure Mechanisms of Modern Power Modules[END_REF]- [START_REF] Ramminger | Crack Mechanism in Wire Bonding Joints[END_REF],

solder joints delamination [START_REF] Woirgard | Contribution à l'étude d'une intégration hybride adaptée à l'électronique automobile[END_REF]- [START_REF] Rizzi | Contribution à l'étude de la fiabilité des modules de puissance pour application automobile[END_REF],

encapsulation gel dielectric breakdown [START_REF] Fabian | Analysis of Insulation Failure Modes in High Power IGBT Modules[END_REF]- [START_REF] Lebey | Contribution to the Dimensioning of The Different Insulating Materials Used in High Voltage Power Electronics Module Manufacturing[END_REF],

substrate dielectric breakdown [START_REF] Malec | Aging of Ceramic Materials in Power Electronic Substrates[END_REF] [START_REF] Decup | Impact of a Surface Laser Treatment on the Dielectric Strength of α-Alumina[END_REF].

The design and the choice of the substrate depend on several parameters, such as the mechanical properties of the metal-ceramic assembly, the thermal properties of the ceramic and its electrical insulating capability. As a basic principle, if the ceramic material is capable of reaching very high operating voltages, its thickness can be reduced, hence decreasing the thermal resistance along the heat evacuation path. The dielectric Figure 1. Schematic drawing of a power module used to package semiconductor dies.

strength is then a critical dimensioning parameter for the ceramic materials used in these power modules. Despite this, there is little available information concerning the relationship among the processing parameters, the ceramic microstructure and the dielectric breakdown field.

For ceramics, as for the other solid materials, the origin of the dielectric breakdown will highly depend on the thermodynamic conditions at which the dielectric breakdown occurs. For example, alumina and other ceramics used at very high temperatures (>600˚C), will show a thermal breakdown due to its high electrical conductivity [START_REF] Yoshimura | Electrical Breakdown Strength of Alumina at High Temperatures[END_REF].

Nevertheless, if we consider much lower temperatures of currently used modules (i.e., 135˚C for Si dies), the substrates are confronted to lower operating temperatures. At these temperatures, thermal breakdown is already excluded as a possible mechanism.

Previous works pointed out that the electric breakdown of ceramics at room temperature could be strongly linked to their mechanical properties [START_REF] Carabajar | Relations between Electrical Breakdown Field and Mechanical Properties of Ceramics[END_REF]. It is assumed that the electric field induces the propagation of the pre-existing cracks within the ceramic material. This assumption has been studied under a high electric field (short term breakdown) [START_REF] Fothergill | Filamentary Electromechanical Breakdown[END_REF] or under a moderate field (aging studies) [START_REF] Zeller | Electrofracture Mechanics of Dielectric Aging[END_REF].

This hypothesis is described in the literature as an "electromechanical dielectric breakdown", based on the analogy of the theory with mechanical breakdown in brittle materials. In such a model, the pre-existing cracks or defects are propagated under the influence of the applied electric field until the final breakdown (Figure 2).

One of the main challenges to verify this hypothesis is to control the size of the initial defect or crack in the ceramic. Previous attempts included the use of ceramics sintered at different temperatures, with different microstructures or with different compositions [17] [22] [23]. This could induce a change in the intrinsic dielectric breakdown strength of the studied material, potentially biasing the results.

In this work, we will employ a mechanical pre-stress in order to produce different defect sizes inside the ceramic, without modifying the chemical composition of its microstructure. Alumina (Al 2 O 3 ) substrates, largely used in power modules, have been selected as the representative samples. Results will be confronted with the electromechanical breakdown model. Finally, our experiments will be analyzed to extract the recommendations on the design rules that could be taken into account in the ceramic substrates manufacturing process.

Origin of the Dielectric Breakdown

The purpose of this work is to understand the mechanisms of the electric breakdown of the ceramics used in the elaboration of substrates for power electronics. The dielectric breakdown of solids can be classified into three types, depending on the underlying physical phenomena: 1) of electronic origin (electronic avalanche or field or intrinsic effect), 2) of thermal origin, which dissipates the heat caused by the current of conduction or dielectric losses and 3) electromechanical origin explaining the inability of the material to withstand mechanical stresses induced by the electric field.

Generally speaking, these mechanisms are not independent and may interact until the final breakdown. It is therefore often difficult to determine and isolate one physical phenomenon responsible for the dielectric breakdown.

Experimental Procedure

Ceramic Material under Study

The chosen ceramic materials are industrially available alumina (Al 2 O 3 ) plates (thickness 635 µm), typically used DBC ceramic substrates for power modules. The purity of the alumina is reported as 96%; with 4% of other sintering additives and impurities (mainly Si, B, Mg, Ca). The alumina has an average grain size of 4 µm observed by scanning electron microscopy (Figure 3).

In order to limit the dispersion due to the processing conditions, we paid special attention to make sure that the alumina plates, originated from the same manufacturer, were extracted from the same batch number. Random sampling (18 tests) showed that the dielectric breakdown of the analyzed batch was 21 kV/mm (at room temperature).

Breakdown experiments are detailed in section 2.3. Alumina plates were then cut out from the as-received dimension of 115 × 115 mm to smaller plates of 38 × 38 mm by laser cutting. After laser cutting, the dielectric breakdown field in the central part of the alumina remained unchanged.

Mechanical Prestress

In the electromechanical breakdown model, the dielectric breakdown strength is limited due to the presence of defects, following the analogy with the Griffith's proposal Figure 3. Microstructure of the alumina under study.

for mechanical breakdown. In its most basic form, these defects could be cracks, crevices or any mechanical discontinuities inside the bulk of the material. As stated previously, in order to avoid any chemical or microstructural modification, we applied a mechanical pre-stress to our alumina samples. Mechanical stresses can be applied by the compression, traction or flexion. Considering the high compressive strength of alumina, we decided to apply a flexion force. One of the simplest configurations to apply a flexural stress is by a 3-point bending test (Figure 4a). We assume that in this configuration cracks will be created at the maximum tensile stress region of the 3-point bending configuration, i.e. in the center region where the upper ram applies the force (Figure 4b). In these conditions, initial cracks will be propagated or new ones will be formed. The large difference in the compressive (2000 GPa) and tensile (200 GPa) strength of alumina will favor the formation of the cracks on the tensile region of the cross-section of the ceramic material (opposite to the contact point of the upper ram).

Total applied mechanical stress is a function of the sample dimensions and test bench geometry as follows: where M is the moment bending, l is the length of the support (outer) span, F is the load (force) at the fracture point, B is the base of the sample, h is the effective thickness of the sample.

Since the formation of a crack or a network of cracks reduces the effective section of the alumina plates, a constant force is applied, the real stress will increase as the crack(s) propagate(s) along the thickness of the material. The final crack dimension will be then related to its initial size, the applied force and the time spent under such a force. A series of tests were performed to determine the mechanical prestress conditions (force and application time) prior to electrical breakdown. The characteristic of the time-to-breakdown (mechanical) vs. the applied force is shown in Figure 5.

The maximum force possible to apply is 180 N, which corresponds to the mechanical flexural strength of the studied alumina. An attempt was made to study the full range of applicable forces. However, samples needed to be transported from the mechanical test bench to the dielectric breakdown test bench, and for high values of the force (>100 N) even a short time of prestress could break the samples during their manipulation. For this study, samples were prestressed from 20 N to 80 N, applying the force for 1 minute.

This allowed the samples to be handled without the risk of breaking them. Samples were marked with a felt-tip pen to identify the region where the force ram was in contact.

Dielectric Breakdown of Prestressed Alumina

The dielectric strength measurements were performed in a plane-tip configuration. The A scheme of the followed procedure is shown in Figure 6.

Results and Discussion

Figure 7 shows the ceramic plates after the dielectric breakdown was performed. The size of alumina plates allowed for 3 to 4 dielectric breakdown tests per sample.

When plotting the results in a Weibull chart we can easily identify that the dispersion of the results obtained for the prestressed samples are significantly higher than the ones obtained for the un-stressed sample (0 N) (Figure 8).

To better observe the impact of the mechanical prestress, only the α parameter (cumulative probability of 63.2%) of the dielectric breakdown extracted from a two-coefficient Weibull distribution was plotted (Figure 9). We should consider that for application purposes, the alumina (or any other ceramic mentioned in section 1) is a part of a complex assembly with other materials, with a metal being in immediate ontact with the substrate assembly. While there are some studies that highlight the importance of the initial mechanical properties of the ceramic within the substrate assembly [START_REF] Dupont | Effects of Metallization Thickness of Ceramic Substrates on the Reliability of Power Assemblies under High Temperature Cycling[END_REF], the data presented here indicates that the other processing parameters of the alumina could be relevant. This could indicate the need to take into account the mechanical history of the ceramic as having an impact on the dielectric strength value.

Conclusion

The dielectric breakdown of the studied alumina depends on the value of the mechanical prestress created by a 3-point bending test. Despite the large dispersion of the results, a global tendency satisfies the currently used model of electromechanical breakdown, where the defect size impacts the short-term dielectric strength of alumina.

Higher values of mechanical prestress induce a decrease in the E br values. This could imply that not only the original dielectric properties of the ceramics should be taken into account for dimensioning purposes, but also their full mechanical history.

Figure 2 .

 2 Figure 2. General model for an electromechanical breakdown (inspired by Griffith's proposal).

Figure 4 .

 4 Figure 4. (a) 3-point bending test bench used for applying a mechanical prestress on alumina plates, (b) tensile stress is created at the bottom part of the alumina plate.

Figure 5 .

 5 Figure 5. Time-to-breakdown (mechanical) vs. applied mechanical force for 635 µm alumina plates under study.

Figure 6 .

 6 Figure 6. (a) Reference sample with an initial crack length l 0 , (b) mechanical prestress that increases the crack length inside the ceramic, total length is l 0 + l m , (c) dielectric breakdown of the prestressed ceramic.

Figure 7 .

 7 Figure 7. Dielectrically broken down alumina sample. Top surface was submitted to compression (black line, felt-tip pen) and bottom suface under tensile stress. Red circles indicate the rupture channel exit points on the surface of the alumina.

Figure 8 .

 8 Figure 8. Full results for the dielectric breakdown of prestressed alumina plates (0 to 80 N, 1 minute).

Figure 9 .

 9 Figure 9. α values for the dielectric breakdown of prestressed alumina plates (0 to 80 N, for t appl = 1 minute), error bars correspond to the 90% confidence intervals.

Acknowledgements

The authors acknowledge the assistance from E. Duhayon for providing access to the mechanical test bench. Z. Jouini thanks the French Ministry of Research and Higher Education for its financial support through a doctoral studies grant.


Mario Spirito

Niels Boulevard

Bohr

Modeling of Isometric Muscle Properties via Controllable Nonlinear Spring, and Hybrid Model of Proprioceptive Receptors

Keywords: Skeletal muscles isometric force modeling, Nonlinear variable spring, Hybrid dynamical modeling, Spindle-Ia & -II, Oculomotor system

This work presents the macroscopic behavior of skeletal muscles seen from a system theoretic standpoint. Then, according to the data available in the literature, we propose a first Evaluation model for the muscle isometric (constant muscle length) force generation depending on the muscle length and neuronal excitation frequency. Such a model is instrumental in providing a more physics-inspired model of such an isometric force by exploiting a nonlinear spring description with controllable characteristics, such as stiffness and rest length. To conclude, we propose a hybrid dynamical filter model to describe the components of the sensory system. This is in charge of feeding back to the Central Nervous System the measurements of muscle length and its rate of change.

Introduction

Skeletal muscles, or simply muscles, are organs of the muscular system whose cells have the characteristic to produce force and movements when electrically stimulated. They are mostly attached to the bones of the skeleton by tendons, and thanks to their contractile and elastic properties they allow the relative motion of the bones of the involved articulations.

Each skeletal muscle consists of thousands of muscle fibers wrapped together by connective tissue sheaths, see any human anatomy book, e.g. [1][Ch.9]. The individual bundles of muscle fibers in a skeletal muscle are known as fasciculi. Each muscle fiber is comprised of several myofibrils containing many myofilaments. When bundled together, all the myofibrils get arranged in a unique striated pattern forming sarcomeres which are the fundamental contractile unit of a skeletal muscle. The two most significant myofilaments are actin and myosin filaments which are arranged distinctively to form various bands on the skeletal muscle [START_REF] Lieber | Skeletal muscle structure, function, and plasticity[END_REF][Ch.1], [START_REF] Frontera | Skeletal muscle: a brief review of structure and function[END_REF] and [START_REF] Sweeney | Muscle contraction[END_REF]. Skeletal muscles show up with neuronal innervations from the sensory-motor system. A set of neurons, the α-motoneurons, constitutes the neuromuscular junction. They are responsible, via their excitation frequency here labeled as α (measured in pulses per second), for the sarcomeres contractions, and hence they allow the generation of the muscle active force. On the other hand, another set of neurons, constituting the proprioceptive receptors, carries sensorial information from the muscle to the Central Nervous System. In particular, the transferred signals are the values of muscle length and rate of change (Spindles -Ia and -II), and the measurement of the force value at the Email address: spirito.mario@gmail.com (Mario Spirito) tendon level (Ib-Spindle or Golgi Tendon Organs (GTO)). We can see the muscle structure, from a system theoretic standpoint, as a multi-input multi-output system with the actuator (α-motoneurons), the plant (the muscles cells, whose contraction produces the exerted force on the joint), and the sensors (Spindles -Ia and -II, and GTO), and the CNS as a feedback controller. A schematic of the involved neurological connections is available in [START_REF] Boyd | The isolated mammalian muscle spindle[END_REF]. terior horn of the grey matter of the spinal cord.

For about half its length the intrafusal muscle bundle is contained within a thick connective tissue capsule which expands in its central 2 mm to form the fluid space organ. Tendon organs lie 'in series' wi extrafusal muscle and, consequently, di charge during a twitch of the muscle.

If the sensory terminals within the spi dle are to send information about musc length to the brain continuously durin The modeling of how muscles exert force started with the work of [START_REF] Hill | The heat of shortening and the dynamic constants of muscle[END_REF] in which the muscle description via a spring-like behavior is presented. In particular, the relationship takes into account a contractile element that produces muscle force and it is put in series to a passive elastic element representing the tendon model between the muscle fibers and the connected bones or tissues. This model has then been exploited to generate more detailed models [START_REF] Huxley | Muscle structure and theories of contraction[END_REF] and [START_REF] Gordon | The variation in isometric tension with sarcomere length in vertebrate muscle fibres[END_REF] in which the authors describe the microscopic behavior of contraction behavior generated by the effects of sarcomeres and cross-bridges. In [START_REF] Mcruer | A neuromuscular actuation system model[END_REF], Hill's model has been enriched with a 'damping element'. In [START_REF] Winters | Analysis of fundamental human movement patterns through the use of in-depth antagonistic muscle models[END_REF], [START_REF] Winters | An improved muscle-reflex actuator for use in large-scale neuromusculoskeletal models[END_REF], and [START_REF] Brown | Measured and modeled properties of mammalian skeletal muscle. ii. the effectsof stimulus frequency on forcelength and force-velocity relationships[END_REF], experimental data has been exploited to construct an analytic expression of the total force generated by the muscles. In these works, the total force has been split into a length-dependent force F L and a velocity-dependent force F V , both depending on the α-motoneuron excitation frequency. The F L term is then split, similarly to Hill's model, into an active plus passive term. The active force describes the part of fibers contraction when excited and produces movements, while the passive force only refers to the elastic mechanical behavior of the muscle fibers.

Experimental data and models of isometric (fixed muscle length) and isokinetic (fixed muscle velocity) forces are available and presented in [START_REF] Winters | Analysis of fundamental human movement patterns through the use of in-depth antagonistic muscle models[END_REF], [START_REF] Winters | An improved muscle-reflex actuator for use in large-scale neuromusculoskeletal models[END_REF], [START_REF] Brown | Measured and modeled properties of mammalian skeletal muscle. ii. the effectsof stimulus frequency on forcelength and force-velocity relationships[END_REF], [START_REF] Joyce | The mechanical properties of cat soleus muscle during controlled lengthening and shortening movements[END_REF], and [START_REF] Rack | The effects of length and stimulus rate on tension in the isometric cat soleus muscle[END_REF]. A system theoretic perspective of the neuromuscular system has been presented in [START_REF] Mcruer | A neuromuscular actuation system model[END_REF][START_REF] Winters | Hill-based muscle models: a systems engineering perspective[END_REF][START_REF] Van Der Helm | Musculoskeletal systems with intrinsic and proprioceptive feedback[END_REF] in which they give a controloriented interpretation of the involved parts. In [START_REF] Hogan | Adaptive control of mechanical impedance by coactivation of antagonist muscles[END_REF], instead, they consider the optimal co-activation problem while minimizing the mechanical impedance of an antagonist-actuated biological joint (e.g. an elbow).

In recent years, with the development of the contraction theory in the control field, the contractile properties of the neuromuscular system (involving Hill's model equations) have been studied in [START_REF] Mcintyre | Does the brain make waves to improve stability?[END_REF]. In [START_REF] Wakeling | A muscle's force depends on the recruitment patterns of its fibers[END_REF] an analytic model of complete and incomplete tetanus contraction has been proposed, and, according to the authors' knowledge, it gives the closest to-reality dynamical model of the muscle behavior. The sensory part involved in the muscle system is related to the measurement of the plant variables, such as force at the tendon level (GTO), and muscle/fibers length and velocity (Spindles -Ia and -II). They have been studied in depth [START_REF] Boyd | The isolated mammalian muscle spindle[END_REF] and [START_REF] Boyd | Analysis of primary and secondary afferent responses to stretch during activation of the dynamic bag 1 fibre or the static bag 2 fibre in cat muscle spindles[END_REF] and the many works have been collected in [START_REF] Boyd | The muscle spindle[END_REF], see [START_REF] Windhorst | Muscle proprioceptive feedback and spinal networks[END_REF] for a more recent review on the topic. From a system engineering perspective, the GTOs dynamical model has been proposed in [START_REF] Mileusnic | Mathematical models of proprioceptors. ii. structure and function of the golgi tendon organ[END_REF] by exploiting a proper transfer function between the force exerted on the tendons and the GTO neurons' firing frequency. For Spindles -Ia and -II, in [START_REF] Maltenfort | Spindle model responsive to mixed fusimotor inputs and testable predictions of β feedback effects[END_REF], a static model has been proposed to describe the relationship among the muscle/fibers length, velocity, and the γ-motoneurons excitation to the spindles firing rate. Such γ-motoneurons are in charge of modulating the firing responses of the connected Spindle system.

To the best of the authors' knowledge, there is no dynamic model of the Spindles organs available in the literature.

The sensory part plays an important role in the movement control architecture. In particular, the presence of internal loops between the muscles and the spinal cord allows the possibility to actuate the reflexes [START_REF] Van Der Helm | Musculoskeletal systems with intrinsic and proprioceptive feedback[END_REF], [START_REF] Schomburg | Spinal sensorimotor systems and their supraspinal control[END_REF], and [START_REF] Matthews | The human stretch reflex and the motor cortex[END_REF], such as stretch re-flexes, reciprocal inhibition, and others. Finally, another relative element is the spinal cord. It connects the CNS and αmotoneurons via the interneurons and it is unavoidably involved in movement control [START_REF] Jankowska | A neuronal system of movement control via muscle spindle secondaries[END_REF][START_REF] Gassel | A critical review of evidence concerning long-loop reflexes excited by muscle afferents in man[END_REF][START_REF] Eccles | Supraspinal control of interneurones mediating spinal reflexes[END_REF].

The paper is organized as follows. We conclude the introduction with some motivations for the work, then in Section 2 we give some authors' considerations about the available data, and in Section 3 we describe models available in the literature that has been obtained exploiting these data. While in Section 4 we describe the developed models for the isometric muscle force. In Section 5 we provide a hybrid modeling of the sensory system available in the muscles. We then provide an application of muscle modeling here proposed to the description of the oculomotor system, in Section 6 and describe the challenges related to it. Section 7 concludes the work by summarizing the results and describing possible future works and advances.

Motivations

This work aims to describe neuromuscular behavior in a control theory paradigm. In particular, we would like to better understand the purposes of the neuromuscular loops and how they are exploited in controlling limbs movements, see for example [START_REF] Iqbal | Stabilizing pid controllers for a single-link biomechanical model with position, velocity, and force feedback[END_REF] for the application of a stabilizing PID action on a single-link biomechanical model or [START_REF] Kistemaker | Control of position and movement is simplified by combined muscle spindle and golgi tendon organ feedback[END_REF] to see how this available measurement feedback can be exploited in position and movement control. The reader can furthermore see [START_REF] Mcintyre | Does the brain make waves to improve stability?[END_REF] to have an insight into how waveform signals can be used to control delayed closed-loop systems, thus mimicking the spiking behavior of biological nervous systems. We strongly believe that this modeling idea would be a good starting point to eventually deeply comprehend the purposes of the brain parts involved in learning and control of dexterous movements, see for example [START_REF] Mtui | Fitzgerald's Clinical Neuroanatomy and Neuroscience E-Book[END_REF], [START_REF] Nolte | The human brain[END_REF], [START_REF] Wolpert | Internal models in the cerebellum[END_REF], [START_REF] Barak | neurophysiological basis of movement. human kinetics, champaign[END_REF], [START_REF] Kiernan | Barr's the human nervous system: an anatomical viewpoint[END_REF], and [START_REF] Patestas | A textbook of neuroanatomy[END_REF], and thus, e.g., describe and simulate pathologies effects affecting movements control. The final objective will then be to have a complete picture and understanding of the structure used by the brain in controlling human dexterous movements and thus mimic such a strategy in the field of nonlinear control systems such as nonlinear output regulation as in [START_REF] Huang | Internal models in control, biology and neuroscience[END_REF]. The study interleaves with the concept of inverse and forward models introduced in neuroscience in the works [START_REF] Wolpert | Multiple paired forward and inverse models for motor control[END_REF], [START_REF] Kawato | Internal models for motor control and trajectory planning[END_REF]. See also [START_REF] Mischiati | Internal models direct dragonfly interception steering[END_REF] for another application of these models. Thus an intermediate step of the study will be to better understand how these models are necessarily formed, trained, and updated in the brain and how are they involved in the control strategy.

Preliminaries on available experimental data in literature

In the muscle force, we can distinguish among three main components: a passive force F p that is always present due to the natural behavior of tissues and only depends on the muscle length, then we have the active static (isometric) force F L only depending on the muscle contraction (length L) and an active dynamic (eccentric and concentric) force F V also depending on the rate of change of muscle contraction L. The last two components are active terms and indeed they depend on the activation frequency α of the α-motoneuron.

Data of isometric force available in the literature are mainly those in the figures of [START_REF] Winters | An improved muscle-reflex actuator for use in large-scale neuromusculoskeletal models[END_REF], [START_REF] Brown | Measured and modeled properties of mammalian skeletal muscle. ii. the effectsof stimulus frequency on forcelength and force-velocity relationships[END_REF], and [START_REF] Hogan | Adaptive control of mechanical impedance by coactivation of antagonist muscles[END_REF], reported in this section for completeness. It is worth noticing that such data only Figure 3: Isometric Force. Relationship between the static Force and the muscle length at different activation frequencies. This plot is half useful due to the unknown length, but we can notice that there are different rest lengths (length at which the force is null) depending on the activation frequency. Plot taken from [START_REF] Hogan | Adaptive control of mechanical impedance by coactivation of antagonist muscles[END_REF].

span the activation frequency in the interval [0, 120] pps (pulse per second), because 120 pps is the maximum excitation frequency the α-motoneuron can 'apply' on the muscle fiber.

At α = 120pps the active force F L (i.e. the total force without passive F p and dynamics F V terms) has a maximum value indexed in literature by the symbol F o . This maximum value F o is associated with a specific length value L o (which is not the maximum length value, but the length L corresponding to the maximum force value F o at the maximum excitation frequency α = 120pps). Almost all available data and proposed models in the literature are normalized with respect to these two values. Indeed we find force and length data of the form F/F o and L/L o , respectively. Due to this fact, we also consider a normalized model for such normalized values.

A crucial point in our model development is the presence in the force characteristic of rest length values (i.e. values of L for which the exerted force is null) that depends on excitation rates Figure 4: Normalized Isometric Force model with correction term given by a negative passive force in [START_REF] Winters | An improved muscle-reflex actuator for use in large-scale neuromusculoskeletal models[END_REF]. Here the activation frequency is normalized with respect to its maximum value α = 120 pps α, see Fig. 3. And this fact leads us to consider a rest length function L 0 (which is not L o as above) depending only on the excitation frequency α and this might be used to describe the contracting behavior of the muscle. A major problem is that most of the data available in the literature lack these rest-length data which is fundamental for our model development. We thus exploit the same technique introduced in [START_REF] Winters | An improved muscle-reflex actuator for use in large-scale neuromusculoskeletal models[END_REF]. In practice, they first consider a model without the presence of rest length (based on Gaussian terms, as it will be clearer at the end of section 3), and then they propose to add a negative passive term to compensate for the lack of rest lengths values in the model and obtain the rest lengths points saturating the force at zero. The approach is very clear in Fig. 4.

In the next section, we provide an overview of the two models available in the literature.

Literature models comparison

In our opinion, the most relevant and complete models in the literature are those available in [START_REF] Brown | Measured and modeled properties of mammalian skeletal muscle. ii. the effectsof stimulus frequency on forcelength and force-velocity relationships[END_REF] and [START_REF] Winters | An improved muscle-reflex actuator for use in large-scale neuromusculoskeletal models[END_REF]. Here we give a short description of such models for the sake of completeness.

Brown's et al model

In [START_REF] Brown | Measured and modeled properties of mammalian skeletal muscle. ii. the effectsof stimulus frequency on forcelength and force-velocity relationships[END_REF], as the Hill's model [START_REF] Hill | The heat of shortening and the dynamic constants of muscle[END_REF], the authors consider both the active isometric and the dynamic forces generators gathered inside a contractile element F CE along which passive term F PE describes the passive force (previously indicated as F p , but we decided to stay faithful to the original nomenclature). The resulting total force F T OT presents as

F T OT = F CE + F PE (1)
and F CE is expressed as

F CE (α, L, L) = R(α)A(α, L)F L (L)F V (L, L) (2)
where R is a fiber recruitment percentage and describes the amount of fiber involved in the force generation depending on the activation frequency, A is the term relating the muscle activation to the excitation frequency α and it would have values between 0 and 1 (the same values range has been considered for R), F L and F V describe respectively the force-length relationship and the force-velocity, and their product provides the force expression of the muscle depending on the current length L, the velocity of contraction L, and excitation frequency α.

In their work, a simplified version of the force function is provided fixing R(α) = 1 and making assumptions of quasistatics (close to isometric) conditions during motion and of a starting length equal to the rest length, hence they move to a force expression of the type:

F CE (α, L, L) = A F (α, L) F L (L)F V (L, L) + F p (L) (3)
where

A f (α, L) = 1 -exp - α 0.56N f (L) N f (L) = 2.11 + 4.16 1 L - 1
F L (L) = exp - L 1.93 -1 1.03 1.87 F p (L) = -0.02e 13.8-18.7L F V (L, L) = -5.72- L -5.72+(1.38+2.09L) L , L ≤ 0 2.5+4.21L+2.67L 2 0.62+ L L, L > 0 (4)
In this model, one can notice that the length F L and velocity F V dependent force terms are multiplied one each other. This might lead to the interpretation that the velocity terms F V are only acting as scaling functions on the 'supporting' isometric force F L . We rather prefer to split the isometric F L and dynamics F V terms to analyze their contribution individually in a spring-damper fashion. Moreover, notice that when L → 0 + in the equation we have F V 1, thus the model loses the continuity with the static isometric properties of the muscle after a lengthening phase.

Winters' Model

In [START_REF] Winters | An improved muscle-reflex actuator for use in large-scale neuromusculoskeletal models[END_REF], Winters's work goes deep inside the basic muscle elements. Indeed in its work, we can find a microscopic description of the muscle components such as intrafusal and extrafusal elements with their relative damping and stiffness terms, and these are connected in parallel and series fashion. Of their work, we only focus on the isometric force that has been modeled as

F ISO (α, L CE ) = F L (α, L CE ) + F p (L CE )
where, for the isometric force model, they exploit a Gaussian function with variable mean value

F L (α, L CE ) = αe - L CE L m0 -1.05-(1-α)0.2 0.4 2 ( 5
)
where L m0 , in Winter's notation, is the optimal muscle length value L o (in this work notation),

L CE = L mt -L t0 -x SE (L m0 + L t0
) and x SE is the (dimensionless) extension of the extrafusal series of spring-element, L mt is the overall musculotendinous length and L t o is the tendon rest length. The value of α is normalized with respect to its maximum value. The passive force is instead given by the expression

F passive =          e 3 0.6 L CE L m0 -1 -1 e 3 -1 , L CE L m0 ≥ 1 -e 6 0.7 1- L CE L m0 -1 e 6 -1 , L CE L m0 < 1 (6)
Notice that, as touched above, the passive force has a negative term for L CE L m0 < 1. We strongly believe that, for our study purposes, the functions exploited in the models presented in this section are very complicated to manipulate and use in control analysis. Moreover, the meaning of the parameters is very difficult to interpret from a physical point of view.

Proposed Models

Our main objective is to describe the isometric force (F L + F p ) as a nonlinear controlled spring with variable stiffness K and rest length L 0

F ISO (α, L) = K(α, L) (L -L 0 (α)) (7)
where α is the α-motoneuron excitation frequency and L is the muscle length. In order to obtain such a final model, we first designed an intermediate/evaluation model following Winters' considerations. Indeed, we exploited the negative passive term introduced in [START_REF] Winters | An improved muscle-reflex actuator for use in large-scale neuromusculoskeletal models[END_REF] to obtain the rest length values for the different excitation frequencies. This step allows us to then construct and fit the parameters of the final/control model (7).

Evaluation Model

In order to construct the intermediate/evaluation model we consider the isometric force as divided in active F L and passive F p terms, as exploited in [START_REF] Winters | An improved muscle-reflex actuator for use in large-scale neuromusculoskeletal models[END_REF]:

F ISO (α, L) = F L (α, L) + F p (L).
Inspired, then, by the modeling approach in [START_REF] Brown | Measured and modeled properties of mammalian skeletal muscle. ii. the effectsof stimulus frequency on forcelength and force-velocity relationships[END_REF] we describe the active force curves as Gaussian functions with mean value µ and standard deviation σ and amplitude A as a function of the excitation frequency α, as shown in (4.1). By running a curve fitting algorithm we are able to fit the active isometric force with a very good order of approximation, i.e. the standard Euclidean norm of the residual vector is |r| 2 = 0.051, where r is the stack of the differences between 0.8 the available force data and the related value of the proposed model. In Fig. 5 we plotted the active force comparison between the available data from [START_REF] Brown | Measured and modeled properties of mammalian skeletal muscle. ii. the effectsof stimulus frequency on forcelength and force-velocity relationships[END_REF] and the intermediate model estimation. The identified coefficients are reported in Tab.1.

F active = Ae -   L -µ σ   2 A = c 11 α 2 + c 12 α + c 13
Exploiting the data available in [START_REF] Winters | An improved muscle-reflex actuator for use in large-scale neuromusculoskeletal models[END_REF] and their negative passive terms considerations we propose a unified continuous model of the passive force term exploiting a sigmoidal and exponential function

F p (L) = f 1 + f 2 -f 1 1 + 10 f 3 ( f 4 -L) + f 5 e f 6 L (8)
were the coefficient values f i , for i = 0, . . . , 6, can be found in Tab.2. The sigmoidal behavior in the passive term, see Fig. 9, does not have an experimental data correspondence but provides a plausible description of the physical limit of every elastic material, i.e. beyond the 'boundary' length value the solicited fibers start to break. Then the total isometric force will result from the saturation at 0 of all negative terms, as in

F ISO (α, L) = F passive (L) + F L (α, L) 0, F Iso ≤ 0 (9)

Control Model

Due to the high complexity of the intermediate/evaluation model [START_REF] Mcruer | A neuromuscular actuation system model[END_REF], and in order to make the force relationship more in- tuitive from a physical standpoint, we exploited the intermediate model to find the rest length points and then the equivalent stiffness values. These data were fundamental to then fit the rest length L 0 (α) and stiffness K(α, L) functions' parameters in terms of the α and L, to construct the final/control model. Thus, the isometric force function is given by [START_REF] Huxley | Muscle structure and theories of contraction[END_REF], reported here for completeness,

Passive force

Available Data and Models Comparison

Passive data Evaluation Model Control Model

F ISO (α, L) = K(α, L) (L -L 0 (α)) (10)
where both L 0 and K are taken as polynomial functions in α and

L L 0 (α) = 1 α 2 + 2 α + 3 K(α, L) = k 1 (L)α 3 + k 2 (L)α 2 + k 3 (L)α + k 4 (L) (11)
where the k i (•), i = 1 . . . 4, a polynomial functions in α

k i (α) = k i1 α 6 + k i2 α 5 + k i3 α 4 + k i4 α 3 + k i5 α 2 + k i6 α + k i7 . ( 12
)
The relative coefficient values are reported in Tables 3 and4, respectively. 6.794 • 10 -5 -0.01271 1

We also report in Fig. 7 a comparison between the stiffness surface K(α, L) and the constructed point from the evaluation model with respect to the α and L values. The model provides a good approximation of the Evaluation Model and, as a consequence, of the experimental data, as it is shown in Fig. 8. Moreover, the new model of the total force

Sensors modelling

In this section, we talk about the dynamic modeling of the organs available as sensors. They measure the muscle length and velocity, along with the applied force at the tendon level. Respectively, these organs are spindle II and Ia, and Ib (also called Golgi Tendon Organs). For a model of the last ones, we refer the reader to the works [START_REF] Mileusnic | Mathematical models of proprioceptors. i. control and transduction in the muscle spindle[END_REF] and [START_REF] Mileusnic | Mathematical models of proprioceptors. ii. structure and function of the golgi tendon organ[END_REF].

Spindles Ia and II

Ia-and II-spindles' neuronal activities are correlated to the length and rate of change of the muscle fibers' contractions. Along with each of the Spindles, we find active elements, i.e. the γ-motoneurons whose excitation frequency modifies the spindles' characteristics. In particular, such modulation can be interpreted as an increase in sensing precision obtained by the Spindles, despite the noise amplification drawback. We can generally distinguish between static, γ s , and dynamic, γ d , γmotoneurons. In the modeling, we only consider the effects of To construct the model of the sensors, we first find the static behavior of the sensors' outputs, we make an analysis on the slopes as similarly done in [START_REF] Maltenfort | Spindle model responsive to mixed fusimotor inputs and testable predictions of β feedback effects[END_REF].

We now first show the hybrid dynamical model found for the Spindle-Ia and then the one of the Spindle-II.

Spindle-Ia

We considered an Hybrid dynamical model and is given by

C Ia :        ẋIa = p Ia K L Ia L -x Ia ṗIa = 0 KL Ia = 0 y Ia = sat 0 x Ia + K γ Ia γ Ia + K LIa L D Ia :        (∀X Ia ∈ D Ia1 ) K + L Ia = 0, p + Ia = pIa (∀X Ia ∈ D Ia2 ) K + L Ia = KL Ia , p + Ia = pIa (∀X Ia ∈ D Ia3 ) K + L Ia = 0, p + Ia = p Ia (∀X Ia ∈ D Ia ) x + Ia = x Ia , y + Ia = y Ia (13)
where x Ia is the filter state, y Ia is its output, the function sat 0 is the function that saturates at zero negative values, p Ia is the filter pole, K L Ia , K V Ia and K γ Ia are the static gains for the muscle In Fig. 10 we compare the output of the dynamical system [START_REF] Joyce | The mechanical properties of cat soleus muscle during controlled lengthening and shortening movements[END_REF] with the available experimental data. The identified parameters are all gathered in Tab.5. where x II is the filter state, y II is system output, the function sat 0 is the function that saturates at zero any negative argument, p II is the filter pole, K L II , K V II and K γ II are the static gains for the muscle length, velocity and the excitation state of γ d -motoneuron; in this case, L, V and γ II are the muscle length, velocity and the excitation state of γ d -motoneuron and are considered as system inputs. Then X II is the vector collecting the state, inputs, and output of the system. The flow set

C II :            ẋII = p II ξ L II K L II L -x II ṗII = 0 ξL II = 0 y II = sat 0 x II + K γ II γ II + K LII L+ 1 -ξ L II K L II L -L D II :        (∀X II ∈ D II1 ) ξ + L II = 0, p + II = pII (∀X II ∈ D II2 ) ξ + L II = 1, p + II = pII (∀X II ∈ D II3 ) ξ + L II = 0, p + II = p II (∀X II ∈ D II ) x + II = x II , y + II = y II (14)
C II = R 3 \D II while the jump set is D II = D II1 ∪ D II2 ∪ D II3
, where

D II1 = {X II : L < 0 & (ξ L II 0||p II pII )}, D II2 = {X II : L > 0 & (ξ L II 1||p II pII )}, D II3 = {X II : L = 0 & (ξ L II 0||p II p II )}.
In Fig. 11 we make a comparison between the output of the dynamical system [START_REF] Rack | The effects of length and stimulus rate on tension in the isometric cat soleus muscle[END_REF] with the experimental available. The identified parameters are all gathered in Tab.5. In this section, we provide an application of the proposed muscle modeling in the framework of the Oculomotor System. In particular, this application has been motivated by the work [START_REF] Broucke | Model of the oculomotor system based on adaptive internal models[END_REF], [START_REF] Broucke | Adaptive internal model theory of the oculomotor system and the cerebellum[END_REF], [START_REF] Battle | Adaptive internal models in the optokinetic system[END_REF], in which a description of the oculomotor system has been provided with the assumption that the actuation is not given by muscular torque but rather given directly as an acceleration, see equation (11b) in [START_REF] Battle | Adaptive internal models in the optokinetic system[END_REF] for example. This, in our opinion, oversimplifies the tackled problem implying that a more complex analysis must be pursued when considering an adaptive control description of the system.

We believe thus that a more insightful model has to take into account the skeletal muscle activation due to the firing rate of the α-neurons α i , i = 1, 2, and that their exert force then generates an equivalent torque on the eyeball. Thus, according to the labels in Fig. 12, we propose the following model for the horizontal motion of an eyeball θ = ω

J ω = -β ω + RK 1 (α 1 , L 1 ) L 1 -L 01 (α 1 ) -RK 2 (α 2 , L 2 ) L 2 -L 02 (α 2 ) ( 15
)
where θ is the angular displacement with respect to the sagittal (or longitudinal) plane, ω is its rate of change, J is the eyeball inertia with respect to the vertical axis, R is its radius, and α i , K 1 , L i , L 0i , for i = 1, 2, are the α-motoneurons firing rates, the muscles variable stiffness, current and rest lengths, respectively. Moreover, the length of each muscle is related to the angle θ , by considering the following relationship, i = 1, 2,

L i = L i0 + (-1) i Rθ
where L i0 are the initial muscle length. And so the model can be written as a function of the system state variable θ = ω

J ω = -β ω + RK 1 (α 1 , L 10 -Rθ ) L 10 -Rθ -L 01 (α 1 ) -RK 2 (α 2 , L 20 + Rθ ) L 20 + Rθ -L 02 (α 2 ) (16)
With the proposed model, we show that such an assumption does not hold any more without additional assumptions on the muscle force generation, e.g., a part of the brain is involved in providing an inversion map between the desired acceleration torque. And it is worth to notice that such an inverse map is not possible to get because the forward mapping from the firing rates α 1 and α 2 to the torque τ

= RK 1 (α 1 , L 1 ) L 1 -L 01 (α 1 ) - RK 2 (α 2 , L 2 ) L 2 -L 02 (α 2 )
is not injective but only surjective. Indeed, in particular, one can exploit the additional degree of freedom to obtain a desired joint stiffness for example.

Moreover, the oculomotor system involves also the sensory part described in sec.5 (one for each muscle). In particular, assuming that L 10 = L 20 the scaled difference between L 2 and L 1 provides the angle θ , while the difference between L2 and L1 provides its rate of change. These signals can be obtained from the response of the spindle Ia and II, of each muscle. Moreover, this information is in synergy with that coming from the visual cortex, when, for example pursuing an object, we also have an additional error signal reconstructed by the visual cortex itself. However, the information provided to the brain-implemented controller is asynchronous, in the sense that the time delay affecting the feedback of these signals is not the same, in general, for the two types of information, i.e., that coming from the sensory system and that coming from the visual system.

The challenge is then to obtain a robust (to the uncertainties of the parameters) control strategy that reconstructs the system state from the Hybrid model of the Spindle Ia and II and guides the firing rates of the two muscles to achieve the desired behavior, such as object pursuing or gaze fixation to name two.

Conclusions and future works

In this work, we dealt with the modeling of muscle force considering the passive and active components of the isometric force. We first provided a brief overview of the model available in the literature and then show two models,i.e., an intermediate and a final one. The first one is more complex and precise than the second but it does not add much to the models available in the literature. On the other hand, it allowed us to construct a set of data that have been useful for the identification of the second model parameters. The latter (the spring model) is easier to handle, although less precise than the intermediate one, and gives a physical intuition about the muscle behavior. Indeed, it describes such muscular behavior as a controllable non-linear spring with controllable (variable) stiffness K and rest length L 0 . We furthermore provide a hybrid dynamical model formulation of the spindles organs that act as muscle length and velocity sensors for the Central Nervous System. We then provide an application to the case of the oculomotor system proposing its dynamic modeling, and describing the related challenges.

Future works may concern the development of a muscle force characteristic depending also on the rate of shortening and lengthening of the muscle fibers' length, to provide a more complete model. Such a model can then be analyzed in a control context involving (biological) joints with antagonist actuation. Moreover, a system theoretic description of the physiological reflexes, such as the stretch reflex, can be investigated by exploiting the proposed model.

~

  Group la~! ~("~'-'~,, i

Fig. 1 .

 1 Fig. 1. The ingoing (afferent) and outgoing (efJerent) pathways connecting the muscle spindle and the spinal corThe pathways carrying information from muscle spindles from the spinal cord to the brain are not shown.

Figure 1 :

 1 Figure 1: Muscle and Central Nervous System schematics of the connections. Picture taken from [5].

Figure 2 :

 2 Figure 2: Normalized Isometric Force. Relationship between the static Force and the muscle length at different activation frequencies. Plot taken from [12].

α 2 +

 2 c 14 α + c 15 µ = c 21 e c 22 α + c 23 e c 24 α σ = c 31 + c 32 cos(c 34 α) + c 33 sin(c 34 α)

Figure 6 :

 6 Figure 6: Passive force comparisons among the available data in green, the Evaluation model in blue, and the control model in magenta. The residual values for the Evaluation model has norm |r EM | 2 = 0.0077 and for the Control model the norm is |r CM | 2 = 0.1706

Figure 7 :

 7 Figure 7: Fitting of the stiffness K(α, L) surface with respect to α and L

Figure 8 :

 8 Figure 8: Data comparison with Spring Model. r indicates the residuals vector.

Figure 9 :

 9 Figure 9: Comparison between the Evaluation and the Control Model for different excitation frequencies, i.e. α ∈ 0, 30, 60, 120 pps.

  (a) Spindle Ia simulation with hybrid dynamical system (b) Experimental Data of Spindle Ia from [21][Ch.19]

Figure 10 :

 10 Figure 10: Spindle Ia Model output with experimental data in [21][Ch.19] comparison

  (a) Spindle II simulation with hybrid dynamical system (b) Experimental Data of Spindle II in [21][Ch.19]

Figure 11 :

 11 Figure 11: Spindle II Model output with experimental data in [21][Ch.19] comparison

Figure 12 :

 12 Figure12: Schematic description of the oculomotor system in the horizontal plane evolution. The oculomotor system presents two antagonist muscles whose combined force produces a torque on the eyeball.

Table 1 :

 1 Parameters of Evaluation Model active force.

		Available Data -E.M. comparison: |r| 2 =0.050687
		1
					Data15pps
		0.8			Data30pps
	Active force	0.4 0.6			Data40pps Data60pps Data120pps EM15pps EM30pps
					EM40pps
		0.2			EM60pps
					EM120pps
		0
			1	1.2	1.4	1.6
				Length
	i	c i1	c i2	c i3	c i4	c i5
	1	1.134	0.723	0.0504	11.3297 643.259
	2	0.574	-0.0211	0.8477	0.0008	n.a.
	3 0.4518 -0.01704 -0.1235 0.04898	n.a.

Figure 5: Data comparison between the available data in [12] and the evaluation/intermediate model. r indicates the residuals vector.

Table 2 :

 2 Parameters of Evaluation Model passive force.

	f 1	f 2	f 3	f 4	f 5	f 6
	0.00381 2.16582 8.16634 1.45388 -7.32 -7.4

Table 3 :

 3 Parameters of Control Model rest legth.

	1	2	3

Table 4 :

 4 Parameters of Control Model stiffness.

	i	k i1	k i2	k i3	k i4	k i5	k i6	k i7
	1 -0.00041 0.0025	-0.0061	0.00763	-0.0052	0.0018	-0.00026
	2	0.0586	-0.352	0.8444	-1.0277	0.6697	-0.2257	0.0319
	3 -1.67504	9.829	-22.692	25.896	-15.268	4.5087	-0.5351
	4	-40.977 244.666 -569.9335 664.1912 -405.645 120.853	-13.045

Table 5 :

 5 Parameters of spindle Ia model.

	pIa	p Ia	KL Ia K LIa	K γ Ia
	85 3.858 7.15 5.15 0.3514
	5.1.2. Spindle-II
	Also for this sensing system a Hybrid dynamical model and
	is given, similarly, by

Table 6 :

 6 Parameters of spindle II model.

	pII	p II	KL II	K LII	K γ II	L
	85 2.858 4.334 0.9184 0.343 2.5


Michael Ovando

Jean-Bernard Blaisot

Françoise Baillot

Françoise Baillot Nonlinear

Franc ¸oise Baillot

Nonlinear acoustic theoretical model considering acoustic transmissivity in the interior of a real fluid object

à 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

The context of this study relies on the thermoacoustic instabilities occuring in Liquid Rocket Engines (LRE). The instabilities results from the coupling between heat release rate fluctuations, q ′ (t) and pressure fluctuations, p ′ (t) which usually tune to the chamber's acoustic eigenmodes (Oefelein and Yang, [START_REF] Oefelein | Comprehensive review of liquid propellant combustion instabilities in f-1 engines[END_REF]). To understand how these instabilities modify the flow, an approach based on the analysis of the response of a fluid flow to acoustic solicitations is considered. In experiments, this is done by using an acoustic forcing acting on the fluid flow that models injection in engine conditions. In this paper we present an extension of a semianalytic approach to model the interaction between a liquid or a surpercritical fluid with acoustics.

A theoretical model capable of describing the effects of the acoustic radiation pressure, considering a non-perfectly reflective object (as opposed to how it is seen in the works of Ficuciello [START_REF] Ficuciello | [END_REF] and Herrera [START_REF] Herrera | [END_REF]) is presented. It is based on the approach of Yosioka and Kawasima [START_REF] Yosioka | Acoustic radiation pressure on a compressible sphere[END_REF] that expressed the acoustic radiation force on spheres by considering the effects of acoustic transmission inside the object medium.

The next section is devoted to the presentation of the mathematical model. Then follows a presentation of the results when considering a practical injection condition typical of supercritical conditions. A conclusion ends this paper.

Acoustic radiation pressure

A potential acoustic velocity field is considered:

- → u = - - → ∇ Φ ext .
With this convention (Yosioka and Kawasima [START_REF] Yosioka | Acoustic radiation pressure on a compressible sphere[END_REF]), the pressure variation in a compressible fluid of density ρ 0 and speed of sound c 0 , where acoustic waves are propagated, can be expressed as given by * Corresponding author: blaisot@coria.fr Proceedings of the European Combustion Meeting 2023 eq. 1 following a development similar to the one of Ficuciello [START_REF] Ficuciello | [END_REF].

p ext = ρ 0 Φext - 1 2 ρ 0 ( ---→ ∇Φ ext ) 2 + ρ 0 2c 2 0 Φ2 ext (1)
The incident field Φ i is scattered in the presence of an object leading to the scattered field Φ s . The computation of the resulting field Φ ext = Φ i + Φ s , is required to estimate the radiation pressure field acting on the surface of the object and the resulting radiation force. Conventionally, the acoustic radiation pressure is expressed as the time average of equation 1 over a time period of the incident acoustic field.

p rad = 1 T T 0 (ρ 0 Φext - 1 2 ρ 0 ( ---→ ∇Φ ext ) 2 + ρ 0 2c 2 0 Φ2 ext )dt = p ζ + p Φ + p q (2)
where p ζ is the contribution due to the motion of the object, p Φ is the time-average volumetric kinetic energy density, p q is the time-average volumetric potential energy density, and Φ ext is the total velocity potential of the surrounding media Since we consider the object to be free to move under the influence of the acoustic field, the derivative of the velocity potential referred to the moving origin is given by:

Φext = d U Φ ext dt - - → U • - → ∇ Φ ext (3)
where -→ U is the barycentric velocity of the object. This velocity can be obtained by applying Newton's second law to the object. The external forces, reduced here to the effect of acoustics, are derived from integration over the object boundary, i.e.:

m - → U = - S p - → n dS = S -Φext ρ 0 - → n dS
Leading to:

- → U = ρ 0 m S -Φ ext - → n dS
In the case of a spherical object of radius r 1 and due to the symmetry properties, this results in :

- → U = 3η 2r 1 π 0 -Φ ext cos(θ ) sin(θ )dθ - → e x ( 4
)
where η = ρ 0 ρ 1 is the ratio between the density of the media and the density of the object. On this study, the index 1 indicates parameters concerning the object, whilst the index 0 indicates the fluid outside of the object. Doing the time average of equation 3, the term d U Φ ext dt reduces to zero and:

Φext = - - → U • - → ∇ Φ ext (5)
Equation 5 is used in the calculations of p ζ as:

p ζ = ρ 0 Φext = -ρ 0 - → U • - → ∇ Φ ext (6)
For p q :

p q = ρ 0 2c 2 Φ2 ext (7)
And finally, the term p φ :

p φ = 1 2 ρ 0 ( - → ∇ Φ ext ) 2 (8)
The combination of equations 6, 7 and 8 allows to estimate the acoustic radiation pressure p rad at the surface of the object.

The approach of Yosioka and Kawasima is used here to evaluate the effect of acoustics under trans-critical or supercritical conditions. Indeed, in such conditions, the object interface can no longer be considered as perfectly reflective and the transmission of the acoustical wave inside the object must be taken into account. This is done through the boundary conditions used to compute p rad .

Modelling of velocity potential fields and boundary conditions

We consider a sphere of a given radius r 1 placed at the origin of the coordinate system. The geometry of the problem is reported in figure 1 Figure 1: Fluid object (1) with a rigid boundary in a surrounding fluid (0), submitted to an acoustic field. In the case of this study, we consider a standing field The incident velocity potential Φ i of the plane standing waves is expressed by:

Φ i = φ 0 e iωt 2 e -ik 0 (x+h) + e ik 0 (x+h) (9)
where x = r cos(θ ) and φ 0 = p a ρ 0 ω is the velocity potential field amplitude, with p a the amplitude of the associated incident acoustic pressure field, and ω the angular frequency. The velocity potential outside the sphere is written as:

Φ ext = Φ i + Φ s
where Φ s is the velocity potential of the scattered waves that is expressed as:

Φ s = φ 0 e iωt ∞ ∑ n=0 (2n + 1) (-i) n A n δ n H (2) n (k 0 r)P n (cosθ ) (10) being H (2)
n the spherical Hankel function of the second kind and P n is the Legendre polynomials. A n are constants to be determined from the boundary conditions. The term δ n = (-1) n e ikh + e -ikh relates to the position of the object in the acoustic field. Yosioka and Kawasima proposed the following expression for the velocity potential Φ int inside a sphere:

Φ int = φ 0 e iωt ∞ ∑ n=0 (2n + 1) (-i) n B n δ n J n (k 1 r)P n (cosθ ) (11)
where k 1 = ω/c 1 is the wave number inside the sphere. B n are other constants determined as A n from boundary conditions, i.e. by considering the continuity of the radial velocity at the object surface (Eq. 12) as well as the continuity of the pressure (Eq. 13).

∂ Φ ext ∂ r = ∂ Φ int ∂ r at r = r 1 ( 12
)
ρ 0 Φext = ρ 1 Φint at r = r 1 ( 13
)
As velocity potential fields are expressed in a base of Legendre polynomials, equations 12 and 13 leads to the following expressions for A n and B n :

A n = ηk 1 J n (k 0 r 1 ) J ′ n (k 1 r 1 ) -k 0 J ′ n (k 0 r 1 ) J n (k 1 r 1 ) k 0 J n (k 1 r 1 ) H ′ (2) n (k 0 r 1 ) -ηk 1 J ′ n (k 1 r 1 ) H (2) n (k 0 r 1 ) (14)
B n = ik 0 / (k 0 r 1 ) 2 k 1 J ′ n (k 1 r 1 ) H (2) n (k 0 r 1 ) -λ k 0 J n (k 1 r 1 ) H (2) ′ n (k 0 r 1 ) (15)
Python implementation All the expressions are implemented in a library acrad.py. Since the expressions are given by infinite sums, the convergence is studied to determine the minimum number of terms N min to be considered in the computation. We introduce the reduced acoustic radiation pressure p * rad (θ ) = p rad (θ ) P , where P = P ac ρ 0 c 2 0 is the mean acoustic energy per unit volume carried by a standing wave in an environment free from any obstacles. The expression used to calculate the convergence is:

error(n) = max 0≤θ ≤2π |p * n+1 (θ ) -p * n (θ )|
The evolution of this error with respect to n is shown in figure 2. The error has been calculated at a velocity antinode (VAN) position in the standing wave acoustic field. We can see a decrease of the error as N min increases. When N min = 3, the error is of the order of 10 -4 . The contribution of the terms over N min > 8 is under the floating point precision of the Python code (about 10 -15 ), therefore not contributing to the sum, but just adding computational time to the code. With these results, the truncated number of terms considered for the calculations presented further is n = 5, since the error calculated for this value is in the order of 10 -8 , which is a good compromise between precision and computational cost.

Studying case

The object is placed at a velocity antinode (VAN) where it was previously shown that maximum deformations occurred under the effect of acoustics [START_REF] Ficuciello | [END_REF][START_REF] Herrera | [END_REF] As indicated in table 1, different equations of state are considered for the calculations of the fluid properties, depending on the fluid considered. For the case of nitrogen, the ideal gas law is used. For the case of ethane, since the studying case is just above its critical point (P c = 4.87 MPa, T c = 305.5 K), the Peng-Robinson equation of state is used:

P = ρRT M -bρ - aαρ 2 M 2 + buρM + vρ 2 (16) with α = 1 + 0.48 + 1.574Ω -0.176Ω 2 T T c 2 , a = 0.42748 R 2 T 2 c P c , b = 0.08664 RT c P c u = 2, v = -ρ 2 ,
and Ω the acentric factor. The velocity potential fields Φ ext and Φ int are given in figure 3 for the thermodynamic conditions given in table 1 and at a time t = t 0 corresponding to the maximum of the instantaneous acoustic velocity at the velocity antinode (VAN). These values represent the amplitude of the velocity potential fields with respect to the maximum amplitude of the incident field when there is no object. On the region of observation, the variation of Φ ext is about 10% of the amplitude of Φ i . This representation is used to visualize all the involved fields, and following this to calculate each one of the terms of p rad . On figure 4 we can see the representation of the squared normalized gradient modulus of the external and internal fields at time t = t 0 . This magnitude is used to represent the volumetric kinetic energy of the system (see equation 8), and since it does not involve the term Φext , it can be spatially represented. The terms which involve Φext (see equations 6 and 7) can only be calculated on the boundary of the object.

Comparison with prior models

The actual model is compared to the one developed by Ficuciello for perfectly reflective objects. On figure 5, the comparison of p * ζ , p * Φ , p * q and p * rad as functions of θ is shown for our model (expressed with "trans" for "transmissive") and for Ficuciello's model ("reflect", for "reflective") along the θ axis, for the same physical conditions. For both models, p q equals zero or about. The other two components p φ and p ζ , have similar distributions over the object surface but the levels are different between the two models. Globally, the levels are higher for the reflective model with p φ negative and p ζ positive, resulting in a distribution for p rad varying around zero with positive maximums and negative minimums. In the case of the transmissive model, p φ is also negative and p ζ positive, but as the levels are lower for both component, it results in negative values for p rad all over the object surface. Also, the amplitude of the variation of p rad over the object surface is larger for the transmissive model. This indicates that the stretching effects that lead the object deformation (not considered here) could be more pronounced with the transmissive model than expected from the reflective model.

Conclusions

The study considered an object at a fixed position in a standing wave acoustic field (velocity antinode, VAN). This position was chosen for its interest relative to potential deformation effects reported in the literature. The implementation of transmission of acoustics in the interior of an object yields a noticeable modification of the prediction of the time-average acoustic radiation pressure distribution over the object surface when compared to reflective models.

A studying case related to engine application was considered, where two fluids: ethane (as the object) and nitrogen (as the surrounding media) are submitted to a standing plane acoustic field at high pressure and temperature conditions. The acoustic radiation pressure was calculated and compared for the case of a perfectly reflective object. Results show a shift of the time-average acoustic radiation pressure with only negative values, and also an increase of the amplitude of the distribution of this acoustic radiation pressure over the object surface. These effects come mainly from the p ζ term, which is obtained by the barycentric velocity approach we took into consideration. Further studies involving density and speed of sound dependance of the equations, as well as different spatial positioning of the object is currently being studied.

On the future, deformation of the object is expected to be implemented, by considering different boundary conditions.

Figure 2 :

 2 Figure 2: Maximum error on the infinite sums for the calculation of the acoustic radiation pressure evolving in the number N min of terms

Figure 3 :

 3 Figure 3: Normalized velocity potential field Φ/φ 0 on both the interior (ethane) and exterior (nitrogen) of a sphere at a velocity antinode (VAN), η = 0.4227, T = 310 K and P = 50 atm. Time t = t 0 corresponds to the maximum instantaneous acoustic velocity at VAN.

Figure 4 :

 4 Figure 4: Squared normalized gradient modulus of the exter-

Figure 5 :

 5 Figure 5: Representation of the normalized acoustic radiation pressure p * rad and its components p * q , p * Φ and p * ζ for η = 0.4227 by two models: perfectively acoustically reflective object and acoustically transmissive object. The bold curves indicate the acoustic radiation pressure for the perfectly reflective (in black) and transmissive (in red) object assumption.

Table 1 :

 1 . The physical parameters used in the calculations are resumed in table 1. Physical parameters used in the simulations in this chapter

		General properties
	T	310 K	P	5 MPa
	r 1	0.003 m	f 0	1000 Hz
	Acoustic position	VAN (λ /4)	P	1000 Pa
	Infinite sum terms n	10	Canvas size	250 x 100 pixels
		Fluid properties
	Ethane		Nitrogen
	Equation of state	Peng-Robinson Equation of state	Ideal gas
	c 0	358.82 m/s	c 1	108.44 m/s
	ρ 0	54.34 kg/m 3	ρ 1	128.53 kg/m 3
	k 0	17.51 m -1	k 1	57.94 m -1

Acknowledgements

The authors would like to acknowledge the Centre National des Études Spatiales (CNES) and the Normandy Region for their funding support.


Joseph Larmarange

Accès à la prophylaxie préexposition (PrEP) en Afrique de l'Ouest : nouvelles questions et enjeux

CONTEXTE

➔ Critères de « comportement à risque » minimum › Ipergay:

deux rapports anaux non protégés dans les 6 derniers mois › Recommandations du groupe d'experts en France :

- ??

ÉVOLUTION DES NORMES ET PRATIQUES ?

» Quelle appropriation de cet outil dans les discours et pratiques ?

› par les utilisateurs et utilisatrices de PrEP ? celleux ne la prenant pas ?

ÉVOLUTION DES NORMES ET PRATIQUES ?

» Quelle appropriation de cet outil dans les discours et pratiques ?

» Que signifie « se protéger » dans un contexte de prévention combinée ?

› Quelle complémentarité entre les différents outils ?

› Comment mesurer dans nos études les comportements préventifs ?

› Prendre en compte les trajectoires individuelles et les parcours de santé : recours alternatif aux outils selon les contextes, les partenaires et les moments de la vie ÉVOLUTION DES NORMES ET PRATIQUES ?

» Quelle appropriation de cet outil dans les discours et pratiques ?

» Que signifie « se protéger » dans un contexte de prévention combinée ?

» La PrEP peut-elle remodeler les relations sociales autour de la sexualité ?

› Outil asymétrique de prévention, permets une plus grande autonomie dans sa gestion personnelle des prises de risque › Influe en même temps les modes de négociation entre partenaires sexuels › Parallèle à faire avec les travaux sur la contraception, notamment la pilule ÉVOLUTION DES NORMES ET PRATIQUES ?

» Quelle appropriation de cet outil dans les discours et pratiques ?

» Que signifie « se protéger » dans un contexte de prévention combinée ?

» La PrEP peut-elle remodeler les relations sociales autour de la sexualité ?

» Quel impact (ou non) sur le marché sexuel ?

› Prise en compte de ce critère dans le choix de ses partenaires sexuels ?

› Évolution des « tarifs » du travail du sexe ?

ÉVOLUTION DES NORMES ET PRATIQUES ?

» Quelle appropriation de cet outil dans les discours et pratiques ?

» Que signifie « se protéger » dans un contexte de prévention combinée ?

» La PrEP peut-elle remodeler les relations sociales autour de la sexualité ?

» Quel impact (ou non) sur le marché sexuel ?

» Bio-médicalisation croissante de la sexualité

»

  Dernière décennie : biomédicalisation accrue de la prévention › Circoncision masculine › Traitement comme moyen de prévention (TasP) › Prophylaxie préexposition (PrEP) › Recherches sur les microbicides, la vaccination… » fondée en partie sur des concepts déjà connus › Rôle des ARVs pour réduire la transmission mère-enfant par exemple › Prophylaxie du paludisme › …. » en lien avec le sentiment de certains acteurs d'un « échec » des approches comportementales › mise en avant d'une efficacité limitée du préservatif AVANT 2014 : DES RÉSULTATS EN DEMI-TEINTE UNE EFFICACITÉ FONCTION DE L'ADHÉRENCE

  plusieurs indications combinant rapports non protégés et partenaires d'un groupe à prévalence élevée ou facteurs physiques d'augmentation du risque » Et pour celles et ceux souhaitant changer d'outil de prévention ? › En soit, la démarche de faire une demande de PrEP est un indicateur d'exposition L'accès doit-il être fondé sur un échec d'autres outils de prévention plutôt que sur l' « exposition » au VIH ? La PrEP nécessite de repenser le suivi et le parcours de soins des VIH-Le dépistage est crucial point d'entrée d'une offre diversifiée de santé une majorité des programmes sont focalisés sur le seul suivi des VIH+ dépistage centré sur l'identification des « nouveaux cas » QUELLE OFFRE DE PREP ? » Risque : des programmes PrEP avec services associés › Quid des non utilisateurs de PrEP ? Peuvent-ils avoir accès aux autres services ? › Ne pas interférer avec ce qui marche déjà If it's not broken, don't fix it! Ne pas conditionner l'accès aux services à la mise sous PrEP SE RECENTRER SUR LES BESOINS DES PERSONNES ? » Privilégier une offre globale en santé sexuelle, dont la PrEP ne serait que l'un des outils mis à disposition › Offrir du choix (one size doesn't fit all) › Un comptoir unique » Simplifier l'accès aux services › D'autant plus pour des personnes qui ne sont pas malades et/ou qui ont d'autres priorités » Adapter les services aux populations activités d'outreach et en clinique Jusqu'où peut-on aller dans les services proposés en mobile ?Un suivi régulier trimestriel est-il possible ? Peut-on atteindre des personnes qui n'iraient pas en clinique ? Sans être trop intrusif ? Coût acceptable ? Peut-on convaincre les bailleurs ?

›

  Tendance déjà existante mais amplifiée par les nouveaux outils de prévention › Intrusion croissante de la santé publique dans la sphère privée › Enjeux éthiques UN ENTHOUSIASME À NUANCER » Un nouvel outil efficace, mais pas forcément adapté à toutes les populations / tous les contextes » Pas si simple d'emploi : › Suivi médical régulier › Nécessite une bonne adhérence » Outil de prévention du VIH, mais ne protège › ni des autres infections sexuellement transmissibles › ni des grossesses non désirées » Au-delà des essais d'efficacité, nécessité des documenter les changements sociaux induits. La PrEP n'est pas la « solution miracle » pour éliminer la transmission du VIH.UNE OPPORTUNITÉ À SAISIR ?La mise en oeuvre de la PrEP est une occasion de faire évoluer les systèmes de santé, tant en matière de › services proposés › modèles de soins › que de politiques publiques à condition de › ne penser uniquement PrEP ou VIH › recentrer l'offre de santé sur l'ensemble des besoins de chaque population › décompartimenter les activités et les financements MERCI joseph.larmarange@ceped.org

incidence > 3% › Demande de recherches opérationnelles pour guider la mise en oeuvre effective de cette recommandation POUR QUI ?

	Partners PrEP
	81% adhérence /
		75% efficacité
		TDF2
		79% adhérence /
	iPrEx	62% efficacité
	51% adhérence / 44% efficacité	Bangkok 67% adhérence /
		49% efficacité
	» Situation hétérogène d'un contexte épidémiologique à un autre :
	› Consommateurs de drogues, migrants, hommes…
	» Les populations que l'on cherche à atteindre sont-elles celles qui
	fréquentent les structures déjà existantes ?

Adhérence (%) Protection contre le VIH FIN 2014 : IPERGAY ET PROUD » 2 essais au Nord auprès d'homosexuels masculins › Efficacité de 86 % (réduction incidence VIH) › PrEP continue (PROUD) ou à la demande (Ipergay) › Résultats annoncés fin 2014 et présentés à la CROI en février 2015 » Septembre 2015 : nouvelles recommandations OMS › Recommandation de la PrEP orale (à base de Tenofovir) › Populations à « risque substantiel », définies « provisoirement » comme

Homosexuels masculins Travailleuses du sexe Jeunes femmes Catégories épidémiologiques d'exposition ≠ groupes sociaux constitués

	QUELS CRITÈRES D'ÉLIGIBILITÉ ?
	» Ce critère de >3% d'incidence pose question
	› Critère populationnel et non individuel
	› Sous-groupes « à risque » parmi les populations « clés » ?


François Dabis

Marie-Louise Newell

Deenan Pillay

Collins Iwuji

Joanna Orne-Gliemann

Till Bärnighausen

Eric Balestre

Sylvie Boyer

Alexandra Calmy

Vincent Calvez

Marie-Laure Chaix

Rosemary Dray- Spira

El Kamal

Kenneth Farouki

Kobus Freedberg

John Herbst

Sophie Imrie

Joseph Karcher

France Larmarange

Richard Lert

Thembisa Lessells

Anne-Geneviève Makowa

Laura Marcelin

Nuala March

Kevi Mcgrath

Colin Naidu

Nonhlanhla Newell

Tulio Okesola

Melanie De Oliveira

Tamsen Plazy

Bruno Rochat

Frank Spire

Rodolphe Tanser

Johannes Thiébaut

Thembelile Viljoen

Zuma

 Trial participants  Africa Centre staff  Traditional Authority

HIV testing of all adult members of a community, followed by immediate ART initiation of all, or nearly all, HIV-infected participants regardless of immunological or clinical staging



  i.e. Universal Test and Treat will prevent onward transmission and reduce HIV incidence in this population TasP trial area Sub-district : Hlabisa Region : KwaZulu-Natal  Zulu speaking people  HIV prevalence is ~30%  Unemployment is ~80%  3 DoH clinics located in the trial area 4 clusters opened in 6 clusters opened in 12 clusters opened in TasP trial design  TasP is a cluster randomized trial.  Each cluster has a population of approx. 1250 adults (16+ years).  The TasP intervention has 2 components: "universal" repeat testing (all clusters) + early treatment (intervention cluster)  In each cluster, rounds of home-based HIV testing repeated every ~ 6 months  All HIV+ identified participants are referred to local TasP clinics (at least one clinic per cluster) Home-based HIV testing  Clinic-based HIV testing for individuals not wanted to test at home  Clinic-based ART  ART counselling  In June 2013, added:  Phone call / Home visit if not linked to TasP clinics within 3 months  Phone call by nurse in case of a missed appointment in TasP clinic  A 'tracking team' for patients not reached by the nurse or requiring additional support



  Started in January 2015 (as in DoH clinics)  Affects mainly the clusters opened in 2014  Conditional power estimates done in November 2014 testing several scenarios (taking into account the new guideline):  power to detect an impact on incidence range from 40 to 95% (depends mainly on the baseline incidence in the 12 clusters opened in 2014)  DSMB concluded that there is no reason to stop this study for futility at this point  Recommendation of implementing new interventions to improve the cascade of HIV care and treatment  May 2015:  Interim analysis among the ten first clusters  DSMB will look at it next week Ongoing discussion  Implementation of a new interventions package?  Community level: improving community engagement and mobilization package  HIV testing: combination of repeat home-based testing, mobile testing and clinic testing  Linkage to care: SMS reminders, counselling and motivational support (phone calls, face-to-face visits), escort to clinics, health navigators, clinical assessment at home, option between home-based or clinic-based ART initiation  Retention in care: SMS reminders, phone calls and home visits in case of missed appointment, simplified ART care for stable patients, additional health services in clinics Ongoing discussion  How to implement this new interventions package?  Contingent upon May 12 DSMB conclusions  Starting during the 2 nd semester of 2015?  Extending follow-up until mid-2018? Conclusions  Real world is complex  Impact of interventions on the cascade, and therefore on HIV incidence takes time  Can we reach 90-90-90 at population level in real world?  Which services/tools to propose in a comprehensive HIV interventions package?  What would be the impact on HIV incidence?  At what cost?

Proportion of all HIV+ on ART in end 2011 39% ART coverage at the beginning of the trial

	Model parameters versus observed estimates The 90-90-90 target
				67% [63-71] Preliminary results 38.1%
	Test acceptance among those offered Linkage to care upon diagnosis among those	80% 70%	HIV ascertainment rate/CR Entry into care within 6 months among  10 clusters 77% [74-80] [44-52] 48% (27.5%  As of 31 st May 2014 without TasP)
	accepting the test		individuals not in care	 Taking into account
					39% both TasP and DoH
					[36-42] clinics
	HIV prevalence in end-2011 (15 years +)	24%	HIV prevalence (first DBS)	30% [29-31]  Including non-observed HIV+ 72.9%
	HIV incidence in end 2011	2.4 / 100 PY Observed HIV incidence	2.35 / 100 PY living in the area
	(15 years +)				[1.40-3.31]

Presented by Iwuji et al. at Melbourne (July 2014) Overall Implementation the 500 CD4 guideline

Acknowledgements