Abstract
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 [1–4], diagnosis and surgery planning [5–7], 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 [33–36].
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.
Here, the subscript i and k have values 1, 2, or 3 and ui represents the three Cartesian components of u.
where μ1,μ2, and K1 are material parameters.
Finite Element Formulation for Direct Deformation Analysis of Hyper-Elastic Materials
where 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.
where m is the normal to the boundary in the reference configuration.
where represents the undeformed shape of in the reference configuration.
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).
where and 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.
where and 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.
where k is the step number, is the step length, and is the search direction. According to Ref. [44], the step length 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 dimension of the matrix L is , 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.
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 [46–48].
The Sensitivity Analysis
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.
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.
Comparing Eqs. (38) and (41), it is clear that they have the same tangent matrix . Therefore, the evaluation of the sensitivity coefficients requires only the evaluation of and can be explicitly determined. is often called the pseudo-load vector , since it appears on the right-hand side of the sensitivity analysis equation.
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.
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.
0% meas. error | 2% meas. error | 5% meas. error | Exact value | |
---|---|---|---|---|
μ1 (Pa) | 0.995 | 0.998 | 1.023 | 1 |
K1 | 9.47 | 8.52 | 8.99 | 10 |
No. of load cases | 1 | 3 | 3 | |
No. of iterations | 9 | 12 | 13 |
0% meas. error | 2% meas. error | 5% meas. error | Exact value | |
---|---|---|---|---|
μ1 (Pa) | 0.995 | 0.998 | 1.023 | 1 |
K1 | 9.47 | 8.52 | 8.99 | 10 |
No. of load cases | 1 | 3 | 3 | |
No. of iterations | 9 | 12 | 13 |
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.
0% meas. error | 2% meas. error | 5% meas. error | 10% meas. error | Exact value | |
---|---|---|---|---|---|
μ1 (Pa) | 0.9998 | 1.009 | 1.025 | 1.052 | 1 |
K1 | 10.09 | 9.69 | 9.70 | 9.46 | 10 |
No. of load cases | 1 | 3 | 3 | 3 | |
No. of iterations | 17 | 21 | 16 | 23 |
0% meas. error | 2% meas. error | 5% meas. error | 10% meas. error | Exact value | |
---|---|---|---|---|---|
μ1 (Pa) | 0.9998 | 1.009 | 1.025 | 1.052 | 1 |
K1 | 10.09 | 9.69 | 9.70 | 9.46 | 10 |
No. of load cases | 1 | 3 | 3 | 3 | |
No. of iterations | 17 | 21 | 16 | 23 |
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.
Measurement error | 0% | 2% | 5% | 10% | ||
---|---|---|---|---|---|---|
Error of the estimated parameters | 6 Sampling points | μ1 (Pa) | 0.5% | 0.2% | 2.3% | — |
K1 | 5.3% | 14.8% | 10.1% | — | ||
8 Sampling points | μ1 (Pa) | 0.02% | 0.9% | 2.5% | 5.2% | |
K1 | 0.9% | 3.1% | 3% | 5.4% |
Measurement error | 0% | 2% | 5% | 10% | ||
---|---|---|---|---|---|---|
Error of the estimated parameters | 6 Sampling points | μ1 (Pa) | 0.5% | 0.2% | 2.3% | — |
K1 | 5.3% | 14.8% | 10.1% | — | ||
8 Sampling points | μ1 (Pa) | 0.02% | 0.9% | 2.5% | 5.2% | |
K1 | 0.9% | 3.1% | 3% | 5.4% |
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).
0% meas. error | 2% meas. error | 5% meas. error | 10% meas. error | Exact value | |
---|---|---|---|---|---|
μ1 (Pa) | 1.001 | 1.009 | 1.025 | 1.051 | 1 |
K1 | 10.05 | 10.09 | 9.85 | 9.51 | 10 |
No. of load cases | 1 | 3 | 3 | 3 | |
No. of iterations | 11 | 17 | 17 | 16 |
0% meas. error | 2% meas. error | 5% meas. error | 10% meas. error | Exact value | |
---|---|---|---|---|---|
μ1 (Pa) | 1.001 | 1.009 | 1.025 | 1.051 | 1 |
K1 | 10.05 | 10.09 | 9.85 | 9.51 | 10 |
No. of load cases | 1 | 3 | 3 | 3 | |
No. of iterations | 11 | 17 | 17 | 16 |
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.
Initial guess | No. of iterations | μ1 (Pa) (error %) | K1 (error %) |
30% of real values | 13 | 23,829 (2.71) | 170,576 (1.49) |
50% of real values | 11 | 23,805 (2.61) | 169,026 (0.57) |
80% of real values | 7 | 23,774 (2.47) | 165,367 (1.61) |
Initial guess | No. of iterations | μ1 (Pa) (error %) | K1 (error %) |
30% of real values | 13 | 23,829 (2.71) | 170,576 (1.49) |
50% of real values | 11 | 23,805 (2.61) | 169,026 (0.57) |
80% of real values | 7 | 23,774 (2.47) | 165,367 (1.61) |
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.
0% meas. error | 2% meas. error | 5% meas. error | 10% meas. error | Exact value | |
---|---|---|---|---|---|
μ1 (Pa) | 0.999 | 1.008 | 1.025 | 1.045 | 1 |
K1 | 9.87 | 9.75 | 9.33 | 9.13 | 10 |
No. of Iterations | 21 | 10 | 14 | 12 |
0% meas. error | 2% meas. error | 5% meas. error | 10% meas. error | Exact value | |
---|---|---|---|---|---|
μ1 (Pa) | 0.999 | 1.008 | 1.025 | 1.045 | 1 |
K1 | 9.87 | 9.75 | 9.33 | 9.13 | 10 |
No. of Iterations | 21 | 10 | 14 | 12 |
Measurement error | 0% | 2% | 5% | 10% | |
---|---|---|---|---|---|
Error of the estimated parameter | μ1 (Pa) | 0.1% | 0.8% | 2.5% | 4.5% |
K1 | 1.3% | 2.5% | 6.7% | 8.7% |
Measurement error | 0% | 2% | 5% | 10% | |
---|---|---|---|---|---|
Error of the estimated parameter | μ1 (Pa) | 0.1% | 0.8% | 2.5% | 4.5% |
K1 | 1.3% | 2.5% | 6.7% | 8.7% |
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.
0% meas. error | 2% meas. error | 5% meas. error | 10% meas. error | Exact value | |
---|---|---|---|---|---|
μ1 (Pa) | 72.06 | 80.13 | 89.20 | 92.35 | 80 |
μ2 (Pa) | 29.56 | 20.77 | 12.25 | 11.00 | 20 |
μ1 + μ2 | 101.62 | 100.90 | 101.45 | 103.35 | 100 |
K1 | 1900 | 1930 | 1880 | 1920 | 2000 |
No. of iterations | 16 | 14 | 12 | 13 |
0% meas. error | 2% meas. error | 5% meas. error | 10% meas. error | Exact value | |
---|---|---|---|---|---|
μ1 (Pa) | 72.06 | 80.13 | 89.20 | 92.35 | 80 |
μ2 (Pa) | 29.56 | 20.77 | 12.25 | 11.00 | 20 |
μ1 + μ2 | 101.62 | 100.90 | 101.45 | 103.35 | 100 |
K1 | 1900 | 1930 | 1880 | 1920 | 2000 |
No. of iterations | 16 | 14 | 12 | 13 |
Measurement error | 0% | 2% | 5% | 10% | |
---|---|---|---|---|---|
Error of the estimated parameter | μ1 (Pa) | 9.93% | 016% | 11.5% | 15.44% |
μ2 (Pa) | 47.8% | 3.85% | 38.75% | 45% | |
μ1+μ2 | 1.62% | 0.9% | 1.45% | 3.35% | |
K1 | 5% | 3.5% | 6% | 4% |
Measurement error | 0% | 2% | 5% | 10% | |
---|---|---|---|---|---|
Error of the estimated parameter | μ1 (Pa) | 9.93% | 016% | 11.5% | 15.44% |
μ2 (Pa) | 47.8% | 3.85% | 38.75% | 45% | |
μ1+μ2 | 1.62% | 0.9% | 1.45% | 3.35% | |
K1 | 5% | 3.5% | 6% | 4% |
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