M C Barick

Y Gaillard

A Lejeune

F Amiot

F Richard

On the uniqueness of intrinsic viscoelastic properties of materials extracted from nanoindentation using FEMU

Keywords: Viscoelasticity, Nanoindentation, Polymers, Identifiability, Uniqueness, FEMU

Instrumented nanoindentation is widely used to extract the material properties from the measured force-displacement curves. In this work, the uniqueness/non-uniqueness of the intrinsic viscoelastic properties of materials determined by nanoindentation during load-unload tests is investigated. A four-parameter viscoelastic law with constant Poisson's ratio is used to model the mechanical behavior of a polymer material and a 2D-axisymetric Finite Element Model (FEM) is used to simulate the nanoindentation test. Firstly, a nanoindentation experimental triangular load-unload test is performed on a bulk sample of polypropylene (PP) with a Berkovich indenter tip at a depth rate of 1000 nm/min. The values of the four material parameters are estimated by the Finite Element Model Updating (FEMU). The numerical results can accurately fit the experimental data. However, several quasi-solutions are shown to exist. These load-unload data allow to identify only three viscoelastic parameters if the Poisson's ratio is known. Secondly, the effect of nanoindentation depth rate, loading type (triangular, trapezoidal, exponential, sinusoidal) and apex angle is numerically investigated using an identifiability index based on the conditioning of the inverse problem. We show a correlation between the identifiability index and the energy dissipated by the material during the tests. The extraction of all material parameters remains impossible using a single test. Finally, some combinations of several nanoindentation triangular tests and indenter tip angles are also investigated. We show that a dual nanoindentation technique (cube corner and Berkovich tips) with triangular load-unload tests is an interesting combination to reliably extract all the viscoelastic parameters, provided that plasticity is taken into account. This result illustrates the interest of using this numerical identifiability index to design nanoindentation experiments to ensure the robustness of the intrinsic viscoelastic properties extraction.

Introduction

Nanoindentation is a powerful technique for mechanical characterization of polymeric materials at small-scale [START_REF] Cheng | Scaling, dimensional analysis, and indentation measurements[END_REF][START_REF] Vanlandingham | Viscoelastic characterization of polymers using instrumented indentation. I. Quasistatic testing[END_REF][START_REF] Zhang | Characterization of mechanical properties of polymers by nanoindentation tests[END_REF]. During the nanoindentation test, a rigid indenter tip is pushed into the surface of a solid and then removed. The force 𝑃(𝑡) and the displacement ℎ(𝑡) are continuously recorded during loading and unloading segments. The resulting force-displacement curve is associated to the material properties, the geometry of the sample and the geometry of the indenter tip. Several models have been proposed to retrieve the elastic modulus and the hardness from such experimental data [START_REF] Doerner | A method for interpreting the data from depth-sensing indentation instruments[END_REF][START_REF] Oliver | An improved technique for determining hardness and elastic modulus using load and displacement sensing indentation experiments[END_REF]. These methods generally assume that the material features a purely elastic behavior during the unloading part and does not exhibit any loading rate dependence [START_REF] Oliver | An improved technique for determining hardness and elastic modulus using load and displacement sensing indentation experiments[END_REF]. However, most of polymers exhibit a significant viscoelastic behavior [START_REF] Tang | Accurate measurement of tip-sample contact size during nanoindentation of viscoelastic materials[END_REF]. It is therefore important to develop a method to retrieve the intrinsic viscoelastic properties of a material from the temporal data of such nanoindentation test. Two approaches have thus been proposed in the literature instead of the classical Oliver & Pharr method [START_REF] Chen | Identification of finite viscoelasticity and adhesion effects in nanoindentation of a soft polymer by inverse method[END_REF]. The first approach is based on the viscoelastic contact theory [START_REF] Lee | Stress analysis in visco-elastic bodies[END_REF][START_REF] Radok | Visco-elastic stress analysis[END_REF][START_REF] Lee | The contact problem for viscoelastic bodies[END_REF][START_REF] Graham | The contact problem in the linear theory of viscoelasticity[END_REF]. The dimensional analysis or the Laplace transform method are used to extend the elastic solution to viscoelastic phenomena. The models yield closed-form solutions, which are used to analyze the nanoindentation test. The parameters are obtained by fitting the experimental forcedisplacement data [START_REF] Cheng | Flat-punch indentation of viscoelastic material[END_REF][START_REF] Vandamme | Viscoelastic solutions for conical indentation[END_REF][START_REF] Liu | Load-displacement relations for nanoindentation of viscoelastic materials[END_REF][START_REF] Oyen | Analytical techniques for indentation of viscoelastic materials[END_REF]. These parameters, such as compliance constants and retardation times describe a mechanical system behavior in which the material is involved, but these are not intrinsic to the material. Indeed, these models generally include correction factors, which mix intrinsic material properties with geometrical consideration, to tune the contact conditions. A good example is the 𝛽 factor which corrects the Sneddon relation of elasticity for non-axisymmetric indentation [START_REF] Oliver | An improved technique for determining hardness and elastic modulus using load and displacement sensing indentation experiments[END_REF]. The vast majority of reported results make use of 𝛽 = 1.034 [START_REF] King | Elastic analysis of some punch problems for a layered medium[END_REF]. It is however known that this factor strongly depends on the tested material [START_REF] Oliver | Measurement of hardness and elastic modulus by instrumented indentation: Advances in understanding and refinements to methodology[END_REF]. This furthermore makes these models often over-parameterized, thus leading to multiple solutions for the fitting procedure [START_REF] Menčík | Characterization of viscoelastic-plastic properties of solid polymers by instrumented indentation[END_REF]. Converting the properties estimated by this approach should into inputs for structural engineering problems based on the finite element analysis is therefore a difficult and non-obvious task. The second approach is carried out by combining Finite Element Method (FEM) and numerical optimization. In this method, the objective function, which is a norm of difference between the numerical nanoindentation force and/or displacement and experimental data, is minimized using optimization techniques [START_REF] Qasmi | Viscoelastic mechanical properties determined by nanoindentation tests and its numerical modelling of polypropylene modified by He+ particle implantation and e-irradiation[END_REF][START_REF] Guessasma | Viscoelasticity properties of biopolymer composite materials determined using finite element calculation and nanoindentation[END_REF]. The parameters of the model are determined as the minimizer of the objective function. However, the uniqueness of this minimizer is generally not assessed in the literature, but it remains a fundamental question, particularly in instrumented nanoindentation. In fact, in the case of elastoplastic behavior, numerous works have shown that a group of materials with distinct elastoplastic properties may yield almost the same conical indentation 𝑃-ℎ curve [START_REF] Cheng | Can stress-strain relationships be obtained from indentation curves using conical and pyramidal indenters?[END_REF][START_REF] Capehart | Determining constitutive models from conical indentation: Sensitivity analysis[END_REF][START_REF] Alkorta | Absence of one-to-one correspondence between elastoplastic properties and sharp-indentation loadpenetration data[END_REF]. It implies that the material properties cannot be uniquely determined by using a single sharp indenter tip. In order to address this problem in the case of elastoplastic behavior, dual or multiple indentation techniques have been proposed by several authors [START_REF] Le | A computational study on the instrumented sharp indentations with dual indenters[END_REF][START_REF] Heinrich | Determination of material properties using nanoindentation and multiple indenter tips[END_REF][START_REF] Le | Improved reverse analysis for material characterization with dual sharp indenters[END_REF]. However, the existence of mystical materials that give almost similar 𝑃-ℎ curves for different indenter tips with half angles ranging from 60° to 80° has also been shown [START_REF] Chen | On the uniqueness of measuring elastoplastic properties from indentation: the indistinguishable mystical materials[END_REF]. Recently, this problem of nonuniqueness of the elastoplastic properties was investigated by [START_REF] Phadikar | On the uniqueness and sensitivity of indentation testing of isotropic materials[END_REF]. They found that non-uniqueness of the solution is caused by a high sensitivity to the experimental errors. They also demonstrated that dual nanoindentation techniques are reliable when the experimental error is within ±1%. This question is poorly addressed in the presence of viscous phenomena (viscoelastic and/or viscoplastic). [START_REF] Constantinescu | On the identification of elastoviscoplastic constitutive laws from indentation tests[END_REF] highlighted this difficulty in the case of Maxwell and Norton-Hoff behaviors.

In this paper, the uniqueness/non-uniqueness of the viscoelastic properties of materials determined by the Finite Element Model Updating (FEMU) of the nanoindentation test is studied. A four-parameter viscoelastic behavior law with constant Poisson's ratio, which has been implemented in a 2D axisymmetric finite element model of the nanoindentation test, is first described in section 2. Then, the FEMU method used to adjust the force-displacement curves on the experimental data is presented. As an experimental reference, a nanoindentation triangular test at constant depth rate on a PP sample is performed. It is shown in the section 3 that depending on the starting point used in the identification algorithm, multiple solutions for the values of the four viscoelastic parameters (𝐸, 𝑐 1 , 𝜈, 𝜂) have been obtained. The quality of this estimation is discussed using an identifiability analysis. In this aim, an identifiability index (I-index) is used to quantify the richness of the nanoindentation data taken into account in the cost function [START_REF] Richard | Viscoelastic modeling and quantitative experimental characterization of normal and osteoarthritic human articular cartilage using indentation[END_REF]. It is also used to define the more adequate experimental procedure and define what would be the experiment(s) to add to a single triangular test in order to assess a unique solution. The identifiability analysis is first performed considering a single nanoindentation triangular test. The effect of nanoindentation depth rate using triangular, trapezoidal, exponential, sinusoidal loading and indenter tip angle is then investigated. The comparison between the I-index results and the dissipation (loss factor) is performed. Finally, the combinations of several nanoindentation triangular tests and indenter tip angles are carried out. A dual nanoindentation technique is shown to yield a unique solution for the inverse problem. To this end, two numerical nanoindentation triangular tests are carried out at a constant depth rate using equivalent cone apex angles of cube corner (42.28°) and Berkovich (70.3°) indenter tips to obtain pseudo-experimental data. The identifiability results for this technique show that a unique solution, which is robust with respect to the noise, can be obtained. However, a poor agreement is found when this solution is compared to the one obtained experimentally. Supported by residual imprint observed on the surface of the sample, plasticity seems to be the cause of this gap. Indeed, a visco-elastic plastic (VEP) model is finally used to describe the material behavior. This VEP model allows a much better agreement between experiment and simulation even if the combination of its five obtained parameters resulting from the updating process is not unique. Finally, it is shown that knowing the yield strength and Poisson's ratio of the material, the identification of the three other parameters is robust, i.e plasticity knowledge helps to identify visco-elastic parameters. These results are supported by a validation performed with tensile response of the material.

Material and methods

Material and viscoelastic law

In this study, a copolymer PP specimen with dimensions of about 10 × 10 × 0.5 mm 3 produced by the Goodfellow company has been used (PP301440). An isotropic linear viscoelastic law with constant Poisson's ratio 𝜈 is chosen to model the behavior of this material. The Helmholtz free energy 𝜓 [START_REF] Lemaitre | Mechanics of solid materials[END_REF] can be written as:

𝜓 = 1 2𝜌 (𝛆 𝑒 : 𝐂: 𝛆 𝑒 + 𝛂 𝑎𝑛 : 𝐂 𝑎𝑛 : 𝛂 𝑎𝑛 ) (1)
where 𝛆 𝑒 is the elastic strains tensor, 𝛂 𝑎𝑛 is the internal variables tensor representing the anelastic phenomena, 𝜌 is the density, 𝐂(𝐸, 𝜈) and 𝐂 𝑎𝑛 (𝑐 1 , 𝜈) are the elastic and anelastic fourth-order stiffness tensors. 𝐸 and 𝑐 1 are the instantaneous modulus, and the anelastic modulus, respectively; the symbol ":" stands for the tensor inner product. The state laws derive from this energy:

𝛔 = 𝜌 ∂𝜓 ∂𝛆 𝑒 and 𝐗 𝑎𝑛 = 𝜌 ∂𝜓 ∂𝛂 𝑎𝑛 (2)
where 𝛔 is the Cauchy stress tensor and 𝐗 𝑎𝑛 is the anelastic stress tensor.

The dissipation potential 𝛺 is defined as:

𝛺 = 𝐸 2𝜂 (𝛔 -𝐗 𝑎𝑛 ): 𝐒: (𝛔 -𝐗 𝑎𝑛 ) ( 3
)
where 𝜂 is the viscosity coefficient and 𝐒 the elastic compliance (fourth-order) tensor such as 𝐒: 𝐂 = 𝐈 (identity tensor). The derivatives of this potential 𝛺 give the internal variables evolutions:

𝛆̇𝑎 𝑛 = ∂𝛺 ∂𝛔 and 𝛂̇𝑎 𝑛 = - 𝛛𝛺 𝛛𝐗 𝑎𝑛 (4)
where 𝛆 𝑎𝑛 is the anelastic strain which is defined as the difference between the total 𝛆 and elastic 𝛆 𝑒 strains.

𝛆 𝑎𝑛 = 𝛆 -𝛆 𝑒 (5)
This linear viscoelastic behavior law with constant Poisson's ratio is controlled by four material parameters, which define the parameter set 𝛉 = (𝜃 1 , 𝜃 2 , 𝜃 3 , 𝜃 4 ) = (𝐸, 𝑐 1 , 𝜈, 𝜂).

From a rheological point of view, the chosen behavior is, for elastic contribution, a linear spring whose stiffness is given by the Young's modulus 𝐸 and Poisson's ratio 𝜈, and for anelastic contribution (delayed elasticity), a classical Kelvin-Voigt model which consists of a linear viscous damper of viscosity 𝜂 and a linear spring of stiffness 𝑐 1 with internal stress X an associated in parallel. One assumes herein a constant Poisson's ratio, so that 𝐸 𝐺 ⁄ = 𝑐 1 𝐺 1 ⁄ = 𝜂 𝜂 𝐺 ⁄ = 2(1 + 𝑣) and 𝐸 𝐾 ⁄ = 𝑐 1 𝐾 1 ⁄ = 𝜂 𝜂 𝐾 ⁄ = 3(1 -2𝑣) (Fig. 1).

Numerical model of the nanoindentation tests

In this paper, a parametric finite element model is built using the ANSYS software, (Ansys 16.0, 2016) (Fig. 1). This 2D-axisymetric model allows the simulation of the nanoindentation test of a viscoelastic material through the ANSYS viscoelastic law. Five conical indenter tips are used, with equivalent half angles equal to 42.28°, 57°, 60°, 65 and 70.3° half angle and are assumed to be rigid. The indenter with 70.3° corresponds to the axisymmetric equivalent cone Berkovich indenter tip used experimentally [START_REF] Chen | On the uniqueness of measuring elastoplastic properties from indentation: the indistinguishable mystical materials[END_REF][START_REF] Fischer-Cripps | Nanoindentation[END_REF]. The Coulomb's friction law is used between contact surfaces and the friction coefficient has been fixed to 0.2. Note that for this indenter geometry the contact friction coefficient does not have any significant effect on the numerical results [START_REF] Bucaille | Determination of plastic properties of metals by instrumented indentation using different sharp indenters[END_REF]. Linear quadrangular elements with 4 nodes (Q4 PLANE182) are used. The size of the modeled sample is 60 times greater than the maximum nanoindentation depth ℎ 𝑚𝑎𝑥 in order to obtain realistic boundary conditions. The mesh size in the area right below the indenter is made finer than in the rest of the sample over a length 8 times greater than the ℎ 𝑚𝑎𝑥 , which makes possible to model the contact and to increase the precision on the result of the simulation. The mesh is progressively coarser when moving away from the indented area, making it possible to reduce the number of elements and thus reduce the computation time down to 30 minutes (on the sequential CPU mode, 4 processors, 2.5 GHz). The model has about 20,000 finite elements. The nodes belonging to the lower surface of the part of the modeled sample are clamped. For the convergence study, several simulations with refined meshes and time increments are performed with all indenter tips. It is found that from 6 elements below the cube corner indenter (48.28°), 9 elements for the indenters 57°, 60°, and 10 elements for the indenter 65° and the Berkovich indenter (70.3°) at ℎ 𝑚𝑎𝑥 , the mesh size has no influence on the obtained results. The nanoindentation test is simulated by two subsequent parts: loading and unloading. During the loading part, the indenter tip penetrates the specimen up to ℎ 𝑚𝑎𝑥 and during the unloading part, the indenter returns to the initial position. In each simulation, 1 to 5 iterations are required at each nanoindentation depth increment to reach an equilibrium configuration. At ℎ 𝑚𝑎𝑥 , about nine elements are in contact with the indenter tip.

Experimental test

A nanoindentation experimental test was performed at room temperature and relative humidity using an Ultra-Nanoindenter (UNHT) with a Berkovich indenter tip. This device allows to impose a force up to 50 mN (±0,1 𝜇N) and a displacement up to 40 𝜇m (±0,1 nm). The test has been carried out in displacement-controlled mode to a maximum value of ℎ 𝑚𝑎𝑥 ≈ 550 nm.

Repeating the same experiment 5 times, a standard deviation of 6µN is observed experimentally at ℎ 𝑚𝑎𝑥 . It consists of a load-unload cycle. Once the contact established, the indenter penetrates the sample at a quasi-constant nanoindentation depth rate ℎ ̇= 1000 nm/min until the maximum displacement ℎ 𝑚𝑎𝑥 . Then during the unloading phase the indenter is removed at the same rate until ℎ = 0. The force 𝑃 𝑒𝑥𝑝 (𝑡) and displacement ℎ 𝑒𝑥𝑝 (𝑡) are recorded during loading and unloading segments for a time 𝑡 𝑚𝑎𝑥 .

Numerical tests

In order to numerically investigate the parameters identifiability with different loading types and determine the better identification procedure, several nanoindentation tests at a depth of ℎ 𝑚𝑎𝑥 = 500 nm and time 𝑡 𝑚𝑎𝑥 are numerically carried out at eight nanoindentation depth rates ℎ ̇= 50, 100, 500, 1000, 2500, 5000, 10000 and 20000 nm/min for triangular and trapezoidal loading, and eight strain rates (ℎ ̇ℎ ⁄ = 0.0104, 0.0207, 0.1036, 0.2072, 0.5179, 1.0359, 2.0722 and 4.1458 s -1 ) for the exponential loading. Normalized time-displacement curves are given in Fig. 2. Other nanoindentation tests at a depth of ℎ 𝑚𝑎𝑥 = 500 nm and maximum load of 𝑃 𝑚𝑎𝑥 = 0.6 mN are numerically carried out in force-controlled mode with eight nanoindentation depth rates (60, 120, 600, 1200, 3000, 6000, 12000 and 24000 µN/min) for triangular loading (Fig. 3a). A monotonic loading test is then performed with superimposed sinusoidal loading to a maximum depth of ℎ 𝑚𝑎𝑥 = 510 nm over a loading time of 30 s and frequency of 4 Hz (Fig. 3b). In the following, the computer-generated load-displacement curves obtained from these different conditions will be referred to as pseudo-experimental data.

Inverse method for the material parameters estimation

The updating process of the numerical model based on the experimental (or pseudoexperimental) data allows to estimate one or more parameters values 𝛉 ̂ which minimize the difference between the force 𝑃(𝑡; 𝛉) resulting from the finite element simulation and the experimental data 𝑃 𝑒𝑥𝑝 (𝑡). The inverse problem is recast as the minimization problem of an objective function ω, which quantifies the difference between the numerical model and the experiment:

𝛉 ̂= argmin 𝛉 ω [𝑃(𝑡; 𝛉), 𝑃 𝑒𝑥𝑝 (𝑡)] (6)
The objective function ω is minimized by a local numerical optimization technique based on the Levenberg-Marquardt algorithm [START_REF] Levenberg | A method for the solution of certain non-linear problems in least squares[END_REF][START_REF] Marquardt | An algorithm for least-squares estimation of nonlinear parameters[END_REF] implemented in MIC2M software (F. [START_REF] Richard | MIC2M Software: Modélisation et Identification du Comportement Mécanique des Matériaux/Modeling and Identification of the Mechanical Behavior of Materials[END_REF]. A starting point is required in this algorithm. The influence of this point will be investigated. The objective function is defined in displacementcontrolled mode as [START_REF] Qasmi | Viscoelastic mechanical properties determined by nanoindentation tests and its numerical modelling of polypropylene modified by He+ particle implantation and e-irradiation[END_REF]:

ω(𝛉) = 1 2𝑇 ∑ [ 𝑃 𝑘 (𝛉) -𝑃 𝑘 𝑒𝑥𝑝 𝑃 𝑚𝑎𝑥 𝑒𝑥𝑝 ] 2 𝑇 𝑘=1 (7)
𝑇 = 1000 is the number of data points for each nanoindentation test, i.e. number of measured force values 𝑃 𝑘 (𝛉) = 𝑃(𝑡 𝑘 ; 𝛉) and 𝑃 𝑘 𝑒𝑥𝑝 = 𝑃 𝑒𝑥𝑝 (𝑡 𝑘 ), and 𝑃 𝑚𝑎𝑥 𝑒𝑥𝑝 is the maximum of the experimental nanoindentation force. 𝑇 is sufficiently large so that it does not influence the reported results.

The uncertainty ∆𝜃 𝑗 on the estimated value of 𝜃 𝑗 after the updating process can be obtained from the following equation in 𝜃 𝑗 = 𝜃 ̂𝑗 :

∆𝜃 𝑗 𝜃 𝑗 = √2ω[𝐇 ̅ -1 ] 𝑗𝑗 (8)
where 𝐇 ̅ is a dimensionless pseudo-hessian matrix computed by forward finite difference method [START_REF] Pac | Microstructural and elasto-plastic material parameters identification by inverse finite elements method of Ti (1-x) AlxN (0< x< 1) sputtered thin films from Berkovich nano-indentation experiments[END_REF][START_REF] Richard | Viscoelastic modeling and quantitative experimental characterization of normal and osteoarthritic human articular cartilage using indentation[END_REF]. The components of 𝐇 ̅ are given as:

𝐻 ̅ 𝑖𝑗 = 1 𝑇 𝜃 𝑖 𝜃 𝑗 𝑃 𝑚𝑎𝑥 2 ∑ 𝜕𝑃 𝑘 (𝛉) 𝜕𝜃 𝑖 𝜕𝑃 𝑘 (𝛉) 𝜕𝜃 𝑗 𝑇 𝑘=1 𝑖, 𝑗 = 1, … ,4 (9)
where 𝑃 𝑚𝑎𝑥 is the maximum of the numerical nanoindentation force.

Note that if the test is force-controlled, the objective function is formulated using the displacement response ℎ(𝑡; 𝛉), instead of 𝑃(𝑡; 𝛉).

When a combination of several nanoindentation tests 𝑛 is used in the updating process, the total objective function is given by the sum of the objective functions of all tests.

ω(𝛉) = ∑ [ 1 2𝑇 ∑ ( 𝑃 𝑘 (𝑒) (𝛉) -𝑃 𝑘 𝑒𝑥𝑝 (𝑒) 𝑃 𝑚𝑎𝑥 𝑒𝑥𝑝 (𝑒) ) 2 𝑇 𝑘=1 ] 𝑛 𝑒=1 (10)
The uncertainty becomes:

∆𝜃 𝑗 𝜃 𝑗 = √ 2 𝑛 ω[𝐇 ̅ -1 ] 𝑗𝑗 (11)
and the pseudo-hessian matrix 𝐇 ̅ is calculated from the following equation:

𝐻 ̅ 𝑖𝑗 = ∑ [ 1 𝑇 𝜃 𝑖 𝜃 𝑗 (𝑃 𝑚𝑎𝑥 (𝑒) ) 2 ∑ 𝜕𝑃 𝑘 (𝑒) (𝛉) 𝜕𝜃 𝑖 𝜕𝑃 𝑘 (𝑒) (𝛉) 𝜕𝜃 𝑗 𝑇 𝑘=1 ] 𝑛 𝑒=1 𝑖, 𝑗 = 1, … ,4 (12)

Parametric identifiability

The sensitivity analysis is performed to evaluate how, and to which extent, variations of the model input data influence the output data. In this analysis, the value of each parameter is changed by 0.1% with respect to its initial value in the study. Using a finite difference scheme, the sensitivity vectors are given by (for a single test):

𝑆 𝑘𝑗 = 𝜃 𝑗 𝑃 𝑚𝑎𝑥 𝜕𝑃 𝑘 𝜕𝜃 𝑗 ( 13
)
The sensitivity of the nanoindentation force to the parameter 𝜃 𝑗 can be computed as:

𝛿 𝑗 = 𝜃 𝑗 𝑃 𝑚𝑎𝑥 √ 1 𝑇 ∑ ( 𝜕𝑃 𝑘 𝜕𝜃 𝑗 ) 2 𝑇 𝑘=1 (14)
A parametric identifiability analysis is used to quantify the reliability of the estimated parameters. The completeness of data contained in the nanoindentation force is quantified by an 𝐼-index [START_REF] Richard | Viscoelastic modeling and quantitative experimental characterization of normal and osteoarthritic human articular cartilage using indentation[END_REF]. This index appears to be convenient to explore and investigate what are the optimal loading conditions to determine the four parameters of the material constitutive law. The analysis can be done before and after the updating process and therefore does not necessarily require the experimental measurements (only pseudo-experimental loading). The 𝐼-index is a measure of the conditioning of the inverse problem and is defined as [START_REF] Pac | Microstructural and elasto-plastic material parameters identification by inverse finite elements method of Ti (1-x) AlxN (0< x< 1) sputtered thin films from Berkovich nano-indentation experiments[END_REF]:

𝐼 = log 10 ( 𝜆 𝑚𝑎𝑥 𝜆 𝑚𝑖𝑛 ) > 0 (15)
where 𝜆 𝑚𝑎𝑥 and 𝜆 𝑚𝑖𝑛 are the maximum and minimum eigenvalue of the matrix 𝐇 ̅ at the considered calculation point 𝛉, respectively. The 𝐼-index is calculated as function of the sensitivity of parameters taken into account in the combination. In the case of combination of nanoindentation tests, I min and I max are the minimum and maximum of the 𝐼-index for the better and worse combination, respectively.

The lower the 𝐼-index, the better conditioned is the matrix, which means its inverse can be calculated with great accuracy. At the opposite, if the 𝐼-index is large, the matrix is considered as ill-conditioned. Some 𝐼-index values defining practical limits can be found in the literature [START_REF] Gujarati | Basic econometrics[END_REF]. For a chosen set of parameters 𝛉 and loading type, this 𝐼-index is calculated for the 11 parameter combinations (6 couples, 4 triplets and 1 quadruplet) and for various time subsets (load-unload, load only, unload only). This procedure allows to distinguish the potentially identifiable combinations (𝐼 ≤ 2) of material parameters from those which are not (𝐼 > 3).

Results and discussions

In this section, the results of the updating process of the FEM are presented in terms of estimated values of the parameters and associated uncertainties. The uniqueness/non-uniqueness of the viscoelastic properties of PP obtained from a single nanoindentation experimental triangular load-unload test is studied. The identifiability analysis is numerically performed using several loading types (trapezoidal, exponential, sinusoidal) and angles apex in order to propose an identification procedure for the four viscoelastic material parameters (𝐸, 𝑐 1 , 𝜈, 𝜂). A link between the identifiability results and the dissipated energy is shown. The addition of the plasticity in behavior law is carried out for a complete description of the material behavior.

Non-uniqueness of the solutions

Updating process from experimental data

Experimental data from the nanoindentation test realized on PP at 1000 nm/min using the triangular loading are used in equation (Eq.6). Three starting points 𝛉 (01) , 𝛉 (02) and 𝛉 (03) are chosen to initialize the minimization algorithm (Table 1) and solve the minimization problem. The evolution of the four parameters during the minimization process are shown in Fig. 4. The parameters 𝐸, 𝑐 1 and 𝜂 tend towards the same values whatever the starting point with acceptable uncertainties (about 15% for the viscosity 𝜂). However the obtained values for Poisson's ratio 𝜈 are multiple with an enormous uncertainty (about 140%). The value of the objective function ω is almost identical for all three cases and remains very low. Therefore there is non-uniqueness of the solution of the inverse problem. Fig. 5 illustrates that the obtained solutions generate almost the same P-h curve as the one obtained experimentally. Objective function ω (3) 5.18 × 10 -3 1.24 × 10 -5 Fig. 4. Evolution of the four parameters (𝐸, 𝑐 1 , 𝜈, 𝜂) during the updating process using three starting points 𝜽 (01) , 𝜽 (02) and 𝜽 (03) . 𝛉 ( 01) 𝛉 ( 02) 𝛉 ( 03) 𝛉 ( 01) 𝛉 ( 02) 𝛉 ( 03) 𝛉 ( 01) 𝛉 ( 02) 𝛉 ( 03) 𝛉 ( 01) 𝛉 ( 02) 𝛉 (03

) 𝛉 ̂(1) 𝛉 ̂(2) 𝛉 ̂(3)
Furthermore, the Poisson's ratio is badly estimated comparing the obtained value with those given in the literature (𝜈~0.3 to 0.42) [START_REF] Gao | Characterisation of interphase nanoscale property variations in glass fibre reinforced polypropylene and epoxy resin composites[END_REF][START_REF] Jakes | Experimental method to account for structural compliance in nanoindentation measurements[END_REF]. Updating processes can also be performed by imposing one of the four parameters. Two particular cases are considered. In the first case A, the Young's modulus 𝐸 is set to 1.50 GPa and the three parameters 𝜈, 𝑐 1 and 𝜂 are estimated (starting point: 𝛉 (04) = (𝐸 = 1.50 GPa (imposed), 𝑐 1 = 12.25 GPa, 𝜈 = 0.4, 𝜂 = 65 GPa. s)). In the second case B, the value of the Poisson's ratio 𝜈 is set to 0.4 and 𝐸, 𝑐 1 and 𝜂 are estimated (starting point: 𝛉 (05) = (𝐸 = 1.50 GPa, 𝑐 1 = 12.25 GPa, 𝜈 = 0.4 (imposed), 𝜂 = 65 GPa. s)). Fig. 6 shows the evolution of the four parameters during the minimization process. The obtained values are multiple except for the parameter c 1 . The parameters estimated resulting from the updating process are presented in Table 2. It can be noted that imposing a parameter increases considerably the uncertainties even on 𝐸 and c 1 .The objective function values are not identical for the two solutions but they are the same order in magnitude. The five solutions 𝛉 ̂(1) , 𝛉 ̂(2) , 𝛉 ̂(3) , 𝛉 ̂(4) and 𝛉 ̂( 5) are summarized in Table 3. All of the P-h curves are in good agreement with the experimental data (Fig. 7). It is interesting to note that Young's modulus 𝐸 = 1.65 GPa obtained from the nanoindentation experimental test at 1000 nm/min using the Oliver-Pharr method [START_REF] Oliver | An improved technique for determining hardness and elastic modulus using load and displacement sensing indentation experiments[END_REF]) using a Poisson's ratio 𝜈 = 0.4 is approximately 12% higher than this estimated by FEMU (case B). 𝛉 (01) 𝛉 ( 01) 𝛉 ( 01)

𝛉 (01) (a) (b) (c) (d) 𝛉 ̂(1) 𝛉 ̂(4) 𝛉 ̂(5)

Sensitivity analysis

To probe the uniqueness of a solution resulting from an updating process may appear to be very time consuming. In fact, an updating process for four parameters requires about 40 FEM computations for each starting point and for each experimental loading. Another way is to investigate, a priori, the identifiability of the parameters. In this way, it is needed to estimate the sensitivity of the nanoindentation force to the material parameters. A sensitivity analysis is performed for the five solutions 𝛉 ̂(1) , 𝛉 ̂(2) , 𝛉 ̂(3) , 𝛉 ̂(4) and 𝛉 ̂(5) (2.5 hours for each solution) using the numerical nanoindentation triangular test realized at 1000 nm/min in displacementcontrolled mode. The norm of sensitivity vectors calculated using equation (Eq.14) are displayed in Fig. 8. The Young's modulus 𝐸 is the most sensitive parameter, which means that it is the most influential to the nanoindentation curve. The sensitivities appear to be well balanced for each solution. An exception can be noted for the Poisson's ratio. In fact, there is a lack of sensitivity to 𝜈 for the solutions 𝛉 ̂(1) , 𝛉 ̂(2) , 𝛉 ̂(3) and 𝛉 ̂(4) where the sensitivity to this coefficient is close to zero (Fig. 8). It is noted that the sensitivity of the nanoindentation force to the friction ratio can be neglected. 3) using triangular test at 1000 nm/min.

𝛉 ̂(1) 𝛉 ̂(2) 𝛉 ̂(3) 𝛉 ̂(4) 𝛉 ̂(5)
The sensitivity vectors of the nanoindentation force have thus been calculated using equation (Eq.13) for the solution 𝛉 ̂(5) . Fig. 9a shows that these vectors are very similar (up to a multiplicative factor) during the loading phase. The identification of these four parameters from the sole loading phase is thus expected to be very difficult. Focusing on the sensitivity to 𝐸 we note that during the unloading part the proportionality with 𝑐 1 and 𝜂 is lost. It seems to indicate that the unloading segment is more suitable to distinguish solution parameters and so to identify them. It also appears that the sensitivity vectors to 𝐸 and 𝜈 are almost collinear during the loading segment (OA) in (Fig. 9b), thus their identification using a single nanoindentation loading data is impossible. This is understood as a consequence of the sensitivity of the Boussinesq's problem to the sole 𝐸 (1 -𝜈 2 ⁄ ) parameter [START_REF] Boussinesq | Application des potentiels à l'étude de l'équilibre et du mouvement des solides élastiques: principalement au calcul des déformations et des pressions que produisent, dans ces solides, des efforts quelconques exercés sur une petite partie de leur surface ou de leur intérieur: mémoire suivi de notes étendues sur divers points de physique[END_REF].

A posteriori identifiability analysis

To quantify the completeness of the data used in the updating process, the 𝐼-index is determined from equation (Eq.15) using the solution 𝛉 ̂(5) with the numerical nanoindentation triangular test realized at 1000 nm/min for the 11 possible combinations of parameters as a function of the considered data subsets (load-unload, load and unload) and summarized in Table 4. The identification of all combinations of two parameters is possible from load-unload test and only unloading segment. It seems possible to identify three parameters (𝐸, 𝑐 1 , 𝜂) when the Poisson's ratio is known. The value of the index 𝐼(𝐸, 𝑐 1 , 𝜂) is greater than 3 only when the sole loading segment is considered and less than 2 as soon as the unloading phase is taken into account. For this combination of parameters, the identifiability is better if only unload is considered (𝐼 = 1.5) than if both load and unload are considered (𝐼 = 1.9). It can be observed that the relevant information is therefore contained in the unloading part, confirming the conclusions drawn from the sensitivity vectors analysis. The 𝐼-index is higher than 3 for the combination of four parameters from single nanoindentation load-unload test. The identification of the full set of parameters is considered to be impossible because of the inverse problem is too ill-posed.

A priori identifiability analysis of experimental approaches based on a single nanoindentation test

In this section, we focus on the possibility to identify the viscoelastic parameters from a single nanoindentation test. The effects of nanoindentation depth rate, loading type (triangular, trapezoidal, exponential and sinusoidal) and the indenter tip angle are investigated.

Effect of the nanoindentation depth rate

The 𝐼-index is calculated for the other numerical nanoindentation triangular load-unload tests (50, 100, 500, 2500, 5000, 10000, 20000 nm/min) with the solution 𝛉 ̂(5) . Table 5 shows that it is possible to identify three parameters when using a nanoindentation depth rate between 500 and 1000 nm/min (𝐼 ≤ 2).The nanoindentation depth rate does not improve the four material parameter 𝐼-index (𝐼 > 3). It can be concluded that whatever the nanoindentation depth rate, in the considered range the identification of the four material parameter from single nanoindentation triangular test is impossible.

Table 5. 𝐼-index for all combinations of parameters for all nanoindentation triangular load-unload tests using the solution 𝜽 ̂(5) . 𝐼 ≤ 2 (green, potentially identifiable), 2 < 𝐼 ≤ 3 (difficult to identify), 𝐼 > 3 (red, not identifiable).

Combination 50 nm/min 100 nm/min 500 nm/min 1000 nm/min 2500 nm/min 5000 nm/min 10000 nm/min 20000 nm/min

Effect of the loading type

The identifiability analysis is also performed using the trapezoidal, exponential and sinusoidal numerical tests (Table 6). Three parameters (𝐸, 𝑐 1 , 𝜂) are potentially identifiable (𝐼 ≤ 2) when the Poisson's ratio is known using trapezoidal loading. Comparing the results from single nanoindentation test in terms of 𝐼-index, the addition of a plateau just after the loading phase does not appear very helpful. It is observed that the better identifiability of the three parameters (𝐸, 𝑐 1 , 𝜂) is obtained using the exponential loading. Finally, the addition of a sinusoidal signal at the chosen frequency to the triangular loading phase is no more interesting for this material and all combinations of parameters are in the best case difficult to identify. Comparing the values of the 𝐼-index, it is found that the angle of the indenter does not have a great influence on the identifiability. The 𝐼-index values are almost the same and the four material parameters are, in the best case difficult to identify.

Table 7. I-index for all combinations of parameters using the nanoindentation load-unload test of 500 nm/min with five indenter angles tips for the solution 𝜽 ̂(5) . 𝐼 ≤ 2 (green, potentially identifiable), 2 < 𝐼 ≤ 3 (difficult to identify), 𝐼 > 3 (red, not identifiable).

Link between dissipation and identifiability

The loss factor (intrinsic damping) tan(𝛿), which is used to measure a viscoelastic material property in the case of tensile and harmonic loading can be approached for any of the considered loading signals considered herein, by discording any signal harmonic so that:

tan(𝛿) = 2𝜋𝑓𝐸𝜂 𝑐 1 (𝐸 + 𝑐 1 ) + (2𝜋𝑓) 2 𝜂 2 (16)
where 𝑓 is the fundamental frequency. In the case of the nanoindentation test 𝑓 = 1 𝑡 𝑚𝑎𝑥 ⁄ .

The 𝐼-index is calculated for the numerical nanoindentation triangular tests with equivalent cube corner and Berkovich indenter tips using displacement-controlled and force-controlled modes to investigate the effect of nanoindentation rate.

Fig. 10 presents the results for three combinations of parameters calculated using the solution 𝛉 ̂(5) with the equivalent cube corner indenter tip. It can be observed that the 𝐼-index values of the four parameters fall for the three slowest rates (60, 120 and 600 µN/min) (Fig. 10b). This numerical problem disappears if the Poisson's ratio is known. The better identifiability of the material parameters is obtained with a nanoindentation depth rate between 500 nm/min (𝑡 𝑚𝑎𝑥 = 120 𝑠) and 1000 nm/min (𝑡 𝑚𝑎𝑥 = 60 𝑠) in displacement-controlled mode and a nanoindentation load rate 600 µN/min (𝑡 𝑚𝑎𝑥 = 120 𝑠) and 1200 µN/min (𝑡 𝑚𝑎𝑥 = 60 𝑠) in force-controlled mode. It also shows that the 𝐼-index is well correlated to the loss factor (Fig. 10). In particular, we note that when the loss factor is maximum, (i.e when the dissipation is maximum) the 𝐼-index is minimum, (i.e the identifiability is the maximum). For the equivalent Berkovich indenter tip (Fig. 11), it can be seen that the better identifiability of the material parameters is obtained with a nanoindentation depth rate between 500 nm/min (𝑡 𝑚𝑎𝑥 = 120 𝑠) and 1000 nm/min (𝑡 𝑚𝑎𝑥 = 60 𝑠) in displacement-controlled mode and a nanoindentation load rate about 1200 µN/min (𝑡 𝑚𝑎𝑥 = 60 𝑠) in force-controlled mode. It also shows that the 𝐼-index is well correlated to the loss factor (Fig. 11), similarly to equivalent cube corner indenter tip. The link between tan(𝛿) and the 𝐼-index is also investigated for trapezoidal and exponential loading types using equation (Eq.16). Fig. 12 shows the 𝐼-index results for 3 combinations of material parameters with equivalent Berkovich indenter tip as function of the loss factor using the solution 𝛉 ̂(5) . It can be seen that, whatever the loading type, the better identifiability corresponds to conditions which maximize the loss factor. The link between tan(𝛿) and the 𝐼index depends on the nanoindentation rate to solicit the dissipative phenomena, which is different from a loading to the other.

Combining tests for a well-posed inverse problem

In order to determine the better identifiability of the four material parameters, the combination of numerical nanoindentation triangular tests and indenter angle tips is investigated. The combination of the equivalent 2D cube corner ( = 42.28°) and Berkovich ( = 70.3°) indenter tips is considered to determine a unique solution for the inverse problem.

Dual nanoindentation

In the following, the combination of the indenter tip angles is used in order to determine the better approach to identify the four viscoelastic material properties. Table 9 presents the 𝐼-index values for all dual nanoindentation using triangular numerical test at 500 nm/min in displacement-controlled mode with the solution 𝛉 ̂(5) . Whatever the combination, it is possible to identify two parameters, however, for three parameters, the indenter tips have to be carefully chosen. The value of the 𝐼-index for the combination of four parameters from the combination of equivalent cube corner ( = 42.28°) and Berkovich ( = 70.3°) indenter tips is equal to 2.

The identification of all combinations of parameters using loading and unloading phases may thus provide a unique solution. The identifiability is difficult from unloading phases and impossible if only loading phases are considered. Fig. 13 displays the simulated P-h curves for the two indenter tips at 500 nm/min for the solution 𝛉 ̂(5) .

Fig. 13. Nanoindentation pseudo-experimental (P-h) curves for equivalent cube corner and Berkovich indenter tips at 500 nm/min obtained using the solution 𝜽 ̂(5) = (𝐸 = 1.47 GPa, 𝑐 1 = 0.94 GPa, 𝜈 = 0.4, 𝜂 = 17.08 GPa.s).

Table 9. 𝐼-index for all combinations of parameters for all dual nanoindentation data at 500 nm/min with the solution 𝜽 ̂(5) . 𝐼 ≤ 2 (green, potentially identifiable), 2 < 𝐼 ≤ 3 (difficult to identify), 𝐼 > 3 (red, not identifiable). The results of the identifiability analysis suggest that dual nanoindentation technique may provide a unique solution for the full set of four unknown parameters. Three updating processes are thus performed using equivalent cube corner and Berkovich indenter tips with the pseudoexperimental tests at 500 nm/min (Fig. 13). These updating processes make use of the three starting points 𝛉 (01) , 𝛉 (02) and 𝛉 (03) . As expected from 𝐼-index results, the obtained solutions are the same for the three starting points (Table 10) and almost equal to the solution 𝛉 ̂(5) . The evolution of the four parameters during the algorithm iterations is illustrated in Fig. 14. Objective function ω (3) 1.17 × 10 -1 2.79 × 10 -5 Fig. 14. Evolution of the 4 parameters (𝐸, 𝑐 1 , 𝜈, 𝜂) during the updating process using the three starting points 𝜽 (01) , 𝜽 (02) and 𝜽 (03) with dual nanoindentation technique.

Combination

The I-index results using two nanoindentation numerical tests performed at nanoindentation depth rate of 500 nm/min with equivalent cube corner and Berkovich indenter tips indicate that it is possible to identify the four parameters (𝐼(𝐸, 𝑐 1 , 𝜈, 𝜂) = 2) (Table 9). This value of the Iindex indicates that the inverse problem is correctly conditioned (i.e not very sensitive to noise).

To illustrate the effect of noise on the updating process, a white Gaussian noise is used to disrupt the nanoindentation forces 𝑃 from the numerical tests at nanoindentation depth rate of 500 nm/min using equivalent cube corner and Berkovich indenter tips. For a single test, the disrupted force 𝑃 ̃ is given by: 𝑃 ̃(𝑡) = 𝑃(𝑡) + 𝒩(0, 𝑠)𝑃 𝑚𝑎𝑥 (17) where 𝒩(0, 𝑠) is a normal (Gaussian) distribution with zero mean and standard deviation 𝑠. 𝑃 𝑚𝑎𝑥 is the maximum value of the force 𝑃(𝑡).

𝛉 ( 01) 𝛉 ( 02) 𝛉 ( 03) 𝛉 ( 01) 𝛉 ( 02) 𝛉 ( 03) 𝛉 ( 01) 𝛉 ( 02) 𝛉 ( 03) 𝛉 ( 01) 𝛉 ( 02) 𝛉 ( 03 Both the nanoindentation forces obtained with 500 nm/min using equivalent cube corner and Berkovich indenter tips are corrupted by noise according to equation (Eq.17) (Fig. 15). The estimated solutions (Eq.6) and the uncertainties (Eq.11) for three noise standard deviations are presented in Table 11. The four material parameters tend towards the same values whatever the starting point (Fig. 16) and these values are close to the reference solution (lower than 5%). The solution is not very sensitive to this type of noise, thereby proving the proposed procedure is adequate to retrieve a unique set of viscoelastic parameters.

Table 11. Estimated solutions for the three levels of measurement noise using dual nanoindentation with the starting point 𝜽 (02) . Reference 𝜽 ̂(5) = (𝐸 = 1.47 GPa, 𝑐 1 = 0.94 GPa, 𝜈 = 0.4, 𝜂 = 17.08 GPa.s). Fig. 16. Evolution of the 4 parameters (𝐸, 𝑐 1 , 𝜈, 𝜂) during the updating process for the three starting points of the minimization algorithm 𝜽 (01) , 𝜽 (02) and 𝜽 (03) using noisy force values of dual nanoindentation data.

Noise

Estimation of the viscoelastic law parameters from dual nanoindentation experimental data

In this section, two experimental nanoindentation tests (𝑛 = 2) carried out using the cube corner and Berkovich indenter tips at depth rate of 500 nm/min are considered. For the updating process, the starting point 𝛉 ( 01) is chosen to initialize the minimization algorithm. The identified viscoelastic (VE) parameters are summarized in Table 12. Unlike the updating process using single test for the viscoelastic behavior (Table 1), the Poisson's ratio tends to the limit value (𝜈 = 0.5). The nanoindentation experimental P-h curves and the numerical results are plotted in Fig. 17. The results show a poor agreement between the experimental data and the updating process results. As it can be seen in the optical images displayed in Fig. 18, residual imprints can be observed. This is a clear indication that plasticity occurs during the indentation of PP. It seems necessary to assess the impact of plastic deformation on the identified VE parameters.

Viscoelastic-plastic law

For a complete description of the material behavior, plasticity should be taken into account. In order to identify the mechanical behavior of PP, viscoelastic-plastic (VEP) behavior has been considered. In this case, the dissipation potential 𝛺 is defined as follows:

𝛺 = 𝐸 2𝜂 (𝛔 -𝐗 𝑎𝑛 ): 𝐒: (𝛔 -𝐗 𝑎𝑛 ) + 1 2𝐾 〈𝑓〉 2 (18)
where 𝐾 is the viscosity coefficient in the viscoplastic domain. In this case, 𝐾 is set to 1 MPa.s to enable the plastic deformation.

The function 𝑓 is given by:

𝑓(𝛔) = σ ̅ -𝜎 𝑦 (19)
where 𝜎 𝑦 is the yield strength and σ ̅ is the von Mises stress. The symbols 〈 〉 denote Macaulay's brackets such as 〈𝑓〉 = 0 if 𝑓 < 0 and 〈𝑓〉 = 𝑓 if 𝑓 ≥ 0.

This viscoelastic-plastic behavior law is controlled by five material parameters, which define the parameter set 𝛉 = (𝜃 1 , 𝜃 2 , 𝜃 3 , 𝜃 4 , 𝜃 5 ) = (𝐸, 𝑐 1 , 𝜈, 𝜂, 𝜎 𝑦 ). The updating process is presented in the following paragraphs.

Results of the updating process for the viscoelastic-plastic law

In this case, the same dual nanoindentation tests (cube corner, Berkovich) are used and three starting points (𝛉 (06) , 𝛉 (07) , 𝛉 (08) ) are considered (Table 13). The starting point 𝛉 (06) is built by including the yield strength of the PP [START_REF] Ashby | Materials selection in mechanical design[END_REF] with the estimated solution for the viscoelastic behavior 𝛉 ̂(5) (𝛉 (06) = (𝐸 = 1.47 GPa, 𝑐 1 = 0.94 GPa, 𝜈 = 0.4 (imposed), 𝜂 = 17.08 GPa.s, 𝜎 𝑦 = 60 MPa). The updating process lasts about 23 days of computation for each starting point. The minimization starting from 𝛉 (06) has been performed under the constrain 𝜈 = 0.4; whereas no constrain was imposed for the starting points 𝛉 (07) and 𝛉 (07) . It can be observed that the Poisson's ratio tends to the 𝜈 = 0.5 when 𝜈 is free. The parameters 12), one can conclude that the viscoelastic-plastic law is more adequate to the description of the material behavior. Then, it is obvious that the numerical and experimental P-h curves are very close for both indenter tips (Fig. 19). Objective function ω (8) 5.39 × 10 -2 1.88 × 10 -3 Fig. 19. Experimental (ℎ ̇= 500 nm/min) and simulated nanoindentation P-h curves using cube corner and Berkovich indenter tips for the three viscoelastic-plastic solutions (Table 13).

Another updating process is also performed by imposing the yield strength 𝜎 𝑦 = 60 MPa and the Poisson's ratio 𝜈 = 0.4. The three viscoelastic parameters 𝐸, 𝑐 1 and 𝜂 are estimated using the starting point: (𝛉 (09) = (𝐸 = 1.47 GPa, 𝑐 1 = 0.94 GPa, 𝜈 = 0.4 (imposed), 𝜂 = 17.08 GPa.s, 𝜎 𝑦 = 60 MPa (imposed)). It can be observed that imposing the yield strength increases the values of the instantaneous 𝐸 and anelastic 𝑐 1 moduli and decreases the viscosity coefficient 𝜂 (Table 14). The obtained P-h curves are in good agreement with the experimental data compared to those in viscoelastic case (Fig. 20). This means that the addition of a yield strength, even when is set to a wrong value, allows a better prediction of the viscoelastic part of the PP behavior. The identifiability analysis shows that when the yield strength σ y and Poisson's ratio 𝜈 are known, the identifiability of the three parameters 𝐸, 𝑐 1 and 𝜂 is possible 𝐼(𝐸, 𝑐 1 , 𝜂) = 2.0. The identification of the five viscoelastic-plastic parameters from this dual nanoindentation is difficult (𝐼(𝐸, 𝑐 1 , ν, 𝜂, 𝜎 𝑦 ) = 2.8). Fig. 20. Experimental (ℎ̇ = 500 nm/min) and simulated nanoindentation (P-h) curves with cube corner and Berkovich indenter tips for the viscoelastic and viscoelastic-plastic laws.

Assessment of the viscoelastic-plastic behavior

In order to validate the identified behavior, tensile tests have been performed using PP dumbbell-shaped specimens with gauge section (20 × 4 × 0.5 mm 3 ). Fig. 21 shows the repetitive progressive loading test carried out at 2 N/s. The full-field strains have been measured by global digital image correlation (DIC) assuming an affine displacement field over the region of interest in both longitudinal (3) and transversal (2) directions. The comparison between the simulation results obtained using the solution 𝛉 ̂(6) = (𝐸 = 1.77 GPa, 𝑐 1 = 1.43 GPa, 𝜈 = 0.4 (imposed), 𝜂 = 15.75 GPa.s, 𝜎 𝑦 = 76.70 MPa) and the tensile test data illustrates that the identified viscoelastic-plastic law predicts the PP behavior better than the viscoelastic one using the solution 𝛉 ̂(5) (Fig. 21).

Conclusion

In this study, the uniqueness\non-uniqueness of the viscoelastic properties of polypropylene determined by the FEMU of the nanoindentation test is studied. A four-parameter viscoelastic behavior law has been implemented in a 2D axisymmetric finite element model. The FEMU process of nanoindentation test illustrates that a single nanoindentation experimental triangular load-unload test realized at constant nanoindentation depth rate (~1000 nm/min) is not sufficient to uniquely determine the four viscoelastic properties of the material. The updating process of the FE 2D-axisymmetric model on these experimental data leads to multiple solutions for the values of the four parameters (𝐸, 𝑐 1 , 𝜈, 𝜂) of the viscoelastic behavior law. The identification of the four parameters from a single nanoindentation triangular test is thus considered as impossible. The identifiability analysis allows the quantification of the ill-posed character of the inverse problem by a scalar 𝐼-index and shows that it is possible to identify three parameters (𝐸, 𝑐 1 , 𝜂) even if only the unloading phase is taken in account, which proves that the relevant information is in the unloading phase.

The effect of nanoindentation depth rate, loading type (triangular, trapezoidal, exponential, sinusoidal) and indenter tip angle on the identifiability is numerically investigated. The comparison between the results from the different loading types shows that the identification of the four material parameters is not possible. The included half angle of the indenter tip does not have a significant influence on the identifiability results.

The comparison between the loss factor and the 𝐼-index results from the nanoindentation triangular tests indicates that the better identifiability of the material parameters is obtained at the maximum loss factor, which corresponds to the maximum of the dissipated energy. It is also observed that, whatever the loading type, the best identifiability is obtained if the loss factor is maximum.

The combination of several triangular load-unload tests improves the identification robustness and does not lead to better 𝐼-index for the four material parameters compared to a single but properly chosen one. We show that the combination of two numerical nanoindentation triangular tests carried out at a constant nanoindentation depth rate using equivalent cone apex angles of cube corner (42.28°) and Berkovich (70.3°) indenter tips allows for the retrieval of a unique solution of the inverse problem, which is robust with respect to the noise. The four material parameters are potentially identifiable using this experimental protocol if the material behavior is viscoelastic. These results are agree with those obtained for elastoplastic behavior.

The viscoelastic-plastic law with five parameters (𝐸, 𝑐 1 , 𝜈, 𝜂, 𝜎 𝑦 ) improves the description of the material behavior. The identifiability results show that if the yield strength and Poisson's ratio 𝜈 are known, the identification of the three viscoelastic parameters (𝐸, 𝑐 1 , 𝜂) from dual nanoindentation experimental data (cube corner, Berkovich) is possible (𝐼 = 2.0). The comparison between the simulation results of the identified viscoelastic and viscoelastic-plastic laws and the tensile test data shows that viscoelastic-plastic law is more adequate to the description of the PP behavior.

This work paves the way for the design of robust experimental protocols based on the 𝐼-index in order to identify more complex behaviors, such as the viscoelastic with non-constant Poisson's ratio, viscoelastic-plastic or viscoelastic-viscoplastic. Indeed, the 𝐼-index can be used to numerically design the nanoindentation tests which allow to activate the dissipative phenomena as much as possible, thus to identify intrinsic and reliable properties.

Fig. 1 .

 1 Fig. 1. 2D-axisymmetric model of the nanoindentation test and the rheological model decomposed into: a deviatoric and a spherical part (𝐸 𝐺 ⁄ = 𝑐 1 𝐺 1 ⁄ = 𝜂 𝜂 𝐺 ⁄ = 2(1 + 𝑣)) and (𝐸 𝐾 ⁄ = 𝑐 1 𝐾 1 ⁄ = 𝜂 𝜂 𝐾 ⁄ = 3(1 -2𝑣)).

Fig. 2 .

 2 Fig. 2. Normalized time-displacement curves, (a) trapezoidal loading. (b) exponential loading.

Fig. 3 .

 3 Fig. 3. Normalized time-force curves, (a) triangular loading. Normalized time-displacement (b) sinusoidal loading.

Fig. 5 .

 5 Fig. 5. Experimental (ℎ ̇= 1000 nm/min) and simulated nanoindentation curves for the three solutions 𝜽 ̂(1) = (𝐸 = 1.63 GPa, 𝑐 1 = 1.05 GPa, 𝜈 = 0.13, 𝜂 = 18.56 GPa.s), 𝜽 ̂(2) = (𝐸 = 1.63 GPa, 𝑐 1 = 1.05 GPa, 𝜈 = 0.03, 𝜂 = 18.48 GPa.s) and 𝜽 ̂(3) = (𝐸 = 1.63 GPa, 𝑐 1 = 1.06 GPa, 𝜈 = 0.03, 𝜂 = 18.50 GPa.s).

Fig. 6 .Fig. 7 .

 67 Fig. 6. Evolution of the four parameters (𝐸, 𝑐 1 , 𝜈, 𝜂) during the updating process for the two particular cases A and B. Comparison with a case with 4 free parameters for the starting point (𝜽 (01) = 𝐸 = 1.5 GPa, 𝑐 1 = 12.25 GPa, 𝜈 = 0.4, 𝜂 = 65 GPa.s).

Fig. 8 .

 8 Fig.8. Sensitivity of the nanoindentation force 𝑃 to 𝜃 𝑗 for the five solutions 𝜽 ̂(1) , 𝜽 ̂(2) , 𝜽 ̂(3) , 𝜽 ̂(4) and 𝜽 ̂(5) (Table3) using triangular test at 1000 nm/min.

Fig. 9 .

 9 Fig. 9. (a) Sensitivity vectors of the nanoindentation force to the material parameters 𝜃 𝑗 during loading and unloading. (b) collinearity between sensitivity vectors 𝑆 𝑘1 and 𝑆 𝑘3 (sensitivities to 𝐸 and 𝜈, respectively).

Fig. 10 .

 10 Fig. 10. 𝐼-index for the nanoindentation triangular tests and the loss factor using the solution 𝜽 ̂(5) with equivalent cube corner indenter tip, (a) displacement-controlled mode. (b) force-controlled mode.

Fig. 11 .

 11 Fig. 11. 𝐼-index for the nanoindentation triangular tests and the loss factor using the solution 𝜽 ̂(5) with equivalent Berkovich indenter tip, (a) displacement-controlled mode. (b) force-controlled mode.

Fig. 12 .

 12 Fig. 12. 𝐼-index versus 𝑡𝑎𝑛(𝛿) using three loading types with the solution 𝜽 ̂(5) for three combinations of parameters, (a) 𝐼(𝐸, 𝜈). (b) 𝐼(𝐸, 𝑐 1 , 𝜂). (c) 𝐼(𝐸, 𝑐 1 , 𝜈, 𝜂).

)Fig. 15 .

 15 Fig. 15. Disrupted pseudo-experimental nanoindentation P-h curves (ℎ ̇= 500 nm/min) for three levels of noise with equivalent cube corner and Berkovich indenter tips obtained using 𝜽 ̂(5) .

Fig. 17 .Fig. 18 .

 1718 Fig. 17. Nanoindentation (P-h) experimental and simulated curves for the viscoelastic behavior from dual nanoindentation with cube corner and Berkovich indenter tips at 500 nm/min.

  1 , 𝜂 and 𝜎 𝑦 tend towards different values with almost the same objective function ω for all three cases. Comparing the objective function values to that obtained for the viscoelastic law (Table

Fig. 21 .

 21 Fig. 21. Experimental data of the tensile test and simulation response using the viscoelastic and viscoelastic-plastic laws.

Table 1 .

 1 Estimated parameters set 𝜽 ̂ (Eq. 6) using three starting points and uncertainties (Eq.8).

		Parameter	Starting value	Estimated value Uncertainty
		𝑗 𝜃 𝑗	𝜃 𝑗 (0)	𝜃 ̂𝑗	∆𝜃 𝑗 /𝜃 𝑗 (%)
		1 𝐸 (GPa)	1.50	1.63	7.0
	Starting point 1: 𝛉 (01)	2 𝑐 1 (GPa) 3 𝜈	12.25 0.4	1.05 0.13	6.0 90
		4 𝜂 (GPa. s)	65	18.56	15
	Objective function ω (1)		4.59 × 10 -2	1.24 × 10 -5
		1 𝐸 (GPa)	1.40	1.63	6.0
	Starting point 2: 𝛉 (02)	2 𝑐 1 (GPa) 3 𝜈	8.0 0.3	1.05 0.03	6.0 137
		4 𝜂 (GPa. s)	45	18.48	15
	Objective function ω (2)		2.18 × 10 -2	1.24 × 10 -5
		1 𝐸 (GPa)	1.20	1.63	7.0
	Starting point 3: 𝛉 (03)	2 𝑐 1 (GPa) 3 𝜈	4.0 0.2	1.06 0.03	6.0 138
		4 𝜂 (GPa. s)	30	18.50	15

Table 2 .

 2 Estimated parameters for the two particular cases A and B.

		Parameter	Starting value	Estimated value Uncertainty
		𝑗 𝜃 𝑗	𝜃 𝑗 (0)	𝜃 ̂𝑗	∆𝜃 𝑗 /𝜃 𝑗 (%)
		1 𝐸 (GPa)	1.50 (imposed) 1.50 (imposed)	10
	Starting point 4: 𝛉 (04)	2 𝑐 1 (GPa)	12.25	0.96	24
	(Case A)	3 𝜈	0.4	0.02	291
		4 𝜂 (GPa. s)	65	22.87	29
	Objective function ω (4)		4.59 × 10 -2	3.27 × 10 -5
		1 𝐸 (GPa)	1.50	1.47	44
	Starting point 5: 𝛉 (05)	2 𝑐 1 (GPa)	12.25	0.94	48
	(Case B)	3 𝜈	0.4 (imposed)	0.4 (imposed)	145
		4 𝜂 (GPa. s)	65	17.08	41
	Objective function ω (5)		4.59 × 10 -2	1.31 × 10 -5

Table 3 .

 3 Five solutions of the inverse problem.

		𝑗 1	2	3	4
	Solution	𝐸 (GPa)	𝑐 1 (GPa)	𝜈	𝜂 (GPa. s)	ω
	𝛉 ̂(1)	1.63	1.05	0.13	18.56	1.24 × 10 -5
	𝛉 ̂(2)	1.63	1.05	0.03	18.48	1.24 × 10 -5
	𝛉 ̂(3)	1.63	1.06	0.03	18.50	1.24 × 10 -5
	𝛉 ̂(4)	1.50 (imposed)	0.96	0.02	22.87	3.27 × 10 -5
	𝛉 ̂(5)	1.47	0.94	0.4 (imposed)	17.08	1.31 × 10 -5

Table 4 .

 4 𝐼-index for all combinations of parameters using triangular test at 1000 nm/min with the solution 𝜽 ̂(5) . 𝐼 ≤ 2 (green, potentially identifiable), 2 < 𝐼 ≤ 3 (difficult to identify), 𝐼 > 3 (red, not identifiable).

	Combination	Load-unload	Load	Unload
	𝐸, 𝑐 1	0.8	2.5	0.2
	𝑐 1 , 𝜂	0.8	2.1	0.8
	𝑐 1 , 𝜈	0.8	2.1	1.0
	𝐸, 𝜂	1.4	2.8	0.8
	𝐸, 𝜈	1.6	3.1	0.9
	𝜈, 𝜂	1.9	3.0	1.8
	𝐸, 𝑐 1 , 𝜂	1.9	3.5	1.5
	𝑐 1 , 𝜈, 𝜂	2.1	3.2	2.2
	𝐸, 𝜈, 𝜂	2.5	3.8	2.1
	𝐸, 𝑐 1 , 𝜈	3.0	3.7	2.6
	𝐸, 𝑐 1 , 𝜈, 𝜂	3.7	3.9	3.5

Table 6 .

 6 𝐼

	Combination Triangular	Trapezoidal	Exponential	Triangular	Sinusoidal
		load-unload	load-unload	load-unload	load	load
	𝐸, 𝑐 1	[0.6; 4.2]	[0.6; 4.2]	[0.6; 5.7]	2.5	2.5
	𝑐 1 , 𝜂	[0.6; 2.3]	[0.7; 2.5]	[0.7; 2.6]	2.1	2.1

-index for different loading types ([𝐼 𝑚𝑖𝑛 ; 𝐼 𝑚𝑎𝑥 ]): triangular, trapezoidal, exponential and a sinusoidal loading with max rate (1020 nm/min). 𝐼 ≤ 2 (green, potentially identifiable), 2 < 𝐼 ≤ 3 (difficult to identify), 𝐼 > 3 (red, not identifiable).

Table 10 .

 10 Estimated solutions for the three starting points using dual nanoindentation.

		Parameter	Starting value	Estimated value
		𝑗 𝜃 𝑗	𝜃 𝑗 (0)	𝜃 ̂𝑗
		1 𝐸 (GPa)	1.50	1.47
	Starting point 5: 𝛉 (05)	2 𝑐 1 (GPa) 3 𝜈	12.25 0.4 (imposed)	0.94 0.4 (imposed)
		4 𝜂 (GPa. s)	65	17.08
	Objective function ω (5)		4.59 × 10 -2	1.31 × 10 -5
		1 𝐸 (GPa)	1.50	1.45
	Starting point 1: 𝛉 (01)	2 𝑐 1 (GPa) 3 𝜈	12.25 0.4	0.95 0.39
		4 𝜂 (GPa. s)	65	17.40
	Objective function ω (1)		8.56 × 10 -1	2.79 × 10 -5
		1 𝐸 (GPa)	1.40	1.45
	Starting point 2: 𝛉 (02)	2 𝑐 1 (GPa) 3 𝜈	8.0 0.3	0.95 0.39
		4 𝜂 (GPa. s)	45	17.40
	Objective function ω (2)		4.40 × 10 -1	2.79 × 10 -5
		1 𝐸 (GPa)	1.20	1.45
	Starting point 3: 𝛉 (03)	2 𝑐 1 (GPa) 3 𝜈	4.0 0.2	0.95 0.39
		4 𝜂 (GPa. s)	30	17.40

Table 12 .

 12 Estimated parameters of the viscoelastic behavior using the dual nanoindentation (cube corner and Berkovich indenter tips).

		Parameter	Starting value	Estimated value
		𝑗 𝜃 𝑗	𝜃 𝑗 (0)	𝜃 ̂𝑗
		1 𝐸 (GPa)	1.50	1.03
	Starting point: 𝛉 (01)	2 𝑐 1 (GPa) 3 𝜈	12.25 0.4	0.17 0.5
		4 𝜂 (GPa. s)	65	49.17
	Objective function ω		2.24	6.71 × 10 -2

Table 13 .

 13 Estimated parameters set 𝜽 ̂ (Eq. 6) for the viscoelastic-plastic law using three starting points.

		Parameter	Starting value	Estimated value
		𝑗 𝜃 𝑗	𝜃 𝑗 (0)	𝜃 ̂𝑗
		1 𝐸 (GPa)	1.47	1.77
		2 𝑐 1 (GPa)	0.94	1.43
	Starting point 6: 𝛉 (06)	3 𝜈	0.4 (imposed)	0.4 (imposed)
		4 𝜂 (GPa. s)	17.08	15.75
		5 𝜎 𝑦 (MPa)	60	76.70
	Objective function ω (6)		9.78 × 10 -2	2.07 × 10 -3
		1 𝐸 (GPa)	4.0	1.61
		2 𝑐 1 (GPa)	3.0	1.31
	Starting point 7: 𝛉 (07)	3 𝜈	0.4	0.5
		4 𝜂 (GPa. s)	60	14.41
		5 𝜎 𝑦 (MPa)	90	69.96
	Objective function ω (7)		1.44	1.86 × 10 -3
		1 𝐸 (GPa)	2.0	1.57
		2 𝑐 1 (GPa)	1.5	1.32
	Starting point 8: 𝛉 (08)	3 𝜈	0.3	0.5
		4 𝜂 (GPa. s)	40	15.25
		5 𝜎 𝑦 (MPa)	70	70.28

Table 15 .

 15 Estimated solution for the viscoelastic-plastic law with imposing yield strength.

		Parameter	Starting value	Estimated value
		𝑗 𝜃 𝑗	𝜃 𝑗 (0)	𝜃 ̂𝑗
		1 𝐸 (GPa)	1.47	1.90
		2 𝑐 1 (GPa)	0.94	1.91
	Starting point 9: 𝛉 (09)	3 𝜈	0.4 (imposed)	0.4 (imposed)
		4 𝜂 (GPa. s)	17.08	10.80
		5 𝜎 𝑦 (MPa)	60 (imposed)	60 (imposed)
	Objective function ω (9)		9.78 × 10 -2	5.36 × 10 -3

Acknowledgements

This work was supported by the French Ministry of Higher Education and Research. We thank P. Delobelle for his valuable scientific advice in the field of nanoindentation.

Several nanoindentation triangular tests

The identifiability analysis is carried out using a set of 2 to 8 nanoindentation triangular numerical tests (50, 100, 500, 1000, 2500, 5000, 10000, 20000 nm/min) for the solution 𝛉 ̂(5) . Table 8 summarizes the 𝐼-index values for all subsets of nanoindentation tests. The combination of nanoindentation triangular tests decreases the maximum 𝐼-index values, and therefore improves the identification robustness. All combinations of two parameters can be identified from the subset of two nanoindentation tests. The subset of multiple nanoindentation tests allows to identify 9 combinations of parameters. The identification of the four material parameters is still impossible despite the additional information. The value of the 𝐼-index for the four material parameters is never less than 3.5. Comparing the 𝐼-index results from single nanoindentation test (Table 4) and the subset of several tests (Table 8), it is observed that the set of eight nanoindentation rates does not necessarily lead to better 𝐼-index compared to a single but properly chosen one.


T T T Nguyễn

T Doanh

A Le Bot

D Dalmas

Investigating the dynamic instabilities of model granular materials in isotropic consolidation and triaxial drained compression

This paper reports the dynamic instabilities of model granular materials under isotropic consolidation and triaxial drained axisymmetric compression. Loose and fully saturated samples of mono sized glass beads exhibit in isotropic drained compression a series of local collapses under undetermined stress and even liquefaction in some rare cases. Stick-slip phenomenon occurs in drained compression, and even rare liquefaction happens for the first slip. These dynamic instabilities (collapse, liquefaction and stick-slip) of loose granular assembly can share the same physical driving mechanisms with strong links to geometrical features, i.e. the collapse of the structural metastable honeycombed macropores or the propagation of the local failures of the force chains; even if the unambiguous identification of the triggering mechanisms is still unknown. The experimental data eliminates the excess pore fluid as the primary cause.

INTRODUCTION

Despite a large number of intensive theoretical and experimental studies, the mechanical behaviour of granular materials is currently far from fully being understood [10,2]. One common idea is to use model granular materials consisting of perfectly rounded particles to overcome the complexity of real ones in laboratory as well as in discrete numerical modelling.

The stick-slip phenomenon in compression shear on glass beads was rapidly observed in geotechnical and geophysical laboratories [1,3,4,5]; however, the precise triggering mechanisms were still missing. The isotropic drained compression is largely unstudied, leading to recent surprising discovery of isotropic liquefaction [6], the third kind of liquefaction in soil mechanics [7]. The unnecessary presence of external shearing stress for this liquefaction has thickened the mystery.

The aim of this paper is to explore some possible experimental links between these unexpected dynamic instabilities (collapse, liquefaction and stick-slip) of fully saturated model granular assembly. The temporal evolution of pore fluid pressure during the dynamic collapse and slip phase in triaxial testing is of particular interest. These links can contribute to a better understanding of the devastating liquefaction phenomenon during earthquake, still awaiting explanation.

EXPERIMENTAL METHOD

Measurement principle. In granular mechanics, the mechanical behaviour of fully saturated granular material is often studied in the versatile triaxial machine in two-step testing on a cylindrical sample. The isotropic drained compression (consolidation) is performed first, up to a desired effective isotropic stress σ r to get the compressibility characteristics; then the triaxial drained shear in compression to get the frictional angle and the dilatancy characteristics under constant σ r . The drained character of these two steps is characterised by a nil interstitial pore pressure U throughout the test duration, satisfying σ = σ, according to the effective stress principle σ = σ + U , σ is the total stress [9]. These simple compressibility and shear characteristics are widely accepted to be the minimum ones in any theoretical constitutive equation on granular material.

Measurement setup. Fig. 1 presents the classical experimental set-up for triaxial testing on a short cylindrical specimen of initial height H 0 and diameter D 0 inside a triaxial cell. The granular assembly is enclosed inside a cylindrical and open-ended latex membrane of 0.3 mm thickness. The back pressure U 0 needed for a full saturation is applied at the bottom of the sample. The pore-water pressure U is recorded on the top cap of the granular sample to ensure an homogeneous stress state.

Compressed air under constant stress rate is used in the first isotropic step up to a desired cell pressure σ r . If the sample is still geometrically stable at σ r , it will be vertically sheared under constant σ r and constant axial strain rate in the second shearing step. This second step adds an additional deviatoric shear stress q = σ a -σ r to the first purely isotropic step; σ a is the axial stress. The excess pore pressure ∆U = U -U 0 is kept small to respect the requirement of full drainage for these two steps (isotropic consolidation and drained compression). The global axial strain ε a = ∆h/H 0 is estimated with the measurement of axial displacement ∆h of the top platen inside the triaxial cell. The global volumetric strain is deduced from ε v = ∆v/V 0 , where V 0 is the initial sample volume, with the measurement of water volume ∆v expelled from or moved into the sample [8]. The volumetric strain is usually translated into void ratio e in the first isotropic step.

Materials and preparation.

Virgin monodisperse and spherical soda-lime glass beads (CVP Sil-glass) of mean diameter of 0.723 mm is used. Granular sample, 70 mm in height (H 0 ) and 70 mm in diameter (D 0 ), is prepared using a modified moist tamping and under compaction method. The combination of water-soluble gas CO 2 together with de-aired distilled water, and a high constant back pressure U 0 of 200 kPa is used to ensure a quasi-saturated state through high Skempton's coefficient B ≥ 0.95.

RESULTS AND ANALYSIS

Isotropic consolidation. Fig. 2a shows the typical results of isotropic consolidation (the first part) with a very different compressibility behaviour. Instead of expected smooth and continuous increase of density with increasing σ , three events with large unexpected drops ∆e in void ratio can be seen at undetermined triggering isotropic stress σ trig . Between two events, continuous decrease of e is observed. Event C with Drained triaxial compression. Fig. 2b andc give the typical results of drained compression shear. Again, instead of expected smooth stress-strain behaviour and continuous increase of q with increasing ε a for an initially loose granular material, stick-slip phenomenon was observed with numerous and quite irregular stick-slip events. Each event consists of a very abrupt deviatoric stress drop ∆q (slip-phase), followed by a gradual hardening (stick-phase) up approximatively to the previous stress level. Upon ignoring these stress drops, a nearly smooth envelope curve is retrieved, conforming to the usual behaviour of loose granular materials characterised by a monotonic increase of q to a horizontal critical state plateau at q max . The stick-slips occur throughout the studied axial strain range, up to 25 %, with five particularly events of large deviatoric stress drop, noted as SS4 to SS29. Rapid succession of precursors with small stress drop can be seen, as before SS29. The largest drop towards the isotropic stress state at 27 % of ε a is a simple unloading-reloading cycle with no stick-slip in unloading.

The corresponding volumetric strain exhibits strongly correlated variations with ∆q for the same values of ε a . Globally, this volumetric strain confirms the loose character previously noticed on the stress-strain relationship with a classical continuous volumetric decrease towards its constant critical state at 2 %. The sawtooth volumetric profile seems to oscillate around this critical value. The slip-phase is always associated with a sudden volumetric compaction and the stick-phase a gradual dilation. Five large compactions are obviously related to the five large deviatoric stress drops. Unexpected excess pore pressure. A closer look reveals that each event, isotropic (collapse and liquefaction) or deviatoric (stick-slip and liquefaction on first slip), commences with a sudden surge of ∆U , followed by a rapid dissipation within a few seconds towards the equilibrium state of constant back-pressure U 0 . It means an instantaneous reduction of effective stress σ , both axial and radial, according to the effective stress principle, followed by its gradual recover. It also means a very brief nonhomogenous stress state. In all events, the excess pore pressure in Fig. 3 vibrates like an oscillating underdamped system with a dominant frequency. ∆U was normalized by the brief stable value ∆U stable . Only selected events are shown to avoid clumsy figure. ∆U was measured simultaneously by two pore pressure transducers, one static (piezoresistive) and one dynamic (piezoelectric) with practically no differences.

The pore pressure evolution can be reasonably decomposed in three phases. The time origin is shifted to the beginning of the transient phase ± 0.1 ms, which is the current time resolution. First a fast transient phase I occurred within 200 ms at nearly constant volume. Then appears a second phase II valid only for liquefaction event, characterized by a stabilized value of ∆U stable for more than two seconds; and finally a last and longest dissipation phase III. The inferred liquefaction levels (solid circles) on the rising edge of ∆U of isotropic instabilities give the level of null effective stress. On liquefaction event it is exactly equal unity, meaning a null σ for the whole phase II, hence the terms liquefaction phase. Liquefaction event is obtained on separated sample.

Fig. 4 shows typical behaviour of isotropic collapses (event C) and compression stick-slip (event SS14), grouping together the detailed development of ∆U , ε v and ε a . The deviatoric stress is added in case of stick-slip. The liquefaction case with very large axial and volumetric deformation is left aside to emphasize the most common instabilities with the notable absence of stabilised phase of ∆U .

With the exception of sudden deviatoric stress drop ∆q to a constant level in less than 0.001 seconds, the temporal evolution of ∆U , ε v and ε a during the slip phase are quite similar to those of spontaneous collapse under isotropic consolidation. Instead of usual jumps of ε v , ε a and even ∆U in geotechnical and geophysical literature, this time-resolved study indicates a gradual dynamic modification of ε v and ε a from one to another steady level in less than 3 seconds, even less than 0.5 seconds for axial strain. Note the occurring of more than half of ε a during the fast transient phase I, the final contraction jump of ε a and the unexpected large extension state in between. These fast axial and volumetric strain rates suggest a dynamic regime for all instabilities.

Driving mechanisms. Fig. 4 also indicates one important insight into the dynamic slip. Since ∆U presents systematically an observed delay of about 3 ms comparing to the beginning of ∆q, it cannot be the first cause of the slip instability. Data from stick-slip experiment are giving preliminary indications for plausible driving mechanisms of the observed instabilities. The brevity of ∆q strongly suggests a possible failure of contact force chains, and ∆U only as a consequence thereof. However, it still plays a significant contributor to the amplification of these dynamic instabilities of fully saturated media by reducing the effective stress. Dynamic instabilities still affect unsaturated and dry granular media, with lesser amplitude.

The dynamic incremental anisotropic coefficient i dyn = ∆ε v /∆ε a and the spherical form of granular particles are highly suspected to be two major ingredients in creating these dynamic instabilities. Fig. 5 indicates the i dyn (black numbers) for selected events, estimated between the beginning and the end of each instability, once the excess pore pressure dissipated, regardless of the loading condition (isotropic or compression shear). Nearly parallel lines with constant i dyn in the dissipated phase III are obtained with different mean values.

These experimental observations point to a plausible conjectural scenario in the structural instability of metastable microstructure within the granular assembly, supplemented by some experimental evidences of cascading collapses of the force chain. The sudden breakage or buckling of the preferential force chains can induce local dynamic shear deformation, and can result in a fast pore pressure buildup due to the incompressibility of water in a fully saturated media. The rolling feature of spherical grain shape can greatly facilitate these chain failures. Then the temporal decrease of the effective stress would promote further deformation. In some rare cases, it can lead to catastrophic liquefaction if the null effective stress can be maintained for a prolonged period of more than one second [6]. Unfortunately, the precise physical reasons are still unknown.

CONCLUSIONS

The mechanical behaviour of fully water-saturated monodisperse glass beads, initially assembled in loose state by moist-tamping technique is investigated. Drained isotropic compression and deviatoric shear are affected by dynamic instabilities: isotropic collapse and liquefaction, deviatoric stick-slip and possibly liquefaction on the first slip. Each event, regardless of the loading condition (isotropic or compression shear) and depositional methods, is characterised by three phases. The first dynamic transient phase I has a fast drop of effective stress due to an excess pore pressure development which vibrates like an oscillating underdamped system. Local liquefaction with null σ is attained. The second phase II of stabilised excess pore pressure can lead to liquefaction failure, if prolonged for more than one second. The quasi-static dissipated phase III follows with gradual increase of axial and volumetric strain and a full recovery of σ to the previous steady state before event.

These unexpected dynamic instabilities probably share the same underlying triggering mechanisms and the excess pore pressure is just the amplification factor on a fully saturated system, not the cause. The speculated structural instability of metastable microstructure is one of the plausible conjectural scenarios. However, the precise driv-

Figure 1 .

 1 Figure 1. Experimental setup for isotropic consolidation and drained triaxial compression.

Figure 2 .

 2 Figure 2. Collapses in isotropic consolidation from 20 to 500 kPa (top) and stickslips in triaxial drained compression at 500 kPa of confining pressure (bottom).

Figure 3 .

 3 Figure 3. Three phases of pore water pressure development. Phase II is absent in isotropic collapses and compression stick-slips. largest volumetric compaction occurs at 476 kPa, just before attaining the required 500 kPa of confining pressure of the drained triaxial compression (the second part). Small red arrows indicate the time direction.

Figure 4 .

 4 Figure 4. Temporal behaviour of isotropic collapse (top, collapse C) and compression stick-slip (bottom, stick-slip 14) with axial extension in transient phase: evolution of deviatoric stress, excess pore pressure, axial and volumetric strain.

Figure 5 .

 5 Figure 5. Dynamic incremental anisotropic coefficient i dyn in the transient phase.

ing mechanisms remain unknown.


Yazan Barazi
email: barazi@laplace.univ-tlse.fr

François Boige
email: francois.boige@gmail.com

Nicolas C Rouger
email: rouger@laplace.univ-tlse.fr

Jean-Marc Blaquiere
email: blaquiere@laplace.univ-tlse.fr

Frédéric Richardeau
email: frederic.richardeau@laplace.univ-tlse.fr

'ig,vgs' Monitoring

Keywords: «Silicon Carbide (SiC)», «signal processing», «Robust control», «intelligent drive», «Protection device

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

The Silicon Carbide (SiC) MOSFETs offers several advantages, low switching losses, higher switching frequencies and high temperature stability [START_REF] Baliga | Silicon Carbide Power Devices[END_REF]. However, today, the short-circuit (SC) withstanding time (T SC ) of SiC MOSFETs is lower than silicon devices one [START_REF] Wang | Design and Performance Evaluation of Overcurrent Protection Schemes for Silicon Carbide (SiC) Power MOSFETs[END_REF], within , [START_REF] Romano | Influence of design parameters on the short-circuit ruggedness of SiC power MOSFETs[END_REF] instead of for Si IGBTs. Power converters designed with these components are therefore less robust. This weakness must be compensated by a protection whose delay must not exceed Tsc/2, as IGBT's standard today, or even less if repeated short-circuit robustness is required.

Power semiconductors devices are exposed to several types of short circuits. In summary, internal Hard switching fault (HSF, SC type1) and external fault under load (FUL, SC type2) are the two main fault behaviors. HSF occurs when the switch turns on with a permanent full bus voltage across itself typically caused by faulted control signals. FUL occurs when the device is already in on-state and an external short-circuit occurs. HSF and FUL can be distinguished by the presence or not of a dv/dt on the drain source of the MOSFETs. In all cases, MOSFETs saturates at high current density leading to a thermal dynamic source of heat in a few microseconds for the least robust components. Many kinds of SC detection methods have been demonstrated for Silicon IGBTs [START_REF] Chen | A Smart IGBT Gate Driver IC with Temperature Compensated Collector Current Sensing[END_REF] and SiC MOSFETs [START_REF] Sadik | Short-Circuit Protection Circuits for Silicon-Carbide Power Transistors[END_REF][START_REF] Awwad | Short-circuit evaluation and overcurrent protection for SiC power MOSFETs[END_REF]. Most of the proposed techniques rely on the drain-source voltage or current sensing. Those techniques are more suitable for bipolar transistors such as IGBTs, and on the other hand they are limited due to the use of a high-voltage rating diode or an additional current sensor. A relatively long blanking time is required in the case of the drain-source voltage detection method in order to achieve a stable state after the turn-on sequence [START_REF] Bertelshofer | Design Rules to Adapt the Desaturation Detection for SiC MOSFET Modules[END_REF]. The 2D diagnosis [START_REF] Bakran | Integrating the New 2D -Short circuit detection method into a power module with a power supply fed by the gate voltage[END_REF] is the best candidate to minimize the timing. Indeed, this family of methods applies during the turn-on sequence and not after, and is time independent. This study presents experimental results of two original detection methods of SiC MOSFETs under short circuit. Both methods are 2D diagnosis based-on low voltage / signals waveforms (I g , V gs ). The first detection method is developed, and depends on the unique behavior of SiC MOSFETs by the means of its gate leakage current runaway at high temperature at the on-state sequence. Where the second one is an adapted detection method, which depends on the gate charge variation on the switching cycle, between normal operation and short-mode. Those methods detect as early as possible the SC and turn softly the device off. In the full paper, the first method will be presented in general terms while the second one will be more detailed and will consist of the core of the article.

Behavior of SiC MOSFETs under short-circuit

Under short circuit fault, the MOSFET SiC presents a very high saturation current caused by the strong electrical field imposed in the channel region to minimize the R DSON in normal operation as reported in [START_REF] Baliga | Silicon Carbide Power Devices[END_REF]. The saturation current combined with the high electrical field in the depletion region dissipates an enormous amount of energy in few microseconds. The released heat increased, especially, the gate oxide temperature leading to hot electron injection inside the gate also called Schottky emission [START_REF] Boige | Physical origin of the gate current surge during short-circuit operation of SiC MOSFET[END_REF]. The MOSFET SiC behavior is displayed in Simplified illustration between normal-turn on NTO and SC-HSF operations. DUT : 1.2kV, 80m@25°C, R gext = 47, T case = 25°C

Detection using the dynamic gate leakage current

The dynamic gate leakage current can be observed in most MOSFET SiC available today and at about half the devices short-circuit withstanding time (Tsc/2) [START_REF] Boige | Gate leakage-current analysis and modelling of planar and trench power SiC MOSFET devices in extreme short-circuit operation[END_REF], Fig. 1-part A. Thus, detecting this current is a SC marker and, with the appropriate electronics, can lead to a fast shutdown to avoid failure.

However, detecting a current on the gate can be the normal switching order of device or a current injection through the parasitic Miller capacitance caused by the switching of another device as crosstalk phenomena in inverter leg. The proposed detection method relies on a low current detection and a dedicated logical circuit to sort out the normal operating condition from the short-circuit behavior. The detection circuit principle is presented in fig. 2(a). The current is measured using the external gate resistor as a shunt with an excellent SNR (Signal Noise Ratio) enabling a high robust detection method. The logical circuit uses the PWM switching signal as a reference starting a blanking time of 2 µs to avoid true false positive that is equivalent to Tsc/4 for the device considered. The proposed detection method has been implemented in an integrated way using SMD components.

The logical circuit is performed with a CPLD. The experimental results show the dynamic gate current detection at 10mA (100mV across 10Ω) and the device shut down in less than 150ns. The waveforms are presented in fig. 2(b). This protection can also handle type II "FUL" in much less than Tsc/2, making it a fast and attractive operational gate-driver. This type of fault operation is not described in this article.

Detection using gate charge method under Hard Switch Fault

Distinctively from the previous method which is not ultra-fast in type I HSF mode but robust in term of SNR criterion, the second method studied in this article is fast and can support high SNR. The method is not based on the gate leakage current monitoring, but still in the gate topology. The studied method depends on the fast integration of the gate current, the gate charge sequence at turn-on. The gate charge method already proved interesting results for IGBTs [START_REF] Oberdieck | Short circuit detection using the gate charge characteristic for Trench/Fieldstop-IGBTs[END_REF] but unclear results for SiC MOSFETs in terms of performances and practical PCB integration [START_REF] Horiguchi | A fast short-circuit protection method using gate charge characteristics of SiC MOSFETs[END_REF]. This study validates experimentally the gate charge method and fault-management for SiC MOSFETs using SMD (Surface-mount technology) components and embedded digital circuit in real full-voltage operation.

In the literature there are different architectures proposed for the gate charge method. In this study the circuit is quite the same in the concept. This gate charge method does not require high voltage diode for sensing or setting a detection period as the desaturation method. Moreover, this method has also, potentially, a high SNR. In the principle, this method requires a resettable integrator to estimate the amount of gate charge Q g ; and two comparators. The first comparator is to distinguish the difference between the charge amount under short circuit Q g-SC and normal turn-on (NTO) Q g-NTO; the second comparator comes to create a reading flag using only a v gs signal threshold crossing [START_REF] Barazi | Comparison between ig integration and vgs derivation methods dedicated to fast short circuit 2D diagnosis for wide band gap power devices[END_REF] and not a blanking time is required, fig. 3.

Fig. 3: Basic schematic of the detection and protection circuit using Q gate monitoring.

The gate charge in function of v GS under NTO and SC-HSF was studied by using fundamental equations and by Simulation using LT-Spice models, then compared to the datasheet, fig. 4(a). This 2D diagram is well known on the study of the gate charge method, offers two important detection keys: the detection zone, and the ratio S= Q g-SC /Q g-NTO .

In order to validate the gate charge method, an experiment was conducted for the gate charge monitoring method. The circuit in fig. 4 The circuit was validated under different V bus values between (0-400V). In fig. 6 & fig. 7 the experiment was under V bus =400V with a R g =10 Ω. Fig. 6 presents oscilloscope waveforms under NTO, fig. 6(a) gate and drain NTO waveforms. Under NTO the amount of the gate charge is higher than under SC, which activate Q g comparator, as shown in fig. 6(b), this activation will be scanned at the arrival of the V gs comparator, therefore, the flag remains low. In the other hand; fig. 7(a) shows the short-circuit effect (no drain-source dv/dt) on the drain and gate waveforms, the circuit is under SC, hence the amount of the gate charge is lower than NTO, the Q g comparator remains low, which at the arrival of the V gs comparator activate the short-circuit detection flag. Both figures were taken using the initial porotype board, on odder to present internal logic signals. The final board presented in Fig. 5 does not include internal logic signal. For a selective detection between NTO vs SC and a better ratio S, the integration can be controlled by adding another threshold level to start the integration. Unfortunately, this reduces the integrator output level, gives less robustness range for other power components and a trade-off must be reached.

C2M0080120D C3M0065090D Q G-NTO 4.2V 2.3V Q G-SC-HSF 1.8V 1V Q Ref 3V 1.7V IB= |Q Ref -Q G |/Q Ref NTO 0.4 0.35 SC 0.4 0.41 S % =Q G-SC-HSF /Q G-NTO (V/V) 43% 42%
S DS% (nC/nC) 34/71=> 48% 17.5/35 => 50% Fig. 9. "V GS image" under 400V NTO and SC operation for two power components, (a).C2M0080120 (b).C3M0065090 @R G_ext .=10.

The detection flag has been activated after 133ns including the buffer propagation time 10 ns with the initial prototype and 118ns with the final board. The detection time does not depend only on how fast our system is, the detection time depends obviously on the value of R g , the input capacitance of the power device, slew rate dv gs /dt and how low the reference level is defined (>V gsM ), etc.. . The slew rate between NTO vs SC change as presented in fig. 67, the V gs comparator under NTO arrives after SC, which is a positive thing to detect faster. Same as the reference levels, the threshold gate level was defined at a robust level 18V, with a lower level the detection will be earlier but with a trade-off on the SNR.

Protection of the device

After the detection of the fault, the SiC MOSFET must be turned off safely, to stop I ds from increasing and to avoid a dangerous over-voltage. Therefore, to protect the circuit safely, a soft shut down SSD system is included in the SMD. If the protection is not well designed it can lead to breakdown voltage due to high negative dI ds /dt caused by the total stray inductance, the protection should be soft to prevent this phenomenon. After activation of the D flip-flop (Detection flag), the 3-state fast driver output is turned under high impedance (HZ), to avoid being short circuited by the SSD circuit. Then, after a delay at least equivalent to the driver propagation time, the SSD is activated, including a serial resistance R SSD =75Ω, to turn off softly the device through discharging the gate input capacitance. The delay has been chosen

to be set at 50ns, for safety, fig. 10. The observed plateau in the gate-source voltage V GS under protection is due to the inductive effect through the kelvin inductance L SK (parasitic inductance of package and terminal pin on PCB in the presence of a di/dt < 0). Fig. 10: Logic Output signal, representing detection and protection time, C2M0080120 DUT, @400V-R g_ext =10.

Fig. 11 presents the waveforms of the device under short circuit, including protection. After the detection flag, the buffer is putted under HZ and the SSD is activated, the V GS starts to softly decrease for 530ns to reach V Drv-(-5V or -4V), limiting I ds to go higher. The V GS comparator goes back to low but the flag remains high until the PWM of the buffer goes low, allowing the system to detect at the next PMW pulse. The high speed detection allows safety turned off at a current level well below the saturation current of the channel. Indeed, the turn-off is done here at 175A while the saturation current is 290A. The SSD could therefore be accelerated in order to further reduce the energy stress at turn-off detection for the SiC MOSFET and to preserve its capability of endurance to the short-circuit cycles that the device could have to undergo throughout its life in accidental conditions of use. Fig. 12 presents the protection behavior under 2 pulse. The behavior is observed under 400V, for both operation, protection not included and included. In order to fully cover the robustness of this detection method, the circuit is validated under different bus voltages (50, 150, 250, 400V) which is lower than the typical 1.2kV (or 0.9kV) operational voltage, which makes the detection between NTO vs SC critical, Fig. 13.The one can see that with the final prototype the drain SC current is lower and stopped earlier than the previous prototype. The current is limited around 130A. is limited to 135A against nearly 290A with no protection. In such a short time and at such a low current value, the dynamic temperature of the chip remains lower than the melting temperature of the aluminum to metal of the die, which would allow a large number of repeated short circuits without ageing effect [START_REF] Fayyaz | Aging and failure mechanisms of SiC Power MOSFETs under repetitive short-circuit pulses of different duration[END_REF]. This last point is one of the future properties to be highlighted in the continuity of this work.

Fig. 1 .

 1 Two short-circuit electrical mechanism characteristics are especially interesting in this study located in part A and B. Part A of the figure might behave differently regarding the fault type (I or II) respectively Hard Switch Fault HSF or Fault Under Load FUL. This study is focused in short-circuit under HSF, whose fault dynamic is faster than FUL one. On other hand part B, is a behavior unique to SiC MOSFETs. Therefore, two main characteristics are studied in this paper related to the gate of the power transistor is the absence of miller plateau (part A), and the current injection in the gate also called dynamic gate leakage current (part B). In order to simplify the short-circuit under HSF, Fig.1.(b) presents an illustration of the fig.1 part-A.

Fig. 1 :

 1 Fig. 1: (a). MOSFET SiC experimental waveforms at two different drain biases under SC, (b). Simplified illustration between normal-turn on NTO and SC-HSF operations. DUT : 1.2kV, 80m@25°C, R gext = 47, T case = 25°C

Fig. 2 :

 2 Fig. 2: (a). Simplified detection method using the dynamic gate leakage current and protection circuit. (b). Oscilloscope waveforms of the DUT under protection

  (b) was proposed using on the low side a 1.2kV-80m, C2M0080120D SiC MOSFET transistor, and on the high side a SiC 600V SDB (Shottky barrier diode) with a load coreless inductor of 270 µH for the normal operation. A copper short-strap across the high side is used in order to get the HSF mode, fig. 4(b) is used. The gate driver used is a 3-state fast buffer, allows a high impedance configuration. The SMD circuit was implemented in the circuit and validated, Fig. 5. The power test bench used has for potential reference the kelvin source of the DUT. Low voltage signals (as V GS and V Drv ) are mostly measured with self-compensated voltage probes (300 V and 1 GHz bandwidth (BP) Tektronix TPP1000). High voltages (as V DS ) is typically measured with a probe (1000 V and 800 MHz BP Tektronix TPP0850). The current is measured through a 2GHz coaxial shunt.

Fig. 4 :

 4 Fig. 4: (a) 2D diagram of the gate charge method. (b) Circuit under test.

Fig. 5 :

 5 Fig.5: The double-pulse / short-circuit test bench and the SMD board based on the gate charge detection method: 6-layer 19mm² PCB as a safe gate-driver plugging-option.

Fig. 6 :Fig. 7 :

 67 Fig. 6: Oscilloscope waveforms under NTO (Initial prototype) @ 400V -R G-Ext =10  -V Supply-Drv =-5/20V C2M0080120D

Fig. 8 .

 8 Fig.8. Integrator output "Q G " under 400V NTO and SC operation for two power components, (a).C2M0080120D (b).C3M0065090D @R g_ext =10.

Fig. 9 .

 9 Fig.9. presents an image of V GS , which represents the voltage V GS adapted to 5V supply functions, therefore this voltage is compared with a reference level V Ref higher than the Miller plateau. This configuration set of the reference will activate the detection. Since the gate switch voltage between both transistor of different generations (G2 & G3) differ, respectively -5/20V -4/15V, the reference level should adapt for each generation (and not each transistor).

Fig. 11 :

 11 Fig.11: Oscilloscope waveforms, Circuit under SSD protection (SMD first prototype) C2M0080120 DUT @400V-R g_ext =10 -R SSD =75Ω.

Fig. 12 :Fig. 13 :

 1213 Fig. 12: DUT before and after protection, C2M0080120 DUT @ 400V R g_ext =10. (Initial prototype)

Table I : Gate charge Immunity Band and Ratio S Component

 I

This research work has been supported financially by the French Ministry in higher education, research and innovation.