Towards a non intrusive method for solving coupled direct retrograde problem in transient dynamics

– The objective of this work is to solve within a standard ﬁnite element software, the direct-retrograde coupled time hyperbolic problems that arise when dealing with the identiﬁcation of material parameters in dynamic in presence of corrupted measurements by means of the adjoint state approach. This question has given rise among others to Riccati methods and shooting methods. Those methods have several drawbacks, namely their numerical complexity and the fact that their implementation in a commercial software would require to entirely re-design the implementation, which is nearly impossible. Therefore we are seeking an iterative method which could be easily implemented. The proposed method solves alternatively a direct problem and a retrograde one. In order to ensure its convergence a relaxation scheme is applied and optimized. Comparisons between the proposed approach and Riccati and shooting methods in terms of complexity, numerical cost and robustness are given in the case of an elastic bar.


Introduction
The identification of physical structural parameters is a difficult task especially in the context of transient dynamics. Therefore this domain has given rise to a lot of proposals and dedicated methods [1][2][3][4][5][6]. A difficulty is that those methods often need to design new specialized softwares. This is also the case of the method we proposed in order to deal with the identification of the dynamic behavior of laminated composites up to failure [3,[7][8][9]. We are thus seeking a method which could be used in a standard finite element code with as few implementation as possible. This paper is a first step in that direction.
The type of inverse problem we are focused on concerns the identification of material parameters in dynamics in the case of corrupted measurements. The identification strategy used is based on the concept of the modified error in the constitutive relation [10].
The minimization of the functional under constraints associated to the formulation leads to coupled directretrograde wave propagation problem which exhibit exponential solutions, thus making the problem highly illconditioned. This problem has previously been solved by a Corresponding author: nouisri@lmt.ens-cachan.fr combining two technics: the resolution of an algebraic equation of Riccati [2,11], which in the non-linear case was adressed using the LATIN method to solve iteratively the problem globally on the whole time interval [7][8][9]. The implantation of those two methods within a commercial software is beyond what can be done in a research team.
The main step of the identification task is about finding a costless method to solve the coupled directretrograde wave propagation problem. Such method can not be one of the methods based on the construction of a global system over space and time because it leads to the resolution of a largest linear system. It is therefore preferable to seek on methods local in time. Among methods local in time, we find the method based on the transition matrix which is well only for short time domain tests. The method based on the algebraic Riccati equation appears relevant for long time domain tests [12]. A combined version of these two methods [7] may be helpfull to evaluante the solution for any arbitrary duration. However these methods became memory-intensive mainly when addressing dynamic problems with large numbers of degrees of freedom for which time and space discretization are coupled. In order to reduce the algorithmic complexity arising from solving the basic problem associated to the identifiation problem, the basic idea of the proposed method is to solve alternatively, starting from an initial guess, a sequence of direct and retrograde differential equations using a time integration scheme (e.g. a Newmark scheme). However, due to the presence of direct and retrograde exponentially increasing functions of time, such method diverges. A waveform relaxation is needed to improve the convergence by mitigating the resonance phenomenon, which is due to the similarity of the eigenvalues of the two coupled equations, especially for wide time range. A numerical example shows the feasibility of the method and comparisons with Riccati and shooting methods are given.

The reference model
In this paper we have restricted the application of the method to a homogeneous isotropic linear elastodynamic structure whose characteristics (Young's modulus E ref , Poisson's ratio ν ref and density ρ ref ) are known. However the method can be applied to more complicated cases. Nevertheless this simplified development of the strategy based on the concept of the modified error in the constitutive relation in presence of corrupted measurements in transient dynamics [10] here allows to test the capabilities of the method. Moreover, work to extend the method to the nonlinear case is in progress. The reference model consists in an elastic bar of length L and cross-section area S (see Fig. 1) that satisfies the following boundary conditions: on the border x = 0, the structure is embedded; on the border x = L, a half-sine chap load The reference problem.

Experimental data
Direct calculation is obtained by means of an implicit θ-method; experimental displacement u and experimental load σ are then obtained by solving the following problem: equilibrium Equation (1): constitutive relation (2): boundary conditions (4): It's important to note that getting experimental measurements of the force in the clamped section is considered a difficult task. The access to such location needs the installation of a sensor between the structure and the frame and that's going to change the kind of the mechanical link. In this case, the measured data will no longer representative of a real embedment. Otherwise an automated photoelastic equipment can be used to determine the stress distribution in a material until embeded area.

Formulation of the inverse problem
Our aim is to find the suitable space-time fields u and σ depending on the Young's modulus E that at the same time satisfies the constitutive relation and minimizes the distance to experimental datasũ d andf d . The question is about finding a compromise between these equations considered as uncertain under the constraints of reliable ones (e.g. equilibrium and time initial conditions). The inverse problem associated to the identification strategy based on the modified error in the constitutive relation [10,13] can then be written as follows: Find the fields u, σ, u d , f d minimizing (5): Under the constraints (6):

Adjoint problem and derivation
The inverse problem involves the minimization of the cost function (5) under constraints (6). The searched fields would be the ones for which the saddle point of the Lagrangian (7) is reached: u * is the Lagrange's mutliplier associated to the model error; λ is the Lagrange's mutliplier for relaxation of the displacements boundary conditions; μ is the Lagrange's mutliplier for relaxation of the loads boundary conditions. This leads to the system (8) which is given in the strong form over ]0, L[ by: The solution must satisfy the time initial and final conditions (9): and the space boundary conditions (10): The system (8) is then projected into a classical finite element space. This leads to the differential equations (11) in which the searched direct and retrograde fields are respectively the unknown nodal vectors U and U * : U and U * must also satisfy the following initial and final conditions in time (12): -M and K are respectively the mass and the stiffness matrices; Classical incremental methods can not be applied due to the presence of coupling between the forward and backward equations in the system (11), (12). Various methods have been proposed in the literature [1][2][3][4]9] to address this type of system.

Towards a non intrusive method
One can transform the systems (11) and (12) into the first-order state-space system (13): with : and being the state vectors; The solution of (13) must satisfy the time initial and final conditions (14): The application of the method of successive approximations [14][15][16] to the system (13) leads to the following process (15): with the initial and final conditions (16):  This iterative procedure consists in solving alternatively the subsystem with initial time values forward from t = 0 to t = T , then the other subsystem with final time values backward from t = T to t = 0. Starting with a given initial guess field Z (0) , the calculation procedure is continued until convergence. During the resolution of a subsystem, the coupled fields are considered as data, obtained from the previous calculation over the range [0, T ]. However, the redundancy of the boundary conditions and the instability due to the presence of an exponentiel phenomena in the solutions of the homogenous equation associated to the reference problem [15] are the main reasons why no convergence arises from such an iterative scheme.

Relaxation method
To help the convergence, the idea of the previous scheme is maintained, but a relaxation scheme is applied. The evaluation of the forward field Y (t) t∈[0, T ] at the (k + 1) th iteration is obtained from a weighted sum of the fields Z(t) t∈[0, T ] calculated at the (k + 1) th iteration and that calculated at the k th iteration (Fig. 2).
Following this line, Equations (15), (16) become a recurrence process : 1. Initialization is performed by an initial guess Z (0) of the retrograde field, followed by an initial evaluation Y (0) of the direct field by solving (17): 2. At step k, the direct and retrograde fields Y (k) and Z (k) are known. Then the fields Z (k+1) and Y (k+1) are determined respectively by means of (18) the retrograde equation and its final condition in time and (19) the direct equation and its initial condition in time:

Numerical tests
In our experiments, the basic solver used is the implicit θ-method with θ = 1 2 (trapezoidal rule) and the timestep is 0.36 μs for a given finite time interval [0, 0.63 ms]. The direct calculation is used to create measurements on a bar (Fig. 1)  To deal with the relaxation scheme, the initial guess field is set to zero and the error tolerance is 10 −4 . The error used here is the holy norm of the differential Equations (11). Plots in Figure 3 gives an idea of the effect of noise on the convergence of the method. Although there is a wide choice of relaxation coefficients for which a few relaxation-iterations are spent until convergence, the nearest values to 1 are considered the best ones; otherwise we do not have a real relaxation.
We report on the right side (see Fig. 4) errors versus the number of the relaxation iterations arising from  the costless values of the relaxation coefficients obtained from the left side one. Besides, for some number of iterations, we note that the error is small with more noise. This is can be explained by the fact that the convergence is not yet established. Solutions at these iterations are resulting from worst quantifications of incompatibilities of the model's parameters towards the redundant boundary conditions.
Algorithm complexities of the proposed relaxation method, the Riccati method and the shooting method are evaluated by calculating the number of Multiply-And-Accumulate (MAC) operations needed to solve the above coupled problem (see Tab. 1). According to this numeri-cal test, the relaxation method looks more suitable than other ones: indeed its complexity is O(42 k mn 2 ) while the Riccati method is O(1 mn 4 ) and the shooting method is O(320 mn 3 ), with k, m and n denoting respectively the number of relaxation iterations, the number of time steps and the number of d.o.f.
An evaluation of the algorithm complexity for various number of d.o.f. is given in Figure 5: it confirms not only the efficiency of the relaxation method but also its reliability compared to the other methods. Robustness of the method to noise is also guaranteed: according to Figure

Conclusions
For a system of time coupled direct-retrograde second order differential equations here derived in the case of an elastic bar, we proposed an iterative method that requires a reduced implementation time, and for which convergence is guaranteed by a relaxation approach. Numerical tests presented here tend to show the better efficiency and robustness of the proposed method compared to Riccati or shooting methods. The main drawback of this method is about finding the suitable relaxation coefficient for which less iterations are spent to identify the parameters. Next step would be about finding a mathematical criterion to replace the numerical way used to obtain the range of convergence.
Although only the one dimensional elastic problem has been studied here, the same approach can directly be used for other linear two-point boundary value problems. The case of nonlinear materials is now under consideration in order to lookfor the methodology in real-word engineering systems incorporating damping and/or in the presence of nonlinearities.