Numerical identification of the elasticity tensor of heterogeneous materials made of Silicon Carbide and Titanium by the Finite Element Model Updating (FEMU)

. This work is subjected to the development of a method to identify the elasticity tensor of homogeneous and heterogeneous materials. The materials are created in the form of checkerboards. We solved the direct problem to obtain the strain field using the finite element method, after obtaining this strain field, we created synthetic experimental displacement data by simulation. A re-calibration of the created experimental and simulated data is done based on the principle of the finite element model updating (FEMU), used in almost all domains, in the inverse problem. The minimization of the cost function obtained by FEMU is done by Levenberg–Marquardt algorithm which is very fast and elegant algorithm. A complete study has been done by studying the sensitivity of the identified values with respect to the refinement of the mesh and with respect to the level of disturbance


Introduction
Today we are concerned with the performance of composite materials and the determination of their elasticity properties which are the coefficients of the elasticity tensor. From a mechanical point of view, industrialists in the manufacture of machines such as automobiles, airplanes, helicopters and in the manufacture of prostheses, implants in the field of medicine are concerned with improving of the rigidity of materials. The improvement of the stiffness of a material cannot be done without its mechanical characterization. Indeed each material is governed by constitutive material parameters that compose the elastic constitutive law. Determining the coefficients of the elastic constitutive law of a material has been made possible by the methods of non-contact field measurements such as Digital Image Correlation (DIC) [1,2] which is the one most used in an inverse process to identify constitutive material parameters.
The principal identification methods known today are: the virtual fields method developed by Grédiac et al. [3] or the equilibrium equation gap method by Claire et al. [4], the constitutive equation gap method by Ladevèze et * e-mail: arnaud.germaneau@univ-poitiers.fr al. [5], the finite element model updating by Kavanagh et al. [6] and the reciprocity gap method by Ikehata [7]. An overview of these methods and their principles have been done by Avril et al. [8]. In [9] the authors make the link between all the methods and compare them to one another.
These methods were most often used in cases of global identification. This type of identification considers that the material is homogeneous which means here that the specimen is made of the same material. With this in mind, Leclerc et al. [10] have used the Finite Element Model Updating (FEMU) which have been integrated to a DIC to identify parameters of a material in a crucifix form. The advantage of this integration was to reduce the effects of the noise. Pottier et al. [11] modified the method and then used it to identify thermoplastic properties of a material. Maček et al. [12] recently used the FEMU considering different geometries and concluded on which test is appropriate to identify a specific property of an orthotropic material. The rest of the methods has been also widely used in the same assumption [13][14][15]. But in this last decade, researchers have been interested in a local identification. For this consideration the material is considered as being heterogeneous or made of at least two materials, which makes even more account composite materials. It is in this perspective Florentin et al. [16] made an identification on a domain considered to be filled with heterogeneous and isotropic fields. The specimen used was a planar specimen of size 10mm × 10mm consisting of 10 × 10 subdomains in one hand and 20 × 20 subdomains in other hand. The identification was done by the constitutive equation gap method (CEGM). A comparison was made with the equilibrium gap method. The relative errors obtained were good and the authors concluded that the CEGM is better than the equilibrium gap method. The sampling of a cruciform specimen was sampled by zones of orthotropic material and The equilibrium gap method was used for identification per orthotropic element [17]. A modification of the virtual fields method by Fourier series has been applied to the case of local identification while providing solutions to the unknown boundary conditions problem [18]. The identification of the first Lamé parameter was done locally in 2D and 3D by Guchhait et al. [19]. Madani et al. [20] used the CEGM to make a local identification and studied the sensitivity of the identified parameters with respect to the initial input values. The speed of convergence was not affected too much. The sensitivity of the identified values with respect to the refinement of the mesh was also studied; the errors decrease considerably but, as expected, the computation time becomes longer and longer. Finally, the authors have studied the sensitivity of the parameters to measurement noise. They concluded, as other authors, that the accuracies become worse as the level of disturbance increases.The FEMU method which is based on the comparison of the experimental and calculated data has been used by Franquet et al. [21] to identify the mechanical properties of an atherosclerotic plaque which is an heterogeneous domain. The authors of [22,23] proposed two types of objective functions, one geometric and the other based on kinematic data. Applied to the identification of heterogeneous material parameters, the geometric cost function is less sensitive to measurement errors. Gras et al. [24] identified the properties of heterogeneous materials and studied the effect of the noise on the identified parameters and to reduce the errors, they have used a regularization method. In 2019, this approach has been used by Pétureau et al. [25] to identify the local parameters of an isotropic heterogeneous material behavior law. Two cost functions have been formulated, one based on experimental and numerical displacements and the other one based on strain fields. The authors concluded that the one based on the strain fields gives better results than the other one. A genetic algorithm was used to identify a large number of parameters with good accuracy. Volume identification of elastic parameters has been performed by the same authors with the same methods and considerations [26]. Ogierman et al. [27] used the finite element model updating method for identification on a two-phase composite material. Each phase is characterized by an isotropic and homogeneous material. The sensitivity of the identified parameters to the refinement of the mesh was studied and they conclude that the errors decrease when the mesh becomes finer.
The objective of this work is to develop a method to identify linear elasticity parameters of heterogeneous materials. The identification method which will be used is the Finite Element Model Updating (FEMU), because of the fact that we can use it in almost all the domains [28]. We have taken up the main steps of this approach such as the formulation of the objective function and its minimisation. For this we use the Levenberg-Marquardt minimization algorithm which is a very fast [29] gradient based algorithm. The identification analysis is conducted based on a comparison of the identified parameters stability with respect to the choice of the initial parameter set, to the refinement of the mesh and on the robustness of the algorithm to handle experimental noisy data on heterogeneous materials.

Direct problem
The physical problem consists here to solve the linear elastic problem in the case of isotropic heterogeneous material. We consider a domain Ω ⊂ R 2 and assume that Ω is bounded by a piecewise smooth boundary ∂Ω such that ∂Ω = ∂ T Ω ∪ ∂ u Ω and ∂ T Ω ∩ ∂ u Ω = ∅ (Fig. 1). We impose a load T on the boundary ∂ T Ω and displacements u 0 on the boundary ∂ u Ω.
Since the volumetric forces are assumed to be zero, the equilibrium equations are given by where σ (M) denotes the Cauchy stress tensor and n (M) the outward unit normal vector at a point M ∈ ∂ T Ω. The strain tensor ε (M) is related to the displacement gradient tensor ∇u by the kinematic relations: (2) The Cauchy stress tensor σ (M) is related to the strain tensor ε (M) through the constitutive law, namely: The components of C, can, either be constant (homogeneous material), or depend on the points of the space (heterogeneous material). For isotropic and homogeneous materials, the constitutive law in the hypothesis of plane stress becomes: where E and ν are respectively Young's modulus and the Poisson's ratio and they are the parameters that we are going to identify. Tr (ε) is the trace of ε and I is the identity tensor. A direct mechanics problem is a classical problem for which the domain Ω governed by the Cauchy stress tensor C and subjected to the stress T, allows us to determine the unknown u knowing the Dirichlet boundary conditions. The determination of the unknowns u is possible by solving the equations (1), (2) and (3) using a numerical finite element method [30] in our case, whith Q4 elements and implemented in FEniCs [31,32] which is an open-source project for solving partial differential equations (PDEs) with the Finite Element Method (FEM) under python. This choice has been made because it is suitable for the use of the finite element method to solve partial differential equations.

Inverse problem 3.1 Finite element model updating (FEMU)
After obtaining the solution of the direct problem, we generate the synthetic displacement fields which takes the place of an experimental field u exp . The finite element model updating method consists in comparing the simulated displacements with the experimental displacements. To do this, we formulate a 2D cost function F (x) as shown in relation (5).
where m is the number of measurement points, u FE i∈{1,··· ,m} are the simulated displacements obtained by the finite element model and u exp i are the experimental displacements by a DIC method. Since the measurement points do not necessarily correspond to the mesh nodes, the interpolation is performed in the element where the measurement is located and compared to the experimental value. x represents the vector of parameters or properties of the specimen considered. The minimization of the functional F (x) by a well chosen algorithm allows to obtain the unknowns x of the problem. The minimization problem is expressed as in (6).
The flowchart of the FEMU method is shown in the Figure 2. In this work we have chosen an isotropic and heterogeneous model. For this, the unknown x will therefore be expressed as x = [E 1 , ν 1 , E 2 , ν 2 ] with E k , ν k the Young's modulus and the Poisson's ratio of the two homogeneous materials which will constituted the heterogeneous material with k = 1, 2.

Minimization algorithm using Levenberg-Marquardt algorithm
The cost function of the equation (5) can be written in a quadratic way as shown by (7): m represents the measurement points and r (x) is the vector of the differences between the simulated and the experimental displacements vector. J ∈ R m×n is the Jacobean matrix and n is the number of unknown parameters, a matrix that contains the derivatives of the function r (x). The The equation (9) can also be written as Levenberg (1944) and later Marquardt (1963) [33][34][35][36] have used an approximation of Gauss-Newton algorithm by using a Lagrange multiplier because the Gauss-Newton algorithm poses the problem of conditioning the matrix J t J. The Levenberg-Marquardt step h lm is defined by: J t J + ζI h lm = −g with g = J t r and ζ ≥ 0 (11) Here J = J (x), r = r (x), and ζ is the Lagrange multiplier. The way we will choose the values of ζ and the stopping criterion is detailed in [37].

Sensitivity analysis
In the process of minimizing the objective function, equation (5), we need to calculate the sensitivity matrix J. For this, we will use the finite difference method to evaluate each coefficient of this matrix. We will use the finite difference scheme called forward which is expressed as follow: ∆x j is a disturbance of the parameter x j and is well explained by the works of Lauwagie [38] and is expressed  by the equation (13).
According to the same author the stable interval where we can choose δ j is 10 −5 , 10 −2 . In our work all the values of δ j in this interval have given good results but δ j = 10 −5 was the best.

Design study Geometry
The specimen used for the different studies considered is a square specimen with checkerboards (Fig. 3). L = 20 mm is the size of the square and l = 5 mm is the size of a checkerboard. p 1 is the couple of parameters of the first material in gray color and p 2 is the couple of parameters of the second material in white color. E k and ν k with k = 1, 2 are respectively Young's modulus and Poisson's ratio and are constant in each checkerboard. On the top face of the specimen y = L, we have imposed a traction force T = 5 × 10 6 N/m. On the bottom face y = 0 a displacement u 0 = 0 has been imposed. The other sides are free of any load. The specimen thus constituted will be divided into quadrilateral (Q4) finite elements.

Table of parameters
Our study will focus on the study of the titanium (material made with checkerboards which have the same properties) to validate our algorithm and the heterogeneous material made of the titanium matrix and silicon carbide fibers (Tab. 1). We have chosen this heterogeneous material because it is used as an innovative material that has a great potential to increase the performance of gas turbine engines in aircraft [39]. We also have a reduction in weight by using this composite material [40]. Since we will use the exact properties mentioned in Table 1 to obtain synthetic data, we can calculate the identification errors using the following formulas: Where E i−imp and ν i−imp are the imposed or reference values that are in Table 1 and E i−est and ν i−est are the identified values on Young's modulus and Poison's ratio. e −i and η −i are the relative errors.
The different statistics of e −i and η −i can be obtained using Monte Carlo simulations. Generating a sufficient number of draws (denoted N MC ) of identified parameters using corrupted data, we can calculate the mean values of the identification relative errors: We can then define the standard deviation of the identification relative errors:

Noise modeling
To obtain synthetic experimental data closed to the real ones, we have chosen to corrupt our data obtained from the direct problem by adding a centered white Gaussian noise. Indeed the main sources of errors in DIC measurements are due to the camera noise and the matching errors during the correlation process. This latter can be expressed by a random error due to the speckle pattern properties (density, gradient, size of speckle), the algorithm, the subset size and the gray-level interpolation used for the DIC minimization process [41,42]. The uncertainty distribution has a null average and Gaussian profile, thus for a standard deviation σ of a white Gaussian noise the perturbations are modeled by equation (17).
with b (%) the percentage of noise, σ the standard deviation and |U | max the maximum displacement in absolute value.

Refinement problem
The sensitivity to the mesh refinement of the parameters which are going to be identified will be studied. For this we have generated a large mesh much more refined so that the nodes of the meshes in the direct problem are not necessarily in the same place as in the mesh when solving the inverse problem. We have then obtained other meshes with N x and N y divisions in the two directions of the space. In doing so, we take into account the errors due to the discretization.

Sensitivity of the procedure to initial values
The iterative Levenberg-Marquardt algorithm requires initial values to start the process. For this initiative, we have chosen 7 sets of initial parameters (see Tab. 2), while considering that the material is heterogeneous. For each case the properties of the material are identified and the results are shown in Figure 4. This figure shows that for each set, all values converge to an optimal solution minimizing the functional defined above. This is reflected by the fact that the objective function is separately convex and always positive. We also notice that the more the initial values are far from the reference values, the greater the number of iterations.

Sensitivity of the identified parameters to the mesh refinement
In this subsection, we study the sensitivity of the identified parameters to the mesh refinement. The number of elements in each case of mesh used here is in Table 3 and the noise level which we use is 1% and N MC = 50 which is the total number of draws. Then we represented the standard relative errors due to each parameter as a function of the number of divisions. Figure 5 shows that as the number of divisions increases, the standard deviations of the relative errors decrease. We lose in precision by passing from the homogeneous material to the heterogeneous one but the parameters are well identified for both homogeneous and heterogeneous material. We remark that when the space step is equal to 0.3 mm (Tab. 3) the relative errors are very small for the two materials. This means that all the parameters are well estimated and the direct problem has then been well achieved. The standard deviations of relative errors on Poisson's ratios are greater than those of Young's moduli and to emphasize this we have calculated the coefficients of variation C v (Tab. 4) of each parameter identified. This reflects the reality compared to what have been found in the literature by Pierron et al. [43]. We observe these large coefficients of variation C v on Poisson's ratios because they influence the real field less than the Young's moduli. Furthermore we assumed to have a homogeneous titanium material organized in a checkerboards so the Young's moduli are the same as for the Poisson's Ratios but after the identification we can remark that there are some errors of estimation as we can see it in Figure 5a. This is due to the errors on the synthetic experimental data and the effect of discontinuities which are here supposed to be negligible.

Sensitivity of the identified parameters to the noise level
After having studied the sensitivity of the parameters to the refinement of the mesh, we have chosen the case where the number of divisions is N x = N y = 64 and N MC = 50 as in the previous paragraph to study the sensitivity of the parameters to the noise level. This choice has been done because we have obtained good accuracies with this     refinement of the mesh. Thus we have studied the relative errors as a function of noise level. As the noise level increases, the standard deviations of the relative errors increase (Fig. 6). Analyzing the figure we remark that all the parameters, both in the homogeneous and the heterogeneous materials, are well identified. The parameters are almost insensitive to noise in the case of the homogeneous material as the Figure 6a shows it. The Young's moduli remain the same at this level of noise but the Poisson's ratios are slightly different when b = 5%. We can then say that a homogeneous material in this configuration can be a source of errors in an identification process. Considering the heterogeneous materials the large standard deviations are observed on the Poisson ratios because of what we earlier said in the previous paragraph.

Reconstruction of the displacement fields
The reference displacements are the displacements calculated using the reference parameters or properties of the materials which are prescribed in Table 1. Those displacements are here noted Ux and Uy respectively along each axis. The reconstructed displacements are calculated using the identified properties of the materials with 3% noise level and N MC = 1 using the identification procedure as explained in the previous paragraphs. Those displacements are noted U x − Rec and U y − Rec respectively along each axis of the space. We represent them as maps as we can see it in the Figures 7 and 8. The Figures 7a and 7b are for the material made of the same checkerboards (homogeneous material) with the titanium properties. Figures 8a  and 8b are the heterogeneous material.
Our objective in this paragraph is to analyze the identified parameters comparing the reference and the reconstructed displacements. Therefore we have calculated  The analysis of the Figures 7c, 7d, 8c and 8d shows that the differences are very small compared to the reference and the reconstructed displacements values. We therefore conclude that the reference and the reconstructed displacements values are very closed and that the parameters are well identified. The non zero differences in the case of the homogeneous (Figs. 7c and 7d) material are due to the fact that the identified and reference parameters are not exactly the same. There are some identification errors due to the errors on the experimental data. Even if the titanium material is homogeneous, it has been geometrically organized in checkerboards in the direct problem in this work. Because of this geometrical organization, we have not exactly the same periodic pair of parameters (E, ν) in each checkerboard after the identification. In the case of the heterogeneous material the discrepancies (Figs. 8c and 8d) are due to the causes that we have mentioned in the case of the homogeneous material.
The analysis of the reference and the reconstructed displacements has permitted us to show that the parameters are well identified in the two cases of materials (homogeneous and heterogeneous).

Conclusion
In this work, we have focused our studies on the identification of the mechanical properties of two types of specimens using a classical FEMU procedure: homogeneous material made of titanium and heterogeneous material made of two different homogeneous materials. The specimens were meshed by four-node quadrilateral elements and identification of the properties of those specimens has been carried. The sensitivity of the identified parameters to the number of divisions, to the choice of the initial parameter set and to the noise levels has been studied. As the number of divisions increases, the standard deviations of the relative errors decrease and as the noise level of the synthetic experimental data increases, these standard deviations increase. Considering the two types of material the parameters were well identified. We can also conclude that the algorithm is robust. A reconstruction of the displacement fields has been done after we have identified the parameters and we calculated the gaps between the reference and the obtained displacements. These gaps were very small compared to the reference and the reconstructed displacements values. This permits us to conclude that the identification process has been very well achieved.

Perspectives
The present study has been performed numerically, therefore our future studies will focus on experimental work, but the geometrical configuration will be a titanium matrix with circular inclusions to be more representative of the reality and 3D studies using Digital Volume Correlation (DVC) will be carried. In this work we have used the Levenberg-Marquardt algorithm which is a gradientbased algorithm limited by the number of parameters to be identified, however it is very efficient to identify parameters with good precisions and it converges very quickly compared to gradient-free algorithms [44]. To identify a lot of parameters of a material which properties are spatially distributed per node or element using finite element method, we have to use other algorithms. So we are planning to use a deep learning approach for identification on heterogeneous materials properties in the future.