S Andrieu

F Schmidt

Y Le Maoult

Infrared heating modelling of thermoplastic sheets in thermoforming process

The thermoforming process involves a heating step before forming, the thermoplastic sheet being warmed up using infrared heaters. The forming temperature is above glass transition temperature for amorphous polymers and slightly below crystalline temperature for semi-crystalline polymers. Today, industrial thermoforming machines generally feature an infrared oven with ceramics heaters that emit infrared long waves. Other electric radiant heaters like quartz tubes or halogen lamps are also used, but not as much. The heating step of the thermoforming process could be improved to obtain best sheet's temperature distribution which governs the thermoformed part thickness distribution. Many parameters influence the heating step such as geometry and number of heaters, heaters temperature, position of the heaters, distance between heaters and sheet… An experimental infrared heating device is developed in order to measure the effect of the different parameters. Software called PLASTIRAD based on the volume method is used to simulate the infrared heating of polymers in order to optimise the infrared oven of an experimental set-up.

1-Introduction

The principle of thermoforming process is sketched in Fig. 1. Thermoforming of thermoplastic sheets requires a heating step before forming. This heating step induces temperature gradients inside the thermoplastic sheet which in turn determines the thickness distribution and consequently the quality of the final part. Due to low polymer thermal conductivity (0.29 W/m/K for PET), this step is performed using infrared heaters which allow rapid heating of the polymers. Many polymers are frequently used in the thermoforming process such as Polystyrene, Polypropylene, ABS, PET, HIPS… [START_REF] Throne | Technology of thermoforming[END_REF]. In the case of PET (crystallisable), the forming temperature is in the range between 130 to 145 °C [START_REF] Illig | Pratique du thermoformage[END_REF]. Generally thermoforming machines are equipped with ceramic heaters rather than halogen lamps which are currently used in injection stretch-blow moulding machines. The industrial problematic is thus to evaluate the capability of halogen lamps in thermoforming process.

Heating with halogen lamps presents different advantages such as fast heat-up and cool down, fast temperature-change responsiveness, zoning capability. These aspects must be considered to figure out the right heating process tuning that permit an optimised temperature distribution and consequently improve the thermoformed parts quality. A lot of different process parameters should be investigated. Therefore, we have developed an experimental infrared heating set-up to study the influence of these parameters. Numerical simulations are performed in order to compute temperature distribution in the thermoplastic sheet versus position of lamps. Computations have been processed using the software called PLASTIRAD based on a control-volume method and taking into account halogen lamps and thermoplastic sheet spectral properties.

2-Characterisation of the infrared heaters

An experimental set-up shown in Fig. 2 has been developed in order to measure with a thermopile the spatial directivity of the infrared lamps used in industrial applications [START_REF] Monteix | Profiles infrared radiative heating in blow molding and thermoforming[END_REF].The thermopile's bandwidth extends between 0.6 and 25 µm. The thermopile moves around the lamp on half circle. It intercepts the radiation in each direction from 0 up to 180° with an increment of 10°. The measured intensity is normalised in relation to the maximum intensity observed in all directions. A comparative directivity study was performed on ceramic heaters (Elstein 650W), quartz tube (TQS 600W) and halogen lamps. These electric radiant heaters used in thermoforming process present a maximum emission in different range of wavelength (long, average and short respectively). In Fig. 3, for each kind of radiant heater, the normalised intensity is plotted versus the angle in comparison with a Lambertian source. The halogen lamp exhibits the closest behaviour to lambertian source compared to ceramic heater and quartz tube. Numerical simulations could be performed using the PLASTIRAD software that computes energy emitted using a view factor calculation.

3-Experimental infrared heating set-up

An experimental infrared heating set-up (see Fig. 4, 5 and Table 1) is developed in order to measure the influence of the process parameters on the heating step of the thermoforming process. This original set-up is designed like a manual thermoplastic machine. It is constituted of an aluminium structure (170*91*130 cm 3 ) that supports a metallic square frame clamping the thermoplastic sheet and a movable infrared oven. Three dimensions of sheet's effective surfaces (500*500, 350*350 and 150*150 mm 2 ) are available to study the effect of the dimension on the sheet's sag for example. In the oven's frame two different modulus could be adapted. One modulus is constituted of a diffusive flat aluminium reflector and eighteen 1000W-Philips halogen lamps (13195Z/98-235V). In the second modulus, individual diffusive aluminium reflectors (IR3-Philips) are adapted on six 1000W (13713Z-235V) and six 500W Philips lamps (13169Z-235V).

This device is designed to control many parameters. The distance between the thermoplastic sheet and the lamps is adjustable. The lamps position and reflectors could be changed. Each lamp is powerregulated. Infrared oven of industrial thermoforming machine (750*750 mm 2 ) could be adapted in this device to study the sheet's temperature distribution that induces. This set-up is equipped with a CDD camera to measure sheet's sag versus heating time and an infrared camera (880 LW AGEMA) to measure surface temperature during heating using a polished aluminium mirror. Thermocouples can be added to measure air temperature between oven and thermoplastic sheet.

4-Control-volume model

The PLASTIRAD software has been developed in order to compute heat transfer between polymer and heaters during an infrared heating step. It is based on a control-volume method [START_REF]Numerical Heat Transfer and Fluid Flow[END_REF]. Previous comparisons performed between numerical simulations and experiments have validated the model [START_REF] Monteix | Experimental study and numerical simulation of preform or sheet exposed to infrared radiative heating[END_REF].

The thermoplastic sheet is meshed (solid mesh) using FEM quadratic elements adapted in control volumes. The halogen lamps used are composed of a coiled tungsten wire enclosed in a quartz tube (diameter: 10 mm) filled with a mixing of halogen and inert gas. Spectral properties of tungsten filament, filling gas and quartz tube have been measured in previous works [START_REF] Monteix | Experimental study and numerical simulation of preform or sheet exposed to infrared radiative heating[END_REF][START_REF] Le Maoult | Measurements and calculation of preform infrared heating[END_REF]. The tungsten filament and the quartz tube are modelled as half-cylinders. Each half-cylinder is meshed (shell mesh) using linear elements.

The temperature balance equation that takes into account radiative and conductive transfers is integrated over each control volume and over time.

∫ ∫ ∫ ∫ ∫ ∫ ∆ ∆ Γ ∆ Γ Γ - Γ - = ∂ ∂ t V t t r c p dt d n q dt d n q dVdt t T C r r r r . . ρ (1)
Where V is the control volume, q c and q r the conductive and radiative fluxes and Γ the surface. In the model, a uniform source temperature is assumed for the tungsten wire and the quartz tube. Temperatures are computed at the cell centres of each control-volume. The amount of incident radiation that reaches each elements of the sheet is computed using a contour integration method based upon Stoke's theorem [START_REF] Rammohan | Efficient evaluation of diffuse view factors for radiation[END_REF]. This method permits good accuracy and is lower time-consuming than other methods, e.g. the Monte-Carlo method.

∫ ∫ Γ Γ = i j j i H H ij r d r rd S F r r ln 2 1 π (2)
Where subscript H refers to the heater i and j to the elementary heater and sheet surfaces. S H is the total heater's surface. Γ i and Γ j are the elementary surfaces coming from the mesh. r refers to the distance between the heater and the sheet surface. The radiative flux received by each volume control surface is computed using Eq. 3:

( ) ∫ ∫ ∆ ∆ + = = Φ λ λ λ λ λ λ λ λ π τ λ π ε τ d T L S S F d T L S S F y Qi Q j Qi Q ij Wi Q j Wi W ij ij ) ( ) ( 0 0 0 (3)
Where F ij W and F ij Q are the view factors referring to tungsten and quartz respectively, S is the surface, τ Qλ the quartz transmittivity, ε λ is the spectral tungsten emissivity, L λ 0 is the blackbody intensity and λ is a given wavelength. The thermoplastic sheet temperature (# 140°C) is lower than the halogen lamp temperature (# 2000°C) so the sheet could be considered as a cold material. Within the sheet thickness, the Beer-Lambert's law gives the transmitted flux:

y k ij ij e y y λ - = Φ = Φ ) 0 ( ) ( (4)
k λ is the spectral absorption coefficient.

5-Numerical simulations 5.1 Process parameters

Simulations are performed to conceive the original configuration of the infrared oven in the experimental heating set-up. The sheet's temperature depending on the position of the lamps in the oven was first analysed. We proceeded to heat a PET sheet (500*1.5*500 mm) with eighteen 1100W(rated power)-Philips halogen lamps (14105Z/98) supplied with 1000W. The distance between lamps and sheet is 50 mm. The spectral properties of PET have been obtained from previous works [START_REF] Monteix | Recent issues in preform radiative heating modelling[END_REF]. Tungsten wire and quartz tube are modelling as half cylinders. Their temperatures, necessary for all numerical simulations are consigned in Table 2.

Table 2 -Operating conditions

Tungsten wire

Quartz Tube Temperature (K) 2539 856

The sheet is meshed (solid-mesh) using 900 elements (15*4*15, length*thickness*width). Quartz tube and tungsten wire of each halogen lamp are meshed with 25 elements (5*5, length*circumference). Three different intervals between lamps 60, 70 and 80 mm are compared with the zigzag configuration described in Fig. 6, 7 and 8 respectively. For the three different configurations, the thermoplastic sheet requires a different heating time (see Table 3) in order to reach the forming temperature. After this heating time, the thermoplastic sheet undergoes the air convection cooling. In Fig. 9, the sheet centre temperature (x = 250 mm and z = 250 mm) is plotted versus time for the surface exposed to radiation (front surface) and the back surface. In the section plane (z = 250 mm), temperature is plotted versus x for the front surface and the back

Results

Fig. 9 demonstrates the temperature evolution versus increased intervals between lamps axis. The heating time to reach the forming temperature increases with the distance between lamps axis. On the other hand, the distribution temperature versus x is better with the 80 mm interval. Actually, with 60 mm interval, the sheet sides are not correctly heated (weakest temperature). A comparison of the different temperature profiles shows that the x direction is privileged. It is due to the lamp geometry, as the heating is better in the x direction perpendicular to the lamps axis. The jump around z=300 mm in Fig. 12 and 13 is due to the presence of a lamp right above the middle line x=250mm from z=250 mm up to z=500 mm. A 5°C gap between front and back surface temperatures is observed. According to all these representations, in our numerical simulations, 80 mm interval simulation is the configuration which gives the best temperature distribution. The main idea is to use energy criteria to compare different simulations independently of the dimensions. Heat balance equation (Eq. 1) could be rewritten (see Eq. 5) in order to deduce the source term S i corresponding to the energy received by each control volume.

6-Energy criteria 6.1 Distribution energy and overall energy criteria

(

) 5
Where S i is the source term. The PLASTIRAD software computes the source term for each control volume. The distribution energy criterion R i presents the discrepancy between the source term of each element and the maximum source term found over the entire thermoplastic sheet.

R i varies from 0 to 1. Therefore, an arrangement that gives R i close to 1 is the eventual goal. This distribution energy criterion is presented for 60 mm, 70mm and 80 mm intervals in Fig. 14,15 and 16 respectively and for the front surface. R i discrepancy for the 60 mm interval is of 50%, for the other configurations is about 30%. The configuration that gives the best uniformity is 80 mm interval except a energy concentration in the sheet centre. Moreover, an overall energy criterion is introduced as a complementary tool in order to compare the infrared radiation efficiency for different configurations. The criterion used E is defined by the ratio energy received by energy emitted:

∑∑ ∫ ∑∑

∫ ∑ ∫         Γ - +         Γ ∇ = ∂ ∂ Σ Σ i S j r conduction i j i V p i ij ij i d n q d n T k dV t T C
In our simulations, we use eighteen 1000W halogen lamp without reflectors so ΣP L =18*1000/2=9000W. Considered that, the calculations of the overall energy criterion are compiled in Table 4. To optimize the infrared heating, a power regulation is planned for introduction in the numerical simulations as well as in the experimental infrared heating set-up. An analysis of Fig. 16 that is the representation of the best configuration (80 mm interval), induces the idea of decreasing the power of the lamps which are located in the zone surrounding the R i maximum position: the lamps that are thus power-regulated are indicated in Fig. 17. The new R i distribution (see Fig. 18) is more homogenous and the overall energy criterion for this particular setting is 13%. The choice of the final configuration depends on what is expected, in our case, it is in a first stage to obtain a temperature uniformity over the all sheet. Consequently, the chosen configuration is the 80 mm interval.

Conclusion

Numerical simulations of PET sheet heating with halogen lamps have been performed to optimize an infrared oven. Due to lamps geometry, a direction is privileged by the choice of the orientation of the lamps. In our case, this direction is perpendicular to the lamp axis. Except the sides of the sheet, a satisfactory temperature distribution is obtained. Energy criteria have been used to compare the different simulations. The distribution energy criterion is useful to avoid a spatial energy distribution. The overall energy criterion is a tool to evaluate efficiency of the heating step. So temperature computations and energy criteria allow to configurate lamps position in the infrared oven. Further numerical model developments will permit to integrate the real direction beam into the sheet and the efficiency of reflectors. Moreover, the experimental infrared heating device will permit to compare experimental measurements and numerical simulations and finally to validate the model developments

Figure 1 -

 1 Figure 1 -The different steps of thermoforming

Figure 2 -

 2 Figure 2 -Experimental set-up used in order to measure spatial directivity

Figure 3 -Figure 4 - 5 -

 345 Figure 3 -Infrared heater's directivity

Figure 6 -Figure 7 -Figure 8 -

 678 Figure 6 -60 mm interval

Figure 10 -Figure 11 -Figure 12 -Figure 13 -

 10111213 Figure 10 -Front surface temperature versus x

Figure 15 -

 15 Figure 15 -70 mm interval distribution energy criterion

  Figure 18 -80 mm interval distribution energy criterion with power regulation

Table 1 -

 1 Designation

	Number	Designation
	1	Aluminum structure
	2	Metallic frame
	3	Thermoplastic sheet
	4	Infrared oven
	5	Electric regulation
	6	CCD camera
	7	Infrared camera
	8	Aluminium mirror
	9	Fixing screws

Table 3 -

 3 Heating times

		60 mm interval	70 mm interval	80 mm interval
	Heating time (s)	50	60	74

Table 4 60

 4

		mm interval	70 mm interval	80 mm interval
	E (%)	17	15	14

Acknowledgements

The authors would like to thank our industrial partners Philips Lighting and EDF-Renardières Companies for technical support and funding.


Harold Berjamin

Bruno Lombard
email: lombard@lma.cnrs-mrs.fr

Guillaume Chiavassa

Nicolas Favrie

Plane-strain waves in nonlinear elastic solids with softening

Keywords: Nonlinear elastodynamics, softening, solitary waves, numerical methods PACS: 43.25.Dc, 02.70.Bf

Propagation of elastic waves in damaged media (concrete, rocks) is studied theoretically and numerically. Such materials exhibit a nonlinear behavior, with long-time softening and recovery processes (slow dynamics). A constitutive model combining Murnaghan hyperelasticity with the slow dynamics is considered, where the softening is represented by the evolution of a scalar variable. The equations of motion in the Lagrangian framework are detailed. These equations are rewritten as a nonlinear hyperbolic system of balance laws, which is solved numerically using a finite-volume method with flux limiters. Numerical examples illustrate specific features of nonlinear elastic waves, as well as the effect of the material's softening. In particular, the generation of solitary waves in a periodic layered medium is illustrated numerically.

Introduction

Geomaterials such as rocks and concrete exhibit nonlinear features at small strains ε ∼ 10 -6 . In particular, longitudinal resonance experiments reveal the generation of higher-order harmonics. Besides this phenomenon, the material exhibits also a softening with increasing strain amplitudes [START_REF] Johnson | Resonance and elastic nonlinear phenomena in rock[END_REF][START_REF] Remillieux | Decoupling nonclassical nonlinear behavior of elastic wave types[END_REF]. After the forcing is stopped, the material recovers gradually its initial stiffness. This softening/recovery process is not instantaneous. The transient regime, named "slow dynamics" in the corresponding literature [START_REF] Ten Cate | Slow dynamics in the nonlinear elastic response of berea sandstone[END_REF][START_REF] Tencate | Slow dynamics of earth materials: An experimental overview[END_REF][START_REF] Rivière | Pump and probe waves in dynamic acousto-elasticity: Comprehensive description and comparison with nonlinear elastic theories[END_REF][START_REF] Haupert | Dynamic acousto-elastic testing[END_REF], has a characteristic time much larger than the period of the forcing.

To describe nonlinear elastic behavior, the finite-strain theory is a self-consistent framework. Various constitutive laws of isotropic hyperelastic material express the stress as a function of a strain tensor, and are compatible with Hooke's law in the infinitesimal strain limit [START_REF] Ogden | Non-linear Elastic Deformations[END_REF][START_REF] Holzapfel | Nonlinear Solid Mechanics: A Continuum Approach for Engineering[END_REF][START_REF] Norris | Finite-amplitude waves in solids[END_REF]. Among them, the Murnaghan law [START_REF] Murnaghan | Finite deformations of an elastic solid[END_REF] is frequently used to describe acoustic nonlinearity in rocks and concrete, and values of the parameters can be found in the literature [START_REF] Johnson | Nonlinear elasticity and stress-induced anisotropy in rock[END_REF][START_REF] Winkler | Measurements of third-order elastic constants in rocks[END_REF][START_REF] Payan | Determination of third order elastic constants in a complex solid applying coda wave interferometry[END_REF]. However, hyperelasticity does not account for the long-time relaxation of elastic constants.

Several models from the literature describe the softening of the material by a dependence of the elastic constants on a scalar variable g [START_REF] Lyakhovsky | Distributed damage, faulting, and friction[END_REF][START_REF] Vakhnenko | Strain-induced kinetics of intergrain defects as the mechanism of slow dynamics in the nonlinear resonant response of humid sandstone bars[END_REF][START_REF] Berjamin | Nonlinear waves in solids with slow dynamics: an internalvariable model[END_REF]. In the present study, the internal-variable model [START_REF] Berjamin | Nonlinear waves in solids with slow dynamics: an internalvariable model[END_REF] is used, where the stress issued from Murnaghan's law is multiplied by 1g. Hence, g ≡ 0 corresponds to the hyperelastic case, and the material softens as g increases. The evolution of g is governed by a first-order differential equation compatible with the principles of thermodynamics. Both softening processes ġ 0 and recovery processes ġ 0 are possible [START_REF] Berjamin | Nonlinear waves in solids with slow dynamics: an internalvariable model[END_REF], accordingly to the experimental observations.

In recent works [START_REF] Berjamin | A finite-volume approach to 1D nonlinear elastic waves: Application to slow dynamics[END_REF], the corresponding one-dimensional equations of motion have been solved numerically using a finite-volume method with flux limiters [START_REF] Leveque | Finite Volume Methods for Hyperbolic Problems[END_REF]. The numerical method has been validated with respect to reference solutions, and qualitative agreement with experimental observations has been obtained. The present article describes a similar method in the two-dimensional plane-strain case. The numerical method is well-suited to the computation of nonlinear waves in the Lagrangian framework, and it can be used for various hyperelastic material models (cf. the related study [START_REF] Bou Matar | A nodal discontinuous Galerkin finite element method for nonlinear elastic wave propagation[END_REF] and references therein).

The article is organized as follows. In Section 2, the equations of motion are detailed, which are rewritten as a nonlinear hyperbolic system of balance laws in two space variables. The numerical method is presented in Section 3. Section 4 shows 2D numerical results illustrating the nonlinear wave propagation in elastic solids with softening. The case of a homogeneous medium is considered, as well as the case of a periodic layered medium with solitary wave solutions.

Governing equations

Lagrangian hyperelasticity with softening

We consider an homogeneous continuum. A particle initially located at some position x 0 of the reference configuration moves to a position x t of the deformed configuration. The deformation gradient is a second-order tensor defined by (see e.g. [START_REF] Ogden | Non-linear Elastic Deformations[END_REF][START_REF] Holzapfel | Nonlinear Solid Mechanics: A Continuum Approach for Engineering[END_REF][START_REF] Norris | Finite-amplitude waves in solids[END_REF][START_REF] Drumheller | Introduction to Wave Propagation in Nonlinear Fluids and Solids[END_REF])

F = grad x t = I + grad u , (1)
where u = x tx 0 denotes the displacement field and grad is the gradient with respect to the material coordinates x 0 (Lagrangian gradient). In the reference configuration, the deformation gradient (1) is equal to the metric tensor I. Here, the Euclidean space is described by an orthonormal basis (e 1 , e 2 , e 3 ) and a Cartesian coordinate system (O, x, y, z), so that the matrix of the coordinates of I is the identity matrix. The Lagrangian representation of motion is used. Hence, the material time derivative Ḟ = ∂ t F of the deformation gradient satisfies

Ḟ = grad v , (2)
where v = u is the velocity field. The conservation of mass implies ρ 0 /ρ = det F, where ρ denotes the mass density in the deformed configuration, and ρ 0 denotes the mass density in the reference configuration. Self-gravitation and heat conduction are neglected, so that the motion is driven by the conservation of momentum

ρ 0 v = div P + f v , where P = det(F) σ • F -⊤ (3)
is the first Piola-Kirchhoff stress tensor, and div denotes the divergence with respect to the material coordinates x 0 . The Cauchy stress tensor σ = (det F) -1 P • F ⊤ = σ ⊤ is detailed later on through a specification of P. The term f v is an external volume force applied to the material.

In hyperelasticity, the only variables of state are the specific entropy and a strain tensor. Here, the Green-Lagrange strain tensor

E = 1 2 (F ⊤ • F -I) is used, i.e. E = 1 2 grad u + grad ⊤ u + grad ⊤ u • grad u . (4)
An internal variable g ∈ [0, 1[ accounting for the softening of the material is added to the previous list of variables of state. We define the internal energy by unit of reference volume as [START_REF] Berjamin | Nonlinear waves in solids with slow dynamics: an internalvariable model[END_REF] ρ

0 e = (1 -g) W(E) + Φ(g) , ( 5
)
where e is the specific internal energy. In [START_REF] Rivière | Pump and probe waves in dynamic acousto-elasticity: Comprehensive description and comparison with nonlinear elastic theories[END_REF], the potential W(E) is the strain energy function of Murnaghan's law [START_REF] Murnaghan | Finite deformations of an elastic solid[END_REF] W

(E) = λ + 2µ 2 E I 2 -2µE II + l + 2m 3 E I 3 -2mE I E II + nE III , (6)
which is expressed as a function of the strain invariants

E I = tr E, E II = 1 2 (tr E) 2 -tr(E 2 )
, and E III = det E. The constants λ, µ are the Lamé parameters and the constants l, m, n are the Murnaghan coefficients (third-order elastic constants). In the case where the Murnaghan coefficients in (6) equal zero, the strain energy of the Saint Venant-Kirchhoff model is recovered. The potential Φ(g) represents a storage energy. Basic requirements are the convexity of Φ and Φ ′ (0) = 0, where Φ ′ is the derivative of Φ (see [START_REF] Berjamin | Nonlinear waves in solids with slow dynamics: an internalvariable model[END_REF]). Moreover, an asymptote at g = 1 is introduced to avoid the destruction of the material. A suitable expression is

Φ(g) = - 1 2 γ ln(1 -g 2 ) , (7)
where γ > 0 is an energy per unit volume. The energy ( 7) is quadratic Φ(g) ≃ 1 2 γg 2 in the limit g → 0.

Under the previous assumptions, the simplest set of constitutive equations with softening which is thermodynamically admissible reads [START_REF] Berjamin | Nonlinear waves in solids with slow dynamics: an internalvariable model[END_REF]

P = (1 -g) F • ∂W ∂E , and
τ 1 ġ = W(E) -Φ ′ (g) , (8)
where τ 1 > 0 in J m -3 s is a material parameter. The first equation in [START_REF] Holzapfel | Nonlinear Solid Mechanics: A Continuum Approach for Engineering[END_REF] is the mechanical constitutive law, which reduces to the case of hyperelasticity if g ≡ 0. According to the second equation in [START_REF] Holzapfel | Nonlinear Solid Mechanics: A Continuum Approach for Engineering[END_REF] which governs the evolution of the variable g, the classical theory of nonlinear elastodynamics is recovered if τ 1 → +∞. Similarly to [START_REF] Berjamin | Nonlinear waves in solids with slow dynamics: an internalvariable model[END_REF][START_REF] Berjamin | A finite-volume approach to 1D nonlinear elastic waves: Application to slow dynamics[END_REF], we rewrite the tensor derivative ∂W/∂E of the strain energy using the invariants' tensor derivatives:

∂W ∂E = α 0 I + α 1 E + α 2 E 2 , (9)
where

α 0 = ∂W ∂E I + E I ∂W ∂E II + E II ∂W ∂E III = λE I + lE I 2 -(2m -n)E II , α 1 = - ∂W ∂E II -E I ∂W ∂E III = 2µ + (2m -n)E I , α 2 = ∂W ∂E III = n . (10)
In the next subsection, we detail the case of plane strain. The case of uniaxial strain is addressed in [START_REF] Berjamin | A finite-volume approach to 1D nonlinear elastic waves: Application to slow dynamics[END_REF].

The plane-strain assumption

The displacement field u is independent of z, and its component u 3 along e 3 is zero. In the basis of unit tensors (e i ⊗e j ) 1 i, j 3 , the matrix of coordinates of the displacement gradient grad u is therefore

u i, j =             u 1,1 u 1,2 0 u 2,1 u 2,2 0 0 0 0             . ( 11
)
Using the Einstein notation with indices in {1, 2}, the coordinates of the Green-Lagrange tensor (4) write E i j = 1 2 u i, j + u j,i + u p,i u p, j . Its invariants in ( 6) are

E I = E nn , E II = 1 2 E I 2 -E i j E i j = ǫ i j E 1i E 2 j
and E III = 0, where ǫ i j is the Levi-Civita symbol of R 2 . The Cayley-Hamilton theorem applied to the restriction of E to R 2 × R 2 reads E im E m j -E I E i j + E II δ i j = 0, where δ i j is the Kronecker delta. Hence, the expression (9) of ∂W/∂E becomes

∂W ∂E i j = α0 δ i j + α1 E i j ( 12
)
in the basis of unit tensors (e i ⊗ e j ) 1 i, j 2 , where

α0 = α 0 -α 2 E II = λE I + lE I 2 -2mE II , α1 = α 1 + α 2 E I = 2 (µ + mE I ) . (13)
The components P i j of the first Piola-Kirchhoff stress tensor P in (8) are therefore

P i j = (1 -g) δ im + u i,m α0 δ m j + α1 E m j ( 14
)
under the plane strain assumption, which does not depend upon the third Murnaghan coefficient n.

When the geometric nonlinearities are negligible, the Green-Lagrange strain tensor (4) is linearized with respect to grad u, i.e. E ≃ 1 2 grad u + grad ⊤ u = ε reduces to the infinitesimal strain tensor. The coordinates E i j of E are replaced by the coordinates ε i j = 1 2 u i, j + u j,i of ε. Moreover, the first Piola-Kirchhoff stress tensor P is linearized with respect to grad u too, i.e. F • ∂W/∂E ≃ ∂W/∂ε. Hence, the equation P i j = (1g) α0 δ i j + α1 ε i j replaces [START_REF] Lyakhovsky | Distributed damage, faulting, and friction[END_REF]. Under this assumption, linear elastodynamics is recovered if g ≡ 0 (i.e., τ 1 → +∞ in (8)), and if the Murnaghan coefficients l, m in (13) are zero (see Appendix A for details).

System of balance laws

Under the plane-strain assumption, the equations of motion ( 2)-( 3)-( 8) are rewritten as a two-dimensional system of balance laws

∂ t q + ∂ x f(q) + ∂ y g(q) = r(q) + s , (15)
where only three strips are displayed (everywhere else, the coefficients in the matrix are zero).

The expression of the coefficients Q i jkℓ in (17) defined by ρ 0 Q i jkℓ = ∂P i j /∂u k,ℓ is detailed in the Appendix A, as well as the expression of the coefficients G i j defined by ρ 0 G i j = ∂P i j /∂g. A similar Jacobian matrix g ′ (q) is obtained for the flux component g along the y-axis. These matrices are diagonalized in the Appendix A. The spectrum of both matrices has the form {-c P , c P , -c S , c S , 0, 0, 0}. In the case of Murnaghan hyperelasticity, the eigenvalues c P , c S can be complex [START_REF] Berjamin | A finite-volume approach to 1D nonlinear elastic waves: Application to slow dynamics[END_REF], so that the system (15)-( 16) is not unconditionally hyperbolic (see e.g. [START_REF] Ndanou | Criterion of hyperbolicity in hyperelasticity in the case of the stored energy in separable form[END_REF] for discussions on hyperbolicity in hyperelasticity). Here, we restrict ourselves to configurations where the eigenvalues c P > c S > 0 are real. Thus, c P and c S correspond to the velocities of compression waves and shear waves, respectively.

Plane waves. We assume furthermore that the displacement field is invariant along a direction, say e 2 , so that u does not depend on y. In this case, the vector of unknown reduces to q =

(u 1,1 , u 2,1 , v 1 , v 2 , g) ⊤ -
the second and fourth rows of ( 15)-( 16) are zero -and the flux g along y is zero. The Jacobian matrix f ′ (q) is obtained from ( 17) by removing the second and fourth rows, as well as the second and fourth columns. Doing so, two zero eigenvalues are removed from the spectrum, which reduces to {-c P , c P , -c S , c S , 0}. We consider the case of Murnaghan material g ≡ 0, with the parameters in Table 1. The latter, found in [START_REF] Payan | Determination of third order elastic constants in a complex solid applying coda wave interferometry[END_REF], have been measured on concrete. Let us introduce the relative variation ∆c/c of the sound velocities c P and c S with respect to the case of Hooke's law, where c P = (λ + 2µ)/ρ 0 ≈ 4458 m/s and c S = µ/ρ 0 ≈ 2700 m/s. Fig. 1 displays the evolution of ∆c/c with respect to the compression strain u 1,1 , when the shear strain u 2,1 is set to zero. One observes that the variations The variation ∆c/c is represented with respect to the compression strain u 1,1 , and the shear strain is u 2,1 = 0. The dotted lines mark Taylor series approximations [START_REF] Leveque | Finite Volume Methods for Hyperbolic Problems[END_REF].

ρ 0 (kg m -3 ) λ (GPa) µ (GPa) l (GPa) m (GPa) γ (J m -3 ) τ 1 (J m
of c P with respect to u 1,1 are much larger than the variations of c S . This is confirmed by the Taylor series approximations

(∆c/c) P = 3 2 + l + 2m λ + 2µ u 1,1 + O(u 1,1 2 ) + O(u 2,1 2 ) ≈ -157 u 1,1 , (∆c/c) S = λ + 2µ 2µ + m 2µ u 1,1 + O(u 1,1 2 ) + O(u 2,1 2 ) ≈ -63.8 u 1,1 , (18)
represented as dotted lines in Fig. 1 (the magnitude of u 1,1 in the figure has been chosen for graphical reasons). These approximations show also that the shear strain u 2,1 has much less influence than the compression strain u 1,1 on the variations of the sound velocities.

Numerical resolution

Numerical strategy

In the examples presented later on, the physical domain is assumed unbounded. We consider a finite numerical domain [0, L x ] × [0, L y ]. It is discretized using a regular grid in space with mesh size ∆x in the x direction, and ∆y in the y direction. The coordinates of the nodes are (x i , y j ) = (i ∆x, j ∆y), where 0 i N x and 0 j N y . The total number of nodes is (N x + 1) × (N y + 1), where N x = L x /∆x and N y = L y /∆y. A variable time step ∆t = t n+1t n is introduced. Therefore, q(x i , y j , t n ) denotes the solution to (15) at the grid node (i, j) and at the nth time step. Numerical approximations of the solution are denoted by q n i, j ≃ q(x i , y j , t n ). The non-homogeneous system of balance laws [START_REF] Vakhnenko | Strain-induced kinetics of intergrain defects as the mechanism of slow dynamics in the nonlinear resonant response of humid sandstone bars[END_REF] is integrated explicitly in time:

q n+1 i, j = q n i, j + ∆q n FV + ∆t r(q n i, j ) + s n i, j , (19)
where the approximation s n i, j of the source term s is specified later on. The increment ∆q n FV is deduced from the integration of ∂ t q + ∂ x f(q) + ∂ y g(q) = 0 over one time step. Usually, one has

∆q n FV = (H x + H y -2) q n i, j
, where the discrete operators

H x q n i, j = q n i, j - ∆t ∆x f n i+1/2, j -f n i-1/2, j , H y q n i, j = q n i, j - ∆t ∆y g n i, j+1/2 -g n i, j-1/2 , (20)
involve the fluxes f n i+1/2, j , g n i+1/2, j of a 2D finite-volume scheme [START_REF] Leveque | Finite Volume Methods for Hyperbolic Problems[END_REF]. Here, a second-order symmetric dimensional splitting [START_REF] Strang | Accurate partial difference methods I: Linear Cauchy problems[END_REF] is used instead. That is to say, H x and H y correspond to the integration of ∂ t q + ∂ x f(q) = 0 and ∂ t q + ∂ y g(q) = 0 over one time step, so that (20) involves the fluxes f n i+1/2, j , g n i+1/2, j of a 1D finite-volume scheme. The increment ∆q n FV is computed according to

∆q n FV = 1 2 H x H y + H y H x -2 q n i, j , (21)
where H x H y denotes the composition of the operators H x and H y .

The numerical fluxes in ( 20)-( 21) are computed according to the flux-limiter method [START_REF] Berjamin | A finite-volume approach to 1D nonlinear elastic waves: Application to slow dynamics[END_REF][START_REF] Leveque | Finite Volume Methods for Hyperbolic Problems[END_REF] described in the next subsection. This finite-volume scheme is well-suited for nonlinear wave propagation and second-order accurate. The operators H x and H y are stable under the Courant-Friedrichs-Lewy (CFL) condition

Co = max 0 i N x 0 j N y max ̺ f ′ (q n i, j ) ∆t ∆x , ̺ g ′ (q n i, j ) ∆t ∆y 1 , (22)
where Co is the maximum Courant number in the x and y directions. The spectral radius ̺ f ′ (q) of f ′ (q) corresponds to c P (expression detailed in the Appendix A), ditto the spectral radius ̺ g ′ (q) of g ′ (q). The stability of the scheme ( 19) is also restricted by the spectral radius of the Jacobian matrix r ′ (q). As in 1D [START_REF] Berjamin | A finite-volume approach to 1D nonlinear elastic waves: Application to slow dynamics[END_REF], the stability limits imply that the scheme ( 19) is stable under the classical CFL condition [START_REF] Strang | Accurate partial difference methods I: Linear Cauchy problems[END_REF]. Hence, given a spatial discretization and a Courant number Co 1, the value of the time step ∆t is imposed by [START_REF] Strang | Accurate partial difference methods I: Linear Cauchy problems[END_REF].

Flux limiter

We describe now the flux-limiter scheme [START_REF] Berjamin | A finite-volume approach to 1D nonlinear elastic waves: Application to slow dynamics[END_REF][START_REF] Leveque | Finite Volume Methods for Hyperbolic Problems[END_REF]. Since the computation of the numerical fluxes in the x and y directions is similar, only the numerical flux f n i+1/2, j in the x direction is detailed here. To do so, we introduce the Jacobian matrix

A i+1/2, j = f ′ 1 2 (q n i, j + q n i+1, j ) (23)
at the arithmetic mean of the grid node values in the x direction. The jump of the numerical solution q n i+1, j -q n i, j along x is decomposed in the basis of right eigenvectors

{p k i+1/2, j , k = 1, . . . , 7} of A i+1/2, j , q n i+1, j -q n i, j = 7 k=1 α k i+1/2, j p k i+1/2, j = 7 k=1 W k i+1/2, j , (24)
which correspond to the eigenvalues {-c P , c P , -c S , c S , 0, 0, 0} (cf. detailed expressions in the Appendix A).

The numerical flux in [START_REF] Drumheller | Introduction to Wave Propagation in Nonlinear Fluids and Solids[END_REF] is the sum of a first-order flux and a second-order limited correction,

f n i+1/2, j = f L i+1/2, j + f H i+1/2, j
, where

f L i+1/2, j = 1 2 f(q n i, j ) + f(q n i+1, j ) - 1 2 c P W 1 i+1/2, j + W 2 i+1/2, j - 1 2 c S W 3 i+1/2, j + W 4 i+1/2, j , f H i+1/2, j = 1 2 c P 1 - ∆t ∆x c P φ(θ 1 i+1/2, j )W 1 i+1/2, j + φ(θ 2 i+1/2, j )W 2 i+1/2, j + 1 2 c S 1 - ∆t ∆x c S φ(θ 3 i+1/2, j )W 3 i+1/2, j + φ(θ 4 i+1/2, j )W 4 i+1/2, j . (25)
The coefficients θ k i+1/2, j where k = 1, . . . , 4 express the upwind variation of the jump [START_REF] Remillieux | Propagation of a finite-amplitude elastic pulse in a bar of berea sandstone: A detailed look at the mechanisms of classical nonlinearity, hysteresis, and nonequilibrium dynamics[END_REF] in the kth characteristic field,

θ 1,3 i+1/2, j = W 1,3 i+3/2, j • W 1,3 i+1/2, j W 1,3 i+1/2, j • W 1,3 i+1/2, j , θ 2,4 i+1/2, j = W 2,4 i-1/2, j • W 2,4 i+1/2, j W 2,4 i+1/2, j • W 2,4 i+1/2, j , (26)
and φ denotes the minmod limiter function φ(θ) = max{0, min{1, θ}}. As such, the weights φ(θ k i+1/2, j ) are designed to avoid spurious oscillations in the numerical solution. Since the eigenvalues indexed by k = 5, . . . , 7 in the decomposition of the jump [START_REF] Remillieux | Propagation of a finite-amplitude elastic pulse in a bar of berea sandstone: A detailed look at the mechanisms of classical nonlinearity, hysteresis, and nonequilibrium dynamics[END_REF] are zero, the corresponding terms W k i+1/2, j do not appear in the numerical flux [START_REF] Fogarty | High-resolution finite-volume methods for acoustic waves in periodic and random media[END_REF]. To carry out one iteration in time ( 19)-( 21) at some grid node (i, j), the numerical values of q at the grid nodes (i-2, . . . , i+2) ×( j-2, . . . , j+2) are required [START_REF] Fogarty | High-resolution finite-volume methods for acoustic waves in periodic and random media[END_REF]. Therefore, two columns and two rows of "ghost cells" are added on the left, the right, the top, and the bottom of the numerical domain. If not specified differently, a zero-order extrapolation of the numerical values is used to update the ghost cell values at each step of [START_REF] Ndanou | Criterion of hyperbolicity in hyperelasticity in the case of the stored energy in separable form[END_REF]. This procedure is detailed in Section 21.8 of [START_REF] Leveque | Finite Volume Methods for Hyperbolic Problems[END_REF], and is used here to simulate outflow boundary conditions (i.e., an infinite physical domain).

Numerical experiments

In the following numerical examples, the Courant number ( 22) is set to Co = 0.9. If not specified otherwise, the physical parameters are given in Table 1. The parameters γ, τ 1 have been chosen so as to obtain significant effects of the softening at the scale of the simulation. The numerical domain is defined by L x = L y = 0.4 m, and is discretized using N x = N y = 800 points in each direction.

Murnaghan hyperelasticity

The first example focuses on nonlinear elastodynamics, i.e., no softening occurs in the material. In ( 15)-( 16), the source term and the relaxation function are removed (s = 0, τ 1 → +∞). We consider a Riemann problem with initial data q(x, y, 0), where the material is initially undeformed and opposite transverse velocities with amplitude V are applied:

q(x, y, 0) =        V 0, 0, 0, 0, sin ϕ, -cos ϕ, 0 ⊤ if x ϕ < 0 , V 0, 0, 0, 0, -sin ϕ, cos ϕ, 0 ⊤ if x ϕ > 0 . ( 27
)
The variable x ϕ = (xx s ) cos ϕ + (yy s ) sin ϕ is the x-abscissa of a new coordinate system, corresponding to a rotation by an angle ϕ and a translation by (x s , y s ) of the original one. Here, 8 the origin is set at (x s , y s ) = (L x , L y )/2, the rotation angle is ϕ = 15 • , and the velocity amplitude is V = 0.1 m/s. To reduce discretization artifacts due to the oblique discontinuity, the average value of ( 27) over the cell [x i-1/2 , x i+1/2 ] × [y j-1/2 , y j+1/2 ] is initially set at the grid node (i, j).

Figures 2 and3 illustrate the coupling between plane shear waves and plane compression waves in hyperelasticity [START_REF] Norris | Finite-amplitude waves in solids[END_REF][START_REF] Favrie | A thermodynamically compatible splitting procedure in hyperelasticity[END_REF], contrary to linear elasticity where both types of waves are decoupled. Fig. 2 displays a map of W 1/8 at t = 0.015 ms, where W is the strain energy (6) obtained numerically with the above method. Fig. 3 display the evolution of the rotated longitudinal velocity (v 1 ) ϕ = v 1 cos ϕ + v 2 sin ϕ and the rotated transverse velocity (v 2 ) ϕ = -v 1 sin ϕ + v 2 cos ϕ, along the solid line displayed in Fig. 2. In the case of Hooke's law of linear elasticity, the solution to the initial-value problem [START_REF] Leveque | Solitary waves in layered nonlinear media[END_REF] writes

q(x, y, t) =                    V 0, 0, 0, 0, sin ϕ, -cos ϕ, 0 ⊤ if x ϕ < -c S t , V c S -sin ϕ cos ϕ, -sin 2 ϕ, cos 2 ϕ, sin ϕ cos ϕ, 0, 0, 0 ⊤ if -c S t < x ϕ < c S t , V 0, 0, 0, 0, -sin ϕ, cos ϕ, 0 ⊤ if c S t < x ϕ , (28)
with c S = µ/ρ 0 , and only shear waves propagate (solid line in Fig. 3). In the hyperelastic case, Fig. 3b shows that shear waves are generated from the initial data ( 27), but faster compression waves are also generated. This amplitude-dependent nonlinear effect is better observed at large amplitudes. The small-amplitude perturbations in the numerical solution are caused by the discretization of the oblique discontinuity [START_REF] Leveque | Solitary waves in layered nonlinear media[END_REF].

(a)

0.1 0.15 0.2 0.25

-1 0 1 x (m) (v 2 ) ϕ /V Murnaghan num. Hooke theo. t = 0 (b) 0.1 0.15 0.2 0.25 0.3 -1 0 1 ×10 -2 x (m) (v 1 ) ϕ /V V = 0.1 m/s V = 0.2 m/s V = 1.0 m/s

Softening

In this second example, we consider the full system ( 15)-( 16). The material is initially undeformed and at rest, q(x, y, 0) = 0, and a volume force f v is used for the forcing s. The volume force is an acoustic point source along x with expression

f v = A v sin(2π f c t) δ(x -x s ) δ(y -y s ) e 1 ,
where δ is the Dirac delta function, A v is the amplitude, and f c is the characteristic frequency. Usually, the increment s n i, j in ( 19) is obtained by averaging the source term s| t=t n of ( 15)-( 16) over the cell [x i-1/2 , x i+1/2 ] × [y j-1/2 , y j+1/2 ]. Here, we approximate the Dirac deltas by a truncated Gaussian function to avoid strain concentration at the source. Thus,

s n i, j = A v ρ 0 sin(2π f c t n ) exp -(d i j /σ c ) 2 πσ c 2 1 -exp -(R/σ c ) 2 1 d i j R (0, 0, 0, 0, 1, 0, 0) ⊤ , ( 29
)
where d i j is the distance between (x i , y j ) and (x s , y s ). Denoted by the indicator function 1 d i j R , the function's support is a disk with radius R = c P /(7.5 f c ), where c P = (λ + 2µ)/ρ 0 is the speed of linear compression waves. The width parameter of the Gaussian function is chosen such that σ c = R/2. The point load has amplitude A v = 0.5 kN/m, frequency f c = 100 kHz, and it is located at the nearest grid node of the domain's center: (x s , y s ) ≃ (L x , L y )/2. The source [START_REF] Samsonov | Strain Solitons in Solids and How to Construct Them[END_REF] is switched on at t = 0, and switched off at t = 0.04 ms. Two receivers R1-R2 record the numerical solution during the simulation. R1 is located at (x r , y r ) = (0.2, 0.22) m, and R2 is located at (x r , y r ) = (0.2, 0.27) m. Figure 4 illustrates the effect of the softening on the wave propagation. Fig. 4a displays a map of the strain energy (1g) W at the time t = 0.04 ms, which shows the propagation of cylindrical waves. Denoted by a bullet point and by a square, the receivers R1-R2 are located in a region of the plane where mainly shear waves propagate. Figs. 4b-4c show the effect of the softening at the position of the receivers for several forcing amplitudes A v . One observes that g increases while the wave passes by the receiver, and that it relaxes towards zero afterwards. This softening/recovery process is all the more important as the forcing amplitude is large. The characteristic time of the slow dynamics τ 1 /γ ≈ 0.05 ms corresponds to the characteristic time of the (b)

0 2 4 6 8
×10 -2 0 -3 3 ×10 -3 v 1 (m/s) • R1 0 2 4 6 8
×10 -2 0 2 4 ×10 -2 t (ms) g (c) 0 2 4 6 8
×10 -2 0 -3 3 ×10 -3 v 1 (m/s) R2 0 2 4 6 8
×10 -2 0 2 4 ×10 -2 t (ms) g A v = 0.5 kN/m A v = 0.7 kN/m A v = 1.0 kN/m
Π = 3 2 Π Π = 1 2 Π. . . d L x v 1 (0, y, t) e 1 e 2
Figure 5: Map of the periodic medium's linear material properties Π ∈ {ρ 0 , λ, µ} over the numerical domain, where Π ∈ { ρ0 , λ, μ} denotes the reference values from Table 1.

recovery [START_REF] Berjamin | Nonlinear waves in solids with slow dynamics: an internalvariable model[END_REF][START_REF] Berjamin | A finite-volume approach to 1D nonlinear elastic waves: Application to slow dynamics[END_REF]. In Figs. 4b-4c, one observes the distortion of the velocity signal during propagation, and its delay due to the increase of g. The recorded signals are similar to experimental ones obtained in a longitudinal configuration [START_REF] Remillieux | Propagation of a finite-amplitude elastic pulse in a bar of berea sandstone: A detailed look at the mechanisms of classical nonlinearity, hysteresis, and nonequilibrium dynamics[END_REF].

Periodic layered medium

Similarly to [START_REF] Fogarty | High-resolution finite-volume methods for acoustic waves in periodic and random media[END_REF][START_REF] Leveque | Finite-volume methods for non-linear elasticity in heterogeneous media[END_REF][START_REF] Leveque | Solitary waves in layered nonlinear media[END_REF], a periodic layered medium is considered. The basis vector e 1 which orientates the x-axis is normal to the interfaces. The layers have the same thickness d = 1 cm and perfect contact is assumed between neighbor layers (continuity of the displacement and of the stress). As illustrated in Fig. 5, the linear material properties Π ∈ {ρ 0 , λ, µ} vary in space according to

Π(x, y) =        3 2 Π if ⌊x/d + 0.5⌋ ≡ 0 (mod 2), 1 2 Π otherwise, (30)
where Π ∈ { ρ0 , λ, μ} denotes the reference parameters given in Table 1. The nonlinear material properties l, m, γ, τ 1 do not vary in space. A velocity pulse with frequency f c = 10 kHz and amplitude V = 1 m/s is imposed at the boundary x = 0:

v 1 (0, y, t) = V sin(π f c t) 2 1 0 t 1/ f c . (31)
Since the configuration is invariant along y, the problem is one-dimensional. The wavelength of compression waves c eff P / f c ≈ 39 cm ≫ d is deduced from the effective sound speed

c eff P = √ 3 2 λ + 2 μ ρ0 ≈ 3861 m/s (32)
in the linear layered medium. The numerical domain defined by L x = 1 m is discretized using N x = 3000 points, i.e. each layer is represented by 30 points. The boundary condition [START_REF] Shu | High order weighted essentially nonoscillatory schemes for convection dominated problems[END_REF] is implemented accordingly to Sec. 7.3.4 of [START_REF] Leveque | Finite Volume Methods for Hyperbolic Problems[END_REF] up to t = 1/ f c . For t > 1/ f c , periodic boundary conditions are applied (see Sec. 7.1 of [START_REF] Leveque | Finite Volume Methods for Hyperbolic Problems[END_REF]). A modification of the numerical method described in Sec. 3 is introduced to account for spatially-varying coefficients (Appendix B). Figure 6 illustrates the effect of spatially-varying coefficients on the numerical solution up to t = 2.5 ms. Figs. 6a-6b display seismograms obtained by unwrapping the velocity signals v 1 (0, y, t) recorded at the abscissa x = 0. The "offset" corresponds to the effective propagation distance from the abscissa x = 0 where the source (31) is imposed. Fig. 6a is obtained with the Murnaghan parameters l, m in Table 1 while the softening variable g is equal to zero (no softening: τ 1 → +∞). Wavefront steepening is observed, leading to a series of oscillations. Fig. 6b shows the same output when the softening variable g evolves according to [START_REF] Holzapfel | Nonlinear Solid Mechanics: A Continuum Approach for Engineering[END_REF], with the parameters γ = 10 5 J m -3 and τ 1 = 5.0 J m -3 s. The parameters γ, τ 1 are chosen such that the variable g reaches values of 3% during the simulation, while the characteristic time τ 1 /γ = 5 × 10 -5 s is the same as in Table 1. With respect to the case without softening (Fig. 6a), a longer propagation distance is needed until the wave separates into a series of oscillations, and the amplitudes of oscillations are modified.

The speed of the largest oscillations over the range of propagation distances [START_REF] Haupert | Dynamic acousto-elastic testing[END_REF][START_REF] Holzapfel | Nonlinear Solid Mechanics: A Continuum Approach for Engineering[END_REF] m is deduced from Figs. 6a-6b and reported in Fig. 6c, where the abscissa is the amplitude of each oscillation at the propagation distance 7 m. This procedure is repeated with a smaller input amplitude V = 0.9 m/s. A nearly linear evolution of the propagation speed with respect to the amplitude is observed, as is the case for solitary waves. At small amplitudes, the propagation speed is close to the effective speed of sound c eff P (32). Similar observations are reported in [START_REF] Leveque | Solitary waves in layered nonlinear media[END_REF] where a different constitutive law is used. Numerically, Fig. 6c shows that the softening does not suppress the solitary waves, but it modifies the relationship between their amplitude and their propagation speed.

Conclusion

Within the Lagrangian finite-strain theory, the constitutive model used in this study expresses the stress as a function of the Green-Lagrange strain tensor and a softening variable g. Also, an evolution equation for g is provided. The system of partial differential equations so-obtained writes as a nonlinear hyperbolic system of balance laws, so that finite-volume methods can be applied. If the softening is neglected, then the material follows Murnaghan's law, where a shear wave excitation induces the propagation of smaller-amplitude compression waves. Otherwise, the propagation of a perturbation is responsible for the softening of the material, which recovers gradually its initial stiffness after excitation has stopped. Several numerical examples are considered along the paper. The latter have been chosen to serve modeling purposes rather than experimental ones, related to the observation of hysteresis and long-time relaxation in geophysics and non-destructive testing. In a periodic layered medium, solitary waves are numerically observed, as known in the case of nonlinear elastodynamics. Now, let us mention possible future works. In the Lagrangian plane-strain case, the numerical method presented here can be used for various hyperelastic constitutive models, such as neo-Hookean, Mooney-Rivlin, Ogden, etc. [START_REF] Ogden | Non-linear Elastic Deformations[END_REF][START_REF] Holzapfel | Nonlinear Solid Mechanics: A Continuum Approach for Engineering[END_REF]. For such models, it would be interesting to solve analytically the symmetric shear impact problem in Fig. 2 for validation purposes. Variations of the evolution equation of g in ( 8) can be considered as well. The influence of the softening on the smoothness of solutions could be investigated numerically and theoretically. To go further towards realistic configurations, viscoelastic attenuation should be accounted for, e.g. in a similar fashion to [START_REF] Berjamin | Modeling longitudinal wave propagation in nonlinear viscoelastic solids with softening[END_REF]. Concerning the nonlinear layered medium, it would be interesting to derive the corresponding nonlinear dispersive wave equations by homogenization, as done in [START_REF] Leveque | Solitary waves in layered nonlinear media[END_REF][START_REF] Samsonov | Strain Solitons in Solids and How to Construct Them[END_REF][START_REF] Andrianov | Dynamic homogenization and wave propagation in a nonlinear 1D composite material[END_REF] and related works. For the development of higher-order methods such as ENO or WENO schemes [START_REF] Shu | High order weighted essentially nonoscillatory schemes for convection dominated problems[END_REF], the eigendecomposition of the fluxes provided in the Appendix A is a useful result.

the coordinates of the strain tensor (4) are written E i j = ε i j + 1 2 Θ u p,i u p, j , where ε i j = 1 2 (u i, j + u j,i ) are the coordinates of the infinitesimal strain tensor. Moreover, the components of the Piola-Kirchhoff stress tensor ( 14) are written P i j = (1g) δ im + Θ u i,m α0 δ m j + α1 E m j , i.e.

P i j = (1 -g) α0 δ i j + α1 E i j + Θ α0 u i, j + α1 u i,m E m j , (A.1)
where α0 , α1 are given in [START_REF] Payan | Determination of third order elastic constants in a complex solid applying coda wave interferometry[END_REF]. The coefficients [START_REF] Berjamin | A finite-volume approach to 1D nonlinear elastic waves: Application to slow dynamics[END_REF] of the Jacobian matrices f ′ (q) and g ′ (q) satisfy ρ 0 Q i jkℓ = ∂P i j /∂u k,ℓ and ρ 0 G i j = ∂P i j /∂g. In the present case of Murnaghan material with softening (A.1), one has

ρ 0 Q i jkℓ = (1 -g) δ i j ∂ α0 ∂u k,ℓ + E i j ∂ α1 ∂u k,ℓ + α1 ∂E i j ∂u k,ℓ + Θ (1 -g) u i, j ∂ α0 ∂u k,ℓ + α0 δ ik δ jℓ + u i,m E m j ∂ α1 ∂u k,ℓ + α1 δ ik E jℓ + u i,m ∂E m j ∂u k,ℓ , ρ 0 G i j = -α0 δ i j + α1 E i j + Θ α0 u i, j + α1 u i,m E m j , (A.2)
where

∂E i j ∂u k,ℓ = 1 2 δ ik δ jℓ + δ jk δ iℓ + 1 2 Θ u k,i δ jℓ + u k, j δ iℓ , ∂E nn ∂u k,ℓ = δ kℓ + Θ u k,ℓ , ∂ α0 ∂u k,ℓ = (λ + 2(l -m)E mm ) ∂E nn ∂u k,ℓ + 2mE i j ∂E i j ∂u k,ℓ , ∂ α1 ∂u k,ℓ = 2m ∂E nn ∂u k,ℓ .
The case of Hookean solids is recovered if g ≡ 0, geometric nonlinearity is neglected (Θ = 0), and the Murnaghan coefficients l, m are zero. In this case, Eq. (A.2) gives ρ 0 Q i jkℓ = λδ i j δ kℓ + µ(δ ik δ jℓ + δ jk δ iℓ ).

Appendix A.2. Eigendecomposition

We provide an eigendecomposition of the Jacobian matrices f ′ (q) and g ′ (q) of the fluxes. The hyperelastic case without softening is recovered by removing the last row and the last column of each matrix in the following paragraphs.

Flux along the x-axis. The Jacobian matrix (17) of f at the linear average ( 23) is diagonalized. Let us write A i+1/2, j = PΛP -1 where P is an invertible real matrix, and Λ is a diagonal real matrix. The matrix of eigenvalues Λ = diag(-c P , c P , -c S , c S , 0, 0, 0) satisfies

c P,S = 1 √ 2 Q 1111 + Q 2121 ± (Q 1111 -Q 2121 ) 2 + 4 Q 1121 Q 2111 , (A.3)
where the plus sign gives the expression of c P (compressional waves), and the minus sign gives the expression of c S (shear waves). The first four right eigenvectors p k i+1/2, j of A i+1/2, j used in ( 24)-( 25) are the first four columns of P, where

P =                                      p
0 0 0 0 0 0 0 1      , (A.4)
with the coefficients

p 11 = 1/c P , p 31 = Q 2111 /c P (c P ) 2 -Q 2121 , p 13 = - (c P ) 2 Q 1121 /c S (c P ) 2 -Q 2121 Q 1111 + Q 1121 Q 2111 , p 33 = 1/c S , p 15 = Q 1112 Q 2121 -Q 1121 Q 2112 Q 1111 Q 2112 -Q 1112 Q 2111 , p 25 = Q 1121 Q 2111 -Q 1111 Q 2121 Q 1111 Q 2112 -Q 1112 Q 2111 , p 16 = Q 1112 Q 2122 -Q 1122 Q 2112 Q 1111 Q 2112 -Q 1112 Q 2111 , p 26 = Q 1122 Q 2111 -Q 1111 Q 2122 Q 1111 Q 2112 -Q 1112 Q 2111 , p 17 = G 21 Q 1112 -G 11 Q 2112 Q 1111 Q 2112 -Q 1112 Q 2111 , p 27 = G 11 Q 2111 -G 21 Q 1111 Q 1111 Q 2112 -Q 1112 Q 2111 .
The matrix P is invertible provided that its determinant is nonzero, i.e.

Q 1121 Q 2111 Q 1111 Q 2121 and Q 1121 Q 2111 -1 4 (Q 1111 -Q 2121 ) 2 .
Let us consider each equality case:

• if Q 1121 Q 2111 = Q 1111 Q 2121 ,
then the eigenvalues of f ′ (q) satisfy c S = 0. Therefore, the reduced system of conservation laws for plane waves propagating along x is not strictly hyperbolic (eigenvalues {-c P , c P , -c S , c S , 0});      q 11 q 12 q 13 q 14 q 15 q 16 q 17 -q 11 -q 12 -q 13 -q 14 q 15 q 16 -q 17 -q 41 -q 42 -q 43 -q 44 q 45 q 46 -q 47 q 41 q 42 q 43 q 44 q 45 q 46 q 47 0 1/p 25 0 -p 26 /p 25 0 0 -p 27 /p 25 0 0 0 1 0 0 0 0 0 0 0 0 0 1 The coefficients α k i+1/2, j in ( 24) are equal to the scalar products α k i+1/2, j = l k i+1/2, j • q n i+1, jq n i, j .

q 12 = 1
Flux along the y-axis. Similarly to [START_REF] Favrie | A thermodynamically compatible splitting procedure in hyperelasticity[END_REF], we introduce the Jacobian matrix g ′ (q) of g at the linear average B i, j+1/2 = g ′ 1 2 (q n i, j + q n i, j+1 ) , and provide an eigendecomposition B i, j+1/2 = PΛP -1 . The matrix of eigenvalues Λ = diag(-c P , c P , -c S , c S , 0, 0, 0) satisfies

c P,S = 1 √ 2 Q 2222 + Q 1212 ± (Q 2222 -Q 1212 ) 2 + 4 Q 2212 Q 1222 , (A.6)
where the plus and minus signs give the expressions of c P and c S , respectively. With similar notations as (A.4), we have

P =                                      0
p 15 = Q 1212 Q 2222 -Q 1222 Q 2212 Q 1211 Q 2212 -Q 1212 Q 2211 , p 25 = Q 1222 Q 2211 -Q 1211 Q 2222 Q 1211 Q 2212 -Q 1212 Q 2211 , p 16 = Q 1212 Q 2221 -Q 1221 Q 2212 Q 1211 Q 2212 -Q 1212 Q 2211 , p 26 = Q 1221 Q 2211 -Q 1211 Q 2221 Q 1211 Q 2212 -Q 1212 Q 2211 , p 17 = G 22 Q 1212 -G 12 Q 2212 Q 1211 Q 2212 -Q 1212 Q 2211 , p 27 = G 12 Q 2211 -G 22 Q 1211 Q 1211 Q 2212 -Q 1212 Q 2211 .      q 11
q 12 q 13 q 14 q 15 q 16 q 17 -q 11 -q 12 -q 13 -q 14 q 15 q 16 -q 17 -q 41 -q 42 -q 43 -q 44 q 45 q 46 -q 47 q 41 q 42 q 43 q 44 q 45 q 46 q 47 1/p 25 0 -p 16 /p 15 0 0 0 -p 17 /p 15 0 0 1 0 0 0 0 0 0 0 0 0 0 1

Figure 1 :

 1 Figure 1: Murnaghan hyperelasticity (g ≡ 0). Relative variation ∆c/c of the speeds c P and c S of compressional and shear waves with respect to the case of Hooke's law. The variation ∆c/c is represented with respect to the compression strain u 1,1 , and the shear strain is u 2,1 = 0. The dotted lines mark Taylor series approximations[START_REF] Leveque | Finite Volume Methods for Hyperbolic Problems[END_REF].

Figure 2 :

 2 Figure 2: Generation of compression waves from pure shear initial data in Murnaghan hyperelasticity (g ≡ 0). Map of W 1/8 at t = 0.015 ms, where W is the strain energy density (J/m 3 ) deduced from the numerical solution. The velocity amplitude of the impact problem (27) is V = 0.1 m/s.

Figure 3 :

 3 Figure 3: Rotated longitudinal and transverse velocities (v 1 ) ϕ , (v 2 ) ϕ along the line y = 0.2 m in Fig. 2. (a) Normalized velocity (v 2 ) ϕ /V in the linear and nonlinear cases for V = 0.1 m/s. (b) Normalized velocity (v 1 ) ϕ /V in the nonlinear case for various amplitudes V.

Figure 4 :

 4 Figure 4: Softening induced by an acoustic point source (29). (a) Map of the strain energy (1g) W in J/m 3 at t = 0.04 ms, where the forcing amplitude is A v = 0.5 kN/m. (b)-(c) Time histories of the velocity component v 1 at the position of the receivers R1-R2 (bullet point and square on the map) for several forcing amplitudes (top); Same for the softening variable g (bottom).

Figure 6 :

 6 Figure6: (a) Seismograms of the particle velocity v 1 for a propagating sinusoidal pulse in a periodic layered Murnaghan material (g ≡ 0). (b) Same output in the case of a periodic layered material with softening (g 0). (c) In each case, we represent the observed speed of the largest oscillations over the range of propagation distances[START_REF] Haupert | Dynamic acousto-elastic testing[END_REF][START_REF] Holzapfel | Nonlinear Solid Mechanics: A Continuum Approach for Engineering[END_REF] m with respect to their amplitude at the propagation distance 7 m.

Table 1 :

 1 Physical parameters of concrete.

  2 p 13p 15 p 33 p 25 (p 11 p 33p 13 p 31 ) , q 42 = 1 2 p 11p 15 p 31 p 25 (p 11 p 33p 13 p 31 ) , 26p 16 p 25 ) p 33p 13 p 26 p 25 (p 11 p 33p 13 p 31 ) , q 44 = 1 2 (p 15 p 26p 16 p 25 ) p 31p 11 p 26 p 25 (p 11 p 33p 13 p 31 ) , q 15 = 1 2 p 11 p 33 p 11 p 33p 13 p 31 , q 45 = -1 2 p 31 p 33 p 11 p 33p 13 p 31 , q 16 = -1 2 p 13 p 11 p 11 p 33p 13 p 31 , q 46 = 1 2 p 11 p 33 p 11 p 33p 13 p 31 , q 17 = 1 2 (p 15 p 27p 17 p 25 ) p 33p 13 p 27 p 25 (p 11 p 33p 13 p 31 ) , q 47 = 1 2 (p 15 p 27p 17 p 25 ) p 31p 11 p 27 p 25 (p 11 p 33p 13 p 31 ) .

	q 13 = -2 q 14 = 1	1 2 (p 15 p p 11 p 33 -p 13 p 31 p 13	,	q 43 = -	1 2	p 11 p 11 p 33 -p 13 p 31	,

  S ) 2 Q 1222 /c P (c P ) 2 -Q 1212 Q 1212 -Q 1222 Q 2212 , p 41 = 1/c P , p 23 = 1/c S , p 43 = -Q 2212 /c S (c P ) 2 -Q 1212 ,

		0	0	0	p 15 p 16 p 17
	p 21 0	-p 21 0	p 23 0	-p 23 0	p 25 p 26 p 27 0 1 0
	p 41 p 21 /p 41 p 21 /p 41 -p 41	p 43 1	1 -p 43	0 1	0 0	0 0	,	(A.7)
	1 0	1 0	p 43 /p 23 p 43 /p 23 0 0 0 0	0 0	0 1	    
	with the coefficients
	p 21 =	(c

Acknowledgments

This work was supported by the interdisciplinary mission of CNRS (INFINITI). The project leading to this publication has received funding from Excellence Initiative of Aix-Marseille University -A*MIDEX, a French "Investissements d'Avenir" programme. It has been carried out in the framework of the Labex MEC. Opération réalisée avec le concours du Programme d'investissements d'avenir du Gouvernement franc ¸ais dont la gestion a été confiée à l'Andra.

Appendix A. Jacobian matrices of the flux Appendix A. [START_REF] Johnson | Resonance and elastic nonlinear phenomena in rock[END_REF]

. Expression of the coefficients

We use the Einstein summation convention with indices in {1, 2}. To encompass both cases with and without geometric nonlinearity in a single equation, we introduce a parameter Θ ∈ {0, 1} such that Θ = 1 corresponds to finite strain and Θ = 0 corresponds to infinitesimal strain. Hence, (A.8) with the coefficients

Appendix B. Spatially-varying coefficients

The problem described in Fig. 5 is one-dimensional. Thus, the variable y and the index j do not appear in the present description of the numerical method, which is a modification of (25) to account for spatially-varying coefficients [START_REF] Fogarty | High-resolution finite-volume methods for acoustic waves in periodic and random media[END_REF][START_REF] Leveque | Finite-volume methods for non-linear elasticity in heterogeneous media[END_REF][START_REF] Leveque | Solitary waves in layered nonlinear media[END_REF] (cf. Chap. 9 of [START_REF] Leveque | Finite Volume Methods for Hyperbolic Problems[END_REF]). The conservation laws [START_REF] Vakhnenko | Strain-induced kinetics of intergrain defects as the mechanism of slow dynamics in the nonlinear resonant response of humid sandstone bars[END_REF] are rewritten in terms of q = (u 1,1 , ρ 0 v 1 , g) ⊤ . The modified method amounts to choosing the eigenvalues {-c P , c P , 0} and corresponding eigenvectors p k i+1/2, j , k ∈ {1, 3} of A i+1/2 in a downwind fashion. Hence, c P is computed at the abscissa x i for k = 1 and at the abscissa x i+1 for k = 2. Since both left-going and right-going waves do not have equal absolute speeds at the cell interface x i+1/2 , writing the low-order part of the algorithm (25) as a flux difference is no longer possible: the "wave-propagation form" is used instead. To guarantee the continuity of the velocity and of the stress at the material interfaces, a modified expression of the coefficients α k i+1/2 defining the vectors W k i+1/2 is used. In particular, the sum of these vectors does not equal the jump q i+1q i anymore, as was the case in a homogeneous medium [START_REF] Remillieux | Propagation of a finite-amplitude elastic pulse in a bar of berea sandstone: A detailed look at the mechanisms of classical nonlinearity, hysteresis, and nonequilibrium dynamics[END_REF]. To avoid numerical instability, a transmission-based limiter is implemented, which modifies the expression (26) of the coefficients θ k i+1/2 .


Ibrahim Busuladzic

Fabrice M Schmidt

Pierre G Lafleur

Distributive Mixing in a Single Screw Extruder Using Boundary Element Method

The relationship between the mixing performance in extrusion and the process conditions is of great importance but is never easy to establish. Although experimental prototyping of mixing equipment provides insight into the mixing process, it is too time consuming and cost prohibitive to fully optimize the geometry and the process itself. Alternatively, numerical methods can be used to provide more detailed information on the mixing flows to alleviate the problems with experimental modeling.

The Boundary Element Method is an alternative technique which offers the possibility of modeling complex systems since it requires only the discretisation of the model boundaries. In this paper, we examine the use of a Boundary Element Method (BEM) model to predict the dispersive mixing performance of Maddock-like mixing heads on single screw extruders. We use this method to calculate the flow number on the middle plan of the studied geometries and compare it between different conditions of extrusion.

The experiments were conducted with a classical Maddock mixing head (in-house design). Two materials were extruded: a polypropylene/ calcium carbonate blend and a polyvinyl chloride/ red pigment blend. The numerical simulations were compared to the observed mixing behavior.

1-Introduction

Mixing is a very important process when we need final products with uniform properties. The state of dispersion (relative size of agglomerates) of filler minerals can be directly correlated to the dynamic and mechanical properties of composites (1)(2). Such as in every polymer process, the mixing is a vital operation in single screw extrusion. Two types of mixing exist: Dispersive mixing and Distributive mixing.

Dispersive mixing involves the rupture of clumps and agglomerates of solid particles such as calcium carbonate or pigments. This rupture will happen only when the forces applied to the agglomerate or droplet exceed the forces keeping it together. In single screw extrusion, we use some devices to improve the mixing performance, like mixing heads, because dispersive mixing depends on the strain rates applied by the mixing devices. To realize our experiments we used Maddock like mixing heads. Classical Maddock mixing heads are primarily dispersive mixers. The material traveling through it is forced over the narrow flight gap before exiting the extruder. The other type of mixing, distributive mixing, spreads the dispersed or secondary phase throughout the primary phase or the mixture.

Beside the experimental studies to characterize mixing numerical analysis can be conducted. Various numerical methods have been used to simulate viscous creeping flows. In this numerical study we use The Boundary Element Method (BEM). The main attractiveness of BEM is that, for many applications, it provides the ability to reduce the dimensionality of the problem by one since only the boundary unknowns calculated in order to get the whole solution. Given the solution on the boundary, the unknowns inside the domain can be obtained.

2.1-Experimental study

The experiments were carried out on a 45 mm extruder. Two types of screw configuration are used: standard single screw and single screw with a Maddock mixing head Figure 1. The first test was to incorporate 5% of calcium carbonate into polypropylene. Then an extruded blown film was produced. The extruded film is transparent and can be observed with a microscope, see Figure 2. The dispersion is evaluated on 20 µm thick film. For the same conditions of extrusion we observed 30 samples from different parts of the film. In order to characterize the dispersion we used the model developed by Bories (3) to calculate the dispersion index f. f can be greater then unity. Indeed, agglomerates are constituted of mineral, entrapped air and polymer, so their density is lower then the mineral's density. The important extrusion parameters are: rotational speed, pressure and screw geometry. The average barrel temperature is 195°C. Rotational speeds are 15, 20 and 25 rpm. The results are shown in Figure 3. It can be noted that the effect of the mixing head is more important when we increase the rotational speed. Also, applying higher pressure leads to better dispersion states. We note an important effect of the mixing head on the quality of the product ( i.e. dispersion of the red pigments).

2.2-Numerical study

The Boundary Element Method (BEM) offers the possibility of modeling complex systems since it requires only discretization of the model boundaries [START_REF] Brebbia | Boundary elements An introductory course[END_REF]. BEM is also a convenient technique for modeling the rotating boundary problems such as the extrusion process because the domain only needs to be discretized once. The main disadvantage of BEM is the treatment of the nonlinear problems. Where the method loses its principal's advantages versus other numerical methods. Some techniques to deal with this problem exist such as Cell integration method or Dual reciprocity method [START_REF] Davis | Investigation of non-linear flows in polymer mixing using the boundary integral method[END_REF].

In this study the polymer is assumed to be a Newtonian fluid and the flow process is considered isothermal. The momentum balance and continuity equation can be written as:

0    k jk x  (1) 0    k k x v (2)
Where jk  is the total stress tensor and k u is the velocity vector. The stress tensor can be written as:

jk jk jk v p    2    (3)
Where p is the pressure, jk  is the Kronecker symbol,  is Newtonian viscosity and jk v is the rate of strain tensor. Finally we obtain the classical boundary integral equation.

Γ Γ Γ Γ d v T d v T v c jk k k jk i k i (4)
Where the multiplier i c depends on the location of i and takes the value: T are the velocity and traction components of fundamental solution respectively. The fundamental solutions correspond to the solutions in the j direction at a field point resulting from an infinite force applied in the k direction to the source point. In the case of viscous fluids for 3D problems we can write:

  k j jk jk r r r u , ,      8 1 (5)   n k j jk r r r r T , , , 2 4 1 
  (6) r r r j j  , ( 7
)
Where r is the distance between the source points and the points of calculation, i r , . In this study we meshed the boundary using two dimensional quadrilateral isoparametric quadratic elements. We proceed to the BEM model assessment with the known analytical solution for a Couette flow. Rotational speed applied on the internal cylinder is 40rad/min. For this simulation we used 400 quadratic elements. CPU time is 1002 s. The total error comparing those two solutions for velocity are less then 0.25% and the same for velocity gradient are less then 3%, Figure 5. Numerical evaluation of mixing was done calculating the Flow number [START_REF] Gramann | Evaluating mixing of polymers in complex systems using the boundary integral method[END_REF][START_REF] Cheng | [END_REF]. The Flow number is a non-dimensional quantity that could be used to asses the flow type. The flow number is defined as:

        (8)
Where   is the magnitude of the strain rate,  is the magnitude of the vorticity tensor. The maximum value for the flow number is 1 which represents an ideal elongational flow shown in Figure 6. Figure 7 shows the first simulation for the case of extrusion using single screw without mixing devices. The boundary conditions applied were simplified. The barrel is rotating and the screw is stationary (8). 802 quadratic elements were used for this simulation. We placed the internal points at the middle plane of channels where we calculate the flow numbers. CPU time is 140 min on IBM RS6000 station. In the case when we apply the rotating velocity of 20 rpm and the pressure difference of 300 kPa, one can see that the value of flow number is not more then 0.62. The average value of flow number is 0.54. The geometry of a classical Maddock mixing head shown in Figure 1. can be divided into three identical section. Because the sections are identical the simulation is performed on two channels, as shown in Figure 8. We apply the following boundary conditions. The barrel is rotating and the mixing head is stationary. Because the flow rate is constant, we impose the known velocity of the polymer at the entrance and in the exit of the mixing head. The polymer is entering by the first channel, forced through the narrow gap and exiting by the second channel, Figure 8b. For this simulation we used 791 quadratic elements. 126 internal points were used and placed at the middle plane inside the geometry to calculate the flow number. For those simulations we worked with a constant flow rate and we used different rotational velocities to study the influence of this parameter on the flow number. CPU time is 125 min on IBM RS6000 station.

For the first simulation we worked with flow rate of 210g/min and we were applying the following rotational velocities a) -15 rpm, b) -40 rpm and c) -25 rpm. The results are shown in Figure 9. We note that the flow number increases on the points entering the narrow gap. Inside this zone, the mixing head produces an elongational type of flow such as shown in Figure 6. The elongational deformation is more important then in others parts of mixing head. By increasing the rotational velocity, the flow number is increasing in the zone entering the narrow gap. Figure 10 shows the flow number distribution where we used a constant flow rate of 400g/min. and different rotational velocities. In this case the factor of pressure difference is more important. With regard to the previous case, the average flow number is greater but is slightly weaker in the region of narrow entrance. Increasing rotational velocity yields higher flow numbers as in the previous case.

3-Conclusions

Experimental study on filler dispersions show that mixing heads are more efficient when we increase the rotational speed. Also increasing the pressure yields to better dispersions. For single screw without mixing head, increasing the rotational speed will result in a decrease of the quality of dispersion. A pressure decrease inside of extruder resulted in worse dispersion. PVC's experiments show the necessities of using a mixing head. The rotational speed is a very important factor to get better dispersions, especially for the screw without mixing head.

Even though in this study the polymer is considered Newtonian, Boundary element method give us very quickly some useful information about the quality of mixing. For this study we have evaluated the dispersion calculating the flow number. Simulation for mixing heads shows that the presence of the elongational mode of the flow is important when the polymer is in the zone of the entrance in the narrow gap. Increasing the rotational speed will increase weakly the flow number. Also the effect of the flow rate caused by pressure difference is evident, producing higher flow number inside the overall geometry.

The numerical results highlight the dependence of the dispersion to pressure difference and to the rotational speed and are in god agreement with experimental results. These results are preliminary and a more elaborate dispersive model will be used to try to obtain quantitative comparison between the state of dispersion and the simulation.

Figure 1 .

 1 Figure 1. Geometry of spiral Madock head

Figure 2

 2 Figure 2 -Sample of

n

  is the number of particles of diameter i d . s A is the total surface observed. f W is the film weight per unit surface, f  is the film density and v  is the mineral volume fraction. This dispersion index f represents the volume ratio of agglomerates larger then 15μm and all secondary phase.

Figure3.

  Figure3. Measured dispersion index for different conditions of extrusion

Figure 4 .

 4 Figure 4. Experimental results quantifying dispersion of red pigment

  r in x, y and z direction. Knowing the external solutions we can obtain the solutions in every point of the domain (internal solutions) applying1  i c

Figure 5 .

 5 Figure 5.Validation of the BEM simulation for Couette flow calculating the angular velocity error [%]

Figure 6 .

 6 Figure 6. Elongational flow

Figure 7 -

 7 Figure 7 -Results for the flow number for the case of extrusion with single screw.

Figure 8 .

 8 Figure 8.-a) Geometry of mixing head, b) Part of the geometry to model and the boundary conditions.

Figure 9 -Figure 10 -

 910 Figure 9 -Flow number distribution with flow rate of 210 g/min and different rotational velocities: a) -30 rpm, b) -40 rpm, c) -50 rpm.

      