Proposal and validation of polyconvex strain-energy function for biological soft tissues

Abstract

BACKGROUND:

Mechanical simulations for biological tissues are effective technology for development of medical equipment, because it can be used to evaluate mechanical influences on the tissues. For such simulations, mechanical properties of biological tissues are required. For most biological soft tissues, stress tends to increase monotonically as strain increases.

OBJECTIVE:

Proposal of a strain-energy function that can guarantee monotonically increasing trend of biological soft tissue stress-strain relationships and applicability confirmation of the proposed function for biological soft tissues.

METHOD:

Based on convexity of invariants, a polyconvex strain-energy function that can reproduce monotonically increasing trend was derived. In addition, to confirm its applicability, curve-fitting of the function to stress-strain relationships of several biological soft tissues was performed.

RESULTS:

A function depending on the first invariant alone was derived. The derived function does not provide such inappropriate negative stress in the tensile region provided by several conventional strain-energy functions.

CONCLUSIONS:

The derived function can reproduce the monotonically increasing trend and is proposed as an appropriate function for biological soft tissues. In addition, as is well-known for functions depending the first invariant alone, uniaxial-compression and equibiaxial-tension of several biological soft tissues can be approximated by curve-fitting to uniaxial-tension alone using the proposed function.

Keywords

Biological tissue biomechanical simulation finite element method hyperelasticity tissue modelling

1. Introduction

Biomechanical simulation—i.e. mechanical simulation for biological tissues—can be used to verify performances of medical equipment, and is thus an effective development tool. The biomechanical simulation requires material properties that numerically represent relationships between forces applied to biological tissues and their deformations. Material properties of biological soft tissues, e.g., skin, muscle, and adipose, must be based on the finite deformation theory because of their nonlinearity and their greater deformation than that of hard tissues, e.g., bone and tooth. Ideally, nonlinearity, anisotropy, and viscoelasticity should be taken into consideration when attempting to reproduce deformation behaviors of these soft tissues. Furthermore, a coupled structural and fluid analysis, e.g., an analysis of blood vessels, requires complicated model [1]. Therefore, in this study, we propose a strain-energy function that focuses on the nonlinearity of biological soft tissues. When applying biomechanical simulation to biological soft tissues, hyperelastic models assuming incompressibility and isotropy—particularly a polynomial strain-energy function formed by Rivlin [2]—have been widely used [3–5]. The polynomial strain-energy function is implemented in many of the finite element method (FEM) systems, e.g., Abaqus, ANSYS, and Marc, and thus can be easily utilized. Therefore, we focused on this incompressible isotropic polynomial strain-energy function.

In material tests on biological soft tissues, such as uniaxial-tension tests, we can see that stress increases monotonically with increasing strain [6,7]. However, from the study of Schröder and Neff [8] showing that the second invariant is non-convex, we judged that the conventionally used polynomial strain-energy function does not guarantee monotonically increasing stress–strain relationships. Additionally, ethical issues have restricted the implementation of tests on biological soft tissues in recent years. Even if they can be done, most of them involve uniaxial-tension. Therefore, a strain-energy function used for biological soft tissues is desirable to be one that can approximate other test modes, such as uniaxial-compression and equibiaxial-tension, from uniaxial-tension alone.

In this study, we proposed a polyconvex polynomial strain-energy function, of which stress increases monotonically with increasing strain. Furthermore, to confirm an applicability of our proposed function to biological soft tissues, we used the function for curve-fitting to uniaxial-tension/compression test of porcine-liver by Chui et al. [9] and uniaxial/equibiaxial-tension test of porcine-coronary-artery by Lally et al. [10]. Compatibility of the curve-fitting results with the experimental data from the original compression and equibiaxial tests was then evaluated.

2. Theoretical considerations

2.1. Strain-energy functions

Biological soft tissues deform greatly, like rubber materials. Therefore, in order to reproduce deformation behaviors of the biological soft tissues using FEM, the strain-energy function formed by Rivlin [2] shown in Eq. ((1)) have been frequently used: $\begin{eqnarray}\displaystyle W=\displaystyle \mathop{\sum }_{i=0,j=0}^{n}c_{ij}(I_{1}-3)^{i}(I_{2}-3)^{j},\quad c_{00}=0 & & \displaystyle\end{eqnarray}$ (1) where W is the strain-energy function of hyperelastic materials, and c _ij represents substance-specific values. I ₁ and I ₂ are the first and the second invariant of the right Cauchy–Green tensor C , respectively, and are defined as follows: $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle I_{1}=\text{tr}∼\boldsymbol{C}={\lambda}_{1}^{2}+{\lambda}_{2}^{2}+{\lambda}_{3}^{2}\end{eqnarray}$ (2) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle I_{2}={\displaystyle \frac{1}{2}}\{(\text{tr}∼\boldsymbol{C})^{2}-\text{tr}∼(\boldsymbol{C}^{2})\}={\lambda}_{1}^{2}{\lambda}_{2}^{2}+{\lambda}_{2}^{2}{\lambda}_{3}^{2}+{\lambda}_{3}^{2}{\lambda}_{1}^{2},\end{eqnarray}$ (3) where 𝜆₁, 𝜆₂, and 𝜆₃, are principal stretches. The deformation gradient tensor F and the right Cauchy-Green tensor C have the following relationship: $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle \boldsymbol{F}=\left[\begin{array}{@{}ccc@{}}{\lambda}_{1} & 0 & 0\\ 0 & {\lambda}_{2} & 0\\ 0 & 0 & {\lambda}_{3}\end{array}\right]\end{eqnarray}$ (4) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle \boldsymbol{C}=\boldsymbol{F}^{T}\boldsymbol{F}=\left[\begin{array}{@{}ccc@{}}{\lambda}_{1}^{2} & 0 & 0\\ 0 & {\lambda}_{2}^{2} & 0\\ 0 & 0 & {\lambda}_{3}^{2}\end{array}\right].\end{eqnarray}$ (5) The strain-energy function W is defined as a scalar that gives the second Piola–Kirchhoff stress tensor S by partial differentiation with the Green–Lagrange strain tensor E or the right Cauchy–Green tensor C as shown in Eq. ((6)) [11]: $\begin{eqnarray}\displaystyle \boldsymbol{S}={\displaystyle \frac{\partial W}{\partial \boldsymbol{E}}}=2{\displaystyle \frac{\partial W}{\partial \boldsymbol{C}}}. & & \displaystyle\end{eqnarray}$ (6) In an incompressible condition, Eq. ((6)) is modified as shown in Eq. ((7)): $\begin{eqnarray}\displaystyle \boldsymbol{S}=-p\boldsymbol{C}^{-1}+2{\displaystyle \frac{\partial W}{\partial \boldsymbol{C}}}. & & \displaystyle\end{eqnarray}$ (7) Equation ((1)) has a problem in which the initial value of the indefinite hydrostatic pressure p in Eq. ((7)) is not zero. In order to avoid this problem, the reduced invariants J _i defined by the following equations are introduced, instead of the invariants I _i [12,13]: $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle J_{1}={\displaystyle \frac{I_{1}}{I_{3}^{1/3}}}\end{eqnarray}$ (8) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle J_{2}={\displaystyle \frac{I_{2}}{I_{3}^{2/3}}}.\end{eqnarray}$ (9) Therefore, the following equation derived from Eqs ((1)), ((8)), and ((9)) is examined: $\begin{eqnarray}\displaystyle W=\displaystyle \mathop{\sum }_{i=0,j=0}^{n}c_{ij}(J_{1}-3)^{i}(J_{2}-3)^{j},\quad c_{00}=0. & & \displaystyle\end{eqnarray}$ (10) Additionally, general-purpose FEM systems implementing the polynomial strain-energy function, e.g., ANSYS, often use the following nine-parameter function, i.e., n = 3 in Eq. ((10)): $\begin{eqnarray}\displaystyle W\text{} & =\text{} & \displaystyle c_{10}(J_{1}-3)+c_{01}(J_{2}-3)\nonumber\\ \displaystyle \text{} & \text{} & \displaystyle \quad +∼c_{20}(J_{1}-3)^{2}+c_{11}(J_{1}-3)(J_{2}-3)+c_{02}(J_{2}-3)^{2}\nonumber\\ \displaystyle \text{} & \text{} & \displaystyle \quad +∼c_{30}(J_{1}-3)^{3}+c_{21}(J_{1}-3)^{2}(J_{2}-3)+c_{12}(J_{1}-3)(J_{2}-3)^{2}+c_{03}(J_{2}-3)^{3}.\end{eqnarray}$ (11) Since Eq. ((11)) can express various well-known strain-energy functions by restrictions of its coefficients, utilization of Eq. ((11)) is convenient when considering various strain-energy functions. For instance, when c ₀₁, c ₂₀, c ₁₁, c ₀₂, c ₃₀, c ₂₁, c ₁₂, c ₀₃ = 0 or c ₂₀, c ₁₁, c ₀₂, c ₃₀, c ₂₁, c ₁₂, c ₀₃ = 0, Eq. ((11)) becomes neo-Hookean model [14,15] or Mooney–Rivlin model [16–19], respectively.

The first Piola–Kirchhoff stress, i.e., the nominal stress, must be determined based on Eq. (7). The second Piola–Kirchhoff stress tensor S and the first Piola–Kirchhoff stress tensor P have the following relationship via deformation gradient tensor F : $\begin{eqnarray}\displaystyle \boldsymbol{P}=\boldsymbol{F}\boldsymbol{S}. & & \displaystyle\end{eqnarray}$ (12) Therefore, we will determine the second Piola–Kirchhoff stress tensor S . Equation (7) is expressed by the partial differential of W by J _i and the partial differential of J _i by C as follows: $\begin{eqnarray}\displaystyle \boldsymbol{S}=-p\boldsymbol{C}^{-1}+2\left({\displaystyle \frac{\partial W}{\partial J_{1}}}{\displaystyle \frac{\partial J_{1}}{\partial \boldsymbol{C}}}+\frac{\partial W}{\partial J_{2}}{\displaystyle \frac{\partial J_{2}}{\partial \boldsymbol{C}}}\right). & & \displaystyle\end{eqnarray}$ (13) Specific ∂W∕J ₁, ∂J ₁∕ C , ∂W∕J ₂, and ∂J ₂∕ C give the second Piola–Kirchhoff stress.

2.2. Uniaxial-tension/compression

For biological soft tissues, a strain-energy function assuming incompressibility and isotropy is frequently used. Under these assumptions, the first Piola–Kirchhoff stress in the case of uniaxial-tension/compression is determined. When an incompressible isotropic object is subjected to uniaxial-tension/compression, the strain within it is uniform. The deformation gradient tensor F , the right Cauchy–Green tensor C , and its inverse tensor C ⁻¹ are given by the following equations ( F ^T is the transpose tensor of F ), with a tension/compression direction as i = 1: $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle \boldsymbol{F}=\left[\begin{array}{@{}ccc@{}}{\lambda} & 0 & 0\\ 0 & {\lambda}^{-1/2} & 0\\ 0 & 0 & {\lambda}^{-1/2}\end{array}\right]\end{eqnarray}$ (14) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle \boldsymbol{C}=\boldsymbol{F}^{T}\boldsymbol{F}=\left[\begin{array}{@{}ccc@{}}{\lambda}^{2} & 0 & 0\\ 0 & {\lambda}^{-1} & 0\\ 0 & 0 & {\lambda}^{-1}\end{array}\right]\end{eqnarray}$ (15) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle \boldsymbol{C}^{-1}=\left[\begin{array}{@{}ccc@{}}{\lambda}^{-2} & 0 & 0\\ 0 & {\lambda} & 0\\ 0 & 0 & {\lambda}\end{array}\right].\end{eqnarray}$ (16) Thus, the (reduced) invariants are given by the following equations: $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle I_{3}=1\end{eqnarray}$ (17) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle I_{1}=J_{1}={\lambda}^{2}+{\displaystyle \frac{2}{{\lambda}}}\end{eqnarray}$ (18) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle I_{2}=J_{2}=2{\lambda}+{\displaystyle \frac{1}{{\lambda}^{2}}}.\end{eqnarray}$ (19) These are substituted into Eq. ((13)), and further from the relationship S ₂₂ = S ₃₃ = 0, p in Eq. ((13)) is obtained. Thus, S ₁₁ is obtained explicitly and is expressed by the following equation: $\begin{eqnarray}\displaystyle S_{11}=2\left(1-{\displaystyle \frac{1}{{\lambda}^{3}}}\right)\left({\displaystyle \frac{\partial W}{\partial J_{1}}}+{\displaystyle \frac{1}{{\lambda}}}{\displaystyle \frac{\partial W}{\partial J_{2}}}\right). & & \displaystyle\end{eqnarray}$ (20) From Eqs ((12)), ((14)), and ((20)), the first Piola–Kirchhoff stress of uniaxial-tension/compression is expressed by the following equation: $\begin{eqnarray}\displaystyle P_{11}={\lambda}S_{11}=2\left(1-{\displaystyle \frac{1}{{\lambda}^{3}}}\right)\left({\lambda}{\displaystyle \frac{\partial W}{\partial J_{1}}}+{\displaystyle \frac{\partial W}{\partial J_{2}}}\right). & & \displaystyle\end{eqnarray}$ (21)

2.3. Equibiaxial-tension

The first Piola–Kirchhoff stress of isotropic equibiaxial-tension is determined in the same manner as the uniaxial-tension/compression. When an incompressible isotropic object is subjected to equibiaxial-tension, the strain within it is uniform. In this case, the deformation gradient tensor F , the right Cauchy-Green tensor C , and C ⁻¹ are expressed by the following equations, with tensile directions as i = 2,3: $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle \boldsymbol{F}=\left[\begin{array}{@{}ccc@{}}1/{\lambda}^{2} & 0 & 0\\ 0 & {\lambda} & 0\\ 0 & 0 & {\lambda}\end{array}\right]\end{eqnarray}$ (22) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle \boldsymbol{C}=\boldsymbol{F}^{T}\boldsymbol{F}=\left[\begin{array}{@{}ccc@{}}1/{\lambda}^{4} & 0 & 0\\ 0 & {\lambda}^{2} & 0\\ 0 & 0 & {\lambda}^{2}\end{array}\right]\end{eqnarray}$ (23) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle \boldsymbol{C}^{-1}=\left[\begin{array}{@{}ccc@{}}{\lambda}^{4} & 0 & 0\\ 0 & 1/{\lambda}^{2} & 0\\ 0 & 0 & 1/{\lambda}^{2}\end{array}\right].\end{eqnarray}$ (24) Equation ((17)) holds according to an incompressibility condition, and the (reduced) invariants are expressed by the following equations: $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle I_{1}=J_{1}=2{\lambda}^{2}+{\displaystyle \frac{1}{{\lambda}^{4}}}\end{eqnarray}$ (25) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle I_{2}=J_{2}={\lambda}^{4}+{\displaystyle \frac{2}{{\lambda}^{2}}}.\end{eqnarray}$ (26) These are substituted into Eq. ((13)), and further from S ₁₁ = 0 and S ₂₂ = S ₃₃, specific p in Eq. ((13)) is obtained. Thus, S ₂₂ (= S ₃₃) is obtained explicitly and is expressed by the following equation: $\begin{eqnarray}\displaystyle S_{22}=S_{33}=2\left({\lambda}-{\displaystyle \frac{1}{{\lambda}^{5}}}\right)\left({\displaystyle \frac{1}{{\lambda}}}{\displaystyle \frac{\partial W}{\partial J_{1}}}+{\lambda}{\displaystyle \frac{\partial W}{\partial J_{2}}}\right). & & \displaystyle\end{eqnarray}$ (27) From Eqs ((12)), ((22)), and ((27)), the first Piola–Kirchhoff stress in the case of equibiaxial-tension is given by the following equation: $\begin{eqnarray}\displaystyle P_{22}=P_{33}={\lambda}S_{22}={\lambda}S_{33}=2\left({\lambda}-{\displaystyle \frac{1}{{\lambda}^{5}}}\right)\left({\displaystyle \frac{\partial W}{\partial J_{1}}}+{\lambda}^{2}{\displaystyle \frac{\partial W}{\partial J_{2}}}\right). & & \displaystyle\end{eqnarray}$ (28)

3. Polyconvex strain-energy function

In uniaxial-tension behaviors reported by Yamada [6] and Abe et al. [7], stress of most biological soft tissues increases monotonically with increasing strain. Therefore, it is desirable for a strain-energy function to have a monotonically increasing stress–strain relationship. For a polyconvex function, a second derivative is positive. Because a second-order differential of a strain-energy function corresponds to a first-order differential of a stress–strain relationship, when a strain-energy function is polyconvex, stress increases monotonically with increasing strain.

Schröder and Neff [8] obtained the following findings with respect to polyconvexity of the reduced invariants J _i in Eq. (10) and a term including them:

•
The first reduced invariant J ₁ is polyconvex.
•
Term (J ₁ − 3)ⁿ is also polyconvex in the range of n ≥ 1.
•
The second reduced invariant J ₂ is not polyconvex.

From these characteristics, if a term having the second reduced invariant is included, the function is not polyconvex; therefore, Eqs ((10)) and ((11)) are not always polyconvex. For the polynomial strain-energy function to be always polyconvex, the coefficients of the terms including J ₂ should be set to zero, i.e., c _ij (j ≠ 0) = 0, to eliminate those terms. Additionally, the terms including J ₁ should be further defined as a positive linear combination, i.e., c _ij (j = 0) ≥ 0. That is, Eq. ((10)) is rewritten as in Eq. ((29)), namely a polyconvex strain-energy function, which is expressed as a reduced polynomial form with a coefficients constraint added. In addition, Eq. ((11)) is rewritten as in Eq. ((30)). $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle W=\displaystyle \mathop{\sum }_{i=0}^{n}c_{i0}(J_{1}-3)^{i},\quad c_{00}=0,c_{i0}\geq 0\end{eqnarray}$ (29) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle W=c_{10}(J_{1}-3)+c_{20}(J_{1}-3)^{2}+c_{30}(J_{1}-3)^{3},\quad c_{10},c_{20},c_{30}\geq 0.\end{eqnarray}$ (30) Because Eqs ((29)) and ((30)) are merely special cases of Eq. ((10)) and ((11)), it can be utilized in many general-purpose FEM systems implementing the polynomial strain-energy function.
4. Applicability of polyconvex strain-energy function

The polyconvex strain-energy function proposed in this study is a function based on the first (reduced) invariant alone. As is well-known, such functions can approximate other test modes from a uniaxial-tension test in rubber [20–23]. In order to verify applicability of the proposed function, it was examined whether other test modes such as uniaxial-compression and equibiaxial-tension can be approximated from uniaxial-tension alone. That is, curve-fitting only to the uniaxial-tension extracted from the data of Chui et al. [9] on experimental uniaxial-tension/compression of porcine-liver was performed with Eq. (30) via Eq. (21); the compression derived from the curve-fitting was then compared with the experimental finding. In addition, a similar comparison for equibiaxial-tension was performed by using the data of Lally et al. [10] on experimental uniaxial/equibiaxial-tension of porcine-coronary-artery. Since, for porcine-liver, Chui et al. [24] suggested that the nine-parameter polynomial strain-energy function, i.e., Eq. (11), can provide a good result, curve-fitting using it was also performed. Additionally, curve-fitting using Eqs (31) and (32) used for porcine-coronary-artery by Lally et al. [10] was then performed, and the results were compared with those of the proposed function: $\begin{eqnarray}\displaystyle W\text{} & =\text{} & \displaystyle c_{10}(J_{1}-3)+c_{01}(J_{2}-3)\nonumber\\ \displaystyle \text{} & \text{} & \displaystyle \quad +∼c_{20}(J_{1}-3)+c_{11}(J_{1}-3)(J_{2}-3)\nonumber\\ \displaystyle \text{} & \text{} & \displaystyle \quad +c_{30}(J_{1}-3)^{3}\end{eqnarray}$ (31) $\begin{eqnarray}\displaystyle W=c_{10}(J_{1}-3)+c_{01}(J_{2}-3)+c_{20}(J_{1}-3)^{2}. & & \displaystyle\end{eqnarray}$ (32) The residuals after curve-fitting were evaluated by the root mean square (RMS) shown in Eq. (33), where N is a number of data points: $\begin{eqnarray}\displaystyle RMS=\sqrt{{\displaystyle \frac{1}{N}}\displaystyle \mathop{\sum }_{i=1}^{N}\{\mathit{data}∼\mathit{calculated(i)}-\mathit{data}∼\mathit{measured}(i)\}^{2}}. & & \displaystyle\end{eqnarray}$ (33)

4.1. Uniaxial-tension/compression

The uniaxial-tension/compression of porcine-liver by Chui et al. [9] was converted to a numerical form by using UnGraph software (Biosoft, Cambridge, UK). The experimental tension/compression has five or six points in both the tension and compression sides. All points were converted to numerical data, and the data only on tension side were submitted to curve-fitting, using Eqs ((11)), ((30)), ((31)), and ((32)). Curve-fitting was performed by minimizing RMS in Eq. ((33)), using a solver function in Microsoft Excel 2013. The curve-fitting results are shown in Fig. 1, and the values of coefficients obtained and the RMS values derived for the tension region and the compression region are shown in Table 1.

Fig. 1.

Results of curve-fitting to uniaxial-tension alone extracted from uniaxial-tension/compression data for porcine-liver.

Table 1

Results of curve-fitting to uniaxial-tension alone of porcine-liver

	Eq. ((11))	Eq. ((30)): Polyconvex	Eq. ((31))	Eq. ((32))
c ₁₀	−1.552 × 10⁻¹	0	1.765 × 10⁻¹	−9.774 ×10⁻²
c ₀₁	1.582 × 10⁻¹	—	−2.062 × 10⁻¹	1.134 × 10⁻¹
c ₂₀	−3.227 × 10⁻¹	0	−3.468 × 10⁻¹	2.486 × 10⁻²
c ₁₁	1.847 × 10⁻¹	—	3.512 × 10⁻¹	—
c ₀₁	3.532 × 10⁻¹	—	—	—
c ₃₀	8.798 × 10⁻²	6.106 × 10⁻³	8.015 × 10⁻²	—
c ₂₁	9.654 × 10⁻³	—	—	—
c ₁₂	−2.357 × 10⁻²	—	—	—
c ₀₃	−3.203 × 10⁻²	—	—	—
RMS for tension	4.647 × 10⁻⁵	3.468 × 10⁻⁴	7.148 × 10⁻⁵	2.180 × 10⁻⁴
RMS for compression	8.929 × 10⁻¹	1.095 × 10⁻³	5.889 × 10⁻²	1.556 × 10⁻¹

[MPa]

4.2. Equibiaxial-tension

The data of Lally et al. [10] on the uniaxial/equibiaxial-tension of porcine-coronary-artery were converted to a numerical form by using UnGraph software (Biosoft, Cambridge, UK). In the report by Lally et al. [10], the experimental uniaxial-tension was shown as a solid line. Therefore, the uniaxial-tension line was converted to points that were as densely located as possible. Using the numerical data only on tension, curve-fitting calculations were performed using Eqs ((11)), ((30)), ((31)), and ((32)), as described in Section 4.1. Figures 2, 3, and Table 2 show the curve-fitting results.

Table 2
Results of curve-fitting to uniaxial-tension alone of porcine-coronary-artery

Eq. (11) Eq. (30): Polyconvex Eq. (31) Eq. (32)

c ₁₀ −1.154 ×10 1.097 × 10⁻¹ −1.365 ×10 3.457

c ₀₁ 1.219 × 10 — 1.436 × 10 −3.897

c ₂₀ 2.887 1.403 −2.549 ×10 7.349 × 10⁻¹

c ₁₁ 3.877 — 4.001 × 10 —

c ₀₁ 3.656 — — —

c ₃₀ 4.442 0 2.122 —

c ₂₁ −3.499 — — —

c ₁₂ −6.097 — — —

c ₀₃ −5.742 — — —

RMS for uniaxial-tension 2.364 × 10⁻² 1.058 × 10⁻¹ 2.323 × 10⁻² 6.918 × 10⁻²

RMS for equibiaxial-tension 5.349 2.703 × 10⁻¹ 1.095 × 10 1.964

	Eq. (11)	Eq. (30): Polyconvex	Eq. (31)	Eq. (32)
c ₁₀	−1.154 ×10	1.097 × 10⁻¹	−1.365 ×10	3.457
c ₀₁	1.219 × 10	—	1.436 × 10	−3.897
c ₂₀	2.887	1.403	−2.549 ×10	7.349 × 10⁻¹
c ₁₁	3.877	—	4.001 × 10	—
c ₀₁	3.656	—	—	—
c ₃₀	4.442	0	2.122	—
c ₂₁	−3.499	—	—	—
c ₁₂	−6.097	—	—	—
c ₀₃	−5.742	—	—	—
RMS for uniaxial-tension	2.364 × 10⁻²	1.058 × 10⁻¹	2.323 × 10⁻²	6.918 × 10⁻²
RMS for equibiaxial-tension	5.349	2.703 × 10⁻¹	1.095 × 10	1.964

[MPa]

5. Discussion

In accordance with the work of Hisada [25], if nonlinear problems in FEM use any strain-energy functions other than polyconvex ones, it is not guaranteed that the functions will have a unique solution for a given boundary condition. Multiple solutions for unique boundary conditions will possibly halt calculation in FEM, so that the use of such functions is not appropriate. When a strain-energy function is polyconvex, uniqueness and stability of solutions in FEM are guaranteed [25]. Therefore, because the proposed strain-energy function is polyconvex, stable execution of FEM can be expected.

The strain-energy function proposed here takes the same form as in the Yeoh model [20], which was proposed for carbon-black-filled rubber vulcanizates. Terms including I ₂ were eliminated in the Yeoh model on the basis of the experiments by Kawabata and Kawai [26] and Seki et al. [27], who found that sensitivity of strain-energy functions to changes of I ₂ was smaller than that to changes of I ₁. On the other hand, we eliminated terms including J ₂ on the basis of the polyconvexity of the reduced invariants and made the remaining terms positive linear combination. It is therefore interesting that the strain-energy functions having the same form were derived from completely different grounds. Yeoh showed that c ₂₀ was definitely negative when evaluating its applicability to carbon-black-filled rubber vulcanizates [20]. This characteristic is probably associated with the deformation behavior of carbon-black-filled rubber vulcanizates. In contrast, in the proposed function, c ₂₀ needs to be zero or positive. Thus, if the characteristic described above represents general deformation behaviors of rubber, it would be difficult to use our proposed function to reproduce them.

Next, we discuss the applicability of the proposed function to the uniaxial-tension/compression of porcine-liver shown in Fig. 1. According to the interpolation of the experimental points by Chui et al. [9], a stress–strain curve should be tangential to an abscissa in the vicinity of strain = 0, and stress should increase monotonically with increasing strain. However, for curve-fitting using Eqs ((11)), ((31)), and ((32)), no monotonic increase is found in the strain range from 0 to 0.28. In addition, in Eqs ((11)) and ((31)), there are inappropriate negative stress, despite in the tensile state. This is because there is no fitting object in the long interval of 0.28 in strain from zero up to the first data point, and the functions are not polyconvex. In contrast, our proposed function does not show inappropriate stress over the whole tension region and provides a monotonically increasing stress–strain relationship, because the function is polyconvex. Consideration of the compression region in Fig. 1 reveals that, because curve-fitting was performed only in the tension region, Eqs ((11)), ((31)), and ((32)) do not match the experimental compression at all, and Eq. ((31)) has inappropriate positive stress in the compression region. In contrast, the proposed function appropriately represents the experimental compression trend, i.e., monotonically increasing, and does not have inappropriate positive stress in the compression region. These results are consistent with the RMS values in Table 1. As in Yeoh’s report [20], the following equation is obtained by combining the equations for the first Piola–Kirchhoff stress in uniaxial-tension/compression, i.e., Eq. ((21)), and the proposed function, i.e., Eq. ((30)): $\begin{eqnarray}\displaystyle {\displaystyle \frac{P_{11}}{{\lambda}-{\lambda}^{-2}}}=2c_{10}+4c_{20}(J_{1}-3)+6c_{30}(J_{1}-3)^{2}. & & \displaystyle\end{eqnarray}$ (34) Figure 4 shows the relationship of (J ₁ − 3) versus P ₁₁∕(𝜆 − 𝜆⁻²) in Eq. ((34)) for the experimental tension/compression of Chui et al. [9] and the curve-fitting result to the tension alone. Because the experimental tension and compression show very similar trends, the curve-fitting result only to the uniaxial-tension approximates the trend of the compression. Considering this, for porcine-liver, the proposed function is suggested that applying to the uniaxial-tension alone can approximate the uniaxial-compression.

Fig. 2.

Approximation results for uniaxial-tension based on curve-fitting to uniaxial-tension alone extracted from experimental data on uniaxial/equibiaxial-tension of porcine-coronary-artery.

Fig. 3.

Approximation results for equibiaxial-tension based on curve-fitting to uniaxial-tension alone extracted from experimental data on uniaxial/equibiaxial-tension of porcine-coronary-artery.

Fig. 4.

Relationships of P ₁₁∕(𝜆 − 𝜆⁻²) to (J ₁ − 1) for uniaxial-tension/compression in the experimental results and the fitting curve for tension alone with respect to Eq. ((34)).

Fig. 5.

Relationships of P ₁₁∕(𝜆 − 𝜆⁻²) to (J ₁ − 1) and P ₂₂(= P ₃₃)∕(𝜆 − 𝜆⁻⁵) to (J ₁ − 1) for uniaxial/equibiaxial-tension in the experimental results and the fitting curve for tension alone with respect to Eqs ((34)) and ((35)).

We then examine the applicability of the proposed function to the uniaxial/equibiaxial-tension behaviors of porcine-coronary-artery by Lally et al. [10]. According to the study of Lally et al. [10], the stress–strain relationships for both uniaxial-tension and equibiaxial-tension of porcine-coronary-artery show a monotonic increasing trend. As shown in Fig. 2, all the functions are similar in uniaxial-tension. However, for Eq. ((32)), there is inappropriate negative stress, despite in the tensile state. This is because Eq. ((32)) is not guaranteed its monotonic increase. In contrast, our proposed function does not show inappropriate stress over the whole tension region and provides a monotonically increasing stress–strain relationship. Although Eqs ((11)) and ((31)) are not guaranteed its monotonic increase, there are not inappropriate values. Here, further consideration is given to the equibiaxial-tension in Fig. 3. Because curve-fitting was performed only in the uniaxial-tension, Eqs ((11)), ((31)), and ((32)) do not match the experimental equibiaxial-tension at all, and Eq. ((32)) has inappropriate negative stress in the equibiaxial-tension. In contrast, the proposed function well represents the feature of the experimental equibiaxial-tension, i.e., monotonically increasing, and does not have inappropriate negative stress in the equibiaxial-tension. These results are consistent with the RMS values in Table 2. As is the case with the former examination for the tension and compression, the following equation is obtained by combining the equation for the first Piola–Kirchhoff stress in the case of equibiaxial-tension, i.e., Eq. ((28)), and the proposed function, i.e., Eq. ((30)): $\begin{eqnarray}\displaystyle {\displaystyle \frac{P_{22}}{{\lambda}-{\lambda}^{-5}}}=\frac{P_{33}}{{\lambda}-{\lambda}^{-5}}=2c_{10}+4c_{20}(J_{1}-3)+6c_{30}(J_{1}-3)^{2}. & & \displaystyle\end{eqnarray}$ (35) Figure 5 shows the relationships of (J ₁ − 3) versus P ₁₁∕(𝜆 − 𝜆⁻²) and (J ₁ − 3) versus P ₂₂(= P ₃₃)∕(𝜆 − 𝜆⁻⁵) in Eqs ((34)) and ((35)) for the experimental uniaxial/equibiaxial-tension of Lally et al. [10] and the curve-fitting result to the tension alone. The trends in both uniaxial-tension and equibiaxial-tension are similar, and the curve-fitting result to the uniaxial-tension alone approximates the experimental trend in the equibiaxial-tension. Considering this, for porcine-coronary-artery, the proposed function is also suggested that applying to the uniaxial-tension alone can approximate the equibiaxial-tension.

From these, we judged that the proposed model can reproduce the deformation behavior of porcine-liver and porcine-coronary-artery more appropriately than the conventionally used strain-energy functions, i.e., Eqs (11), (31), and (32), from the viewpoints of monotonically increasing trend and uniaxial-compression/equibiaxial-tension obtained from the curve-fitting to uniaxial-tension alone.

The scope of application of our proposed function is the static deformation behavior of materials that can be assumed to be incompressible and have the monotonically increasing relationship between stress and strain, e.g., biological soft tissues. Also in consideration of the features as “stability of calculation based on the uniqueness of solution” and “reproducibility of other modes’ testing from uniaxial test alone” as described above, the proposed function can be effective used for biomechanical simulation involving various test modes, e.g., compression and equibiaxial-tension.

Many biomechanical simulations with various deformation to verify the effect of mattresses on the human body and the stress of intervertebral discs have been studied in recent years [28–44]. In such biomechanical simulation, the intervertebral discs are subjected to tension and compression with the deformation behavior of the waist bending. When only tensile tests are applied to the conventional strain-energy function, tensile stresses are calculated for compressive deformations, and appropriate solutions may not be derived. This is not the case with the proposed function, which shows compressive stress for compressive deformation. We believe that the proposed function will effectively provide the success of such static biomechanical simulations from uniaxial test alone.

6. Conclusions

In this study, we proposed a polyconvex strain-energy function. We evaluated the proposed function by applying it to uniaxial-tension/compression and uniaxial/equibiaxial-tension of biological soft tissues. The results obtained are summarized as follows:

The proposed function can reproduce a monotonically increasing trend of stress–strain relationships.

The proposed function guarantees unique solutions in FEM.

The uniaxial-compression behavior of porcine-liver can be approximated by curve-fitting to the uniaxial-tension alone using the proposed function.

The equibiaxial-tension behavior of porcine-coronary-artery can be approximated by curve-fitting to the uniaxial-tension alone using the proposed function.

From these facts, the proposed function can provide stable static biomechanical simulations including multiple deformation behaviors, e.g., tension, compression, and equibiaxial-tension, from uniaxial test alone. For example, the proposed model can be used for stress analysis of the intervertebral discs due to the bending motion of the waist that occurs when a person sleeps on a mattress. Since the proposed function can be realized by the coefficient restriction to the function implemented in general-purpose FEM systems, we believe that the proposed function greatly contributes to implementation of biomechanical simulation for biological soft tissues.

Footnotes

Acknowledgements

T. Funai is grateful to Dr. Osamu Sugiyama of the Industrial Research Institute of Shizuoka Prefecture for his help in checking the grammar of this manuscript.

Conflict of interest

None to report.

References

Nishikawa

Yokobori

A.T.J.

Ohmi

and Saiki

, Analysis of mechanical behavior of blood vessel wall on pulsatile flow of blood vessel with aneurysm, Proc Conf Tohoku Branch 49 (2014), 129–130.

Rivlin

R.S.

, Large elastic deformations, in: Rheology, Academic Press, New York, 1956, pp. 351–385.

Aranda-Iglesias

Ramón-Lozano

and Rodríguez-Martínez

J.A.

, Nonlinear resonances of an idealized saccular aneurysm, Int J Eng Sci 121 (2017), 154–166.

Hassan

M.A.

Hamdi

and Noma

, The nonlinear elastic and viscoelastic passive properties of left ventricular papillary muscle of a Guinea pig heart, J Mech Behav Biomed Mater 5(1) (2012), 99–109.

Uhlir

Mayo

Lin

P.H.

Chen

Lee

Y.-T.

Hershey

, Biomechanical characterization of the periodontal ligament: Orthodontic tooth movement, Angle Orthod 87(2) (2017), 183–192.

Yamada

, Strength of Biological Materials Evans

F.G.

(ed.), Lippincott Williams & Wilkins, US, Baltimore, 1970.

Abe

Hayashi

and Sato

, Data Book on Mechanical Properties of Living Cells, Tissues, and Organs, Springer, Tokyo, 1996.

Schröder

and Neff

, Invariant formulation of hyperelastic transverse isotropy based on polyconvex free energy functions, Int J Solids Struct 40(2) (2003), 401–445.

Chui

Kobayashi

Chen

Hisada

and Sakuma

, Transversely isotropic properties of porcine liver tissue: Experiments and constitutive modelling, Med Biol Eng Comput 45(1) (2007), 99–106.

10.

Lally

Reid

A.J.

and Prendergast

P.J.

, Elastic behavior of porcine coronary artery tissue under uniaxial and equibiaxial tension, Ann Biomed Eng 32(10) (2004), 1355–1364.

11.

Holzapfel

G.A.

, Nonlinear Solid Mechanics: A Continuum Approach for Engineering, 1st edn, Wiley, Chichester, New York, 2000.

12.

Simo

J.C.

Taylor

R.L.

and Pister

K.S.

, Variational and projection methods for the volume constraint in finite deformation elasto-plasticity, Comput Methods Appl Mech Eng 51(1–3) (1985), 177–208.

13.

Watanabe

and Hisada

, Mixed finite element analysis of incompressible hyperelastic materials, Trans Jpn Soc Mech Eng Ser A 62(595) (1996), 745–752.

14.

Treloar

L.R.G.

, The elasticity of a network of long-chain molecules—II, Trans Faraday Soc 39(0) (1943), 241–246.

15.

Treloar

L.R.G.

, The elasticity of a network of long-chain molecules. I, Trans Faraday Soc 39 (1943), 36–41.

16.

Mooney

, A theory of large elastic deformation, J Appl Phys 11(9) (1940), 582–592.

17.

Rivlin

R.S.

, Large elastic deformations of isotropic materials. IV. Further developments of the general theory, Philos Trans R Soc Math Phys Eng Sci 241(835) (1948), 379–397.

18.

Rivlin

R.S.

, Large elastic deformations of isotropic materials. V. The problem of flexure, Proc R Soc Math Phys Eng Sci 195(1043) (1949), 463–473.

19.

Rivlin

R.S.

, Large Elastic Deformations of Isotropic Materials. VI. Further Results in the Theory of Torsion, Shear and Flexure, Philos Trans R Soc Math Phys Eng Sci 242(845) (1949), 173–195.

20.

Yeoh

O.H.

, Characterization of elastic properties of carbon-black-filled rubber vulcanizates, Rubber Chem Technol 63(5) (1990), 792–805.

21.

Gough

Gregory

I.H.

and Muhr

A.H.

, Determination of constitutive equations for vulcanized rubber, in: Finite Element Analysis of Elastomers Boast

and Coveney

V.A.

(eds), Professional Engineering Pub, London, 1999, pp. 5–26.

22.

Marlow

R.S.

, A general first-invariant hyperelastic constitutive model, in: Constitutive Models for Rubber III: Proceedings of the Third European Conference on Constitutive Models for Rubber, 15–17 September 2003, London, UK, A.A. Balkema, 2003, pp. 157–160.

23.

Wex

Arndt

Stoll

Bruns

and Kupriyanova

, Isotropic incompressible hyperelastic models for modelling the mechanical behaviour of biological tissues: A review, Biomed Eng Biomed Tech 60(6) (2015).

24.

Chui

Kobayashi

Chen

Hisada

and Sakuma

, Combined compression and elongation experiments and non-linear modelling of liver tissue for surgical simulation, Med Biol Eng Comput 42(6) (2004), 787–798.

25.

Hisada

, Hisenkei yugenyousohou no tame no tensorukaiseki no kiso (Fundamentals of tensor analysis for nonlinear finite element method), Maruzen Publishing Co., Ltd., Tokyo, 1992.

26.

Kawabata

and Kawai

, Strain energy density functions of rubber vulcanizates from biaxial extension, in: Molecular Properties, Springer, Berlin, Heidelberg, 1977, pp. 89–124.

27.

Seki

Fukahori

Iseda

and Matsunaga

, A large-deformation finite-element analysis for multilayer elastomeric bearings, Rubber Chem Technol 60(5) (1987), 856–869.

28.

Funai

Suzuki

Nakamura

Yokota

Hakamata

and Maejima

, Estimating mattress deformation with a person in the supine position using biomechanical simulation, Proc Asian Pac Conf Biomech Emerg Sci Technol Biomech 8(0) (2015), 264.

29.

Lee

Won

B.H.

and Cho

S.W.

, Finite element modeling for predicting the contact pressure between a foam mattress and the human body in a supine position, Comput Methods Biomech Biomed Engin 20(1) (2017), 104–117.

30.

Yuan

and Li

, A novel method for comfort assessment in a supine sleep position using three-dimensional scanning technology, Int J Ind Ergon 67 (2018), 104–113.

31.

Rohlmann

Zander

Schmidt

Wilke

H.-J.

and Bergmann

, Analysis of the influence of disc degeneration on the mechanical behaviour of a lumbar motion segment using the finite element method, J Biomech 39(13) (2006), 2484–2490.

32.

and Wang

, Intervertebral disc biomechanical analysis using the finite element modeling based on medical images, Comput Med Imaging Graph 30(6) (2006), 363–370.

33.

Williams

J.R.

Natarajan

R.N.

and Andersson

G.B.J.

, Inclusion of regional poroelastic material properties better predicts biomechanical behavior of lumbar discs subjected to dynamic loading, J Biomech 40(9) (2007), 1981–1987.

34.

Hussain

Gay

R.E.

and An

K.-N.

, Reduction in disk and fiber stresses by axial distraction is higher in cervical disk with fibers oriented toward the vertical rather than horizontal plane: A finite element model analysis, J Manipulative Physiol Ther 33(4) (2010), 252–260.

35.

Kuo

C.-S.

H.-T.

Lin

R.-M.

Huang

K.-Y.

Lin

P.-C.

Zhong

Z.-C.

, Biomechanical analysis of the lumbar spine on facet joint force and intradiscal pressure-a finite element study, Cancel Bone (2010), 13.

36.

Chagnon

Aubin

C.-É.

and Villemure

, Biomechanical influence of disk properties on the load transfer of healthy and degenerated disks using a poroelastic finite element model, J Biomech Eng 132(11) (2010), 111006.

37.

Kamińska

Roman-Liu

Zagrajek

and Borkowski

, Differences in lumbar spine load due to posture and upper limb external load, Int J Occup Saf Ergon 16(4) (2010), 421–430.

38.

Hussain

Natarajan

R.N.

Chaudhary

H.S.

and Andersson

G.B.

, Simulation of inhomogeneous rather than homogeneous poroelastic tissue material properties within disc annulus and nucleus better predicts cervical spine response: A C3–T1 finite element model analysis under compression and moment loadings, Spine 36(4) (2011), E245.

39.

Jin

Y.J.

Kim

Y.E.

Seo

J.H.

Choi

H.W.

and Jahng

T.-A.

, Effects of rod stiffness and fusion mass on the adjacent segments after floating mono-segmental fusion: A study using finite element analysis, Eur Spine J 22(5) (2013), 1066–1077.

40.

Toumanidou

and Noailly

, Musculoskeletal modeling of the lumbar spine to explore functional interactions between back muscle loads and intervertebral disk multiphysics, Front Bioeng Biotechnol (2015), 3.

41.

Akbarian

Rouhi

Mashhadi

M.M.

and Herzog

, Biomechanics of cervical spine following implantation of a semi-constrained artificial disc with upward center of rotation: A finite element investigation, J Mech Med Biol 15(4) (2015), 1550063.

42.

Lipscomb

K.E.

Sarigul-Klijn

Klineberg

and Mohan

, Biomechanical effects of human lumbar discography: In vitro experiments and their finite element validation, Clin Spine Surg 30(3) (2017), E219.

43.

Shen

and Li

Y.-K.

, A finite element analysis of stress distribution and disk displacement in response to lumbar rotation manipulation in the sitting and side-lying positions, J Manipulative Physiol Ther 40(8) (2017), 580–586.

44.

Mills

M.J.

and Sarigul-Klijn

, Validation of an in vivo medical image-based young human lumbar spine finite element model, J Biomech Eng 141(3) (2019), 031003.