A novel method for the identification of the unloaded configuration of a deformed hyperelastic body

ABSTRACT

The unloaded configuration of a body without residual stresses is generally used as the starting point to compute the displacement field but may be unknown for some applications. In this paper, we propose a novel inverse formulation to identify the unloaded configuration of a deformed hyperelastic body using finite element based discretization schemes. The inverse problem consists of finding a domain that is geometrically defined by its boundary such that it deforms in a way consistent with the measurements taken on its deformed configuration. The coordinates of boundary nodes in the deformed configuration can be measured and are assumed to be known for the inverse problem. Since only the coordinates of boundary nodes in the unloaded configuration are updated in each step of the inverse analysis, re-meshing is essential to continue the optimization process. The Akin’s method is employed to map the mesh used for the unloaded configuration in the previous step for the updated one. Both numerical and experimental data sets are used in 2D/3D to demonstrate the effectiveness of the proposed inverse technique. This work could have practical applications in the design of elastomeric parts or to recover the unloaded configuration of soft tissues in vivo.

KEYWORDS:

• Unloaded configuration
• inverse problems
• finite element methods
• hyperelasticity
• Akin’s method

SUBJECT CLASSIFICATIONS CODE:

• 7405
• 74A05
• 74A10
• 74S05

1. Introduction

Given the unloaded and free of initial stresses configuration of a body, it is straightforward to determine the deformed configuration of a hyperelastic model for given material properties and boundary conditions, by solving a boundary value problem. This can be done analytically for simple problems or computationally using a numerical technique such as the finite element method (FEM) for general problems. However, determining the undeformed configuration for given current configuration and boundary conditions (loading) is a daunting task. This could be useful in designing the unloaded configuration of elastomeric parts, given the desired deformed configuration subject to loading.

In past works, Govindjee and Mihalic [1] proposed a formulation for compressible materials that employs re-parameterization of the equilibrium equations to compute the undeformed geometry of the body where the deformed shape and boundary condition are known. Govindjee and Mihalic extended their work to nearly incompressible materials using a mixed displacement-pressure formulation [2]. Further, they utilized their approach to design a gasket in a wedge channel to prevent leakage.

In the field of soft tissue biomechanics, determining the unloaded configuration of tissue has numerous applications. The unloaded configuration for tissues such as the breast is not known and has inherent deformations due to its own weight. Determining the unloaded configuration of breast tissue is essential to locate tumours during a biopsy and cosmetic surgeries [3]. The most common screening methods for tumours in breast tissue are magnetic resonance imaging (MRI) and mammography. During MRI imaging, the breast will undergo deformations due to gravity, while for mammography the breast is compressed tightly between two plates and undergoes large deformations. To compare findings from the MRI and mammogram images properly, these images could be mapped into their undeformed configurations. Further, during surgical removal of the tumour or biopsy, the patient lies in the supine position and the breast deforms during the procedure, aggravating the surgeon to locate the precise extent of the tumour. Some researchers developed finite element (FE) models of the breast to predict its deformed shape during diagnostic imaging or surgery. However, most of these models do not accurately take into account the undeformed configuration that is essential to correctly model breast deformations.

Various techniques have been employed to determine the undeformed configuration of the breast. Rajagopal et al. [4] immersed volunteers’ breasts in water in prone positions to obtain an approximate undeformed state of the breast. Lee et al. [5] simply inverted the direction of gravity to determine the unloaded breast configuration. However, since the biomechanical response of breast tissues is nonlinear, inverting gravitational forces on breast tissues is not sufficient to compute the unloaded configuration [6]. Rajagopal et al. [7,8] employed a hyperelastic material model to determine the unloaded configuration using an iterative method. The unloaded state was estimated at each iteration by minimizing residual forces. Pathmanathan et al. [9] solved a finite elasticity problem for inverted gravitational forces on breast tissue to obtain the unloaded configuration in prone position. The undeformed state was then used to simulate mammography and to predict the tumour location. Carter et al. [10] employed an iterative finite element (FE) method to compute the unloaded configuration of breast tissues from their deformed configuration in supine position. Based on this method, at first, gravity is applied in the reverse direction to obtain an initial guess for the unloaded state. Then, gravity is applied to the calculated unloaded state in its actual direction. If the locations of the nodes of the computed deformed state are sufficiently coincident with the locations of the nodes of the supine model, the estimated reference state is accepted, otherwise, the nodes of the reference state are adjusted and the process is repeated. Eiben et al. [11,12] and Eder et al. [13] employed simple iterative techniques to compute the unloaded configuration of the breast as well.

Identification of the unloaded configuration for various other tissues was also pursued by some researchers. Raghavan et al. [14] obtained the undeformed geometry of arterial aneurysms using a zero-order iterative method. The undeformed geometry was then used to determine the critical loading on the arterial wall that could lead to rupture. Bols et al. [15] employed a simple iterative FEM to obtain the zero-pressure undeformed configuration of blood vessels and utilized this configuration to calculate the stress distribution through the vessel wall under pressure loading. Riverous et al. [16] employed the same method to determine the zero-pressure geometry of abdominal aortic aneurysms using an anisotropic hyperelastic material model. Vavourakis et al. [17] developed an FE displacement-pressure formulation that re-parameterized the balance equations in the deformed state to calculate the unloaded configuration of soft tissues. Nikou et al. [18] studied the effects of employing the unloaded configuration on the accuracy of the computed in vivo diastolic properties of the heart. They presented an optimization technique to minimize the difference between deformed configurations obtained by the FEM and MRI measurements. Sellier [19] proposed a simple and successful gradient-free method to identify the unloaded configuration of a body. Rausch et al. [20] used Sellier’s method and proposed a method to determine the stress-free configuration of soft tissues.

In all of the works discussed above, the positions of the nodes within the domain of the body and on its boundary in the unloaded configuration are considered as unknowns of the inverse problem. In this study, a novel gradient-based inverse FEM is proposed that only considers positions of boundary nodes of a body with arbitrary shape as unknowns of the inverse problem. This has the advantage that the number of unknowns significantly reduces and the inverse analysis can be performed more efficiently. However, it should be noted that for each direct analysis during the inverse problem solution, the displacement components of boundary and bulk nodes are unknown. The following strategy is pursued in this paper. Suppose that the coordinates of the boundary points of the deformed configuration for a hyperelastic body are given. Next, a finite element mesh is generated. An initial guess for the unloaded configuration (pre-calculated unloaded configuration) is obtained by inversion of the applied load on the hyperelastic body. The load in its actual direction is then applied to the estimated unloaded configuration to obtain the calculated deformed configuration. An objective function is defined as the discrepancy between the measured deformed configuration and calculated deformed configuration. The damped Gauss–Newton method is used to minimize the objective function iteratively. At each minimization step, the boundary node coordinates in the unloaded configuration are updated. However, to solve the direct problem with the FEM at each step, the updated nodal coordinates of the FE mesh must be known in the entire domain. This is accomplished by utilizing a mapping approach between two consecutive steps of the inverse analysis.

2. Constitutive Equation and FE formulation

Various forms of strain energy density functions have been proposed over the decades to model the mechanical behaviour of hyperelastic materials. In this study, we have chosen the neo-Hookean and Mooney-Rivlin material models. These material models are relatively simple and widely used, but the proposed inverse method is not limited to these models and can be performed with any alternative hyperelastic material model.

The constitutive equations of a hyperelastic material can be formulated in terms of the left Cauchy Green deformation tensor $\mathbf{B}$, which is defined as follows [21]: (1) $\mathbf{B}=\mathbf{F}{\mathbf{F}}^T$(1) where $\mathbf{F}=\mathrm{\partial }\mathbf{x}/\mathrm{\partial }\mathbf{X}$ is the Green-Lagrange deformation gradient tensor. $\mathbf{X}$ and $\mathbf{x}$ represent the position of a particle in the undeformed and deformed configurations, respectively. We rewrite the left Cauchy Green deformation tensor in the form $\overline{\mathbf{B}}=\mathbf{B}/J^{2/3}$, where $J=\text{det(}\mathbf{F})$ represents the Jacobian of the deformation gradient. $\overline{\mathbf{B}}$ is more suitable for modelling nearly incompressible hyperelastic materials. The invariants of $\overline{\mathbf{B}}$ are defined as follows: (2) ${\overline{I}}_1={\overline{B}}_{kk},{\overline{I}}_2=\frac{1}{2}(\overline{I}{}_1^2-{\overline{B}}_{ik}{\overline{B}}_{ki})$(2)

The strain energy density function for the neo-Hookean hyperelastic material model is given as follows [22]: (3) $\overline{W}=\frac{{\mu }_1}{\text{2}}({\overline{I}}_1-3)+\frac{K_1}{2}(J-1{)}^2$(3) where µ₁ and K₁ are the model parameters. The strain energy density function for the Mooney-Rivlin hyperelastic model is defined as follows [22]: (4) $\overline{W}=\frac{{\mu }_1}{\text{2}}({\overline{I}}_1-3)+\frac{{\mu }_2}{\text{2}}({\overline{I}}_2-3)+\frac{K_1}{2}(J-1{)}^2$(4) where ${\mu }_1$, ${\mu }_2$ and $K_1$ are the model parameters.

Following standard procedures, the principle of virtual work is employed to obtain the weak form that can be discretized to derive the FE formulation. The weak form is expressed as follows: (5) $\underset{V}{\int }\mathit{\sigma }:\frac{\mathrm{\partial }\delta \mathbf{u}}{\mathrm{\partial }\mathbf{x}}dV-\underset{V}{\int }\rho \mathbf{b}\cdot \delta \mathbf{u}dV-\underset{\mathrm{\partial }{V}_2}{\int }{\mathbf{t}}^{\ast }\cdot \delta \mathbf{u}dA=0$(5) where V is the volume of the body in the deformed configuration, and $\mathrm{\partial }{V}_2$ is the part of the domain boundary with prescribed traction boundary condition. σ denotes the Cauchy stress tensor, b is the body force vector, ρ is the density and t* is the traction on $\mathrm{\partial }{V}_2$. In Equation (5), the dot product and tensor inner product are represented by ‘.’ and ‘:’, respectively. The first term in Equation (5) represents the virtual work of internal forces. The second and third terms in Equation (5) represent the virtual work of body and surface forces, respectively.

Next, we discretize the displacement field and virtual displacement field as follows: (6) $\mathbf{u}(\mathbf{x})=\sum _{a=1}^n{M}^a(\mathbf{x}){\mathbf{u}}^a$(6) (7) $\delta \mathbf{u}(\mathbf{x})=\sum _{a=1}^n{M}^a(\mathbf{x})\delta {\mathbf{u}}^a$(7) where $M^a(\mathbf{x})$ denotes the shape function of the a^th node and $n$ represents the total number of nodes. In this work, 4-node quadrilateral and 8-node brick elements are used for 2D and 3D problems, respectively. By Substituting Equations (6) and (7) into Equation (5), performing some mathematical manipulations, and employing the Newton–Raphson method, the following matrix equation is obtained to compute the increments of the global displacement vector $(\mathrm{\Delta }\mathbf{u})$ [22,23]: (8) $\mathbf{K}\mathrm{\Delta }\mathbf{u}+\mathbf{R}-\mathbf{F}=0$(8) where $\mathbf{K}$ represents the stiffness matrix, $\mathbf{R}$ and $\mathbf{F}$ are vectors and are expressed as follows [22,23]: (9) $\begin{array}{rl}{K}_{ik}^{ab}& =\underset{V_0}{\int }\frac{\mathrm{\partial }{\tau }_{ij}}{\mathrm{\partial }{F}_{kl}}\frac{\mathrm{\partial }{M}^b}{\mathrm{\partial }{X}_l}\frac{\mathrm{\partial }{M}^a}{\mathrm{\partial }{X}_m}{F}_{mj}^{-1}d{V}_0-\underset{V_0}{\int }{\tau }_{ij}\frac{\mathrm{\partial }{M}^a}{\mathrm{\partial }{X}_m}{F}_{mk}^{-1}\frac{\mathrm{\partial }{M}^b}{\mathrm{\partial }{X}_p}{F}_{pj}^{-1}d{V}_0\\ & -\underset{{\mathrm{\partial }}_2{V}_0}{\int }{t}_i^{\ast }{M}^a\frac{\mathrm{\partial }\eta }{\mathrm{\partial }{u}_k^b}d{A}_0\end{array}$(9) (10) $\begin{array}{rl}{R}_i^a& =\underset{V_0}{\int }{\tau }_{ij}\frac{\mathrm{\partial }{M}^a}{\mathrm{\partial }{X}_m}{F}_{mj}^{-1}d{V}_0\end{array}$(10) (11) $\begin{array}{rl}{F}_i^a& =\underset{V_0}{\int }{\rho }_0{b}_i{M}^{a}d{V}_0+\underset{{\mathrm{\partial }}_2{V}_0}{\int }{t}_i^{\ast }{M}^a\eta d{A}_0\end{array}$(11) where ${\tau }_{ij}=J{\sigma }_{ij}$ represents the Kirchhoff stress tensor. The reference and deformed area and volume of the elements are related through the relations $\eta =dA/d{A}_0=J\sqrt{m_i{F}_{ik}^{-1}{F}_{jk}^{-1}{m}_j}$ and $dV=Jd{V}_0$, respectively [24]. $m_i$ are the components of the unit normal vector in the reference configuration that is normal to the boundary of the body.

3. Inverse problem formulation

In this section, we present the inverse formulation to determine the unloaded configuration of a hyperelastic body. It is emphasized that by unloaded configuration we mean that the body is unloaded and free of initial/residual stresses. More precisely, we will solve for the boundary coordinates of the unloaded problem domain, while the material properties, the geometry of the boundary in the deformed configuration, and boundary conditions are given (see Figure 1). The unknown boundary coordinates of the unloaded configuration are defined in vector representation as follows: (12) $\mathbf{Y}=\text{[}{X}_1^{(1)}{X}_2^{(1)}{X}_3^{(1)}{X}_1^{(2)}{X}_2^{(2)}{X}_3^{(2)}\dots X_1^{(P)}{X}_2^{(P)}{X}_3^{(P)}{\text{]}}^T$(12) where $X_1^{(j)}$, $X_2^{(j)}$ and $X_3^{(j)}$ are components of the position vector of the jth node on the boundary of the body in three-dimensional space.

Figure 1. The hyperelastic body in its unloaded and deformed configurations.

The known or measured position vectors at $P$ boundary nodes are denoted by (13) $\mathbf{M}=\text{[}{\overline{x}}_1^{(1)}{\overline{x}}_2^{(1)}{\overline{x}}_3^{(1)}{\overline{x}}_1^{(2)}{\overline{x}}_2^{(2)}{\overline{x}}_3^{(2)}\dots {\overline{x}}_1^{(P)}{\overline{x}}_2^{(P)}{\overline{x}}_3^{(P)}{\text{]}}^T$(13) where ${\overline{x}}_1^{(j)}$, ${\overline{x}}_2^{(j)}$and ${\overline{x}}_3^{(j)}$ are components of the position vector of the jth node on the boundary of the deformed body.

In the presented inverse method, manual mesh generation is performed just once for the given deformed configuration. For this mesh generation, a mesh study can be performed to make sure that the quality of the mesh is in accordance with the standards of the FEM. After discretizing the deformed body by a sufficient number of finite elements, the load in its reverse direction is applied to the deformed body to obtain an initial guess of the unloaded configuration (pre-calculated unloaded configuration). In this analysis, the deformed body is assumed to be stress-free because we have no information about the actual stress distribution in the deformed body and therefore, the pre-calculated unloaded configuration is just an estimation of the real one. The deformed mesh obtained by applying the reverse load can be used for the next analysis (the mesh for the initial guess). The quality of this mesh can also be studied and if it does not satisfy the standards of the FEM, the mesh can be manually modified. In succeeding steps of the optimization procedure, changes in the estimates of the unloaded configuration between two consecutive steps are small. Therefore, the mapping procedure can be utilized to generate the mesh for the domain in the next step. As a result, the quality of the mesh is preserved and no other controls are needed. The load in its actual direction is then applied to the pre-calculated unloaded configuration (initial guess in the first step of the inverse analysis) to obtain the position of the material points in the deformed configuration. The computed position vector for $P$ boundary nodes are represented with the following notation: (14) $\mathbf{C}=\text{[}{x}_1^{(1)}{x}_2^{(1)}{x}_3^{(1)}{x}_1^{(2)}{x}_2^{(2)}{x}_3^{(2)}\dots x_1^{(P)}{x}_2^{(P)}{x}_3^{(P)}{\text{]}}^T$(14) where $x_1^{(j)}$, $x_2^{(j)}$ and $x_3^{(j)}$ are the components of the position vector at the j^th point on the boundary of the body. Next, we define a cost function as the difference between measured and computed boundary position vectors as follows: (15) $\mathrm{\Phi }=(\mathbf{M}-\mathbf{C}{)}^T(\mathbf{M}-\mathbf{C})$(15)

We minimize Equation (15) for the unknown position vector of boundary nodes in the unloaded configuration. Different optimization techniques such as the Levenberg-Marquardt or Gauss–Newton methods can be employed to minimize the objective function in Equation (15). In the Gauss–Newton method, the calculated search direction may be too large and the inverse algorithm may not converge for some cases [25]. Therefore, by adding a step factor to the conventional Gauss–Newton method, the damped Gauss–Newton method [26,27] is obtained in which by controlling the value of the step factor, the objective function in the $(k+1\text{)th}$ step is forced to be smaller than the objective function in the Kth step [26]. The damped Gauss–Newton method is employed in this study to minimize the objective function. This method is very efficient for solving inverse problems; see for example [28,29]. Hereby, the vector of unknowns is updated at each minimization step as follows: (16) ${\mathbf{Y}}^{k+1}={\mathbf{Y}}^k+{\gamma }^k{\mathbf{D}}^k$(16) where ${\gamma }^k$ is the step factor, ${\mathbf{D}}^k$ is the search direction and $k$ is the step number. According to [27], the first choice for the step factor ${\gamma }^k$ is one. If this value increases the objective function, the step factor is divided by two and the step is replicated. This procedure is repeated until the value of the cost function in the current step becomes less than its value in the previous step.

The search direction in Equation (16) is calculated as follows: (17) ${\mathbf{D}}^k=\text{[(}{\mathbf{S}}^k{)}^T{\mathbf{S}}^k{\text{]}}^{-1}\text{[(}{\mathbf{S}}^k{)}^T(\mathbf{M}-{\mathbf{C}}^k\text{)]}$(17) where $\mathbf{S}$ is the sensitivity matrix. The components of the sensitivity matrix are computed as follows: (18) $S_{ij}=\frac{\mathrm{\partial }{C}_i}{\mathrm{\partial }{Y}_j}$(18)

Various techniques based on direct differentiation and finite difference approximation can be employed to calculate the components of the sensitivity matrix. In this study, the components of the vectors $\mathbf{C}$ and $\mathbf{Y}$ are the deformed $(\mathbf{x})$ and undeformed $(\mathbf{X})$ position vectors of the boundary nodes, respectively. The sensitivity coefficients in Equation (18) for each node is equal to the components of the deformation gradient $(\mathbf{F}=\mathrm{\partial }\mathbf{x}/\mathrm{\partial }\mathbf{X})$ at that node. The deformation gradient is calculated at the element level. In the process of the problem solution, in order to compute element integrals, the components of the deformation gradient tensor are calculated at integration points. The deformation gradient tensor at nodes of the FE mesh is computed by an interpolation method using the values at the integration points. Therefore, during the solution of the direct problem in each step, the deformation gradient for boundary points obtained at the element level are stored to form the sensitivity matrix. For more clarification, suppose that we want to calculate the deformation gradient at the FE node A shown in Figure 2. This node belongs to four elements. During solution of the direct problem, the deformation gradients for all of the integration points within these four elements are stored. Therefore, by interpolating the deformation gradient tensor between 3rd integration point of element one, 4th integration point of element 2, 1st integration point of element 3 and 2nd integration point of element 4, the deformation gradient at node A is obtained. The algorithms of the present paper have been implemented in an in-house written MATLAB code.

Figure 2. Interpolating the deformation gradient tensor at a node of the FE mesh.

The optimization process for determination of the unknown unloaded configuration is repeated until a specified tolerance, $e$, e.g. 0.1 percent of the largest measured displacement, is achieved for successive solutions of the position vectors according to (19) $||{\mathbf{Y}}^{k+1}-{\mathbf{Y}}^k||\le e$(19) At each minimization step, a new estimate of the boundary node locations for the unloaded configuration is obtained and consequently, generation of interior mesh is required. For this purpose, the mesh used in the previous step is mapped for the newly calculated unloaded configuration. To this end, the Akin’s mapping method [30,31] will be used, which is based on the solution of the Laplace equation. In the following, we will briefly review this method. Let us define two consecutive unloaded configurations shown in Figure 3. Domain $A$ is obtained in the kth step and we wish to generate a similar mesh for domain $A^{\mathrm{\prime }}$ in the $(k+1\text{)th}$ step. Now, suppose that each node with coordinates $X$, $Y$ and $Z$ in domain$A$ is mapped into a node with coordinates $X^{\mathrm{\prime }}$, $Y^{\mathrm{\prime }}$ and $Z^{\mathrm{\prime }}$ in domain $A^{\mathrm{\prime }}$ with$X^{\mathrm{\prime }}=X+\mathrm{\Delta }X$, $Y^{\mathrm{\prime }}=Y+\mathrm{\Delta }Y$ and $Z^{\mathrm{\prime }}=Z+\mathrm{\Delta }Z$. Here, $\mathrm{\Delta }X$, $\mathrm{\Delta }Y$ and $\mathrm{\Delta }Z$ are known everywhere on the boundary of domain $A$. To determine $\mathrm{\Delta }X$, $\mathrm{\Delta }Y$ and $\mathrm{\Delta }Z$ for interior nodes in domain $A$, we solve Laplace’s equations given by: (20) $\frac{{\mathrm{\partial }}^2(\mathrm{\Delta }X)}{\mathrm{\partial }{X}^2}+\frac{{\mathrm{\partial }}^2(\mathrm{\Delta }X)}{\mathrm{\partial }{Y}^2}+\frac{{\mathrm{\partial }}^2(\mathrm{\Delta }X)}{\mathrm{\partial }{Z}^2}=0\text{with}\mathrm{\Delta }X\text{prescribed on the boundary of}A$(20) (21) $\frac{{\mathrm{\partial }}^2(\mathrm{\Delta }Y)}{\mathrm{\partial }{X}^2}+\frac{{\mathrm{\partial }}^2(\mathrm{\Delta }Y)}{\mathrm{\partial }{Y}^2}+\frac{{\mathrm{\partial }}^2(\mathrm{\Delta }Y)}{\mathrm{\partial }{Z}^2}=0\text{with}\mathrm{\Delta }Y\text{prescribed on the boundary of}A$(21) (22) $\frac{{\mathrm{\partial }}^2(\mathrm{\Delta }Z)}{\mathrm{\partial }{X}^2}+\frac{{\mathrm{\partial }}^2(\mathrm{\Delta }Z)}{\mathrm{\partial }{Y}^2}+\frac{{\mathrm{\partial }}^2(\mathrm{\Delta }Z)}{\mathrm{\partial }{Z}^2}=0\text{with}\mathrm{\Delta }Z\text{prescribed on the boundary of}A$(22)

Figure 3. Two consecutive calculated unloaded configurations.

Equations (20)–(22) are solved using the standard FEM [32]. Finally, we can use $\mathrm{\Delta }X$, $\mathrm{\Delta }Y$ and $\mathrm{\Delta }Z$to update the mesh in domain $A^{\mathrm{\prime }}$. We emphasize that the solution to Laplace’s equation only serves the purpose to update the mesh at each iteration and it is not for updating the deformation of the body. In the procedure of the optimization problem, updating the deformation of the body is neither performed by solving the Laplace equation nor by solving the linear elasticity problem. But, in each step of the optimization process, after performing the sensitivity analysis, the deformation of the body is updated by solving a hyperelastic problem. The algorithms of the presented method can be described as follows:

• Step 1: Mesh the given deformed configuration of the body.

• Step 2: Obtain the initial guess (pre-calculated unloaded configuration) by applying the load to the body in its reverse direction.

• Step 3: Apply the load in its actual direction to the pre-calculated unloaded configuration.

• Step 4: Compute the sensitivity matrix and update the position of the boundary nodes using the damped Gauss–Newton method.

• Step 5: Re-mesh the domain of the current estimate for the unloaded configuration.

• Step 6: Check to see if the specified tolerance is acquired. If not go to step 3.

4. Results and discussions

In this section, we present numerical examples in 2D and 3D as well as experimental examples to verify and validate the proposed method. All of the examples in this paper involve solely gravity loading, and we do not include external loading conditions that are applied actively. The main reason for this is that in many cases, actively applied loading conditions can be removed to obtain the unloaded configuration, while the same does not hold for gravity loading. Furthermore, we will investigate the effects of measurement and modelling errors on identifying the computed unloaded configuration. For the experimental part to validate this method, we chose a sample made of silicone material that is deformed under its own weight. We note that whenever displacement boundary conditions are imposed (here mainly fixed boundary conditions), the undeformed coordinates at these nodes are known, and thus, are not introduced as unknowns.

4.1. Numerical study in 2D

A hyperelastic solid with dimensions and boundary conditions as depicted in Figure 4 in its unloaded configuration is considered. The left edge of the body is fixed and the body is deformed under its own weight with $\rho =1\text{Kg}/{\text{m}}^3$. A neo-Hookean material with model parameters ${\mu }_1=2400\text{Pa}$ and $K=24000\text{Pa}$ has been chosen to represent the mechanical behaviour of the body.

Figure 4. Dimensions of the hyperelastic body with fixed boundary conditions in its unloaded configuration.

No experiments are performed for this example, and we solely focus on the verification of the proposed method with synthetic data. In this case, the result of the inverse problem (undeformed configuration of the body) is exactly known. To this end, we create synthetic data by solving a boundary value problem. From the solution of this direct problem, we extract the deformed boundary and use it as the ‘measured’ data in the inverse problem.

The mesh consists of 220 four-node quadrilateral elements and 55 nodes are defined on the boundary. Therefore, the vector of unknown boundary positions in the unloaded configuration is given by: (23) $\mathbf{Y}=\text{[}{X}_1^{(1)}{X}_2^{(1)}{X}_1^{(2)}{X}_2^{(2)}\dots X_1^{(55)}{X}_2^{(55)}{\text{]}}^T$(23) where $X_1^{(j)}$ and $X_2^{(j)}$ are components of the position vector of the j^th node on the boundary of the body. Gravity is then applied in reverse direction to the deformed configuration to obtain a pre-calculated configuration, i.e. initial guess for the inverse analysis. The unloaded (exact), deformed, and pre-calculated configurations for plane stress and plane strain conditions are depicted in Figures 5 and 6, respectively.

Figure 5. The unloaded, deformed and pre-calculated configurations of the body in the plane stress condition.

Figure 6. The unloaded, deformed and pre-calculated configurations of the body in plane strain condition.

The recovered unloaded configurations in comparison with the exact ones for plane stress and plane strain conditions are depicted in Figure 7. For both plane stress and plane strain conditions, the inverse problem recovers the unloaded configuration very accurately.

Figure 7. The recovered unloaded configuration in comparison with the exact one in (a) the plane stress condition, (b) the plane strain condition.

Knowing the exact solution of the inverse problem in this example, we can define an error term as the difference between exact and calculated unloaded configurations as follows: (24) $\sqrt{\sum _{i=1}^P{(Y_i-Y_i^0)}^2/\sum _{i=1}^P{(Y_i^0)}^2}\times 100\mathrm{\%}$(24) where $Y_i^0$ are the components of the position vector for boundary node $i$ of the exact unloaded configuration. The error calculations for both, plane stress and plane strain conditions are provided in Table 1.

Table 1. The inverse problem results for plane stress and plane strain conditions.

	No. of iterations	Error (%)
Plane stress, present method	8	0.027
Plane strain, present method	6	0.031
Plane strain, Sellier’s method [19]	45	0.039

It is observed that the computational cost for (1) the direct solver, (2) the sensitivity analysis and inverse solver, and (3) Laplace equation solver are 93.8%, 4%, and 2.2%, respectively. It should be mentioned that most of the matrices which are used in the sensitivity analysis are the same as those in the direct solver and they are not computed again.

Next, we compare the efficiency of the proposed method with the Sellier’s method [19], which is a zero-order approach. Thus, we also solve the above inverse problem for plane strain with the Sellier’s method. The Sellier’s method starts with the given deformed configuration as the initial guess for the unknown unloaded configuration. In each step of the inverse analysis, the load in its actual direction is applied to the current estimate of the unloaded configuration, resulting in an approximate deformed configuration. The difference between the calculated deformed configuration and the actual one is computed and is subtracted from the approximate unloaded configuration to update it. The procedure is continued until convergence is achieved.

The results of the inverse problem for this case are also provided in Table 1 and reveal that the proposed inverse algorithm recovers the unloaded configuration more efficiently using a gradient-based minimization approach. The computational cost of the proposed method was 6.2 times less than the Sellier’s method. It is noted that since Sellier’s method is based on a zero-order optimization method, it does not need sensitivity analysis and can be used with any direct solver code.

It is also noted that in all of the presented examples, a load, which produces a significant deformation, is considered. If a problem with a smaller deformation is considered, the problem approaches a linear problem, for which the required number of iterations by the zero-order method (Sellier’s method) approaches that of the first-order method (proposed method). On the other hand, for problems with large deformations, the nonlinearity of the problem increases and the number of iterations of Sellier’s method increases more than that of the proposed method.

Next, we ask the question of how the computational time would be affected if our method includes the coordinates of interior nodes in addition to the boundary nodes as unknowns. Thus, we solve the same inverse problem for the plane strain case with all the nodes in the FE mesh as unknowns instead of using only boundary nodes. The number of iterations and the solution error for this case are 7 and 0.028%, respectively. Our approach to optimize only for boundary nodes turns out to be 3.5 faster than optimizing for all nodes, which shows the efficiency of the proposed method.

In actual experiments, error in measurements is unavoidable. Additionally, the model used for describing the material response may not precisely represent the actual material behaviour. We will add these error sources by varying the model parameter ${\mu }_1$ of the neo-Hookean hyperelastic model. Changing ${\mu }_1$ from 2400 to 2300 Pa leads to a maximum error of 3% in the displacement field in the deformed configuration. The results for the plane strain condition are provided in Table 2. It should be noted that the exact solution corresponds to ${\mu }_1=2400\text{ Pa}$, while the inverse problem has been solved with ${\mu }_1=2300\text{Pa}$.

Table 2. Results of the inverse problem for the plane strain case with different levels of added errors in the input data.

Error in input data	No. of iterations	Solution error
3%	6	1.32%
6%	5	2.34%

Next, we assume that the actual model parameter is ${\mu }_1=2600\text{ Pa}$, while we solve the inverse problem with an assumed model parameter ${\mu }_1=2300\text{ Pa}$. This creates a maximum error of 6% in the displacement field. In Table 2, the results of the inverse problem for this case are summarized.

The unloaded configuration (exact), computed unloaded configuration, deformed configuration of the body with 6% error and without measurement error are depicted in Figure 8. From this figure it can be seen that when the amount of error in the displacements of deformed configuration is increased, the unloaded configuration is recovered with less accuracy; however, the result appears to be reasonable with a solution error of 2.34%.

Figure 8. The unloaded configuration and calculated unloaded configuration together with deformed configurations of the body with 6% error and without error in the input data for the plane strain condition.

4.2. Validation with experimental data

In this section, we will collect measurements from actual experiments and solve the inverse problem. A sample made of silicone with dimensions shown in Figure 9 was fabricated. The neo-Hookean model parameters for the silicone sample were known with ${\mu }_1=12000\text{Pa}$ and $K_1=120000\text{Pa}$. Since the thickness of the specimen was relatively small (8 mm), the problem has been solved under the plane stress condition. Further, we fixed the silicone sample at two edges as shown in Figure 10 and placed it on a graph paper with resolution of 1 mm. At first, the specimen was placed in horizontal position and leaning on the plate underneath, such that gravity did not deform it much. This state was considered to be the unloaded configuration of the body. The specimen is then placed vertically to allow gravity to bend it under its own weight as shown in Figure 11. The position of the boundary nodes in this configuration are the measured data for the inverse problem.

Figure 9. Dimensions of the silicone body.

Figure 10. The silicone body in the horizontal state.

Figure 11. The silicone body in the vertical position.

The deformed configuration of the specimen with a manually generated mesh and the pre-calculated configuration with the mapped mesh are depicted in Figure 12.

Figure 12. The deformed and pre-calculated configurations of the silicone body.

After 7 iterations, the unloaded configuration of the silicone sample was recovered as shown in Figure 13. We observe that the presented method has calculated the unloaded configuration of the silicone body with good accuracy.

Figure 13. The unloaded, deformed, and calculated unloaded configurations of the silicone body.

4.3. Numerical study in 3D

In this section, we test the proposed method with a hyperelastic material in 3D, having a relatively complicated geometry as depicted in Figure 14 in its unloaded configuration. The posterior surface of the body is fixed and the body is deformed under its own weight with $\rho =1000\text{Kg}/{\text{m}}^3$. The Mooney-Rivlin hyperelastic material model with constants ${\mu }_1=7000\text{Pa}$, ${\mu }_2=2000\text{Pa}$ and $K_1=90\text{kPa}$ is employed to simulate the mechanical behaviour of the body. First, the direct problem is solved and the positions of boundary (here a surface) nodes in the deformed configuration are extracted and utilized as measured data for the inverse problem.

Figure 14. Dimensions and boundary conditions of the 3D hyperelastic body in its unloaded configuration, (a) side view, (b) 3D view.

A mesh of 1419 eight-node brick elements is generated within the domain of the deformed body. This mesh has 487 nodes on the surface (boundary) of the body. Therefore, the vector of unknowns is given by: (25) $\mathbf{Y}=\text{[}{X}_1^{(1)}{X}_2^{(1)}{X}_3^{(1)}{X}_1^{(2)}{X}_2^{(2)}{X}_3^{(2)}\dots X_1^{(487)}{X}_2^{(487)}{X}_3^{(487)}{\text{]}}^T$(25) where $X_1^{(j)}$, $X_2^{(j)}$ and $X_3^{(j)}$ are components of the position vector of the j^th node on the surface of the body. Gravitational load is then applied in reverse direction to the deformed configuration of the body to obtain an initial guess for the inverse analysis. The unloaded, deformed and pre-calculated configurations of the body are depicted in Figure 15.

Figure 15. The unloaded, deformed and pre-calculated configurations of the body.

The inverse problem recovers the unloaded configuration of the body after 7 iterations very accurately. The recovered unloaded configuration in the x-y plane in comparison to the exact one is depicted in Figure 16.

Figure 16. The recovered unloaded configuration in x-y plane in comparison to the exact unloaded configuration.

5. Conclusions

In this paper, a gradient-based inverse method was presented that employed boundary-only variables to identify the unloaded configuration of a hyperelastic body. The Akin’s mesh mapping method was used to avoid manual mesh generation in each step of the inverse analysis. We verified the effectiveness of the method with 2D problems under plane stress and plane strain conditions in a numerical study. It was observed that the proposed gradient-based algorithm performs significantly better than the existing zero-order method. The effect of errors on identification of the unloaded configuration was also investigated. It was observed that when the error in measured data and material model increased, the unloaded configuration was obtained with less accuracy; however, it was still acceptable. To validate this approach, a silicone-based specimen was created and the proposed inverse method was employed to retrieve its unloaded configuration after being deformed under its weight. Despite the presence of inherent error in the experimental example, the unloaded configuration was recovered accurately. The efficiency of the proposed algorithm was also examined with a 3D numerical example having a relatively complicated geometry.

We note that, in this work the type of loading on the body was gravity. However, the presented method can be employed to identify the unloaded configuration of the body with any other type of loading as well. We also emphasize that the neo-Hookean and Mooney-Rivlin hyperelastic models were employed to describe mechanical behaviour of the body in this study, but, the presented inverse method can be used with other hyperelastic models as well.

Disclosure statement

No potential conflict of interest was reported by the authors.

ORCID

M. R. Hematiyan http://orcid.org/0000-0001-8926-7163