Identifying the stored energy of a hyperelastic structure by using an attenuated Landweber method

# Identifying the stored energy of a hyperelastic structure by using an attenuated Landweber method

## Abstract

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.

s

tored energy function, Cauchy’s equation of motion, hyperelasticity, conic combination, attenuated Landweber method, local tangential cone condition

{AMS}

35L70, 65M32, 74B20

## 1 Introduction

The starting point of our inverse problem is Cauchy’s equation of motion for an hyperelastic material

 ρ(x)¨u(t,x)−∇⋅∇YC(x,Ju(t,x))=f(t,x),

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

 ρ(x)¨u(t,x)−∇⋅∇YC(x,Ju(t,x))=f(t,x) (1)

for and with

 u(0,⋅) = u0∈H2(Ω,R3), (2) ˙u(0,⋅) = u1∈H1(Ω,R3) (3)

and

 u(t,x)=0,t∈[0,T],x∈∂Ω. (4)

Equation (1) 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 [15]. 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 (1). 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 [8]. In [9] Bourgeois and others have applied and implemented the linear sampling method, introduced by Colton and Kirsch in [12] 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 [10] 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 [18]. 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 [2] 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 [29]. The reconstruction of an anisotropic elasticity tensor from a finite number of displacement fields for the linear, stationary elasticity equation is represented in [3]. Lechleiter and Schlasche considered in [22] 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 [21]. A semi-smooth Newton iteration for the same problem was implemented in [7] 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 [6] 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 [19] and following the authors of [27, 31] we assume to have a dictionary consisting of appropriate functions , , given such that

 C(x,Y)=N∑K=1αKCK(x,Y) (5)

with nonnegative constants , . In that way the mentioned dictionary should consist of physically meaningful elements such as polyconvex functions, see [4, 17]. The wave equation then is then given as

 ρ(x)¨u(t,x)−N∑K=1αK∇⋅[∇YCK(x,Ju(t,x))]=f(t,x) (6)

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) (6), (2)–(4).
Denoting by the forward operator, which maps a vector to the unique solution of the IBVP (6), (2)–(4) for fixed, then (IP I) is represented by the nonlinear operator equation

 T(α)=~u. (7)

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,

 QT(α)=~y, (8)

where is the observation operator, which maps the solution of (6) to measured data. Thus the observation operator includes the measurement modalities to our mathematical model. Following the article [6] 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

 Q[u](t)=(∫∂Ω⟨gk,u(t)⟩R3dξ)k=1,...,l,

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 [6]. The second inverse problem is defined as follows.
(IP II) Given as well as data , compute the coefficients , such that satisfies (8).
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 (6), (2)–(4) (Theorem 2) and results of the article [28]. 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 [31] and [28]. The first one is a uniqueness result for the IBVP (6), (2)–(4) from [31].

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

 κ[0]K∥Y∥2F≤CK(x,Y)≤μ[0]K∥Y∥2F (9)

and

 κ[1]K∥H∥2F≤⟨⟨H|∇Y∇YCK(x,Y)H⟩⟩≤μ[1]K∥H∥2F (10)

hold for all and for almost everywhere. Let be the Frobenius norm induced by the inner product of matrices

 ⟨⟨A|B⟩⟩:=tr(A⊤B)for A,B∈R3×3.

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

 ∥∂Ypq∂Yij∂YklCK∥L∞(Ω×R3×3)≤μ[2]K (11)
 ∥∂Yab∂Ypq∂Yij∂YklCK∥L∞(Ω×R3×3)≤μ[3]K (12)
 ∥∂l∂YklCK∥L∞(Ω×R3×3)≤μ[4]K (13)
 ∥∂Yij∂l∂YklCK∥L∞(Ω×R3×3)≤μ[5]K (14)
 ∥∂l∂Yij∂YklCK∥L∞(Ω×R3×3)≤μ[6]K (15)
 ∥∂Ypq∂l∂Yij∂YklCK∥L∞(Ω×R3×3)≤μ[7]K (16)

for and . Additionally we require

 ∂Yij∂l∂YklC(x,Y)=∂l∂Yij∂YklC(x,Y) (17)

for all , which holds true if e.g. . We suppose that the mapping is three times continuously differentiable for almost everywhere.
Furthermore, the set of admissible coefficient vectors of the conic combination (5) is supposed to be restricted by assuming

 α∈C((κ[a])a=1,2,(μ[b])b=1,...,7) := {α∈(0,∞)N:∑NK=1αKκ[a]K≥κ[a],∑NK=1αKμ[b]K≤μ[b] for alla=1,2 and b=1,...,7}.

It is easy to see that this set is coupled to the nonlinearity conditions of (9)–(16) via the constants for .
After all we define the following set of admissible solutions of IBVP (1)–(4). Let be given the constants , . Then we set

 A := A(M0,M1,M2,M3,M4) (18) := ⎧⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪⎩u∈L∞((0,T)×Ω,R3)∩W1,∞((0,T),H1(Ω,R3)):∥∂l∂ju∥L∞((0,T),L2(Ω,R3))≤M0,∥∂l˙uk∥L∞((0,T)×Ω)≤M1,∥∂l∂j˙uk∥L∞((0,T)×Ω)≤M2,∥∂l∂juk∥L∞((0,T)×Ω)≤M3,∥∂luk∥L∞((0,T)×Ω)≤M4for all% i,j,k,l=1,2,3⎫⎪ ⎪ ⎪ ⎪ ⎪⎬⎪ ⎪ ⎪ ⎪ ⎪⎭.

Let us notice that this set is only a subset of the set of admissible solutions mentioned in [31] and [28] 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 (6), (2)–(4) for given , which has been presented in [31].

{theorem}

[[31, Theorem 2.1]] Let , be two solutions to the initial boundary value problem (6), (2)–(4) corresponding to the parameters, initial values and right-hand sides and , respectively. Furthermore, assume that . If, in addition, the condition

 78μ<κ<98μ (19)

is satisfied for

 κ:=N∑K=1αKκ[1]Kandμ:=N∑K=1αKμ[1]K (20)

and if there are constants and , so that

 κ≥κ(α)>0andμ≤μ(α), (21)

then there exist constants , and , such that the stability estimate

 [ρ∥(˙u−˙¯u)(t,⋅)∥2L2(Ω,R3)+κ(α)∥(Ju−J¯u)(t,⋅)∥2L2(Ω,R3×3) +ρ∥(¨u−¨¯u)(t,⋅)∥2L2(Ω,R3)+κ(α)∥(J˙u−J˙¯u)(t,⋅)∥2L2(Ω,R3×3) +∥(u−¯u)(t,⋅)∥2H2(Ω,R3)]12 ≤¯C0[μ(α)∥u0−¯u0∥2H2(Ω,R3)+∥u1−¯u1∥2H1(Ω,R3)]12 +¯C1∥f−¯f∥W1,1((0,T),L2(Ω,R3))+¯C2∥α−¯α∥∞

is valid for all . Thereby, the constants , and only depend on , , , , ,

 ¯C(α):=N∑K=1αKμ[2]K(N∑K=1αKκ[1]K)−1 (22)

and

 ^C(α):=^K1−√1−ϵN∑K=1αKμ[1]K(N∑K=1αKκ[1]K)−2, (23)

where is a constant, whose existence is ensured by inequality (19). The constant is defined by the continuity of the embedding ,

 ∥g∥H2(Ω,R3)≤^K∥g∥H2,10(Ω,R3)=^K(3∑k=1∫Ω3∑l=13∑j=1(∂i∂jgk(x))2dx)12

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 :

 0<ζ:=min1≤K≤Nμ[2]Kmax1≤K≤Nκ[1]K≤¯C(α)≤max1≤K≤Nμ[2]Kmin1≤K≤Nκ[1]K=:η<∞. (24)

Next we prove an estimate being necessary for the proof of Theorem 4 and following from Theorem 2. {corollary} Let be two solutions of the problem (1) - (4). Then there is for and all

 3∑i,j=1∫Ω|∂j~ui(t,x)|4dx≤(M4)2∥~u(t,⋅)∥2H10(Ω,R3). (25)
{proof}

Using the condition it follows by the triangle inequality

 ∥∂j~ui∥L∞((0,T)×Ω)=∥∂jui−∂j¯ui∥L∞((0,T)×Ω)≤∥∂jui∥L∞((0,T)×Ω)+∥∂j¯ui∥L∞((0,T)×Ω)≤2M4

for all . Hence we have and therefore

 Missing or unrecognized delimiter for \bigg

for all and . Finally we obtain the estimate

 Missing or unrecognized delimiter for \bigg ≤ 3∑i,j=1(M4)4∫Ω(|∂j~ui(t,x)|M4)2dx=(M4)2∥~u(t,⋅)∥2H10(Ω,R3)

and hence the assertion of the corollary.

Next we define for the remainder of the article the following spaces

 V:=H1(Ω,R3) and H:=L2(Ω,R3)

and identify with its dual space . Then we obtain the Gelfand triple

 V⊂H=H′⊂V′

with dense, continuous embeddings. Furthermore let be

 U:=H10(Ω,R3)

with

 ∥u∥U:=∥Ju∥L2(Ω,R3×3)

for all and thereby

 U⊂V⊂H=H′⊂V′⊂U′

with dense, continuous embeddings and . Then the Poincaré inequality yields that there is a constant with

 ∥u∥V≤√1+CΩ∥u∥U (26)

for all . Besides it follows directly from the definition of the norms in and

 ∥u∥H≤∥u∥V (27)

for all .
For the proof of the local tangential cone condition we need Gronwall’s lemma.

{lemma}

[Gronwall’s lemma] Let and be non-negative functions. If satisfies

 ψ(τ)≤a+τ∫0b(t)ψ(t)dt+τ∫0k(t)ψ(t)pdt

for all with constants and , then

 ψ(τ)≤exp(τ∫0b(t)dt)[a1−p+(1−p)τ∫0k(t)exp((p−1)t∫0b(σ)dσ)dt]1/(1−p) (28)

is valid for all .

A proof of this version is given in [1]. The following results in connection with the Gâteaux/Fréchet derivative of the forward operator and its adjoint operator are proven in [28]. 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

 D(T):={α∈RN+:the IBVP (% ???), (???)--(???) has a unique % solution u∈A}. (29)

At first we characterize the Gâteaux derivative of .

{lemma}

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

 ρ¨v(t,x)−∇⋅[∇Y∇YCα(x,Ju(t,x)):Jv(t,x)]=∇⋅[∇YCh(x,Ju(t,x))] (30)

for and ,

 v(0,x)=˙v(0,x)=0for x∈Ω (31)

and

 v(t,x)=0for x∈∂Ω. (32)

Here, we used the notations

 Cα=N∑K=1αKCK respectively Ch=N∑K=1hKCK.

It was proven in [28] that the IBVP (30)–(32) has a unique solution and hence the Gâteaux derivative is well defined.

{theorem}

The IBVP (30)–(32) has a unique, weak solution in .

The aim of [28] 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 (30) is linear in and . The continuity was subject of the following theorem.

{theorem}

[[28, Prop. 3.4]] Adopt the assumption of Lemma 2. 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}

[[28, Th. 3.8]] Adopt the assumptions of Lemma 2 and let . There is a constant depending only on , and such that

 ∥u(α+h)−u(α)−v∥L2(0,T;V)≤L2∥h∥32∞for ∥h∥∞→0.\vspace2mm (33)

We continue by recapitulating the representation of the adjoint operator of the Fréchet derivative , which was also proven in [28]. 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

 Bv=f (34)

for , if we define the mapping by

 Bv:=ρ¨v−∇⋅[∇Y∇YCα(x,Ju):Jv] (35)

for all . Then we have and the mapping is bijective since (34) with (31) and (32) is uniquely solvable according to Theorem 2. This yields , , where solves (34) with (31) and (32). Finally endowed with the norm turns into a Hilbert space, which is a closed subspace of . The following lemma states that the embedding even is continuous.

{lemma}

The embedding is continuous, i.e. there is a constant not depending on such that

 ∥v∥L2(0,T;U)∩H1(0,T;H)≤C∥v∥X,v∈X.

A representation of the adjoint operator for fixed is presented in the following theorem.

{theorem}

Let be fixed and . The adjoint operator of the Fréchet derivative is given by

 T′(α)∗w=[−T∫0∫Ω∇YCK(x,Ju(t,x)):J(B−1)∗w(t,x)dxdt]K=1,...,N∈RN, (36)

where is the weak solution of the hyperbolic, backward IBVP

 ρ¨p(t,x)−∇⋅[∇Y∇YCα(x,Ju(t,x)):Jp(t,x)]=w(t,x), (37) p(T,x)=˙p(T,x)=0,x∈Ω, (38) p(t,ξ)=0,(t,ξ)∈[0,T]×∂Ω. (39)

We have now all ingredients together for the proof of the convergence result in Section 4 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

 Q[u](t)=(∫∂Ω⟨gk,u(t)⟩R3dξ)k=1,...,l=(⟨gk,Γ(u(t))⟩L2(∂Ω,R3))k=1,...,l∈Rl, (40)

where is the trace operator and given weight functions (see Section 1). Then represents the data space. Because of (40) it is easy to see that is linear in . Furthermore it is proven in [6] that is continuous in and that the adjoint operator has for all the representation

 Q∗a=l∑k=1akΓ∗gk. (41)

## 3 Numerical scheme

In this section we outline a numerical solution scheme for computing a solution of (7) with input data . We want to use the attenuated Landweber method (see for example [13], [25] or [20]). This iteration is defined for solving (IP I) via

 α(k+1)=α(k)−ωT′(α(k))∗(T(α(k))−y),k=0,1,2,... (42)

and for solving (IP II) via

 α(k+1)=α(k)−ωT′(α(k))∗Q∗(QT(α(k))−y),k=0,1,2,... (43)

with initial value and relaxation parameter

 ω∈(0,1C2) (44)

with

 C:=sup{∥T′(α)∥:α∈Bρ(α(0))}. (45)

We always denote by the closed ball centered about with radius . In (42) respectively (43) we can see that we have to solve in every iteration step of the attenuated Landweber method respectively once the forward problem (6), (2)–(4) and the adjoint problem (37), (38)–(39). 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

 {˙u(t,x)−r(t,x)=0ρ˙r(t,x)−∇⋅∇YC(x,Ju(t,x))=f(t,x) (46)

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

 ⎧⎨⎩uj−uj−1k=θrj+(1−θ)rj−1ρrj−rj−1k=∇⋅∇YC(x,θJuj+(1−θ)Juj−1)+θfj+(1−θ)fj−1.

After some reformulations we get the system

 ⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩uj=uj−1+krj−1+k2θρ∇⋅∇YC(x,θJuj+(1−θ)Juj−1)+k2θ2ρfj+k2θ(1−θ)ρfj−1rj=rj−1+kρ∇⋅∇YC(x,θJuj+(1−θ)Juj−1)+kθρfj+k(1−θ)ρfj−1. (47)

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

 F(ujl) = ujl