Kỷ yếu hội nghị khoa học và công nghệ toàn quốc về cơ khí - Lần thứ IV

STABLE AND PHYSICAL RESPONSE OF CONVEX HYPERELASTIC
MODELS FOR FIBRE-REINFORCED MATERIALS
ỨNG XỬ TỰ NHIÊN VÀ ỔN ĐỊNH CỦA CÁC MÔ HÌNH SIÊU ĐÀN HỒI LỒI
MÔ TẢ CÁC VẬT LIỆU COMPOSIT CỐT SỢI
Minh Tuan Duong1a, Thi Thanh Hai Tran1b
Hanoi University of Science and Technology, Hanoi, Vietnam
a
;
1

ABSTRACT
Fibre-reinforced hyperelastic models characterising rubbers, fibre-reinforced
elastomers, and biological soft tissues (e.g. arteries), can be unstable and lead to unphysical
behaviour. This problem is due to the ill-conditioned constitutive matrix of the model. Thus,
the paper analyzes this numerical instability by investigating a hyperelastic model in which
the fibres characterised by a 2D Fung-type potential and an isotropic term representing the
matrix, namely the modified Fung-type one (MFH). The unrealistic response of the convex
MFH might take place in physiological range of soft tissues or in interested range of
hyperelastic materials. Several tension tests are conducted with the material constants of the
MFH obtained from curve fitting to experimental data of porcine arterial samples by G. A.
Holzapfel, which show poor stress solutions associated with unrealistic thickness thickening
of the tested specimens. Consequently, the paper proposes a novel Fung-type potential which
solves this numerical instability and can ensure physical response in the physiological
deformation range of the artery.
Keywords: Unphysical response, Ill-conditioning, FEM, Hyperelastic models, Fibrereinforced materials.
TÓM TẮT
Ứng xử phi tự nhiên có thể xảy ra đối với các mô hình siêu đàn hồi dùng cho mô tả cao
su, vật liệu dẻo cốt sợi hay mô sinh học. Vấn đề này xảy ra là do ma trận cấu thành của mô
hình vật liệu ở dạng điều kiện yếu. Vì vậy bài báo phân tích sự mất ổn định số này bằng việc

this unrealistic phenomenon also takes place in an equi-biaxial test using the MFH model and
the Holzapfel model in which the fibre directions are unchanged under deformation, see [6].
Therefore, it is needed to analyze this numerical problem to ensure the correct simulations.
This paper presents further solutions to this problem, especially for the MFH model. To this
end, we propose a novel model comprising a new strain energy component which is capable
of solving the unphysical response of the MFH model in the interested deformation range.
2. HYPERELASTIC MATERIAL MODEL
The stress tensors of homogeneous hyperelastic materials are directly derived (by means
of the Clausius-Planck form of the second law of thermodynamics) from a given scalar-valued
strain-energy function

S=2

∂W (C) ∂W (E)
=
,
∂C
∂E

(1)

where E is the Cauchy-Green tensor, S is the second Piola-Kirchhoff stress tensor and
the left Cauchy tensor C = FT F with the deformation gradient tensor F . The Cauchy stress
tensor σ = J −1FSFT , in which J is the determinant of the deformation gradient (volume
ratio). These types of equations are known as constitutive equations. The constitutive
equations for composite materials are expressed in a finite element model as in the following.
2.1. Finite element formulation
The total Lagrangian formulation is used since considered tensors in the previous
section are described in the reference configuration. A variational theorem for finite elasticity
may be written using hyperelastic materials as

δ uT ⋅  ∫ BT SdV −
Ω
 0


T
T
−
N
b
dV
N
T
dS
0,
=
0
∫
∫

Ω0
Γ0


(3)

where B and N are the displacement deformation matrix and the shape function
matrix, respectively. In general, solving nonlinear problems needs to adopt load increments
through a load factor λ. The principle of virtual work with the external load factor λ reads as


Ω0

∫N

Ω0

T

b0 dV − ∫ N T T dS

(6)

Γ0

The system in (5) is nonlinear discrete governing equation because of both
geometrically and physically nonlinearity from large deformation and the constitutive relation
of the material, respectively.
As a result, the total classical tangent stiffness matrix K t must be defined as
K t ∆u = R ,

Kt =

∂Fint (u)
∂F (u)
− λ ext
.
∂u
∂u

(7)

reinforcing fibres embedded in a soft material. Numerically, this problem is due to the illconditioned constitutive matrix of the model [6].

Figure 1: Fibre families of an artery (without intima) and a strip for tension test
When a 2D Fung potential models fibres in a plane Oθ z then the model with the fibre
 
orientation vectors a 0 , g 0 in Figure 1, the MFH model [3] is written as

) µ ( I1 − 3) + C exp {Q fibre (E)} − 1 ;
WMFH ( E=


2
892

Q fibre (E=
) c1 Eθθ2 + c2 Ezz2 + c3 Eθθ Ezz (9)


Kỷ yếu hội nghị khoa học và công nghệ toàn quốc về cơ khí - Lần thứ IV
where Err , Eθθ , and Ezz are principal components of the Green deformation strain E
with respect to reference configuration in Figure 1; the invariant I1 of the right Cauchy-Green
tensor C : I1 = tr C ; µ equivalent to the small strain shear modulus; c1 , c2 , and c3 are
dimensionless parameters so that the angle between two fibre families along the hoop
2c3
direction is defined by tan(2α ) =
, tan(2α l )(c1 − c2 ) ≤ 2c3 ≤ tan(2α u )(c1 − c2 )
c1 − c2
where α l and α u are the lower and upper limits for the fibre angle, respectively. The
component Waniso contributes mechanical properties only in the fibre plane Oθ z as an
exponential function, whereas the polynomial function W=


373.14

c3

α [°]

245.05

55.75


Kỷ yếu hội nghị khoa học và công nghệ toàn quốc về cơ khí - Lần thứ IV
Simulation results of the uniaxial tension test using the material constants in Table 1 with
the shear modulus µ * = 22.72 [kPa] are shown in Figure 3, the unphysical behavior occurs at
very low stretch, leading to an inaccurate Cauchy stress σ z curve of the MFH in Figure 2.
Three values of the shear modulus (remaining constants were unchanged) were
employed in the simulation, illustrated in Figure 3; the red one is for µ = 0.1µ * , the blue one
is simulated with the shear modulus µ = µ * and the last one (violet) uses µ = 10 µ * . It is
obvious that the largest shear modulus accounts for the best radial stretch which gradually
decreases due to the incompressibility constraint=
J det
=
F 1 . Subsequently, this value
makes the model more stable and the result is, of course, closer to the analytical. The
numerical stability of the model strongly depends on the isotropic strain-energy function (neoHookean one). The unphysical behavior (the radial stretch increasing) is clearly exhibited as
the shear modulus decreases.
This problem is observed as the constitutive matrix is ill-conditioned and is caused by
the large difference between the isotropic energy and the anisotropic energy [6]. Thus, we can
control properly the balance between the isotropic and the anisotropic strain energy functions.

If the exponent Qcoup is convex then the proposed model is logarithmically convex [9].
The exponent Qcoup is a quadratic function of E therefore Acoup is the Hessian matrix of

Qcoup . Thus, if Acoup is positive definite then the exponent Qcoup is convex. Specifically,
Acoup is positive definite when constraints on 8 material constants are expressed as
0 ≤ ν 1 ≤ 1, ( µ , C , c1 , c2 , c3 , c4 , c5 ) > 0 , (1 −ν 1 ){(c12 − c22 )(1 −ν 1 ) + c1c3ν 1} > 0 and
c1 (1 −ν 1 )

c2 (1 −ν 1 )

c2 (1 −ν 1 )

c2 (1 −ν 1 ) c1 (1 −ν 1 ) + c3ν 1 c2 (1 −ν 1 ) + c5ν 1 > 0.
c2 (1 −ν 1 ) c2 (1 −ν 1 ) + c5ν 1 c1 (1 −ν 1 ) + c4ν 1
894

(11)


Kỷ yếu hội nghị khoa học và công nghệ toàn quốc về cơ khí - Lần thứ IV
The constraints in (11) were checked after the material constants of the model were
obtained from the curve fitting process. Only the satisfying parameters were utilized in the
numerical example section. The material coefficients of the model become more meaningful
compared to those of the general Fung-type strain-energy function [10]. Moreover, the
weighting coefficient ν 1 , which might be considered as a fibre dispersion factor, plays an
important role to balance the difference between the isotropic strain energy and the
anisotropic one. If ν 1 = 0 then the whole model becomes fully isotropic. In contrast, ν 1 = 1 ,
the term Wcoup totally turns into an anisotropic one, and the overall model becomes the MFH.
This material law can be comparable to others taking the fibre dispersion into account [5].
3.2

1
2

(FT F - I) =

1
2

(C - I) .

(13)

Employing the Voigt’s notation, the modified Green-Lagrange strain is expressed in the
matrix form as E = [ E11 , E 22 , E 33 , 2 E12 , 2 E 23 , 2 E 31 ]T . The novel Fung-type strain-energy
function in a form of nearly incompressibility is

=
Wcoup C[exp( E T ⋅ A ⋅ E ) − 1].

(14)

The second Piola-Kirchhoff stress vector, which is calculated from Wcoup is
=
Scoup

∂Wcoup ∂E ∂Wcoup
∂E
( AE ),
=
= 2CeQ

∂ 2Wcoup ∂E
∂Q ∂E
∂ ∂E
[2CeQ ( A + AE ⊗
)
( )].
=
+ 2CeQ ( AE ) ⊗
2
∂E
∂E
∂E ∂E
∂E ∂E

Derivatives in (17) are calculated as
∂E ∂C
∂C
∂C
∂C
∂E
J 2/3 I=
2
2 J 2/3 I 6×6 ;
=
= =
=
=
6×6 ;
−2/3
−2/3

895

(17)


Kỷ yếu hội nghị khoa học và công nghệ toàn quốc về cơ khí - Lần thứ IV
∂ ∂E
∂E ∂J −2/3 1 −2/3 ∂C −1 ∂C −1
1
− J [
C +
C] .
( )=( I 6×6 − C ⊗ C −1 )2
∂E ∂E
∂E ∂C
∂E
∂E
3
3

4. NUNERICAL RESULTS
The MF2 can solve the unphysical response of the MFH. By utilizing simple tension
tests we can investigate the performance of the novel model in detail. In these tests, Table 2
show the material constants for the MF2 and the MFH that were obtained from the curve
fitting to the adventitial data [8], see Figure 1 and Figure 4. To prevent the ill-conditioning
problem of the stiffness matrix at high values of Poisson's ratio (0.4996 < ν < 0.5) , the Uzawa
algorithm for an augmented Lagrangian method is used for all numerical tests [11]. The MF2
is also capable of representing compressive state as shown in Figure 5, in which the neoHookean, the MFH and the Holzapfel models cannot describe the stress state.

Figure 4. Curve fitting for MF2 in uniaxial tests (left) and compression tests (right)


0.06

200.00

400.00

225.35

546.44

c5

ν1

-

-

200.00 0.75

In the uniaxial tension test, the numerical results employing the MF2 are more accurate
than those of the MFH, see Figure 5. The stretch and Cauchy stress computed for the MF2 are
not influenced by the ill-conditioning problem and are in good agreement with the analytical
ones, whereas the response of the MFH are severely affected by the ill-conditioned
constitutive matrix. Moreover, the numerical results for equi-biaxial tests of the MF2 are also
stable and are in good agreement with the analytical, see also [6].
896



[9] Conway, J. B., Functions of One Complex Variable. 2nd ed. New York: Springer Verlag, 2001.
[10] Holzapfel, G. A., Nonlinear solid mechanics, A Continuum Approach for Engineering:
Chichester John Wiley & Sons, 2000.
[11] Taylor, R. L., A Finite
www.ce.berkeley.edu/projects/feap/

Element

Analysis

Program,

2011,

from

[12] Fung, Y. C., Fronek, K., & Patitucci, P., Pseudoelasticity of arteries and the choice of its
mathematical expression, Am J Physiol., 1979. vol. 237, no. 5, p. H620–31.
AUTHOR’S INFORMATION
1.

TS. Minh Tuan DUONG, Hanoi University of Science and Technology, Hanoi, Vietnam
Tel: +84947036686, +84438680073

2.

TS. Thi Thanh Hai TRAN, Hanoi University of Science and Technology, Hanoi, Vietnam
Tel: +84978263926, +84438680073

897