M Collet

M Ouisse

M Ichchou

M Ruzzene

Semi-active optimization of 2D wave's dispersion into mechanical systems by the mean of periodically distributed shunted piezoelectric patches: a new class of adaptive metamaterials

Keywords: Distributed control, 2D Waves Dispersion, Bloch Theorem, Shunted piezoelectric, Mid-Frequency Optimization

In this paper, we present an application of the Floquet-Bloch theorem in the context of electrodynamics for vibroacoustique power flow optimization by mean of distributed and shunted piezoelectric patches. The main purpose of this work is first to propose a dedicated numerical approach able to compute the multi-modal wave dispersions curves into the whole first Brillouin zone for periodically distributed 2D shunted piezo-mechanical systems. By using a specific indicator evaluating the evanescent part Bloch's waves, we optimize, in a second time, the piezoelectric shunting electrical impedance for controlling energy diffusion into the proposed semi-active distributed set of cells. A 3D modeling of semi-distributed distribution of the optimal smart metamaterial is used for validating the obtained cell design.

INTRODUCTION

Tailoring the dynamical behavior of wave-guide structures can provide an efficient and physically elegant approach for optimizing mechanical components with regards to vibration and acoustic criteria, among others. However, achieving this objective may lead to different outcomes depending on the context of the optimization. In the preliminary stages of a product's development, one mainly needs optimization tools capable of rapidly providing global design directions. Such optimization will also depend on the frequency range of interest. One usually discriminates between the low frequency (LF) range and the medium frequency (MF) range, especially if vibration and noise are considered. However, it should be noted that LF optimization of vibration is more common in the literature than MF optimization. For example, piezoelectric materials and other adaptive and smart systems are employed to improve the vibroacoustic quality of structural components, especially in the LF range [START_REF] Preumont | Vibration control of structures : An introduction[END_REF][START_REF] Nelson | Active Control of Sound[END_REF][START_REF] Banks | Smart material structures Modeling Estimation and Control[END_REF] even if ditributed transducers are used. [START_REF] Batra | Multimode vibration suppression with passive twoterminal distributed network incorporating piezoceramic transducers[END_REF][START_REF] Tzou | A study of segmentation of distributed piezoelectric sensors and actuators, part i: Theorical analysis[END_REF] Recently, much effort has been spent on developing new multi-functional structures integrating electro-mechanical systems in order to optimize their vibroacoustic behavior over a larger frequency band of interest. [START_REF] Collet | Wave Motion Optimization in Periodically Distributed Shunted Piezocomposite Beam Structures[END_REF][START_REF] Thorp | Attenuation and localization of wave propagation in rods with periodic shunted piezoelectric patches[END_REF][START_REF] Beck | Broadband vibration suppression assessment of negative impedance shunts[END_REF][START_REF] Beck | Experimental analysis of a cantilever beam with a shunted piezoelectric periodic array[END_REF][START_REF] Bartoli | Modeling wave propagation in damped waveguides of arbitrary cross-section[END_REF][START_REF] Foda | Vibration suppression of composite laminated beams using distributed piezoelectric patches[END_REF][START_REF] Dawood | Three-dimensional static shape control analysis of composite plates using distributed piezoelectric actuators[END_REF][START_REF] Jiang | Decentralized robust vibration control of smart structures with parameter uncertainties[END_REF][START_REF] Ren | Finite difference based vibration simulation analysis of a segmented distributed piezoelectric structronic plate system[END_REF] However, there is still a lack of studies in the literature for MF optimization of structural vibration. To that end, the focus of this study is to provide a suitable numerical tool for computing wave dispersion in 2D periodic systems incorporating controlling electronics devices. The main final aim is to allow their optimization in order to optimize vibroacoustic diffusion in 2D wave guides. Two numerical approaches can be distinguished for computing that dispersion: the semi-analytical finite element method (SAFE) and the wave finite element (WFE) method. The main disadvantage of the SAFE method is that FE used are not standard so they must be specifically defined for each application. Nevertheless, a large amount of FE has been developed since 1975 to compute dispersion curves of rails, 15 laminated composite plates [START_REF] Mace | Modelling wave propagation in two-dimensional structures using finite element analysis[END_REF][START_REF] Gonella | Analysis of in-plane wave propagation in hexagonal and re-entrant lattices[END_REF] and viscoelastic laminated composite plates. [START_REF] Akrout | Comportement dynamique dterministe et large bande des structures guides[END_REF] To avoid development of specific FE, the WFE method considers the structures as periodic in order to model, with standard FE, a period of the structure. By using the periodic structure theory (PST) introduced by Mead, [START_REF] Mead | A general theory of harmonic wave propagation in linear periodic systems with multiple coupling[END_REF] an eigenvalue problem can be formulated from the stiffness and mass matrices of the FE model to find wave numbers of all the propagating waves. Contrary to SAFE method, the displacement field is now approximated in the direction of propagation. Thus, some numerical issues can arise when the size of FE are too coarse. As recommended by Mace and Manconi, 16 a minimum of six elements per wavelength is a good rule of thumb to ensure a reliable analysis. The WFE method has been successfully used to deal with wave propagation in two dimensional structures. [START_REF] Manconi | The Wave Finite Element Method for 2-dimensional Structures[END_REF][START_REF] Berthaut | Multi-mode wave propagation in ribbed plates : Part i k-space characteristics[END_REF] One of the main problem all these approaches is the difficulty to compute the damped wave numbers in the whole Brillouin domain necessary for optimizing vibroacoustic behavior of smart periodic structures. After recalling the Floquet-Bloch theorems, we introduce a new numerical formulation for computing the multimodal damped wave numbers dispersion in the whole first Brillouin domain of a periodical smart structure made of periodically distributed shunted piezoelectric patches. Based on this wave modeling, optimization of the electrical impedance of the shunted circuit is made in order to decrease group velocity of flexural waves. The obtained optimal impedance is also tested in controlling the HF response of a semi-distributed system.

PIEZO-ELASTO-DYNAMICAL APPLICATION OF THE FLOQUET-BLOCH THEOREM

In this section the application of the celebrated Floquet-Bloch theorem is presented for piezo-elastodynamic problems. Based on the well known results obtained by Floquet 21 in one-dimensional and later rediscovered by Bloch 22 in multidimensional problems, we propose an original application to bi-dimensional piezo-elastodynamical problem leading to very general numerical implementation for computing waves dispersion for periodically smart distributed mechanical systems incorporating electronic components and damping effects. [START_REF] Collet | Numerical tools for semi-active optimization of 2d waves dispersion into mechanical system[END_REF]

The Bloch Theorem

The Bloch theorem gives the form of homogeneous states of Schrödinger equation with periodic potential. This theorem can be considered as a multidimensionnal application of the Floquet theorem. [START_REF] Joannopoulos | Photonic Crystals: Molding the Flow of Light[END_REF] The periodic medium (or potential) properties satisfy M (x + R.m) = M (x), m ∈ Z 3 where R = [r 1 , r 2 , r 3 ] ∈ R 3×3 is a matrix grouping the three lattice's basis vectors (in 3D). We can also define the primitive cell as a convex polyhedron of R 3 called Ω x . The reciprocal unit cell is denoted by Ω k limited by the reciprocal lattice vector defined by the three vectors g j so that: r i .g j = 2πδ ij (δ i,j the Kronecker index). We note G = [g 1 , g 2 , g 3 ] the reciprocal lattice matrix in the later. If Ω x is the irreductible primitive cell, Ω k corresponds to the first Brillouin zone of the lattice. One can see [START_REF] Kittel | Introduction to Solid State Physics[END_REF] for details.

The Bloch Theorem stipulates that any functions u(x) ∈ L 2 (R 3 , C n ) can be expressed as

u(x) = Ω k e ikx ũ(x, k)dk (1)
where the Bloch amplitude ũ(x, k) is Ω x -periodic and has the representations

ũ(x, k) = n∈Z 3 û(k + Gn)e iGn.x , u(x) = |Ω x | (2π) 3 n∈Z 3 u(x + Rn)e ik(x+Rn) (2)
where û(k) stands for the Fourier transform of u(x). One can also demonstrate that the mean value of the Bloch amplitude is the Fourier amplitude of u(x) for the corresponding wave vector: ũ(., k) Ωx = û(k). Using the Bloch theorem to represent the solutions of periodical partial derivative equations implies that all derivatives are shifted by k in the sense given by the used spatial operator.

Based on that theorem one can define the expansion functions v m (x, k), called the Bloch eigen modes, such that they can be used to represent the Bloch amplitudes of any solution of the corresponding partial derivative equation as

ũ(x, k) = m u m (k)v m (x, k) ( 3 )
and at the same time diagonalize the partial derivative equations. One notes that the expansion coefficients u m (k) depend on the applied disturbance and also on the induced wave vector (see 26 for details). Let us consider a piezo-elastodynamic problem made of infinite periodic distribution of unitary cell described in figure 1. The harmonic homogeneous dynamical equilibrium of system is driven by the following partial derivative equation :

ρ ẅ(x) -∇σ(x) = 0 ∀x ∈ Ω x -∇D(x) = 0 ∀x ∈ Ω x ( 4
)
where w(x) ∈ R 3 (Ω x ) is the displacement vector, σ represents the Cauchy stress tensor, = ∇ sym w = 1 2 (∇w T (x) + w(x)∇ T ) the Green strain tensor, D(x) the electric displacement. The linear constitutive material behavior relationships can be written as

σ = C E (x) -e T (x)E (5) D = e(x) + ε S (x)E (6)
where E = -∇V the electric field vector (V the voltage), C E the elasticity tensor at constant electrical field, e T the piezoelectric coupling tensor and ε S the dielectric permittivity at constant strain. We add to this set of equilibrium equations an output expression

q o = - St D.ndS (7)
allowing the introduction of the charge measurement on the piezoelectric's top electrode and hence the dual counterpart of the imposed electrical Dirichlet boundary condition for applying the shunt impedance operator. The equations above are consistent for each kind of material to the extent that null piezoelectric and permittivity tensors can be used when passive materials are considered. All of these tensors also depend on the spatial location vector x. The piezo-elastodynamic equilibrium can also be written as :

ρω 2 w(x) + ∇C∇ sym (w(x)) + ∇e T (x)∇V (x) = 0 ∀x ∈ Ω x (8) -∇e(x)∇ sym (w(x)) + ∇ε S (x)∇V (x) = 0 ∀x ∈ Ω ( 9 )
As the problem is 2D infinitely periodic, only electrostatic boundary conditions have to be considered on each cell: (10) where S b is the grounded bottom electrode of the piezoelectric layer, S t is the top electrode connected to the external shunt and s l the lateral electrode less boundary. The top electrode applied feedback voltage V o depends on the shunt characteristic and on the collected charges q o (7) and can be expressed in the Fourier space by:

⎧ ⎨ ⎩ V = 0 ∀x ∈ S b V = V o ∀x ∈ S t D.n = 0 ∀x ∈ S l
V o (iω) = -Z(iω)q o (iω) (11)
By considering a primitive cell of the periodic problem Ω x and by using the Bloch theorem, we can compute the associated Bloch eigenmodes (3) and the dispersion functions by searching the eigen solutions of the homogeneous problem ( 8) and ( 9) as:

u(x) = w(x) V (x) = u n,k (x)e ik.x (12) with u n,k (x) = w n,k (x) V n,k (x)
, Ω x periodic functions. By introducing expression (12) in the piezo-elastidynamic equations ( 8), ( 9), one can demonstrate that w n,k (x), V n,k (x) and ω n (k) are solutions of the generalized eigenvalues problem:

ρω 2 n (k)w n,k (x) + ∇C∇ sym (w n,k (x)) + ik {(C∇ sym (w n,k (x))).Φ + ∇(CΞ n,k (x))} -k 2 (CΞ n,k (x)).Φ +∇e T ∇V n,k (x) + ik (∇e T V n,k (x)).Φ + (e T ∇V n,k (x)).Φ -k 2 V n,k (x)(e T Φ).Φ = 0 ∀x ∈ Ω x ( 13
)
-∇e∇ sym (w n,k (x)) -ik {∇(eΞ n,k (x)) + (e∇ sym (w n,k (x))).Φ} + k 2 (eΞ n,k (x)).Φ +∇ε S ∇V n,k (x) + ik {(∇ε S V n,k (x)).Φ + (ε S ∇V n,k )(x).Φ} -k 2 (ε S ΦV n,k (x)).Φ = 0 ∀x ∈ Ω x ( 14
)
with the associated bounadry conditions :

⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ w n,k (x -R.m) = w n,k (x) ∀x ∈ S r m ∈ Z 2 V n,k (x) = 0 ∀x ∈ S b V n,k (x)V = -Z(iω)q o n,k ∀x ∈ S t D.n = 0 ∀x ∈ S l
where k = k ⎡ ⎣ cos(φ) sin(φ) 0 ⎤ ⎦ = kΦ where φ represent the direction angles into the reciprocal 2D lattice domain and

Ξ n,k (x) = 1 2 (w n,k (x).Φ T + Φ.w T n,k (x)
) the symmetric dyadic tensor or the dyadic product of the displacement w n,k (x) and direction vector Φ. S r are the interface of the cells continuum, and R the matrix grouping the two lattice's basis vectors (in 2D in the considered problem). In the electrical boundary conditions, q t n,k is given by :

q o n,k = St [-e(∇ sym (w n,k (x)) + ik∇eΞ n,k (x)) + ε S (∇V n,k (x) + ikV n,k (x)Φ].ndS ( 15
)
where n is the outpointing unitary normal vector. The proposed formulation is also based on the computation of the Floquet vectors (equation ( 13), ( 14)), instead of computing the Floquet propagators commonly used for elastodynamic applications. Our approach allows to obtain the full 2D waves dispersions functions and to clearly introduce damping and electrical impedance into the piezo-elastodynamic operator. The adopted methodology allows the computation of the complete complex map of the dispersion curves incorporating computation of evanescent waves and allowing the introduction of damping and shunt operator if any.

Weak Formulation and computation of waves dispersion functions in periodical piezo-composite lattice

Let us consider the partial derivative equations ( 13), ( 14) on a unit cell Ω. It stands for a generalized eigenvalue problem leading to compute the dispersion functions ω n (k) and the corresponding Floquet eigenvectors u n,k (x).

For computing the 2D dispersions curves, we need to introduce a suitable weak formulation. If u n,k (x) is a solution of equations ( 13), ( 14), also

∀ wn,k (x) ∈ H 1 (Ω, C 3 )/ wn,k (x -Rm) = wn,k (x) ∀x ∈ S r and Ṽn,k (x) ∈ H 1 (Ω, C)/ Ṽn,k (x) = 0 ∀x ∈ S b and Ṽn,k (x) = V o ∀x ∈ S t we have : Ω ρω 2 n (k) wn,k (x)w n,k (x) -(ε n,k (x) -ik Ξn,k (x))C(ε n,k (x) + ikΞ n,k (x)) +(ε n,k (x) -ik Ξn,k (x))e T (∇V n,k (x) + ikV n,k (x)Φ) -(∇ Ṽn,k (x) -ik Ṽn,k (x)Φ)e(ε n,k (x) + ikΞ n,k (x)) +(∇ Ṽn,k (x) -ik Ṽn,k (x)Φ)ε S (∇V n,k (x) + ikV n,k (x)Φ)dΩ - Ṽ t n,k V t n,k Z(iω) = 0 ( 16
)
This weak formulation is simply obtained by integrating equation ( 13), ( 14) projected onto any test function ũn,k (x). The boundary integral vanishes as the test functions are chosen so that wn,k (x -Rm) = wn,k (x) on S r . For a polyhedron cell, each boundary is generally a polyhedral plane sub-domain that can be associated with a parallel opposite one. The symmetry conditions called wn,k (x -Rm) = wn,k (x) explicitly link these associated surfaces.

Numerical Computation of the Bloch's waves

The numerical implementation is obtained by using a standard finite elements method to discretize the weak formulation (16). The assembled matrix equation is given by:

(K(Z(iω n (λ, φ)) + λL(φ, Z(iω n (λ, φ))) -λ 2 H(φ, Z(iω n (λ, φ))) -ω 2 n (λ, φ)M )u n,k (φ) = 0, (17)
where λ = ik, M and K(Z(iω n (λ, φ))) are respectively the standard symmetric semi-definite mass and stiffness matrices (the mass matrix is semi definite because elastostatic equation are condensed into the equation), L(φ, Z(iω n (λ, φ))) is a skew-symmetric matrix and H(φ, Z(iω n (λ, φ))) is a symmetric semi-definite positive matrix.

When k and φ are fixed and Z does not depend on ω the system ( 17) is a linear eigen value problem allowing us to compute the dispersion functions ω 2 n (k, φ) and the associated Bloch eigenvector u n,k (φ).

This approach has been widely used for developing homogenization techniques and spectral asymptotic analysis like in the work of. [START_REF] Allaire | Bloch waves homogenization and spectral asymptotic analysis[END_REF] It can also be applied for computing wave's dispersion even if Floquet propagators is preferred for 1D or quasi 1D computation, as indicated in. [START_REF] Ichchou | Guided waves group and energy velocities via finite elements[END_REF][START_REF] Houillon | Wave motion in thin-walled structures[END_REF][START_REF] Mencik | Multi-mode propagation and diffusion in structures through finite elements[END_REF] Nevertheless these approaches have been only developed for undamped mechanical systems that is to say represented by a set of real matrices. In this case, most of the previously published works present techniques based on the mesh of a real k-space (i.e k or λ and φ) inside the first Brillouin zone for obtaining the corresponding frequency dispersion and the associated Floquet vectors. For undamped system only propagative or evanescent waves exist corresponding to a family of eigen solutions purely real or imaginary. Discrimination between each class of waves is easy. If a damped system is considered (K, L, H are complex frequency dependent) or frequency dependence of the electrical shunt impedance is considered, the obtained eigenvalue problem is not quadratic and a complex specific numerical methodology has to be implemented. Furthermore, evanescent part of propagating waves appear as the imaginary part of ω 2 n (λ, Φ). It then becomes very difficult to distinguish the propagative and evanescent waves but also to compute the corresponding physical wave's movements by applying spatial deconvolution.

Another much more suitable possibility for computing damped system, dedicated for time/space deconvolution and for computation of diffusion properties as defined by, [START_REF] Collet | Wave Motion Optimization in Periodically Distributed Shunted Piezocomposite Beam Structures[END_REF][START_REF] Mencik | Multi-mode propagation and diffusion in structures through finite elements[END_REF] is to consider the following generalized eigen value problem:

(K(Z(ω) -ω 2 M + λ n (ω, φ)L(φ, Z(ω)) -λ 2 n (ω, φ)H(φ, Z(ω)))u n (ω, φ) = 0. ( 18
)
In this problem, the pulsation ω is a real parameter corresponding to the harmonic frequency. Wave's numbers and Floquet vectors are then computed. An inverse Fourier transformation in the k-space domain can lead us to evaluate the physical wave's displacements and energy diffusion operator when the periodic distribution is connected to another system as in. [START_REF] Collet | Wave Motion Optimization in Periodically Distributed Shunted Piezocomposite Beam Structures[END_REF] Another temporal inverse Fourier transformation can furnish a way to access spatio-temporal response for non-homogeneous initial conditions. As L is skew-symmetric, the obtained eigen values are quadruple (λ, λ, -λ, -λ) collapsing into real or imaginary pairs (or a single zero) when all matrices are real (i.e. for an undamped system). In this case a real pair of eigen values correspond to evanescent modes oriented in two opposite directions on the k-space and imaginary values to two traveling waves propagating in opposite direction.

As previously mentioned, the real part of k = kΦ vector is restricted to stand inside the first Brillouin zone. In the quadratic eigen value problem (18) nothing restricts computation to only find eigen values satisfying this condition. For direction vector Φ orthogonal to the lattice facelets (i.e. for Φ p1 = [1, 0] T and Φ p2 = [0, 1] T in bi-dimensional rectangular cell), the periodical conditions expressed for one dimensional wave guide are still valid: if λ j (ω, Φ p ) is an eigen value associated to w j (ω, Φ p ) then ∀m ∈ Z 3 , λ + i.Φ T p (G.m) is also an eigen value associated to w j (ω, Φ p )e -i.Φ T p (G.m)x . Thus, for undamped systems, all obtained eigenvalues are periodically distributed in the k-space along its principal directions.

Computation of the group velocity and evanescence criterion

The main aim of this paper is to provide a numerical methodology for optimizing the piezoelectric shunt impedance Z(ω) for controlling energy flow into the periodically distributed piezo-composite structure. For doing this, we need to define a suitable criterion. The waves group velocities indicate how energy is transported into the considered system and allow to distinguish the 'propagative' and 'evanescent' waves. If one Bloch eigen solution (i.e u n (ω, φ), k n (ω)) is considered, the associated group velocity vector 31 is given by :

C g n (ω, φ) = ∇ k ω = S e tot = I E tot ( 19
)
where : is the spatial and time average respectivelly on one cell and one period, S is the density of energy flux defined in, 31 I the mean intensity and e tot , E tot the total piezomechanical energy and its time average on a period (see [START_REF] Maysenhölder | Körperschall-energie Grundlagen zur Berechnung von Energiedichten und Intensitäten[END_REF] for details). In this problem, we only consider mechanical energy transportation as the electrostatic coupling is decentralized and can be condensed as a mechanical interface as proved in [START_REF] Hagood | Damping of structural vibrations with piezoelectric materials and passive electrical networks[END_REF] and generally computed in. [START_REF] Collet | Modal Synthesis and Dynamical Condensation Methods for Accurate Piezoelectric Systems Impedance Computation[END_REF] So we also compute the intensity vector I by :

I n = - ω 2 Re( Ωx C(ε n (x, ω, φ) + ikΞ n (x, ω, φ)).(w * n (x, ω, φ))dΩ/V ol ( 20
)
where . * is the complex conjugate and V ol the domain volume. As the spatio-temporal average of the system Lagragian is null (see [START_REF] Maysenhölder | Körperschall-energie Grundlagen zur Berechnung von Energiedichten und Intensitäten[END_REF] ), the total energy average is approximated by only computing the kinetic energy average:

E tot = 1 2 ( Ωx ρω 2 w n (x, ω, φ).w * n (x, ω, φ)dΩ/V ol ( 21
)
The group velocity vectors C g n (ω, φ) is computed for all wave numbers at each frequency. In order to focus our analysis on only flexural modes (S and SH ones) we introduce an indicator allowing to select them by computing the ratio of kinetic energy average on out of plane displacement as: (22) with wz n (x, ω, φ) being the (Oz) component of vector w n (x, ω, φ).

Ind(n, ω, φ) = 1 2 ( Ωx ρω 2 wz n (x, ω, φ)wz * n (x, ω, φ)dΩ/V ol E tot

Piezoelectric Material

OPTIMISATION OF THE FLEXURAL ENERGY FLOW INSIDE THE SHUNTED PERIODIC PIEZO-COMPOSITE

The considered piezo-composite cell is presented in figure 1. The supporting plate material is standard aluminum with 0.1 % of hysteretic damping ratio and the piezoelectric material characteristics is given in table 1.

The used methodology for optimizing the shunt impedance Z(iω) is based on the minimization of the maximal group velocity collinear to the wave number vector (19) for waves having a ratio of transported flexural kinetic energy (22) greater than 0.8. The used criterion can also be written as:

Crit(Z(iω), φ) = max n/Ind(n,ω,φ)>0.8 (C g n (ω, φ).Φ) (23)
The used numerical optimization of the criteria is based on a multidimensional unconstrained nonlinear minimization (Nelder-Mead).

Case 1 : Z(iω) ∈ C

In a first test, we optimize the criterion by considering any frequency dependent complex impedance. We present in figure 2(a) and 2(b) the obtained real parts and imaginary parts of the wave number kx n (iω) along (Ox) axis. The red circles mark the dispersions curves for Z = 0 and the blue crossed the optimal dispersion. The corresponding group velocities along (Ox) are presented in figure 3 while the real and imaginary parts of the optimal impedance are plotted in figure 4.

We immediately observe that the optimization of the shunt impedance leads to greatly modify the group velocity of the A o mode in the first part of the spectrum (i.e before the first bending band gap between 22 and 27 kHz). The bending waves also propagate energy with a lower velocity and can even been very low as evanescent waves before 6 kHz. The corresponding real and imaginary parts of the dispersion curves are slightly modified. We notice an increase of the imaginary values indicating an increase of the spatial decay rates. The optimal impedance values are almost real, and correspond to those obtained if a constant negative capacitance is used. The corresponding average value is -150.05 pC.V -1 . Some imaginary parts of the optimal impedance are negative which indicate that the optimization leads to provide energy to the system for controlling mechanical damping effect introduced with hysteretic damping ratios into the model, and, also, obtain a fully conservative system.

Case 2 : Z(iω) ∈ R

As the obtained negative impedance imaginary part does not appear realistic for physical implementation, the second case considers optimization by considering frequency dependent real impedance. We present in figure 5(a) and 5(b) the obtained real parts and imaginary parts of the wave number kx n (iω) along (Ox) axis. The red circles mark the dispersions curves for Z = 0 and the blue crossed the optimal dispersion. The corresponding group velocities along (Ox) is presented in figure 6 while the real and imaginary parts of the optimal impedance are plotted in figure 7. We immediately observe that the obtained results are globally similar to those provided by the complete complex optimization presented before, even if the equivalent constant negative capacitance is now -149.88 pC.V -1 . The numerical cost is lower and the optimization faster. In order to study the translation of the wave properties obtained by optimizing the previously described criterion on the dynamical response of a finite dimension system, we apply the optimal impedance on a finite set of shunted piezo-composite cells semi-distributed onto a plate system as described in figure 8. The harmonic response of this system is also computed at different frequency when optimal impedance is connected or not to each patches. We plot in figure 9 the obtained results in term of kinetic energy density.

Validation on a periodically semi-distributed set of adaptive cells

These numerical results clearly show up a strong influence of the modifications of the wave dispersion (i.e on the group velocities) in the standing wave responses plotted in figure 9. We observe at 5 kHz an increase in the dynamic response of the system when optimal shunt is connected. A large part of the system energy remains in its left part where the applied forces is located. The semi-distributed interface also change the system admittance and filter wave diffusion by increasing is reflexibility property. At 22 kHz, the energy diffusion is clearly condensed into the left part of the system with a largely decreased amplitude compared to this obtained with open circuit. The structural dynamical admittance has been decreased by connecting the shunt circuits. Finally, at 60 kHz, the energy is concentrated on the right part of the system between the adaptive interface and the clamped boundary condition. We observe something similar to a wave trap effect. The average value is, one more time, largely decrease by using the shunt circuits.

CONCLUSIONS

This paper presents a numerical procedure able to compute the damped wave's dispersion functions in the whole first Brillouin domain of multi dimensionnal piezo-elastodynamical wave guides. The method was applied for determining the optimal impedance allowing to minimize the group velocities of the flexural waves. Based on this approach, some numerical test on a finite dimension system incorporating a semi-distributed set of shunted piezo-composite cells has been performed. We underline a strong influence of the designed shunt circuits in the dynamical response of the system. Even if the link between the obtained wave properties are not clearly established, we also demonstrated that our developed numerical procedures can be used for optimizing the energy diffusion operator of such adaptive mechanical interface. To do so, additional work has to be done for optimizing the complete interface scattering and for controlling the evanescent waves playing an important role in the finite system dynamical response. The proposed methodology can also be used for studying particular dissipation phenomenon such as those induced by complex shunted piezoelectric patches as proposed by [START_REF] Beck | Broadband vibration suppression assessment of negative impedance shunts[END_REF] and, [START_REF] Casadei | Vibration control of plates featuring periodic arrays of hybrid shunted piezoelectric patches[END_REF] or even foams or complex polymers behaviors. The proposed method furnishes an efficient tool for future optimization of distributed smart cells as proposed in the case of 1D wave guide by.

Figure 1 .

 1 Figure 1. Generic 3D piezocomposite periodic cell

Figure 2 .Figure 3 .4

 23 Figure 2. Real parts (a) and imaginary parts (b) of the wave number kxn(iω) along (Ox), the red circles mark values obtained for Z = 0 and the blue crossed the optimal ones

Figure 4 .Figure 5 .Figure 6 .1

 456 Figure 4. Real and imaginary parts of the optimal impedance

Figure 7 .

 7 Figure 7. Real and imaginary parts of the optimal impedance

Figure 8 .

 8 Figure 8. 3D piezocomposite periodic semi-distributed Cells

6

 6

Figure 9 .

 9 Figure 9. Kinetic Energy density for different source frequency ((a) 5 kHz, (b) 22 kHz, (c) 60 kHz). On the left with optimal impedance connected to the piezo-patches and on the right with open circuit

Table 1 .

 1 Piezoelectric patch characteristics

	Symbol	Value	Property
	s E 11 = s E 22 = s33 E 11.6 × 10 -12 P a -1 s E 12 = s E 13 = s E -3.33 × 10 -12 P a -1 12, 13 and 23 compliance matrix coefficient 11, 22 and 33 compliance matrix coefficients 23 s E 44 = s E 55 = s E 66 45.0 × 10 -11 P a -1 44, 55 and 66 compliance matrix coefficients
	η	0.1 %	Hysteretic Damping ratio
	d 31 = d 32	-6 × 10 -11 C/N	31 and 32 piezoelectric matrix coefficients
	d 33	15.2 × 10 -11 C/N	33 piezoelectric matrix coefficient
	d 24 = d 15	730 × 10 -12 C/N	24 and 15 piezoelectric matrix coefficients
	ρ	7600kg/m 3	Density
	ε T 11 = ε T 22 ε T 33	504.1 ε o C/V/m 270 ε o C/V/m	Dielectric Permittivity Dielectric Permittivity

ACKNOWLEDGMENTS

This work was carried out with a grant of French agency ANR number N T 09-617542. We gratefully acknowledge the French ANR and CNRS for supporting this program.


Chaos and scale

As rightly emphasized by Ruxton, the relative timescales of reproduction and dispersal are critical for the dynamic of a metapopulation when local dynamics are taken into account. Two mechanisms of dispersal can be distinguished as extremes along a continuum -the steppingstone mechanism whereby the dispersers can reach only the neighbouring patches, and the 'common pool' of dispersers where the distance between patches is unimportant for the 1 dispersers. This general scheme is further complicated by habitat selection when habitat quality varies among patches, and by the triggering mechanism of dispersal. The 'common pool' mechanism can have a strong stabilizing influence both on population and metapopulation dynamics even when local dynamics are strongly chaotic 9,10 . This is achieved providing the dispersal rate is sufficiently low. When dispersal rate increases, it 'homogenizes' the whole system (as seasonality can do 11 ) so that chaos occurs. Stepping-stone dispersal contributes critically to metapopulation persistence 3,4 ; however, its stabilizing influence on global dynamics is yet to be demonstrated under realistic biological assumptions 12 . Both dispersal mechanisms could be important for the same species at different spatial or temporal scales.

Paradis

  Institut des Sciences de l'Evolution, CC 64, Université Montpellier II, Place Eugène Bataillon, F-34095 Montpellier cedex 5, France In a recent issue of TREE, Ruxton l points out Csilling et al.'s results 2 concerning the dynamic of a theoretical metapopulation. They observed that a set of linked chaotic populations exhibits neither chaotic temporal variations nor spatial structuring. In other studies, the opposite result was observed 3-6 . Ruxton argued that this contrast could be explained by the fact that the timescales of reproduction and dispersal are decoupled in Csilling et al.'s model. I think the absence of chaos in Csilling et al.'s study could be explained by a particular feature of their model. The modelled metapopulation occupies a landscape composed of a finite number of patches. All local populations have the same chaotic dynamic. Dispersal occurs when the population density exceeds a fixed threshold according to a stepping-stone rule. A fixed fraction of the population emigrates to neighbouring patches. The dispersal process is reiterated until all populations are below the threshold, then reproduction starts again. Since a given population may receive dispersers from other patches, numerous dispersal steps may occur in what Csilling et al. call a 'migration avalanche'. The individuals dispersing out of the landscape are lost, so that all excess individuals are eliminated. Such a dispersal process is an efficient density-dependent regulating mechanism. It should be interpreted as an inhibition of reproduction in local populations as long as at least one population density exceeds the threshold. Such a mechanism does not seem biologically realistic, since it assumes that reproduction in a population is not only directly influenced by its own density but also by densities in the other populations. May and his collaborators 7,8 pointed out a long time ago that chaos occurs when there is a discrepancy between the timescales of reproduction and its density-dependent inhibition. In Csilling et al.'s model, the local density-dependent mechanisms are outweighed by the dispersal mechanism.


Clémentine Cottineau
email: c.cottineau@ucl.ac.uk

Juste Raimbault

Marion Le Texier

Florent Le Néchet

Romain Reuillon

Initial spatial conditions in simulation models: the missing leg of sensitivity analyses?

Keywords: Space, Initial conditions, Sensitivity, ABM

Although simulation models of geographical systems in general and agent-based models in particular represent a fantastic opportunity to explore socio-spatial behaviours and to test a variety of scenarios for public policy, the validity of generative models is uncertain until their results are proven robust. Sensitivity analysis usually include the analysis of the effect of stochasticity on the variability of results, as well as the effects of small parameter changes. However, initial spatial conditions are usually taken for granted in geographical models, thus leaving completely unexplored the effect of spatial arrangements on the interaction of agents and of their interactions with the environment. In this contribution, we present a method to assess the effect of initial spatial conditions on simulation models, using a systematic generator controlled by metaparameter to create density grids used in spatial simulation models. We show, with the example of two very classical agent-based models (Schelling's models of segregation and Sugarscape) that the effect of space in simulation is significant, and sometimes even larger than parameters themselves. We do so using high performance computing in a very simple and straightforward open-source workflow.

Introduction

Simulation has been recognised and increasingly used by geographers to explore various geographical processes and problems within virtual laboratories. It appears as a very fruitful way to overcome the difficulty of analytic resolution of many spatial models developed in the past, and to explore the possible trajectories of social and ecological systems in space and time. The specificity of geographical models compared with models from other social sciences is usually regarded as the way geographers consider space and spatial interactions, driven by an explicit interest in the way space influences the outcomes of the model. Geographers are indeed concerned about understanding and modelling how space plays a role in social interactions and environmental processes, and whether its action is placed-based or place-neutral. think simulation can become a very good tool to achieve this, provided that models include relevant spatial description and modelling, behavioural rules that take space into account, and provided model evaluation stresses the sensitivity of output variations to the way space is modelled. This communication focuses on a way to do so: by using a sensitivity method to test the model's outcomes to initial spatial conditions.

2 Why spatial patterns are expected to impact social simulations 2.1 Because actual cities are not regular grids Although this may seem obvious, cities are not regular grids of isotropic densities. Modelling may be about abstracting features to highlight processes in a way which is simpler than reality. However, the uniform grid which represents space in most simulation models is potentially not enough to represent urban processes because density and accessibility have environmental, economic and social consequences. An intermediate and more meaningful way of abstracting space might thus be to consider, not the peculiarities of every city, but their broad density structures. In Europe for example, one can find broad types of density distributions (Le [START_REF] Néchet | De la forme urbaine à la structure métropolitaine: une typologie de la configuration interne des densités pour les principales méropoles européennes de l'Audit Urbain[END_REF]). At initialisation, we then expect agents to be distributed in different patterns, which influence results in the long run. A first approach could be to evaluate the effects of the same set of social mechanisms, but modelled in different types of cities, for example monocentric, polycentric and discontinuous.

Because actual agents are not random walkers

We expect spatial patterns to influence model outcomes also because the agents' rule of action itself might depend on the spatial structure of the environment. Indeed, mechanisms of surrounding sensing will be impacted by different distributions of density. For example, agents tend to create buffer zones if they are allowed network-based rather than isotropic movements [START_REF] Banos | Network effects in Schelling's model of segregation: new evidences from agent-based simulation[END_REF]), while the vision / sensing mechanism is sensitive to the scale of modelled environments [START_REF] Laurie | Role of'vision'in neighbourhood racial segregation: a variant of the Schelling segregation model[END_REF]; [START_REF] Fossett | Effects of city size, shape, and form, and neighborhood size and shape in agent-based models of residential segregation: are Schelling-style preference effects robust?[END_REF]). Finally, households can have different preferences with respect to the built-environment [START_REF] Spielman | The Co-evolution of Residential Segregation and the Built Environment at the Turn of the 20th Century: A Schelling Model[END_REF]), thus creating a sensitivity to the initial form of the city modelled.

Method

In this paper, we suggest tackling sensitivity to spatial conditions by generating a variety of density grids to feed into simulation models at initialisation, and exploring the sensitivity of the model outcomes to the variation of meta-parameters (i.e. parameters used to generate initial spatial conditions). The purpose is two-fold: 1/ to test the robustness of simulation results to small variations of meta-parameters within a typical category of space (a monocentric case for example) and 2/ to study the non-trivial effects of typical categories of spatial distribution (monocentric vs. polycentric for example) on the results of a given model (cf. figure 1).

1: General workflow of our method

Constructing a spatial generator

Our spatial generator applies an urban morphogenesis model [START_REF] Batty | Cities and complexity: understanding cities with cellular automata, agent-based models, and fractals[END_REF]), which has been generalised, explored and calibrated [START_REF] Raimbault | A hybrid network/grid model of urban morphogenesis and optimization[END_REF]). In a nutshell, grids are generated through an iterative process which adds a quantity N (population) at each time step, allocating it through preferential attachment characterised by its strength of attraction α. This first growth process is then smoothed n times using a diffusion process of strength β. Grids are thus generated from the combination of the values of these four meta-parameters α, β, n and N , in addition to the random seed. To ease our exploration, only the distribution of density is allowed to vary rather than the size of the grid, which we fix to a 50x50 square environment of 100,000 units (cf. figure 2).

Including a (typical) variety of initial conditions in the sensitivity analysis

In order to generate density grids which correspond to empirical density distributions, we select among the generated grids using an objective function which matches the point cloud of 110 metropolitan areas in Europe described by four dimensions: their concentration index, hierarchy index, centrality index and continuity index (cf. Le [START_REF] Néchet | De la forme urbaine à la structure métropolitaine: une typologie de la configuration interne des densités pour les principales méropoles européennes de l'Audit Urbain[END_REF]). A stochastic exploration of a Latin Hypercube Sampling of 2000 points in the 4-dimensional space of parameters α, β, n, N gives a subset of 170 interesting grids matching empirical densities, which constituted our set of different initial spatial conditions. These are further clustered into three classes of morphology: compact (e.g. Vienna), polycentric (Liege) and discontinuous (Augsburg) in order to evaluate the non-trivial effects of urban form on simulation results. We select 15 grids of each type to work with.

Figure 2 :

 2 Figure 2: Four examples of grids produced by the spatial generator

Application cases

We use our of initial density grids in two spatial simulation models: Schelling and Sugarscape.

Schelling's model of residential segregation

Schelling's model consists in a abstract urban housing market where agents of different nature sense their environment, evaluate their satisfaction in terms of neighbourhood composition, and relocate if unsatisfied. It has been shown by [START_REF] Aziz-Alaoui | Models of segregation[END_REF] that even tolerant agents tended to produce segregated patterns because of the complexity of their local interactions. The main parameters of this model are the tolerance level (% of agents similar to ego), the scope of sensing and the percentage of vacant spaces in the housing market. In addition, we are interested in testing the impact of the spatial distribution of housing capacity in this project, using the generated grids.

Sugarscape model of resource extraction and population settlement

Sugarscape is a model of resource extraction which simulates the unequal distribution of wealth within a heterogenous population [START_REF] Epstein | Growing artificial societies: Social science from the bottom up (complex adaptive systems)[END_REF]). Agents of different vision scopes and different metabolisms harvest a self-regenerating resource available heterogeneously in the initial landscape, they settle and collect this resource, which leads some of them to survive and others to perish. The main parameters of this model are the number of agents, their minimal and maximal resource. In addition, we are interested in testing the impact of the spatial distribution of the resource in this project, using the generated grids. The outcome of the model is measured as a phase diagram of an index of inequality for ressource distribution (Gini index).

Results

We proceed to 4,500,000 simulation runs of the Schelling model (1000 parameter combinations x 45 density grids x 100 replications), using OpenMOLE to distribute the computation, and apply segregation measures to characterise the results. We find that the different urban morphologies impact the parameter interaction patterns, and that polycentric and discontinuous cities appear systematically more segregated than compact cities in terms of dissimilarity and entropy index. For Sugarscape, 2,500,000 simulations (1000 parameter points x 50 density grids x 50 replications) allow us to show that the model is more sensitive to space than to its other parameters, both qualitatively and quantitatively: the amplitude of variations across density grids is larger than the amplitude in each phase diagram, and the behavior of phase diagram is qualitatively different in different regions of the morphological space.

Acknowledgements

The authors acknowledge the funding of their institutions and the EPSRC project number EP/M023583/1.


Jean-Marc Meunier
email: jmeunier@univ-paris8.fr

Samuel Szoniecky

AN ONTOLOGY-BASED ENRICHED EBOOK FOR TEACHING STATISTICS

With the development of Internet and the use of digital pedagogy, the production of open educational resources grows very quickly. Sharing these resources has become a major challenge both from an economic point of view, but also for the dissemination and sharing of practices among teachers. The problematics are psychologic (help for learning), ergonomic (usability, acceptability) but also documentary (reuse resources) and editorial (enhancement of digital archives). Based on an ontology developed by University of Oxford and translated in French, we developed a digital archive of open educational resources for teaching statistics. We present and discuss e-book enrichment scenarios in undergraduate statistics courses. We present how this ontology-based approach allows thinking enriched e-books as an adaptive and collaborative learning device.

INTRODUCTION

Teaching statistics requires the coordination of a functioning logic (the theoretical basis of statistics) and a use logic (what test to choose and how to implement it concretely). The first presupposes the acquisition of a very abstract knowledge, the second presupposes the very fine particularization of this knowledge [START_REF] Richard | Logique du fonctionnement et logique de l'utilisation[END_REF]. This coordination is based on the structure of knowledge, in particular the organization of concepts and their properties. These properties can be structural or functional (procedural-related). For example, knowing which properties a sample must have in order to be able to apply a parametric test is the first logical aspect. Knowing how to structure a data table according to the software used, or how to adapt the procedure depending on whether you are working on individual data or a distribution, is a matter of the second. This articulation of knowledge and know-how through the structure of knowledge offers a formalization of the notion of competencies [START_REF] Paquette | An Ontology and a Software Framework for Competency Modeling and Management[END_REF] and at the same time offers ways to develop a logical chain of learning. Above all, it makes it possible to design new types of teaching materials that make it possible to break away from strict linearity of reading and thus offer its adaptability to the diversity of learners and situations in which teaching resources are consulted. It is this possibility that we propose to explore in this article.

FORMALIZE THE KNOWLEDGE STRUCTURE TO INSTRUMENTALIZE IT

This knowledge structure can be formalized in a domain ontology. It is the description of an area of knowledge using a set of hierarchical and interrelated concepts. It is based on a class taxonomy and a set of inference rules. Within the context of a project called Ontostats (http://www.ontostats.univ-paris8.fr/), supported by the French Ministry of Higher Education, Research and Innovation, we translated into French the ontology of STATO statistical concepts developed by the University of Oxford [START_REF] Gonzalez | Stato -Statistics Ontology[END_REF].

The interest of such an ontology is to allow the sharing of knowledge, not only about concepts but also about the relationships between them. Ontology is a vocabulary used to label concepts and the relationships between them. The separation of labels and concepts allows for a multilingual description of the domain of knowledge. Relationships are class inclusion relationships, but also properties or attributes that can take multiple values. They allow us to deduce what is usually implicit in the discourse and contribute to the development of decision support systems. Ontologies also have many advantages for indexing educational resources compared to using a norm [START_REF] Bouzeghoub | Ontology Mapping for Learning Objects Repositories Interoperability[END_REF]. However, ontologies are promising especially for building learning aids. They provide a conceptual model for the study of learning disabilities. An ontology can be used to design training paths by adapting the structure of teaching sequences [START_REF] Guo | An Ontology Infrastructure for an E-Learning Scenario[END_REF], to design resource recommendation systems based on relationships between concepts [START_REF] Shen | Ontology-based learning content recommendation (English)[END_REF], and also when ontology is associated with learning traces analysis, course personalization [START_REF] Castro | Learning Objects and Ontologies to Perform Educational Data Mining[END_REF] and learner evaluation [START_REF] Romero | Pedagogically-Driven Ontology Network for Conceptualizing the e-Learning Assessment Domain[END_REF].

The STATO ontology has been used to index a set of open educational resources, notably those developed by French thematic digital universities. These resources have been collected in a digital archive using Omeka-S. In addition, the archive offers each identified teacher the opportunity to work on his or her own collection of resources and to propose a personalized editorialization. It is on this last feature that our enriched digital book project is based.

TOWARDS A NEW CONCEPTION OF A DIGITAL TEXTBOOK

The book is a medium through which a teacher can express and structure knowledge at some time in order to enable a learner to acquire it. The appropriation of this knowledge through the reading of a textbook depends, of course, on the objectives of reading, the level of knowledge, the activities motivating reading and the knowledge related to the use of the reading device itself (use of glossary, search engine, navigation function). The teacher, when developing the pedagogical design and editorializing the educational resources, makes choices and places them in an original production that takes into account the context, domain and audience. [START_REF] Armatte | Le rôle de l'histoire dans l'enseignement de la statistique[END_REF] distinguishes between four approaches to the teaching of statistics: (i) a formal, hypothetical-deductive approach (focused on the teaching of axioms); (ii) an inductive or experimental approach, centred on simulations and data mining; (iii) a problem-oriented approach based on the logic of the social field or the problems of daily life; (iv) the historical approach consisting in presenting to students the consequences of problematization that led to the development of a particular method. We can add a fifth, procedural one, which consists in focusing on the when and how to do, therefore very dependent on the analysis tool, particularly software, used during teaching. The expertise coordinates its five approaches in a more or less important way, but it is rare to have the time, within the same course, to present them all. Depending on the discipline and the role that statistics play in it, a teacher may prefer one to another. In the same way, a student, depending on his or her previous background, may be more or less sensitive to one or the other. Giving students a bridge between these different points of view is one of the challenges and not the least of them to allow access to these different points of view.

However, the current digital versions of the statistics textbooks because, for the vast majority of them, they are transpositions of the paper version offer very few advantages. For many actors in the field, e-books are nothing more than "the combination of a classical book structure, or rather the familiar concept of books, with characteristics related to the electronic environment" [START_REF] Landoni | Electronic Books[END_REF]. It does not offer the possibility of varying points of view and linking them to a more general body of knowledge, so that the book is enclosed within itself. To overcome this enclosure, it is necessary to imagine a new form of writing which according to [START_REF] Clément | Le e-book est-il le futur du livre ? In Les savoirs déroutés[END_REF] is characterized by three essential aspects: hypertextuality, its distribution on the network and finally its dynamic and multimedia character.

The first characteristic makes it possible to overcome the linearity of the discourse. If this is a valuable help for beginners, it can be a brake on learning as skills increase. It can even be blocking when the book is used in problem solving or synthesis activities on a given theme [START_REF] Roy | Statistiques mieux « dévorées » à l'aide d'un ebook ?[END_REF]. With hypertexts, it becomes possible to browse a document in a non-linear way. However, we must be careful. Many studies in cognitive psychology show that a badly controlled hypertextual structure can be a major barrier to learning, but well used, an opportunity [START_REF] Vörös | Effect of high-level content organizers on hypertext learning[END_REF]. We propose to rely on an ontology to allow the transition from a thematic structure, the table of contents, to a structure oriented towards the relationship between concepts and procedures, to see them at different levels while maintaining an overall coherence. Ontology becomes the guide to navigating and working on course content.

Figure 2 Principles of Enrichment of an e-book from an ontology

The second characteristic is distribution on the network. This is a reality that all teachers now know. A course can no longer be limited to itself. Students as needed complement it or compare it with other resources on the Internet. We may regret it and develop defensive strategies, or take advantage of the opportunity for diversity to enrich our course or offer students a personalized learning experience. This involves rethinking the way a course is written to take full advantage of computers by inserting the text into a larger set of documents that, if correctly indexed, make it possible to associate with a book additional or alternative open source and editorial resources that correspond to the teaching objective. The concept of document discussed here is very general as it concerns texts or multimedia documents, exercises, datasets, notes or synthesis of students' lectures, as well as forum or frequently asked questions or information stored in the ontology itself (definitions, multilingual lexicon and relations with other concepts). From this point of view, the digital manual is not just one of several documents on a learning management system (LMS), it is the LMS itself. Such a conception of the textbook makes it a dynamic object that can be adapted according to the needs of the learners, or even adapted to each of them to offer them a personalized reading and learning experience. Figure 2 summarizes the relationships between ontology, resource repository and manual development.

CONCLUSION AND FUTURE WORK

The platform will be in production in the autumn of 2018, at the same time as the realization of a first prototype manual. At the same time, we have initiated two studies. The results of the first one on the use of open and editorial educational resources and electronic books should allow us to refine the specifications for such an electronic book. The second one aims to identify the profiles and needs of the learners in order to personalize the learning paths. This includes identifying a minimum number of questions to assess the learner's state in terms of expert knowledge and providing him/her with the most relevant resources and activities to help him/her progress.

Figure 1

 1 Figure 1 An extract from the ontology of statistical concepts (Meunier, Szoniecky, & Lamolle, 2017)


I. INTRODUCTION

The main goal of this project was to build a basic musculo-skeletal model of the human body in order to be able of analyzing and simulating human motion. To this end, an OpenSim [START_REF] Delp | Opensim: Open-source software to create and analyze dynamic simulations of movement[END_REF] musculo-skeletal model was modified and augmented based on previous works of the literature. The resulting model is also compatible with Pinocchio [START_REF] Carpentier | Pinocchio: fast forward and inverse dynamics for poly-articulated systems[END_REF], an open source C++ library for efficient rigid multi-body dynamics computations based on revisited Roy Featherstone's algorithms. This allows users for creating powerful simulations of whole-body human-based motions with anthropomorphic systems. Using Pinocchio with the model permits also to use state-of-the-art libraries written in C++ or python (e.g. from robotics, machine learning, ...).

The first version of the whole-body model was modified and built based on the following specific criteria:

• The model should express the main elementary motions of humans performing parkour [START_REF] Maldonado | Angular momentum regulation strategies for highly dynamic landing in parkour[END_REF], [START_REF] Maldonado | Strategies of parkour practitioners for executing soft precision landings[END_REF].

• The model has to be used for the analysis and generation of anthropomorphic motion [START_REF] Maldonado | From biomechanics to robotics[END_REF]. This implies that the model should be more complicated than a humanoid robot model for instance, but not as a complex as a complex human model. • The model should be whole-body and three dimensional as sportive practicers require whole-body coordination for accomplishing complex motions.

• The model has to include surface joint stabilizer muscles that favor the study of muscular activity in highly dynamic motion. For example the agonist-antagonist muscles that perform abduction-adduction, flexion-extension or plantar flexiondorsiflexion. In the sequel we provide a technical report of the skeletal and muscular modeling. This model is freely available for download at https://github.com/GaloMALDONADO/ospi along with all the libraries to parse the model to use it with Pinocchio. Some examples of utilization and motion data are also provided.

II. METHODS

In this section we aim at providing the description of the anthropomorphic model, of the joints modeling, angle conventions and ranges of motions.

A. Skeletal model

The model includes 42 degrees of freedom (Table I), 18 joints and 19 segments. The pose (position and orientation) of a reference frame attached to the pelvis is expressed with respect an intertial world frame (Fig 1). The skeleton is modeled as a friction-less articulated multi-body dynamic system as shown in Fig 2 . The pelvis, thighs, shanks and feet segments mass and inertia were set according to the model of [START_REF] Frank | A Dynamic Optimization Solution for Vertical Jumping in Three Dimensions[END_REF] where data was obtained based on the studies of [START_REF] Mcconville | Anthropometric relationships of body and body segment moments of inertia[END_REF]. In this model, the length of the segments was taken from [START_REF] Delp | An interactive graphics-based model of the lower extremity to study orthopaedic surgical procedures[END_REF] model. The same model was used in [START_REF] Samuel | Muscle contributions to propulsion and support during running[END_REF] for simulating running motions. Upper limb anthropometry are based on the running model of [START_REF] Samuel | Muscle contributions to propulsion and support during running[END_REF] which was estimated based on the regression equations of [START_REF] De | Adjustements to Zatsiorsky-Seluyanov's segment inertia paremeters[END_REF] and [START_REF] Holzbaur | A Model of the Upper Extremity for Simulating Musculoskeletal Surgery and Analyzing Neuromuscular Control[END_REF]. Antrhopometry of the torso and head segments (including the neck) were estimated from the regression equations of [START_REF] Dumas | Adjustments to mcconville et al. and young et al. body segment inertial parameters[END_REF], [START_REF] Dumas | Corrigendum to "Adjustments to McConville et al. and Young et al. body segment inertial parameters[END_REF].

B. Joints modeling

The characteristics of the joint modeling are presented in Table I. The joint modeling adapted from the literature is described below:

• Each lower extremity has 7 DoF. The hip is modeled as a ball-and-socket joint, the knee is modeled as a hinge joint, the ankle is modeled as 2D hinge joints (flexion-extension and inversion-eversion), and the toes are modeled with one hinge joint at the metatarsals. Lower extremity joint definitions and the pelvis ones are based on [START_REF] Delp | An interactive graphics-based model of the lower extremity to study orthopaedic surgical procedures[END_REF] corrected in [START_REF] Arnold | A model of the lower limb for analysis of human movement[END_REF]. • The pelvis joint is modeled as a free-flyer joint to permit the model to translate and rotate in the 3D space. This 6D joint is attached to the free-floating base (root frame) of the under-actuated systems. The lumbar motion is modeled as a ball-and-socket joint and the neck joint is also modeled as a ball-and-socket joint. Joint definitions are based on [START_REF] Frank | A Dynamic Optimization Solution for Vertical Jumping in Three Dimensions[END_REF]. • Each arm is modeled with 8 DoF. The shoulder is modeled as a ball-and-socket joint, the elbow and forearm rotations are modeled with hinge joints to represent flexion-extension and pronation-supination [START_REF] Holzbaur | A Model of the Upper Extremity for Simulating Musculoskeletal Surgery and Analyzing Neuromuscular Control[END_REF], the wrist flexion-extension and radial-ulnar deviations are modeled with hinge joints, and the hand fingers are modeled with one hinge joint for all fingers. Joint definitions are based on [START_REF] Samuel | Muscle contributions to propulsion and support during running[END_REF]. For further details about the location of each coordinate system with respect to bony landmarks, see [START_REF] Frank | A Dynamic Optimization Solution for Vertical Jumping in Three Dimensions[END_REF]. In the sequel we present some illustrations showing the ranges of motion of the model and the respective angle conventions.

III. MARKER SET

The model includes a whole-body marker set with 48 markers placed on anatomical bony landmarks (13) selected based on the recommendations of the International Society of Biomechanics (ISB) [START_REF] Wu | ISB recommendation on definitions of joint coordinate system of various joints for the reporting of human joint motion-Part I: ankle, hip, and spine[END_REF], [START_REF] Wu | ISB recommendation on definitions of joint coordinate systems of various joints for the reporting of human joint motion-Part II: shoulder, elbow, wrist and hand[END_REF]. The musculoskeletal model (Fig. 14) has been assembled to perform whole-body analysis of 3d human motion. The model contains 48 superficial muscles (Fig. 14) that represent the more relevant joint stabilizer muscles to study sport movements. The Muscles are modeled using hill-type models based on the works described in [START_REF] Schutte | Improving the efficacy of electrical stimulation-induced leg cycle ergometry: an analysis based on a dynamic musculoskeletal model[END_REF], [START_REF] Darryl G Thelen | Adjustment of muscle mechanics model parameters to simulate dynamic contractions in older adults[END_REF] which are used in OpenSim [START_REF] Delp | Opensim: Open-source software to create and analyze dynamic simulations of movement[END_REF]. The muscles and their principal mechanical actions in the execution of sport movements, are presented in Table III for the lower limb muscles, in Table IV for the upper limb muscles and in Table V for the trunk and head muscles.

A. Lower limbs muscles

Twenty three superficial muscles that stabilize the hip, knee and ankle joints were selected based on a functional analysis of the muscles mechanical action in the execution of jumps and landing techniques (Table III). Four superficial muscles that stabilize the lumbar and neck joints were selected based on a functional analysis of the muscles mechanical action in the execution of the sport techniques V. The sternocledomastoid and the trapecius ascendus muscles were used for studying the motion of the head. Nevertheless, the trapecius ascendus muscle is also important for the scapulae movement in the execution of climbing techniques.

Fig. 1 :

 1 Fig. 1: Kinematics chain of the model. In order to model hinge joints (2 DoF), extra segments were added (talus and calcaneus, ulna and radius). The visual elements are based on the running model of Hammer [8].

Fig. 2 :

 2 Fig. 2: Kinematics chain of the model. The numbers inside the colored balls indicate the number of degrees of freedom between segments. In order to model hinge joints (2 DoF) with Pinocchio, extra segments were added (talus and calcaneus, ulna and radius) to the model.

  Fig 2 show how joints are linked from parent to child segments in the kinematic chain. Euler XYZ body-fixed rotation angles are used to express the orientation using OpenSim.

Fig. 3 :

 3 Fig. 3: Shoulder joint conventions and range of motions. (a) Maximum shoulder extension (-) equal to 90 deg, anatomical pose and maximum shoulder flexion (+) equal to 180 deg. (b) Maximum shoulder abduction (-) equal to 180 deg, anatomical pose and maximum shoulder aduction (+) equal to 35 deg. (c) Maximum shoulder internal rotation (+) equal to 100 deg, anatomical pose and maximum shoulder external rotation (-) equal to 100 deg.

Fig. 4 :

 4 Fig. 4: Elbow joint conventions and range of motions. (a) Maximum forearm supination equal to 110 deg, anatomical pose and maximum forearm pronation equal to 110 deg. (b) Maximum elbow extension which is the same as the anatomical pose and maximum elbow flexion equal to 170 deg.

Fig. 5 :Fig. 6 :

 56 Fig. 5: Wrist joint conventions and range of motions. (a) Maximum radial deviation equal to 25 deg, anatomical pose and maximum ulnar deviation equal to 35 deg. (b) Maximum wrist extension equal to 90 deg, anatomical pose and maximum wrist flexion equal to 90 deg.

Fig. 7 :

 7 Fig. 7: Neck joint conventions and range of motions. (a) Maximum neck extension (-) equal to 60 deg, the anatomical pose and maximum neck flexion (-) equal to 60 deg. (b) Maximum neck right bending (+) equal to 45 deg, the anatomical pose and maximum neck left bending (-) equal to 45 deg. (c) Maximum neck right rotation (-) equal to 80 deg, the anatomical pose and maximum neck left rotation (+) equal to 80 deg.

Fig. 8 :

 8 Fig. 8: Lumbar joint conventions and range of motions. (a) Maximum lumbar extension (+) equal to 30 deg, anatomical pose and maximum lumbar flexion (-) equal to 90 deg. (b) Maximum lumbar right bending (+) equal to 30 deg, anatomical pose and maximum lumbar left bending (-) equal to 30 deg. (c) Maximum lumbar right rotation (-) equal to 30 deg, anatomical pose and maximum lumbar left rotation (+) equal to 30 deg.

Fig. 9 :

 9 Fig. 9: Hip joint conventions and range of motions. (a) Maximum hip extension (-) equal to 30 deg, anatomical pose and maximum hip flexion (+) equal to 120 deg. (b) Maximum hip abduction (-) equal to 50 deg, the anatomical pose and maximum hip adduction (+) equal to 30 deg. (c) Maximum hip external rotation (-) equal to 45 deg, the anatomical pose and maximum hip internal rotation (+) equal to 35 deg.

Fig. 10 :

 10 Fig. 10: Knee joint conventions and range of motions. (a) Maximum knee flexion (-) equal to 140 deg, anatomical pose and maximum knee extension (+) equal to 10 deg.

Fig. 11 :

 11 Fig. 11: Ankle joint conventions and range of motions. (a) Maximum ankle dorsiflexion (+) equal to 20 deg, anatomical pose, maximum ankle plantarflexion (-) equal to 50 deg. (b) Maximum ankle eversion (-) equal to 50 deg, anatomical pose and maximum ankle inversion (+) equal to 30 deg.

Fig. 12 :

 12 Fig. 12: Toes joint conventions and range of motions. (a) Maximum toes extension (+) equal to 90 deg, anatomical pose and maximum toes flexion (-) equal to 45 deg.

Fig. 13 :

 13 Fig. 13: Marker set depicted on the whole body skeletal model

Fig. 14 :

 14 Fig. 14: Musculoskeletal model. The visual elements are based on the running model of Hammer [8].

  medius Abduct an rotate medially the thigh at the hip. (2) Sartorius Flex, abduct, and rotate laterally the thigh at the hip. Flex the leg at the knee. (3) Tensor fasciae latae Extend the thigh at the hip. (4) Biceps femoris Flex and rotate the leg laterally at the knee. Extend the thigh at the hip. (5) Semitendinoseus Flex and rotate the leg medially at the knee. Extend the thigh at the hip. (6) Vastus lateralis Extend the leg at the knee. (7) Vastus medialis Extend the leg at the knee. (8) Lateralis gastrocnemius Plantar flex the foot at the ankle. (9) Medialis gastrocnemius Plantar flex the foot at the ankle and flex the leg at the knee. (10) Soleus Plantar flex the foot at the ankle. (11) Tibialis posterior Plantar flex the foot at the ankle. Slight inversion of the foot at the ankle. (12) Tibialis anterior Dorsiflex the foot at the ankle. Slight inversion of the foot at the ankle.

  Erector spinae longissimus Extension of the lumbar spine at the lumbar joint (25) Rectus abdominis Flexion of the lumbar spine at the lumbar joint (26) Sternocledomastoid Oblique rotation and flexion of the head at the neck joint (27) Trapecius ascendus Scapulae motion and extension of the head at the neck joint

TABLE I :

 I Synovial joints and free flyer joint of the whole body model.

TABLE II :

 II Marker set IV. MUSCLES MODELING

TABLE III :

 III Lower body muscles of the whole body model and their mechanical actionB. Upper limbs musclesTwenty one superficial muscles that stabilize the shoulder, elbow and wrist joints were selected based on a functional analysis of the muscles mechanical action in the execution of climbing techniques (TableIV). Triceps brachii lateralis Extension of the forearm at the elbow joint.(20) Biceps brachii longus Supination and flexion of the forearm at the elbow joint. (21) Brachioradialis Flexion of the forearm at the elbow joint (22) Carpi radialis flexor Flexion and radial abduction of the hand at the wrist joint (23) Carpialis extensor Extension and radial abduction of the hand at the wrist joint.

		Upper Limbs
	Muscle	Main Functions
	(13) Anterior deltoideus	Abduction, flexion and medial rotation of the forearm at the shoulder joint
	(14) Posterior deltoideus	Abduction, extension and lateral rotation of the upper arm at the shoulder joint
	(15) Medial deltoideus	Abduction of the upper arm at the shoulder joint
	(16) Pectoralis major	Flexion, adduction and medial rotation of the upper arm at the shoulder joint
	(17) Latisimus dorsi	Extension, adduction, horizontal abduction, flexion from an extended position, and medial rotation
		of the upper arm at the shoulder joint
	(18) Triceps brachii longus	Extension of the forearm at the elbow joint. Also adducts and may assist in extension of the
		shoulder joint.
	(19)

TABLE IV :

 IV Upper body muscles of the whole body model and their mechanical action C. Trunk and head muscles

TABLE V :

 V Trunk muscles of the whole body model and their mechanical action


Martial Sauceau
email: martial.sauceau@enstimac.fr

Clémence Nikitine

Élisabeth Rodier

Jacques Fages

E Rodier

Effect of supercritical carbon dioxide on polystyrene extrusion

Keywords: Polystyrene, Polymer, Foam, Extrusion, Supercritical carbon dioxide

A study on the extrusion of polystyrene was carried out using supercritical carbon dioxide (scCO 2 ) as foaming agent. scCO 2 modifies the rheological properties of the material in the barrel of the extruder and acts as a blowing agent during the relaxation at the passage through the die.

For experiments, a single-screw extruder was modified to be able to inject scCO 2 within the extruded material. The effect of operating parameters on material porosity was studied. Samples were characterized by using water-pycnometry, mercury-porosimetry and scanning electron microscopy. Polystyrene with expansion rate about 15-25% was manufactured. A rapid cooling just downstream the die is important to solidify the structure. The die temperature allows the control of the porosity structure. CO 2 concentration shows no significant influence.

Introduction

Extrusion is a process of converting a raw material into a product of uniform shape and density by forcing it through a die under controlled conditions [START_REF] Rodier | Polymer Extrusion[END_REF]. Industrial applications of the extrusion process dates back to the 1930s. It is extensively applied in the plastic and rubber industries, where it is one of the most important manufacturing processes. A particular application concerns the built up of polymeric foams. Conventional foamed products are produced using either chemical or physical blowing agents. Various chemical blowing agents, which are generally low molecular weight organic compounds, are mixed into a polymer matrix and decompose when heated beyond a threshold temperature. This results in the release of a gas, and thus the apparition of bubbles. This implies however the presence of residues in material and the need for an additional stage to eliminate them. Techniques using physical agents include the mixing of a gas with melted polymer. Then, a thermodynamic instability is induced through a pressure or temperature change to nucleate bubbles of the blowing agent. After nucleation, the bubble growth is controlled by reducing the temperature below the glass transition temperature of the polymer.

A supercritical fluid (SCF) is defined as a substance for which both pressure and temperature are above the critical value. In general, SCFs offer mass transfer advantages over conventional organic solvents because of their gas-like diffusivity, liquid-like density, low viscosity and surface tension [START_REF] Eckert | Supercritical fluids as solvents for chemical and materials processing[END_REF]. In particular, supercritical carbon dioxide (scCO 2 ) has emerged as an important SCF due to its many desirable attributes such as low cost, non-toxicity, non-flammability, chemical inertia and easily accessible supercritical conditions (T c = 31.1 • C, P c = 7.38 MPa). When CO 2 is raised above its critical point, its physicochemical properties can be continuously tuned between vapour-like and liquid-like limits by varying the system pressure and/or temperature. Moreover, scCO 2 has been widely used in polymer processing [START_REF] Kazarian | Polymer processing with supercritical fluids[END_REF][START_REF] Nalawade | Supercritical carbon dioxide as a green solvent for processing polymer melts: processing aspects and applications[END_REF] as it is possible to dissolve large quantities of CO 2 in a variety of polymers [START_REF] Aubert | Solubility of carbon dioxide in polymers by the quartz crystal microbalance technique[END_REF][START_REF] Sato | Solubilities and diffusion coefficients of carbon dioxide and nitrogen in polypropylene, high-density polyethylene, and polystyrene under high pressures and temperatures[END_REF]. This involves in particular the swelling and the plasticization of the polymer, inducing changes in mechanical and physical properties [START_REF] Kazarian | Polymer processing with supercritical fluids[END_REF]. For instance, scCO 2 decreases the glass transition temperature, the interfacial tension or the viscosity of various polymers, without changing the viscoelastic behaviour [START_REF] Gourgouillon | Simultaneous viscosity and density measurement of supercritical CO 2 -saturated PEG 400[END_REF]. With polystyrene (PS), the glass transition temperature decreases from 105 • C to 98 • C and to 46.4 • C with a CO 2 mass fraction of 1 and 5.9%, respectively [START_REF] Condo | Glass transition behavior including retrograde vitrification of polymers with compressed fluid diluents[END_REF] and the interfacial tension decreases linearly from 28 to 17 mJ m -2 on a CO 2 pressure range from 1 to 10 MPa at 200 • C [START_REF] Li | Effect of carbon dioxide on the interfacial tension of polymer melts[END_REF][START_REF] Park | Surface tension measurement of polystyrene melts in supercritical carbon dioxide[END_REF].

The presence of scCO 2 in extrusion process modifies rheological properties of the polymer in the barrel of the extruder and acts as a blowing agent during the relaxation at the passage through the die. Thus, its high solubilization in the polymer results in extensive expansion at the die. The reduction of viscosity decreases the mechanical constraints and the operating temperature within the extruder.

In classical theory, the homogeneous nucleation rate (i.e. the number of cells created per units of time and volume) can be expressed as a function of the Gibbs free energy barrier for the formation of a critical bubble nucleus and temperature [START_REF] Han | Bubble nucleation in polymeric liquids. II. Theoretical considerations[END_REF]. The theory suggests that the energy barrier and the interfacial tension decrease as the supersaturation pressure increases. Consequently, the homogeneous nucleation rate increases and a large number of small cells is obtained. In fact, controlling supercritical CO 2 quantity allows the modification of the cell density and the cell size in microcellular foaming.

The objective of this work is to develop a process of extrusion coupled with the use of scCO 2 . In a preliminary work, an extruder was modified in order to be able to inject scCO 2 within the extruded material and some experiments were performed on polystyrene (PS) [START_REF] Sauceau | Développement de procédés d'extrusion assistés par un fluide supercritique[END_REF]. A pressure decrease within the extruder was observed upon CO 2 injection and samples of porous PS were obtained, with a decrease of their apparent density from 30 to 40%. However, the porous structure was irregular, with bubble diameters varying from 100 m to 1 mm. This could be explained considering the quality of the mixture in the extruder. Indeed, an insufficient homogeneity of mixture PS-CO 2 could explain the irregularity of porosity. Therefore, the experimental device has been modified in order to improve the mixture between the two components. In this work, this novel device has been tested with PS and the influence of cooling, die temperature and CO 2 content on the material properties has been studied.

Experimental

Materials and equipment

Polystyrene (choc, BP Chemicals) was used in experiments and CO 2 (N45, Air Liquide) was applied as the physical foaming agent.

Fig. 1 shows the flow diagram of the foaming system. The extruder used in this study is a single-screw machine with a 20 mm screw diameter and a length to diameter ratio (L/D) of 20 (Rhéoscam, SCAMEX). Carbon dioxide is pumped from a cylinder by a syringe pump (260D, ISCO), which allows the delivery of small flow rates (i.e. 0.1 l mn -1 ). The pump is run in constant volumetric flow rate mode and CO 2 is introduced at the same pressure as the pressure prevailing in the extruder. The injection position is located at about 12D from hopper, in a zone where the screw diameter is constant. A checking valve is inserted in the line to prevent pollution of CO 2 . The polymer-CO 2 mixture is then passed through a mixer (type 'pineapple'), composed of diamond spike. This element, inserted before the die, provides a dynamic distributive mixing between the two components [START_REF] Vergnes | Extrusion monovis (partie 1)[END_REF]. The pressure drop through the die can be adjusted by modifying restriction opening by means of a screw placed in the polymer flow.

The pressure, the temperature and the volumetric CO 2 flow rate are measured within the syringe pump. CO 2 density is obtained by Peng Robinson's equation of state [START_REF] Peng | A new two-constant equation of state[END_REF]. In operating conditions used, CO 2 concentration is always smaller than saturation concentration. Consequently, polystyrene-CO 2 blend is in a single phase at thermodynamic equilibrium [START_REF] Sato | Solubilities of carbon dioxide and nitrogen in polystyrene under high temperature and pressure[END_REF][START_REF] Sato | Solubilities and diffusion coefficients of carbon dioxide in poly(vinyl acetate) and polystyrene[END_REF].

The temperature inside the barrel is regulated at four locations: T 1 , T 2 before and after the CO 2 injection, T 3 before the mixing element and T 4 before the die. Only the temperature of the die T 4 is changed at each experiment. Three other temperatures and screw speed are kept constant at 180 • C and 50 rpm, respectively. At each new condition, temperature T mat and pressure P mat of material are measured just before entry within the die. The flow rate of polymer Q mat , of the order of 0.6 kg h -1 , is obtained by weighing the extrudate at the end of the experiment. CO 2 concentration C CO 2 in the polystyrene can thus be expressed as weight flow rate ratio.

A large decrease of the pressure into the extruded matter was observed when scCO 2 was injected, mainly due to the decrease of melt viscosity [START_REF] Kwag | Rheology of molten polystyrene with dissolved supercritical and near-critical gases[END_REF]. For all experiments, 3-5 mm outer diameter PS cylinders were obtained. Once stable conditions have been established, samples were collected and submitted to two kinds of cooling. The first one consists in leaving the polymer at ambient conditions at about 22 • C. In the second one, material is plunged into a water bath at ambient temperature immediately after its exit from the die.

To study the influence of operating parameters on porous structure, several experiments have been carried out at 12 MPa varying cooling mode, die temperature and CO 2 concentration, other parameters being constant. Experiments are time and material consuming and we have thus chosen to study die temperature and CO 2 content at three levels to observe the main tendencies of the influence of these operating parameters.

Characterization

Porosity ε is defined as the ratio of void volume to the total volume of the sample and can be calculated by the following equation:

ε = V e V t = 1 - V m V t with V t = V e + V m (1)
where V t is the total volume of the sample, V e the void volume and V m is the solid volume. Two methods are used to characterize porosity: the water-pycnometry and the Hg-porosimetry (Autopore III, Micrometerics).

The water-pycnometry allows the measurement of the apparent density ρ app of expanded samples and thus to calculate the total porosity ε:

ε = 1 - ρ app ρ PS ( 2
)
where ρ ps is the bulk density of raw polystyrene, which corresponds to polystyrene extruded without CO 2 injection. This density, also measured by water-pycnometry, is equal to 1 g cm -3 . In Hg-porosimetry, mercury is used as penetrating fluid in a calibrated pycnometer. It is forced to enter into the pores by means of an increasing pressure related to the diameter of the filled pore. Each pressure level required to fill the sample provides the diameter D p and the volume V p of corresponding pores. In this method, only the opened porosity ε eff is filled by mercury and can thus be calculated with the following equation:

ε eff = V p V t (3)
To quantify the fraction of the total porosity which is opened, an open porosity rate X eff can be defined as the ratio of the opened porosity to the total porosity:

X eff = ε eff ε (4)
It has to be noted that a minimal pressure is required to begin a measurement in Hg-porosimetry. This lowest usable pressure corresponds to a maximal value of pore diameter equal to about 200 m.

To complete the characterization of the porosity structure, samples were examined by scanning electron microscopy (XL30 ESEM, FEG, Philips).

Results and discussion

Influence of the cooling mode

Fig. 2 compares three PS samples processed at different die temperatures. All of them have an external skin with a roughly constant thickness. It is probably formed by an effect of hardening during the cooling. The differences appear in the centre of the cylinders. Indeed, with increasing temperature, pores become larger and more irregular. Finally, for highest temperatures, large pores occupy all the central part. This effect is more important in the rod centre, which is more slowly cooled. Moreover, at lower temperature, the porous structure is improved by the addition of the mixer [START_REF] Sauceau | Développement de procédés d'extrusion assistés par un fluide supercritique[END_REF].

Table 1 shows the evolution of porosity ε as a function of temperature. An important increase of ε occurs from 175 to 185 • C, probably due to growth of the pores which conserve a 3 compares two samples processed in the same operating conditions but cooled by using the two different ways previously described. At this temperature, the influence of the cooling is obvious: the cooling in water allows smaller and homogeneous pores and a porosity reduced from 30 to 19%. Growth and coalescence are attenuated and it is thus important to solidify rapidly the structure to carry out relevant observations on the influence of operating parameters. Arora et al. [START_REF] Arora | Preparation and characterization of microcellular polystyrene foams processed in supercritical carbon dioxide[END_REF] have used a batch process to foam PS by using CO 2 between 40 and 120 • C. They have also observed fewer and larger pores at higher temperatures. This can be explained by the decrease of the viscosity of the substrate material causing the restrictive force restricting cell growth to decrease and the increase of the diffusivity of CO 2 within the substrate. With a continuous process, Park et al. [START_REF] Park | Low density microcellular foam processing in extrusion using CO 2[END_REF] have shown that the improvement of pore structure is obtained at the price of the decrease of expansion. They have thus observed an increase of cell coalescence with increasing die temperature. More recently using also a continuous process, Han et al. [START_REF] Han | Effect of die temperature on the morphology of microcellular foams[END_REF] have done the same observation, with however coupled effects of pressure and temperature. Indeed, the authors have worked with constant die geometry and, a temperature decrease has involved a pressure drop increase because of the viscosity increase, which may also influence the structure.

Effect of die temperature

The relationship between die temperature and homogenous nucleation rate is complex because many physical properties depend on this temperature. Indeed, by using classical theory, it can be shown that increasing the temperature enhances the homogenous nucleation rate [START_REF] Han | Effect of die temperature on the morphology of microcellular foams[END_REF]. Thus, as previously observed in the case of air cooling, it seems that the observed evolution cannot be explained by nucleation theory but rather by an increase in the growth and coalescence phenomena, involving fewer and larger pores [START_REF] Arora | Preparation and characterization of microcellular polystyrene foams processed in supercritical carbon dioxide[END_REF]. Fig. 5 shows the evolution of porosity with temperature. At higher CO 2 concentration, total porosity ε increases from 180 to 200 • C and then remains almost constant from 200 to 220 • C. At the lowest CO 2 concentration, the temperature exhibits no more influence on porosity. At the same time, opened porosity ε eff decreases slowly with increasing temperature (Fig. 6). This leads to the decrease of the opened porosity rate ε eff from 68% at 180 • C to 37% at 220 • C (Fig. 7). Growth is therefore the predominant effect at lower temperatures, while cell coalescence occurs at higher temperatures, involving the closing of these cells.

The pore size distribution obtained by Hg-porosimetry is represented at different temperatures on Figs. 8 and9. At 180 • C, diameter of pores is between 5 and 20 m. When the temperature increases, pore size increases to values higher than 200 m, which is the limit of the apparatus. This evolution confirms the former conclusion that coalescence and growth phenomena occur when temperature increase and lead to larger porosity.

Effect of CO 2 weight ratio

On SEM pictures, the effect of CO 2 concentration is less obvious than that of the die temperature (Fig. 4). On Figs. 5 and6, it can be seen that total porosity ε and opened porosity ε eff decrease slightly with increasing CO 2 concentration, with however a higher effect at lower temperature. This finally leads to an opened porosity rate ε eff almost constant (Fig. 7). Pore size distributions obtained by Hg-porosimetry confirm the lim- ited influence of CO 2 concentration, as for instance at 200 • C (Fig. 10).

It is clear that increasing CO 2 concentration increases nucleation rate [START_REF] Han | Effect of die temperature on the morphology of microcellular foams[END_REF]. However, CO 2 concentration increase involves also the modification of other key physical variables, such as for instance properties like viscosity or surface tension [START_REF] Gourgouillon | Simultaneous viscosity and density measurement of supercritical CO 2 -saturated PEG 400[END_REF][START_REF] Condo | Glass transition behavior including retrograde vitrification of polymers with compressed fluid diluents[END_REF][START_REF] Li | Effect of carbon dioxide on the interfacial tension of polymer melts[END_REF][START_REF] Park | Surface tension measurement of polystyrene melts in supercritical carbon dioxide[END_REF]. This difficulty is well illustrated by Han et al. [START_REF] Han | Continuous microcellular polystyrene foam extrusion with supercritical CO 2[END_REF] who have also studied the effect of CO 2 concentration. By means of a computational fluid dynamics program, these authors have represented profiles of pressure and saturation CO 2 fraction along the nozzle. Because of their constant die geometry, the increase of CO 2 concentration was coupled with a pressure decrease. They have then found that the more CO 2 in the melt, the earlier the bubbles begin to nucleate and grow. In other words, a high CO 2 content gives more nuclei, but a longer time for bubble growth. However, during their experiments, it was noted that both the pressure and the pressure drop rate are independent of CO 2 content. Finally, they observed that, at low CO 2 concentrations, cell size decreases and cell density increases with increasing CO 2 content. While at high CO 2 concentrations (above 4%), cell size tends to be independent of CO 2 content, probably because of extra gas, which is not used to nucleate cells. If we translate these experimental observations in terms of total porosity, no clear trend appears. Indeed, cell size increase and cell density decrease show opposite effects on total porosity.

In this work, no significant influence of CO 2 concentration has been observed. However, PS-CO 2 blend is theoretically in a single phase at thermodynamic equilibrium, as CO 2 weight solubility was estimated to be at about 5% at 180 • C and 12 MPa [START_REF] Sato | Solubilities of carbon dioxide and nitrogen in polystyrene under high temperature and pressure[END_REF]. Thus, as observed by Han et al. [START_REF] Han | Continuous microcellular polystyrene foam extrusion with supercritical CO 2[END_REF], it seems that a part of injected CO 2 remains as extra gas and is not used to nucleate cells. Furthermore, higher CO 2 content than studied range usually causes unstable polymer flow with oscillating pressure, what would tend to confirm this assumption.

Conclusion

A process based on extrusion coupled with the use of supercritical carbon dioxide (scCO 2 ) was implemented. scCO 2 modifies the rheological properties of the material in the barrel of the extruder and acts as a blowing agent during the relaxation at the passage through the die. A novel experimental device based on a single-screw extruder was developed. It allows the injection of scCO 2 into extruded melt, the mixing of both compounds and the creation of porosity into the extrudated polymer. In this work, it was applied to polystyrene (PS). Operating temperatures and mechanical constraints within the extruded material were reduced by the injection of CO 2 . The effects of cooling mode, die temperature and CO 2 concentration on material porosity were studied. PS with porosity from 15 to 25% was manufactured. A rapid cooling is important to solidify the structure. Increasing die temperature increases growth and coalescence phenomena, and thus increases total porosity and pore diameter. Finally, this parameter allows the control of the porosity structure. CO 2 concentration shows no significant influence. It was explained by the fact that a part of injected CO 2 remains as extra gas and is not used to nucleate cells.

Fig. 1 .

 1 Fig. 1. Flow diagram of the extruder.

Fig. 2 .

 2 Fig. 2. SEM pictures of samples cooled in air processed at P mat = (12 ± 1) MPa, C CO 2 = (5 ± 0.2)% and (a) T mat = 175 • C; (b) T mat = 185 • C; (c) T mat = 200 • C.

Fig. 4

 4 Fig.4shows SEM pictures of samples obtained at different die temperatures. Increasing the temperature increases the pore size and decreases the pore number. It appears a radial distribution of pore size, with larger pores in the cylinder centre. The shape evolves also with temperature, the pores becoming more regular and finally almost spherical at 220 • C.Arora et al.[START_REF] Arora | Preparation and characterization of microcellular polystyrene foams processed in supercritical carbon dioxide[END_REF] have used a batch process to foam PS by using CO 2 between 40 and 120 • C. They have also observed fewer and larger pores at higher temperatures. This can be explained by the decrease of the viscosity of the substrate material causing the restrictive force restricting cell growth to decrease and the increase of the diffusivity of CO 2 within the substrate. With a continuous process, Park et al.[START_REF] Park | Low density microcellular foam processing in extrusion using CO 2[END_REF] have shown that the improvement of pore structure is obtained at the price of the decrease of expansion. They have thus observed an increase of cell coalescence with increasing die temperature. More recently using also a continuous process, Han et al.[START_REF] Han | Effect of die temperature on the morphology of microcellular foams[END_REF] have done the same observation, with however coupled effects of pressure and temperature.

Fig. 3 .

 3 Fig. 3. SEM pictures of samples processed at P mat = (12 ± 1) MPa, T mat = 200 • C, C CO 2 = (5 ± 0.2)% and cooled (a) in air; (b) in water.

Fig. 4 .

 4 Fig. 4. SEM pictures of samples cooled in water processed at P mat = (12 ± 1) MPa.

Fig. 5 .

 5 Fig. 5. Total porosity ε as function of CO 2 concentration at different temperatures.

Fig. 6 .Fig. 7 .

 67 Fig. 6. Opened porosity ε eff as function of CO 2 concentration at different temperatures.

Fig. 8 .

 8 Fig. 8. Pore size distribution with 2.2% of CO 2 at different temperatures.

Fig. 9 .

 9 Fig. 9. Pore size distribution with 3.8% of CO 2 at different temperatures.

Fig. 10 .

 10 Fig. 10. Pore size distribution at 200 • C and different CO 2 concentrations.

Table 1

 1 Porosity ε at different temperatures T mat of samples processed at P mat = (12 ± 1) MPa, C CO 2 = (5 ± 0.2)% and cooling in air

	T mat ( • C)	ε (%)
	175	16
	185	31
	200	30
	regular shape. Then, porosity ε remains constant from 185 to
	200 • C. It seems that growth is attenuated and that coalescence phenomena are at the origin of the structure observed at 200 • C.
	Fig.