Modified constitutive relation error identification strategy for transient dynamics with corrupted data: The elastic case

Comput. Methods Appl. Mech. Engrg. 196 (2007) 1968–1983 www.elsevier.com/locate/cma

Modified constitutive relation error identification strategy for transient dynamics with corrupted data: The elastic case P. Feissel a

a,*

, O. Allix

b

Laboratoire Roberval de Me´canique, Universite´ de Technologie de Compie`gne, BP 20529, rue Personne de Roberval, 60205 Compie`gne Cedex, France b LMT-Cachan, ENS de Cachan/CNRS/Universite´ Paris VI, 61, avenue du Pre´sident Wilson, 94235 Cachan Cedex, France Received 10 July 2006; received in revised form 30 October 2006; accepted 31 October 2006

Abstract In this paper, we present an identification strategy for transient dynamic tests with corrupted data and its first development for the identification of elastic properties. This method is based on the modified constitutive relation error principle, which was initially introduced for model updating in vibration analysis. The governing principle of the method consists in splitting the equations of the problem into two groups, a ‘‘reliable’’ group and an ‘‘uncertain’’ group, and defining mechanical fields which verify the equations of the reliable group exactly and those of the uncertain group only approximately. Then, these fields are used to estimate the quality of the model to be identified. We tested our method on various examples and found it to be very robust with respect to high levels of perturbations. Furthermore, we show that this good performance is due to the fact that the method takes into account all the experimental data in a single calculation and introduces a distance between the measurements and the calculation. We also propose an alternative strategy dedicated to cases where there is no a priori knowledge of the perturbations.  2006 Elsevier B.V. All rights reserved. Keywords: Dynamics; Identification; Inverse problem; Uncertain measurements

1. Introduction This work was undertaken in order to identify a dynamic damage model for the prediction of the rupture of composite crash absorbers. The identification tests whose results were available had been performed on a system derived from a Split Hopkinson Pressure Bar (SHPB) apparatus. Due to the fracture phenomena associated with damage and localization, the data obtained from those tests, namely the velocities and tractions at the ends of the specimen, were unreproducible and highly scattered. Therefore, we sought an identification strategy which would be as unaffected by this scatter as possible.

*

Corresponding author. Tel.: +33 3 4423 4604; fax: +33 1 4740 2785. E-mail address: [email protected] (P. Feissel). URLs: http://www.utc.fr/lrm (P. Feissel), http://www.lmt.ens-cachan. fr (O. Allix). 0045-7825/$ - see front matter  2006 Elsevier B.V. All rights reserved. doi:10.1016/j.cma.2006.10.005

In this context, a natural choice would be to use a Kalman filtering approach. Kalman filters, which were introduced specifically for such uncertain frameworks in control, were first developed for state evaluation in the linear case [12], then applied to nonlinear models (through extended Kalman filters) and to parameter identification [16]. When addressing problems with large numbers of degrees of freedom in the state vector, Kalman filtering approaches appear to create stability and convergence problems [10]; this was confirmed when we tried to apply extended Kalman filters to our case [3]. Therefore, these approaches cannot be applied directly. One way to improve the Kalman approaches relies on ‘‘Unscented Kalman filters’’ [20,9], whose development seems promising but in which some questions, such as the choice of the parameters of the algorithm, remain unanswered. We decided to explore another option by developing a variational inverse approach. Coupling such an approach

P. Feissel, O. Allix / Comput. Methods Appl. Mech. Engrg. 196 (2007) 1968–1983

with a Kalman filtering approach should be an interesting extension to our work. In the framework of variational approaches, Tikhonov regularization techniques [19] enable one to improve the performance of the methods with respect to perturbations of the data [4]. However, the question of the choice of the weighting factor leading to a relevant solution remains and requires some a priori knowledge of the perturbations [7]. In [17], the authors propose a strategy for choosing the weighting factor based on the minimization of a criterion. What also makes this method interesting is that it provides a balance between the experiments and the theory which guarantees its robustness with respect to the perturbations. In [8], Tikhonov regularization coupled with an iterative strategy for Cauchy inverse problems in statics dispenses the measured values from being enforced strictly when calculating the inverse step. Both the idea of a balance between a model error term and a measurement error term and the idea of relaxing the measured values remain in the methods based on the Constitutive Relation Error (CRE) principle [14]. This approach was developed mainly for model updating in vibration analysis [11]. Its governing principle consists in dividing the relations which characterize the problem into a reliable set and an uncertain set. The latter is dealt with only approximately throughout the strategy. Subsequently, this approach was extended to transient dynamics with corrupted measurements, first in the linear elasticity framework in order to test its capabilities. The first developments, in [2], showed promising results. [1] described the numerical strategy associated with the resolution of the non-standard mechanical problems involved in the CRE method. In the present paper, our goal is, first, to confirm the potential of the method and its robustness on various examples, then to improve the understanding of its performance so that some improvements can be proposed. In the first part of the paper, we present the identification framework and describe the various identification strategies which have been tested (i.e. the Rota-based method, the basic CRE method and alternative formulations derived from the CRE method). The results presented in the following parts are based on short-duration tests only, but other test conditions would not affect our conclusions about the method. First, we apply the standard CRE method to the simple example of a homogeneous bar up to high perturbation levels and confirm its robustness. Then, we show that this method still gives good results on more complex examples, such as the identification of property fields. The last part of the paper focuses on understanding the robustness of the method, which is related to the fact that it takes into account all the experimental data in one calculation and to the introduction of a distance between the simulation and the experimental measurements. Such an understanding enables us to improve the effectiveness of the method and to propose an interesting

1969

alternative for cases with no a priori knowledge of the perturbations. 2. Identification framework: elastic properties in transient dynamics In this paper, we chose to present the proposed strategy in the simplified case of the identification of the properties of an elastic rod, but the method can also be applied to more complicated cases, such as field measurements. Furthermore, work is in progress to extend the method to the nonlinear case. Nevertheless, this simplified framework is sufficient to test the capabilities of the method. The presented work will only deal with synthetic datas created by Finite Element simulations. The corresponding direct problem consists in loading a one-dimensional elastic specimen of length L over a period of time [0, T]. The measured quantities, similar to what could be obtained from an SHPB test, are the displacements and tractions at both ends of the rod. Therefore, the direct problem is described by the following equations and data, where u is the scalar displacement and r the scalar stress: • • • • •

equilibrium: q€uðx; tÞ  r;x ðx; tÞ ¼ 0; constitutive relation: r = Eu,x; _ 0Þ ¼ u_ 0 ; initial conditions: u(x, 0) = u0 and uðx; measured displacements: ~ud ðtÞ at both x = 0 and x = L; measured tractions: f~ d ðtÞ at both x = 0 and x = L. Let us note that the domain is a unidimensional bar and its section is arbitrary and plays no role here. Therefore, the boundary condition in terms of stresses can be expressed: rn = fd. f~ d is hence chosen to be a traction per unit surface.

Let us observe that if all these equations were verified this would lead to an ill-posed direct problem, the boundary values being overspecified. In fact, it is this redundancy which makes the identification possible. Now, the objective is to identify the Young’s modulus properties of the rod from the measurements. We considered various cases: • The rod is homogeneous and one seeks its Young’s modulus E. • The rod is heterogeneous and composed of n parts, and the spatial distribution of the properties is known. One seeks the values of the Young’s modulus in each part (E1, . . . , En). • The rod is heterogeneous and one seeks E(x) in the absence of any a priori information. • The rod is heterogeneous and the stiffnesses of the different parts are known. One seeks the size of the different zones. These cases will be treated in the example sections. The various identification strategies were applied on numerical examples in which a standard direct problem

1970

P. Feissel, O. Allix / Comput. Methods Appl. Mech. Engrg. 196 (2007) 1968–1983

• uniform or Gaussian white noise sampled at each time step; • deterministic perturbations, such as high-frequency sine functions (with a different frequency for each measurement); • superposition of a white noise with a deterministic signal (e.g. non-zero-mean white noise).

ex Fig. 1. Boundary conditions uex d ðtÞ and fd ðtÞ from the reference calculation.

was first solved numerically using an explicit time scheme. This problem consisted in a rod clamped at its right hand side and subjected to a half-sine traction at the other end over a period of time of about 1.5 wavepaths. The rod’s elastic properties varied depending on the case being treated. The typical case was calculated using unit quantities to describe the problem: E0 = 1 Pa, q = 1 kg m3 and L = 1 m. One can always choose different values for these three quantities and deduce all the other quantities of the problem through dimensional considerations. Furthermore, all the examples were based on the same space and time discretizations with 100 linear spatial finite elements and 600 time steps. Since the time scheme is an explicit scheme, the Courant stability condition implies a relation between the time step and the length of the spatial elements. Here, the time step is chosen equal to one fourth of the critical time step, therefore this leads to a typical time period of T = 1.5 s. With these remarks, since the unit quantities are meaningless, some results will be presented as functions of the time step and element length. From this first calculation, the boundary fields were recovered in terms of displacements uex d ðtÞ and tractions fdex ðtÞ at both ends, as shown in Fig. 1. These can be considered to be the unperturbed measurements. In order to test the robustness of the method, perturbations in terms of both displacements and tractions were added to the boundary conditions. Thus, the measurements became: f~ d ðtÞ ¼ fdex ðtÞ þ df ðtÞ and ~ ud ðtÞ ¼ uex d ðtÞ þ duðtÞ. Several types of perturbations can be added, as illustrated in Fig. 2 for one of the traction boundary conditions:

The magnitude of the perturbation with respect to the exact signal, unless otherwise specified, was defined as the ratio of the maximum value of the noise to the maximum value of the signal. Finally, let us note that our framework was deterministic since each identification was carried out on one set of measurements.

3. Identification strategies The identification framework presented above is characterized by the transient dynamics and the corrupted measurements. In this section, we first present an identification strategy dedicated to transient dynamics and more particularly to SHPB tests [18]. Since the performance of this strategy in the case of highly perturbed measurements is not robust enough, another strategy based on the CRE principle [1] was developed. Two variations of the initial formulation of the latter strategy are also proposed. All these methods will be compared in Sections 4 and 5. In a variational inverse approach where one only seeks the model’s parameters, one evaluates the quality of the model’s parameters with respect to the experimental data using a cost function, then finds the best parameters by minimizing this cost function. Thus, the first step consists in confronting the model with the experiment, which is usually done through simulation. Then, the cost function is derived from the solution fields of this simulation. 3.1. An inverse approach for dynamic tests [18] The first method studied was inspired by [18]. It consists in defining, using the experimental data, two auxiliary problems: one with prescribed displacements and the other with prescribed tractions. These are well-posed problems

Fig. 2. Traction boundary conditions at one end for various types of perturbations. (a) Sine perturbation – 40%. (b) Uniform white noise perturbation – 20%. (c) Uniform white noise perturbation, non-zero mean (0.2) – 20%.

P. Feissel, O. Allix / Comput. Methods Appl. Mech. Engrg. 196 (2007) 1968–1983

which can be solved numerically. They yield some solution fields uKA and uDA which depend on the Young’s modulus E chosen. KA (for Kinematically Admissible) corresponds to the problem with prescribed displacements, while DA (for Dynamically Admissible) corresponds to the problem with prescribed tractions. Then, one defines a distance between the two simulations based on their solution fields. This distance is used to quantify the consistency of E with the boundary conditions. If the two calculations yield the same solution fields, this distance is equal to zero (perfect consistency). The identification consists in selecting E which minimizes the distance. Among the possible definitions of the distance, one can choose the distance based on the constitutive relation (1): Z TZ L eðuKA ; uDA Þ ¼ ðrDA  rKA ÞðDA  KA Þdx dt; ð1Þ 0

1971

Table 1 Reliable and unreliable relations in the case of 1D elasticity Reliable Equilibrium: Initial conditions:

Unreliable q€u  r;x ¼ 0 u(x, 0) = u0 _ 0Þ ¼ u_ 0 uðx;

Constitutive relation: Boundary conditions:

r = E u,x u~d and f~ d

basic problem yields some displacement, stress and boundary condition mechanical fields. In the second step, a cost function of the elastic parameters which must be identified is defined based on these solution fields and is minimized. The basic problem one must solve becomes:

0

which, in elasticity, is always positive and where (rDA, DA) are the stress and strain of the solution of the problem with prescribed tractions and, similarly (rKA, KA) correspond to the problem with prescribed displacements. Let us note that such a method can be implemented quite easily since it makes use only of standard direct calculations. However, it possesses no special mechanism to deal with perturbations on the measurements. In order to do that, a new approach is proposed in the next section.

Find the fields u(x, t), r(x, t), ud(t), fd(t), with (x, t) 2 [0, L] · [0, T], which minimize: Z Z L 1 T E1 ðr  Eu;x Þ2 dx J A ðr; u; ud ; fd Þ ¼ 2 0 0  2 L 2 L ~ þ Ajðud  ~ud Þ j0 þ Bjðfd  f d Þ j0 dt under the constraints: u 2 UAd ðud Þ;

r 2 DAd ðfd ; uÞ;

ð2Þ

3.2. Towards a robust identification strategy: formulation L

As will be shown in Section 4.1.1, the method proposed in the previous section appears to be limited with respect to experimental perturbations. There are, in fact, two questionable points in this method: there are several ways of dividing the experimental data into two sets, and the uncertain measurements are strictly prescribed in the calculation. Through our attempts to work around these difficulties, we came up with the following method. In order to confront the model with the experimental data, one has to derive a new problem from the ill-posed direct problem described in Section 2. In order to do that, one must relax some of the constraints of this direct problem. Therefore, during the identification process, our guiding principle, directly inspired by model updating techniques in vibration problems [13,11], was to concentrate on the exact verification of the properties which are considered to be reliable. Then, the uncertain relations can be taken into account by minimizing a modified constitutive relation error [15], leading to an optimal control problem. In our case, let us divide the relations into two groups, as shown in Table 1 for the case of 1D elasticity. The formulation is presented for a one-dimensional domain [0, L] over a time period [0, T], which corresponds to the problem presented in Section 2. The identification strategy proceeds in two steps. In the first step, one defines a basic problem confronting the model to the measurements by applying the CRE principle described above. This

where for all f, jf ðxÞj0 ¼ f ð0Þ þ f ðLÞ. The spaces involved in (2) are defined as follows: • UAd ðud Þ ¼ fu 2 H 1 ð0; L½Þ=ujx2f0;Lg ¼ ud and ujt¼0 ¼ _ t¼0 ¼ u_ 0 g; such a field u is called KA (Kinematiu0 ; uj cally Admissible); let us note that the initial positions and velocities are considered to be perfectly known quantities. • DAd ðfd ; uÞ ¼ fr 2 H 1 ð0; L½Þ=rnjx2f0;Lg ¼ fd ; q€ u  r;x ¼ 0 on ½0; Lg; where n(0) = 1 and n(L) = 1. Such a field r is called DA (Dynamically Admissible). A and B are two parameters which, depending on the relative confidence one has in those terms, enable one to balance the two terms defining the error (i.e. the term associated with the constitutive relation and the term associated with the measurements). One choice, which also allows the functional to be dimensionally correct, consists in dividing each term by its order of magnitude. This choice appears to be satisfactory in terms of the identification results. Once the basic problem has been defined, the solution fields are used to evaluate the quality of the model’s parameters through the introduction of a cost function. In the classic form of the method, the same functional is also used to define the cost function. This, however, is not mandatory and an interesting alternative is proposed in Section 3.3 whereby the model’s identification is performed on the model’s error term alone while the term corresponding

1972

P. Feissel, O. Allix / Comput. Methods Appl. Mech. Engrg. 196 (2007) 1968–1983

to the distance to the measurements remains in the basic problem. Hence, the identification step is the following: min J A ðrðEÞ; uðEÞ; ud ðEÞ; fd ðEÞ; EÞ; E

Find the fields u, r which minimize:

ð3Þ J B ðr; uÞ ¼

Z 0

T

1 2

Z

L

ðr  Eu;x ÞE1 ðr  Eu;x Þdx dt 0

where Functional JA is defined in (2). The following notation is also used:

under the constraints:

gA ðEÞ ¼ J A ðrðEÞ; uðEÞ; ud ðEÞ; fd ðEÞ; EÞ;

u 2 UAd ð~ud Þ;

ð4Þ

where (r(E), u(E), ud(E), fd(E)) are solutions of (2). In practice, the identification is performed in two steps as described here. First, the functional to be minimized is evaluated for a given Young’s modulus E, leading to the determination of the unknown boundary conditions and inner fields, which depend implicitly on E. Then, the minimization of gA with respect to the value of E is performed. It should be noted that the identification cost function gA (which is a functional of the Young Modulus only) is equal to the value of the functional of the basic problem JA only when its argument is taken as the solution of the basic problem.

r 2 DAd ðf~ d ; uÞ:

ð5Þ

In this case, the measurements are strongly enforced in the calculation. This basic problem would correspond to the basic problem (2) where the weighting factors A and B would tend to infinity. In the second step, one must define a cost function based on the solution fields of the basic problem. Since the first term of the functional of the basic problem is a model error term, one can use this term to estimate the quality of the model’s parameters. Thus, the identification step becomes: min gB ðEÞ ¼ min J B ðrðEÞ; uðEÞÞ; E

ð6Þ

E

3.3. Alternative identification strategies The two steps which constitute the above method can be summarized as follows: • in the first step, for given parameters, some mechanical fields are built from the confrontation of the model with the experimental data; • in the second step, a cost function is defined from these mechanical fields. This cost function is an estimate of the model’s quality whose minimization constitutes the identification. In fact, these two steps are relatively independent of one another, and it is not mandatory to use the same functional in both. Let us now address the question of what functional can be introduced in each step. In the first step, the previous method (Section 3.2) has two key features: • The whole information from the experiment and from the model leads to only a single calculation. This is possible by relaxing the constitutive relation. • A distance between the measured boundary conditions and the solution fields is introduced. This distance is expected to regularize the solution fields and filter perturbations. Since this distance involves both the displacement and the tractions, there is no excitations on the constraints (state equations) of the basic problem and the system is solicited only through the minimization of the cost function. The second feature is not mandatory in order for the basic problem to be well-posed and, therefore, one can define the basic problem as follows:

where (r(E), u(E)) are the solutions of (2) or (5). Consequently, we now have three identification strategies based on the CRE principle depending on the different terms which are introduced, see Table 2. These three methods will be compared in Section 5 through some examples. 3.4. The basic problem related to the CRE method 3.4.1. The equations to be solved This section presents the resolution of (2). The constraints of the basic problem are relaxed by assigning a Lagrange multiplier to each, which leads to the following functional: Lðr; u; ud ; fd ; u ; k; l; EÞ

Z

T L

¼ J A ðr; u; ud ; fd ; EÞ  jðu  ud Þkj0 dt 0 Z TZ L Z T L   ðq€u  r;x Þu dt  jðr:n  fd Þlj0 dt; 0

0

ð7Þ

0

where u*, k, l are the three Lagrange multipliers associated with the constraints. The system which must be verified by the solution is obtained by enforcing the stationarity of L: dL ¼ 0;

ð8Þ

Table 2 The three possible methods CRE method

CRE-A

CRE-B

CRE-C

Functional of the basic problem Identification cost function

JA gA

JB gB

JA gB

P. Feissel, O. Allix / Comput. Methods Appl. Mech. Engrg. 196 (2007) 1968–1983

where dL has the following expression:  Z T Z L L 1   L dL ¼ drE ðr  Eðu þ u Þ;x Þdx þ jdrnu j0  jdrnlj0 dt 0 0  Z T Z L   du;x ðr  Eu;x Þ þ d€ uqu dx þ jdukjL0 dt  0 0 Z TZ L du ðq€u  r;x Þdxdt  0 0 Z T   jdfd ðl þ Bðfd  f~ d ÞÞjL0 þ jdud ðk þ Aðud  ~ ud ÞÞjL0 dt þ Z0 T   ð9Þ  jdlðrn  fd ÞjL0 þ jdkðu  ud ÞjL0 dt: 0

Taking (9) into account in (8) and partially integrating the term d€ u with respect to t, the above expression leads to the following strong form of the associated equations for all t 2 [0, T]: 8 over 0; L½; r ¼ Eðu þ u Þ;x ; > > > > > q€ u  ðEðu þ u Þ;x Þ;x ¼ 0; > > > > > q€ u  ðEu;x Þ;x ¼ 0; < ð10Þ for x 2 f0; Lg; l ¼ u ; > > > k þ ððr  Eu;x ÞÞn ¼ 0; > > > > > ud Þ ¼ 0; k þ Aðud  ~ > > : l þ Bðfd  f~ d Þ ¼ 0: In the above set of equations, one can express all other quantities as functions of u and u*. Thus, the basic problem for a given Young’s modulus E becomes: Find the fields u(x,t) and u*(x,t) in H1([0, L] · [0, T]) such that 8 over 0; L½; r ¼ Eðu;x þ u;x Þ; > > > > q€ u  ðEu;x þ Eu;x Þ;x ¼ 0; > > > > < q€ u  ðEu;x Þ;x ¼ 0; ð11Þ > over f0; Lg; Eðu þ u Þn þ 1 u ¼ f~ ; > ;x d > ;x > B > > > E  > : ~d : u  u;x n ¼ u A with the initial and final conditions:  _ 0Þ ¼ u_ 0 ; over ½0; L; uðx; 0Þ ¼ u0 ; uðx; ð12Þ u ðx; T Þ ¼ 0; u_  ðx; T Þ ¼ 0: Thus, the space weak formulation associated with this problem is for any time: 8Z L   L   u >    >  ~ > ðq€ udu þ Eðu;x þ u;x Þdu;x Þdx   f d  du  ¼ 0; > > B < 0 0 Z L L > ðq€ u du þ Eu;x du;x Þdx  jAðu  ~ ud Þduj0 ¼ 0; > > > 0 > :  r ¼ Eðu;x þ u;x Þ ð13Þ with the same initial and final conditions in time for all x 2 [0, L]. The presence of initial and final conditions when solving the basic problem poses a major difficulty in relation to the time discretization, even though in dynamics time and

1973

space are coupled. Regarding the space discretization, a classical finite element approach is used. In order to deal with the initial and final conditions, several methods have been applied. Such problems are actually encountered in the framework of optimal control. The details of the methods applied to our problem can be found in [1]. There are two families of methods: global methods and local methods in time. Global methods, such as the time finite element method, add one dimension to the problem, which leads to the resolution of a large system. In local approaches, one cannot solve the problem taking into account both the initial and final conditions simultaneously. One has therefore to deal with the unstable nature of the problem when it is not controlled by its final conditions. As shown in [2], the methods which appear to be the most effective are based on the Riccati approach. Here, for efficiency reasons a method which consists of guessing the missing initial conditions thanks to the transition matrix associated with the problem was used and appears to be only relevant for linear problems and short-duration tests, which is the framework we placed ourselves here. Besides, it is interesting to note that, when the same functional is used in the definition of the basic problem and in the identification step, the cost function and its gradient can be deduced from the solution fields of the basic problem with no additional calculation. The gradient is given by DgA ðEÞdE ¼

oL dE oE

where L was introduced in ð7Þ:

ð14Þ

3.4.2. Solution fields of the basic problem The first step of the method consists in calculating the solution fields of the basic problem (2), which will be used to estimate the cost function of the material parameters. We will present these fields for various cases depending on the value of the Young’s modulus and the perturbations on the measurements. The measurements were created as presented in Section 2, from a first direct calculation with a given Young’s modulus E0. Fig. 3 shows the displacement and stress fields deduced from the calculation providing the measurements. In the x · t plane, one can observe the stress wave propagating from the boundary x = 0 and reflecting on the clamped end at x = L. In the case of exact measurements and for a Young’s modulus equal to the reference value E0, the solution of the basic problem corresponds to the experimental fields (Fig. 3). Indeed, these fields satisfy the constraints exactly and the error is zero, therefore minimum. Furthermore, since r = E(u + u*), u* is identically equal to zero. In the case of exact measurements, but with a Young’s modulus twice the reference value E0, the boundary conditions are no longer compatible with the constitutive relation and the change in wave velocity leads to incompatibilities along the wave paths. The solution fields u and u*,

1974

P. Feissel, O. Allix / Comput. Methods Appl. Mech. Engrg. 196 (2007) 1968–1983

Fig. 3. Displacement u and stress r in the space · time plane – reference calculation.

Fig. 4. Solution fields in the space · time plane, E = 2E0.

shown in Fig. 4, do not verify either the constitutive relation or the boundary conditions. From the solution fields, it is possible to study the map of the error contribution, shown in Fig. 5. The field E(u*)2 represents the density of the error in constitutive relation, from which one can extract the following contributions:

• Instantaneous error contribution of Element k at time step ti: Z tiþ1Z 1 Eðu Þ2 dx dt: ð15Þ Eti ;k ¼ Dtle ti k • Error contribution of Element k integrated over time: 1 Ek ¼ Tle

Z TZ 0

2

Eðu Þ dx dt:

ð16Þ

k

where Dt is the time step and le the length of the spatial elements.

Fig. 5. Map of the error contribution Eti ;k in the space · time plane.

Fig. 5 shows the field Eti ;k in the space · time plane in the case E = 2E0 with no perturbation on the measurements. The constitutive relation appears to have been relaxed along the wave path due to the incompatible boundary conditions. Fig. 6 shows the error contribution Ek of each element for various types of perturbations. One can observe that for E = E0 and a white noise perturbation of magnitude 10% the error contribution is localized near the ends of the bar. In Section 5.2.1, we will show that this is related to the

P. Feissel, O. Allix / Comput. Methods Appl. Mech. Engrg. 196 (2007) 1968–1983

exact measurements, E = 2E0 perturbed measurements (sine), E = 2E0 perturbed measurements (sine), E = E0 perturbed measurements (white noise), E = E0 3 2.5

εk

2 1.5 1 0.5 0 0

10

20

30

40

50

60

70

80

90

100

element number

1975

4.1.1. Identification based on the standard method [18] This identification was first carried out using the method presented in Section 3.1 which is our reference for identification in transient dynamics (and SHPB tests in particular). Fig. 7 shows the cost function to be minimized with respect to the relative Young’s modulus (E/E0). In the absence of perturbations, the minimum is obtained for E = E0 and the method leads to proper identification. However, for large perturbations, the curve has no minimum in the interval being considered and would yield highly erroneous identification results. This is due to the fact that the two calculations are very different, not because of the material parameters, but because of the perturbations on the boundary conditions which introduce artificial inconsistencies between the Young’s modulus and the boundary conditions.

Fig. 6. Error contribution Ek of each element integrated over time.

frequency range of the perturbations. Furthermore, the solution fields differ from the experimental fields, i.e. the error is non-zero. For E = 2E0, we compare the error contribution Ek in the case of unperturbed and perturbed measurements. One can see that the error contribution due to the perturbations is small compared to the modeling errors. Even if, in a qualitative sense, as shown in Fig. 6, the total error can be viewed as the superposition of the error due to the Young’s modulus and the error due to the perturbed measurement, this superposition is not strictly additive: the curve for perturbed measurements and a wrong parameter is not the sum of the curve with perturbed measurements only and the curve with a wrong parameter only. This is due to the fact that the solution fields are not linear with respect to the Young’s modulus.

4.1.2. Identification based on the CRE-A method The same example was treated with the CRE-A method. The basic problem was solved for a range of material parameters, which enabled us to estimate the cost function (which depends on only one parameter). Then, the minimum of this cost function could be easily found, leading to the identified value of E. Fig. 8(a) shows the identification curves for the same cases as for the previous method (Fig. 7). The identified Young’s modulus corresponds to the minimum of each curve. In all three cases, the identification achieved by the method was good (less than 3% error) and the sensitivity to moderate perturbations very small. The method appears to be very robust with respect to perturbations of the measurements.

4. Identification performance of the CRE-A method

4.2. Example of a heterogeneous 1D bar

In this section, we first compare the CRE-A method, as defined in Table 2, which was the standard CRE method used in model updating until now [11], with the method from [18] on an example with various levels of perturbations. Then, we apply the CRE-A method to more complex examples in which the rod is assumed to be heterogeneous, so that the identification is carried out with several parameters.

The second example deals with the identification of the material properties of a heterogeneous elastic bar subjected to dynamic loading. The identification framework remains

4.1. Example of a homogeneous 1D rod – comparison with [18] This is the same example as the one we described in Section 2, in which we assume the rod to be homogeneous. Thus, the objective is to identify only a single parameter: the Young’s modulus E of the rod. The measurements were created by performing an initial calculation with a reference value E0 and adding perturbations in the form of high-frequency sines of various magnitudes (but the conclusions would be the same for other types of perturbations).

Fig. 7. The cost function at various perturbation levels (based on [18]).

1976

P. Feissel, O. Allix / Comput. Methods Appl. Mech. Engrg. 196 (2007) 1968–1983

Fig. 8. The cost function for various levels of perturbations. (a) Total error, gA(E). (b) Constitutive relation error term, gB(E).

the same as in the previous example. Two cases were treated: • Case A: The spatial distribution of the properties is known and only the values of these properties are sought. This radical assumption eases the identification process on a reduced number of parameters. • Case B: No assumption is made on the distribution of the properties; therefore, a field E(x) is sought. This can be an ill-posed inverse problem. 4.2.1. Identification procedure In both cases, the identification strategy was the CRE-A method proposed in Section 3.2 with two terms in the cost function. For the identification step, one has to minimize the modified constitutive relation error with respect to the set of parameters to be identified. If one denotes E(x) the field of the parameters to be identified, depending on the case being treated, E is piecewise constant or function of the space variable. In order to perform this minimization, we first used the simple algorithm presented in Fig. 9 (constant-step gradient method), based on the direct estimation (14) of the gradient from the solution fields of the basic problem. Although inefficient, this first algorithm led to the solution. Now that the capability of the method has been proven, work is in progress in order to improve the algorithm by using an optimum a(n) coupled with a BFGS method. The procedure is described in Fig. 9. With the cost function of the one-dimensional case, the gradient is oL oJ ðdEÞ ¼ ðdEÞ oE oE   Z TZ L o E1 2 ¼ ðr  E  u;x Þ dE dx dt 2 0 0 oE Z TZ L  u;x  ðu þ 2u;x ÞdE dx dt: ¼ 2 ;x 0 0

DgA ðEÞðdEÞ ¼

Fig. 9. The identification algorithm.

4.2.2. Identification: Case A When the spatial distribution of the properties is known, one must calculate the gradient in each zone, given by X Z TZ ðu;x Þ ðu;x þ 2u;x Þdx dt dEi DgA ðEÞðdEÞ ¼ 2 Xi 0 i |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} ¼

X

rgi ðEÞ

rgi ðEÞdEi ;

ð18Þ

i

where Xi = {x 2 [0, L]/E(x) = Ei} with Ei constant. The method was then applied to the example of a beam consisting of four parts of the same size, but with different Young’s moduli, as shown in Fig. 10. The experimental measurements were created by an initial direct calculation of this bar, as presented in Section 2, clamped at the E4-

ð17Þ Fig. 10. Piecewise homogeneous elastic rod.

P. Feissel, O. Allix / Comput. Methods Appl. Mech. Engrg. 196 (2007) 1968–1983

4.2.3. Identification: Case B The method was then applied to the identification of a property field E(x) with no assumption on its spatial distribution. Therefore, the gradient is given by Z T ðu;x Þ  rgA ðEÞðxÞ ¼ ð19Þ ðu;x þ 2u;x Þdt: 2 0 In order to compute the identification of E(x), we decided to discretize it as a piecewise constant field and we choose it to be constant in each space element. This choice actually simplifies the problem (which corresponds now to the case A), but the great number of parts raises some of the difficulties encountered in the case of the identification of a continuous field E(x). The loading was the same and the minimization proceeded in the same way after being initialized with a homogeneous field such that the mean velocity in the bar be the same as in the exact bar. The property field at various iterations is shown in Fig. 13. By the end of the identification process, the identified distribution was not exact, but the trend was good. Recon-

3.5 3

E

4

Young's modulus

2.5 2 E2

1.5 1

E1

0.5

E3

0 0

200

400

600

800

1000

1200

1400

iterations

Fig. 12. Identification by zones, 5% perturbed measurements.

3

iteration 1 iteration 10 iteration 50 iteration 500 iteration 1800 reference

2.5

Young's modulus

modulus end with a half-sine traction applied at the other end. The results of the identification are presented in Fig. 11. Fig. 11(a) shows the values of the four Young’s moduli throughout the iterations. The initialization was done with a constant Young’s modulus such that the mean wave velocity in the bar be the same as the exact value. Such a value was determined by identification under the assumption that the rod was homogeneous. There are two convergence phases corresponding to two different gradient steps. As can be seen in Fig. 11(b) showing the descent in the E3 · E4 subspace, the descent step was reduced because due to the step being too large the convergence rate was becoming too small. The method led to the correct values of the parameters. We treated a second example in which the measurements were perturbed by a uniform white noise sample with a magnitude of 5%. At the end of the identification process, the errors in the Young’s moduli were small, as can be seen in Fig. 12.

1977

2

1.5

1

0.5

0

0

10

20

30

40

50

60

70

80

90

100

element number

Fig. 13. Modulus field at various time steps, unperturbed measurements.

structing the complementary boundary conditions through a direct simulation of the bar subjected to one of the measurements at each end and with the identified material properties, we obtained boundary fields very close to the

Fig. 11. Identification by constant zones, unperturbed measurements. (a) Young’s moduli of the four zones throughout the iterations. (b) Descent in the E3 · E4 space.

1978

P. Feissel, O. Allix / Comput. Methods Appl. Mech. Engrg. 196 (2007) 1968–1983

5.1. Taking into account the whole information alone: the CRE-B method

force boundary conditions 2

force

1 0 –1 –2 –3

0

measurements, x = 0 measurements, x = L x = 0, iteration 1800 x = L, iteration 1800 x = 0, iteration 1 x = L, iteration 1 0.5

time (s)

1

1.5

displacement boundary conditions 0.4

displacement

0.3 0.2 0.1 0

–0.1

0

0.5

1

1.5

time (s)

Fig. 14. Modulus field identification, reconstructed boundary conditions.

other two measurements (see Fig. 14). Therefore, one may consider that convergence had been reached. This means that the identification problem has several solutions, or at least that the inverse problem is ill-posed due to the discontinuity of E(x) with respect to the measurements. In order to deal with such a situation, one can either enrich the experimental data, for example by extending the observation time, or add some a priori information on the property field as we did in Case A. These initial examples showed that the identification strategy is robust with respect to perturbations of the measurements and seems usable for complex problems. Furthermore, the resolution of the basic problem gave a direct estimate of the cost function and of its gradient, which was used in the heterogeneous case. 5. Understanding and improving the proposed method The objective of this section is to identify the key points of the method and to improve its effectiveness, given that the two main features of the method are that it takes into account the whole information in one calculation and that it introduces a distance between the measurements and the calculation to define the mechanical fields. In order to do that, we will discuss the performances of the various CRE methods (CRE-A, CRE-B and CRE-C) on several examples. First, by looking at the CRE-B method, we will see that taking into account the whole information in one calculation improves the accuracy, but not sufficiently. Then, we will illustrate the role of the distance to the measurements on some examples, which will lead, in particular, to an improvement of the CRE-A method in the case where one has some a priori knowledge of the perturbations. Finally, we will show the effectiveness of the CRE-C method when this knowledge is lacking.

In order to show the role played by the distance-to-measurements term, let us study the CRE-B method. In this case, there is no distance-to-measurements term and, therefore, the measurements are strictly prescribed in the basic problem. 5.1.1. Pollution of the modeling error Fig. 15 shows the contribution of each element to the error over time in various cases, for perturbations created with a high-frequency sine function with a magnitude 20% of the maximum value of each perturbed signal. In the case with the exact Young’s modulus, the error lies primarily in the vicinity of the boundary conditions, where the perturbations are. However, with a stiffer Young’s modulus, the error field is more perturbed than in the case of the modified constitutive relation error (the CRE-A method). We can conclude that the CRE-B method results in a greater pollution of the inner fields by the boundary perturbations. 5.1.2. Cost function The consequence of such behavior can be seen in the effect of perturbations on the cost function. This function is obtained through the solution fields of (5) and is denoted gB as defined by (6). The identification curves are presented in Fig. 16 for various levels of perturbations. For small perturbations, the method remains accurate, but when the level of perturbations becomes high the cost function is greatly affected by the perturbations and its mean value is about ten times higher than in the absence of perturbations. The minimum is less clearly defined, even though the method still makes an identification possible, which was not the case of the first method tested in Section 4.1.1. This example shows that taking into account all experimental data in a single calculation is a key point of

Fig. 15. Contribution of each element to the error over time.

P. Feissel, O. Allix / Comput. Methods Appl. Mech. Engrg. 196 (2007) 1968–1983

1979

Fig. 16. Cost function gB for CRE-B method as a function of the Young’s modulus.

the method, and also that the introduction of a distance between the calculation and the measurements improves the effectiveness of the identification. 5.2. Distance to the measurements and choice of the weighting factor The previous example pointed out that the distance-tomeasurements term improves the robustness of the identification strategy. Here, we focus on the effects of this term on the cost function and on the solution fields of the basic problem. 5.2.1. Error terms The second term of the modified constitutive relation error was introduced in order to deal with discrepancies in the measurements. Its value is expected to reflect the quality of the measurements, whereas the first term ought to give an evaluation of the quality of the model alone. Since the decomposition of the error is not that rigorous, one can discuss the pollution in the modeling error term due to the perturbations (and also the pollution in the distance-to-measurements term due to the modeling error). This pollution is expected to be relatively small. Fig. 17 shows the two terms of the modified error in the case of white noise at various levels of perturbation. One can see that the perturbations have very

Fig. 18. Spatial distribution of the error Ek , effect of the frequency, E = E0.

little effect on the modeling error term, as was stated in Fig. 6 showing the effectiveness of the measurement error term. Nevertheless, Fig. 18 shows the effect of the frequency range of the perturbation on the pollution in the modeling error term. More precisely, it shows the contribution of each element to the error Ek for sine perturbations with a magnitude of 20% of the maximum value at various frequencies. The higher the frequency, the smaller the pollution level. Altogether, the pollution level remains reasonable. These remarks explain the discrepancy between the pollution in the cost function observed for a high-frequency sine perturbation and for a white noise (Figs. 8(b) and 17). 5.2.2. Boundary conditions The comparison of the boundary conditions between the simulation and the measurements shows the regularizing effect of the second term of the modified constitutive relation error. Fig. 19 shows the measured boundary conditions and those deduced from a calculation with a zero-mean white

Fig. 17. Cost function – perturbation: uniform white noise sample. (a) Constitutive relation error term. (b) Distance to the displacements.

1980

P. Feissel, O. Allix / Comput. Methods Appl. Mech. Engrg. 196 (2007) 1968–1983

Fig. 19. Displacement: measurements and boundary conditions, white noise sample with zero mean.

noise perturbation for E = E0. The boundary fields were indeed regularized. In the case of a non-zero-mean perturbation as illustrated in Fig. 20, the calculated boundary conditions are smoother, but they also differ from the measurements in their mean values. In Fig. 20, these are compared with the exact boundary conditions (reference simulation before perturbation). This actually means that the reconstructed fields are not fully satisfactory. Despites this drawback, the identification based on the cost function, which is the key point, leads to good results. 5.2.3. Choice of the weighting factor The distance-to-measurements term can be viewed as a regularization term, even though a zero weight on this term would lead to a trivial solution of the basic problem and, therefore, would prevent any identification. Thus, the choice of the relative weights of the two terms depends on the reliability of the measurements, a higher weight yielding a better fit to the measurements. The first step in the choice of the weighting factor consists in choosing a factor such that the two terms of the error have the same order of magnitude [6]. This is done by dividing each term by the norm used in the term applied on the measured quantity (the constitutive relation error term is divided by the strain energy of the structure integrated over time).

Fig. 21. Distribution of the identified E (1200 samples of Gaussian white noise with a 20% standard deviation): choice of the weighting factor.

Table 3 Statistical data on the identified Young’s modulus Weighting factor

1

1/10

Mean value Standard deviation

0.987 2.6 · 102

0.998 1.3 · 102

Based on this initial weighting factor, one can choose a more accurate weighting factor when there is some a priori knowledge of the quality of the model as well as of the measurements. Some techniques based on the Tikhonov regularization scheme can then be applied [7,5]. In [6], the choice was based on a probabilistic approach and the study led to the conclusion that for model updating in vibration problems the initial non-dimensional weighting factor yields good results. In order to illustrate this point, a study based on a Monte Carlo approach for white Gaussian noise perturbations with a 20% standard deviation was conducted. For each noise sample, an identification is conducted and its result stored. The stochastic output of the identification process is therefore a distribution of identified Young’s modulus and is shown in Fig. 21 using this white noise for two weighting factors: the non-dimensional factor and the same divided by ten. Table 3 gives the characteristics of the identified Young’s modulus distribution and emphasizes the improvement due to a good weighting choice. In both cases, the identification remained accurate. However, in our context, there is no a priori knowledge of the perturbations, and the non-dimensional weighting factor is the more reasonable choice. 5.3. CRE-C: an alternative method in the absence of knowledge of the perturbations

Fig. 20. Displacement: measurements and boundary conditions, white noise perturbation with non-zero mean.

In the CRE-C method, the basic problem is defined with the modified constitutive relation error, that is with the distance to the measurements term, and the cost function is based on the modeling error term (i.e. the constitutive relation error) alone. The objective of this section is to study this mixed method (CRE-C) and compare it with the first method (CRE-A) through two examples.

P. Feissel, O. Allix / Comput. Methods Appl. Mech. Engrg. 196 (2007) 1968–1983

5.3.1. Comparison of the stochastic responses of the methods The first comparison of the CRE-A and CRE-C methods was performed by comparing the results of the two methods (i.e. the distributions of the identified Young’s moduli obtained by a Monte Carlo technique) under a Gaussian white noise perturbation on the measurements. The identification methods were applied to the same clamped homogeneous bar problem as in the previous examples, in which the perturbations were created in the form of a Gaussian white noise with a standard deviation equal to 20% of the maximum value of the signal for each boundary condition and for each time step, see Fig. 22. Since this is a large perturbation, the question of the weighting of the terms of the functional was raised as discussed in Section 5.2.3. Therefore, for each perturbation sample, the identification was performed with both the CRE-A and CRE-C methods and with two different weighting factors, the non-dimensional factor and onetenth of the same. Fig. 23(a) and (b) show the distribution of the identified Young’s moduli for 1200 perturbation samples. Table 4 gives the mean values and standard deviations

Fig. 22. Perturbed displacement boundary conditions.

1981

Table 4 Statistical data on the identified Young’s moduli Weighting factor

1

Method Mean value Standard deviation

CRE-A 0.987 2.6 · 102

1/10 CRE-C 0.997 8.6 · 103

CRE-A 0.998 1.3 · 102

CRE-C 0.999 9.6 · 103

of these distributions, which confirm the robustness of both methods, since the mean values are close to the reference value and, in the worst case, the standard deviation is 2.8%. Furthermore, the comparison of the two methods leads to the conclusion that the mixed method yields better results, even though a good choice of the weighting factor improves the CRE-A methods. The CRE-C method appears to be a good option in the absence of a priori knowledge of the perturbations because it is affected only minimally by the weighting factor in this example. 5.3.2. Identification of the extent of a heterogeneity The second example on which we propose to compare the two methods concerns an elastic bar with Young’s modulus E0 containing a heterogeneous section with modulus E0/2 in the middle. The objective is to identify the extent of the heterogeneity. The entire bar was meshed with 100 elements. An initial calculation was performed in order to create the measurements for various extents of the heterogeneity (between 1 and 20 elements). The boundary fields with or without added perturbations were used as the measurements for the identification process. Both methods were applied with the various heterogeneity extents. Fig. 24 shows the identified extent of the heterogeneity as a function of the extent used to create the measurements. Without perturbations, both methods yield the right extent in all cases. When the measurements were perturbed by a 15% white noise sample, there were errors in the identification results in some cases, but the CRE-C methods gave much better results.

Fig. 23. Distribution of identified E (1200 samples of Gaussian white noise with a 20% standard deviation). (a) Methods A and C with non-dimensional weighting. (b) Methods A and C with 1/10 non-dimensional weighting.

1982

P. Feissel, O. Allix / Comput. Methods Appl. Mech. Engrg. 196 (2007) 1968–1983

Fig. 24. Identification of the extent of the heterogeneity.

5.3.3. Identification strategy regarding the knowledge of the perturbations These two examples confirmed the robustness of the method and, in the case of the CRE-A methods, emphasized the importance of the weighting factor. They also showed that the CRE-C method seems to be more accurate and should be effective in the absence of a priori information on the perturbations, even though the change of cost function leads to the lack of a direct evaluation of its gradient from the solution fields of the basic problem. These results will have to be confirmed on more examples.

6. Conclusion We sought a solution to the problem of the identification of material parameters from dynamic tests with high measurement uncertainties. The method we proposed is based on the use of the modified constitutive relation error for identification purposes in the framework of transient dynamics. By adhering strictly to the basic idea which consists of partitioning relations into a reliable set and an uncertain set, we ended up with an optimal control problem leading to a two points boundary value problem requiring the use of suitable methods, involving a rather large computation time. Despite this drawback, the approach has the great advantage of being robust with respect to measurement uncertainties, at least on the examples we considered. The study presented here enabled us to show that the robustness of the method relies mainly on two principles: taking into account the whole experimental data in a single calculation and introducing a distance between the measurements and the calculation. This distance-to-measurements term is necessary in order to regularize the boundary conditions and the mechanical fields, but useless for the identification itself. Therefore, we proposed an alternative method which consists in introducing this distance in the functional of the basic problem alone, then carrying out the identification with a cost function based on the modeling error term. This approach seems particularly suitable in the absence of a priori knowledge of the perturbation level. The next steps of our work will

focus on the adaptation of this method to the case of damage, then to more complex structures. References [1] O. Allix, P. Feissel, H.M. Nguyen, Identification strategy in the presence of corrupted measurements, Engrg. Comput. 22 (5/6) (2005) 487–504. [2] O. Allix, P. Feissel, A delay damage mesomodel of laminates under dynamic loading: basic aspects and identification issues, Comput. Struct. 81 (12) (2003) 1177–1191. [3] O. Allix, H.M. Nguyen, P. Feissel, An identification strategy based on the CRE in the presence of highly scattered measurements and its comparison with an approach based on the extended Kalman filter, in: 8th European Mechanics of Materials Conference, Cachan, France, 2005. [4] S. Andrieux, F. Voldoire, Stress identification in steam generator tubes from profile measurements, Nucl. Engrg. Des. 158 (2–3) (1995) 417–427. [5] H.W. Engl, M. Hanke, A. Neubauer, Regularization of Inverse Problems, Kluwer Academic Publishers, 1996. [6] D. Barthe, A. Deraemaeker, P. Ladeve`ze, Validation of dynamics models with uncertainties based on the constitutive relation error, in: Complas VII – 7th International Conference on Computational Plasticity, Barcelona, Spain, CIMNE, 2003, pp. cd-rom. [7] H.D. Bui, Inverse Problems in the Mechanics of Materials: An Introduction, CRC Press, 1994. [8] A. Cimetie`re, F. Delvare, F. Pons, Une me´thode inverse a` re´gularisation e´vanescente, C. R. Acad. Sci. – Series IIB – Mech. 328 (9) (2000) 639–644. [9] A. Corigliano, A. Ghisiy, S. Mariani, Parameter identification of nonlinear constitutive laws by Unscented Kalman Filters, in: E. On˜ate D.R.J. Owen (Eds.), Complas VIII – 8th International Conference on Computational Plasticity, Barcelona, Spain, 2005. [10] A. Corigliano, S. Mariani, Parameter identification in explicit structural dynamics: performance of the extended kalman filter, Comput. Methods Appl. Mech. Engrg. 193 (2004) 3807–3835. [11] A. Deraemaeker, P. Ladeve`ze, Ph. Leconte, Reduced bases for model updating in structural dynamics based on constitutive relation error, Comput. Methods Appl. Mech. Engrg. 191 (21–22) (2002) 2427–2444. [12] R.E. Kalman, A new approach to linear filtering and prediction problems, Trans. ASME J. Basic Engrg. 82 (1960) 35–45. [13] P. Ladeve`ze, A. Chouaki, Application of a posteriori error estimation for structural model updating, Inverse Problems 15 (1999) 49– 58. [14] P. Ladeve`ze, M. Reynier, N.M.M. Maia, Error on the constitutive relation in dynamics, in: M. Tanaka, al. H.D. Bui (Eds.), Inverse problems in Engineering, 1994, pp. 251–256.

P. Feissel, O. Allix / Comput. Methods Appl. Mech. Engrg. 196 (2007) 1968–1983 [15] P. Ladeve`ze, A modelling error estimator for dynamical structural model updating, in: P. Ladeve`ze, J.T. Oden (Eds.), Advances in Adaptive Computational Methods in Mechanics, Elsevier, 1998. [16] L. Ljung, Asymptotic behavior of the extended Kalman filter as a parameter estimator for linear systems, IEEE Trans. Automat. Control AC 24 (1979) 36–50. [17] J. Orkisz, A. Skrzat, Reconstruction of residual stresses in railroad vehicle wheels based on enhanced saw cut measurements, formulation and benchmark tests, Wear 191 (1–2) (1996) 188–198.

1983

[18] L. Rota, An inverse approach for identification of dynamic constitutive equations, in: Tanaka Bui (Ed.), ISIP 94: Second International Symposium on Inverse Problems in Engineering Mechanics, Paris, 2–4 November 1994, Balkema, 1994, pp. 157–164. [19] A.N. Tikhonov, V.Y. Arsenin, in: Solutions of Ill-Posed Problems, Winston and Sons, 1977. [20] E.A. Wan, R. van der Merwe, The unscented kalman filter, in: S. Haykin (Ed.), Kalman Filtering and Neural Networks, John Wiley and Sons, Inc., New York, 2001, pp. 221–280.