Linfeng Ding
email: dingli@uni-mainz.de

Manuel Thieme

Sylvie Demouchy

Clemens Kunisch

Boris J P Kaus

Effect of pressure and temperature on viscosity of a borosilicate glass

Keywords: Borosilicate glass, high pressure, deformation experiment, viscosity, Maxwell model

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

Introduction

Hardness measurements [START_REF] Peter | Densification and flow phenomena of glass in indentation experiments[END_REF][START_REF] Shang | Creep behavior of soda-lime glass in the 100-500 K temperature range by indentation creep test[END_REF] have reported that glass flows viscously at room temperature as well as at high temperature (e.g., up to 630 °C for soda-lime-silica float glass [START_REF] Wilantewicz | Vickers indentation behavior of several commercial glasses at high temperatures[END_REF] ). Nevertheless, during industrial production processes, scribing with defined tools and breaking are used to separate glass sheets. If scribing velocity or scribing force is not well adjusted, it can lead to poor edge quality with non-desirable cracks at the edges or cracks are not deep enough (or sufficiently well developed) for further breaking. [START_REF] Peter | Densification and flow phenomena of glass in indentation experiments[END_REF][START_REF] Dick | New experiments on the microplasticity of glass[END_REF][START_REF] An | Experimental investigation on cutting characteristics in nanometric plunge-cutting of BK7 and fused silica glasses[END_REF] The propagation of these cracks, desired or not, depends on the stress state underneath the indenter during scribing. [START_REF] Hagan | The origin of median and lateral cracks around plastic indents in brittle materials[END_REF][START_REF] Lawn | Kinetics of shear-activated indentation crack initiation in soda-lime glass[END_REF] The actual distribution of stress components in the glass during scribing remains, to date, poorly quantified, and thus continues to be challenging to model numerically. [START_REF] Lee | The contact problem for viscoelastic bodies[END_REF][START_REF] Alkorta | Absence of one-to-one correspondence between elastoplastic properties and sharp-indentation load-penetration data[END_REF] To better quantify the viscous contribution to the deformation [START_REF] Rouxel | The brittle to ductile transition in a soda-lime-silica glass[END_REF] underneath an indenter or scribing tool, the effect of pressure on the viscosity of an industrial optical glass, needs to be quantified experimentally by in situ measurements.

The temperature dependence of the viscosity of glasses at ambient pressure has already been studied experimentally leading to numerous temperature-dependent models of viscosity for glasses. Notably, the Vogel-Tammann-Fulcher (VTF) model [START_REF] Fulcher | Analysis of recent measurements of the viscosity of glasses[END_REF] which is similar to the Williams-Landel-Ferry (WLF) model, [START_REF] Williams | The temperature dependence of relaxation mechanisms in amorphous polymers and other glass-forming liquids[END_REF] the free-volume model, [START_REF] Cohen | Molecular transport in liquids and glasses[END_REF][START_REF] Greet | Test of adam-gibbs liquid viscosity model with O-Terphenyl specific-heat data[END_REF][START_REF] Turnbull | On the free-volume model of the liquid-glass transition[END_REF] the Adam-Gibbs (AG) model [START_REF] Adam | On the temperature dependence of cooperative relaxation properties in glass-forming liquids[END_REF] and the Mauro-Yue-Ellison-Gupta-Allan(MYEGA) model. [START_REF] Mauro | Viscosity of glass-forming liquids[END_REF] The VTF model can accurately fit the temperature dependence of viscosity of SCHOTT N-BK7

 glass in the range from 1 to 13 log 10 (Pa•s).

Several experimental studies have quantified the pressure dependence of the viscosity of silicate melts, with a special focus on melts relevant to geological and volcanic systems, such as molten olivine, molten jadeite, dacite and basalt at very high temperatures (> 1700 °C) [START_REF] Mysen | The structure of silicate melts: Implications for chemical and physical properties of natural magma[END_REF][START_REF] Wang | Atomistic insight into viscosity and density of silicate melts under pressure[END_REF][START_REF] Kushiro | Changes in viscosity and structure of melt of NaAlSi 2 O 6 composition at high-pressures[END_REF][START_REF] Kushiro | Viscosity and structural-changes of albite (NaAlSi 3 O 8 ) melt at highpressures[END_REF] and pressure up to 8 GPa. Different behaviors were observed with increasing pressure as a function of the degree of polymerization. [START_REF] Wang | Atomistic insight into viscosity and density of silicate melts under pressure[END_REF][START_REF] Scarfe | Pressure dependence of the viscosity of silicate melts. Magmatic processes: Physicochemical principles[END_REF] The polymerization parameter is best described by the ratio of non-bridging oxygen (NBO) and tetrahedrally coordinated cations (T, e.g., Si 4+ , Al 4+ , Ti 4+ , Fe 3+ ). The ratio of NBO/T can be simply calculated by the ratio of O/T at ambient pressure with satisfying accuracy. [START_REF] Wang | Atomistic insight into viscosity and density of silicate melts under pressure[END_REF] However, under high pressure, highly polymerized silicate melts (NBO/T<1, e.g., obsidian, basalt, dacite) have a lower tetrahedral connectivity (less oxygen atoms in TO 4 will bond to adjacent tetrahedron and thus forming less 3D networks) than at ambient pressure. These highly polymerized silicate melts [START_REF] Wang | Atomistic insight into viscosity and density of silicate melts under pressure[END_REF] display a viscosity decrease with an increase of pressure. While the viscosity of depolymerized silicate melts (NBO/T≥2, molten peridotite, molten olivine, Al-poor melt), with higher tetrahedral connectivity, show a viscosity increase with increasing pressure. [START_REF] Wang | Atomistic insight into viscosity and density of silicate melts under pressure[END_REF] For borosilicate glasses, in addition to the NBO/T ratio, the structure is greatly influenced by the BO 3 /BO 4 ratio. [START_REF] Limbach | Plasticity, crack initiation and defect resistance in alkali-borosilicate glasses: From normal to anomalous behavior[END_REF] Indeed, the viscosity model presented by Jantzen [START_REF] Jantzen | Using polymerization, glass structure, and quasicrystalline theory to produce high level radioactive borosilicate glass remotely: A 20+ year legacy[END_REF] showed an irregular behavior of viscosity versus NBO in radio-active borosilicate glasses, making the prediction of the pressure effect on borosilicate problematic.

Previous experimental studies on the effect of pressure on viscosity and density of glasses have mostly reported a linear relationship for pressure dependence. For example, the pressure dependence of the viscosity of B 2 O 3 , studied by Sperry and Makenzie To date, the latest high-pressure dataset available for SCHOTT N-BK7  is based on relatively fast-cooled glass under high pressure (0.5 GPa to 1.5 GPa using a piston cylinder). Combined with the post-mortem density and volume recovery measurements, [START_REF] Ding | Pressure dependence of density and structural relaxation of glass near the glass transition region[END_REF] which also indicate a linear viscosity increase with increasing pressure (∆/∆P =0.05 log 10 (Pa•s)/100 MPa 28 ).

Several phenomenological models [START_REF] Tribone | Pressure-jump volume-relaxation studies of polystyrene in the glass-transition region[END_REF][START_REF] Hodge | Enthalpy relaxation and recovery in amorphous materials[END_REF][START_REF] Simon | Enthalpy recovery of a glass-forming liquid constrained in a nanoporous matrix: Negative pressure effects[END_REF][START_REF] Grassia | Modeling volume relaxation of amorphous polymers: Modification of the equation for the relaxation time in the kahr model[END_REF] describe the time dependence of glass relaxation as a function of both temperature and pressure history. A recent relaxation model was developed by Ding et al., [START_REF] Ding | Pressure dependence of density and structural relaxation of glass near the glass transition region[END_REF] following the approach of Gupta [START_REF] Gupta | Fictive pressure effects in structural relaxation[END_REF] to use two internal parameters, named fictive temperature 34 ( ) and fictive pressure 33 ( ). This type of phenomenological models can predict the pressure dependence of the viscosity of glass forming materials through the Maxwell relation:

(

where is viscosity, is the Maxwell relaxation time, and is the unrelaxed elastic shear modulus. [START_REF] Dingwell | Structural relaxation in silicate melts and non-newtonian melt rheology in geologic processes[END_REF] For instance, the two-internal-parameter relaxation model presented by Ding et al., [START_REF] Ding | Pressure dependence of density and structural relaxation of glass near the glass transition region[END_REF] predicts a positive linear dependence of pressure on viscosity for pressures up to 1.5 GPa based on the density and volume recovery measurements.

On the other hand, some studies directly express the effect of pressure on the viscosity model. For example, the model proposed by Avramov [START_REF] Avramov | Pressure dependence of viscosity of glassforming melts[END_REF] describes the pressure dependence of viscosity by relating the viscosity to the entropy of the glass forming melts, using a master equation . The equation derived by Gupta, [START_REF] Gupta | Negative-pressure dependence of viscosity[END_REF] which follows the Adam-Gibbs model, [START_REF] Adam | On the temperature dependence of cooperative relaxation properties in glass-forming liquids[END_REF] predicting a negative pressure dependence of viscosity.

SCHOTT N-BK7

 is available in large formats with high homogeneity. Internally at SCHOTT, it is used as a reference material. It was already used in several studies, [START_REF] Mosaddegh | Friction measurement in precision glass molding: An experimental study[END_REF][START_REF] Dyer | 157-nm laser micromachining of N-BK7 glass and replication for microcontact printing[END_REF][START_REF] Koontz | Prony series spectra of structural relaxation in N-BK7 for finite element modeling[END_REF] , 5 but the rheological properties were never quantified experimentally in situ at pressures above 1-atm. To quantify the viscosity of SCHOTT N-BK7  glass, compared to (28) we use a totally different experiment set up by performing in situ deformation experiments in a gas-medium pressure vessel.

Experimental method

Sample description and preparation

The borosilicate glass SCHOTT N-BK7  (chemical composition: 70.0 % SiO 2 , 11.5 % B 2 O 3 , 9.5 % Na 2 O, 7.5 % K 2 O and 1.5 % BaO) is an optical glass widely used for lenses and prisms, due to its high homogeneity and high light yield. The density of SCHOTT N-BK7  at room temperature and ambient pressure is 2.508  0.001 g/cm [START_REF] Wilantewicz | Vickers indentation behavior of several commercial glasses at high temperatures[END_REF] .

The glass transition temperature T g = 561 °C is defined as the temperature at which the viscosity is 12 log 10 (Pa•s). Cylinders of SCHOTT N-BK7  with a diameter of 9.95  0.05 mm and a length of 17.00  0.10 mm were drilled from large glass blocks, which had been slowly cooled at a rate of -0.4 °C/h. Both ends of the glass cylinder were mirror polished, and parallelism of both extremities was ensured.

Deformation experiments

Uniaxial deformation experiments were performed with a high temperature and high pressure in situ deformation gas-medium apparatus also called a Paterson Press. [START_REF] Paterson | Rock deformation experimentation[END_REF] The apparatus is used with a routine confining pressure of 300 ± 1 MPa, but the confining pressure can range from 50 to 350 MPa. Using Argon as pressure-medium prevents chemical reactions with the samples. Diffusion experiments at high pressure and high temperature have reported an Argon diffusivity of about -3.3 × 10 -20 cm 2 s -1 at 600 °C in silicate glasses. [START_REF] Carroll | Argon solubility and diffusion in silica glass -implications for the solution behavior of molecular gases[END_REF] Thus, for the short duration of our experiments (1h), the incorporation of Argon is considered negligible. The force is measured using an internal load cell (2 kN to 75 kN), resulting in an uncertainty in net axial stress (σ load ) of less than 20 MPa. In this study, uniaxial compressive deformation experiments on SCHOTT

N-BK7

 were performed at constant displacement rates (~1 × 10 -5 to ~5 × 10 -5 mm/s) resulting in almost constant strain rates.

For the deformation experiment, the glass samples were encapsulated with zirconia and alumina pistons and alumina spacers inside a copper jacket as shown in Fig. 1, following well-established experimental protocols. [START_REF] Karato | Rheology of synthetic olivine aggregatesinfluence of grain-size and water[END_REF][START_REF] Mei | Influence of water on plastic deformation of olivine aggregates 1. Diffusion creep regime[END_REF][START_REF] Mei | Influence of water on plastic deformation of olivine aggregates 2. Dislocation creep regime[END_REF][START_REF] Demouchy | Experimental deformation of olivine single crystals at lithospheric temperatures[END_REF][START_REF] Demouchy | Low strength of earth's uppermost mantle inferred from tri-axial deformation experiments on dry olivine crystals[END_REF] A self-made internal furnace is used to anneal the sample at high temperature, with a low thermal gradient (±5 °C) along the sample zone. The thermal gradient is established from separate calibration experiments before deformation experiments. The temperature was increased at a rate of ~15 °C/min to the target temperature and maintained stable for 20 to 30 minutes to insure sufficient relaxation of the glass structure before applying an axial load to the sample. The temperature was kept constant using a proportional-integral-derivative (PID) controller during the mechanical experiment as a small variation in temperature could result in significant variations of viscosity. A data point is recorded each 1 or 2 s. After the end of the deformation, the piston was maintained in position, and the temperature was decreased at a rate of 70 °C/min by decreasing power. Afterward, the pressure was slowly released at a rate of ~10 MPa/min to ambient pressure.

Data treatment

The mechanical data must be corrected to account for several experimental approximations adequately. Firstly, the displacement has to be corrected for the elastic distortion [START_REF] Paterson | A high-pressure, high-temperature apparatus for rock deformation[END_REF] of the apparatus (stiffness = 82.5 kN/mm). Secondly, to determine the contribution of the copper jacket to the measured force, the power flow law for copper from [START_REF] Frost | Deformation mechanism maps: The plasticity and creep of metals and ceramics[END_REF] is used. Thirdly, values of stresses and strain rates were corrected for instantaneous changes in the cross-sectional area of the copper jacket and the glass cylinder during deformation assuming that samples maintained a cylindrical shape (i.e., negligible barreling) and constant volume. These corrections of the mechanical data are identical to corrections used in numerous published experimental studies [START_REF] Karato | Rheology of synthetic olivine aggregatesinfluence of grain-size and water[END_REF][START_REF] Mei | Influence of water on plastic deformation of olivine aggregates 1. Diffusion creep regime[END_REF][START_REF] Mei | Influence of water on plastic deformation of olivine aggregates 2. Dislocation creep regime[END_REF][START_REF] Demouchy | Experimental deformation of olivine single crystals at lithospheric temperatures[END_REF][START_REF] Demouchy | Low strength of earth's uppermost mantle inferred from tri-axial deformation experiments on dry olivine crystals[END_REF][START_REF] Demouchy | Viscoplasticity of polycrystalline olivine experimentally deformed at high pressure and 900 degrees C[END_REF] on polycrystalline silicates using the same gas-medium deformation apparatus, guaranteeing that this source of uncertainty was properly corrected.

The mean stress (total pressure) on the glass sample during deformation is defined as:

σ mean = (σ load + 2σ Ar ) / 3 (2)
where σ load is the net uniaxial stress [START_REF] Paterson | Experimental rock deformation-the brittle field[END_REF] measured by an internal load cell and σ Ar is the confining pressure of Argon.

Results

Experimental conditions and measured mechanical data obtained from the in situ deformation experiments are reported in Table 1. Repeating the experiments permits to check the reproducibility of the mechanical results and the number of the repeated experiments at identical experimental conditions are reported in Table 1 and shown in

Mechanical data

The stress-strain curves for all deformation experiments are shown in Fig. 2. The deformation experiments at constant displacement rates yield strain rates ranging from 9.70 × 10 -6 s -1 to 4.98 × 10 -5 s -1 and finite strains ranging from 1.4 % to 8.9 %. To find the proper temperature and strain rate to study the pressure dependence of the viscosity of our samples, we ran experiments at different temperatures. The temperature dependence is shown in Fig. 2a for experiments performed at a confining pressure of 300 MPa, at a constant displacement rate of ~5 × 10 -5 s -1 and for temperatures from

Fig. 2 Net uniaxial stress versus strain (a) at a constant confining pressure and constant strain rate for increasing temperatures; (b) at a constant temperature and constant strain rate for increasing confining pressures; (c) at a constant temperature and constant confining pressure for increasing strain rates. For all curves, 1 data point out of every 50 points is shown for clarity.

To determine the pressure dependence, we conducted experiments at a constant temperature of 575 °C, at a constant strain rate of ~5 × 10 -5 s -1 and confining pressures ranging from 100 MPa to 300 MPa as shown in Fig. 2b. The σ load increases with increasing pressure from 340, to 430, and to 560 MPa for a confining pressure of 100, 200 and 300 MPa, respectively (see Table 1). The sample at 100 MPa confining pressure reached steady-state deformation after approximately 4% of strain, which is significantly earlier than for the experiments at 200 MPa and 300 MPa. The stress of the experiment at 100 MPa slightly decreases after reaching steady state, indicating some limited weakening.

To quantify the effect of strain rate, experiments at a constant temperature of 575 °C, a constant pressure of 300 MPa and axial deformation at different strain rates ranging from 9.7 × 10 -6 s -1 to 4.82 × 10 -5 s -1 , as shown in Fig. 2c. The maximal stress increases with increasing strain rate. Indeed, the σ load increases by ~91 MPa from ~1 × 10 -5 s -1 to ~3 × 10 -5 s -1 , while the σ load increases by ~325 MPa from ~3 × 10 -5 s -1 to ~5 × 10 -5 s -1 (Fig. 2c).

Maxwell model

We use both a single Maxwell and a two-Maxwell viscoelastic model to fit the experimental data combined with an optimization algorithm (MATLAB scripts are provided as supplementary materials).

The single Maxwell viscoelastic model (1-Maxwell) is defined by [START_REF] Kaus | Effects of elasticity on the rayleigh-taylor instability: Implications for large-scale geodynamics[END_REF] (3)

where is the total strain rate, is the strain rate of elastic deformation, is the strain rate of viscous deformation, is the shear stress converted from the net uniaxial stress by assuming a constant Poisson's ratio of 0.206 for SCHOTT N-BK7  , is the viscosity and G is the elastic shear modulus. To keep the model simple, we neglect the contributions from the volumetric bulk deformation strain. Under a constantly applied strain rate , the stress can be solved by integrating over time as:

(

where is the Maxwell relaxation time can be calculated by Eq. ( 1), and is the initial stress (in the model, we assume ). For the elastic shear modulus, we assume G as a free parameter, determined during the inversion process. In the limited temperature range of the present study, the dependence of viscosity on temperature can be described by an Arrhenius law. To account for the influence of pressure on viscosity, the Arrhenius equation is extended to Eq. ( 6) by introducing an activation volume as has been done in several studies. [START_REF] Ding | Pressure dependence of density and structural relaxation of glass near the glass transition region[END_REF][START_REF] Jin | Pressure effect on the structural relaxation and glass transition in metallic glasses[END_REF]

Combining Eq. ( 6) with the idea from Narayanaswamy [START_REF] Narayanaswamy | A model of structural relaxation in glass[END_REF] , we assume that the viscosity depends on temperature and pressure as (7) where = 10 5 Pa and =561 °C are reference material parameters, and the initial viscosity , activation energy and activation volume are free parameters [START_REF] Ding | Pressure dependence of density and structural relaxation of glass near the glass transition region[END_REF] .

The 1-Maxwell model can fit each experiment separately very well by allowing a change of elastic shear modulus G between the different experiments. However, the elastic shear modulus is known to be nearly constant in both glass science and geoscience applications [START_REF] Narayanaswamy | A model of structural relaxation in glass[END_REF][START_REF] Scherer | Theories of relaxation[END_REF][START_REF] Webb | Non-newtonian rheology of igneous melts at high stresses and strain rates -experimental results for rhyolite, andesite, basalt, and nephelinite[END_REF][START_REF] Zheng | Viscosity of glass-forming systems[END_REF] for small pressure changes (<2 GPa). Therefore, we tested fitting a single set of material parameters to all experimental mechanical data using a 1-Maxwell model. In order to fit the model to the experiments, we employ a multidimensional unconstrained nonlinear minimization optimization algorithm (using the Nelder-Mead Simplex search algorithm as incorporated in the MATLAB function fminsearch [START_REF] Lagarias | Convergence properties of the nelder-mead simplex method in low dimensions[END_REF] ), to find the best-fit parameters for the model. The best-fit model we obtained (after experimenting with a range of starting values , GPa, ∆H/R = 7.4×10 4 K, = 8.129×10 -6 m 3 /mol [START_REF] Jain | Finite element modelling of stress relaxation in glass lens moulding using measured, temperature-dependent elastic modulus and viscosity data of glass[END_REF] ) is shown in Fig. 3a (best-fit yields G=24.9 GPa ∆H/R= 5.4×10 4 K which are different than the reference value of G ref =34

GPa [START_REF] Jain | Finite element modelling of stress relaxation in glass lens moulding using measured, temperature-dependent elastic modulus and viscosity data of glass[END_REF] and (∆H/R) ref = 7.4×10 4 K at 1-atm), which is rather poor. For many glass types, more than one Maxwell element is needed to describe the experimental data. [START_REF] Koontz | Prony series spectra of structural relaxation in N-BK7 for finite element modeling[END_REF] To improve the fitting, we applied a generalized Maxwell model with 2 Maxwell elements [START_REF] Jain | Finite element modelling of stress relaxation in glass lens moulding using measured, temperature-dependent elastic modulus and viscosity data of glass[END_REF] (2-Maxwell model applied, MATLAB script available in supplemental materials), which models the rheology of glass as two parallel Maxwell bodies, which was previously found to approximate the uniaxial stress relaxation of glass successfully. [START_REF] Scherer | Viscoelastic-elastic composites: I, general theory[END_REF] For the stress revolution, we introduce a weighting parameter w for 2-Maxwell model and a free parameter fac to separate the viscosities, 54 (8) where the stresses are (9) (10) and the viscosities are

(11) (12) (13) (14) (15) (16)
The best overall fit parameters and curves of the 2-Maxwell model are shown in Fig. 3b. We note that the best-fit is not unique, and it is slightly influenced by the initial guess of the free parameters. We have performed a series of tests with different initial parameters and found an initial guess of (best-fit of 1-Maxwell model), GPa, ∆H/R = 7.4×10 4 K, = 8.129×10 -6 m 3 /mol 59 , w=0.7 and fac = 10

to give an overall lowest misfit. The best-fit G=34.2 GPa is fit with the reference value of 34 GPa, [START_REF] Jain | Finite element modelling of stress relaxation in glass lens moulding using measured, temperature-dependent elastic modulus and viscosity data of glass[END_REF] while the best-fit activation energy ∆H/R= 6.1×10 4 K (482 kJ•mol -1 ) is lower than the value obtained at 1-atm of ∆H/R= 7.4×10 4 K (615 kJ•mol -1 ). [START_REF] Jain | Finite element modelling of stress relaxation in glass lens moulding using measured, temperature-dependent elastic modulus and viscosity data of glass[END_REF] The fitting is improved compared to the 1-Maxwell model. Still, the fit is not perfect, especially for the data measured at 1×10 -5 s -1 and 3×10 -5 s -1 . A generalized Maxwell model with 3-Maxwell elements was also tested, but the results show no significant fitting improvement compared to the 2-Maxwell model. Therefore, the outcome of a 3-Maxwell model is not further discussed here.

Viscosity

The glass sample is in a uniaxial stress state during the deformation experiment, and the viscosity is then calculated as follows: 61

where is the net uniaxial stress, and is the instantaneous strain rate. The uncertainty of the applied force results in a stress uncertainty of ± 6 MPa. Together with a 1 % error of the strain rate, this yields an uncertainty of the viscosity of ± 0.01 The effect of pressure on the viscosity of SCHOTT N-BK7  is rather small.

Viscosity-time curves for 100 MPa and 200 MPa of confining pressure, shown in Fig. 4b, almost overlap each other for finite strains < 1 %. Increasing the strain rate to ~ 5 × 10 -5 s -1 results in a viscosity offset of ± 0.2 log 10 (Pa•s) (Fig. 4c), which is close to the uncertainty of the machine and also consistent with the results from the 2-Maxwell model. Besides, samples deformed with different strain rates achieve steady state (Fig. 4c) at the similar strain of approximately 2 %, indicating that the viscosity of SCHOTT N-BK7  does not significantly depend on strain for this strain level.

Discussion

Since the viscosity increases with increasing pressure, SCHOTT N-BK7

 glass behaves as a depolymerized glass. The pressure yields a positive effect on forming the 3D networks (oxygen atoms in TO 4 connect to adjacent tetrahedron) in SCHOTT N-BK7  glass.

The present in situ experiments confirm the previous overall temperature dependence results fitted by an Arrhenius equation. However, the viscosity dependence on pressure is rather small and sensitive to uncertainties in temperature. Indeed, ± 5 °C can yield an uncertainty as high as ± 0.22 log 10 (Pa•s) on the viscosity at 300 MPa.

Glass (i.e., non-crystalline materials) typically behaves like a Newtonian fluid [START_REF] Simmons | Non-newtonian viscous-flow in glass[END_REF] under applied stress, yet non-Newtonian viscous flow has also been observed in sodalime silicate glass [START_REF] Simmons | Non-newtonian viscous-flow in glass[END_REF] at strain rates around 10 -4 s -1 . Here, mechanical data as a function of strain rate indicate that the SCHOTT N-BK7  glass at 300 MPa and 575 °C might have a non-Newtonian behavior (Fig. 2c), which partly affects the overall fit of the 2-Maxwell model in Fig. 3b. During the industrial glass production, the process parameters are chosen to avoid fracture, that is, to be slow enough to remain perfectly Newtonian. Thus, the influence of a large strain rate range on SCHOTT N-BK7  viscosity will need further experiments to properly quantify. On the other hand, the nonlinear viscosity change with increasing strain rates might also be due to the uncertainty in temperature, such as a local heterogeneous temperature distribution (e.g., radial gradient or the heat distribution difference between the deformation assembly containing the glass cylinder and the furnace calibration assembly containing an alumina cylinder).

For the large absolute viscosity discrepancy between experimental data and the latest two-internal-parameter relaxation model (based on static high pressure and volume recovery measurement), we do not have an unambiguous explanation. It may be due to the following aspects: (i) the deformation apparatus does not permit experiments at 1-atm and high temperature and thus the current data set cannot be confirmed by experimental data at 1-atm and 575 °C. (ii) The temperature calibration was done by using a hollow alumina cylinder instead of the SCHOTT N-BK7  glass cylinder, which may have affected the temperature distribution. We have tried to calibrate the furnace using a specially designed hollow SCHOTT N-BK7  glass cylinder, but at high temperature and pressure the hole for the thermocouple (which is at room pressure) shrunk as a result of viscous flow under high confining pressure, preventing any successful temperature calibration with a glass specimen. The different thermal conductivity between the SCHOTT N-BK7  glass hollow cylinder and the alumina hollow cylinder could have caused a slight change in temperature distribution. (iii) the surface tension [START_REF] Parikh | Effect of atmosphere on surface tension of glass[END_REF] and friction [START_REF] Liu | Sliding friction of copper[END_REF] between the glass sample and the copper jacket were not taken into account during the stress determination; or (iv) some other fundamental process is missing in the two-internal parameter relaxation model.

Conclusions

Deformation experiments in uniaxial compression on SCHOTT N-BK7  glass were performed at temperatures near the glass transition temperature (561 °C) and at confining pressures of 100 to 300 MPa. The main results of this study are: [START_REF] Ding | Pressure dependence of density and structural relaxation of glass near the glass transition region[END_REF] The reason for this discrepancy has not been identified clearly.

(d)

The mechanical data can be overall fitted by a Maxwell viscoelastic model with two parallel elements.

Fig. 1

 1 Fig. 1 (a) Photograph of the high-pressure assembly for a deformation experiment; 1-Zirconia piston; 2-Alumina piston; 3-Alumina spacers; 4-SCHOTT N-BK7  glass cylinder; 5-Fitted copper jacket; 6-Bottom steel piston; 7-Top steel piston.

Figure S1 (

 S1 Figure S1 (supplementary material).A few samples showed a variation in diameter along the deformed cylinder post-deformation (i.e., classic barreling), which could indicate that the temperature distribution on the sample is not perfectly homogenous. Or more likely that, the viscous flow continued to occur post-data acquisition and during quenching, when the temperature distribution has become heterogeneous (i.e., the bottom cylinder being hotter than top).

  550 °C up to 595 °C. The glass sample showed an approximately constant stress increase at 550 °C, with the stress-strain curve displaying an almost pure elastic behavior. Between 565 and 575 °C, the samples behave viscoelastic. Steady state (viscous) flow is finally achieved at 595 °C. The stress at 595 °C (81 MPa, sample LD18) is rather low with an applied force of approximately 7 kN.

Fig. 3

 3 Fig. 3 (a) Best-fit from the 1-Maxwell model with the fitting parameters: =4.9×10 12 Pa•s, G=24.9 GPa, = 5.4×10 4 K, and =2.0×10 -6 m 3 /mol; (b) Best-fit from the 2-Maxwell model with the fitting parameters: =5.4×10 12 Pa•s, G=34.2 GPa, = 6.1×10 4 K, =2.9×10 -6 m 3 /mol, w=0.83, and fac=15.6. The uncertainty of the shear stress is calculated from the standard deviation of steady-state viscosity and strain rate.

log 10 (

 10 Pa•s). The uncertainty of the temperature distribution along the sample (maximal ± 5 °C) would yield more important uncertainties on viscosity, up to ± 0.22 log 10 (Pa•s), based on the Arrhenius equation at 1-atm. The latter values are reported in Table 1. Our results show that to change the viscosity of SCHOTT N-BK7  glass by 0.1 log 10 (Pa•s) has to change the temperature by 2.5 °C (at 300 MPa) or the pressure by 100 MPa. As shown in Fig. 4a, the viscosity decreases by ~0.4 log 10 (Pa•s)/10 °C from 575 °C to 595 °C. Furthermore, the viscosity curve at 595 °C furthermore takes less time to achieve steady state than deformation experiments at 575 °C.

Fig. 4 Fig. 5 .

 45 Fig. 4 Time dependence of viscosity (a) at a constant confining pressure and constantstrain rate for increasing temperatures; (b) at a constant temperature and constant strain rate for increasing confining pressures; (c) at a constant temperature and constant confining pressure for increasing strain rates. The time dependence of strain is also shown in (c). For all curves, 1 data point out of every 50 points is shown for clarity.

Fig. 5

 5 Fig. 5 Temperature dependence of the viscosity of SCHOTT N-BK7  at high pressure compared to an Arrhenius model at 1-atm, a 2-Maxwell model at 300 MPa and the latest two-internal-parameter relaxation model at 300 MPa (Ding et al., 2018).

Fig. 6

 6 Fig.6 Viscosity as a function of the mean stress of SCHOTT N-BK7  for the experimental data at 575 °C compared to the values calculated by using a 2-Maxwell model as well as the latest two-internal-parameter relaxation model(Ding et al., 2018).

1 log 10 (

 110 viscosity of SCHOTT N-BK7  glass increases linearly with increasing pressure at ~0.1 log 10 (Pa•s)/100 MPa in the range of confining pressure investigated. (b) The pressure dependence of viscosity measured in this study (~0.Pa•s)/100 MPa) is consistent with the latest two-internal-parameter relaxation model. (c) At the same temperature and the same pressure, the viscosity measured in this study is approximately +1.0 log 10 (Pa•s) higher than in Ding et al. (2018).

  [START_REF] Sperry | Pressure dependence of viscosity of B 2 O 3[END_REF] at 380-465 °C and from 0.1 to 30 MPa, demonstrated a linear increase of viscosity with increasing pressure (∆/∆P = 0.28 to 0.59 log 10 (Pa•s)/100 MPa[START_REF] Sperry | Pressure dependence of viscosity of B 2 O 3[END_REF] ). A special parallelplate viscometer was developed by Schulze et al.,[START_REF] Schulze | Determination of the influence of pressure and dissolved water on the viscosity of highly viscous melts: Application of a new parallelplate viscometer[END_REF] using an internally heated pressure vessel (IHPV) to deform the standard melt DDG1, Di 100 (Di=CaMgSi 2 O 6 ) and Ab 55 Di45

(Ab=NaAlSi 3 O 8 ) silicate melt. And their results also suggest a linear change with increasing pressure up to 350 MPa (∆/∆P = -0.12 to 0.23 log 10 (Pa•s)/100 MPa 26 ).

Deformation of float glass with high water content

[START_REF] Gaudio | Viscosity and glass transition temperature of hydrous float glass[END_REF]

likewise reported a linear viscosity increase with increasing pressure from 100 MPa to 400 MPa, and a negative correlation with increasing water content (0.03-4.87 wt.%), (∆/∆P < 0.22 log 10 (Pa•s)/100 MPa 27 ).

Table 4 - 1

 41 Experimental conditions for deformation experiments and mechanical results

	Samples	Length [mm]	Diameter [mm]	Confining Pressure [MPa]	Temperature [°C]	Displacement Rate [mm/s]	Strain Rate [s -1 ]	Finite Strain [%]	Axial Stress [MPa]	Mean Stress [MPa]	Viscosity* Log 10 [Pa•s]	Repeat times
				Constant pressure and strain rate, variable temperature
	LD6	17.10	9.95	300	550 ± 5	5 × 10 -5	1.50 × 10 -5	1.4 802 ± 6 467 ± 7	-
	LD9	17.07	9.95	303	565 ± 5	5 × 10 -5	3.89 × 10 -5	3.9 792 ± 6 466 ± 7	-
	LD15 17.05	9.95	302	575 ± 5	5 × 10 -5	4.82 × 10 -5	6.5 560 ± 6 388 ± 7 12.6 ± 0.2
	LD17 17.10	9.95	303	585 ± 5	5 × 10 -5	4.77 × 10 -5	8.7 219 ± 6 275 ± 7 12.2 ± 0.2
	LD18 17.10	9.95	300	595 ± 5	5 × 10 -5	4.98 × 10 -5	8.4	81 ± 6 227 ± 7 11.7 ± 0.2
				Constant temperature and strain rate, variable pressure
	LD13 17.05	9.95	105	575 ± 5	5 × 10 -5	4.67 × 10 -5	5.7 340 ± 6 183 ± 7 12.4 ± 0.2
	LD16 17.05	9.95	203	575 ± 5	5 × 10 -5	4.86 × 10 -5	6.5 430 ± 6 279 ± 7 12.5 ± 0.2
	LD15 17.05	9.95	302	575 ± 5	5 × 10 -5	4.82 × 10 -5	6.5 560 ± 6 388 ± 7 12.6 ± 0.2
				Constant temperature and pressure, variable strain rate
	LD20 17.00	9.95	301	575 ± 5	1 × 10 -5	9.70 × 10 -6	7.6 144 ± 6 249 ± 7 12.7 ± 0.2
	LD19 17.02	9.95	302	575 ± 5	3 × 10 -5	2.88 × 10 -5	8.9 235 ± 6 280 ± 7 12.4 ± 0.2
	LD15 17.05	9.95	302	575 ± 5	5 × 10 -5	4.82 × 10 -5	6.5 560 ± 6 388 ± 7 12.6 ± 0.2
	* Viscosity is calculated according to Eq. 4.1

Acknowledgments S.D. thanks D. Mainprice and M.A. Bouhifd for informative discussions during this study. The authors thank the technicians and engineers at the SCHOTT machine-shop for the sample preparation. This project has received funding from the European Union's Horizon 2020 research and innovation program under the Marie Sklodowska-Curie grant agreement N° 642029.

The model developed by Ding et al. [START_REF] Ding | Pressure dependence of density and structural relaxation of glass near the glass transition region[END_REF] can be transformed into a viscosity model (18) and gives a positive linear dependence of viscosity on the pressure of SCHOTT N-

BK7

 , with a slope of 0.05 log 10 (Pa•s)/100 MPa. The pressure dependent viscosity data from this study, together with the viscosity predicted by the 2-Maxwell model and the latest two-internal-parameter relaxation model [START_REF] Ding | Pressure dependence of density and structural relaxation of glass near the glass transition region[END_REF] are shown in Fig. 6. We recall that the mean stress used in Fig. 6 equals to σ mean = (σ load + 2σ Ar ) / 3 (see Eq. ( 2), [START_REF] Paterson | Experimental rock deformation-the brittle field[END_REF] ) and that the viscosity is calculated using Eq. ( 17


Ema Falomir

Serge Chaumette
email: serge.chaumette@labri.fr

Gilles Guerrini
email: gilles.guerrini@fr.thalesgroup.com

Mobility Strategies based on Virtual Forces for Swarms of Autonomous UAVs in Constrained Environments

Keywords: UAV, Swarm, Autonomous Swarm of UAVs, Compactness, Virtual Forces, Distributed Robotic Systems

The usage of autonomous unmanned aerial vehicles (UAVs) has recently become a major question. For wide area surveillance missions, a swarm of UAVs can be much more efficient than a single vehicle. In this case, several aircrafts cooperate in order to fulfill a mission while avoiding collisions between each other and with obstacles. This article proposes original distributed mobility strategies for autonomous swarms of UAVs, the goal of which is to fulfill a surveillance mission. Our work is based on virtual forces and our approach allows narrow areas crossing that require a compact formation of the autonomous swarm.

INTRODUCTION

In the past few years, the usage of unmanned aerial vehicles (UAVs) has been widely studied. In this paper, we consider rotor wings, which allow low speed maneuvering and hovering.

A set of several UAVs deployed on the same area, independent of any infrastructure and able to establish peer-to-peer connections form a flying adhoc network (FANET) [START_REF] Bekmezci | Flying Ad-Hoc Networks (FANETs): A survey[END_REF]. A FANET offers several advantages in comparison to a single UAV. Additionally to a larger coverage, if one aircraft of a FANET encounters a failure, the assigned mission can still be achieved. Furthermore, each UAV can embed a different sensor, which could not fit inside a single aircraft.

Several studies have already been carried out on surveillance missions lead by FANETs. The UAVs can all be multirotors [START_REF] Chaumette | CARUS, an operational retasking application for a swarm of autonomous UAVs: First return on experience[END_REF] but multirotors also can be associated with fixed-wings [START_REF] Jaimes | An approach to surveillance an area using swarm of fixed wing and quad-rotor unmanned aerial vehicles UAV(s)[END_REF], [START_REF] Bouvry | Using heterogeneous multilevel swarms of UAVs and high-level data fusion to support situation management in surveillance scenarios[END_REF].

In this paper, we will consider a FANET composed of multirotors with the same kinematics characteristics, able to fulfill a mission autonomously. We will refer to this particular kind of FANET as "autonomous swarm". In this configuration, an operator only deals with a single entity contrary to a multiplatform system.

To evolve in partially unknown environment autonomously, the UAVs have to calculate their flight plans online. Here, not to be dependent of a single UAV, each one calculates its own flight plan.

Methods based on artificial potential fields (introduced in [START_REF] Khatib | Real-Time Obstacle Avoidance for Manipulators and Mobile Robots[END_REF]) are often chosen for online and decentralized applications thanks to their limited needs in terms of computation and because they are distributed by nature. They can be used to maximize area coverage [START_REF] Howard | Mobile Sensor Network Deployment using Potential Fields: A Distributed, Scalable Solution to the Area Coverage Problem[END_REF], [START_REF] Gobel | A Model of Autonomous Motion in Ad Hoc Networks to Maximise Area Coverage[END_REF] or to obtain specific connectivity characteristics [START_REF] Casteigts | Biconnecting a Network of Mobile Robots using Virtual Angular Forces[END_REF]. Some studies also aim at reaching a target with robot cars thanks to artificial forces approaches [START_REF] Boonyarak | Real-Time Obstacle Avoidance for Car Robot Using Potential Field and Local Incremental Planning[END_REF], [START_REF] Jin | A Stable Switched-System Approach to Obstacle Avoidance for Mobile Robots in SE(2)[END_REF]) but these methods have not been widely extended to swarms.

In this paper, we consider a cooperative and autonomous swarm of multirotors in the context of a surveillance mission of two areas of interest (AoIs), as represented on the Figure 1. The two AoIs are separated by a narrow passage (a tunnel), possibly containing static obstacles of unknown position and size, requiring a compact formation of the swarm. We propose a new distributed mobility strategy for autonomous swarms of UAVs, based on virtual forces in order to avoid collisions between the UAVs and with obstacles, and maximize the communication links along a surveillance mission.

The remainder of this paper is organized as follows. Section 2 presents the scenario and the associated hypothesis. Section 3 introduces the forces used in the mobility strategies thereafter presented in section 4. Section 5 introduces the metrics used to evaluate the efficiency of our mobility strategies. The results are shown in section 6. Finally, section 7 is dedicated to the conclusion of this work and sketches future research directions.

SCENARIO AND HYPOTHESIS Scenario

Our mission aims at the successive surveillance, by an autonomous swarm of UAVs, of two AoIs separated by a narrow passage possibly containing obstacles unknown prior to the mission. The UAVs are multirotors equipped with limited range sensors supporting the surveillance mission and allowing obstacle detection and with communication systems.

Hypothesis

We suppose that there are no obstacles in the AoIs but only in the narrow passage which is the phase of the scenario that we want to study in details. We also assume consider that the UAVs know the coordinates of the AoIs corners and of the passage entrance. Additionally, it is supposed that no UAV is added or removed from the system.

As we are interested in mobility strategies, we do not consider neither "low-level" details in this first study (communication protocols, localization, etc.) nor the UAVs flight dynamics. Moreover, it is supposed that the UAVs have enough energy to fulfill the complete mission.

CONTRIBUTION 1: VIRTUAL FORCES

To achieve a surveillance mission and to avoid collisions, mobility strategies are proposed based on virtual forces and on distances of interest which are presented in the following subsections.

For this study, we have chosen to define the forces with a magnitude comprised between -1 and 1. Attractive forces are positive, repulsive forces are negative. This permits to easily compare the influence of several forces in a system which can be complex. Moreover, as this is a preliminary study, forces are defined on the basis of linear functions.

Safety Distance

In order to avoid collisions, a minimal safety distance, D has to be maintained between each UAV and other objects. For this preliminary study, we define this distance as follows:

1.5

(1)

Inter-UAV Wanted Distance

During the mission, two goals must be achieved. First, the maintenance of the communication links, and second the coverage of the AoI with little redundancy. Consequently, we define an inter-UAV wanted distance which takes into account the communication and the surveillance sensor ranges:

min 2 , (2)
where k ∈ 0; 1 .

It is important to note that D must be greater than (or equal to) D to use the strategy described in the following sections. In practice, this condition is not limiting as UAVs always embed sensors with range larger than their own size.

Attraction-repulsion between UAVs

To fulfill their mission, the UAVs should remain spaced of D . In this aim, an attractive force is applied when they are further and a repulsive one when they are closer, depending on their distance d . As the UAVs safety is more important than the mission accomplishment, the repulsive forces have to be stronger than the attractive ones. The attractive-repulsive virtual forces between UAVs are directed from the considered UAV to its neighbor and of magnitude F (see Figure 2) defined in (3). Note: These notations will be used all along the paper.

1 1 if 0.5 if 2 1 2 1
(3)

Attraction towards a Goal

To fulfill the mission, the UAVs have to deploy on the AoIs, fly over them, and cross the passage. To achieve these goals, each UAV calculates temporary targets all along the mission and is subject to an attractive force towards them.

We have chosen to give a constant magnitude to the attractive force towards a goal till a certain distance. In the area around the target, the force decreases linearly with respect to the robot-goal distance d , in order to avoid oscillations around the goal. This decrease near the target also favors the other forces in this area.

The attraction to a goal is directed towards the goal and is of magnitude F (represented on Figure 3) defined as follows: Note: This notation will be used all along this paper.

0.5 if 2 else (4)

Repulsion Due to Obstacles

The repulsion due to obstacles is based on the safety distance D . In this first study, we admitted that if a UAV is further away than three times its safety distance from an obstacle, it should not consider the obstacle. Nonetheless, in future work this repulsive force should be based on other characteristics. The repulsive force due to an obstacle is directed at the opposite of the obstacle. A first magnitude has been tested, F , defined as follows:

3 1 if 3 0 e l s e (5)
where d is the UAV-obstacle distance. Using this first approach, the repulsion due to obstacles was too high compared to the attraction-repulsion between UAVs. As a result, the swarm did not remain connected in the tunnel, and collisions between UAVs occurred. Indeed, between 2D and 3D , the obstacle had a too high influence. So we modified the magnitude to F , represented in the Figure 4, defined as follows:

0 i f 3 9 20 1 if 2 1 10 3 else (6)

Weighted Average

At each step, each UAV is subject to several virtual forces as illustrated on the Figure 5. UAVs calculate a weighted average of all the virtual forces applied to them and their next direction is resulting of this weighted average. The weight of each force is the exponential of its magnitude, as proposed for spring forces in a work on network biconnection [START_REF] Casteigts | Biconnecting a Network of Mobile Robots using Virtual Angular Forces[END_REF].

CONTRIBUTION 2: MOBILITY STRATEGIES

As introduced in section 2, our scenario is composed of three main phases: survey the AoIs, make-up of the compact formation and cross the tunnel. In the following subsections we detail the mobility strategies for each of these three steps.

Surveillance Mobility Strategy

During the surveillance phase, each UAV should have, at each moment, a vision of the whole AoI as recent as possible. To be able to build this view in a cooperative manner, the UAVs have to be deployed all over the area and maintain communication links. An adapted mobility strategy consists in creating an "S-shaped" travel over the AoI, as illustrated in Figure 6. The mesh to support the movements of the UAVs is calculated at the very beginning of the mission and is composed of cells of size D D . In this mobility strategy, the targets of the UAVs are located at the extremity of the lines. Once a UAV reaches its temporary goal, it calculates its next target, depending on its identifier and of the number of vehicles in the swarm. Its move towards this target starts when its 1-and 2-hop neighbors have reached their own goals. In this way, the connectivity of the swarm is maintained. During the surveillance, each UAV is subject to attractive-repulsive force with its neighbors and to attractive force towards its successive targets.

Compact Swarm Formation Setup

Once the first AoI has been covered, the swarm flies toward the passage in which it can encounter obstacles. In this study, the UAVs do not share the obstacles location, but in a further work this should be implemented. To facilitate this, a biconnectivity of the swarm is sought. In case of a very narrow passage (i.e. of a width smaller than 3D ), UAVs can only get through the tunnel in a queue (Figure 7 a), so biconnectivitv cannot be guaranteed. Else, the swarm can organize itself in a staggered rows formation (see Figure 7 b) andc)). The 2-row staggered formation can be setup when the passage width is greater than 2D √3/2D . To allow 3 rows, the tunnel width has to be greater than 2D √3 D . These values can easily be geometrically calculated. From the passage width, the number of vehicles in the swarm and its identifier, each UAV calculates a target close to the tunnel entrance, in a deterministic manner. The connectivity of the swarm during the surveillance ensures that all the UAVs know when they have to go towards the tunnel entrance. The positions of the UAVs in compact formation are separated of D , so there will be some communication links between the UAVs. So, only the repulsive forces between the UAVs and thee attractive force to their goal are necessary when the swarm travels towards the tunnel entrance. The swarm starts to get through the passage only when all the UAVs are at their position or when UAVs are at their correct location for a defined duration.

Passage Crossing Mobility Strategy

The swarm has to cross the tunnel to reach the second AoI.

When the UAVs get through the tunnel, they are subject to the virtual attractive-repulsive forces between each other, the repulsive force due to the obstacles and the attractive force towards a goal. As the UAVs do not know the shape of the passage, they calculate at each step a new temporary target in the direction of the second AoI. This target is located at a small distance (we arbitrarily chose D /4) and so produces a weak attraction, the goal being to favor collision avoidance and to maintain of the communication links.

If an obstacle is located between the UAV and its target, the UAV stops being attracted towards this point and it favors the attraction towards its neighbors so as not to stay in a local minima.

We consider that the UAVs can communicate even if an obstacle is located between them. Additionally, if a UAV shares a communication link with another one in the opposite direction of its target and if an obstacle is between them, the first UAV will not take into account the attractive force towards this neighbor (see Figure 8). A further simplification is made in this first study concerning the repulsion due to obstacles. Each UAV decomposes its environment in four equal areas: upper, below, backwards and ahead. They take into account a maximum of one repulsive force due to obstacles in each quadrant. If several obstacles are in a same quadrant, the closest one only is taken into account. This is enough to avoid obstacles. Nonetheless, limitations exist. Indeed, there is no definition of a new target if there is an obstacle between it and the UAV, and favoring the attraction towards other UAVs is not always sufficient as illustrated on Figure 9.

EVALUATION CRITERIA

The presents the Measures of Performance (MoPs) for the three parts of the scenario. For each MoP, an ideal result is given, not always achievable.

SIMULATION AND RESULTS

The following subsections present the simulation tool used in this study, the parameters of the simulations and the results for the different steps, according to the MOPs defined in section 5.

JBotSim Simulation Tool

To simulate the mobility strategies defined in the previous section, we chose the open source simulation library JBotSim [START_REF] Casteigts | JBotSim: a Tool for Fast Prototyping of Distributed Algorithms in Dynamic Networks[END_REF]. It is a tool for distributed algorithms fast prototyping in dynamic networks. Contrary to other well-known simulators such as NS3 [START_REF] Riley | The ns-3 Network Simulator[END_REF] or OMNet++ [START_REF] Varga | The OMNeT++ discrete event simulation system[END_REF], JBotSim does not implement real-world networking protocols, which was not required in this preliminary study.

Simulation Scenario and Parameters

We have conducted a series of simulation to validate our mobility strategies. For all the simulations the k factor used in the D definition is fixed at 0.75.

We chose to simulate UAVs of circular shapes of 20cm diameter, equipped with a communication system of 1.10m range and with a sensor used for both surveillance and obstacle detection of 50cm range.

The first AoI is 9m long and 7.50m large. At each simulation step, the UAVs travel at a random distance comprised between 5cm and 10cm in the resulting direction of the virtual forces average. Simulations have been performed with swarms composed of 1 to 6 UAVs and with 5 tunnels (see Figure 7, Figure 10 and Figure 11). For each parameter set we have run 30 simulations. At the entrance, shape a tunnel is 95cm large, tunnels b, d and e are 151cm large and finally tunnel c is 226cm large.

Simulation Results, Evaluation and Analysis

Overall Mission

During the whole mission, the priority is to safeguard the UAVs. To achieve this goal, the distance between a UAV and another object has to remain above the safety distance defined previously.

The minimal distance between the UAVs as well as between the UAVs and an obstacle, remained above 1.3 times D , for all the simulation ran. Thus, the main objective of the project is reached.

Surveillance

As expected, the simulations have shown that the more UAVs compose the swarm, the quicker the complete surveillance is performed (see Table 2).

The very small difference of duration between 4 and 5 UAVs is due to the number of rows of the mesh supporting the UAVs movements in the AoI. Indeed, a 4-UAV swarm has to cross the 12-row area three times to fully cover it, as a 5-UAV swarm.

Furthermore, a UAV moves towards its next target only when its neighbors have reached their own one. Indeed the more UAVs compose the swarm, the higher is the probability to have to wait for other UAVs. This is why the speed-up, defined as follows, is not linear:

SpeedUp n UAVs duration 1 UAV duration n UAVs (7)
Even if the objective of measuring a linear speed-up is not reached, we can observe a clear decrease in duration when the number of UAVs increases. Finally, as the mobility strategy for the surveillance phase is almost deterministic (random speed of the UAVs), the standard deviation is low and so an estimation of the surveillance duration can be done by running a configuration only once.

Furthermore, we studied the creation and loss of links during the surveillance step (see the values in Table 3). We can see that there are on average a few loss of links, and as many creation as loss. The objective for this MOP is thus reached.

Compact Formation Setup

The Figure 12 shows the duration between the end of the surveillance and the moment when the compact formation is set up, normalized by the necessary time to go from the further UAV position to the tunnel entrance. As UAV target is not exactly at the tunnel entrance, the normalization can be smaller than one. Width of tunnel shape b, d and e are the same, so we wait for an identical duration of compact formation setting up. Nevertheless, differences appear between the 3 curves, due to the zoom effect.

In the case of a single UAV, the duration is the smallest for all the tunnels, because there is no repulsive force with neighbors, and so the UAV is not subject to forces opposed to the target attraction.

For all the configurations, the normalized compact formation setup duration is smaller than 1.4, which is a good result. This means that even for a swarm composed of numerous UAVs, and independently of the compact formation shape, the formation setup is rather quick. We note that we tried to simulate a thinner width for Tunnel shape a (of width 2D instead of 3D ). In this case, only multi-UAVs swarms entered: for a single UAV, the repulsive forces at the entrance of the tunnel were greater than the attractive force towards the target. In the case of several UAVs, the attraction-repulsion between them was in the same direction as the attraction towards their target, which allowed the swarm to enter in the tunnel. Nevertheless it occurred that a 6-UAV swarm reached a local equilibrium at the passage entrance (see the progression of the swarm on Figure 13). In the case of such equilibrium, the closest UAVs to the tunnel were slowed down because of the repulsion due to the walls, while the other UAVs had a strong attraction towards their target and so did not remain behind the others.

Tunnel Crossing

The duration of tunnel crossing is the duration between the entrance of the first UAV and the exit of the last one. Results are represented in Figure 14.

The tunnel crossing duration for the thinner tunnel does not vary in function of the number of UAVs. This is because the UAVs have to be queued in this case and so the attractive-repulsive forces between them are in the direction of their target. As expected, tunnels of shape b and c show the same trend even if the largest one is crossed quicker.

We note that in the case of several UAVs and tunnel of shape b to e, the crossing duration is longer than for a single UAV. This is due to the attractiverepulsive forces between the UAVs which are not in the direction of their target.

For the tunnel shape d, the crossing duration follows the trend of tunnel shapes b and c until 4 UAVs and skyrockets for 5 UAVs. This is due to the creation of additional links between the UAVs during the travel through the linear parts which have to be broken to pass over the corner.

In the case of the tunnel shape e, we can notice a strong increase of the crossing duration from 4 UAVs, due to the two obstacles which are too close to allow a biconnectivity. For that reason, some links between the UAVs have to be broken, while an attractive force between them is still applied. From 4 UAVs, two links have to be broken and only one in the case of 3 UAVs.

Moreover, in all the simulations, the connectivity of the swarm was maintained in the tunnel.

CONCLUSIONS AND FUTURE WORK

In this paper we propose original mobility strategies based on virtual forces for a swarm of autonomous UAVs. Following these strategies, the UAVs, all having equivalent roles, can autonomously fulfill a surveillance mission of two AoIs separated by a narrow passage. To travel from the initial AoI to the other, the swarm organizes itself as a compact formation favoring communication links between the UAVs, which travel through the tunnel while avoiding obstacles unknown prior to the mission. We have run many simulations and evaluated them using a number of criteria. Our results show that our approach gives very good results.

Nevertheless, numerous topics have to be further explored. First of all, the UAV safety distance should depend on the range of the embedded sensor used to detect the obstacles, as well as on the speed of the UAVs. Furthermore, each UAV that detects an obstacle should share its location with the other aircrafts of the swarm. Thanks to this information, the UAVs could calculate a new target points taking the obstacle into consideration. The presence of obstacles on the AoI could also be considered. A wide subject of study could be the self-organization of the swarm in order to make up a compact formation depending on the tunnel width, instead of the deterministic configuration presented in this work. Finally, magnitudes of the three forces could be defined by other functions, as polynomial or exponential. This could improve the model strategies by speeding up surveillance or tunnel crossing.

Figure 1 :

 1 Figure 1: Representation of the two Areas of Interest separated by a narrow passage (a tunnel).

Figure 2 :

 2 Figure 2: Graphical representation of the attractiverepulsive force between UAVs. On the schemes, UAVs are represented by discs, attractive and repulsive forces are represented by open and triangle arrows respectively. The lines linking the UAVs represent communication links. Note: These notations will be used all along the paper.

Figure 3 :

 3 Figure 3: Graphical representation of the attractive force to a goal. The UAV target is represented by nested circles. Note: This notation will be used all along this paper.

Figure 4 :

 4 Figure 4: Graphical representation of the repulsive force due to obstacles. The obstacle is represented by a line.

Figure 5 :

 5 Figure 5: Scheme of the repulsive and attractive forces applied to a UAV in a swarm, in a narrow passage.

Figure 6 :

 6 Figure 6: Movements of a swarm composed of two UAVs during the surveillance mission over an AoI.

Figure 7

 7 Figure 7: a) Too narrow passage for a biconnectivity. This will be referred to as Tunnel shape a. b) 2-row staggered biconnected formation. (Tunnel shape b.) c) Biconnected formation for larger passages. (Tunnel shape c.).

Figure 8 :

 8 Figure 8: Case when an attractive force is not taken into account (dotted open arrow). The forces are plotted for the UAV of interest.

Figure 9 :

 9 Figure 9: A limitation of the model. The forces are plotted for only two UAVs for the sake of visibility. The two upper UAVs do not go backward for joining the others.

  of one UAV in the tunnel by the shortest path None Equal to the number of disconnections A view of the AoIs, tunnel and UAVs extracted from JBotSim is shown on the Figure 10.

Figure 10 :

 10 Figure 10: View of the two AoIs and of the passage referred to as Tunnel d.

Figure 11 :

 11 Figure 11: Passage referred as Tunnel e.

Figure 12 :

 12 Figure 12: Mean compact formation setup duration for each simulated configuration, normalized.

Figure 13 :

 13 Figure 13: Movements of a 6-UAV swarm aiming at entering in the tunnel. At the time t , the swarm is in a local equilibrium at the tunnel entrance.

Figure 14 :

 14 Figure 14: Tunnel crossing duration for each simulated configuration, normalized by the duration of tunnel crossing by a single UAV by the shortest path at maximum speed.

Table 1 :

 1 Evaluation criteria and ideal result for each step of the mission.

		Operational Objective	MoPs	Ideal Result
	Overall mission	Safeguard of the UAV	Distance between a UAV and another object	Greater than a given safety distance
	Surveillance mission	links Deploy over the AoI with communication	Number of connections Duration of an AoI coverage Number of disconnections	Inversely proportional to the number of UAVs None

Table 2 :

 2 Mean and standard deviation (SD) of the first AoI surveillance duration for various numbers of UAVs in the swarm, calculated on 150 runs for each number of UAVs.

	Number of UAVs	Duration of surveillance (number of steps) Mean SD	Speed Up
	1	1210	7	1
	2	764	14	1.6
	3	555	5	2.2
	4	425	8	2.8
	5	416	4	2.9
	6	233	11	5.2

Table 3 :

 3 Mean and standard deviation (SD) of the number of links lost and created during the surveillance phase, calculated on 150 runs for each number of UAVs.

	Number of	Lost Links	Created Links
	UAVs	Mean	SD	Mean	SD
	2	0.05	0.26	0.05	0.26
	3	0.09	0.29	0.09	0.30
	4	0.77	2.33	0.76	2.32
	5	2.03	4.80	2.02	4.77
	6	3.51	5.41	4.49	5.32


Ema Falomir

Serge Chaumette

Gilles Guerrini

Gilles Guerrini Mobility

Mobility Strategies for Swarms of Unmanned Aerial Vehicles using Artificial Potential Fields and Global Path Planning

come

Mobility Strategies for Swarms of Unmanned Aerial Vehicles using Artificial Potential Fields and Global Path Planning

Ema Falomir, Serge Chaumette, Gilles Guerrini

Motion Decision

Mission Each UAV has a mission represented by its own objective map, evolving all along the mission.

The UAVs collaborate to create a shared constraints map, updated all along the mission. Our Approach -The environment is discretized into square cells.

-The UAVs trajectories depend on the cells width.

-Addition of a potential field for the obstacle avoidance anticipation to the repulsive field around the obstacles (see figure below).

-Calculation of a global path for several iterations.

-Speed proportional to distance of next waypoint.

Objectives Whole project

Allow autonomous Unmanned Aerial Vehicles (UAVs) to quickly perform collaborative tasks, such as wide area surveillance. The UAVs communicate between each other and have similar characteristics: they make a swarm. They are considered as a unique entity as seen by an operator permitting to decrease workload of the operator.

PhD goal

Develop & validate a distributed mobility model for a swarm of autonomous UAVs which takes into account the embedded sensors capacities to allow the realization of a mission.

Artificial potential field resulting from our proposal, showing an anticipation of collision avoidance.

  @u-bordeaux.fr Pr. Serge Chaumette, Director schaumette@labri.fr Dr. Gilles Guerrini, industrial coordinator gilles.guerrini@fr.thalesgroup.com Publication & Patent E. Falomir, S. Chaumette and G. Guerrini. Mobility Strategies based on Virtual Forces for Swarms of Autonomous UAVs in Constrained Environments. 14th International Conference on Informatics in Control, Automation and Robotics, July 2017. E. Falomir, G. Guerrini, P. Garrec, Essaim constitué d'une pluralité de drones volants légers. Trajectories of 3 UAVs composing a swarm in an unknown environment.

  moves within an APF and goes towards the lowest potential.-The APF is a combination of an attraction towards a target and a repulsion from obstacles.


Philippe Nouillot

Françoise Natta
email: francoise.natta@insep.fr

Influence of velocity on the human global postural strategies during the movement leading up to the vertical upside-down position

Keywords: Anticipatory postural adjustments, Balance, Biomechanical parameters, Motor strategy

Influence of velocity on the human global postural strategies during the movement leading up to the vertical upside-down position

Philippe Nouillot, Françoise

When a voluntary movement is carried out, it is well established that postural adjustments generally appear before its initiation. In segmental movements performed in the upright stance, the local (segmental) and global (acceleration of center of mass (CoM), ground reaction, center of pressure (CoP)) anticipatory postural adjustments (APA) vary with the parameters of intentional movement, such as the localization, inertia and velocity [START_REF] Bouisset | Anticipatory postural movements and dynamic asymmetry[END_REF][START_REF] Horak | The effects of movement velocity, mass displaced and task certainty on associated postural adjustments made by normal and hemiplegic individuals[END_REF].APA also depend on the stability of equilibrium. In the erect postural stance, when the postural base is more unstable, either because of the reduced support area [START_REF] Aruin | Anticipatory postural adjustments in conditions of postural instability[END_REF][START_REF] Zattara | Posturo-kinetic organization during upper limb movement performed from varied foot positions[END_REF] or owing to a movable support surface [START_REF] Cordo | Properties of postural movements related to a voluntary movement[END_REF], the APA are increased. According to the hypothesis proposed by Bouisset and Zattara [START_REF] Bouisset | Segmental movement as a perturbation to balance? Facts and concepts[END_REF] that the dynamic of the intentional movement is a perturbation to postural equilibrium, the APA, which are dynamic phenomena, serve to minimize the subsequent postural destabilization. In the upper segmental ('focal chain') movements, the initial postural base, characterized by the support area and the 'postural chain', remains the same.

In movements with body progression the APA serve moreover to create appropriate global dynamics, allowing the displacement of the CoM, such as in step or gait initiation [START_REF] Brenière | Are dynamic phenomena prior to stepping essential to walking?[END_REF], and in lower limb elevation or flexion [START_REF] Do | Is balance or posture at the end of voluntary movement programmed?[END_REF][START_REF] Mouchnino | Séquence posturo-cinétique pendant l'élévation d'une jambe: timing des événements et effets de l'entraînement[END_REF][START_REF] Rogers | Dynamic transitions in stance support accompanying leg flexion movements in man[END_REF]. In these movements the duration and amplitude of anticipatory global postural dynamics increase with the speed of movement. However, in certain cases: [START_REF] Aruin | Anticipatory postural adjustments in conditions of postural instability[END_REF] when the movement is performed under global mechanical constraints, such as the CoP magnitude in the initiation of step, although the speed varies, APA can remain the same [START_REF] Ito | Anticipatory control in the initiation of a single step under biomechanical constraints in humans[END_REF]; (2) when the equilibrium state is too precarious, the APA can be absent [START_REF] Pai | Center of mass velocity-position predictions for balance control[END_REF]. In this last case, it was hypothesized that this absence of APA can correspond to another strategic response of the central nervous system (CNS) [START_REF] Aruin | Anticipatory postural adjustments in conditions of postural instability[END_REF]. The global dynamic strategy can also change with speed. For example, in lower leg flexion, when the movement is performed very slowly, the direction of the horizontal ground reaction force under the stance leg is reversed [START_REF] Rogers | Dynamic transitions in stance support accompanying leg flexion movements in man[END_REF]. In the movements with body progression, when the postural base is continuously changing, as in locomotion, after the loss of balance, the CoM must be continuously in motion, whereas in lower limb movements, in which the postural base is transitorily changing, the CoM's speed must be stopped and balance restored. In these latter movements, horizontal speed can be considered as a factor responsible for the greater difficulty in controlling the final balance. Indeed, in erect stance in the quasi static equilibrium on the floor, the gravitational line must be projected inside the support area, according to the concept of 'CoM position limits'. On the other hand, when the CoM is in horizontal movement before the final stance, the concept of 'CoM velocity-position limits', proposed by Pai and Patton [START_REF] Rogers | Dynamic transitions in stance support accompanying leg flexion movements in man[END_REF], must then be used. This concept is grounded on the 'functional' aspect of the support area, which is reducing conversely with the speed of movement. Therefore, in movement with body progression, when the speed of CoM is increasing, it can considered that the final balance is more precarious.

In acrobatic sports, there are numerous tasks with transient postural bases in which the final balance is more precarious than in the upright stance. This is the case when a gymnast is standing on his hands on the floor after a reverse movement. Compared to the normal standing posture, the supporting plane is smaller, and the forearm muscles are likely less strong. In this kind of movement, the final postural balance can be compared to an inversed pendulum [START_REF] Clément | Motor behavior underlying the control of an upside-down vertical posture[END_REF][START_REF] Pai | Center of mass velocity-position predictions for balance control[END_REF].

The aim of this study was to answer the following question: which global motor strategies are involved in movement in which final balance is particularly precarious?

Seven subjects, all expert gymnasts, five men and two women from 17 to 22 years (mean age 19.2 ± 2.2 years), participated in the experiment with their full consent. They performed an acrobatic movement from an initial posture with both hands and feet on the ground (Fig. 1).

The hands were placed on the lateral axis, flat on the ground. Feet were placed along the anteroposterior axis, one in front of the other so that the rear lower limb was straight and the front one flexed. The corresponding task was to place the body in an upside-down posture. The final posture, in equilibrium on the hands, had to be held for 1-2 s. The aim of the task was, after an auditory signal, to raise up the body on the hands by projecting the straight lower limb, and pushing on the ground with the flexed limb. This projection was considered as the onset of intentional movement (t0) and was determined from a mono-axial accelerometer placed on the heel. For each test the subject repositioned the feet and hands in the same place. Two experimental conditions were considered, distinguished by the velocity: (1) natural velocity (normal condition: NV); (2) fast velocity (FV), i.e. as fast as possible. About ten tests per condition had to be performed. Between each trial, the subjects had a short rest. For safety reasons, an operator was positioned in front of the subject so as to retain him in case there was a dorsal imbalance, leading to a fall or a hand movement. The interest in this kind of movement concerns the initial stand base, particularly wide (in quadruple stance), and the final stand base, particularly narrow (manual stance).

As in a previous paper [START_REF] Nouillot | Caractérisation dynamique d'un mouvement acrobatique: le renversement en appui tendu sur les mains[END_REF] the method of analysis was based on data exploitation from a force plate of large dimensions, so that the subjects could have the four supports on it without the force plate being a constraint. From the force plate, which displayed the CoP and the ground reaction forces, the global dynamics applied to the CoM were obtained, according to Newton's laws. The biomechanical variables taken into account in our results were accelerations of the CoM and speeds obtained by simple integration of accelerations (Euler's method) on the antero-posterior (x″, x′) and vertical (z″, z′) axes. The acquisition of data was realized with a period of sampling rate of 750 Hz. The statistic methods used were t-test and paired t-test. The significant differences were considered when P < 0:05. Concerning the performance, i.e. the speed of CoM, according to the order given to subjects the peaks of speed showed greater values in FV than NV, but only on the vertical axis (Table 1). Indeed, whereas speed values were about 16% greater on the vertical axis, they were similar for the two conditions on the antero-posterior axis.

Concerning the description of the task from mechanical traces, the first variations of traces (t-1, Figs. 1 and2) occurred systematically before the onset of intentional movement (t0), which was at the time of the break-down observed on the traces, generally on the anteroposterior axis (see Table 2). This corresponded to the onset of a purely postural phase, i.e. to a global APA (APA phase). Before t 2 1 there was no variation of the traces, corresponding to a period of motor latency when the subject remained motionless (motionless phase). Between t0 an t1 the subject was on three supports, the two hands and one foot, the rear foot being up from the floor. The body was turning upside down. Therefore, during this phase as the CoM was only forwards accelerated, and simultaneously the displacement of the CoP backwards directed, we can consider that the subject was pushing backwards on the floor with the flexed limb (pushing phase). From t-1 to t1, which is the sum of the APA and pushing phases, on the antero-posterior axis (x) the variation of the acceleration was forwards directed, and simultaneously the CoP backwards directed (amplitudes of displacement: between 57 ± 16 and 359 ± 19 mm). This period corresponds to a propulsive phase. In the same time, the vertical acceleration (z) was generally upwards directed (systematically for five subjects). However, during the APA phase vertical acceleration was slightly downwards directed only for two subjects. For the anteroposterior and vertical acceleration traces there were one or two peaks, depending on the subject.

The individual variations on the two axes have been more precisely described in a preliminary study [START_REF] Nouillot | Caractérisation dynamique d'un mouvement acrobatique: le renversement en appui tendu sur les mains[END_REF]. At t1, the CoM was beginning to be backwards and downwards accelerated, which corresponds to an onset of a breaking phase along both the axes. Consequently, from this time the speeds of the CoM were reducing. This time can be considered as the take-off of the front foot. The subject was then only on his hands. The period from t1 to the end of CoM displacement corresponds to a tilting phase.

Concerning the durations, the comparison between the two experimental conditions showed that they were significantly shorter in FV for all phases (Table 2). The differences were about 23% in the APA phase on the antero-posterior axis and about 28%on the vertical axis, about 27% in the body turning phase on the antero-posterior axis and about 37% on the vertical axis, and about 26% in the body propulsive phase on both the antero-posterior and vertical axes.

Concerning the peak magnitudes of the acceleration, phases must be separately considered (Table 2). For the APA phase, on the antero-posterior axis, peaks were significantly greater in FV than in NV. In this axis five subjects presented positive peaks, and two subjects presented negative peaks. However, for these two categories of subjects the comparison of the peaks between the two experimental conditions was also significantly different. For the group of five subjects see Table 2. For the two other subjects, differences were P = 0.001 for the first subject and P = 0.02 for the second subject (t-test). For the pushing phase on the antero-posterior axis, there was no difference between FV and NV, whereas on the vertical axis peaks were greater in FV than NV.

The purpose of this study was to examine the influence of the velocity on the global postural organization. In the movement of standing up into an upside-down posture, we were interested by the dynamic phases preceding the initiation of the intentional movement, i.e. to APA, as well as the phase of execution, which precedes the tilting movement towards the final balance position of the body. Concerning the speed of movement, subjects had to perform the movement with a spontaneous velocity, and as fast as possible. In the FV condition, the speed of CoM was significantly increased compared to the NV. What is particularly interesting is that speed varied only on the vertical axis. The equilibrium problems must therefore be considered separately on each axis.

On the vertical axis, it had been hypothesized in a preliminary study of a lower limb flexion performed in erect stance that it was not necessary to preserve equilibrium [START_REF] Nouillot | Are there anticipatory segmental adjustments associated with lower limb flexions when balance is poor in humans?[END_REF]. For example, the body segments are stacked up and all the vertical displacements of the CoM remain above the support area, i.e. they do not change the gravity line projection contrary to an anteroposterior displacement.

On the antero-posterior axis, the control of body balance should depend upon the posturokinetic capacity (PKC) of the body system to minimize the postural destabilization [START_REF] Bouisset | Posturo-kinetic capacity assessed in paraplegics and parkinsonians[END_REF].

However, the PKC is lower on this axis compared to the vertical axis, as has been shown in normal erect stance. Displacements are limited by the support area, the muscle forces, and the possible joint mobilizations [START_REF] Kuo | A biomechanical analysis of muscle strength as a limiting factor in standing posture[END_REF][START_REF] Mccollum | Form and exploration of mechanical stability limits in erect stance[END_REF]. Results of studies have shown that dynamic variations depend on the axis. In a study of a forearm flexion performed in erect stance, it was reported by Friedli et al. [START_REF] Friedli | Postural adjustments associated with rapid voluntary arm movements. II. Biomechanical analysis[END_REF] that on the antero-posterior axis the acceleration created by the voluntary movement at the elbow level was considerably reduced at the feet stance level, whereas on the vertical axis the acceleration was totally transmitted. Thus, according to the hypothesis that a greater speed of theCoM reduces the 'functional' support area, we can assume that the CNS 'chooses' to conserve a spontaneous speed which perhaps corresponds to an 'optimal' and/or a maximal value, i.e. in particular compatible with the PKC of the neuromuscular system.

Concerning both the APA duration and amplitude, whilst amplitudes of mechanical traces were increasing with the speed, the durations were decreasing. This result is likely to depend on the specificity of tasks. In segmental movements, the aim of the APA is to counteract in advance the perturbation due to the forthcoming voluntary movement in order to maintain the initial postural equilibrium [START_REF] Bouisset | Anticipatory postural movements and dynamic asymmetry[END_REF]. In the initiation of gait the dynamic of movement is created using gravity. The greater the speed of CoM, the longer the time the gravity action acts, because gravity is a constant value. It has been demonstrated by Brenière et al. [START_REF] Brenière | Are dynamic phenomena prior to stepping essential to walking?[END_REF] that the horizontal speed of CoM was a programmed parameter. In the sprint start, in which initial posture is similar to the task of raising up on the hands, the increasing APA durations are linked firstly to the forthcoming antero-posterior speed and secondly eventually to the stand up on the feet, i.e. against gravity [START_REF] Natta | Le départ de sprint en starting-blocks: analyse des mouvements anticipateurs de la course[END_REF]. Although in our study the movement was also performed against gravity, the antero-posterior speed must be stopped, contrary to the sprint start. As the antero-posterior speed remained similar and the vertical speed was greater, the reduction of durations is therefore not linked to the reduction of the speed. This reduction of durations does not depend on the instability of initial balance which is obviously very stable, or on the initial global mechanical constraints [START_REF] Ito | Anticipatory control in the initiation of a single step under biomechanical constraints in humans[END_REF]. The APA reduction seems to be linked to the precariousness of the final balance. Indeed, according to the impulsion law (∑F∆t = m ∆V, with initial velocity = 0, because the subject was initially motionless), when the magnitude of the antero-posterior reaction force was increasing, the reduction of duration allowed the conservation of the 'optimal' speed. This dynamical strategy occurs during the pushing phase, which follows the APA phase. Therefore, as suggested by Bouisset et al. [START_REF] Bouisset | Are amplitude and duration of anticipatory postural adjustments identically scaled to focal movement parameters in human?[END_REF], the duration seems a parameter particularly sensitive to the postural condition.

During the pushing phase, a braking activity assumed by the upper limbs (which are postural limbs) was probably exerted. This assumption is relayed on the tilting phase traces, which follows the pushing phase, and on results from the dynamic of the postural lower limb in the study of lower limb flexion [START_REF] Zattara | Posturo-kinetic organization during upper limb movement performed from varied foot positions[END_REF]. However, this upper limb breaking activity should be more important during the pushing phase in FV, as global magnitude values remained similar (Table 3).

To sum up, according to a previous study, in which Do et al. [START_REF] Do | Is balance or posture at the end of voluntary movement programmed?[END_REF] had shown that final balance was a programmed parameter, from our results we can confirm that the greater or lesser stability of the final balance was a determining factor of the motor program. In our study, a specific global dynamic strategy was to increase the body speed only on the vertical axis, on which final balance, in terms of CoM position limits, does not exist. This strategy was characterized by a reduction of the APA duration, contrary to other studies in which APA were increasing with speed of movement. The global dynamic strategy has been allowed on the grounds of the postural stance characteristics during the body propulsive phase.

Table 1

Speed values of the CoM NV (mean ± SD) FV (mean ± SD) P value (t-test, n = 7)

x' (ms -1 ) 0.74 ± 0.12 0.77 ± 0.14 NS z' (ms -1 ) 0.89 ± 0.16 1.06 ± 0.17 <0.001

Values correspond to the maximum speed magnitude (m/s) of the CoM along the anteroposterior (x') and vertical (z') axes. NV, FV, movement performed respectively in natural and in fast velocity. The comparison between the two conditions is performed with the t-test. The values correspond to the durations of the different phases. APA phase (1): duration between the first variations of mechanical traces (t-1) and the onset of intentional movement (t0). The negative sign indicates that durations are before t0. Pushing phase (2): duration between the onset of intentional movement (t0) and the take-off of the front foot from the floor (t1). Propulsive phase: duration between the first variations of mechanical traces (t 2 1) and the take-off of the front foot from the ground (t1). This last phase corresponds to the sum of the absolute values (1 + 2). x″, z″: forwards and upwards accelerations of CoM, respectively. NV, FV, movement performed in natural and in fast velocity, respectively. Mean (of subjects) and standard deviations (SD) are expressed in milliseconds. APA phase: duration between the first variations of mechanical traces (t-1) and the onset of intentional movement (t0). Pushing phase: duration between the onset of intentional movement (t0) and the take-off of the front foot from the ground. x″ , z″, forward and upward accelerations of CoM, respectively. The comparisons between the two conditions were performed with the paired t-test ( a n = 7, b n = 5).

Figure 1 .

 1 Figure 1.

Figure 2 .

 2 Figure 2.

Table 2

 2

	Kinematic durations
	Phases	NV (mean ± SD)	FV (mean ± SD)	P value (t-test, n = 7)
	APA phase (1)
	x″ (ms)	-333 ± 53	-257 ± 54	<0.03
	z″ (ms)	-318 ± 67	-230 ± 61	<0.03
	Pushing phase (2)
	x″ (ms)	336 ± 75	246 ± 67	<0.04
	z″ (ms)	330 ± 61	208 ± 105	<0.03
	Propulsive phase (1 + 2)
	x″ (ms)	665 ± 116	490 ± 86	<0.01
	z″ (ms)	627 ± 127	467 ± 72	<0.02

Table 3 Differences of kinematic magnitudes between NV and FV

 3

	Phases	NV (mean ± SD)	FV (mean ± SD)	FV -NV (mean ± SD)	P value (paired t-test)
	APA phase
	x″ (ms -2 )	1.45 ± 0.99	2.72 ± 1.40	1.19 ± 0.64	<0.01 a
	z″ (ms -2 )	1.16 ± 0.54	2.51 ± 0.84	1.35 ± 0.39	<0.01 b
	Pushing phase
	x″ (ms -2 )	2.49 ± 1.07	2.60 ± 1.63	0.11 ± 0.61	NS a
	z″ (ms -2 )	3.66 ± 1.05	5.87 ± 1.87	2.21 ± 1.28	<0.01 a