Inverse problems in elasticity

TL;DR: In this paper, a review is devoted to some inverse problems arising in the context of linear elasticity, namely the identification of distributions of elastic moduli, model parameters or buried objects such as cracks.

Abstract: This review is devoted to some inverse problems arising in the context of linear elasticity, namely the identification of distributions of elastic moduli, model parameters or buried objects such as cracks. These inverse problems are considered mainly for three-dimensional elastic media under equilibrium or dynamical conditions, and also for thin elastic plates. The main goal is to overview some recent results, in an effort to bridge the gap between studies of a mathematical nature and problems defined from engineering practice. Accordingly, emphasis is given to formulations and solution techniques which are well suited to general-purpose numerical methods for solving elasticity problems on complex configurations, in particular the finite element method and the boundary element method. An underlying thread of the discussion is the fact that useful tools for the formulation, analysis and solution of inverse problems arising in linear elasticity, namely the reciprocity gap and the error in constitutive equation, stem from variational and virtual work principles, i.e., fundamental principles governing the mechanics of deformable solid continua. In addition, the virtual work principle is shown to be instrumental for establishing computationally efficient formulae for parameter or geometrical sensitivity, based on the adjoint solution method. Sensitivity formulae are presented for various situations, especially in connection with contact mechanics, cavity and crack shape perturbations, thus enriching the already extensive known repertoire of such results. Finally, the concept of topological derivative and its implementation for the identification of cavities or inclusions are expounded.

Summary

1. Introduction

• Elasticity theory describes the reversible deformation of solid bodies subjected to excitations of various physical natures: mechanical, thermal, electromagnetical etc.
• Such excitations, applied as distributions over the body (e.g. gravitation, Lorentz forces, thermal expansion) or over the boundary (pressure, contact forces), generate strains (i.e. local deformations) and stresses (i.e. local forces) in the material.
• Models of complex engineering structures often feature local parameters that are not known with sufficient accuracy, and therefore need to be corrected by exploiting experimental information on the mechanical response of the structure.
• Constitutive parameters of some common isotropic elastic materials: E Young modulus, ν Poisson ratio, ρ mass density; σ Y and ε Y define the elastic limit, i.e. are the stress and strain levels beyond which the material is no longer elastic.

2.1. Fundamental field equations for three-dimensional elasticity

• The linearized elasticity theory [74] is established on the assumption of small strains, namely |∇u(x, t)| ≪ 1.
• The material is characterized by two constitutive parameters: its mass density distribution ρ(x), associated with the kinetic energy T (u) = 1 2 Ω ρ| u| 2 dV (where the dot denotes time differentiation) and the fourth-order tensor of elastic moduli C(x), hereafter referred to as the elasticity tensor, associated with the elastic strain energy EQUATION.
• The elasticity tensor C defines a positive definite quadratic form over the 6-dimensional space of symmetric second-order tensors.
• Therefore, C has the following symmetries: EQUATION and hence has at most 21 independent coefficients.
• C can be expressed in terms of the Lamé coefficients (λ, µ): EQUATION.

2.2. Direct problems

• The typical direct problem of elastodynamics is the initial-boundary value problem for the unknown displacement field u(x, t) defined by the field equation ( 9), the initial conditions (11) and the boundary conditions ( 12) and ( 13).
• The direct elastic equilibrium problem has a similar structure, except that all field variables are time-independent, so that the field equation ( 9) reduces to the analog for linear elasticity of the Laplace equation, i.e. EQUATION and initial conditions such as (11) are not needed.
• The direct equilibrium problem for thin plates is defined in Section 5.3.1.

2.3. The principle of virtual work

• The principle of virtual work [129] is a fundamental principle of the dynamics of continua, from which the balance equation (7) can then be deduced.
• Equation (18) is general in that it does not refer to specific constitutive properties (e.g. elasticity).
• The following identity holds for elastic equilibrium problems in the absence of body forces: EQUATION.

3. The virtual work principle as an observation equation

• The distribution of displacements and in-plane strains can be measured over the surface of a body by optical means.
• Therefore, displacement and in-plane strain fields on the surface are experimentally available.
• Inversion strategies based on the assumed availability of continuous data have therefore a practical relevance.
• In such situations, the virtual work principle allows to formulate observation equations, i.e. mathematical relationships between the observations and the unknown quantities, in a direct and effective way.
• The virtual fields method will also be briefly described for completeness (section 3.4).

3.1. The reciprocity gap for elastic moduli identification

• The reciprocity gap concept is presented in connection with the identification of an unknown distribution of elastic moduli from complete non-intrusive experiments conducted in quasi-static conditions.
• For general anisotropic elastic media, the following argument suggests that the answer is negative.
• Definition (30) implies that L(y) possesses the major symmetry L ijkℓ =.
• Therefore, one can find fourth-order tensors C having the symmetry properties (3) such that the fourth-order tensor L linked to C by ( 30) is an elasticity tensor.
• A technical error was subsequently found by these authors [117] and Eskin and Ralston [60] .

3.1.2. The reciprocity gap functional.

• Taking into account the symmetry properties (3) of C and C ⋆ and the overdetermined boundary data (27) , one arrives at the identity EQUATION which defines the reciprocity gap R(C ⋆ −C; u ⋆ , w).
• The terminology comes from the fact that the identity (32) can also be established from the Betti reciprocity theorem [129] applied to the two states w and u ⋆ .
• The Betti theorem being based on the assumption that both states refer to the same elasticity tensors, a "reciprocity gap" occurs, here in the form of the domain integral in (32) when that assumption is not correct.
• For any selection of the elastostatic virtual field w (sometimes also termed "adjoint field"), the left-hand side of ( 32) is known, while the right-hand side depends on the unknown moduli C ⋆ .
• It is therefore natural to consider its linearized counterpart, the linearized reciprocity gap.

3.2. Identification of elastic moduli perturbations using the linearized reciprocity gap.

• The linearized reciprocity gap (LRG) identity ( 33) is useful for either theoretical investigations of the linearized inverse problem or numerical inversion of data.
• The reference elasticity tensor C is assumed isotropic, i.e. of the form (4).

3.3. Identification of cracks using the reciprocity gap

• Crack identification problems consist in identifying a crack (or a set of cracks) from a set of overdetermined force-displacement boundary measurements.
• Such problems can be formulated within different physical contexts such as electrostatics, elastomagnetism, acoustics or elastodynamics.
• In particular, the elastostatic RG identity is obtained from ( 49) by removing the integrations in time and suppressing the last term involving time derivatives.
• This remark permits the definition of several numerical techniques to identify the position of a crack.
• The colours of the incoming wavefronts denote the value of the instantaneous RG and the plot is stopped when the instantaneous RG becomes nonzero.

3.4. Direct applications of the virtual work principle

• In some situations, measurements of field variables may be available not only on the boundary but also over the whole sample body.
• This is in particular the case for measurements made on thin samples (for which plane stress conditions may be assumed).
• In addition, experimental techniques based on X-ray tomography allow the non-invasive measurement of displacement fields inside three-dimensional deformable solids [63] .
• In such situations, the virtual work principle provides a direct link between the experimental data and the unknown quantities, which are usually related to constitutive properties, defects or damage.

4.1. Least-squares functionals

• In a variety of practical cases one can assume that the quantity to be identified (e.g. elastic moduli, boundary tractions of displacements, cracks.
• Moreover, the available experimental data is also often discrete in nature due to a variety of practical limitations.
• These measurements result from the application of a known excitation φ (e.g. a distribution of forces over the boundary).
• Regularized least-squares cost functions of the form EQUATION are widely used.
• Approach (iii), which replaces the solution of B derivative problems with that of one adjoint problem, is the most efficient when the direct problem is linear, while being applicable to other situations as well.

4.2. Adjoint solution for the identification of elastic moduli

• As an example, the adjoint state approach is formulated in this section for the problem of identifying distributions of elastic moduli from elastostatic boundary measurements.
• Let u = u[C] denote the displacement created in the body endowed with a trial elastic tensor C(x) by the application of known forces φ over S. Note that the format (55) allows for incomplete data, and even pointwise measurements, whereas RG functionals such as (32) require complete data.
• The parameter sensitivity formula (59) treats C(x) as the main unknown.
• Then, the constrained minimization of J (C) is recast as a saddle point search for L(u, w, C), characterized by the stationarity equations EQUATION For an assumed value of C, the first two stationarity equations are readily seen to define the direct problem (54) and the adjoint problem (61) .
• It is especially convenient when all equations (64) are linear, which may happen in other types of inverse problems amenable to the same type of formulation, as the inverse problem can then be solved by means of one calculation.

4.3. Indentation test and adjoint state for unilateral contact

• The indentation test is used for the determination of various types of mechanical parameters, in particular those associated with constitutive properties of materials.
• Even for linear elastic materials, the indentation curve is nonlinear due to the unilateral contact conditions and (depending on the properties of the sample surface) the frictional forces.
• The practical usefulness of the test lies in its simplicity and the fact that no specimens have to be prepared.
• Due to the non-penetration inequality constraint, the actual contact surface Γ C (t) is not known a priori.
• On differentiating the first-order optimality conditions ( 66) and (67c), the variation of the solution to the direct problem ( 66)-( 67) is found to be governed by the weak formulation EQUATION where, as in Section 4.2, δC = (∂ Ĉ/∂b).δb, and provided strict complementarity holds (i.e. that the inequalities are strict almost everywhere) in (67) .

5. Formulations based on the error in constitutive equation

• As mentioned earlier, many inverse problems in solid mechanics, and in other areas of physics as well, revolve around the identification of constitutive parameters, i.e. parameters appearing in the constitutive relations modelling the physical behaviour of materials.
• Among the key issues in that area is the identification of distributed parameters.
• The sets of all admissible displacements and stresses are defined e.g. by ( 22) and ( 23) with reference to equilibrium equations and boundary conditions, and admissible strains are then obtained from the compatibility equation (1) .
• ECE-based cost functions are involved in e.g. [141] (identification of permeability distributions) or [89] (identification of electrostatic conductivity distributions).
• Independently, the concept of ECE has first been introduced in elasticity by Ladevèze and Leguillon [101] in connection with error estimation in finite element computations.

5.1. Error in constitutive equation in elastostatics

• The variational principles of elasticity are very useful in defining cost functions based on the error in constitutive equation (ECE), which are well suited to the identification of distributed parameters.
• The main properties of such cost functions are now discussed in more detail and illustrated, in connection with elastostatics (section 5.1.1) and elastic vibrations (section 5.5).

5.1.1. ECE functionals.

• Let us again consider the problem of reconstructing the unknown distribution of elastic moduli C ⋆ from a series of overspecified boundary data pairs of displacements and traction (ξ, φ).
• One can introduce and motivate ECE functionals by considering the sum E(v, s, C) of the potential energy (20) and the complementary energy (21) for an assumed elasticity tensor C and a well-posed (i.e. not overspecified) set of boundary conditions of the form (12-13): EQUATION.
• In particular, in view of ( 78), the statement (25) means that an admissible pair (u, σ) solves an elastic equilibrium problem if and only if they are related through the elastic constitutive equation ( 8), as expected.
• These remarks suggest to introduce the following definition for E(C), the error in constitutive equation (ECE) functional for the identification of a distribution of elastic moduli: EQUATION in which an important difference with (76) is that overspecified boundary data of the form ( 27) is now involved.
• One can then easily show through integration by parts and algebraic manipulations that E(C) can be expressed in any of the following equivalent forms: EQUATION Note that E(C) expressed by either of the above formulae depends on C explicitly, but also implicitly through u D and u N .

5.2. Identification algorithm based on alternating directions

• The lack of mathematical results for the properties of E in elasticity notwithstanding, this functional has been used to define an alternating directions algorithm for the numerical reconstruction of distributed elastic moduli [48] , along the lines proposed in [89] for isotropic conductivities.
• The minimization of step (ii) with respect to the eigenmoduli c k is straightforward and yields EQUATION.
• It is therefore applicable only within certain a priori assumptions regarding material symmetries, so that the ζ k are known beforehand.

5.3.1. Fundamental equations for thin elastic plates. A thin plate is a planar solid occupying a domain of the form

• Only flexural motions of plates are considered here, leaving aside stretching motions (i.e. in-plane deformations) which are governed by equations similar to those of section 2.1.
• In the most general case, the bending of anisotropic elastic plates is therefore described by six independent moduli.

5.3.2. ECE functional.

• D on ∂ω where the above definition of C and S, the sets of kinematically and statically admissible fields, is consistent with the overdetermined nature of the boundary data and are the counterparts for plates of equations ( 22) and ( 23).
• Identifiability and uniqueness issues have been investigated by Ikehata [82, 83] .
• For anisotropic plates, he has shown [83] that flexural rigidity tensors belong to either of two classes, uniqueness being valid for one but not for the other.
• The general-purpose FEM code CAST3M has, as with the plane-strain example of Section 5.2, been used as a basis for the implementation.
• The synthetic measured displacement data consists in the deflection u at the location of the concentrated load and the deflection distribution on ∂ω.

5.4. Other ECE-based functionals in elastostatics 5.4.1. Modified ECE functional.

• There are many other ways of incorporating the error in constitutive equation into cost functions for the purposes of parameter identification.
• One natural extension of (79) consists of relaxing the satisfaction of overdetermined boundary data (ξ, φ) in the definition of spaces of admissible fields.
• Adjusting them by means of the L-curve method [79] has not been done so far but is certainly worthy of investigation.
• Simulated data for a square defect occupying the four finite elements located in rows (10 ,11) and columns (11 ,12) , and whose elastic moduli are µ def = 2µ and ν = 0.4, are computed for six parabolic distributions of normal loads applied at six different locations.
• The modified ECE functional appears to have better convexity in the vicinity of the "true" defect than J (C).

5.4.2. ECE as a penalty term.

• The major drawback of the least squares functionals are generally bad stability properties, in the sense that small data errors induce large errors in the solution.
• The classical technique for restoring stability is regularization, a subject covered by a vast literature (see e.g. [58, 138] ).
• The value of the minimum W C(b) (u) is problem-dependent, making this approach somewhat impractical.
• Accordingly, the primal-dual formulation [39] of the identification problem is then defined in terms of the cost function EQUATION EQUATION Formulations ( 99) and ( 95) are clearly linked, with coefficients (A, B) in ( 95) playing the same role as the penalty parameter η in (99).

5.5. Error in constitutive equation in dynamics

• Theoretical studies [12] [13] [14] show in particular that the knowledge of all eigenfrequencies for one set of boundary conditions is usually insufficient for a reconstruction of such characteristics, even assuming that the latter are in the form of a small unknown perturbation of a known reference value.
• On the other hand, mechanical structures of engineering interest are too complex for the reconstruction inverse problem to be amenable to the same kind of mathematical analysis.
• Structural FE model updating using measured vibrational data has been the subject of many investigations during the recent years.

5.5.1. ECE-based functional.

• Consider for definiteness the free vibrations of a linearly elastic solid occupying the domain Ω and endowed with the elasticity tensor C ⋆ (x) and mass density ρ ⋆ (x) distributions.
• Assume that imperfect values ρ(x), C(x) of the true distributions of mass density ρ ⋆ (x) and elasticity tensor C ⋆ (x) are known.
• Thus, a primal-dual formulation similar to (99) could conceivably be set up, with γ then becoming a penalty parameter allowed to assume arbitrarily large values.
• The first-order variation of L then reads: EQUATION where, for each measured frequency: EQUATION EQUATION.

5.6. Prior localization of defects

• This useful feature is supported by many numerical experiments [27, 41, 57, 102, 126] , exploiting either synthetic or real experimental data, and similar results are available for static problems [30] .
• This observation can, again, be given a partial formal explanation, as follows.
• This in particular corroborates the fact, generally observed in numerical experiments, that the geometrical localization using ECE density is better for stiffness defects than for mass defects.
• The a priori localization properties are in practice good even with few and/or inexact eigendisplacement measurements.
• Two illustrative example of the a priori localization capabilities of the ECE density are now presented.

Example 1 (synthetic data).

• The eigendisplacement field ξ associated with a vibrational mode has been computed, together with the distributed ECE, under the assumption that a complete measurement of ξ is available.
• Figure 6 presents the FE mesh and (through the distribution of Young moduli) the simulated defect, while the distribution of ECE over the mesh is depicted on Figure 7 .
• Thus, it is apparent that in the most favorable situation of a complete and exact measurement of the field ξ (disregarding the FE discretization error), the distribution of ECE provides a very accurate estimation of the geometrical support of the defect.

Example 2 (experimental data

• Eigendisplacements in the direction normal to the shell are measured at 48 sensor locations (4 rings of 12 equally spaced sensors located on the cylindrical part of the structure, respectively 22,5cm, 45cm, 67,5cm and 90cm above the disk-shaped bottom).
• A finite element model of the structure has been set up using the general-purpose code CAST3M, in which the cylindrical part is made of 4 identical rings of 24 eight-noded shell elements.
• The sensor locations coincide with mesh nodes.
• The figure 9 depicts the values of the ECE integrated over each finite element (normalized to max=1) obtained for the "small" and "large" defects (respectively located on elements 61-62 and 61-62-63-64), using experimental data for 13 eigenmodes (i.e. 13 frequencies and 13 × 48 displacements).

6. Geometrical sensitivity techniques for defect identification

• Defect identificaton problems are often formulated in terms of the minimization of a cost function featuring the experimental data and prior information.
• Traditional minimization methods are usually preferred to global search techniques such as evolutionary algorithms, the latter being in most cases infeasible due to the prohibitive computational cost of large numbers of direct elastic scattering solutions.
• To perform optimally, gradient-based optimization techniques are used in conjunction with shape sensitivity formulations.
• Geometrical sensitivity techniques are presented in connection with the problem of determining the shape and position of an unknown object (cavity or crack) embedded in the elastic body from elastodynamic experimental data.
• Time-harmonic conditions are assumed for simplicity, although a treatment in the time domain could have been presented as well.

Denoting by

• To establish shape sensitivity formulae for cost functions of the form (122) , one may either use an optimal control approach whereby a Lagrangian is introduced so as to incorporate the direct elastodynamic problem as an equality constraint, or exploit a reciprocity identity together with a direct differentiation of (122) with respect to the defect shape.
• One involves domain integrals and is therefore wellsuited to domain discretization method such as the finite element method (section 6.1.1), while the other features boundary integrals and is better suited to solution techniques based on boundary integral equation methods, which are frequently used for this kind of inverse problem (section 6.1.2).

6.1.1. Sensitivity formula in domain integral form.

• When applied to the computation of the gradient of J (Γ) it entails the solution of one derivative problem (126) for each geometrical parameter involved.
• The adjoint solution is defined as the particular elastodynamic state that fulfills the well-posed set of boundary conditions EQUATION.
• On substituting these boundary conditions into (127) and adding the resulting equation to equation (124) , one obtains the following sensitivity formula: EQUATION.

6.1.2. Sensitivity formula in boundary integral form.

• To ensure that formula (131) can actually be evaluated using only the boundary traces of the direct and adjoint solutions, the bilinear form σ[w] : ∇u must be expressed in terms of the surface gradients ∇.
• The evaluation of shape sensitivities using ( 131)-( 132) is very useful in cases where the direct and adjoint elastodynamic solutions are computed by means of the boundary element method.
• Such implementation is presented for three-dimensional elastodynamics in [78, 119].

6.2. Adjoint formulations for crack shape sensitivity

• Defect identification is often concerned with finding cracks.
• One is therefore led to enquire whether the sensitivity formulae ( 129) and ( 131) remain valid for crack identification problems and, if not, to establish valid substitute formulae.
• Besides, an examination of the steps leading to (129) reveals that all domain integrals involved in the derivation are integrable.
• This contradicts e.g. well-known results of fracture mechanics [65] showing that the mechanical potential energy of a cracked solid is affected by in-plane crack extensions.
• The functions K[w] I,II,III (s) are known as the stress intensity factors (SIFs).

6.3. Shape sensitivity formulae in the time domain

• Without going again through a detailed derivation, sensitivity formulae similar to ( 129), ( 131) and ( 135) can be established in connection with elastodynamics in the time domain.
• Then, the counterparts of the sensitivity formulae ( 129), ( 131) and ( 135) are respectively EQUATION EQUATION EQUATION where, in particular, the SIFs K[u Γ ] I etc. are now functions of both s and t.
• Note that the last integral in (122) , which is unaffected by whether time-harmonic or transient conditions are assumed, has been omitted in (136) for the sake of brevity.

6.4. Topological derivative

• These pitfalls can be circumvented to some extent by prior information, which, if in good agreement with the physical reality, may lead to improvements for both the choice of initial guess and the cost function properties.
• In [76] , this idea was considered in the context of inverse elastic scattering pertaining to semi-infinite and infinite domains, where the availability of suitable fundamental solutions made it possible to establish explicit expressions for T (x o ), while developments along similar lines for 2D elastostatics are presented in [66, 85] .
• In what follows, a procedure for the computation of the topological derivative applicable to arbitrary elastic bodies is presented.
• One may note that this definition is not restricted to spherical infinitesimal cavities (for which B is the unit ball, S the unit sphere and |B| = 4π/3).

6.4.3. Numerical examples in linear acoustics.

• To illustrate the foregoing notions, sample results from numerical experiments, for frequency-domain acoustics, are now presented.
• Similar experiments are currently being conducted for the preliminary imaging of cavities and elastic inclusions using elastodynamic data.
• Again, the values obtained for T (x o ) in a series of horizontal planes are consistent with the actual location of the "true" inclusion, despite the fact that the latter is of finite size while the asymptotic formula (153) only holds in the limit ε → 0.

7. Conclusion and further developments

• Several types of inverse problems arising in the linear theory of elasticity have been considered, revolving around the identification of material parameters, distributions of elastic moduli, and geometrical objects such as cracks and inclusions.
• In addition, the virtual work principle appears to be instrumental in establishing computationally efficient formulae for parameter or geometrical sensitivity, in connection to the adjoint solution method.
• Finally, the concept of topological derivative and its implementation for the identification of cavities or inclusions have been expounded.
• In particular, recent optical techniques now allow to measure displacements fields or in-plane strain fields on the boundary, and microtomographic techniques allow to measure displacements inside a body.
• Real-world mechanical systems feature complex geometries and limited data, whereas rigorous analyses of fundamental issues such as existence or uniqueness of solutions, or convergence properties of inversion algorithms are usually available only for idealized situations.

Figures

Figure 2. Identificatication of a distribution of Poisson ratio: “true” distribution used to create the synthetic data (left); identified distribution after 5 iterations (center) and after 32 iterations (right).

Figure 3. Identificatication of a distribution of Poisson ratio: “true” distribution used to create the synthetic data (left); identified distribution after 2 iterations for noisy data (center) and after 14 iterations for noise-free data (right).

Figure 1. Identification of the convex hull of a set of two straight cracks using the instantaneous reciprocity gap.

Figure 5. Contour maps of J (C) (top), H(C) with (A,B) = (50, 50) (middle) and H(C) with (A,B) = (1, 50) (bottom), plotted as functions of the location of the center of a trial four-element square defect. The graphs in the left and right columns correspond to no noise and a simulated relative noise of amplitude 1% on the synthetic data, respectively. The white square indicates the location of the “true” four-element defect.

Figure 6. Error in constitutive equation and free vibrations (synthetic example): FE mesh and distribution of Young modulus showing the simulated defect.

Figure 7. Error in constitutive equation (ECE) and free vibrations (synthetic example): FE mesh and distributed ECE computed from one perfectly known modal displacement.

Figure 4. Indentification of a distribution of the bending stiffness in an isotropic plate: “true” distribution used to create the synthetic data (left) and identified distribution after 20 iterations using the second eigenmodule (right).

Figure 12. Identification of a penetrable obstacle of ellipsoidal shape (center (−1, 2,−3), semiaxes 1, 0.6, 0.4 aligned with the global Cartesian frame, wave velocity contrast γ = 0.7, mass density contrast β = 0.45) embedded in an acoustic half-space: distribution of T in nine regularly spaced horizontal planes x3 = a with a = −1, . . . ,−5. The horizontal contour of the “true” penetrable obstacle is shown (white ellipse). The same color scale is used in all nine graphics. Case k = 4.

Figure 8. Description of the test structure and FEM model.

Figure 9. ECE distributions (normalized to max=1) for the “small” (elements 61-62, left) and “large” (elements 61-62-63-64, right) defects, using real data for 13 measured eigenmodes.

Figure 11. Identification of a penetrable obstacle of ellipsoidal shape (center (−1, 2,−3), semiaxes 1, 0.6, 0.4 aligned with the global Cartesian frame, wave velocity contrast γ = 0.7, mass density contrast β = 0.45) embedded in an acoustic half-space: distribution of T in nine regularly spaced horizontal planes x3 = a with a = −1, . . . ,−5. The horizontal contour of the “true” penetrable obstacle is shown (white ellipse). The same color scale is used in all nine graphics. Case k = 1.

Table 1. Constitutive parameters of some common isotropic elastic materials: E Young modulus, ν Poisson ratio, ρ mass density; σY and εY define the elastic limit, i.e. are the stress and strain levels beyond which the material is no longer elastic.

Figure 10. Identification of a set of three hard spherical objects (with respective centers (−3, 2,−3), (2, 3,−3) and (0,−4,−4) and radii 1, 1 and 2) in an acoustic half-space: distribution of T in the x3 = −3 (horizontal) plane. The horizontal contours of the three objects are shown (white circles).

Full text

Inverse problems in elasticity
Marc BONNET, Andrei CONSTANTINESCU
Laboratoire de M
´
ecanique des Solides (UMR CNRS 7649), Ecole Polytechnique, 91128
Palaiseau cedex
E-mail: bonnet@lms.polytechnique.fr,
andrei.constantinescu@lms.polytechnique.fr
Abstract. This article is devoted to some inverse problems arising in the context of linear
elasticity, namely the identification of distributions of elastic moduli, model parameters, or
buried objects such as cracks. These inverse problems are considered mainly for three-
dimensional elastic media under equilibrium or dynamical conditions, and also for thin
elastic plates. The main goal is to overview some recent results, in an effort to bridge
the gap between studies of a mathematical nature and problems defined from engineering
practice. Accordingly, emphasis is given to formulations and solution techniques which are
well suited to general-purpose numerical methods for solving elasticity problems on complex
configurations, in particular the finite element method and the boundary element method. An
underlying thread of the discussion is the fact that useful tools for the formulation, analysis and
solution of inverse problems arising in linear elasticity, namely the reciprocity gap and the error
in constitutive equation, stem from variational and virtual work principles, i.e. fundamental
principles governing the mechanics of deformable solid continua. In addition, the virtual
work principle is shown to be instrumental for establishing computationally efficient formulae
for parameter or geometrical sensitivity, based on the adjoint solution method. Sensitivity
formulae are presented for various situations, especially in connection with contact mechanics,
cavity and crack shape perturbations, thus enriching the already extensive known repertoire of
such results. Finally, the concept of topological derivative and its implementation for the
identification of cavities or inclusions are expounded.
Inverse Problems 21:R1–R50 (2005)

Inverse problems in elasticity 2
1. Introduction
Elasticity theory describes the reversible deformation of solid bodies subjected to excitations
of various physical natures: mechanical, thermal, electromagnetical etc. Such excitations,
applied as distributions over the body (e.g. gravitation, Lorentz forces, thermal expansion)
or over the boundary (pressure, contact forces), generate strains (i.e. local deformations) and
stresses (i.e. local forces) in the material. Elasticity is a mechanical constitutive property of
materials whereby (i) a one-to-one relationship between instantaneous strains and stresses on
the current deformed configuration is assumed, and (ii) the material reverts to its initial state
if the sollicitation history is reversed.
Almost all natural or manufactured solid materials have a deformation range within
which their mechanical behaviour can be modelled by elasticity theory. For sufficiently small
strains, the elastic behaviour is considered as linear, i.e. strains and stresses are assumed to
be proportional to each other. A vast body of engineering experience shows that the theory
of linear elasticity allows an accurate modelling of many man-made or natural objects: civil
engineering structures, transportation vehicles, machines, the earth mantle (to list just a few),
and provides an essential tool for analysis and design. In addition to the basic theory for three-
dimensional solid media, specialized approches have been developed for cases featuring two
or more dissimilar length scales: composite media, slender structures (beams, plates, shells).
The theory of linearized elasticity has developed into one of the now classical areas of
mathemematical physics. Equilibrium problems are governed by elliptic partial differential
equations, similar to those of electrostatics but more complex in that physical quantities of
interest are described by tensor fields rather than vector fields. Closed-form solutions are
available only for simple geometries (usually corresponding to separable coordinate systems),
so that most real-life modelling studies are based on numerical solution methods. Dynamic
conditions give rise to hyperbolic partial differential equations.
The main types of inverse problems that arise in the context of linear elasticity, and more
generally of the mechanics of deformable solids, are similar to those encountered in other
areas of physics involving continuous media and distributed physical quantities, e.g. acous-
tics, electrostatics and electromagnetism. They are usually motivated by the desire or need
to overcome a lack of information concerning the properties of the system (a deformable
solid body or structure). Mathematical and numerical techniques for the reconstruction of
buried objects of a geometrical nature, such as cracks, cavities or inclusions, are the subject of
many investigations, e.g. [3, 7, 8, 15, 26, 31, 32, 44–46, 61, 76, 95, 96, 118, 120]. Mechanical
waves, such as ultrasonic or Lamb waves, are also frequently used in practical non-destructive
testing of structures, see e.g. [103–105, 112, 131] and the references provided therein.
The identification of distributed parameters [5, 12–14, 29, 39, 48, 49, 60, 70, 81–83, 116] (e.g.
elastic moduli, mass density, wave velocity) arises in connection with e.g. medical imaging
of tissues [11] or seismic exploration [98, 114, 123, 130, 136, 140]. The reconstruction of
residual stresses [9, 67, 106, 127] is a related topic with important engineering implications.
Models of complex engineering structures often feature local parameters that are not known
with sufficient accuracy, and therefore need to be corrected by exploiting experimental

Inverse problems in elasticity 3
information on the mechanical response of the structure. Model updating is often treated
as an inverse problem [17, 27, 41, 57, 100, 107, 126, 142], in particular because corrections
affect distributed parameters over a limited region of space which is not known beforehand.
Identification of sources or inaccessible boundary values (i.e. Cauchy problems in elasticity)
are also encountered [42, 50, 55, 94, 109–111]. Finally, the identification of homogeneous
constitutive properties is increasingly often made on the basis of measurements taken on
structures, for which simplifying assumptions such as constant states of strain or stress are
invalid, and inverse techniques are then developed for that purpose [51, 64, 69, 108, 137].
This article is devoted to some of the above-mentioned inverse problems, namely the
identification of (i) distributions of elastic moduli, (ii) model parameters, (iii) buried cracks
or other geometrical objects. These inverse problems will be considered mainly for three-
dimensional elastic media under equilibrium or dynamical conditions, and also for thin elastic
plates. The main goal is to overview some recent results, in an effort to bridge the gap
between studies of a mathematical nature and problems defined from engineering practice.
Accordingly, emphasis will be given to formulations and solution techniques which are well
suited to general-purpose numerical methods for solving elasticity problems on complex
configurations, in particular the finite element method [28, 125] and the boundary element
method [23, 40]. The important role of the variational principles of elasticity, which in
particular provide the foundations of the above-mentioned numerical solution methods, will
be highlighted throughout this article. Additionnally, investigations of a more mathematical
nature will also be reviewed.
The article is organised as follows. An overview of basic theory and equations of linear
elasticity under the small strain hypothesis, including the virtual work principle and varia-
tional formulations of equilibrium problems, is presented in Section 2. Then, Section 3 des-
cribes strategies for the identification of distributions of elastic moduli or cracks exploiting
the virtual work principle as an observation equation, with emphasis on the reciprocity gap
concept. The virtual work principle is also an effective tool for setting up parameter sensitivity
analyses and computing gradients of cost functions associated to identification problems, as
shown in Section 4. Then, Section 5 is devoted to cost functions and parameter identification
techniques based on the error in constitutive equation (ECE). Formulations for both three-
dimensional bodies and plates are presented, and the ability of the energy density function
associated to the ECE to outline the geometrical support of defects is discussed. Finally,
Section 6 is devoted to geometrical sensitivity tools, based on adjoint solutions, for defect
identification. A concise formulation for cavity shape sensitivity is followed by more recent
results concerning crack shape sensitivity and a presentation of the topological derivative
associated with wave scattering in the limit of vanishingly small objects.
A brief review of the typical orders of magnitude involved in elastic solid bodies will
close this introductory section. Values of elastic constitutive parameters for common isotropic
materials are given in table 1, where the Young modulus E and the Poisson ratio ν are
defined in terms of the basic experiment performed by applying a traction force F to both
extremities of a cylindrical bar. Under this experiment, the axial length stretches from −ℓ
0
to ℓ while the cross-section shrinks from S
0
to S, and one sets E = (F/S
0
)/(ℓ/ℓ
0
−1) and

Inverse problems in elasticity 4
E (GPa) ν ρ (kg/m
3
) σ
Y
(GPa) ε
Y
aluminium 71 0.34 2.6 150 – 400 .002 – .006
steel 210 0.29 7.8 200 – 1600 .002 – .007
titanium 105 0.34 4.5 700 – 900 .006 – .009
marble 26 0.3 2.8 10 .0004
glass 60 0.2 - 0.3 2.5 – 2.9 1200 .02
Table 1. Constitutive parameters of some common isotropic elastic materials: E Young
modulus, ν Poisson ratio, ρ mass density; σ
Y
and ε
Y
define the elastic limit, i.e. are the
stress and strain levels beyond which the material is no longer elastic.
ν = −(
p
S/S
0
−1)/(ℓ/ℓ
0
−1). Linear elasticity is usually valid when strains are small (typical
magnitudes are 10
−3
or less), and is also restricted to stress levels below a certain threshold
σ
Y
beyond which irreversible constitutive properties, e.g. plasticity, set in (see Table 1).
Strains and displacements can be measured directly, e.g. using strain gages, whereas
stresses can only be measured indirectly. Classical strain gages are reliable up to 10
−6
but
offer only a “pointwise” measurement, in practice over an area of a few square millimeters.
Modern technology based on laser interferometry or image correlation techniques [21] are
reliable up to 10
−5
but allow measurements of practically continuous fields over extended
areas. Such experimental techniques yield rich experimental data and are therefore well-suited
to identification problems. The importance of the latter techniques is increasing as inversion
techniques specifically exploiting availability of field quantities (either on the boundary or
over part of the domain itself) become accessible.
2. Review of governing equations
2.1. Fundamental field equations for three-dimensional elasticity
The deformation of an elastic body, occupying in its undeformed state the region Ω ⊂ R
3
bounded by the surface S, is usually described in terms of a vector displacement field u(x, t)
(x ∈ Ω) which is such that the deformation process moves a small material element lying
at x to its new position x + u(x, t). The linearized elasticity theory [74] is established on
the assumption of small strains, namely |∇u(x, t)| ≪ 1. In that case, the changes in metric
induced by the deformation are described by the linearized strain tensor ε(x, t), defined as a
differential operator on u by:
ε[u](x, t) = (∇u(x, t) + ∇u
T
(x, t))/2 (1)
This equation is often referred to as the compatibility equation for small deformations. The
strain ε(x, t) is a symmetric second-order tensor.
The material is characterized by two constitutive parameters: its mass density distribution
ρ(x), associated with the kinetic energy
T (u) =
1
2
Z
Ω
ρ|
˙
u|
2
dV

Inverse problems in elasticity 5
(where the dot denotes time differentiation) and the fourth-order tensor of elastic moduli
C(x), hereafter referred to as the elasticity tensor, associated with the elastic strain energy
E(u) =
1
2
Z
Ω
ε[u]:C :ε[u] dV (2)
The elasticity tensor C defines a positive definite quadratic form over the 6-dimensional space
of symmetric second-order tensors. Therefore, C has the following symmetries:
C
ijkℓ
= C
kℓij
= C
jikℓ
(1 ≤ i, j, k, ℓ ≤ 3) (3)
and hence has at most 21 independent coefficients. In the simplest situation of isotropic
elasticity, C depends on only two independent moduli. For instance, C can be expressed in
terms of the Lam
´
e coefficients (λ, µ):
C
ijkℓ
= λδ
ij
δ
kℓ
+ µ(δ
ik
δ
jℓ
+ δ
jk
δ
iℓ
) (4)
Other commonly used elastic parameters are the Young modulus E and the Poisson ratio ν,
which are related to the Lam
´
e constants by
E =
µ(3λ + 2µ)
λ + µ
ν =
λ
2(λ + µ)
(5)
The stress tensor σ describes internal forces: the traction vector
p[n](x, t) = σ(x, t).n(x) (6)
is such that p[n] dS is the elementary force applied on a infinitesimal surface patch dS
of unit normal n located at x ∈ Ω. The fundamental balance equation of the dynamics of
deformable bodies (an extension of Newton’s second law to a small material element) is then:
div σ(x, t) + f(x, t) − ρ(x)
¨
u(x, t) = 0 (7)
where f (x, t) is a given distribution of body forces. The constitutive assumption of linearized
elasticity, adopted in this article, postulates that the stress tensor σ(x, t) depends linearly on
the linearized strain tensor, i.e.:
σ(x, t) = C(x):ε[u](x, t) (8)
On combining the three field equations (1), (7) and (8) and eliminating ε and σ, the
displacement field is found to be governed by the partial differential equation
[A
C
u](x, t) + f (x, t) − ρ(x)
¨
u(x, t) = 0 (9)
with the elasticity operator A
C
defined by:
A
C
u := div (C : ε[u]) = div (C : ∇u) (10)
(where the last equality stems from the constitutive symmetries (3)). Equation (9) is the
analog for linear elasticity of the hyperbolic linear wave equation. Besides, a well-posed
elastodynamic problem features initial conditions
u(x, 0) = u
0
(x)
˙
u(x, 0) = v
0
(x) (x in Ω) (11)

Frequently asked questions

Q1. What are the contributions mentioned in the paper "Inverse problems in elasticity" ?

This article is devoted to some inverse problems arising in the context of linear elasticity, namely the identification of distributions of elastic moduli, model parameters, or buried objects such as cracks. 

Q2. What future works have the authors mentioned in the paper "Inverse problems in elasticity" ?

Some directions for further work directly related to topics presented in this article include ( i ) a more systematic development and testing of the reciprocity gap approach as a method for the identification of cracks, ( ii ) the study of convexity properties of functionals based on the error in constitutive equations, and develop similar functional in connection with the identification of nonlinear constitutive properties and ( iii ) the development of topological sensitivity techniques for time-domain formulations and in refined forms ( in particular based on higher-order expansions with respect to defect size ) and their integration in defect identification strategies. 

Q3. What is the class of inverse problem of engineering interest?

The reconstruction of distributed parameters (such as the flexural stiffness or the mass density) of mechanical structures from vibrational data, i.e. measured values of eigenfrequencies and eigenmodal displacements, is a class of inverse problem of engineering interest, especially in connection with updating finite element (FE) models of mechanical structures, i.e. correcting FE models so that they agree best with measurements on the real structure. 

Q4. What is the main drawback of the least squares functionals?

The major drawback of the least squares functionals are generally bad stability properties, in the sense that small data errors induce large errors in the solution. 

Q5. What is the importance of the minimization of the shape of using such methods?

The minimization of J with respect to Γ using such methods needs in turn, for efficiency, the ability to evaluate the gradient of the functional J with respect to perturbations of the shape of Γ, in addition to J (Γ) itself. 

Q6. How many cases of applied forces were generated?

Synthetic data was generated by solving the plate bending problem for the true distribution of D for nine cases of applied forces, in this case concentrated loads applied at nine different locations. 

Q7. What is the topological derivative approach for inverse elastostatics?

In cases featuring only one characteristic length (one scatterer embedded in an unbounded medium an illuminated by a plane wave), the topological derivative approach essentially provides (up to a scaling factor) the lowest-order moment of the normalized scattering amplitude in the theory of low-frequency direct and inverse scattering [37, 53, 54]. 

Q8. What are some directions for further work related to the topics presented in this article?

Some directions for further work directly related to topics presented in this article include(i) a more systematic development and testing of the reciprocity gap approach as a method for the identification of cracks, (ii) the study of convexity properties of functionals based on the error in constitutive equations, and develop similar functional in connection with the identification of nonlinear constitutive properties and (iii) the development of topological sensitivity techniques for time-domain formulations and in refined forms (in particular based on higher-order expansions with respect to defect size) and their integration in defect identification strategies. 

Q9. What is the fundamental balance equation of the dynamics of deformable bodies?

The fundamental balance equation of the dynamics of deformable bodies (an extension of Newton’s second law to a small material element) is then:div σ(x, t) + f(x, t) − ρ(x)ü(x, t) = 0 (7)where f(x, t) is a given distribution of body forces. 

Q10. What is the need to keep the number of direct computations as low as possible?

The need to keep the number of direct computations as low as possible suggests instead to stick with classical gradient-based optimization algorithm. 

Q11. What is the role of the virtual work principle in determining the quantity of displacement fields?

In such situations, the virtual work principle provides a direct link between theexperimental data and the unknown quantities, which are usually related to constitutive properties, defects or damage. 

Q12. How many independent combinations of the elastic coefficients Cijkl(x) can be?

using ϕ[z] = ϕ(1)[z] in (41) for all y and all ϑ allows to find at most five independent linear combinations of the elastic coefficients δCijkℓ(x). 

Q13. What are the types of inverse problems in the linear theory of elasticity?

In this article, several types of inverse problems arising in the linear theory of elasticity have been considered, revolving around the identification of material parameters, distributions of elastic moduli, and geometrical objects such as cracks and inclusions. 

Q14. What is the topological derivative for acoustic scattering?

Numerical experiments based on the topological derivative in elastodynamics recently appeared in connnection with identification of cavities [26, 76] and of penetrable elastic inclusions [77].6.4.2. Topological derivative for acoustic scattering. 

Q15. What is the corresponding integral in (134)?

Note that the last integral in (122), which is unaffected by whether time-harmonic or transient conditions are assumed, has been omitted in (136) for the sake of brevity. 

Q16. What is the topological derivative for the case of a hard obstacle?

For the simplest case where B is the unit sphere, one hasDij(B, β, γ) = 3(1 − β)2 + β δijThe limiting situation β = 0 in (154) yields the expression of the topological derivative for the case of a hard (i.e. rigid) obstacle of vanishing size ε.6.4.3. 

Q17. What are the techniques used to improve the effectiveness of the virtual field method?

To enhance the effectiveness and efficiency of this technique, a variety of techniques have been proposed for the construction of families of virtual fields w(x) tailored for specific classes of constitutive parameter identification problems.