T Guo

-M C Ho

Ba Tho

Biomechanical factors effects on cementless Total Hip Arthoplasty

Keywords: bone material property, bone remodelling, FEA, THA, femur MOTS-CLÉS: Propriété matérielle osseuse, remodelling osseuse, FEA, THA, fémur

The objective of this study is to simulate the dynamic bone adaptation, and to apply this method on patients to predict the effect of different biomechanical factors on cementless Total Hip Arthroplasty (THA) results. The results show that the factors considered and their combination had obvious effects on the bone remodelling results. RÉSUMÉ: L'objectif de cette étude est de simuler l'adaptation dynamique osseuse et d'appliquer cette méthode sur des patients pour prédire l'effet de facteurs biomécaniques sur les arthroplasties de hanche sans ciment. Les résultats montrent que les facteurs considérés et leur combinaison ont des effets évidents sur le remodelage osseux.

Introduction

Bone has the ability to change its structure to adapt the alteration of mechanical environment. It is known as bone remodelling. For the cementless THA, bone remodelling happens after the operation. But for the different individual, long-term operation results are quite different. Some researchers concentrated on mechanical characteristics on bone-prosthesis system, but they did not consider the bone adaptation process [START_REF] Mann | Coulomb frictional interfaces in modeling centered total hip replacements: a more realistic model[END_REF][START_REF] Couteau | Validation of a three dimensional finite element model of a femur with a customized hip implant[END_REF][START_REF] Viceconti | Large-sliding contact elements accurately predict levels of bone-implant micromotion relevant to osseo[END_REF]; some researchers described bone as living material, but suppose the interface between bone and prosthesis can only transfer the common press stress but not shear stress [START_REF] Van Rietbergen | The Mechanism of bone remodeling and resorption around press-fitted titanium stems[END_REF]; some considered both bone remodelling and friction between bone and prosthesis [START_REF] Fernandes | A contact model with ingrowth control for bone remodelling around cementless stems[END_REF], but they applied mathematic optimised method to optimise the bone structure, this method can work, but without physiological meaning behind bone remodelling rule; others simulated bone remodelling process with the consideration of bone material properties from the literature (Bisaskos et al., 2005).

Biomechanical factors are the main controlling factors to effect the structure of bone (Li et al., 2000). we chose four main biomechanical factors, and studied the effect of these factors on the long term THA bone remodelling results.

Our approach was to combine non linear FEA with physiological bone remodelling equations, simulate the bone remodelling on post operation process by applying real patients database.

Materials and methods

The different steps of the methodology are described below:

1) Specific geometry model was built, which was developed from clinical X-ray image, illustrated in figure 1(a), provided by Polyclinique St. Côme.

2) Develop mesh and input material properties in Patran (MSC.Software). Our FEA model includes 565 2D quadrilateral elements. The initial density distribution after surgery (14 days) is illustrated in Figure 1 (a), which means in this simulation, we considered different material properties of cortical bone, spongy bone and layer bone.

3) Boundary conditions, loads were applied and non linear analysis was performed using Marc. Analysis (MSC.Software). We fixed all the nodes at the bottom of the model, muscle forces were applied and four areas of interest are pointed (Figure 1). We will observe bone density distribution of these areas after long-term operation. (Li et al., 2000).

⎟ ⎠ ⎞ ⎜ ⎝ ⎛ - ρ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ⋅ = ρ ∑ = n n i i K / U n B dt d 1 1 [1] 2 0 0 2 1 2 1 n n n E K ε ρ ε σ ρ ⋅ = ⋅ ⋅ = [2]
Here, B is the bone remodelling speed ratio, n ε represents four bone remodelling thresholds (n = 1,2,3,4), K n represents four specific bone physiological status, which correspond with four strain threshold (Figure 2). With these bone physiological bone remodelling equations, bone adaptation behaviour under all kinds of external loads was expressed.

Biomechanical factors considered here include: weight, bone material properties, bone remodelling speed ratio and muscle forces. Each factor was classified into three degrees: a) Weight: for which 3 classes (40kg, 70kg, 100kg) were considered; b) Bone material properties: the values are mentioned in table 1 (Ho Ba Tho et al., 1992). c) Bone remodelling speed ratio, for which 3 cases (0.002, 0.005, 0.008 kgN -1 m -1 kgm -3 s -1 ) were considered; d) Muscle forces: include the forces came from Gluteus maximus, Gluteus medius, Gluteus minimus, psoas, illiacus and poriformis. Each muscle force was divided into 3 classes. The middle class was from the literature [START_REF] Bitsakos | The effect of muscle loading on the simulation of bone remodelling in the proximal femur[END_REF], and the other two classes were calculated on scale according different weight classes.

Functional analyses were allowed to study the effect of the factors simultaneously, three levels for each of four factors, so we had 3 4 kinds of analyses.

Results

In our 81 results, 5 ideal results, 30 bad results, 46 average results were obtained. Almost all the typical clinical operation results are included in our simulated results.

Here we present two typical cases of results: 1) Ideal results represent the best operation results, that means after 14 years, the material properties of femur has nearly no change, bone is nearly the same as the beginning (figure 3); 2) Bad results represent the worst operation results (except fractured prosthesis), that means obvious remodelling or resorption happened in two or more of these four areas (figure 4).

It is found that in 5 ideal results, 4 of them had very good bone properties (low weight and high density) from the beginning or had harmony combination among the biomechanical factors (four factors are in the closed level), all of them had low bone remodelling speed ratio; and in 30 bad results, 27 of them had bad bone properties from the beginning or had improper combination among the biomechanical factors (four factors are not in the closed level); and 13 of them had high bone remodelling speed ratio.

Discussion

Our simulated results have good qualitative agreement with the long-term clinical results. From the comparison between our simulation results and clinical operation results after 14 years, we could find: each of the biomechanical factors had obvious effects to the simulation results, bone material properties are essential to obtain an ideal operation results, muscle force and bone remodelling speed ratio were important, too. Too much value of muscle force and bone remodelling speed ratio would cause abnormal bone remodelling, and proper these factors could help to obtain the ideal results. Not only the biomechanical factors that we considered, but also the combination of these factors are very important to the long-term operation results. Acknowledgements: The financial support of CNRS and Dr. Roux from Polyclinique St. Côme are gratefully acknowledged.

Figure 1 :

 1 Figure 1: The patient's X-Ray of 14 days after the operation (a) and FEA model (b), areas pointed by rows represent four areas of interest.

Figure 2 .

 2 Figure2. Bone Biomechanical adjusting and controlling system(Li et al., 2000).

Figure 3 .

 3 Figure 3. Comparison of the ideal patient's X-ray results 14 years after surgery (a) and visualization of simulated ideal result (b).

Figure 4 .

 4 Figure 4. Comparison of the bad patient's X-ray results 14 years after surgery (a) and visualization of simulated bad result (b).

Table 1 :

 1 3 classes of bone material properties

	;

Giens. Volume X -No. X/2005


Jamil Kanfoud -M A Hamdi
email: jamil.kanfoud@utc.frhamdi@utc.fr

Preliminary study of acoustic protection of Ariane V fairings using layers of porous materials

Keywords: BIOT equations, Porous material, 1D model, Absorption

During lift-off, space launchers are submitted to high level acoustic loads which may damage the payloads. High cost acoustic absorbers have been previously developed by to protect the payload in Ariane V. CNES launched a program in order to reduce the cost by using a layer of porous materials protected by a fabric.

The purpose of this paper is to propose a 1D model of porous layers based on modified BIOT equations. This model has been validated by impedance tube measurements (absorption coefficient) performed by RIETER ) and the ENTPE.

Based on this model, a parametric study of the material mechanical and acoustic parameters has been carried out with respect of aerospace constraints. The study shows the relevance of the resistivity and the mechanical parameters of the skeleton for the absorption. The effect of the mounting conditions has also been studied experimentally and theoretically. Results reveal the importance of the constraints on the absorption mainly at the low frequency range.

INTRODUCTION

The purpose of this paper is to describe wave propagation in elastic porous layers in order to predict sound absorption inside an impedance tube. The theory used here is based on BIOT theory [1] which has been adapted to acoustic absorbers by Allard & Al [2]. studied transmission through structures lined with elastic porous layers; he described porous materials in terms of displacements of both skeleton and interstitial fluid. Here, we use modified BIOT equations described in terms of skeleton displacement and interstitial fluid pressure. 3 waves propagate inside a laterally infinite porous layer, two compression and one shear wave.

The present work is organised as following. In the first section, we give the solution of the modified BIOT equations propagating inside a laterally infinite porous layer excited by a plane wave under oblique angle. Boundary conditions are explored; the coupling between two porous layers and between a porous layer and a fluid cavity are presented. The case of a septum separating layers (typically used in automotive industry) is studied. In the second section, the case of the novel absorber to be used in Ariane V launcher is studied. The absorber consists of a porous layer covered by a fabric; the latter is simulated as a porous layer as well. A good agreement has been obtained between measurements and simulations. A comparative study of the relevant parameters of the porous layer as well as the effect of the mounting conditions has been carried out.

MODIFIED BIOT EQUATIONS

The modified BIOT equations (with e -iωt sign convention) written in terms of skeleton displacement and interstitial fluid pressure are

( ) 2 . s s U p ρ ω σ αφ δ β φ + ∇ - + ∇ = ( ) 0 p (2.1) 2 1 . ( ) . f p p U U R φ φ β α ρ ω ⎛ ⎞ ∇ ∇ - + + ∇ ⎜ ⎟ ⎜ ⎟ ⎝ ⎠ 0 = (2.2)
The equations (2.1) and (2.2) relate respectively the skeleton and interstitial fluid dynamic equilibrium. They look like elasto-dynamic equilibrium equation (2.1) and Helmholtz equation (2.2) coupled by elastic coupling factor α and inertial coupling factor β. The coupling is more complex as the effective mass densities ρ s and ρ f are derived from the physical densities, the porosity φ, the tortuosity τ and the viscous coupling factor b. Hence, the effective mass densities reflect the inertial coupling between the two phases as well as the viscous coupling which is frequency dependant.

WAVE PROPAGATION IN POROUS MEDIA

The problem investigated in this paper involves plane waves propagation inside a finite thickness layer of porous material with infinite lateral dimensions. The explicit variables expressions (pressure, skeleton and fluid displacements, stress…) are given in terms of 3 waves ; 2 compression waves and a shear wave. Each wave is made of a forward and a backward moving waves toward the thickness direction. The six corresponding wave amplitudes are determined by the appropriate boundary conditions.

The skeleton displacement could be split into a rotational and potential displacements U = U l + U r with ∇^U l = 0 and ∇.U r = 0. Substituting respectively the rotational and the potential wave form in the modified BIOT equations give wave numbers by (3.1) for the rotational wave and (3.2) for the potential wave.

2 s r ρ γ ω μ = - (3.1) 4 2 2 4 2 2 ( 2 ) ( ) s f s f f a c b l l R R ρ ρ ω ρ λ μ λ μ ρ ω β α ρ γ γ ⎛ ⎞ + + + + + - + ⎜ ⎟ ⎜ ⎟ ⎝ ⎠ 0 = (3.2)

Boundary conditions -Coupling fluid layer / porous layer

Coupling a porous and a fluid layers requires 4 conditions. The fluid layer pressure and displacement are noted p 0 and U 0 .

Direct coupling 0 2 2 2 2 ( ) ( f z z z z i U U i U U ω ω φ - - =- -) 0 2 p p = 0 2 z p σ -= 2 0 xz τ = Coupling through impervious film 0 2 z z U U = 2 2 2 ( ) f z z i U U ωφ 0 - - = 0 2 z p σ -= 2 0 xz τ =
The direct coupling happens when no facings separate the fluid and the porous layers. The 1 st equation is cinematic, it expresses the fluid transfer between the porous layer and the fluid layer. Hence, we write the continuity of the air flow on the interface where φ is the surface porosity. In fact, surface porosity is immeasurable, and volume porosity is always used instead for simulations. The 2

( f n n q i U U ωφ = - -)
nd condition expresses the pressure continuity between both layers as the porous material is open. The 3 rd and 4 th equations are dynamic, they express the stress continuity between layers. When an impervious weightless film separates the layers, there are two cinematic conditions. The 1 st equation expresses skeleton and external fluid normal displacements continuity. The 2 nd equation relates air flow nullity as the porous layer is henceforth closed. The dynamic coupling is the same as the direct coupling.

-Coupling porous layer / porous layer

Coupling two porous layers requires 6 conditions. We assume here that there is no relative sliding between the porous layers.

Direct coupling

1 2 z z U U = 1 2 x x U U = 1 2 z z σ σ = 1 2 xz x z τ τ = 1 2 p p = 1 1 1 2 2 2 ( ) ( f f z z z z i U U i U U ωφ ωφ - - =- -)
Coupling through impervious film

1 2 z z U U = 1 2 x x U U = 1 2 z z σ σ = 1 2 xz x z τ τ = 1 1 1 ( ) f z z i U U ωφ 0 - - = 2 2 2 ( ) f z z i U U ωφ 0 - - =
We notice here that 4 equations are the same for direct and film coupling. These equations are similar to the coupling between elastic layers. 1 st and 2 nd equations express cinematic coupling of skeleton displacements. 3 rd and 4 th equations express dynamic coupling of the total stresses (stress on skeleton and interstitial fluid). The two latter equations differ with the considered configuration. For direct coupling, we have pressure and air flow conservation across interface. For coupling through impervious interface, air flow vanishes for both layers.

NORMAL INCIDENCE ABSORPTION COEFFICIENT

We compared measurement and simulations of a 10cm porous layer made of foam. This layer will be integrated in Ariane V fairings as an acoustic protection as a replacement of the high cost fap. Measurements of the absorption coefficient of the naked foam and the foam covered by a fabric have been carried out in a 10cm diameter impedance tube . Results are obtained in the frequency range from 100Hz to 4000Hz. Figure .1 shows a good agreement between measurements and simulations, particularly with the naked foam. The fabric, which aim is the protection of the payload from foam dust, affects highly the absorption coefficient. The absorption coefficient is lowered at the high frequency range (above 500 Hz). Below, fabric enhances the absorption efficiency of the porous complex.

Porous material behaviour is very sensitive to the boundary conditions. In this paper we investigate the sensitivity of the porous material to the lateral and normal boundary conditions. For the normal mounting conditions, two configurations are studied : porous material stuck to the background (S) or laid against the background (L). The former case (S) is simulated by imposing null normal displacement of the porous structure and fluid. The latter (S) is simulated by adding a thin fluid layer between the porous layer and the background (1mm thickness). The main difference between the two configurations is that the frequency of the structure resonance is displaced due to difference in limit conditions. Moreover, the structure resonance is less pronounced with the L configuration. As expected, toward the high frequency range the two configurations tend to converge (Figure 2.1). At the lower frequency range, no clear tendency emerges. For example, with the simulations carried out, the L configuration is better from 150 to 500 Hz. Below, the S configuration ensures better absorption. Lateral boundary conditions have been studied experimentally from 40Hz to 200Hz range. A set of measurement of the foam covered by a fabric has been carried out in a 46 mm diameter impedance tube. Figure 2.2 shows that increasing constraints, at low frequencies (below 250Hz), don't have steady effect. Indeed, with growing frequency, more constraints result first in increasing absorption, then in decreasing it. The combinations of analytical and experimental results stress the importance of boundary conditions on the porous material behaviour mainly at the low frequency range.

A parametric study of the porous materials parameters has been carried out. It shows that the main influencing parameter of the porous material is its resistivity. At high frequency range, a low resistivity enhances absorption. At the low frequency range, no steady effect of the resistivity is observed. Hence, at a given frequency, there's an optimum value of resistivity which gives a maximum of absorption. It is worth noting that the optimum resistivity value depends closely of the mechanical parameters of the porous material and the mounting conditions.

CONCLUSION

In this paper, the BIOT theory has been used to predict the absorption coefficient of a complex absorber made of a porous foam covered by a high resistif fabric. The comparison of the theory and the measurements for the porous complex shows excellent correlation, particularly with naked foam. Through simulations and measurements, it has been proved that mounting conditions (lateral and normal) are critical with regard to the porous material absorption efficiency. It has also been shown, by parametric studies, that resistivity is the critical acoustic parameter for absorption. At high frequencies, growing resistivity enhances absorption. At low frequencies, a critical value of resistivity gives optimum absorption. Besides, low frequency absorption coefficient is very sensitive to the elastic parameters of the skeleton. Since the model is semi-analytical, an optimisation has been carried out and showed that a modification of the material critical parameters enhances considerably its absorption efficiency (up to 25% at the 63Hz 1/3octave ).

  Fig 1. confrontation simulation/measurements (S configuration)


Marcela Cruchaga
email: mcruchag@lauca.usach.clr

Diego Celentano

Piotr Breitkopf

Alain Rassineux

Pierre Villon

Remaillage des fronts lagrangiens dans les écoulements bi-fluides

Keywords: MOTS-CLÉS :Interfaces fluide-fluide, remaillage, Approximation Diffuse Fluid-fluid interfaces, remeshing, Diffuse Approximation

La technique MLIT (Moving Lagrangian Interface Technique) utilise un maillage fixe pour la résolution du problème de dynamique des fluides, l'interface entre les deux fluides étant modélisée par un ensemble de marqueurs mobiles. Toutefois, comme un schéma lagrangien est employé pour avancer la position de l'interface, la distribution des marqueurs varie au cours du calcul et peut nécessiter des remaillages. Nous proposons ici une technique de placement des nouveaux marqueurs, basée sur l'approximation diffuse de la courbure. La technique proposée intègre la contrainte de préservation du volume et tient compte des restrictions géométriques. Cette approche est illustrée par un exemple de mouvement propre d'un liquide dans un réservoir.

ABSTRACT. In the framework of fixed-mesh finite element approaches for the computation of the fluid dynamics of the problem, a moving lagrangian interface technique (MLIT) was proposed to describe the interface. Nevertheless, as a lagrangian scheme is used to move the interface, its updated position could be represented by a highly distorted distribution of markers. To avoid this problem, we propose in this work a new remeshing technique applied to the interface. The remeshing is performed according to a curvature-based criterion by using a diffuse approximation technique. Moreover, the new distribution of points needs to be volume preserving. This condition is included as a constraint in the remeshing algorithm and it is the novel contribution of this work. Geometric restrictions, size distribution quality and physical aspects are all considered in the proposed formulation. Particular features of their capabilities are evaluated in simple test: a sloshing problem.

Introduction

Many engineering applications require the analysis of time-dependent incompressible flows with moving two-liquid interfaces or free surfaces. Filling processes or open channels flows give an idea of the wide range of problems involving moving boundaries and interfaces. Many investigations have been devoted to numerical developments of formulations able to deal with the treatment of this kind of problems. Either moving or fixed domain discretizations are adopted to describe problems involving interfaces. All of them provide alternative solutions to update the material front and to overcome severe difficulties present in the numerical simulations of moving interfaces. In addition to the inconveniences related to the numerical schemes used to solve the fluid dynamics equations usually coming from the incompressibility constraint and relevant convection terms, there are many others due to the presence of two fluids and to the algorithm used to move the interface. In particular, global mass preserving and discontinuity in material properties are two of the most important aspects to be properly described in order to avoid numerical oscillations and to produce accurate results.

The mass preserving algorithm for the remeshed interface

A mass conservation algorithm applied to the remeshed interface is presented in this section. A remeshed interface position given by (x rem ,y rem ) is not, in general, a mass preserving configuration. Then, a further correction in the front coordinates, denoted as ( vp x rem , vp y rem ), needs to be performed in order to attain a known given volume Ω 1 (Figure 1 The coordinates of each node j are (x rem (j),y rem (j)) for j=1,…,n rem , where n rem is the number of nodes. In these conditions, an 1D isoparametric representation of the interface within the interval [j,j+1] is given by: , y

) ( N ) ( y 1,2 k k rem k rem ∑ = = θ θ k rem k rem k 1,2 x ( ) N ( )x θ θ = = ∑ (1)
where -1≤θ≤1 the isoparametric coordinate associated to the shape function

matrix N ( ) 1 ( 5 . 0 N 1 θ - = and ) 1 ( 5 . 0 N 2 θ + =
) such that x 1 rem = x rem (j), x 2 rem = x rem (j+1) and y 1 rem = y rem (j), y 2 rem = y rem (j+1). The computation of the grey area in Figure 1.a is performed through the Green-Riemann formula rem

1 n 1 e 1 1 T rem , rem rem QY X d ) ( x ) ( y Area rem ∑ ∫ - = - = = θ θ θ (2)
where , and the matrix Q is computed as the assembly of the (n

)) n ( x ,
θ d N N Q 1 1 ' j i e ij ∫ - = (3)
resulting, for the adopted linear shape functions, in the following expression:

⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ - - - - = 2 / 1 2 /
0 0 ... 2 / 1 0 2 / 1 0 0 0 ... 0 2 / 1 2 / 1 Q (4)
of a tridiagonal matrix whose components are independent of the distance between markers.

The volume to preserve can be written, according to equation (2), as:

1 rem vp T rem vp Y Q X Ω = (5)
considering that the volume preserving new coordinates are given by: δX new and δY new being the admissible variations of the coordinates and β is a scalar parameter. These variations are assumed to be described by:

[ ] t u(j) (j) (j) (j) x rem δ δ n V ⋅ - = sign , [ ] t v(j) (j) (j) (j) y rem δ δ n V ⋅ - = sign (7)
where V=(u,v). The expressions (7) are proposed on the basis that the admissible motion of the markers is in the direction of their velocities and, besides, they must expected to take place. The computations are performed with two different time steps of sizes 0.5 s and 0.1 s.

The time history of the relative wave height at the left side of the container (i.e., (h A -0.3)/0.3, h A being the instantaneous wave height at point A) obtained using different numerical approaches is presented in Figure 3. The results computed with both methods using the finer mesh and the smaller time step are practically identical. These results also adjust to those computed in [7,8] that are considered as the reference ones. Moreover, when the coarse mesh and the smaller time step are considered, the MLIRT provides slightly better results than those corresponding to ETILT. Nevertheless, erroneous results are obtained when using the ETILT with greater time step for both meshes. On the other hand, the relative wave height at point A predicted by the MLIRT exhibits a good response also for the greater time step. This satisfactory behavior in relatively coarse both time and space discretizations is an attractive aspect of the proposed methodology. Further investigation on this subject needs to be done in more severe problems.

Conclusions

In this work, a volume preserving remeshing technique based on Diffuse Approximation of curvatures aimed at redefining the markers' distribution on the interface has been described. The remeshing scheme is embedded in the two-fluid dynamic computation. The numerical predictions of the interface position during time have been found to be nearly independent of the number of initial markers and their distribution.

Figure 1 .

 1 Figure 1. Mass preserving algorithm: a) interface description, b) move directions.

Figure 2 .Figure 3 .Figure 4 .

 234 Figure 2. Two-fluid interface problem: a) problem description and meshes used in the analysis: b) fineand c) coarse (respectively identified as "fm" and "cm").

Acknowledgements

The support provided by the Chilean Council for Research and Technology CONICYT (FONDECYT Project Nº 1020029), the Cooperation Project CONICYT -CNRS, the Department of Technological and Scientific Research at the Universidad de Santiago de Chile and the Université de Téchnologie de Compiègne (Visitor Professor Program) are all gratefully acknowledged.

/2 ⎤ ⎦ 2 compensate the loss (or gain) of mass flow through the interface. It is important to remark that with this procedure a physically-based correction is attained. Then, by introducing equations ( 6) and ( 7) into (5), we have the following quadratic equation on the β parameter:

where the values of a, b and c are:

The solution of equation ( 8) results in two real roots (the variations (7) were found to always guarantee the condition b 2 -4ac ≥ 0). In practice, β is chosen as the root with a minimum absolute value. With this parameter, equation (6) gives a mass preserving interface position. This approach provides a remeshed interface position with zero error in the volume. No iterations are needed to achieve convergence. This mass conservation corrector algorithm can be applied between two different configurations with different number of markers. The technique can be also extended to 3D where in this case the expression equivalent to (8) is of third degree on β.

Numerical example: Two-liquid interface problem

A closed container (Figure 2.a) with dimensions 0.8 m x 0.6 m is filled with two fluids with the lighter one being on top of the heavier one. The initial, inclined interface is linear with an average height of 0.3 m and slope 1:4. The fluid properties [START_REF] Tezduyar | A new strategy for finite element computations involving moving boundaries and interfaces -the deforming-spatial-domain/space-time procedure, part I and II[END_REF] used in the present analysis are ρ 1 = 2 kg/m 3 , ρ 2 = 1 kg/m 3 , µ 1 = µ 2 = 0.001 kg/ms with g=0.294 m/s 2 (the effect of different set of properties has been assessed in [49]). The boundary conditions consist of no-slip conditions at the upper and lower walls, and slip conditions at the side walls. Both fluids are initially at rest. Two different meshes are used (Figures 2.b and 2.c) to test the performance of the improved formulation when, in particular, coarse meshes are used. Note that the meshes in Figure 2 are refined in the zone where the movement of the interface is


Jérôme Fehrenbach

Davuth Ia

Mohamed Masmoudi

Rémi Souchon

Philippe Trompette

Medical imaging using elastography

Keywords: élastographie, imagerie médicale, optimisation topologique, problème inverse, élasticité, échographie elastography, medical imaging, topological optimization, inverse problems, elasticity, echography

Le tenseur de Hooke d'un matériau élastique (linéaire isotrope) dépend de deux paramètres : le module d'Young ¤ et le coefficient de Poisson ¥ . Un problème intéressant -du point de vue médical -est de détecter des inhomogénéités où le module d'Young prend une valeur différente du reste.

On propose ici une méthode pour reconstruire une approximation du module d'Young, en supposant qu'il est constant sauf dans certaines régions où il prend une valeur différente. L'outil principal est une méthode générale pour les problèmes inverses : il s'agit d'une implémentation économe en mémoire et en calculs de la méthode de Gauss-Newton, basée sur l'utilisation de la dérivation algorithmique en mode direct et en mode adjoint. Cette méthode est validée par des résultats expérimentaux : la propriété de régularisation de l'algorithme de Gauss-Newton permet de localiser les plus grosses inhomogénéités.

Introduction

Prostate and breast tumors can have a Young modulus much higher than the surrounding safe tissue. Their detection by clinical palpation requires that the hard nodulus have to be big or near enough from the skin. Elastography was introduced by J. Ophir et al. [OPH 91] to provide an elasticity image by comparison of two sonograms before and after a small external compression. The main drawback of this technique is that elastograms are noisy and difficult to read. A recent review article on elastography is [OPH 02].

The inverse problem framework [KIR 96] is used to improve elastograms. Gauss Newton algorithm is a regularization method for inverse problems, it chooses a solution of minimal norm. Indeed, it is shown here that the Gauss Newton method can be implemented thanks to both forward and reverse modes of algorithmic differentiation. This implementation of the Gauss Newton algorithm adds data from the physical model to go through the lack of measured data, it has the advantage of being no more expensive than the gradient method and gives a better convergence.

In this paper, the direct problem for elasticity is recalled (section 2). We explain then (section 3) the strategy of the inverse problem : reconstructing the Young modulus from the radial displacement under known boundary conditions. The implementation of Gauss-Newton algorithm using forward and reverse modes of algorithmic differentiation is then splitted in short algorithms. In the last section, experimental results are given.

Direct problem

Consider a smooth domain ¦ in the plane. The boundary of ¦ is divided into two parts : § ¨¦ © . The domain ¦ is filled with an elastic material, subject to a displacement ! " $# &% '¦ )( 10 # 32 . There are no volume forces, there is a known displa- cement 4 on and known forces 5 on . The Lamé coefficients of the material 6 © 6 % 87 2 © 6 @9 BA % 87 2 and C D© EC F% 87 2 © EC 9 3A % G7 2 depend on the space variable (this is equivalent to assuming the Poisson ratio H is constant, and the Young modulus A % G7 2 depends on the space variable).

The boundary values problem is to find a displacement field and a stress field I defined in ¦ by :

% GP RQ 2 S T T U T T V I © 6 XW 1Y a`% 8 2 cb ed gf C `% G 2 ¦ h pi Xq sr I © t ¦ I vu w © 5 © 4
Perform a degree one finite elements method : the finite elements version of the mixed boundary problem % GP RQ 2 is equivalent to the linear system :

x y © (

[1]
where y ! 0 # is the displacement vector, g! 0 # 1 is defined by the boundary conditions and

x

is the stiffness matrix.

The stiffness matrix

x depends on the Young modulus distribution A , this matrix is denoted

x % A 2 . Note that A x % A 2 is an affine map. Let x be its linear part.

Inverse problem

We look for inhomogeneities having Young modulus A # in a material having Young modulus

A )

. The known measure is the radial displacement under known boundary conditions, and the unknown is the location of the inhomogeneities.

The radial displacement at a point is given by ¨ © " pu ( where " is a linear operator. The w D f w matrix of the discretized version of this operator can be easily formed, this matrix is also denoted " .

We want to minimize the following quantity :

d % A 2 © fe f Fg " u vh h ji lk nm o Gp rq ns g # © fe f Fg ut % A 2 g # ( where t % A 2 © " u h h i lk nm 1o 8p rq ns , with h solution of x % A 2 h © v .
The idea is to apply one step of the Gauss-Newton algorithm, starting from an homogeneous distribution The vector 4 will be searched by the conjugate gradient method. Indeed, the conjugate gradient method needs only to know the product of the given matrix by a vector, avoiding thus to compute the whole jacobian matrix, a matrix-matrix multiplication and the use of memory to store a (non sparse) matrix.

A xw yA

For the computation of z t { z t 7 we use the algorithmic differentiation rules, and proceed as follows :

1) Computing © z t 7

It is the directional derivative of a vectorial function : use forward differentiation.

2 x % A 2 © h x % 87 2 h , compute © " pu (see proposition 1), solve for h the equation x { h © " { ~, the components % G 2 of are given by © y%

x % ' B 2 h h 2 (see pro- position f ).

Experimental results

We have applied the algorithm above to experimental data obtained by R.Souchon [SOU 04] : inside a phantom of gelatin including six inhomogeneities (of different sizes), a ball is inflated (see fig. 1) : the measures are the radial displacements between before and after inflation. The imaging system is based on an ultrasound scanner (Combison 311, Kretz, Austria) equipped with a transrectal probe (IRW 77AK, Kretz, Austria).

The data are treated as follows :

1) the boundary conditions are not exactly known, but reconstructed using an efficient implementation of Gauss-Newton algorithm similar to the one explained above (but simpler) : we look for Dirichlet conditions on the inner boundary, and Neumann conditions on the outer boundary.

In order to show that the gradient of the cost functional gives no relevant information, we show in fig. 2 the value of the gradient. The inhomogeneities are located inside the gelatin, not on the boundary.

2) once the boundary conditions are estimated, algorithm 1 is applied to estimate the Young modulus distribution. Indeed, a Tikhonov regularization term is added to the cost functional : F%

A 2 © # t % A 2 s # d D A s # .
The results are shown in fig. 3. We can distinguish four among the inclusions, the other two inclusions are smaller and do not appear on our results.

Figure 1 .

 1 Figure 1. The phantom experiment Figure 2. The gradient

Figure 3 .

 3 Figure 3. Estimation of the contrast


Bernard Rousselet

Gérard

G Vanderborck

U N S A Laboratoire

J A Dieudonné

U M R C N R S 6621

Non destructive testing using non linear vibroacoustic

Keywords: Non destructive testing, non linear vibroacoustic MOTS-CLÉS : Essais non destructifs, vibroacoustique non linà c aire

Some previous papers

-In Ekimov-Didenkulov-Kasakov (1999), [A.E 99],

-In Zaitsev-Sas (1999), [V.Z 99],

-Other results may be found in Sedunov-Tsionsky-Donskoy(2002) [D.D 02],Sutin-Donskoy (1998), [A.S 98], Moussatov-Castagnede-Gusev(2002), [GUS 02] ...

-GDR 2501 (Etude de la propagation ultrasonore en milieux inhomogènes en vue du controle non destructif) -In Vanderborck-Lagier-Groby (2003) [P.D 03], "a vibro-acoustic method, based on frequency modulation, is developed in order to detect defects on aluminum and concrete beams.

Very recent experiments

-have been performed on a real bridge by G. Vanderborck with four prestressed cables: two undamaged cables, a damaged cable and a safe one but damaged at the anchor; -See Vanderborck-Lagier(2004) [VAN ] for a presentation of the results of the experiment with a new post processing graphic presentation of experimental results.

Difficulties of the experiments:

-non linearities of the shakers (including piezoelectric actuators) -Natural non linearities: supports, links of complex multi structures as air planes, bridges etc

Simplest mechanical example

in which we can exhibit intermodulations. We consider a 1 d.o.f example of a spring mass system with a non linear spring. ¡ ¢ ¤£ ¦¥ ¨ § ¨© ¢ ¥ ¢ "! $# &% (' 0) 21 43 65 with § ¨© £ 7

(2.1) with initial conditions ¢ ) (8 95 £ A@ CB ED GF¢ H£ AI PB (2.2)

We are going to solve this equation symbolically with an asymptotic expansion with respect to : ¢ ¤£ ¢ B Q © ¢ © ¢ SR &R TR ; then numerically...

The linear case

The first term is solution of:

¡ ¢ B £ ¥ ¨ § ¨© ¢ B 4 !
# &% (' ) 21 U3 65 with ¢ ) (8 95 £ @ CB ED GF¢ ¤£ I VB which gives:

(2.3) with W £ ! ) 21 © ¥ X § © 5 D ¢ B £ @ B 6Y T`# 9) § 3 65 ) I VB § 1 § W 5 A# &% (' ) § 3 65 ¥ W # &% (' 0) 21 43 65 (2.4)
Remarque 1. If we set @ B £ 8 D aI B £ 8 , then the term of pulsation § has magnitude bc d times the magnitude of the term of pulsation 1 ; this is not a good choice for the non linear case in which bc d is of order e 8 f8 ; it seems good choice @ CB £ 8 D I PB £ g § ) ¥ bc d e 5 h i qp b Pr ts c r d Eu £ ) ¥ 1 § 5 W

Other terms

The term ¢ © is zero but the third term satisfies:

¡ ¢ £ ¥ ¨ § v© ¢ B ¢ B (2.5)
to simplify, we assume @ CB £ 8 and set W £ h i qp b Pr ts c r d &u D Aw £ yx c d b Pc 9d so that:

¢ B £ ) w # &% ' 0) § 3 65 ¥ W # &% (' 0) 21 43 65 A5 £ (2.6) ¥ e E w ¨ # &% (' ) ( § 3 A5 w ) e t w © W © 5 A# &% (' ) § 3 65 (2.7) ¥ t w W © # &% (' 0) A) § 1 05 3 65 ¥ w W © # &% (' ) A) § ¥ 1 5 3 65 (2.8) t w © W # &% (' 0) A) 21 § 5 3 65 w © W # &% ' 0) A) 21 ¥ § 5 3 A5 (2.9) ¥ f ) w © W e t W 5 A# &% (' ) 21 U3 65 e E 9W # &% (' 0) ( 1 43 A5 (2.10)
If we go on in the expansion, we get terms of pulsation 1 § D 1 § etc General tendency:

-The pick of § is much larger than the pick in 1 ¤ § which are the most natural picks in the experiments; -it is delicate to find datas such that the secondary picks at 1 § actually appear when the differential equation is solved numerically.

-Question: algorithm and software for detecting the secondary picks? -then find (by optimization) datas such that the secondary picks are important: criteria for damage.

Transverse vibrations: vibrating masses on streched cables in large displacement

Work performed with Theissen (doctoral student of U. Muenster); Erasmus students N. Goris and I. Altrogge worked on this topic during their stay in UNSA (2004UNSA ( -2005)). We consider n masses attached to horizontal springs (or cables) which are in tension D at rest ; the tension is positive when the cable is in traction which is assumed; at rest the mass is submited to the force the masses are moving (vertically) transversely to the springs; we denote by uper case letters quantities in the rest position and lower case in the current configuration.

Masses in vertical displacement

Here we assume that the masses can move only vertically.

A possible damage of a cable

is breakage of several fibers, this will cause decrease of rigidity B say for cable 1. -Let us start with undamaged cables of same rigidity . If we note , the common length of the unstressed cables, and L their common stressed lenth, their tension is £ ) ¥ 5 ; -now, after damage, B d £ © , cable 1 becomes longer and cable 2 shorter, B fe © , the tension goes down to A £ B ) B ¥ 5 £ © ) © ¥ 5 ; -note the limit case of cable 1 broken is B £ 8 so that the cable 2 gets lenth but the system is no longer working properly! -Before such a breakdown, if the change of tension is substantial, this causes a substantial change of the fundamental frequency; indeed, this is the routine monitoring of cable bridges! -The nonlinear vibroacoustic testing aims at monitoring the cables before such a substantial change.

A non linear string model

A model of non linear string has been introduced first by Kirchoff in 1877 and rederived by Carrier in 1945.

¢ 9g hg 4¥ ) i j ¢ C© k 5 ¢ k Ek £ ml (4.1)
Several mathematical studies of this type of equations have been performed recently (Medeiros(1994), Clark-Lima (1997).

n n n n n n n n n n n n n n n n o o o o o o o o o o o o
p p p p p p p p p p p p p p p p p p p p q q q q q q q q q q q q q q q r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s

L 1 L 3 l 1 X Y 2 L m 2 θ 3 m 1 θ 2 l 2 l 3
θ1

Two masses on stretched cables (cable 2 damaged) moving freely (5.2)

Masses free to move in a plane

Bar models with defects

Bar models with longitudinal waves (dynamical traction and compression) are considered.

0 © 3 © ¥ ' t £ ml ) t D 3 A5 (6.1)
With a non linear stress-strain law:

' £ ) ( t A a f ( ) t 5 5

(6.2)
We could as well assume that the applied load is of order epsilon without any assumption on the nonlinearity. Assuming to be small an approximate solution is searched for with the following "ansatz": £ A B $R &R &R Theoretical justification of the expansions: the situation is complex in full generality: non linear hyperbolic equations exhibit a singularity after a finite time! But: the experiments are performed during a short time interval and the Fourier transforms are computed on these time intervals! Following a suggestion of Guy Metivier we are addressing the problem during a small initial time interval in which the solution is smooth: plan to use an approximation of the equation with a fixed point method proposed in Majda. In any case we should smooth the characteristic function (the material is changing smoothly)!

Conclusion

-Some simple models governed by ODE pr PDE show intermodulations;

-But what is the relative level of secondary peaks for a given set of datas deserves investigations: indeed it is also the difficulty of the real experiments -Need to include other behaviors: shocks, friction -Need of more precise models: non linear beams including tractional, flexural, torsional effects -Mixture of local models for the defect and global models for the undamaged structure to obtain precise results at low computational cost.

  Equation of the dynamics (t u (v w£ ¥ Y T`# 9) (x 5 zy B Y T`# 9) (x zy B 5 l t % £ e R TR &R ' (5.1) { 2¢ f v £ ¥ # &% (' ) (x 5 |y B # &% (' 0) (x zy B 5 ~} % £ e R &R &R '