Three-dimensional transient elastodynamic inversion using the modified error in constitutive relation

This work is concerned with large-scale three-dimensional inversion under transient elastodynamic conditions by means of the modified error in constitutive relation (MECR), an energy-based cost functional. A peculiarity of time-domain MECR formulations is that each evaluation involves the computation of two elastodynamic states (one forward, one backward) which moreover are coupled. This coupling creates a major computational bottleneck, making MECR-based inversion difficult for spatially 2D or 3D configurations. To overcome this obstacle, we propose an approach whose main ingredients are (a) setting the entire computational procedure in a consistent time-discrete framework that incorporates the chosen time-stepping algorithm, and (b) using an iterative SOR-like method for the resulting stationarity equations. The resulting MECR-based inversion algorithm is demonstrated on a 3D transient elastodynamic example involving over 500,000 unknown elastic moduli.

This work aims at formulating and demonstrating a MECR-based approach suitable for large-scale three-dimensional inversion under transient elastodynamic conditions. The forwardbackward coupling bottleneck is overcome by (a) setting up and exploiting the stationarity conditions in a time-discrete framework that incorporates the chosen time-stepping algorithm, and (b) an iterative SOR-like treatment for the resulting stationarity equations.
where u is the displacement field, ρ is the mass density, σ denotes the stress tensor, f is a given body force density, g and Γ N ⊂ Γ are the given surface force density (traction) and its support, and n is the unit normal vector pointing outward from Ω; (b) the kinematic compatibility equations where ε[u] denotes the linearized strain tensor associated with the displacement u and Γ D = Γ\Γ N is the constrained part of Γ; (c) the homogeneous initial conditions and (d) the constitutive (linear elastic) relation, involving the fourth-order elasticity tensor C: Inverse problem. Supplementing the basic equations (1)-(4), additional information in the form of measured transient displacements u in Ω m ×[0, T ] (where Ω m ⊂ Ω is a measurement region) is assumed to be available. The inverse problem then consists in finding the spatial distribution of the elastic moduli C from the data u, subject to equations (1)-(4) with given excitations f , g.
Modified error in constitutive relation approach. Measurement noise is to be expected, in which case the exact enforcement of experimental data is normally not desirable. Accordingly, a modified error in constitutive relation (MECR) functional [2,6] is defined by treating the discrepancy between measured quantities and their computed counterparts as a penalty term: where (f , g) stands for the L 2 (Ω) scalar product of f and g, ξ > 0 is a tunable dimensionless weight and f a fixed constant ensuring consistent units. In this framework, the inverse problem is cast as the optimization problem whose solution (u , σ , C ) satisfies equations (1) to (3) (i.e. balance, compatibility and initial conditions) while achieving a compromise between (i) satisfying the constitutive relation (4) and (ii) matching the measurements u. When using MECR functionals of the form (5), the mismatch between an assumed material property C and measurements u is quantified through the modified error in constitutive relation E(C) := U (u C , σ C , C), where the mechanical fields (u C , σ C ) are determined by the partial minimization of U (u, σ, C) for given C, and in particular depend on C: A natural approach for solving the minimization (6), previously used for frequency-domain problems in [3], then consists in an alternating-directions strategy whereby the transition from a with (u C , σ C ) solving the partial minimization (7), and (ii) finding the material update C by solving the partial minimization problem where (u C , σ C ) are the outcome of (7) and are kept fixed. This approach has been implemented in this work and is employed for the numerical experiments.
Discrete formulation of MECR-based inversion. We now consider a time-discrete formulation of the transient dynamics involved in the definition and exploitation of the MECR functional, based on a Newmark time-stepping scheme. Introducing a sequence t 0 = 0, t 1 = h, . . . , t k = kh, . . . , t n = N h = T of discrete time instants (with a constant time step h = T /N used for simplicity), we let f k denote the value at t = t k of a generic time-dependent quantity f . The Newmark update relations linking the displacements u k , the velocities v k :=u k and the accelerations a k :=ü k , which introduce two additional algorithmic parameters β, γ and read (whereβ := 1 2 −β andγ := 1−γ, and with β = 1/4, γ = 1/2 used in the numerical experiments). In this setting, introducing for convenience the shorthand notations u, v, a, σ for denoting mechanical field histories (u k , v k , a k , σ k ) 0≤k≤N , the inverse problem consists in reconstructing the spatial distribution C of the elasticity tensor, and additionally in finding the mechanical field histories σ and (u , v , a ), such that where the admissibility spaces U and S(a) respectively incorporate kinematic constraints (2), (3) and balance equations (1), and the time-discrete counterpart U N of the MECR functional (5) is defined by Lagrangian and optimality conditions. As problem (10) involves kinematical and dynamical admissibility constraints embedded in spaces U and S, it is now reformulated in terms of its first-order Karush-Kuhn-Tucker necessary optimality conditions. The latter are derived from the following Lagrangian functional: which incorporates the MECR functional (11) together with the constraints resulting from (i) the initial balance equation, with multiplier fieldū 0 ; (ii) the current balance equations, with multiplier fieldsū k+1 and (iii) the Newmark update equations (9), with multiplier fieldsā k+1 andv k+1 . Moreover, the quantitiesū,v,ā are time-discrete sequences of Lagrange multiplier fields, (·, ·) denotes L 2 (Ω) scalar products, and F k are linear forms arising from the prescribed time-dependent excitations f , g featured in (1). The (first-order) necessary optimality conditions are found by setting to zero the first-order derivatives of the Lagrangian (12) with respect to (a) the stress history σ, (b) the kinematic histories (u, v, a, σ), (c) the multiplier histories (ū,v,ā) and (d) the elasticity tensor C. (a) Differentiating L with respect to σ, the stress history in Ω is found to be given in terms of the kinematic histories by while the multiplier history is found to satisfy the boundary condition (b) Setting to zero the derivatives of L with respect to the kinematic histories (u, v, a), using (13) to eliminate stresses, and performing a finite element discretization, the following forward time-stepping equations are obtained for the space-discrete histories (u, v, a) and (ū,v,ā): (i) Initial conditions: (ii) Forward time-stepping equations (0 ≤ k ≤ N − 1): Ku k+1 + Ma k+1 = Kū k+1 + F k+1 , and Newmark update equations (9). (15b) (c) Setting to zero the derivatives of L with respect to the multiplier histories (ū,v,ā) simply restores all constraints introduced in L. Upon elimination of stresses using (13) and finite element discretization, the following set of backward time-stepping equations is obtained.
(i) Final conditions: (iii) Backward transition equation (last): Steps (a), (b) and (c) achieve the partial minimization (7), where C is kept fixed. K, M, F k respectively denote the stiffness and mass matrices arising from FE discretization, and the nodal excitation at time t k , while the matrix D corresponds to L 2 scalar products on Ω m . (d) This step, constituting the second part of the alternating-direction approach, implements the partial minimization (8). Enforcing the first-order necessary optimality condition ∂ C U (u C , σ C , C)[Ĉ] = 0 ∀Ĉ for the latter problem and assuming isotropic linear elasticity, the explicit pointwise updating formulas are found, with B and G denoting (spatially-dependent) bulk and shear moduli, respectively. Since B , G > 0 by construction, any such update C of C is positive definite.  where the definition of matrices A, C, D and vectors F, G can be easily obtained by identification from (15a,b) and (16a-d). In particular, the matrix A is associated to the forward Newmark scheme, and the backward time-stepping can be shown to be represented by the adjoint matrix A T . This implies in particular that the forward and backward time stepping schemes have identical stability conditions. Then, the iterative scheme is defined by introducing auxiliary unknowns W,W through so that one iteration of the SOR algorithm consists, for given (X (i) ,X (i) ) in solving equations (where 0 < η < 2 is the relaxation parameter) and obtain (X (i+1) ,X (i+1) ) from (19) with the converged values of W,W . The SOR iterations (20) are quite simple to implement using existing FE dynamical analysis codes, as they involve the usual stiffness and mass structural matrices. We used the finite element code SIERRA/SDA [13], which has massively parallel capabilities including efficient multilevel domain decomposition strategies for solving linear systems.
3D numerical example. A unit cube containing a cylindrical inclusion with an ellipsoidal cross section (Figure 1), its bottom surface being fixed, is loaded on its other faces with a normal uniform time-harmonic force density, whose magnitude is one unit on the top surface and 0.5 units on the side surfaces (load frequency 1 Hz, total load duration 1 s, with N = 100 time steps of h = 0.01 s). The true shear and bulk moduli of the background material are taken as G = 2 and B = 4, while those of the inclusion are G = 6 and B = 16. The 3D domain was discretized using 414, 000 tetrahedral finite elements (75, 000 nodes). The data u is synthetically generated by (i) computing all displacement components at all nodes of the fine mesh for 100 time steps, (ii) interpolating them onto a coarser uniform mesh (i.e. cylindrical inclusion not included), consisting of 275, 000 tetrahedral elements (50, 000 nodes), and (iii) polluting it with Gaussian noise with relative noise level δ = 0.01. The total number of unknown moduli for this example is 550, 000. We note that reconstruction of moduli from interior data is the main focus of the active research area of biomechanical imaging (see e.g. [12] and references therein) .The value of ξ is selected according to δ via Morozov's discrepancy principle, i.e. choosing the largest value of ξ meeting the condition (wherein (u k ) 0≤k≤N satisfy the optimality system (18)): Figure 1 shows threshold plots of the reconstructed shear and bulk moduli, respectively, with the lower threshold set to 5 and 12, respectively, for plotting G and B. Either plot shows the inclusion location and geometry to be correctly identified (with some aberrations occurring near the top surface, which may result from measurements being less sensitive to moduli in those areas). Figure 2 shows cuts perpendicular to the x, y and z axes, going through the cube center.