Inverse identification of hyperelastic parameters by metaheuristic optimization algorithm

: In the present study, a numerical method based on a metaheuristic parametric algorithm has been developed to identify the constitutive parameters of hyperelastic models, by using FE simulations and full kinematic ﬁ eld measurements. The full kinematic ﬁ eld was measured at the surface of a cruciform specimen submitted to equibiaxial tension. The test was simulated by using the ﬁ nite element method (FEM). The constitutive parameters used in the numerical model were modi ﬁ ed through the optimization process, for the predicted kinematic ﬁ eld to ﬁ t with the experimental one. The cost function was formulated as the minimization of the difference between these two kinematic ﬁ elds. The optimization algorithm is an adaptation of the Particle Swarm Optimization algorithm, based on the PageRank algorithm used by the famous search engine Google.


INTRODUCTION
The constitutive parameters of hyperelastic models are generally identified from three homogeneous tests, basically the uniaxial tension, the pure shear and the equibiaxial tension.From about 10 years, an alternative methodology has been developed (Guélon et al. 2009, Promma et al. 2009, Johlitz et al. 2011, Seibert et al. 2014), and consists in performing only one hete-rogeneous test as long as the field is sufficiently hete-rogeneous.This is typically the case when a multiaxial loading is applied to a 3 branch (Guélon et al. 2009) or a 4-branch (Promma et al. 2009) cruciform specimen, which induces a large number of strain states at the specimen surface.The Digital Image Correlation (DIC) technique is generally used to characterize the full kinematic field at the specimen surface.The induced heterogeneity is analysed through the distribution of the biaxiality ratio and the maximal eigen value of the strain.Thus, a large number of experimental data is provided for identifying the constitutive parameters of the considered model.
Several methods have been recently developed to identify parameters from experimental full field measurements, typically the finite element updating method (FEMU), the constitutive equation gap (CEGM), the virtual fields method (VFM), the equilibrium gap method (EGM) and the reciprocity gap method (RGM).These methods are fully reviewed in (Avril et al. 2008).
In the present study a new methodology is proposed in order to minimize the cost function in the FEMU approach.The optimization algorithm used is based on the Particle Swarm Optimization (PSO) algorithm and the artificially smart PageRank algorithm used by the famous search engine Google.This algorithm enables us to minimize the full kinematic field differences by modifying the constitutive parameters, while minimizing the CPU calculation time.Even though the final objective is the identification of complex constitutive models, i.e. a large number of constitutive parameters, the two-parameter Mooney's model (Mooney 1940) is presented in this paper to illustrate the methodology.

EXPERIMENTAL SETUP
The material considered here is a carbon black filled natural rubber.The specimen geometry is given in Figure 1.It is a 2 mm thick cruciform specimen.Figure 2 presents an overview of the experimental setup composed of a home-made biaxial testing machine and an optical camera.The machine is composed of four independent RCP4-RA6C-I-56P-4-300-P3-M (IAI) electrical actuators.They were driven by a PCON-CA-56P-I-PLP-2-0 controller and four PCON-CA (IAI) position controllers.The actuators were controlled by an in-house LabVIEW program.Two load cells, whose capacity is equal to 1094 N, store the force variation in the two perpendicular directions.In the present study, one equibiaxial loading was carried out in such a way that the specimen'scentre was motionless for the displacement measurement to be easier.The displacement and loading rate were set at 70 mm and 150 mm/min respectively for the four independent actuators.
Images of the specimen surface at increasing stretches were stored at a frequency equal to 5 Hz with a IDS camera equipped with a 55 mm telecentric objective.The charge-coupled device (CCD) of the camera has 1920 × 1200 joined pixels.The Digital Image Correlation (DIC) technique is used to determine the displacement field at the specimen surface.The software used for the correlation process was SeptD (Vacher et al. 1999), and a uniform cold lighting was ensured by a home-made LED lamp.The spatial resolution, defined as the smallest distance between two independent points was equal to 4 pixels corresponding to 0.343mm.The Region Of Interest (ROI) and the grid used for the digital correlation is represented in Figure 3.

NUMERICAL MODEL
A finite element calculation is performed by assuming plane stress state and material incompressibility.The four-node PLANE182 ANSYS element is used.The mesh is made of 9600 nodes and 9353 elements.A biaxial traction load is obtained by prescribing the same displacement of 70 mm on the four branches of the specimen.The two-parameters hyperelastic Mooney model is used for the calculation.The values of the constitutive parameters are changing at each iteration of the optimization process, as described in the next section.
As the spatial resolution between the numerical and the measured kinematic fields was different, the experimental kinematic field was fitted by a polynomialbased function.In this way, the numerical kinematic field was compared, for each node, with the experimental field at the same position in the specimen.With the fitting method applied, the difference between the experimental field and the polynomial-based function was less than 0.2 mm for every point considered.

METAHEURISTIC OPTIMIZATION STRATEGY
The aim of the optimization process was to find the constitutive parameters for the numerical kinematic field to fit the experimental one.The cost function f was the squared difference between the experimental and the numerical fields, considering the force too, as follows: where N is the number of nodes in the numerical ROI, U x,exp is the polynomial-based experimental displacement, U x,num is the numerical horizontal displacement, F exp is the experimental horizontal force, and F num is the numerical horizontal force.
The optimization algorithm used is an adaptation of the classical Particle Swarm Optimization algorithm.In this version, all the particles are influenced by all the others, by considering this influence to be adapted as a function of the respective performance of the particles.The population of PSO particles is then considered as a Markov chain, in which the particles are the nodes, and the links between them are considered as the transition probabilities.For each particle, the Page-Rank valuethat is the steady state of the considered Markov chainis given by equation ( 2).In this way, the PageRank value of each particle is deduced from its performance compared to the best one.
It is possible to deduce the transition connectivity matrix C giving the influence of all the particles on all the others by using a pseudo-random process.The classical equations of PSO are modified, weighing the influence of all the particles by using the components of C, as follow: where V tþ1 i is the speed of the i th particle at iteration t+1, c 1 and c 2 are confident parameters, ω is the inertia weight, V tþ1 i is the position of particle i at iteration t+1, r 1 and r 2 are random numbers given in [0,1], V tþ1 j;best is the personal best position of particle i at iteration t+1, and C is the transition connectivity matrix of the considered Markov chain.This Inverse-Page-Rank-PSO algorithm is fully described in (Di Cesare et al. 2015).

RESULTS AND DISCUSSION
As the particles are initially randomly defined, the optimization was launched 10 times, to compare the obtained solutions and be sure that the global minimum of the cost function was reached.The convergence curves of the 10 launched optimization calculations are represented in Figure 4.The obtained values of the cost function and constitutive parameters are given in Table 1.
The validation of the optimized results is checked by comparing the displacements and efforts measured in the sample with the numerical optimized one.In the final numerical model, the values of the constitutive parameters have been set to the mean of the obtained optimized values found in the 10 different calculations launched.Figure 5 shows the difference between the experimental polynomial-based This difference is presented in a quantitative way in Figure 6 showing the difference between the two fields for every point of the ROI.One can note that the difference is always less than 6.5 % of the experimental kinematic field.For the force, the experimental value was 176.02 N, while the numerical value obtained with the optimized values of the constitutive parameters was 176.37 N, which leads to a difference up to 0.4 %.

CONCLUSIONS
This work is proposing a new numerical method for inverse identification of hyperelastic parameters.It is based on the use of a PSO-based parametric optimization algorithm.Experimental and numerical kinematic fields are compared to finally be fitted through the optimization process, while the constitutive parameters are smartly modified.This process is able to find the constitutive parameters reproducing the mechanical response of the specimen and the

Figure 4 .
Figure 4. Convergence curve of the 10 optimization calculations launched

Figure 5 .
Figure 5.Comparison between the two kinematic fields after the optimization process

Table 1 .
Obtained results