I Murrieta-Pazos

C Gaiani

L Galet

B Cuq

S Desobry

J Scher

Comparative study of particle structure evolution during water sorption: Skim and whole milk powders

Keywords: Dairy powder Sorption isotherm Microstructure Surface composition AFM, atomic force microscopy, aw, water activity, BET, Brunauer-Emmett-Tellet, DVS, dynamic vapor sorption, EDX, energy dispersive X-ray, Ra, mean roughness, RH, relative humidity, SEM, scanning electron microscopy, SMP, skim milk powder, WMP, whole milk powder, XPS, X-ray photoelectron spectroscopy

Surface composition of dairy powders influences significantly a quantity of functional properties such as rehydration, caking, agglomeration. Nevertheless, the kinetic of water uptake by the powders was never directly related to the structure and the composition of the surface. In this work, the effect of relative humidity on the structural reorganization of two types of dairy powder was studied. The water-powder interaction for industrial whole milk powder, and skim milk powder was studied using dynamic vapor sorption. The water sorption isotherms were fitted with a Brunner-Emmet-Teller model and each stage of the sorption curve was analyzed with a Fickian diffusion. The water content in the monolayer predicted for each powder and the moisture diffusivity calculated were discussed and compared. Concurrently, powders microstructure and powders surface under variable relative humidity were assessed by X-ray photoelectron spectroscopy, scanning electron microscopy coupled with energy dispersive X-ray and atomic force microscopy. A correlation between the data obtained from the sorption isotherms and the modifications of structure allowed us to conclude that powder microstructure and chemical state of the components could play an important role in determining the water diffusivity.

Introduction

One of the most important physical phenomena occurring in dairy powder during milk reconstitution [START_REF] Gaiani | Use of a turbidity sensor to characterize micellar casein powder rehydration: Influence of some technological effects[END_REF][START_REF] Foster | The prediction of moisture sorption isotherms for dairy powders[END_REF] or humid storage [START_REF] Gaiani | Native phosphocaseinate powder during storage: lipids released onto the surface[END_REF][START_REF] Haque | Protein conformational modifications and kinetics of water-protein interactions in milk protein concentrate powder upon aging: effect on solubility[END_REF][START_REF] Gonzales | Maillard reaction kinetics in milk powder: effect of water activity at mild temperatures[END_REF][START_REF] Silalai | Roles of water and solids composition in the control of glass transition and stickiness of milk powders[END_REF] is water transfer. During reconstitution, water transfer optimization is fundamental to obtain a quick process. This depends on the succession of different steps: the wettability, the sinkability, the dispersibility; and finally the solubility, which is the ability of components to dissolve in water [START_REF] Freudig | Dispersion of powders in liquids in a stirred vessel[END_REF][START_REF] Gaiani | The dissolution behaviour of native phosphocaseinate as a function of concentration and temperature using a rheological approach[END_REF][START_REF] Gaiani | Dairy powder rehydration: Influence of protein state, incorporation mode, and agglomeration[END_REF][START_REF] Baldwin | Insolubility of milk powder products -A minireview[END_REF][START_REF] Baldwin | Factors affecting the reconstitution properties of whole milk powder[END_REF]. During storage, diffusion of water inside the powder particles was often studied with the determination of the water sorption isotherms. The case of a complex matrix (e.g. dairy powders) in which the structure is heterogeneous was often investigated [START_REF] Gaiani | Native phosphocaseinate powder during storage: lipids released onto the surface[END_REF][START_REF] Mathlouthi | Water vapour sorption isotherms and the caking of food powders[END_REF][START_REF] Hogan | Influence of milk proteins on the development of lactose-induced stickiness in dairy powders[END_REF][START_REF] Sun | Comparison of the mechanical properties of cellulose and starch films[END_REF][START_REF] Shrestha | Water sorption and glass transition properties of spray dried lactose hydrolysed skim milk powder[END_REF]. Water diffusivity was already studied in detail for complex food structures such as cakes [START_REF] Guillard | Moisture diffusivity in sponge cake as related to porous structure evaluation and moisture content[END_REF], bread crust [START_REF] Meinders | Modeling water sorption dynamics of cellular solid food systems using free volume theory[END_REF][START_REF] Van Nieuwenhuijzen | Water uptake mechanism in crispy bread crust[END_REF], or starches [START_REF] Chivrac | How does water diffuse in starch/montmorillonite nano-biocomposite materials? Carbohydr[END_REF][START_REF] Yu | Determination of the bulk moisture diffusion coefficient for corn starch using an automated water sorption instrument[END_REF][START_REF] Enrione | Sorption and diffusional studies of extruded waxy maize starch-glycerol systems[END_REF]. Model food systems were also investigated [START_REF] Roca | Predicting moisture transfer and shelf-life of multidomain food products[END_REF][START_REF] Roca | Effective moisture diffusivity modelling versus food structure and hygroscopicity[END_REF]. Furthermore, little is known about how the dairy powder structure is modified when the moisture content of the material increases. Microstructure and diffusivity relationships in a porous food system were investigated [START_REF] Guillard | Moisture diffusivity in sponge cake as related to porous structure evaluation and moisture content[END_REF]. It was observed that the level of porosity in the structure was directly related to the water diffusion. Nevertheless, deep studies coupling sorption isotherms and water diffusivity are lacking in the dairy powder field.

Recently, some authors have demonstrated that dairy powder surface influences some functional properties [START_REF] Gaiani | Native phosphocaseinate powder during storage: lipids released onto the surface[END_REF][START_REF] Faldt | Spray-dried whey protein/lactose/soybean oil emulsions. 2. Redispersability, wettability and particle structure[END_REF][START_REF] Gaiani | Surface composition of dairy powders observed by X-ray photoelectron spectroscopy and effects on their rehydration properties[END_REF][START_REF] Gaiani | How surface composition of high milk proteins powders is influenced by spray-drying temperature[END_REF][START_REF] Gaiani | Relationships between dairy powder surface composition and wetting properties during storage: importance of residual lipids[END_REF][START_REF] Kim | Surface characterization of four industrial spray-dried dairy powders in relation to chemical composition, structure and wetting property[END_REF][START_REF] Kim | Effect of surface composition on the flowability of industrial spray-dried dairy powders[END_REF][START_REF] Nijdam | The effect of surface composition on the functional properties of milk powders[END_REF][START_REF] Vega | Invited review: spray-dried dairy and dairy-like emulsions -Compositional considerations[END_REF]. The application of X-ray photoelectron spectroscopy (XPS) was extensively developed by these authors. From the C, O and N percentages, surface contents in protein, lactose and fat were calculated. These surface contents are representative of the extreme surface of the particle (around 5 nm) [START_REF] Briggs | Practical Surface Analysis, Auger and X-ray Photoelectron Spectroscopy[END_REF]. Another approach is also possible with the development of scanning electron microscopy (SEM) coupled with energy dispersive X-ray (EDX). This complementary technique rends possible the surface analysis at a different depth (up to few micrometers for EDX) [START_REF] Kaciulis | Chemical composition of superconducting SmFeAsO doped with fluorine[END_REF]. Nevertheless, this technique was never applied to food powders.

Applicable conclusions on the influence of processes or storage parameters on dairy powder structures can only be evaluated while considering the supra-molecular structure of the material building the particle [START_REF] Palzer | The relation between material properties and supra-molecular structure of water-soluble food solids[END_REF]. For this purpose, industrial skim and whole milk powders (containing 1% and 28% of lipids respectively) were used as different materials. Weights changes at different equilibrating relative humidities (RH) were followed in continue with a dynamic vapor sorption (DVS) equipment. For each relative humidity, diffusivity value was determined using the second Fick's law. Concurrently, the structure was characterized by SEM, XPS, EDX and Atomic Force Microscopy (AFM). Finally, the variations in water sorption and diffusivity were related to milk powders surface and microstructure modifications.

The aim of this paper was to study milk powder/water interactions considering sorption and diffusion in humid air and integrating swelling and crystallization occurring during powder water vapor hydration.

Materials and methods

Materials

Two industrial milk powders obtained from Lactalis (Laval, France) were studied. The following compositions (% w/w) were measured respectively for skim milk powder (SMP): water (3.8), proteins (37.1), lactose (52.0), lipids (1.4) and ashes (5.7) and whole milk powder (WMP): water (3.9), proteins (26.2), lactose (36.7), lipids (27.3) and ashes (5.9). These powders were sieved at 75 m to present a similar particle size. The sieving did not affect the chemical composition of the powders (data not shown).

Samples were equilibrated at 25 Images were recorded using an MFP3D-BIO instrument (Asylum Research Technology, Atomic Force F&E GmbH, Mannheim, Germany). Silicon nitride cantilevers of conical shape were purchased from Veeco (MLCT-AUNM, Veeco Instruments SAS, Dourdan, France), and their spring constants were determined using the thermal calibration method [START_REF] Lévy | Measuring the spring constant of atomic force microscope cantilevers: thermal fluctuations and other methods[END_REF], providing values of ∼10.4 ± 1.7 pN nm -1 . Prior to experiment, the geometry of the tip was systematically controlled using a commercial grid for 3-D visualization (TGT1, NT-MTD Compagny, Moscow, Russia). Curvature of the tip in its extremity was found to lay in the range ∼20-50 nm. Experiments were performed in air at room temperature. For powder immobilization a piece of wax (Tempfix mounting adhesive, Spi Supplies, USA) was placed on a metal disk support and was warmed up until the wax softens (around 70 • C). After the wax softens, the metal was removed from the heater. The powder was gently spread when the surface of the wax just starts to solidify. Finally, the sample was ready for imaging when the thermal wax becomes totally solid To evaluate the surface topology, the mean roughness (Ra) of the powders was calculated according to Eq. (1):

R = 1 N x N y Ny i=1 Ny j=1 z ij -z mean (1)
where z ij is the height of a given pixel, z mean is the average height of the pixel, and N x and N y are the numbers of pixels in the x and y directions.

Scanning electron microscopy (SEM)

The powders were examined with a Hitachi SEM S2500 instrument operating at 10 kV. They were spread at the surface of a sticky plastic circle fixed on a support. Then the samples were covered with gold particles by sputtering (Bio-Rad type SC 502).

Particles size

Analyses were performed using a laser light diffusion granulometer at a wavelength of 632.8 nm (Mastersizer S, Malvern instruments, UK). The Malvern apparatus is equipped with a sample dispersion module which is suited for powders in suspension. For each measure, 0.25 g of the powder was dispersed in 75 ml of ethanol (to obtain a good obscuration and to avoid particles super-position). A monomodal Gaussian distribution was obtained. As a consequence, three measurements of the d 50 were conducted with each powder sample.

Surface composition 2.4.1. X-ray photoelectron spectroscopy

XPS is a quantitative spectroscopic technique that measures the elemental composition up to 5 nm depth [START_REF]Auger and X-ray Photoelectron Spectroscopy[END_REF]. It provides elemental and chemical state data in solid samples. XPS is performed in ultra high vacuum (10 -8 kPa) which may limit the technique in the food fields (liquids). Nevertheless, it is not an inconvenient in the case of food powders. [START_REF] James | Surface structure and composition of fresh and bloomed chocolate analysed using X-ray photoelectron spectroscopy, cryo-scanning electron microscopy and environmental scanning electron microscopy[END_REF] The use of XPS was extended for the last 10 years to the surface composition of dairy powders [START_REF] Gaiani | Surface composition of dairy powders observed by X-ray photoelectron spectroscopy and effects on their rehydration properties[END_REF][START_REF] Faldt | The surface coverage of fat on food powders analyzed by ESCA (electron spectroscopy for chemical analysis)[END_REF]. From the C, O and N percentages, the contents in protein, lactose and fat at surface were calculated with a matrix formula, the elemental composition in the sample is assumed to be a linear combination of pure components conforming the sample. Milk powders are generally composed of fat, proteins and lactose but also vitamins and traces of mineral elements. By using the precedent matrix, only lactose, fat and protein are take into account; others components are neglected [START_REF] Faldt | The surface coverage of fat on food powders analyzed by ESCA (electron spectroscopy for chemical analysis)[END_REF][START_REF] Faldt | The surface composition of spray-dried protein-lactose powders[END_REF] The XPS analyses were carried out with a Kratos Axis Ultra (Kratos Analytical, Manchester, UK) photoelectron spectrometer. The instrument uses a monochromatic Al K␣ X-ray source. The powder samples were attached to the sample holder with a double side conductive tape. The analyzed area was currently about 700 × 300 m. Spectra were analyzed using the Vision software from Kratos (Vision 2.2.0). Quantification was performed using the photoemission cross-sections and the transmission coefficients given in the Vision package.

Energy dispersive X-ray analysis

EDX was carried out on a Philips CM20 microscope. This microscope works at 200 kV and has an Oxford EDX detector with an atmospheric thin window and a link analytical system. An average elemental composition of the samples was obtained by a data collection at 3 different mm 2 sized windows on the powder surface

Dynamic vapor sorption

Sorption isotherms of powders were obtained with a dynamic vapor sorption analyzer DVS (Surface Measurement Systems Ltd., London, UK) equipped with a Cahn microbalance. The experiments were carried out at a constant temperature (20 • C) and different relative humidity values ranging between 0% and 98%. Approximately 15 mg of powder was loaded onto the quartz sample pan. First, the samples were dehydrated in the DVS chamber (RH = 0%) for 200 min. Then, the samples were submitted to a program allowing hydration with 10% RH steps. The samples were considered to be at equilibrium when the value dm/dt (slope of the changing in mass with time) was set to be <0.002 mg min -1 or equilibration time exceeded 1000 min. Experiments were performed in triplicate for each sample.

The Brunauer-Emmett-Tellet (BET) equation has been often used to model water sorption data from 0% to 50% RH [START_REF] Gaiani | Native phosphocaseinate powder during storage: lipids released onto the surface[END_REF][START_REF] Shrestha | Water sorption and glass transition properties of spray dried lactose hydrolysed skim milk powder[END_REF][START_REF] Haque | Water sorption and plasticization behavior of spraydried lactose/protein mixtures[END_REF][START_REF] Haque | Water plasticization and crystallization of lactose in spray-dried lactose/protein mixtures[END_REF][START_REF] Saad | Impact of re-grinding on hydration properties and surface composition of wheat flour[END_REF]. The BET isotherm model is given by:

X = X m C BET a w (1 -a w ) (1 -a w + C BET a w ) (2)
In Eq. ( 2), a w is the water activity, X is the moisture content (g 100 g -1 ), X m is the monolayer moisture capacity (g 100 g -1 ) and C BET is a constant related to the water binding energy. The fitting of BET model was checked by calculating the R 2 with the Origin 6.1 software (OriginLab Corporation, Northampton, USA).

Results and discussion

Surface and size characterization of powders freshly manufactured

Whole milk and skim milk powders were observed by SEM. Due to the particle sieving, the particle size was relatively homogeneous (Figs. 1A and2A). As observed in Figs. 1B and 2B the average size of each sample was 50 m. Nevertheless, the surface microstructure of the powders was totally different depending on the powder composition (Figs. 1C and2C). WMP particles presented a regular surface whereas SMP particles were characterized by a "brain type" surface with some deep and shallow folds. The lactose was in an amorphous state as no sharp crystals were noticed for each powder (and confirmed by DSC, data not shown). All these features were typical for fresh whole and skim milk powders [START_REF] Kim | Surface characterization of four industrial spray-dried dairy powders in relation to chemical composition, structure and wetting property[END_REF].

Concurrently, AFM was performed on the same powders to better analyze the particles surfaces (Fig. 2). Up to now, AFM imaging in tap-mode has never been realized on dairy powders. The AFM images obtained by 3-D projection were in total agreement with those obtained by SEM. In addition, the average surface rough-ness (Ra) was determined and was respectively found around 306 nm for SMP and 146 nm for WMP for a definite surface area (10 × 10 m). The smoother visual aspect of WMP surface in comparison with SMP surface (a "brain" type surface) was associated with the lower Ra value. The surface composition was determined by XPS with the matrix formula develop in detail for milk powders [START_REF] Gaiani | Surface composition of dairy powders observed by X-ray photoelectron spectroscopy and effects on their rehydration properties[END_REF][START_REF] Gaiani | How surface composition of high milk proteins powders is influenced by spray-drying temperature[END_REF][START_REF] Kim | Surface characterization of four industrial spray-dried dairy powders in relation to chemical composition, structure and wetting property[END_REF][START_REF] Faldt | The surface composition of spray-dried protein-lactose powders[END_REF][START_REF] Jayasundera | Surface modification of spray dried food and emulsion powders with surface-active proteins: a review[END_REF][START_REF] Shrestha | Effect of protein concentration on the surface composition, water sorption and glass transition temperature of spray-dried skim milk powders[END_REF]. It was measured that the surface of WMP was composed at 91.8% of fat, 7.2% of proteins and 0.8% of lactose. The SMP surface presented more proteins and lactose at the surface (45.8% and 30.8% respectively) and less fat (22.3%). Trace minerals were measured at the surface of each powder. These surface compositions were in total agreement with others studies done by XPS on milk powders [START_REF] Kim | Surface characterization of four industrial spray-dried dairy powders in relation to chemical composition, structure and wetting property[END_REF][START_REF] Kim | Surface composition of industrial spray-dried milk powders. 3. Changes in the surface composition during long-term storage[END_REF]. The over-representation of fat at the surface in comparison with the bulk composition of the powders is now generally accepted and well understood [START_REF] Gaiani | Native phosphocaseinate powder during storage: lipids released onto the surface[END_REF][START_REF] Gaiani | How surface composition of high milk proteins powders is influenced by spray-drying temperature[END_REF][START_REF] Kentish | Milk skin formation during drying[END_REF][START_REF] Kim | Surface composition of industrial spray-dried milk powders. 2. Effects of spray drying conditions on the surface composition[END_REF]. It is postulated that, during the spray drying of fat-containing materials (e.g. whole milk) larger fat globules are preferentially present at the surface of droplets and so fat appears in high concentration on the powder surface [START_REF] Kim | Surface composition of industrial spray-dried milk powders. 2. Effects of spray drying conditions on the surface composition[END_REF]. In addition, during spray-drying, there is a further increase in the amount of fat probably due to segregation of components within the drying [START_REF] Gaiani | How surface composition of high milk proteins powders is influenced by spray-drying temperature[END_REF]. It is assumed that, even when fat is present at low concentrations (e.g. skim milk) in the milk droplet, residual fat may be preferentially present at the air/droplet interface, and thus may appears at relatively high concentration on the powder surface after complete drying [START_REF] Gaiani | Native phosphocaseinate powder during storage: lipids released onto the surface[END_REF][START_REF] Nijdam | The effect of surface composition on the functional properties of milk powders[END_REF][START_REF] Kim | On the mechanisms of surface formation and the surface compositions of industrial milk powders[END_REF]. This phenomenon may explain the presence of 22% fat at the surface of SMP whereas the particle bulk was composed of only 1.5%. In the presence of surface active components (such as proteins), it appears that this component accumulates at the surface of a milk droplet/particle at the expense of lactose during drying [START_REF] Kim | Surface characterization of four industrial spray-dried dairy powders in relation to chemical composition, structure and wetting property[END_REF][START_REF] Jayasundera | Surface modification of spray dried food and emulsion powders with surface-active proteins: a review[END_REF][START_REF] Shrestha | Effect of protein concentration on the surface composition, water sorption and glass transition temperature of spray-dried skim milk powders[END_REF][START_REF] Özkan | Characterization of stickiness and cake formation in whole and skim milk powders[END_REF]. Indeed, only 30% of lactose was observed at the surface of SMP instead of 54% in the bulk. In complement to XPS, EDX analyses were performed. These two techniques were for the first time coupled to analyze the surface of dairy powders at different deepness. Therefore, EDX analyses were performed to evaluate the surface at 1 m depth whereas XPS was operating at 5-10 nm depth. It was observed in Table 1 that under the surface established by XPS, the WMP presented proteins (25.8%), lactose (13.7%), lipids (58.6%) and minerals (1.9%). In comparison with the first 5 nm, EDX measurements suggested that the surface of WMP may be externally covered by a thin layer of fat (observed by XPS) with a matrix of proteins, lactose and minerals under this layer. For SMP, EDX measures gave the following percentages: proteins (34.0%), lactose (16.8%), lipids (46.6%) and minerals (2.8%). Contrary to WMP, the lipids were more present in the first micrometer layer in comparison with the five first nanometers. These results may be explained by the dominance of proteins at the surface (more tensioactive); more fat being located just under. For each powder, minerals were more located in the first 5 nm.

3.2. Effect of water uptake on powder structure and particle size 3.2.1. Powder structure Each powder was introduced in saturated atmospheres during 4 days allowing water exchanges between the powder of low water activity (a w ) and the relative humidity of the surrounding air. SEM was performed on powders equilibrated at 11 a w ranging between 0.11 and 0.97. For WMP, Fig. 3 indicates that some changes in surface structure with relative humidity occurred. Between a w = 0.11 and 0.33, the surface was relatively smooth and regular. From a w = 0.43, an irregular surface appears corresponding to lactose crystals formation. These crystals were assumed to be localized just under the surface. As reported previously [START_REF] Özkan | Characterization of stickiness and cake formation in whole and skim milk powders[END_REF][START_REF] Faldt | Changes in surface composition of spray-dried food powders due to lactose crystallization[END_REF] and in agreement with the XPS results (Table 1), it was not possible to observe sharp lactose crystals at the surface of WMP certainly due to the fact that the surface was mainly covered with a layer of fat. At higher a w values, lactose crystals were more present but still localized under the surface. Surface evolution of SMP with the relative humidity was totally different (Fig. 4). A substantial amount of needle like lactose crystals was observed on the surface, which confirmed that the lactose was transformed from the rubbery state to the crystalline state around a relative humidity of 54%. Then, the number of dents at the particle surface increased regularly from a w = 0.54 until 0.75. The relationship between the extent of crystallization and relative humidity was found parabolic with the maximum crystallization at a relative humidity around 0.7 [START_REF] Jouppila | Glass transition, water plasticization, and lactose crystallization in skim milk powder[END_REF]. Nevertheless, these assumptions are not obvious on our SEM images.

In freshly spray-dried milk powders (Fig. 1) and at low relative humidity (Figs. 3 and4) lactose was in the metastable amorphous state. In this state, lactose is hygroscopic and forms a continuous matrix containing embedded fat globules, protein, minerals and air vacuoles. In contact to a high relative humidity and/or an increased temperature, the amorphous state present in the powder may progress into a stable crystalline state. Moreover, previous research has shown that crystallization of lactose may be delayed in the presence of proteins [START_REF] Hogan | Influence of milk proteins on the development of lactose-induced stickiness in dairy powders[END_REF][START_REF] Nasirpour | Lactose crystallization delay in model infant foods made with lactose, [beta]-lactoglobulin, and starch[END_REF]. Indeed, preferential sorption of water by non-amorphous constituents was found to delay the rate at which lactose endured a change from the glassy to the rubbery form [START_REF] Hogan | Influence of milk proteins on the development of lactose-induced stickiness in dairy powders[END_REF]. It may also be the case with WMP and SMP (Figs. 3 and4). In SMP, the total protein content was higher (Table 1) and the presence of crystal was observed for slightly higher relative humidity than WMP. Hypothesis for this delay was also supported by the fact that proteins may significantly lower the availability of water for plasticization and crystallization of lactose [START_REF] Hogan | Influence of milk proteins on the development of lactose-induced stickiness in dairy powders[END_REF][START_REF] Shrestha | Effect of protein concentration on the surface composition, water sorption and glass transition temperature of spray-dried skim milk powders[END_REF].

Particle size

Particle size was followed during water uptake for WMP and SMP (Fig. 5). First, a significant size increase was observed for each powder between 11% and 43% RH (Fig. 5). Up to 43% RH, the d 50 increases may be attributed to particle swelling due to water migration into the amorphous molecular matrix. Then, a particle size decrease was found around 54% RH and may be explain by the collapse of the amorphous matrix (lactose/proteins/lipids) triggered by crystallization of lactose. This relative humidity was well correlated with the apparition of crystals on microscopy images (Figs. 3 and4). Finally, the particle size increases strongly from 54% RH. This phase, related to powder caking was already strongly documented [START_REF] Silalai | Roles of water and solids composition in the control of glass transition and stickiness of milk powders[END_REF][START_REF] Mathlouthi | Water vapour sorption isotherms and the caking of food powders[END_REF][START_REF] Hogan | Influence of milk proteins on the development of lactose-induced stickiness in dairy powders[END_REF][START_REF] Shrestha | Water sorption and glass transition properties of spray dried lactose hydrolysed skim milk powder[END_REF][START_REF] Nijdam | The effect of surface composition on the functional properties of milk powders[END_REF][START_REF] Fitzpatrick | Glass transition and the flowability and caking of powders containing amorphous lactose[END_REF]. For WMP, the caking may be mainly attributed to the presence of fat at the surface even if significant caking has been observed principally in powders containing a total fat content of at least 40% [START_REF] Foster | The contribution of milk fat towards the caking of dairy powders[END_REF]. For SMP, lactose may be responsible of the caking phenomena even if lactose under its crystalline form is less sensitive to caking than amorphous lactose [START_REF] Palzer | The relation between material properties and supra-molecular structure of water-soluble food solids[END_REF].

Water sorption isotherms

Kinetic data obtained during water adsorption of SMP and WMP at 20 • C are given in Fig. 6. Extend of hydration was the lowest for whole milk powder particularly at high RH. Indeed, fat in whole milk powder is not a water-absorbing component, consequently fat containing samples exhibit little moisture sorption [START_REF] Palzer | The relation between material properties and supra-molecular structure of water-soluble food solids[END_REF]. If the moisture content was calculated on a non-fat basis, the isotherms of SMP and WMP are almost similar which was also observed by others studies [START_REF] Foster | The prediction of moisture sorption isotherms for dairy powders[END_REF][START_REF] Jouppila | Water sorption and time-dependent phenomena of milk powders[END_REF][START_REF] Lin | Desorption isotherm of milk powders at elevated temperatures and over a wide range of relative humidity[END_REF]. In agreement with electron microscopy images, lactose crystallization was observed. As crystalline form has a lower water capacity than amorphous form [START_REF] Palzer | The relation between material properties and supra-molecular structure of water-soluble food solids[END_REF], water is liberated and is characterized by a decrease in the sorption curve at around 50% RH. The BET monolayer moisture contents for SMP and WMP was calculated at 4.45 and 2.44 g 100 g -1 dry solid respectively (Fig. 7). These values were in agreement with the literature and are dependent of the powder chemical, composition, the temperature [START_REF] Palzer | The relation between material properties and supra-molecular structure of water-soluble food solids[END_REF].

Relationships between moisture diffusivity, powder structure and surfaces

Apparent moisture diffusion (called diffusivity in the continuation of the paper) as a function of relative humidity at 20 • C is presented Fig. 8. This representation is based on the assumption that moisture sorption kinetic is limited by the diffusion. Water diffusivity in whole milk and skim milk powders was calculated for each RH level from the water sorption kinetic curve (Fig. 6) measured using the DVS and the slope method developed by others [START_REF] Yu | Determination of the bulk moisture diffusion coefficient for corn starch using an automated water sorption instrument[END_REF]. For this purpose, different hypothesis were formulated. First, the diffusivity was assumed to be constant for each RH level and could change from one stage to another. Second, particles size is assumed to be homogeneous (average diameter apply to all particles) when placed in the DVS pan. Third, in this work, to improve the method, swelling during adsorption was considered. As indicated by Fig. 5, the powder size was variable upon water sorption and has been taken into account for a more accurate diffusivity determination [START_REF] Chivrac | How does water diffuse in starch/montmorillonite nano-biocomposite materials? Carbohydr[END_REF]. The resultant diffusivity values were plotted as a function of the average water activity of level for all the successive RH levels investigated (Fig. 8). A bell like curve was obtained for each powder in agreement with a great number of studies [START_REF] Chivrac | How does water diffuse in starch/montmorillonite nano-biocomposite materials? Carbohydr[END_REF][START_REF] Yu | Determination of the bulk moisture diffusion coefficient for corn starch using an automated water sorption instrument[END_REF][START_REF] Enrione | Sorption and diffusional studies of extruded waxy maize starch-glycerol systems[END_REF][START_REF] Roca | Predicting moisture transfer and shelf-life of multidomain food products[END_REF][START_REF] Roca | Effective moisture diffusivity modelling versus food structure and hygroscopicity[END_REF]. The maximum diffusivity value was observed at a w = 0.41 for skim milk powder and a w between 0.41 and 0.53 for whole milk powder. D max was respectively at 3.1 × 10 -9 and 1.7 × 10 -9 m 2 s -1 for SMP and WMP. These calculated values were in the range of those reported in literature for food products [START_REF] Guillard | Moisture diffusivity in sponge cake as related to porous structure evaluation and moisture content[END_REF][START_REF] Labuza | Moisture migration and control in multi-domain foods[END_REF]. Comparison with dairy powders was impossible due to lack of data in this field up to now. As observed by others [START_REF] Guillard | Moisture diffusivity in sponge cake as related to porous structure evaluation and moisture content[END_REF], D max values were not related to the moisture content corresponding to the monolayer value (4.45 and 2.44 g 100 g -1 for SMP and WMP respectively) but were slightly higher. Nevertheless, the amorphous-crystallized transition of lac-tose was estimated to be at about a w = 0.43 for WMP and 0.54 for SMP from MEB observation. This corresponds well to the point where the maximum in diffusion rate is observed in Fig. 8. The D max was dependent on the state of the sample [START_REF] Van Nieuwenhuijzen | Water uptake mechanism in crispy bread crust[END_REF] and also on differences in the method used to determined the moisture diffusivity [START_REF] Guillard | Moisture diffusivity in sponge cake as related to porous structure evaluation and moisture content[END_REF][START_REF] Labuza | Moisture migration and control in multi-domain foods[END_REF].

The bell like curves obtained may be explained by different processes occurring during water uptake. Comparison between curves obtained for WMP and SMP (Fig. 8) reveals that at low a w , the moisture diffusivity was lower for WMP. During this first phase, water adsorbs to the surface of the particles. As expected, powders presenting important surface fat coverage (WMP) were less prone to adsorb water and the rate of sorption was lower than lactose containing powders (SMP). Water diffusion is known to be accelerated in polar hydrophilic matrices [START_REF] Palzer | The relation between material properties and supra-molecular structure of water-soluble food solids[END_REF]. In addition, the average surface roughness of SMP particles was significantly higher (Figs. 1 and 2) increasing the effective area. Next to the maximum, the diffusivity decrease may be attributed to change in lactose state. It is well established that the permeability of crystalline and amorphous substance differs significantly. As a consequence, crystalline lactose has a lower water capacity and a lower water diffusion than amorphous lactose [START_REF] Palzer | The relation between material properties and supra-molecular structure of water-soluble food solids[END_REF]. This can be explained by differences in the free volume of these two types of matrices. The low free volume and the small hole size in crystalline structure may be responsible of a lower diffusion. As a consequence, when the quantity of crystallized lactose increases, the diffusivity decreases. Others incriminated the caking of the macroscopic structure [START_REF] Roca | Controlling moisture transport in a cereal porous product by modification of structural or formulation parameters[END_REF] which in turn reduce the effective area exposed to the solvent and therefore decrease the sorption rate or changes in particle porosity [START_REF] Guillard | Moisture diffusivity in sponge cake as related to porous structure evaluation and moisture content[END_REF]. Nevertheless, in this work, the evolution of particle size (swelling, collapse and caking) was taken into account in the mathematical calculation of the moisture diffusivity. Another possible phenomenon may be that at these humidities the transport of water vapor through the air may be a limiting factor. An important and unexpected decrease of the diffusivity value was observed for high relative humidity. Indeed, an additional phenomenon may occur and the Fickian model used to obtain the diffusivity value may be not entirely appropriate to represent this different diffusional mechanism. Indeed, water that migrates into the dry matrix (low RH) will probably do that at a different rate as compared to water migrating into the already hydrated matrix [START_REF] Van Nieuwenhuijzen | Water uptake mechanism in crispy bread crust[END_REF].

Conclusion

Water diffusion in dairy powders was studied coupled with structure and surface investigations. A relationship between lactose state and water diffusivity was highlighted. The differences observed were also related to the surface composition (more or less lipids). In the future, knowledge of these factors can be used to better control powder ageing and rehydration. Moreover, the combination of XPS and EDX, to estimate the surface composition of powders at different depth, with SEM and AFM, to investigate the structure of the powder particles, may be a promising method to better understand the repartition of components under the surface and their structure. An in-depth study is actually under development (with pure milk components) to validate the use of XPS and EDX in the field of dairy powders.

Fig. 1 .

 1 Fig. 1. Scanning electron microscopy of skim milk powder (1) and whole milk powder (2). A: ×50, B: ×1000, C: ×2500.

Fig. 2 .

 2 Fig. 2. Atomic force microscopy images (10 × 10 m) of skim milk powder (1) and whole milk powder (2) at high resolution (A) and the 3D projection (B).

Fig. 3 .

 3 Fig. 3. Surface characterization of whole milk particles observed by scanning electron microscopy at variable relative humidity (from 0.11 to 0.97).

Fig. 4 .

 4 Fig. 4. Surface characterization of skim milk particles observed by scanning electron microscopy at variable relative humidity (from 0.11 to 0.97).

Fig. 5 .

 5 Fig. 5. Evolution of particle size (d50, m) at 11 different aw (from 0.11 to 0.97) for skim and whole milk powders determined with a Malvern apparatus. Mean of three independent analyses.

Fig. 6 .

 6 Fig. 6. Example of kinetic data obtained during water adsorption of milk powders at 20 • C and measured by DVS automatic sorption analyzer (powders sieved at 50 m).

Fig. 7 .

 7 Fig.7. Experimental and predicted water vapor sorption isotherms (mean of three independent analyses). Determination of the monolayer moisture content (Xm) and constant C with the BET model for relative humidity ranging between 0% and 50%.

Fig. 8 .

 8 Fig. 8. Moisture diffusivity of skim milk and whole milk powders as a function of aw. Mean of three independent analyses.

  Water sorption isotherms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 3.4. Relationships between moisture diffusivity, powder structure and surfaces . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 4. Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

	3.3.

Table 1

 1 Composition (%) of milk powders: bulk composition, surface composition calculated from XPS analyses (≈5 nm) and surface composition calculated from EDX (≈1 m). Mean of at least two independent analyses.

		Proteins	Lactose	Lipids	Minerals	Water
	Bulk composition (%)
	Skim milk powder Whole milk powder First 5 nm composition from XPS (%)	34.1 ± 0.2 25.2 ± 0.3	53.9 ± 0.3 38.8 ± 0.1	1.5 ± 0.0 26.2 ± 0.1	6.5 ± 0.0 5.0 ± 0.5	4.0 ± 0.1 3.8 ± 0.1
	Skim milk powder Whole milk powder First 1 m composition from EDX (%)	45.8 ± 1.4 7.2 ± 0.5	30.8 ± 1.2 0.8 ± 0.1	22.3 ± 1.1 91.8 ± 2.1	1.1 ± 0.1 0.2 ± 0.0	--
	Skim milk powder Whole milk powder	34.0 ± 1.8 25.8 ± 0.9	16.8 ± 0.8 13.7 ± 1.3	46.4 ± 2.3 58.6 ± 2.8	2.8 ± 0.6 1.9 ± 0.4	--

Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 Acknowledgments ANR funding from program "Reactive Powder" are gratefully thanks. LCPME laboratory (Laboratoire de Chimie Physique et de Microbiologie pour l'Environnement, Nancy) is thanks for performing XPS and AFM analyses.


H Berjamin

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

G Chiavassa

N Favrie

AFinite-Volume Approachto1DNonlinear Elastic Waves: Application to SlowDynamics

Keywords: no, 43, 25, Dc, 02, 70, Bf

Anumerical method for longitudinal wave propagation in nonlinear elastic solids is presented. Here, we consider polynomial stress-strain relationships, which are widely used in nondestructive evaluation. The large-strain and infinitesimal-strain constitutive laws deduced from Murnaghan'sl aw are detailed, and polynomial expressions are obtained. The Lagrangian equations of motion yield ahyperbolic system of conservation laws. The latter is solved numerically using afi nite-volume method with flux limiters based on Roe linearization. The method is tested on the Riemann problem, which yields nonsmooth solutions. The method is then applied to acontinuum model with one scalar internal variable, accounting for the softening of the material (slowdynamics).

Introduction

Rocks and concrete are known to behave nonlinearly when vibrating longitudinally,evenatvery lowamplitudes [START_REF] Guyer | Nonlinear mesoscopic elasticity: evidence for an ew class of materials[END_REF][START_REF] Tencate | Slowdynamics of earth materials: an experimental overview[END_REF]. Firstly,d ynamic acoustoelastic testing (DAET) [START_REF] Renaud | Revealing highly complexelastic nonlinear (anelastic)behavior of Earth materials applying anew probe: Dynamic acoustoelastic testing[END_REF][START_REF] Rivière | Pump and probe wavesindynamic acoustoelasticity: comprehensive description and comparison with nonlinear elastic theories[END_REF] reveals that the speed of sound measured locally decreases with time, and recovers its initial value after the excitation is stopped. This softening occurs overatime scale larger than the period of the dynamic loading, which highlights the phenomenon of slow dynamics.S econdly,t he evolution of this speed with respect to the strain presents an hysteretical behavior.L astly,a ll these phenomena are enhanced when the forcing amplitude is increased. Such observations are not compatible with linear elastodynamics, where the speed of sound is aconstant.

Polynomial nonlinear stress-strain relationships are widely used in nondestructive testing [START_REF] Va Nden Abeele | Elastic pulsed wave propagation in media with second-or higher-order nonlinearity.P art I. Theoretical framework[END_REF][START_REF] Payan | Applying nonlinear resonant ultrasound spectroscopytoimproving thermal damage assessment in concrete[END_REF]. However, nonlinear elasticity is not sufficient to represent the softening and the hysteresis revealed by DAET.Several models can be found in the literature to reproduce these phenomena (see e.g. [START_REF] Meurer | Wave propagation in nonlinear and hysteretic media-a numerical study[END_REF][START_REF] Li | Pseudo-spectral simulation of 1D nonlinear propagation in heterogeneous elastic media[END_REF]). The soft-ratchet model by Va khnenkoe t al. [START_REF] Va Khnenko | Strain-induced kinetics of intergrain defects as the mechanism of slowdynamics in the nonlinear resonant response of humid sandstone bars[END_REF][START_REF] Favrie | Fast and slowdynamics in anonlinear elastic bar excited by longitudinal vibrations[END_REF] consists in introducing as calar variable g to describe the softening of the material. As imilar model with refinements wasproposed by Lyakhovskyand coauthors in aseries of papers [START_REF] Lyakhovsky | Distributed damage, faulting, and friction[END_REF][START_REF] Lyakhovsky | Nonlinear elastic behaviour of damaged rocks[END_REF]. Recently,a3Dmodel of continuum has been developed by the authors [START_REF] Berjamin | Nonlinear wavesi ns olids with slowd ynamics: an internalvariable model[END_REF] in the framework of continuum thermodynamics with internal variables of state [START_REF] Maugin | Thermodynamics with internal variables[END_REF][START_REF] Maugin | The sagaofinternal variables of state in continuum thermo-mechanics (1893-2013[END_REF]. This model, which generalizes the soft-ratchet model to 3D geometries and fixes thermodynamical issues, will be used along the present paper. Regardless the fact that slowd ynamics is taken into account or not, the equations of the model appear as a nonlinear hyperbolic system of conservation laws with relaxation, the solutions of which may be discontinuous (shock waves).The numerical computation of nonsmooth solutions requires aparticular care. Here, dedicated finitevolume methods [START_REF] Godlewski | Numerical Approximation of Hyperbolic Systems of Conservation Laws[END_REF][START_REF] Leveque | Finite-Volume Methods for Hyperbolic Problems[END_REF][START_REF] Ve Lasco-Segura | Afinite volume approach for the simulation of nonlinear dissipative acoustic wave propagation[END_REF] are used to compute such solutions in an on-oscillatory manner.A1D scheme based on the Roe linearization and flux limiters is adapted to the present system of equations. Due to the separation of timescales between the wave propagation and the softening of the material, the coupling between nonlinear elastodynamics and the slowdynamics is straightforward.

The article is organized as follows. Forpedagogical purposes, the case of nonlinear elastic solids without slowdynamics is first considered (sections 2and 3).Section 2derivest he equations of nonlinear elastodynamics, and the case of Murnaghan hyperelastic material [START_REF] Murnaghan | Finite deformations of an elastic solid[END_REF] is tackled. The use of this lawi sd iscussed, in particular the wayi t relates to polynomial stress-strain relationships [START_REF] Ostrovsky | Dynamic nonlinear elasticity in geomaterials[END_REF][START_REF] Payan | Determination of third order elastic constants in ac omplexs olid applying coda wave interferometry[END_REF]. The numerical method is presented in section 3, including the construction of aR oe matrix. In section 3.4, the numerical method is validated with the analytical solution of the Riemann problem of nonlinear elastodynamics [START_REF] Berjamin | Analytical solution to 1D nonlinear elastodynamics with general constitutive laws[END_REF]. Section 4introduces the modifications of the equations to account for the slowd ynamics, as well as the modifications of the numerical method. Then, the propagation of a sinusoidal wave in the material is addressed (section 4.4). The results are in qualitative agreement with experimental observations.

Nonlinear elastodynamics

Governing equations

Lagrangianhyperelasticity. Let us consider an homogeneous continuum in which no heat transfer occurs. Furthermore, self-gravitation is neglected. Aparticle initially located at some position x 0 of the reference configuration movest oap osition x t of the current configuration. The deformation gradient is as econd-order tensor defined by (see e.g. [START_REF] Ogden | Non-Linear Elastic Deformations[END_REF][START_REF] Holzapfel | Nonlinear Solid Mechanics: AContinuum Approach for Engineering[END_REF][START_REF] Drumheller | Introduction to wave propagation in nonlinear fluids and solids[END_REF][START_REF] Norris | Finite-amplitude wavesinsolids[END_REF])

F = grad x t = I + grad u, (1)
where u = x t -x 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 )a nd aC artesian coordinate system (O, x, y, z). In this case, the matrix of the coordinates of I is the identity matrix.

Here, the Lagrangian representation of motion is used. The material derivative of the deformation gradient satisfies

.

F = grad v, (2)
where v(x 0 ,t)i st he 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. The motion is also drivenbythe conservation of momentum

ρ 0 . v = div det(F ) σ • F -+ f v , ( 3
)
where div denotes the divergence with respect to the material coordinates. The tensor σ = σ is the Cauchystress tensor,a nd f v is an external volume force applied to the material.

In hyperelasticity,the only variables of state are the specifice ntropy η and as train tensor.M oreover, the dissipation in the material is zero, i.e. the thermodynamic process is isentropic. Here, the deformation of the material is represented by the Green-Lagrange strain tensor,but other choices are possible. The Green-Lagrange strain tensor is defined as E = 1 2 (C -I) where C = F • F is the right Cauchy-Green strain tensor,o re quivalently as af unction of the displacement gradient tensor,

E = 1 2 grad u + grad u + grad u • grad u .(4)
The internal energy per unit of reference volume is W (E), where W is the strain energy density function. Under these assumptions, the expression of the Cauchystress tensor is

σ = 1 det(F ) F • ∂W ∂E • F . ( 5
)
Thus, the first Piola-Kirchhoff tensor in (3) satisfies

det(F ) σ • F -= F • ∂W ∂E . ( 6
)
In the isotropic case, the strain energy is afunction of the invariants of E defined by

E I = tr(E), E II = 1 2 tr(E) 2 -tr(E 2 ) , E III = det(E). (7)
The derivative of the strain energy with respect to the strain tensor writes

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

α 0 = ∂W ∂E I + E I ∂W ∂E II + E II ∂W ∂E III , α 1 = - ∂W ∂E II -E I ∂W ∂E III , α 2 = ∂W ∂E III , ( 9
)
are functions of the invariants (E I ,E II ,E III ).

Longitudinal plane waves. We makethe assumption that the displacement field u has no component along e 2 and e 3 .M oreover, its component u along e 1 is independent on y and z.T herefore, the displacement gradient writes grad u = ε (e 1 ⊗ e 1 ), where ε = ∂ x u>-1i st he axial component of the displacement gradient. The invariants [START_REF] Meurer | Wave propagation in nonlinear and hysteretic media-a numerical study[END_REF] of the Green-Lagrange strain tensor are E I = ε + 1 2 ε 2 and E II = 0 = E III .T hus, the strain energy W is now af unction of ε only.I nt he longitudinal case, the 11coordinate of the first Piola-Kirchhoff tensor (6) is equal to the coordinate σ of the Cauchystress tensor.Equations ( 8)-( 9) and the expression of the invariants provide the stress-strain relationship

σ = (1 + ε) α 0 + E I α 1 + E I 2 α 2 = ∂E I ∂ε ∂W ∂E I = W (ε), (10)
where the denotes the total derivative with respect to ε.

Constitutive laws. The Murnaghan model of hyperelasticity [START_REF] Murnaghan | Finite deformations of an elastic solid[END_REF] is widely used in the communities of geophysics and nondestructive testing [START_REF] Mccall | Theoretical study of nonlinear elastic wave propagation[END_REF][START_REF] Johnson | Nonlinear elasticity and stress-induced anisotropyi nr ock[END_REF][START_REF] Ostrovsky | Dynamic nonlinear elasticity in geomaterials[END_REF][START_REF] Payan | Determination of third order elastic constants in ac omplexs olid applying coda wave interferometry[END_REF]. Its strain energy density function is

W = λ + 2µ 2 E I 2 -2µE II + l + 2m 3 E I 3 -2m E I E II + n E III , (11)
where (λ, µ)are the Lamé parameters and (l, m, n)are the Murnaghan coefficients. Sometimes, Landau'sl aw with parameters (A, B, C)isused instead, and the relationship with Murnaghan'sl aw [START_REF] Lyakhovsky | Distributed damage, faulting, and friction[END_REF] is specified in [START_REF] Norris | Finite-amplitude wavesinsolids[END_REF]. With the above expression of the strain energy [START_REF] Lyakhovsky | Distributed damage, faulting, and friction[END_REF],the axial component of the Cauchystress [START_REF] Favrie | Fast and slowdynamics in anonlinear elastic bar excited by longitudinal vibrations[END_REF] is

σ = M 0 ε 1 + 3 2 + ϑ ε + 1 2 + 2ϑ ε 2 + 5ϑ 4 ε 3 + ϑ 4 ε 4 , ( 12
)
where M 0 = λ + 2µ>0i st he elastic modulus, and ϑ = (l+2m)/(λ+2µ). If the Murnaghan coefficients equal zero (ϑ = 0),o nly geometric nonlinearities remain, and the Saint Ve nant-Kirchhoff model of hyperelasticity is recovered. The later reduces to the classical Hooke'slaw in the case of infinitesimal strain ε 0. When geometric nonlinearities are neglected (e.g., when the Murnaghan coefficients are very large), the Green-Lagrange strain tensor is linearized with respect to the components of grad u,sothat the strain tensor ( 4) is replaced by the infinitesimal strain tensor: E ε.D oing this, the first invariant of the strain tensor is E I ε in the longitudinal case. Moreover, the first Piola-Kirchhoff stress is linearized with respect to the components of grad u as well, i.e. F • ∂W/∂E ∂W/∂ε in ( 6).D oing this, the longitudinal constitutive law( 10)r educes to σ = α 0 + α 1 ε + α 2 ε 2 ,and ( 12) becomes

σ = M 0 ε 1 + ϑε . ( 13
)
This constitutive lawc orresponds to aq uadratic polyno-

mial σ = M 0 ε (1 -βε),with β = -ϑ.
Acubic polynomial constitutive law [START_REF] Meurer | Wave propagation in nonlinear and hysteretic media-a numerical study[END_REF][START_REF] Li | Pseudo-spectral simulation of 1D nonlinear propagation in heterogeneous elastic media[END_REF][START_REF] Favrie | Fast and slowdynamics in anonlinear elastic bar excited by longitudinal vibrations[END_REF]]

σ = M 0 ε 1 -βε -δε 2 (14)
is also widely used to describe nonlinear elasticity in solids. If the geometric nonlinearities are negligible [START_REF] Berjamin | Nonlinear wavesi ns olids with slowd ynamics: an internalvariable model[END_REF], then the parameters β = -ϑ and δ = 0correspond to Murnaghan'slaw.Ifthe geometric nonlinearities are taken into account [START_REF] Lyakhovsky | Nonlinear elastic behaviour of damaged rocks[END_REF],then the choice -β = 3 2 + ϑ and -δ = 1 2 + 2ϑ makes the twomodels coincide up to the third order.Contrary to [START_REF] Lyakhovsky | Nonlinear elastic behaviour of damaged rocks[END_REF],orders 4and 5are not taken into account in [START_REF] Maugin | Thermodynamics with internal variables[END_REF],which shows that both models are not equivalent.

The major difference between the stress-strain relationships ( 12) and ( 14) is the number of independent parameters. In ( 14),the cubic term can be set independently of the quadratic term. In [START_REF] Lyakhovsky | Nonlinear elastic behaviour of damaged rocks[END_REF],the cubic term is not independent on the quadratic term: if ( 12) and ( 14) are assimilated, then δ = 5 2 + 2β.However,e xperimental evidence shows that δ is larger than β by several orders of magnitude, so that aconstitutive lawofthe Murnaghan type is not sufficient to represent accurately elastic nonlinearity in geomaterials [START_REF] Renaud | Revealing highly complexelastic nonlinear (anelastic)behavior of Earth materials applying anew probe: Dynamic acoustoelastic testing[END_REF]. From nowon, the cubic polynomial law( 14)is used.

Hyperbolic system of conservation laws

The equations of motion ( 2)-( 3) write as an on-homogeneous system of conservation laws with respect to the variables q = (ε, v) ,where v is the particle velocity component along e 1 ,

∂ t q + ∂ x f (q) = s. ( 15
)
The flux function is f (q) = -(v, σ(ε)/ρ 0 ) ,and the forcing is s = (0,f v /ρ 0 ) ,where f v is the volume force component along e 1 .The Jacobian matrix of the flux is

f (q) = - 01 σ ( ε ) /ρ 0 0 , ( 16
)
where σ = W is the derivative of σ with respect to ε. The eigenvalues of f (q)are {-c(q),c(q)},where c(q) = σ (ε)/ρ 0 is the speed of sound. The system of conservation laws is strictly hyperbolic if σ (ε) > 0, i.e. overadomain where the strain energy W is astrictly convex function of ε.Some properties of the system (15) are listed belowwithout proof. Interested readers are referred to standard textbooks for more details about hyperbolic systems [START_REF] Godlewski | Numerical Approximation of Hyperbolic Systems of Conservation Laws[END_REF][START_REF] Leveque | Finite-Volume Methods for Hyperbolic Problems[END_REF].

In the case of the cubic law( 14), strict hyperbolicity is ensured if [START_REF] Berjamin | Analytical solution to 1D nonlinear elastodynamics with general constitutive laws[END_REF] ε ∈ 1

β -β 2 + 3δ , 1 β + β 2 + 3δ . ( 17
)
If β = δ = 0, then the polynomial law( 14)a mounts to Hooke'sl aw,a nd the characteristic fields are linearly degenerate. If β = 0a nd δ = 0, then the polynomial law is quadratic, and the stress is either as trictly convex or a strictly concave function of the strain. The characteristic fields are genuinely nonlinear.I fδ =0, then the polynomial law [START_REF] Maugin | Thermodynamics with internal variables[END_REF]isneither convex nor concave.Indeed, an inflexion point is located at ε 0 = -β/3δ.The characteristic fields are neither genuinely nonlinear nor linearly degenerate (i.e., theyare nongenuinely nonlinear). In the case of Murnaghan'slaw [START_REF] Lyakhovsky | Nonlinear elastic behaviour of damaged rocks[END_REF],asimilar analysis is carried out in the appendix.

Finite volumes with flux limiters

Conservative scheme

In the examples presented later on, the physical domain is unbounded. We consider afi nite numerical domain [x 0 ,x N ].I ti sd iscretized using ar egular grid in space with step Δx = (x N -x 0 )/N .A lso, av ariable time step Δt = t n+1 -t n is introduced. Therefore, q(x i ,t n )d enotes the solution to (15) at the abscissa x i = x 0 + i Δx and at the time t n .The volume force f v is assumed to be apoint load f v = δ s (x)ϕ(t), where δ s (x) = δ(x -x s )isthe Dirac delta located at the abscissa x = x s ,and ϕ(t)isthe source signal.

The non-homogeneous system of conservation laws (15) is integrated explicitly:

q n+1 i = q n i - Δt Δx (f n i+1/2 -f n i-1/2 ) + Δt Δx s n i , ( 18
)
where q n i q(x i ,t n )approximates the solution at the grid nodes, and s n i = s(x i ,t n ). The numerical flux f n i+1/2 of the flux-limiter method is specified later on. This scheme is stable under the classical CFL condition

κ = Δt Δx c n max 1, ( 19
)
where κ is the Courant number,and c n max is the maximum sound speed that is encountered at time t n .Ifσis convex or concave,then the maximum sound speed at time t n is

c n max = max 0 i N c(q n i ). ( 20
)
If the constitutive lawisneither convex nor concave,sound speeds larger than (20) may be reached (see e.g. section 16.1 in [START_REF] Leveque | Finite-Volume Methods for Hyperbolic Problems[END_REF]). The local maximum sound speed is obtained by maximizing c(q)for q between q n i and q n i+1 .Then, one obtains

c n max = max 0 i<N 1 ρ 0 max ε∈D i σ (ε) , ( 21
)
where D i is the interval with bounds ε n i and ε n i+1 .F inally,the method has av ariable time step satisfying Δt = κ Δx/c n max . Since the flux function f is nonlinear with respect to q, an initial-value problem (orC auchyp roblem)o ft he homogeneous system defined by the data q(x, 0) at the time t = 0can have several weak solutions (i.e. solutions of an integral form of ( 15)). It is not straightforward for an umerical method to convergetowards the correct weak solution. In particular,the numerical fluxes f n i+1/2 in (44) must be computed carefully.G iven that the eigenvalues ±c(q) of ( 16) have constant sign, no transsonic rarefaction can occur.Moreover, no slow-moving shock can occur either, due to the monotonicity of σ overt he domain of hyperbolicity [START_REF] Leveque | Finite-Volume Methods for Hyperbolic Problems[END_REF].T herefore, Roe linearization can be used to construct an accurate numerical scheme.

Roe linearization

The nonlinear flux f (q)i sa pproximated locally by the linear flux A i+1/2 q.ARoe matrix A i+1/2 approximates the Jacobian f (q n i+1/2 )a tt he midpoint of [x i ,x i+1 ]a nd the time t n ,insuch away that

1. A i+1/2 is diagonalizable with real eigenvalues; 2. A i+1/2 q n i+1 -q n i = f (q n i+1 ) -f (q n i ); 3. lim q n i+1 →q n i A i+1/2 = f (q n i ).
According to section 15.3.2 of [START_REF] Leveque | Finite-Volume Methods for Hyperbolic Problems[END_REF], such amatrix may be expressed by the formula

A i+1/2 = 1 0 f q n i + ξ (q n i+1 -q n i ) dξ, (22)
which ensures that the properties 2. and 3. are satisfied.

In the case of the system [START_REF] Maugin | The sagaofinternal variables of state in continuum thermo-mechanics (1893-2013[END_REF],w hich corresponds to exercise 15.1.(a) p. 349 of [START_REF] Leveque | Finite-Volume Methods for Hyperbolic Problems[END_REF], the following matrix is obtained:

A i+1/2 = 0 -1 -a i+1/2 0 , ( 23
)
where

a i+1/2 = W (ε n i+1 ) -W (ε n i ) ρ 0 ε n i+1 -ε n i . ( 24
)
To avoid divisions by zero when ε n i = ε n i+1 ,one computes

a i+1/2 = W (ε n i ) ρ 0 , ( 25
)
in this particular case. The eigenvalues of the matrix ( 23) are {-s i+1/2 ,s i+1/2 },w here s i+1/2 = a i+1/2 .S ince the strain energy function W is convex overthe hyperbolicity domain, the coefficient a i+1/2 is positive.T herefore, the eigenvalue s i+1/2 is real, and the property 1. is satisfied.

The matrix (47) is aRoe matrix.

Forlater use, we introduce the decomposition of q n i+1q n i in the basis of right eigenvectors of the Roe matrix

p 1 i+1/2 = 1,s i+1/2 , p 2 i+1/2 = 1,-s i+1/2 . ( 26
)
The jump of the solution can be expanded as

q n i+1 -q n i = 2 k=1 α k i+1/2 p k i+1/2 , = 2 k=1 W k i+1/2 , ( 27
)
with the coefficients

α 1 i+1/2 = 1 2 ε n i+1 -ε n i + v n i+1 -v n i s i+1/2 , α 2 i+1/2 = 1 2 ε n i+1 -ε n i - v n i+1 -v n i s i+1/2 .
(28)

Flux limiter

We describe nowt he flux-limiter scheme. The numerical flux in [START_REF] Ve Lasco-Segura | Afinite volume approach for the simulation of nonlinear dissipative acoustic wave propagation[END_REF] takes the form (section 15.4 in [START_REF] Leveque | Finite-Volume Methods for Hyperbolic Problems[END_REF])

f n i+1/2 = f L i+1/2 + f H i+1/2 , ( 29
)
where f L i+1/2 is the flux of Godunov'sm ethod or one of its approximated versions, and f H i+1/2 is ah igher-order correction deduced from the Lax-Wendroff method. Here, the Godunovflux is approximated linearly by the classical Roe flux, and the Roe matrix (47) is used. Thus, where W k i+1/2 is defined in [START_REF] Mccall | Theoretical study of nonlinear elastic wave propagation[END_REF].The higher-order correction of the flux-limiter method writes

f L i+1/2 = 1 2 f (q n i ) + f (q n i+1 ) - 1 2 s i+1/2 W 1 i+1/2 + W 2 i+1/2 , (30)
f H i+1/2 = 1 2 s i+1/2 1 - Δt Δx s i+1/2 • φ(θ 1 i+1/2 )W 1 i+1/2 + φ(θ 2 i+1/2 )W 2 i+1/2 , ( 31
)
where

θ 1 i+1/2 = W 1 i+3/2 • W 1 i+1/2 W 1 i+1/2 • W 1 i+1/2 , θ 2 i+1/2 = W 2 i-1/2 • W 2 i+1/2 W 2 i+1/2 • W 2 i+1/2 , ( 32
)
and φ is alimiter function. Here, the minmod limiter

φ(θ) = max{0, min{1,θ}} (33)
is used. The weights φ(θ k i+1/2 )a re designed to avoid spurious oscillations in the numerical solution. If the weights φ(θ k i+1/2 )equal one in [START_REF] Berjamin | RiemannElasto1D Toolbox[END_REF],then the Lax-Wendroff method based on the Roe matrix (47) is recovered. In the case of nonconvex flux functions, the minmod limiter (33) has shown better convergence properties than other limiter functions. Similar observations are reported in [START_REF] Kurganov | Adaptive semidiscrete central-upwind schemes for nonconvex hyperbolic conservation laws[END_REF][START_REF] Vo Ss | Exact Riemann solution for the Euler equations with nonconvex and nonsmooth equation of state[END_REF].

To carry out one iteration in time at some grid node i, the numerical values of q at the grid nodes i -2,...,i+2 are required. Therefore, two" ghost cells" must be added on the left and on the right of the numerical domain. Here, az ero-order extrapolation of the numerical solution can be used to achieve outflowboundary conditions. Thus, we simply set

q n -2 = q n -1 = q n 0 , q n N+2 = q n N+1 = q n N , (34)
at each time step [START_REF] Leveque | Finite-Volume Methods for Hyperbolic Problems[END_REF].

Numerical test case

The orders of magnitude of the elastic parameters for Berea sandstone ρ 0 , M 0 , β, δ in table Ihavebeen taken in Table Iand Figure 5a of [START_REF] Renaud | Revealing highly complexelastic nonlinear (anelastic)behavior of Earth materials applying anew probe: Dynamic acoustoelastic testing[END_REF]. From [START_REF] Leveque | Finite-Volume Methods for Hyperbolic Problems[END_REF],itfollows that the model ( 15) is hyperbolic if |ε| < 5.7 • 10 -5 .The numerical domain is [x 0 ,x N ] = [-0.5,0.5] m. To avoid instability, the Courant number ( 19) is set to κ = 0.95. This test is carried out to validate the ability of the numerical scheme to represent accurately nonsmooth solutions. Here, no volume force is applied in the material. We consider piecewise constant initial data with as ingle discontinuity at the abscissa x = 0, i.e. aRiemann problem. Here, the initial data is ε = 2 • 10 -5 , v = 0ifx<0, and ε = -2.6 • 10 -5 , v = 0.097 m/s if x>0. The analytical solution to this initial-value problem is displayed in Figure 1a. Details about its computation can be found in [START_REF] Berjamin | Analytical solution to 1D nonlinear elastodynamics with general constitutive laws[END_REF]. Moreover, an interactive application and aM atlab toolbox [START_REF] Berjamin | RiemannElasto1D Toolbox[END_REF] can be found at http://gchiavassa.perso.centralemarseille.fr/RiemannElasto/.

The numerical solution (44) is computed up to t = 0.19 ms on agrid with N points. Figures 1b and1c compare the analytical solution with the numerical solution. On Figure 1b, the left-going discontinuity is represented (shock wave). On Figure 1c, the right-going compound wave with ac ontinuous part followed by ad iscontinuous part is represented (rarefaction shock wave). Both waves are well-located and accurately computed by the numerical scheme. In particular,n os purious oscillations arise, and no convergence failure is noticed as N increases from 1000 to 2000.

Coupling with slowd ynamics

As specified in the introduction, nonlinear elastodynamics presented in Section 2i sn ot sufficient to describe the softening of the material. In [START_REF] Berjamin | Nonlinear wavesi ns olids with slowd ynamics: an internalvariable model[END_REF], as imple approach was proposed to account for such phenomena. The model was presented in 3D, and analytical computations were carried out in the longitudinal case. Here, we showh ow the previous numerical method is adapted to this model of slow dynamics.

Modified equations of motion

Internal-variable model. An internal-variable g in [0, 1[ is introduced, to represent the softening of the material. The internal energy per unit of reference volume ρ 0 e = W (ε)becomes [START_REF] Berjamin | Nonlinear wavesi ns olids with slowd ynamics: an internalvariable model[END_REF]

ρ 0 e = (1 -g) W (ε) +Φ(g), ( 35
)
where e is the specificinternal energy.In( 35), Φ(g)r epresents astorage energy with expression

Φ(g) = - 1 2 γ ln 1 -g 2 , ( 36
)
where γ>0i sa ne nergy per unit volume, buto ther choices are suitable. As discussed in [START_REF] Berjamin | Nonlinear wavesi ns olids with slowd ynamics: an internalvariable model[END_REF], the storage energy ( 36) is chosen such that equilibrium points of the model are unique, g = 0i sa ne quilibrium point, and g is bounded by 1. With such achoice, one observes that the internal energy per unit of reference volume (35) is equal to the strain energy W when g is equal to zero. The second principle of thermodynamics yields the axial component of the Cauchystress [START_REF] Berjamin | Nonlinear wavesi ns olids with slowd ynamics: an internalvariable model[END_REF] [START_REF] Favrie | Fast and slowdynamics in anonlinear elastic bar excited by longitudinal vibrations[END_REF]. The simplest thermodynamically admissible choice of evolution equation is [START_REF] Berjamin | Nonlinear wavesi ns olids with slowd ynamics: an internalvariable model[END_REF]

σ = 1 -g W (ε)( 37) instead of W (ε)
τ 1 . g = W (ε) -Φ (g), ( 38
)
where τ 1 = τ • 1J/m 3 and τ>0isatime constant. With the laws (37) and (38),the dissipation per unit of reference volume is [START_REF] Berjamin | Nonlinear wavesi ns olids with slowd ynamics: an internalvariable model[END_REF]

ρ 0 T . η = -ρ 0 ∂e ∂g . g = τ 1 . g 2 0, ( 39
)
where T>0isthe absolute temperature, and η is the specificentropy.

Conservation laws with relaxation. The variable g is added to the vector of unknowns q,sothat q = (ε, v, g) . The system of conservation laws (15) becomes

∂ t q + ∂ x f (q) = r(q) + s, (40)
where

f (q) = -v, (1 -g) W (ε)/ρ 0 , 0 , r(q) = 0, 0, W (ε) -Φ (g) /τ 1 , s = (0,f v /ρ 0 , 0) . ( 41
)
The Jacobian matrix f (q)o ffhas the eigenvalues {-c(q),c(q),0},where

c(q) = (1 -g) W (ε) ρ 0 (42)
is the speed of sound. The speed of sound ( 42) is real and nonzero -ino ther words, the system (40) is strictly hyperbolic -provided that g<1and W (ε) > 0. This second condition wasalready required in the elastic case [START_REF] Maugin | The sagaofinternal variables of state in continuum thermo-mechanics (1893-2013[END_REF]. Nowlet us examine the spectrum of the relaxation function in (40).The Jacobian matrix r (q)ofrhas the eigenvalues {0, 0, -Φ (g)/τ 1 }.T he expression (36) ensures that r (q)i sn egative semi-definite. Its spectral radius is

Υ = γ τ 1 1 + g 2 (1 -g 2 ) 2 , (43)
which involves the relaxation time τ 1 /γ characteristic of the slowdynamics [START_REF] Berjamin | Nonlinear wavesi ns olids with slowd ynamics: an internalvariable model[END_REF].

Modified numerical method

Conservative scheme. Applying the same explicit discretization as [START_REF] Ve Lasco-Segura | Afinite volume approach for the simulation of nonlinear dissipative acoustic wave propagation[END_REF] to (40) yields

q n+1 i =q n i - Δt Δx (f n i+1/2 -f n i-1/2 ) +Δtr(q n i )+ Δt Δx s n i , (44)
where the numerical flux f n i+1/2 is specified later on. According to the time-marching formula (44),explicit singlestep time integration is used.

Numerical stability imposes abound of the form

Δt min Δx c n max , 2 Υ n max , ( 45
)
where c n max is the maximum sound velocity (42) that is encountered at time t n ,a nd Υ n max is the maximum spectral radius of the relaxation function (43).Since the relaxation time of the slowdynamics is much larger than the period of exciting signals, the condition (45) reduces to the classical CFL condition [START_REF] Murnaghan | Finite deformations of an elastic solid[END_REF].Ifthe constitutive law ε → W (ε) is convex or concave,t hen the maximum sound speed at time t n is givenby [START_REF] Ostrovsky | Dynamic nonlinear elasticity in geomaterials[END_REF], with the sound speed (42).Ifthe constitutive lawi sn either convex nor concave,t hen (21) becomes

c n max = max 0 i<N 1 -min{g n i ,g n i+1 } ρ 0 max ε∈D i W (ε) , ( 46
)
where D i is the interval with bounds ε n i and ε n i+1 .

Roe linearization. In the case of the system (40),the following matrix is obtained from [START_REF] Berjamin | Analytical solution to 1D nonlinear elastodynamics with general constitutive laws[END_REF]:

A i+1/2 =   0 -10 -a i + 1 / 2 0 b i + 1 / 2 00 0   , (47)
where

a i+1/2 = σ n i+1 -σ n i ρ 0 ε n i+1 -ε n i + g n i+1 -g n i ε n i+1 -ε n i b i+1/2 , b i+1/2 = W (ε n i+1 ) -W (ε n i ) ρ 0 ε n i+1 -ε n i , (48)
and the stresses (37) are

σ n i = (1 -g n i ) W (ε n i ).
To avoid divisions by zero when ε n i = ε n i+1 ,one computes

a i+1/2 = 1 - g n i + g n i+1 2 W (ε n i ) ρ 0 , b i+1/2 = W (ε n i ) ρ 0 , (49)
in this particular case. The eigenvalues of the matrix (47) with the coefficients (48)-(49)a re {-s i+1/2 ,s i+1/2 ,0},w here the expression s i+1/2 = a i+1/2 is unchanged. We rewrite the coefficient a i+1/2 from (48) as

a i+1/2 = W (ε n i+1 ) -W (ε n i ) ρ 0 ε n i+1 -ε n i -g n i W (ε n i+1 ) -W (ε n i ) -W (ε n i ) ε n i+1 -ε n i ρ 0 ε n i+1 -ε n i 2 -g n i+1 W (ε n i ) -W (ε n i+1 ) -W (ε n i+1 ) ε n i -ε n i+1 ρ 0 ε n i+1 -ε n i 2 . ( 50
)
One can note that the first term in (50) corresponds to the elastic case [START_REF] Holzapfel | Nonlinear Solid Mechanics: AContinuum Approach for Engineering[END_REF].M oreover, the next terms vanish when g n i = 0 = g n i+1 ,i .e. when no softening occurs. Since the strain energy function W is convex and g<1o vert he hyperbolicity domain, the coefficient a i+1/2 in (50) is positive.F inally,t he eigenvalue s i+1/2 is real, and the property 1. is satisfied. The matrix (47) is aR oe matrix for (40). Now, the jump q n i+1 -q n i is decomposed in the basis of right eigenvectors of the Roe matrix (47)

p 1 i+1/2 = 1,s i+1/2 ,0 , p 2 i+1/2 = 1,-s i+1/2 ,0 , p 3 i+1/2 = 1,0, a i+1/2 b i+1/2 . ( 51
)
The jump of the solution can be expanded as as um [START_REF] Mccall | Theoretical study of nonlinear elastic wave propagation[END_REF] of three W k i+1/2 = α k i+1/2 p k i+1/2 ,where the coefficients are

α 1 i+1/2 = 1 2 v n i+1 -v n i s i+1/2 + σ n i+1 -σ n i ρ 0 s i+1/2 2 , α 2 i+1/2 = 1 2 - v n i+1 -v n i s i+1/2 + σ n i+1 -σ n i ρ 0 s i+1/2 2 , α 3 i+1/2 = ε n i+1 -ε n i - σ n i+1 -σ n i ρ 0 s i+1/2 2 . ( 52
)
Since the eigenvalue corresponding to k = 3inthe decomposition of the jump is zero, W 3 i+1/2 does not appear in the numerical flux [START_REF] Kurganov | Adaptive semidiscrete central-upwind schemes for nonconvex hyperbolic conservation laws[END_REF] of the flux-limiter method. Therefore, the formulas ( 29) to [START_REF] Shu | High order weighted essentially nonoscillatory schemes for convection dominated problems[END_REF] can be applied without modification with the flux function (41),t he eigenvalue s i+1/2 deduced from (48)-( 49), and the decomposition of the jump (51)-(52).

Pulse propagation

Figure 3of [START_REF] Renaud | Revealing highly complexelastic nonlinear (anelastic)behavior of Earth materials applying anew probe: Dynamic acoustoelastic testing[END_REF] provides an order of magnitude of the characteristic time of the slowdynamics, and an order of magnitude of the average softening ΔM/M for agiven strain amplitude V .According to [START_REF] Berjamin | Nonlinear wavesi ns olids with slowd ynamics: an internalvariable model[END_REF], the characteristic time of the slowdynamics is τ 1 /γ,and the average softening satisfies ΔM/M -1 4 ( M 0 /γ + 6δ) V 2 .Ifthis scaling rule is combined with the values in [START_REF] Renaud | Revealing highly complexelastic nonlinear (anelastic)behavior of Earth materials applying anew probe: Dynamic acoustoelastic testing[END_REF], then an egative value of γ is obtained, which is not reliable. Thus, the parameters γ and τ in table Iare not chosen to reach quantitative agreement with [START_REF] Renaud | Revealing highly complexelastic nonlinear (anelastic)behavior of Earth materials applying anew probe: Dynamic acoustoelastic testing[END_REF], butrather to obtain comparable orders of magnitude for the characteristics of softening.

The numerical domain and the Courant number are the same as in section 3.4, buth ere, the softening is taken into account. The initial data is zero, and ap oint load f v = δ s (x)ϕ(t)l ocated at the abscissa x s = 0g enerates as mooth pulse with angular frequency ω c = 2πf c ,w ith f c = 6kHz. The source signal in (40) satisfies

ϕ(t) = 2M 0 V 4 m=1 a m sin 2 m-1 ω c t , (53)
where V = 6.64 • 10 -6 , a 1 = 1, a 2 = -21/32, a 3 = 63/768, and a 4 = -1/512. It is turned on from t = 0t o t=1 /f c ,w hich corresponds to one fundamental period. Chosen for the smoothness of its time-evolution (53),the point source generates left-going and right-going waves with strain amplitude 1.507 V ≈ 10 -5 . The reference solution is an oversampled numerical solution (44) computed at t = 0.19 ms on ag rid with N ref = 2 15 points (Figure 2a). The pulse injected at x s = 0 has propagated towards both increasing and decreasing x. By symmetry,only the right-going part is displayed here. No shock wave is observed: the waveform is slightly distorted butk eeps smooth. Figure 2b illustrates the convergence of the numerical method. Foragiven coarse spatial discretization where N N ref ,the numerical solution is computed up to t ≈ 0.19 ms, and is compared to the reference numerical solution at the same final time. The evolution of the L 2 global error between both strain waveforms is represented in Figure 2b with respect to Δx = 1/N .One can observethat the order of convergence is between one and two.

Dynamic acousto-elasticity

The setup is the same as in the previous section, buthere, the point load generates asinusoidal strain with amplitude V ≈ 10 -6 and angular frequency ω c = 2πf c ,w ith f c = 5kHz:

ϕ(t) = 2M 0 V sin ω c t . ( 54
)
The source (54) is turned on from t = 0t ot = 40 ms, which corresponds to 200 periods of signal. Ar eceiver records the numerical solution at the abscissa x r = 0.3m. The numerical solution (44) is computed up to t = 80 ms on ag rid with N = 80 points, which corresponds until the source is stopped. Simultaneously,the strain signals are smooth sinusoids, so that N = 80 is sufficient. After the source is stopped, the elastic modulus recovers slowly its initial value, while the strain is equal to zero.

Figure 3b focuses on the steady-state solution. Here, the last 80 numerical values before t = 40 ms are used, which corresponds to twop eriods of signal at the frequency f c . When ΔM/M from (55) is represented with respect to the strain recorded at the position of the receiver x r ,a hysteresis curvei so btained. The orders of magnitude of these phenomena -duration of the transients, magnitude of the softening, size and shape of the hysteresis curvesare very similar to those reported in [START_REF] Renaud | Revealing highly complexelastic nonlinear (anelastic)behavior of Earth materials applying anew probe: Dynamic acoustoelastic testing[END_REF] for Berea sandstone. In particular,o ne can note that the average softening ΔM/M is proportional to V 2 ,a sp redicted in [START_REF] Berjamin | Nonlinear wavesi ns olids with slowd ynamics: an internalvariable model[END_REF]. However, only qualitative agreement with Figure 5a of [START_REF] Renaud | Revealing highly complexelastic nonlinear (anelastic)behavior of Earth materials applying anew probe: Dynamic acoustoelastic testing[END_REF] is obtained.

Conclusion

Afi nite-volume method with flux-limiters for nonlinear longitudinal elastodynamics is implemented, with various polynomial constitutive laws. Based on Roe linearization, this numerical method is well-suited to the present system of partial differential equations, and has been validated with the Riemann problem of nonlinear elastodynamics. Nonsmooth solutions such as shock wavesa re well-captured, even in the case of nonconvex stress-strain relationships. The method is adapted to acase with an additional scalar evolution equation, which is deduced from the internal-variable model [START_REF] Berjamin | Nonlinear wavesi ns olids with slowd ynamics: an internalvariable model[END_REF] of slowdynamics. The numerical model reproduces qualitatively experimental observations related to dynamic acoustoelasticity.

We mention here afew improvements to be introduced. Boundary conditions such as free edges and oscillating walls need to be implemented so as to reproduce real configurations in amore realistic way(section 7of [START_REF] Leveque | Finite-Volume Methods for Hyperbolic Problems[END_REF]). Due to resonance, viscoelastic attenuation needs then to be incorporated in the model, which has been carried out simi-larly to [START_REF] Favrie | Fast and slowdynamics in anonlinear elastic bar excited by longitudinal vibrations[END_REF] in recent works [START_REF] Berjamin | Modeling longitudinal wave propagation in nonlinear viscoelastic solids with softening[END_REF]. Currently,higher-order shock-capturing methods [START_REF] Shu | High order weighted essentially nonoscillatory schemes for convection dominated problems[END_REF] are developed, as well as similar methods in multiple space dimensions.

  Received18October 2017, accepted 14 May 2018.

Figure 1 .

 1 Figure 1. (a) Analytical solution to aR iemann problem. Comparison with the numerical solution (b) on the left-going strain wave,(c) and on the right-going strain wave.

Figure 2 .

 2 Figure 2. Assessment of the numerical method'sconvergence on the propagation of asmooth pulse. (a) Reference solution at t = 0.19 ms: strain, particle velocity,and softening variable. (b) Error measurement.

Figure 3 .Figure 4 .

 34 Figure 3. (Colour online)( a) Strain signal and softening of the material, as recorded by the receiveratthe abscissa x r = 0.3m. (b) Hysteresis curves in steady-state. The arrowi ndicates how time increases along the curve.to 40 points per wavelength at the angular frequency ω c . The computations are performed in C++.Each simulation lasts around 1.5 s, when arecent desktop computer is used (Intel Core i5-4690, 3.5 GHz,16 Go, 2015). Figure3adisplays the strain ε and the variation of the elastic modulusΔM M = ρ 0 c 2 -M 0 M 0 (55)deduced from the sound speed (42),which are recorded at the position x r of the receiver. Aslowdecrease of the elastic modulus combined with fast oscillations is observed

Table I .

 I Physical parameters of Berea sandstone.

	ρ 0 (kg/m 3 ) M 0 (GPa)	βδγ (J/m 3 ) τ (ms)
	2.2 • 10 3	14	50 10 8	1.0	10

Acknowledgments

This work wassupported by the interdisciplinary mission of CNRS (INFINITI). The project leading to this publication has receivedf unding from Excellence Initiative of Aix-Marseille University -A *MIDEX, aF rench "Investissements d'Avenir" programme. It has been carried out in the framework of the LabexM EC. Opération réalisée avec le concours du Programme d'investissements d'avenir du Gouvernement français dont la gestion aé té confiée àl'Andra.

Appendix

In the case of Murnaghan'sl aw [START_REF] Lyakhovsky | Nonlinear elastic behaviour of damaged rocks[END_REF],s trict hyperbolicity of ( 15) is ensured if

Therefore, one must have

where

The constitutive law(12)isconvex if σ (ε) 0, i.e.

3

Hence, Murnaghan'slaw is locally concave at small strains (ε 0) provided that ϑ -3/2, otherwise it is locally convex at small strains. The inflexion point

is represented in Figure 4, as well as the hyperbolicity domains (A2) of Murnaghan'slaw.