J Monnier

A M Benselama

I Cotoi

Flow Patterns in the Vicinity of Triple Line Dynamics arising from a Local Surface Tension Model

We model and simulate numerically a droplet impact onto a solid substrate. The triple line dynamics modelling is implicit (as opposed to classical explicit mobility relations), it is based on the Shikhmurzaev equations. These equations include generalized Navier slip type boundary conditions with extra local surface tension gradient terms. Numerical results when spreading are presented. A particular attention is paid to flow patterns near the contact line.

Introduction

One of the main difficulties to simulate present droplet flows is the contact line dynamics modelling. This problem has been widely studied and still remains an open problem. Let us cite macroscopic models with slip boundary conditions, see e.g. [START_REF] Dussan | On the motion of a fluid-fluid interface along a solid surface[END_REF], [START_REF] Hocking | A moving fluid interface. Part II. The removal of the force singularity by a slip flow[END_REF], [START_REF] De Gennes | Wettings: statics and dynamics[END_REF], mesoscopic models with diffuse interface, see e.g. [START_REF] Jacqmin | Contact-line dynamics of a diffuse fluid interface[END_REF] and also [START_REF] Pomeau | Recent progress in the moving contact line problem: a review[END_REF], and connection with molecular dynamics, see e.g. [START_REF] Qian | Generalized Navier boundary condition for the moving contact line[END_REF]. The most frequent contact line model is an explicit mobility relation giving the contact line velocity in function of the wetting angle value at any time (Tanner type laws). Nevertheless, it has been showed, [START_REF] Blake | Experimental evidence of nonlocal hydrodynamic influence on the dynamic contact angle[END_REF], that the dynamic wetting angle cannot be determined inherently by such a mobility relation.

In the present paper, we consider an "implicit" model based on the Shikhmurzaev theory, [START_REF] Shikhmurzaev | The moving contact line problem on a smooth solid surface[END_REF], [START_REF] Blake | Experimental evidence of nonlocal hydrodynamic influence on the dynamic contact angle[END_REF]. This model lead to generalized Navier slip type boundary conditions with local surface tension gradient terms. It does not impose the contact line velocity nor the wetting angle since they are a response of the full model. In Section 2, we present the equations: 2D axi-symmetric Navier-Stokes equations with an ALE formulation. In Section 3, we discretize the equations using a finite element scheme. We introduce a contact line algorithm which imitates caterpillar motion observed in experiments. Curvature is computed using a local Bezier least-square approximation. In Section 4, first we consider a Tanner type law and we test the robustness of our algorithm. Then, we consider both spreading phase and recoiling phase, using the present "implicit" model based on the Shikhmurzaev theory. The extra local terms appearing in this model are set by preliminary results presented in [START_REF] Monnier | Modélisation Numérique de la Ligne Triple[END_REF] and [START_REF] Monnier | Analysis of a Local Hydrodynamic Model with Marangoni Effect[END_REF]. We show different flow patterns generated by these terms near the triple point. Finally, we obtain encouraging results for both phases.

Mathematical Model

The equations

The droplet dynamics is modelled by the2D axi-symmetric unsteady incompressible Navier-Stokes equations. We denote by u = (u r , u z ) T the fluid velocity, p its pressure, Σ the stress tensor, D the deformation 1

tensor and Re the Reynolds number. We denote by ( τ , n) the unit tangential and external normal vectors such that it is direct. We set: Σ n = Σ. n; Σ n = Σ n n + Σ τ τ .

Free surface dynamics and ALE formulation The free surface Γ f r is transported by the velocity field u. It is described by an implicit function: φ(t, x(t), y(t), z(t)) = 0, and we have:

dφ dt = ∂φ ∂t + u • ∇φ = 0 in (0, T ) × Ω (1)
with given initial conditions.

We present the ALE formulation we use as well as the characteristics method. We refer to [START_REF] Hughes | Lagrangian-Eulerian Finite Element Formulation for incompressible Viscous Flows[END_REF], [START_REF] Maury | Characteristics ALE method for unsteady free surface flows with surface tension Z[END_REF].

The principle of ALE method is to define an equivalent velocity field γ (in the sense γ. n = u. n on ∂Ω)

such that it preserves the mesh inside the fluid domain Ω. To this end, at any time t, we define the deformation field γ on the free surface Γ f r , then we extend it all over the domain by solving a linear elasticity system.

Let t → C(x, τ ; t) be the characteristics lines associated to the velocity field γ, then one can define the ALE variables:

u τ = u • C and p τ = p • C.
In the ALE space, fluid particles have a velocity ( u -γ) at first order in time (t -τ ). And for t = τ one has:

∂ u ∂t = ∂ u τ ∂t -γ τ ∇ u τ , with γ τ = γ • C.
If we denote by t → X(x, τ ; t) the Lagrangian characteristic lines (i.e. those associated to the velocity field u), ALE method consists to define a regular field γ such that:

Γ t f r = C(Γ 0 f r , t 0 ; t) = X(Γ 0 f r , t 0 ; t).
The r-momentum equation with ALE formulation is:

Re ∂u τ r ∂t + ( u τ -γ) • ∇u τ r = - ∂p τ ∂r + 2 1 r ∂(r D rr ( u τ )) ∂r + ∂(D rz ( u τ )) ∂z - u τ r r 2 (2)
The z-momentum equation is similar. The continuity equation is:

1 r ∂(ru τ r ) ∂r + ∂u τ z ∂z = 0
with initial conditions. Boundary conditions on the free surface (liquid-gas) are:

Σ n = (-p ext + κ τ Ca ) n + h τ in (0, T ) × Γ f ( 3
)
where Ca is the Capillary number, κ is the mean curvature and p ext is the external pressure. The extra term h is given (see Shikhmurzaev's theory below). We have classical boundary conditions on the symmetry axis. We decompose the liquid-solid interface into two parts: Γ ad and Γ sl . Γ sl denotes a "small" part of the liquid-solid interface near the triple point. We consider a generalized Navier slip type boundary condition (local slipping):

u τ . n = 0 in (0, T ) × Γ sl Σ τ = -[β u τ + g]. τ in (0, T ) × Γ sl ( 4
)
where β is a given sliding coefficient. The extra term g is given, it models either a surface tension gradient in the Shikhmurzaev's theory or a uncompensated Young stress in the Qian-Wang-Sheng theory, see below.

On Γ ad , we impose adherence boundary conditions: u = 0.

The solution ( u τ , p τ ) • C is an approximation of ( u, p) at first order in time, and curvatures κ τ approximates κ at first order in time too, see e.g. [START_REF] Maury | Characteristics ALE method for unsteady free surface flows with surface tension Z[END_REF].

In all the sequel we simplify notations by dropping the arrows and the superscript τ . 2

Contact line dynamics modelling

We consider two different types of model for the contact line dynamics: an explicit model (Tanner type law) and an "implicit" one deriving from Shikhmurzaev theory.

Mobility relation (Tanner type Law). The contact line velocity and the wetting angle are related by:

U CL (t) = k (θ(t) -θ eq ) γ (θ eq ) γ for t ∈]0, T [ (5)
where U CL is the contact line velocity, θ eq is the wetting angle at equilibrium (Young's law), θ is the (dynamic) wetting angle. k and γ are parameters usually determined using experimental results.

Local flow modelling and Shikhmurzaev's theory. The "implicit" model does not impose the contact line velocity and the wetting angle but consider them as a response of the model. This model is based on the Shikhmurzaev's theory, [START_REF] Shikhmurzaev | The moving contact line problem on a smooth solid surface[END_REF], [START_REF] Blake | Experimental evidence of nonlocal hydrodynamic influence on the dynamic contact angle[END_REF], which introduces the generalized Navier slip condition (4) and the condition (3). These conditions are local since the extra terms g and h vanish except in a vicinity of the triple point. The basic idea of this theory is to consider that the rolling motion observed in experiments, [START_REF] Dussan | On the motion of a fluid-fluid interface along a solid surface[END_REF], implies that particles of the liquid-gas interface, whom become an element of the solid-liquid, lose their properties in a finite time. Then the surface tension value associated to this particle must change to its new equilibrium value relative to the solid-liquid interface. This process would gives rise to a surface tension gradient in a small vicinity of the advancing contact line (hence a local Marangoni effect). In other respect, the Young equation would remain valid at any time. In this theory,

g = - 1 2Ca ∇σ and h = 1 Ca ∇σ LG . τ (6)
where σ and σ LG are the liquid-solid and the liquid-gas surface tension coefficient respectively. In [START_REF] Monnier | Modélisation Numérique de la Ligne Triple[END_REF] and [START_REF] Monnier | Analysis of a Local Hydrodynamic Model with Marangoni Effect[END_REF], a mathematical and numerical study presents some qualitative behaviors of g and h arising from Shikhmurzaev's theory. The mesoscopic Local Surface Model Briefly, the so-called mesoscopic LSM as it is established in [START_REF] Shikhmurzaev | The moving contact line problem on a smooth solid surface[END_REF] is as follows. The interfaces are described by surface densities ρ s (a variable that describes unambiguously the surface state). These surface densities are solution of surface continuity equations.

A linear state equation gives the relation between ρ s and the surface tension coefficients σ.

We denote: by ρ s i , i = 1, 2, the surface density on Γ f r (i=1) and on Γ solid (i=2), by σ 1 = σ and σ 2 = σ SL . The surface tension is related to the excess density through the linear state equation:

σ i = γ(ρ s 0 -ρ s i ) i = 1, 2. ( 7
)
where γ and ρ s 0 are given constants. According to [START_REF] Billingham | A reappraisal of a model for the motion of a contact line on a smooth solid surface[END_REF], γ ≈ 2.10 6 and ρ s 0 ≈ 10 -7 . We have the surface continuity equation:

∂ρ s i ∂t + div(ρ s i v s i ) + 1 τ * (ρ s i -ρ eq i ) = 0 i = 1, 2 (8)
where τ * is the relaxation time relative to the rolling motion, τ * ≈ 10 -3 -see [3]-, v s i is a mean velocity inside the layer and ρ eq i is its density at equilibrium, [START_REF] Shikhmurzaev | The moving contact line problem on a smooth solid surface[END_REF]:

σ i (ρ eq i ) = σ eq i , i = 1, 2. The velocity v s 1 (respectively v s 2 ) is related to ρ s 1 (respectively ρ s 2 )
and to the fluid velocity u. We have the following Darcy laws type, [START_REF] Shikhmurzaev | The moving contact line problem on a smooth solid surface[END_REF]:

(1 + 4α 1 α 2 ) ∇σ f r = 4α 2 ( v s 1 -u) and v s 2 = α 1 ∇σ solid + 1 2 u (9)
where α i , i = 1, 2, are given constants characterizing the viscous properties of the interface. According to [START_REF] Billingham | A reappraisal of a model for the motion of a contact line on a smooth solid surface[END_REF], α 1 ≈ h µ , h being the surface layer thick, h ≈ 10 -10 .

At the triple line, the surface flux continuity is imposed:

(ρ s 1 v s 1 ) e f = (ρ s 2 v s 2 ) e g ( 10
)
where e f and e g are unit vectors normal to the contact line and tangential to the gas-liquid and gas-solid interface respectively. Let us notice that: cos(θ) = -e f . e g .

Remark.

The most important feature of the Shikhmurzaev model is supposed to be the following.

It does not impose neither the contact line velocity nor the wetting angle since both are a full response of the model. Nevertheless the LSM is constituted by the two surface equations ( 8) and the boundary condition [START_REF] Maury | Characteristics ALE method for unsteady free surface flows with surface tension Z[END_REF] at contact line. Then, as it is noticed in ( [START_REF] Billingham | A reappraisal of a model for the motion of a contact line on a smooth solid surface[END_REF]), the condition (10) at contact line is not enough to close the system. To this end the problem misses an extra condition at contact line. This can be very easily showed for some linear second order differential equations with only a flux continuity condition at junction point.

The missing condition cannot be the Young equation otherwise it would imply that the wetting angle must be given, and it would not be a response of the model anymore. The author of [START_REF] Billingham | A reappraisal of a model for the motion of a contact line on a smooth solid surface[END_REF] proposes an extra condition of the following type:

(ρ s 1 v s 1 ) e f = F (ρ s 1 , ρ s 2 ) (11)
where F (., .) is a bilinear given function related to a chemical potential law.

A connection with Qian-Wang-Sheng theory. From molecular dynamics simulations on immiscible fluids, the authors of [START_REF] Qian | Generalized Navier boundary condition for the moving contact line[END_REF], [START_REF] Qian | Molecular hydrodynamics of the moving contact line in two-phase immiscible flows[END_REF] show that the relative slipping between the fluid and the solid wall follows a generalized Navier slip boundary condition similar to [START_REF] De Gennes | Wettings: statics and dynamics[END_REF]. In this theory, the extra term g in ( 4)

would model the contribution of the tangential stress due to non-shear stress that appears because of the anisotropy of the pressure (thus its non-shear origin) across the fluid-fluid interface; normal pressure term is no longer equal to tangential pressure term across this interface, as this usually occurs in the bulk. Hence this term, which is present as far as two(fluid)-phase region is involved, is connected to the fluid-fluid interfacial tension. The extra tangential stress g, called interfacial uncompensated (or unbalanced) Young stress, is defined as follows: Γint g dy = σ(cosθ -cosθ eq ), where Γint dy denotes the integral across the interface Γ f r .

From a macroscopic point of view, it is obvious that it arises from the deviation of the fluid-fluid interface from its static configuration and is narrowly distributed in the fluid-fluid interface region. Therefore, both theories (Shikhmurzaev's and molecular ones) lead to such a generalized Navier slip boundary condition (equation ( 4)). The derivation and the interpretation of the extra term g is different for each theory but in a mathematical modelling point of view, the boundary condition is similar. Moreover, the authors of [START_REF] Qian | Generalized Navier boundary condition for the moving contact line[END_REF], [START_REF] Qian | Molecular hydrodynamics of the moving contact line in two-phase immiscible flows[END_REF] point out the fact that, independently of the nearness to the fluid-fluid interface, the generalized Navier slip condition is a more general and suitable condition than the classical non-slip condition, provided slip coefficient β and g are well adjusted.

3 Discretization and Algorithms

Finite element scheme

We discretize the full model in time using the Euler implicit scheme and a first order characteristic method. Let {t 0 = 0, t 1 , ..., t M = T } be a discretization of ]0, T [ with a constant time step ∆t = T M . Semi-discrete r-momentum equation is (similar discretisation of z-momentum equation is done):

Re u m+1 r -u m r • χ m ∆t = - ∂p m+1 ∂r + 2 1 r ∂(r D rr ( u m+1 )) ∂r + ∂(D rz ( u m+1 )) ∂z - u m+1 r r 2 (12)
and

Re u m+1 z -u m z • χ m ∆t = - ∂p m+1 ∂z + 2 1 r ∂(r D zr ( u m+1 )) ∂r + ∂(D zz ( u m+1 )) ∂z (13)
where the superscript m represent the evaluation of the quantity at t = t m and χ m denotes the characteristic associated to the velocity field ( u -γ), defined by:

dχ m dt (t) = ( u -γ)(χ m (t), t) t ∈]t m , t m+1 [ χ m (t m+1 ) = x (14)
Then, χ m (x) is approximated by: χ

m (x) ≈ x -∆t( u m -γ m )(x).
The semi-discrete z-momentum equation is similar.

We set:

u ax = ( u r r 2 , 0) T
and we denote by

div rz (•) = 1 r ∂ ∂r (r (•)) + ∂ ∂z (•).
Then we can write ( 12)-( 13) as follows:

λ u m+1 -2div rz (D( u m+1 )) + ∇p m+1 + 2 u m+1 ax = λ ( u m • χ m ) ( 15
)
with λ = Re dt .

Let us recall that ( u, p) does not denote the Eulerian velocity-pressure but its ALE counterpart.

Equation ( 1) is semi-discretized using the Euler implicit scheme and the first order characteristic method.

This gives:

φ m+1 = φ m • X m in Ω, m = 1, 2... ( 16
)
where X m denotes the characteristic associated to the velocity field u at time t m . It is approximated as

follows: X m (x) ≈ x -∆t u m (x).
The space discretization is done using the second order Taylor-Hood finite element. The scheme is implemented using a public C++ finite element library Rheolef, [START_REF] Saramito | The finite element system[END_REF].

The ALE Free Surface Algorithm

The free surface problem is solved using an ALE formulation and a characteristics method. The algorithm is the following. At time step n, given Ω n and ( u n , p n ) in Ω n .

Step 1. Compute the new free surface position Γ n f following the Lagrangian characteristic lines (each point of the boundary is translated by u n ∆t). Then, compute an equivalent deformation field γ n preserving the mesh by solving a linear elasticity system.

Step 1bis. Compute the new solution (

u n+1 τ , p n+1 τ ) of Navier-Stokes in Ω n .
Step 2. Update the domain Ω n+1 = γ n (Ω n ) (mesh transport by γ n ).

Step 3. Obtain the new solution ( u n+1 , p n+1 ) in Ω n+1 by setting:

u n+1 | Ω n+1 = u n+1 τ | Ω n and p n+1 | Ω n+1 = p n+1 τ | Ω n .

Spreading and caterpillar motion

When spreading our algorithm imitates the caterpillar motion observed in experiments, [START_REF] Dussan | On the motion of a fluid-fluid interface along a solid surface[END_REF]. Fig. 3 represents typical successive time steps. If at time step n, the first point (or several points) on the free surface are projected on the solid substrate, all the points on the solid wall but the last ones become no-slip points. The last two points remain slip point.

Recoiling motion

When de-wetting (recoiling), the triple point and the slip points are allowed to move horizontally, also vertically if their vertical velocity is positive, see Fig. 4. If the triple point moves away from the solid wall, the first point to the left becomes the triple point, hence becoming a free point. The second point to the left becomes a slip point.

Curvature Computation

We seek to compute the curvature of the droplet surface. Computing accurately the curvature is a difficult task since its interface is a piecewise linear curve, hence not C 2 differentiable. In addition, points defining this piecewise linear curve result from the full free surface algorithm, hence may comprises some nonnegligible numerical errors. Then, we seek to estimate the curvature of an underlying smooth surface.

Computing a discrete surface curvature is a classical (and difficult) problem. Usually in the CAGD context, surfaces are 3D and triangularized, and the objectives are to smooth the mesh, to simplify it, but not to quantify a local variation of curvature, see e.g. [START_REF] Meyer | Discrete differential-geometry operators for triangulated 2-manifolds[END_REF].

The present problem is easier since the curve is 2D only but we are facing the following dilemma: we seek to get rid of numerical errors connected to the points defining the curve while we seek to detect as accurate as possible local variation of curvature, in order to obtain the resulting surface tension forces accurately.

We do not consider a direct computation by a finite difference method since it is very sensitive to data error. We do not consider neither a polynomial reconstruction of the underlying smooth surface then evaluate its curvature, since it leads to inaccurate results and unexpected behavior. Following [START_REF] Hamann | Curvature approximation for triangulated[END_REF], we consider a local least-square approximation then we evaluate the curvature. In the present algorithm, we consider a second order local Bezier approximation. As the numerical tests show (see below), this method filters noise quite reasonably while it allows to detect local and rapid variation of curvature.

The Algorithm

Given N points X i = (r i , z i ) T , i = 1..N , defining the liquid-gas interface, the main idea is to approximate these data using a local least-square approximation by a Bezier's curve.

The Bezier's curve C(t) is given by:

C(t) = (r(t), z(t)) T = M j=1 P j B M -1 j-1 (t) for t ∈ [0, 1]
where P j = (α j , β j ) T ∈ IR 2 are the control points and {B m j (t)} 0≤j≤M -1 is the classical Bernstein's basis, B m j ∈ P m , B m j (t) = C m j (1 -t) m-j t j , C m j being the binomial coefficients.

We set M = 3, hence we consider 3 points of control P j and second degree curves.

For an inner point X i , see Fig. 5, we compute the least square approximation of the 5 points {X i-2 , .., X i+2 } by a Bezier's curve as follows. We minimize:

J(P 1 , .., P M ) = i+2 l=i-2 M j=1 P j B M -1 j-1 (t l ) -X l 2
where {t i-2 = 0, .., t i+2 = 1} is an uniform subdivision of [0, 1]. The unique minimum is computed by solving the corresponding normal equations.

For the extremal point X 1 , we consider a Bezier's curve approximating the points X i for i = 1, .., 4.

For X 2 , we consider a Bezier's curve approximating the points X i for i = 1, .., 5. For the extremal points X N -1 and X N , the principle is similar.

z X i X i-1 X i-2
Curvature expression. Once a Beziers's curve C(t) = (r(t), z(t)) T is computed for each point X i , we evaluate the (plane) curvature as follows:

κ i = κ(t i ) = r z" -r"z (r 2 + z 2 ) 3 2 (t i )
where (r , z ) and (r , z ) are computed using de Casteljau's algorithm, t i being the parameter value corresponding to X i .

Sensitivity to random noise. Since the control points defining the (optimal) droplet shape are resulting from the full free surface algorithm, they may be perturbed by some non-negligible numerical errors.

Hence, we test the robustness of our algorithm to data perturbation (see below).

If we denote: N (r, z) = (r z" -r"z ) and D(r, z) = (r 2 + z 2 ) 3 2 , then: κ(r, z)(t) = N (r,z) D(r,z) (t). Let δz be a perturbation on z-coordinate of data X i , i = 1..N , then we have:

∂κ ∂z (r, z).δz = N (r, δz) D(r, z) -3 κ(r, z) (r 2 + z 2 ) z δz
This formulae expresses the curvature sensitivity to perturbation on z-coordinate.

Noise introduced in the numerical tests below is a random perturbation (on the z-coordinate of data

X i , i = 1..N ).
It is a normal distribution with zero mean value and unit variance.

Numerical tests

We present some numerical tests using two different methods:

Method a. Standard second order finite difference scheme directly applied to the N data

X i = (r i , z i ) T , i = 1..N , Method b.
Local approximation by a second order Bezier's curve then curvature computation.

In order to compare these two methods, we consider the following three curves (with noise introduced): As mentioned above, noise is defined as a perturbation on the z-coordinate of data X i , i = 1..N .

The circle is defined by N points as follows:

r(i) = Rcos(η i ), z(i) = Rsin(η i ), with η i = ( (i-1)π 2(N -1) ), i = 1..N .
The so-called "oscillating-curve" is defined by N points as follows:

r(s) = (R + cos(a.s)) cos( π 2 s), z(s) = (R + cos(a.s)) sin( π 2
s) with s ∈ [0, 1], s discretized by N points similarly to η and = R 10 , a = 10. The exact curvature of the "oscillating-circle" is straightforwardly obtained. This curve is interesting since it presents smooth variations of curvature with change of sign.

The so-called "double-circle" is defined by N points as follows:

For i = 1..M , r(i) = Rcos(η i ), z(i) = Rsin(η i ); for i = M +1..N , r(i) = a+ R 2 cos(η i ), a = (1-R)cos(η M ) and z(i) = b + R 2 sin(η i ), b = (1 -R)sin(η M
). This curve is interesting since it presents a sharp variation of curvature (non-continuous in fact).

We set R = 1, M = 8 and N = 30.

Without noise (i.e. if considering the exact curves), the two methods allow to retrieve perfectly the curvature values of the three curves.

In figures 6, 7 and 8, we also compare the two methods a. and b. by computing the curvature of curves 1, 2 and 3 respectively when some noise is introduced. In such a case, all our computational tests we performed showed that method b. gives much better results than method a.

In other respects, we compared methods a. and b. with a "global" method: we interpolate the given points using a parametric cubic spline then we compute the curvature of the interpolating spline. This last method is much less accurate and much more sensitive to noise than methods a. and b.

Spreading phase using a Tanner type law

We start by testing the efficiency of the algorithm. To this end, we consider the Tanner type law as triple line dynamic modelling, and we make fit this law with available experimental data, [START_REF] Allaman | Impact and spreading of micro-drops onto solid substrates[END_REF], related to a spreading phase. The triple point position is imposed at each time. We focus on the volume conservation, the height and diameter of the spreading (splat radius) and the deformation of the mesh (number of remeshing necessary), see Fig. 9. With a time step dt = 5.10 -5 s and a coarse mesh (≈ 600 elements), volume lost after 650 iterations is roughly 3%. The loss occurs mainly at the very beginning of spreading.

Mesh transport is efficient hence re-meshing occurs only when a projection occurs (see previous section).

From a qualitative point of view, numerical results are similar to experimental data, [START_REF] Allaman | Impact and spreading of micro-drops onto solid substrates[END_REF]. We present some influence of terms specific to the Shikmurzaev's model by considering the following steady-state Stokes system: Find ( u, p) such that:

2 Ω D(u) : D(v)dx - Ω p div(v)dx + 2 Ω u a vdx + β Γ slip u r v r ds = Γ f ree κ Ca v nds + Γ f h v τ ds - Γ sl g r v r ds, ∀v (17)
Ω div(u) q dx = 0, ∀q

We try to find the main effects of parameters β, g, h from numerical simulations. To do so, we fix one of the three parameters and we simulate the Stokes system (17)- [START_REF] Qian | Molecular hydrodynamics of the moving contact line in two-phase immiscible flows[END_REF] above with different values for the remaining ones. The set of absolute values used for the parameters is {0, 10 2 , 10 4 } with signs taken as previously discribed. According to Fig. 10 the effect of high β values is similar to a no-slip boundary condition on Γ sl . For the g-term, it is clear that a high |g| (g < 0) is equivalent to accelerating the spreading; in order to get the fewer constraint, the tangential velocity have to be of the order of -g β .

The h-term creates an upwind force which should be responsible for recoiling, see Fig. 11 and12, for β = 0, 10 4 and different values of g and h. Also, because of continuity, Fig. 10 and 13 suggest the that caterpillar-like motion is driven by solid-liquid surface properties as well as by those of the liquid-gas interface. In addition, according to Fig. 12 and 13, -g β acts differently on the free interface whether the contact angle is greater or fewer than π 2 .

Spreading phase using the Shikhmurzaev model

We consider the spreading phase using the Shikhmurzaev's model and the algorithm presented previously.

As mentioned above, g acts like a control on the position of the triple point. A decrease for g < 0 translates in a faster spreading. If g is small, this influence is likely not to be seen in the first part of the spreading, when the inertial forces dominate all other forces in presence. By increasing |g| over a threshold, which depends on β, it accelerates the spreading and modifies the flow nature. For |g| large, large interface distortions occurs and mesh generator fails to re-mesh. In order to observe the influence of the g-term,

we fix β = 1000, and we perform 1000 iterations with a step size of dt = 5.0 × 10 -6 for g ∈ {-10, -100}, Fig. 14. One can observe that after the inertial phase, the triple point position is farther to the right with an increasing |g| as shown more clearly by Fig. 15.

Fig. 16 shows a droplet subjected to large distortion (similar to breakup) when |g| is getting larger.

Simulation can go no further because of the change of topology it implies is not supported by the ALE method we are using.

We try now to see the influence of the sign of parameter g on the droplet behavior. Since we have seen that parameter g monitor the spreading, it is also supposed to do so with recoiling. In Fig. 17, simulation is started with initially spread droplet. Also here, the ratio -g β significantly influences the behavior of the droplet while parameter h influences more particulary on the droplet curvature than on the triple-point position.

Conclusion

We have discretized and implemented macroscopic part of the Shikhmurzaev model, [START_REF] Shikhmurzaev | The moving contact line problem on a smooth solid surface[END_REF], [START_REF] Blake | Experimental evidence of nonlocal hydrodynamic influence on the dynamic contact angle[END_REF], using an ALE formulation and a finite element scheme. Also, a connection with Qian-Wang-Sheng's results, [START_REF] Qian | Generalized Navier boundary condition for the moving contact line[END_REF],

has

SOLIDFigure 1 :

 1 Figure 1: Advancing angle and Shikhmurzaev's model

Figure 2 :

 2 Figure 2: Algorithm of resolution

Figure 3 :

 3 Figure 3: Spreading process. Circles represent no-slip nodes, triangles represent slip nodes uz = 0, squares represent "free" nodes, diamond is the triple point.

Figure 4 :

 4 Figure 4: Recoiling process. Circles represent no-slip nodes u = 0, triangles represent slip nodes uz = 0, squares represent "free" nodes, diamond is the triple point.

Figure 5 :

 5 Figure 5: Inner point X i . Local least square approximation using Bezier's curve.

Curve 1 . 6 , 7 ,

 167 a circle, Fig. Curve 2. an "oscillating-curve", Fig. Curve 3. a "double circle', Fig. 8.

Figure 6 :

 6 Figure 6: Curve 1. with noise introduced and computed curvature value

Figure 7 :Figure 8 :

 78 Figure 7: Curve 2. with noise introduced and computed curvature value

Figure 9 :

 9 Figure 9: Spreading: droplet profiles and volume in function of time steps

  been presented. The present model comprises generalized Navier slip boundary conditions with local surface tension gradient terms in the vicinity of the contact line. The local term values have been chosen upon a previous mathematical and numerical study done in[START_REF] Monnier | Analysis of a Local Hydrodynamic Model with Marangoni Effect[END_REF]. To test the efficiency of the present algorithm, first we considered a classical mobility relation and we obtained spreading droplet profiles qualitatively comparable to experimental results. Then we showed how the introduction of these new local terms in the vicinity of the contact line allow to model the spreading phase. Concerning the recoiling phase, further investigation to implement properly our algorithm is under progress. Also, a full coupling between the macroscopic part and the microscopic part of the Shikhmurzaev model is under progress.

Figure 13 : 18 0Figure 14 :

 131814 Figure 13: Case 3: influence of β, g and h, h = 2g. Velocity in the vicinity of the triple point for:(β, g) = {(0, 0), (10 4 , 0)} (up), {(0, -10 2 ), (10 4 , -10 2 )} (middle), {(0, -10 4 ), (10 4 , -10 4 )} (down). Scale factors are 0.05, 0.05, 0.005, 0.1, 5 × 10 -5 , 0.005, respectively. θ > π 2

Figure 15 :

 15 Figure 15: Triple-point porition for (β, g) = {(10 3 , -10), (10 3 , -10 2 )}

Figure 16 :Figure 17 :

 1617 Figure 16: Large interface distortion, similar to a droplet breakup g -10 4 }

Acknowledgements. The authors would like to thank A. Soucemarianadin and J.-L. Achard from Laboratory of Geophysical and Industrial Fluid Flows -LEGI-Grenoble, France, for their numerous and fruitful discussions.


D Jalocha
email: jalocha@lms.polytechnique.fr

• A Constantinescu

• R Neviere

A Constantinescu
email: andrei.constantinescu@lms.polytechnique.fr

Prestrain-dependent viscosity of a highly filled elastomer: experiments and modeling

Keywords: Prestrain, Time spectrum of relaxation, Nonlinear viscoelasticity, Constitutive behavior B D. Jalocha

Highly filled elastomers exhibit a complex microstructure made up of rigid fillers bounded by a thin layer polymeric matrix. The interactions between the fillers and the binder amplify locally the applied strains and induce a nonlinear viscoelastic behavior. The aim here is to analyze the influence of prestrain on the viscoelastic behavior. This paper proposes a prestrain-dependent viscoelastic constitutive model. The model is a superposition of three relaxation spectra, each corresponding to a family of polymer chains, and can be regarded in either its continuous or discrete expression. More specifically, one of these relaxation spectra is modified to assure the prestrain sensitivity. The parameters of the discrete model are identified from relaxation and DMA experiments performed on a solid propellant, and the obtained predictions match closely the experiments. The novelty of the analysis proposed in this paper is threefold. On the one hand, we report a new series of experimental measures, performed for a large range of frequencies for the DMA experiment and relaxation times for the relaxation experiment, and, on the other hand, we propose a constitutive law compatible with the principles of thermodynamics, which predicts closely the measurements. Finally, the analysis is performed comparing both relaxation and DMA experiments using the spectrum of relaxation times. A peculiarity of the present discussion is the novel identification method used, which identifies directly the relaxation times. This technique leads to models with smaller and optimum numbers of parameters than classical methods based on a logarithmic distribution of relaxation times.

Introduction

Composite materials are employed in various applications, in which they exploit their exceptional properties. Carbon fiber composites appreciated for the excellent stiffness-to-weight ratio are increasingly applied in the aerospace industry [START_REF] Botelho | Mechanical behavior of carbon fiber reinforced polyamide composites[END_REF]. Reinforced rubbers are solicited in large elastic strain regime with damping and are recommended for the manufacturing of tires, sealing, etc. [START_REF] Roeder | Elastomeric bearing: State of the art[END_REF]. Solid propellants, another class of reinforced rubbers, are retained for their energetic properties in solid boosters [START_REF] Ozupek | Constitutive Modeling of Hight Elongation Solid Propellants[END_REF]. The mechanical behavior of these different composites is generally viscoelastic and tailored with respect to the final application. Therefore, the modeling of the material behavior at different scales took an important place in the last decades (see, e.g., [START_REF] Huber | Finite deformation viscoelasticity laws[END_REF][START_REF] Christensen | A nonlinear theory of viscoelasticity for application to elastomers[END_REF][START_REF] Swanson | A constitutive formulation for high elongation propellants[END_REF].

The models describing the mechanical viscoelastic behavior of reinforced elastomers or rubbers are of different nature. Some studies propose hyperelastic constitutive models for filled elastomer [START_REF] Arruda | A three-dimensional model for the large stretch behavior of rubber elastic materials[END_REF]Oscar 2010), and others reproduce the timedependent behavior of the material based on phenomenological consideration [START_REF] Ozupek | Constitutive modeling of high-elongation solid propellants[END_REF][START_REF] Simo | On a fully three-dimensional finite-strain viscoelastic damage model: Formulation and computational aspects[END_REF][START_REF] Steinmann | Hyperelastic models for rubber-like materials: Consistent tangent operators and suitability for Treloar's data[END_REF] or on homogenization theory [START_REF] Xu | Constitutive modeling of solid propellant materials with evolving microstructural damage[END_REF]. Specifically, nonlinear viscoelastic models have recently been explored by combining the physics at the scale of the microstructure with the mechanical observation at the macroscopic scale of the applications (Azoug et al. 2014a(Azoug et al. , 2014b(Azoug et al. , 2014c)). Other studies [START_REF] Lion | The Payne effect in finite viscoelasticity: Constitutive modelling based on fractional derivatives and intrinsic time scales[END_REF] proposed complex constitutive models in the large-strain regime taking into account the nonlinear effects due to the presence of high-volume fraction of fillers.

Reinforced polymers and rubbers will cover a large range of the volumetric filler fraction, lying between 20 % and 70 %. Small filler fractions are used in rubbers, whereas high filler fractions are found in propellants. We remark that generally the nonlinearity of the mechanical behavior increases with increasing filler fraction.

The aim of this paper is to explore the particular nonlinearity induced by a prestrain on the viscoelastic properties of a highly filled elastomer. This effect was observed during an experimental campaign on solid propellants, performed by [START_REF] Azoug | Micromécanismes et comportement macroscopique d'un élastomère fortement chargement[END_REF] and reported in several articles. In [START_REF] Thorin | Influence of prestrain on mechanical properties of highly-filled elastomers: Measurements and modeling[END_REF] the influence of prestrain on viscoelastic properties for different HTPB composite propellants was measured using dynamic mechanical analysis (DMA). The prestrain-dependent behavior was modeled using a modified generalized Maxwell model, where the stiffnesses of each Maxwell element was linearly dependent on prestrain. The paper (Azoug et al. 2013b) was devoted to a biaxial experiment having a prestrained DMA on one axis and a fixed prestrain on the other. The experiments showed that the storage and loss moduli increase with the prestrain under both loadings and that the nonlinear behavior is quantitatively modified by adding an orthogonal prestrain. Moreover, the modification of the behavior under a horizontal prestrain was canceled out by an increase of the vertical prestrain, a phenomenon that may be explained by the alignment of the fillers with the prestrain. The experimental results of a complete biaxial prestrained DMA experiment have been recently reported in Jalocha et al. (2015b). The experiments on propellant specimen confirmed the previous results. Moreover, it has been shown that the dependence of prestrain in biaxial experiments can be characterized using the second invariant of the prestrain. The limitation of the reported analysis is the unique frequencies where the viscoelastic nonlinearity has been measured. The novelty of the analysis proposed in this paper is threefold. On the one hand, we report a new series of experimental measures, performed for a large range of frequencies for the DMA experiment and large range of times for the relaxation experiment. On the other hand, we propose a constitutive law, compatible with the principles of thermodynamics, which predicts closely the measurements. Finally, the analysis is performed comparing both relaxation and DMA experiments.

The analysis of the data starts from the idea that relaxation spectrum and generalized Maxwell models represent continuous and approximate discrete viscoelastic models, respectively. As a consequence, the variations of the experimental observations as functions of the prestrain are considered as variations of the measured complex and relaxation moduli. In this paper, we propose a prestrain-dependent viscoelastic constitutive model. The model is a superposition of three relaxation spectra, each corresponding to a family of polymer chains, and can be regarded in either its continuous or discrete expression. More specifically, one of these relaxation spectra is modified to assure the prestrain sensitivity. The parameters of the discrete model are identified from relaxation and DMA experiments performed on a solid propellant, and the obtained predictions match closely the experimental. A particularity of the present discussion is the novel identification method used, which identifies directly the relaxation times. This technique leads to models with smaller and optimum numbers of parameters than classical methods based on a logarithmic distribution of relaxation times. A detailed mathematical description of this identification algorithm and a discussion of its robustness are reported in Jalocha et al. (2015a).

The paper is organized as follows. The first section presents the material and experiments. The second and third sections discuss the continuous and discrete models, respectively. The final sections present the results and conclusions.

Material and experiments

Material

The material studied here is a solid propellant, used for the propulsion of rockets and classified as a highly filled elastomer. Fillers represent 70-80 % of the total volume and are bonded by a viscoelastic matrix [START_REF] Azoug | Micromécanismes et comportement macroscopique d'un élastomère fortement chargement[END_REF], see Fig. 1(a). The performance of the rocket engine will depend on both energetic capacity of the filler and the mechanical properties of the propellant in order to guarantee the performance and integrity of the engine. The mechanism of the engine can be summarized as follows. The oxidation-reduction reaction decomposes the fillers and produces highly pressurized gases in the combustion chamber, and the gases are ejected through the nozzle of the engine and propel the rocket. Integrity of the engine imposes that the evolution of the propellant during the burning process stays stable and in accordance with the initial design in spite of the important mechanical loadings of the propulsion.

The high-volume fraction of fillers demands exceptional mechanical properties from the matrix that is a network of thin bridges surrounding the rigid fillers. The characteristic of the polymer matrix assured by its microscopic chemical and physical properties. The manufacturing process will generally create an over cross linked elastomer matrix. At the microscopic scale, one will encounter a multiphase binder as already discussed in [START_REF] Azoug | Micromécanismes et comportement macroscopique d'un élastomère fortement chargement[END_REF] and composed of: (i) a principal cross linked network polymer that bonds the solid particles together and (ii) nonlinked polymer chains floating freely through the principal network. An idealized structure is displayed in a schematic diagram in Fig. 1(b).

Experimental procedure

The viscoelastic behavior of the propellant at the macroscopic scale of the structure can be characterized using standard relaxation or DMA experiments. Let us next introduce some The viscoelastic behavior of a material is completely represented by the continuous spectrum of relaxation time H (τ ) as exposed, for example, in [START_REF] Findley | Creep and Relaxation of Nonlinear Viscoelastic Materials[END_REF]. The spectrum H (τ ) is related to the relaxation modulus E(t) over the time domain and to the dynamical modulus E * (ω) = E (ω) + iE (ω) over the frequency domain by the following equations:

E(t) = E 0 + ∞ -∞ H (τ ) e -t τ d ln(τ ), E * (ω) = E 0 + ∞ -∞ H (τ ) iωτ 1 + iωτ d ln(τ ). (1)
According to [START_REF] Markovitz | Boltzmann and the beginnings of rheology[END_REF], the relaxation and complex moduli permit to establish the classical viscoelastic constitutive equations, expressed in the time domain as

σ (t) = t -∞ E(t -τ ) dε(τ ) dτ dτ (2)
and in the frequency domain as

σ * (ω) = E * (ω)ε * (ω), (3)
where σ and ε denote the stress and strain fields, respectively. The representations in (1) and (3) are equivalent, and the passage between the time and frequency domains is obtained through a direct or inverse Fourier transform. As a consequence, σ * is the direct Fourier transform of the stress σ , ε * is the direct Fourier transform of the strain ε, and the complex modulus E * is the direct Fourier transform of the relaxation modulus E.

Let us further recall that the spectrum of relaxation times H (τ ) cannot be measured directly. It is only recovered from either the relaxation modulus E(t) or the complex modulus E * (ω), which are the direct outputs of the relaxation and DMA experiments; for a detailed explanation, see [START_REF] Knauss | Handbook of Experimental Solid Mechanics[END_REF]. As a consequence, we can conclude that the relaxation and the DMA experiments are dual to each other, one representing the time domain, and the other the frequency domain.

The time or frequency domain, covered by the relaxation or the complex modulus, respectively, tends to cover several orders of magnitude. Hence, it is difficult to cover experimentally the complete viscoelastic domain, and based on the time-temperature superposition principle, the so-called master curves of relaxation and complex modulus are constructed. The method is explained, for example, in [START_REF] Brinson | Polymer Engineering Science and Viscoelasticity. An Introduction[END_REF], [START_REF] Williams | The temperature dependence of relaxation mechanisms in amorphous polymers and other glass-forming liquids[END_REF].

The objective here is to propose a viscoelastic constitutive model with prestrain sensitivity capable of predicting data from relaxation and DMA experiments. The main features of the experimental procedure are summarized in Fig. 2. In both cases, strains are imposed, stresses are measured, and the moduli are computed from the measured data. The data will be first collected in the format of a master curve encompassing several orders of magnitudes in the time or frequency range and used after for identification of the model. The prestrain is applied in both experiments as an initial load step as represented schematically in Fig. 3. The mathematical duality of the measured data in the sense of the Fourier transform discussed implies that the effect of the applied prestrain should be represented by the same phenomena to cover both loading cases.

In the time domain, the relaxation modulus E(t) is measured in the following way: a tensile prestrain ε s is imposed, and then, after a short delay of 10 min considered large enough for the stress to reach a stabilized state, the relaxation test is started. The history of the imposed strain is displayed in Fig. 3(a). The relaxation test is performed imposing a tensile strain step superimposed with the prestrain: ε(t) = ε s + ε 0 h(t) (where h(t) is the Heaviside function), and the decreasing stress σ (t) is recording in function of the time. The relaxation modulus is computed from the data as the secant modulus defined by E(t) = σ (t) ε 0 , where σ (t) is the difference between the relaxed stress and the residual stress as a function of the imposed prestrain.

In the frequency domain, the complex modulus E * (ω) is measured in the following steps: a tensile prestrain ε s is imposed, and, after a short delay of 10 min, the DMA test is started. The history of the imposed strain is displayed in Fig. 3(b). The DMA test is performed imposing sinusoidal strain around the given prestrain, ε(t) = ε s + ε d sin(ωt), and the stress σ (t) is recorded. The complex modulus is computed as the secant modulus using the formula

E * (ω) = σ d ε d (cos δ + i sin δ)
, where σ d is the amplitude of the sinusoidal stress, ε d is the amplitude of the sinusoidal strain, and δ is the time delay between the strain and stress signals measured during the stabilized period.

As mentioned before, both moduli are computed as secant moduli. This guarantees the consistency and correlation in the data processing between the two experiments.

In experiments involving viscoelastic filled materials under relaxation or DMA tests, one encounters the Payne effect [START_REF] Payne | The dynamic properties of carbon black-loaded natural rubber vulcanizates. Part I[END_REF] or, close to the large-strain region, the Mullins effect [START_REF] Mullins | Softening of rubber by deformation[END_REF]). In the time domain, the Mullins effect characterizes an instantaneous softening of the stress-strain curve that occurs whenever the load increases beyond its prior all-time maximum value. For the propellant under scrutiny, we remark that a small softening of the Mullins type is generally observed. However, the effect is completely erased if a sufficiently long time lapse occurs between the two loadings. In the frequency domain, the Payne effect characterizes a softening of the norm of the complex modulus, which increases the applied dynamic strain amplitude. The effect is present in the material and was the reason to chose only small dynamic strain amplitudes. However, the chosen amplitude is smaller than that encountered in practical applications. Let us further remark that the two effects are generally considered to be physically closely related and could be at the origin of the small difference in viscoelastic properties measured from the relaxation and the DMA experiments. Except this fact, the process history of the material does not change the experimental results.

Mathematically, it is straightforward that the essential features of the model are expressed as the spectrum of relaxation times H (τ ), and the aim is to identify the spectrum from both relaxation and DMA experiments. Therefore, DMA results are transposed into the time domain using the inverse Fourier transform based on the empirical relation

E(t) = E * ( t 2π ) .

Experimental results

The relaxation tests have been performed for the following prestrain values: ε s = 0; 1.2; 2.5; 3.7; 5; 6; 7.2; 8.4 % at temperatures between +20 °C and +40 °C. The choice of the maximum prestrain value of 8.4 % ensures that the experiments are performed in the incompressible domain of the material and that damage (i.e., filler debonding, matrix fracture, etc.) does not occur. The temperature range is coherent with the application constraints of the model. Moreover, results given at smaller temperature than the lower limit of +20 °C can later be replaced with DMA experiments, which provide data of better accuracy. The upper limit +40 °C covers, after the time-temperature superimposition, the application of the model.

The master curves obtained at a reference temperature of T ref = 20 °C for different prestrain values are presented in Fig. 4. Panel (a) presents the time evolution of the relaxation modulus E(t), and panel (b) displays the instantaneous values of the relaxation modulus at t = 0.03 s. Both representations show that the relaxation modulus increases nonlinearly as a function of the prestrain ε s .

The DMA tests are performed for the values of prestrain ε s = 0; 3; 4; 6 % at temperatures between -90 °C to +90 °C with an increment of 5 °C. The imposed strain amplitude ε d for the DMA test was equal to 0.02 %, and, as a consequence, the ratio between the cyclic varying strain and the static strain load was in minimum 2 × 10 -3 .

The master curves obtained from the DMA data for the reference temperature of T ref = 20 °C are presented in Fig. 5. Panel (a) represents the evolution of the norm of the amplitude of the complex modulus E * (ω) for increasing pulsation ω and different prestrains. One can remark that the influence of the prestrains occurs essentially at small pulsation. Moreover, a detailed analysis of the evolution of E * (ω) with increasing prestrains, at fixed pulsation ω = 31.4 rad/s, displayed in panel (b), exhibits a twofold increase of the amplitude of the complex moduli accompanied with a 30 % decrease of tan(δ) over the tested prestrain range. The discrepancy for low frequencies, smaller than 10 5 rad/s, is the direct effect of the prestrain. The different curves show that the amplitude of the complex modulus increases with respect to the prestrain value. As in the case of the relaxation modulus, the complex modulus amplitude increases as a function of prestrain ε s , but the phase angle δ decreases nonlinearly as a function of the same prestrain; see Fig. 5

(b).

A complete master curve for the relaxation modulus relaxation in the absence of prestrain ε s = 0 % is exhibited in Fig. 6. The curve is an overlap of the relaxation data, which cover the large times, and the inverse Fourier transform of the DMA data, which cover small times. It will cover a total time range between 10 -15 s and 10 8 s for a reference temperature of 20 °C. In panel (a) of Fig. 6, we displayed the relaxation curve, and one can observe a good match between the two data sets, in spite of a small slope change captured by the DMA relaxation data. Panel (b) of Fig. 6 exhibits a comparison between DMA and relaxation data in the evolution of the relaxation modulus with prestrain at t = 0.03 s. The data sets have a similar nonlinear evolution with a relative error of ≈10 % for small prestrains, which grows to ≈15 % for large prestrains.

Continuous model

Material without prestrain

Physically, for polymer material, the relaxation spectrum H (τ ) is related to the statistical distribution of the molecular mobility of polymer chains [START_REF] Doi | Molecular theory of the viscoelastic properties of concentrated polymer solutions[END_REF]. At the microscopic scale, the molecular mobility can now be measured by nuclear magnetic resonance (NMR) spectrometry to investigate the spin relaxation times of the chains. Moreover, the results from NMR published in [START_REF] Mowery | Sensitivity of proton NMR relaxation times in a HTPB based polyurethane elastomer to thermo-oxidative aging[END_REF], [START_REF] Litvinov | EPDM-carbon black interactions and the reinforcement mechanisms, as studied by low-resolution 1H NMR[END_REF] prove the existence of a multiphase matrix in reinforced polymers and explore the influence of the binder matrix interactions on the mechanical properties.

At the macroscopic scale, this spectrum can be computed directly from the complete master curve of the relaxation modulus using the inverse Laplace transform; see Eq. (1). Following [START_REF] Baumgaertel | Interrelation between continuous and discrete relaxation time spectra[END_REF], a linear relaxation modulus curve can be generally expressed by the following analytical expression:

E(t) = E 0 + E g -E 0 (1 + t t 0 ) β , ( 4
)
where E 0 , E g , t 0 , and β are adjustable parameters. Starting from the two remarks and ap-

plying the change of variable f (z) = H ( 1 z ) z
in (1), the analytical expression of the spectrum H (τ ) is given by the inverse Laplace transform of (4):

H (τ ) = (E g -E 0 ) Γ (β) τ t 0 -β e -t 0 τ (5)
with Γ (β) the Gamma function of β. See [START_REF] Smith | Empirical equations for representing viscoelastic functions and for deriving spectra[END_REF] and [START_REF] Widder | The Laplace Transform[END_REF] for a complete explanation of the operations performed.

An inspection of the curve of relaxation modulus of solid propellant without prestrain displayed in Fig. 6(a) shows that the relaxation modulus is not linear. Based on NMR measurements, Azoug et al. (2014c[START_REF] Azoug | Microstructure and deformation mechanisms of a solid propellant using 1H NMR spectroscopy[END_REF] showed that propellant has a multiphased matrix due to his manufacturing process. Moreover, they were able to bring up three phases in the matrix of the propellant associated to families of polymer chains and mechanisms of mobility to explain the observed nonlinearity in the relaxation. [START_REF] Zimm | Reptation of a polymer chain in an irregular matrix: Diffusion and electrophoresis[END_REF] proposed to associate mechanism of molecular mobility to phases of the polymer matrix. As a consequence of the previous observation, we propose to associate to each a phase of the matrix of the propellant and correspondingly to each chain mobility mechanism a relaxation modulus of the form (4) and to sum up the contributions in the relaxation modulus:

E(t) = E 1 (t) + E 2 (t) + E 3 (t). ( 6
)
As a consequence, the modulus is approximated by the following analytic expression:

E(t) = E g1 (1 + t t 1 ) β 1 + E g2 (1 + t t 2 ) β 2 + E g3 (1 + t t 3 ) β 3 . ( 7
)
From a physical point of view, as explained in Azoug et al. (2013a), this decomposition is coherent with the assumption of the existence of three families of different molecular mobilities in the binder:

• The first family (subscript "1") represents the molecular mobility of the principal reticulated network polymer of the matrix. • The second family (subscript "2") represents the molecular mobility of the free polymer chains through the principal network of the matrix. • The third family (subscript "3") represents the long-time polymer chain mobility.

As the Laplace transform is a linear operator, from (7) it follows that the relaxation spectrum is the sum of the contributions for each phase H (τ ) = H 1 (τ ) + H 2 (τ ) + H 3 (τ ), or explicitly:

H (τ ) = E g1 Γ (β 1 ) τ t 1 -β 1 e -t 1 τ + E g2 Γ (β 2 ) τ t 2 -β 2 e -t 2 τ + E g3 Γ (β 3 ) τ t 3 -β 3 e -t 3 τ . ( 8
)
3.2 Prestrained material [START_REF] Simhambhatla | On the rheological modeling of filled polymers with particle matrix interactions[END_REF] assumed that the prestrain influences only the mobility of the free polymer chain, and the change induced into the fillers network will only affect the elastic part of the behavior, not the viscous part. Physically, one can visualize the effect of prestrain as a trapping mechanism of the fillers and the reticulated polymer network on the free chains, which have a growing restriction in the movement with increasing prestrain. Moreover, the NMR measurements on prestrained samples, reported in [START_REF] Azoug | Microstructure and deformation mechanisms of a solid propellant using 1H NMR spectroscopy[END_REF], show that the prestrain influences only the second family corresponding to free polymer chains, confirming the assumptions. This means that the prestrain affects the slope of the second contribution to the relaxation modulus, which implies that β 2 becomes a function of the prestrain ε s . The prestrain-modified spectrum H takes the following expression: 7) master curve of the relaxation modulus and the evolution of the relaxation moduli of the three phases

H (τ, ε s ) = E g1 Γ (β 1 ) τ t 1 -β 1 e -t 1 τ + E g2 Γ (β 2 (ε s )) τ t 2 -β 2 (εs ) e -t 2 τ + E g3 Γ (β 3 ) τ t 3 -β 3 e -t 3 τ . (9)

Identification

The parameters of the three families, that is, (E gi , t i , β i ) 1≤i≤3 , are now identified from the experimental relaxation modulus curves in two steps: (i) identification of the parameters in the absence of prestress ε s = 0, that is, the parameter β 2 is considered constant, and (ii) identification of the prestress-dependent parameter β 2 (ε s ). Both steps are performed by minimizing the least square error between the model using the NonLinearModelFit operator of MATHEMATICA ® (Mathematica).

The identified parameter values in the absence of prestress ε s = 0 are summarized in Table 1, and the comparison between prediction and measurements is plotted in Fig. 7. The model matches measurements over the total range of time, with the exception of the first decades, where a large error is observed. A closer inspection of plot and of the identified parameters shows that viscoelasticity is predominately driven by the molecular mobility of the first phase, representing the reticulated network. Nevertheless, in the time interval 10 -4 ≤ t ≤ 10 s, one can remark that the second phase of free polymer chains has a significant influence, even if its amplitude represents only 1 % of the amplitude of the principal reticulated network ( E g2 E g1 ≈ 0.01). The contribution of the third phase is very small when compared with the two others. However, it describes the viscoelastic behavior at long times.

The identified values of the model are now used to display in Fig. 8 the continuous spectrum of relaxation times and the contribution of the three phases of the polymer matrix. The phases representing molecular mobility will each contribute to a local maximum situation times τ of 1 × 10 -11 s, 1 × 10 -4 s, and 1.8 × 1 s, respectively. As expected, one can remark a similar trend with the spectra of spin-spin relaxation times of the chains obtained by NMR spectrometry.

The prestrain dependence of the complete behavior is embodied only in the parameter β 2 as mentioned earlier. The identification is performed using experimental data exhibited in Fig. 8 Identified spectrum of relaxation times of solid propellant H (τ ) (see ( 8)) with a focus on the time between 10 -6 ≤ τ ≤ 10 6 in order to observe the second local maximum.

β 2 (ε s ) = β 0 2 (1 + ( εs εr ) q ) . ( 10
)
The nonlinear regression leads to β 0 2 = β 2 (0) = 0.312 (see Table 1), ε r = 0.045, and q = 2.1. Figure 9 presents the comparison of the prediction of the model and experimental data, represented as a continuous line and dots, respectively.

The complete identification of the prestrain-dependent model permits to inspect the evolution of the spectrum of relaxation as a function of time and prestrain, that is, H (τ, ε s ). Figure 10 exhibits the time evolution of H (τ, ε s ) for the prestrain values ε s = 0, 2, 4, 6 %. The plot shows a finite support in time and a decreasing of the second protuberance with increasing prestrain. This protuberance represents the contribution of the second prestraindependent molecular mobility of free polymer chains. From the physical point of view, the decreasing contribution is the signature of the increasing trapping of the free chains with increasing prestrain.

Generalized Maxwell model with prestrain effect

Let us consider a viscoelastic body occupying the domain Ω. Using the standard notation in continuum mechanics, we shall denote the current point in the reference and the actual Fig. 10 Evolution of the spectrum of relaxation time between 10 -6 ≤ t ≤ 10 s for prestrain ε s = 0, 2, 4, 6 %. The second protuberance is decreasing with increasing prestrain configuration as X and x, respectively. The transformation gradient F is defined as F = ∂x ∂X . Under the assumption of small strain, the transformation gradient becomes F = I+Grad(u), where u = x -X stands for the field of displacements. The small strain tensor ε is computed as ε = 1 2 (Grad(u) + Grad T (u)). The viscoelastic constitutive law is presented in the framework proposed initially in [START_REF] Halphen | Sur les matériaux standards généralisés[END_REF], [START_REF] Holzapfel | Non Linear Solid Mechanics: A Continuum Approach for Engineering[END_REF] and denoted as generalized standard materials. This formal framework is thermodynamical and assures the complete consistency with the principles of thermodynamics and conducts naturally to efficient and robust numerical integration algorithms. The constitutive law is defined in terms of the free energy ψ and the dissipation potential φ. For the viscoelastic bodies studied here, we consider a decoupled representation of the free energy ψ:

ψ ε, α 1 , . . . , α m = ψ ∞ (ε) + n i=1 ψ i ε, α i , ( 11
)
where ψ ∞ and ψ i with 1 ≤ i ≤ n represent the elastic and viscoelastic parts of the free energy, and α i denote a series of internal variables; their physical significance will be clarified later in the presentation.

In this framework, according to [START_REF] Coleman | Thermodynamics with internal state variables[END_REF], the total stress σ is defined as

σ = ∂ψ ∂ε = ∂ψ ∞ ∂ε + n i=1 ∂ψ i ∂ε . ( 12
)
The formal construction assures that for any isothermal process, the Clausius-Duhem inequality

0 ≤ σ : ε - ψ is satisfied if ∂ψ i ∂α i + ∂φ i ∂ αi = 0 (13)
with φ i the dissipation potentials. A detailed discussion of this point for the case of viscoelasticity is given, for example, in [START_REF] Letallec | Numerical models of steady rolling for non linear viscoelastic structures in finite deformations[END_REF].

In order to specify detailed expressions of the free energy and the dissipation potential, we make the following choice of a quadratic expression for the elastic part for the free energy:

ψ ∞ (ε) = 1 2 E ∞ ε : ε, ψ i ε, α i = 1 2 E i ε -α i : ε -α i , (14)
with 0 ≤ E ∞ and 0 ≤ E i . This choice leads to a classical model in linear viscoelasticity and is also presented, for example, in [START_REF] Simo | On a fully three-dimensional finite-strain viscoelastic damage model: Formulation and computational aspects[END_REF].

If we further assume that the internal variables α i represent the strain in the viscous dampers, then the fictitious stress q i acting on each damper will be associated to the internal variable α i . Similarly to the reasoning in [START_REF] Holzapfel | Non Linear Solid Mechanics: A Continuum Approach for Engineering[END_REF], by analogy with (12) we postulate that q i = -∂ψ ∂α i . Therefore, applying ( 14) leads to

q i = E i ε -α i , 1 ≤ i ≤ n. ( 15
)
The dissipation potentials are equally chosen as quadratic potentials:

φ i = 1 2 η i αi : αi , 1 ≤ i ≤ n, (16)
where the viscosity η i is a positive real number. The construction of the constitutive model does not restrict so far a strain-dependent viscosity parameter. A similar construction and extended discussion of the subject are presented in [START_REF] Amin | Nonlinear dependence of viscosity in modeling the rate-dependent response of natural and high damping rubbers in compression and shear: Experimental identification and numerical verification[END_REF]. As a consequence, we define the constant relaxation time τ i as τ i = η i E i with η i the viscosity value at the origin of the process.

A series of simple algebraic transformations of Eqs. ( 12) and ( 14) lead to the following expression of stresses:

σ = E ∞ ε + n i=1 E i ε -α i . ( 17
)
Moreover, ( 13) and ( 16) lead to the following system of evolution equations for the internal variables α i :

η i αi -E i ε -α i = 0. ( 18
)
Complete equations of the constitutive model are obtained after a final inference, which starts from the assumption q ∞ = E ∞ ε, then uses the definition of the fictitious stress q i , and finally replaces the expression into ( 17) and ( 18). The final expression of the constitutive equations in terms of stresses and internal variables is

σ = q ∞ + n i=1 q i , ( 19
)
η i ε = τ i qi + q i . ( 20
)
The viscoelastic constitutive model described previously is denoted as a generalized Maxwell model. Moreover, it is consistent with the principles of thermodynamics by construction. It can also be described in terms of rheological elements by the so-called Prony The n viscoelastic elements are completely described by the characteristic parameters (τ i , η i ) 1≤i≤n , where τ i is the discrete relaxation time, and η i is the associated viscosity. It is generally schematically represented by simple rheological elements using springs and dampers; see, for example, [START_REF] Findley | Creep and Relaxation of Nonlinear Viscoelastic Materials[END_REF][START_REF] Simo | Computational Inelasticity[END_REF] and the graphic representation in Fig. 11(a).

For a homogeneous one-dimensional specimen, the model can be easily integrated numerically. The backward Euler method leads to the following time stepping algorithm, where the stress σ m at the time instant t m is computed as a function of the stress σ m-1 at the preceding time instant t m-1 and a load increment expressed as the total strain increment in the following steps:

(i) t = t m -t m-1 ; (ii) q i m = η i (εm-ε m-1 )+τ i q i m-1 τ i + t ; (iii) q ∞ m = E ∞ ε m ; (iv) σ m = q ∞ m + n i=1 q i m .
This algorithm can then be generalized into full 3D evolution equations and implemented in standard finite element codes.

Results and discussion

The viscoelastic relaxation spectrum H (τ, ε s ) is a continuous function of the time τ already discussed. Moreover, the continuous spectrum accepts a discrete approximation in the form of the prestrain ε s dependent generalized Maxwell model, represented by the so-called Prony series. The identification and comparison of its predictions with experimental data will be discussed next.

Let us start with the continuous spectrum in the absence of prestrain H (τ, 0). The identification method of such a spectrum is performed in two steps: (i) the choice of the characteristic times of the viscoelastic elements and (ii) the determination of the viscosities.

The classical choice of characteristic times is an equidistant distribution on the logarithm scale as suggested, for example, in [START_REF] Baumgaertel | Interrelation between continuous and discrete relaxation time spectra[END_REF]. Here we use an optimization method minimizing the distance between the continuous spectrum and its discrete approximation. The mathematical background of the method and its application has extensively been discussed in Jalocha et al. (2015a). In order to cover the features in this paper, we propose a short review in the Appendix.

If we consider the support of the continuous spectrum defined by a cut of 0.005 of its maximal value, 0.005 H max ≤ H (τ ) ≤ H max , then the corresponding time interval is defined for τ as 10 -14 s ≤ τ ≤ 10 20 s. The demanded accuracy in the optimization process for a relaxation spectrum H with three protuberances, one for each phase of the matrix, has led to n = 48 discrete characteristic relaxation times τ i . The Prony series is completed by the values of the viscosity η i obtained by the formula

η i = H (τ i , 0)τ i ln(r i ) ( 21
)
with r i = τ i+1 τ i-1 . The n elastic moduli E i are computed as the ratio of viscosities and characteristic times:

E i = η i τ i . ( 22
)
To continue the identification method, we now introduce the prestrain dependence of the continuous spectrum H (τ, ε s ). The prestrain dependence is only carried by the viscosity η i of each element, and therefore, the prestrain-dependent viscosity naturally becomes

η i (ε s ) = H (τ i , ε s )τ i ln(r i ). ( 23
)
Let us finally remark that the construction preserves prestrain-independent characteristic times, whereas the viscosity is prestrain dependent. The expression of the thermodynamic potentials preserve a quadratic free energy in terms of strains and is coherent with the principles.

A schematic view of the model is presented in panel (a) of the Fig. 11, and panel (b) illustrates the identified prestrain viscosity evolution of the 17th and 28th elements of the identified model. One can remark that the elements have different evolutions, which is also a signature of different nonlinearities in different viscoelastic elements.

The predictions of the identified prestrain-dependent generalized Maxwell model will be compared next with experimental data from the relaxation or the DMA tests.

The first comparison concerns the relaxation test starting from several values of prestrain, with an imposed strain history as displayed in Fig. 3(a). The value of the strain step was ε 0 = 1 %.

The comparisons of model predictions and experimental data in terms of evolution of the relaxation modulus as a function of time are plotted in Fig. 12(a). The figure represents the time interval 0.03 ≤ t ≤ 1000 s, where the nonlinearity is strongest. The model matches both the evolution of the relaxation modulus as a function of the time and the order of magnitude imposed by the presence of different prestrains.

The second comparison concerns the DMA test. We only represent the evolution of the norm of the complex modulus and of the loss modulus as a function of prestrain at a given angular frequency for several prestrain values, with an imposed strain history as displayed in Fig. 3(b). The chosen pulsation and strain step values were ω = 31 rad/s and ε d = 0.02 %, respectively. The viscosity η i of the prestrain-dependent model is considered to follow the . 23). We further assumed that the prestrain ε s can be reasonably approximated by the total strain ε.

A comparison, in terms of the norm of the complex modulus as a function of the prestrain, between numerical predictions and the measured experimental data is plotted in Fig. 12(b). One can easily observe that the model predictions match the experimental data. Moreover, the experimental increase of the amplitude of the complex modulus as a function of the prestrain is closely described by the model. small panel in the same picture exhibits the capacity of the same model to correctly predict the decrease of the shift factor δ as a function of the prestrain. A close inspection of the values shows a small overestimation of the numerical prediction when compared with the experimental data.

Conclusion

We presented a series of relaxation and DMA experiments involving prestrain and proposed a thermodynamical coherent constitutive model to predict the measured data. The proposed model is a prestrain-dependent generalized Maxwell model. The prestrain dependence was carried by the viscosities of each viscoelastic element, whereas the characteristic relaxation times and elastic moduli stayed prestrain independent. The parameters of the model were then identified from experimental data with a robust identification method, and numerical predictions using a simple backward Euler integration algorithm allowed us to confront the model to the experiments.

The material under discussion was a high filled elastomer in the form of a solid propellant. It exhibited a high-volume fraction of the filler with respect to the matrix and, as a consequence, manifested a complex viscoelastic behavior.

The results show that the nonlinearity of the macroscopic material behavior exhibited during relaxation and DMA experiments could be matched by the predictions of the experiment. The nonlinearity discussed in this paper was the evolution of the viscoelastic properties as a function of the prestrain. The results presented here complete both in terms of experiments and modeling the initial results studied in [START_REF] Azoug | Micromécanismes et comportement macroscopique d'un élastomère fortement chargement[END_REF], Azoug et al. (2013b), [START_REF] Thorin | Influence of prestrain on mechanical properties of highly-filled elastomers: Measurements and modeling[END_REF].

The nonlinear evolution of viscoelastic behavior as a function of the prestrain is a direct physical signature of the mobility of the molecular chains in the polymeric binder. The experimental results discussed here confirm the NMR result obtained at a microscopic scale and reported in [START_REF] Azoug | Microstructure and deformation mechanisms of a solid propellant using 1H NMR spectroscopy[END_REF]. Moreover, they also show that the prestrain influences the relaxation modulus up to time t = 10 -4 s.

The proposed prestrain-dependent generalized Maxwell model predicts closely the experimentally observed nonlinear viscoelastic behavior in terms of (i) complex and loss moduli or (ii) relaxation moduli. The thermodynamic consistent construction of the model is compatible with two or three spatial dimensions and can therefore easily be extended for multi-axial loading configurations. An interesting perspective is the exploration of the prestrain-dependent biaxial experimental. A further challenge is the programming of a complete numerical integration algorithm.

The influence of the prestrain on the second phase, characterizing the molecular mobility of the free polymer chains, seems to be a characteristic of highly filled elastomer. Preliminary finite element computations on 3D microstructures have shown a tenfold amplification of the strain inside the matrix when compared with the applied macroscopic strain. This strain amplification diminishes the space between the fillers and therefore restrains the movements of the free polymer chains. A work in progress analyzes the evolution of this effect as a function of the filler volume fraction.

This papers opens new perspectives both in terms of experiments and theory for the nonlinear viscoelastic behavior of solid propellant and, more generally, for highly filled elastomer under multiaxial prestrain. The optimal distribution of the relaxation times τ i is defined in terms of the endpoints of each unit box function τ i = t i-1 +t i 2

, i = 1, . . . , n, for given n. Moreover, it will minimize the residual function R. The end points of the time series t 0 and t n are supposed to be a priori chosen. See Fig. 13 for a schematic view.

Fig. 1

 1 Fig. 1 Microstructure of solid propellant: (a) optic observation under a scale factor of 10 and (b) schematic representation of a cross section

Fig. 2

 2 Fig. 2 Schematic view of the experimental procedure

Fig. 3

 3 Fig. 3 Experimental setup to measure the relaxation modulus (a) and the complex modulus (b) under prestrain ε s

Fig. 4

 4 Fig. 4 (a) Relaxation modulus master curves under prestrain between 0 % ≤ ε s ≤ 8.4 % for a relaxation strain step of ε 0 = 1% and for a reference temperature of T ref = 20 °C. (b) Nonlinear evolution of the relaxation modulus in function of the prestrain ε s for a given time t = 0.03 s

Fig. 6

 6 Fig. 6 Complete classical master curve of relaxation modulus using relaxation and DMA data. (a) Evolution of the relaxation modulus for ε s = 0 % as a function of time from 10 -15 s to 10 8 s. (b) Evolution of the relaxation modulus as a function of prestrain ε s

Fig. 9

 9 Fig. 9 Evolution of β 2 in function of the prestrain: dots are identified from experimental results, a solid line is the proposed model (10) to catch this evolution

.

Fig. 11

 11 Fig. 11 (a) Schematic view of a generalized Maxwell model with variable, prestrain dependent viscosity. (b) Evolution of viscosity with prestrain in the 17th and 28th element of the identified model

Fig. 12

 12 Fig. 12 Comparison between the predictions of the prestrain dependent generalized Maxwell model and experimental results. (a) Evolution of the relaxation modulus as a function of time for different values of prestrain. (b) Evolution of the norm of complex modulus as a function of prestrain for ω = 31 rad/s and ε d = 0.02 %

Fig. 13

 13 Fig. 13 An example of continuous spectrum of relaxation time H (τ ) and a staircase function S(τ ). Computation principle of τ i

Table 1

 1 Identified parameter of the relaxation modulus in the absence of prestrain

		E gi [MPa]	t i [s]	β i [adim]
	E 1	12191	3.4 × 10 -12	0.333
	E 2	91	3.1 × 10 -5	0.312
	E 3	6	1 .27	0.07
	Fig. 7 A comparison of
	measured and predicted (

Herakles, Centre de Recherche du Bouchet, 9 rue Lavoisier, 91710 Vert Le Petit, France

Acknowledgements

The authors would like to thank the Direction Générale de l'Armement, especially Laurent Munier (DGA), and HERAKLES-SAFRAN for financial support for this study.

Appendix: Discretization of the continuous spectrum of relaxation time

The generalized Maxwell model, represented as a Prony series is a discrete approximation of the continuous spectrum of relaxation H (τ ). The standard discretization method starts with the discretization of the time axis into a finite set of instants (τ i ), separated by equidistant intervals on the logarithmic scale. Then, the associated elastic modulus is computed using the formula

) and leads to the complete Prony series (τ i , E i ), i = 1, . . . , n, with n the number of viscoelastic elements of the model. The associated viscosity is defined as the ratio of elastic moduli and characteristic time, that is, η i = E i /τ i . A complete description of the method is given, for example, in [START_REF] Baumgaertel | Interrelation between continuous and discrete relaxation time spectra[END_REF].

The method discussed next is based on the general remark that a continuous function can be approximated by a staircase function of variables widths (see [START_REF] Kurtz | Theory of Integration[END_REF]). The method will compute the set relaxation times τ i based on a nonlinear optimization method minimizing the integral distance between the continuous relaxation function and its staircase approximation. A detailed presentation of the method is given in Jalocha et al. (2015a).

The discretized spectrum, represented as a staircase function S(τ ) with n elements (see Fig. 13), is expressed as

H (τ i )U (t i-1 ,t i ) (τ ) (24)