Identifying the stored energy of a hyperelastic structure by using an attenuated Landweber method
We consider the nonlinear, inverse problem of identifying the stored energy function of a hyperelastic material from full knowledge of the displacement field as well as from surface sensor measurements. The displacement field is represented as a solution of Cauchy’s equation of motion, which is a nonlinear, elastic wave equation. Hyperelasticity means that the first Piola-Kirchhoff stress tensor is given as the gradient of the stored energy function. We assume that a dictionary of suitable functions is available and the aim is to recover the stored energy with respect to this dictionary. The considered inverse problem is of vital interest for the development of structural health monitoring systems which are constructed to detect defects in elastic materials from boundary measurements of the displacement field, since the stored energy encodes the mechanical peroperties of the underlying structure. In this article we develope a numerical solver for both settings using the attenuated Landweber method. We show that the parameter-to-solution map satisfies the local tangential cone condition. This result can be used to prove local convergence of the attenuated Landweber method in case that the full displacement field is measured. In our numerical experiments we demonstrate how to construct an appropriate dictionary and show that our algorithm is well suited to localize damages in various situations.
Key words. stored energy function, Cauchy’s equation of motion, hyperelasticity, conic combination, attenuated Landweber method, local tangential cone condition
AMS subject classifications. 35L70, 65M32, 74B20
The starting point of our inverse problem is Cauchy’s equation of motion for an hyperelastic material
where denotes time and a point in a bounded, open domain in . Furthermore, denotes the mass density in and is an external body force. The function , with , is called stored energy function and encodes the mechanical properties of the material. The derivative is to be understood componentwise and represents the displacement gradient in and . We refer to standard textbooks like [11, 17, 23] for detailed introductions and derivations of Cauchy’s equation.
Given initial values and and assuming that we have homogeneous Dirichlet boundary values we end up with
for and with
Equation (LABEL:cauchy-hyper) describes the behavior of hyperelastic materials. An example for hyperelastic materials are carbon fibre reinforced composites (CFC). This is why the inverse problem of identifying from measurements of the displacement field can be very important for the detection of defects in such structures, the so called structural health monitoring, see . Structural health monitoring systems consist of a number of actors and sensors which are applied to the structure’s surface. The idea is that the actors generate guided waves propagating through the structure which interact with defects and are subsequently acquired by the sensors. Mathematically this can be described as an inverse problem of determining material properties from boundary data of the displacement field which is a solution of (LABEL:cauchy-hyper). In this article we investigate the reconstruction problem of the stored energy function where we consider two different scenarios for the data acquisition. On the one hand we assume to have as data the full displacement field available and on the other hand we suppose that the data are measured on parts of the boundary. The second scenario is used as mathematical model for sensor measurements.
Because there are numerous publications on inverse identification problems in elastic media for different settings, we only summarize here articles which are associated with the scope of this article. A comprehensive overview of various inverse problems in the field of elasticity offers the article . In  Bourgeois and others have applied and implemented the linear sampling method, introduced by Colton and Kirsch in  for detection of reverberant scatterers, for the isotropic Navier Lamé equation. Because the method represents a possibility to detect defects in isotropic materials and damages, which are described by such a scatterer, it was used in  for the identification of cracks. The inverse problem on determining the spatial component of the source term in a hyperbolic equation with time-dependent principal part is investigated in . For solving this problem numerically the authors adopt the classical Tikhonov regularization to transform the inverse problem into an output least-squares minimation that can be solved by the iterative thresholding algorithm. In  the authors developed an algorithm for the quantitative reconstruction of constitutive parameters, namely the two eigenvalues of the elasticity tensor, in isotropic linear elasticity from noisy full-field measurements. An algorithm, that guarantees the conservation of the total energy as well as the conservation of momentum and angular momentum is found in . The reconstruction of an anisotropic elasticity tensor from a finite number of displacement fields for the linear, stationary elasticity equation is represented in . Lechleiter and Schlasche considered in  the identification of the Lamé parameters of the second order elastic wave equation from time-dependent elastic wave measurements at the boundary. For this reason they dispose an inexact Newton iteration validating the Fréchet differentiability of the parameter-to-solution map in terms of . A semi-smooth Newton iteration for the same problem was implemented in  also delivering an expression for the Fréchet derivative. The topic of our considerations is the identification of spatially variable stored energy functions from time-dependent boundary data which has not been investigated so far, to the best of our knowledge. In  an algorithm for defect localization in fibre-reinforced composites from surface sensor measurements was proposed using the equations of linear elastodynamics as mathematical model. The key idea of the method is the interpretation of defects as if they were induced by an external volume force. In some sense this article can be seen as the foundation of the method presented here.
We want to specify the inverse problems to be investigated in this article. Inspired by Kaltenbacher and Lorenzi  and following the authors of [27, 31] we assume to have a dictionary consisting of appropriate functions , , given such that
for and . We consider at first the following inverse problem.
(IP I) Given as well as the displacement field for and , compute the coefficients , such that satisfies the initial boundary value problem (IBVP) (LABEL:cauchy-hyper-comb), (LABEL:anfangswert1)–(LABEL:randwert).
Denoting by the forward operator, which maps a vector to the unique solution of the IBVP (LABEL:cauchy-hyper-comb), (LABEL:anfangswert1)–(LABEL:randwert) for fixed, then (IP I) is represented by the nonlinear operator equation
Here, denotes the domain of to be specified in Section 3. In that case we assume to have the full knowledge of the displacemt field. We solve this problem numerically. Since in practical applications measurements usually are only acquired at the structure’s surface we consider a further inverse problem,
where is the observation operator, which maps the solution of (LABEL:cauchy-hyper-comb) to measured data. Thus the observation operator includes the measurement modalities to our mathematical model. Following the article  we want to model having sensors in mind that average the displacement field on small parts of the boundary of the structure . For this reason we define by
where is the number of sensors and , , are weight functions that display the localization of the particular sensors. We assume that the functions for all have small support on . For more information about this definition of the observation operator see .
The second inverse problem is defined as follows.
(IP II) Given as well as data , compute the coefficients , such that satisfies (LABEL:allgGlgQ).
Outline. Section 2 provides all mathematical ingredients and tools which are necessary to deduce the results of the article. In particular we summarize an existing uniqueness result for the solution of the IBVP (LABEL:cauchy-hyper-comb), (LABEL:anfangswert1)–(LABEL:randwert) (Theorem LABEL:theorem21imp) and results of the article . Section 3 describes the implementation of a numerical solver the inverse problems (IP I), (IP II) using the attenuated Landweber method. In Section 4 we prove the local tangential cone condition for (IP I) and relying on that a convergence result for the attenuated Landweber iteration when applied to (IP I). Numerical results for (IP I) and (IP II) are subject of Section 5. Section 6 concludes the article.
2 Setting the stage
For the investigations and proofs in this article we need at first some results the most of which are proven in  and . The first one is a uniqueness result for the IBVP (LABEL:cauchy-hyper-comb), (LABEL:anfangswert1)–(LABEL:randwert) from .
In the following we assume that the conditions and are valid for the function for all and . Furthermore we restrict the nonlinearity of all functions and hence of by supposing, that there exist positive constants for , such that
hold for all and for almost everywhere. Let be the Frobenius norm induced by the inner product of matrices
In addition we require the existence and boundedness of higher derivatives of with respect to . More precisely we assume, that there are constants for with
for and . Additionally we require
for all , which holds true if e.g. . We suppose that the mapping is three times continuously differentiable for
Furthermore, the set of admissible coefficient vectors of the conic combination (LABEL:con-comb) is supposed to be restricted by assuming
It is easy to see that this set is coupled to the nonlinearity conditions of (LABEL:bed1)–(LABEL:beschr6) via the constants for .
After all we define the following set of admissible solutions of IBVP (LABEL:cauchy-hyper)–(LABEL:randwert). Let be given the constants , . Then we set
Let us notice that this set is only a subset of the set of admissible solutions mentioned in  and  because of the additional condition
for all . holds true, if e.g. , , , and are sufficiently smooth.
Now all necessary conditions are mentioned to prove the following uniqueness result for the solution of the IBVP (LABEL:cauchy-hyper-comb), (LABEL:anfangswert1)–(LABEL:randwert) for given , which has been presented in .
Theorem 2.1 ([31, Theorem 2.1])
Let , be two solutions to the initial boundary value problem (LABEL:cauchy-hyper-comb), (LABEL:anfangswert1)–(LABEL:randwert) corresponding to the parameters, initial values and right-hand sides and , respectively. Furthermore, assume that . If, in addition, the condition
is satisfied for
and if there are constants and , so that
then there exist constants , and , such that the stability estimate
is valid for all . Thereby, the constants , and only depend on , , , , ,
where is a constant, whose existence is ensured by inequality (LABEL:dim). The constant is defined by the continuity of the embedding ,
for all . Moreover, the constants , and are uniformly bounded, if we take with bounded.
The function is positive and bounded in the following way because of the non-negativity of the coefficients :
Next we prove an estimate being necessary for the proof of Theorem LABEL:kegelbed and following from Theorem LABEL:theorem21imp.
Let be two solutions of the problem (LABEL:cauchy-hyper) - (LABEL:randwert). Then there is for and all
Proof. Using the condition it follows by the triangle inequality
for all . Hence we have and therefore
for all and . Finally we obtain the estimate
and hence the assertion of the corollary.
Next we define for the remainder of the article the following spaces
and identify with its dual space . Then we obtain the Gelfand triple
with dense, continuous embeddings. Furthermore let be
for all and thereby
with dense, continuous embeddings and . Then the Poincaré inequality yields that there is a constant with
for all . Besides it follows directly from the definition of the norms in and
for all .
For the proof of the local tangential cone condition we need Gronwall’s lemma.
Lemma 2.3 (Gronwall’s lemma)
Let and be non-negative functions. If satisfies
for all with constants and , then
is valid for all .
A proof of this version is given in . The following results in connection with the Gâteaux/Fréchet derivative of the forward operator and its adjoint operator are proven in . That is the reason why there are no proofs in the remainder of this section. For this purpose let be the domain of the forward operator defined by
At first we characterize the Gâteaux derivative of .
Let be an interior point of . The Gâteaux derivative of the solution operator fulfills for the following linear system of differential equations with homogeneous initial and boundary value conditions
for and ,
Here, we used the notations
It was proven in  that the IBVP (LABEL:gateaux-diff)–(LABEL:vRandwert) has a unique solution and hence the Gâteaux derivative is well defined.
The IBVP (LABEL:gateaux-diff)–(LABEL:vRandwert) has a unique, weak solution in .
The aim of  was to show, that even is Fréchet differentiable.
To this end it was proven that the mapping , , is linear and bounded.
It is linear because (LABEL:gateaux-diff) is linear in and . The continuity was subject of the following theorem.
Theorem 2.6 ([28, Prop. 3.4])
Adopt the assumption of Lemma LABEL:frechet-system. The Gâteaux derivative is continuous in for all , i.e. there is a constant with .
Next we state the Fréchet differentiability of .
Theorem 2.7 ([28, Th. 3.8])
Adopt the assumptions of Lemma LABEL:frechet-system and let . There is a constant depending only on , and such that
We continue by recapitulating the representation of the adjoint operator of the Fréchet derivative , which was also proven in . We will see that the adjoint is important when applying iterative solvers as the Landweber method to the inverse problem . Let be the space consisting of all solutions of
for , if we define the mapping by
for all . Then we have and the mapping is bijective
since (LABEL:Bveq) with (LABEL:vAnfangswert) and (LABEL:vRandwert) is uniquely solvable according to Theorem LABEL:gateaux-ex. This yields , , where solves (LABEL:Bveq) with (LABEL:vAnfangswert) and (LABEL:vRandwert).
Finally endowed with the norm turns into a Hilbert space, which is a closed subspace of . The following lemma states that the embedding even
The embedding is continuous, i.e. there is a constant not depending on such that
A representation of the adjoint operator for fixed is presented in the following theorem.
Let be fixed and . The adjoint operator of the Fréchet derivative is given by
where is the weak solution of the hyperbolic, backward IBVP
We have now all ingredients together for the proof of the convergence result in Section LABEL:sec:convergence and the numerical solution of (IP I). The last point to consider in this section is the observation operator, which is needed for (IP II). Let be the observation operator with
where is the trace operator and given weight functions (see Section LABEL:sec:introduction). Then represents the data space. Because of (LABEL:beobachtung) it is easy to see that is linear in . Furthermore it is proven in  that is continuous in and that the adjoint operator has for all the representation
3 Numerical scheme
In this section we outline a numerical solution scheme for computing a solution of (LABEL:allgGlg) with input data . We want to use the attenuated Landweber method (see for example ,  or ). This iteration is defined for solving (IP I) via
and for solving (IP II) via
with initial value and relaxation parameter
We always denote by the closed ball centered about with radius . In (LABEL:Landweberalpha) respectively (LABEL:LandweberalphaQ) we can see that we have to solve in every iteration step of the attenuated Landweber method respectively once the forward problem (LABEL:cauchy-hyper-comb), (LABEL:anfangswert1)–(LABEL:randwert) and the adjoint problem (LABEL:backward-1), (LABEL:backward-2)–(LABEL:backward-3). So we need at first an algorithm for solving the forward problem. Because of the second time derivative in the differential equation we split initially the differential equation and then we get with for all the equivalent formulation
for all . We discretized this equation system at first in time with the -method with . Let be equidistantly partitioned in equal time steps of length and points , . Then we get with for all in every time step for all the following formulation of the time discretized equations
After some reformulations we get the system
It is easy to see that the first equation is not linear in and the second one is linear in . That is the reason why we have to implement a nonlinear solver for the first equation. Then we can discretize in space and solve both equations. For using the Newton method to solve the first equation we define