Parameter identification of two-time-scale nonlinear transient models

The aim of this paper is to study two-time-scale nonlinear transient models and their associated parameter identification. When it is possible to consider two well-separated time scales, and when the fast component of the applied loading is periodic, a periodic time homogenization scheme, similar to what exists in space homogenization, can be used to derive a homogenized model. A parameter identification process for this latter is then proposed, and applied to an academic example, which allows to show the benefits of such a strategy.

In order to give accurate predictions, such time-homogenized models have to be compared with experimental data. The key point is to define an identification strategy able to deal with such models and using a process that remains cheap and efficient. The aim of this paper is to analyze on a specific academic example how a relevant identification strategy can be proposed, and to comment on the different choices that can be made throughout the whole identification process. First, the basic elements of the periodic time homogenization method are described, as well as one particular case of study showing the efficiency of the method; then, the general parameter identification strategy that will be adapted to the case of time-homogenized models is introduced. After this introduction, the main part of the paper focuses on a simple academic example allowing to give the theoretical developments associated with a truly two-time-scale identification strategy, and to analyze the associated numerical results upon which conclusions can be drawn.

Main ingredients
Periodic time homogenization, as it was initially proposed in Ref. [2], can be seen as a transposition to time of the classical periodic space homogenization methods, such as the techniques described in Ref. [3] or in Ref. [4]: it is typically based on the introduction of a two-scale asymptotic expansion for every time-dependent variable sought for in the equations of the reference problem: where t stands for the so-called slow time scale and τ = t/ξ for the fast time scale. Indeed, when the ratio ξ is very small, it is possible to separate the two time scales by considering that they are independent one from the other, implying that any derivative with respect to time has to use the partial derivatives with respect to the two time scales: where d t •, ∂ t • =• and ∂ τ • = • ′ stand for the total time derivative, the partial derivative with respect to the slow time and the partial derivative with respect to the fast time respectively. Moreover, if the applied loading has a component, which is periodic with respect to the fast time τ, it is possible to assume that each variable is quasiperiodic, meaning that it is periodic with respect to τ . Figure 1 gives an illustration of such a quantity. This quasiperiodicity assumption is all the more justified as the ratio ξ is small, and can then be written as: where T f is the period associated with the fast periodic component of the loading.
Using the asymptotic expansion (1) for each quantity in the reference equations, and balancing the different orders of ξ, the time-homogenized equations are determined by averaging over a fast period the different quantities: (1) α(t, τ ) = α 0 (t, τ ) + ξα 1 (t, τ ) + ξ 2 α 2 (t, τ ) This fast-time average actually allows to separate slow-evolving phenomena and fasttime periodic components by using the fact that the quasiperiodicity assumption can actually be rewritten for any time-dependent variable as: In addition to the different time-homogenized quantities α , the residuals associated with this fast-time average are then denoted as α * = α − �α�, and depend on both time scales t and τ a priori. They usually have to verify a very simple problem (typically a linear, elastic one), which will be solved for a time interval corresponding to one fast period T f . Eventually, the obtained time-homogenized equations are solved relatively to the slow time scale only, by introducing the averaged influence of the fast cycles corresponding to the solution of the fast problem at each slow time step, hence allowing to solve for all the zeroth-order time-homogenized variables with a drastic reduction in the computational cost when compared with what would be required to solve the reference problem (for example by using the classical rule of thumb of 20 time steps per fast cycle).

Cases of study
References on the time homogenization method still tend to be quite scarce, as showed a recent tentative review in Ref. [5]. Nevertheless, the possible applications are numerous, as we can see with a significant number of references dealing with topics as diverse as the behavior of viscoelastic materials in Ref. [6], the vibration of preloaded beams in Ref. [7], the ultrasonic welding of composites in Ref. [8], the ultrasonic imaging of tissues in Ref. [9], the cyclic loading of normally consolidated clay in Ref. [10] or the viscoelastic-viscoplastic behavior of polymers in Ref. [11].
Recently, we focused on validating the method for different cases of simulations of structures withstanding fatigue loads with two periodic components: • material fatigue with a viscoplastic law defining two hardening variables in Ref. [12]; • material fatigue with an isotropic damage law in Ref. [13]; • extension of the method to three different time scales in Refs. [5] and [14]; • extension of the method to the case of a resonant excitation in Ref. [5]. In this latter case, a real case example, coming from the European Project PREMECCY on fatigue, related to gas turbine blades, and briefly described in Ref. [15], was presented: indeed, these latter, during operation, withstand a loading with two periodic components whose time evolution is sketched in Figure 2. The PREMECCY project's principle was to experimentally study CCF with several blade-shaped specimens, such as the one depicted in Figure 3. Numerical simulations using the periodic time homogenization method have been run on such a geometry, with a loading implying frequencies  [15].

Figure 3
Numerical results for a specific specimen. Bending II mode at 1,400 Hz and final local longitudinal plastic strain after 100 slow cycles and 1,000,000 fast cycles. F = 0.14 Hz and F /ξ = 1,400 Hz, this latter value corresponding to the specimen's second bending eigenfrequency. A specific result is shown in Figure 3, where it is easily seen that the highest values of the longitudinal plastic strain are located at the specimen's edges, and that these values nearly vanish everywhere else. More details can be found in Ref. [5].

Gradient-based parameter identification using adjoint state formulations
To improve the quality of the model's predictions, especially the fatigue life estimate, it is mandatory to compare computed quantities with experimental data and to define an identification scheme allowing to update the model's parameters judiciously. In this section, a general description of the proposed identification strategy is given, which, as we will see, is relevant for both reference and time-homogenized problems: it will thus be illustrated in this latter case in the "Methods".
The forward problem is considered as an implicit formulation with a vector function F over a time interval [0, T ] : where d t and d 2 t are the first-and second-order time derivatives respectively. U 0 and V 0 stand for the initial conditions of the dynamic problem. Whereas u is the state vector of size N, composed of all the time-dependent degrees of freedom (DOFs) describing the studied problem, p stands for the vector containing the P scalar parameters associated with the differential equation (6).

Formulation of the identification problem
The identification problem consists in finding the parameter vector p opt such that the solution u(t; p opt ) of (6) obtained with the parameters p opt is as close to the available experimental data as possible. These latter are indeed compared with the corresponding quantities Au(t; p), where A is a projection operator allowing to select, for each quantity, the closest DOF to the experimental measurement point. In order to use consistent notations, the corresponding experimental quantity is denoted Au exp (t); however, it does not mean that such a vector u exp (t) actually exists.
The following misfit function is then introduced: it consists of a norm measuring the discrepancy between the quantities predicted with the forward model (6) and experimental data: where u(t; p) satisfies Eq. (6). The L 2 -norm proposed here is completed with a Tikhonov regularization term, allowing to deal with the ill-posedness of the identification problem, by bounding the magnitude of the parameter vector p to be identified: this regularization term uses a vector p 0 containing nominal values corresponding to a priori experience, and a diagonal weighing matrix R. Eventually, the solution of the identification problem can be sought as the parameter vector p opt minimizing the misfit function J (p):

Adjoint state formulation
The determination of this minimum is achieved using gradient-based minimization methods, therefore the question of avoiding local minima by means of an appropriate regularization process should be carefully addressed. In some cases, rather than using the classical Tikhonov regularization term, the a priori experience can be introduced in some specific ways, as in Ref. [16]. Similarly, the fact of using a homogenized model in the parameter identification process can introduce a regularizing effect, just as explained in Ref. [17]. However, we will not address here this specific question, but rather focus on the identification process itself.
To estimate the gradient of the misfit function, we solve here an adjoint state problem. A typical example in mechanical engineering is given in Ref. [18], where the parameters of an elastoplastic material law are identified with indentation tests. In the strategy proposed here, the generic form of the adjoint state problem, which can be obtained by expressing the stationarity of a Lagrangian as in Ref. [19], is as follows: where ∇ u F , ∇ d t u F and ∇ d 2 t u F stand for the directional derivatives of F with respect to u, d t u and d 2 t u respectively. The adjoint state problem is then a time-backward differential equation with two final conditions, and where the first-order sensitivities of the forward problem are concerned.
Once the adjoint state problem (9) is solved, it can be shown that the misfit function's gradient with respect to the parameter vector p can be expressed as: This specific way of estimating the misfit function's gradient can be compared with a classical finite difference formula, such as the central finite difference scheme: in this latter case, when the parameter vector is of size P, the gradient calculation is obtained by evaluating the misfit function in 2P additional 'points' , each couple of points corresponding to two symmetrical perturbations of the misfit function associated with each parameter in the vector p. The resulting computational cost for each gradient evaluation consists of 2P solutions of the forward problem (6) and 2P time integral evaluations. By contrast, when the adjoint state solution is used, only two differential equation solutions are required: one for the forward problem (6) and one for the adjoint problem (9). The associated computational cost for each gradient evaluation is then two differential equation solutions, and P time integral evaluations. The resulting gain is as high as the number of parameters to be identified. Moreover, it is easier to control the accuracy of the gradient estimate with the adjoint state method than with finite difference formulas, for which the choice of the discretization steps has a strong influence on the final estimate.

Methods
We propose to study an academic example in order to discuss the different steps of the parameter identification problem associated with a time-homogenized model. To the best of our knowledge, it is the first time that such a two-time-scale identification strategy using an adjoint state formulation is proposed: only one alternative has been briefly studied in Ref. [20] with the use of a genetic method. First, we will recall how the identification process is formulated on a classical reference problem, i.e. with no time homogenization. Then, we will apply the identification strategy on the time-homogenized version of the problem, and see what connections exist with the identification process associated with the reference problem.
The case of study detailed here consists of a straight bar of length L, clamped at one end, and withstanding at the other end a normal force with two periodic components. The measured displacement at this end is then used to determine the parameter values of the material elastic viscoplastic law.

Forward problem
The dynamic equilibrium of the bar is a scalar equation, with the normal (nonviscous) stress σ (x, t) and the longitudinal displacement u(x, t) defined at each point x ∈ [0, L] and for any t ∈ [0, T ]: where ∂ x is the partial space derivative, ρ is the mass density, and c K is the damping ratio where damping is assumed as proportional to stiffness. Homogeneous initial conditions for the displacement and the velocity are assumed. No load is applied along the bar, and, whereas this latter is clamped at This surface force is biperiodic, with the ratio of the associated low frequency over the fast one equal to 10 −4 .
The constitutive relation links the normal (nonviscous) stress to the longitudinal strain along the bar for any ( where E is the Young's modulus, and ε p (x, t) stands for the longitudinal plastic strain. One assumes that this latter verifies a Norton's viscoplastic evolution law for any (x, t) ∈ [0, L] × [0, T ], which depends on the nonviscous stress only: . K and n are the two constant parameters to be identified. The numerical implementation consists of quadratic 1D finite elements under MAT-LAB for the spatial discretization of the displacement, whereas the longitudinal plastic strain is linearly interpolated using the 'external' nodes of the different elements. The time integration uses the 'ode45' procedure, based on a fourth-order embedded Runge-Kutta formula according to Ref. [21]: depending on the number of finite elements used, the time step to be chosen can verify the classical rule of thumb of 20 time steps per fast cycle, or must be smaller than the maximal time step associated with the Courant-Friedrichs-Lewy condition, making in both cases the computational cost become unaffordable as soon as a high number of cycles is calculated.

Parameter identification process
It is assumed that experimental data consist of the knowledge of the displacement u exp (t) measured at x = L. The proposed misfit function is then: where u(x, t; K , n) is the solution of the forward problem (11)-(15) with parameter values K and n.
In order to derive the adjoint state problem, the following Lagrangian is introduced: where (u, ε p , K , n, z, ) are considered as independent quantities. The minimization of the misfit function (16) under the constraints (11)-(15) is then equivalent to express the stationarity of the Lagrangian (17) with no constraint [19]. The adjoint states z(x, t) and (x, t) are Lagrange multipliers verifying, for any (x, t) ∈ (0, L) × [0, T ], equations coming from the stationarity of the Lagrangian with respect to u: Formally, the latter equations correspond to a dynamic elasto-viscoplastic problem (with negative damping), when one considers the first Lagrange multiplier z as a kind of displacement (even if it is not homogeneous to a length) and the second Lagrange multiplier as a kind of plastic strain (even it is not dimensionless). The discrepancy between the forward model's predictions and experimental data is directly imposed as a Neumann boundary condition at x = L, i.e. where the measurements are available.
One should notice that this problem is time-backward with final conditions (at t = T): once the adjoint state problem has been properly discretized with respect to space, the obtained time differential equations can be solved with classical time integration methods by applying a change in variables θ = T − t in order to deal with initial conditions rather than final ones. Once solved, the adjoint state solution is used to evaluate the misfit function's gradient, whose components are equal to the partial derivatives of the Lagrangian (17) with respect to the parameters: One should notice that, beside the forward state solutions, only the second Lagrange multiplier appears in the previous expressions: this is logical here, because one is interested in the parameters (K , n) of the evolution law (15) only. The first Lagrange multiplier z would have appeared in these expressions if one had tried to identify parameters associated with the dynamic equilibrium equation, such as the Young's modulus E or the damping ratio c K .
These estimates of the misfit function's gradient components are used by the minimization algorithm to determine the best search direction and step size. However, even if the strategy using the adjoint state problem is the cheapest one in terms of computational cost, this latter can be awfully prohibitive if a high number of fast cycles has to (20) ∂J ∂n (K , n) = ∂L ∂n (u(K , n), ε p (K , n), K , n, z(u, ε p ), (u, ε p )) calculated, as it is the case here, since viscoplasticity is a slow-evolving phenomenon, implying that the change in the longitudinal viscoplastic strain is very small for each individual fast cycle.

Time-homogenized problem
In order to obtain a drastic reduction in the computational cost, the periodic time homogenization method is applied. Following the different steps briefly described in the "Background", the time-homogenized equations to be solved are obtained.

Dynamic equilibrium equation
The first equation to be dealt with is the dynamic equilibrium equation, which can be rewritten in a normalized way in order to compare the orders of magnitude of the different terms: where F is the frequency of the slow component of the loading, and the different normalized quantities are denoted as •. Different cases arise, depending on the orders of magnitude of ρL 2 F 2 /E and Fc K . Here, we assume that: where β ≤ O(1) and γ ≤ O(1). Physically, the first assumption (26) is equivalent to having L/ = √ β, where is the wavelength associated with the fast loading frequency: this means that the order of magnitude of should be at least equal to the length of the bar.
Using these two assumptions and introducing an asymptotic expansion of the displacement u = u 0 + ξ u 1 + ξ 2 u 2 + O(ξ 3 ) allows to express the acceleration as follows: It can then be shown [5] that the two time scales are separable and that the homogenization method can be applied: the zeroth-order dynamic equilibrium equation becomes: and the associated zeroth-order boundary condition at x = L reads: Eventually, by taking the fast-time average (4) of the two previous equations, and by using the quasiperiodicity assumption (5), one can obtain the time-homogenized zerothorder dynamic equilibrium equation: for any (x, t) ∈ (0, L) × [0, T ], and the time-homogenized zeroth-order boundary condition: which corresponds to the solution of a quasistatic problem.

Evolution law
In the same way, the evolution law (15) gives, up to order zero: First, the 1/ξ-order has to be considered, giving: ε p 0 ′ = 0. This is equivalent to say that the zeroth-order plastic strain only depends on the slow time variable: The physical interpretation of this latter result is that viscoplasticity is a slow-evolving phenomenon, even if the material withstands a high-frequency loading. This slow variation is described by the zeroth-order evolution law, which can be expressed as: In order to make the contribution of the first-order plastic strain disappear from this relation, the fast-time average (4) of this latter is evaluated. Since, on the one hand, Eq. (34) is equivalent to say that �ε with a zero initial condition.

Zeroth-order time-homogenized problem
The final system consists of the following equations to be solved at the slow time scale only, for any (x, t) ∈ (0, L) × [0, T ]: The resulting system then corresponds to a quasistatic, elastic viscoplastic problem, where the loading consists of the slow component of the surface force. Of course, the fast component does have an effect on these homogenized equations, by means of the evolution law (41): since this latter is nonlinear, it is required to calculate the zeroth-order 'instantaneous' stress field σ 0 instead of directly using the corresponding homogenized quantity σ 0 . To do this, it is mandatory to estimate the residual quantity σ * 0 (x, t, τ ), defined according to the following relation: This residual is then the solution of the following system of residual equations: obtained by subtracting the equations of the time-homogenized problem (37)-(40) to the zeroth-order equations of the reference problem. This residual system consists in solving, for the fast time scale τ, and at each slow time step t k , a dynamic problem for a purely viscoelastic material, with a surface force f * s = f s − �f s �. Since viscoplasticity (which is the only source of nonlinearity here) appears in the zeroth-order time-homogenized equations only, this residual problem is linear, and can be conveniently solved in the frequency domain. If needed, additional orders could be addressed, as detailed in Ref. [5].
The evaluation of the fast-time average appearing in (41) can be achieved by means of a numerical integration formula, such as the trapezoidal rule with N + 1 points, defined as: This choice is not only related to its formal simplicity, but can actually be justified using the conclusions presented in Ref. [22]: the author shows that the trapezoidal rule converges extremely fast when integrating smooth periodic functions, as it is the case here, and that no substantial additional gain would be obtained with the use of more elaborated formulas, such as Simpson's for example.

Associated parameter identification strategy
The identification strategy using adjoint state formulations described in the "Background" is applied now on the time-homogenized problem. First, a consistent misfit function has to be selected, then the adjoint state problem required to estimate the misfit function's gradient is introduced.

Choice of the misfit function
The first step consists in describing, through the misfit function, the discrepancy between the time-homogenized model's predictions and experimental data: indeed, on the one hand, the model has been solved on slow time steps only, whereas, on the other hand, the experimental data can be available on a much finer scale. The most efficient choice in terms of computation cost is to use time-homogenized quantities in the misfit function (16), since it allows to address the time integral on slow time steps only: Whereas �u 0 �(L, t; E, K , n) is the solution of the zeroth-order time-homogenized forward problem (37)-(41), u exp (t) stands for the corresponding experimental quantity, which is obtained by fast-averaging the experimental data for each slow time step t k of the time-homogenized displacement:

Adjoint state problem
The misfit function's gradient is evaluated in the same way as in the section associated with the reference problem, using the solution of an adjoint state problem. Starting from the Lagrangian associated with the misfit function (48) and the time-homogenized forward problem (37)-(41), it is first obtained that the second Lagrange multiplier depends on the slow time scale only: because it is associated with the zeroth-order time-homogenized evolution law (41), which depends on t only. On the contrary, the first Lagrange multiplier is decomposed as �z 0 �(x, t) + z * 0 (x, t, τ ) to deal with the two systems (37)-(40) and (43)-(46) respectively. Then, the two following PDEs corresponding to the adjoint state problem are obtained, for any (x, t) ∈ (0, L) × [0, T ]:  corresponding adjoint state solutions and 0 for the identification process detailed in the "Results and discussion": the adjoint state solution 0 corresponding to the timehomogenized problem is homogeneous along the bar, and is very close to the (homogeneous) fast-time average of the (heterogeneous) adjoint state solution associated with the reference problem.
Indeed, in this case, the previous equations can be decoupled in order to get the PDE verified by the second Lagrange multiplier 0 for any (x, t) ∈ [0, L] × [0, T ]: which is a time-backward ODE with a final condition equal to zero. This equation can be solved using the slow time steps t k only, which allows to derive the solution in a way as efficient as for the time-homogenized forward solution.

Misfit function's gradient
The misfit function's gradient then consists of the two following partial derivatives: Once again, these relations correspond to the zeroth-order time-homogenized estimates of the two misfit function's gradient components (23)-(24) obtained for the identification problem associated with the reference problem. Actually, this is a result that already occurs in periodic space homogenization, as shown for example in Ref. [23].

Validation of the time homogenization method
The studied bar, which is 1 m long, is discretized in ten quadratic finite elements. The surface force at x = L is a combined cycle load, defined as the sum of two sines of frequencies F and F /ξ. For this case of study, the assumption (26) is verified, since the two loading frequencies are F = 0.05 Hz and F /ξ = 500 Hz (cyclic loads with respective The associated results are compared in Figure 6 presenting how evolves, for the first slow period, the zeroth-order longitudinal plastic strain at x = 0 and x = L: the estimates are in excellent agreement, as confirmed by Table 1, and tend to indicate that the predictions can remain good even if a high number of slow cycles is applied. Figure 7 shows the previous time evolutions (at x = 0) between two close time steps, and allows us to see that the zeroth-order longitudinal plastic strain evolves smoothly, while the reference strain increases step by step for each fast loading period. All these results allow to validate with this typical example the periodic time homogenization method, and to get an estimate of the computational gain it allows through the drastic reduction in the

Table 1 Validation of the periodic time homogenization method
Reference and zeroth-order time-homogenized estimates, at t = 20 s, of the longitudinal plastic strain at x = 0 and x = L. number of required time steps for the numerical simulation. This number is indeed the main criterion concerning the computational cost of the time homogenization method, since, when compared with the reference problem, the only additional calculations consist in solving the fast residual system (once and for all) and in estimating the fast-time averages with the trapezoidal rule (47).

Parameter identification
In order to evaluate how the misfit function performs, synthetic data u exp (t) are created by solving the Eqs. (11)-(15) of the reference forward problem, using K exp = 100 × 10 6 SI Units and n exp = 10 as parameter values. As previously, the bar withstands a combined cycle loading with F = 0.05 Hz and F /ξ = 500 Hz. K 0 = 50 × 10 6 SI Units and n 0 = 5 are chosen as initial parameter values for the identification process, which is based on an interior-reflective Newton method [24] with a BFGS formula in order to minimize the misfit function J 0 . Since we do not want to address specifically the question of regularization, no regularization term is added to the misfit function (48). Eventually, given the fact that the two parameters have very different orders of magnitude, a normalization is introduced, using the values of K 0 and n 0 , in order to prevent numerical inaccuracies in the gradient estimate. The logarithm of the misfit function is depicted in Figure 8, showing a well-defined global minimum despite a crescent-shaped valley. Subsequent results are listed in Table 2: the parameter values identified are very close to the ones used to create the synthetic data. Figure 9 shows the comparison between the identified model and the synthetic reference, more precisely the variations of the longitudinal plastic strain, which is not directly observable.
In addition, to demonstrate the robustness of the process, centered Gaussian noise is added to the synthetic data used in the identification strategy: different levels (corresponding to a standard deviation proportional to the mean of the measured displacement) are proposed, and the associated identification results are listed in Table 2. They show that the identification works, even for high levels of noise: the choice of the misfit function (48) with the application of the fast-time average to the noisy data has allowed to filter significantly the introduced noise, hence strongly limiting its impact on the identification results. This is a particularly interesting property, whereas identifying using the reference model and the whole set of noisy data undoubtedly would have led to poor results: indeed, as said above, the fact of using a homogenized model in a parameter identification process generally introduces a regularizing effect, as explained in Ref. [17]; unfortunately, even in this academic case of study, the numerical cost for solving the reference identification problem is too prohibitive to rigorously check this property.
Eventually, the computational cost associated with the identification process is significantly reduced when compared with what is obtained when the inverse problem related

Table 2 Parameter identification process
Results obtained for synthetic data (with K exp = 100 × 10 6 SI Units and n exp = 10) and different noise levels. to the reference problem (11)-(15) is considered, which allows to address the case of slow-evolving phenomena, requiring the simulation over a large time interval.

Conclusions
Here we have proposed a first preliminary study of a two-time-scale parameter identification process, using time-homogenized models: the adaptation of a classical identification strategy based on an adjoint state formulation to estimate the misfit function's gradient can be used in this specific framework: this leads, on the proposed example, to the determination of the time-homogenized counterpart of the adjoint solution associated with the reference identification problem. Despite its simplicity, the academical example studied here showed the relevance of the strategy and its reduced computational cost: these results can be viewed as a first step before dealing with more complex cases of study.
To come back more specifically to CCF life prediction, a possible prospect could be to use the experimental information available (with classical high cycle fatigue tests for instance) to tune correctly the parameters of the considered material laws. In addition, to further improve the reduction in the computational cost, the extension to models with three different time scales, such as those described in Ref. [5], should be straightforward as well.