James E Warner

Manuel I Diaz

Wilkins Aquino

Marc Bonnet
email: mbonnet@ensta.fr

Inverse material identification in coupled acoustic-structure interaction using a modified error in constitutive equation functional

Keywords: Parameter estimation, acoustic-structure interaction, error in constitutive equation, inverse problem

This work focuses on the identification of heterogeneous linear elastic moduli in the context of frequency-domain, coupled acoustic-structure interaction (ASI), using either solid displacement or fluid pressure measurement data. The approach postulates the inverse problem as an optimization problem where the solution is obtained by minimizing a modified error in constitutive equation (MECE) functional. The latter measures the discrepancy in the constitutive equations that connect kinematically admissible strains and dynamically admissible stresses, while incorporating the measurement data as additional quadratic error terms.

We demonstrate two strategies for selecting the MECE weighting coefficient to produce regularized solutions to the ill-posed identification problem: 1) the discrepancy principle of Morozov, and 2) an error-balance approach that selects the weight parameter as the minimizer of another functional involving the ECE and the data misfit. Numerical results demonstrate that the proposed methodology can successfully recover elastic parameters in 2D and 3D ASI systems from response measurements taken in either the solid or fluid subdomains. Furthermore, both regularization strategies are shown to produce accurate reconstructions when the measurement data is polluted with noise. The discrepancy principle is shown to produce nearly optimal solutions, while the error-balance approach, although not optimal, remains effective and does not need a priori information on the noise level.

Introduction

The noninvasive characterization of material properties in a physical system is of great importance in a variety of science and engineering fields. Along these lines, considerable research efforts have been made to formulate and solve inverse problems in which experimental measurements of the mechanical response of a system are used to infer its defining material parameters. Such parameter estimation problems are prevalent in areas like damage detection in structures, geotechnical exploration, biomechanical imaging, etc. [START_REF] Oberai | Evaluation of the adjoint equation based algorithm for elasticity imaging[END_REF][START_REF] Levasseur | Soil parameter identification using a genetic algorithm[END_REF][START_REF] Stull | Model-based structural health monitoring of naval ship hulls[END_REF][START_REF] Banks | Real time computational algorithms for eddy-current-based damage detection[END_REF]. In this work, we narrow our focus on the problem of elasticity imaging in systems that involve coupled acoustic-structure interaction (ASI).

While an inverse problem of this nature could arise in many scenarios, from oceanic oil discovery to the nondestructive evaluation of marine structures, we are motivated mainly by applications in medical imaging. Here, the modeling of interaction between an acoustic fluid and biomechanical structure is necessary for imaging areas like the heart wall, arteries, and bladder that have direct contact with blood and other bodily fluid. Since elastic properties are known to be an indicator for distinguishing diseased from healthy human tissue [START_REF] Greenleaf | Selected methods for imaging elastic properties in biological tissues[END_REF][START_REF] Bernal | Material property estimation for tubes and arteries using ultrasound radiation force and analysis of propagating modes[END_REF][START_REF] Plewes | Visualization and quantification of breast cancer biomechanical properties with magnetic resonance elastography[END_REF][START_REF] Kingwell | Arterial stiffness and prediction of cardiovascular risk[END_REF], elasticity imaging has become an essential tool in detecting the onset and monitoring of the progression of a number of diseases. Depending on the imaging modality used, either the response of the tissue itself [START_REF] Fatemi | Quantifying elasticity and viscosity from measurement of shear wave speed dispersion[END_REF][START_REF] Muthupillai | Magnetic resonance elastography by direct visualization of propogating acoustic strain waves[END_REF] or the acoustic emission in the surrounding fluid [START_REF] Fatemi | Ultrasound-stimulated vibro-acoustic spectrography[END_REF][START_REF] Fatemi | Vibro-acoustography: An imaging modality based on ultrasoundstimulated acoustic emission[END_REF] under a prescribed excitation can be used as data for the inverse identification of the tissue's elastic parameters.

Despite its importance in medical imaging and other applications, there exist relatively few computational approaches for parameter estimation in systems with ASI. In [START_REF] Brigham | Inverse estimation of viscoelastic material properties for solids immersed in fluids using vibroacoustic techniques[END_REF], an approach was developed to inversely estimate the viscoelastic properties of a submerged solid using pressure data from the steadystate dynamic response of the system. Since solving the inverse problem requires many costly evaluations of a forward ASI finite element solver, this work was extended in [START_REF] Brigham | Inverse viscoelastic material characterization using pod reduced-order modeling in acoustic-structure interaction[END_REF] to use reduced-order modeling with proper orthogonal decomposition for the forward problem to reduce computation time. The estimation of material parameters in submerged, orthotropic elastic cylinders was performed in [START_REF] Rosario | Identification of material properties of orthotropic elastic cylinders immersed in fluid using vibroacoustic techniques[END_REF], where it was shown that a surface velocity measurement was sufficient to recover the orthotropic parameters, while an acoustic pressure measurement contained only enough information to recover two of the three unknown parameters. While the previously mentioned works operated under the assumption of homogeneous materials with geometries that were known a priori, the work in [START_REF] Aguilo | An inverse problem approach for elasticity imaging through vibroacoustics[END_REF] sought pointwise reconstructions by representing the spatial variation of the unknown elastic modulus using Gaussian radial basis functions. Inclusions were properly identified within a submerged solid using both fluid pressure measurements and surface velocity measurements on the solid.

A common theme in the above studies was the minimization of a norm (generally L 2 ) of the error between computed and measured system responses using non-gradient based optimization techniques to obtain parameter estimates. In this work, as the main point of departure, we postulate the inverse problem as the minimization of a modified error in constitutive equation (MECE) functional. The MECE functional is a combination of the error in constitutive equation (ECE), introduced in [START_REF] Ladeveze | Error estimate procedure in the finite element method and applications[END_REF] as a measure of the discrepancy in the constitutive equations that connect kinematically admissible strains and dynamically admissible stresses, and quadratic error terms that incorporate the measurement data. MECE-based approaches for identification initially appeared in the context of model updating from vibrational data [START_REF] Ladevèze | Updating of finite element models using vibration tests[END_REF][START_REF] Barthe | Validation and updating of industrial models based on the constitutive relation error[END_REF]. The underlying principle was to split the equations into a reliable set (containing, for example, the equilibrium equations, initial conditions, and boundary conditions), to be enforced strictly as constraints, and an unreliable set (that included measured data and the unknown constitutive properties), contributing terms in the MECE functional to be minimized. More recent extensions of the MECE approach to time-domain formulations [START_REF] Allix | Identification strategy in the presence of corrupted measurements[END_REF][START_REF] Feissel | Modified constitutive relation error identification strategy for transient dynamics with corrupted data[END_REF] were shown to be very robust in the presence of high levels of noise, while also providing an inherent error estimate through the ECE. The MECE method was extended to large scale identification problems in [START_REF] Banerjee | Large scale parameter estimation problems in frequency-domain elastodynamics using an error in constitutive equation functional[END_REF], where the authors also showed that this approach displayed fast convergence and accuracy as compared to conventional least-squares minimization approaches.

In this work, we extend the MECE approach for the inverse estimation of elastic material parameters to the context of frequency-domain, coupled ASI systems. The MECE functional is modified to include an additional error term for fluid pressure data, allowing for a general formulation where measurements of either the displacement response in the solid or the acoustic emission in the surrounding fluid can be used as data for the inverse problem. In doing so, the governing equations of the solid system, the acoustic fluid system, and the coupling conditions between them are enforced as constraints, while the MECE functional is minimized to obtain an estimate for the unknown material parameters. The weight parameter that multiplies the misfit data term in the MECE functional in essence behaves as a regularization parameter. Two different strategies for selecting this weight parameter are demonstrated and compared in numerical examples: 1) the discrepancy principle of Morozov [START_REF] Isakov | Inverse Problems for Partial Differential Equations[END_REF][START_REF] Colton | Inverse Acoustic and Electromagnetic Scattering Theory[END_REF], and 2) an error-balance approach where the sum of the squares of the ECE and data mismatch terms that enter the definition of the MECE functional is minimized.

The article is organized as follows. The following section formulates both the forward and inverse problems for a coupled ASI system and then details the MECE approach for solving the inverse problem. Section 3 then comments on some of the practical aspects of the MECE algorithm, including the two regularization methods that are demonstrated in this work. The performance of the numerical method is then demonstrated in Section 4 through results for the recovery of elastic parameters given noisy measurements taken in either the solid or the fluid, on 2D and 3D examples. Concluding remarks are given in Section 5.

Formulation

Forward Model

In this section, we present the governing equations for the steady-state acoustic-structure interaction problem. We consider a linear elastic body Ω s immersed in a semi-infinite fluid domain Ω f with a separating fluid-structure interface denoted as Γ fs . In this work, we assume both small strains and deformations in the solid domain as well as negligible flow and small pressure amplitudes in the fluid domain. Furthermore, body forces throughout each subdomain are taken to be negligible. We consider a constant mass density within the solid domain and a constant mass density and bulk modulus in the fluid domain. The remaining quantities can have spatial and/or frequency dependence, which is suppressed in our notation for simplicity.

The linear elastic solid domain Ω s undergoing time-harmonic motion is governed by (a) the balance equations

∇•σ = -ρ s ω 2 u in Ω s , (1a) σ•n s = t on Γ t , (1b)
σ•n s = -pn s + g on Γ fs , ( 1c
)
where σ is the stress tensor, ρ s is the solid mass density, ω is the angular frequency, u is the displacement field, n s is the unit normal vector pointing outward from Ω s , t and Γ t ⊂ ∂Ω s are the specified traction and its support, p is the acoustic pressure; (b) the kinematic compatibility equations

ε[u] = 1 2 ∇u + ∇u T in Ω s , (2a)
u = 0 on Γ u , ( 2b
)
where ε is the linearized strain tensor and Γ u = ∂Ω s \ Γ t is the constrained part of the boundary; (c) the constitutive (linear elastic) equation

σ = C : ε in Ω s , (3)
where C is the heterogeneous fourth-order elasticity tensor. The boundary condition in Equation (1c) prescribes balance between forces arising from the solid and fluid domains along the wet interface Γ fs , with g denoting a specified traction that may be present over Γ fs . The latter traction would arise, for instance, in biomedical imaging problems where the radiation force of ultrasound is used to excite part of the wet interface [START_REF] Hasegawa | A general theory of rayleigh and langevin radiation pressures[END_REF][START_REF] Zhang | Estimation of tissues elasticity with surface wave speed[END_REF].

The governing equations for the acoustic fluid system can be written in the frequency domain as

∇ 2 p + k 2 p = 0 in Ω f , (4a) p = 0 on Γ p , (4b)
∂ n f p = -f (ik + β) p on Ψ R , (4c)
∂ n f p = ρ f ω 2 u•n f on Γ fs (4d)
where k = ω ρ f /K is the wave number with ρ f and K denoting the fluid mass density and bulk modulus, respectively. The vector n f is the unit normal pointing outward from Ω f . Here, we introduce the artificial boundary Ψ R in order to truncate the semi-infinite fluid domain for computational purposes. Equation (4c) is the non-reflecting radiation condition applied to this boundary as a first order approximation to the Sommerfield radiation condition, where f and β are geometry-specific constants. While this simple treatment suffices in this work, we note that more sophisticated techniques like higher order absorbing conditions or perfectly matched layers [START_REF] Berenger | A perfectly matched layer for the absorption of electromagnetic waves[END_REF] could be substituted here. Equation (4d) arises from the continuity in displacements of fluid and solid particles normal to Γ fs . This interface condition along with Equation (1c) provide the coupling between the two sets of governing equations for the solid and fluid systems.

Since the constitutive tensor C is the main unknown of the inverse problem, the constitutive equation (3) will (in remaining consistent with earlier works based on MECE formulations) be included in the MECE functional to be minimized, while equations (1), ( 2) and (4) will be enforced exactly, in weak form. For the latter purpose, let the spaces of trial and test solutions for displacements and acoustic pressures, respectively, be defined as

U = u| u ∈ [H 1 (Ω s )] d , u = 0 on Γ u (5a) P = p| p ∈ H 1 (Ω f ), p = 0 on Γ p ( 5b
)
where d is the spatial dimension. Furthermore, the space of dynamically-admissible stresses in the solid is defined by

S(u) := σ| σ ∈ H div (Ω s ), ∇•σ = -ρ s ω 2 u in Ω s , σ•n s = t on Γ t (6)
With these definitions, the weak formulation for equations (1) and (2a) is

B(σ, u, w) = F s (p, w) + F(w) ∀w ∈ U (7)
having set

B(σ, u, w) := (σ, ε [w]) Ωs -ρ s ω 2 (u, w) Ωs , (8a)
F s (p, w) := (pn s , w) Γ fs , (8b)
F(w) := (t, w) Γt + (g, w) Γ fs (8c)
and where (•, •) Φ denotes the L 2 (Φ) inner product. Specifically, for two second-order tensor fields x, y, one has (x, y) Φ :=

Φ x : ȳ dΦ = Φ x ij ȳij dΦ (9)
where the over-bar ȳ denotes complex conjugation and the repeated indices indicate summation over the components of x and y.

For the set of equations ( 4) governing the fluid system, the corresponding weak formulation is given as

A f (p, v) = F f (u, v) ∀v ∈ P (10)
where

A f (p, v) := (∇p, ∇v) Ω f -k 2 (p, v) Ω f + f (ik + β) (p, v) Ψ R (11a) F f (u, v) := ρ f ω 2 (u•n f , v) Γ fs (11b)

Modified error in constitutive equation approach for inverse ASI problems

The inverse problem associated with the coupled ASI forward problem consists of estimating the spatial distribution of elastic moduli in Ω s that define the constitutive tensor C, given measured acoustic pressures p m ( xf ), xf ∈ Ω m f ⊆ Ω f and/or measured solid displacements u m ( xs ), xs ∈ Ω m s ⊆ Ω s , obtained at one or more frequencies. Additionally, we assume that the properties of the fluid are known.

In this work, the inverse problem is cast as an optimization problem in which the unknown constitutive tensor is estimated by minimizing a MECE functional. The MECE functional combines two types of error terms: 1) an error in constitutive equation (ECE) functional [START_REF] Ladeveze | Error estimate procedure in the finite element method and applications[END_REF] that measures the discrepancy in the constitutive equations that connect kinematically admissible strains and dynamically admissible stresses, and 2) quadratic error terms quantifying the mismatch between measurement data and computed fields for given material parameter values. The MECE functional in the present ASI context is expressed as

Λ(u, σ, p; C) = U (u, σ; C) + κ s 2 u -u m 2 L 2 (Ω m s ) + κ f 2 p -p m 2 L 2 (Ω m f ) (12)
where

U (u, σ; C) := 1 2 Ωs (σ -C : ε[u]) : C -1 : (σ -C : ε[u]) dΩ s ( 13
)
is the ECE functional for linear elastic materials and κ s and κ f are weighting parameters for the solid displacement and acoustic pressure data, respectively, that control the relative importance of these terms in the inverse problem. The ECE term [START_REF] Brigham | Inverse estimation of viscoelastic material properties for solids immersed in fluids using vibroacoustic techniques[END_REF] quantifies the discrepancy between kinematically admissible displacements and dynamically admissible stresses for a given C that bears the important properties

U (u, σ; C) ≥ 0 ∀C (14)
U (u, σ; C) = 0 ⇐⇒ σ = C : ε[u] (15)
Now, the solution of an inverse problem cast using a MECE framework for an ASI system is given by

(u ⋆ , σ ⋆ , p ⋆ , C ⋆ ) = arg min u∈U , σ∈S(u), p∈P, C∈C Λ(u, σ, p; C), ( 16
)
where C is the search space for the unknown constitutive tensor comprised of all fourth-order tensor fields that are symmetric, positive definite, and bounded, while the remaining spaces U , S, P are defined by (5a), (5b) and [START_REF] Bernal | Material property estimation for tubes and arteries using ultrasound radiation force and analysis of propagating modes[END_REF]. On noting that (16) defines a PDE-constrained optimization problem, we define a Lagrangian functional L :

U × S × P × U × P × C → R as L(u, σ, p, w, v; C) := Λ(u, σ, p; C) -Re [B(σ, u, w) -F s (p, w) -F(w)] -Re [A f (p, v) -F f (u, v)] ( 17
)
where B, A f , F s , F f and F were defined in Section 2.1 and Re[•] denotes the real part of a complex number. Note that the test functions w ∈ U and v ∈ P used in the variational forms of the elastic and acoustic systems, respectively, act as Lagrange multipliers in ( 17) [START_REF] Hinze | Optimization with PDE Constraints[END_REF]. The remainder of this section is devoted to the derivation (and solution strategy) of the first-order optimality conditions for the minimization problem (16).

Derivation of the first-order optimality conditions

We now derive the first-order optimality conditions for the MECE inverse problem ( 16) by taking directional derivatives of the Lagrangian ( 17) with respect to u, σ, p, w, v, C and setting each to zero. The unknown fields u, σ, p, C that satisfy these conditions will be those that minimize the MECE functional [START_REF] Fatemi | Vibro-acoustography: An imaging modality based on ultrasoundstimulated acoustic emission[END_REF] while also satisfying the forward ASI problem. We denote, for example, the directional derivative of the Lagrangian functional L with respect to the stress field σ as L ′ σ , which is given as

L ′ σ , σ = d dθ [L(u, σ + θ σ, p, w, v; C)] θ=0 (18)
for all variations σ ∈ H div (Ω s ). By carrying out the calculation in [START_REF] Ladevèze | Updating of finite element models using vibration tests[END_REF] and setting the result equal to zero, we obtain the following expression for the stress field

L ′ σ , σ = Re ( σ, C -1 : σ) Ωs -( σ, ε[u]) Ωs -( σ, ε[w]) Ωs = Re ( σ, C -1 : σ -ε[u] -ε[w]) Ωs ∀ σ ∈ S (19)
σ = C : ε[u + w] (20)
Proceeding in the same fashion, the partial derivatives of the Lagrangian functional with respect to the remaining mechanical fields, the Lagrange multipliers, and the constitutive tensor are given by

L ′ w , ŵ = Re (σ, ε[ ŵ]) Ωs -ω 2 (ρ s u, ŵ) Ωs -(t, ŵ) Γt -(g, ŵ) Γ fs -(pn s , ŵ) Γ fs ( ŵ ∈ U) L ′ v , v = Re -(∇p, ∇v) Ω f + k 2 (p, v) Ω f -f (ik + β)(p, v) Ψ R + ρ f ω 2 (u•n f , v) Γ fs (v ∈ P) L ′ u , û = Re (C : ε[u] -σ, ε[ û]) Ωs + κ s (u -u m , û) Ωm + ρ s ω 2 (w, û) Ωs + ρ f ω 2 ( û•n f , v) Γ fs ( û ∈ U) L ′ p , p = Re -(∇p, ∇v) Ω f + k 2 (p, v) Ω f -f (ik + β)(p, v) Ψ R -(pn s , w) Γ fs + κ f (p -p m , p) Ω m f (p ∈ P) L ′ C , Ĉ = Ĉ , ε[u] ⊗ ε[ ū] -(C -1 : σ) ⊗ (C -1 : σ) Ωs ( Ĉ ∈ C)
By setting these partial derivatives equal to zero, we obtain the following set of coupled variational equations

A s (u, ŵ) + (C : ε[w], ε[ ŵ]) Ωs -F s (p, ŵ) -F( ŵ) = 0 ∀ ŵ ∈ U (21a) A f (p, v) -F f (u, v) = 0 ∀v ∈ P (21b) A s (w, û) -ρ f ω 2 (v, û•n f ) Γ fs + κ s (u -u m , û) Ω m s = 0 ∀ û ∈ U (21c) (w, pn s ) Γ fs + A f (v, p) + κ f (p m -p, p) Ω m f = 0 ∀p ∈ P (21d) Ĉ , ε[u] ⊗ ε[ ū] -(C -1 : σ) ⊗ (C -1 : σ) Ωs = 0 ∀ Ĉ ∈ C (21e)
To simplify, we have used Equation [START_REF] Allix | Identification strategy in the presence of corrupted measurements[END_REF] in Equations (21a)-(21d) and introduced the bilinear forms A f , given by (11a), and A s , defined by

A s (u, w) := B s (C : ε [u] , u, w) = (C : ε [u] , ε [w]) Ωs -ρ s ω 2 (u, w) Ωs . (22)
The set of variational equations ( 21), along with [START_REF] Allix | Identification strategy in the presence of corrupted measurements[END_REF], represent the first-order optimality conditions for the MECE minimization problem [START_REF] Aguilo | An inverse problem approach for elasticity imaging through vibroacoustics[END_REF]. We would like to point that we do not substitute for σ in Equation (21e) as this form of the equation is conveniently used in the block Gauss-Seidel strategy explained in the next section.

Solution of the first-order optimality conditions

The system of equations ( 21) can be subdivided into (a) a set of four linear equations (21a)-(21d), and (b) the non-linear equation (21e). Like in earlier applications of MECE to parameter identification, e.g. [START_REF] Barthe | Validation and updating of industrial models based on the constitutive relation error[END_REF][START_REF] Banerjee | Large scale parameter estimation problems in frequency-domain elastodynamics using an error in constitutive equation functional[END_REF], this subdivision makes it natural to adopt an iterative alternating strategy of block Gauss-Seidel type whereby each solution iteration for [START_REF] Feissel | Modified constitutive relation error identification strategy for transient dynamics with corrupted data[END_REF] consists in (i) solving the set of linear equations (21a)-(21d) for the mechanical fields and Lagrange multipliers, with C kept fixed, and (ii) updating C via equation (21e) with mechanical fields and multipliers fixed. Here, this method will be seen to entail solving a 4 × 4 complex symmetric block linear system followed (when considering isotropic materials) by simple and explicit update formulas for the elastic moduli, similar to those used in [START_REF] Banerjee | Large scale parameter estimation problems in frequency-domain elastodynamics using an error in constitutive equation functional[END_REF].

Step (i): field and multiplier update The finite element method [START_REF] Bathe | Finite Element Procedures[END_REF] is used to discretize and transform the coupled weak formulation (21a)-(21d) into a linear system of equations. Using standard Voigt notation, the displacement fields, pressure field, and corresponding test functions are replaced with the following discrete approximations

u h = [N s ] {u} , ûh = [N s ] {û} , ε[u h ] = [B s ] {u} , w h = [N s ] {w} , ŵh = [N s ] { ŵ} , ε[w h ] = [B s ] {w} , p h = [N f ] {p} , ph = [N f ] {p} , ∇p h = [B f ] {p} , v h = [N f ] {v} , vh = [N f ] {v} , ∇v h = [B f ] {v} ,
where [N s(f ) ] and [B s(f ) ] represent matrices of finite element shape functions and their derivatives with respect to spatial coordinates, respectively, for the solid (fluid) system. The terms, {•}, represent vectors of nodal quantities in the finite element mesh. Inserting these approximations into the variational problem (21a)-(21d) and simplifying, we arrive at the following 4 × 4 block system of equations:

     [H s ] -[S] T [K s ] [0] -ρ f ω 2 [S] [H f ] [0] [0] -κ s [Q s ] [0] [H s ] -ρ f ω 2 [S] T [0] -κ f [Q f ] -[S] [H f ]]               {u} {p} {w} {v}          =          {P s } {0} -κ s {R s } -κ f {R f }          (23)
where the sub-matrices and vectors in the system are defined as

[H s ] = [K s ] -ω 2 [M s ] , [H f ] = [K f ] + iω [C f ] -ω 2 [M f ] [K s ] = elements Ω e s [B s ] T [C] [B s ] dΩ, [K f ] = elements Ω e f [B f ] T [B f ] dΩ [M s ] = elements Ω e s ρ s [N s ] T [N s ] dΩ, [M f ] = elements Ω e f ρ f K [N f ] T [N f ] dΩ [Q s ] = elements Ω e s [N s ] T [N s ] dΩ, [Q f ] = elements Ω e f [N f ] T [N f ] dΩ [S] = elements Γ e fs [N f ] T n T f [N s ] dS, [C f ] = (iω) -1 elements ∂Ψ e R f (ik + β) [N f ] T [N f ]dS {R s } = elements Ω e s [N s ] T u m dΩ s , {R f } = elements Ω e f [N f ] T p m dΩ s {P s } = elements Γ e t [N s ] T t dS + elements Γ e f s [N s ] T g dS
where [C] is the linear elastic constitutive matrix in Voigt notation. Here, the expressions for

[Q s ], [Q f ], {R s },
and {R s } assume that the full fields are available for the measurement data u m and p m . In the case where only sparse data is available, the matrices [Q s ] and [Q f ] are replaced with diagonal Boolean matrices whose non-zero entries correspond to degrees of freedom that have been measured. Likewise, the vectors {R s } and {R s } are modified to contain the measurement values at those degrees of freedom and zeros elsewhere. It is noted that the equations in [START_REF] Isakov | Inverse Problems for Partial Differential Equations[END_REF] can easily be rearranged to yield a complex symmetric matrix, but the system remains indefinite. In this study, we employ the parallel, direct linear solver PARDISO [START_REF] Schenk | Solving unsymmetric sparse systems of linear equations with pardiso[END_REF][START_REF] Schenk | On fast factorization pivoting methods for symmetric indefinite systems[END_REF] which provides efficient solutions to such systems.

Step (ii): constitutive update This step consists in solving equation (21e) for C, with u, w, p, v set to the solution of step (i). As shown in [START_REF] Banerjee | Large scale parameter estimation problems in frequency-domain elastodynamics using an error in constitutive equation functional[END_REF], this reduces to explicit update formulas when considering isotropic linear elastic materials, for which the elasticity tensor can be expressed in terms of the bulk modulus B and the shear modulus G as

C = B - 2 3 (I ⊗ I) + 2GI (24)
(with I and I denoting the second and fourth order identity tensor, respectively). The update formulae for B and G are obtained by first decoupling the stress and strain tensors into deviatoric and volumetric components

σ = σ d + qI, ε[u] = ε d [u] + 1 3 e u I
where σ d and ε d are the deviatoric stress and strain tensors, respectively, q = 1 3 tr(σ) is the mean stress and e u = tr(ε[u]) is the volumetric strain. Then, using [START_REF] Colton | Inverse Acoustic and Electromagnetic Scattering Theory[END_REF]

in (21e), yields Ĝ 2(ε d [u], ε d [u]) - (σ d , σ d ) 2G 2 + B (e u , e u ) - (q, q) B 2 = 0 ∀ B, Ĝ (25)
From which, given u, w at iteration n, we obtain at iteration n + 1

B n+1 = (q, q) 1/2 (e u , e u ) 1/2
= B n (e u + e w , e u + e w ) 1/2 (e u , e u ) 1/2 (26a)

G n+1 = (σ d , σ d ) 1/2 (ε d [u], ε d [u]) 1/2 = G n (ε d [u + w], ε d [u + w]) 1/2 (ε d [u], ε d [u]) 1/2 (26b)
where the second equality results from the dependence of σ on both u and w as in Equation [START_REF] Allix | Identification strategy in the presence of corrupted measurements[END_REF]. The domain of definition for the inner products appearing in ( 26) can be interpreted, depending on the problem at hand, as finite elements, points, or, in general, parts of the body over which material properties remain constant. The latter is useful, for example, in inverse problems involving homogeneous materials with known geometries where one seeks only the elastic parameters defining each region. In summary, the MECE inverse problem for ASI in Eq. ( 16) is solved by first forming the Lagrangian functional in [START_REF] Ladeveze | Error estimate procedure in the finite element method and applications[END_REF]. Then, by taking directional derivatives of the Lagrangian with respect to the unknown mechanical fields, the Lagrange multipliers, and the constitutive tensor, we arrive at the coupled set of equations in ( 20) and ( 21), representing the first-order optimality conditions. A block Gauss-Seidel solution strategy is adopted to solve this set of equations by alternating between the solution of Equations (21a) -(21d) with the current value of C and then Equation (21e) with the updated values of u, p, w, v, and σ. We reiterate that in practice, one iteration of this Gauss-Seidel algorithm amounts to solving the block system in [START_REF] Isakov | Inverse Problems for Partial Differential Equations[END_REF], evaluating the stress field using [START_REF] Allix | Identification strategy in the presence of corrupted measurements[END_REF], and then updating the elastic moduli using [START_REF] Zhang | Estimation of tissues elasticity with surface wave speed[END_REF].

Remark 1

The MECE minimization strategy is easily adaptable to situations where measurement data is acquired at multiple frequencies by defining the MECE functional as

Λ(u 1 , . . . , u N , σ 1 , . . . , σ N , p 1 , . . . , p N , C) = N i=1 Λ i (u i , σ i , p i ; C) ( 27
)
where data is obtained at the frequencies {ω i } N i=1 and Λ i depends on ω i through the data at that frequency. A MECE iteration hence consists of solving equations (21a)-(21d) for each frequency (with Λ = Λ i , yielding fields u i , w i , p i , v i ), followed by a constitutive update where (for the isotropic case) the righthand sides in [START_REF] Zhang | Estimation of tissues elasticity with surface wave speed[END_REF] are summed over the relevant frequencies.

Remark 2

The block Gauss-Seidel strategy described above can be instead interpreted as an alternating minimization method [START_REF] Banerjee | Large scale parameter estimation problems in frequency-domain elastodynamics using an error in constitutive equation functional[END_REF], since steps (i) and (ii) are the stationarity equations for the partial minimization of the MECE functional with moduli C fixed (to the previous iterate) and with fields u, w, p, v fixed (to their latest value), respectively.

Remark 3 Methods based on Newton or quasi-newton approaches [START_REF] Epanomeritakis | A newton-cg method for large-scale three-dimensional elastic full-waveform seismic inversion[END_REF] may also be used to solve the MECE optimization problem [START_REF] Aguilo | An inverse problem approach for elasticity imaging through vibroacoustics[END_REF]. These approaches have received little or no attention in the current literature related to ECE-type methods and present an interesting future direction to investigate.

Weighting parameter selection and regularization in MECE

In this section, we discuss the role and adjustment of the weighting parameters, κ s and κ f , in the MECE functional [START_REF] Fatemi | Vibro-acoustography: An imaging modality based on ultrasoundstimulated acoustic emission[END_REF]. These parameters, which define the balance between minimizing the ECE ( 13) and matching the experimental data, are of great importance in the quality of the reconstruction of C. In that respect, the MECE functional ( 12) is analogous to a regularized cost functional [START_REF] Isakov | Inverse Problems for Partial Differential Equations[END_REF][START_REF] Calvetti | Tikhonov regularization and the l-curve for large discrete ill-posed problems[END_REF], with regularization provided by the ECE term. In this analogy, κ s and κ f act as reciprocals of regularization parameters in the usual sense.

We now introduce some additional notation to facilitate the discussion to follow. First, to ensure consistent units (i.e. energy) and proper scaling among the components of the MECE functional [START_REF] Fatemi | Vibro-acoustography: An imaging modality based on ultrasoundstimulated acoustic emission[END_REF], the coefficients κ s and κ f are recast in the following form, as proposed in [START_REF] Banerjee | Large scale parameter estimation problems in frequency-domain elastodynamics using an error in constitutive equation functional[END_REF]:

κ s = α s A s , κ f = α f A f , with A s := U (C 0 ) u m 2 L 2 (Ω m s ) , A f := U (C 0 ) p m 2 L 2 (Ω m f ) (28)
where U (C 0 ) is the strain energy in the solid system of the solution of the forward problem ( 1)-( 4) with the initial guess C = C 0 and α s , α f are non-dimensional weighting parameters. Moreover, let u α , σ α , p α and C α denote the fields and elastic moduli obtained upon convergence of the minimization problem [START_REF] Aguilo | An inverse problem approach for elasticity imaging through vibroacoustics[END_REF] for given α := {α s α f }. The converged values of the ECE, displacement misfit and pressure misfit components of the MECE functional (12) are, respectively,

U (α) := U (u α , σ α , C α ), D s (α) := 1 2 A s u α -u m 2 L 2 (Ω m s ) , D f (α) := 1 2 A f p α -p m 2 L 2 (Ω m f ) , (29)
(with A s , A f as defined in [START_REF] Hinze | Optimization with PDE Constraints[END_REF]) so that the converged MECE value Λ(α) is given by

Λ(α) = U (α) + α s D s (α) + α f D f (α) (30)
Adjusting the MECE weighting coefficients is now a matter of finding an effective strategy for setting the parameters α. As the parameters in α increase, D s (α) and D f (α) decrease. That is, greater emphasis is put on reducing the L 2 data discrepancies and the measured data is approximated more closely. However, the measurement noise is also emphasized as these parameters increase. Alternatively, decreasing α yields reconstructions that better reduce the ECE, i.e. emphasize the satisfaction of the constitutive relation (3), while eventually losing important information contained in the measured data. Since α ultimately decides the tradeoff between over-smoothing a solution and over-fitting noisy input data, selecting appropriate values is akin to regularization with the MECE algorithm. We present two strategies for the selection of α. The first is based on the discrepancy principle of Morozov [START_REF] Isakov | Inverse Problems for Partial Differential Equations[END_REF][START_REF] Colton | Inverse Acoustic and Electromagnetic Scattering Theory[END_REF], a well-established approach, which assumes the level of noise in the data to be known a priori. In the second approach, applicable to cases with unknown noise, α is selected by minimizing an error-balance function that we propose on the basis of heuristic considerations

Adjustment method 1: Morozov discrepancy principle. This approach, which assumes that the level of noise δ in the measurement data is known, exploits the discrepancy principle of Morozov [START_REF] Isakov | Inverse Problems for Partial Differential Equations[END_REF][START_REF] Colton | Inverse Acoustic and Electromagnetic Scattering Theory[END_REF]. Here, the parameters in α are chosen to be the smallest positive numbers such that the final discrepancy between the computed and measured system response is at the level of the noise. This amounts (considering displacement data for definiteness) to choosing α such that

u α -u m L 2 (Ω m s ) u m L 2 (Ω m s ) = cδ, (31)
where c is a constant (taken as c = 1 in this work). This strategy is implemented by solving the inverse problem for increasing values of α (starting from a very low value) until condition [START_REF] Schenk | On fast factorization pivoting methods for symmetric indefinite systems[END_REF] is met within a chosen tolerance.

Adjustment method 2: error balance. For the situation when the noise level δ on the data is unknown, we now propose an approach for weighting parameter selection that is based on an errorbalance between the terms in the converged MECE functional [START_REF] Schenk | Solving unsymmetric sparse systems of linear equations with pardiso[END_REF]. The underlying, heuristic, idea is to adjust α so as to strike a balance between the minimization of the different components of the MECE functional, hoping to achieve a satisfactory tradeoff between over-smoothing the solution (over-emphasis on ECE minimization) and over-fitting the data (over-emphasis on data discrepancy minimization). In practice, we propose to choose α so as to minimize the following error-balance function J (α):

J (α) := U 2 (α) + D 2 s (α) + D 2 f (α) (32)
Without loss of generality and for the sake of simplicity, we will limit the focus of this work to cases in which either pressure or displacements are measured, but not both. Hence, α will be treated as scalarvalued henceforth. Herein, we adopt a simple approach to approximate that value of α. The optimization problem ( 16) is solved for multiple values of α in a predetermined range, yielding corresponding values for U (α), D s (α), D f (α). The optimal value of J (α) is then obtained by solving a one-dimensional minimization problem. Although this approach seems very computationally expensive, in practice solutions for different values of α are independent and can be obtained simultaneously through parallel solutions of optimization problems ( 16). An adaptive approach for determining more efficiently an optimal value of α is highly desirable and will be pursued in future work.

Numerical results

In this section, we apply the MECE algorithm to estimate unknown material properties in ASI systems for three different examples. We consider nearly-incompressible materials assuming a known bulk modulus equal to that of the surrounding fluid medium and seek to recover an unknown shear modulus field throughout the solid domain. The mean-dilatation approach [START_REF] Hughes | The finite element method: linear static and dynamic finite element anaysis[END_REF] is used to handle near-incompressibility. We examine the case where the full displacement field is measured as well as the case where displacement data is only available in one direction. We also test the performance of the algorithm in recovering elastic parameters given pressure measurements in the fluid domain. In all examples, we add random noise to the inverse problem data to simulate the measurement errors that are inherent in practice.

Denoting the synthetic (computed and interpolated) displacement at a node i as ûi , the corresponding noisy measurement u m i is given as

u m i = ûi (1 + δr i ), ( 33
)
where r i is a normal random variable with zero mean and unit variance and the parameter δ is a prescribed relative noise level. In our examples, we use δ = 0.01, 0.03, and 0.05.

In Example 1, we consider a two-dimensional problem where the shear modulus field is estimated using incomplete and noisy displacement data in the solid domain. Example 2 is another two-dimensional problem in which the shear modulus of three different materials with geometries that are known a priori are recovered using sparse pressure measurements from sensors within the fluid domain. In Example 3, we image a shear modulus field in a three-dimensional acoustic-structure system using a full displacement field.

Remark 4

The frequencies in the examples presented herein were chosen so as to maintain a low shear wavenumber in the solid. This choice was made without loss of generality and in the interest of avoiding excessively fine meshes. Notice that because of the large differences in the bulk and shear moduli, the wavelengths in the fluid are very large as compared to the solid domain. This fact is not a problem when pressure is used to identify shear moduli as long as the pressure field remains sensitive to changes in the material.

Example 1: 2D imaging with displacement data

We consider a two-dimensional acoustic-structure system in which a square solid domain under plane strain conditions is submerged in an infinite fluid medium. We seek to use the noisy and incomplete (i.e. one-directional) displacement response of the solid under an applied traction to recover the unknown shear modulus field. We study the quality of the reconstruction for varying levels of data noise and illustrate the performance of both the discrepancy principle and the error-balance technique for choosing the weight α s .

A sketch of the problem domain for this example can be seen in Figure 1(a). The fluid medium has been truncated using a circle of radius R = 0.1 m on which the non-reflecting radiation condition (4c) has been specified. The solid domain is a square with sides of length L = 0.1 m containing a centered, circular inclusion with radius r = 0.02 m. A uniform traction g = -ge y , with g = 10 3 N/m, was applied in the ydirection along the top surface of the solid. As was mentioned in Section 2.1, this type of traction may be generated, for instance, through the radiation force of ultrasound. The coupling conditions (1c) and (4d) are present along the solid-fluid interface. The mass density and bulk modulus are assumed to be constant and the same for both the fluid and solid domains: ρ f = ρ s = 1000 kg/m 3 and K = B = 2.2 × 10 9 Pa. We consider an inclusion shear modulus (G i ) that is four times stiffer than that of the background material (G b ): G i = 4G b = 6.0 × 10 4 Pa.

! ! !! x y L r R A B Ω s Ω f G i G b (a) ! ! ! ! ! ! ! ! ! ! ! ! ! !!!!!!! ! ! ! ! ! ! !! R Ω s Ω f Sensors r 3 r 2 G 3 G 2 G 1 (b)
The displacement data used for the inverse problem is generated artificially by solving the coupled ASI problem (Equations ( 1) to ( 4)) with the true shear moduli on an unstructured finite element mesh with 70, 000 nodes. For the solution of the inverse problem with MECE, however, a coarser mesh with 12, 000 nodes is used, where the solid domain is a uniform 58 × 58 grid (i.e. the geometry of the inclusion is not explicity meshed). Both meshes are made of four-node, bilinear elements. The displacement within the solid domain is interpolated from the finer mesh to the coarser mesh, which avoids committing the "inverse crime" and deliberately introduces a modeling error. Only those displacements in the y-direction are then kept as data to recover the shear modulus. The data is collected for two frequencies (f = 40.0, 50.0 Hz), and thus the whole input data for the problem consists of 6, 738 displacement measurements.

Remark 5

The reconstruction of elastic moduli using full or one-directional interior data, such as the ones presented in the examples herein, is relevant to many important practical areas. For instance, reconstruction of moduli from interior data is the main focus of the active research area of elasticity imaging (see e.g. [START_REF] Sinkus | Imaging anisotropic and viscous properties of breast tissue by magnetic resonance-elastography[END_REF][START_REF] Kwon | Shear modulus decomposition algorithm in magnetic resonance elastography[END_REF] and references therein).

Shear reconstructions using the discrepancy principle

The inverse problem ( 16) is first solved using the discrepancy principle of Morozov, i.e. selecting α s such that equation ( 31) is approximately satisfied for c = 1 when the MECE algorithm has converged (note that for this example α f = 0). We note that in this study, the algorithm is deemed to have converged when the relative change in the functional [START_REF] Fatemi | Vibro-acoustography: An imaging modality based on ultrasoundstimulated acoustic emission[END_REF] between two successive iterations is below 1%. The latter criterion was confirmed to be adequate for the examples presented herein as it was verified (from extensive numerical testing) that the relative change in mechanical fields and moduli was negligible for this level of change in the functional. For the following results, an initial guess of G 0 = G b = 1.5 × 10 4 Pa is used for the shear modulus throughout the entire solid domain. The effect of measurement error on the recovered shear modulus is studied by considering solutions to the inverse problem in the presence of three different levels of noise (1%, 3%, and 5%).

The selected value of α s for each noise level considered is given in Table 1, together with the final value of the relative discrepancy between the measured and computed displacements (left hand side of (31)) when the algorithm converged and the number of MECE iterations until convergence. As expected, the selected value of α s increases as the noise level δ decreases. This is due to the fact that higher weight (larger α s ) is needed on the discrepancy term of the MECE functional [START_REF] Fatemi | Vibro-acoustography: An imaging modality based on ultrasoundstimulated acoustic emission[END_REF] for the minimization to reach lower thresholds in the discrepancy equation [START_REF] Schenk | On fast factorization pivoting methods for symmetric indefinite systems[END_REF]. Note also that the higher value of α s for the case of 1% caused faster convergence (5 iterations versus 21 for 3% and 5% noise), an observation that is consistent with earlier work [START_REF] Banerjee | Large scale parameter estimation problems in frequency-domain elastodynamics using an error in constitutive equation functional[END_REF].

The shear modulus field estimates for each noise level, with α s obtained using the discrepancy principle, are displayed in Figure 2. Figure 2(a) shows the solutions for each noise level plotted across line AB of Figure 1(a) compared with the reference solution. Figures 2(b,c,d) show the two-dimensional recovered shear modulus field for 1%, 3%, and 5% noise levels, respectively. It can be seen that the level of accuracy for the case of 3% and 5% noise is comparable, while the solution for 1% noise is more accurate, as expected. Indeed, Figure 2 Table 1: The value of MECE weight α s selected for each noise level using the discrepancy principle [START_REF] Schenk | On fast factorization pivoting methods for symmetric indefinite systems[END_REF] along with the final value of the relative discrepancy (l.h.s. of ( 31)) and the number of iterations it took for the MECE algorithm to converge slightly overestimated for the higher noise levels. Additionally, the two-dimensional plots show that the boundary of the inclusion is more distinctly identified for 1% noise than for the higher noise levels. Generally speaking, however, the solutions for each case provide satisfactory estimates of the background and inclusion shear modulus and accurately resolve the boundary between the materials.

Shear reconstructions using error balance

We now demonstrate the use of the error balance approach for selecting α s . For comparison purposes, the problem setup described in Section 4.1 remains unchanged. To carry out the approach, the inverse problem ( 16) is solved for weighting coefficients sampling the range 0.1 ≤ α s ≤ 100, and postprocessing is done to find the value of α s that minimizes the error balance function J (α s ) defined by [START_REF] Epanomeritakis | A newton-cg method for large-scale three-dimensional elastic full-waveform seismic inversion[END_REF]. The resulting plot of J versus α s , displayed in Figure 3 for the different noise levels, shows that the method selects α s to be 5, 10, and 10 for 1%, 3%, and 5% noise, respectively. The corresponding shear field reconstructions for each noise level using the error balance approach are shown in Figure 4. As in the previous section, the solutions are compared with the reference along line AB in Figure 4(a), while Figures 4(b,c,d) show the two-dimensional estimated shear modulus fields for 1%, 3%, and 5% noise levels, respectively. Again, the relative accuracy of the recovered fields are satisfactory for each noise level in terms of the estimated magnitudes for each material as well as the resolution of the boundary between them. The line plots in Figure 4(a) show that the higher noise levels tend to overestimate the inclusion shear modulus more so than the case of 1% noise. However, it can be seen in this plot, as well as in Figures 4(b,c,d), that the inclusion geometry is more accurately identified for 3% and 5% noise than for 1% noise. The latter observation seems to be a reflection of the error-balance approach producing a weight parameter value that is too small for the lowest noise level case, further emphasizing the non-optimal nature of the approach. The relative performance of using the discrepancy principle versus the error balance approach for selecting the weighting coefficient is now discussed. To facilitate a comparison of the accuracy between the two methods, the solution error is computed for a wide range of α s values for each noise level. The error metric used is a relative ℓ 2 norm over the elements in the solid domain, given as Figure 5 shows ÊG (α s ) plotted against α s for each noise level. In addition, the values for α s selected using the discrepancy principle are indicated by circle markers, while those selected using the error balance approach are marked by triangles. For each noise level, it can be seen that the discrepancy principle produces a value of α s that is closer to the optimal value α ⋆ s (such that ÊG (α ⋆ s ) is the lowest possible solution error for a given data set) than the error balance approach . Indeed, using the discrepancy principle leads to values of α s that are very close to optimal for 3% and 5% noise, and although it is farther off for 1% noise, the resulting error ÊG is still close to its optimal value ÊG (α ⋆ s ). On the other hand, the error balance method yields choices of α s that are also near optimal for 3% and 5% noise but significantly far from the optimal value for 1% noise. Furthermore, the error balance selection of α s appears to be relatively insensitive to the noise level, which is a drawback of the approach. The results from this example suggest that using the discrepancy principle to select the MECE weighting coefficient outperforms the error balance approach. The solutions for all noise levels using each method appear to be satisfactory in Figures 2 and4, but the rigorous error comparison in Figure 5 shows that, in terms of accuracy, the discrepancy principle outperforms the error-balance approach. Perhaps more significant than the difference in accuracy is the trending behavior of each method for varying noise levels. The discrepancy principle yields values of α s that increase consistently with decreasing noise, following the same expected trend as the optimal value α ⋆ s . The error balance approach shows little variance for different noise levels and follows the opposite trend, selecting a slightly smaller weighting coefficient for 1% noise. Of course, employing the discrepancy principle in practice requires a priori knowledge of the noise level, while the error balance approach does not. Hence, we reiterate that the main appeal of the latter method is its generality and argue that in the absence of information about noise level, it represents a viable approach for selecting the weighting parameter.

ÊG (α s ) = NE i=1 G (i) true -G (i) est (α s ) 2 NE i=1 G (i) true 2 (34)

Example 2: 2D modulus estimation with pressure data

We again consider a two-dimensional acoustic-structure system that consists of a square solid domain under the plane strain assumption submerged in an infinite fluid medium, see Figure 1(b). In this case, however, we use sparse and noisy measurements of the acoustic pressure in the fluid to estimate the elastic properties of the solid. The latter is made of two inclusions of different size and material embedded in the background matrix and thus involves three homogeneous materials whose spatial distribution is assumed to be known. Taking the center of the square as the origin, material two is centered at (-0.025 m,-0.025 m) with radius r 2 = 0.015 m and material three is centered at (0.02 m,0.02 m) with radius r 3 = 0.02 m (Fig. 1(b)).

The boundary conditions and dimensions of the solid and fluid domains, as well as the material properties of the fluid and the bulk modulus of the solid, are as in Example 1, while the shear modulus values are: G 1 = 1.5 × 10 4 Pa, G 2 = 3.0 × 10 4 Pa, and G 3 = 6.0 × 10 4 Pa. The acoustic pressure is measured at thirty-two sensors spaced evenly on a ring in the fluid domain. The system is excited at a The inverse problem reduces to estimating the values G 1 , G 2 , G 3 of the shear moduli in the background and inclusions. The relevant MECE weight α f is selected using both the discrepancy principle and the error balance approach. The pressure data for the inverse problem is generated by solving the forward problem with true moduli values, storing the pressure values at sensor locations, and adding artificial noise according to [START_REF] Calvetti | Tikhonov regularization and the l-curve for large discrete ill-posed problems[END_REF] with u m replaced by p m . The forward problem is solved on a finite element mesh with 54, 000 nodes while a coarser mesh with 17, 000 nodes is used for the inverse problem. In this case, the inclusions are meshed in the grid used for the inversion since a known spatial distribution of moduli is assumed.

The unknown moduli G = [G 1 G 2 G 3 ] T are estimated using a uniform initial guess G 0 = 1.5 × 10 4 Pa with pressure data polluted with 1%, 3%, and 5% noise. The relative solution accuracy is defined as

ĒG (α f ) = G true -G est (α f ) G true , (35)
The value of ĒG (α f ) is plotted versus α f for each noise level in Figure 6, with the values of α f selected using the discrepancy principle (1.5, 0.2, and 0.1 for 1%, 3%, and 5% noise, respectively) and the error balance approach (2.5 for each noise level), respectively, indicated by circles and triangles. Selecting α f using the discrepancy principle produces significantly more accurate solutions for 3% and 5% noise levels while the error balance solution is slightly more accurate for 1% noise. As in Example 1, α f increases as the noise level δ decreases when selected using the discrepancy principle, but is insensitive to δ when selected using error balance. Reconstruction errors ĒG (α f ) increase with δ for both weight adjustment methods, as expected. The overall identification accuracy is nevertheless satisfactory for all considered cases, with all relative errors on individual moduli found to be below 10% and most of them well below 5%. We mention in passing that we also studied the effect of different initial guesses of shear moduli on the accuracy of the resulting reconstructions for the current example. We found that the behavior of the proposed algorithm was insensitive to the initial guess (at least for the current example). Furthermore, the results reported above for the error balance and discrepancy principle approaches remained the same. We note that the pointwise reconstruction of the shear modulus field without assuming known geometry a priori may not be feasible with the current problem setup. This is because having pressure data is akin to having only surface displacement measurements, and therefore less information is available about the unknown shear moduli than when interior solid data is available. A pointwise reconstruction of the shear field might be feasible using pressure data collected for a range of frequencies and different loading configurations; e.g. when using approaches similar to those found in seismic imaging. However, this is beyond the scope of this work.

Example 3: 3D imaging with displacement data

The MECE algorithm is now applied to a larger scale, three-dimensional, inverse problem. Data consisting of the full displacement field in all directions is used to recover an unknown shear modulus in a solid cube The displacement data used for the inverse problem is generated artificially by solving the coupled ASI problem (1)-( 4) with the true shear moduli on a mesh with 43, 000. A coarser mesh with 27, 000 nodes (13, 000 nodes in the solid) is then used for the inverse problem. Quadratic (10-node) tetahedral elements were used in both meshes. The inverse problem was solved for noise levels of 1%, 3%, and 5% with an initial guess equal to the background shear modulus. In this example, only the discrepancy principle was used to select α s , which yielded values of 50.0, 3.0, and 1.25 for 1%, 3%, and 5% noise, respectively.

Figure 7 shows the reconstructed shear modulus on its horizontal plane of symmetry, for the case of 3% noise. The inclusion is easily distinguished from the background matrix, with the estimated values of the shear modulus in each region in good agreement with their true couinterparts (G i = 6.0 × well with the true inclusion, whose meshed shape is shown for comparison. Finally, a comparison of the estimated shear modulus distribution for different noise levels is shown in Figure 9. The solutions are plotted across the x-axis through the center of the solid domain and compared with the reference solution. As expected, the reconstruction accuracy somewhat declines, and the discontinuity between the two materials becomes in particular less well resolved, as the noise level is increased. The solution for each noise level nevertheless clearly identifies the embedded inclusion and provides satisfactory estimates of the shear moduli. Moreover, a low number of iterations of the MECE algorithm was found sufficient to achieve convergence to the solution of the optimality system for all values of α s . For instance, the solutions for 1%, 3%, and 5% noise were obtained in 11, 19, and 17 iterations, respectively.

Conclusions

In this study, we proposed a framework for the inverse identification of material properties in frequencydomain, coupled acoustic-structure interaction (ASI) systems using a modified error in constitutive equations (MECE) approach. Our formulation allows for both measurements of displacements in the solid and the acoustic pressure in the neighboring fluid to estimate the linear elastic parameters that define the solid. Furthermore, we demonstrated two different methods to select the MECE weighting coefficient: (i) the discrepancy principle of Morozov and (ii) an error-balance approach. The latter, which did not produce optimal results and was outperformed by the former from an accuracy standpoint, has the advantage of remaining effective (for practical purposes) even when the level of noise on the input data is unknown a priori. Overall, the numerical results in this paper showed the effectiveness of each approach in generating satisfactory reconstructions in the presence of noisy and incomplete data.

Figure 1 :

 1 Figure 1: Diagrams of the problem domains in (a) Example 1 and (b) Example 2.

  (a) shows that the true magnitude of the inclusion shear modulus is Noise

Figure 2 :

 2 Figure 2: Shear reconstruction, MECE weight α s selected using the discrepancy principle. (a) The solution for each noise level, plotted along line segment AB of Fig. 1(a). (b) Recovered shear field for 1% noise. (c) 3% noise. (d) 5% noise. Units: Pa, m

Figure 3 :Figure 4 :

 34 Figure 3: The error-balance function J (Eq. (32)) versus α s , with minimum values denoted with an "×".

Figure 5 :

 5 Figure 5: The shear reconstruction error (34) versus α s for different noise levels in Example 1. The circle and triangle markers denote the values of α s selected using the discrepancy principle and the error balance technique, respectively

Figure 6 :

 6 Figure 6: The relative error (35) in the estimated shear moduli versus α f for different noise levels in Example 2. The circle and triangle markers denote the values of α f selected using the discrepancy principle and the error balance technique, respectively.

Figure 7 :

 7 Figure 7: A clip plane contour plot of the recovered shear modulus in Example 3 for 3% noisy data. Units: Pa

10 4 Figure 9 :

 49 Figure 8: A threshold plot of the recovered shear modulus in Example 3 for 3% noisy data. The mesh outline of the true inclusion is shown for comparison. Units: Pa.

Acknowledgments

This work was partially supported by NIH Grants #EB002640 and #EB002167.