D Ambard

F Cherblanc

Mechanical behavior of annulus fibrosus: a microstructural model of fibers reorientation

Keywords: intervertebral disc, annulus fibrosus, fibers reorientation, microstructural model, mechanical behavior

Experimental uniaxial tensile tests have been carried out on annulus tissue samples harvested on pig and lamb lumbar intervertebral discs. When subjecting the samples to loading cycles, the stress-strain curves exhibit strong non-linearities and hysteresis. This particular behavior results from the anisotropic microstructure of annulus tissue composed of woven oriented collagen fibers embedded in the extracellular matrix. During uniaxial tension, the collagen fibers reorient toward the loading direction increasing its global stiffness. To describe this behavior, we propose an heuristic two-dimensional rheological model based on three mechanical and one geometrical characteristics. The latter one is the fibers orientation angle becoming the key parameter that govern the macroscopic mechanical behavior. The experimental results are used to identify the physical properties associated with the rheological model, leading to an accurate representation of the stress-strain curve over a complete loading cycle. In this framework, the fibers reorientation can solely account for the rigidity increase while the hysteresis is associated with liquid viscous flows through the matrix. Based on this representation, unusual coupling effects between strains and fluid flows can be observed, that would significantly affect the cell nutrients transport mechanisms.

Introduction

The intervertebral disc is a highly-specialized element of the spine that provides flexibility and adsorbing capacities. When mechanical loads are transmitted along the spine, the intervertebral disc mainly supports compression and flexion stresses. This results in a hydrostatic excessive pressure in the nucleus pulposus and generates circumferential tensile stresses in the surrounding annulus fibrosus. Maximal physiologic strains observed along the annulus superficial tissue are about 4% when the disc is in compression or torsion, and about 6% when the disc is in flexion or extension [START_REF] Bass | Biaxial testing of human annulus fibrosus and its implications for a constitutive formulation[END_REF][START_REF] Costi | Direct measurement of intervertebral disc maximum shear strain in six degrees of freedom: Motions that place disc tissue at risk of injury[END_REF] , but never exceed 10% [START_REF] Shirazi-Adl | Strain in fibers of a lumbar disc. analysis of the role of lifting in producing disc prolapse[END_REF][START_REF] Costi | Direct measurement of intervertebral disc maximum shear strain in six degrees of freedom: Motions that place disc tissue at risk of injury[END_REF] . To hold these large circumferential strains, the annulus tissue is composed of a woven oriented structure of collagen fibers embedded in a highly hydrated matrix.

This structured layout provides some interesting mechanical properties, indeed, an anisotropic and non-linear behavior is observed when subjected to uniaxial tensile tests [START_REF] Fujita | Radial tensile properties of the lumbar annulus fibrosus are site and degeneration dependent[END_REF][START_REF] Klisch | Application of a fiber-reinforced continuum theory to multiple deformations of the annulus fibrosus[END_REF][START_REF] Elliott | Anisotropic and inhomogeneous tensile behavior of the human anulus fibrosus: experimental measurement and material model predictions[END_REF][START_REF] Tower | Fiber alignment imaging during mechanical testing of soft tissues[END_REF][START_REF] Wagner | Theoretical model and experimental results for the nonlinear elastic behavior of human annulus fibrosus[END_REF][START_REF] Holzapfel | Single lamellar mechanics of the human lumbar annulus fibrosus[END_REF][START_REF] Iatridis | Mechanical damage to the intervertebral disc annulus fibrosus subjected to tensile loading[END_REF][START_REF] Guerin | Quantifying the Contributions of Structure to Annulus Fibrosus Mechanical Function Using a Nonlinear, Anisotropic, Hyperelastic Model[END_REF][START_REF] Smith | Elastic Fibers Enhance the Mechanical Integrity of the Human Lumbar Anulus Fibrosus in the Radial Direction[END_REF] .

Using fiber-reinforced continuum theory, hyperelastic formulations have been developed to account for this non-linear behavior. This leads to some exponential strain energy functions able to curvefit the experimental stress-strain behaviors from various configurations: tension or compression on axial, radial or circumferential samples [START_REF] Klisch | Application of a fiber-reinforced continuum theory to multiple deformations of the annulus fibrosus[END_REF] . Some attempts to investigate the contributions of individual structural components (matrix, fibers, fiber-matrix or fiber-fiber interactions) have been done [START_REF] Guerin | Quantifying the Contributions of Structure to Annulus Fibrosus Mechanical Function Using a Nonlinear, Anisotropic, Hyperelastic Model[END_REF] . Nevertheless, with these homogenized approaches, the geometrical modifications of the micro-structure are hidden in the strain energy functions, what make tough to identify the micro-structural origin of the annulus fibrosus complex behavior. Actually, when subjected to uniaxial tension, the collagen fibers reorient toward the loading direction. With techniques of digital image processing, this reorientation has been quantified showing significant geometrical modifications of the micro-structure [START_REF] Tower | Fiber alignment imaging during mechanical testing of soft tissues[END_REF][START_REF] Guerin | Degeneration affects the fiber reorientation of human annulus fibrosus under tensile load[END_REF] . These experimental observations imply that the fibers orientation variations could be large enough to induce a macroscopic stiffness increase. This point of view is confirmed by biaxial tests that point out higher stiffnesses but less non-linearities [START_REF] Bass | Biaxial testing of human annulus fibrosus and its implications for a constitutive formulation[END_REF][START_REF] Bruehlmann | In situ intercellular mechanics of the bovine outer annulus fibrosus subjected to biaxial strains[END_REF] .

Therefore, the fibers orientation could be the key parameter that governs the mechanical behavior of annulus and this geometrical characteristic could be sufficient to account for the material properties differences between outer and inner annulus tissue.

Concerning the biological processes, the intervertebral disc is an avascular organ; it means that the cell nutrients are transferred from the upper and lower vertebrae through the disc by advective or diffusive mechanisms [START_REF] Grunhagen | Nutrient supply and intervertebral disc metabolism[END_REF][START_REF] Yao | Three-dimensional inhomogeneous triphasic finite-element analysis of physical signals and solute transport in human intervertebral disc under axial compression[END_REF] . Resulting from daily loading cycles on the disc, extracellular fluid displacements increase mass exchange between outer and inner disc regions and participate to the nutrients transport. These liquid flows could be modified by tissue mechanical strains [START_REF] Ferguson | Fluid flow and convective transport of solutes within the intervertebral disc[END_REF][START_REF] Soukane | Computation of coupled diffusion of oxygen, glucose and lactic acid in an intervertebral disc[END_REF] or flow characteristics of the cartilaginous end-plates [START_REF] Accadbled | A measurement technique to evaluate the macroscopic permeability of the vertebral end-plate[END_REF] . In this framework, the contribution played by the mechanical behavior of the annulus tissue has not been clearly established. Since the annulus fibrosus is highly hydrated, the liquid phase plays a significant role in its mechanical behavior. Some viscous effects are observed resulting from fluid filtration through the matrix and leads to hysteretic behaviors. Despite the fact that the intervertebral disc is subjected to daily loading cycles, the annulus response to cycling tensile test has not been much investigated.

In saturated porous media, the coupling effects between strains and fluid transfers are usually taken into account using a poro-elastic approach [START_REF] Best | Compressive mechanical properties of the human annulus fibrosus and their relationship to biochemical composition[END_REF][START_REF] Argoubi | Poroelastic creep response analysis of a lumbar motion segment in compression[END_REF][START_REF] Magnier | Nutrient distribution and metabolism in the intervertebral disc in the unloaded state: A parametric study[END_REF] . In such formulations, the skeleton Poisson's coefficient that governs the total volume change rate is of the highest importance. When dealing with lengthy fibers composite materials, an anisotropic behavior is observed. By relying on an homogeneous material point of view, the measured macroscopic Poisson's coefficients can be greater than the classical upper limit given by 0.5. In such cases, the total volume variation resulting from uniaxial tensile stresses can be negative depending on the fibers orientation. If the annulus tissue behaves as such composite materials, some non-usual coupling effects between mechanical behavior and fluid transfers can be observed. Its analysis and description is a major issue in order to develop a comprehensive model of the intervertebral disc.

The goal of this study is to develop a simple heuristic model able to describe the homogenized mechanical behavior (non-linearity, hysteresis) and to retrieve micro-structural characteristics (fibers reorientation, coupling with the fluid phase). This model is based on simple rheological elements such as linear springs or dampers laid out in a geometrical configuration similar to the woven micro-structure of annulus. Therefore, macroscopic non-linearities are generated by geometrical reorientation of fibers during stretching. This description underlines that the nonlinear mechanical behavior of annulus tissue mainly results from its particular woven architecture rather than from micro-scale nonlinearities. The characteristics of this model are identified using some experimental uniaxial tensile tests realized on pig and lamb annulus samples. Eventually, the mechanical behavior is discussed in terms of cyclic effects.

Microstructure model

The microstructure of annulus fibrosus is composed of collagen fibers embedded in the extracellular matrix. This matrix is mainly composed of water and elastin fibers which perform a junction between collagen fibers [START_REF] Yu | Microfibrils, elastin fibres and collagen fibres in the human intervertebral disc and bovine tail disc[END_REF] . The collagen fibers have a relatively large tensile ultimate strength while being flexible, they confer the softness and resistance to this biological tissue. The secondary network of elastin fibers allows to preserve the lateral micro-displacement of collagen fibers. The collagen fibers are oriented with an angle θ 0 with respect to the horizontal plane. This angle, θ 0 , generally ranges from 25 • to 45 • 6 . The annulus fibrosus is structured with concentric lamellae around the nucleus pulposus alternately oriented with an angle of +θ 0 and -θ 0 .

When focusing on an experimental sample (Fig. 1a), the microstructure of the annulus tissue can be represented by a regular layout of collagen fibers embedded in an homogeneous matrix (Fig. 1b). Using the symmetry properties of this structure, a representative elementary volume (REV) is described by one collagen fiber interacting with a triangular zone of matrix (Fig. 1c).

The collagen fiber is oriented with an initial angle, θ 0 . The fiber is supposed to have a linear elastic behavior characterized by an elastic modulus, E f , which is a fair assumption when subjected to low strains [START_REF] Holzapfel | Single lamellar mechanics of the human lumbar annulus fibrosus[END_REF] . The matrix composed of an isotropic mixture of elastin fibers saturated with physiological liquid is modeled through a visco-elastic behavior characterized by an elastic modulus, E m , and a viscosity, η m . Actually, the complete rheological model presented in Fig. (1c) should be assimilated to a mechanism since the orientation, θ, changes during macroscopic stretching. Indeed, since the rheological behaviors introduced are linear (elastic and visco-elastic), macroscopic non-linearities will directly results from geometrical modifications imposed by the circumferential strain, ε θ . Thereby, the mechanical behavior of annulus fibrosus depends on 4 parameters: 3 mechanical parameters (E m , η m , E f ) and 1 geometrical parameter θ 0 . This rheological model should be able to simulate both non-linear mechanical behavior and hysteresis over a loading cycle.

Geometrical relations

The geometrical orientation of collagen fibers, θ, is related to the initial angle, θ 0 , through the mechanical strains (Eqs. 1 and 2). The first relation (Eq. 1) gives the reorientation as a function of the collagen fibers strain, ε f , and the sample circumferential strain, ε θ , while the second one (Eq. 2) shows the dependence on the tansversal strain, ε z .

cos(θ) = cos(θ 0 ) 1 + ε θ 1 + ε f (1) sin(θ) = sin(θ 0 ) 1 + ε z 1 + ε f (2)
Experimental tensile tests are carried out at a constant velocity rate, V , for the length of the samples h 0 so that the circumferential strain, ε θ , is given by:

ε θ = V h 0 t (3)
where h 0 is the initial length of an experimental sample (≈ 10 mm).

Behavior laws

Concerning the mechanical behavior constitutive laws, a linear elastic model is used for the collagen fibers (Eq. 4) and a visco-elastic model for the matrix (Eqs. 5 and 6). These relations depends on the 3 mechanical parameters (E m , η m , E f ) introduced in the previous section.

σ f = E f ε f (4)
σ θm = E m ε θ + η m εθ (5)
σ zm = E m ε z + η m εz (6)

Equilibrium equations

By relying on the symmetry conditions associated with the rheological model described in Fig. (1), the corresponding mechanical equilibrium equations are given below:

σ θ = σ θm + 2σ f cos(θ 0 ) 1 + ε θ 1 + ε f (7) 0 = σ zm -2σ f sin(θ 0 ) 1 + ε z 1 + ε f (8)

Numerical simulation and parameters identification

This set of first-order non-linear differential equations (Eqs. 4 to 8) is solved using an implicit Euler's scheme associated with a Newton-Raphson's algorithm. The numerical time step is identical than for experimental tests, dt = 0.01 s. The convergence of Newton-Raphson's procedure is fixed to 10 -7 Mpa with respect to the equilibrium conditions.

Materials and method

Specimen preparation

Two lumbar discs (L3-L4) were harvested from cadaver of lambs. These lambs were 6 months old and obtained from the veterinary school of Lyon, France. To avoid a blood-clot formation into the cartilage end-plates, a subcutaneous injection of heparin (Innohep 4500U I, tinzaparine, Leo Pharma) was administered the day before the sacrifice. Under general anesthesia and appropriate analgesia, the spinal segment (L3-L4) was removed. This spinal segment was stored into a 10% DMSO solution (Dimethyl sulfoxide, a cryoprotectant) in a -12 • C freezer. Associated to this lumbar disc, four lumbar discs (L1-L5) were harvested from the cadaver of a domestic pig (160 Kg) obtained in a local slaughter house. This lumbar spine was excised and the specimens were also stored in a -12 • C freezer.

The discs were separated from the vertebral bodies by blunt dissection. Next, the posterior, lateral and anterior quadrants were separated from each disc using a surgical knife. In brief, all specimens were excised in anatomical cut parallel to the horizontal plane along the circumferential direction e θ . From each quadrant (Fig. 2a), a subsample of biological tissue was extracted to determine the solid mass fraction. It is evaluated by differential weighing with a high-precision balance (Sartorius) after 24 hours drying at 105 • . Some significant differences in composition between the pig and lamb specimens have been observed, revealed by the statistical t-test (p = 0.0027) (Fig. 3b). Another subsample of each quadrant has been analyzed by means of transmission electron microscopy (Fig. 3a). Using the same magnification (x3000), the microscopic analysis has shown a higher matrix fraction with lamb specimens than with pig ones. This confirms the global measures obtained by weighing technique (Fig. 3b). This microscopic analysis has also underlined a significant difference (p = 2 × 10 -7 ) between the collagen fibers mean size of lamb and pig (Fig. 3c).

Specimen preparation for mechanical testing

From each quadrant (Fig. 2a), one plane-parallel specimen was extracted using a specific tool.

This tool was realized with two razor blades separated with a 2 mm-spacer. The length of each specimen h 0 was limited using a surgical knife at approximately 10 mm. Then, to get stuck into the testing device grips, both specimen ends were glued into aluminum rings using cyanoacrylate adhesive. Specimen images were taken with digital video camera disposed in an optical microscope (ZEISS). The image resolution was 4µm/pixel. The specimen geometry was measured using a custom imaging analysis software. A multi-threshold technique was associated with a mean square method to obtained the main size in the radial direction, L r , and in the vertical direction,

L z . In

Mechanical testing

Prior to testing, the specimens were stored in a 0.15 M ol/l N aCl solution at T = 37 • C. The testing device was composed of a Texture Analyzer (TA-XT2, Stable Micro Systems, UK) with a 50 N load cell. Each sample was placed into the testing device grips and stretched during 5 loading cycles from 0 to 1 mm lengthening in order to reach a maximum strain of ε max = 10%.

The extension rate was 0.01 mm/s and the tension load was recorded every 0.01 s.

Data and statistical analysis

All statistical analysis were realized with the GRETL software [START_REF] Yalta | Gretl 1.6.0 and its numerical accuracy[END_REF] . The significance p-value is obtained with a Student paired t-test. A value of p < 0.05 indicates a significant difference while a value of p > 0.05 means a negligible difference.

Results

Stress-strain curves from uniaxial tensile tests are presented first, then the ability of the rheological model to represent the experimental behavior is addressed and finally, the geometrical implications of this model is analyzed.

Cyclic uniaxial tensile tests

The whole set of annulus samples were subjected to the experimental protocol described in subsection Mechanical testing. Concerning two annulus samples, the evolution of the circumferential tensile stress as a function of strain is plotted in Fig. (5). Only the first two cycles are presented here. These results exhibit a strong non-linear behavior with a rigidity increase as the stress increases. Similar observations have been done with all investigated specimens. This particular behavior leads to the stiffening of the intervertebral disc tissue for large strains. Concerning the loading cycle, a significant hysteresis is observed. This features an energy dissipation measured by the area inside one cycle. This agrees with the essential dissipative function of intervertebral discs. In the second cycle, the dissipated energy is systematically lower. The successive loading cycles generally converge towards a single curve as already shown by Tower et al. [START_REF] Tower | Fiber alignment imaging during mechanical testing of soft tissues[END_REF] .

Mechanical behaviors

From these experimental results, the physical properties associated with the rheological model, i.e., the initial fibers orientation, θ 0 , and the mechanical characteristics, (E f , E m and η m ), are computed using a Levenberg-Marquardt identification algorithm. Concerning two samples, the rheological model is compared to the first experimental loading cycle in Fig. (6). A fairly good agreement is obtained since our rheological model can accurately describe the stress-strain curve along the whole cycle. The physical characteristics identified with all specimens are summarized in Tab. (1). The fitting quality is evaluated based on the correlation factor, R 2 , showing a very good accuracy in each cases. Regarding to the mechanical parameters, the confidence intervals reveal high variabilities of the identified properties. This can partly be attributed to some sample geometry discrepancies and is discussed in the last section. However, when focusing on the initial fibers orientation, high-confidence values have been obtained without significant differences between pig and lamb samples (p = 0.38).

Model results

To estimate the efficiency of this rheological model, the mechanical behavior of two annulus samples have been simulated using the physical characteristics identified in the previous section. The characteristics of the pig specimen are (E f =6.6 Mpa, E m =0.0191 Mpa, η m =0.61 Mpa.s, θ 0 =22.01 • ), while the lamb specimen ones are (E f =12.3 Mpa, E m =0.023 Mpa, η m =0.59 Mpa.s, θ 0 =22.89 • ).

In as a function of the circumferential strain, ε θ . This shows that the collagen fibers extension is relatively low with a maximal extension of 2%. In compression, the orientation angle is large and the fibers are weakly stressed.

Discussion

The microstructural informations gathered in this work agree with previous researches concerning the global water content [START_REF] Johnstone | The fluid content of the human intervertebral disc. comparison between fluid content and swelling pressure profiles of discs removed at surgery and those taken postmortem[END_REF][START_REF] Iatridis | Measurements of proteoglycan and water content distribution in human lumbar intervertebral discs[END_REF] and the collagen fibers size [START_REF] Skaggs | Regional variation in tensile properties and biochemical composition of the human lumbar annulus fibrosus[END_REF][START_REF] Iatridis | Mechanisms for mechanical damage in the intervertebral disc annulus fibrosus[END_REF] . The significant differences observed between pig and lamb samples mainly results from the ages difference (lamb: 6 months, pig: 12 months). This is confirmed by the stronger mechanical characteristics observed with pig samples than those obtained with lamb ones.

The experimental unidimensional tensile tests reveal classical features of the annulus mechanical behavior, i.e., a non-linear behavior with stiffness increasing, an hysteretic evolution under loading cycle and a mechanical adaptation beyond the second loading cycle. These hysteretic cycles are associated with the visco-elastic nature of annulus tissue [START_REF] Iatridis | Shear mechanical properties of human lumbar annulus fibrosus[END_REF] . It gives prominence to the coupling between the global mechanical behavior and the viscous flow through the matrix that characterizes of poro-elastic media. Indeed, the visco-elastic properties of annulus tissue should be related to its liquid permeability [START_REF] Iatridis | The viscoelastic behavior of the non-degenerate human lumbar nucleus pulposus in shear[END_REF][START_REF] Riches | The internal mechanics of the intervertebral disc under cyclic loading[END_REF][START_REF] Yao | Effects of Swelling Pressure and Hydraulic Permeability on Dynamic Compressive Behavior of Lumbar Annulus Fibrosus[END_REF] .

Between the different specimens tested, some discrepancies, revealed by the large confidence intervals given in Tab. (1), are observed. Some possible reasons are enumerated below.

The low number of specimens under investigation do not allow to identify the rheological parameters with confidence. This aspect is emphasized in the case of lamb where only 6 samples were tested, while the 12 pig samples lead to a better estimation.

Moreover, the mechanical behavior of annulus tissue highly depends on the localization inside the disc [START_REF] Best | Compressive mechanical properties of the human annulus fibrosus and their relationship to biochemical composition[END_REF][START_REF] Yin | A homogenization model of the annulus fibrosus[END_REF] . Indeed, the macroscopic mechanical characteristics are governed by the fibers orientation while it ranges from ≈ 25 • in outer annulus lamellae to ≈ 45 • close to the nucleus [START_REF] Cassidy | Hierarchical structure of the intervertebral disc[END_REF][START_REF] Marchand | Investigation of the laminate structure of lumbar disc annulus fibrosus[END_REF][START_REF] Wagner | Theoretical model and experimental results for the nonlinear elastic behavior of human annulus fibrosus[END_REF] .

Meanwhile, the collagen fibers density increases in external annulus tissue [START_REF] Marchand | Investigation of the laminate structure of lumbar disc annulus fibrosus[END_REF] . These microstructural and composition variations affect significantly the reproducibility of experimental tests since the sampling is not sufficiently precise with respect to the localization.

Eventually, annulus samples were realized with a 2×2 mm cross-section. In the radial direction, e r , a sample includes several lamellae (Fig. 4b). Since the viscous flow through the matrix influences the sample global stiffness, the number of lamellae acts as a characteristic size in the visco-elastic behavior by governing the liquid pressure relaxation time in the radial direction.

Therefore, unixial tensile tests do not bring up the intrinsic mechanical behavior of one lamella, but rather provide the macroscopic behavior of an heterogeneous anisotropic structure [START_REF] Holzapfel | Single lamellar mechanics of the human lumbar annulus fibrosus[END_REF] . Since the number of lamellae differs from one sample to another, this aspect could explain the observed discrepancies.

With a low number of parameters, 3 mechanical and 1 geometrical characteristics, this rheological model is able to represent several features of the mechanical behavior of annulus fibrosus tissue, i.e., non-linear stiffness, hysteretic visco-elastic behavior, reorientation of collagen fibers while loading and anisotropy in the ( e θ ; e z ) plan. It should be noticed that, with this model, the strong non-linear behavior is obtained only for a limited range of the initial fibers angle, θ 0 ≈ 23 • . This is confirmed by the small confidence interval computed when identifying the initial fibers angle from experimental measures (Tab. 1). This initial angle is close to the fibers angle directly measured on outer annulus tissue, θ 0 ≈ 25 • 6 . With a larger fibers orientation angle, as observed with inner annulus tissue, the mechanical behavior is more linear [START_REF] Bass | Biaxial testing of human annulus fibrosus and its implications for a constitutive formulation[END_REF][START_REF] Wagner | Theoretical model and experimental results for the nonlinear elastic behavior of human annulus fibrosus[END_REF] . However, our tensile tests have been performed without preconditioning as usually done in previous experimental works [START_REF] Fujita | Radial tensile properties of the lumbar annulus fibrosus are site and degeneration dependent[END_REF][START_REF] Bass | Biaxial testing of human annulus fibrosus and its implications for a constitutive formulation[END_REF][START_REF] Guerin | Degeneration affects the fiber reorientation of human annulus fibrosus under tensile load[END_REF] .

As reported by Tower et al. [START_REF] Tower | Fiber alignment imaging during mechanical testing of soft tissues[END_REF] , a transition occurs when the collagen fibers become uncrimped.

By imaging techniques, they have shown that a significant fibers reorientation takes place during the preconditioning and depends on the choice of the reference configuration. Thus, in our case, not taking into account this transition could affect the quality of identified properties and leads to under-estimate the fibers initial angle. Anyway, this heuristic rheological model brings out that the fibers reorientation can account solely for the strong non-linear behavior without relying on the assumption that the stiffening results from a transition between crimped and uncrimped fibers [START_REF] Holzapfel | A constitutive model for fibrous tissues considering collagen fiber crimp[END_REF] .

Therefore, we assess that the non-linear behavior mainly results from the collagen fibers reorientation. This point of view relies on the strong microstructural assumption that the fibers reorientation is governed by the macroscopic strain, i.e., there is no relative displacement between the fibers and the matrix. This leads to a fibers reorientation rate ranging between -1.7 • /% and -2.0 • /% (Fig. b). This rate is about 2 times larger than direct measures based on imaging techniques proposed by Guerin and Elliott [START_REF] Guerin | Degeneration affects the fiber reorientation of human annulus fibrosus under tensile load[END_REF] , which seems to indicate that this microstructural assumption is not completely efficient.

Concerning Young's modulus of the collagen fibers, the average estimated value, E f =8.27

Mpa, is in agreement with already published results for low fibers strain, ε f ≤ 2% [START_REF] Holzapfel | Single lamellar mechanics of the human lumbar annulus fibrosus[END_REF][START_REF] Smith | Elastic Fibers Enhance the Mechanical Integrity of the Human Lumbar Anulus Fibrosus in the Radial Direction[END_REF] . In these conditions, the elastic linear assumption is therefore a fair representation of fibers mechanical behavior [START_REF] Eppell | Nano measurements with microdevices: mechanical properties of hydrated collagen fibrils[END_REF] . Regarding to the matrix stiffness, E m , the values presented in this work are about 20 times less than previous experimental results [START_REF] Iatridis | Mechanisms for mechanical damage in the intervertebral disc annulus fibrosus[END_REF][START_REF] Yin | A homogenization model of the annulus fibrosus[END_REF][START_REF] Smith | Elastic Fibers Enhance the Mechanical Integrity of the Human Lumbar Anulus Fibrosus in the Radial Direction[END_REF] . However, these characteristics should not be compared to each other since they do not have the same physical meaning. So far, the matrix has been described using only one stiffness while we proposed a visco-elastic behavior based on two properties, a stiffness, E m , and a viscosity, η m . Thus, matrix stresses are partly supported by the viscous flow leading to a lower matrix stiffness. To check this proposition, tensile tests must be carried out with different extension velocity (Eq. 3) and this is the scope of future work.

Relying on such description of the poro-mechanical behavior will have strong implications on the coupling between strains and fluid transfers. In particular, this kind of woven oriented structure can exhibit, at the macroscopic scale, values of Poisson's coefficients that lie outside of the classical range valid for homogeneous materials, 0 < ν < 0.5. This unusual behavior mainly depends on the fibers orientation angle, θ. With angles around 25 • , as observed with annulus tissue, the macroscopic transverse shrinkage is larger than the longitudinal extension. In the case of uniaxial tension, this leads to a negative total volume variation that cannot be recovered with classical poro-elastic formulations usually used [START_REF] Best | Compressive mechanical properties of the human annulus fibrosus and their relationship to biochemical composition[END_REF][START_REF] Argoubi | Poroelastic creep response analysis of a lumbar motion segment in compression[END_REF][START_REF] Magnier | Nutrient distribution and metabolism in the intervertebral disc in the unloaded state: A parametric study[END_REF] . Instead of a global swelling, the matrix shrinkage generates an additional fluid pressure and an outward liquid flow.

Therefore, when focusing on the cell nutrients transfer mechanisms from the vertebrae toward the nucleus under cyclic loading, these particular coupling effects significantly influence liquid flows. It is therefore a major issue to develop an accurate poro-mechanical formulation in order to propose a comprehensive model of the cell nutrition in the intervertebral disc.

Figure 1 :

 1 Figure 1: a) Representation of one sample and loading direction. b) Description of internal microstructure of annulus fibrosus. c) Rheological model of annulus fibrosus microstructure.

Figure 2 :

 2 Figure 2: a) Lumbar disc with the specimens locations. b) Specific tool used for cutting rectangular specimen. c) Annulus specimen with aluminum ring.

Figure 3 :

 3 Figure 3: a) Transmission electronic microscopy images of a lamb and a pig specimen. White regions correspond to the matrix and black ones to the collagen fibers. b) Solid mass fractions of lamb and pig specimens. c) Collagen fibers size for lamb and pig specimens measured with ImageJ.

Figure 4 :

 4 Figure 4: Numerical threshold technique. a) Specimen size evaluated in the plane ( e θ , e z ). b) Specimen size evaluated in the plane ( e θ , e r ).

Fig. ( 4 )

 4 , one can clearly observes the fibers orientation on the left image and the successive lamellae on the right one. Finally, 12 porcine and 6 lamb specimens were prepared to analyze the unidimensional mechanical behavior of annulus along the circumferential direction e θ .

Figure 5 :

 5 Figure 5: Experimental curves of the stress versus strain: a) pig specimen, b) lamb specimen.

Figure 6 :

 6 Figure 6: First stress-strain cycle curve of a) a pig specimen, b) a lamb specimen. Comparison between experimental measures and computed behavior.

  Fig. (), three plots are presented. The first one represents the mechanical behavior of both samples during the first loading cycle (Fig. a), emphasizing the strong non-linear behavior. The tension / compression asymmetrical behavior must be noticed since the computed stiffness is lower in compression than in traction as already shown by Wagner and Lotz[START_REF] Wagner | Theoretical model and experimental results for the nonlinear elastic behavior of human annulus fibrosus[END_REF] . This particular behavior mainly results from the reorientation of collagen fibers along the traction direction. This is represented by the evolution of the fibers angle, θ, as a function of strain, ε θ , in Fig.(b). The reorientation rate, defined by the average slope of this plot, is roughly constant and ranges between -1.7 • /% and -2.0 • /%. Indeed, as the fibers get oriented along the loading direction, e θ , they participate more in the sample rigidity that leads to a global stiffness increase. The major role of fibers reorientation is confirmed when focusing on the fibers strain, ε f , plotted in Fig.(c)

Figure 7 :

 7 Figure 7: Computed mechanical behaviors of two annulus samples subjected to one loading traction-compression cycle: a) stress-strain curve, b) fibers orientation, c) fibers strain.

Table 1 :

 1 Mechanical properties of pig and lamb specimens identified from experimental data (average[confidence interval at 95%]).

		E f (Mpa)	Em(Mpa)	ηm(Mpa.s)	θ 0 ( • )	R 2
	Lamb specimen
	(n=6)	12.3 [-7.1;31.8]	0.023 [-0.01;0.056]	0.59 [-1.49;2.68]	22.89 [21.58;24.20]	0.995 [0.99;0.998]
	Pig specimen
	(n=12)	6.6 [3.3;9.8]	0.0191 [0.009;0.028]	0.61 [0.25;0.97]	22.01 [20.69;23.33]	0.995 [0.993;0.997]
	Total (n=18)	8.3 [3.4;13]	0.0196 [0.01;0.028]	0.6 [0.3;0.91]	22.27 [21.32;23.22]	0.995 [0.993;0.996]

Acknowledgments

This research was supported by the C.N.R.S (Projets Exploratoires Pluridisciplinaires) and the University of Montpellier II.


Fabien Cherblanc

Jérôme Boscus

Jean-Claude Benet

Jean-Claude Bénet

Electro-osmosis in gel -Application to Agar-Agar

come

Introduction

Many geological and biological media (clays, cartilages, ...) exhibit swelling-shrinkage behaviour when submitted to hydraulic, electric or osmotic stresses. Since similar theoretical approaches are involved, the coupling between the mechanical behavior and transport phenomena is of great interest in geomechanical and biomedical engineering [START_REF] Huyghe | Preface on physicochemical and electromechanical interactions in porous media[END_REF]. The strains observed stem from liquid phase transport phenomena that take place in porous media (osmosis, electro-osmosis). These coupling effects observed at the macroscale generally results from electrical interactions at the microscale [START_REF] Mitchell | Fundamentals of Soil Behavior[END_REF]. The link between these two scales of description and the physical properties associated can de derived theoretically using an upscaling procedure [START_REF] Moyne | Electro-chemo-mechanical couplings in swelling clays derived from a micro/macro-homogenization procedure[END_REF], but also relies on a macroscopic experimental characterization since the microscopic quantities are hardly measurable.

When a porous medium saturated by an ionic solution exhibits some fixed electrical charges at the surface of the solid phase, the spatial distribution of ions and cations in the liquid phase is ruled by the electrical double layer theory [START_REF] Mitchell | Fundamentals of Soil Behavior[END_REF]. If the fixed charges are negatives, excess mobile cations are required in the solution to ensure electroneutrality. When submitted to an external electric field, these excess mobile cations exert an additional momentum to the pore fluid that drives the liquid phase toward the cathode. This net movement of fluid relative to the solid structure under an imposed electric potential gradient is called electro-osmosis [START_REF] Casagrande | Electro-osmosis in soils[END_REF]. In practice, the electro-osmotic velocity is described at a macroscopic scale by

v L = -k E ∇Φ (1)
where v L is an averaged velocity also called Darcy's velocity or filtration velocity, and Φ is the electrical potential. The electro-osmotic permeability, k E , characterizes the porous medium by relating the flux v L to the force ∇Φ. This relation can also be derived theoretically from the linear thermodynamics of irreversible processes (TIP) ideas [START_REF] Kuiken | Thermodynamics of Irreversible Processes[END_REF]. Electro-osmotic permeability of soils has been widely investigated in the civil engineering area since some potential applications in hazardous waste management are based on electrokinetic phenomena [START_REF] Alshawabkeh | Electrokinetic Remediation. II: Theoretical Model[END_REF][START_REF] Virkutyte | Electrokinetic soil remediation -Critical overview[END_REF]. For a large panel of soil, from clay to sand, the practical range of electro-osmotic permeability lies between 1 × 10 -9 to 1 × 10 -8 m 2 .V -1 .s -1 , while the hydraulic permeability can range over 6 decades [START_REF] Mitchell | Fundamentals of Soil Behavior[END_REF][START_REF] Yeung | Datla Fundamentals Formulation of Electrokinetic Extraction of Contaminants from Soil[END_REF]. This makes it possible some interesting application in clays, since a pressure-induced liquid flow can hardly be obtained in such low-permeable medium. While geological porous media have received a great deal of attention, little results are available about gels. Thereby, this contribution presents an experimental investigation dealing with the electro-osmotic characterization of Agar-Agar gel.

Material

Gels are usually assumed to be a cross-linked polymer-based network saturated with a liquid phase. The solid matrix ensures the gels holding excellent mechanical properties while the large water content make it possible the transport of diffusive ions. Gels are able to convert the chemical, electrical and thermal energy to the mechanical energy easily, what make them a very competitive candidate in the development of artificial muscles, sensors/actuators [START_REF] Chen | Modelling deformation behaviour of polyelectrolyte gels under chemo-electromechanical coupling effects[END_REF]. Their physical properties similar to articular cartilages account for their large use in tissue engineering [START_REF] Gu | Transport of multi-electrolytes in charged hydrated biological soft tissues[END_REF].

Agar-Agar is typically a strong gelling component widely used in the food industry as a thickening agent (E406), a vegetarian gelatin substitute, or a clarifying agent in brewing. Agar is an unbranched polysaccharides extracted from the cell walls of some species of red algae or seaweed. It can be considered as an heterogeneous mixture of two classes of polysaccharides: agaropectin and agarose. Although both polysaccharides classes share the same galactose-based backbone, agaropectin is heavily modified with acidic side-groups, such as sulfate and pyruvate [START_REF] Labropoulos | Dynamic rheology of agar gels: theory and experiments. Part II. gelation behavior of agar sols and fitting of a theoretical rheological model[END_REF]. The neutral charge and lower degree of chemical complexity of agarose make it less likely to interact with biomolecules, such as proteins. This confers to agarose the ability of being used as a culture medium for microbiological work and electrophoretic identification of proteins in DNA sampling.

Gel samples are prepared by dissolving dehydrated agar powder in hot water with a 7% mass ratio, and heating the solution at 90 • C to have an homogeneous and random distribution of agar molecules in the solution. Gelation occurs by cooling the agar solution to a temperature below the gelation point estimated around 45 • C [START_REF] Labropoulos | Dynamic rheology of agar gels: theory and experiments. Part II. gelation behavior of agar sols and fitting of a theoretical rheological model[END_REF]. In the widely accepted double helix model, agar chains are cross-linked to each other by forming a double helix of 1.9nm pitch. The resulting 3D polymeric network has a microscopic characteristic length about 0.2 µm [START_REF] Mrani | Two dimensional simulation of dehydration of highly deformable gel : moisture content, stress and strain fields[END_REF][START_REF] Pernodet | Pore size of agarose gels by atomic force microscopy[END_REF]. The final water content, measured by differential weighing after 48-hours drying at 105 • C, is 1300%. The solid structure exhibits some electrical negative charges resulting from the ionization of acidic groups carried by agaropectin chains; Agar-Agar is generally described as a polyelectrolyte gel. Consequently, Agar-Agar gel should be sensitive to electrokinetic phenomena, what is the scope of the next section dealing with the experimental characterization of the electro-osmotic permeability.

Experimental Device

The complete experimental device is given in Fig. 1. The Agar-Agar gel sample (a) (64 mm diameter and 26 mm thick) is held between two parallels perforated Plexiglas plates (b) (64 mm diameter and 5 mm thick). These plates applied a slight axial compression stress on the gel sample to ensure the no-flow condition on the periphery. The whole system is saturated by pure distilled water.

The capillary tube on the left (c) is used to determine the water flow through the sample. A 4 mm diameter tube allows us to measure very low fluxes (Fig. 1) while the water evaporation is prevented by disposing an oil bubble on the left edge (d). The suction resulting from capillary effects at the water-oil interface is fairly negligible. For one experimental measure, the electroosmotic flux in the sample causes an horizontal displacement, d, of water in the capillary tube less than 15 cm. This corresponds to a water rise in the right tank about 0.1 mm, which is not sufficient to create some pressure gradient discrepancies.

To impose an electrical potential gradient through the gel sample, two electrodes (anode on the left (e) and cathode on the right (f)) are immersed in each water tank. They are connected to a controlled voltage generator (g) via an amperometer (h). These electrodes are made of carbon cloth to take advantage of the inert characteristic of carbon with respect to water. A second electrical circuit is used to measure the effective electrical gradient. It is composed of two platinum electrodes (i), each of them is stuck between the porous plate and the sample. These electrodes are connected to a high impedance voltmeter (j).

The time required for one experimental measure is less than 1 hour, while the electrical intensity is always lower than 1 mA. Thus, the gas volume of oxygen produced by electrolysis at the anode is negligible and do not leads to some discrepancies while estimating the liquid flow rate through the sample.

The experimental procedure is relatively straightforward. Imposing an electrical potential difference between both electrodes (e)-(f) generates an electrical gradient inside the gel sample. The electro-osmotic flow produced is evaluated through the measure of the oil bubble displacement. Some experimental results are presented in the next section.

Results and Discussion

With a 10 V potential difference applied between electrodes, the experimental measures are given below. The oil bubble displacement, d, is given in Fig. 2a. Using a centered finite difference scheme, the corresponding water flow rate is computed in Fig. 2b. The derivation operator amplifies the measurement errors, what explains the observed fluctuations. Anyway, a smoothed trend curve could easily be plotted inside the confidence intervals.

While the electrical potential difference applied between the carbon electrodes (e)-(f) is regulated at 10 V , the electrical potential difference effectively applied at the boundaries of the sample is measured using the platinum electrodes (i) (Fig. 3a). A slight increase was always observed. Indeed, due to the water electrolysis, water in the right tank gets acid, while it gets basic in the left tank; this increases the water electrical conductivity in both tanks. Thus, since the electrical dropdown gets lower in the tanks, the electrical difference applied between both sides of the sample is increasing. The electrical intensity, resulting from ionic transport in the liquid phase, is given in Fig. 3b. A slight decrease is generally observed in early time. This corresponds to the establishment of the diffusive layer around the carbon electrodes.

time [s] time [s] displacement [mm] electro-osmotic flow [m /s]
Using Eq. 1, the electro-osmotic permeability k E is calculated at each time step (Fig. 4). Except for the early times (t < 1000 s), the electro-osmotic permeability measured is fairly constant and can be approximated by k E = 1, 52×10 -9 m 2 .s -1 .V -1 . This result validates the ability of this experimental procedure to determine the electro-osmotic permeability in gel. This procedure has been repeated by applying increasing electrical potential difference 10, 20 and 30 V on different gel sample. All the results are plotted in Fig. 5. On can notes that the water velocity is a superficial velocity, since it is the flow rate divided by the sample cross-area. The electrical gradient corresponds to the electrical potential difference measured between platinum electrodes divided by the sample thickness. By averaging all the experimental measures, the resulting electro-osmotic permeability of Agar-Agar gel is

k E = 1, 29 × 10 -9 m 2 .s -1 .V -1 (2)
The discrepancies observed between results are partly attributed to the differences of gel characteristics between samples. Indeed, the physical properties of Agar-Agar gel are very sensitive to the gelification process (temperature, heating time, cooling velocity, ...) [START_REF] Aymard | Influence of thermal history on the structural and mechanical properties of agarose gels[END_REF], and the elaboration repeatability can hardly be achieved.

As noted in the introduction, experimental values of the electro-osmotic conductivity is between 1 × 10 -9 to 1 × 10 -8 m 2 .V -1 .s -1 for a broad range of soils, from clay to sand [START_REF] Mitchell | Fundamentals of Soil Behavior[END_REF][START_REF] Yeung | Datla Fundamentals Formulation of Electrokinetic Extraction of Contaminants from Soil[END_REF]. The electroosmotic conductivity of Agar-Agar gel, k E (Eq. 2), also lies in this range, even if this material does not share any micro-structural characteristics with soils. It should be recalled that the liquid phase saturating the gel is distilled water, i.e., the only ionic solute results from the water dissociation ([OH -] = [H + ] = 10 -7 mol.l -1 ). It seems that this range of values would have an universal meaning which do not depend on the kind of porous media, the micro-structural geometry or the ionic concentration of the solution.

Based on the Helmholtz-Smoluchowski theory [START_REF] Vane | Effect of aqueous phase properties on clay particle zeta potential and electro-osmotic permeability: Implications for electro-kinetic soil remediation processes[END_REF], the zeta-potential, ζ, can be deduced from the electro-osmotic permeability according to where µ = 0.89 × 10 -3 kg.m -1 .s -1 is the dynamic viscosity of water, τ = 1, 3 is the tortuosity estimated from Lakatos et al., 2004 [16], φ = 0.96 is the porosity, ε = 6, 90 × 10 -10 F.m -1 is the permittivity of water. In the case of Agar-Agar gel, this leads to ζ = -2, 93 × 10 -3 V , which is about one order lower than the values measured with soils.

ζ = - µτ 2 φε k E (3)

Conclusion

In food-and bio-engineering research, Agar-Agar gel is commonly used as reference material to model the swelling/shrinking behaviour induced by hydraulic, electric or osmotic stresses. An experimental investigation dealing with the electro-osmotic characterization of Agar-Agar gel has been presented. The electro-osmotic phenomenon is described by a linear relation, widely used in geotechnical applications, which defines the electro-osmotic permeability. The experimental results show that the electro-osmotic permeability of Agar-Agar gel lies in the limited range usually observed on a large panel of soils, even if the micro-structural characteristics are really different.

Figure 1 :

 1 Figure 1: Experimental device for the determination of the electro-osmotic permeability of Agar-Agar gel sample

Figure 2 :Figure 3 :Figure 4 :

 234 Figure 2: Experimental results with a 26mm thick sample submitted to a 10V electrical potential difference: (a) oil bubble displacement as a function of time -(b) electro-osmotic flowrate through the sample as a function of time.

1 Figure 5 :

 15 Figure 5: Synthesis of all experimental results: water velocity through the gel samples versus the electrical potential gradient inside the samples. The equation of the trend line is given with its correlation factor.


Jocelyn Étienne
email: etienne@ujf-grenoble.fr

J Ér

Ôme Loh Éac

Pierre Saramito
email: pierre.saramito@imag.fr

A LAGRANGE-MULTIPLIER APPROACH FOR THE NUMERICAL SIMULATION OF AN INEXTENSIBLE MEMBRANE OR THREAD IMMERSED IN A FLUID *

Keywords: Fluid-structure interactions, surface divergence, saddle-point problem, boundary Lagrange multiplier AMS subject classifications. 65M60, 74K05, 74K15, 74F10, 76D05, 76D07

The inextensibility constraint is encountered in many physical problems involving thin solids interacting with a fluid. It is generally imposed in numerical simulations by means of a penalty method. Here, we propose a novel saddle-point approach allowing to impose it through a Lagrange multiplier defined on the thin structure, the tension. The functional analysis of the problem allows to determine which boundary conditions are needed for this problem. The forces originating from the structure appear as a boundary condition for the fluid problem, defined on a moving boundary which represents the structure. The problem is discretised with mixed finite elements. The mesh of the thin solid is included in the mesh of the bulk, and is advected by its velocity in the course of time iterations. The appropriate choice of the finite element spaces for this mixed approach is discussed and it is shown that boundary conditions on the thin structure edges impact on this choice. Numerical tests are performed which demonstrate the convergence and robustness of the method. These examples include the simulation of a closed, membrane bound object (a vesicle) and of a filament or flag in a large Reynolds-number flow.

1. Introduction. From sailing to micro-organism swimming, there are innumerable examples of dynamics involving the interactions between a fluid and a thin elastic sheet or a slender filament. Because they are thin or slender, the geometry of these elastic structures introduces a small parameter: the thickness, or the radius, compared to the length. When this parameter is taken to its limit, this problem is a generalisation of interfacial dynamics problems, where the governing equations in the fluid domain are supplemented with a jump boundary condition for the stress at the interface.

An elastic thread is fully characterised by its bending and extension-compression elasticity, while an elastic membrane additionally has a shear elasticity [START_REF] Helfrich | Elastic properties of lipid bilayers: Theory and possible experiments[END_REF]. Here we are interested in the case when the extension-compression modulus is extremely high, that is, of nearly inextensible (and incompressible) threads and membranes. This case is relevant to many physical materials, including most fabric structures such as threads, sails, flags; metal shells such as rods or wings; and biological membranes and flagellae. Some form slender structures (one-dimensional), such as threads, others form thin sheets (two-dimensional), such as flags, leather balls or biological membranes. Another classification is whether the structures are 'closed' such as balls and living cells, or 'open' such as threads, flags, and flagellae.

A very high extension modulus makes the problem very ill-conditionned, and not practical for numerical resolution. Two options are available : either to reduce artificially the value of the modulus, but keeping it high, this is the penalty approach, or to consider that this modulus is infinite, making the thread totally inextensible, which leads to a constrained problem. The penalty approach introduces a parameter, the penalty extension modulus, which needs to be adjusted in simulations to achieve convergence. The numerical result will depend on the penalty modulus. Also, since external forces may vary locally, the actual, non-physical extension of the thread will be higher in regions of strong sollicitation. On the other hand, the constrained problem approach introduces no additional free parameter and imposes inextensibility uniformly. It is thus physically more attractive for approaching very high extension modulus cases, but requires more involved numerical methods.

In this paper we develop a numerical method for the resolution of this constrained problem. The augmented Lagrangian approach of Fortin and Glowinski [START_REF] Fortin | Augmented Lagrangian methods, applications to the numerical solution of boundary value problems[END_REF] is widely used for solving problems subject to a constraint, such as volume incompressibility in Stokes and Navier-Stokes equations. Its application to the case of an inextensible thread or membrane immersed in a fluid is described in section 2. A saddle-point approach allows us to characterise the solution in a weak formulation, which is discretised using mixed finite elements in section 3. The choice of finite-element discrete spaces is discussed. Numerical experiments validate the method in section 4, where simulations of both 'closed' membrane-bound objects and 'open' threads with a 'free' endpoint are presented, both in Stokes and high-Reynolds number Navier-Stokes flows.

2. Governing equations and saddle-point approach. Consider a two dimensional domain Ω and a curve Γ 0 = {γ 0 (s), s ∈ (0, 1)} ⊂ Ω. This is the initial geometry of the problem, the curve Γ(t) = {γ(t, s), s ∈ (0, 1)} ⊂ Ω being one of its unknowns, with γ(0, s) = γ 0 (s), and represents the thread immersed in a twodimensional fluid flow in Ω (alternatively, using plane symmetry it can represent a membrane of infinite width in a three-dimensional fluid flow in Ω × (-∞, +∞), or, using axial symmetry, a membrane ring in Ω × (0, 2π)). We denote by A(t) = γ(t, 0) and B(t) = γ(t, 1) the endpoints of Γ(t), and s is the curvilinear coordinate along Γ(t). We introduce t = ∂γ/∂s the tangent to Γ(t), and n the normal such that t × n = 1.

Let us assume that this thread or membrane is inextensible and has no resistance to bending, and that the fluid is incompressible and inertia-less. (The inertial case will be considered in the end of this section.) The fluid is set into motion by Dirichlet boundary conditions, u = u D , imposed on Γ D = ∂Ω. The extension to Neumann boundary conditions is straightforward. The question of the necessity of prescribing additional boundary conditions to determine the behaviour of the thread or membrane will be addressed in this section as we expose the governing equations of the problem.

Over a time interval (0, T ), the problem writes:

Find γ ∈ C 0 ([0, T ], H 3/2 (0, 1)), u ∈ L ∞ 0, T, V div (u D , t) , such that: J(u(t)) = min v∈V div (u D ,t) J(v) a.e. t ∈ [0, T ], γ(t, s) = γ 0 (s) + t 0 u(τ, γ(τ, s)) dτ, ∀(t, s) ∈ [0, T ] × [0, 1], with Γ(t) = {γ(t, s), s ∈ (0, 1)} , V div (u D , t) = v ∈ H 1 (Ω) 2 , v Γ(t) ∈ H 3/2 (Γ(t)) 2 , v = u D on Γ D , div v = 0 in L 2 (Ω), div s v = 0 in Z ′ (t) , and
J(v) = Ω |D(v)| 2 dV,
where

D(v) = 1 2 grad v + (grad v) T
is the rate-of-strain tensor, and the integration variable is dV = dx dy in plane symmetry, and dV = r dr dz in axial symmetry. (In that case, the functional spaces are understood as weighted Sobolev spaces, see [START_REF] Kufner | Weighted Sobolev Spaces[END_REF]). The surface divergence operator is defined as div s v = div v -((grad v)n) • n. The functional space Z ′ (t) is defined on Γ(t) and will be characterised in the sequel. Note that in this particular case, no non-dimensional group arises when deriving these equations.

When Γ(t) is known, this problem is equivalent to finding a saddle point (u;

p, ζ) ∈ V (u D ) × Q × Z(t) of the Lagrangian function [2]: L(v; q, ξ) = Ω |D(v)| 2 dV - Ω q div v dV - Γ(t) ξ div s v dA, (2.1)
with dA = ds or dA = rds respectively in plane and axial symmetries. We have introduced the spaces :

V (v D ) = v ∈ H 1 (Ω) 2 , v = v D on Γ D , Q = q ∈ L 2 0 (Ω) ,
and Z(t) ⊂ H 1/2 (Γ(t)), the dual of Z ′ (t), may have to include boundary conditions for ζ. These boundary conditions will be detailed in the sequel. The Lagrange multipliers p and ζ are respectively the pressure in the fluid and the tension in the membrane.

A saddle-point is characterised by the weak formulation:

Ω 2D(u) : D(v) dV - Ω p div v dV - Γ(t) ζ div s v dA = 0 ∀v ∈ V (0) - Ω q div u dV = 0 ∀q ∈ Q - Γ(t) ξ div s u dA = 0 ∀ξ ∈ Z(t)
A generalisation of Green's formula for curved boundaries (e.g. [START_REF] Laadhari | On the equilibrium equation for a generalized biological membrane energy by using a shape optimization approach[END_REF]) can be formally applied to the boundary integrals above, expliciting the operator div s . Note that this introduces a source term when Γ(t) is curved (denoting κ its curvature):

Γ(t) ξ div s v dA = - Γ(t) ∂ξ ∂s t + κξn • v dA + [ξv • t] B(t) A(t) . (2.2)
The above equation does not generally make sense for ζ ∈ H1/2 (Γ) and v ∈ H 1/2 (Γ) [6, 7, lecture 33]. However, the physics of the problem dictates boundary conditions at points {A, B} = ∂Γ(t) for these functions, which are, either, a fixed velocity or a fixed tension. For the sake of simplicity, we will only consider homogeneous cases of zero velocity or zero tension.

Zero-tension boundary conditions thus imply the choice

Z(t) = H 1/2 00 (Γ(t)). Then Eq. 2.2 makes sense if v ∈ L ∞ (Γ(t)) 2 , and κn = ∂ 2 γ/∂s 2 is in the dual of H 1/2 00 (Γ(t)).
A sufficient condition for these two requirements is that u and v be in

H 3/2 (Γ(t)) 2 .
The case of zero-velocity conditions can also be treated, taking

Z(t) = H 1/2 (Γ(t)) and u ∈ H 3/2 0 (Γ(t)) 2 ⊂ H 1/2 00 (Γ(t)) 2 . Indeed, the product u ξ is in H 1/2
00 (Γ(t)) (this can be verified using Hardy's inequality, see [START_REF] Tartar | An introduction to Sobolev spaces and interpolation spaces[END_REF], and the fact that H 1/2 (Γ) ֒→ L 8 (Γ)). Finally, one can apply a fixed velocity at one end of the moving boundary and a fixed tension at the other using the same ideas.

Note that the tension, unlike the pressure in Dirichlet boundary problems, is not defined up to a constant. The absolute value of the tension is either fixed by a boundary condition at a 'free' endpoint (one not in Γ D ) or by the term κξn if κ is not everywhere zero 1 . The value of the tension at a 'free' endpoint should be zero [8, p. 130], unless there is an external force other than hydrodynamic forces applied on that end.

We define:

a(u, v) = Ω 2D(u) : D(v) dV b(v, q) = - Ω p div v dV c Γ(t) (v, ξ) = Γ(t) ∂ξ ∂s t + κξn • u dA
Using these notations, we rewrite the characterisation of the saddle-point

(u; p, ζ) ∈ V (u D ) × Q × Z(t) defined in Eq. 2.1: a(u, v) + b(v, p) + c(v, ζ) = 0 ∀v ∈ V (0), (2.3a) b(u, q) = 0 ∀q ∈ Q, (2.3b) c Γ(t) (u, ξ) = 0 ∀ξ ∈ Z(t), (2.3c)
For a sufficiently regular solution, the saddle point problem is thus equivalent to the following strong formulation in terms of its momentum equation and the constraints:

grad p -div 2D(u) = 0 in Ω \ Γ(t) div u = 0 in Ω div s u = 0 on Γ(t) Γ(t) = γ(t, s) = γ 0 (s) + t 0 u(τ, γ(τ )) dτ, ∀s ∈ (0, 1)
subject to the boundary condition:

-[p]n + [2D(u)]n = ∂ζ ∂s t + κζn on Γ(t) u = u D on Γ D ζ(P ) = 0 for P ∈ ∂Γ(t) \ Γ D
where the notation [•] refers to the jump of a quantity accross the interface Γ(t).

In the case when there is inertia, the energy variation term J needs to include kinetic energy of both the fluid and the thread, and becomes :

J(v) = Re d dt Ω v 2 dV + Re Γ d dt Γ(t) v 2 dA + Ω |D(v)| 2 dV,
where Re = ρU L/η and Re Γ = λU/η are two Reynolds numbers associated with the fluid of surfacic density ρ and viscosity η, and the thread of lineic density λ and length L, for a typical velocity U . An initial condition for u must then be provided. The same procedure can be applied to impose the inextensibility of the thread. Finally, let us note that any other volume or surface force can be easily taken into account in this setting, such as a bending rigidity [START_REF] Helfrich | Elastic properties of lipid bilayers: Theory and possible experiments[END_REF].

3. Finite element discretisation. In order to have explicit meshes of both Ω and Γ(t), and no interpolation between the spaces defined on these meshes, we first introduce an approximation of Γ 0 :

Γ 0 h = e 0 i = [x 0 i , x 0 i+1 ], x 0 i ∈ Γ(t), |e i | < h, ∀i = 0, ..., N -1, x 0 0 = A 0 , x 0 N = B 0
The finite element mesh T 0 h of Ω, composed of triangles of diameter smaller than h, is then defined such that all the edges e 0 i are among the sides of the triangles. Being given an approximation u n h of the velocity field u at time t n , a sequence of approximations of the moving interface Γ(t) at time t n+1 = t n + ∆t can be constructed as:

Γ n+1 h = e n+1 i = e n i + ∆t [u n h (x n i ), u n h (x n i+1 )], ∀i = 0, ..., N -1 . (3.1)
and hence a triangulation T n+1 h can be obtained either by remeshing or by advecting interior node with an ad-hoc mesh velocity c n h , this latter option being known as the arbitrary Lagrangian-Eulerian technique [START_REF] Hirt | An arbitrary Lagrangian-Eulerian computing method for all flow speeds[END_REF]. Inextensibility allows to derive an upper bound in O(h(1 + ∆t) n ) for the length of the edges e n+1 i , periodic remeshing of Ω n+1 h and occasionnally of Γ n+1 h are nevertheless necessary to ensure mesh quality. We define

c n h ∈ X n h (u n h ) = V (0) ∩ {v h ∈ P 1 (T n h ), v h (x n i ) = u n h (x n i ), ∀x n i ∈ Γ n h , i = 0, ..., N -1} as the solution of the problem: Ω D(c n h ) : D(v h ) dV = 0 ∀v h ∈ X n h (0). (3.2) Fig. 3.1.
The choice of finite element spaces for a thread with no 'free' endpoint. Velocity u is continuous, quadratic by piece; pressure p linear by piece and continuous except accross the thread; and tension ζ is continuous, linear by piece (defined only on the thread). Velocity discretisation is shown for two elements only.

We then introduce the finite element spaces P k (T h ) of piecewise polynomial functions of degree k. We define

V n h (v D,h ) = V (v D,h ) ∩ P 2 (T n h ) ∩ C 0 (Ω), and Q n h = Q ∩ P 1 (T n h ) ∩ C 0 (Ω \ Γ n h ).
The pressure approximation is thus taken continuous except accross the (approximation of the) interface. Because this discontinuity is thus localised, this does not break the Babuška-Brezzi condition [START_REF] Brezzi | Mixed and hybrid finite elements methods[END_REF] as far as the velocitypressure pair is concerned.

Further, we introduce

Z n m,h = ξ h ∈ P m (Γ n h ), ξ h (P ) = 0 for P ∈ ∂Γ n h \ Γ D ∩ C 0 (Γ n h ),
whose degree m has to be determined such that the Babuška-Brezzi condition holds for the composite Lagrange multiplier formed of pressure and tension. It turns out empirically that, when the contiuous space is Z(t) = H 1/2 (Γ(t)), using m = 1 (which correponds to the discretisation pictured in Fig. 3.1) leads to a stable scheme. This scheme is shown to converge to the analytical solution on an example in section 4.1.

Thus this seems to be the appropriate choice when the interface has no 'free' endpoint (i.e., when ∂Γ(t) \ Γ D = ∅). However, in the case of an interface with one 'free' endpoint, where a boundary condition must be imposed on the tension (e.g., when Z(t) = H 1/2 00 (Γ(t))), oscillations appear in the velocity field on the interface, close to this 'free' endpoint. From a discrete viewpoint, these correspond to modes in the kernel of the inextensibilty operator, due to the deletion of the degree of freedom corresponding to the Dirichlet boundary condition on the tension. For Γ 0 with κ = 0 (straight segment), it is in fact quite easy to construct elements ũh of V 0 h such that div s ũh = 0 but c Γ n h (ũ h , ξ h ) = 0 for all ξ h ∈ Z 0 1,h . The interpretation is thus that the P 1 (Γ h ) test-space is too small to enforce inextensibility close to the endpoint of the moving boundary. Adding an extra degree of freedom for the tension locally in the last edge element does not help this behaviour, but only reproduces the same problem in the penultimate edge element. Using m = 2 in turn, that is, quadratic-by-piece tension, allows to enforce inextensibility, but results in a spurious oscillation mode at the mid-edge values of the tension. It is thus not obvious empirically, whether an adequate choice of spaces exists that would verify the Babuška-Brezzi condition. However, in the case m = 2, oscillations in tension can be filtered a posteriori by projection to Z n 1,h [START_REF] Brezzi | Mixed and hybrid finite elements methods[END_REF]. This strategy was employed in sections 4. [START_REF] Fortin | Augmented Lagrangian methods, applications to the numerical solution of boundary value problems[END_REF] Step b is the problem considered by Fortin ang Glowinski [START_REF] Fortin | Augmented Lagrangian methods, applications to the numerical solution of boundary value problems[END_REF], and is solved using an augmented Lagrangian technique with Uzawa algorithm.

) ∈ V n h (u D,h ) × Q n h × Z n m,h , such that, for all (v h , q h , ξ h ) ∈ V n h (0) × Q n h × Z n m,h : a(u n h , v h ) + b(v h , p n h ) + c Γ n h (v h , ζ n h ) = 0, (3.3a) b(u n h , q h ) = 0, (3.3b) c Γ n h (u n h , ξ h ) = 0. (3.3c) c) If m = 2, find ζ n h ∈ Z n 1,h such that Γ n h (ζ n h -ζn h )ξ h dA = 0 ∀ξ h ∈ Z n 1,h , else set ζ n h = ζn h . d)
When Re > 0, we make use of the method of characteristics to calculate the material derivative, as described for both the Lagrangian and ALE methods in [START_REF] Étienne | A Lagrangian-Eulerian approach for the numerical simulation of free-surface flow of a viscoelastic material[END_REF]. An extra term then appears in the right hand side of Eq. 3.3a (and also in the right hand side of Eq. 3.3c in the case Re Γ > 0).

Numerical tests.

Swollen vesicle in Stokes flow.

A spherical object, filled with an incompressible fluid and bounded with an inextensible membrane is rendered undeformable by the combination of these two constraints. This can easily be seen from the fact that any volume-preserving deformation of a sphere increases its area. As a result, such a vesicle must behave as a hard sphere, possibly with surface flow at the boundary since the membrane is fluid. This surface flow may be quenched by choosing a setup with appropriate symmetries, so that the vesicle will behave as a hard sphere with no-slip boundary conditions.

Thus it appears that a demanding test of our method is to calculate the permanent flow generated by the fall or rise of a buoyant vesicle with a fluid membrane in a fluid at rest, see Fig. 4.1. In the frame of the vesicle, this problem writes, Find u, p and ζ such that: subject to the boundary conditions: ρ being the buoyancy, set to one in Ω d and zero elsewhere. No boundary condition is necessary for the tension, since the endpoints of Γ are on the axis Γ axis and have zero velocity in the direction e r , which is tangent to the membrane at these points. Here the problem is written as static, because we will see that it results in u = 0 on the membrane Γ provided that u 0 is chosen such that gravity and drag forces balance.

grad p -div 2D(u) = -ρe z in (R + × R) \ Γ, div u = 0 in R + × R, div s u = 0 on Γ,
-[p]n + [2D(u)]n = ∂ζ ∂s t + κζn on Γ, u → u 0 e z as r 2 + z 2 → ∞, u • e r = 0 ∂u • e z ∂r = 0 on Γ axis ,
Indeed, because of the axial symmetry of the problem, the membrane is actually at rest in the frame of the mass centre of the vesicle (no tangential flow), and thus the vesicle behaves as a solid sphere. Because in addition there is no shear source in the fluid inside the vesicle, it is at rest and the problem reduces to the flow generated by a buoyant inelastic solid sphere. This problem consists of the problem above stated in Ω e only and with boundary condition u = 0 on Γ and has an analytical solution [START_REF] Lamb | Hydrodynamics[END_REF],

u • e R = u R = u 0 sin Ψ 1 - 3 2R + 1 2R 3 , u • e Ψ = u Ψ = u 0 cos Ψ 1 - 3 4R - 1 4R 3 ,
written in spherical coordinates with R = √ r 2 + z 2 and Ψ = tan -1 z r , e R and e Ψ being the corresponding basis vectors. This also yields the drag force, and hence the value of u 0 balancing the gravity force, u 0 = 2/9.

Using this analytical solution, we can prescribe a boundary condition at a finite distance from the sphere (Ω = (0, 10)×(-10, 10)). We simulate the problem with four meshes of different refinement, Ω h with h ∈ {2h 0 , h 0 , h 0 /2, h 0 /4}, which are obtained by refining or coarsening the original mesh Ω h0 of approximately 2500 elements. The H 1 error u e -u h 1 is found to decrease as O(h 1.5 ) (Fig. 4.2), which corresponds to the expected rate of convergence. The sub-optimality of this rate compared to O(h 2 ) convergence of Stokes flow in a polygonal domain is due to the absence of an isoparametric implementation of the finite elements method. The method is seen to be nearly as accurate as the direct implementation of a solid sphere.

4.2.

Flow past a plate. We now turn to the case when the curve Γ(t) has one 'free' end, but still in conditions in which there will be no time evolution of Γ(t) because inextensibility forces balance hydrodynamic forces. This is the case of a 'flag', or a piece of thread in a two-dimensional flow such as a soap film [START_REF] Zhang | Flexible filaments in a flowing soap film as a model for one-dimensional flags in a two-dimensional wind[END_REF] if the thread is initially aligned with the flow direction. This setup simulates the flow past a plate, a

ζ ∼ a √ 1 -x, u x ∼ a 4 p max{x -1, 0} and ∂u x ∂y (x, 0 + ) = -∂ζ ∂x ∼ a 2 √ 1-x .
Parameter a was identified by best fit of ζ h over the range x ∈ [0.9, 0.999] and thence used for all the comparisons. The friction stress is multivalued at B, however its limit for x = (1ε, 0) is close to the asymptotic value. The tension should tend to 0 at A, this is a boundary condition for the P 2 field ζh , but is not verified pointwise by its projection ζ h on P 1 (see insert). However, the value is correct for degrees of freedom at x < 1.

problem equivalent to the Stokes 'stick-slip' problem, for which an asymptotic solution is known [START_REF] Richardson | A 'stick-slip' problem related to the motion of a free jet at low Reynolds numbers[END_REF][START_REF] Georgiou | A singular finite element for Stokes flow: the stick-slip problem[END_REF]. Fig. 4.3 describes the geometry and boundary conditions, with Ω = [-5, 5] × [-2 × 2] and ∂Ω = Γ out ∪ Γ D . The initial condition is γ 0 (s) = (s, 0) T , and our analysis states that at later times, γ(t, s) = γ 0 (s). We test the numerical implementation with respect to this result and the asymptotic solution of the Stokes 'stick-slip' problem.

We find empirically that the discrete inf-sup condition is violated for the choice (u h , p h , ζ h ) ∈ V h (u D,h )×Q h ×Z 1,h : the velocity on Γ h oscillates in a neighbourhood of point B. Using a Lagrange multiplier ζh ∈ Z 2,h these velocity oscillations disappear. Oscillations however appear on ζh , which can be filtered thanks to the projection method described in Section 3 to obtain the tension ζ h ∈ Z 1,h . In Fig. 4.4 we compare the numerical solution obtained with our scheme to the asymptotic solution around point B. It is seen that the numerical solution verifies closely this asymptotic solution.

In particular, the velocity on Γ h is zero up to floating point operations error (10 -11 ), and increases sharply in the neighbourhood of the flag (see Fig. 4.5). The numerical scheme is thus successful at describing a truly inextensible thread in conditions where the stress diverges at the endpoint.

Flag in a high

Reynolds-number flow. With the same initial conditions as in the previous section, we find that the flow can be unstable for a finite value of the Reynolds number. Such an instability is found experimentally for a filament in a twodimensional fluid flow [START_REF] Zhang | Flexible filaments in a flowing soap film as a model for one-dimensional flags in a two-dimensional wind[END_REF], although it is generally believed that it is due to the interplay of fluid inertia and (bending) elasticity of the filament [START_REF] Argentina | Fluid-flow-induced flutter of a flag[END_REF]. Here there is no In this example, we choose Re = Re Γ = 10 3 . The evolution of Γ h (t) is shown on Fig. 4.6. It is found that perturbations of the order of the numerical error are sufficient to trigger oscillations that grow to first order amplitude. These oscillations are the topic of a forthcoming paper.

Conclusions.

We develop and demonstrate a novel numerical method for the simulation of problems involving a constraint on a moving boundary, in the present case, the inextensibility constraint. This is done by making use of a surface Lagrange Fig. 4.6. Time evolution of the filament Γ h (t) for Re = Re Γ = 10 3 . The y-axis is magnified differently in each plot, there is no magnification in the last one. Different gray levels correspond to different successive instants, from dark to light gray. At t ≃ 0.5, small random perturbations are seen, which amplify and result in travelling waves of growing amplitudes for t 2. For t ≃ 6, the amplitude of these waves is of the same order as the filament length. multiplier (called tension), which allows to enforce accurately the inextensibility constraint locally on the thread or membrane. A variational formulation of this constrained problem is given, using a saddle-point approach in an analogous manner as for the incompressibility constraint. Here, however, the tension force associated with the Lagrange multiplier is a surfacic force applied on the inextensible boundary. Because this boundary may be curved, one finds that in addition to the tangential surfacic tension force, which is the surface gradient of the tension, a normal component arises which involves the local value of the tension and of the curvature. This differs with the case of bulk incompressibility, where only the gradient of the corresponding Lagrange multiplier (pressure) appears. The boundary conditions specific to the inextensible moving boundary are explicited in the context of the variational formulation of the saddle-point problem, they are in agreement with the physical consideration that at the endpoints of the boundary, either the tension (corresponding to a point force) or the tangential velocity must be imposed.

To take the best advantage of this method, which allows to enforce that the velocity on the moving boundary verifies the inextensibility constraint up to roundoff error, it is better not to interpolate the value of velocity over a non-matching bulk mesh, as would be done in an immersed boundary or a level-set method. Using a Lagrangian or an arbitrary Lagrangian-Eulerian (ALE) method allows to have a bulk mesh that includes the moving boundary's mesh in the set of its edges. Hence, there is no need for interpolation and the inextensibility force is actually treated as a boundary condition rather than approximated by a bulk force. This allows our implementation to yield sharp velocity profiles as in the example in Fig. 4.5. It is shown in this and other examples that the combination of the Lagrange multiplier approach and an ALE method based on finite elements yields a very accurate and highly robust numerical method, allowing to focus on stiff problems of fluid-structure interactions: in the swollen vesicle example, the combination of incompressiblity and inextensibility constraints maintain the shape of the vesicle like an inflated ball; in the 'flag' example, the structure interacts with a high-Reynolds number fluid flow. Some questions arise however which are not addressed within the scope of this paper. Indeed, it is well known that mixed methods require a compatibility condition between the finite element space of the Lagrange multiplier and the one of the original variable. Here, it is found empirically that, if the moving boundary forms a closed loop, or has both endpoints fixed by means of a Dirichlet boundary condition on the velocity, the surfacic Lagrange multiplier can be taken continuous linear by piece (P 1 (Γ h )) when the velocity is continuous quadratic by piece (P 2 (Ω h )). However, this is not the case anymore when at least one of the endpoints of the moving boundary is free to move, that is, when the boundary condition is imposed on tension rather than velocity. In that case, we find that a P 1 (Γ h ) test-space is too small to enforce inextensibility close to the endpoint of the moving boundary, as there exist discrete velocity fields which do not verify the (continuous) inextensibility constraint but are in the kernel of the corresponding discrete surface divergence operator. Employing a P 2 (Γ h ) test-space and filtering the resulting Lagrange multiplier to obtain the tension did circumvent this problem. It remains an open question whether in this case there exists an adequate choice of spaces that would verify the Babuška-Brezzi condition.

Finally, we would like to stress that there are many physical situations of interest in which the strong points of this method could prove decisive. Indeed, any penalty method implies to choose globally an extensibility modulus. Thus, regions of strong sollicitation will extend proportionately more in a penalty method, which will not be the case using a Lagrange multiplier approach. Yet, a physical situation as common and simple as a filament extended by a flow field, fixed by one of its ends, exhibits a sollicitation that goes to infinity at its free end. It is thus crucial to pay attention to this issue.

Fig. 4 . 1 .Fig. 4 . 2 .

 4142 Fig. 4.1. Boundary conditions for the swollen vesicle in Stokes flow problem. The flow imposed through the boundary condition counterbalances exactly the gravity force acting on the object.Through the saddle-point approach, pressure and surface forces are calculated that enforce the undeformability of the subdomain Ω d .

Fig. 4 . 3 .

 43 Fig. 4.3. Geometry and boundary conditions for the flow past a flag problem. The flow past a plate can be simulated when the curve Γ(t) is a straight segment aligned with ex.

Fig. 4 . 4 .

 44 Fig. 4.4. Comparison of the numerical solution for the tension ζ h (x, 0), velocity u x h (x, 0) and friction force ∂u x h ∂y (x, 0 + ) on the line y = 0 around point B = (1, 0) T , with the asymptotic behaviours ζ ∼ a √ 1x, u x ∼ a

Fig. 4 . 5 .

 45 Fig. 4.5. Approximation u x h of the velocity in the neighbourhood of the flag. Visualisation uses P 1 -iso -P 2 approximation.

  2 , for two (related) reasons: functions of H 1/2 (Γ) are not continuous, and the bilinear form (ξ, v t ) → Γ (∂ξ/∂s)v t dA is not continuous on H 1/2 (Γ) 2[5, page 32]. This is related to the fact that the derivatives of functions in H 1/2 (Γ) are in the dual of H

	1/2 00 (Γ) rather than in the dual of H	1/2 0 (Γ)

  Calculate the mesh advection velocity c n h using Eq. 3.2.

	e) Calculate mesh Ω n+1 h necessary.	by advection of Ω n h by ∆t c n h , remesh if

The case of a straight thread with both endpoints in Γ D being either trivial, or ill-posed if a normal force is applied on the thread.

Acknowledgments. J. E. wishes to thank Richard Michel for helpful discussions on some aspects of the functional analysis of this problem.

supported by ANR project MOSICOB "Modélisation et Simulation de Fluides Complexes Biomimétiques".