F Trivaudey

V Guicheret-Retel

L Boubakar
email: lamine.boubakar@univ-fcomte.fr

Effect of frictional sliding on the unilateral damaged behaviour for laminate composite

Keywords: Closure/opening effect, anisotropic damage, composites

. Some numerical results are presented to illustrate the different situations.

Introduction

For several decades, many models have been proposed to describe composite behaviour, and especially the asymmetry in traction and compression loads [START_REF] Ladeveze | XVIth int. cong. of theoretical and applied mechanics[END_REF][START_REF] Chaboche | [END_REF]7,8]. A previous model of unilateral damage has been developed by the authors [START_REF] Trivaudey | Unilateral Damage Model for Laminate Analysis[END_REF] which was able to report matrix micro-cracking taking account for the cracks states. When cracks were opened, a damaged behaviour was used, when cracks were closed and stuck the virgin behaviour was recovered and when cracks were closed and sliding, virgin behaviour was recovered in compression and damaged shearing behaviour was present. Actually, when cracks were closed with a slight negative transverse stress component and loaded with a large shearing component, no damage increased contrary to expected. The new version of the model is able to let damage increase even with closed cracks.

Unilateral damaged elastic behaviour

Damage variable definition

The damage model presented here is proposed for laminates made of layers of polymer matrix reinforced with long glass fibres. Using a meso-macro approach, the material behaviour is only needed at the level of one ply. It represents the degradation of the elastic properties due to the initiation and the growth of cracks in the matrix along the fibres direction (See fig. 1).

Although the damage is anisotropic, it can be modelled by a single damage variable D differently acting, according to the cracks opening modes through a damage tensor ) (D H (See fig. 2). This tensor increases some components of the virgin compliance tensor S when damage is acting.

The compliance of the damaged elastic material is called S ˆand the elastic strain is obtained as follows: This formulation corresponds to the one developed by Boubakar et al. in [START_REF] Boubakar | A meso-macro finite element modelling of laminate structures. PartI: time-independent behaviour[END_REF]. The evolution law of D is classically obtained from the principle of maximal dissipation. A yield damage function is chosen in the following form:

) ( p c D qD Y Y f + - = (4) p q Y c , ,
being material constants.

Cracks closure

As soon as the cracks close ( give a blocked strain due to the cracks closure. The first part will let the compliance to be either the virgin one or a partly damaged one according to the sign of the sliding function s f 12 . The new coefficient and the sliding function will be discussed further.

(M2) (M1)                     = 66 44 22 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 H H H H 2 x r 22 22 1 S D D H - = 22 11 44 1 S S D D H - = 22 66 1 S D D H - = A 2 x r 1 x r (M1) B 3 x r 2 x r (M2) C 1 x r fibre 2 x r micro-crack (M1) micro-crack (M2) (M2) matrix                     = *

Damage kinetic with closed cracks

When cracks are closed ( )

0 22 < σ
, the virgin elastic behaviour is recovered in compression thanks to 22 S which is kept equal to the value of the virgin material all along the closure of the cracks. At the same time, the shear behaviour is supposed to depend on the contact between the cracks lips. If the compression stress is much more important than the shearing one, cracks are assumed to be stuck and the virgin shearing behaviour should then be recovered through the nullity of the sliding function s ij f . If compression is much weaker than shearing, cracks lips are assumed to be sliding one with respect to the other, leading to a more important strain than the previous case.

The virgin compliance is then increased through ζ and s ij f . In order to record as closely as possible the sliding of the cracks lips, the sliding function is defined as follows:

22 σ σ f f ij s ij - = (6)
where f is a friction coefficient. The strain increment is now split in two parts:

0 γ ε ε + = e (7)
The first term is able to report the virgin/partly damaged behaviour in case of stuck/sliding cracks; the second is constant and obtained at the closure of the cracks. When 22 σ is negative, the free specific energy is chosen as follows:

( ) ( ) ( )       + = - e e e D H S D ε ζ ε ρ ε ψ : : 2 1 , 1 * (8)
where ζ is a parameter which deals with the effect of damage on closed sliding cracks. s ij f is a criterion function to detect whether cracks will be able to slide when they are closed. The damage tensor will act or not according of to the sign of the criterion.

                    = 66 23 44 12 * 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 H f H f H s s where      < ≥ = 0 0 0 1 s ij s ij s ij f if f if f (9)
The state law, the damage force and the intrinsic dissipation becomes now:

( ) e e H S ε ζ ε ψ ρ σ : 1 * - + = ∂ ∂ = with σ σ ζ : ' : 2 1 * H Y = (10)
The damage dissipation is considered in the same way than opened cracks. The ζ coefficient must be positive but smaller than one to record an increase of compliance if cracks slide but not more than when cracks are opened. The damage driving force will be smaller than with opened cracks and as a result, damage increase will be harder with closed sliding cracks. The virgin shear behaviour is recovered as soon as the sliding function is negative, what means when cracks are stuck. When they can slide, a partly damaged behaviour is obtained through the ζ coefficient. The important point here is that damage can increase even with closed cracks if the damage force becomes high enough.

In order to insure the strain continuity when cracks slide, ( )

12
σ ε σ ε + = + = (11)
Dop is the damage value at cracks opening. This could result in a jump of strain due to the blocked part of the strain which is subtly released by the cracks opening.

Numerical results

Different simulations will be presented here in order to show the ability of the model to easily describe different behaviours of a representative volume element of plate according to the cracks states. The results obtained with only tensile-compression stress component in [START_REF] Trivaudey | Unilateral Damage Model for Laminate Analysis[END_REF] are not affected by this new formulation and will not be related here. The only shearing stress is 12 σ . The transverse normal stress 22 σ is chosen negative to close the cracks and weak enough to let them slide. Three criteria have to be taking account for to develop the model. They are written in terms of stress, so as soon as the stress state is known the behaviour is identified. The loading is applied step by step and the strain increment is deduced from the behaviour laws corresponding to the current damage and cracks state.

The numerical results presented here have been obtained with the following fictive material characteristics:

E 1 =45680 MPa; E 2 =16470 MPa; G 12 =76920 MPa, 12 ν =0.3, 23 ν =0.3; C Y =0.0043MPa; q=1.37; p=0.96
The material coefficients related to the sliding strains are studied below.

Effect of the sliding coefficient ζ

To point out the possibilities of the current model, shearing tests are presented here with alternatively closed and opened cracks. The compressive stress value is small enough to let the cracks slide when they are closed. Two kinds of test are carried out: a) shear loading is applied with opened cracks and unloading is performed with closed cracks; b) shear loading is applied with close cracks and unloading with opened cracks.

Two different values ( 15 . 0 = ζ and 0.75) have been chosen for the sliding coefficient. If ζ was equal to 0, the material would recover virgin shear behaviour as soon as cracks would be closed because no damage would be taken into account in the stress-strain law, and no stress would appear in the damage driving force. If ζ was equal to 1, the material would damage in the same way in shearing than with opened cracks as soon as cracks slide. The effect of this coefficient is very sensitive to the value of damage.

During the first unloading of figure 3.a), no difference is visible between the two curves, but during the two last unloading when damage is important, the slope of the plain curve remains high because material is not supposed to be altered a lot by damage with a small value of ζ when cracks are closed, which is different with high value of ζ where the decreasing slope is closed to the one with opened cracks.

In figure 3.b) shear loading is applied with closed cracks which are then opened at different value of the shearing stress. The strain jump records the quick release of the blocked strain and the damage increase as soon as compressive stress disappears. At the cracks opening, both curves reach the curve corresponding to opened cracks.

When ζ is small, the damage driving force can not be high enough to activate damage increase during the loading and as a result, the plain slope is straight up to the opening. On the other hand, when ζ is high enough, it could even make the damage increase as we can see around 120MPa. The unloading curve is common to the two curves for instance in figure5). As soon as cracks start sliding, a slight (f=0.1) or high (f=0.4) jump is observed on the stress/strain curve corresponding to a jump in damage. When the friction coefficient is high (over 0.4), damage brusquely increases when the sliding function becomes positive because the damage force is much larger than the damage criterion and results in a jump of damage. This jump induces a loss of stiffness and results in a jump of shearing strain. Compared to the value of damage reached with opened cracks, the curve of damage with the sliding cracks is always below. The lower value of ζ , the lower the curve is. As far as unloading is concerned, cracks are closed before unload. For f=0.1, cracks are sliding down to 60 MPa following a partly damaged slope. For medium value of f, ( f=0.4) cracks are first sliding (down to 120 MPa) and become stuck when the shearing stress goes on decreasing. When cracks are stuck, the virgin behaviour is recovered. At 12 σ reaches 60MPa, cracks are opened and the damaged behaviour is observed.

Strain driven numerical integration

As soon as the behaviour has to be implemented in a finite element code, the numerical integration is strain driven [START_REF] Trivaudey | Opening/Closure effect on laminates behaviour[END_REF]. In this model, three criteria written in terms of stress have to be used. The point is: knowing the stress state n σ , damage value n D at the end of the previous step and the new strain increment ε ∆ , find the stress increment σ ∆ which let the three criteria be verified. This results in a serial of predictor-corrector schemes to evaluate the stress increment related to the given strain increment. ) and, according to the criterion, a damage increment is computed if necessary and a new stress increment is obtained.

If the new computed stress state verifies the cracks state criterion ( )

0 22 < σ
, the stress increment is kept, otherwise, cracks have been opened during the step. As a result, the step is split in two parts to treat the first part as above (closed cracks) and the second part as opened cracks. Just before computing the second part, the new shear stress is evaluated through the damaged behaviour with the value of damage at the cracks opening.

The figure 6 shows results when the integration is strain or stress driven. The applied strain state, when the computation is strain driven, is chosen in a way to give the same stress state than the stress driven computation. The two stress-strain curves are slightly different. When it is stress driven, the shear stress is kept constant during the opening of the cracks, the opening results then in a jump of strain because material becomes more compliant. When it is strain driven, shear strain is kept constant when cracks open and as a result, shear stress useful to keep the shear strain constant decreases because of the increasing compliance.

Conclusion

A new damage model accounting for the opening/closure effect is presented. This work follows a previously developed model whose weaknesses have been removed by letting damage increase even with closed cracks. It was already able to describe unsymmetrical behaviour in tensile or compression loading. A first criterion based on the sign of the transverse normal stress let decide if cracks are opened or closed. Using a new formulation and adding two new material coefficients and one new criterion, it is now able to report different behaviour according to the status of the closed cracks. The new criterion introduced here compares the absolute value of the compressive stress and the shearing stress and let decide if cracks are sliding or stuck. When cracks are supposed to be stuck, virgin shear behaviour is recovered. When cracks are supposed to be sliding, damage increased is possible as soon as the damage driving force is high enough to violate the damage criterion.

Figure 1 :Figure 2 :

 12 Figure 1 : Micro-cracks orientation in the matrix

  behaviour is modified by the presence of the cracks which are closed but whose lips are still able to slide. A new damaged compliance * S is defined from the virgin compliance modifying only the terms possibly affected by damage:

DS

  is the damage value at the closure. The new compliance can be written in the following way with the sliding function s f 23 ). c D is the damage value at the cracks closure. The second part of * 44 S will

Figure 3 :

 3 Figure 3 : Shearing stress-strain curve. a) loading with opened cracks, b) loading with closed cracks The figure 4 shows how damage is affected by the sliding coefficient. The damage increase with closed cracks can be seen especially for the high value of ζ . The jump of damage can also be pointed out on this figure.

Figure 4 :

 4 Figure 4 : Damage evolution according to the cracks state 3.2 Effect of the friction coefficient f This coefficient let the sliding criterion be more (or not) easily reached to let the cracks slide. If it is high, it can delay the cracks sliding and prevent damage from increasing ( 7 . 0 = f

Figure 5 :

 5 Figure 5 : Stress/strain curve (a) and damage evolution (b) according to the friction coefficient 3.3 Strain driven numerical integrationAs soon as the behaviour has to be implemented in a finite element code, the numerical integration is strain driven[START_REF] Trivaudey | Opening/Closure effect on laminates behaviour[END_REF]. In this model, three criteria written in terms of stress have to be used. The point is: knowing the stress state n

Figure 6 :

 6 Figure 6 : Stress/strain curve(a) and damage evolution (b) according to integration scheme If cracks were opened at the end of the previous increment, we predict they will remain open A new prediction is made with no damage increase If damage criteria is verified ( ) 0 ≤ d f

  or to stick. This explains the presence of the sliding function in the *

					S	44	+	ζ	H	* 44	σ	12	=	S	44	σ	must be verify when cracks
	start to slide H tensor.
	2.5 Cracks opening As soon as the cracks open (	σ	22 ≥	0	), the damaged behaviour is recovered ( ) ( Dop opening H S 23 ) ( 66 66 23 ) opening opening Dop opening H S 12 ) ( 44 44 12

Acknowledgements


Joseph Lardiès
email: joseph.lardies@univ-fcomte.fr

Olivia Arbey

Marc Berthillier

Analysis of the pull-in voltage in capacitive mechanical sensors

MEMS based capacitive type sensors offer advantages due to their small size, relative high sensitivity, batch fabrication capability, low power consumption and low noise features. A MEMS capacitive sensor is basically an electrostatic transducer and an analytical approach is used to model a MEMS-based capacitive-type sensor.

A closed-form model to evaluate the pull-in voltage associated with a clamped square diaphragm and a circular diaphragm subject to electrostatic forces due to a bias voltage is developed. The approach is based on a linearized uniform approximation of the nonlinear electrostatic force due to the bias voltage and the use of a 2D load deflection model for MEMS based capacitive acoustical sensor. The spring hardening effect associated with nonlinear stretching of the central region of a clamped diagram is also considered. The resulting electrostatic pressure on the diaphragm, the pull-in voltage, and the deflection of the midpoint of the diaphragm for different bias voltage are studied. A comparison of the results obtained using the developed analytical model of the communication with the results obtained by Senturia is discussed. Numerical results are presented showing the effectiveness of the method in nonlinear identification problems.

Introduction

Any micro-electro-mechanical system (MEMS) includes input transducers, mechanical resonators and output transducers. The input transducers convert the input electrical signals into an electrical force. The force acting on mechanical resonators sufficiently isolated from their surroundings causes the corresponding vibration of the resonators. The output transducers will sense the motion of the mechanical resonator and generate the corresponding electrical signals. The central component in micro-electro-mechanical systems is the mechanical resonator which constitutes a capacitive transducer and is formed with two plates: a fixed plate and a movable plate. Due to the electrical force, when the gap between the two plates becomes two thirds of the initial gap, the movable plate is not stable, we have a ''push-down'' phenomenon and the MEMS fails. From the dynamics point of view, the system loses its stability. The gap being equal to two-thirds of the initial gap is termed the minimum gap in MEMS [START_REF] Senturia | Microsystem Design[END_REF]. In fact, stoppers are used in the design of MEMS for avoiding such a failure of the capacitor. The main function of MEMS is to transfer input signals to output signals with a specific frequency relationship through a mechanical resonator. If the exact frequency of MEMS cannot be determined, the vibration of mechanical resonator in MEMS cannot be easily induced by input signals with certain frequencies. Furthermore, the output signals cannot be sensed out, which leads to the failure of MEMS. Due to the enormous advance in microelectronics and signal processing, small size, bandwidth, high sensibility, low power consumption and batch fabrication, MEMS capacitive sensors, usually called CMUT (capacitive micromachined ultrasonic transducer) have been an attractive alternative to conventional piezoelectric transducers in the generation and detection of ultrasound in air and water. The CMUT is constructed directly on a silicon substrate by means of CMOS technology and provides the advantage of increased bandwidth with comparable sensitivity to conventional piezoelectric transducers as well as ease of fabrication and electronics integration. The CMUT consists of a lot of very thin square or circular silicon nitride membranes in parallel suspended above silicon substrate. Basically, the single cell of a CMUT is an electrostatic transducer like the well known condenser microphone [START_REF] Senturia | Microsystem Design[END_REF]; nevertheless, this kind of capacitive transducer is able not only to detect but also to generate acoustic waves in fluids. This capacitive-type sensor is then basically an electrostatic transducer that depends on electrical energy in terms of constant voltage (voltage drive) to facilitate monitoring of capacitance change due to an external mechanical excitation such as force or acoustical pressure. The electrostatic force associated with the voltage is non linear due to its inverse square relationship with the airgap thickness between the capacitor electrodes. This gives rise to a phenomenon known as ''pull-in'' that causes the movable structure (membrane) to collapse if the bias voltage exceeds the pull-in limit and limits the effective range of deformation of the structure. Accurate determination of the pull-in voltage, or the collapse voltage, is critical in the design process to determine the sensitivity, harmonic distortion and the dynamic range of a MEMS-based capacitive transducer with a square or a circular diaphragm. Pull-in instability is fundamental to the understating of many MEMS devices and is one of the basic parameters of the design of many electrostatic MEMS devices. Lack of an accurate model for predicting pull-in voltage for basic beam structures such as fixed-fixed beams, fixed-free beams or for basic membranes necessitates the clear need for a closed form expression for the pull-in voltage. Attempts have been made by several authors in the past [START_REF] Senturia | Microsystem Design[END_REF][START_REF] Puers | Electrostatic forces and their effects on capacitive mechanical sensors[END_REF][START_REF] Chowdhury | A new analytical model for the pull-in voltage of clamped diaphragms subject to electrostatic force[END_REF][START_REF] Chowdhury | A closed form model the pull-in voltage of electrostatically actuated cantilever beams[END_REF] to derive a closed-form expression for the pull-in voltage. In this communication an analytical solution is presented to calculate the pull-in voltage and diaphragm deflection for a clamped square diaphragm and a circular clamped diaphragm under electrostatic actuation. The method includes the nonlinearities of the electrostatic force and the large deflection model for a clamped square and a clamped circular diaphragm deflection.

The model

The basic structure of a MEMS capacitive microphone is shown in Figure 1 : a metalized membrane stretched by a tensile force is put in front of a fixed conducting backplate by means of a surrounding border which assures a separation distance. In fact, the structure can be viewed as a parallel plate capacitor consisting of a top diaphragm and a bottom backplate separated by a small airgap acting as the dielectric material. When an acoustical sound wave is incident on the diaphragm, it causes the diaphragm to deflect and the gap between the diaphragm and the backplate decreases causing an increase in the capacitance between them. As the diaphragm vibrates in accordance with the frequency of the acoustical wave, the capacitance between the electrodes changes accordingly due to a variable airgap. If a battery is connected across the diaphragm and the backplate, following the principle of energy conservation, electrical charge will flow to and away from the battery in accordance with diaphragm vibrations. By connecting a suitable charge flow sensing electrical circuit to the system, a usable voltage signal representation of the incident acoustical wave can be obtained. The fixed plate of the capacitor with area A is connected with a constant supply voltage V. The other plate of the capacitor with mass m and area A is movable and rigid. The support of the moving plate is modeled through an equivalent spring with stiffness k. Without any electrostatic force, the gap between two plates of the capacitor in the MEMS is d 0 . The coordinate x is shown in Figure 2 and the origin is at equilibrium without any voltage. By applying a voltage across the plates, an electrostatic attractive force F e (x) is induced which leads to a decrease of the gap spacing, thereby stretching the spring. This results in an increase of the mechanical elastic force (or spring force) F m (x) which counteracts the electrostatic force. Pull-in instability occurs as a result of the fact that the electrostatic force increases non-linearly with decreasing gap spacing, whereas the mechanical elastic force is a linear function of the change in the gap spacing. In simple terms, the pull-in voltage can be defined as the voltage at which the restoring spring force can no longer balance the attractive electrostatic force. Neglecting any damping within the system, the equation of motion of the movable plate due to an electrostatic attraction force F e (x) caused by a constant supply voltage V is:

m 2 2 t d x d +kx = F e (x) (1)
The mechanical elastic force is F m (x) = kx and the capacitance of the movable-plate parallel-plate capacitor is given by :

C(x) = x d A 0 0 - ε (2)
where 0 ε is the permittivity of the free space. The electrostatic attraction force F e (x) between the plates due to the charges on the plates can be found by differentiating the stored energy of the capacitor with respect to the position of the movable plate and is expressed as :

F e (x) =       ) x ( C V 2 1 x d d 2 = 2 0 2 0 ) x d ( 2 V A - ε (3)
At the static equilibrium the mechanical elastic force equals the electrostatic attraction force and the relationship between the voltage V and displacement of the movable plate is :

V=(d 0 -x) ) A ε ( / x k 2 0 (4)
The maximum of the voltage is obtained for dV/dx=0 and from this equation we obtain the distance where the pull-in occurs:

dV/dx = - ) A ε ( / x k 2 0 +(d 0 -x)k ) A ε x k 2 /( 1 0 (5)
The distance where the pull-in occurs is x pi = d 0 /3 and the pull-in gap is d pi = 2d 0 /3. The pull-in voltage for this ideal parallel plate structure is then :

V pi = ) A ε 27 ( / d k 8 0 3 0 (6)
and the spring constant of the movable plate is given by:

k = 27 0 ε AV pi 3 / (8 d 0 3 ) (7)
A graph of the normalized voltage as a function of the normalized displacement of the movable plate is shown in Figure 3. A stable displacement occurs for a normalized gap greater than 2/3 (or a normalized displacement of the movable plate smaller than 1/3).

Figure 3. Variations of the normalized voltage for the electrostatic actuator

If the applied voltage is increased beyond the pull-in voltage, the resulting electrostatic force will overcome the elastic restoring force and will cause the movable plate to collapse on the fixed plate and the capacitor will be short circuited.

To obtain stiffness due to the electrostatic force we expand equation (1) using a Taylor series approximation about the nominal distance x 0 :

m 2 2 t d x d +         - - 3 0 0 2 0 ) x (d V A ε k x =         ∑ - - - - - = - - N 3 n 1 n 0 0 1 n 0 0 0 0 2 0 0 2 0 ) x d ( ) x x ( n ) x d ( x 2 1 ) x d ( V A ε 2 1 ( 8
)
The electrostatic attraction force effectively modifies the spring constant k of the movable plate and the effective spring constant at a specified voltage V is :

k effective =         - - 3 0 0 2 0 ) x (d V A ε k (9)
The resonant frequency of the structure is shifted from m / k res = ω to m / k ( effective res = ω .

3 Deflection of a square diaphragm Due to the presence of residual stress and a significantly large deflection of the diaphragm compared to its thickness, the developed strain energy in the middle of the diaphragm causes a stretch of the diaphragm middle surface. Figure 4 shows the principle of deflection of a diaphragm. The deflection of the diaphragm middle surface corresponds to a nonlinear behavior of a rigidly clamped diaphragm and is known as spring hardening. Tabata et al [START_REF] Tabata | Mechanical property measurements of thin films using load-deflection of composite rectangular membranes[END_REF] developed an analytical solution for the loaddeflection of membranes. They found a relationship between the external load and the membrane deflection to determine the residual stress and Young's modulus of thin films. Pan et al [START_REF] Pan | Verification of FEM analysis of load-deflection methods for measuring mechanical properties of thin films[END_REF] compared the analytical solution with FEM analysis. They found that the functional form of the analytical results is correct, but some constants have to be corrected. Pan et al also found that the analytical forms of the membrane's bending lines do not describe the real behavior very accurately. Maier-Schneider et al [START_REF] Maier-Schneider | A new analytical solution for the load-deflection of square membranes[END_REF] found an analytical solution for the loaddeflection behavior of a membrane by minimization of the total potential energy and assuming a functional form of the deflected membrane shape which extracted form [START_REF] Timoshenko | Theory of Plates and Shells[END_REF]. The deflection of the midpoint of a square diaphragm under a uniform pressure load P can be expressed as [START_REF] Timoshenko | Theory of Plates and Shells[END_REF] :

P(h 0 ) = C 1 0 2 h a e σ + C 2 ( ν ) 3 0 4 h a E e ( 10
)
where P is the applied uniform pressure, h 0 the deflection of the diaphragm midpoint, e the diaphragm thickness, 2a the diaphragm side length, E the Young's modulus, ν the in-plane Poisson's ration and σ the residual stress.

The dimensionless constants C 1 and C 2 are numerical parameters which are obtained from [START_REF] Timoshenko | Theory of Plates and Shells[END_REF]: C 1 = 3,45 and C 2 ( ν ) = 1,994(1-0,271 ν )/(1ν ). The non-linear spring constant of the diaphragm is then :

k nl = 0 h A ) 0 h ( P = (C 1 2 a e σ + C 2 ( ν ) 2 0 4 h a E ê ) A (11)
where E ˆ=E/(1-ν 2 ) is the effective Young's modulus. The deflection-dependent nonlinearity due to spring hardening appears in equation ( 11) where the square of the midpoint deflection variable h 0 is considered. For a test device we consider the parameters given in table 1. A plot of the variations of the non-linear spring as a function of the deflection of the diaphragm midpoint is shown in figure 5. From this figure we can obtain the value of the non-linear spring. Note that this spring hardening has been obtained under the hypothesis of an applied uniform pressure.

Figure 5. Spring constant and diaphragm deflection

The deflection of the diaphragm from mid-side to mid-side can be calculated from midpoint deflection h 0 as [START_REF] Timoshenko | Theory of Plates and Shells[END_REF] :

h(y,0)= h 0               + 2 a y 401 , 0 1 cos         π a 2 y ( 12
)
For a parallel plate configuration, the nonlinear electrostatic force is always uniform. However, for a rigidly clamped diaphragm, the electrostatic force becomes non-uniform due to a hemispherical deformation profile of the diaphragm (Figure 6). Thus, to evaluate the deflection of a rigidly clamped diaphragm under an electrostatic force, it is necessary to obtain a uniform linear model of the electrostatic force that can be applied in loaddeflection equation (10). A uniform linearized model of the electrostatic force can be obtained from ( 8) by linearizing the electrostatic force about the zero deflection point x=0. We obtain then:

        + ε 3 0 2 0 2 0 d x d 2 1 V A = k x (13)
The left-hand of this equation equals an approximate uniform electrostatic force and the uniform electrostatic pressure on the diaphragm is :

0 0.5 1 1.5 2
P unif =         + ε 3 0 2 0 2 0 d x d 2 1 V ( 14
)
This equation is general and can be applied to any sort of diaphragms. The pull-in deflection is obtained for x=d 0 /3=h 0 and the effective pull-in pressure is :

P(h 0 ) = C 1 2 0 a 3 d e σ + C 2 ( ν ) 4 3 0 a 27 d E ê = 2 0 2 pi 0 d 6 V 5 ε (15)
and the pull-in voltage for the square clamped diaphragm is given by :

V pi =         + 4 3 0 2 0 1 0 2 0 a 27 d E ê ) (ν 2 C a 3 d σ e C ε 5 d 6 ( 16
)
Figure 7 shows the variations of electrostatic and elastic forces as a function of the displacement x, obtained from (13). The spring constant k is obtained using equation (11) or Figure 5. At a distance of one-third of the original airgap, the elastic force graph intersects the 17.45 volts constant voltage graph and thus gives the voltage where the pull-in occurs. The pull-in voltage can be calculated using equation ( 16). If we use equation ( 6) we obtain V pi = 15.02 volts, a value which is 2.43 volts smaller than the value obtained with the method proposed in this work. The mid-side to mid-side deflection profiles of the diaphragm for different voltages is shown in figure 8. The analytical method developed in the preceding section to determine the pull-in voltage for a rigidly clamped square diaphragm can be adopted to evaluate the pull-in voltage for a circular diaphragm under electrostatic load. Figure 10 shows a clamped circular diaphragm of radius a and thickness e under a uniform transverse load P.

1 0 V 1 0 V 1 5 V 1 7 .
Figure 10. Schematic of a clamped circular diaphragm For a rigidly clamped circular diaphragm the deflection of the midpoint of the diaphragm under a uniform pressure P can be expressed as [START_REF] Sheplak | A MEMS microphone for aeroacoutics measurements[END_REF] :

h(r)= ( )         - + - ν - 2 2 2 3 4 2 a 4 r a ) T ( 1 I T 2 ) T ( 0 I a / r T 0 I e E 2 T a P ) 1 ( 12 ( 17
)
where r is the radial position, h(r) the defection at a radial position r, T the tension parameter, I 0 and I 1 are the modified function Bessel functions of first kind of zeroth and first order respectively. The tension parameter is defined as:

T = E ) 1 ( 12 e a 2 σ ν - (18)
The deflection of the center of the diaphragm is obtained from equation (17): and we obtain the applied pressure P(h 0 ) in the format of equation (15) :

P(h 0 ) = 1 ) T ( I T 2 ) T ( I 4 1 a ) 1 ( 12 0 h e E T 1 0 4 2 3 2 - - ν -       (20)
Using the same method as the square diaphragm case, and considering h 0 in the pull-in deflection (h 0 =d 0 /3) we obtain then : 22) provides a model that includes spring-hardening phenomena and a linearized approximation to the nonlinear electrostatic force, resulting in a more accurate determination of the pull-in voltage for circular diaphragms.

Conclusion

A new relatively simple closed-form model to evaluate the pull-in voltage associated with a rigidly clamped square or circular diaphragm subject to an electrostatic force has been developed. The method incorporates the nonlinear and non uniform nature of the electrostatic force associated with a clamped diaphragm deformation. However, the approach is based on a linearized uniform approximation of the nonlinear electrostatic force due to a voltage an the use of a 2D load deflection model. The spring hardening effect associated with deflection of a clamped diaphragm has been evaluated. The method can be extended to determine the pull-in voltage for other microstructures such as cantilever beams under electrostatic excitation.

Figure 1 .

 1 Figure 1. Cross-section of a MEMS-based capacitive type acoustical sensor backplate

Figure 2 .

 2 Figure 2. A simplified mechanical model of the MEMS with two parallel plates

Figure 4 .

 4 Figure 4. A schematic of electrostatic actuation

Figure 6 .

 6 Figure 6. Non-uniformity of the pressure profile (top) Unifomity of the pressure profile (bottom)

Figure 7 .

 7 Figure 7. Electrostatic and elastic forces

Figure 8 . 4 D

 84 Figure 8. Mid-side to mid-side deflection profiles of a clamped square diaphragm

Figure 9

 9 Figure9shows the pull-in voltage as a function of the airgap thickness. For a value of d 0 = 3.5 m µ we obtain the pull-in voltage of 17.45 volts.

Figure 9 .

 9 Figure 9. Tension of collapse for a clamped square diaphragm

  Figure11shows the pull-in voltage for a range of diaphragm radius.

Figure 11 .

 11 Figure 11. Pull-in voltage for clamped circular diaphragms Equation (22) provides a model that includes spring-hardening phenomena and a linearized approximation to the nonlinear electrostatic force, resulting in a more accurate determination of the pull-in voltage for circular diaphragms.

Table 1 .

 1 Model parameters

	parameter	e	a	d 0	E	ν	σ
	value	0,8	1,2	3,5	169	0,28	20
	unity	m µ	m µ	m µ	GPa	-	MPa


J Lardiès

H Ma

M Berthillier

Localization of radiating sources by an acoustical array

Arrays of sensors are used in many fields to detect signals, to resolve closely spaced targets, to estimate the bearing, the position, the strength and other properties of radiating sources whose signals arrive from different directions. The purpose of the work described in the communication is to provide estimators which are used to localize acoustical sources by an acoustical array of sensors. These estimators are based on acoustic processing algorithms: the conventional beamforming, the Minimum Variance and the constrained Capon algorithm. It is shown that these adaptive algorithms can detect radiating sources which are not detectable by classical weighted arrays of sensors. Numerical results are presented.

Introduction

Arrays of receivers are used in many fields : in radio, radar, sonar, seismic exploration, chemical detection, ultrasonic diagnosis, and so forth, to detect weak signals, to resolve closely spaced targets and to estimate the bearing and other properties of a signal source. The receiving sensors may be any transducers that convert the receiving energy to electrical signals. The type of sensors used to detect these signals differ accordingly: microphones for the acoustic signals, electromagnetic antennas for radio waves, accelerometers/seismometers for the detection of earthquakes, ultrasonic probes and X-ray detectors in medical imaging, containers with membranes or biosensors for gas and vapour detection and so forth. In all these highly diverse applications of array signal processing, the sensors are designed with one basic objective in mind: to provide an interface between the environment in which the array is embedded and the signal processing part of the system and the physical manner in which this interface is established depends on the application of interest. In these applications, the goal is to determinate the distribution of the emitted energy in the medium (air, water, rock, etc) that surrounds the array. For example, in industrial environments to localize complex noise fields, in underwater surveillance with sonar systems, in communications to separate speakers, in seismology for the monitoring and analysis of global earthquakes, in all these cases signals received by sensor arrays are processed to obtain estimates of their strength and direction of arrival. In this communication we are concerned in industrial noises localization which have significant effects on labor's health and community living quality. It is highly desired to develop methods that are capable of locating noise sources in an accurate and systematic manner before any noise control measure can be applied. In all that follows, we consider only idealized arrays, those in which there are no cross-couplings between individual receivers. The receivers are assumed to sample the spatial field without distorting it. We further assume that the spatial field in the vicinity of the array is homogenous, an assumption that is often valid when all signal sources are distant from the array, and that all receivers have the same sensitivity. In this communication we are concerned in industrial noise localization which has significant effects on labor's health and community living quality. It is highly desired to develop methods that are capable of locating noise sources in an accurate and systematic manner before any noise control measure can be applied. Conventional ways of noise source identifications include, for example, sound pressure measurement, sound intensity measurement and acoustic holography. These methods suffer from the drawbacks of being either inaccurate or being restrictive in only small areas or short distances when applied to complex noise fields in industrial environment. An alternative is the use of an array of sensors whose outputs are processed. The performances of these methods in the localization of closely spaced sources are presented. The approach taken here is to assume that the signal field at the array is comprised of K independent plane wave arrivals from unknown directions and the problem reduces to estimating the K directions in a background noise environment. One main of the work described in the paper is to provide estimators which are simple to implement on line. These estimators are based on acoustic processing algorithms: the conventional beamforming estimator, the Minimum Variance (MV) estimator and the robust Capon estimator.

Beamforming

Conventional beamformer

The principles of conventional beamforming were established many decades ago. Consider a linear array as shown in figure 1. A plane wavefront arriving from a direction normal to the array (broadside array) arrives at all the receivers in phase, and so if summed the signal outputs from the receivers will add in phase and reinforce one another; signals from other directions will not be in phase, and so will not be reinforced. If it is desired to receive narrowband signals from some other direction than broadside, we simple shift the phases of the receiver outputs by appropriate amounts to bring all the receiver outputs into phase before summing. More generally, if the signals are broadband, we time delay the receiver outputs appropriately before summation. This process is called beamsteering. It can readily be seen that the same principles can be extended to arrays of arbitrary geometry in two or three dimensions: we simple insert the appropriate time delays.

We consider a linear array of N equally spaced elements, located along the horizontal axis, with an element spacing of d. The beamformer output usually called the space or array factor is given by [START_REF] Elliot | Antenna theory[END_REF]:

y = ) , ( i i N 1 i 0 θ θ, k v a ∑ = = ))] 0 θ cos - θ (cos( d k 1) i ( j exp[ a - ∑ = N 1 i i ( 1
)
where a i are weighting coefficients used to control the shape of the directivity diagram, 0 θ determines the direction of the maximum output (the direction of main beam) and k the wavenumber of interest. Note a and v the columns vectors whose entries are { } i a and { } i v . Equation (1) may be written as :

y = a T v = v T a (2) Source direction d θ N 3 2 1
Figure 1. A linear array of N elements where T is the transpose vector. In the case of a uniform weighting array we have a i = 1 , i=1, 2,…,N and it can be shown, that in order to avoid grating lobes, the condition for a scanned beam becomes :

) cos(θ 1 1 d/λ 0 + ≤ (3)
If we plot the average power output as a function of signal direction, we obtain the polar diagram or polar response. By convention, the polar response is normalized so its response in the beamsteered direction 0 θ is unity. The beamwidth, that is the width of the main lobe, often defined as the angle between the points at which the power response falls by half, indicates the ability to resolve closely-spaced signal sources; the level of the sidelobes is a measure of susceptibility to interference by unwanted signals in directions away from the main beam. The directivity pattern, also known as the polar diagram is obtained by taking the normalized magnitude of y as a function of the polar angle θ for a given frequency. For a uniform weighting array we have the normalized magnitude:

D = ) 2 / Ψ sin( N ) 2 / Ψ N sin( (4) whereΨ = ) θ cos θ cos ( d k 0 -
. For broadside arrays, where the beam maximum occurs in the direction perpendicular to the array axis, we have 0 θ = π/2 . For end-fire arrays where the beam maximum occurs in the direction of the array axis, we have 0 It is instructive to examine the behaviour of the directivity diagram. In particular, it can be manipulated to yield information about locations of major maxima or beams, locations of nulls, angular widths of major maxima or beamwidths and locations of minor maxima or sidelobes.

θ = 0. Figures 2 (a

Dolph-Tchebychev and binomial beamformer

Shading is a well-known technique used to modify the polar diagram of an array and low sidelobes are usually desirable to reduce the effect of interfering signals. Dolph [START_REF] Elliot | Antenna theory[END_REF] addressed that problem by noting the behaviour of the Tchebychev polynomials and devising suitable transformations of variables to link the behaviour of the polynomials to array side lobe levels. With the Dolph-Tchebychev shading the sidelobes all have the same peak level. For an even number of sensors the beamformer output is:

y = 2 ] 2 / Ψ ) 1 i 2 cos[( i 2 1 i - ∑ = N/ a (5)
We consider t = cos( Ψ /2) and use the Moivre formula to write each cos(p Ψ /2) function as a polynomial of degree p, called Tchebychev polynomial and noted T p (t). With this Tchebychev polynomial equation ( 5) becomes

y = 2 ) ( 1 i 2 T i 2 / 1 i t a N - ∑ = = T n-1 (t Z 0 ) (6)
where Z 0 is associated to the sidelobes level. By identification we obtain the weighting coefficients a i . Figure 3(a) shows the end-fire polar diagram of the Dolph-Tchebychev shading with six elements, d/ λ = 0,4 and a constant sidelobe level of 20 dB. The sidelobes all have the same peak level which has been selected to be much lower than the levels for the conventional beamformer.

In order to suppress the sidelobes we use weighting coefficients p N 1 -C extracted of the binomial expansion of

(1 + x) N-1 .
For an even number of elements the beamformer output is:

y= ] 2 / Ψ ) 1 i 2 cos[( C 2 / 1 i i - /2 1 - - ∑ = N N N = ) 2 / Ψ ( cos 2 1 1 - - N N (7)
The fundamental property of the binomial weighting concerns the absence of sidelobes in the polar diagram; this is because the directivity function has a general form in cos n-1 ( Ψ /2) and this function does not oscillate between 0 and π/2 . Figure 2

Superdirectivity and optimal beamforming

It is sometimes desirable to use a given array at low frequencies. Conventional beamforming then yields a broad beam. However, by appropriate choice of the weight vector a, it is possible to generate narrow beamwidths and the array is called superdirective. But there is a penalty : as the frequency is reduced, the magnitude of the weight required becomes large, making the practical implementation of the superdirectivity difficult. Schelkunoff [START_REF] Schelkunoff | A mathematical theory of linear arrays[END_REF] showed that every linear array with commensurable separations between the elements can be represented by a polynomial and every polynomial can be interpreted as a linear array. The analytical representation of an output array is accomplished with the aid of the following transformation: z=exp(j Ψ ) and we have :

y = ∑ = 1 0 i i i z a - N (8)
This polynomial can be written in terms of its (n-1) roots (z 1 , z 2 , …, z N-1 ) and Schelkunoff has shown that if these roots are regularly spaced between 0 et - The beamwidth, or angular space between the half-power points on each side of the main lobe, is narrower than that of a uniform, Dolph-Tchebychev or binomial weighted end-fire array. We address now the problem of selecting a set of weights which maximize the directivity index of a passive acoustical array. The directivity index is defined by:

D (f) = ∫ Ω π ) ( dΩ 2 ) Ω , ( 4
1 2 ) 0 Ω ( f y f, y = a ) W(f a a V(f a T * *T ) 0 , Ω (9)
where

2 ) 0 Ω ( f, y = a ) 0 Ω V(f, a *T
is the power out of the beamformer steered in direction implied Ω 0 , or in the direction of the main beam, and

∫ Ω π ) ( dΩ 2 ) Ω , ( 4
1 f y = *T
a Wa the average output power. In these expressions we have :

V(f, Ω 0 )=v * ( f, Ω 0 )v T ( f, Ω 0 ) and W ij = ∫ Ω Ω ) ( d π 4 1 j v * i v
. The denominator of (9) is the output power when the acoustical array is in a far field of plane waves arriving from all directions and having the same level. This is also the definition of the omnidirectional noise. In fact, the directivity index represents also the output signal-to-noise ratio which may be defined as the ratio of the power received per unit solid angle in the direction of the signal to the average noise power received per unit solid angle. The purpose is to determine an optimum set of weights { } i a in amplitude and phase to yield the maximum directivity index with respect to a prespecified number of elements, spacing and frequency. The index directivity is expressed as the ration of two quadratic forms and it is shown [START_REF] Gantmacher | The theory of matrices[END_REF] that maximisation of ( 9) is obtained when: a = W -1 v and the directivity index However, these weighted arrays have a fundamental resolution problem in the localisation of different sources.

is : D Max (f) = v *T W -1 v.
Figure 5 shows the performances of the optimal end-fire array in the localisation of two sources at 0° and 20° relative to the antenna axis. An important beamwidth and important sidelobes are clearly present in this spectra and it is impossible to obtain the direction of arrival of these two sources. 3 Adaptive array

Spatial correlation matrix

The receiving array has N omnidirectional sensors and is immersed in an acoustic noise field which consists of K independent discrete sources. Because of the geometric positions of the sensors, the total signal power incident on each sensor is the same, but the phase information is different on each receiver.

x(t) = ∑ = K 1 i b( i θ ) s i (t) + n(t) ( 10
)
where s i (t) is the complex amplitude of the i th source. It is a zero-mean complex random variable : E[s i (t)] = 0. The signal power p i of the i th source which we wish to localize is represented by its variance

p i = Var[s i (t)]=E[s i (t)s i (t) * ].
Here E[ ] denotes an ensemble average and the superscript* represents the complex conjugate. The direction of arrival of the i th signal is represented by the N-dimensional vector b )

( i θ =[ T ] ) i θ b .... ) i (θ b ) i (θ b ( N 2 1
, often called the array manifold vector or the steering vector or the directional mode vector and n(t) is the additive noise. The noise is assumed to be spatially white (independent or uncorrelated from sensor to sensor) and the same power level of noise is present in each receiver. With these assumptions, the cross-spectral matrix for the noise alone is R N = E[n(t)n(t) H ] = p N I where p N is the noise power, I is the (NXN) identity matrix and the superscript H denotes the complex-conjugate transpose operation. Equation (10) may we rewritten in the matrix form

x(t) = B s(t) + n (t) t ∈ {t 1 , t 2 ,…t N } (11)
The (NXK) matrix B where each column is a source direction vector, is the so-called array manifold matrix. For any single plane wave arrival, the outputs from the N individual receivers will differ in phase by an amount determined by the geometry of the array and the arrival direction. In other words, the elements B qr of the matrix B are functions of the signal arrival angles and the array elements locations. Thus, one has B qr = exp(j φ qr ) where φ qr is the phase of the signal at the q th receiver from the r th source, measured relative to some arbitrary reference point. That is, B qr depends on the q th array element, its position relative to the origin of the coordinate system, and its response to a signal incident from the direction of the r th source. s(t) is the K-dimensional vector, the components of which are the complex amplitudes of the sources. It can readily be seen that the output signal from the q th sensor may be written as :

x q (t) = ∑ = K 1 r B qr ( θ ) s r (t) + n q (t) (12)
Since the K arrivals are by assumption independent, the correlation matrix between the different signal sources is

R S =E[s(t)s(t) H ] = diag(p 1 , p 2 ,…,p K ) (13)
and at the operating frequency, the spatial correlation matrix (or covariance matrix) of the receiver outputs may be expressed, for signals uncorrelated of each other and of noise, as

R = E[x(t)x(t) H ] = BR S B H + R N (14)
In practice, the spatial correlation matrix is estimated by a finite number of time domain samples (snapshots) and the following estimated form is used :

R ˆ= ∑

= M 1 i M 1 x(t i )x(t i ) H (15)
Where x(t i ) is the array signal vector sampled at time t i and M is the number of such samples. The caret (^) denotes an estimated value.

We can now derive estimators for source localization using an acoustical array.

The minimum variance estimator

The minimum variance (MV) estimator was originally proposed by Capon who conducted a frequency wavenumber analysis on earthquake data analysis. The conventional beamformer can be considered as a kind of linear spatial filter with data-independent coefficients. In contrast, the minimum variance method (or Capon's method) can be considered as a kind of data-dependent spatial filter, in which the coefficients are chosen such that the filter has constant gain at a particular direction while its output power is minimized. The underlying principle of the method amounts to finding an optimal steering vector w opt such that the array output power is minimized while maintaining the gain along the look direction to be constant, say unity [START_REF] Stoica | Introduction to spectral analysis[END_REF]. That is,

min w H R w (16) w subject w H b( θ ) = 1 (17)
Minimizing the resulting beam energy reduces the contributions to this energy from sources and noise not propagating in the direction of look. The solution of this constrained optimization problem occurs often in the derivation of adaptive array processing algorithms. The solution technique is to use a Lagrange multiplier λ and a cost function

H(w) = w H R w + 2 λ ( w H b( θ ) -1) (18)
The gradient with respect to w is

grad H(w) = 2 R w + 2 λ b( θ ) (19)
and the minimum of the cost function is obtained when

w opt = -λ R -1 b( θ ) (20)
But w opt must satisfy the constraint, one then has

λ = (-b H ( θ ) R -1 b( θ )) -1 ( 21
)
w opt = R -1 b( θ )(b H ( θ ) R -1 b( θ )) -1 (22)
and the corresponding array output power is

p MV ( θ ) = E[ 2 ) t ( y ] = (b H ( θ ) R -1 b( θ )) -1 (23)
The goal of the minimum variance estimator is that the contributions of the signals from directions other than θ to the array output are minimized while the signal at direction θ passes through without any distortion. Equation ( 23) is also known as the angular power spectrum of the minimum variance estimator and the signal directions are found by the locations of the spectrum peaks. The peaks level of the spectrum give a good estimate of the true targets power and the spectrum also has uniformly low sidelobes. Simulations have shown that this estimator gives satisfactory resolution if the number of snapshots is high. The algorithm does not require any knowledge of the number of sources present and can also be used with irregular arrays. It is expected that this estimator performs better than the classical beamforming and has super-resolution provided that the SNR is moderately high, the sources are not strongly correlated and the number of snapshots is sufficient.

The common advantage of this estimator is that it does not assume anything about the statistical properties of the data and, therefore, it can be used in situations where we lack information about these properties. Figure 7 shows the power spectral density of the adaptive array with six sensors and d/ λ = 0,4. Two sources at 0° and 20° relative to the antenna axis are present are localized with this estimator. The resolution of the antenna array can be improved using a robust adaptive array.

The robust minimum variance estimator

The minimum variance estimator adaptively selects the weight vector to minimize the array output power subject to the linear constraint that the signal of interest does not suffer from any distortion. This beamformer has a better resolution and much better interference rejection capability than the standard beamformer, provided that the array steering vector corresponding to the signal of interest is accurately known. However, the knowledge of the signal of interest steering vector can be imprecise, which is often the case in practice due to the differences between the assumed signal arrival angle and the true arrival angle. Whenever this happens, the minimum variance estimator can give poor results. A robust minimum variance estimator is then proposed based on a diagonal loading [START_REF] Li | On robust Capon beamforming and diagonal loading[END_REF] where a scaled identity matrix is added to the covariance matrix given in (14). The main idea of diagonal loading is to replace the covariance matrix R by the matrix R+rI where the diagonal-loading factor r is a user-selected parameter. The weight vector so obtained is given by : w = (R+ rI) -1 b( θ )[b *T ( θ )(R+rI) -1 b( θ )] -1

(24)

The robust minimum variance power estimate of array output is then :

P RMV ( θ )

Conclusion

Beamforming techniques are useful for focusing sound waves but are not adapted for the localization of multiple sources. Two methods using adaptive arrays have been proposed and are based on the minimum variance and the robust minimum variance estimators. These algorithms have been tested by computer simulation studies and they perform satisfactory. In the future the resolution capability of these algorithms and their limitations in a broadband environment will be studied.

Figure 2 .

 2 Figure 2. Directivity diagrams for broadside (a) and end-fire (b) arrays

Figure 3 .

 3 Figure 3. Directivity diagrams for Dolph-Tchebychev (a) and binomial (b) weighting

,

  the visible range of z, we have a superdirective array.

Figure 4 (

 4 Figure 4(a) shows the end-fire directivity pattern with a Schelkunoff distribution, six elements and d/ λ = 0,4.

Figure 4 (Figure 4 .

 44 Figure 4. Directivity diagrams for superdirective (a) and optimal (b) weighting A comparison of different weighting techniques developed in the communication are presented in table 1 for a linear end-fire array of six elements and d/ λ = 0,4.

Figure 5 .

 5 Figure 5. Direction of arrival of two sources using the optimal beamforming

Figure 6 .

 6 Figure 6. The setup of the source location problemThe purpose of any estimator is to use the phase information in some way to infer which signals reached the receiving array and the goal of sensor array source localization is to find the locations of sources of wavefields that impinge on the acoustical array. The available information is the geometry of the array, the parameters of the medium where wavefields propagate and the time measurements or outputs of the sensors. For purposes of exposition, we first focus on the narrowband scenario. For a set of K sources, the signals observed at the outputs of the N sensors array are represented by the N-dimensional vector[START_REF] Stoica | Introduction to spectral analysis[END_REF][START_REF] Lardiès | Power estimation of multiple sources by an acoustical array[END_REF][START_REF] Lardiès | Performance of high-resolution sensor array processing algorithms in the localisation of acoustic sources[END_REF]

Figure 7 .

 7 Figure 7. Direction of arrival of two sources using the adaptive estimator

  Figure8shows the power spectral density of the robust adaptive array with six sensors and d/ λ = 0,4. Two sources at 0° and 20° relative to the antenna axis are present are localized with this estimator which has a better resolution than the minimum variance estimator.

Figure 8 .

 8 Figure 8. Direction of arrival of two sources using the robust adaptive estimator

Table 1 .

 1 Comparison of weighting methods

	weighting	uniform	Dolph-Tchebychev	binomial	Schelkunoff	optimal
	beamwidth	84 °	90 °	126 °	78 °	58 °
	sidelobes	important	low	absent	low	low


Mohamed Bellaredj

Marc Berthillier

Joseph Lardies

Brahim Belgacem

Fabrication and experimental characterization of capacitive micromachined ultrasonic transducers made with anodic-bonding technology

Fabrication and experimental characterization of capacitive micromachined ultrasonic transducers made with anodic-bonding technology

M.L.F.Bellaredj, M. Berthillier, J. Lardiès, B.Belgacem, G.Bourbon FEMTO-ST, 24 Chemin de l'Epitaphe, 25000 Besançon France Département de Mécanique Appliquée (Mohamed.Bellaredj@femto-st.fr)

Capacitive micromachined ultrasonic transducers CMUTs introduced about two decades ago have been shown to be a good alternative to conventional piezoelectric transducers in various aspects, such as sensitivity, transduction efficiency, bandwidth etc… CMUT arrays are commonly fabricated by means of surface micromachining or wafer bonding techniques. Surface micromachined cMUTs have a common drawback, i.e., their electrode gap may not be made too large.

Otherwise the deposition period of the sacrificial layer will be unacceptably long. Small electrode gap itself does not constitute problem to transducer operation, but it does influence its fabrication, because the vibrating membranes are usually released by wet etching process, whereas in small distance range, capillary force plays an important role. Another problem is its inferior surface flatness. Surface micromachining consists of successive deposition and patterning. During this procedure, the wafer's surface will repetitively be modified. Unless extra planarization steps are introduced, the surface topology can never be flat. Such unevenness is also disadvantageous to the uniformity of the transducer [START_REF] Jia | Silicon Micromachined Ultrasonic Transducers[END_REF]. Direct bonding places stringent requirements on the cleanliness and flatness of the surfaces to be bonded and it requires a high temperature annealing to establish high strength bonding [START_REF] Thomas | Detailed characterization of anodic bonding process between glass and thin-film coated silicon substrates[END_REF]. To realize reliable bonding, the contact areas roughness should not exceed 1~2nm, and the wafer flatness should be smaller than 2μm. Furthermore, to avoid bonding voids, no particles should be present on the wafer's surface [START_REF] Jia | Silicon Micromachined Ultrasonic Transducers[END_REF]. Moreover, the use of a silicon wafer as mechanical support of the transducer is not optimal for the device operation, as it generates ringing in the pulse-echo signal because of its high acoustic impedance [START_REF] Coppa | Building CMUTs for imaging applications from top to bottom[END_REF] To solve these drawbacks, a CMUT transducer has been designed and realized using a new fabrication process called anodic bonding introduced in [START_REF] Clatot | Concéption,fabriacation et caractérisation de transducteurs capacitifs micro-usinés pour l'imagerie ultrasonore[END_REF]. It's a particularly promising technique for CMUTs technology because it does not require either a postheat treatment at high temperature or hydrophilic cleaning. Moreover, it's easy and reliable to control and requires no expensive installation. Glass-to-silicon stacks are commonly applied because of the ease of bonding with high alignment accuracy due to glass transparency. Anodic bonding offers advantages of high bond strength at moderate or low processing temperature, usually lower than the glass softening point and the melting points of the materials selected for this purpose (200°C-400°C). Thus, materials with large differences in the thermal expansion coefficient can be bonded together with better reliability. The use of low temperature for bonding can avoid damaging pre-fabricated devices and integrated circuitries. In comparison with direct bonding, the surface quality requirement of anodic bonding is much tolerant and does not require an ultraclean environment. In fact, even if a 1~2μm high step exists on the interface, realization of anodic bonding is still possible [START_REF] Jia | Silicon Micromachined Ultrasonic Transducers[END_REF]. Other advantages of glass to silicon bonding are the benefits of the electrical, thermal and optical properties of the glass. The glass provides dielectric isolation properties that significantly reduce parasitic capacitance in capacitance transducers (compared with silicon/silicon bonding). In manufacturing, the glass optical transparency enables backside and internal state visual inspections and optical diagnostics of the sealed microstructures [START_REF] Van Le | Anodic Bonding and the Crack Opening Method Experiment[END_REF] Here, we introduce a new process to the CMUTs technology that uses the anodic-bonding technique to fabricate membranes of silicon on a glass wafer for acoustic applications. This new design, based on silicon-on-insulator SOI anodic bonding technology utilizes a structured glass as well as a highly doped and, therefore, conductive membrane of SOI silicon. This new technology enables the fabrication of large membranes with large gaps, and can expand the frequency span of the CMUTs to a low frequency operation range. Moreover, anodic-bondingbased CMUT fabrication as its wafer-bonding counterpart provides design flexibility in terms of choice of materials and fabrication of membranes with different thickness profiles.

After introducing the basic operating principles of the anodic bonding technique, we will provide a brief description of the design process as an alternate method of CMUTs technology that utilizes a glass wafer as a

10th International Workshop on Micromachined Ultrasonic Transducers, Salerno, Italy June 16-17 2011

cMUT-ROME 2011 frame and SOI silicon as a membrane. The designed cells have a diameter varying from 50µm to 150µm and a cell to cell pitch decreasing from 150µm to 20µm. The element to element pitch is defined to be less than half of the wavelength of the acoustic ultrasound wave for each designed cell. The ratio between the bottom electrode radius and the cavity radius is 80%, which represents a surface fill factor of 64%. The performance of the CMUT transducers is evaluated by demonstrating basic characterization results and their relevance to the device performance in order to determine the yield and the uniformity of the anodic bonding process. It involves laservibrometry as well as electrical impedance measurements

The vibration behaviour of CMUT individual cells and arrays was investigated using a laser vibrometer coupled to a displacement decoder. The impulsional response of each cell was measured (under AC+DC excitations) and its spectral response deduced, allowing us to extract its band width directly as well as its fundamental resonance frequency and first harmonics and compare them with theoretical predictions based on circular membrane models from FEM and analytical simulations. The measured spectral response was confirmed by making a frequency sweep excitation from 30KHz to 20MHz for every tested cell. To eliminate the effect of acousto-optic interaction and calibrate the displacement decoder, the laser beam was initially focused on the posts of the cells where the actual displacement would be negligible in comparison to that of the cells as well as on a static flat polished silicon wafer. Assuming the posts and the flat wafer to be fixed, the interferometer data gathered was later used to correct for the actual membrane displacements [START_REF] Bayram | Characterization of crosscoupling in capacitive micromachined ultrasonic transducers[END_REF]. CMUT devices with different frequency spans were fabricated using this technology, and tested. These devices have been applied to a displacement measurement in air. So far, cMUTs single cells and arrays have been produced with resonance frequencies in the range from 290Khz to 2.4MHz. Peak displacement levels varied approximately between 1nm to 35nm for the 50µm cells and between 5 to 35nm for the 150µm cells. Several over-responsive and under responsive elements could be seen as well as a non-responsive element. The locations of the under-and over-responsive elements were random in nature. The most likely explanation for these measured differences is due to variations in the CMUT membrane characteristics (residual stress, thickness, bonding quality to support).

The visualization, analysis of surface topography (static and dynamic) and structural vibration patterns (dynamic) of the CMUT individual cells/arrays have been investigated using the MS-500 Polytec MEMS analyzer. The resulting vibration maps gave us important informations about the static and dynamic mechanical behavior of cMUTs. In particularly, these helped us in quantifying the anodic bonding quality regarding the resulting curvature (tensile-compressive strain) and the initial deflection of the membrane due to the vacuum bonding.

The CMUT cells and arrays were also electrically characterised in air using a capacitor meter and vector network analyzer in order to obtain the electrical resonance frequency and bandwidth. Input electrical impedance in air was measured as a function of frequency. The electro-mechanical coupling factor evaluation was based on capacitance determination and, the resulting resonance and antiresonance peaks.

Initial investigation results showed that, this newly developed CMUT has the advantages of relatively uniform cavity definition and membrane bonding with reliable vibrational behavior for acoustic applications


J Lardies
email: joseph.lardies@univ-fcomte.fr

STRUCTURAL HEALTH MONITORING OF STAY CABLES BY THE SCRUTON NUMBER

Keywords: Vibrations, Stay cables, Cable tension, Scruton number

Identification of modal parameters: eigenfrequencies, damping ratios and mode shapes from empirical data is of fundamental engineering importance in the dynamical analysis of mechanical structures and particularly in the analysis of stay cables. These modal parameters will serve to perform structural health monitoring, damage detection and safety evaluation. The identification is based on the knowledge of output-only responses, without using excitation information. We propose in this communication a time domain method based on a subspace algorithm to extract the modal parameters of vibrating stay cables from outputonly measurements. We can then use these identified parameters for structural health monitoring by the computation of the Scruton number and the cable tension.

INTRODUCTION

Identification of modal parameters from empirical data is very important in mechanical vibrations and in the dynamical analysis of mechanical structures excited by external forces. An application constitutes the vibration analysis of stay cables in the civil engineering domain. The identified modal parameters will serve to perform structural health monitoring, damage detection and safety evaluation of such structures. The identification is based on the knowledge of output-only responses, without using excitation information, and is known as Operational Modal Analysis (OMA), also named as ambient or natural excitation or output-only modal analysis [START_REF] Peeters | Reference based stochastic subspace identification for out-put only modal analysis[END_REF][START_REF] Lardies | A stochastic realization algorithm with application to modal parameter estimation[END_REF][START_REF] Petsounis | Parametric time-domain method for the identification of vibrating structures-a critical comparison and assessment[END_REF][START_REF] Bodeux | Modal identification and damage detection using the datadriven stochastic and ARMAV methods[END_REF][START_REF] Zhang | Modal parameter identification using response data only[END_REF][START_REF] Zhang | Modal indicators for operational modal identification[END_REF]. The principal advantages of OMA are the ease of use and the fast to conduct since no artificial excitation equipment is required. We propose in this communication a time domain method, based on a subspace algorithm to extract the modal parameters of vibrating stay cables from output-only measurements. We can then use these identified parameters for structural health monitoring.

The subspace identification method is a time domain method and uses the covariance matrices between output data. The subspace method is based on the observability and controllability properties of linear systems and we use shifted properties of the controllability matrix to identify the modal parameters. Unfortunately spurious modes, which are inherent to the subspace method, appear [START_REF] Zhang | Modal indicators for operational modal identification[END_REF][START_REF] Juang | Applied system identification[END_REF]. To eliminate such spurious poles a criterion using modal coherence indicators is used. This criterion describes the coherence between each mode of the state space model and the modes directly inferred from the measured signal. The criterion serves as a model quality measure. Modal parameters are then obtained from stability diagrams. Experimental results of a single tower, double row cable-stayed bridge supported by 112 stay cables are presented, where ambient vibrations of each stay cable are carried out using accelerometers. The subspace method is then applied to determine the eigenfrequencies and damping coefficients of different cables. An estimate of the fundamental frequency f i for each cable is then obtained. We can then compute the tension force for each cable. For each natural frequency identified a respective damping coefficient is assigned. Thus, a Scruton number for each cable and for each mode of vibration is determined [START_REF] Ren | Experimental modal analysis of stay cables in cable-stayed bridges[END_REF][START_REF] Geier | Cable-stayed bridges and their dynamic response[END_REF]. High values of the Scruton number tend to suppress the oscillation and bring up to the start of instability at high wind speeds. For example, a Scruton number superior to 10 is sufficient to prevent rain, wind and traffic induced vibrations. In a continuous monitoring and modal analysis process, the tension forces and Scruton numbers could be used to assess the health of stay cables in cablestayed bridges.

VIBRATION ANALYSIS OF A CABLE AND THE SUBSPACE IDENTIFICATION METHOD

Vibration analysis of a cable

The free vibration equation of cable motion is given by [10]

EI 4 x t) w(x, 4   + ρ S 2 t t) w(x, 2   -P 2 x t) w(x, 2   = = f(x,t) (1)
where E is the Young's modulus, I the moment of inertia of the cable cross-section, ρ the specific mass, S the cross-sectional area, P the axial load, f(x,t) the external load, w(x,t) the transversal displacement which is assumed to be small, x the position along the cable and t the time variable. The solution of the free vibration equation can be obtained using the method of separation variables as w(x,t)=q(t) W(x) (

Assuming q(t) harmonic with natural frequency ω in rad.s -1 and substitution of equation (2) into equation [START_REF] Peeters | Reference based stochastic subspace identification for out-put only modal analysis[END_REF] gives

EI 4 x W 4 d  -P 2 x W 2 d  -ρ S ω 2 W = 0 (3)
By assuming the solution W(x) to be

W(x) = C e s x ( 4
)
in equation ( 3), the auxiliary equation can be obtained

s 4 -I E P s 2 - I E 2 ω A ρ = 0 (5)
and the roots of this equation are

s 1 2 , s 2 2 = I E 2 P  1/2 I E 2 ω A ρ 2 I 2 E 4 2 P            (6)
The solution of W(x) can be expressed as

W(x) = C 1 cosh s 1 x + C 2 sinh s 1 x + C 3 cos s 2 x + C 4 sin s 2 x ( 7
)
where the constants C 1 to C 4 are to be determined from the boundary conditions:

W (0) = 0 ; 2 x d W(0) 2 d = 0 (8) W(L) = 0 ; 2 x d W(L) 2 d = 0 (9)
Equations ( 8) require that C 1 = 0 and C 3 = 0, and so

W(x) = C 2 sinh s 1 x + C 4 sin s 2 x ( 10
)
The application of ( 9) to ( 10) leads to sinh s 1 L sin s 2 L = 0 (11) Since sinhs 1 L > 0 for all values of s 1 L  0 , the only roots to this equation are

s 2 L = k π , k = 0, 1, 2,….. ( 12
) Thus equations ( 12) and ( 6) give the natural frequencies of vibration:

f k = 1/2 I E 2 π 2 L P 2 k 4 k 1/2 S ρ I E 2 2L π                  (13)
Further, note that the smallest Euler buckling load for a simply supported beam is given by P cri = π 2 E I / L 2 . Thus equation (13) can be written as

f k = 1/2 cri P P 2 k 4 k 1/2 S ρ I E 2 2L π                    (14)
We would like to approximate the desired partial differentiate equation solution w(x,t) in a separable form as series expansion of time varying coefficients q i (t) and spatially varying basis functions W i (x) :

w ~(x,t) =   n' 1 i i i (x) W (t) q (15)
The Galerkin procedure can then be used to specify the equations of motion for the coefficients q i (t) and for an n' degrees of freedom system with viscous damping, the motion equation is

M q  + C q  + K q = f (16)
where M, C and K are respectively the mass, damping and stiffness matrices. In [START_REF] Barbieri | Dynamic analysis of transmission line cables[END_REF] Barbieri et al. propose a procedure based on experimental and simulated data to identify damping of transmission line cables. The experimental data are collected through accelerometers and the simulated data are obtained using the finite element method. In this communication, we identify eigenfrequencies and damping ratios of line cables from output data only using the subspace method.

The subspace identification method

The subspace identification method assumes that the dynamic behaviour of a structure excited by ambient forces can be described by a discrete time stochastic state space model [START_REF] Peeters | Reference based stochastic subspace identification for out-put only modal analysis[END_REF][START_REF] Lardies | A stochastic realization algorithm with application to modal parameter estimation[END_REF] :

z k+1 = A z k + w k state equation ( 17
)
y k = C z k + v k observation equation (18)
where z k is the unobserved state vector of dimension n=2n'; y k is the (mx1) vector of observations or measured output vector at discrete time instant k ; w k , v k are white noise terms representing process noise and measurement noise together with the unknown inputs, it is assumed that the excitation effect appears in the disturbances w k and v k , since the system input can not be measured; A is the (nxn) transition matrix describing the dynamics of the system and C is the (mxn) output or observation matrix, translating the internal state of the system into observations. The stochastic identification problem deals with the determination of the transition matrix A using output-only data.

The modal parameters of a vibrating system are obtained by applying the eigenvalue decomposition of the transition matrix

A A =  Λ  -1 (19)
where

() i diag   
, i=1,2,…,n, is the diagonal matrix containing the complex eigenvalues and  contains the eigenvectors of A as columns. The eigenfrequencies f i and damping ratios i  are obtained from the eigenvalues which are complex conjugate pair:

  2 2 1 ln 4 cos 4 2 ii i i i ii f Ar t                          (20)     2 2 2 ln ln 4 cos 2 ii i ii ii ii Ar                         (21)
with Δt the sampling period of analyzed signals.

The mode shapes evaluated at the sensor locations are the columns of the matrix Φ obtained by multiplying the output matrix C with the matrix of eigenvectors  : Φ = C  (22) Our purpose is to determine the transition matrix A in order to indentify the modal parameters of the vibrating system. Define the (mfx1) and (mpx1) future and past data vectors as y + k = [y T k , y T k+1 , . . ., y T k+f-1 ] T and y - k = [y T k , y T k+1 , . . ., y T k-p+1 ] T . The (mfxmp) covariance matrix between the future and the past is given by [START_REF] Lardies | A stochastic realization algorithm with application to modal parameter estimation[END_REF] :

H = E[y + k y - k T ] =                     1 p f 1 f f 1 p 3 2 p 2 1 R . R R . . . . R . R R R . R R (23)
where E denotes the expectation operator and the superscript T the transpose operation. H is the block Hankel matrix formed with the (mxm) individual theoretical auto-covariance matrices

R i = E[y k+i y T k ] = CA i-1 G, with G=E[x k+1 y T k ].
In practice, the auto-covariance matrices are estimated from N data points and are computed by

Ri = N -1  1 - N 1 = k T k i + k y y
; i = 0,1, . . ., p+f and with these estimated auto-covariance we form the block Hankel matrix H. In order to identify the transition matrix A two matrix factorizations of H are employed. The first factorization uses the singular value decomposition (SVD) of

H H = U Σ V T = U Σ 1/2 Σ 1/2 V T (24)
with U T U and V T V identity matrices and Σ a diagonal matrix of singular values. The second factorization of the block Hankel matrix H considers its (mfxn) observability and (nxmp) controllability matrices, O and K, as

H =                 G CA . G CA G CA . . . . G CA . G CA CAG G CA . CAG CG 2 - p + f f 1 - f p 2 1 - p =                 1 - f CA . CA C [G AG. . . A p-1 G] = O K (25)
By identification we obtain the observability matrix

O =                 1 - f CA . CA C (26)
and the controllability matrix

K = [G AG. . . A p-1 G] (27)
The two factorizations of the block Hankel matrix are equated to give

H = U Σ 1/2 Σ 1/2 V T = OK ( 28
)
implying O = U Σ 1/2 and K = Σ 1/2 V T . To determine the transition matrix A we use properties of the controllability matrix (properties of the observability matrix can also be used). We introduce the two following matrices:

 K = [ AG A 2 G.....A p-1 G ] and  K = [ G AG A 2 G.....A p-2 G ] ( 29
)
where  K is the matrix obtained by deleting the first block column of K and  K is the matrix obtained by deleting the last block column K. We obtain then:

 K = A  K or  ) V ( T 1/2 Σ = A  ) V Σ ( T 1/2 (30)
We use properties of shifting columns operators ''  '' and ''  '' : let Φ and Ψ be two (axb) and (bxc) matrices, we have the following properties:

( Φ Ψ )  =( Φ Ψ  ) and ( Φ Ψ )  =( Φ Ψ  ), consequently we have ) V Σ ( T 1/2  = A ) V Σ ( T 1/2  (31)
The transition matrix is

A = ) V Σ ( T 1/2  ) V Σ ( T 1/2  +
(32) where ( ) + represents the pseudo inverse of a matrix, and the eigenvalues of the transition matrix are given by

) ( A  = λ [ ) V ( T  ) V ( T  + ] (33)
This approach constitutes a subspace modal identification method and in the case of noisy data a problem of model order determination occurs: when extracting mechanical (or physical) modes this algorithm can generate spurious modes. For these reasons, the assumed number of modes, or model order, is incremented over a wide range of values and we plot the stability diagram. The stability diagram involves tracking the estimates of eigenfrequencies and damping ratios as a function of model order. As the model order is increased, more and more modal frequencies and damping ratios are estimated, hopefully, the estimates of the physical modal parameters stabilize using a criterion based on the modal coherence of measured and identified modes. Using this criterion we can separate mechanical modes and spurious modes.

APPLICATION OF THE SUBSPACE IDENTIFICATION METHOD

The subspace identification method is applied to the analysis of stay cables of the Jinma cable-stayed bridge (Figure 1), that connects Guangzhou and Zhaoqing in Guangdong Province, China. It is a single tower, double row cable-stayed bridge, supported by 28*4 =112 stay cables. Inputs could evidently not be measured, so only acceleration data from accelerometers in contact with a cable are available. For the ambient vibration measurement of each stay cable an accelerometer was mounted securely to the cables, the sample frequency is 40 Hz and the recording time is 140.8 second, which results in total 5632 data points. Cables 1, 56, 57 and 112 are the longest and cables 28, 29, 84 and 85 are the shortest. A full description of the test can be found in [START_REF] Ren | Experimental modal analysis of stay cables in cable-stayed bridges[END_REF]. Figures 2 and3 show the stabilization diagrams on eigenfrequencies of cable 1 and cable 25 using the subspace method and the modal coherence indicator. These diagrams show remarkable stable eigenfrequencies and from these plots we determine the eigenfrequencies of the two stay cables. The subspace identification procedure is applied to determine the eigenfrequencies of 112 cables of the bridge and we obtain the Figure 4. It can be seen that the fundamental frequencies of both bridge sides are almost identical and the fundamental frequency distribution is symmetric with respect to the single tower. The fundamental frequencies vary between 0.533 Hz for the longest cable and 2.703 Hz for the shortest cable. The cable tension can be estimated by the expression T = 4 ρL 2 f 0 2 , where ρ is the linear density of the cable (ρ = 66.94 kg/m), L is the length of the cable and f 0 the fundamental frequency. The maximum and minimum cable forces for the Jinma bridge are then : T max = 5052 kN (cable number 57), T min = 2490 kN (cable number 84). These cable forces can be considered as reference tensions and used as indicators in the field of health monitoring process.

The instability of a cable is assessed by the Scruton number [START_REF] Geier | Cable-stayed bridges and their dynamic response[END_REF] which is defined for each mode of vibration by ,

2

.

i ci a S D   
, where i  is the damping ratio for each mode, ρ a is the density of the air and is the cable diameter. High values of the Scruton number tend to suppress the oscillation and bring up the start of instability at high wind speeds. Considering ρ a = 1.2 kg/m 3 and D = 0.203m, the Scruton number for each mode of cable 25 is presented in Table 1. In the Scruton number damping is a main factor and most vibration problems in stay cables are subject to low damping values.

CONCLUSION

Cable vibrations induced by wind, rain and traffic are observed in a number of cablestayed bridges and degradation of these cables may have catastrophic effects. The determination of tension and Scruton number of these cables from output only measurements can be used as structural health monitoring indicators. Operational modal analysis applied to the dynamic data of stay cables provides useful information to determine the current condition of stay cables accurately.

Figure 1 :

 1 Figure 1: View and schematic of the Jimma cable-stayed bridge

Figure 2 :

 2 Figure 2: Ambient time response of cable n°1 and its stabilization diagram on eigenfrequencies

Figure 3 :

 3 Figure 3: Ambient time response of cable n°25 and its stabilization diagram on eigenfrequencies

Figure 4 :

 4 Figure 4: Fundamental frequency of each cable on upstream side and on downstream side

Table 1 :

 1 Natural frequencies, damping ratios and Scruton numbers for cable 25

	Modes	1	2	3	4	5	6
	f k (Hz)	1.960	3.938	5.915	7.893	9.856	11.884
	k  (%)	0.314	0.046	0.178	0.078	0.081	0.094
	Sc k	4.250	0.622	2.410	1.056	1.096	1.271


Lukas Capek

Zdenek Lochman

L Dzan

Emmanuelle Jacquet

Biaxial

Biaxial extensometer for measuring of the human skin anisotropy in vivo

I. INTRODUCTION

ESEARCH into the biomechanics of the human skin has attracted a lot attention among dermatologists. Skin performs a wide variety of functions resulting from chemical and physical reactions inside its components. The major function of skin is to act as a barrier to the exterior environment. From the mechanical point of view human skin is a heterogeneous material composed of collagen and elastin fibres in a proteoglycan matrix [START_REF] Robinson | Surgery of the Skin: Procedural Dermatology[END_REF]. Acting together, these components are responsible for the complex mechanical behaviour of skin as manifested in stress, strain curves. These curves are characterized by a low-stiffness region at small strains followed by a dramatic increase in stiffness as the strain becomes large [START_REF] Fung | Biomechanics[END_REF], [START_REF] Agache | Measuring the skin[END_REF]. Moreover, its behaviour depends on direction.

Wijn et al. conducted uniaxial experiments on the leg in parallel and perpendicular directions to Langer´s lines and discovered that the Young´s moduli was 3.8 times greater in parallel direction [START_REF] Pff. Wijn | The alinear viscoelastic properties of the human skin in vivo related to sex and age, Lancaster: Marks and Payne PA Editions Bioengineering and the skin[END_REF]. Khatyr et al. showed that the anisotropy ratio on the forearm was equal to 5 [START_REF] Khatyr | Model of the viscoelastic behaviour of skin in vivo and study of anisotropy‖[END_REF]. Manschot et al. showed that on the leg 76 % of elastic fibres were oriented along Langer´s lines and only 5.1 % were oriented perpendicularly [START_REF] Manschotte | Measurement and modeling of the mechanical properties of the human skin in vivo‖[END_REF].

There are several non-invasive and invasive techniques for the evaluation of mechanical properties of human skin. The most commonly used methods are based on the measurement of suction, torsion and traction [START_REF] Agache | Measuring the skin[END_REF], [START_REF] Evans | Controlled quasi-static testing of human skin in vivo‖[END_REF], [START_REF] Gunner | The mechanical properties of skin in vivo-a portable hand held extensometer‖[END_REF],

Manuscript received August 10, 2010. This work was supported in part by the Regional hospital of Liberec under Grant VR 090303 L. Capek, Department of Applied Mechanics, Technical University of Liberec, Czech Republic (e-mail:lukas.capek@tul.cz)corresponding author.

Z. Lochman, Department of Applied Mechanics, Technical University of Liberec, Czech Republic (e-mail: zdenek.lochman@tul.cz) L. Dzan Department of Maxillo-Facial Surgery, Regional Hospital of Liberec, Czech Republic (e-mail: ladislav.dzan@nemlib.cz).

E. Jacquet, FEMTO-ST, University of Franche-Comte, France (e-mail: emmanuelle.jacquet@univ-fcomte.fr) [START_REF] Vescovo | Validation of measuring technique for identifying of Young´s moduli of human skin[END_REF], [START_REF] Lim | New extensometer to measure in vivo uniaxial mechanical properties of human skin‖[END_REF]. Other experimental methods used to measure the mechanical behavior of the skin elastic wave propagation and indentation [START_REF] Pailler-Matte´i | Analysis of adhesive behaviour of human skin in vivo by an indentation test‖[END_REF]. However, the data obtained by these methods are mainly descriptive and often very different, depending on the experimental conditions. The most common test is the uniaxial one. It works by applying displacements to the skin using two extensible pads that are attached to the skin. The force needed for deformation is measured by a load cell as the skin deforms.

The aim of this study was to construct a prototype of a biaxial extensometer allowing the measurement of the directionally dependence of biomechanical properties of human skin in vivo.

II. MATERIAL AND METHODS

In spite of the fact that there are lot of uniaxial in vivo extensometers, there is no device, to our knowledge, enabling to load the human skin in two directions simultaneously in vivo. The design of the extensometer was done in the CAD software Pro/Engineer (PTC software ltd., Czech Republic). The device consists of two miniature load cells, four actuators, fours pads and a movable frame, Fig. 1.

The input parameters: low mass, high durability, variability and low cost. For optimizing of these parameters an finite element analysis was performed. The movement of two opposite actuators is tied with maximal speed of 600 steps/second and maximal displacement of 26 mm. The movement from the actuators to pads is transfered by a lever system. The pads are joined to the skin surface by doubled

Biaxial Extensometer for Measuring of the Human Skin Anisotropy in Vivo

L. Capek, Z. Lochman., L. Dzan and E. Jacquet R sided adhesive tape. The distance between two opposite pads is adjustable, but the minimum distance is 30 mm. The final design of the real biaxial extensometer can be seen on the Fig. 2. All components were made from duralumin, which meets the low mass with high durability criteria. We have conducted three main test to validate our new biaxial extensometer. All of them were performed at the region of human forearm:  The first one was done in axis x (parallel to Langer´s lines),  the second one in y axis (perpendicular to Langer´s lines)  the last one in x and y axis together.

III. RESULTS

All gained stress and time curves exhibit the initial, lowstiffness region, where the collagen fibres align themselves parallel to the maximum stretch direction while either the elastin and/or the proteoglycan matrix provides resistance to deformation. Once the collagen fibres are aligned, further extension of the skin requires extension of the collagen fibbers, resulting in a significant increase in skin stiffness until failure.

The typical curves are presented below. In the first test carried out in x direction the maximum recorded force was 2.45 N, Fig. 3. In the second test carried out in y direction the maximum recorded force was 6.5 N, Fig. 4. In the last test the maximum recorded force in x direction was 6.6 N and in perpendicular direction 7.5 N, Fig. 5.

IV. CONCLUSIONS

We assume that the tensile tests are, from a mechanical point of view, one of the best possible load applications to identify skin anisotropy. Human skin exhibits nonlinear anisotropic properties. To identify human skin anisotropy biaxial tests must be performed. Nowadays the most common tests used to identify human skin anisotropy are uniaxial tests, which are done separately in different directions. From the mechanical point of view this doesn't provide us with all the information about human skin anisotropy. According to literature, there is no biaxial device enabling tensile tests on human skin in vivo.

The proposed biaxial device is the first prototype of its kind. All gained stress and time curves exhibit the nonlinear behaviour with values close known from the literature [START_REF] Fung | Biomechanics[END_REF], [START_REF] Agache | Measuring the skin[END_REF]. The most important data can be seen Fig. 6, where two stress strain curves were gained in the same time. Human skin anisotropy is well-known property of the human skin. Nevertheless, if on looks at the literature in the field it has not received much attention. It has mostly been dealt with within the context of uniaxial tests.

It should be noted, that only preliminary tests were done and results are not valid from the statistical point of view. In the future the results will be confronted by a correlation system Dantec Q-400 (Dantec Dynamics LtD., Germany) and the statistical results must be gained. The aim of this paper was to present our new biaxial system constructed for measuring human skin in vivo.

Fig. 1 .

 1 Fig. 1. The CAD model of the biaxial system.

Fig. 2 .

 2 Fig. 2. The real biaxial system.

Fig. 3 .

 3 Fig. 3. The force displacement curve from uniaxial in vivo test along the Langer´s lines on the forearm.

Fig. 4 .Fig. 5 .

 45 Fig. 4. The force displacement curve from uniaxial in vivo test perpendicular to the Langer´s lines on the forearm.

ACKNOWLEDGMENTS

We would like to thank Dr. Aleš Lufinka for technical help during this project.