Leonardo Sanches
email: leonardo.sanches@isae.fr

Guilhem Michon
email: guilhem.michon@isae.fr

INSTABILITY ZONES IDENTIFICATION IN MULTIBLADE ROTOR DYNAMICS

Keywords: Ground Resonance, Nonlinear Dynamics, Floquet Method, Multiple Scales Method

Helicopter ground resonance is an unstable dynamic phenomenon which can lead to the total destruction of the aircraft during take-off or landing phases. Studies have been developed by researchers which considered a simplified mathematical model. With the goal of further comprehend the phenomenon, predictions of unstable motions are done and compared. First, Floquet's theory is applied to solve the linear equations of motion including parametric and periodic terms. Then, the multiple scales method is applied on the nonlinear model. The analyses highlight that, by keeping the nonlinear terms in the equations, other instability zones are identified.

INTRODUCTION

The ground resonance phenomenon in helicopters consists of a self-excited oscillation caused by the interaction of the rotor blades' lagging motion with other helicopter modes of motion. Several accidents have been recorded and show that, under some conditions, it can be very violent and can lead to the total destruction of the aircraft. [START_REF] Coleman | Theory of self-excited mechanical oscillations of helicopter rotors with hinged blades[END_REF] performed some of the earliest research on the phenomenon, and laid the foundation for all of the work that was to follow. Between those, [START_REF] Donham | Ground and air resonance characteristics of a soft inplane rigid rotor system[END_REF], and [START_REF] Lytwyn | Airborne and ground resonance of hingeless rotors[END_REF] studied this phenomenon taking into account also the air resonance effect. Major contributions in understanding the ground resonance phenomenon in hingeless and bearingless rotors were done by Army researchers, as for example, [START_REF] Dawson | An experimental investigation of the stability of a bearingless model rotor in hover[END_REF] and [START_REF] Hodges | An aeromechanical stability analysis for bearingless rotor helicopters[END_REF]. Recently, [START_REF] Kunz | Nonlinear Analysis of helicopter ground resonance[END_REF] analysed the influence of nonlinear springs and dampers (elastomeric elements) on predicting the rotor's instability zone. [START_REF] Byers | Rotor blade with radial orber[END_REF] explored the passive control of the problem.

However, in order to establish the criteria to avoid the unstable oscillations, the equations of motion were linearized and then simplified by eliminating the periodic and parametric termes. This last process, applicable only for isotropic systems, is identified as the Coleman Transformation and more generally called the multi blade coordinate transformation.

A recent study given by [START_REF] Skjoldan | On the similarity of the Coleman and Lyapunov-Floquet transformations for modal analysis of bladed rotor structures[END_REF] kept the periodic and parametric terms in the equations of motion of a wind turbine. The instability zone predicted by using the Floquet's theory over a certain range of rotation speeds is similar to that predicted by Coleman and Feingold, for a isotropic rotor.

With the objective of further understanding the ground resonance phenomenon, an analysis of the nonlinear equations of motion is envisioned. This paper sets up the nonlinear equations of motion by considering all terms. Multiple scales method is used to treat the equations of motion.

MECHANICAL MODEL

Similar to that proposed by Coleman and Feingold, the present mechanical model is developed to characterize the dynamical behavior of a helicopter with a hinged rotor. In other words, it consists of figuring out the relation between the longitudinal -x(t) -and lateral displacement -y(t) -of the fuselage, and the k th blade lag anglek ϕ (t) -in terms of the rotor speed Ω and time t. Figure 1 illustrates a general schema of the system.

The fuselage is considered to be a rigid body with its center of gravity at point O. At the initial time, the origin of an inertial coordinate system (X 0 ,Y 0 ,Z 0 ) is coincident at this point. The body is connected to springs which represent the flexibility of the landing skid.

The rotor head system consists of an assembly of one rigid rotor hub with Nb blades. Each blade is represented by a concentrated mass located at a distance b from the lag articulation (point B) and, on each articulation, a torsional spring is present. The origin of a mobile coordinate system (x, y, z), parallel to the inertial one, is located at the geometric center of the rotor hub (point A). Both, body and rotor head, are joined by a rigid shaft and the aerodynamical forces on the blades are neglected.

Equation of motion

The equations of motion are obtained by applying Lagrange's equation on the kinetic and potential energy expressions of the system (body and rotor head). To reach this goal, some conditions have been considered, which are:

• The fuselage has mass m f and the spring constants connected to it are K f-X and K f-Y through x and y directions, respectively; • The rotor is composed of Nb = 4 blades and each blade k has an azimuth angle of ( )

2 1 k k Nb θ π = -
with the x -axis; • Each blade has the same mass m b and moment of inertia around the z -axis is to Iz b ;

• The angular spring constant for each k blade is K b-k ; • The k blade position projected in the inertial coordinate system is: cos( ) cos( (t)) (t) x a t b t x b_k k k k Ω θ Ω θ ϕ = + + + + + (1) sin( ) sin( (t)) (t) y a t b t y b_k k k k Ω θ Ω θ ϕ = + + + + + (2)
where a is the hinge offset.

Considering the general Laplace variable q as

{ } { } { } 1..6 1 2 3 4 1..6 (t) (t) (t) (t) (t) (t) q u T T x(t) y(t) ϕ ϕ ϕ ϕ = = (3)
where u is a vector made up of the degree of freedom variables of the system. Equations (4 -6), are nonlinear equations in the function of u. Equations ( 4) and ( 5) correspond to the motion equation in the x and y directions, respectively. Equation ( 6) represents the dynamical equations of blades 1 to 4, changing n from 3 to 6 in that order.

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 2 2 1 2 2 2 2 2 d (t) (t) dt d d cos sin (t) (t) sin cos (t) (t) dt dt d d 2 cos cos (t) (t) 2 sin sin (t) (t) dt dt sin si 1 1 m n 2 n n m n 2 n n m n 2 n n m n 2 n n m n 2 u u r b Ω t u u r b Ω t u u Ω r b Ω t u u Ω r b Ω t u u r b Ω t ω θ θ θ θ θ - - - - -   +           - + - +                 - + + +     +     + + ( ) ( ) ( ) ( ) ( ) ( ) ( ) 2 2 2 2 0 d d n (t) (t) cos cos (t) (t) dt dt cos cos (t) sin sin (t) 6 n 3 n n m n 2 n n m n 2 n m n 2 n u u r b Ω t u u Ω r b Ω t u Ω r b Ω t u θ θ θ = - - -             =           - +             - + + +     ∑ (4) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 2 2 2 2 2 2 2 2 d (t) (t) dt d d cos cos (t) (t) sin sin (t) (t) dt dt d d 2 cos sin (t) (t) 2 sin cos (t) (t) dt dt sin cos 2 2 m n 2 n n m n 2 n n n 2 n n m n 2 n n m n 2 u u r b Ω t u u r b Ω t u u Ω α b Ω t u u Ω r b Ω t u u r b Ω t ω θ θ θ θ θ - - - - -   +           + + - +                 - + - +     +     - + ( ) ( ) ( ) ( ) ( ) ( ) ( ) 3 2 2 2 2 0 d d (t) (t) cos sin (t) (t) dt dt cos sin (t) sin cos (t) 6 n n n m n 2 n n m n 2 n m n 2 n u u r b Ω t u u Ω r b Ω t u Ω r b Ω t u θ θ θ = - - -             =           - +             - + - +     ∑ (5) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 2 2 2 1 2 2 2 2 1 2 2 2 2 2 2 d d (t) (t) sin (t) sin cos (t) (t) dt dt d d cos sin (t) (t) cos (t) (t) dt dt d sin sin (t) (t) 0 dt 2 n n n b n b n 2 n b n 2 n b n 2 n b n 2 n u u r a Ω u r Ω t u u r Ω t u u r cos Ω t u u r Ω t u u ; ω θ θ θ θ - - - -     + + - +                 - + + +               - + =       3 6 n ... = (6) where _ _ _1..4 2 2 2 1 2 3..6 2 1 , , , , f X f Y b b m b f b f b f b b b K K K m r r b m Nb m m Nb m m Nb m b m Iz ω ω ω = = = = = + + + +
The terms 1 ω and 2 ω are the resonance frequency of the fuselage in the x and y directions, respectively. Moreover, 3..6 ω are the resonance frequencies of blades 1 to 4. The terms r m and r b are smaller than one.

FLOQUET'S METHOD

In this section a simplified mathematical model is expanded. A Taylor's series of first order is developed for the terms sin(u n ) and cos(u n ) in the equations (4-6) and, furthermore, the nonlinear terms are neglected. Consequently, a new dynamical system is reached, and is represented as below.

ext M u G u K u F + + = ɺɺ ɺ ( 7
)
where M is the inertia matrix, G is the damping matrix, K is the stiffness matrix, and ext F is the external vector force.

Due to the fact that parametric terms exist, M becomes a symmetric and non-diagonal matrix, G and K are non- symmetric and non-diagonal matrices. Also, because of the rotor head's symmetry, ext F is equal to zero. Then, Floquet's theory is used to solve differential equations with periodic coefficients. Considering S(t), the statespace matrix with period T of Eq. ( 7), and its state-space representation as

[ ] 0 0 0 ( ) ( ) ( ), ( ) , ( ) ( ) ( ) v S v v v v u u T t t t t t t t t t = > = = ɺ ɺ (8)
By Floquet's Theory, a transition matrix Φ , which relates ( ) v t and 0 ( ) v t , is defined as

( ) ( ) ( ) 0 0 0 , , Q Φ P t t t t t t e - = (9)
The monodromy matrix R and the constant matrix Q are defined, then:

( ) 0 0 , R Φ t T t = + and
( )

1 log Q R T = (10) 1 1 0 1 0 1 0 0 0 0 1 0 1 0 0 1 0 0 2 2 0 1 0 1 1 1 0 1 3 2 2 0 0 1 0 0 1 0 1 2 2 2 2 0 0 1 2 3 2 0 0 1 0 D , , 1 D , D , 2 , 1 1 D , 2 2 1 D , D 2 R R R R R I Ω T T I Ω T T n n n R n n I T T I T T n n n R I T T n n n R n u T T u T T I b u T T e I b u T T e u T T I b u T T e I b e I b u T T e I b u T σ σ Ω σ Ω σ Ω σ ω α α Ω α α Ω α α Ω + + - - + + - - + - - + =   +       + +     - + + ( ) ( ) ( ) ( ) ( ) ( ) [ ] 0 1 0 1 0 6 3 0 1 2 3 0 0 1 , . 1 D , 2 R R I T T n I T T n n T e c c I b u T T e Ω σ Ω σ α + = + -                       +         +         ∑ (14) D , , 1 D , D , 2 , 1 1 D , 2 2 1 D , D 2 R R R R R I T T I T T n n n R n n I T T I T T n n n R I T T n n n R n u T T u T T I b u T T e I b u T T e u T T I b u T T e I b e I b u T T e I b u Ω σ Ω σ Ω σ Ω σ Ω σ ω α α Ω α α Ω α α Ω + + - - + + - - + - - + =   - +       + +     - - + ( ) ( ) ( ) ( ) ( ) ( ) [
]

0 1 0 1 0 6 3 0 1 2 2 0 0 1 , . 1 D , 2 R R I T T n I T T n n T T e c c I b u T T e Ω σ Ω σ α + = + -                       +         -         ∑ (15) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) [
]

1 1 0 1 0 1 0 0 0 0 1 0 1 0 0 2 2 0 0 1 0 1 3 2 2 2 0 1 0 1 0 2 0 1 0 1 2 3 2 2 0 1 0 1 0 2 0 1 D , , 1 1 D , D , , 2 2 1 1 D , D , . 3..6 2 2 R R R R n n n I T T I T T n n n I T T I T T n n u T T u T T I u T T e I u T T e u T T I u T T e I u T T e c c n Ω σ Ω σ Ω σ Ω σ ω β β β β + + - - + + - - + =   - +     - - + = (16)
where [c.c] corresponds to the complex conjugate of the previous terms.

As observed, the right-hand side of the above equations is all dependent on the steady-state responses 0 n u .

Secular terms -possible instable regions

The developing of Eq. (14-16), after Eq. ( 13) have been substituted into them, lead to combinations of n ω and R Ω multiplied by T 0 as exponent to each equation of motion. Therefore, depending on the applied specific rotor speed R Ω , secular terms may exist. This fact implies that the response(s) 1 n u grows without limit as t increases and, consequently, it does not provide a small correction to 0 n u .

Based on that, a set of possible rotor speed values of R Ω for each equation is listed and it indicates instability zones of the system. These zones are well defined through the analysis of the amplitude responses calculated around a critical value for R Ω by shifting the detuning frequency parameter.

Body's Equation -Eq. ( 14):

1 1 1 1 0 2 3 3..6 n n n n ω ω ω ω ω ω ω   ± ± ± =   (17)
Body's Equation -Eq. ( 15):

2 2 2 2 0 2 3 3..6 n n n n ω ω ω ω ω ω ω   ± ± ± =   (18)
Blade's Equations -Eq. ( 16):

1,2 1,2 1,2 2 3..6 n n n ω ω ω ω ω   - - =   (19)

RESULTS

The goal of this section is to predict and collect the speed rotor value R Ω for which the system becomes instable from both mathematical models developed previously. Later, the results are compared. The numerical values used for the system's inputs are defined in Tab. 1.

f = 2902.9 [kg] m b = 31.9 [kg] a = 0.2 [m] b = 2.5 [m] K f-X = K f-Y = 1.077e3 [kN/m] K b-1..4 = 40.716 [kNm/rad] I zb = 259 [kg m 2 ]
Applying the numerical values, the resonance frequency, in Hz, of the fuselage 1 ω and 2 ω are equal to 3.00 and, for the blades 3 ω , 4 ω , 5 ω and 6 ω are equal to 1.50. By varying the rotor speed from 0 to 7 Hz, Fig. 2 illustrates the instability zone for 4.358< Ω <5.187, in Hz, by using Floquet's Method. It is defined when the real part of the characteristic multipliers are positive and nonzero.

CONCLUSION

As for the comprehension and analysis, this paper envisioned to verify the influence of the nonlinear terms in the equations of motion on the prediction of the rotor speed which the ground resonance phenomenon may occur. Two different models were developed and compared: one keeps the linear, parametric and periodic terms and another considers all terms, including the nonlinear ones.

The numerical results obtained in the first model, by applying the Floquet's Theory, show a well defined instability zone for the rotor speed between 4.358 and 5.187 Hz. Time response analysis for Ω =4.77 Hz evidences the divergent behaviour of the fuselage displacement.

However, the analysis of the secular terms in the nonlinear model, by using multiple scales method, verifies many others possibilities of instability zones. They are defined as a combination of the resonance frequency of the fuselage and the blades', by considering the linear system of motion as uncoupled. It evidences the importance of considering the nonlinear terms into the model in order to predict the instability zones.

The perspectives for this work consist in studying the amplitude responses in function of the detuning frequency parameter. It is envisioned too the analysis of the sensibility of each unstable region and the amplitude response in terms of internal detuning frequency parameters.

A passive control and an experimental validation of the results, through an experimental set up, are envisioned.

Figure 1 .

 1 Figure 1. General schema of the mechanical model

Figure 2 .

 2 Figure 2. Campbell Diagram by Floquet Method.However, by using Multiple Scales Method, the rotor speed values in Hz, at which unstable oscillations for the fuselage and the blades may occurs, are ordered as follow:

Table 1 .

 1 Numerical values of the fuselage and rotor head inputs

	Fuselage	Rotor
	m

The dynamical system Eq. ( 7) is stable if the eigenvalues (µ) of Q are negative or, similarly, if the norm of the eigenvalues of R is less than one.

MULTIPLE SCALES METHOD

The perturbation method used to obtain the solution for these nonlinear equations is the multiple scales method [START_REF] Nayfeh | Nonlinear Oscillations[END_REF]. It consists of assuming that in any function dependent of time, the independent variable t is a function of multiple scales of the time at the first order. Then, the perturbation substitutions created by introducing a bookkeeping parameter ε are:

where, Ω R is the nominal rotor speed, σ is the detuning frequency parameter of the rotation speed and the terms T 0 and T 1 correspondent to t and ε t, respectively.

The terms sin(u n ) and cos(u n ) are expanded by a Taylor's series of first order in u and, after that, the perturbation substitutions are done on the equations of motion. Consequently, several sets of equations are obtained by grouping them by power of ε.

Order ε 0 equations

The set of order ε 0 equations of the system is collected. They represent the steady-state response of the uncoupled rotor and fuselage system.

where D 0 is the derivative with respect to T 0 .

The solutions of these equations are trivial and are in the form of

, since the steady-state responses are zero. The term A n (T 1 ) is complex and defined as

, where C n and B n are the amplitude and phase angle of each response, respectively.

Order ε 1 equations

The following, Eq. (14-16), is the set equations of order ε 1 . (

) ( )


Muhammad Ilyas
email: muhammad.ilyas@isae.fr

Christine Espinosa
email: christine.espinosa@isae.fr

Frédéric Lachaud
email: frederic.lachaud@isae.fr

Laurent Michel
email: laurent.michel@isae.fr

Michel Salaün
email: michel.salaun@isae.fr

Modelling Aeronautical Composite Laminates Behaviour under Impact using a Saturation Damage and Delamination Continuous Material Model

Keywords: unidirectional composite, damage mechanics, impact, strain rate behavior, numerical simulations

We show that the behavior of T700/M21s and T800/M21s composite panels are affected by the influence of strain rates together with local shear and crush punch or global flexural strengths of the structure. A deterministic continuous composite material model has been developed as a LS-DYNA user defined material model for unidirectional composites on the basis of the Matzenmiller model widely used for woven composites. Initiation and evolution up to saturation and fracture are implemented for various and coupled damage mechanisms including delamination. Quasi-static and dynamic characterization tests laminates have been carried out on balanced angle ply [±θ] and used for calibration of numerical values. Impact induced damage from experiment's measures and numerical predictions are compared for T800/M21S aeronautical samples impacted at 15J.

Introduction

Where usual rules and design admissible criteria are well known for metallic materials in airplane certification procedures, the global residual strength of composites in case of damage is not well known. Catastrophic cases are the ruin of a structure due to uncontrolled propagation of inner damages. In order to better understand what happens during an impact, and to aid the design of composite structures, it is our aim to experimentally investigate the behavior of T700/M21s and T800/M21s composite panels and to propose a behavior model that could be used to predict the amount of damage after impacts for a wide range of energies or velocities. We present here some part of the work we have done on experimental characterization, theoretical behavior modeling and development. To investigate strain rate effects, 2 real experiments are compared to simulations predictions: one at a low mass and medium velocity, one at a medium mass and a low velocity, both experiment giving the same impact energy of about 15J.

Experimental characterization T700/T800

Quasi-Static damage saturation. 150mmx20mm coupons of T700/M21S and T800/M21S [±45°] 2s were fabricated and tested under quasi-static axial length traction at 2 mm/min that is about 2x10 -4 s -1 (Fig. 1a). The general linear elastic behavior of the ply is defined up to 1.2% deformation and a related yield stress of 50 MPa for both materials. Following destructive examination of the tested samples, and considering the different dimension of the T700 and the T800 fiber diameters, it is suggested that the curved part of each curve between 50 MPa and 60 MPa is related to a coupling effect between the inner ply tension and the interface behavior. Above 60 MPa, the T800 based lamina is governed by the interface behavior whereas in the T700 lamina, there is still a competition between the inner ply and the interface irreversible behaviors. It is concluded that the inner plies of both materials can not suffer perpendicular matrix damage higher than a saturation value.

Dynamic damage and rupture. Dynamic compression tests have been conducted using a Split Hopkinson Pressure Bar test system at ISAE. Bars are stainless steel 2m long, and 20 mm diameter rods. Different angles were tested at different strain rates in order to identify in one hand the effect of rates on the global behavior, and on a second hand, the contribution of the pure to coupled loadings between inner ply fibers and matrix, and between plies interfaces. Square 20mmx20mm samples of [±θ] 3s were tested at angles of 15°, 30°, 45° (Fig. 1b), 60°, 75°.

The dynamic behaviour of [±45°] 3s is a 4 parts stress-strain curve: 1st linear (up to 1.2% strain), 2nd a curved part passing a maximum stress (at about 2% strain), 3rd a rather linear decreasing part with a strain rate dependent slope (from 3% strain), 4th rupture at a failure strain increasing with increasing the strain rate (always ≤ 10%). Other angles do not exhibit the third part of the curves. Apparent rigidities are all about 20GPa (±2 GPa depending on strain rate). As for QS testing, using destructive analysis and previous studies observations [START_REF] Ilyas | Dynamic delamination using cohesive finite elements[END_REF], it is suggested that the2nd part is the coupling of inner and inter ply irreversible behavior (onset of damage) and the 3rd part is related to interfaces and parallel matrix cracking only present at [±45°] interfaces. It is suggested that the fourth part of the curves is related to damage saturation in the ply of T800/M21s samples. Since the deformation yield is always the same, whatever the strain rate, it is supposed that the saturation damage is also the same. The growth of damage and the consequent maximum stress level in the 2nd part of the curves are strain rate dependent above a reference level of 750s -1 . Strength versus strain rate has been interpolated (Fig. 1c) using Eq. 1.

max 0 . 1 log ref S S C ε ε     = +            
with c=4,7 and S 0 =120 MPa.

(1) Material model Starting model. The basis material model is the non linear elastic anisotropic continuum damage mechanics one developed by Matzenmiller &al. [START_REF] Matzenmiller | A constitutive model for anisotropic damage in fiber composites[END_REF]. As done "classically", damage is modeled through its effect on the elastic rigidity loss in further loading or unloading until the damage reaches the value of 1 meaning rupture. The model has been generalized by Xiao &al. [START_REF] Xiao | Progressive damage and delamination in plain weave S-2 glass/SC-15 composites under quasi-static punch-shear loading[END_REF] for 3D composites and essentially used for woven composites with success for high energy impacts. The model distinguishes 6 damage variables {ω i }, i=1,6 and 5 ruin modes {r j }, j=1,5. The onset and growth of the damage variables can be expressed as a linear combination of the thresholds evolution or using maximum values with the following expressions where q ij is a coupling tensor between the individual damage modes ω i and the thresholds r j and m is a strain softening parameter that can be adjusted separately for each mode (Eq. 7).

New criteria for UD. New failure criteria are proposed hereafter. Note that when ω 2 reaches saturation then ω 4 and ω 5 are saturated too, leading to a failed E 22 G 12 and G 23 . The summation method is chosen with a coupling tensor given in Eq. 7.

2 2 2 11 2 12 13 1 1 2 ( , , ) 0 T fs f r r X S σ σ σ σ ω     + = + -=           Tension/Shear . ( 2
) 2 11 22 33 2 2 2 2 ( , , ) 0 2 C f r r X σ σ σ σ ω   - + -- = -=       Compression. (3) a) b) c) 2 11 22 33 2 3 3 ( , , ) 0 3 C f r r Z σ σ σ σ ω  -- -  = -=     Crush. ( 4
C f r r X Y S S σ σ σ σ σ ω σ ϕ σ ϕ        -  = + + + -=         + - + -         Perp. matrix. ( 5
T f r r Z S S σ σ σ σ ω σ ϕ σ ϕ       = + + + -=       + - + -       Parallel matrix . (6)
[ ]

1 1 1 0 0 0 0 1 1 0 0 0 1 0 1 1 1 1 1 1 0 0 1 1 1 1 1 1 0 1 q         =             5 1 ( , , ) . ( , , ) i ij j j q ω σ ω ε φ σ ω ε = = ∑ ( ) 1

Application to impact on T800

Impact tests and the corresponding numerical simulations have been performed on simply supported laminates [-45/+45/0/90/0/0/-45/+45/0]s 150mmx100mm T800/M21S panels at about 15J, with two different couples of mass and velocities for the projectile: drop tower (2.368kg, 3,47m/s), canon (17g, 40,82m/s). Results under consideration are the global behavior of the impacted plates given by the force versus time and external craters, and strain rate effects on damage predictions. Numerical model. The numerical model is composed of 8 nodes solid elements, one element per ply in the thickness. The support is rigid as is the projectile and contact are introduced with the plate.

Global behavior. Models presented on Fig. 2a are a drop test curve (Test; 7400N/3.17ms) and numerical curves using tensile E 11 Young modulus with elastic plies (simu_elastic), or with our model and a constant quasi static strength of 50 MPa (sigr_qst), or with our model and the strain rate dependent strength (sigr_dyn). Fig. 2b shows the drop test curve and numerical curves using the flexural modulus Ef with a cohesive interface (int_coh_flex) previous model [6], or with our model and the strain rate dependent strength (sigr_dyn_flex; 8300N/3.05ms). Errors are 12%/4% in force/time. On Fig. 2a and 2b, the first linear elastic part is well recovered by all the models (up to 0,7 ms), even though it is obvious that the global rigidity is better retrieved with the flex-modulus. Between t=0,7ms and t= 1,5 ms, the global frequency of vibration is correctly recovered with the strain rate dependent models, especially our dyn_flex model even if slightly delayed in time. After a small plateau 1,5-1,7 ms our model and the experimental curve decrease identically down to 3 ms.

Numerical indentation (relative displacement of the upper and bottom faces) at max force time (Fig. 2c) show a promising numerical prediction with our model also compared with a cohesive model we previously developed [6]. Typical values are 0.16 mm for test, 0.22 mm for our dyn_flex model, 0.32 mm for the cohesive model. Note that the indent has been measured 48h after the test while the numerical depth does not take into account any dynamic relaxation. Strain rate effects. The spatial localization of elements failed by saturation in the dyn_flex model (Fig. 3b) gives a prediction of the delaminated area (0°_26,5mm, 90°_24mm) similar to what can be observed on C-Scan after a 15J Drop Tower test (0°_25,5mm, 90°_29mm). Our previous cohesive model is more diffusive (0°_33,5mm, 90°_39mm), as is also the qst model.

Conclusion

In this study a composite damage material model has been presented that is able to cover a wide range of strain rate loadings, while calibrated up to now with strain rates lower or equal to 10 3 s -1 .

On the experiments analysis, we have proposed and presented our model that consists in proposing new damage criteria and coupling for UD aeronautical composites in the basis model of Matzenmiller &al. and Xiao & al.. A saturation effect was observed for both T700 and T800/M21s samples and taken into account in the dynamic model. The experimentally determined data set has been used to predict impact induced damage with an acceptable ability of the model to represent both the global flexural and vibrating behavior and the local outer indent and inner damages suffered by the composite plates. Further work consist in doing tests at higher strain rates, realizing numerical design simulations to better identify the limitation of our model, improve the cohesive model by adding rate effects, and better use it to establish links between the outer indent and the inner damages at various impact levels.

Figure 1 .

 1 Figure 1. a) Quasi-static tension of T700/M21S and T800/M21S [±45°] 2s ; b) Dynamic compression of T800/M21S [±45°] 3s ; c) Strain rate effects on T800/M21S [±45°] 2s lamina strength.

  Calibration of material properties. Comparisons between QS and dynamic tests and simulations have been done for T800/M21S determine the following values, and the flexural Young modulus following previous work recommendations[START_REF] Ilyas | Dynamic delamination using cohesive finite elements[END_REF]. E 11 = 165 GPa E 22 = 7.64 GPa E 33 = 7.64 GPa E f = 112 GPa X T = 2.2 GPa X C = 1.2 GPa ν 21 = 0.016162 ν 31 = 0.016162 ν 32 = 0.4 m i =10 Y T = 45 MPa Y C =280MPa G 12 = 5.61 GPa G 23 = 2.75 GPa G 13 = 5.61 GPa Z T = 45 MPa Z C = 0.7 GPa S ffc = 0.5 GPa S 12 = 0.05 GPa S 23 = 0.05 GPa S 31 = 0.05 GPa S fs = 1.5 GPa ω max =0.87 ϕ = 10°

Figures 2 .

 2 Figures 2. a) and b) Force [N] versus time [ms] history curves; c) Indentations [mm]. a) b) c)

Figure 3 .

 3 Figure 3. a) Delaminated cohesive or b) saturated elements model, and c) DT 15J C-Scan.

Table 1 .

 1 Data set for our model.

Acknowledgement

The authors are grateful of Higher Education Commission of Pakistan for partial funding. Special thanks are also extended to IMPETUS Afea France and all the students and technical staff of ISAE for their valuable input towards numerical and experimental aspects, respectively.


Rony T Chung
email: ronychung.br@gmail.com

Flávio Luiz

S Bussamra

Pierre Selva

Joseph Morlier

NUMERICAL AND MECHANICAL MODELING OF THE INNER EAR VESTIBULAR APPARATUS

Keywords: Inner ear, Vestibular System, Finite Elements Modeling, Computational Fluid Dynamics, Bioengineering

As you read this text, tilt your head left and right. Note then that the text stands still. Now draw a clear line onto the ground, close your eyes and walk over it. You most likely succeeded. The aforementioned facts or abilities are remarkable skills given by a small inner ear organ set, the vestibular system. Usually operating in the background, this structure infrequently comes to the attention of our conscious mind. In this work, three finite element models (in CATIA V5R18 and COMSOL Multiphysics 3.5a environments) are developed to provide a mechanical understanding of this formidable nature's gift. In these numerical models, a fluid structure interaction modeling is presented to show the system behavior during a horizontal head spin. Finally, a physical augmented scale model is proposed, including the choice of materials and a similitude study.

INTRODUCTION

Sense of balance or equilibrioception can be defined as one's physiological capacity of perception of the external environment, granting movements, displacements and gestures with reasonable efficacy. The simple task of standing up right or walking on a straight line is only possible for the spatial equilibrium sense exists. From this perspective, a clear permanent mental representation of one's body in the external ambiance, conscious or unconscious, along with a prompt repertoire of motor reactions adapted to specific situations, becomes critical. The balance notion, for its turn, is a plurimodal sensorial nucleus [START_REF] Roman | Encyclopédie Médico-Chirurgicale[END_REF] whose origins are:

-Visual, providing plastic cues for localization in space and relative velocities -Vestibular, informing the linear and angular accelerations and the head position and inclination in relation to the gravitational axis -Proprioceptives, describing the relative disposition of different parts of the body, their tension and accelerations. This function is remarkably effectuated by the muscles acting in the upright position and profound neck muscles output nerves.

The vestibular system becomes critical in low visibility situations, such as night or cloud non-instrumental aircraft flight. Furthermore, it provides fair interpretation of high frequency movements and a vital capacity to human quotidian life, the Vestibular-Ocular-Reflex (VOR), which refers to a reflex eye movement that allows humans to fixate the gaze at one desired object while in motion. The image is stabilized on the retina through a series of rapid eye movements in the opposite direction of the body motion. For slight head movements are present all the time, the VOR is crucial. For instance, patients with impaired VOR find it difficult to read due to normal involuntary head tremors. Pilot illusions are frequent issues in aeronautics and many of them are related to the vestibular mechanism.

For the above mentioned reasons, it may be fair enough to classify the vestibular system as acutely important for good quality of life.

ANATOMY AND PHISIOLOGY OF THE SEMICIRCULAR CANALS

Most animals posses a pair of inner ears disposed symmetrically in relation to the Sagittal plane (vertical plane between the eyes). These structures are relatively small, of the magnitude of millimeters. The membranous labyrinth is filled up with a Newtonian incompressible potassium rich fluid, the endolymph, crucial for the orientation process [START_REF] Selva | A Matlab/Simulink model of the inner ear angular accelerometers sensors[END_REF]. The labyrinth of the inner year is constituted of two major branches, the cochlea, a peripheral part of the auditory system, and the vestibular system, a critical equilibrioception structure. For its part, the latter is subdivided in two functional sensors: the semicircular channels and the otolith, capable of measuring angular and linear rotations, respectively.

As shown in Figure 1, there are three orthogonal semicircular canals in each inner ear: the horizontal also known as lateral or external, the posterior, also named superior, and the anterior one. In the regular head position, the horizontal canal rests horizontally in space, whilst the others remain vertical.

Specifically, only the vertical canals present a common branch linking them to the vestibule, though all canals end up in the vestibule by a dilated orifice, the ampulla, and a simple orifice. The ampulla is a swelling which contains its own hair cell sensory epithelium, the crista. The bundles of crista hair cells protrude into the ampulla interior involved by a gelatinous acellular membrane, the cupula, as represented in Figure 2. This diaphragm plays a major role as a vestibular angular acceleration detector. [START_REF] Roman | Encyclopédie Médico-Chirurgicale[END_REF] Figure 2 : Cupula Membrane [START_REF] Corwin | The Vestibular System[END_REF] When there is a rotation in the horizontal plane, the horizontal semicircular canal follows. The endolymph within, however, lags behind due to its inertia, originating a relative movement between the fluid and the canal inner walls. Since the cupula is coupled to the canal, it is also in relative movement with the endolymph, which, for its part, exerts a force on the membrane, finally deflecting it. Grosso modo, each semicircular canal is predominantly excited by a pure rotation in its correspondent plan by the intersection point of the three central axes of the canals. However, during a head rotation, usually all three canals are excited because the inner year is not located exactly in the center of the head and the canals are not perfectly oriented and orthogonal. The Central Nervous System (CNS) is charged of interpreting and calibrating this experience with the perceived reality.

For both left and right canals, a deflection towards the rostral end, as shown in Figure 3, depolarizes the hair cells. Thus, during a counter clockwise rotation, a depolarization of the hair cells occurs on the left canal and a hyperpolarization takes place on the right one. Since the canals are mirror-opposite, they always operate antagonistically as push-pull pairs, or, in simple terms, when one is hyperpolarized, the other one is depolarized. Figure 3 illustrates the scheme. These signals are transmitted to the neurons that innervate the canals and are then delivered to the CNS for post-processing. Figure 4 illustrates the deformation and transduction process. A similar line can be applied to explain the posterior and anterior semicircular canals transducers. It is important to remark, though, that once the acceleration ceases, the endolymph tends to spin along with the canal walls due to its viscosity and the non slipping condition. Moreover, the cupula gradually returns to its unstrained state due to its elastic properties, describing more precisely an exponential curve [START_REF] Selva | A Matlab/Simulink model of the inner ear angular accelerometers sensors[END_REF].

NUMERICAL MODEL

Fluid structure interaction

The problem must be resolved with a multiphysical algorithm, since pressure and viscous forces exerted by the fluid on the structure cause it to deform, and this deformation may significantly alter the original flow geometry. The arbitrary Lagrangian-Eulerian (ALE) method handles the geometry deformation and fluid changing boundaries with a moving grid. In this approach, the so called mesh smoothing takes place, with the computation of new mesh coordinates based on the movement of structure boundaries. The Navier-Stokes equations that solve the flow are employed in this set of moving coordinates. For its turn, the structural mechanics part of the model does not require the ALE method and is thus solved in a fixed coordinate system as usual. Nonetheless, its strains serve as feedback information for computing the deformed coordinates with ALE.

The fluid flow is derived from the incompressible Navier-Stokes equations for the velocity field 𝒖 = (𝑢, 𝑣, 𝑤). Being 𝜌 the fluid density, η the dynamic viscosity, 𝑝 the pressure and 𝒖 𝒎 the coordinate system velocity,

𝜌 𝜕𝒖 𝜕𝑡 -∇ . [-𝑝𝑰 + η(∇𝐮 + (∇𝐮) T )] + 𝜌�(𝒖 -𝒖 𝑚 ). ∇�𝒖 = 𝑭 (1) -∇. u = 0 (2)
where 𝑰 denotes the unit diagonal matrix and 𝑭 the volume force on the fluid. These equations are solved referring to boundary settings such as the no-slip condition (𝑢 = 𝑣 = 𝑤 = 0) in the contact with a fixed wall, a reference pressure and velocities equal to deformation rates in the deformable solid and fluid interface. The structural strains are calculated through an elastic formulation and a nonlinear geometry formulation for the larger deformations. It is then necessary to solve

𝜌 𝑠 𝜕 2 𝒖 𝑠 𝜕𝑡 2 = ∇. 𝝈 𝒔 (3)
where 𝒖 𝑠 = (𝑢 𝑠 , 𝑣 𝑠 , 𝑤 𝑠 ) is the displacement vector, 𝜌 𝑠 is the material density and 𝝈 𝒔 is the stress tensor. Considering an elastic isotropic material, the generalized Hookes law gives

𝝈 𝒔 = 2𝜅𝝐 + 𝜆𝑡𝑟(𝝐)𝑰 (4)
The coefficients 𝜇 and 𝜆 are the Lame's coefficients and can be related to the Young Modulus 𝐸 and the Poisson's ratio 𝜈 through

𝜅 = 𝐸 2(1 + 𝜈) (5) 𝜆 = 𝜈𝐸 (1 + 𝜈)(1 -2𝜈) (6)
As boundary conditions, the solid on the fluid solid interface is subject to a fluid pressure 𝑭 𝑻 . Also, as aforementioned, the velocity field on the fluid solid interface is obtained with the deformation rate

𝑭 𝑇 = -𝒏. �-𝑝𝑰 + η(∇𝐮 + (∇𝐮) T )� (7) 𝒖 = 𝑑𝒖 𝒔 𝑑𝑡 (8)
where 𝒏 denotes the boundary normal vector.The moving mesh can be computed with algorithms such as the Winslow smoothing. In this case, the solved equation system is

𝜕 2 𝑋 𝜕𝑥 2 + 𝜕 2 𝑋 𝜕𝑦 2 = 0 (9) 𝜕 2 𝑌 𝜕𝑥 2 + 𝜕 2 𝑌 𝜕𝑦 2 = 0 (10)
where 𝑥 and 𝑦 are the spatial coordinates in the frame to be the defined and 𝑋 and 𝑌 are are the reference coordinates of the reference frame. At the boundaries of the solid obstacle, the mesh displacement is the same as the structural deformation. At the exterior boundaries of the flow domain, this deformation is null in all directions.

In COMSOL Multiphysics (COMSOL, 2008), which is employed in the numerical analysis in this work, the resolution is supported by three application modes which interact to solve the multiphysic problem: the incompressible Navier-Stokes for the fluid flow, the Plane Strain for the structure mechanical segment and the Moving Mesh (ALE) for computing the mesh deformation.

The computational model has been developed in COMSOL Multiphysics with the following simplifications: the model comprises the utricule and the semicircular canals with their ampullas; the saccule and other labyrinth components influence is assumed neglectful; the geometry of the semicircular canals are deemed similar and with constant elliptic section and, for geometry simplicity, the cupula is modeled in such a manner that it entirely blocks the ampulla cross-section. The model geometry was created with computer aided design (CAD) in CATIA.

Geometry

In this section, the CAD design took part in CATIA V5R18. In general terms, three models have been designed: the horizontal (single), horizontal/posterior (double) and horizontal/posterior/anterior (triple) semicircular canals models, all of which containing the utricule and ampulla cupula segment.

Firstly, the geometry of the single horizontal vestibular canal -along with the ampulla cupula -has been designed according to the dimensions shown in Figure 5. In simplified terms, the utricule cupula part has been constructed with the aid of cross sections, connected with the Multi-sections Solid tool. For the horizontal canal, the nervure tool has been used with a semicircular guide curve and an elliptic cross section. A connection root has been extruded on the utricule surface for pressure input and future link with the other canals. Also, four defined cross sections have been delimited for the cupula modeling in COMSOL.

Figure 5: Horizontal Canal Geometry [START_REF] Curthoys | Dimensions of the horizontal semicircular duct, ampulla and utricule in the human[END_REF] Should the similarity of the three semicircular canals be taken, the mentioned structure could be rotated by 90 degrees around the x and y axes to build up the anterior and posterior canals, with slight geometrical differences. The newly formed canals are connected to the utricule by a common branch.

As the main objective was the exportation to COMSOL Multiphysics, an effort has been made to construct continuous soft surface CAD geometries which the latter could eventually recognize. For this same reason, the Boolean partial relimitation tool has been used in the anterior and posterior canals connections to the utricule and for the common branch bifurcation. The multiphisical model is obtained from the importation of the CATIA stp format file with the COMSOL CAD Import Module.

Model Study Case

The multiphysical model consists of angular velocity imposition on the canal walls and utricule endolymph, which for its part originates pressure on the cupula diaphragm deflecting it. The walls are set to be rigid, except for the cupula adjacent regions, whose mesh is set to move freely (ALE), and the cupula (Solid Subdomain), set to deform due to endolymph pressure (Fluid Subdomain). A temporal analysis is utilized and in the post processing phase, after a horizontal head input angular velocity, the fluid velocity field and the cupula deflection can be observed.

The tested input movement is a pure turn around the horizontal canal central normal axis. The Δ 𝑥 and Δ 𝑦 displacement and 𝑣 𝑥 and 𝑣 𝑦 velocities are given by Δ 𝑥 = [𝑋. cos(𝜔. 𝑡) + 𝑌. 𝑠𝑒𝑛(𝜔. 𝑡) -𝑋] * 𝑓𝑙𝑐2ℎ𝑠(𝑡 -𝑡 0 , Δ𝑡) (11) Δ 𝑦 = [𝑋. cos(𝜔. 𝑡) + 𝑌. 𝑠𝑒𝑛(𝜔. 𝑡) -𝑋] * 𝑓𝑙𝑐2ℎ𝑠(𝑡 -𝑡 0 , Δ𝑡) ( 12)

𝑣 𝑥 = 𝑑Δ 𝑥 𝑑𝑡 ( 13
)
𝑣 𝑦 = 𝑑Δ 𝑦 𝑑𝑡 ( 14
)
where 𝜔 = 0.5 𝑟𝑎𝑑/𝑠 is the steady state angular speed and the function 𝑓𝑙𝑐2ℎ𝑠 is used to soften the step input in angular velocity in the beginning of the simulation, avoiding problems in the derivatives 𝑣 𝑥 and 𝑣 𝑦 calculations. 𝑡 0 -Δ𝑡 and at 𝑡 0 + Δ𝑡 correspond, respectively, to the beginning and the end of the step input. In the single canal problem, 𝑡 0 = 0.13 𝑠 and Δ𝑡 = 0.12 𝑠. For the double and triple canal models, 𝑡 0 = 0.18 𝑠 and Δ𝑡 = 0.16 𝑠. Figure 7 and Figure 8 show the imposed angular velocities 𝜔.

Material Properties

The Cupula is modeled as a membrane (Solid Subdomain) and the Endolymph as a fluid (Fluid Subdomain). Their pertinent utilized properties are shown below in Table 1. Those properties are applied to their respective elements after the meshing process.

Boundary Conditions

As the CATIA models are rigid walls closed organs, an outside reference pressure is needed as a boundary condition. For that reason, a built in entrance with a free flow condition is imposed in all three models. The fluid no-slip condition in relation to the canal walls is modeled by making the fluid boundaries spin along with the walls at the same rate. Thus, the canal walls obey Eq.(11) to Eq.( 14). In the fluid solid interface, the endolymph moves at a velocity according to Eq.( 8) and the cupula boundaries are finally set to face the pressure field given by Eq.( 7).

Meshing

The meshing was conducted primarily with COMSOL Automatic Mesher. As the existence of small and inverted elements is hazardous for the solver convergence and calculation velocity, regions where they were present have been manually remeshed. For the

𝜔 [rad/s] time [s]
ampulla and cupula sections, a final additional refinement was implemented, not only because of the cupula membrane, but also due to the high variation in the canal cross-section and medium line, establishing special flow conditions. Tetrahedral elements have been utilized with approximately 16000, 20000 and 25000 elements for the single, double and triple canal models (see Figure 6). Later refinements were limited due to limited solving memory capacity (8 GB RAM).

Post-processing

The velocity field in the canal reference frame and the cupula deformation with time can be seen from Figure 7 to Figure 13. for each canal in the following pages. From Figure 7, it can be seen that an endolymph counterclockwise recirculation flow is established on the left side of the cupula. The Reynolds number being low due to low speed and canal dimensions, the flow is totally laminar. For 𝑡 = 0.13 𝑠, the maximum input angular velocity is imposed, and thus the maximum flow speed is reached on top left side. A consequence of this flow velocity field is a higher pressure on the right side of the cupula, which results in its leftwards displacement shown in Figure 8. The maximum strain is developed in a time lag of Δ𝑡 1 = 0.04 𝑠 at 𝑡 1 = 0.17 𝑠.

For the double canal model (Figure 9 and Figure 10), the obtained results do not appear to fit well the real system and the achieved simulation time was short. The velocity at the free flow boundary is highest among all, edging the single canal model maximum flow speed by ten times. The originated swirl is then clockwise, probably also due to the introduction of the posterior canal, lowering locally the pressure. The flow in this case exits the semicircular canal, describes a clockwise swirl in the utricule and exits by the open boundary. This model presented some singularities in the velocity field close to the cupula free adjacent regions. The horizontal canal velocity reaches a maximum at the end of simulation 𝑡 = 0.26𝑠 in the semicircular canal, and the maximum strain is developed at 𝑡 1 = 0.23 𝑠.

In the triple canal model the highest velocities and strains (Figure 11) are reached in the anterior canal due to the fact that this canal's ampulla is inclined in comparison with the horizontal plane. Thus the effect of gravity pulls down the cupula and originates a higher upper flow pressure on the membrane, deflecting it downwards . On the other hand, the flow in the posterior canal is not as intense, as one would expect from a pure horizontal rotation stimulus. For better viewing the horizontal canal, Figure 12 and Figure 13. omit the anterior canal results. It can be seen that, for 𝑡 = 0.198 𝑠, the maximum horizontal canal velocity is obtained and after Δ𝑡 1 = 0.202 𝑠, or 𝑡 1 = 2.200 𝑠, the maximum displacement is obtained. It is interesting to note that, because of gravity, the maximum velocity is now situated on the outfall of the posterior canal into the horizontal one, but is in magnitude similar to the one canal model. Also, this flow results in a final horizontal cupula in one order of magnitude higher.

SCALE MODEL

The objective of this part was constructing a physically similar augmented model of one semicircular canal and the utricule for both quantitative and qualitative testing. Price, construction simplicity and material availability have been determinant variables during this project. The project requirements are listed: size of the model large enough to be easily observed and photographed during the experiments. This restraint could be translated as a minimum scale factor of 20, considering the small dimensions of the vestibular system (around 7mm long). On the other hand, its scale should also be limited superiorly by 100; transparence of the fluid representing the endolymph, due to the visual experiments; low manufacturing complexity; reasonable price; cupula Material Young Modulus non inferior to 𝐸′ 𝑚𝑖𝑛𝑖𝑚𝑢𝑚 = 0.005 𝐺𝑃𝑎 (silicone rubber). Equally, for the tests: maximum Scale Model angular velocity limited by 𝜔′ 𝑚𝑎𝑥 = 250°/s, due to equipment restraints; minimum Simulated Angular Velocity of 𝜔 𝑚𝑖𝑛 = 10°/s.

Similitude Analysis

The fabrication of the scale model was preceded by a Similitude Study, with the use of the Buckingham Pi theorem [START_REF] Dym | Principles of Mathematical Modeling[END_REF], in order to verify the materials required properties and thus their feasibility in the project. The dimensionless numbers theory shall then be derived. Given the fundamental units 𝑀(mass), 𝐿 (length) and 𝑇 (time), the physical parameter 𝐺 𝑖 may be written as

𝐺 𝑖 = 𝐿 𝑎 1 𝑀 𝑏 1 𝑇 𝑐 1 (10)
where α 1 , α 2 and α 3 are integer coefficients. More specifically, Fluid density:

𝐺 1 = 𝜌 = 𝑀 𝐿 3
Angular velocity:

𝐺 4 = 𝜔 = 1 𝑇 Cross section radius: 𝐺 2 = 𝑟 = 𝐿 Cupula Young Modulus: 𝐺 5 = 𝐸 = 𝑀 𝐿𝑇 2
Dynamic viscosity:

𝐺 3 = 𝜇 = 𝑀 𝐿𝑇 Canal length 2𝜋𝑅: 𝐺 6 = 𝑙 = 𝐿 Table 2: Physical parameters dimensional analysis
The fluid velocity is function of

𝑦̇= 𝑓(𝜌, 𝑟, 𝜇, 𝐸, 𝜔, 𝑙) = 𝐿 𝑥 𝑀 𝑦 𝑇 𝑧 = 𝐿 𝑇 (11) 𝑥 = 1; 𝑦 = 0; 𝑧 = -1 (12)
In terms of the physical properties, the velocity can be expressed as

𝑦̇= 𝐺 1 𝛼 1 𝐺 2 𝛼 2 𝐺 3 𝛼 3 𝐺 4 𝛼 4 𝐺 5 𝛼 5 𝐺 6 𝛼 6 (13)
Substituting the expressions in Table 2 in Eq.( 13) is derived from the correspondent linear system:

� -3𝛼 1 + 𝛼 2 -𝛼 3 -𝑎 5 𝛼 5 + 𝑎 6 𝛼 6 = 1 𝛼 1 +𝛼 3 + 𝑎 5 = 0 -𝛼 3 -𝛼 4 -2𝑎 5 = -1 (14) Finally, 𝑦̇= 𝜇 𝜌𝑟 � 𝜌𝑟 2 𝐸 𝜇 2 � 𝛼 5 � 𝜌𝑟 2 𝜔 𝜇 � 𝛼 4 � 𝑙 𝑟 � 𝛼 6 (17) 𝑦̇= 𝜇 𝜌𝑟 𝑓 �� 𝜌𝑟 2 𝐸 𝜇 2 � , � 𝜌𝑟 2 𝜔 𝜇 � , � 𝑙 𝑟 �� (18)
It is therefore required to keep constant the three dimensionless numbers in Table 3, in order to pursuit physical similitude:

Dimensionless Number Comments � 𝑙 𝑟 �
This geometrical relation can be easily maintained by adjusting the correct canal length to a chosen cross-section radius.

� 𝜌𝑟 2 𝐸 𝜇 2 �
The main restraint of the scale model problem resides in maintaining this second dimensionless number, because it requires an extremely low Young Modulus for the model Cupula Material.

𝑅𝑒 = � 𝜌𝑟 2 𝜔 𝜇 �
The third dimensionless number is the Reynolds number, in this case expressed in terms of angular velocity. Let the scale model variables are referred with an apostrophe. Thus, if 𝑟 means the inner ear cross-section radius, 𝑟 ′ represents the scale model cross-section radius. The scale factor 𝐴 exists so that 𝑟 ′ = 𝐴𝑟 is the new cross section radius. For conserving the second dimensionless number, it may be written

� 𝜌𝑟 2 𝐸 𝜇 2 � = � 𝜌′𝑟 ′ 2 𝐸′ 𝜇′ 2 � (19)
From which

𝐸 ′ = � 𝜇 ′ 𝐴𝜇 � 2 � 𝜌 𝜌 ′ � 𝐸 (20)
Given that

𝜈 = 𝜇 𝜌 (21)
The Young Modulus expression assumes the form

𝐸 ′ = � 𝜈 ′ 𝐴𝜈 � 2 � 𝜌′ 𝜌 � 𝐸 (22)
By naming the linear coefficient

𝑘 = � 𝜈 ′ 𝐴𝜈 � 2 � 𝜌′ 𝜌 � (23)
We have

𝐸 ′ = 𝑘. 𝐸 (24)
The main limitation of the problem is that 𝐸 is considerably small, namely 𝐸~5 𝑃𝑎 . For this reason, the coefficient 𝑘 has to be large enough so that 𝐸 ′ = 𝑘. 𝐸 ≥ 𝐸 𝑚𝑖𝑛𝑖𝑚𝑢𝑚 . In a practical approach, the density 𝜌′ and the scale factor 𝐴 play a minor role in influencing 𝑘, as they don't vary much in magnitude as 𝜈 ′ does (𝜌′ values are generally the same order of magnitude of water density and 20 ≤ 𝐴 ≤ 100 is a project requirement).𝜈′ values, for their part, vary in range length superior to 10 8 units. For solving this Cupula Material limitation, then, it is essential that the model fluid kinematic viscosity 𝜈 ′ be high, so that 𝐸 ′ does surpass 𝐸′ 𝑚𝑖𝑛𝑖𝑚𝑢𝑚 .

As long as 𝑅𝑒 = 𝑅𝑒′, by choosing a large enough value for 𝜈 ′ , 𝜔 is relatively larger than 𝜔' in a scale model experiment when all other variables are kept constant. In other words, the simulation of a chosen head rotation requires an even larger scale model rotation. This fact can be shown by

𝑅𝑒 = 𝑅𝑒′ (25) � 𝜌𝑟 2 𝜔 𝜇 � = � 𝜌′𝑟 ′ 2 𝜔′ 𝜇′ � (26)
Obtaining the simulated angular velocity

𝜔 = �𝐴 2 𝜇𝜌 ′ 𝜇 ′ 𝜌 � 𝜔′ (27)
Or

𝜔 = �𝐴 2 𝜈 𝜈 ′ � 𝜔′ (28)
As ω ′ is superiorly limited by equipment restraints, Eq.( 28) shows that ν ′ cannot be as large as desired, otherwise the project criterion ω ≥ ω min would not be met. As aforementioned, ν ′ should be, though, large enough due to the restraint E ′ ≥ E′ minimum . Concluding, a balance value for ν′ must be established so that both specifications are obeyed. A larger value for A aids less substantially in the rotation problem, but shall be considered in a second analysis.

Materials Choice

For modeling the endolymph, the Polydimethylsiloxane Silicone Oil has been chosen for representing high enough Cupula elasticity modulus, providing experiments with higher rotational rates. Equally, the commercial easiness to find it renders it an attractive choice. The silicone rubber presents other interesting features, to mention: clear, colorless, odorless, fluid; extremely high kinematic viscosities, covering a wide range; little viscosity change with temperature; high temperature stability and high oxidation resistance. It shall also be pointed that the mentioned fluid covers a wide range of possible kinematic viscosities (from less than 0.65 till 20 million cSt), whilst maintaining practically constant density (978 kg/m 3 ).

For the Cupula, low elasticity modulus silicone rubbers, with 𝐸′ ~ 𝐸 𝑚𝑖𝑛𝑖𝑚𝑢𝑚 = 0.005 𝐺𝑃𝑎 have been considered, for highly elastic materials are required for high viscosity fluids. By choosing the material's lowest possible Young modulus, lower part of high kinematic viscosities 𝜈 ′ can be utilized, minimizing the rotation ratio limitation.

The vestibular canal walls can be modeled with transparent flexible plastic hoses, a simple and low-cost solution.

Sample Proposed Simulations

For the largest scale factor, 𝐴 = 100, the expression 𝜔 = �𝐴 2 𝜈 𝜈 ′ � 𝜔′ allows the construction of Table 4, with sample construction possibilities and their respective simulated speeds. A set of experiments could also be determined, depending on the availability of the materials in the market and their price, and as a result, the chosen fluid cinematic viscosity was 𝜈 ′ = 60,000 𝑐𝑆𝑡 = 0,06 𝑚 2 /𝑠 (Wacker® AK 60,000 cSt Silicone Oil). The suggested cupula silicone rubber can be the Smooth-On Ecoflex 00-10, a mold making material whose Young Modulus is close to the required. It is composed of a base and a catalyst that must be mixed and cured to mold the silicone rubber in desired shape.

CONCLUSION

The numerical finite element model provides a visual understanding of the cupula response to external stimulus. It has been seen that the flows highest speeds take place in a flow recirculation region illustrating the pressure mechanism with which the endolymph exerts a force due to the fluid pressure difference on both sides of the membrane. If the single canal model functioned well, the double canal model presented singularities that may render its results dubious. For the triple canal model, the anterior semicircular canal presented a great difference in the model conception for being greatly influenced by gravity.

For future works regarding the vestibular modeling problem, one may leave as a suggestion for different FEM models the elaboration of triple canal geometry where the ampulla/cupula region of the anterior canal be in an horizontal position. This would present better results for analyzing the deflection of this membrane due to a pure horizontal rotation. Furthermore, a model can have an improved realism by modeling the cupulas as not entirely obstructive in the ampullas' endolymph flow. It would indeed provide interesting results by eliminating one boundary fixed constraint for the diaphragm structure and allowing a communication between the left and right side of the flow.

The parametric study for choosing the best materials for the augmented scale model was successfully conducted and provided satisfactory results regarding simplicity, prices and physical similitude. It will be of great value for future visual experiments.

Figure 1 :

 1 Figure 1 : Semicircular canals relative disposition (Roman and Thomassin, 2000)Figure2: Cupula Membrane[START_REF] Corwin | The Vestibular System[END_REF]

Figure 3 :

 3 Figure 3: Canal Excitation In Both Ears[START_REF] Corwin | The Vestibular System[END_REF]

Figure 7 :

 7 Figure 7: One Canal Input Velocity with Time

Figure 8 :

 8 Figure 8: Double and Triple Canal Input Angular Velocity with Time

Figure 6 :

 6 Figure 6: Model Meshed geometry configurations

Figure 7 :

 7 Figure 7: Single canal velocity field in organ reference frame

Table 1 :

 1 Material Properties

	Property [unit]	Symbol Cupula Endolymph
	Young Modulus [Pa] Kinematic Viscosity [m 2 .s] 𝜈 𝐸 Density [kg/m 3 ] 𝜌	5.00 -1020	-0.33 1020

Table 3 :

 3 Dimensionless numbers description

Table 4 :

 4 Possible experiments. Highlighted cells indicate cases where the minimum project simulated velocity ω min = 250°/s is met