G Bastos

L Sales

N Di

A Tayeb

J.-B Le Cam

N Di Cesare

Inverse-PageRank-Particle Swarm Optimisation for inverse identification of hyperelastic models: a feasibility study

Keywords: Inverse identification, hyperelasticity, FEMU, Particle Swarm Optimisation, Rubber, Heterogeneous test, Digital image correlation

In this study, the Finite Element Model Updating (FEMU) technique is used in order to identify hyperelastic parameters from only one heterogeneous test. A residual considering measured and identified stretches as well as the global reaction force of the specimen is built. The originality of this paper is to investigate the feasibility of the resolution of this minimisation problem by using the Inverse-PageRank-Particle Swarm Optimisation (PSO) for identifying hyperelastic parameters. For that purpose, the so-called PSO technique has been enriched with a PageRank algorithm to adapt iteratively the PSO parameters. As the paper examines whether Inverse-PageRank-PSO is adapted or not to the minimisation of the objective function in the present case, only two basic hyperelastic models have been considered.

Introduction

Hyperelasticity is widely used to predict the mechanical response of soft materials under large strains, for instance tissues [START_REF] Masson | Finite dynamic deformations of a hyperelastic, anisotropic, incompressible and prestressed tube. applications to in vivo arteries[END_REF][START_REF] Chagnon | Hyperelastic energy densities for soft biological tissues: A review[END_REF][START_REF] Chaimoon | An anisotropic hyperelastic model with an application to soft tissues[END_REF], laminates [START_REF] Li | Microscopic instabilities and elastic wave propagation in finitely deformed laminates with compressible hyperelastic phases[END_REF][START_REF] Aguiar | Analytical and numerical investigation of failure of ellipticity for a class of hyperelastic laminates[END_REF] and membranes [START_REF] Patil | A Finite inflation of an initially stretched hyperelastic circular membrane[END_REF]. Hyperelastic models are generally identified from several homogeneous tests, see [START_REF] Bazkiaei | A framework for model base hyper-elastic material simulation[END_REF] and [START_REF] Said | Efficiency of rubber-pad cushion in bending process of a thin aluminum sheet[END_REF] for instance, since the values of their constitutive parameters strongly depend on the strain state [START_REF] Treloar | Stress-strain data for vulcanised rubber under various types of deformation[END_REF]. Three homogeneous tests are classically considered, namely the uniaxial tension (UT), the pure shear (PS) and the equibiaxial tension (EQT), to completely describe the domain of possible loading paths [START_REF] Treloar | The Physics of Rubber Elasticity[END_REF][START_REF] Ward | An Introduction to the Mechanical Properties of Solid Polymers[END_REF]. A trade-off between the sets of values obtained with the different tests has therefore to be found to obtain parameters that can reasonably be considered as intrinsic to the mechanical behaviour of the material. As explained in [START_REF] Charlès | Inverse identification from heat source fields: a local approach applied to hyperelasticity[END_REF], such identification approach has many disadvantages: ⁻ each of the tests carried out are assumed to induce a homogeneous strain state, which is a strong assumption in case of PS and EQT tests, -several specimen geometries are required, -several testing devices are needed to apply these loading conditions, -dispersion obtained for each test requires testing several specimens for each loading condition, -the comparison between the constitutive parameters identified from different loadings is a matter of debate, -the elaboration process may differ from one specimen geometry to another one (typically compression moulding versus injection moulding).

This can have a significant effect on the values of the identified constitutive parameters and therefore on the predicted mechanical response. An alternative method consists in performing only one heterogeneous test as long as the strain/stress fields are sufficiently heterogeneous. This is typically the case when a multiaxial loading is applied to a 3-branch [START_REF] Guélon | A new characterization method for rubbers[END_REF] or a 4-branch (cruciform) [START_REF] Promma | Application of the virtual fields method to mechanical characterization of elastomeric materials[END_REF] specimen, which induces a large number of strain/stress states at the specimen's surface. This approach was also explored in Johlitz and Diebels [START_REF] Johlitz | Characterisation of a polymer using biaxial tension tests. Part I: Hyperelasticity[END_REF], Sasso et al. [START_REF] Sasso | Visco-hyper-pseudo-elastic characterization of a fluoro-silicone rubber[END_REF] and Seibert et al. [START_REF] Seibert | Biaxial testing of elastomers -experimental setup, measurement and experimental optimisation of specimen's shape[END_REF]. Such tests are all the more interesting that the full kinematic field can be measured and used to enrich the identification process. As no analytical relationship is available between measurements and parameters to be identified, an inverse identification procedure has to be used.

In the present paper, the identification is carried out with the Finite Element Model Updating (FEMU) technique. A residual considering measured and identified stretches as well as the global reaction force of the specimen is built. As this type of heterogeneous tests can be used to identify more complex behaviours, also governed by viscosity, permanent set and softening, to name a few, the number of parameters to be identified can increase drastically and methods such as the Particle Swarm Optimisation (PSO) [START_REF] Kennedy | Particle swarm optimization[END_REF][START_REF] Dreo | Métaheuristiques pour l'optimisation difficile[END_REF][START_REF] Clerc | The particle swarm -explosion, stability, and convergence in a multidimensional complex space[END_REF][START_REF] Sengupta | Particle swarm optimization: A survey of historical and recent developments with hybridization perspectives[END_REF][START_REF] Bonyadi | A theoretical guideline for designing an effective adaptive particle swarm[END_REF][START_REF] Pan | Hybrid particle swarm optimization with simulated annealing[END_REF][START_REF] Uriarte | A Characterization of recycled rubber using particle swarm optimization techniques[END_REF] can be envisaged for the resolution of this minimisation problem.

The main goal of this paper is to investigate the feasibility of using the Inverse-PageRank-Particle Swarm Optimisation (I-PR-PSO) and has therefore to be considered as a feasibility study. For that purpose, the so-called PSO technique has been enriched with a PageRank algorithm to adapt iteratively the PSO parameters to increase the convergence capabilities of the particles as the optimisation calculation is going on. As the paper examines whether Inverse-PageRank-PSO is adapted or not to the minimisation of the objective function in the present case, only two basic hyperelastic models have been considered.

The paper is organised as follows. First, the main features of the PSO and the PageRank algorithm used are presented. Secondly, a numerical study illustrates the relevance of the approach with a special emphasis on the influence of measurement noise on identified parameters. Third, the experimental test corresponding to an equibiaxial tensile test is described and the identification results are discussed. Concluding remarks close the paper.

Inverse identification method

In this section, the inverse identification of the constitutive parameters is presented. First, the finite element (FE) model is described. Secondly, the optimisation strategy carried out with a metaheuristic population based algorithm (the PSO algorithm) involving artificial intelligence (the PageRank algorithm) is precisely detailed. Third, the methodology is validated with a purely numerical approach. The section closes with a discussion on the effects of the measurement noise on the identified parameters.

FE model

The geometry chosen is presented in Fig. 1. It is a 105 mm long and 2 mm thick cruciform specimen. Such a geometry leads to a strong heterogeneity of the strain fields in terms of both the strain states and the distribution of the strain levels for a given strain state (see [START_REF] Promma | Application of the virtual fields method to mechanical characterization of elastomeric materials[END_REF] and [START_REF] Guélon | A new characterization method for rubbers[END_REF] for further details). Fig. 1 Specimen geometry (dimensions in mm) and FE mesh. It should be noted that with FEMU, a problem is often encountered, which is to determine the correct force distribution along the boundary of the specimen. Indeed, slippage in the grips inducing heterogeneities in the force distribution cannot be really measured. The problem is resolved by adding small cylinders at the end of the four branches of the specimen and by using suitable grips.

They prevent slipping into the machine's grips

The FE code used is ANSYS APDL. The finite element calculation is performed by assuming plane stress state and material incompressibility. For that purpose, the four-noded PLANE182 element is used. The mesh is made of 9600 nodes and 9353 elements. It is illustrated in Fig. 1.

The equibiaxial tensile loading is prescribed by applying the same maximum displacement of 70 mm to the four specimen branch ends. The time increment is chosen in such a way that the predicted and measured strain fields can be compared at the same prescribed displacements. Seven displacement levels were considered: 11, 21, 31, 41, 51, 61 and 70 mm. Two hyperelastic models were chosen. The first one is the Mooney model [START_REF] Mooney | A theory of large elastic deformation[END_REF] and is given by the following form of the strain energy density:

𝑊 ! = 𝐶 " (𝐼 " -3) + 𝐶 # (𝐼 # -3) (1)
In Eq. 1, 𝐼 " and 𝐼 # are respectively the first and second invariants of the right Cauchy-Green strain tensor C.

They are calculated by considering the material as incompressible (detF = 1). This model predicts quite well the nonlinear strain-stress relationship up to a moderate strain, i.e. until a stress hardening leading to an inflection in the mechanical response curve is observed at large strains.

In order to account for the stress hardening effect, the Yeoh model has also been considered [START_REF] Yeoh | Some forms of the strain energy function for rubber[END_REF]. The corresponding strain energy density is written as follows:

𝑊 $ = 𝐶 "% (𝐼 " -3) + 𝐶 #% (𝐼 " -3) # + 𝐶 &% (𝐼 " -3) & (2)
The values of the constitutive parameters to be identified were evolving at each iteration of the optimisation process. The value of the incompressibility parameter (𝐾 '" ) was set to 10 '( 𝑀𝑃𝑎 '" for all the FE calculations proceeded, which is low enough to consider the material incompressibility [START_REF] Ansys I | ANSYS Mechanical APDL Material reference[END_REF]. Note that only one and four (equally-spaced) displacement levels are sufficient for the identification of the Mooney and Yeoh models, respectively. This has been found from both experimental and generated data considering several combinations of the displacement levels considered in the identification.

Metaheuristic optimisation strategy

Definition of the objective function

The optimisation process aims at determining the constitutive parameters for the predicted data to fit the experimental ones. Experimental data considered here are the kinematic fields at different displacement levels in the specimen's branches, as well as the reaction force in the branches. Indeed, as the hyperelastic models chosen involve a linear relationship between the stress and the constitutive parameters, the kinematic fields can be the same for two different sets of constitutive parameters. This is the reason why the force also needs to be considered in the objective function calculation. The objective function to be minimised is thus defined as the squared relative difference between the experimental (exp) and the numerical (num) data: + $%& (() '+ ,-.

𝑚𝑖𝑛 ∑ ∑ 6
(()

+ $%& (() 8 # 9 , ) ,-" . /-" ( 3
)
where N is the number of nodes and k corresponds to the number of strain states (displacement levels applied) considered in the optimisation process. 𝜆 012 (/,,) and 𝜆 0/6 (/,,) are the maximum and the minimum in-plane principal stretches. 𝐹 (/) is the horizontal force measured at the machine horizontal grip. Only one strain field (𝑘 = 1) is sufficient to identify the two constitutive parameters of the Mooney model as soon as a sufficient distribution in the strain level and the strain state is induced. Since the stress hardening cannot be predicted by the Mooney model, the corresponding displacement fields were not used within the identification procedure. For the Yeoh model, all the strain fields are used, obtained at different displacement levels, to activate the three constitutive parameters and therefore to be able to take the inflection point into account, corresponding to a stress hardening in the stressstrain curve.

Classical Particle Swarm Optimisation algorithm

The classical Particle Swarm Optimisation (PSO) algorithm is considering particles (initially belonging to a bird flock or fish swarm, see for instance [START_REF] Kennedy | Particle swarm optimization[END_REF] for further details) that are each a potential solution to the optimisation problem to solve [START_REF] Schutte | Parallel global optimization with the particle swarm algorithm[END_REF][START_REF] Parpinelli | A comparison of swarm intelligence algorithms for structural engineering optimization[END_REF]. These particles are then "flying" through the solution domain, whatever its dimension, to converge together to the global minimum by smartly following each other. For that purpose, the particles positions, denoted 𝑿, and velocity2 , denoted 𝑽, have to be calculated for every particle at each iteration of the optimisation process. Their movement in the solution domain is obtained by making them follow the best one denoted 𝑮 789: and the best one from their neighbourhood denoted 𝑷 789: . The speed and position of particle m at iteration 𝑡 + 1 are given as: C 𝑽 0 :;" = 𝜔 × 𝑽 0 : + 𝑐 " × 𝑟𝑎𝑛𝑑 " × I𝑷 0,789:

:

-𝑿 0 : J + 𝑐 # × 𝑟𝑎𝑛𝑑 # × (𝑮 789: : -𝑿 0 : ) 𝑿 0 :;" = 𝑿 0 : + 𝑽 0 :;" , (4)
𝑷 0,789: is the best particle in particle 𝑚'𝑠 neighborhood, 𝑐 " and 𝑐 # are confident parameters that weight the importance of the neighbourhoods' memory and the global swarms' memory, respectively. The inertia weight 𝜔 is weighing the influence of the last iteration's speed on the calculation of the new one. 𝜔 is then considered as an inertia imposed to the particle as the calculation is going on. A proper value of 𝜔 can balance the particles' ability to (i) explore the solution domain when 𝜔 is large (around 1:2), and (ii) exploit interesting areas of the solution domain when 𝜔 is small (around 0:4).

The particle's speed and position have to be constrained for the calculation to converge. The extremal speed of the particles is then defined in [-𝑉 012 ; 𝑉 012 ] where 𝑉 012 is given as a function of the solution domain, such as 𝑉 012 = (𝑋 012 -𝑋 0/6 ) = 5 where 𝑋 012 and 𝑋 0/6 are the boundaries of the solution domain. The neighbourhood of every particle can be defined as a link to the other particles [START_REF] Kennedy | Population structure and particle swarm performance[END_REF][START_REF] Mendes | Watch the neighbor or how the swarm can learn from its environment[END_REF][START_REF] Mendes | The fully informed particle swarm: Simpler, maybe better[END_REF][START_REF] Arvind | Neighborhood re-structuring in particle swarm optimization[END_REF], i.e. the way the particles are influencing each other. Some static neighbourhoods have been proposed in the literature [START_REF] Kennedy | Population structure and particle swarm performance[END_REF], such as the Global Best (GBEST) topology, in which every particle is influenced by all the others, and the Local Best (LBEST) topology, in which every particle is linked and influenced by two of its peers. Moreover, some dynamic particles topologies have been proposed in the literature [START_REF] Akbari | A rank-based particle swarm optimization algorithm with dynamic adaptation[END_REF][START_REF] Janson | A hierarchical particle swarm optimizer[END_REF][START_REF] Pasupuleti | The gregarious particle swarm optimizer (G-PSO)[END_REF][START_REF] Suganthan | Particle Swarm Optimiser with neighbourhood operator[END_REF][START_REF] Hong | Analysis of parallel sub-swarm PSO with the same total particle numbers[END_REF], in which the way the particles are linked is evolving as the calculation is going on.

Inverse-PageRank-PSO

The Inverse-PageRank-PSO algorithm used in this work is based on both the PSO and the PageRank algorithm used by the world-famous search engine Google [START_REF] Brin | The anatomy of a large-scale hypertextual web search engine[END_REF][START_REF] Fercoq | Optimization of perron eigenvectors and applications: from web ranking to chronotherapeutics[END_REF][START_REF] Langville | Google's PageRank and beyond: The science of search engine rankings[END_REF]. In this strategy the population's topology is based on the GBEST one, but the links between all the particles are weighed regarding their fitness i.e. their instant capability to converge to the global minimum of the considered objective function. This means that the closest particles to the littlest minimum found so far have a greater influence on the others than the distant ones. As demonstrated in [START_REF] Cesare | A new hybrid PSO algorithm based on a stochastic markov chain model[END_REF], the particles topology can be seen as an oriented graph, in which the particles are the nodes, and the links between them represent the influence of the particles on each other. As these weights are not depending on the past, this oriented graph can be seen as a Markov chain, in which the transition probabilities are the probabilities for a particle to follow one or the other of its peers. To calculate the weighted influence of every particle on the others, an inverse version of the PageRank algorithm is used. The PageRank algorithm is used by the search engine Google to rank the webpages when a search is formulated: the more the links going in a webpage, the higher the rank of this webpage. The aim here is to do the opposite, which is to deduce the probability transitions of the Markov chain by considering the PageRank, i.e. the steady-state of the Markov chain, as known. Each particle performance fitness (𝑷 0 ) is compared to the best one 𝑮 789: and transformed into a transition probability in the PageRank vector π as follows:

𝝅 :1<=8: > (1, 𝑚) = U 𝑓𝑖𝑡𝑛𝑒𝑠𝑠(𝑮 789: ) × 100 𝑓𝑖𝑡𝑛𝑒𝑠𝑠(𝑮 789: ) -𝑓𝑖𝑡𝑛𝑒𝑠𝑠(𝑷 0 ) + 𝜖 U ∀ 𝑚 ∈ [1, 𝑝] ( 5
)
where p is the number of particles in the swarm and 𝑓𝑖𝑡𝑛𝑒𝑠𝑠(𝑿) is the objective function value calculated for particle 𝑿, previously defined in Eq. 3. Then, a pseudo-random process called inverse PageRank is used to calculate the transition probability matrix 𝑪 corresponding to the previous 𝝅 vector. In the final converged connectivity matrix, the elements of row 𝑟 are relative to the links going out of node 𝑟, while the elements of column 𝑐 are relative to the links going into node 𝑐. The components of matrix 𝑪 are then considered as the influence of all the particles on the others, that is the probability for each of them to follow every other one. In this way, the best particles -that is the ones obtaining the littlest values of the objective function -obtain the greatest values in 𝝅 and in 𝑪 as well. This inverse PageRank process is launched at every iteration of the optimisation process, for the particles topology to be smartly evolving and adapting to the swarm's performance as the calculation is going on. The neighbourhood part of Eq. 4 is then modified, weighing the influence of all the particles by using the components of 𝑪, as follows: ] 𝑽 / :;" = 𝜔 × 𝑽 / : + 𝑐 " × 𝑟𝑎𝑛𝑑 " × I𝑷 /,789:

:

-𝑿 / : J + 𝑐 # × 𝑟𝑎𝑛𝑑 # × ^𝑪/, 6 ,-"
× _𝑷 /,789: :;" -𝑿 / : X/ :;" = 𝑿 / : + 𝑽 / :;"

. ( 6
)
In this work, every particle will be a vector containing as much variables as the model contains parameters to be identified, and the objective function defined in Eq. 3 will be calculated for each of them. In the present study, the authors have made their own algorithm, coded with Matlab software [START_REF]MATLAB: version 9.4 (R2018a[END_REF].

Numerical validation

In this section, the inverse identification methodology proposed is applied to data obtained from FE simulations of the biaxial test previously defined in 2.1 with a given set of constitutive parameters for both the Mooney and Yeoh models. The constitutive parameters used are reported in Table 1. The research domains for each constitutive parameter to be retrieved are reported in Table 2. Note that 𝐶 #% is imposed to be negative for the Yeoh model, in order to fairly predict the shear modulus for all ranges of strain as explained in [START_REF] Yeoh | Some forms of the strain energy function for rubber[END_REF]. The stability criterion defined by Drucker [14] is used within the FE code ANSYS APDL in order to ensure the behaviour law stability. If the energy density is unstable, then the FE simulation will not converge. In this case, a large value of the objective function is affected (10 ? typically). Therefore, this research domain area will not be considered as a promising zone to be exploited by the PSO particles in the next iterations of the optimisation process. The objective here is to retrieve the values of the constitutive parameters for the two hyperelastic models, which predict the mechanical response of the cruciform specimen. As above mentioned, the kinematic fields used in the identification procedure are the maximum and minimum in-plane principal stretches. The distribution of the maximal principal stretch for a displacement of 70 mm is reported in Fig. 2a. It lies between 1.8 for equibiaxial tension state and 3.2 for uniaxial tension. The heterogeneity is enough for the identification to be carried out. This is in a good agreement with the conclusion previously drawn for such a test in [START_REF] Promma | Application of the virtual fields method to mechanical characterization of elastomeric materials[END_REF].

In order to account for experimental measurements, the stability of the identification procedure to measurement uncertainties is examined by adding three amplitude noise levels 𝑁𝑜𝑖𝑠𝑒 @0A of 0.01, 0.05 and 0.1 to the numerical principle stretches fields, which are dimensionless quantities. This noise was generated by using Matlab function 𝑟𝑎𝑛𝑑𝑛. It should be noted that noise in real camera sensors is signal-dependent (or heteroscedastic) and that a uniform (or homoscedastic) noise is considered here instead for the sake of simplicity, as in similar studies dealing with identification [START_REF] Promma | Application of the virtual fields method to mechanical characterization of elastomeric materials[END_REF].

The noised numerical data are reported in Figs 2b, 2c and 2d for 0.01, 0.05 and 0.1 noise amplitudes, respectively. It should be noted that the maximum error caused by the noise introduced to the numerical data can exceed 11% from its original value. This error is greater than the correlation error calculated in Section 3.1. The identification methodology was applied to identify the constitutive parameters for both Mooney and Yeoh models by considering the three noise levels for each of them. Given that the Mooney model is not able to describe the hardening phenomenon, a maximum displacement of 51 𝑚𝑚 is considered for the identification procedure. This value corresponds to the appearance of the inflection point in the experimental force-displacement curve. Only one strain state (𝑘 = 1 in Eq. 3) corresponding to the maximum displacement applied is considered for the identification. It should be noted that no post-processing is needed for the numerical data since the kinematic fields were obtained from the same Ansys FE model. The identified parameters from reference and noised numerical data are reported in Table 3. The first parameter 𝐶 " is fairly identified since the relative error is inferior to 3 % for the reference and noised data. However, the second parameter was not fairly retrieved from the noised data. This is due to its greater sensitivity to noise as it is not sufficiently activated at this displacement level, which was also seen in [START_REF] Promma | Application of the virtual fields method to mechanical characterization of elastomeric materials[END_REF]. Moreover, following the formulation of the Mooney model, this error is not as significant as it seems because the first parameter has a greater contribution to the strain energy density of Eq. 1 than the second one. This is due to its association with the second invariant 𝐼 # , which needs a larger EQT zone to enhance its contribution to the global response of the material. The overall veracity of the identification procedure is seen in the values of the force error and the objective function values. Fig. 3 illustrates the accuracy of the identified parameters from the force response of the cruciform specimen during the biaxial test.

Identification results for the Yeoh model

For the Yeoh model, four strain states corresponding to a prescribed displacement of 11 mm, 31 mm, 51 mm and 70 mm (𝑘 = 4 in Eq. 3) were used in the identification of the constitutive parameters. These several strain states are necessary to account for the inflection point in the force-displacement curve 3 . The identified parameters from the reference and noised numerical data are reported in Table 4. These parameters show the capacity of the identification procedure to retrieve the constitutive parameters even from noised data. This is also shown in terms of the value of the final force error, which lies between 10 '? and 10 '# and the final objective function values, which lies from 10 'C to 10 '" for all the analysed data. Furthermore, the evolution of the force during the biaxial test is well predicted, especially the inflection point in the force-displacement curve 2 . This is shown in Fig. 4. As a conclusion, for the two models under consideration, the identification procedure showed its relevancy even when the kinematic fields are strongly noised. This methodology is applied to experimental data in the next section.

Application to experimental data

In the previous part, the methodology has been successfully applied to numerical data. It was shown how much noise effect influences the identification methodology. The result obtained clearly shows that the methodology is well suitable for identifying constitutive parameters from experimental, i.e. noisy, kinematic fields, which is the aim of the present section. Fig. 5 Experimental setup Fig. 5 presents an overview of the experimental setup. It consists in a home-made biaxial testing machine and an optical camera. The machine is composed of four independent electrical actuators, controlled by an inhouse LabVIEW program. Two load cells, whose capacity is equal to 1000 𝑁, measure the force in the two perpendicular directions. In the present study, four equibiaxial load-unload cycles were prescribed in terms of grip displacement. The displacement and the loading rate were set at 70 mm and 150 mm/min respectively for each of them.

Full kinematic field measurement

Images of the specimen surface were stored at increasing displacement of the grips at a frequency equal to 5 Hz, which is the lowest frequency allowed by the current experimental setup, with an IDS camera equipped with a 55 mm telecentric objective. The charge-coupled device (CCD) of the camera has 1.920 × 1.200 joined pixels. The displacement field at the specimen surface was determined by using the DIC technique. It consists in correlating the grey levels between two different images of a given zone, each image corresponding to a different strain state. DIC hardware and analysis parameters are reported in Tables 5 and6, respectively, such as recommended in the guideline [START_REF] Jones | A good practices guide for digital image correlation[END_REF]. In order to improve the image contrast, a white paint was sprayed on the specimen surface before the test, which led to a random speckle with a feature pattern of about 6 pixels approximatively. A uniform cold lighting at the specimen surface was ensured by a home-made LED lighting system. DIC was performed using SeptD [START_REF] Vacher | Bidimensional strain measurement using digital images[END_REF]. The parameters used for the DIC are reported in Table 6. Namely, the smallest distance between two independent points, also referred to as the step size, was equal to 4 pixels (here the size of the zones of interest (ZOIs)) corresponding to 289.6 µ𝑚. The Region Of Interest (ROI) used to compute the displacement field with the DIC technique is depicted in Fig. 6. It corresponds to an area of 384 × 680 pixels. The displacements and strains noisefloors reported in Table 6 were obtained by performing the DIC using static images i.e. before application of the mechanical loading. Note that these values present the standard deviation of the displacement and strain fields over the ROI. These are the typical values for a good DIC performance, see [START_REF] Jones | A good practices guide for digital image correlation[END_REF] for further details. Given that the applied displacement and strains are large, the strain formulation used in this work is the logarithmic one, which is recommended for the finite strain framework [START_REF] Hencky | Uber die form des elastizitatsgesetzes bei ideal elastischen stoffen[END_REF].

It should finally be noted that the strain measurement used here is the maximum and minimum in-plane principal stretches of the transformation gradient tensor F, denoted 𝜆 012 and 𝜆 0/6 respectively.

Camera Image Resolution

Lens

Kinematic fields processing

As the data point coordinates of experimental and numerical fields did not coincide, the experimental kinematic fields were fitted locally by polynomial functions. When the stretch fields obtained by the DIC technique exhibited uncorrelated zones, the values in these zones were interpolated from surrounding ZOIs and then fitted by the polynomial functions. This post-processing step reduces the noise effect, without altering the signal, as no strong strain gradients are induced under such loading conditions. The polynomial fit was carried out at the local scale, i.e. for every point in the ROI, by considering the 6 closest points to calculate the 6 coefficients of the secondorder polynomial used. The experimental data were replaced by the polynomial functions in the objective function calculation, for the errors to be calculated at exactly the same coordinates.

Results and discussion

In this section, the experimental data are first given; the mechanical response, the kinematic field measurement and post-processing and the results of the identification methodology proposed are highlighted.

Mechanical response

Fig. 7 gives the experimental mechanical response obtained in terms of the force in the horizontal direction versus the displacement applied for each specimen branch. The force used in the identification procedure is the one obtained during the fourth load (plotted in red colour Fig. 7). This choice was made to identify the stabilised behaviour of the specimen, since the material accommodates during the first cycles [START_REF] Mullins | Effect of stretching on the properties of rubber[END_REF]. It should be noted that the accommodation level is different from one point to another in the specimen. Therefore, the identification is done in the case of an accommodated structure. As above-mentioned, the curve exhibits a strong non-linearity that is amplified at the largest displacement applied due to the hardening effect. An inflection point is indeed observed at a prescribed displacement of 50 mm. Fig. 10 provides the same fields as in Fig. 9, but for a displacement equal to 70 mm, which was the maximum displacement applied during the test. In order to evaluate the relevancy of the experimental data postprocessing, the differences between the kinematic fields obtained with the DIC technique and the polynomialbased functions are calculated for every point in the ROI. The calculations were performed for displacements equal to 51 𝑚𝑚 and 70 𝑚𝑚. For a displacement equal to 51 𝑚𝑚, 98.51 % of the points in the ROI exhibit an error inferior to 1% for the maximal stretch and 97.52 % of the points in the ROI exhibit an error inferior to 1 % for the minimal stretch. The results are better for the largest displacement (70 𝑚𝑚): all the points exhibit an error inferior to 1 %. 8. The inertia of the particles has been chosen to be linearly decreasing as the calculation is going on, to encourage the particles to explore the design domain at the beginning of the calculation, and to focus on interesting areas of the domain at the end of the calculation. Fig. 11 gives the convergence curves of the optimisation calculations for the two models. The values of the optimised objective function and design variables are given in Table 9. The errors between the experimental and numerical values of the principal stretches for every point of the ROI are depicted in Fig. 12 and 13 for the Mooney and the Yeoh model, respectively 4 . We recall here that several strain states were necessary to identify the parameters for the Yeoh model. Contrarily to what was done in the numerical validation section, the strain states used were not homogeneously distributed with respect to the prescribed displacement. As explained previously, the mechanical response has not exactly the same shape as the one given by the Yeoh model. It should be noted that several possibilities have been tested beforehand. 8 Optimisation parameters For the Mooney model, the maximal stretch field was fairly retrieved since 92.25 % of the points in the ROI exhibit an error inferior to 3%, whereas for the minimal stretch, over 90 % of the points have an error inferior to 5 %. But this deviation has no significant effect on the global response as it is seen in Fig. 14. For the Yeoh model, both maximal and minimal stretches were fairly retrieved since the error was inferior to 5 % for the whole ROI. The validation of the optimised values of the constitutive parameters is performed by comparing the experimental and numerical forces -displacement curves. The mechanical responses measured and predicted are given in Fig. 14. For the Mooney model, a good agreement is found for a maximum displacement equal to 51 mm. For the Yeoh model, the result was satisfactory for the whole displacement range investigated here, i.e. including the inflection in the curve, meaning that the hardening was predicted quite well. This result shows that only one deformed state is sufficient to identify the Mooney constitutive parameters and only 4 deformed states for the Yeoh model. It should be noted that several combinations were tested for the identification of the Yeoh model.

Conclusion

In this study, an identification methodology based on the coupling between PSO and FE simulations is proposed for determining hyperelastic model parameters from only one heterogeneous tests. A PageRank algorithm is used to adapt iteratively the PSO parameters. The mechanical test considered consists in stretching a cross-shaped specimen, which induces simultaneously a wide range of strain states, especially the uniaxial tension, the pure shear, the equibiaxial tension and many intermediary strain states. The data used for the identification is the kinematic fields and the reaction force in one of the specimen's branches. The identification methodology is first benchmarked with a purely numerical approach. Then, it has been successfully applied to experimental data. The kinematic field was characterised with the DIC technique. The identified parameters enabled us to satisfactorily model the material behaviour. This study provides numerous perspectives as the number of unknown constitutive parameters can increase significantly if other phenomena have to be taken into account in the mechanical response of rubbers: hysteresis, permanent set, strain-induced anisotropy, strain-induced crystallisation, accommodation (the Mullins effect [START_REF] Mullins | Effect of stretching on the properties of rubber[END_REF]), to name a few.

  (a) Numerical data (b) 𝑁𝑜𝑖𝑠𝑒 @0A = 0.01 (c) 𝑁𝑜𝑖𝑠𝑒 @0A = 0.05 (d) 𝑁𝑜𝑖𝑠𝑒 @0A = 0.1 Fig. 2 Maximal principal stretch fields for a displacement equal to 70 𝑚𝑚: numerical and noised data Identification results for the Mooney model

Fig 3

 3 Fig 3 Force response versus displacement for parameters ofTable 3 (Mooney model)

Fig. 4

 4 Fig. 4 Force response versus displacement for parameters ofTable 4 (Yeoh model)

ApertureField

  Fig. 6 Region of Interest (ROI)

Figure 7 :

 7 Figure 7: Mechanical response in terms of the horizontal force versus the displacement applied 1.2. Kinematic field measurement and post-processing Fig.8illustrates the distribution of the images considered to determine the strain fields. 8 images are used to determine 7 strain fields at increasing displacement. The first image is the reference one, which corresponds to the undeformed state. The DIC process was carried out iteratively. The correlation is first done between two successive images and the strain field is then built with respect to the undeformed state. In the case of the Mooney model, the inflection in the mechanical response curve is not predicted. Only one image before the inflection point (Image #6, corresponding to a displacement equal to 51 𝑚𝑚 applied to each branch) was required for the identification of the two constitutive parameters. In the case of the Yeoh model, several images are considered to account for the hardening up to a displacement equal to 70 𝑚𝑚 at the branch end. Fig.9gives the fields of maximum and minimum (Figs. 9a and 9c respectively) principal in-plane stretches for a displacement of 51 mm. In these fields, the values in uncorrelated zones (white ZOIs in the measured fields) are first interpolated from surrounding ZOIs and then fitted by the polynomial functions as previously described in Section 3.2. The obtained fields are given in Figs.9b and 9drespectively.

Fig 8 8

 8 Fig 8 8 images are considered to determine the 7 strain fields during the test

Fig. 11

 11 Fig. 11 Convergence curves of the optimisation processes Mooney Yeoh Iterations 50 30 Particles 15 20 𝜔 1.2 -0.4 1.2 -0.4Table 8 Optimisation parameters

Fig. 12

 12 Fig. 12 Mooney model: final results

Fig 14

 14 Fig 14 Comparison between optimised numerical results and experimental data

Table 2

 2 Research domains for the constitutive parameters

		Mooney	Yeoh
	0.4 0.04 ---Table 1 Reference Constitutive parameters -𝐶 " MPa -𝐶 # MPa 0.5 𝐶 "% MPa -0.02 𝐶 #% MPa 0.005 𝐶 &% MPa
		Min	Max
	𝐶 " MPa 𝐶 # MPa 𝐶 "% MPa 𝐶 #% MPa 𝐶 &% MPa	0.1 0.01 0.1 -0.04 0.001	1 0.1 1 0 0.01

Table 3 (

 3 Mooney model)

	𝑁𝑜𝑖𝑠𝑒 @0A 𝐶 " MPa Error 𝐶 # MPa Error Force error (N)	0 0.4014 0.35 % 0.0377 -5.7 % 3.4 10	0.01 0.3902 -2.4 % 0.0563 40.75 %	0.05 0.3892 -2.7 % 0.0579 44.75 %	0.1 0.3944 -1.4 % 0.0491 22.75 %

'?

5.

[START_REF] Sasso | Visco-hyper-pseudo-elastic characterization of a fluoro-silicone rubber[END_REF]

10 '( 8.4 10 '( 1.4 10 '( Objective function 1.16 10 'B 2.9 10 'B 7.1 10 'B 2 10 '"%

Table 3

 3 Identification results for Mooney model

Table 4 (

 4 Yeoh model)

	𝑁𝑜𝑖𝑠𝑒 @0A 𝐶 "% MPa Error 𝐶 #% MPa Error 𝐶 &% MPa Error Force error (N)	0 0.502 0.4 % -0.0201 -0.5 % 0.00496 -0.8 % 2.6 10	0.01 0.4997 -0.06 % -0.0195 2.5 % 0.00494 -1.2 %	0.05 0.5004 0.08 % -0.0198 1 % 0.00497 -0.6 %	0.1 0.5034 0.68 % -0.0203 -1.5 % 0.0048 -4 %

'& 9.2 10 '? 1.14 10 '& 1 10 '# Objective function 3.37 10 '( 5.4 10 'C 6 10 'C 4.67 10 '"

Table 4

 4 Identification results for Yeoh model

Table 6

 6 DIC Analysis parameters

Table 7 .

 7

				Mooney
		Research domain				Speed limits
		min	max			min	max
	𝐶 " MPa 𝐶 # MPa	0 0	1 0.1	Yeoh	𝐶 " MPa 𝐶 # MPa	-0.5 -0.05	0.5 0.05
		Research domain				Speed limits
		min	max			min	max
	𝐶 "% MPa 𝐶 #% MPa 𝐶 &% MPa	0.1 -0.035 0.001	1 0 0.01		𝐶 "% MPa 𝐶 #% MPa 𝐶 &% MPa	-0.105 -0.003 -0.00025	0.105 0.003 0.00025

Table 7

 7 Research domainsOptimisation algorithm parameters are given in Table

Table 9

 9 Results of the optimisation calculation

	Mooney	Yeoh

called "the velocity" in the literature, this parameter actually represents the particles displacement in the solution domain

Note that in this case, the mechanical behavior is fully described by the Yeoh model, which is not necessarily the case when dealing with experimental force-displacement curve. Thus, other intermediate strain states, i.e. strain states that would not be homogeneously distributed during the test, can be used in the identification procedure.

As the CPU time needed to calculate the objective function strongly depends on the machine used, the CPU time needed to converge is not a reliable indicator to evaluate the performances of the proposed methodology. Nevertheless, a more relevant information could be the number of calls to the OF calculation. It can be seen in Fig.11that the calculation can be considered as converged after 15 iterations, that is after 225 and 300 calls for the objective function calculation, for the Mooney and Yeoh models, respectively.

Acknowledgements

Cooper Standard is acknowledged for providing the specimens. The authors thank Dr Eric Robin for his help on image processing. Authors also thank Dr Mathieu Miroir, Mr Vincent Burgaud and Mr Mickäel Le Fur for having designed the biaxial tensile machine.

Declarations

Conflict of interest: The authors certify that they have NO affiliations with or involvement in any organisation or entity with any financial interest (such as honoraria; educational grants; participation in speakers' bureaus; membership, employment, consultancies, stock ownership, or other equity interest; and expert testimony or patentlicensing arrangements), or non-financial interest (such as personal or professional relationships, affiliations, knowledge or beliefs) in the subject matter or materials discussed in this manuscript.