Identification in transient dynamics using a geometry-based cost function in Finite Element Model Updating method

Finite Elements in Analysis and Design 122 (2016) 49–60

Contents lists available at ScienceDirect

Finite Elements in Analysis and Design journal homepage: www.elsevier.com/locate/finel

Identification in transient dynamics using a geometry-based cost function in Finite Element Model Updating method Clément Touzeau a, Benoit Magnain a, Bruno Emile b, Hélène Laurent a, Eric Florentin a,n a b

INSA-CVL / Université d'Orléans – PRISME, 88, boulevard Lahitolle, F-18020 Bourges, France Université d'Orléans / INSA-CVL – PRISME, 2 avenue François Mitterrand, F-36000 Châteauroux, France

art ic l e i nf o

a b s t r a c t

Article history: Received 1 April 2016 Received in revised form 18 July 2016 Accepted 6 September 2016

In this article, we introduce a new numerical method for identifying mechanical parameters of hyperelastic materials in a dynamic framework. Using a Finite Element Model Updating (FEMU) procedure, we propose a new cost function family. The goal is to avoid the use of a random speckle in association with a Digital Image Correlation (DIC) tool that are both needed when a cost function based on full displacement fields is used. The experimental data consist in a set of images. The developed method then uses a simple segmentation of these images without requiring any DIC procedure and associated random speckle technique. This advantage is made possible through the use of a new cost function based on geometry quantities. Examples based on synthetic data illustrate the performance of the proposed method on transient dynamics problems where the flow of information can be very important. & 2016 Elsevier B.V. All rights reserved.

Keywords: Identification Contactless measurements Finite element Transient dynamics Image segmentation

1. Introduction In a design study, engineers need information on the mechanical characteristics of the used material. To that end, standardized tests are performed on a sample of the considered material. In this article, we are interested in identification methods based on the use of contactless measurements. In that context, a very popular technique is Digital Image Correlation method (DIC) [22] that allows us to obtain the displacement fields: applying a random speckle on the specimen and the different images taken during the test are compared to lead to a full displacement field. Starting from the displacement fields, different methods have been developed in mechanical engineering to identify material parameters. A review is presented in [2]. Among the classical methods, the Constitutive Equation Gap Method (CEGM) [10,5], the Virtual Field Method (VFM) [11,3,12], the Equilibrium Gap Method (EGM) [6], the Reciprocity Gap Method (RGM) [13] can be mentioned. A very large amount of works are also based on the Finite Element Method Updating (FEMU) principle [14,7] applied in the specific field of dynamics [16]. The method used in this work relies on this FEMU principle: the results of a finite element model are compared to those extracted from images taken during the test. n

Corresponding author. E-mail address: eric.fl[email protected] (E. Florentin).

http://dx.doi.org/10.1016/j.finel.2016.09.003 0168-874X/& 2016 Elsevier B.V. All rights reserved.

Major problems are highlighted when random speckle in combination with DIC is used. The first one comes from the realization of the random texture itself whose quality largely depends on the speckle resolution and the resulting contrast. Interesting works are based on goal-oriented filtering technique in which data are combined into new output fields which are strongly correlated with specific quantities of interest [15]. The second problem comes from the DIC and the calculation cost it generates, especially when the number of images to be processed increases sharply. The intrinsic errors related to the DIC can also lead to difficulties, see [1,19,4]. We propose in this article to use geometric quantities (area and second moment of area) as base of the cost function for updating the finite element model. The proposed work thus consists in studying a variant of the FEMU for identify simple behaviors. Tracking random speckle is replaced by a simple detection of the specimen shape in the scene. DIC is therefore replaced by a simple segmentation of the picture easier to implement and much less CPU and memory consuming. In this work, all data are synthetic. This makes it possible to check the accuracy of the developed method. To its quality, the method is compared to a classical displacement based cost function. We show the developed method is an interesting proposal for simple behaviors. The paper is organized as follows: in Section 2, we the problem is defined. Section 3 is dedicated to the presentation of the proposed FEMU based identification method. Sections 4–6 present the

50

C. Touzeau et al. / Finite Elements in Analysis and Design 122 (2016) 49–60

results obtained on numerical tests, and carefully study the effect of various acquisition parameters that could potentially modify the precision.



where Td is the traction imposed to small surface dS and N 0 its normal vector. Dirichlet boundary conditions: u ¼ ud

on

Γu  τ

ð5Þ

2. Problem definition Let us consider θ ¼ ðθ1 ; …; θn Þ the set of parameters introduced in the constitutive relation of the material defined by: S¼

∂WðE; θÞ ∂E

 Unilateral contact conditions: g Z0

rn Z 0

and

g  rn ¼ 0

on

Γc

ð6Þ

ð1Þ

where S is the second Piola–Kirchhoff stress tensor, E is the Green–Lagrange strain tensor. W is the strain energy density function and depends of material parameters θ. The mechanical problem is defined considering large deformations in dynamic framework and possible contacts. Furthermore, the material studied is supposed to be isotropic, homogeneous and elastic. Thus, we define the set of equations used in the proposed method of identification.

where g is the gap between the solid and the obstacle and rn is the contact reaction force acting along the obstacle outward normal. In this study, the contact is considered without friction.

2.1. Equations of the problem

 Saint-Venant Kirchhoff strain energy density function

As shown in Fig. 1, the specimen may be represented as a closed subspace of R2 , denoted Ω, and its boundary ∂Ω on which one defines the partition:

In this first case, we identify 2 parameters, θ ¼ ðE; νÞ, respectively the Young modulus and Poisson's ratio of the considered material. The strain energy density function, WðE; ðE; νÞÞ, is defined as a function of E and θ ¼ ðE; νÞ are parameters:

∂Ω ¼ Γ u [ Γ σ [ Γ c ; ∅ ¼ Γu \ Γσ ¼ Γσ \ Γc ¼ Γu [ Γc ;



on

Γσ  τ

W ðE; ðE; νÞÞ ¼


 Eν E ðtrðE ÞÞ2 þ tr E 2 2ð1 þ νÞ 2ð1 þ νÞð1 2νÞ

ð7Þ

from which we obtain the constitutive relation linking S and E: S¼

Eν E trðE Þ:Id þ E ð1 þ νÞð1  2νÞ ð1 þ νÞ

ð8Þ

where Id is the identity matrix.

 Blatz–Ko strain energy density function ð3Þ

where P is the first Piola–Kirchhoff stress tensor, b represents the body forces, ρ is the density of the material and u€ is the acceleration vector. Neumann boundary conditions: P:N 0 ¼ T d

In this work we use two models of hyper-elastic behavior of which we briefly present the main elements.

ð2Þ

where the displacements are imposed on Γ u , the external forces are imposed on Γ σ and Γ c is a potential area of contact. The study is made on the time interval τ ¼ ½0; t f . In addition of the constitutive law (1) and within the total Lagrangian formulation, the physical problem is modeled by a set of relationships involving spatio-temporal quantities defined on Ω  τ:

 local equilibrium equations: DivP þ b ¼ ρu€ in Ω  τ

2.2. Constitutive relations

ð4Þ

In this second case, one parameter is identified, namely the shear modulus, θ ¼ ðGÞ. For practical reason, the elastic deformation potential of this second case refers to the invariants of right Cauchy–Green deformation tensor. The strain energy density function, WðE; ðE; νÞÞ is defined by as a function of E and θ ¼ ðGÞ are parameters:   pffiffiffiffiffiffi G C1 ð9Þ þ 2 C 3 5 ¼ WðE; GÞ wðC; GÞ ¼ 2 C2 where C 1 ¼ trðC Þ;

C2 ¼


 1
2 C  tr C 2 ; 2 1

C 3 ¼ detC

ð10Þ

Then: S¼

Fig. 1. Continuous problem.

∂WðE; GÞ ∂wðC; GÞ ¼2 ∂E ∂C

ð11Þ

and we classically obtain the constitutive relation: n o SðEÞ ¼ G Jð2E þ IdÞ  1  ð2E þIdÞ  2

ð12Þ

where pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi J ¼ detð2E þ IdÞ

ð13Þ

C. Touzeau et al. / Finite Elements in Analysis and Design 122 (2016) 49–60

2.3. Finite element analysis Field u is chosen in a finite element form: u ¼ NT U

ð14Þ

where N is the matrix containing the finite element shape functions and U the discrete displacement (nodal displacements) at instant t. Finite element method provides displacement fields u from its discrete form U solving: M U€ ðtÞ ¼ F ext ðtÞ  F int ðtÞ þRc ðtÞ

ð15Þ

where

    

51

One can alternatively define the n-order moment centered on G, the geometrical center of the studied solid: Z I n;G ¼ ðx  xG Þn dΩ ð18Þ Ω

Note: of course, in the special case where n¼ 0 the considered geometrical moment represents in dimension 2 a measure of the object surface area which will be denoted I0 hereinafter. The cost function, denoted C ðn;P;xÞ ðθÞ, is then: Z tf cal 2 ½I mes ð19Þ C ðn;P;xÞ ðθÞ ¼ n;P  I n;P ðθÞ dt 0

where and I cal n;P ðθÞ are the n-order moments based respectively on measure and finite element calculation, where P can be O or G as depicted before. We then find parameters values by identifying θid minimizing the function C ðn;P;xÞ ðθÞ: I mes n;P

M is the mass matrix, U€ is the discrete acceleration vector, F ext is the discrete external force vector, F int is the discrete internal force vector, Rc is the discrete contact reaction force vector.

θid ¼ argminC ðn;P;xÞ ðθÞ

To solve Eq. (15), we use an explicit integration scheme in time: the central difference scheme. Regarding the contact, we use the bi-potential method applied to the dynamic as it was developed in [8]. In the following, the displacement u corresponding to a coarse mesh and coarse time step is denoted ucal and is obtained using (14) and (15) from U cal . On the other hand, measurements that should normally be obtained experimentally via a contactless technique resulting in a video made of a succession of Nmes images acquired at different time steps t i A f1‥Ng is also simulated. The field umes is obtained using (14) and (15) from U mes . The simulation is done with a very fine mesh and time step. From this simulation we extract a set of pictures (synthetic data). We note f t1 ðxÞ…f tN ðxÞ the Nmes pictures available at these different time steps. An image is represented by the function f ðxÞ that associates a discrete value (ranging from 0 to 255 in the case of an 8-bit encoding) to each position x.

To perform this minimization, the geometrical moment I mes n;P is required for every instant t. Here we get this information by an image segmentation technique. The principle is to transform each image f ðxÞ in a binary picture f T ðxÞ based on a threshold T: f T ðxÞ ¼ 1

if f ðxÞ Z T

f T ðxÞ ¼ 0

otherwhise

We develop a variant of the FEMU, based on the geometry (Section 3.1) and compare it to a reference, the FEMU based on the full displacement field (Section 3.2). The two methods are compared in the next section (Section 4) for identifying material parameters. To this end, a cost function is used for comparing synthetic experimental results to calculated results. Our study focuses on the choice of the cost function to be minimized. FEMU identification based methods can be summarized (see [2] for more details) in searching θid such that: θ

ð16Þ

where CðθÞ is a cost function. This cost function CðθÞ compares experience (indexed mes) and calculated (indexed cal) quantities. 3.1. Proposed cost function: geometric based function We introduce here a family of cost functions CðθÞ ¼ C ðn;P;xÞ ðθÞ (see Eq. (19)) based on different geometrical quantities. On domain Ω, we define the different n-order geometrical moments as follows: Z xn dΩ ð17Þ I n;O ¼ Ω

where O is a fixed point, the origin of the coordinate system.

ð21Þ

Among various thresholding techniques, this technique called global thresholding is the simplest one. It is the one we used in the following work. 3.2. Reference cost function: kinematic Hereafter we introduce the classical cost function based on kinematic, that will be considered as the reference cost function in the numerical tests. Cost function C u ðθÞ is introduced, based on the following displacement: Z tf Z C u ðθ Þ ¼ ½umes  ucal ðθÞ2 dΩdt ð22Þ 0

3. FEMU based identification method

θid ¼ argminCðθÞ

ð20Þ

θ

Ω

The parameters values are then obtained by identifying θid that minimizes the function C u ðθÞ:

θid ¼ argminC u ðθÞ

ð23Þ

θ

Field umes is required and generally obtained using image correlation technique. The key idea is to minimize the gap between two successive images defined on Ω: f ðxÞ at instant t, and gðxÞ at instant t þdt. The displacement field is chosen in a finite element form at time t such that: Z ð24Þ umes ¼ argminu ¼ N T U ½f ðxÞ  gðx þ uÞ2 dΩ Ω

where N is the matrix containing the finite element shape functions and U the discrete displacement (nodal displacements) at instant t. Further details can be found on the image correlation technique in [17,18]. 3.3. Minimization algorithm We present in this subsection the algorithm we use for minimizing cost functions. It uses the bisection method (Algorithm 1). For simplicity, it is presented in the particular case of Saint-Venant Kirchhoff law with θ ¼ ðE; νÞ.

52

C. Touzeau et al. / Finite Elements in Analysis and Design 122 (2016) 49–60

Algorithm 1. Bisection method adopted to minimize CðθÞ.

C. Touzeau et al. / Finite Elements in Analysis and Design 122 (2016) 49–60

53

Table 1 Characteristics used to define synthetic data Mes1. Characteristics

Test case 1

Test case 2

NFE tf Δt Nmes

1024 2:5E  4 s 1:0E  7 s 250

41,730 4:0E  3 s 1:0E  7 s 400

Table 2 eE and eν errors (%) versus NFE. N FE ¼ 16

N FE ¼ 64

N FE ¼ 256

N FE ¼ 1024

Cu C ð0;P;xÞ C ð2;O;xÞ C ð2;G;xÞ

1.5 5:4E  1 1.64 1.92

0.49 1:3E  1 5:2E  2 2:6E  1

0.31 2:8E  3 4:1E  1 2:8E  3

1:1E  2 5:9E  4 4:5E  5 4:5E  5

Cu C ð0;P;xÞ C ð2;O;xÞ C ð2;G;xÞ

5.86 1.5 9.83 5.92

1.89 4:6E  1 3.32 8:4E  1

0.98 6:5E  2 2.67 1:3E  1

1.17 4:2E  3 7:3E  3 6:6E  2

Characteristics Fig. 2. Test case 1. eE ð%Þ

eν ð%Þ

scenario (external force for the first case and initial speed for the second one) and by the constitutive relation of the material to be identified (respectively Saint-Venant Kirchhoff and Blatz-Ko).

Fig. 3. Applied pressure load.

Fig. 4. Test case 2.

4.1.1. Test case 1 In the first test, a square of side L ¼ 0:1 m is placed on a rigid foundation when a dynamic pressure is applied on its top edge (Fig. 2). The profile of the transient load is shown as a function of time in Fig. 3. We use the Saint-Venant Kirchhoff behavior related to Young modulus E ¼ 69 GPa and Poisson's ratio ν ¼ 0:3. Density is ρ ¼ 700 kg m  3 . The maximum value of the applied pressure is equal to P max ¼ 1:0E10 Pa. The duration of the loading is set to t i ¼ 2:0E  4 s. In this test case, both the elastic modulus and Poisson's ratio are identified. 4.1.2. Test case 2 In the second test (Fig. 4), a disc of R ¼ 2:5E  2 m radius is launched with an initial velocity V 0 ¼  4:43 m s  1 onto a rigid barrier. The initial time is the very first time of contact between the disc and the obstacle. The duration of experience is t f ¼ 4:0E  3 s. The Blatz-Ko behavior relationship has been selected to model the material with a shear modulus initially taken at G ¼ 5:0E5 Pa and a density ρ ¼ 700 kg m  3 . In test case 2, the single shear modulus parameter is identified. 4.2. Synthetic data used

4. Numerical illustration In this section, we present results obtained with two test cases in order to illustrate the performance of the developed methodology. Both test cases are first introduced in details. Then, we present the two types of measures which have been considered. Finally a detailed review of the results is carried out. 4.1. Presentation of test cases Two test cases are used to compare the cost functions C ðn;P;xÞ ðθÞ and C u ðθÞ. These two examples are distinguished both by their

In this work, we use synthetic measured data. In this comparative study, we proceed in two steps to evaluate the performance of the proposed cost function. In the first stage (Section 4.2.1: Mes1), we consider that the measured field is perfect (with respect to the machine accuracy) while the second step (Section 4.2.2: Mes2) introduces the intrinsic defects coming from images processing. 4.2.1. Mes1: perfect measurements The measurements Mes1 result of a converged finite element calculation. Specifically, the calculation is made using an implicit integration scheme with a very small time step Δt. We also use a

54

C. Touzeau et al. / Finite Elements in Analysis and Design 122 (2016) 49–60 10 C C C

1

0.0102

C

0.01 0.0098 0.0096

I0 (m2)

eE (%)

0.1

0.01

0.001

0.0094 0.0092 0.009 0.0088

0.0001

P =1.0E7 P =1.0E8 P =1.0E9 P =1.0E10

0.0086

1e-05 10

100

1000

10000

0.0084

0

5e-05

NFE

C C C

eν (%)

0.1

0.01

1000

10000

Fig. 6. eν error (%) versus NFE.

Table 3 eE and eν errors (%) versus ΔtðsÞ. Characteristics

eE ð%Þ

eν ð%Þ

Cu C ð0;P;xÞ C ð2;O;xÞ C ð2;G;xÞ Cu C ð0;P;xÞ C ð2;O;xÞ C ð2;G;xÞ

Δt ¼ 1:0E  7 s

Δt ¼ 2:0E  7 s

Δt ¼ 1:0E  8 s

0.31 2:8E  3 4:1E  1 2:8E  3 0.98 6:5E  2 2.67 1:3E  1

0.3 8:5E  2 3:7E  1 8:5E  2 1.5 4:5E  1 2.21 3:2E  1

0.32 8:5E  2 4:2E  1 1:9E  1 0.52 3:9E  1 2.8 7:2E  1

Table 4 eE and eν errors (%) versus Nmes. Characteristics

eE ð%Þ

eν ð%Þ

Cu C ð0;P;xÞ C ð2;O;xÞ C ð2;G;xÞ Cu C ð0;P;xÞ C ð2;O;xÞ C ð2;G;xÞ

0.00025

consistent mass matrix and a very fine mesh. For each test case, Table 1 summarizes useful features: NFE is the number of finite elements and Nmes the number of measurement points stored as synthetic data. All the quantities useful for cost functions considered in this work are calculated during the finite elements process and therefore benefit from the machine precision. The advantage of this synthetic perfect measure is to test the quality of cost functions regardless of the quality of image processing methods used to obtain the measured fields (image correlation for the cost function Cu, for example).

C

1

NFE

0.0002

Fig. 7. Influence of P max ðPaÞ on I 0 ðtÞ.

10

100

0.00015 time (s)

Fig. 5. eE error (%) versus NFE.

0.001 10

0.0001

N mes ¼ 25

N mes ¼ 125

N mes ¼ 250

0.3 1:7E  2 4:6E  1 1:7E  2 1.5 1:3E  1 2.99 6:5E  2

0.31 1:7E  2 4:0E  1 1:7E  2 0.98 1:3E  1 2.54 6:5E  2

0.31 2:8E  3 4:1E  1 2:8E  3 0.98 6:5E  2 2.67 1:3E  1

4.2.2. Mes2: digital images To get closer to actual conditions of an identification process, we use the Mes1 measurement and the post-processing software GMSH [9] to build over time a series of Nmes images for each test case. Several output formats are considered in this study (compressed and uncompressed). This sequence of digital images is then post-treated to obtain the temporal evolution of the geometric quantities (I0, I 2;O and I 2;G ). To that end, a basic segmentation method is performed. Based on the analysis of the gray-level histogram, a simple thresholding allows us to convert images into binary ones. The afterwards extracted Mes2 measurements will then be compared to Mes1 ones.

5. Results on Mes1 In this section, we explore the potential of cost functions using geometric moments when the measured field comes from a converged Finite Element Analysis (Mes1). The performance and robustness of the method are studied in detail. The obtained results are particularly highlighted by comparison with those obtained using the conventional cost function with respect to the kinematic fields. 5.1. Mesh influence Based on test case 1 with t f ¼ 2:5E  4 s, Table 2 shows the evolution of relative identification errors eE ¼ j E  Eid j =E and eν ¼ j ν  νid j =ν, for each function and cost depending on the mesh size used. Four meshes are considered from (N FE ¼ 16) to (N FE ¼ 1024).

C. Touzeau et al. / Finite Elements in Analysis and Design 122 (2016) 49–60

55

Table 5 eE and eν errors (%) versus Pmax(Pa). P max ¼ 1:0E7 Pa

P max ¼ 1:0E8 Pa

P max ¼ 1:0E9 Pa

P max ¼ 1:0E10 Pa

Cu C ð0;P;xÞ C ð2;O;xÞ C ð2;G;xÞ

0.16 4:0E  2 1:4E  1 4:0E  2

0.16 3:0E  2 1:5E  1 4:0E  2

0.22 1:1E  2 1:8E  1 1:7E  2

0.31 2:8E  3 4:0E  1 2:8E  3

Cu C ð0;P;xÞ C ð2;O;xÞ C ð2;G;xÞ

0.59 6:5E  2 6:5E  1 6:5E  2

0.59 6:5E  2 6:5E  1 6:5E  2

0.26 6:5E  2 7:8E  1 6:5E  2

0.98 6:5E  2 2.67 1:3E  1

Characteristics

eE ð%Þ

eν ð%Þ

1

10 C C C

C

eG (%)

eE (%)

0.1

1

0.01

C C C

C

0.001 1e+07

1e+08

1e+09

1e+10

Pmax (Pa)

0.1 100000

1e+06

1e+07

1e+08

G(Pa)

Fig. 8. eE error (%) versus Pmax(Pa).

Fig. 10. eG error (%) versus G(Pa).

10

eν (%)

1

0.1

C C C

0.01 1e+07

1e+08

1e+09

C

1e+10

Pmax (Pa)

Fig. 9. eν error (%) versus Pmax(Pa).

It is observed that the geometric cost function C ð0;P;xÞ leads to the best results for coarse meshes. However, other costs functions lead to satisfactory results. Figs. 5 and 6 show that, for equal mesh, geometric cost functions allow better identification or at least an identification of similar quality. In the following, all identification processes are performed on the basis of the third mesh (N FE ¼ 256). 5.2. Time step influence We are now studying the influence of the time step used in the identification process. The explicit integration scheme requires

Fig. 11. Segmentation obtained after thresholding of a JPEG image with T ¼ 254=255.

checking a criterion to be stable (Courant–Friedrichs–Lewy) and thus Δt takes the three values defined in Table 3. The results in Table 3 show the robustness of the discussed methods with respect to time step.

56

C. Touzeau et al. / Finite Elements in Analysis and Design 122 (2016) 49–60 0.00205

1e-11

0.002 1e-12

0.0019

C(2,G,x)

I0 (m2)

0.00195

0.00185 0.0018

1e-14

0.00175 0.0017

1e-13

T=252/255 T=253/255 T=254/255 Mes 1

T=252/255 T=253/255 T=254/255 Mes1

0

0.0005

0.001

0.0015

0.002

0.0025

0.003

0.0035

1e-15 200000

0.004

300000

400000

500000

600000

700000

800000

G(Pa)

time (s)

Fig. 12. Influence of segmentation thresholds on I 0 ðtÞ estimates.

Fig. 15. Influence of segmentation thresholds on C ð2;G;xÞ .

1e-05

0.002

0.00195 1e-06

2

I0 (m )

C(0,P,x)

0.0019

1e-07

0.00185

0.0018 1e-08

1e-09 100000

200000

T=254/255 T=254/255,se=Disk,r=3px T=254/255,se=Disk,r=7px T=254/255,se=Square,s=5px T=254/255,se=Square,s=13px Mes1

0.00175

T=252/255 T=253/255 T=254/255 Mes 1

300000

400000

500000

600000

700000

0.0017

800000

0

0.0005

0.001

0.0015

0.002

0.003

0.0035

0.004

Fig. 16. Influence of morphological operations on I 0 ðtÞ estimates.

Fig. 13. Influence of segmentation thresholds on C ð0;P;xÞ .

1e-11

1e-05

1e-12

1e-06

C(0,P,x)

C(2,O,x)

0.0025

time (s)

G(Pa)

1e-13

1e-14

1e-07

1e-08 T=252/255 T=253/255 T=254/255 Mes 1

1e-15 200000 300000 400000 500000 600000 700000 800000 900000

1e+06

G(Pa)

Fig. 14. Influence of segmentation thresholds on C ð2;O;xÞ .

1e-09 100000

T=254/255 T=254/255,se=Disk,r=3px T=254/255,se=Disk,r=7px T=254/255,se=Square,s=5px T=254/255,se=Square,s=13px Mes 1

200000

300000

400000

500000

600000

700000

800000

G(Pa)

Fig. 17. Influence of morphological operations on C ð0;P;xÞ .

5.3. Influence of the number of records

5.4. Influence of the load intensity

In this part, we study the influence of the number of records available to make the identification. For the same value of tf and Δt, we vary the amount of measurement points as shown in Table 4. Again, the presented methods prove robust.

We now vary the intensity of the load. As shown in Fig. 7, we can thus act on the amplitude of the deformations undergone by the domain Ω. Table 5 shows identification errors for values of Pmax between 1:0E7 Pa and 1:0E10 Pa.

C. Touzeau et al. / Finite Elements in Analysis and Design 122 (2016) 49–60

57

1e-11

C(2,O,x)

1e-12

1e-13

1e-14

1e-15 200000

T=254/255 T=254/255,se=Disk,r=3px T=254/255,se=Disk,r=7px T=254/255,se=Square,s=5px T=254/255,se=Square,s=13px Mes 1

300000

400000

500000

600000

700000

800000

900000

1e+06

G(Pa)

Fig. 18. Influence of morphological operations on C ð2;O;xÞ .

1e-11

C(2,G,x)

1e-12

1e-13

1e-14

1e-15 200000

T=254/255 T=254/255,se=Disk,r=3px T=254/255,se=Disk,r=7px T=254/255,se=Square,s=5px T=254/255,se=Square,s=13px Mes 1

300000

400000

500000

600000

700000

800000

G(Pa)

Fig. 19. Influence of morphological operations on C ð2;G;xÞ .

Table 6 eG errors (%) achieved for various segmentation parameters. Characteristics

C ð0;P;xÞ

C ð2;O;xÞ

C ð2;G;xÞ

T ¼ 252=255 T ¼ 253=255 T ¼ 254=255 T ¼ 254=255;  se ¼ Disk;  r ¼ 3px T ¼ 254=255;  se ¼ Disk;  r ¼ 7px T ¼ 254=255;  se ¼ Square;  s ¼ 5px T ¼ 254=255;  se ¼ Square;  s ¼ 13px Mes1

2.57 1.33 70.0 3.89 1.82 3.89 9:7E  1 6:3E  1

2.10 1.43 26.84 4.36 3.59 4.36 3.04 1:6E  1

1.72 2.18 44.39 6.84 4.68 6.84 3.64 2:6E  1

The results show that, with Mes1, all costs functions do not deteriorate, even for very small changes in the physical quantity used. Fig. 20. Pixelisation levels P1, P3 and P6.

5.5. Test case 2

6. Results using Mes2

A similar study was conducted for test case 2. The conclusions that can be drawn are identical to those described above. In addition, we studied the influence of the stiffness of the material on the quality of identification. Fig. 10 shows the evolution of eG ¼ j G  Gid j =G according to G.

In this second part, reference measurements are built with the post-processing software GMSH, allowing us to obtain, for each situation, a series of images representing the studied test case. As no speckle is available we still consider here geometric cost functions whose interest has been highlighted in the previous

58

C. Touzeau et al. / Finite Elements in Analysis and Design 122 (2016) 49–60 1e-11

Table 7 Levels of pixelisation. Object pixels

Contour pixels

Ratio (%)

P1 P2 P3 P4 P5 P6

735,902 257,867 116,270 33,736 5210 402

7744 4584 3072 1656 644 180

1.05 1.77 2.64 4.90 12.36 44.44

1e-12

C(2,O,x)

Characteristics

1e-13

1e-14

0.00205

Mes1 P1 P2 P3 P4 P5 P6

0.002 1e-15 300000

400000

500000

I0 (m2)

0.00195

600000

700000

800000

900000

1e+06

G(Pa)

Fig. 23. Influence of pixelisation of images on C ð2;O;xÞ . 0.0019

Table 8 eG errors (%) resolutions.

0.00185 Mes1 Mes2 with P1 Mes2 with P2 Mes2 with P3 Mes2 with P4 Mes2 with P5 Mes2 with P6

0.0018

0.00175

0

0.0005

0.001

0.0015

0.002

0.0025

0.003

0.0035

0.004

time (s)

Fig. 21. Influence of pixelisation on calculation of I 0 ðtÞ estimates.

achieved

for

various

image

Characteristics

C ð0;P;xÞ

C ð2;O;xÞ

Mes1 P1 P2 P3 P4 P5 P6

6:3E  1 2:3E  1 1.58 1.97 1:3E  1 1.96 60.0

1:6E  1 6:3E  1 1:6E  1 2.02 1.14 5.1 80.0

1e-05

0.010000

1e-06

0.010000

1e-07

I0 (m2)

C(0,P,x)

0.010000

1e-08

1e-09 200000

Mes1 P1 P2 P3 P4 P5 P6

300000

400000

500000

600000

700000

0.010000 0.009999 0.009999

800000

G(Pa)

Fig. 22. Influence of pixelisation on C ð0;PxÞ .

0.009999 0.009999

Mes1 Mes2

0

5e-05

0.0001

0.00015

0.0002

0.00025

time (s)

section. The objective is then to mobilize usual image processing techniques to extract mechanical characteristics and to study the influence of the chosen image processing parameters on the identification accuracy. To that end, we focused on test case 2 in which the identification of the shear modulus of the Blatz-Ko law is pursued. 6.1. Segmentation influence The obtained images present the launched disk clearly distinguishable against the background. As in many applications where the gray levels of pixels belonging to the object are substantially different from the ones of the background, thresholding constitutes a simple but effective tool to separate objects from the background [20]. A basic thresholding can then be considered to process the segmentation step. In our synthetic frames, white background has been used. The threshold value T to choose must

Fig. 24. Influence of Pmax(Pa) on I 0 ðtÞ estimates (P max ¼ 1:0E7 Pa).

therefore simply allow us to eliminate pure white. We then set T¼ 254/255. If this threshold operates when uncompressed formats are used, it can be defective if JPEG format is considered. The image compression indeed induces some approximations in pixel quantization. This can result, after binarization, in erratic foreground pixels detection as shown in Fig. 11. With T ¼ 254=255, numerous pixels of the background are falsely incorporated into the object mask altering scale and object area estimation and more generally the computation of geometric cost functions, making them less accurate or even ineffective. One way to overcome this problem is to relax constraints by choosing a lower threshold. Figs. 12, 13 and 14 and 15 illustrate this situation by respectively presenting the evolution over time of I0, C ð0;P;xÞ , C ð2;O;xÞ and C ð2;G;xÞ estimates with thresholds T ¼ 253=255 and T ¼ 252=255.

C. Touzeau et al. / Finite Elements in Analysis and Design 122 (2016) 49–60

59

1000

0.01002 0.01001

C C C C C C

100

mes1 mes2 mes1 mes2 mes1 mes2

0.01001

10

eE (%)

I0 (m2)

0.01001 0.01000

1

0.01000

0.1 0.00999

0.01

0.00999 Mes1 Mes2

0.00998

0

5e-05

0.0001

0.00015

0.0002

0.00025

time (s)

0.01

I0 (m2)

0.00998 0.00996 0.00994 0.00992 0.0099 Mes1 Mes2

5e-05

0.0001

1e+09

1e+10

Fig. 28. eE errors (%) versus Pmax(Pa) (comparison between Mes1 and Mes2).

0.01002

0

1e+08

Pmax

Fig. 25. Influence of Pmax(Pa) on I 0 ðtÞ estimates (P max ¼ 1:0E8 Pa).

0.00988

0.001 1e+07

0.00015

0.0002

0.00025

time (s)

Fig. 26. Influence of Pmax(Pa) on I 0 ðtÞ estimates (P max ¼ 1:0E9 Pa).

respectively present the evolution over time of I0, C ð0;P;xÞ , C ð2;O;xÞ and C ð2;G;xÞ estimates with threshold T ¼ 254=255 and the above mentioned morphological parameters. Table 6 summarizes the identification errors achieved for each segmentation process. However the cost function may be defined, the use of raw binarization (T ¼ 254=255, without morphological operations) does not allow a good identification of the shear modulus. By contrast, as soon as the segmentation threshold is relaxed or morphological post-treatment is applied, the identification is quite satisfactory. The proposed method is then quite robust to the segmentation step. Any of the considered geometrical cost functions may operate properly, provided that acquisition conditions present the object on a uniform background. However, best results are obtained considering either the relaxed threshold T ¼ 253=255 or the raw thresholding (T ¼ 254=255) combined with morphological opening. 6.2. Image sampling influence

0.0102 0.01 0.0098

I0 (m2)

0.0096 0.0094 0.0092 0.009 0.0088 0.0086 0.0084

Mes1 Mes2

0

5e-05

0.0001

0.00015

0.0002

0.00025

time (s)

Fig. 27. Influence of Pmax(Pa) on I0 ðtÞ estimates (P max ¼ 1:0E10 Pa).

An alternative to reduce the influence of noise that is due to the use of compressed format and raw binarization is to apply morphological operations [21]. To that end, a morphological opening operator seems to be the best candidate to eliminate isolated black pixels while preserving the essential shape characteristics of the object. As images structures are selectively filtered out, the selection depending on the used structural element (se), we tested several shapes and sizes, namely a disk with 3 or 7 pixels radii (r) and a square with 5 or 13 pixels sides (s). Figs. 16, 17, 18 and 19

In order to study the influence of image resolution on the performances, we simulated various image sampling, thus varying the number of pixels effectively characterizing the object. The image format is uncompressed PNG that offers the best image quality. The study is conducted on I0 and I 2;O with 6 pixelisation levels named P1 to P6. For the PNG format, the threshold is fixed to T ¼ 254=255. As shown in Fig. 20, P1 presents the highest number of pixels constituting the object and consequently the most precise resolution. Table 7 presents, for each level of sampling, the number of pixels within the object, the number of pixels on the object contour and the ratio between perimeter pixels and area pixels. Fig. 21 shows that the image sampling may impact the detection of the deformable body and therefore the computation of C ð0;P;xÞ , when extremely low resolution is employed. Figs. 22 and 23 show the evolution of cost functions C ð0;P;xÞ and C ð2;O;xÞ for different levels of pixelation. As expected, the quality of identification increases with the level of resolution. Whatever the cost function used, P6 level does not allow us to identify the shear modulus. For other resolutions, a good estimate of the shear modulus is still possible even at low resolutions (see Table 8). 6.3. Test case 1 On Mes1, the machine precision allowed a very good identification even if the deformation caused by loading was very low (Figs. 7–9). We also studied the behavior of Mes2 face to such

60

C. Touzeau et al. / Finite Elements in Analysis and Design 122 (2016) 49–60

situations. We used PNG images, with a P3 pixelisation level and a segmentation threshold T ¼ 254=255. Figs. 24–27 present, depending on Pmax(Pa) applied, the I 0 ðtÞ estimates obtained by Mes1 and Mes2. It is clearly shown that when the magnitude of the load is too low, the deformations of the solid can be too small to be detected by image processing. In Fig. 28, a final comparison is proposed. We find:

 a loss of precision when moving from Mes1 to Mes2, for all cost 

functions, a real limitation of the method in the real case (with Mes2) when amplitudes of deformations are very small (image processing measurements do not contain the information).

[3]

[4]

[5]

[6] [7] [8]

[9]

7. Conclusions [10]

In this article, a new method of identification of material parameters is tested in the framework of dynamics. The identification problem based on FEMU method involves a minimization problem of carefully chosen cost functions. The novelty of the present work relies in the choice of a geometrical cost function. The advantage of such a cost function is that it only requires an image segmentation operation to be estimated from the measurements. This measure allows us to filter essential information from a wide amount of available data, which is interesting in transient dynamics. The numerical results presented illustrate the interest of the method (quality, robustness, etc.). In future development, this idea may be used in case of more complex constitutive equations. Indeed, the method seems adapted to simple test cases where the solution is unique, as treated in this work. More complex laws could lead to non unique solution (for example heterogeneities in material). An indicator based on the derivative could be constructed in such cases.

References

[11]

[12]

[13]

[14] [15]

[16] [17]

[18]

[19]

[20] [21]

[1] M. Alfano, G. Lubineau, G.H. Paulino, Global sensitivity analysis in the identification of cohesive models using full-field kinematic data, Int. J. Solids Struct. 55 (2015) 66–78. [2] S. Avril, M. Bonnet, A.S. Bretelle, M. Grédiac, F. Hild, P. Ienny, F. Latourte, D. Lemosse, S. Pagano, E. Pagnacco, F. Pierron, Overview of identification

[22]

methods of mechanical parameters based on full-field measurements, Exp. Mech. 48 (4) (2008) 381–402. S. Avril, J.M. Huntley, F. Pierron, D.D. Steele, 3d heterogeneous stiffness reconstruction using mri and the virtual fields method, Exp. Mech. 48 (4) (2008) 479–494. B Pan BW, G. Lubineau, Comparison of subset-based local and fe-based global digital image correlation: theoretical error analysis and simulation, Opt. Lasers Eng. V. 82 (2016) 148–158. M. Ben Azzouna, P. Feissel, P. Villon, Robust identification of elastic properties using the modified constitutive relation error, Comput. Methods Appl. Mech. Eng. 295 (2015) 196–218. D. Claire, F. Hild, S. Roux, Identification of damage fields using kinematic measurements, C. R. Mécan. 330 (11) (2002) 729–734. J.D. Collins, G.C. Hart, T.K. Haaselman, B. Kennedy, Statistical identification of structures, AIAA J. 12 (2) (1974) 185–190. Z.Q. Feng, B. Magnain, J.M. Cros, Solution of large deformation impact problems with friction between Blatz–Ko hyperelastic bodies, Int. J. Eng. Sci. 44 (2006) 113–126. C. Geuzaine, J.F. Remacle, Gmsh: a three-dimensional finite element mesh generator with built-in pre- and post-processing facilities, Int. J. Numer. Methods Eng. 79 (11) (2009) 1309–1331. G. Geymonat, F. Hild, S. Pagano, Identification of elastic parameters by displacement field measurement, C. R. Mécan. 330 (6) (2002) 403–408. M. Grédiac, Principe des travaux virtuels et identification, C. R. De. l'Acad. Des. Sci. Sér. 2 Mécan. Phys. Chim. Sci. De. l'univers Sci. De. La. Terre 309 (1) (1989) 1–5. M. Grédiac, F. Pierron, S. Avril, E. Toussaint, The virtual fields method for extracting constitutive parameters from full-field measurements: a review, Strain 42 (4) (2006) 233–253. M. Ikehata, Inversion formulas for the linearized problem for an inverse boundary value problem in elastic prospection, SIAM J Appl. Math. 50 (6) (1990) 1635–1644. K.T. Kavanagh, R.W. Clough, Finite element applications in the characterization of elastic solids, Int. J. Solids Struct. 7 (1) (1971) 11–23. G. Lubineau, A goal-oriented field measurement filtering technique for the identification of material model parameters, Comput. Mech. 44 (5) (2009) 591–603. M.I. Friswell JEM, Finite Element Model Updating in Structural Dynamics, Kluwer, London, 1995. J.C. Passieux, J.N. Périé, High resolution digital image correlation using proper generalized decomposition: Pgd-dic, Int. J. Numer. Methods Eng. 92 (6) (2012) 531–550. J.C. Passieux, P. Navarro, J.N. Périé, S. Marguet, J.F. Ferrero, A digital image correlation method for tracking planar motions of rigid spheres: application to medium velocity impacts, Exp. Mech. 54 (8) (2014) 1453–1466. R. Tao, G.L.A. Moussawi, B. Pang, Accurate kinematic measurement at interfaces between dissimilar materials using conforming finite-element-based digital image correlation, Opt. Lasers Eng. 81 (2016) 103–112. M. Sezgin, B. Sankur, Survey over image thresholding techniques and quantitative performance evaluation, J. Electron. Imaging 13 (1) (2004) 146–168. P. Soille, Morphological Image Analysis: Principles and Applications, Springer, Heidelberg, 2013. M. Sutton, W. Wolters, W. Peters, W. Ranson, S. McNeill, Determination of displacements using an improved digital correlation method, Image Vision. Comput. 1 (3) (1983) 133–139.