A Size-Dependent Cost Function to Solve the Inverse Elasticity Problem

Abstract

Characterizing nonhomogeneous elastic property distribution of solids is of great significance in various engineering fields. In this paper, we observe that the solution to the inverse problem utilizing the standard optimization-based inverse approach is sensitive to the sizes of inclusions. The standard optimization-based inverse approach minimizes a cost function, containing the absolute error between the measured and computed displacements in L2 norm. To address this issue, we propose a novel inverse scheme to characterize nonhomogeneous shear modulus distribution of solids. In this novel method, the cost function is modified, and is dependent on the size of the inclusions. A number of simulated experiments are performed, and demonstrate that the proposed approach is capable of improving the shear modulus contrast in inclusions and reducing the size sensitivity. Furthermore, a theoretical analysis is conducted to validate what we have observed in simulated experiments. This theoretical analysis reveals that what we have observed in the simulated experiments is not induced by the numerical issues Instead, the size sensitivity issue is induced by regularization. The findings of this work encourage us to propose new cost functions for the optimization-based inverse approach to improve the quality of the shear modulus reconstruction.

1. Introduction

The advancement of high-resolution imaging modalities like digital imaging correlation (DIC) makes it possible for us to measure full-field displacement fields with high accuracy [1,2,3]. The measured full-fields displacement fields can be utilized to identify the nonhomogeneous elastic distribution of solids. Mapping the nonhomogeneous elastic distribution usually requires solving an inverse problem in elasticity. There are many inverse approaches that have been developed to solve the inverse problem, including virtual fields methods [4,5,6], optimization-based approaches [7,8], etc. Among them, optimization-based inverse approaches have been commonly used, as they are rigorous and highly robust. In the optimization-based inverse approach, the inverse problem in elasticity is posed to be a constrained minimization problem, where a cost function accounting for the discrepancy between measured and computed displacements is usually defined in the L2 norm and minimized iteratively. Since noise in the measured data will affect the quality of the final reconstructed results, regularization is introduced to prevent over-fitting. The regularization term is usually added as an additional term in the cost function (see Equation (3)). By setting an optimal weight for the regularization term, the inverse problem is solved without fully minimizing the displacement discrepancy term.

However, recent studies have revealed that the cost function that minimizes the absolute error in the displacement discrepancy term does not perform well in some practical cases. In [9], the authors found that the solution to the inverse problem utilizing this standard cost function is boundary-sensitive. That is, the reconstructed elastic property distribution of solids is sensitive to the prescribed boundary conditions (see Figure 1).

To address this issue, a modified cost function was proposed and significantly improved the reconstruction results—in particular, the stiffness contrast between different heterogeneous regions (reconstructed shear moduli in inclusions in Figure 1). The modified cost function was spatially weighted by the strain fields.

In this paper, we will investigate another limitation of the standard optimization-based inverse approach. That is, the solution to the inverse problem depends on the size of the inclusions (see Figure 2). To address this issue, a novel cost function that is size-dependent will be proposed. The solid studied in this paper is assumed to be incompressible linear elasticity. We will adopt a number of simulations and a simple theoretical analysis to show why the proposed approach performs better than the standard approach.

The paper is organized as follows. We will first discuss the standard and proposed inverse algorithms in Methods section. We will then solve the inverse problem of elasticity by both methods with the simulated displacement fields acquired by finite element simulations, in the Results section. To understand the performance of the standard and proposed inverse approaches, we will employ a one-dimensional (1-D) model composed of two elastic bars and perform a simple theoretical analysis. Afterwards, we will thoroughly elaborate the reconstruction results in the Discussion section.

2. Methods

2.1. Forward Problem

The forward problem in linear elasticity is defined as follows:

$$\begin{array}{ll}\mathrm{div}\left(\mathsf{\sigma }\right)=0& \mathrm{on} \Omega \\ \mathbf{u}=\overline{\mathbf{u}}& \mathrm{on} {\Gamma }_u\\ \mathsf{\sigma }\cdot \mathbf{n}=\overline{\mathbf{t}}& \mathrm{on} {\Gamma }_{\sigma }\end{array}$$

(1)

where div is the divergence operator; $\mathsf{\sigma }$ and $\mathbf{u}$ are the stress tensor and displacement vector, respectively; $\overline{\mathbf{u}}$ and $\overline{\mathbf{t}}$ are the displacement and traction values on the displacement and traction boundaries ${\Gamma }_u$ and ${\Gamma }_{\sigma }$, respectively; and ${\Gamma }_u\cup {\Gamma }_{\sigma }$ constitutes the entire boundary of the domain $\Omega$. In addition, n is the outward normal vector. The incompressible linear elastic assumption is made in this paper, and the associated stress–strain relation in two-dimensional cases is written as follows:

$$\mathsf{\sigma }=-p\mathbf{I}+2\mu \mathsf{\epsilon }$$

(2)

where $\mu$ is the shear modulus; p is the hydrostatic pressure; and $\mathsf{\epsilon }$ and I are the strain and identity tensors, respectively. To solve Equations (1) and (2), we adopt the finite element approach in this paper.

2.2. Inverse Problem

The inverse problem is solved by an optimization-based algorithm [7], where a cost function is minimized. Typically, a L2 norm correlation term, that is, the absolute error of the measured and computed displacement fields, is used to define the cost function:

$$F=\frac{1}{2}\parallel {\mathbf{u}}^c(\mu )-{\mathbf{u}}^{meas}{\parallel }_0^2+\frac{1}{2}\alpha \sqrt{{\mu }^2+c^2}$$

(3)

where nodally-defined ${\mathbf{u}}^{meas}$ and ${\mathbf{u}}^c$ are the measured and computed displacement fields, respectively. The computed displacement field is acquired by solving the forward problem using the finite element method at the current estimated shear modulus distribution of the domain of interest. The second term n Equation (3) is the total variation diminishing (TVD) regularization term to regularize the inverse problem. In addition, $\alpha$ is the regularization factor determining the significance of regularization term and chosen by smoothness criteria [10], while $c={10}^{-2}$ is a small constant that resolves the singularity issue when computing the gradients of the cost function with respect to nodal shear moduli. The inverse problem is solved by a quasi-Newton approach referred to as the L-BFGS method [11,12]. The gradients of the cost function with respect to the nodal shear moduli are evaluated efficiently by the adjoint method [7].

In this paper, we propose a novel approach to address the issue that the solution to the inverse problem in elasticity is dependent on the size of the inclusions (see Figure 3b,e). The approach is divided into the following steps. First, we adopt the inverse approach with the standard cost function (Equation (3)) to detect the location of different heterogeneous regions. Second, we delimitate and compute the area of each inclusion from the reconstructed shear modulus distribution. Finally, we solve the inverse problem again by optimization-based inverse algorithms, with the following size-dependent cost function:

$$\pi =\frac{1}{2}\parallel \frac{\left({\mathbf{u}}^c(\mu )-{\mathbf{u}}^{meas}\right)}{{\overline{A}}^i}{\parallel }_0^2+\frac{1}{2}\alpha \sqrt{{\mu }^2+c^2}$$

(4)

${\overline{A}}^i$ is a weight function depending on the normalized area of the i-th inclusion. That is, for the nodes within the i-th inclusion, ${\overline{A}}^i=A^i/A_{\min }$ is the area ratio between the i-th inclusion and the smallest inclusion (the areas of the i-th inclusion and the smallest inclusion are denoted by $A^i$ and $A_{\min }$, respectively). For the background of the domain of interest, $A^i=1$. In the next section, we will employ the proposed approach to solve the inverse problem.

3. Results

3.1. Modulus Reconstruction Acquired by Simulated Data

We test the standard and newly proposed inverse algorithms by the simulated data. The displacement fields are acquired by solving forward problems with a known shear modulus distribution (see Figure 2a), assuming that the solid is in the state of the incompressible plane strain. The forward problem is solved by the stabilized finite element method [13]. In the target shear modulus distribution, there are two inclusions embedded in the soft background. The stiffness ratio between the two inclusions and the background is 5:1. The size ratio between the right and left inclusions varies from 1 to 2. We applied 1% uniform compression on the top edge and restricted the vertical motion of the bottom edge. To prevent rigid body motion, we fix the center point of the bottom edge. The problem domain is discretized uniformly by 14,400 bilinear elements (120 elements along the x and y directions, respectively). To mimic the experimental data, we add 3% noise into the simulated displacement data. In addition, we utilize the vertical displacement component as the measured data in solving the inverse problem, since the displacements perpendicular to the direction of the compressive loading are usually measured with low precision in practice [14].

When these two inclusions have the same size (see Figure 2b), we observe that these two inclusions could be detected clearly from the background utilizing the standard inverse scheme. However, the shear modulus values of these two inclusions are roughly half of the target. This is probably due to the fact that the sizes of the inclusions are too small compared to the overall dimension of the problem domain. Thus, the shear modulus reconstruction is highly sensitive to the noise. We also find that the recovered shear modulus of the left is larger than that of the right inclusion, due to the noise effect.

When we increase the radius of the right inclusion (see Figure 3), the reconstructed results utilizing the standard approach demonstrate that the left inclusion is remarkably underestimated (see Figure 3b,e). Moreover, the underestimation became more significant when the size ratio between the two inclusions increases. More specifically, the recovered shear modulus contrast between the right and left inclusions increases from 1.3 to 1.4. The reason for this will be explained in the next section. When we utilize the proposed inverse approach (see Figure 3c,f), it is observed that the stiffness values of the left inclusion improve significantly. That is, the recovered shear modulus contrasts between the right and left inclusions for Figure 3c,f are 0.99 and 0.95, respectively. However, we also observe that the shape of left inclusion is slightly distorted compared to the shear modulus reconstruction by the standard approach. This is because the weight of the nodes in the inclusion is inversely proportional to its area. Thus, the smaller inclusion gain more weight in the new cost function (Equation (4)), leading to the full minimization of displacements of the left inclusion in the new method.

3.2. Modulus Reconstruction by a One-Dimensional Coupled Model

To represent enormous cases with varying sizes and shear moduli of these two inclusions, we employ an 1-D coupled model for analysis, shown in Figure 4. In Figure 4, the black regions in the two elastic bars represent two inclusions in the two-dimensional (2-D) cases presented in Section 3.1, and the rest denotes the background. The target shear modulus value of the background is fixed to be 1. The target shear modulus of the left and right inclusions is $\overline{\mu }$. In this case, the cost function in the standard approach can be simplified to

$$F={\displaystyle \sum _{i=1}^2\left(\frac{1}{2}{\displaystyle {\int }_0^L{\left(u_i^c-u_i^{meas}\right)}^2\mathrm{d}}y+\alpha \left|{\mu }_{in}^i-1\right|\right)}$$

(5)

where i = 1 denotes the variables for the left and right bars, and ${\mu }_{in}^i$ denotes the recovered shear modulus value in the inclusion of the i-th bar. The recovered shear modulus value of the background is assumed to be reconstructed perfectly, thus being equal to 1.

To establish the relationship between ${\mu }_{in}^1$ and ${\mu }_{in}^2$, we take the gradient of the cost function F with respect to ${\mu }_{in}^i$ and set it to zero. Accordingly, we reach the following relation:

$$\frac{\left({\mu }^1-\overline{\mu }\right){\left(L-a^1\right)}^2{a}^1}{{\left(L{\mu }^1+a^1-{\mu }^1{a}^1\right)}^3\left(L\overline{\mu }+a^1-a^1\overline{\mu }\right)}=\frac{\left({\mu }^2-\overline{\mu }\right){\left(L-a^2\right)}^2{a}^2}{{\left(L{\mu }^2+a^2-{\mu }^2{a}^2\right)}^3\left(L\overline{\mu }+a^2-a^2\overline{\mu }\right)}$$

(6)

where L and aⁱ are the total length of each elastic bar and the size of the black region of the i-th bar. Since we have elaborated on a similar derivation in [15,16], we will not discuss herein.

For the proposed approach, the cost function is rewritten as

$$F={\displaystyle \sum _{i=1}^2\left(\frac{1}{2}\left({\displaystyle {\int }_0^{\frac{\left(L-a^i\right)}{2}}{\left(u_i^c-u_i^{meas}\right)}^2\mathrm{d}}y+{\displaystyle {\int }_{\frac{\left(L-a^i\right)}{2}}^{\frac{\left(L+a^i\right)}{2}}{\left(\frac{u_i^c-u_i^{meas}}{a^i}\right)}^2\mathrm{d}}y+{\displaystyle {\int }_{\frac{\left(L+a^i\right)}{2}}^L{\left(u_i^c-u_i^{meas}\right)}^2\mathrm{d}}y\right)+\alpha \left|{\mu }_{in}^i-1\right|\right)}$$

(7)

In this case, the relationship between the recovered shear moduli of these two inclusions is

$$\frac{\left({\mu }^1-\overline{\mu }\right){\left(L-a^1\right)}^2\left(1+\left(1-a^1\right)a^1\right)a^1}{{\left(L{\mu }^1+a^1-{\mu }^1{a}^1\right)}^3\left(L\overline{\mu }+a^1-a^1\overline{\mu }\right)}=\frac{\left({\mu }^2-\overline{\mu }\right){\left(L-a^2\right)}^2\left(1+\left(1-a^2\right)a^2\right)a^2}{{\left(L{\mu }^2+a^2-{\mu }^2{a}^2\right)}^3\left(L\overline{\mu }+a^2-a^2\overline{\mu }\right)}$$

(8)

If the recovered shear modulus ${\mu }^2$ is determined, the recovered shear modulus of the inclusion of the left bar ${\mu }^1$ can be calculated via Equation (6) or Equation (8). Thus, we can quantitatively compare these two methods by the significance of the underestimation of shear modulus ${\mu }^1$. More specifically, we assume that the recovered shear modulus of the right bar ${\mu }^2$ is underestimated by ε. To this end, ${\mu }^2=\left(1-\epsilon \right)\overline{\mu }$, and ${\mu }^1$ should be $\left(1-\epsilon \right)\overline{\mu }$, as the target shear moduli of these two inclusions are the same. We define a relative error as $error=\left({\mu }^1-\left(1-\epsilon \right)\overline{\mu }\right)/\left(\left(1-\epsilon \right)\overline{\mu }\right)\times 100\%$, and compute the relative error with varying sizes of the right inclusion $a^2$, as shown in Figure 5. In Figure 5, the length of these two bars L, the size of the left inclusion $a^1$, and the target shear modulus $\overline{\mu }$ are fixed to be 1, 0.1, and 5, respectively.

Figure 5 demonstrates that when the sizes of these two inclusions are the same, then ${\mu }^1={\mu }^2$, indicating no error in the reconstruction. However, increasing the size of the right inclusion causes underestimation in the shear moduli of the left inclusion, as the relative error is negative in both methods. The underestimation becomes more serious with the increase of the size of the right inclusion. Importantly, employing the proposed approach is capable of significantly reducing the relative error. This result is consistent with what we have observed in the 2-D simulated experiments presented in Section 3.1.

4. Discussion

This paper presents a novel inverse method to characterize the shear modulus distribution of solids. This method modifies the standard cost function that minimizes the absolute error between the measured and computed displacements, and it is capable of improving the shear modulus contrast of inclusions with varying sizes. We have performed a number of 2-D numerical simulations to test the feasibility of this approach, and compared the performance of the standard and proposed approaches. The reconstruction results revealed that the proposed approach performs much better than the standard method, and improves the shear modulus contrast in the inclusions significantly. Additionally, a 1-D theoretical analysis was conducted to validate what we have observed in 2-D simulations.

The simulations presented in Section 3.1 revealed that the solution to the inverse problem utilizing the standard approach is size-dependent. More specifically, the smaller inclusion was underestimated, and the underestimation became more significant with the reduction of the size of the inclusion. To address this size-sensitivity issue, a modified cost function depending on the size of the inclusion was proposed in this paper. The new approach improved the shear modulus contrast in inclusions remarkably. In fact, the limitation of the standard L2 norm-based optimization inverse approaches was also discussed in [16]. The authors in [16] observed that the solution to the inverse problem by minimizing the standard L2 norm cost function is boundary-sensitive (see Figure 1). Therefore, minimizing L2 norm cost function might not be an optimal way to solve the inverse problem in practical engineering to characterize mechanical property distribution of solids.

We also employed an 1-D theoretical analysis to compare the new and standard approaches with a wide range of size ratios of the two inclusions (see Figure 5). The main merit of this 1-D analysis is that a large number of cases can be tested within seconds. This analysis showed that the new approach reduced the relative error in the shear modulus contrast in the inclusions with different sizes. We should note that we utilized the exact simulated displacement data to evaluate Equations (6) and (8) without introducing noise, since the underestimation of the shear modulus contrast in inclusions using the standard inverse approach is induced by the regularization. Though the new approach is capable of improving the shear modulus contrast in inclusions, the recovered shear modulus of the smaller inclusion is still underestimated, particularly in cases with a larger size ratio of the inclusions. In addition, we merely restricts our attention in the size sensitivity issue in this paper, but do not consider the boundary sensitivity issue. Thus, an improved and more general approach should be proposed in future work to completely resolve the issue. Besides, the displacement data can be obtained in high resolution by digital imaging correlation (DIC) [15,17,18,19]. In the future, we will test the proposed approach using experimental data from the DIC system.

5. Conclusions

This paper presents a novel inverse method to reconstruct nonhomogeneous shear modulus distribution of solids. The proposed inverse approach is capable of reducing the sensitivity of the solution of the inverse problem to the size of the inclusions. A number of numerical experiments have been conducted to test the feasibility of the novel approach. From these simulations, we observed that the novel approach significantly improved the shear modulus contrast in the inclusions. In addition, we have investigated a theoretical model and varied the size of the inclusions to validate what we have observed in simulations. This study will strengthen our understanding of the limitation of the standard, optimization-based inverse approach, which minimizes the absolute error between the measured and computed displacements in the L2 norm and triggers us to propose optimal inverse approaches to identify the nonhomogeneous mechanical property distribution of solids with high accuracy.