Hiperelasticidad
Hiperelasticidad
Hiperelasticidad
Facultad de ingeniería
Departamento de Obras civiles
Project 3
Hyperelasticity
Abstract
In this project, we modeled the behavior of a hyperelastic material using both the ANSYS
program and Matlab. The material constants were derived through Curve Fitting based on
available experimental data. Our observations revealed that, although ANSYS produced a curve
that resembled the expected behavior, the Matlab-generated code demonstrated a superior ap-
proximation to the reference data. Notably, the model based on the Mooney-Rivlin formulation
exhibited better performance compared to the Neo-Hookean model.
1. Introduction
The incompressible Mooney-Rivlin hyperelastic model is a classical example of an isotropic hy-
perelastic material model. Its constitutive behaviour is defined by the strain energy function Ψe
as:
µ1 µ2
Ψ (C) = [I1 (C) − 3] − [I2 (C) − 3] (1)
2 2
where I1 (C) and I2 (C) are respectively the first and second invariant, C is the right Cauchy-
Green strain tensor, µ1 and µ2 are material constants with the shear modulus µ = µ1 −µ2 . Departing
from equation (1), the incompressible neo-Hookean hyperelastic model can be obtained from the
Mooney-Rivlin energy function by taking µ = µ1 and µ2 = 0. This leads to:
µ1
Ψe (C) =[I1 (C) − 3] (2)
2
The energy functions presented in (1) and (2) can be extended to the compressible range by
adding and additional function U (J) that depends on the Jacobian of the deformation gradient J,
yielding:
µ1 µ2
Ψe (C) = [I1 (C) − 3] − [I2 (C) − 3] + U (J) (3)
2 2
for the coupled compressible Mooney-Rivlin hyperelastic model, and
µ2
Ψe (C) = [I1 (C) − 3] + U (J) (4)
2
for the coupled compressible neo-Hookean hyperelastic model. The function U (J) can be defined
as:
1
U (J) = λ (J − 1)2 − µln (J), (5) (5)
2
or alternatively as
1
U (J) = λ ln2 (J) − µ ln (J) (6)
2
where λ and µ are the Lammé constants.
An alternative representation of the strain energy function Ψe (C) for compressible hyperelastic
materials can be formulated starting from the multiplicative split of the deformation process. This
procedure involves a multiplicative decomposition of the deformation gradient tensor F into a volume-
preserving isochoric or distortional part and a volume-changing dilatational part. This multiplicative
split yields:
2. Theorical study
2.1. Stress-Strain relation
In this section, expressions for the second Piola-Kirchhoff stress tensor S will be derived in terms
of the right Cauchy-Green strain tensor C for both the coupled and decoupled compressible Mooney-
Rivlin and Neo-Hookean hyperelastic material models. Additionally, an expression for the Cauchy
stress tensor σ in terms of the left Cauchy-Green strain tensor b will be obtained for the coupled
compressible Mooney-Rivlin and Neo-Hookean hyperelastic material models. These derivations con-
sider the functions (3) and (4) for the coupled compressible Mooney-Rivlin and Neo-Hookean models,
taking U (J) into account as defined in equations (6) and (10).
∂Ψ
S=2· (11)
∂C
Where, in this case, the strain energy function for the Mooney-Rivlin model is given by:
µ1 µ2
Ψe = [IC − 3] − [IIC − 3] + U (J) (12)
2 2
And the strain energy function for the Neo-Hookean model is given by:
µ1
Ψe = [IC − 3] + U (J) (13)
2
2.1.1.1. M-R
With U (J) = 1
2 · λ · ln2 (J) − µ · ln(J), S is expressed as:
∂Ψ
S=2· = µ1 · I − µ2 · [tr(C) · I − C] + λ · ln(J) · C−1 − µ · C−1 (14)
∂C
1
With U (J) = 2 · K · (J − 1)2 , S is expressed as:
∂Ψ
S=2· = µ1 · I − µ2 · [tr(C) · I − C] + K · (J − 1) · C−1 (15)
∂C
2.1.1.2. N-H
With U (J) = 1
2 · λ · ln2 (J) − µ · ln(J), S is expressed as:
∂Ψ
S=2· = µ1 · I + λ · ln(J) · C−1 − µ · C−1 (16)
∂C
1
With U (J) = 2 · K · (J − 1)2 , S is expressed as:
∂Ψ
S=2· = µ1 · I + K · (J − 1) · C−1 (17)
∂C
τ (b) = F · S · FT (18)
and that σ = Jτ .
2.1.2.1. M-R
With U (J) = 1
2 · λ · ln2 (J) − µ · ln(J) τ is expressed as:
µ1 · b − µ2 · [tr(b) · b − b2 ] + λ · ln(J) − µ
σ= (20)
J
1
With U (J) = 2 · K · (J − 1)2 τ is expressed as:
µ1 · b − µ2 · [tr(b) · b − b2 ]
σ= + K · (J − 1) (22)
J
2.1.2.2. N-H
With U (J) = 1
2 · λ · ln2 (J) − µ · ln(J) τ is expressed as:
µ · b + λ · ln(J) − µ
σ= (24)
J
1
With U (J) = 2 · K · (J − 1)2 τ is expressed as:
τ (b) = µ1 · b + K · (J − 1) · J (25)
Then:
µ1 · b
σ= + K · (J − 1) (26)
J
2.1.3. Decoupled compressible Mooney-Rivlin and neo-Hookean in terms
of S and C
In this case, the strain energy function for the Mooney-Rivlin model is given by:
µ1 2 µ2 4
Ψe = · [J − 3 · IC − 3] − · [J − 3 · IIC − 3] + U (J) (27)
2 2
And the strain energy function for the Neo-Hookean model is given by:
µ 2
Ψe = · [J − 3 · IC − 3] + U (J) (28)
2
2.1.3.1. M-R
With U (J) = 1
2 · λ · ln2 (J) − µ · ln(J), S is expressed as:
" 2 #
−J − 3 · C−1 −2
− 32 4
S = µ1 · J · I + I C · ( ) −µ2 ·J − 3 · · C−1 · IIC + tr (C) · I − C +λ·ln(J)·C−1 −µ·C−1
3 3
(29)
With U (J) = 12 · K · (J − 1)2 , S is expressed as:
" 2 #
−J − 3 · C−1 −2
− 32 4
S = µ1 · J · I + IC · ( ) −µ2 ·J − 3 · · C−1 · IIC + tr (C) · I − C +K·(J − 1)·J·C−1
3 3
(30)
2.1.3.2. N-H
With U (J) = 1
2 · λ · ln2 (J) − µ · ln(J), S is expressed as:
" 2 #
− 23 −J − 3 · C−1
S = µ1 · J · I + IC · ( ) + λ · ln(J) · C−1 − µ · C−1 (31)
3
1
With U (J) = 2 · K · (J − 1)2 , S is expressed as:
" 2 #
− 23 −J − 3 · C−1
S = µ1 · J · I + IC · ( ) + K · (J − 1) · J · C−1 (32)
3
∂S J
4
C = 2· = 2· µ · C−14 Is C−1 + λ · J 2 C−1 ⊗ C−1 − J 2 C−14 Is C−1 − C−1 ⊗ C−1 + JC−14 Is C−1
∂C 2
(34)
3. Numerical study
In this section, the three-parameter Mooney-Rivlin model and the Neo-Hookean model were im-
plemented in Matlab to analyze the behavior of a hyperelastic material studied by Hidajat Sugihardjo,
Tavio, and Yudha Lesmana. Additionally, experiments were simulated in ANSYS to compare the
results.
First, a curve fitting was conducted with the experimental data to obtain the material constants
C as referenced in (1).
Results:
Figure 5: Deformation
Results:
Figure 8: Deformation
Results:
Results:
4. Conclusions
Satisfactorily, expressions for the second Piola-Kirchhoff stress tensor in terms of C were success-
fully obtained for both the coupled and decoupled cases, as well as an expression for σ in terms of b.
Additionally, the fourth-order material elasticity tensor was derived from the indicated expressions.
Both the Mooney-Rivlin and Neo-Hookean models exhibited a good fit to the reference experi-
mental data, whether implemented in Ansys or Matlab. The similarity between the curves generated
by Ansys and Matlab indicates a robust implementation of the method, given that identical algebraic
expressions were employed.
To enhance the ANSYS model, a mesh convergence study and optimization of stress extraction
from each test could be conducted.
5. References
• Juan Carlos Pina. (2022). “Lecture Notes Tensor Algebra” and “Constitutive Relations".
• Sugihardjo, H., Tavio, T., & Lesmana, Y. (2018). FE model of low grade rubber for modeling
housing’s low-cost rubber base isolators. Civil Engineering Journal, 4(1), 24. https://doi.org/10.28991/cej-
030966
• Holzapfel, G. A. (2000). Nonlinear solid mechanics: A continuum approach for engineering (1a
ed.). John Wiley & Sons.
1. Appendix
1.1. Useful equivalences
The invariants are defined as:
• IC = tr(C) = C : I
1
tr(C)2 − tr(C2 )
• IIC = 2
∂C = ∂C = 2 ∂C
1
2 2tr(C) ∂tr(C)
∂C −
∂C:C
∂C = 1
2 (2tr(C)I − 2C) = tr(C)I − C
3. ∂J
∂C = 12 JC−1
∂I
4. ∂C =0
∂C
5. ∂C =4 Is
∂C−1
6. ∂C = −C−1 ·4 Is · C−1
proof:
∂I ∂C−1 C
∂C = ∂C =0
∂C−1 −1 ∂C
∂C C + C ∂C = 0
∂C−1 −1 ∂C
∂C C = −C ∂C
∂C−1 −1 4 s
∂C = −C · I · C−1 = −C−1 ⊗ C−T
13 /PREP7
14 ! DATOS HIPERELASTICOS
15
16 ! Mat: MOONEY-RIVLIN
17 TB,HYPE,1,1,3,MOON
18 TBTEMP,0
19 TBDATA,,C10,C01,C11,De1,,
20
21 ! ELEMENTO
22 ET,1,SOLID185
23
24 !GEOMETRIA
25 ! CILINDRO COMPLETO
26 CYL4,0,0,0,0,25,360,2
27
28 ! MALLADO
29
38 MSHAPE,0,3D
39 MSHKEY,1
40 VMESH,ALL
41
42 ! COORDENADAS POLARES
43 CSYS,1 ! Cambiar el eje coordenado de cartesiana a polares
44 NROTAT,ALL ! Rota el eje de coordenadas de cada nodo
45
46 ! Solucion
47 /SOLU
48 ANTYPE,0 ! Analisis Estructural
49 NLGEOM,ON ! Se activa la no-linealidad geometrica
50
51 ! Restricciones
52 DA,2,UZ,0 ! Restriccion de movimiento en Z
53 DA,15,UZ,0 ! Restriccion de movimiento en Z
54 DA,21,UZ,0 ! Restriccion de movimiento en Z
55 DA,18,UZ,0 ! Restriccion de movimiento en Z
56
67 ! Tiempo de estudio
68 t0 = 0 ! Tiempo inicial
69 tf = 1 ! Tiempo final
70 dt = 0.01 ! Salto de tiempo
71 steps = (tf-t0)/dt + 1
72 NSUBST,steps,steps,steps ! Define el numero de pasos de analisis
73 OUTRES,ALL,ALL ! Solicitar los resultados
74
75 ALLSEL,ALL
76 SOLVE
14 /PREP7
15 ! DATOS HIPERELASTICOS
16
17 ! Mat: MOONEY-RIVLIN
18 TB,HYPE,1,1,3,MOON
19 TBTEMP,0
20 TBDATA,,C10,C01,C11,De1,,
21
22 ! ELEMENTO
23 ET,1,SOLID185
24
25 !GEOMETRIA
26 ! Block
27 BLOCK,0,50,0,3,0,2,
28
29 ! MALLADO
30 VATT,1,,1,,
31 MSHKEY,1
32 MSHAPE,0,3D
33 ESIZE,1 ! Tamaño de elementos
34 ALLSEL,ALL
35 VMESH,ALL
36
37
38 ! Solucion
39 /SOLU
40 ANTYPE,0 ! Analisis Estructural
41 NLGEOM,ON ! Se activa la no-linealidad geometrica
42
43 !!RESTRICCIONES
44
56 ASEL,S,LOC,Z,0!Fijar caras en Z
57 NSLA,S,1
58 D,ALL,UZ,0 !Restriccion
59 allsel,all
60
61 !Desplazamiento
62 asel,S,LOC,X,50
63 NSLA,S,1
64 D,all,UX,50
65 ALLSELL,ALL
66
67 ! Tiempo de estudio
68 t0 = 0 ! Tiempo inicial
69 tf = 1 ! Tiempo final
70 dt = 0.01 ! Salto de tiempo
71 steps = (tf-t0)/dt + 1
72 NSUBST,steps,steps,steps ! Define el numero de pasos de analisis
73 OUTRES,ALL,ALL ! Solicitar los resultados
74
75 ALLSEL,ALL
76 SOLVE
13 /PREP7
14 ! DATOS HIPERELASTICOS
15
16 ! Mat: MOONEY-RIVLIN
17 TB,HYPE,1,1,3,MOON
18 TBTEMP,0
19 TBDATA,,C10,C01,C11,De1,,
20
21 ! ELEMENTO
22 ET,1,SOLID185
23
24 !GEOMETRIA
25 ! Block
26 BLOCK,0,150,0,15,0,2,
27
28 ! MALLADO
29 VATT,1,,1,,
30 MSHKEY,1
31 MSHAPE,0,3D
32 ESIZE,5,0, ! Tamaño de elementos
33 ALLSEL,ALL
34 VMESH,ALL
35
36
37 ! Solucion
38 /SOLU
39 ANTYPE,0 ! Analisis Estructural
40 NLGEOM,ON ! Se activa la no-linealidad geometrica
41
42 ! Restricciones
43
44 ASEL,S,LOC,Y,0
45 NSLA,S,1
46 D,ALL,UY,0
47 ALLSEL,ALL
48
49 ASEL,S,LOC,Y,0
50 NSLA,S,1
51 D,ALL,UX,0
52 ALLSEL,ALL
53
54 ASEL,S,LOC,Y,0
55 NSLA,S,1
56 D,ALL,UZ,0
57 ALLSEL,ALL
58
59 ASEL,S,LOC,Y,15
60 NSLA,S,1
61 D,ALL,UX,0
62 ALLSEL,ALL
63
64 ASEL,S,LOC,Y,15
65 NSLA,S,1
66 D,ALL,UZ,0
67 ALLSEL,ALL
68
69 !Desplazamiento
70 asel,S,LOC,Y,15
71 NSLA,S,1
72 D,ALL,UY,15
73 ALLSEL,ALL
74
75 ! Tiempo de estudio
76 t0 = 0 ! Tiempo inicial
77 tf = 1 ! Tiempo final
78 dt = 0.01 ! Salto de tiempo
79 steps = (tf-t0)/dt + 1
80 NSUBST,steps,steps,steps ! Define el numero de pasos de analisis
81 OUTRES,ALL,ALL ! Solicitar los resultados
82
83 ALLSEL,ALL
84 SOLVE
13 /PREP7
14 ! DATOS HIPERELASTICOS
15
16 ! Mat: MOONEY-RIVLIN
17 TB,HYPE,1,1,3,MOON
18 TBTEMP,0
19 TBDATA,,C10,C01,C11,De1,,
20
21 ! ELEMENTO
22 ET,1,SOLID185
23
24 !GEOMETRIA
25 ! CILINDRO COMPLETO
26 CYL4,0,0,0,0,3.175,360,10
27
28 ! MALLADO
29
42 ! COORDENADAS POLARES
43 CSYS,1 ! Cambiar el eje coordenado de cartesiana a polares
44 NROTAT,ALL ! Rota el eje de coordenadas de cada nodo
45
46 ! Solucion
47 /SOLU
48 ANTYPE,0 ! Analisis Estructural
49 NLGEOM,ON ! Se activa la no-linealidad geometrica
50
51 ! Restricciones
52 DA,2,UZ,-2 ! Restriccion de movimiento en Z
53 DA,15,UZ,-2 ! Restriccion de movimiento en Z
54 DA,21,UZ,-2 ! Restriccion de movimiento en Z
55 DA,18,UZ,-2 ! Restriccion de movimiento en Z
56
72 ! Tiempo de estudio
73 t0 = 0 ! Tiempo inicial
74 tf = 1 ! Tiempo final
75 dt = 0.01 ! Salto de tiempo
76 steps = (tf-t0)/dt + 1
77 NSUBST,steps,steps,steps ! Define el numero de pasos de analisis
78 OUTRES,ALL,ALL ! Solicitar los resultados
79
80 ALLSEL,ALL
81 SOLVE
For Neo-Hookean, the same codes are used, replacing the following:
1 ! Parametros para Neo-Hookean
2 mu = 0.438899369 ! en MPa
3 De1 = 0.0006057 ! en MPa^-1
4
5 /PREP7
6 ! DATOS HIPERELASTICOS
7
8 ! Mat: Neo-Hookean
9 TB,HYPE,1,1,2,NEO
10 TBTEMP,0
11 TBDATA,,mu,De1,,,,
2 close all;
3 clear;
4
5 fclose(’all’);
6 %% Starting Values
7 unit_tensors; % Initialize unit tensors.
8 nd =3; % Problem dimensions.
9
10 %% Material Properties
11 K = 2/0.0006057 ; % Bulk modulus [MPa^-1]
12 c10 = 0.182353216; % 1st M-R constant. [MPa]
13 c01 = 0.034139345; % 2nd M-R constant. [MPa]
14 c11 = -0.001180849; % 3rd M-R constant. [MPa]
15
16 mat=1;
17 %% Method
18 for ensayo=1:4 %
19 gama = 0; %Amplifier
20 nincr = 51; %N° of Increments
21
22 for incr=1:nincr
23 %% Deformation Gradient Tensor
24 if incr == 1
25 gama = 0;
26 else
27 gama = gama + 1/(nincr-1); % Monotonically increasing loading factor
28 end
29
30 if ensayo == 1
31 gradu0_vol = [1-gama*0.08 0 0; 0 1 0; 0 0 1]; %Vol def - Vol test
32 elseif ensayo == 2
33 gradu0_vol = [1+gama*1 0 0; 0 1/sqrt(1+gama*1) 0; 0 0 1/sqrt(1+gama*1)]; %Vol def -
,→ uniaxial test
34 elseif ensayo == 3
35 gradu0_vol = [1+gama*1 0 0; 0 (1+gama*1)^-1 0; 0 0 1]; %Vol def - pla shear
,→ test
36 else
37 gradu0_vol = [1+gama*1 0 0; 0 1+gama*1 0; 0 0 (1+gama*1)^-2]; %Vol def - pla
,→ shear test
38 end
39
66
77 % Hencky strain
78 [eigVecU,eigValU] = eig(U,’vector’);
79 eigValU = log(eigValU);
80 epsilon = zeros(nd,nd);
81 for i = 1:3
82 epsilon = eigValU(i)* eigVecU(:,i)*eigVecU(:,i)’ + epsilon;
83 end
84
89
90 %% VOLUMETRI TEST
91 if ensayo == 1
92 VOLUMETRIC_s(incr,mat)=2/(2/K)*(J-1);
93 VOLUMETRIC_e(incr,mat) = F(1,1);
94
95 %% UNIAXIAL TEST
96 elseif ensayo == 2
97 true_stress_u=2*(l^2-l^-1)*(c10 + l^-1*c01 + c11*((I2C-3)+l^-1*(I1C-3)));
98 eng_stress_u = true_stress_u/l;
99 UNIAXIAL_s(incr,mat) = eng_stress_u;
100 UNIAXIAL_e(incr,mat) = F(1,1)-1;
101
129
139
To adapt the code for the Neo-Hookean model, the following modification is made in the code:
1 K = 2/0.0006057 ;
2 mu = 0.438899369 ;
3 c10 = mu/2;
4 c01 = 0;
5 c11 = 0;