Large scale parameter estimation problems in frequency-domain elastodynamics using an error in constitutive equation functional

https://doi.org/10.1016/j.cma.2012.08.023Get rights and content

Abstract

This paper presents the formulation and implementation of an Error in Constitutive Equations (ECE) method suitable for large-scale inverse identification of linear elastic material properties in the context of steady-state elastodynamics. In ECE-based methods, the inverse problem is postulated as an optimization problem in which the cost functional measures the discrepancy in the constitutive equations that connect kinematically admissible strains and dynamically admissible stresses. Furthermore, in a more recent modality of this methodology introduced by Feissel and Allix [17], referred to as the Modified ECE (MECE), the measured data is incorporated into the formulation as a quadratic penalty term. We show that a simple and efficient continuation scheme for the penalty term, suggested by the theory of quadratic penalty methods, can significantly accelerate the convergence of the MECE algorithm. Furthermore, a (block) successive over-relaxation (SOR) technique is introduced, enabling the use of existing parallel finite element codes with minimal modification to solve the coupled system of equations that arises from the optimality conditions in MECE methods. Our numerical results demonstrate that the proposed methodology can successfully reconstruct the spatial distribution of elastic material parameters from partial and noisy measurements in as few as ten iterations in a 2D example and fifty in a 3D example. We show (through numerical experiments) that the proposed continuation scheme can improve the rate of convergence of MECE methods by at least an order of magnitude versus the alternative of using a fixed penalty parameter. Furthermore, the proposed block SOR strategy coupled with existing parallel solvers produces a computationally efficient MECE method that can be used for large scale materials identification problems, as demonstrated on a 3D example involving about 400,000 unknown moduli. Finally, our numerical results suggest that the proposed MECE approach can be significantly faster than the conventional approach of L2 minimization using quasi-Newton methods.

Highlights

► We develop an ECE strategy suitable for large scale parameter identification. ► We propose a continuation scheme for accelerating convergence. ► We formulate a block-SOR scheme for large scale problems. ► The proposed MECE scheme showed faster convergence than a conventional L2/BFGS algorithm.

Introduction

The identification of material parameters (e.g. elastic moduli) are of paramount importance in science and engineering. For instance, this type of inverse problem arises in connection with applications such as damage detection and structural health monitoring, seismic exploration, and biomechanical imaging, among other applications. Models of complex engineering structures often feature local or heterogeneous parameters that are unknown and, therefore, need to be identified by exploiting experimental information on the mechanical response of the structure. Considerable research effort has been dedicated to develop algorithms for material identification. A review on these methods in the context of elasticity can be found in [1]. Despite current advances, numerical algorithms are still often challenged by the inherent ill-posedness of inverse problems, especially when, as in this article, the identification of heterogeneous material properties is considered.

Parameter identification problems are often solved using nonlinear optimization schemes. The most common form of such approaches involves minimizing the L2-norm of the error between computed and measured responses (e.g. displacement, strains) [2], [3], [4], [5], [6], [7]. Quasi-Newton methods, which require only the computation of the gradient of the objective function using adjoint methods, are often preferred due to their ease of implementation [3], [8]. Full-Newton methods have also been successfully used for large-scale identification problems [9]. The latter methods converge significantly faster than quasi-Newton methods [8], but are more difficult to implement. In general, gradient-based nonlinear optimization methods have the advantage of allowing the handling of a large number of unknowns and noisy data (using regularization techniques). The main disadvantages of these approaches, when a L2-norm misfit functional is used, is their sensitivity to the initial guess due to the non-convexity of the functional and the large number of iterations required to get an acceptable solution when quasi-Newton methods are used.

In recent years, a new paradigm for inverse material identification has emerged which uses the concept of error in constitutive equations (ECE) for defining cost functionals whose physical meaning is stronger than that of usual L2 functionals and directly related to the material identification problems at hand. The basic premise in the ECE approach is that, given an over-determined set of boundary or internal data (e.g. displacements and tractions), a cost functional is defined based on the least residual (measured in terms of an energy norm) in the constitutive equations that connect stresses and strains that are constrained to be dynamically and kinematically admissible, respectively, with respect to the available experimental data. This type of cost functional has the important property of being zero for the exact constitutive equations and strictly positive otherwise. ECE functionals have been initially introduced for error estimation in FEM computations [10] before being also applied to various material parameter identification problems under linear static [11], nonlinear quasistatic [12], time-harmonic [13], [14], [15] or, more recently, transient [16], [17], [18] conditions. ECE functionals were also found to be useful for solving data completion (Cauchy) problems [19]. Similar energy-error functionals have been introduced for the identification of scalar spatially-varying conductivity coefficients in e.g. [20], [21], with mathematical and numerical issues also discussed in [22], [23].

In particular, a new ECE-based method, named Modified Error in Constitutive Equations (MECE) approach, was recently proposed for time-domain dynamics problems by Feissel and Allix [17] and Nguyen et al. [18]. In MECE, fields and equations are separated into reliable and unreliable sets. The reliable set typically includes the equilibrium equations, initial conditions and known boundary conditions, while the unreliable set includes measured data, constitutive properties, and (when applicable) imperfectly known boundary conditions. The identification problem is then posed as the minimization of a MECE functional, with reliable equations enforced strictly (e.g. as constraints, using Lagrange multipliers) and unreliable equations treated through an ECE-type functional and penalty terms. Allix and Feissel [17] showed that MECE is very robust and can produce reliable identification results in problems with high levels of data noise.

To our best knowledge, MECE-based functionals have not yet been applied to large-scale dynamical identification problems involving unknown heterogeneous material parameters despite their clear potential advantages for this type of problem. With respect to such large-scale applications, the main drawbacks of existing MECE-based identification formulations are twofold. First, they involve a coupled system of partial differential equations (PDEs) that needs to be solved for each iteration. Second, a penalty parameter, which plays a very important role in accuracy of the inverse problem solution, must be provided a priori.

This paper is devoted to the formulation, implementation and validation of the MECE for large-scale heterogeneous material parameter identification problems in the context of frequency domain elastodynamics. Our main contributions in this work are threefold. First, we provide a simple continuation approach for evolving the penalty parameter that appears in MECE, eliminating the need to estimate this parameter a priori; second, we show a Successive Over Relaxation (SOR) scheme that can be suitable for solving large scale problems with MECE; and third, we demonstrate the feasibility of using the MECE method in two-dimensional and three-dimensional time-harmonic elasticity imaging problems with numerical experiments involving up to 400,000 unknown material parameters.

The article is organized as follows. The inverse problem and the MECE-based identification approach are introduced in Section 2. Proposed refinements of the MECE approach are presented in Sections 3 Normalization and a continuation scheme for the penalty term, 4 Successive over-relaxation (SOR) technique for the solution of the coupled problem, leading to the algorithm of Section 5. Section 6 is devoted to 2D numerical examples designed to assess the performance of the main ingredients of the proposed MECE algorithm (evolutive penalty parameter, block-SOR technique for solving the MECE optimality equations), and to compare them against a L2-BFGS minimization approach. Finally, results on a large-scale 3D identification problem are shown in Section 7, before closing the paper with concluding remarks.

Section snippets

Forward problem

In this section, we briefly describe the strong and weak forms of the forward steady-state elastodynamics problem. Let Ω¯=ΩΩRd (1d3) denote a bounded and connected body. The governing equations for a linear elastic body undergoing time-harmonic motion can be written as·σ+b=-ρω2uinΩu=u0onΓuσ·ns=tonΓtσ=C:[u]=12(u+uT)where u is the displacement field, ω is the angular frequency, ρ is the mass density, b is the body force density, denotes the linearized strain tensor, σ denotes the

Normalization and a continuation scheme for the penalty term

In this section, we address a very important issue in the MECE approach: choosing or determining the penalty coefficient κ that appears in the MECE functional (14). As shown in [17], the penalty coefficient plays a fundamental role in the quality of the reconstruction of material parameters, especially in cases where high levels of noise are present in the data. Our goal is to devise a simple strategy for estimating this coefficient adaptively. To this end, we first propose the following form

Successive over-relaxation (SOR) technique for the solution of the coupled problem

One of the main computational challenges that arises in the MECE method is obtaining a solution for the coupled system (32). This can be achieved by either employing direct or iterative solvers on the whole system or using a block strategy. Applying direct linear solvers to (32) for large scale problems may be prohibitive. On the other hand, designing efficient and effective pre-conditioners for the entire coupled system is not a trivial task if iterative solvers were to be used. Hence, our

MECE algorithm

In summary, each iteration in the MECE approach consists in solving the coupled system given in Eq. (32) followed by updating the material parameters using the formulae given in Eq. (31a), (31b). These steps are executed until a maximum number of iterations is reached or a predefined error tolerance is satisfied. The flow of an implementation of the MECE approach is as follows.

  • Set initial values for shear and bulk moduli, fix lower and/or upper bounds and normalization factors.

  • Loop until a

2D numerical examples and results

This section presents a set of 2D numerical examples that were designed to study the performance of the MECE approach described herein for reconstructing spatially varying bulk and shear moduli in the presence of incomplete and noisy data. A 3D example using a massively parallel iterative solver for (32) will then be presented in Section 7. The examples presented herein are motivated from inverse elasticity problems related to biomedical applications [4], [31], [32]. In these problems,

3D example using parallel algorithm

In this section, we consider a three-dimensional example in which the proposed SOR algorithm and the FETI-H parallel linear solver are used to obtain solutions of the coupled block system (32) at each MECE iteration. The domain of the problem consists of a stiff cylinder embedded in a soft cubic matrix as shown in Fig. 11. The “true” Young’s modulus and Poisson’s ratio of the matrix were taken as Ematrix=1Pa and νmatrix=0.3, respectively. The corresponding values for the inclusion were Einc=10Pa

Concluding remarks

We proposed a methodology based on the Modified Error in Constitutive Equations (MECE) approach that is suitable for large scale inverse identification of elastic parameters in frequency-domain dynamics. To this end, we proposed a continuation scheme for the penalty term that appears in MECE. This continuation scheme leads to faster convergence and improved accuracy in the reconstructed solution as compared to the case when a constant penalty coefficient is used. Furthermore, we put forward a

Acknowledgments

The work of Biswanath Banerjee and Wilkins Aquino was partially supported by the United States National Institute for Biomedical Imaging and Bioengineering, Grant# 5R01EB002640-13. The work of Wilkins Aquino was also supported by the CSAR Program at Sandia National Laboratories. Sandia is a multiprogram laboratory operated by Sandia Corporation, a Lockheed Martin Company, for the United States Department of Energy (DE-AC04-94AL85000).

References (34)

  • J. Nocedal et al.

    Numerical Optimization

    (2000)
  • I. Epanomeritakis et al.

    A Newton-CG method for large-scale three-dimensional elastic full-waveform seismic inversion

    Inverse Probl.

    (2008)
  • P. Ladevèze et al.

    Error estimate procedure in the finite element method and applications

    SIAM J. Numer. Anal.

    (1983)
  • G. Geymonat et al.

    Identification of mechanical properties by displacement field measurement: a variational approach

    Meccanica

    (2003)
  • F. Latourte et al.

    Elastoplastic behavior identification for heterogeneous loadings and materials

    Exp. Mech.

    (2008)
  • P. Ladevèze et al.

    Updating of finite element models using vibration tests

    AIAA J.

    (1994)
  • P. Ladevèze et al.

    Application of a posteriori error estimation for structural model updating

    Inverse Probl.

    (1999)
  • Cited by (71)

    • A domain decomposition strategy for mCRE-based model updating in dynamics

      2023, Computer Methods in Applied Mechanics and Engineering
    • Numerical study based on the Constitutive Relation Error for identifying semi-rigid joint parameters between planar structural elements

      2021, Engineering Structures
      Citation Excerpt :

      On its augmented version, the Modified Constitutive Relation Error (MCRE) is characterized by the adding a regularization term to the performance function. This new form broadens the possibilities of including the available information instead of making new assumptions [26,27]. The CRE technique can also be adapted to modern experimental devices.

    View all citing articles on Scopus
    View full text