Identification of material properties of hyper-elastic materials such as soft tissues of the human body or rubber-like materials has been the subject of many works in recent decades. Boundary conditions generally play an important role in solving an inverse problem for material identification, while their knowledge has been taken for granted. In reality, however, boundary conditions may not be available on parts of the problem domain such as for an engineering part, e.g., a polymer that could be modeled as a hyper-elastic material, mounted on a system or an in vivo soft tissue. In these cases, using hypothetical boundary conditions will yield misleading results. In this paper, an inverse algorithm for the characterization of hyper-elastic material properties is developed, which takes into consideration unknown conditions on a part of the boundary. A cost function based on measured and calculated displacements is defined and is minimized using the Gauss–Newton method. A sensitivity analysis is carried out by employing analytic differentiation and using the finite element method (FEM). The effectiveness of the proposed method is demonstrated through numerical and experimental examples. The novel method is tested with a neo–Hookean and a Mooney–Rivlin hyper-elastic material model. In the experimental example, the material parameters of a silicone based specimen with unknown boundary condition are evaluated. In all the examples, the obtained results are verified and it is observed that the results are satisfactory and reliable.

Introduction

A wide range of materials with relevance in medicine and industry exhibit large elastic deformations upon application of loads. The mechanical behavior of these materials under quasi-static loadings may be described by hyper-elastic material models. Characterizing the mechanical behavior of soft tissues of the human body has great potential for improvements in virtual reality-based surgical simulation systems [14], diagnosis and surgery planning [57], trauma research [8], and detection of breast and liver cancer [9,10]. Further, identification of the constitutive behavior of rubber-like materials is very important for accurate modeling and analysis in different industrial applications, see for example [11,12]. Generally, a hyper-elastic material model is assumed to represent the solid and model parameters are determined to describe the underlying mechanical behavior.

Identification of material parameters using deformation measurements has been the subject of study for many years, and various approaches have been developed over time. A well-known and widely used approach is to fabricate a sample of the material with standard shapes and conduct simple tests, such as uniaxial or biaxial tests, to record the stress versus strain behavior, which can be fitted to a proper choice of a hyper-elastic material model.

In the following, we briefly summarize important approaches to characterize mechanical properties of biological soft tissues and rubber-like materials with various testing methods in which the materials' behavior has been modeled using hyper-elastic models.

Indentation tests along with finite element methods (FEM) have been widely used for identification of local material parameters. Samani and Plewes [13] and O'Hagan and Samani [14] used indentation test to measure hyper-elastic material properties of breast and breast tumor tissue, respectively. Tran et al. [15] modeled the three layers of human skin with a nearly incompressible neo–Hookean material and characterized their model parameters in vivo using indentation tests. Ruggiero et al. [16] and Chen et al. [17] employed an indentation procedure and an inverse method to obtain material parameters of rubber-like materials. Liu et al. [18] developed an algorithm for tissue parameters identification using a rolling indentation test. Zisis et al. [19] used instrumented indentation responses and an inverse method to characterize Mooney–Rivlin mechanical parameters of incompressible hyper-elastic materials. MacManus et al. [20] used a micro indentation test and an inverse method to characterize regional dynamic mechanical properties of brain tissue in vitro and in situ.

Other researchers have employed pipette aspiration experiments to characterize the mechanical behavior of hyper-elastic materials. Hendriks et al. [21] presented an experimental-numerical technique to characterize nonlinear Mooney–Rivlin mechanical properties of human dermis. They performed suction experiments on volar forearm skin of ten cases at various pressures and measured deformation of dermis and fat during the suction. Nava et al. [22] developed a technique to determine the mechanical properties of soft tissues of human body by applying a weak vacuum to the target tissue through a small tube. In the work of Delalleau et al. [23], the homogenous mechanical behavior of skin was identified by applying suction deformation on volar aspect of the forearm of a subject.

Other testing methods for hyper-elastic material parameter identification include propagation of elastic waves in the sample [9], inflation experiments [24], and equi-biaxial tension [25]. Mei et al. [26,27] recovered inclusions representative of tumors using only displacements known at the boundary of samples with no priori assumptions about location of inclusions, shape, and stiffness distribution.

In the work of Gambarotta et al. [28], in vivo experiments were performed during surgery on undermined scalp skin flaps to identify its isotropic hyper-elastic material parameters. Wang et al. [29] identified hyper-elastic material parameters of a domain with a breast-like geometry using undeformed and deformed configurations. Lago et al. [30] identified hyper-elastic material parameters of human cornea in vivo by using noncontact tonometry. In the work of Simón-Allué et al. [31], mechanical behavior of abdominal wall was characterized in vivo for New Zealand rabbits by performing pneumoperitoneum tests. Coordinates of several boundary points on the surface of the abdomen was recorded and used to create a three-dimensional finite element model.

In some identification techniques, such as the methods based on indentation and pipette aspiration tests, where the sample size is much greater than the dimension of the instrument contact size, boundary conditions far from the contact region may have negligible effects. Further, for testing methods based on tracking the propagation of elastic waves in a soft tissue, when the sample dimension is much greater than the wavelength of the elastic wave, effects of boundary conditions will not come into play. However, in other identification methods, which are based on elasto-static tests on the whole material, boundary conditions play an important role. In other words, a small error in boundary conditions may result in significant amplification of error in the identified material parameters [32]. Evaluating the properties of such a material is more difficult and should be carried out by special inverse analysis.

The main focus of this research is on identifying hyper-elastic material properties of a member given incomplete boundary conditions, using displacement measurements from static deformation tests. In reality, conditions at some parts of the boundary of the problem domain may be unknown and solving the inverse problem with hypothetical boundary conditions will likely result in drastic errors in the recovered material parameters. For instance, in the problem of in vivo identification of breast material properties, the boundary condition of the breast-chest interface is unknown. In all of the works discussed previously and to the best of our knowledge, in all of the previous works in this field, boundary conditions have been considered to be known everywhere for the inverse material identification problem. In the present work, an inverse algorithm for the characterization of hyper-elastic material properties is developed that handles a part of the boundary with unknown conditions. The proposed inverse method makes use of displacement measurements for the identification of these unknowns. A cost function based on measured and calculated displacements is defined and is minimized using the Gauss-Newton method. The sensitivity analysis is carried out by an analytic differentiation and using the FEM.

Constitutive Equations for Isotropic Hyper-Elastic Materials

The stress field of a perfectly elastic material does not depend on the deformation path. Hyper-elastic materials are capable of representing the nonlinear stress–strain response of elastic materials when subjected to large deformations. Hyper-elastic material modeling has received much attention to in the past and still continuous in analyzing the mechanical behavior of polymeric materials to soft tissues [3336].

There are many material models that can be used to describe the behavior of hyper-elastic materials. Some of these models such as neo–Hookean, Mooney–Rivlin, Ogden, Yeoh, and Arruda-Boyce have been developed for rubber-like materials [37,38], but have been used for soft tissues as well, see for example Refs. [14], [15], [20], and [21]. However, the nonlinear mechanical behavior of soft tissues at large deformations is usually different from rubber-like materials [39] and special hyper-elastic material models should be used for describing the mechanical behavior of soft tissues. Demiray et al. [40] and Holmes and Mow [41] used material models including exponential terms for isotropic soft tissues. Many anisotropic models have also been developed to capture the direction dependent mechanical behavior of soft tissues [39]. In this paper, neo–Hookean and Mooney–Rivlin models are employed to demonstrate that their material parameters can be recovered without knowing boundary conditions on some part of the problem domain. These material models were used by many other researchers, e.g., see references in the Introduction. However, the formulations presented in this paper can be adopted and extended for any other hyper-elastic material model. We note, however, that the number of model parameters varies with the particular hyper-elastic material model and the inverse and sensitivity analysis may behave differently for alternative hyper-elastic material models, and is beyond the scope of this paper.

Consider a deformable solid at a stress-free condition. As shown in Fig. 1, the stress-free condition is denoted by V0 and referred to as the reference configuration. Traction boundary condition (t*) and displacement boundary condition (u*) are applied on separate parts of the boundary named as V1 and V2, respectively, and induce the displacement field u(X). The deformed state of the solid is denoted by V and referred to as the current configuration. The position of a material particle in its reference and deformed configuration are denoted by X and x, respectively. Therefore, the following relation for displacement field holds: 
ui(Xk)=xiXi
(1)

Here, the subscript i and k have values 1, 2, or 3 and ui represents the three Cartesian components of u.

The deformation gradient and its Jacobian are defined as follows: 
Fij=δij+uiXj
(2)
 
J=det(F)
(3)
The left Cauchy Green deformation tensor is denoted by B and expressed as follows: 
B=FFT
(4)
Invariants of B are defined as follows: 
I1=tr(B)=Bkk
(5)
 
I2=12[(tr(B))2tr(B2)]=12(I12BikBki)
(6)
 
I3=detB=J2
(7)
where tr(B) represents the trace of B. The invariants of B can be expressed in terms of the invariants of B¯=J2/3B. Nearly incompressible hyper-elastic materials are usually modeled using invariants of B¯, which are defined as follows: 
I¯1=I1J2/3=BkkJ2/3
(8)
 
I¯2=I2J4/3=12(I¯12BikBkiJ4/3)
(9)
 
J=detB
(10)
It is assumed that a Helmholtz free energy function Ψ exists for a hyper-elastic material. Ψ is defined per unit reference volume. In a case that Ψ is only a function of strain or deformation gradient tensor, it is referred to as the strain energy function [42]. The strain energy density can be described in terms of three invariants of the left Cauchy–Green deformation tensor as follows: 
Ψ(F)=U(I1,I2,I3)=U¯(I¯1,I¯2,J)
(11)
The strain energy density function for neo–Hookean model is defined as follows [43]: 
U¯=μ12(I¯13)+K12(J1)2
(12)
in which μ1 and K1 are material parameters that can be adjusted to describe different materials.
The strain energy density for Mooney–Rivlin model is expressed as follows [43]: 
U¯=μ12(I¯13)+μ22(I¯23)+K12(J1)2
(13)

where μ1,μ2, and K1 are material parameters.

Finite Element Formulation for Direct Deformation Analysis of Hyper-Elastic Materials

For finite element formulation of the problem, the principal of virtual work is employed. This principle states that the change in potential energy during a virtual displacement from equilibrium is zero. We consider a virtual displacement field, δu(x), that satisfies the displacement boundary conditions. The principle of virtual work is then expressed as follows: 
VσijδuixjdVVρbiδuidVV2ti*δuidA=0
(14)

where σij are the components of the Cauchy stress tensor, ρ is density, and b is the force vector. The first integral is virtual strain energy and the second and third integrals are the virtual work done by body and surface forces, respectively.

The principle of virtual work expressed in Eq. (14) is based on the current configuration, which is unknown. The following relations hold between reference and current configurations for area and volume elements, respectively [42]: 
dAdA0=η=JmiFik1Fjk1mj
(15)
 
dV=JdV0
(16)

where m is the normal to the boundary in the reference configuration.

Using relations (15) and (16), the domain of the integration in Eq. (14) changes to the known reference configuration as follows: 
V0JσijδuixjdV0V0JρbiδuidV0(V0)2ti*δuiηdA0=0
(17)

where (V0)2 represents the undeformed shape of V2 in the reference configuration.

For finite element analysis, the problem domain is discretized with elements connected at nodal points in the element boundaries. Unknown displacement at each element node is denoted by uia. The superscript a ranges from 1 to n, where n is the total number of nodes. Displacement field at an arbitrary point within the solid, X, is obtained by interpolating nodal displacement values as follows: 
ui(X)=a=1nNa(X)uia
(18)
where Na(X) is the shape function of node a in the domain. The virtual displacement field can be interpolated within the solid as follows: 
δui(X)=a=1nNa(X)δuia
(19)
Substituting Eqs. (18) and (19) in Eq. (17) and performing some mathematical operations [40], the virtual work equation is obtained as follows: 
a=1n[V0τij(F)NaXmFmj1dV0V0ρ0biNadV0(V0)2ti*NaηdA0]δuia=0
(20)
where ρ0 is the density of the material in reference configuration and τ is the Kirchhoff stress tensor defined as follows: 
τij=Jσij
(21)
Equation (20) must hold for all δuia. Therefore, the following relation is obtained: 
V0τij(F)NaXmFmj1dV0V0ρ0biNadV0(V0)2ti*NaηdA0=0
(22)
Equation (22) represents 2n nonlinear equations in two-dimensional (2D) problems. The Newton–Raphson method is used in this research to solve these nonlinear equations. Using this method to solve for the increments of displacement, Eq. (22) becomes as follows [43]: 
Kaibkdukb+RiaFia=0
(23)
where Kaibk, Ria and Fia are expressed as follows: 
Kaibk=V0(τijFklNbXlNaXm)Fmj1dV0V0τijNaXmFmk1NbXpFpj1dV0(V0)2ti*NaηukbdA0
(24)
 
Ria=V0τijNaXmFmj1dV0
(25)
 
Fia=V0ρ0biNadV0+(V0)2ti*NaηdA0
(26)

Inverse Analysis

In this section, the inverse formulation to determine material parameters of a member represented by a neo–Hookean and a Mooney–Rivlin hyper-elastic material with unknown boundary conditions on a part of its boundary is presented.

In this study, unknowns of the inverse problem are of two different types. The first type of unknowns is the material parameters of the hyper-elastic model. From Eq. (12), material parameters for the neo–Hookean hyper-elastic model that must be determined are μ1 and K1. Further, from Eq. (13), μ1,μ2, and K1 are material parameters that must be determined for Mooney–Rivlin hyper-elastic model. Unknowns of the second type are boundary conditions, which are unknown displacements over a part of the boundary.

Consider a hyper-elastic material as the one shown in Fig. 2(a). The target material is attached to a different material or the same material at a part of its boundary. The boundary condition, i.e., displacement or traction is unknown on this part of the boundary. To identify the material parameters of such a hyper-elastic material, the whole set is loaded and for inverse analysis the target material is considered in a situation shown in Fig. 2(b). The points A and B on the interface with unknown condition in Fig. 2(b) are surface points and their displacements can be measured easily. However, displacement boundary conditions of all the other points on this interface are unknown. P points on this part of boundary are chosen and named as controlling points in the inverse analysis. It should be mentioned that the number of controlling points on the interface is not necessarily the same as the number of nodes on the interface. To reduce the number of unknowns of the inverse problem, only a few points, less than the number of nodes on the interface are considered as controlling points. The displacements of these controlling points are unknown in the inverse problem. In the inverse analysis, the displacement of other points on this interface is interpolated using the displacements of controlling points. In this study, the problems are solved for plane strain and plane stress conditions and therefore, two unknown displacement components are considered at each controlling point.

Measured displacement data at some points on the boundary or within the domain with known deformation, called sampling points, are used to determine the unknown data in the inverse analysis. In real cases with plane strain condition, the end boundaries (the domain of the problem in its 2D model) may not be available for any measurement, while for plane stress problems, the domain of the problem, i.e., the member's surface is accessible. Therefore, sampling points are assumed to be located on the boundary of the problem for plane strain condition, while for problems with a plane stress condition, sampling points maybe located within the domain or over the boundary. M sampling points are considered as shown in Fig. 2(b).

The vector of unknowns for neo–Hookean and Mooney–Rivlin hyper-elastic materials in the inverse analysis is respectively as follows: 
V=[μ1K1u1c1u2c1u1c2u2c2u1cPu2cP]T
(27)
 
V=[μ1μ2K1u1c1u2c1u1c2u2c2u1cPu2cP]T
(28)

where u1cj and u2cj are the first and second components of the displacement vector at the jth controlling point. Since the nature of components of the unknown vector V is of different types, the numerical values of the components of V may have different orders of magnitude. These differences may produce some difficulties in the further computations. To overcome these difficulties, we can replace the components of the vector V with dimensionless values.

The vector of measured data at the M sampling points is expressed as follows: 
S=[u¯1s1u¯2s1u¯1s2u¯2s2u¯1s3u¯2s3u¯1sMu¯2sM]T
(29)

where u¯1sj and u¯2sj are the first and second components of the displacement vector at the jth sampling point.

In the first step of the inverse analysis procedure, an initial guess for the unknown vector is selected based on past experience, approximated values, or reported values of similar materials.

After each iteration of the inverse analysis procedure, the displacements of sampling points are calculated and the vector of calculated displacements at sampling points is constructed as follows: 
C=[u1s1u2s1u1s2u2s2u1s3u2s3u1sMu2sM]T
(30)
where u1sj and u2sj are the first and second components of the computed displacement at the jth sampling point. An objective function is defined and minimized to obtain the vector of unknowns. The objective function is defined as follows: 
ψ=(SC)T(SC)
(31)
By minimizing this objective function, the differences between measured and computed data at sampling points are minimized. Different methods can be used for minimizing the objective function. In this study, the damped Gauss–Newton optimization technique [44,45] is employed, where the vector of unknowns is updated in each step as follows: 
Vk+1=Vk+αkRk
(32)

where k is the step number, αk is the step length, and Rk is the search direction. According to Ref. [44], the step length αk is first set to be one. If this choice of step length increases the value of the objective function instead of decreasing it, the step length is divided by two and the step is repeated. This process continues until the objective function in (k + 1)th step is less than the objective function in kth step.

The search direction is computed using the following equation: 
Rk=[(Lk)TLk]1[(Lk)T(SCk)]
(33)
where L is the sensitivity matrix and its components are expressed as follows: 
Lij=CiVj
(34)

The dimension of the matrix L is M×Q, where M is the number of sampling points and Q is the number of unknowns.

Different methods based on finite difference approximation and direct differentiation can be used to compute elements of the sensitivity matrix. In this research, direct differentiation formulation is derived and used in the inverse analysis. The details of the sensitivity analysis are described in Sec. 5.

Convergence rule for the optimization method is defined as follows: 
Vk+1Vke
(35)

where e is a specified tolerance.

In order to enrich measured displacement data sets, we utilize intermediate measurements while applying intermediate external loads, i.e., applying 50%, 75%, and 100% of the full load. The unknown material parameters are the same for all load cases; however, the displacement components of the controlling points on the boundary interface are different for each load case. The measured data obtained from all load cases are simultaneously used in the inverse analysis to find the unknowns. By this approach, the ratio of the total number of unknowns to the total number of measured data reduces and it is expected to get more reliable solution from the inverse analysis. This approach is more effective for problems in which the number of sampling points is not sufficiently large [4648].

The Sensitivity Analysis

An important part in solving an inverse problem is the sensitivity analysis. The objective of the sensitivity analysis in this research is to evaluate the influence of perturbed material parameters and boundary conditions (design parameters) on the displacement response. Consider Eq. (22) in an equivalent form as follows: 
Λ(u)=0
(36)
As discussed earlier, this problem is solved using an iterative Newton–Raphson method. If the solution in the current step is ui, the residual in Eq. (36) is updated in the next step as follows: 
Λ(ui+1)Λ(ui)+Λ(ui)udu=0
(37)
du is the incremental response (dukb in Eq. (23)) and Λ(ui)/u is the tangent stiffness matrix (K in Eq. (24)). The incremental response can be obtained from Eq. (37) as follows: 
Λ(ui)udu=Λ(ui)
(38)

The response of a system, also called the design response, depends on design parameters such as material constants, shape parameters, and load. The residual in Eq. (36) depends explicitly and implicitly on the design parameters. In this study, displacements at sampling points represent the design response, while material parameters and unknown displacement boundary conditions represent the design parameters. The sensitivity analysis for each design parameter should be carried out separately.

First, assume that Vα is a material parameter. In other words, Vα can be one of the first two elements of V in Eq. (27) or one of the first three elements of V in Eq. (28). The residual in Eq. (36) can be expressed as follows: 
Λ(u(Vα),Vα)=0
(39)

Assume that we have solved Eq. (36) at the end of an increment and obtained the converged solution u. The sensitivity expression is obtained by differentiating each design response with respect to design parameters.

The sensitivity of the design response (displacement at sampling points) with respect to material parameters (Vα) can be simply extracted from components of u/Vα. In order to determine u/Vα, Eq. (39) is differentiated as follows: 
ΛuuVα+ΛVα=0
(40)
Therefore, u/Vα can be obtained by solving the following equation: 
ΛuuVα=ΛVα
(41)

Comparing Eqs. (38) and (41), it is clear that they have the same tangent matrix K=Λ/u. Therefore, the evaluation of the sensitivity coefficients u/Vα requires only the evaluation of Λ/Vα and can be explicitly determined. Λ/Vα is often called the pseudo-load vector Λ̃, since it appears on the right-hand side of the sensitivity analysis equation.

Equation (41) can be rewritten as follows: 
KuVα=Λ̃
(42)

The stiffness matrix K in Eq. (42) is obtained by solving for the incremental displacement. Therefore, only the pseudo-load vector Λ̃ has to be computed in an element by element manner.

Second, assume that Vβ is a displacement component at controlling points such as given in Eq. (27) or (28). The residual in Eq. (36) can be expressed as follows: 
Λ(u(Vβ),Vβ)=0
(43)
In order to compute u/Vβ, Eq. (43) is differentiated as follows: 
ΛuuVβ=ΛVβ
(44)

Again, by comparing Eqs. (38) and (44), it is clear that they have the same tangent matrix K=Λ/u. The elements of Λ/Vβ can be obtained from the stiffness matrix directly before applying boundary conditions at controlling points, and therefore, no additional calculations are needed.

Results and Discussions

In this research, the examples are solved under plane strain/stress conditions, which have many applications in engineering, such as hyper-elastic materials mounted on structures and systems. Applications of identifying material properties of soft tissues in vivo under plane strain/stress assumptions were made, in past works see for example Refs. [49] and [50]. For three-dimensional problems, displacements of surface points can be measured, but displacement components within the domain cannot be measured using a digital image correlation system. Therefore, it is important to be able to find the unknown material parameters by using measurement data solely from surface points. To take this fact into account, in the presented inverse method, displacements at boundary points is solely assumed to be measured and used to identify unknown material properties under plane strain conditions.

At first, a theoretical/hypothetical example is numerically studied to verify the effectiveness of the presented inverse method. In this example, the effects of measurement error, location of sampling points, and values of initial guesses are investigated under plane stress and plane strain conditions for both neo–Hookean and Mooney–Rivlin hyper-elastic models. In the second example, the material properties of a member made of silicone are estimated by performing a simple tension test utilizing an INSTRON machine.

Numerical Study.

The geometry of the direct problem and its boundary conditions are shown in Fig. 3. A part of the domain of the direct problem, i.e., the rectangle ABCD, is considered as the domain of the inverse problem (Fig. 4). For the inverse problem, boundary conditions on edge AB and the hyper-elastic material parameters are unknown, while displacements at some points on edges AD, CD, and BC are known.

First, the direct problem is solved and the results from the direct problem are used as the measured data for solving an inverse problem. The exact solution of the inverse problem is known by this approach.

The input for the inverse problem are measured displacements of some boundary/domain points. The displacements assumed to be known and used to solve the inverse problem are sampling points and we introduce points on the interface AB with unknown displacement values are called controlling points. If the area with unknown boundary condition is large, more controlling points should be considered on the interface and as a result, the number of unknowns increases, which makes the inverse problem more difficult to be solved.

609 elements were used to solve the direct problem shown in Fig. 3. The rectangle ABCD, i.e., the problem domain for the inverse problem, is discretized by 17 × 15 = 255 elements. Our numerical study showed that this mesh models the problem with a sufficient accuracy. A suitable mesh size in the inverse analysis results in material parameters, which can reconstruct the measured date with sufficient accuracy.

As mentioned earlier, under plane strain condition, some of the boundaries may not be accessible for measurement and as a result, sampling points are assumed to be located only on known boundaries of the problem. This limitation makes the inverse problem more difficult. To overcome this difficulty, multiple load cases should be used under plane strain conditions to collect a rich measured data set.

A neo–Hookean material with constants μ1 = 1 Pa and K1 = 10 is considered under plane strain conditions. First, only one controlling point is considered on the interface AB as shown in Fig. 4 and a quadratic shape is assumed for the deformed shape of AB. The sampling points 1 to 4, shown in Fig. 4, are used in the inverse analysis in which three load cases with the total tensile force of 0.4, 0.7 and 1.1 N are considered. The inverse analysis did not converge for this case, which shows that the number of sampling points is not enough to solve the inverse problem. Therefore, the sampling points 1 to 6 and the same three load cases are used to solve the inverse problem. In Table 1, results from the inverse analysis for the neo–Hookean hyper-elastic material with constants μ1=1 Pa and K1=10, one controlling point and six sampling points are provided. In this case, the initial guesses for material parameters are considered to be half of their real values.

At first, material parameters are obtained using ideal inputs. However, noise is unavoidable in experimental data. To investigate the capability of the proposed method in handling noisy experimental data, different cases with a maximum error of 2%, 5%, and 10% assuming Gaussian distribution in the measured data are considered as well. In this case, the inverse analysis does not converge for measurement errors of more than 5%. Using the reconstructed material parameters reported in Table 1, the deformed shape of the interface AB is reconstructed for different measurement errors as shown in Fig. 5. As it can be seen from Table 1, for the case without measurement error, one load case is sufficient to solve the inverse problem. However, even in this case, the obtained results have some errors and the interface is not accurately reconstructed. This is due to the fact that only one controlling point is used in the inverse analysis.

For better convergence and more accuracy, the sampling points 1–8 and the same three load cases are used to solve the inverse problem. In Table 2, results of the inverse analysis for the neo–Hookean hyper-elastic material with constants μ1 = 1 Pa and K1 = 10, one controlling point and eight sampling points are provided for the plane strain condition. In this case, the initial guesses for material parameters are considered to be half of their real values. Using reconstructed material parameters from Table 2, the deformed shape of the interface AB is reconstructed for different measurement errors as shown in Fig. 6. As it can be seen from Fig. 6, the unknown shape of the interface AB is determined more precisely with eight sampling points.

The errors of the results given in Tables 1 and 2 for compared in Table 3. Since identification of the material parameters is the main concern of the inverse problem, only the error of μ1 and K1 is reported in this table. As it is observed in Table 3, for the case with six sampling points, the error is acceptable when no measurement error is present. However, for 2% and 5% measurement errors, the error in K1 is relatively large compared to the value of the measurement error.

When eight sampling points are used, even with 10% measurement error, the material constants are obtained with acceptable accuracy. The fact that the obtained constants are within an acceptable range of error shows the efficiency of the inverse method in the presence of measurement errors. Increasing the number of sampling points from 4 to 8 provides a richer data set given on other parts of the boundary, and therefore, more accurate results are obtained. An additional increase of the density of sampling points in a specific area will not provide more information and therefore does not improve the results. In summary, no matter how complex the geometry of the domain is, if the sampling points cover the material surface with adequate density, the material constants will be obtained with acceptable accuracy even with noisy measured data.

In order to examine the effectiveness of the presented method that uses measured data from different portions of the full load, the same inverse problem under plane strain conditions with one (1.1 N) and two (0.55 and 1.1 N) load cases is solved. In the presence of measurement error and considering eight sampling points, the case with one loading yields poor results and cannot minimize the objective function sufficiently. The results for the cases with two and three loadings are almost the same. For example, in the presence of 10% measurement error, the computed parameters are μ1 = 1.052 Pa and K1 = 9.46, for two and μ1 = 1.042 Pa and K1 = 9.41 for three load cases, respectively. In the following parts of this section, the results of the inverse problem corresponding to three load cases are reported.

In Table 4, results of the same inverse problem are provided for the plane strain condition when two controlling points are considered on the interface AB (Fig. 7).

Comparing the results of Tables 2 and 4 indicates that considering two controlling points instead of one improves the results of the inverse problem considerably and the unknown interface AB is determined more precisely as shown in Fig. 8. It should be noted that if a large number of controlling points is considered, the number of unknowns of the inverse problem significantly increases and we have to consider more sampling points and more load cases to obtain an accurate solution.

The inverse problem is also solved without considering any unknown parameters for the interface AB. In this case, the exact values of the displacements of the points A and B are used, while the displacements of other points on AB are linearly interpolated. The inverse problem is solved using eight sampling points and the same three load cases. The obtained results for material parameters are μ1 = 0.956 Pa and K1 = 5.56, which shows that ignoring the deformation of the interface results in poor results.

In order to investigate the effects of the initial guess and robustness of the proposed inverse method, the neo–Hookean material parameters obtained in Ref. [51] for muscle with μ1 = 23,200 Pa and K1 = 168,067.23 are considered in the next analysis. 5% measurement error is applied and using different initial guesses and three load cases, the material constants are estimated. The results are provided in Table 5. From Table 5 it is observed that, even in the cases with an initial guess far from its real value, the material constants are successfully obtained. Further, it is clear that with a better initial guess, the material constants are obtained with fewer number of iterations.

For plane stress problems, the problem domain, i.e., the member's surface is accessible; therefore, sampling points can be on the boundaries and within the domain. As a result, using one load case is sufficient to solve the inverse problem.

For Table 6, results of the inverse problem shown in Fig. 9, for the neo–Hookean material with constants μ1 = 1 Pa and K1 = 10, are provided for the plane stress condition. As can be seen from Fig. 9, ten sampling points are used to solve the inverse problem. Two of these points are in the domain. Again, the problem is solved with and without measurement errors. Three controlling points are used for modeling the deformation of the interface AB in the inverse analysis. In Table 7, the errors of the recovered constants reported in Table 6 are provided. From Table 7, it is observed that the errors of the estimated parameters lie within an acceptable range for all measurement errors.

In Table 8, results of the inverse problem are provided for the Mooney–Rivlin material with constants μ1 = 80 Pa, μ2 = 20 Pa, and K1 = 2000. The problem is solved with and without measurement errors. Ten sampling points are used and the problem is solved under plain stress condition. The initial guesses for material parameters are half of their exact values. In Table 9, the errors of the recovered constants reported in Table 8 are provided. It is observed that the error in material constants for some cases is considerably large. However, according to the publication by Rauchs et al. [52], for the Mooney–Rivlin hyper-elastic model, the individual values of μ1 and μ2 are not very important, but their summation (half the shear modulus) plays an important role for their mechanical response. From Table 9, it is observed that the error of μ1 + μ2 falls within an acceptable range.

In Fig. 10, the nominal stress–stretch curves for the Mooney−Rivlin material with constants μ1 = 80 Pa, μ2 = 20 Pa, and K1 = 2000, and the material parameters reported in Table 8 are plotted. As observed in Fig. 10, for the case without measurement error (μ1 = 72.06 Pa, μ2 = 29.56 Pa, and K1 = 1900), the nominal stress–stretch curve lies on real material response (μ1 = 80 Pa, μ2 = 20 Pa, and K1 = 2000). We also observe that for the reconstructed material parameters with 10% noise, the deformation response is almost the same as the exact deformation from the direct analysis. This confirms the fact that for the Mooney–Rivlin hyper-elastic model, the summation of μ1 and μ2 plays the major role in material behavior.

Experimental Study.

Digital image correlation (DIC) is a straightforward imaging method that provides full field displacement data on the surface of the sample under study. DIC has been used in past in combination with inverse methods to estimate material properties of silicone gel and soft tissues [24,51,53,54]. In DIC setup, several cameras are used to obtain sequential images of the sample during deformation. To correlate these images, the surface of the sample needs to be coated with a random gray-scale pattern as shown in Fig. 11. The software that analyzes DIC images traces unique features in this pattern within small facets on the surface of the sample. In order to create these features, black paint was sprayed on the surface of the sample to form a pattern that could be traced with the DIC system.

In our study, silicone samples (Ecoflex 00-10) were made and used for experiments. At first, a rectangular sample with dimensions shown in Fig. 12(a) was created. Uniaxial tension was applied to this sample using an INSTRON machine with loading rate of 10 mm per minute, which is the minimum loading rate of the machine. In order to prevent the sample from slipping between the grips, sandpaper was glued to the sample at both ends. The force-stretch curve of the rectangular sample for uniaxial tension is shown in Fig. 13. This curve is used to verify the results of the inverse problem.

Another sample with dimensions shown in Fig. 14(a) was molded in a customized wooden mold as shown in Fig. 14(b). This sample was used for the inverse problem. A DIC system (DANTEC DYNAMICS, Germany) with two cameras was used to record the deformation of the sample during simple tension up to 30 mm using INSTRON 5567 as shown in Fig. 15. During the deformation, top and bottom parts of the sample, which are in the grips of INSTRON machine, are held between wooden pieces to avoid slip. ISTRA 4D software (DANTEC DYNAMICS, Germany) was then used to analyze the images and calculate displacements of different points on the surface of the sample.

Since our purpose was to identify material properties with partially unknown boundary conditions, only the region of the sample shown in Fig. 16 is used for inverse analysis. In the inverse problem, the material parameters of this rectangular region and the boundary condition on the interface AB are unknowns to be found.

A finite element model of the rectangular domain with 517 elements in the undeformed configuration was constructed. Since the thickness of the sample (0.008 m) is small in comparison with other dimensions, the problem is considered to be under plane stress condition. The load was applied on the upper edge of the model as surface traction and displacement in x-direction at this edge was held at zero. Three controlling points with unknown displacements, i.e., six unknowns, are considered on the interface AB. The neo–Hookean hyper-elastic model was used to predict the behavior of the material.

Since DIC provides full field displacement on the surface of the member, sampling points could be considered both within the domain and on the boundary. Therefore, 12 sampling points in the domain and on the boundary as shown in Fig. 17 were used to solve the inverse problem.

The results of the inverse problem are presented in Table 10. In order to verify the results of the inverse problem, a finite element model of the rectangular sample shown in Fig. 12 was created and using the neo–Hookean material parameters of Table 10, the force-stretch curve of the model in comparison with the experimental curve is shown in Fig. 18. It is observed that the generated force-stretch curve is in good agreement with the experimental one. In other words, the material constants obtained from the inverse analysis predict the material behavior very well.

Further, for verification of the obtained neo–Hookean material parameters of Table 10, the developed finite element model of the problem shown in Fig. 17 was solved with these parameters. Displacement contour of the solution in comparison with real displacement contour obtained from ISTRA 4D software is shown in Fig. 19. From Fig. 19, it is observed that the obtained material parameters can predict material behavior very well. In Fig. 20, the unknown interface AB reconstructed from the inverse analysis is compared with the same interface obtained from the experiment using ISTRA 4D software. As it can be seen form Fig. 20, the unknown interface is reconstructed with a good accuracy.

Conclusions

An inverse method for identification of material properties of a hyper-elastic member with partially unknown boundary conditions was presented. The inverse problem was solved for neo–Hookean and Mooney–Rivlin hyper-elastic models, which are well-known models for hyper-elastic materials. Different levels of error were applied to examine the effectiveness of the method in dealing with measurement errors. The inverse problems were solved in plane stress and plane strain conditions and in both cases the results were promising. The effects of measurement error, sampling points location, and initial guess were also investigated. It was observed that when the initial guess was closer to exact values, fewer number of iterations were required to obtain the final solution. However, even in the cases with initial guess far from exact values, the method converged to a solution with sufficient accuracy.

After a numerical study of the proposed inverse method, the material parameters of a silicone specimen with unknown boundary conditions were determined by employing the presented inverse method and using real measurement data. The results were satisfactory in this case, too.

Acknowledgment

We would like to thank Dr. Terry Creasy for his Instron machine, and Dr. Arun Srinivasa and Ms. Nazanin Afsar Kazerooni for helping with the use of the Instron machine.

Funding Data

  • National Science Foundation (Grant No. CMMI #1663435).

  • The Haythornthwaite Foundation.

Nomenclature

     
  • A =

    area in deformed configuration

  •  
  • A0 =

    area in reference configuration

  •  
  • B =

    left Cauchy Green deformation tensor

  •  
  • F =

    deformation gradient

  •  
  • I1 =

    first invariant of left Cauchy Green deformation tensor

  •  
  • I2 =

    second invariant of left Cauchy Green deformation tensor

  •  
  • I3 =

    third invariant of left Cauchy Green deformation tensor

  •  
  • J =

    Jacobian of deformation gradient

  •  
  • K1 =

    material parameter of Mooney–Rivlin and neo–Hookean hyper-elastic models

  •  
  • u =

    displacement field, m

  •  
  • V =

    volume in deformed configuration

  •  
  • V0 =

    volume in reference configuration

  •  
  • x =

    position of a material particle in its deformed configuration, m

  •  
  • X =

    position of a material particle in its reference configuration, m

  •  
  • σ =

    Cauchy stress tensor, Pa

  •  
  • ρ =

    material density in deformed configuration, kg/m3

  •  
  • ρ0 =

    material density in undeformed configuration

  •  
  • Ψ =

    Helmholtz free energy function

  •  
  • μ1 =

    material parameter of Mooney–Rivlin and neo–Hookean hyper-elastic models

  •  
  • μ2 =

    material parameter of Mooney–Rivlin hyper-elastic model

References

References
1.
Basdogan
,
C.
,
Ho
,
C. H.
, and
Srinivasan
,
M. A.
,
2001
, “
Virtual Environments for Medical Training: Graphical and Haptic Simulation of Laparoscopic Common Bile Duct Exploration
,”
IEEE/ASME Trans. Mechatronics
,
6
(
3
), pp.
269
285
.
2.
Satava
,
R. M.
, and
Jones
,
S. B.
,
1997
, “
Virtual Environments for Medical Training and Education
,”
Presence: Teleoperators Virtual Environ.
,
6
(
2
), pp.
139
146
.
3.
Tendick
,
F.
,
Downes
,
M.
,
Goktekin
,
T.
,
Cavusoglu
,
M. C.
,
Feygin
,
D.
,
Wu
,
X.
,
Eyal
,
R.
,
Hegarty
,
M.
, and
Way
,
L. W.
,
2000
, “
A Virtual Environment Testbed for Training Laparoscopic Surgical Skills
,”
Presence: Teleoperators Virtual Environ.
,
9
(
3
), pp.
236
255
.
4.
Avis
,
N. J.
,
2000
, “
Virtual Environment Technologies
,”
Minimally Invasive Ther. Allied Technol.
,
9
(
5
), pp.
333
339
.
5.
Ayache
,
N.
,
Cotin
,
S.
,
Delingette
,
H.
,
Clement
,
J. M.
,
Russier
,
Y.
, and
Marescaux
,
J.
,
1998
, “
Simulation of Endoscopic Surgery
,”
Minimally Invasive Ther. Allied Technol.
,
7
(
2
), pp.
71
77
.
6.
Cotin
,
S.
,
Delingette
,
H.
, and
Ayache
,
N.
,
2000
, “
A Hybrid Elastic Model for Real-Time Cutting, Deformations, and Force Feedback for Surgery Training and Simulation
,”
Visual Comput.
,
16
(
8
), pp.
437
452
.
7.
Picinbono
,
G.
,
Lombardo
,
J. C.
,
Delingette
,
H.
, and
Ayache
,
N.
,
2002
, “
Improving Realism of a Surgery Simulator: Linear Anisotropic Elasticity, Complex Interactions and Force Extrapolation
,”
Comput. Animation Virtual Worlds
,
13
(
3
), pp.
147
167
.
8.
Snedeker
,
J. G.
,
Bajka
,
M.
,
Hug
,
J. M.
,
Szekely
,
G.
, and
Niederer
,
P.
,
2002
, “
The Creation of a High‐Fidelity Finite Element Model of the Kidney for Use in Trauma Research
,”
Comput. Animation Virtual Worlds
,
13
(
1
), pp.
53
64
.
9.
Jiang
,
Y.
,
Li
,
G. Y.
,
Qian
,
L. X.
,
Hu
,
X. D.
,
Liu
,
D.
,
Liang
,
S.
, and
Cao
,
Y.
,
2015
, “
Characterization of the Nonlinear Elastic Properties of Soft Tissues Using the Supersonic Shear Imaging (SSI) Technique: Inverse Method, Ex Vivo and In Vivo Experiments
,”
Med. Image Anal.
,
20
(
1
), pp.
97
111
.
10.
Goenezen
,
S.
,
Dord
,
J. F.
,
Sink
,
Z.
,
Barbone
,
P. E.
,
Jiang
,
J.
,
Hall
,
T. J.
, and
Oberai
,
A. A.
,
2012
, “
Linear and Nonlinear Elastic Modulus Imaging: An Application to Breast Cancer Diagnosis
,”
IEEE Trans. Med. Imaging
,
31
(
8
), pp.
1628
1637
.
11.
Mazurkiewicz
,
D.
,
2010
, “
Problems of Identification of Strength Properties of Rubber Materials for Purposes of Numerical Analysis: A Review
,”
Arch. Civ. Mech. Eng.
,
10
(
1
), pp.
69
84
.
12.
Le Saux
,
V.
,
Marco
,
Y.
,
Bles
,
G.
,
Calloch
,
S.
,
Moyne
,
S.
,
Plessis
,
S.
, and
Charrier
,
P.
,
2011
, “
Identification of Constitutive Model for Rubber Elasticity From Micro-Indentation Tests on Natural Rubber and Validation by Macroscopic Tests
,”
Mech. Mater.
,
43
(
12
), pp.
775
786
.
13.
Samani
,
A.
, and
Plewes
,
D.
,
2004
, “
A Method to Measure the Hyperelastic Parameters of Ex Vivo Breast Tissue Samples
,”
Phys. Med. Biol.
,
49
(
18
), p.
4395
.
14.
O'Hagan
,
J. J.
, and
Samani
,
A.
,
2009
, “
Measurement of the Hyperelastic Properties of 44 Pathological Ex Vivo Breast Tissue Samples
,”
Phys. Med. Biol.
,
54
(
8
), p.
2557
.
15.
Tran
,
H. V.
,
Charleux
,
F.
,
Rachik
,
M.
,
Ehrlacher
,
A.
, and
Ho Ba Tho
,
M. C.
,
2007
, “
In Vivo Characterization of the Mechanical Properties of Human Skin Derived From MRI and Indentation Techniques
,”
Comput. Methods Biomech. Biomed. Eng.
,
10
(
6
), pp.
401
407
.
16.
Ruggiero
,
L.
,
Sol
,
H.
,
Sahli
,
H.
,
Adriaenssens
,
S.
, and
Adriaenssens
,
N.
,
2011
, “
An Inverse Method to Determine Material Properties of Soft Tissues
,”
Mechanics of Biological Systems and Materials
, Vol. 2, pp.
19
32
.
17.
Chen
,
Z.
,
Scheffer
,
T.
,
Seibert
,
H.
, and
Diebels
,
S.
,
2013
, “
Macroindentation of a Soft Polymer: Identification of Hyperelasticity and Validation by Uni/Biaxial Tensile Tests
,”
Mech. Mater.
,
64
, pp.
111
127
.
18.
Liu
,
H.
,
Sangpradit
,
K.
,
Li
,
M.
,
Dasgupta
,
P.
,
Althoefer
,
K.
, and
Seneviratne
,
L. D.
,
2014
, “
Inverse Finite-Element Modeling for Tissue Parameter Identification Using a Rolling Indentation Probe
,”
Med. Biol. Eng. Comput.
,
52
(
1
), pp.
17
28
.
19.
Zisis
,
T.
,
Zafiropoulou
,
V. I.
, and
Giannakopoulos
,
A. E.
,
2015
, “
Evaluation of Material Properties of Incompressible Hyperelastic Materials Based on Instrumented Indentation of an Equal-Biaxial Prestretched Substrate
,”
Int. J. Solids Struct.
,
64–65
, pp.
132
144
.
20.
MacManus
,
D. B.
,
Pierrat
,
B.
,
Murphy
,
J. G.
, and
Gilchrist
,
M. D.
,
2016
, “
Mechanical Characterization of the P56 Mouse Brain Under Large-Deformation Dynamic Indentation
,”
Sci. Rep.
,
6
(
1
), p.
21569
.
21.
Hendriks
,
F. M.
,
Brokken
,
D. V.
,
Van Eemeren
,
J. T. W. M.
,
Oomens
,
C. W. J.
,
Baaijens
,
F. P. T.
, and
Horsten
,
J. B. A. M.
,
2003
, “
A Numerical‐Experimental Method to Characterize the Non‐Linear Mechanical Behaviour of Human Skin
,”
Skin Res. Technol.
,
9
(
3
), pp.
274
283
.
22.
Nava
,
A.
,
Mazza
,
E.
,
Kleinermann
,
F.
,
Avis
,
N. J.
, and
McClure
,
J.
,
2003
, “
Determination of the Mechanical Properties of Soft Human Tissues Through Aspiration Experiments
,”
International Conference on Medical Image Computing and Computer-Assisted Intervention
(MICCAI), Montreal, QC, Canada, Nov. 15–18, pp.
222
229
.
23.
Delalleau
,
A.
,
Josse
,
G.
,
Lagarde
,
J. M.
,
Zahouani
,
H.
, and
Bergheau
,
J. M.
,
2008
, “
A Nonlinear Elastic Behavior to Identify the Mechanical Parameters of Human Skin In Vivo
,”
Skin Res. Technol.
,
14
(
2
), pp.
152
164
.
24.
Nguyen
,
T. D.
, and
Boyce
,
B. L.
,
2011
, “
An Inverse Finite Element Method for Determining the Anisotropic Properties of the Cornea
,”
Biomech. Model. Mechanobiol.
,
10
(
3
), pp.
323
337
.
25.
Delgadillo
,
J. O. V.
,
Delorme
,
S.
,
Thibault
,
F.
,
DiRaddo
,
R.
, and
Hatzikiriakos
,
S. G.
,
2015
, “
Large Deformation Characterization of Porcine Thoracic Aortas: Inverse Modeling Fitting of Uniaxial and Biaxial Tests
,”
J. Biomed. Sci. Eng.
,
8
(
10
), p.
717
.
26.
Mei
,
Y.
,
Fulmer
,
R.
,
Raja
,
V.
,
Wang
,
S.
, and
Goenezen
,
S.
,
2016
, “
Estimating the Non-Homogeneous Elastic Modulus Distribution From Surface Deformations
,”
Int. J. Solids Struct.
,
83
, pp.
73
80
.
27.
Mei
,
Y.
,
Wang
,
S.
,
Shen
,
X.
,
Rabke
,
S.
, and
Goenezen
,
S.
,
2017
, “
Mechanics Based Tomography: A Preliminary Feasibility Study
,”
Sensors
,
17
(
5
), p.
1075
.
28.
Gambarotta
,
L.
,
Massabo
,
R.
,
Morbiducci
,
R.
,
Raposio
,
E.
, and
Santi
,
P.
,
2005
, “
In Vivo Experimental Testing and Model Identification of Human Scalp Skin
,”
J. Biomech.
,
38
(
11
), pp.
2237
2247
.
29.
Wang
,
Z. G.
,
Liu
,
Y.
,
Wang
,
G.
, and
Sun
,
L. Z.
,
2011
, “
Nonlinear Elasto-Mammography for Characterization of Breast Tissue Properties
,”
Int. J. Biomed. Imaging
,
2011
, p.
5
.
30.
Lago
,
M. A.
,
Rupérez
,
M. J.
,
Martínez-Martínez
,
F.
,
Monserrat
,
C.
,
Larra
,
E.
,
Güell
,
J. L.
, and
Peris-Martínez
,
C.
,
2015
, “
A New Methodology for the In Vivo Estimation of the Elastic Constants That Characterize the Patient-Specific Biomechanical Behavior of the Human Cornea
,”
J. Biomech.
,
48
(
1
), pp.
38
43
.
31.
Simón-Allué
,
R.
,
Calvo
,
B.
,
Oberai
,
A. A.
, and
Barbone
,
P. E.
,
2017
, “
Towards the Mechanical Characterization of Abdominal Wall by Inverse Analysis
,”
J. Mech. Behav. Biomed. Mater.
,
66
, pp.
127
137
.
32.
Roan
,
E.
, and
Vemaganti
,
K.
,
2007
, “
The Nonlinear Material Properties of Liver Tissue Determined From No-Slip Uniaxial Compression Experiments
,”
ASME J. Biomech. Eng.
,
129
(
3
), pp.
450
456
.
33.
Zhang
,
C.
,
Wu
,
J.
,
Hwang
,
K. C.
, and
Huang
,
Y.
,
2016
, “
Postbuckling of Hyperelastic Plates
,”
ASME J. Appl. Mech.
,
83
(
5
), p.
051012
.
34.
Breslavsky
,
I. D.
,
Amabili
,
M.
, and
Legrand
,
M.
,
2016
, “
Static and Dynamic Behavior of Circular Cylindrical Shell Made of Hyperelastic Arterial Material
,”
ASME J. Appl. Mech.
,
83
(
5
), p.
051002
.
35.
Zhang
,
C.
,
Wu
,
J.
, and
Hwang
,
K. C.
,
2015
, “
Hyperelastic Thin Shells: Equilibrium Equations and Boundary Conditions
,”
ASME J. Appl. Mech.
,
82
(
9
), p.
094502
.
36.
Dhavale
,
N. N.
,
Tamadapu
,
G.
, and
DasGupta
,
A.
,
2014
, “
Finite Inflation Analysis of Two Circumferentially Bonded Hyperelastic Circular Flat Membranes
,”
ASME J. Appl. Mech.
,
81
(
9
), p.
091012
.
37.
Boyce
,
M. C.
, and
Arruda
,
E. M.
,
2000
, “
Constitutive Models of Rubber Elasticity: A Review
,”
Rubber Chem. Technol.
,
73
(
3
), pp.
504
523
.
38.
Vahapoğlu
,
V.
, and
Karadeniz
,
S.
,
2006
, “
Constitutive Equations for Isotropic Rubber-Like Materials Using Phenomenological Approach: A Bibliography (1930–2003)
,”
Rubber Chem. Technol.
,
79
(
3
), pp.
489
499
.
39.
Payan
,
Y.
, and
Ohayon
,
J.
,
2017
,
Biomechanics of Living Organs: Hyperelastic Constitutive Laws for Finite Element Modeling
,
World Bank Publications
, London.
40.
Demiray
,
H.
,
Weizsäcker
,
H. W.
,
Pascale
,
K.
, and
Erbay
,
H.
,
1988
, “
A Stress-Strain Relation for a Rat Abdominal Aorta
,”
J. Biomech.
,
21
(
5
), pp.
369
374
.
41.
Holmes
,
M. H.
, and
Mow
,
V. C.
,
1990
, “
The Nonlinear Characteristics of Soft Gels and Hydrated Connective Tissues in Ultrafiltration
,”
J. Biomech.
,
23
(
11
), pp.
1145
1156
.
42.
Holzapfel
,
A. G.
,
2000
,
Nonlinear Solid Mechanics II
, Wiley, West Sussex, UK.
43.
Bower
,
A. F.
,
2009
,
Applied Mechanics of Solids
,
CRC Press
, Boca Raton, FL.
44.
Björck
,
Å.
,
1996
,
Numerical Methods for Least Squares Problems
,
Society for Industrial and Applied Mathematics
, North-Holland, The Netherlands.
45.
Ortega
,
J. M.
, and
Rheinboldt
,
W. C.
,
2000
,
Iterative Solution of Nonlinear Equations in Several Variables
,
Society for Industrial and Applied Mathematics
, Philadelphia, PA.
46.
Hematiyan
,
M. R.
,
Khosravifard
,
A.
,
Shiah
,
Y. C.
, and
Tan
,
C. L.
,
2012
, “
Identification of Material Parameters of Two-Dimensional Anisotropic Bodies Using an Inverse Multi-Loading Boundary Element Technique
,”
Comput. Model. Eng. Sci. (CMES)
,
87
(
1
), pp.
55
76
.
47.
Hematiyan
,
M. R.
,
Khosravifard
,
A.
, and
Shiah
,
Y. C.
,
2015
, “
A Novel Inverse Method for Identification of 3D Thermal Conductivity Coefficients of Anisotropic Media by the Boundary Element Analysis
,”
Int. J. Heat Mass Transfer
,
89
, pp.
685
693
.
48.
Hematiyan
,
M. R.
,
Khosravifard
,
A.
, and
Shiah
,
Y. C.
,
2017
, “
A New Stable Inverse Method for Identification of the Elastic Constants of a Three-Dimensional Generally Anisotropic Solid
,”
Int. J. Solids Struct.
,
106–107
, pp.
240
250
.
49.
Avril
,
S.
,
Bouten
,
L.
,
Dubuis
,
L.
,
Drapier
,
S.
, and
Pouget
,
J. F.
,
2010
, “
Mixed Experimental and Numerical Approach for Characterizing the Biomechanical Response of the Human Leg Under Elastic Compression
,”
ASME J. Biomech. Eng.
,
132
(
3
), p.
031006
.
50.
Franquet
,
A.
,
Avril
,
S.
,
Le Riche
,
R.
, and
Badel
,
P.
,
2012
, “
Identification of Heterogeneous Elastic Properties in Stenosed Arteries: A Numerical Plane Strain Study
,”
Comput. Methods Biomech. Biomed. Eng.
,
15
(
1
), pp.
49
58
.
51.
Affagard
,
J. S.
,
Feissel
,
P.
, and
Bensamoun
,
S. F.
,
2015
, “
Identification of Hyperelastic Properties of Passive Thigh Muscle Under Compression With an Inverse Method From a Displacement Field Measurement
,”
J. Biomech.
,
48
(
15
), pp.
4081
4086
.
52.
Rauchs
,
G.
,
Bardon
,
J.
, and
Georges
,
D.
,
2010
, “
Identification of the Material Parameters of a Viscous Hyperelastic Constitutive Law From Spherical Indentation Tests of Rubber and Validation by Tensile Tests
,”
Mech. Mater.
,
42
(
11
), pp.
961
973
.
53.
Moerman
,
K. M.
,
Holt
,
C. A.
,
Evans
,
S. L.
, and
Simms
,
C. K.
,
2009
, “
Digital Image Correlation and Finite Element Modelling as a Method to Determine Mechanical Properties of Human Soft Tissue In Vivo
,”
J. Biomech.
,
42
(
8
), pp.
1150
1153
.
54.
Genovese
,
K.
,
Casaletto
,
L.
,
Humphrey
,
J. D.
, and
Lu
,
J.
,
2014
, “
Digital Image Correlation-Based Point-Wise Inverse Characterization of Heterogeneous Material Properties of Gallbladder In Vitro
,”
Proc. R. Soc. A
,
470
(
2167
), p.
20140152
.