Aravind Rajan Ayagara
email: aravind-rajan.ayagara@univ-orleans.fr

André Langlet
email: andre.langlet@univ-orleans.fr

Ridha Hambli
email: ridha.hambli@univ-orleans.fr

Finite-Element simulations of Ribs subjected to dynamic three

Keywords: Blunt impact, Porcine ribs, HR-pQCT, Dynamic response, Split Hopkinson Pressure Bar, Finite Element simulations, LS-Dyna

come

Introduction

The bone fragments originating from a blunt thoracic impact pose a risk of penetration thereby causing complex medical conditions, for example, pneumothorax, hemothorax or sometimes even a fail chest leading to eventual death. Therefore, trauma to the chest is one of the main contributors to the development of medical complexity and even mortality among polytrauma patients [START_REF] Horst | Thoracic trauma now and then : A 10 year experience from 16,773 severely injured patients[END_REF]. A recent study on chest trauma proposed by [2] states that a chest trauma can be caused either by a blunt impact or a penetrating impact. As the projectile does not penetrate into the body, an blunt impact cannot be identified as easily as a penetrating impact. Therefore, the blunt trauma (a trauma to the chest by blunt impact) needs more attention and time for recovery [2], [3]. Another aspect of blunt trauma is that the injury mechanism does play an important role [START_REF] Burnside | Blunt thoracic trauma[END_REF], [START_REF] Fallouh | Blunt thoracic trauma[END_REF]. Among which road traffic accidents are the most commonly encountered injury mechanisms, which are important to be known prior to a radiological examination. Therefore, we are interested to understand the behavior of rib under blunt impact environment.

The advantages of FE simulations over expensive and time staking experimental procedures (e.g. human dummies) have drawn attention for years now. Several numerical human thorax models are have been developed : Isolated thorax models [START_REF] Schoell | Age-and sex-specific thorax finite element model development and simulation[END_REF], H-model [START_REF] Haug | Computational Models for the Human Body, Handbook of Numerical Analysis[END_REF], RADIOSS [START_REF] Arnoux | Radioss finite element model of the Thor dummy[END_REF], THUMS [START_REF] Maeno | Development of a Finite Element Model of the Total Human Model for Safety (THUMS) and Application to Car-Pedestrian Impacts[END_REF], and HUMOS [START_REF] Robin | HUman MOdel for Safety -A joint effort towards the development of refined human like car occupant models in[END_REF]. Despite their efficiency and frequent use, these FE models still have some constraints, such as : geometrical representation, use elements, constitutive law, and damage law. To understand the behavior of ribs under such environment, the mechanical properties, the geometrical representation, the constitutive behavior and progressive nature of damage must be given equal importance.

Quite a lot of researches can be found on the mechanical properties of bones and the properties influencing them through regressions. We can notice that some authors consider both cortical and trabecular bone are formed by same morphological tissue i.e. the tissue density of both bones is the same and the difference in mechanical behavior is due to the architectural arrangements. While other authors state that they both are of different morphological tissue. Therefore, leading to the use of different independent variables. These studies can be categorized in terms of the independent variable used. For example, the apparent density ρ app was considered by [START_REF] Carter | Compressive behavior of bone as a two-phase porous structure[END_REF], [START_REF] Gibson | Mechanical behavior of cancellous bones[END_REF] and [START_REF] Keyak | Correlations between orthogonal mechanical properties and density of trabecular bone : use of different densitometric measures[END_REF]. In contrast to ρ app , mineralization based parameters such as ash fraction α or apparent ash density ρ α were also considered as independent variables. Some researches regarding the later are [START_REF] Vose | Bone strength-its relationship to X-ray-determined ash content[END_REF], [START_REF] Keller | Predicting the compressive mechanical behavior of bone[END_REF]. One of the recent studies on the parameters influencing bone properties put forth by [START_REF] Hernandez | Influence of bone volume fraction and ash fraction on bone strength and modulus[END_REF] had considered both volume fraction V f and ash fraction α. They had stated the bone tissue density ρ tissue is influenced by its mineralization and had proposed regression formulae.

Apart from the aforementioned variables, there are variables that influence the mechanical behavior of bones. One of the most important variables influencing the mechanical behavior of bones is the strain rate ε i.e. bone behaves differently under different mechanical loading conditions. An increase in Young's modulus E and an increase in ultimate strength σ ult was reported by [START_REF] Mcelhaney | Dynamic response of bone and muscle tissue[END_REF], [START_REF] Crowninshield | Response of compact bone in tension at various strain rates[END_REF], [START_REF] Currey | The effects of strain rate, reconstruction and mineral content on some mechanical properties of bovine bone[END_REF], [START_REF] Wright | Tensile testing of bone over a wide range strain rates : effects of strain rate, microstructure and density[END_REF] and [START_REF] Hansen | Effects of strain rate on mechanical properties of human cortical bone[END_REF]. A more interesting fact is that the strain at fracture ε f decreases with an increase in strain rate [START_REF] Mcelhaney | Dynamic response of bone and muscle tissue[END_REF], [START_REF] Wood | Dynamic response of human cranial bone[END_REF].

Similar to the mechanical properties, bone had also been studied through FE simulations for years now. The most commonly implemented constitutive law for cortical bone is elastic-plastic [START_REF] Garcia | Elastic-plastic damage laws for cortical bone[END_REF], [START_REF] Pavier | Contribution à la compréhension des phénomènes physiques lors de l'impact d'un corps sur un modèle de structure biologique[END_REF]. We can see the use of elastic-plastic consitutitve law for trabecular bone in [START_REF] Li | Rib fractures under anterior-posterior dynamic loads : Experimental study and Finite Element Simulations[END_REF], [START_REF] Yates | Subject specific modeling of human ribs : Finite Element simulations of rib bending tests, mesh sensitivity, model testing with data from coupon tests[END_REF]. Yet in contrast, there is evidence of using foam like constitutive laws for trabecular bone [START_REF] Gibson | Mechanical behavior of cancellous bones[END_REF], [START_REF] Kelly | Experimental and numerical characterization of the elastic-plastic properties of bovine trabecular bone and a trabecular bone analogue[END_REF], [START_REF] Karkar | Personnalisation morpho-mécanique de la voûte crânienne humaine à différentes vitesses de sollicitations[END_REF] In this study, we propose an elasticviscoplastic constitutive law with isotropic hardening for cortical bone and an elasticviscoplastic constitutive law with kinematic hardening for trabecular bone. The constitutive law is coupled to an external path dependent damage law that is capable of considering non linear damage accumulation, reduction in rigidity through continuum damage mechanics, and strain rate effects on fracture strain.

Material and Methods

Experimental Procedure

The experimental tests were carried using a modified Split Hopkinson Pressure Bar (SHPB) setup for studying the behavior of material under dynamic three point bending. The principle of SHPB tests is to generate a compressive stress pulse through an impact of a striker (preferably of same material as bars) launched at the input bar. This compressive waves propagates (ε inc (t )) along the length of input bar. Once the wave arrives at input bar sample interface, a part of it is reflected and a part of it is transmitted to the output bar(s) through the sample. The amount of wave reflected (ε ref (t )) and transmitted (ε tr1 (t ) and ε tr2 (t )), depends on the difference in the mechanical impedance of bar and sample. Since the signals are acquired at a known distance form the interfaces, a time and space shifting procedure is required to differentiate these three waves. This includes the definition of arrival time for the incident wave pulse at the input bar-sample interface and then expressing all strain signals to the relative time (see Eq. 1) as such they were measured at the interface.

t = t -t 0 (1)
If the stress magnitudes in the bar are below yield limit of the bar material, then the plane wave formula is applicable. This formula leads to the dynamic variables used for describing the response of the sample i.e. Force at interfaces (F inp , F out1 , and F out2 ), velocity at interfaces (v inp , v out1 , and v out2 ), displacements at interface (u inp , u out1 , and u out2 ) expressed from Eq. 2 to Eq. 10

The porcine rib samples were brought from a local butcher and the surrounding muscles and tissues were removed prior to experiments. The experimental setup consists of a striker, an input bar (3 m long) and two output bars (1.75 m long). The bars of SHPB setup were made of nylon in order to adapt for the impedance difference between bars and rib sample. The anterior and posterior extremities of rib were sanded to establish a planar contact at the output bar-sample interface.

F inp = A b E b [ε inc (t) + ε ref (t)] (2) F out1 = A b E b ε tr1 (t) (3) F out2 = A b E b ε tr2 (t) (4) v inp = -c b [ε inc (t) -ε ref (t)] (5) v out1 = -c b ε tr1 (t) (6) v out2 = -c b ε tr2 (t) (7)
u inp = ˆt 0 v inp dt (8)
u out1 = ˆt 0 v out1 dt (9)
u out2 = ˆt 0 v out2 dt (10)

Numerical Rib Geometry

The numerical porcine rib sample was generated from the contours of DICOM images obtained through High Resolution-peripheral Quantitative Computed Tomography (HR-pQCT) scans. The porcine rib samples brought from a local butcher were soaked in saline solution, stored in a refrigerator at and, carefully defatted prior to CT scans. The raw DICOM image data were the imported to ScanIP module of Simpleware (Synopsys) for image processing, visualization and geometry generation of rib segments. A Five-sample calibration phantom with K 2 HPO 4 and H 2 O densities was placed below the ribs in the course of scan. The raw linear attenuation data µ were converted to normalized data GS, expressed in Hounsfield Units using the following relation

GS [HU] = 1000 µ b -µ w µ w -µ a ; µ w = 0; µ a = -1000; (11)
here, µ b , µ w , and µ a are µ values of bone, water, and air respectively.

In order to test the influence of different independent variables on the mechanical properties of bone, three different densities were calibrated :

1. Apparent Dry Density ρ app , determined from the hypothesis presented by [START_REF] Taylor | Determination of orhtotropic bone elastic constants using FEA and modal analysis[END_REF]. Firstly, the effective density ρ ef f (Eq. 12) 1 . Once ρ ef f is obtained, the ρ app was interpolated to GS through a linear relation as in Eq. 13

ρ ef f [kg.m -3 ] = 0.523GS + 1000 (12) ρ app [kg.m -3 ] = 1.207GS (13)
2. Equivalent Mineral Density ρ eqm . The ρ eqm was found from GS values through the following equation ρ eqm [kg.m -3 ] = 0.6452GS + 0.487 [START_REF] Vose | Bone strength-its relationship to X-ray-determined ash content[END_REF] 3. Apparent Ash Density ρ α . Two methods were followed obtain ρ α of porcine rib (a) form ρ eqm , through Eq.

ρ α (ρ eqm ) [kg.m -3 ] = 1.27ρ eqm + 0.0452 (15)
(b) based on the relation of ρ α with ρ app , proposed by [START_REF] Keller | Predicting the compressive mechanical behavior of bone[END_REF] as expressed as in Eq. 16 and Eq. 17

Mechanical properties of bone

As said in introduction, the mechanical properties of rib constituents were determined using either single variable power law regression (Model 1a, 1b and 1c) or two variable power law regression (Model 2). These models predict mechanical properties such as Young modulus (E), Yield strength in compression and tension (σ yC ) and (σ yT ), Ultimate strength in compression and tension (σ uC ) and (σ uT ) respectively.

Single variable models

These models consider only one independent variable either ρ app (model-1a) or ρ α (model-1b and model-1c).

Model-1a

The Model-1a is based on the hypothesis that ρ app is the only variable influencing the mechanical properties of bone. Since the porcine rib has both cortical and trabecular bone, corresponding porosity P o values were calculated as follows

V f app = ρ app ρ tissue (18)
P o app = 1 -V f app × 100 (19)
The ρ tissue used in Eq. 18 was considered as 1850 kg.m -3 (same for both cortical and trabecular bone), that is represented as ρ tissue literature in Fig.

Model-1b

In contrast to Model-1a, this model uses the apparent ash density obtained from [START_REF] Keller | Predicting the compressive mechanical behavior of bone[END_REF] ρ α (ρ app ) as independent variable. The mechanical properties were obtained though the equations given below

Model-1c

This model uses the apparent ash density obtained from ρ eqm as independent variable.The mechanical properties were obtained though the equations given below

E α(ρeqm) [GPa] = 3.284 × 10 -6 ρ 1.92 α (23) σ uC α(ρeqm) [MPa] = 4.399 × 10 -5 ρ 2.27 α (24) σ uT α(ρeqm) [MPa] = 3.079 × 10 -5 ρ 2.27 α (25)

Two variable model

The two variable model used in this study was based on [START_REF] Hernandez | Influence of bone volume fraction and ash fraction on bone strength and modulus[END_REF], that considers both the volume fraction V f and ash fraction α of the porcine rib sample. The assumption proposed by [START_REF] Hernandez | Influence of bone volume fraction and ash fraction on bone strength and modulus[END_REF] is the ash fraction and volume fraction of bone influence the mechanical properties, therefore, ρ tissue in this model was interpolated through Eq.26 as opposed to a constant value of 1850 kg.m -3 in Model-1a.

ρ tissue = ρ app 1.41 + 1.29α (26)
Since ρ tissue is not constant, the V f does not have the same evolution for different GS values as of Model-1a, thereby changing the P o values. The mechanical properties in this model were obtained through the following equations

E α,V f [GPa] = 84.37 V 2.58 f α 2.74 (27) σ uC α,V f [MPa] = 794.33 V 1.92 f α 2.79 (28)
σ uT α,V f [MPa] = 556.00 V 1.92 f α 2.79 (29)

Comparison of regression models

In order to make sure which hypothesis shall be considered for the material properties of bone, the results of single variable model and two variable model were compared and verified with the element distribution in the numerical rib sample.

First of all, the different density measures obtained through all these models were compared. The ash density of model-1b ρ α (ρ app ) is greater than that of model-1c ρ α (ρ eqm ). This is due to the sample origin used for correlation i.e. the data presented by [START_REF] Keller | Predicting the compressive mechanical behavior of bone[END_REF] was obtained from pooling data of human femur and human vertebrae, whereas, ρ α (ρ eqm ) of Model-1c was specific to porcine ribs. This shows the ash fraction varies for different species, thereby explaining the overestimation of density data.

Secondly, the volume fractions and porosity from Model-1a and Model-2 V f app ,P o app and V f (α), P o(α) were compared. This again proved that the ash fraction had a tendency of over estimating the values.

Lastly, estimated mechanical properties common to all 4 power law regression models (Young modulus, ultimate strength in compression and tension) were compared. Regarding the mechanical properties, following conclusions can be drawn

• The young's modulus obtained thorough apparent density and two variable model are of similar range and are largely high to those obtained from apparent ash density based models (see Fig. 2a).

• The ultimate properties are highly influenced by α or ρ α as it can be seen in Fig. 2b. Even though the model 2 considers α, it does produce results similar to those obtained from apparent density.

We can conclude that the ash content in bone has a strong positive correlation on strength of bone and a negative correlation on the initial elastic behavior of bone. The reason for this influence in particular for human ribs or its biological surrogates is yet to be explained. The properties obtained from apparent density are in accordance with that used by previous authors for human rib. The comparison in presented below curve plotted in Fig. 4a.

In this study, the static fracture strain ε f s value is 2 × 10 -2 m/m for cortical and 3 × 10 -2 m/m for trabecular bone. To approach the realistic fracture and propagation behavior using our damage model, the influence of strain rate ε on the fracture strain ε f d was implemented by an equation from [START_REF] Wood | Dynamic response of human cranial bone[END_REF] as follows :

ε f d ε f s = A -B Ln ε (32)
where A = 0.63. Five different values of B = 1.25 × 10 -2 . The corresponding ε f d /ε f s vs ε input curve is plotted in Fig. 4b. Note that the static fracture strain provided by [START_REF] Wood | Dynamic response of human cranial bone[END_REF] is 0.86 × 10 -2 m/m, which corresponds to a quasistatic strain-rate of 3.18 × 10 -3 s -1 . Similarly, in the present study Eq. ( 32) is not applicable for ε lower than 3.00 × 10 -3 s -1 , which are quasi static. For all values 0 s -1 ≤ ε ≤ 3.00 × 10 -3 s -1 , a linear approximation is applied so that

ε f d /ε f s = 1 at ε = 0 s -1 .
The increments in damage ∆D is expressed in terms of plastic strain ∆ε p and triaxility η as follows :

∆D = n ε f d (η, ε) D (1-1/n) ∆ε p (33)
The damage exponent n is to differentiate between linear and non-linear damage accumulation.

The reduction in load carrying capacity due to an increase in damage is taken into account through a modified equation of Lemaitre's effective stress principle :

σ = σ (1 -D m ) (34)
When the damage D reaches 1.0, fracture occurs and the Gaussian integral point of the corresponding element no longer exists in the calculation.

The damage exponent n values for cortical and trabecular bone were 1.25 and 2.0 respectively. Unfortunately, the lack of database on the fading exponent m for bone had led us to use m = 1.0.

When the yield stress is reached damage D starts to accumulate according to Eq. [START_REF] Fillipe | On the prediction of damage and failure in LS-Dyna : a comparison of GISSMO and DIEM[END_REF]. The fracture strain ε f d influences this accumulation. The more ε f d decreases with increasing strain-rate, more drastic will be the increase in damage accumulation ∆D and in ∆ε p .

The ultimate damage value D = 1 is reached through the non-linear coupling between damage D, fracture strain ε f d (η, ε), plastic strain ε p and softening. When D = 1, the corresponding element is removed from the structure, thereby initiating the crack propagation until fracture.

The strain rate influence on the post-yield behavior is also clear on the energy W absorbed by the structure. The more rapid is the decrease ε f d with respect to ε, less is the energy absorbed.

Results and Discussions

The results were analyzed in terms of force-time response and force-displacement response.

Experimental results

The forces in Fig. 5a were obtained through Eq. 2, Eq. 3 and Eq. 4 respectively. This response was classed into the category "Fracture before dynamic equilibrium" i.e. fracture of rib long before a dynamic equilibrium is established between F inp and (F out1 +F out2 ). The time taken by other porcine rib samples to establish dynamic equilibrium was ≈ 6.00 × 10 -4 s. The rupture of rib is characterized by a brutal decrease in force, occurs at ≈ 5.20 × 10 -4 s.

On the other hand, the structural response of rib was analyzed by force-displacement curves. Three different corridors were established with relative interface velocities and the Corridor-3 is shown in Fig. 5b which represent a group of force-displacement curves enclosed in maxima and minima of all curves of corresponding corridor.

Numerical results

The numerical responses were also obtained through the same shifting procedure of the numerical strain signal. The striker was given a velocity of 17.40 m.s -1 to simulate a similar experimental case of fracture before equilibrium. The numerical responses are shown in figure below compared to the corresponding experimental results. Following deductions can be drawn, 1. A keen observation to the ε inc and ε ref , reveals that there is a slight difference between numerical and experimental signals. On the other hand, the difference in the transmitted strain signals can be explained by the difference in geometry. Yet, the numerical interface velocities and displacements are almost similar to that of experimental.

2. The numerical displacement of output interface 2 shows a peculiar response. It is clear that it superposes with experimental signal till t ≈ 500 × 10 -4 s. Beyond this point of time the magnitude of u out2 Num is greater than u out1 Exp, that can be explained by the loss of contact between the sample and output bar interface 2. The increase in displacement beyond this time represents the rigid displacement of output bar. The numerical force-displacement response is shown in figure below. Initially, the structure behaves as an elastic material (till t ≈ 0.75 × 10 -4 s). During this period of time, the response depends on material properties and the stiffness distribution. Beyond this time, the force starts to increase till it reaches the maximum value and decreases till t ≈ 2.05 × 10 -4 s. This is followed by an increase in force, which represents the reflection of waves from output interface. The second decrease in force starts to occur when a crack appears in macro scale thereby proving the reduction in force carrying capacity of the rib until rupture at t ≈ 5.00 × 10 -4 s or U ≈ 5.00 × 10 -3 m.

F 1 -

 1 ρ α (ρ app ) [kg.m -3 ] = 0.522ρ app + 0.007; f or ρ app < 1000 kg.m -3 (16) ρ α (ρ app ) [kg.m -3 ] = 0.779ρ app -0.250; f or ρ app ≥ 1000 kg.m -3 Different density measures for porcine rib

E

  α(ρapp) [GPa] = 3.284 × 10 -6 ρ 1.92 α (20) σ uC α(ρapp) [MPa] = 4.399 × 10 -5 ρ 2.27 α (21) σ uT α(ρapp) [MPa] = 3.079 × 10 -5 ρ 2.27 α (22)

F 2 -

 2 Comparison of mechanical properties for porcine ribUnfortunately, there is a lack of data on the influence of α, ρ α or two variable models on σ y .The constituents of porcine rib, obtained from model 1a and model 2 were compared with the element distribution of numerical rib sample based on porosity values. It is clearly seen that Po app is in strict accordance with the element distribution in contradiction to those obtained by Po α,V f of model 2. This element distribution is shown in figure below and corresponding density and porosity values.

  distribution in porcine rib sample

F 4 -

 4 ε f d /ε f s vs ε Input curves of damage law

F 6 -

 6 Experimental and Numerical strain, interface velocities and displacements

  1. Corresponding power law regression equations available in literature were used for cortical and trabecular bone constituents based on P o app (see Tab. 1)

	Property	Cortical Bone	Trabecular Bone
	Young's Modulus E app , [MPa]	1.109 × 10 -6 ρ 3.09 app	0.5903 ρ 1.20 app
	Compressive Strength σ uC app , [MPa] 1.659 × 10 -4 ρ 1.88 app	8.534 × 10 -4 ρ 1.53 app
	Tensile Strength σ uT app , [MPa]	1.61 × 10 -4 ρ 1.88 app	8.201 × 10 -4 ρ 1.07

app Compressive Yield σ yC app , [MPa] 1.442 × 10 -4 ρ 1.85 app 5.167 × 10 -4 ρ 1.60 app Tensile Yield σ yT app , [MPa] 3.064 × 10 -4 ρ 1.67 app 7.5686 × 10 -3 ρ 1.04 app T 1 -Power law regression equations used for apparent mechanical properties (Model-1a)

density of the sample including the bone marrow and fluids

These data prove that either V f or ρ app based regression models are efficient and reliable for having a decent insight on mechanical properties given that it is correctly correlated to GS.

FE Simulations

The FE simulations presented in this study were carried out at experimental scale. An elastic constitutive model was chosen for bars. The striker was given an initial velocity of 17.40 [m.s -1 ].

A good constitutive law for bone material must consider the strain rate effects. As proven by [START_REF] Zioupos | Microcracking damage and the fracture process in relation to strain rate in human cortical bone tensile fracture[END_REF] that at higher strain rates, the strain rate has more effects on the post yield regime and the use of elasticviscoplastic constitutive law by [START_REF] Li | Rib fractures under anterior-posterior dynamic loads : Experimental study and Finite Element Simulations[END_REF] proves that the elastoviscoplastic constitutive law is adaptable for bone material. Therefore, a classic elastic plastic constitutive law was modified to consider the strain rate effects on post yield behavior through the cowper-symmonds model, that can be expressed as

here, σ yd is the strain sensitive yield stress, σ y0 is the initial yield stress, ε is the strain rate, C and P are cowper-symmonds parameters. The σ y0 is identified as the tensile yield stress presented in Tab. 2. The cowper-symmonds parameters are adapted from [START_REF] Li | Rib fractures under anterior-posterior dynamic loads : Experimental study and Finite Element Simulations[END_REF] and are specified here for the sake of clarity C = 2.5 GPa and P = 7.0. The cortical and trabecular bone were supposed to obey isotropic and kinematic hardening respectively as found in previous studies.

In order to simulate the progressive nature of fracture propagation and damage accumulation, an external damage law was coupled to the constitutive law. This law is incremental and stress-state dependent and is capable of considering the reduction in load carrying capacity and effects of strain rate. This model is based on the works of [START_REF] Frieder | Considering damage history in crash-worthiness simulations[END_REF] and [START_REF] Fillipe | On the prediction of damage and failure in LS-Dyna : a comparison of GISSMO and DIEM[END_REF]. The dependency on stress state is considered through stress triaxiality as follows

where : 1 3 tr σ is the hydrostatic pressure (average of the σ stress tensor trace) and σ eq is the equivalent stress.

The relation between the fracture strain ε f d and triaxiality η -Eq. ( 31) -is given in the form of an input

Conclusions

Dynamic three point bending tests using the split hopkinson pressure bar on rib was successfully implemented to have insights on fracture of ribs.

Different power law regression models were tested. The apparent density was the one that gave reliable results. The power law regression models with apparent ash density or ash fraction had a tendency to over estimate the ultimate properties of porcine rib. Moreover, the average apparent properties of porcine rib constituents are in good accordance with those of human ribs. This proves that the choice of porcine rib as human surrogates is biologically reliable.

The numerical results in this study prove that the FE model presented in this study is pertinent. The strain rate effects in behavior and damage till rupture were successfully taken into account. Even though the FE model results are coherent with that of experiments, optimization of damage law parameters through uni axial tensile test data could further improve the predicting capability of the model.


Aravind Rajan Ayagara

André Langlet

Gregory Girault

Fluid Structure Interaction between plate and blast waves : Numerical simulations

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Fluid Structure Interaction between plat and blast waves : Numerical simulations

Aravind Rajan Ayagara a , André Langlet b ,Grégory Girault c a. Laboratoire de Mécanique Gabiel Lamé, Universtié d'Orléans : aravind-rajan.ayagara@univ-orleans.fr b. Laboratoire de Mécanique Gabiel Lamé, Universtié d'Orléans : andre.langlet@univ-orleans.fr c. Institut de Recherche Dupuy de Lôme, Université de Bretagne Sud, Université de Brest, ENSTA Bretagne, Centre National de la Recherche Scientifique : gregory.girault@univ-ubs.fr

Résumé :

Cette étude présente un modèle numérique de l'interaction fluide-structure, capable de simuler la réponse d'une plaque à une explosion aérienne. Le calcul a été effectué en deux phases. Depuis l'explosion jusqu'à la plaque, la propagation des ondes de chocs est modélisée en 1D par l'analogue de ballon. Après l'impacte d'onde, le calcul est réalisé en 3D par projection des résultats 1D en 3D. Cette méthodologie est capable de reproduire la phénomène de réflexion des ondes, la formation de pied de mach. La plaque d'aluminium est modélisée avec la loi de comportement élastique-plastique. La répartition des contraintes sur la plaque dépasse la limite élastique. Donc, une loi de comportement couplée avec une loi d'endommagement pourrais donner des résultats précises et l'information sur l'évolution d'endommagement de la plaque.

Abstract :

This study, presents a numerical methodology for fluid-structure interaction between blast waves and plate using the LS Dyna explicit code. The numerical simulation was executed in two stages, (i) a one dimensional calculation was carried out using the balloon analogue and, (ii) the results of 1D calculation were remapped to three dimensions. The explosive used in this study was a mixture of propane and oxygen. The 3D model was able to reproduce the complex reflection phenomenon of shock waves, including the mach stem. A more refined mesh for the air domain shall produce accurate representation of this phenomenon. The aluminum plate was modeled using an elastic plastic constitutive law without damage. Results show that the stress in the plate are beyond its yield limit, therefore, an external damage law coupled to the constitutive law will further improve the results by giving us vital information on damage evolution and rupture.

Mots clefs : Air blast, Fluid-Structure Interaction, Mach Stem, MMALE, LS Dyna.

Introduction

The interaction of blast waves with a thin plate is investigated in this study. A blast wave originating from high explosives have a unique effect, known as the primary blast. This is due to the impact of over-pressurized waves on a structure and can be divided into 4 processes as follows 1. the detonation process of explosive 2. propagation of shock waves in the atmosphere around 3. reflection of shock waves from the obstacle (plate)

4. response of the plate and its consequences in the domain below the plate.

According to [START_REF] Kinney | Explosive Shocks in Air[END_REF] and [START_REF] Baker | Similarity methods in engineering[END_REF], the reflection of blast wave from a flat surface can be described accurately by mathematical equations (coalescence of incident and reflected waves at the triple point). One the blast waves impacts the structure, it resembles a fast moving discontinuity, followed by a continuously evolving non-uniformly distributed pressure.

There are two methods for simulation of blast loads (i) using empirical formulas and TNT energy equivalence and (ii) the compressed balloon analogue. The first method uses non linear equation of state such as JWL equation of state for the mechanics of explosion with an equivalent TNT energy and has been successfully implemented in [START_REF] Wang | Three dimensional ghost-fluid large-scale numerical investigation on air explosion[END_REF]. Due to the difference between gaseous and solid explosives, this methodology might lead to discrepancies in the results. Therefore, the compressed balloon analogue was chosen in this study. The principle of balloon analogue method is that the gaseous explosives are confined in balloon of certain radius and are separated from the ambient atmosphere by an fictitious membrane. The detonation is achieved by inducing a spark to the mixture set to the Chapmann-Jouget pressure. The successful application of this analogue can be seen in [START_REF] Benselama | A 1D-3D mixed method for the numerical simulation of blast waves in confined geometries[END_REF] and [START_REF] Langlet | Transient response of a plate-liquid system under an aerial detonation : Simulations and experiments[END_REF].

The explosion process is very complicated since, this process divides the fluid domain into in-homogeneous parts through a thin reactive zone and does induce large deformations. Due to this reason, the numerical simulations of explosion using classic numerical methods is very difficult. For which, the ALE or Eulerian Multi-Material formulations can be used as an alternative approach [START_REF] Langlet | Air blast reflecting on a rigid cylinder : simulation and reduced scale experiments[END_REF].

As said above, the explosive used in this study is a mixture of propane and oxygen (C 3 H 8 + 5O 2 ), which consists of two different elements. Therefore, a Multi-Material Arbitrary Lagrange-Euler (MMALE) formulation was chosen in this study. Thanks to the interface tracking algorithm in MMALE method, it is capable of tracking the interfaces between the explosive gas and the surrounding air. This cannot be achieved by ALE method. In the MMALE formulation, the FE mesh moves independently to the material flow and each element can allow the solver to track different materials. The working principle of MMALE is explained in detail in [START_REF] Aquelet | Euler-Lagrange coupling with dampening effects : Application to slamming problems[END_REF] and [START_REF] Langlet | Air blast reflecting on a rigid cylinder : simulation and reduced scale experiments[END_REF] and the equations are presented here for the sake of ease. Let us consider a variable f . It's time derivative in the referential coordinate system can be expressed as

df (X, t) dt = ∂f (x, t) ∂t + (v -w) • gradf (x, t) (1)
here, X is the Lagrangian coordinate, x is the ALE coordinate, v is the particle velocity vector and w is the velocity vector of the reference coordinate. Then the equations conservation of mass, momentum and energy in ALE coordinates can be written as

∂ρ ∂t + ρdivv + (v -w)gradρ = 0 (2)
ρ ∂v ∂t + ρ(v -w)grad v = div σ + f (3)
ρ ∂e ∂t + ρ(v -w)grad • e = σ grad v + f.v (4)
To solve the equations 2 to 4 in ALE form, there are two ways. The first is to solve fully coupled equations for fluid mechanics and the second is by using Operator-Split method. The first methods cannot handle multiple materials in the fluid domain, therefore, the later was chosen in this study. Each time step calculation in an Operator-Split method is divided into two phases, with an interface tracking algorithm intervention between each phase. In a nutshell, a time step in Operator-Split method consists of a Lagrangian phase, an intervention of interface tracking algorithm and the advection phase or transportation phase. The advection phase in this study was carried out through second order van-Leer algorithm.

Since, the fluid mechanics problem is compressible, we need an Equation of State (EOS) for the fluid domain. Apart from these, for the interaction between solid and fluid domains at the point of contact, we might also need the Constrained-Lagrange-in-Solid card of LS-Dyna.

Numerical Simulations

The model consists of an aluminum plate of 0.5 mm thickness, which is subjected to blast loads originating from an explosion at 283 mm above the plate. The explosive is a mixture of propane and oxygen and is enclosed in a hemisperichal soap balloon of 62.5 mm radius.

The numerical simulations were carried out using the LS Dyna explicit hydrodynamic code. The simulation was divided into two parts for the sake of time saving. The first part was calculation of pressure variation in one dimension and the second was remapping the one dimensional results into three dimensions. In fact, the three dimensional model was carried out on 1/4 th scale of real case.

To simulate explosions and propagation of shock waves, one must need equations of state for both explosive and air. There are different equations of state available in LS Dyna. Since, the compressed balloon approach followed here is applicable for gaseous mixture explosive, we had chosen a linear polynomial equation of state for explosive gas and also air. A linear polynomial equation of state is expressed as

follows p = C 0 + C 1 µ + C 2 µ 2 + C 3 µ 3 + C 4 + C 5 µ + C 6 µ 2 e (5)
The constants C 0 to C 6 are user defined and µ is expressed in terms of relative volume as follows

µ = 1 v r -1 (6)
If the gases are considered as perfect gas, then the constants of Eq.5 have to be considered as follows

C 0 = C 1 = C 2 = C 3 = C 6 = 0 and C 4 = C 5 = γ -1
Therefore, the pressure for an ideal gas using the linear polynomial equation of state can be rewritten in terms of specific heats γ as p = (γ -1)

ρ ρ 0 e (7)
ρ and ρ 0 are the densities at a given iteration and initial density respectively. The equations of state for T 1 -EOS properties for explosive gas and air

1D model

The 1D FE model consists only of explosive balloon and the ambient air of standalone distance, meshed with beam elements. The NULL material model was used for both air and explosive. The reason for using NULL material model is to avoid the calculation of deviatoric stress components. Here, r b is the radius of the balloon which was considered 65.5 mm and d s is the standalone distance 283 mm. Since the material model NULL requires only the initial densities for this problem, the densities of explosive gas and air were considered to be 1.474 kg.m -3 and 1.290 kg.m -3 respectively. Other than these constituents, the 1D model was also provided with the ALE module. A Multi-Material ALE group was declared for explosive and ambient air domains.

3D model

The three dimensional model consists of two domains (i) the ambient air above and below plate (see Fig. 2a) and (ii) an aluminum plate of 0.5 mm thickness (see Fig2b).

Boundary and Initial conditions

Since the model was considered to be axi-symmetric with symmetry boundary conditions along the XY and YZ planes and the other faces of air domain were provided with non-reflecting boundary conditions. Where as the dotted black lines in Fig. 2b, represent the clamped boundary condition of the plate.

The initial condition for the 3D model was the solution of 1D model.

Results

It is interesting to see the variation of pressure with respect to space and time. For example the red and green curves in Fig. 3a are recorded at explosion point and at z = d s /4 which are respectively the results of 1D calculation. The variation in pressure with respect to depth below the plate is shown in Fig. 3b.

The Fig. 4, shows the corresponding pressure contours in the atmosphere above the plate. We can notice that the simulation is capable of reproducing the complex reflection of shock wave including the incident wave (solid black line), the reflected wave (dotted black line) and the mach stem where incident and reflected waves meet (vertical black line) at t = 6.097 × 10 -3 s. Similarly, the pressure contours below the plate are presented in Fig. 5. These are represented in the form of Isosurfaces for the sake of clarity. The range of pressure was fixed so that the wave propagation below the plate could be seen clearly. Note that the white blank space in Fig. 5 is due to the pressure values out of the range presented in the figure. The variation of pressure along the length of plate about X axis is shown in Fig. 6. The curves in solid line represent the pressure above the plate whereas the curves in dots represent the counterpart below the plate. The positive peak pressure in the curves below the plate represent the propagation of load front in the fluid below due to the deflection of plate. Meanwhile, the variation of von-mises stress in the plate is presented in Fig : 7. We were able to notice the reflection of waves in the plate from the fixed (clamped) end. This can clearly be seen between time instances t = 8.097 × 10 -3 s and t = 1.000 × 10 -3 s i.e. till the end of calculation.

Conclusions

We were able to successfully implement the balloon analogue in the LS Dyna explicit code. Consideration of linear polynomial equation state for a perfect gas gives coherent results. The use of ALE1D for one dimensional calculation and remapping the results for three dimensional from 1D results is a time effective method. The three dimensional model with current mesh size was able to reproduce the complex reflection phenomenon of shock wave including the mach stem formation. The current configuration had proved that the plate starts to yield, therefore if the model is coupled with an external damage law will produce the evolution of damage and rupture of plate if possible.

explosive and ambient air are given below Domain C 4 =

 4 C 5 = γ -1 e Explosive 0.2715 16 MJ.m -3 Air 0.4000 0.0 MJ.m -3

F 1 -

 1 1D Model with explosive (red) and ambient air (blue)

F 2 - 2 -

 22 (a) Fluid domain above and below plate (b) Plate submerged in fluid domain Fluid and Solid domains in the 3D model The fluid domain were meshed with 8 node hexahedral solid elements whilst the plate was meshed with shell elements. For the sake of accuracy, fully integrated shell elements with 5 integration points were considered. The plate was modeled using elastic-plastic law with following properties Property value Density ρ Kg.m -3 2780.00 Young's Modulus E [GPa] 73.00 Yield Strength σ y [MPa] Properties of aluminum plate

F 4 -

 4 Contours of pressure above the plate in MbarF 5 -Iso Surfaces of pressure below the plate in Mbar

2 F 6 -

 26 Pressure curves at different X on either side of plate F 7 -von-mises stress contours on plate


Anne-Laure Bernard
email: albernard@recherche.loreal.com.

Marie-Alice Guedeau-Boudeville

Olivier Sandre

Serge Palacin

Jean-Marc Di Meglio

Ludovic Jullien

A.-L Bernard

J.-M Di Meglio

Permeation through Lipid Bilayers by Adhesion of Giant Vesicles on Decorated Surfaces

Introduction

Controlling the transport of water-soluble solutes through lipid bilayers is a major issue in developing new approaches in drug targeting and delivery. Many strategies have already been designed to enhance the transmembrane passage of species as small as ions or as large as proteins. Chemical strategies generally make use of the addition of a suitable agent (detergent, carrier, channel, ...) whose conformation and favorable partition between the aqueous phase and the bilayer is responsible for the increase of permeability of the lipid barrier. In these strategies, transport is expected to occur in a steady state as long as the promoting factor is present in the system. Despite encouraging results obtained by acting upon the permanent endocytotic processes for transporting the largest solutes such as DNA, 1,2 chemical strategies are most often restricted to the transport of small species. [START_REF] Cevc | Phospholipid handbook[END_REF] Physical strategies consist in applying an external perturbation (temperature rise, [START_REF] Lasic | Liposomes: from Physics to Applications[END_REF] pulses of electrical field, [START_REF] Zhelev | [END_REF] ...) which leads to the formation of large permeating defects or pores in the bilayer. Unfortunately, together with often severe practical requirements, the physical approaches are intrinsically limited by the relaxation time of return to the nonpermeable equilibrium state via lipid migration or flip-flop processes for instance. Indeed, this characteristic relaxation time determines a kinetic window for permeation, in terms of both its extent and duration; the permeation occurs under kinetic control in analogy with the kinetic/thermodynamic controls of wide use in chemistry. 6 In the continuity of our previous studies devoted to the setup of chemical 7 and physical 8,9 strategies for inducing solute transport through lipid bilayers, we have been concerned with a new approach making use of a kinetic control that should lead to the transport of the largest solutes. Its basic principle lies in inducing the formation of transient pores in the membrane of giant vesicles 10,11 adsorbed on suitably decorated surfaces interacting with the membrane.

In this work, the interaction between vesicles and substrates originates from the strong and fast 12 adhesion of a charged bilayer onto a surface of opposite charge. Anticipating that a contrast between alternating more or less adhesive elements of surface should facilitate permeation, attention has been focused on decorated surfaces with a surface fraction x S of the more adhesive portions (x s ) 0.5, Figure 1). We have examined the significance of the scale and of the geometry of the surface patterns. In view of the bilayer thickness (a few nanometers) and of the typical radius of giant vesicles (10-100 µm), patterns with periodicities ranging from 1 to 50 µm have been chosen. Hence, the characteristic size of the pattern should be (i) large enough to promote large morphological changes of the vesicle and (ii) small enough to get a permeation that does not depend on the size and position of the vesicle. Stripe, square, triangle, and herringbone geometrical shapes have been retained. In addition to promote permeation, we are also interested in highlighting (i) the dynamics of vesicle spreading onto decorated surfaces with respect to the results obtained on smooth surfaces. 13 Despite the extensive use of vesicles for preparing supported bilayers and the numerous theoretical reports on the equilibrium morphology of vesicles interacting with surfaces, 14 the dynamic aspect has been rarely studied, 15 (ii) the most favorable patterns for promoting (or hindering) vesicle adhesion by changing the characteristic lengths and symmetries with regard to the thermal fluctuations of the vesicles, 16,17 and (iii) any vesicle motility driven by Marangoni effect on suitably decorated surfaces. 18

Experimental Section

Giant vesicles have been prepared by swelling a thin film obtained after evaporation of freshly prepared chloroform/ methanol 9:1 solutions of EPC (egg phosphatidylcholine, Sigma-Aldrich) in 50 mM sucrose aqueous solutions in the presence of an alternating electrical field. 19 To induce the sedimentation onto the surface, the vesicles are transferred at constant osmolarity into a less dense aqueous solution (50 mM glucose). 20 It has been experimentally observed that the trace amounts of negatively charged species resulting from EPC hydrolysis (evaluated at about 0.1% mol/mol 21 ) are sufficient to promote vesicle adhesion to the positively charged surfaces. The decorated surfaces have been obtained using the microcontact printing technique introduced by Whitesides. 22 PDMS (poly(dimethylsiloxane)) patterned stamps impregnated with a 1 mM ethanol solution of dodecanethiol were applied onto ITO (indium tin oxide)-coated glass surfaces covered by a thin optically transparent gold layer. After removal of the stamp, the resulting surface was stepwise dipped for 1 min into a 1 mM water/ethanol 1:9 solution of sodium 3-mercaptopropane sulfonate, then cautiously rinsed with Millipore water and finally immersed for 24 h in a 0.01% aqueous solution of positively charged poly-L-lysine (MW ) 150000-300000, Sigma-Aldrich). The decorated surfaces are then immediately used after extensive washing with Millipore water. Selective etching processes and atomic force microscopy (AFM) imaging in friction mode have been used to characterize the decorated surfaces. They have shown that the patterns have been correctly transferred from the PDMS stamps to the gold surface, producing alternating hydrophilic and hydrophobic areas on the surface. The patterns have been also directly observed by reflection interference contrast microscopy (RICM). 23 The spreading of the giant vesicles onto the decorated surfaces as well as induced bilayer permeability have been examined using a custombuilt microscope working in either phase contrast (to determine the vesicle radius) or evanescent wave induced fluorescence (EWIF). 13,24 In the latter case, a suitable fluorescent water-soluble dye (calcein) was added to the external vesicular pool; it becomes thus possible to visualize the contact zone between the spreading vesicle and the decorated surfaces. EWIF microscopy also allows the determination of the shape of the adsorbed vesicle at the micrometer scale. 13,24 In addition, one can evidence any permeability of the membrane toward the fluorescent probe by analyzing the contrast of the vesicle (dark) with regard to the background (bright). Complementary experiments have been performed with RICM or with epifluorescence microscopy using a lipophilic aminostyryl pyridininium dye (S-467, Molecular Probes). These techniques enable study of the topology of the bilayer with a high resolution.

Results

Figures 2, 3, 4, and 5 display the typical behavior of an EPC giant vesicle spreading onto surfaces patterned with stripes (Figure 2), herringbone (Figure 3), squares (Figure 4) and grids (Figure 5) exhibiting different attraction strengths for EPC vesicles. The initial time is taken as the moment when the vesicle sedimenting toward the decorated surfaces stops its descent. From this time, three different regimes have been evidenced. The longest regime (regime I; characteristic time τ 1 in the minute range) is associated with the formation of a first contact point between the bilayer and the surface. By comparison of the behavior on stripes of different widths (2.5, 5, and 12 µm), τ 1 is shown to increase when the surface fraction of attractive domains x S is reduced. Nevertheless, τ 1 remains in the same range of values as observed for smooth surfaces under comparable conditions. 13 By contrast to the case of smooth surfaces, the second regime reveals by EWIF microscopy a long sequential process of spreading (regime II, characteristic time τ 2 in the range of 10 s to be compared with 50 ms on smooth surfaces). The vesicle first undergoes a local spreading; for instance, it spreads rapidly along a stripe (Figure 2B). Then, one observes spreading along a neighboring attractive stripe (Figure 2C). This growth process occurs repeatedly until part of a homogeneous adhering surface of the vesicle is observed (Figures 2E, 3E, and4B). This process looks analogous to the Haines jump of a liquid on a corrugated surface. 25 The adhered part of the vesicle appears homogeneously dark in EWIF microscopy, and the intensity of the vesicle image does not evolve during the second regime. The vesicle shape remains circular as observed by phase-contrast microscopy. Then, the vesicle pauses for a period equal to a few seconds on square and grids patterns (Figures 4 and5) or to 1 min on striped and herringbone patterns (Figures 2 and3). The last regime (regime III, characteristic time τ 3 in the minute range) begins with the formation of a pore in the lipid bilayer along the contact line. In EWIF microscopy, the pore is evidenced by the emission of a diffuse dark cloud (see arrow in Figures 2F, 3F, 4C, and5E). Finally one notices spreading of the lipid bilayer along the attractive zones at a constant velocity of about 20 µm's -1 . Epifluorescence microscopy has been used to demonstrate that the vesicle parts spreading along attractive lines are not simply adsorbed bilayers but true tubular extensions (pseudopods) of the vesicle undergoing a Rayleigh-type instability (Figure 6A). 26 In view of the difference in luminosity between attractive and nonattractive stripes, EWIF microscopy confirms the preceding observation. Indeed, the tubes lower the fluorescence intensity above the attractive stripes in EWIF microscopy (Figure 6B). EWIF also shows that the vesicle/background optical contrast drops during regime III. At the end of the spreading process, the vesicle adopts a shape that depends on the surface pattern. In the case of the 2.5 µm wide striped pattern, phase contrast microscopy shows that the vesicle is elongated along the stripe axis (Figure 6C).

The spindle adopted morphology is confirmed by RICM (Figure 6D) and epifluorescence microscopy (Figure 6E).

A typical contact angle of the vesicle along the in-plane perpendicular axis to the stripes equal to 15-20°was extracted from the spacing between the interference fringes observed by RICM. Among the different regimes observed, regime III is the only one to exhibit scale-and geometry-dependent features. Within the series of striped patterns, regime III is only present for patterns with the 2.5 µm periodicity. For patterns made of larger stripes at the same value of x S , (i) the spreading phenomenon stops at the end of regime II (regime III does not occur) and (ii) no significant solute permeation can be observed and the vesicular shape as observed by phase-contrast microscopy remains circular. Note that under the experimental conditions investigated, the vesicle never blows up on stripe patterns. In contrast, regime III was always present on the other decorated surfaces examined (triangles, squares, grids, herringbone) and led to solute permeation. Moreover, the spreading vesicles could not be observed anymore by phase contrast, and it is thus suggested that the latter geometries promote vesicle rupture.

No motion of the center of mass of vesicles has been ever noticed during experiments, even when a driving force was expected (large triangle and square patterns for example). But an oscillating behavior has been frequently observed when aged patterned stamps of 2.5 µm stripes were used instead of the new stamps that were used to obtain the results reported in the preceding paragraph. Under such conditions, EPC vesicles hardly spread during regime III but one observes the regular opening and closing of a pore formed on the contact line at the border between zones exhibiting different attractive interactions (Figure 7). At the end of this beating step, one notices a loss of vesicle contrast in EWIF microscopy which is evidence of the internalization of the fluorescent probe.

Discussion

Our observations emphasize the complexity of the spreading behavior of giant vesicles on adhesive surfaces. Similar to the case of smooth surfaces, 13 a sequential dynamic process with three regimes has been observed (Figure 8).

Regime I can be described as the approach regime. After the sedimentation process controlled by buoyant and viscous forces, the vesicle motion considerably slows down as soon as the distance between vesicle and surface becomes small enough (comparable to the membrane fluctuations). The vesicle then appears immobile in phase contrast microscopy. The corresponding scale τ 1 is the longest time encountered during the entire adhesion process. In the case of undecorated surfaces which sufficiently attract the lipid bilayer, 13 spreading begins as soon as a contact between the membrane and the surface is established. In the case of decorated surfaces, one also has to take into account the probability of creating a nucleation point whose stabilizing energy exceeds the energy of thermal fluctuations (k B T). In the range of scales investigated, this probability is governed by the surface fraction of attractive zones. Nevertheless, it is reasonable to anticipate that some more complicated dependence of τ 1 on the pattern geometry of the decorated surface may exist but at a smaller length scale than that investigated in the present study. A crude estimate of the typical length scale at which such effects should manifest themselves can be obtained as follows. In the preceding model used to derive τ 1 , the assumption of a planar bilayer has been used. In fact, a fluctuating vesicle at a distance h of a planar surface can be described as the sum of independent elements whose characteristic lengths are ξ | ) (κ c /kTh) 1/2 and ξ ⊥ ≈ h, where κ c designates the Helfrich bending modulus. 27 If the surface is made of zones in the ξ | scale range (about 600 nm in the present case) and whose symmetry is more or less commensurate to the bilayer excursions, one should be able to observe geometrydependent effects.

Once the nucleation point is created, fast spreading of the vesicle onto the surface indicates the onset of regime II. It occurs at constant surface and volume of the vesicle. In the case of anisotropic patterns such as stripes, the spreading develops faster along the stripe axis than along the corresponding in-plane perpendicular axis. Indeed, although spreading is limited in both cases by the drainage of a thin aqueous film, a zip process takes place along the stripe axis whereas jumps are required along the perpendicular direction. 25,28,29 Despite the expected difference of kinetics and in marked contrast with regime III, the shape of the vesicle remains spherical at the end of regime II. It is also interesting to note that in the time scale investigated in regime II, the absence of optical contrast in EWIF microscopy between the stripes of different attractiveness for the bilayer underlines the absence of Spreading of an EPC vesicle on grid patterns (width of the attractive square grid, 5 µm; side of the square, 50 µm) as a function of time. Vesicle (50 µm diameter) as observed by phase contrast microscopy at t ) 0 (A). Vesicle as observed by EWIF microscopy (the depth of penetration of the evanescent wave is estimated between 600 nm and 1 µm) at t ) 0.3 (B), 2.7 (C), 24.5 (D), 25.3 (E), 28.1 (F), 34.9 (G), 45.9 (H), and 47 s (I). The arrow on picture E emphasizes the formation of a macroscopic pore along the contact line. differences in the thicknesses of the aqueous film located between the surface and the bilayer. Due to the strong adhesion energy, a rather large tension is exerted on the bilayer at the end of regime II. 30 Regime III starts with the release of the corresponding tension by the formation of a macroscopic pore with a radius r that can be estimated as follows. The volume of the nonfluorescent solution expelled from the vesicle between t ) 67 s and t ) 68 s (panels F and G of Figure 2) corresponds to a volumic flow J of about 10 4 µm 3 's -1 (assuming a spherical dispersion of the vesicle contents). The Laplace pressure drop responsible for this flow J is where Σ is the membrane tension and R, the vesicle radius that can be deduced from Figure 2A. If we consider that the flow J through the pore is limited by the viscous dissipation in the aqueous medium (this is justified by an a posteriori Reynolds number that gives Re ≈ 10 -3 ), we have where η is the solution viscosity. We have assumed that the pressure drop ∆P relaxes on a length r. Taking Σ ) 5 × 10 -3 mN/m, 13 R ) 40 µm (Figure 2A) and η ) 10 -3 Pa's, we obtain a pore radius r ≈ 5 µm. This estimate compares with the repeat distance of the decorated pattern and is larger than the recent observed radii of the membrane pores formed under different conditions. 31 Because of the leakage caused by the pore formation, an excess of area becomes available. Buds form and tubular pseudopods now spread. 32 In the case of the 2.5 µm periodic stripe pattern, the spreading stops once enough tension is released. In contrast, for some other patterns (triangles, squares and herringbone), the spreading vesicles blow up before releasing enough tension to close the pores.

To propose a model for the dynamics of vesicle spreading during the regime III on all the decorated surfaces tested experimentally is out of the scope of the present paper. Nevertheless, we bring semiquantitative arguments for the stripped pattern. In ref 33, the authors study the profile of an adhering vesicle on a solid plate in the vicinity of a pinning center. 33 The latter exerts locally a force f 0 on the bilayer. To make the analogy with the adhesive stripes, we write f 0 ) Wa, where W is the adhesion energy per unit area and a is the stripe width. The membrane resists against the deformation with an effective line tension τ given by the relation (see Figure 1) where κ c is the bending modulus, R c is the radius of curvature close to the adhesion area, and λ is the angular correlation length that was first introduced in ref 34. It arises from the requirement for the membrane to merge smoothly with the surface (in contrast with a liquid droplet that makes a sharp wedge). This length is written as λ ) (κ c /Σ) 1/2 and can be measured by RICM. 35 By estimating (30) In the case of smooth surfaces, the adhesion energy W has been evaluated to (1-5) × 10 -3 mJ/m 2 from the consideration of the charge densities of the lipid bilayer and of the surface. Such a value is generally considered to correspond to a strong adhesion. The corresponding surface tension Σ is related to W and to the contact angle of the vesicle on the surface, θ, by the Young-Dupre ´equation (see: Seifert, U.; Lipowsky, R. Phys. Rev. A 1990, 42, 4768. Nardi, J.; Bruinsma, R.; Sackmann, E. Phys. Rev. E. 1998, 58, 6340.): W ) Σ(1 -cos θ). In the present system, the bilayer surface tension at the end of regime II is evaluated to 5 × 10 -3 mN/m from both experimental observations and theoretical predictions. See ∆P ) 2Σ/R

(1)

J ≈ ∆Pr 3 /η (2) τ ) 1 2 κ c R c 2 λ (3) also R c ≈ λ/θ c
, where θ c is the contact angle, one obtains As sketched in Figure 1, the bilayer adopts locally a triangular shape as long as 2τ > f 0 . In contrast, if the bilayer cannot oppose a sufficient line tension to the applied force, an instability was predicted: 33 a tether is pulled from the vesicle. In our experiments, Σ ≈ 5 × 10 -6 N'm -1 , 30 θ c ≈ 15°, W ) Σ(1 -cos θ c ) ≈ 10 -7 N'm -1 , κ c ≈ 4 × 10 -20 J, 36 λ ≈ 10 -7 m, and 2τ ) 2 × 10 -14 N. The latter estimate is 1 order of magnitude smaller than the pulling force f 0 ) 3 × 10 -13 N. Thus we conclude that the growth of cylindrical protrusions from a vesicle can effectively be driven by spreading of the bilayer onto the charged stripes of the substrate. Since the initially floppy vesicle becomes tense and adopts a spherical cap shape at the end of regime II, the growth of the tethers during regime III would lead to dramatically increase the membrane tension Σ due to the corresponding increase of the contact area. Hence if there was no mechanism to release the membrane tension, the tether instability would stop immediately. In fact the process of vesicle through experimentally observed pores keeps long enough a low surface tension to create long protusions.

In ref 31, only vesicles that slowly release a highly viscous solution can bear transient pores as large as those estimated in the present paper (in a low viscosity medium, both the outflow and the tension release occur too fast to let the pores grow). A possible explanation for the present observations is that the pores nucleate at a larger initial radius r c on the decorated surface, where the pattern provides germs above 1 µm. The nucleation radius is also strongly dependent on the line energy of the bilayer (γ ≈ 10 -11 N for an EPC bilayer) since r c ) γ/Σ. 37 We have invoked the presence of minor amounts of some hydrolysis products in our phospholipidic vesicles to account for their strong adhesion on oppositely charged surfaces. 13 Among them, the water-soluble lysolecithins with a single hydrocarbon tail decrease γ by a factor of 2 at concentrations as low as 3 × 10 -4 mg/mL. 38 Thus r c ≈ 1 µm could be expected. Then, the chemical heterogeneities of the solid surface could easily nucleate pores.

Conclusion

The present work demonstrates that the adhesion of giant vesicles onto decorated surfaces can lead to the permeation of large hydrophilic solutes through lipid bilayers. The observations made by EWIF microscopy on giant vesicles directly demonstrate both the internalization of calcein (size ≈ 1 nm) and the creation of a transient macroscopic pore (diameter ≈ 1 µm) as a result of the spreading process. In addition to these new and interesting features, some selectivity for transport might be expected. In fact, the size and the lifetime of the pore could probably be tuned by conveniently tailoring the size and the geometry of the patterns as well as the relative surface fractions x s of the decorated surfaces.

Figure 1 .

 1 Figure 1. Forcing permeation through lipid bilayers by membrane adhesion on decorated attractive surfaces. Striped patterns have been chosen for the illustration with xs ) 0.5.From left to right: (A) nonadhesive regime, the vesicle stays just above the surface, the bilayer is not stretched and exhibits a low permeability toward water-soluble species; (B) dynamic regime of spreading, the vesicle spreads out onto the surface and adapts its shape to decrease the system energy, in this regime the bilayer becomes strained and its permeability toward hydrophilic solutes is increased as a consequence of the creation of microscopic defects or of macroscopic pores; (C) regime of equilibrium. Several mechanisms (lipid migration, lipid flipflop, ...) lead to the partial relaxation of the membrane tension. The defect density drops and the pores close.

( 26 )

 26 Bar-Ziv, R.; Moses, E. Phys. Rev. Lett. 1994,73, 1392.

Figure 2 .

 2 Figure 2. Spreading of an EPC vesicle on striped patterns (width ) 2.5 µm) as a function of time. Vesicle as observed in phase contrast microscopy at t ) 0 (A). Vesicle as observed in EWIF microscopy (the depth of penetration of the evanescent wave is estimated between 600 nm and 1 µm) at t ) 1 (B), 3 (C), 4 (D), 15 (E), 67 (F), 68 (G), 73 (H), 84 (I), and 92 s (J). The arrow on photograph F emphasizes the formation of a macroscopic pore along the contact line.

Figure 3 .

 3 Figure 3. Spreading of an EPC vesicle on herringbone patterns (width ) 50 µm) as a function of time. Patterning stamp as observed by phase contrast microscopy (A). Vesicle as observed in phase contrast microscopy at t ) 0 (B). Vesicle as observed in EWIF microscopy (the depth of penetration of the evanescent wave is estimated between 600 nm and 1 µm) at t ) 1 (C), 6 (D), 31 (E), 31.3 (F), 32 (G), 33 (H), 35 (I), 37 (J), 39 (K), 40 (L), and 43 s (M). The arrow on picture F emphasizes the formation of a macroscopic pore along the contact line.

Figure 4 .

 4 Figure 4. Spreading of an EPC vesicle on square patterns (width of the square grid, 5 µm; side of the attractive square, 50 µm) as a function of time. Vesicle (diameter 50 µm) as observed in phase contrast microscopy at t ) 0 (A). Vesicle as observed by EWIF microscopy (the depth of penetration of the evanescent wave is estimated between 600 nm and 1 µm) at t ) 5 (B), 30 (C), 34 (D), and 35 s (E). The arrow on picture C emphasizes the formation of a macroscopic pore along the contact line.

Figure 5 .

 5 Figure 5.Spreading of an EPC vesicle on grid patterns (width of the attractive square grid, 5 µm; side of the square, 50 µm) as a function of time. Vesicle (50 µm diameter) as observed by phase contrast microscopy at t ) 0 (A). Vesicle as observed by EWIF microscopy (the depth of penetration of the evanescent wave is estimated between 600 nm and 1 µm) at t ) 0.3 (B), 2.7 (C), 24.5 (D), 25.3 (E), 28.1 (F), 34.9 (G), 45.9 (H), and 47 s (I). The arrow on picture E emphasizes the formation of a macroscopic pore along the contact line.

Figure 6 .

 6 Figure 6. Typical shape of an EPC vesicle spreading on striped patterns (2.5 µm width) of contrasted attraction as observed in epifluorescence microscopy (A), EWIF microscopy (B), phase contrast microscopy (C), RICM (D); projection of 30 fluorescence cross sections in epifluorescence microscopy; please note the curvature of the membrane induced by the pattern (E).

( 36 )

 36 Evans, E.; Rawicz, W. Phys. Rev. Lett. 1990, 64 (17), 2094. (37) Taupin, C.; Dvolaitzky, M.; Sauterey, C. Biochemistry. 1975, 14 (21), 4768. (38) Chernomordik, L. V.; Kozlov, M. M.; Melikyan, G. B.; Abidor, I. G.; Markin, V. S.; Chizmadzhev, Y. A. Biochim. Biophys. Acta 1985, 812, 643.

Figure 7 .

 7 Figure 7. EPC vesicle exhibiting a periodic opening-closing process of a pore during spreading on aged 2.5 µm stripe patterns. Vesicle observed in phase contrast microscopy before adhesion (A). Same vesicle observed in EWIF microscopy after adhesion (B, C). The periodic opening-closing process (0.3 s duration, 3 s periodicity) allows the leakage of the fluorescent inner medium of the vesicle (D). Vesicle observed in phase contrast microscopy at the end of the phenomenon (E).

Figure 8 .

 8 Figure 8. Spreading/permeation mechanism of an EPC vesicle on striped patterns (2.5 µm width) as a function of time (from left to right): (1) regime I, approach; (2 and 3) regime II, spreading at constant lipid surface but increased contact area and constant volume, no permeation; (4) regime III, formation of a macroscopic pore, spreading at constant lipid surface, the internal volume decreases; (5) equilibrium vesicular shape, the macroscopic pore is closed. (A) Vesicle image in phase contrast microscopy. (B) Vesicle image by EWIF microscopy. (C) Vesicle profile.

  (13) Bernard, A.-L.; Guedeau-Boudeville, M.-A.; Jullien, L.; di Meglio, J.-M. Strong Adhesion of Giant Vesicles on Surfaces: Dynamics and Permeability. Langmuir 2000, 16, 6809.

(14) Seifert U.; Lipowsky, R. Langmuir 1991, 7, 1867. Seifert U.; Lipowsky, R. Mol. Cryst. Liq. Cryst. 1991, 202, 17. (15) Baumga ¨rtner A.; Bhattacharjee, S. M. J. Chem. Phys. 1997, 107, 4390. (16) Brochard F.; Lennon, J.-F. J. Phys. (Paris) 1975, 36, 1035. (17) Singhvi R.; et al. Science 1994, 264, 696. (18) Durand, I.; Jo ¨nson, P.; Misbah, C.; Valance, A.; Kassner, K. Phys. Rev. E 1997, 4, 56. (19) Angelova, M. I. Mol. Cryst. Liq. Cryst. 1987, 152, 89. (20) Gravity does not here play a significant role on the permeation processes. See for instance: Kraus, M.; Seifert, U.; Lipowsky, R. Europhys. Lett. 1995, 32, 431. (21) Pincet, F.; Cribier, S.; Perez, E. Eur. Phys. J. 1999, 311, 127. (22) Whitesides, G. Angew. Chem., Int. Ed. 1998, 37, 550. (23) Ra ¨dler, J., Sackmann, E. J. Phys. II 1993, 3, 727. (24) Bernard, A.-L.; Guedeau-Boudeville, M.-A.; Jullien, L.; di Meglio, J.-M. Europhys. Lett. 1999, 46 (1), 101. (25) Haines, W. G. J. Agric. Sci. 1930, 20, 97.

Author manuscript of article published in Langmuir 2000, 16, 6801-6808

Langmuir, Vol. 16, No. 17, 2000 Bernard et al.

Langmuir, Vol. 16, No. 17, 2000 Bernard et al.

Langmuir, Vol. 16, No. 17, 2000 Bernard et al.

Langmuir, Vol. 16, No. 17, 2000 Bernard et al.

Acknowledgment. This study has been supported by L'Ore ´al. We greatly acknowledge Professor Sackmann for the access to RICM facilities and Dr. V. Marchi-Artzner for her assistance during preliminary RICM measurements. We thank Dr. Vanlerberghe for fruitful discussions. DOI: 10.1021/LA991342P

Author manuscript of article with DOI 10.1021/la991342p American Chemical Society -Published on Web 07/22/2000


Nicolas Poirier

C Gandini

Spatial dynamics of settlement pattern and cultivated areas from

Keywords: Settlement, agrarian areas, spatial statistics, time-space dynamics, comparative approach

come   L'archive ouverte pluridisciplinaire

Introduction

The French Archaedyn project1 , directed by F. Favory and L. Nuninger, is dedicated to the study of territorial dynamics in the long term, especially regarding settlements and catchment's areas identified by ceramic scatters discovered during fieldwalking campaigns.

The aim of the project is to delve deeper into the analysis of spatial interactions between social groups, the resources that they exploit to insure their subsistence, the manufacturing of objects and trade with other groups over a short, medium and long distance, and the itineraries which govern the flow at different periods. The objective is clearly to model the evolution of spaces which concern these occupations and activities, and to characterize, on the one hand, their rate of human influence and, on the other hand, the evolution of this human influence: modes of the agrarian conquest of space, growth and decline, eventually the abandonment of the occupation, types of activities, the consumption of what kinds of objects. With this aim, the objective is to specify the quantitative and qualitative indicators which enable us to measure territorial dynamics, their attractiveness, their skills, their relation to the social hierarchy as a factor of stimulation and economic development. It will then be possible to establish a typology of territories according to their evolution and in function of different fixed indicators, and to compare the different paths taken by these territorial types [START_REF] Nuninger | From Archaeomedes to ArchaeDyn[END_REF][START_REF] Gandini | Hierarchical Typology and Settlement Patterns Modelling at Interregional Scale[END_REF] The parallel works done on settlements networks [START_REF] Gandini | Hierarchical Typology and Settlement Patterns Modelling at Interregional Scale[END_REF][START_REF] Bertoncello | The spatio-temporal dynamics of settlement patterns from 800 BC to 800 AD in Central and Southern Gaul: models for an interregional comparison over the long term[END_REF]) and cultivated areas [START_REF] Poirier | Measurements of Diachronic Stability of Agrarian Exploitation[END_REF][START_REF] Poirier | Measurements of Diachronic Stability of Agrarian Exploitation[END_REF]) led to a crossed approach, allowing to measure the extent of spatial dynamics that affected the exploitation centres on one hand, and cultivated areas on the other, from the Iron Age to Middle Ages .

Methodological principles

In order to analyze the spatial pattern and the degree of dispersion of the archaeological features over the long term, with a common and simple index, we used the method of the mean centre and standard deviation ellipse.

This technical procedure developed one century ago belongs to the field of descriptive spatial statistics ( [START_REF] Lefever | Measuring geographic concentration by means of the standard deviational ellipse[END_REF]. The principle of the mean centre (or barycentre, or centre of gravity) consists in summarising a seed point by calculating the geographical mean of the spatial entities studied, which can possibly be weighted. The centre of gravity is then "a useful reference for comparing several seed points in a same geographical area or even for comparing the position of the same seed points over time" [START_REF] Pumain | Les interactions spatiales[END_REF]. This method has already been used to this purpose in archaeology [START_REF] Hodder | Spatial analysis in archaeology[END_REF]Gauthier 2004;[START_REF] Gauthier | Recherche d'une méthode d'analyse comparée des flux de métal à l'Âge du Bronze (XVII e -IX e s.) en France orientale et en Transdanubie. Utilisation et adaptation de la méthode des sommes focales[END_REF][START_REF] Mordant | Dynamisme et espaces culturels. De la notion de mobilité au sein des populations, du Bronze moyen à l'étape initiale du Bronze final en France orientale (XV-XIIIe s[END_REF]). In our case, the innovative principle is to compare the mobility vs. stability of the centre of human activity based on several types of remains from these human activities.

After having summarized the available information on a single point, the second aspect of this descriptive spatial analysis consists in measuring the dispersion of the values around this index of central tendency, as the standard deviation shows the dispersion of values around a mean. The standard deviation ellipse applied in geography from the 1970's, is particularly helpful for comparative purposes [START_REF] Raine | Summarizing point patterns with the standard deviational ellipse[END_REF]. It was used in our project to measure the dispersion of the values around the mean centre of each phase. The standard ellipse takes into account the anisotropy of the distribution, in other words, its deformation in one direction or another [START_REF] Zaninetti | Statistiques spatiales, méthodes et applications géomatiques[END_REF]. The long axis of the ellipse indicates the direction of the greatest variability between the values. Naturally, the direction of this ellipse can be influenced by the very pattern of the study area.

In both thematics -settlements and agrarian areas -, the variations in the size of the standard deviation ellipses, from one phase to another, were the object of a simple rate of variation (t = ((Va -Vd) /Vd) * 100 where Va is the end value (or size) and Vd is the start value (or size)).

This standardization of measurements makes it possible to visualize on a single graph the variations of two chosen indices within different study areas [START_REF] Saligny | 7 millenia of territorial dynamics : settlement pattern, production and trades[END_REF].

Settlement patterns

The aim of the workgroup dealing with "Settlement patterns" is to analyze the intensity, the stability and the modalities of rural settlement over the long term (from 800 BC to 800 AD) and at an interregional scale, since 12 study areas are concerned [START_REF] Gandini | Hierarchical Typology and Settlement Patterns Modelling at Interregional Scale[END_REF][START_REF] Bertoncello | The spatio-temporal dynamics of settlement patterns from 800 BC to 800 AD in Central and Southern Gaul: models for an interregional comparison over the long term[END_REF].

Settlement is considered here in a broad sense, such as it is defined in geography by R. Brunet, including as well dwellings as simply agricultural constructions. (Brunet et al., 1992: 229).

Comparison of several spatial distributions, under different historical and geographical contexts, is one of the difficulties faced by archaeological studies of settlement. By using these two indices, we aim at identifying the phases of extension or withdrawal and pinpointing the variations in the localization of the occupied areas.

The method was tested on five study areas, located in central and southern France, and the calculations concerned the settlements occupied between the 5 th c BC and the 7 th c AD, that correspond at the most documented period in all the study areas. So, for each chronological period, the mean centre corresponds to the mean of the X and Y coordinates of the settlements, weighted by the hierarchical level of the settlements, such as it had been defined by an automatic classification [START_REF] Bertoncello | The spatio-temporal dynamics of settlement patterns from 800 BC to 800 AD in Central and Southern Gaul: models for an interregional comparison over the long term[END_REF]. This weighting was determined empirically after several tests. Thus, it doesn't take into account only the quantitative change in the settlement but also its structural transformation. Effectively, the variations of the number of occupied sites do not have the same significance depending on the nature and the function of the settlements: small uninhabited agricultural buildings have different demographic and do not play the same role as the villas or farms.

The position of each mean centre has not interest in itself. That is interpretable and carries historical information is the displacement of the location of the barycentre over the long term. The period by period comparison of the displacements of the barycentre gives a first synthetic indicator of the distribution of habitat in a given space.

In order to compare the spatial dynamics in the different micro-regions, it was necessary to normalize the measurements of the mean centre's displacements. For this, we looked at the rate of variations in the mean centre localization per period, which can then be expressed in percentages and visualize on a graph. As mean centre localization for each period can be influenced by the size and shape of the study area, the measurement was related to the greatest possible distance within the study area. The greatest possible distance is easily obtained by drawing a rectangle on the dataset and by measuring the diagonal. It is the sum (in meters) of displacements of the mean centre, between each century, related to the largest possible distance (in meters as well). This approach is simple to implement and enables quickly to compare the trajectory of the study areas between them.

As shown on figure 1, the settlement evolution varies according to different rhythms and modalities from one region to another. The study area comparison shows that the spatial evolution of settlement in the Languedoc, between 5thC BC and 7 th C AD, is the most stable, while the Berry is the most unstable region. The comparison of the 2 profiles of Eastern Languedoc (located in Transalpine Gaul) and Limagne (located in the Three Gauls) globally shows the same tendency, but with a more unstable evolution in the Limagne area, with important changes of localization from one period to another. On the contrary, in Languedoc, after a phase of mobility of the mean centres between the 5th and the 1st c BC, the localization of the settlement tends to stabilize during the roman period. The standard deviation ellipse expresses the dispersion of the settlements around the mean centre of each period. The more the ellipse is big and distended, the more the spatial and quantitative variability of the settlement pattern is important; on the contrary, a small ellipse contracted around the mean centre evidences the settlements concentration. The variations in the size of the standard deviation ellipse clarify the results of the mean centre analysis for Limagne and Languedoc:

-In Languedoc, between the 5th and 1st c BC, the ellipse alternates between phases of contraction and expansion, before stabilizing during the roman period, -In Limagne, the ellipse is very narrow and stretched between the 5th and the 3rd c BC. It becomes wider in the 2nd c BC, expressing the phase of settlements' creations in the 2nd c BC. This level is maintained until the end of the period, even if the ellipse tends to contract at the end of Antiquity, in relation to a drop in the number of occupied settlements.

These results obviously must be interpreted by combining them with the results of the other analyses, as the evolution curves. But this example already shows the utility of these two basic indices of spatial statistics for apprehending the spatio-temporal dynamics of settlement patterns and for making synthetic regional comparisons, whose aim is to identify common tendencies to local particularities which will have to be explained.

Agrarian manurings

In the context of the workgroup dedicated to "Catchment areas", the calculation of the mean centres and the standard deviation ellipses was used with the goal of measuring the extent of the spatio-temporal mobility which affects the localization of the cultivated areas over time and to model the displacements and the dispersion of agrarian areas at each time step. The location of infield areas is assessed by the finding of off-site material during fieldwalking campaigns. This off-site material is interpreted as remain of manure. Manure collected in stables is supposed to have been mixed with domestic refuses before its spreading in cultivated fields. Interpreting off-site scatters as remains of manure practices is an idea now supported by most researchers interested in the history of landscape and settlements [START_REF] Wilkinson | The definition of ancient manured zones by means of extensive sherd-sampling techniques[END_REF][START_REF] Bintliff | Off-site potery distributions: A regional and interregional perspective[END_REF]. This interpretation is based on several ancient texts which offer proof of this practice from Antiquity to the present [START_REF] Oschinsky | Walter of Henley and other treaties on estate management and accounting[END_REF]. It is also based on recent examples of excavated structures used to prepare manure [START_REF] Puig | Du fumier à l'épandage, enrichir le sol en Rousillon entre les XIe et XIVe siècles ap[END_REF] or manured zones identified by geochemical measurements, such as phosphates [START_REF] Neil Rimmington | Soil geochemistry and artefact scatters in Beotia, Greece[END_REF]. Since the 2000s, many projects have been based on the study of timespace distribution of off-site material in order to approach past landscapes and territories [START_REF] Nuninger | Peuplement et territoires protohistoriques du 8e au 1er s. av. J.-C. en Languedoc oriental[END_REF][START_REF] Nuninger | From Archaeological Sherds to Qualitative Information for Settlement Pattern Studies[END_REF][START_REF] Jones | Signatures in the soil: The use of pottery in manure scatters in the identification of medieval arable farming regimes[END_REF]). Analysing scatters of artefacts over the long term enables us to study the dynamics of arable spaces with precision, in terms of spatial variability and the quantitative investment of ancient societies within these spaces.

The application of mean centres to agrarian manuring was developed in the framework of a PhD thesis [START_REF] Poirier | Un espace rural en Berry dans la longue durée : expérience de micro-analyse des dynamiques spatio-temporelles du paysage et du peuplement dans la région de Sancergues (Cher)[END_REF][START_REF] Poirier | Un espace rural à la loupe : paysage, peuplement et territoires en Berry de la Préhistoire à nos jours[END_REF]). The calculated mean centre is the centre of gravity of a point cloud formed by all of the centroids of the collection units of off-site artefacts at different periods. The localization of the seed point remains, of course, unchanged from one phase to another; what changes is the weighting of each point according to the density of off-sites artefacts collected in each area for each phase under consideration. The mean centre is thus drawn towards the collection units which supplied the highest densities of off-site artefacts, those that are interpreted as having been the most intensely exploited during the chronological phase under study.

A different mean centre corresponds to each time step (or chronological phase); each mean centre is associated with a standard deviation ellipse whose size and orientation are also different (Figure 4). The measurement of the surface taken by each ellipse, and the observation of the variations in size of this ellipse from one phase to another can be interpreted as an index of the extension or the contraction of the intensely exploited areas. In practical terms, the more the ellipse is large and dilated, the more the spatial and quantitative variability of the manurings is important; on the contrary, an ellipse of a small size and restricted around the mean centre is a proof of the concentration of cultivated areas. To compare the dynamics between the two study areas, we also had to normalize the raw values of distance between the different mean centres. We choose to normalize it by the maximum possible distance within the extent of the study-area (i.e. the diagonal of the extent rectangle).

The variations of mean centres of the agrarian manurings within every micro-region draw rather different curves for each area (Figure 6). While the region of Sancergues presents a rather chaotic evolution, resulting from an alternation between phases of big mobility on spaces recently exploited and episodes of stabilization, the mobility curve of the agrarian spaces at Neuvy-le-Roi is much more regular. At Neuvy-le-Roi, the phase of the greatest mobility concerns the transition between the phase [5-1] and the phase [1-3]. The mean centre of manurings records then a movement of about 8 % of the biggest length of the study-area. This mobility decreases then regularly till the end of the Middle Ages (phase [11][12][13][14][15]). The transition between the phases [8-11] and [11][12][13][14][15] presents a mobility only of 2 %, equal besides in that recorded at Sancergues. The transition between late Middle Ages and modern period presents an increase of agrarian spaces' mobility, the mean centre moving of about 4 %. This mobility is however much lower than that recorded in the region of Sancergues where the value of this movement exceeds the 10 %.

Crossed approach: comparing the spatial dynamics of settlements and agrarian areas.

Using the same indices of spatial statistics, the two workgroups opened up the way to cross the observations realized. It is to compare the dynamics of the settlements with that of their associated agrarian areas. In other words, are the rhythms of evolution of settlements patterns the same as those of cultivated areas?

But this ambition is confronted by a certain number of difficulties in particular related to the different resolutions of the processes observed:  First, spatial resolutions which do not always make it possible to compare evolutions readable on a local scale and others on a regional scale. That's why the method is tested only on two study areas: Berry_Sancergues and Touraine_Neuvy-le-Roi, where settlements and agrarian areas were equally documented during fieldwalking campaigns and spatially overlap.  Second, temporal resolutions, determined in a large part by our capabilities to date archaeological objects. So, we will use the broad chronological phases of 3 or 4 centuries that are used to date agrarian manurings. In Berry_Sancergues (Figure 7), the transition between Iron Age and Roman period shows a low mobility of agrarian manurings (4% of greatest distance), while the settlement network shows a bigger mobility (9%).

The transition between the roman times and the late Antiquity and very early middle ages shows a low mobility of the settlement network, because a lot of the occupied sites were created during the previous chronological phase. On the contrary, the agrarian manurings shows a great mobility of more than 10% of the greatest possible distance. This can be explained by the fact that the manurings were very centralized around occupied sites in the previous phase, on the contrary to the Late Antiquity, where manurings are more dispersed all around the study area, and not only in the immediate vicinity of the settlements.

The transition between the Late Antiquity and the Early Middle Ages shows a big mobility, both for settlement networks and agrarian manurings. Most of the occupied sites of this phase are new ones, as for the majority of the cultivated spaces. This transition is the one of the enhancement of new spaces of the study area. Low values are observed for both settlements and manurings mobility during the transition between Early Middle Ages and Classical Middle Ages, essentially because we see a kind of stabilization of the settlements at the places that are still occupied today, while the agrarian space is organized around these settlements.

The greatest values of mobility are observed for the last transition, between Middle Ages and Modern Era, with a displacement of the mean centres of more than 10 and 20 % of the greatest possible distance, both for settlements and manurings. We know that this great mobility is related to the settling of new farms, exploiting completely new spaces in the Far East of the study area (Poirier 2010: 63). The Touraine_NLR study area shows a general trend to the stabilization of both settlement pattern and agrarian spaces, characterized by a regular decrease of the mobility value of the mean centres (Figure 8). This can be explained by the observation of an early polarization of settlement and agrarian activities around the settlements that are still occupied today, proven by the medieval written sources [START_REF] Poirot | Neuvy-le-Roi (Indre-et-Loire), depuis ses origines jusqu'au XIXe siècle[END_REF]. The great moving of the settlement's mean centre between phase 1 and 2 can be explained by a data bias related to the fact that only a single settlement was discovered for the Iron Age while proofs of agrarian activities for this phase were identified in the whole study area. The observation of the distance separating the settlements' mean centre and the manurings' mean centre illustrates the disconnection between the archaeologically observed settlements pattern and spaces actually amended (Figure 9). This disconnection is high enough for phase 1 (Protohistory), decreases for phases 2 and 3 (Antiquity, Late Antiquity, early Middle Ages), and then increases steadily until the modern period (phase 6) where it is very strong in Berry_Sancergues (nearly 25% of the maximum possible distance within the study area).

Discussion

The steady increase of this disconnection is likely to be explained largely by the effect of sources. Indeed, the high values observed for phase 1 illustrate the difficulties in detecting protohistoric sites through their surface image which is often tenuous. Similarly, it should be noted that from the Middle Ages, even the early Middle Ages, the only archaeological image of the settlement pattern is certainly insufficient to accurately reflect all the points of settlement that had really existed. It should be possible to include in the analysis all the settlements that are mentioned in written sources and those that are currently still occupied, and implanted in these periods. It is not possible, because we can not precisely date the creation of these habitats, their first mention in written sources only providing a terminus ante quem.

Finally, spatial statistic indicators appear as interesting tools to develop comparative approaches: It allows a comparison between the dynamics of different kinds of archaeological entities (settlements vs agrarian areas), and comparison of the dynamics over different microregions. The normalization of the raw values of mean centre's displacements permit to put all the areas at the same level, and then confront the rhythms of their development, and enhancement of new spaces. These spatial statistic indicators allow going further in spatial analysis, especially by evaluating biases and sources effects and re-introducing the spatial dimension in spatial analysis.

To go further, our aim will be to compute these spatial statistic indicators in each area under study, both for settlements and cultivated areas, in order to reconstruct the general dynamics of land use in the long term. Those rhythms will then be compared to determine global tendencies and local particularities, considering historical, geographical and natural contexts.

Figure 1 :

 1 Figure 1: Indicator of general tendency of the stability / instability of settlement in 5 study-areas, between the 5thC BC and 7 th C AD

Figure 2 :

 2 Figure 2: Displacements of the settlement mean centre between the 5th c BC and the 7th c AD (data: Languedoc: F. Favory, E. Fovet, L. Nuninger, Cl. Raynaud ; Limagne: D. Coin, B. Dousteyssier, F. Trément, M. Ségard ; calculation: C. Gandini)

Figure 3 :

 3 Figure 3: Changes in the standard deviation ellipse in the Languedoc and Limagne study areas between the 5th C BC and the 7th C AD. (data: Languedoc: F. Favory, E. Fovet, L. Nuninger, Cl. Raynaud ; Limagne: D. Coin, B. Dousteyssier, F. Trément, M. Ségard ; calculation: C. Gandini)

Figure 4 :

 4 Figure 4: Study of agrarian manuring -Berry-Sancergues study area: Mean centre and dispersion (Sources and map: N. Poirier)Two study-areas were integrated in the calculation of these spatial indices: The area around Sancergues (in Berry) and the area around Neuvy-le-Roi (in Touraine) (Figure5).

Figure 5 :

 5 Figure 5: Location of the study-areas of Berry_Sancergues and Touraine_NLR (sources : Poirier 2007 ; Poirot 1998 ; map : N. Poirier)

Figure 6 :

 6 Figure 6: Dynamics of agrarian manurings in Sancergues (Berry, France) and Neuvy-le-Roi (Touraine, France)

Figure 7 :

 7 Figure 7: Mean centres' mobility of settlements and agrarian areas in Sancergues (Berry, France)

Figure 8 :

 8 Figure 8: Mean centres' mobility of settlements and agrarian areas in Neuvy-le-Roi (Touraine, France)

Figure 9 :

 9 Figure 9: Distance between mean centres of settlements and agrarian areas in Berry_Sancergues and Touraine_Neuvy-le-Roi.

Programme "Blanc"

2008, contrat ANR-08-BLAN-0157-01, 2009-2011.