- Email: [email protected]
Material-symmetries congruency in transversely isotropic and orthotropic hyperelastic materials
Accepted Manuscript Material-symmetries congruency in transversely isotropic and orthotropic hyperelastic materials Marcos Latorre, Francisco Javier Montáns PII:
S0997-7538(15)00027-3
DOI:
10.1016/j.euromechsol.2015.03.007
Reference:
EJMSOL 3175
To appear in:
European Journal of Mechanics / A Solids
Received Date: 29 August 2014 Revised Date:
22 January 2015
Accepted Date: 23 March 2015
Please cite this article as: Latorre, M., Montáns, F.J., Material-symmetries congruency in transversely isotropic and orthotropic hyperelastic materials, European Journal of Mechanics / A Solids (2015), doi: 10.1016/j.euromechsol.2015.03.007. This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT
Marcos Latorre, Francisco Javier Mont´ans∗
RI PT
Material-symmetries congruency in transversely isotropic and orthotropic hyperelastic materials
SC
Escuela T´ecnica Superior de Ingenier´ıa Aeron´ autica y del Espacio Universidad Polit´ecnica de Madrid Plaza Cardenal Cisneros, 3, 28040-Madrid, Spain
Abstract
M AN U
In this work we present the material-symmetries congruency property for anisotropic hyperelastic materials. A transversely isotropic or orthotropic material model holding this property guarantees the recovery of the isotropic behavior when, for example, the volume of reinforcement in a composite material vanishes. A material presenting material-symmetries congruency
TE D
must show it not only from an analytical, theoretical point of view, but also from a numerical, practical point of view. The former may be obtained from construction of the stored energy function, whereas the latter must be obtained guaranteeing that the parameter-fitting procedure yields material
EP
parameters which combined with the specific model results in an isotropic behavior prediction for isotropic materials. We show that some anisotropic
AC C
models do not present either analytical or numerical material-symmetries congruency.
Keywords: Hyperelasticity, anisotropy, incompressible materials. ∗
Corresponding author. Tel.:+34 637 908 304. Email addresses: [email protected] (Marcos Latorre), [email protected] (Francisco Javier Mont´ans) Preprint submitted to European Journal of Mechanics - A/Solids
January 22, 2015
ACCEPTED MANUSCRIPT
1
1. Introduction Many material models have been developed for analyzing hyperelastic
3
transversely isotropic materials and for orthotropic materials, see for ex-
4
ample [1],[2],[3],[4],[5],[6],[7],[8],[9],[10],[11],[12],[13],[14],[15],[16],[18],[19],[23],
5
among others. The main objective of those material models is to account for
6
different material behavior encountered in different material directions for
7
some materials like rubber [1],[14],[17], fibre-reinforced composites [15],[16]
8
and biological materials [4],[5],[6],[7],[8],[9],[10],[11],[12],[13],[18],[19].
M AN U
SC
RI PT
2
The different material properties in the different directions are frequently
10
due to a reinforcing constituent which has some preference in alignment.
11
This reinforcement may be laid-out during manufacturing or, as in biological
12
tissues, by nature. A common approach is to test the composite material as
13
a whole, see for example [12], [17], [23], [24], among others.
TE D
9
Hyperelastic material models typically consist on a stored energy function
15
of some deformation-based invariants [3], some of which are derived from
16
structural tensors obtained from the preferred material directions. These
17
energy functions have a pre-defined shape and some material parameters
18
that are obtained fitting experimental values from some tests using a proper
19
optimization algorithm, typically the Levenberg-Marquardt algorithm, see
20
for example [1], [25] and [26]. Some proposed stored energy functions take
21
into account the fact that anisotropy is mainly due to reinforcement and they
22
separate the stored energy into an isotropic (bulk) part and a reinforcement
23
part, see for example [3], [15], [18]. In fact, these models may include specific
24
parameters that specifically take into account the amount of anisotropy or
AC C
EP
14
2
ACCEPTED MANUSCRIPT
25
reinforcement [3], [15], [16]. Other proposals are more phenomenological and
26
consider from the onset the material as a whole, see for example [1],[2],[19]. Furthermore, in anisotropy it is also frequent to have available less exper-
28
imental data than those necessary to fully characterize the material [31] and
29
then the missing information may strongly influence the overall behavior [2].
30
Many of these available constitutive models for transversely isotropic and
31
orthotropic materials intend to represent the behavior of different materi-
32
als with different degrees of reinforcement, and they intend to do so with
33
parameters obtained from some tests over the composite. One would ex-
34
pect that when the reinforcement volume goes to zero in the composite, the
35
material captures the behavior of the isotropic matrix in an isotropic way.
36
Furthermore, in a orthotropic material, if the reinforcement in one direction
37
vanishes, one would also expect that the model captures the behavior of the
38
transversely isotropic material in a transversely isotropic manner. We say of
39
a model behaving this way that the model has material-symmetries congru-
40
ency because it converges to an isotropic material when the data corresponds
41
to that of an isotropic material.
TE D
M AN U
SC
RI PT
27
The material-symmetries congruency is not only an analytical issue on
43
the form of the stored energy function. A hyperelastic constitutive model
44
may present congruency from an analytical standpoint by construction, but
45
still lack that congruency from a numerical point of view, i.e. the model
46
using the parameters obtained from (quasi-)isotropic experimental curves
47
may not behave in a (quasi-)isotropic manner, simply because the parameter
48
fitting algorithm does not distinguish between quasi isotropic, transversely
49
isotropic and orthotropic materials. Many possible minima can be obtained
AC C
EP
42
3
ACCEPTED MANUSCRIPT
during the parameter fitting procedure which strongly influence the predicted
51
behavior even in the isotropic case [26]. In this case, if a model with analyt-
52
ical material-symmetries congruency loses this congruency due to the actual
53
parameters identified from the experiments, we say that whereas the model
54
has analytical material-symmetries congruency, it lacks the numerical one.
55
Material-symmetries congruency is specially important when some experi-
56
mental data are not available, and also in models that use decompositions of
57
the type of Valanis-Landel [33] because this decomposition has been verified
58
for isotropic materials both experimentally and analytically up to high orders
59
of strain [22], but in general remains to be an assumption on the form of the
60
stored energy for anisotropic materials. If material-symmetries congruency is
61
preserved, only the decomposition on the deviation from isotropy remains as
62
an assumption for moderate large strains. Furthermore, this deviation will
63
tend to vanish in the limit of isotropic behavior.
M AN U
SC
RI PT
50
The purpose of this paper is to present the material-symmetries congru-
65
ency for anisotropic hyperelastic materials and to analyze some models as
66
examples, both from an analytical perspective and from a numerical one.
67
2. Theoretical and numerical material-symmetries congruency
EP
Let us define the space L of the proper orthogonal tensors Q
AC C
68
TE D
64
Q ∈ L / QQT = I
4
and
det(Q) = 1
(1)
ACCEPTED MANUSCRIPT
69
and the symmetry group for orthotropic materials Lor ⊂ L containing the
70
set of rotation tensors Qor ∈ Lor = {Qπa0 , Qπb }
RI PT
(2)
0
71
where Qπa0 = a0 ⊗a0 −b0 ⊗b0 −c0 ⊗c0 and Qπb = −a0 ⊗a0 +b0 ⊗b0 −c0 ⊗c0 0
represent two rotations of π radians around two preferred material directions,
73
defined by the orthogonal vectors a0 and b0 . Note that the tensor Qπc0 =
74
Qπa0 Qπb associated to the third preferred direction, c0 = a0 × b0 , is also 0
75
included in Lor .
M AN U
76
SC
72
For an anisotropic hyperelastic material, in general ( ) W (E) ̸= W QEQT
(3)
for a given Q and a Lagrangian strain measure E. If the material is or-
78
thotropic, it has two (subsequently, three) planes of material symmetry, with
79
a0 and b0 defining these two planes and being perpendicular to them. Then,
80
it can be stated that
TE D
77
EP
( ) W (E, a0 , b0 ) = W Qor EQTor , a0 , b0
(4)
where Qor stands now for any arbitrary rotation included in the symmetry
82
group Lor . If Q ∈ / Lor , then W(E, a0 , b0 ) ̸= W(QEQT , a0 , b0 ) in general. If
83
completely arbitrary rotations Q are considered, this more general invariance
AC C
81
5
ACCEPTED MANUSCRIPT
84
relation (i.e. frame independence) is to be fulfilled ( ) W (E, a0 , b0 ) = W QEQT , Qa0 , Qb0
RI PT
(5)
that is, W(E, a0 , b0 ) must be an isotropic function with respect to its three
86
arguments simultaneously. Note that Eq. (4) stands for symmetry consider-
87
ations of the strain energy function while Eq. (5) stands for frame indepen-
88
dence so W(E, a0 , b0 ) must satisfy both of them. A material model satisfying
89
only Eq. (5) and not Eq. (4) has not necessarily orthotropic symmetries and
90
it would not be an acceptable formulation for an orthotropic material.
M AN U
SC
85
Another property that the function W(E, a0 , b0 ) must obey and which
92
seem to have never been suggested in the literature (to the best of the author’s
93
knowledge) is what we call material-symmetries congruency. This require-
94
ment states that for an orthotropic hyperelastic material, the model being
95
employed should tend to an isotropic formulation when slightly orthotropic,
96
nearly isotropic, materials are being characterized. That is, the tendency of
97
the orthotropic strain energy function when the limit of isotropic behavior is
98
considered should be such that
TE D
91
EP
( ) W (E, a0 , b0 ) ≈ W QEQT , a0 , b0
100
for any proper orthogonal tensor Q. Moreover, if an isotropic material is
AC C
99
(6)
characterized using an orthotropic model, then it is required that ) ( W (E, a0 , b0 ) = W QEQT , a0 , b0
6
(7)
ACCEPTED MANUSCRIPT
for any Q, as it effectively occurs for isotropic materials, i.e. W(E) =
102
W(QEQT ). Note the difference between Eq. (4), only forced to be ful-
103
filled for Qor ∈ Lor , and Eqs. (6) and (7), enforced for arbitrary rotations
104
Q ∈ L.
RI PT
101
We find in the literature some formulations satisfying the requirements
106
given in Eqs. (4) and (5), but which do not fulfill Eqs. (6) and (7) either
107
from a theoretical or from a, practical, numerical point of view. For example,
108
consider the incompressible uncoupled orthotropic model
SC
105
M AN U
W (E, a0 , b0 ) = ω11 (E11 )+ω22 (E22 )+ω33 (E33 )+2ω12 (E12 )+2ω23 (E23 )+2ω31 (E31 ) (8)
that we present in Ref. [2]. In Eq. (8), the arguments of the ω functions
110
are the components of the (deviatoric) Hencky strain tensor in the material
111
preferred basis X pr = {e1 , e2 , e3 } = {a0 , b0 , c0 }, that is Eij = ei · Eej , with
112
i, j = {1, 2, 3}. Note that these logarithmic strain components are deviatoric
113
(traceless) due to the incompressibility condition, and that the product of the
114
principal stretches is one because of the same reason. Thus, the invariance
115
relation given in Eq. (5) is automatically satisfied and the model is invariant
116
under simultaneous rotations of E, a0 and b0 . However, note that if any
117
of the shear terms ωij , i ̸= j, is not an even function of its corresponding
118
argument Eij , then the model does not satisfy the symmetry requirement of
119
Eq. (4) and the model would not be physically correct. Therefore, we only
120
consider even functions for the shear terms ωij as we properly address in
121
Ref. [2]. This fact shows that both symmetry and invariance requirements
122
must be satisfied for any orthotropic model and that Eq. (5) is not a suffi-
AC C
EP
TE D
109
7
ACCEPTED MANUSCRIPT
cient condition for a material model to reproduce an orthotropic behavior,
124
in contrast to what one could erroneously infer from, for example, Ref. [3]
125
(cf. page 274).
RI PT
123
It can be readily shown that if the orthotropic model of Eq. (8) is used
127
to characterize an assumed symmetric-Valanis-Landel-type incompressible
128
isotropic material from the corresponding set of six different experimental
129
data curves, then the computed model using the procedure detailed in Ref.
130
[2] reduces to
SC
126
M AN U
W (E, a0 , b0 ) = ω (E11 ) + ω (E22 ) + ω (E33 ) + 2ω (E12 ) + 2ω (E23 ) + 2ω (E31 ) (9)
with ω (E) being an even function of the strain component E. If as a result of
132
the model determination procedure the function ω contains quadratic terms
133
on its argument, it can be shown that the complete invariant E : E is
134
recovered in W multiplied by the proper coefficient. On the contrary, if as a
135
result of the model determination procedure the function ω contains fourth-
136
order terms on its argument, it can be shown that the complete invariant
137
E 4 : I is not recovered in W multiplied by the proper coefficient, so Eq. (4)
138
does not hold. Hence, although an arbitrary strain energy for an orthotropic
139
material can be approximated as closely as desired by a suitable polynomial of
140
the corresponding symmetry-preserving invariants (cf. Ref. [22], pages 211–
141
212), we remark that this “suitable” polynomial should contain the isotropic
142
formulation for the corresponding limit.
AC C
EP
TE D
131
143
Regarding the previous formulation, a possibility to modify the model to
144
obtain material-symmetries congruency simply consists in adding an isotropic
8
ACCEPTED MANUSCRIPT
145
contribution to the orthotropic one W (E, a0 , b0 ) = Wis (E) + Wor (E, a0 , b0 )
RI PT
(10)
146
where Wis (E) is the Valanis-Landel-type isotropic model expressed in terms
147
of principal logarithmic strains defined by Sussman and Bathe in Ref. [20],
148
i.e.
SC
Wis (E) = ω (E1 ) + ω (E2 ) + ω (E3 )
(11)
and Wor (E, a0 , b0 ) is given by Eq. (8). With respect to the modified model
150
provided in Eq. (10), the material-symmetries congruency requirements Eqs.
151
(6) and (7) will be effectively satisfied if both Wor (E, a0 , b0 ) tends to vanish
152
when slightly orthotropic, nearly isotropic, materials are characterized and
153
Wor (E, a0 , b0 ) = 0 in the limit of pure isotropic behavior, respectively. In
154
this case, due to the use of the spline-based methodology and the inversion
155
formula (cf. Refs. [19], [2], [20]), the equilibrium equations of the tests
156
are solved exactly for the isotropic model. This important fact implies that
157
Wis (E) exactly reproduces the prescribed isotropic behavior and that, as a
158
consequence, Wor (E, a0 , b0 ) vanishes for these type of materials.
159
3. Examples
160
3.1. Spline-based model according to Latorre and Mont´ans [19]
AC C
EP
TE D
M AN U
149
161
In Ref. [19] we present a spline-based transversely isotropic model which
162
is capable of reproducing the uniaxial experimental data of Diani et al. [23]
163
over transversely isotropic calendered rubber with a very good (exact in
164
practice) agreement. In that work, we also compare the results obtained with 9
ACCEPTED MANUSCRIPT
our model to those predicted by the model of Itskov and Aksel [1], showing a
166
better accuracy achieved by the spline model. Following similar reasonings to
167
those detailed above, it is readily shown that the spline-based model that we
168
presented in Ref. [19] does not have isotropic material-symmetry congruency.
169
However, as mentioned above, the addition of the isotropic contribution of
170
Eq. (11) circumvents this issue.
RI PT
165
In this example, and the next ones, we study two limit cases in which the
172
two curves presented by Diani et al coalesce to one of both curves. The aim
173
of these examples is to predict with transversely isotropic models two possi-
174
ble isotropic behaviors obtained as isotropic limits of a transversely isotropic
175
material and analyze the numerical material-symmetries congruency of the
176
models. Since the uniaxial tests are performed on an initially assumed trans-
177
versely isotropic material and along its “preferred” material directions, no
178
shear terms are needed and the extended spline-based transversely isotropic
179
model of Ref. [19] reduces to
TE D
M AN U
SC
171
W (E, a0 ) = Wis (E) + Wtr (E, a0 )
EP
= ω1# (E1 ) + ω1# (E2 ) + ω3# (E3 )
(12) (13)
where the simplifications ωn# (E) = ω (E) + ωn (E), n = 1, 3, are considered
181
for practical purposes. As in Ref. [19], the subscripts 1 and 2 are associated
182
to two orthogonal directions within the isotropic plane and the subscript 3
183
represents the assumed preferred direction of the material.
184
AC C
180
For the first case, we prescribe two discrete stress-stretch distributions
10
M AN U
SC
RI PT
ACCEPTED MANUSCRIPT
Figure 1: Nominal stresses P3 and P1 in terms of the stretches λ3 and λ1 , respectively, from both uniaxial tests (black solid dots) performed on a transversely isotropic material in a hypothetical limit of isotropic behavior and predictions obtained with the spline-based model. Original experimental data is obtained from the test in calendering direction from Diani et al [23].
which are identical to the original measured data points in the uniaxial test
186
performed along the calendering direction, i.e. direction 3, as it is shown in
187
Figure 1 with black solid dots. Using these “experimental” distributions as
188
initial data of the procedure detailed in Table 3 of Ref. [19] (where the func-
189
tions ωn are simply substituted by ωn# ), the solution functions are found to be
190
ω1# (E) = ω3# (E), with the transverse strain distribution E2 (E1 ) = −E1 /2
191
(ν12 = 1/2) computed also as a part of the solution. Clearly, the solution
192
converges to the corresponding isotropic solution. That is, considering the
193
relations ωn# (E) = ω (E) + ωn (E), n = 1, 3, we deduce
AC C
EP
TE D
185
11
ACCEPTED MANUSCRIPT
ωn (E) = 0
(14)
RI PT
ωn# (E) = ω (E) ;
These results indicate that the model given in Eqs. (12)–(13) has also
195
isotropic material-symmetry congruency from a numerical point of view, as
196
one would expect. The stress predictions in both directions provided by
197
the computed model are also shown in Figure 1. A strictly exact fit of the
198
experimental data is attained, as it should occur within the isotropic context
199
(cf. Refs. [20] and [19]).
200
3.2. Model of Itskov and Aksel [1]
M AN U
SC
194
Regarding the Itskov-Aksel work of Ref. [1], the orthotropic model pro-
202
vided therein in terms of generalized invariants, namely Eq. (58) of that
203
reference, includes the isotropic generalized Mooney-Rivlin model, Eq. (68)
204
of the same reference, as a particular case. Hence, the proposed strain energy
205
function has theoretical material-symmetries congruency. However, when cer-
206
tain slightly transversely isotropic, nearly isotropic, materials are character-
207
ized with the parameter identification procedure detailed in their work, the
208
tendency to the Mooney-Rivlin formulation may not be observed. Further-
209
more, when two identical stress-strain distributions are used to initialize the
210
parameter identification process, i.e. pure isotropic behavior is prescribed,
211
several specific transversely isotropic solutions are found to be more optimum
212
(in a least-squares sense) than the isotropic solution. This is due to the fact
213
that the optimization process employed to identify the material parameters
214
does not ensure that their limit values for an isotropic material are obtained.
AC C
EP
TE D
201
12
ACCEPTED MANUSCRIPT
This fact manifests that the Itskov-Aksel formulation, as given in their pa-
216
per, does not have material-symmetries congruency from a numerical point
217
of view.
218
219
RI PT
215
In order to show this unexpected result, the transversely isotropic incompressible model by Itskov and Aksel, namely Eq. (95) of Ref. [1]
[( ( )αr ] )βr 3 3 1 ∑ ∑ 1 (r) (r) W= µr wi λ2i −1 + wi λ−2 − 1 i 4 r=1 αr βr i=1 i=1 (r)
SC
3 1∑
where µr , αr , βr and wi
221
relations in Eq. (98) of the same Reference, i.e. (r)
are material parameters, along with the symmetry
(r)
w2 = w3 =
M AN U
220
(15)
) 1( (r) 1 − w1 , r = 1, 2, 3 2
(16)
are employed to characterize the isotropic materials of the preceding example
223
using the minimization procedure detailed in Ref. [1]. A standard least
224
squares curve-fitting algorithm (specifically, the Matlab function lsqcurvefit)
225
has been used to fit the two point distributions shown in Figure 1 using Eqs.
226
(100)–(102) of Ref. [1]. The optimization procedure has been initialized
227
with the set of material parameters given in Ref. [1], i.e. Eq. (104), and
228
the principal stretches in the calendering direction λ1 (j = 1, ..., m) for the
229
uniaxial test in the transversely isotropic direction 2 obtained as solution of
230
their parameter identification procedure (note that in this example we are
231
using the index numbering of Ref. [1], which is different from the one used in
232
the preceding example). All these material parameters and stretches initially
233
correspond to a transversely isotropic behavior. In this first attempt to assess
EP
TE D
222
AC C
(j)
13
ACCEPTED MANUSCRIPT
234
the numerical material-symmetries congruency of the model, the values of
235
the parameters which have been found as “best” solution using these initial
236
values are given in Table 1. The fact that w1 ̸= 1/3, r = 1, 2, 3 indicates
237
that the computed solution has not converged to the associated isotropic
238
solution (cf. Eqs. (67) and (98) in Ref. [1]). In Figure 2, it can be observed
239
that the stresses predicted by the model defined with the parameters given in
240
Table 1 reproduce the tendency of the “experimental” data, although little
241
differences between both curves (due to the lack of isotropy) are observable. (1)
w1 = 0.900 (2) w1 (3) w1
α1 = 4.27 α2 = 6.84 α3 = 1.0
β1 = 9.44 β2 = 7.78 β3 = 1.59
M AN U
µ1 = 3.84 × 10−3 MPa µ2 = 2.03 × 10−3 MPa µ3 = 5.44 MPa
SC
RI PT
(r)
= 0.0
= 0.323
Table 1: Computed material parameters for the transversely isotropic model of Ref. [1] corresponding to the first approach in Figure 2. Initial guess (associated to a transversely isotropic material) obtained from Ref. [1].
Clearly, the foregoing solution for the set of material parameters being
243
investigated corresponds to the nearest local minimum to the initial values
244
given in Eq. (104) of Ref. [1]. Obviously, a local solution of any optimization
245
procedure is strongly affected by the initial numerical values given to the
246
unknowns of the optimization problem. Hence, one may arguably think that
247
the isotropic limit solution could be obtained if another, more appropriate, set
248
of initial values nearer to the isotropic solution is taken. In a second attempt
249
to assess the numerical material-symmetries congruency of the model, we
250
have previously obtained the values of the nine material parameters µr , αr
251
and βr , r = 1, 2, 3, of the isotropic model (with w1 = 1/3, r = 1, 2, 3) which
252
best fit one of the two identical experimental distributions of Figure 1. For
AC C
EP
TE D
242
(r)
14
TE D
M AN U
SC
RI PT
ACCEPTED MANUSCRIPT
AC C
EP
Figure 2: Nominal stresses P3 and P1 in terms of the stretches λ3 and λ1 , respectively, from both uniaxial tests (black solid dots) performed on a transversely isotropic material in a hypothetical limit of isotropic behavior and predictions obtained with the model of Itskov and Aksel, Ref [1], with two different sets of material parameters. Original experimental data is obtained from the test in calendering direction from Diani et al [23]. Index numbering is consistent with Figure 1.
15
ACCEPTED MANUSCRIPT
253
this fist step, the experimental data points have been fitted by means of Eq.
254
(100) of Ref. [1] (with w1 = 1/3). The values of the parameters which have
255
been found as “best” solution for the isotropic model are given in Table 2.
256
These values provide a very good agreement of the isotropic model with the
257
experimental data. α1 = 4.495 α2 = 4.338 α3 = 1.003
β1 = 1.004 β2 = 9.017 β3 = 1.0
SC
µ1 = 2.30 × 10−1 MPa µ2 = 2.46 × 10−3 MPa µ3 = 5.39 MPa
RI PT
(r)
M AN U
Table 2: Material parameters of the isotropic model of Ref. [1] computed by means of the fitting of the experimental data (only one curve) of Figure 2.
258
Subsequently, as before, we have proceeded to fit the two point distribu-
259
tions of Figure 1 using Eqs. (100)–(102) of Ref. [1]. In this second step,
260
the optimization procedure has been initialized with the computed material
261
parameters given in Table 2 together with the values w1 = 1/3, r = 1, 2, 3,
262
and the principal stretches λ1
263
describe an isotropic model. Again, a solution corresponding to a trans-
264
versely isotropic description of the model, i.e. with w1 ̸= 1/3, r = 1, 2, 3, is
265
found to be the “best” solution for the considered initial values. In Figure
266
2, no noticeable differences are observable between the stress distributions
267
predicted by the model in both directions. However, this is just an apparent
268
isotropic response of the model in a specific situation (which has been inten-
269
tionally fitted). We want to emphasize that in other more generic situations
270
an isotropic mechanical response of the material will not be guaranteed if
271
the (transversely isotropic) model defined with the constants in Table 3 is
272
used. Finally, we have shown that this model, even though having theoretical
(r)
= (λ2 )−1/2 (j = 1, ..., m), which initially (j)
(r)
AC C
EP
TE D
(j)
16
ACCEPTED MANUSCRIPT
273
material-symmetries congruency, this congruency is not observed in practice
274
using the parameter identification procedure. (1)
w1 = 0.310 (2) w1 (3) w1
= 0.0 = 0.336
α1 = 4.354 α2 = 4.663 α3 = 1.0
β1 = 1.0 β2 = 8.523 β3 = 1.0
RI PT
µ1 = 2.16 × 10−1 MPa µ2 = 2.79 × 10−3 MPa µ3 = 5.39 MPa
SC
Table 3: Computed material parameters for the transversely isotropic model of Ref. [1] corresponding to the second approach in Figure 2. Initial guess (associated to an isotropic (r) (j) (j) material) given in Table 2 along with w1 = 1/3 (r = 1, 2, 3) and λ1 = (λ2 )−1/2 (j = 1, ..., m).
At this point, we note that, besides to the lack of numerical isotropic
276
material-symmetries congruency of this model, the optimization procedure
277
has resulted to be strongly dependent on the initial values given to the ma-
278
terial parameters. Hence, it is convenient that the material parameters of
279
a model have a clear physical meaning that facilitates the task of assigning
280
proper initial values. This fact is also encountered in isotropic materials [26].
281
Similar results are obtained if the two “experimental” discrete stress-
282
stretch distributions are assumed to be equal to the original measured data
283
points in the uniaxial test performed along the transverse direction, as it is
284
shown in Figure 3.
EP
TE D
M AN U
275
Again, the solution for the spline model is such that ω1# (E) = ω3# (E) =
286
ω (E) and E2 (E1 ) = −E1 /2. That is, the transversely isotropic contribution
287
Wtr (E, a0 ) in Eq. (12) vanishes and the model converges to the correspond-
288
ing isotropic solution. The exact predictions of the prescribed data points
289
provided by the computed isotropic limit model are shown in Figure 3.
AC C
285
290
In that Figure, we also represent the predictions given by the Itskov-
291
Aksel model with two different “best” solutions for the material parameters 17
ACCEPTED MANUSCRIPT
Experimental points
Prediction (Itskov-Aksel model - First approach)
Prediction (Spline model)
Prediction (Itskov-Aksel model - Second approach)
3
3
P (MPa)
P (MPa) 2
2
1
1
1
1.5
2
λ3
0
2.5
1
1.5
2
SC
0
RI PT
1
3
λ1
2.5
M AN U
Figure 3: Nominal stresses P3 and P1 in terms of the stretches λ3 and λ1 , respectively, from both uniaxial tests (black solid dots) performed on a transversely isotropic material in another hypothetical limit of isotropic behavior. Predictions obtained with the splinebased model and the Itskov-Aksel model of Ref. [1] with two different sets of material parameters. Original experimental data is obtained from the test in transverse direction from Diani et al [23]. Index numbering is consistent with previous Figures.
being investigated. The first set of them (first approach) corresponds to the
293
transversely isotropic solution given in Table 4. These material parameters
294
have been computed with the optimization procedure detailed in their work
295
using as initial values the results provided by these authors in Ref. [1],
296
i.e. Eq. (104) and the associated principal stretches λ1 (j = 1, ..., m) for
297
the test in the isotropic direction 2. A second set of transversely isotropic
298
material parameters (second approach) are provided in Table 5. In this case,
299
these material parameters have been obtained using as initial data the values
300
w1 = 1/3, r = 1, 2, 3, the stretches λ1 = (λ2 )−1/2 , j = 1, ..., m, and the
301
associated isotropic material parameters previously computed to fit only one
302
of the two stress-stretch distributions in Figure 3 using Eq. (100) of Ref. [1]
303
along with w1 = 1/3 (see Table 6). Again, we deduce from Tables 4 and
TE D
292
(r)
AC C
EP
(j)
(j)
(r)
18
(j)
ACCEPTED MANUSCRIPT
5 that numerical convergence to the values of the associated isotropic model
305
is not attained in any case, even though apparent isotropic predictions are
306
certainly observable in Figure 3. µ1 = 3.22 × 10−4 MPa µ2 = 1.17 × 10−3 MPa µ3 = 4.77 MPa
(1)
w1 = 0.0 (2) w1 (3) w1
= 0.136 = 0.330
α1 = 6.94 α2 = 7.18 α3 = 1.0
RI PT
304
β1 = 1.0 β2 = 11.1 β3 = 1.0
µ1 = 1.26 × 10−3 MPa µ2 = 1.82 × 10−3 MPa µ3 = 4.75 MPa
(1)
M AN U
SC
Table 4: Computed material parameters for the transversely isotropic model of Ref. [1] corresponding to the first approach in Figure 3. Initial guess (associated to a transversely isotropic material) obtained from Ref. [1].
w1 = 0.0 (2) w1 (3) w1
= 0.237 = 0.334
α1 = 6.35 α2 = 6.37 α3 = 1.0
β1 = 1.0 β2 = 12.4 β3 = 1.0
TE D
Table 5: Computed material parameters for the transversely isotropic model of Ref. [1] corresponding to the second approach in Figure 3. Initial guess (associated to an isotropic (j) (j) (r) material) given in Table 6 along with w1 = 1/3 (r = 1, 2, 3) and λ1 = (λ2 )−1/2 (j = 1, ..., m).
µ1 = 1.15 × 10−3 MPa µ2 = 3.80 × 10−3 MPa µ3 = 4.76 MPa
EP
α1 = 6.73 α2 = 6.72 α3 = 1.0
β1 = 1.01 β2 = 13.1 β3 = 1.0
307
AC C
Table 6: Material parameters of the isotropic model of Ref. [1] computed by means of the fitting of the experimental data (only one curve) of Figure 3.
3.3. Model of Holzapfel et al. [18]
308
Other types of strain energy functions widely encountered in the literature
309
are based on the notion of structural tensors and integrity basis, see Ref. 19
ACCEPTED MANUSCRIPT
[27]. For an incompressible transversely isotropic hyperelastic material, the
311
isochoric strain energy function is then formulated in terms of the following
312
four invariants I1 = C : I
I2 =
] 1[ (C : I)2 − C 2 : I 2
313
I5 = C 2 : A0
(17) (18)
SC
I4 = C : A0
RI PT
310
where C is the isochoric Cauchy-Green deformation tensor (i.e. det(C) =
315
I3 = 1) and A0 = a0 ⊗a0 is the structural tensor associated to the anisotropic
316
direction a0 . Hence, in the most general case W = W (I1 , I2 , I4 , I5 ), where
317
coupling terms between the four invariants may be present in the expression
318
of W. In biomechanics, it is convenient to split the strain energy function
319
into a part Wis associated with isotropic deformations and a part Wanis
320
associated with anisotropic deformations [28]. For the case at hand, we write
TE D
M AN U
314
W (I1 , I2 , I4 , I5 ) = Wis (I1 , I2 ) + Wtr (I4 , I5 )
(19)
Proceeding in that way, note that Eq. (12) has been retrieved, in this case
322
formulated in terms of C, i.e.
EP
321
(20)
AC C
W (C, a0 ) = Wis (C) + Wtr (C, a0 )
323
so material-symmetries congruency is guaranteed for these models from a
324
theoretical viewpoint. Eq. (19) is often simplified as (c.f. Ref. [29]) W (I1 , I4 ) = Wis (I1 ) + Wtr (I4 ) 20
(21)
ACCEPTED MANUSCRIPT
According to Refs. [30],[31] the additive, uncoupled, decomposition given in
326
Eq. (21) can be determined from a specific set of two experimental curves
327
as in Ref. [11]. That is, two independent behavior “curves” are used to
328
determine the constitutive functions Wis (I1 ) and Wtr (I4 ).
RI PT
325
In the following example we show that (reduced) models of the form of
330
Eq. (21) may also lack the isotropic material-symmetries congruency from
331
a numerical standpoint if the parameter fitting procedure is not properly
332
designed. As an example of Eq. (21) we take the transversely isotropic
333
counterpart of the well-known model of Holzapfel et al. [18] (see also Ref.
334
[32]), namely W (I1 , I4 ) =
M AN U
SC
329
( ) ) k1 ( c (I1 − 3) + exp k2 (I4 − 1)2 − 1 2 k2
(22)
This model is specifically intended for describing the mechanical behavior
336
of arterial walls, representing a simple starting point to model this type of
337
materials. However, one would expect that if two equal experimental data
338
points representations (obtained from a material tested in different directions,
339
as in the previous example) are used to determine the model constants c, k1
340
and k2 , then k1 should strictly vanish. This is, indeed, what we define as
341
numerical material-symmetries congruency.
EP
TE D
335
Following an analogous procedure to that of the previous example, we
343
have tried to fit the nominal stresses given by the model Eq. (22) to the
344
modified, purely isotropic, experimental data of Figure 1. For isochoric de-
345
formations λ1 λ2 λ3 = 1. Aside, for a uniaxial test performed in direction 1,
346
λ2 = λ3 if the plane 2 − 3 is isotropic, so λ3 = λ1
347
initial isotropic behavior for the initial guess values of the parameters, i.e.
AC C
342
−1/2
21
. We have assumed an
ACCEPTED MANUSCRIPT
348
for example c = 1 MPa and k1 = k2 = 0
(23) −1/2
along with the principal stretches in the “anisotropic” direction λ3 = λ1
350
for the uniaxial test performed over the isotropic direction (denoted herein
351
as axis 1). We note that the first partial derivative of W (I1 , I4 ) with respect
352
to I4 is given by
SC
( ) ∂W (I1 , I4 ) = 2k1 (I4 − 1) exp k2 (I4 − 1)2 ∂I4
RI PT
349
(24)
so the initial specification k2 = 0 gives no evaluation problems in the derived
354
stress expressions. In order to fit the stress predictions of the model to the
355
two point distributions P1 (λ1 ) and P3 (λ3 ) given in Figure 1, the Matlab
356
function lsqcurvefit has been employed. The objective function has been
357
defined as for the previous example. The solution for the material constants
358
that have been computed in this case are
TE D
M AN U
353
c = 1.659 MPa , k1 = 1.024 × 10−4 MPa and k2 = 0.211
361
EP
360
Moreover, as a result of the parameter identification procedure, the ob(j)
tained prediction for every “experimental” point (j) is also such that λ3 ̸= ( )−1/2 (j) λ1 for the uniaxial test performed over the “isotropic” direction.
AC C
359
(25)
362
From these results, it is apparent that the numerical procedure has not con-
363
verged to the corresponding isotropic solution included in the model given in
364
Eq. (22) because we have obtained Wtr (I4 ) ̸= 0. In Figure 4 it can be ob-
365
served the different responses predicted by the model over the two isotropic
22
ACCEPTED MANUSCRIPT
P (MPa) 4
Experimental points P3 (prediction)
RI PT
P1 (prediction)
3
SC
2
0 1
1.5
M AN U
1
2
λ
2.5
TE D
Figure 4: Nominal stresses P3 and P1 in terms of the stretches λ3 and λ1 , respectively, from both uniaxial tests (black solid dots) performed on a transversely isotropic material in a hypothetical limit of isotropic behavior and predictions obtained with the model of Holzapfel et al., Ref [18], initialized with isotropic material constants. Original experimental data is obtained from the test in calendering direction from Diani et al [23].
directions, respectively. We want to emphasize that the fact that k1 is rather
367
small, in relative terms (when compared to the isotropic material constant
368
c), causes that the isotropic behavior is dominant for moderate stretches in
369
both directions (so it is numerically congruent for moderately large strains).
370
However, for larger values of λ the anisotropic contribution becomes domi-
371
nant for the test in the “anisotropic” direction, which causes the response to
372
be different from the one in the transverse isotropic direction.
AC C
EP
366
373
Additionally, we show in Figure 5 the results obtained for the fitting of
374
this model to the two first Piola-Kirchhoff stress distributions of Figure 3.
23
ACCEPTED MANUSCRIPT
P (MPa) Experimental points
3
P3 (prediction)
RI PT
P1 (prediction) 2
0
1
1.5
2
SC
1
λ
2.5
M AN U
Figure 5: Nominal stresses P3 and P1 in terms of the stretches λ3 and λ1 , respectively, from both uniaxial tests (black solid dots) performed on a transversely isotropic material in another hypothetical limit of isotropic behavior. Predictions obtained with the model of Holzapfel et al., Ref [18], initialized with isotropic material constants. Original experimental data is obtained from the test in transverse direction from Diani et al [23].
Numerical convergence to the exact isotropic response of the model (with
376
k1 = 0) is not attained either, although better results are obtained for the
377
range of deformations under study.
TE D
375
Finally, we note that in order to prevent this undesired issue (specially
379
when slightly anisotropic, nearly isotropic materials are characterized with
380
anisotropic strain energy functions) it is possible to design a constrained
381
parameter fitting procedure so as to attain numerical material-symmetries
382
congruency for analytically congruent models.
383
4. Conclusions
AC C
EP
378
384
In this paper we have presented the material-symmetries congruency
385
property for anisotropic hyperelastic materials. This is a desirable material 24
ACCEPTED MANUSCRIPT
property for several reasons. For example, it guarantees that the isotropic be-
387
havior of the matrix is recovered when the volume of reinforcement vanishes,
388
as one should expect. It also guarantees for orthotropic materials that when
389
the volume of one of the fibers vanishes, the transversely isotropic behavior
390
is recovered.
RI PT
386
The material-symmetries congruency property is not only an analytical
392
issue, it is also a numerical one. The lack of numerical material-symmetries
393
congruency in analytically congruent materials is due to the design of the
394
parameter-fitting procedure inherent to the model from composite material
395
experiments. The use of material parameters with a clear physical inter-
396
pretation in the material models easies the development of procedures that
397
preserve the material-symmetries congruency not only from an analytical but
398
also from a numerical point of view.
TE D
Acknowledgements
M AN U
SC
391
Partial financial support for this work has been given by grant DPI201126635 from the Direcci´on General de Proyectos de Investigaci´on of the Min-
AC C
References
EP
isterio de Econom´ıa y Competitividad of Spain.
[1] M. Itskov, N. Aksel. A class of orthotropic and transversely isotropic hyperelastic constitutive models based on a polyconvex strain energy function. International Journal of Solids and Structures, 41, 3833–48, 2004. 25
ACCEPTED MANUSCRIPT
[2] M. Latorre, F.J. Mont´ans. What-You-Prescribe-Is-What-You-Get orthotropic hyperelasticity. Computational Mechanics, 53(6), 1279–1298,
RI PT
2014. [3] G.A. Holzapfel. Nonlinear Solid Mechanics. A Continuum Approach For Engineering. Wiley, Chichester, 2000.
SC
[4] G.A. Holzapfel, R.W. Ogden. Constitutive modelling of passive myocardium. A structurally-based framework for material characterization.
3445–3475, 2009.
M AN U
Philosophical Transactions of the Royal Society of London, Series A, 367,
[5] S. G¨oktepem S.N.S. Acharya, J. Wong, E. Kuhl. Computational modeling of passive myocardium. International Journal for Numerical Methods in Engineering, 27, 1–12, 2011.
TE D
[6] J. Schr¨oder, P. Neff, D. Balzani. A variational stable anisotropic hyperelasticity. International Journal of Solids and Structures, 42, 4352–4371, 2005.
EP
[7] M.P. Nash, P.J. Hunter. Computational mechanics of the heart. Journal of Elasticity, 61, 113–141, 2000.
AC C
[8] D. Tang, C. Yang, T. Geva, R. Rathod, H. Yamauchi, V. Gooty, A. Tang, M.H. Kural, K.L. Billiar, G. Gaudette, P.J. del Nido. A multiphysics modeling approach to develop right ventricle pulmonary valve replacement surgical procedures with a contracting band to improve ventricle ejection fraction. Computers & Structures, 122, 78–87, 2013.
26
ACCEPTED MANUSCRIPT
[9] C.O. Horgan, J.G. Murphy. Simple shearing of soft biological tissues. Proceedings of the Royal Society of London, Series A, 467, 760–777,
RI PT
2010. [10] K. Costa, J. Holmes, A. McCulloch. Modelling cardiac mechanical prop-
ciety of London, Series A, 359, 1233–1250, 2001.
SC
erties in three dimensions. Philosophical Transactions of the Royal So-
[11] T. Gasser, R.W. Ogden, G.A. Holzapfel. Hyperelastic modelling of arterial layers with distributed collagen fibre orientations. Journal of the
M AN U
Royal Society Interface, 3, 13–35, 2006.
[12] H. Schmid, P. O’Callagan, M.P. Nash, W. Lin, I.J. LeGrice, B.H. Smaill, A.A. Young,P.J. Hunter. Myocarcial material parameter estimation. A non-homogeneous finite element study from simple shear tests. Biome-
TE D
chanics and Modeling in Mechanobiology, 7, 161–173, 2008. [13] S. Rossi, T. Lassila, R. Ruiz-Baier, A. Sequeira, A. Quarteroni. Thermodynamically consistent orthotropic activation model capturing ventricular systolic wall thickning in cardiac electromechanics.
EP
European Journal of Mechanics, A/Solids, 2014, in Press. DOI: 10.1016/j.euromechsol.2013.10.009.
AC C
[14] F. Vogel, S. G¨oktepe, P. Steinmann, E.Kuhl. Modeling and simulation of viscous electro-active polymers. European Journal of Mechanics, A/Solids, 2014, In Press. DOI: 10.1016/j.euromechsol.2014.02.001. [15] J. Merodio, R.W. Ogden. Instabilities and loss of ellipticity in fiber-
27
ACCEPTED MANUSCRIPT
reinforced compressible non-linearly elastic solids under plane deformation. International Journal of Solids and Structures, 40, 4707–4727, 2003.
RI PT
[16] J. Merodio, R.W. Ogden. Material instabilities in fiber-reinforced non-
linearly elastic solids under plane deformation. Archives of Mechanics, 54, 525–552, 2002.
SC
[17] J.E. Bischoff, E.M. Arruda, K. Grosh. Finite element simulations of orthotropic hyperelasticity. Finite Elements in Analysis and Design, 38,
M AN U
983–998, 2002.
[18] G.A. Holzapfel, T. Gasser, R.W.Ogden. A new constitutive framework for arterial wall mechanics and a comparative study of material models. Journal of Elasticity, 61, 1–18, 2000.
[19] M. Latorre, F.J. Mont´ans. Extension of the Sussman–Bathe spline-based
TE D
hyperelastic model to incompressible transversely isotropic materials. Computers & Structures, 122, 13–26, 2013. [20] T. Sussman, K.J. Bathe. A Model of Incompressible Isotropic Hy-
EP
perelastic Material Behavior using Spline Interpolations of TensionCompression Test Data. Communications in Numerical Methods in En-
AC C
gineering, 25, 53–63, 2009.
[21] C. Truesdell, W. Noll. The non-linear field theories of mechanics. Handbuch der Physik, Vol. III/3, Springer, Berlin, New York, 2004. [22] R.W. Ogden. Nonlinear Elastic Deformations. Dover, Mineola, New York, 1997. 28
ACCEPTED MANUSCRIPT
[23] J. Diani, M. Brieu, J.M. Vacherand, A. Rezgui. Directional model for isotropic and anisotropic hyperelastic rubber-like materials. Mechanics
RI PT
of Materials, 36, 313–21, 2004. [24] L.L. Demer, F.C.P. Yin. Passive biaxial mechanical properties of isolated canine myocardium. Journal of Physiology, 339, 615–630, 1983.
[25] E.H. Twizell, R.W. Ogden. Non-linear optimization of the material con-
SC
stants in Ogden’s stress-deformation relation for incompressible isotropic elastic materials. J. Aust. Math. Soc. 24, 424–434, 1983.
M AN U
[26] R.W. Ogden, G. Saccomandi, I. Sgura. Fitting hyperelastic models to experimental data. Computational Mechanics, 34(6), 484–502, 2004. [27] A.J.M. Spencer. Constitutive theory for strongly anisotropic solids. In: A.J.M. Spencer (ed.), Continuum Theory of the Mechanics of FibreReinforced Composites, CISM Courses and Lectures No. 282, Inter-
1984.
TE D
national Centre for Mechanical Sciences, Springer-Verlag, Wien, 1–32,
[28] H.W. Weizs¨acker, G.A. Holzapfel, G.W. Desch, K. Pascale. Strain en-
EP
ergy density function for arteries from different topographical sites. Biomediz. Tech., 40, 139–141, 1995.
AC C
[29] J.D. Humphrey. Cardiovascular Solid Mechanics. Cells, Tissues, and Organs. Springer, New York, 2002. [30] J.D. Humphrey,F.C.P. Yin. On constitutive relations and finite deformations of passive cardiac tissue. Part I: A pseudo-strain energy function. Journal of Biomechanical Engineering, 109, 298–304, 1987. 29
ACCEPTED MANUSCRIPT
[31] G.A. Holzapfel, R.W. Ogden. On planar biaxial tests for anisotropic nonlinearly elastic solids. A continuum mechanical framework. Mathe-
RI PT
matics and Mechanics of Solids, 14, 474–489, 2009. [32] G.A. Holzapfel, R.W. Ogden. Constitutive modelling of arteries. Proc. R. Soc. A, 466, 1551–1597, 2010.
SC
[33] K.C.Valanis, R.F.Landel. The strain-energy function of a hyperelastic
material in terms of the extension ratios. Journal of Applied Physics,
AC C
EP
TE D
M AN U
38, 2997–3002, 1967.
30
ACCEPTED MANUSCRIPT Highlights: • •
AC C
EP
TE D
M AN U
SC
RI PT
•
The property of material-symmetries congruency for orthotropic and transversely isotropic hyperelastic materials is presented An anisotropic material model holding this property recovers an isotropic behavior when the experimental data used to calibrate the model corresponds to that of an isotropic material A model may have analytical material-symmetries congruency and still lack numerical one due to the procedure employed to obtain the material parameters from experimental data
S0997-7538(15)00027-3
DOI:
10.1016/j.euromechsol.2015.03.007
Reference:
EJMSOL 3175
To appear in:
European Journal of Mechanics / A Solids
Received Date: 29 August 2014 Revised Date:
22 January 2015
Accepted Date: 23 March 2015
Please cite this article as: Latorre, M., Montáns, F.J., Material-symmetries congruency in transversely isotropic and orthotropic hyperelastic materials, European Journal of Mechanics / A Solids (2015), doi: 10.1016/j.euromechsol.2015.03.007. This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT
Marcos Latorre, Francisco Javier Mont´ans∗
RI PT
Material-symmetries congruency in transversely isotropic and orthotropic hyperelastic materials
SC
Escuela T´ecnica Superior de Ingenier´ıa Aeron´ autica y del Espacio Universidad Polit´ecnica de Madrid Plaza Cardenal Cisneros, 3, 28040-Madrid, Spain
Abstract
M AN U
In this work we present the material-symmetries congruency property for anisotropic hyperelastic materials. A transversely isotropic or orthotropic material model holding this property guarantees the recovery of the isotropic behavior when, for example, the volume of reinforcement in a composite material vanishes. A material presenting material-symmetries congruency
TE D
must show it not only from an analytical, theoretical point of view, but also from a numerical, practical point of view. The former may be obtained from construction of the stored energy function, whereas the latter must be obtained guaranteeing that the parameter-fitting procedure yields material
EP
parameters which combined with the specific model results in an isotropic behavior prediction for isotropic materials. We show that some anisotropic
AC C
models do not present either analytical or numerical material-symmetries congruency.
Keywords: Hyperelasticity, anisotropy, incompressible materials. ∗
Corresponding author. Tel.:+34 637 908 304. Email addresses: [email protected] (Marcos Latorre), [email protected] (Francisco Javier Mont´ans) Preprint submitted to European Journal of Mechanics - A/Solids
January 22, 2015
ACCEPTED MANUSCRIPT
1
1. Introduction Many material models have been developed for analyzing hyperelastic
3
transversely isotropic materials and for orthotropic materials, see for ex-
4
ample [1],[2],[3],[4],[5],[6],[7],[8],[9],[10],[11],[12],[13],[14],[15],[16],[18],[19],[23],
5
among others. The main objective of those material models is to account for
6
different material behavior encountered in different material directions for
7
some materials like rubber [1],[14],[17], fibre-reinforced composites [15],[16]
8
and biological materials [4],[5],[6],[7],[8],[9],[10],[11],[12],[13],[18],[19].
M AN U
SC
RI PT
2
The different material properties in the different directions are frequently
10
due to a reinforcing constituent which has some preference in alignment.
11
This reinforcement may be laid-out during manufacturing or, as in biological
12
tissues, by nature. A common approach is to test the composite material as
13
a whole, see for example [12], [17], [23], [24], among others.
TE D
9
Hyperelastic material models typically consist on a stored energy function
15
of some deformation-based invariants [3], some of which are derived from
16
structural tensors obtained from the preferred material directions. These
17
energy functions have a pre-defined shape and some material parameters
18
that are obtained fitting experimental values from some tests using a proper
19
optimization algorithm, typically the Levenberg-Marquardt algorithm, see
20
for example [1], [25] and [26]. Some proposed stored energy functions take
21
into account the fact that anisotropy is mainly due to reinforcement and they
22
separate the stored energy into an isotropic (bulk) part and a reinforcement
23
part, see for example [3], [15], [18]. In fact, these models may include specific
24
parameters that specifically take into account the amount of anisotropy or
AC C
EP
14
2
ACCEPTED MANUSCRIPT
25
reinforcement [3], [15], [16]. Other proposals are more phenomenological and
26
consider from the onset the material as a whole, see for example [1],[2],[19]. Furthermore, in anisotropy it is also frequent to have available less exper-
28
imental data than those necessary to fully characterize the material [31] and
29
then the missing information may strongly influence the overall behavior [2].
30
Many of these available constitutive models for transversely isotropic and
31
orthotropic materials intend to represent the behavior of different materi-
32
als with different degrees of reinforcement, and they intend to do so with
33
parameters obtained from some tests over the composite. One would ex-
34
pect that when the reinforcement volume goes to zero in the composite, the
35
material captures the behavior of the isotropic matrix in an isotropic way.
36
Furthermore, in a orthotropic material, if the reinforcement in one direction
37
vanishes, one would also expect that the model captures the behavior of the
38
transversely isotropic material in a transversely isotropic manner. We say of
39
a model behaving this way that the model has material-symmetries congru-
40
ency because it converges to an isotropic material when the data corresponds
41
to that of an isotropic material.
TE D
M AN U
SC
RI PT
27
The material-symmetries congruency is not only an analytical issue on
43
the form of the stored energy function. A hyperelastic constitutive model
44
may present congruency from an analytical standpoint by construction, but
45
still lack that congruency from a numerical point of view, i.e. the model
46
using the parameters obtained from (quasi-)isotropic experimental curves
47
may not behave in a (quasi-)isotropic manner, simply because the parameter
48
fitting algorithm does not distinguish between quasi isotropic, transversely
49
isotropic and orthotropic materials. Many possible minima can be obtained
AC C
EP
42
3
ACCEPTED MANUSCRIPT
during the parameter fitting procedure which strongly influence the predicted
51
behavior even in the isotropic case [26]. In this case, if a model with analyt-
52
ical material-symmetries congruency loses this congruency due to the actual
53
parameters identified from the experiments, we say that whereas the model
54
has analytical material-symmetries congruency, it lacks the numerical one.
55
Material-symmetries congruency is specially important when some experi-
56
mental data are not available, and also in models that use decompositions of
57
the type of Valanis-Landel [33] because this decomposition has been verified
58
for isotropic materials both experimentally and analytically up to high orders
59
of strain [22], but in general remains to be an assumption on the form of the
60
stored energy for anisotropic materials. If material-symmetries congruency is
61
preserved, only the decomposition on the deviation from isotropy remains as
62
an assumption for moderate large strains. Furthermore, this deviation will
63
tend to vanish in the limit of isotropic behavior.
M AN U
SC
RI PT
50
The purpose of this paper is to present the material-symmetries congru-
65
ency for anisotropic hyperelastic materials and to analyze some models as
66
examples, both from an analytical perspective and from a numerical one.
67
2. Theoretical and numerical material-symmetries congruency
EP
Let us define the space L of the proper orthogonal tensors Q
AC C
68
TE D
64
Q ∈ L / QQT = I
4
and
det(Q) = 1
(1)
ACCEPTED MANUSCRIPT
69
and the symmetry group for orthotropic materials Lor ⊂ L containing the
70
set of rotation tensors Qor ∈ Lor = {Qπa0 , Qπb }
RI PT
(2)
0
71
where Qπa0 = a0 ⊗a0 −b0 ⊗b0 −c0 ⊗c0 and Qπb = −a0 ⊗a0 +b0 ⊗b0 −c0 ⊗c0 0
represent two rotations of π radians around two preferred material directions,
73
defined by the orthogonal vectors a0 and b0 . Note that the tensor Qπc0 =
74
Qπa0 Qπb associated to the third preferred direction, c0 = a0 × b0 , is also 0
75
included in Lor .
M AN U
76
SC
72
For an anisotropic hyperelastic material, in general ( ) W (E) ̸= W QEQT
(3)
for a given Q and a Lagrangian strain measure E. If the material is or-
78
thotropic, it has two (subsequently, three) planes of material symmetry, with
79
a0 and b0 defining these two planes and being perpendicular to them. Then,
80
it can be stated that
TE D
77
EP
( ) W (E, a0 , b0 ) = W Qor EQTor , a0 , b0
(4)
where Qor stands now for any arbitrary rotation included in the symmetry
82
group Lor . If Q ∈ / Lor , then W(E, a0 , b0 ) ̸= W(QEQT , a0 , b0 ) in general. If
83
completely arbitrary rotations Q are considered, this more general invariance
AC C
81
5
ACCEPTED MANUSCRIPT
84
relation (i.e. frame independence) is to be fulfilled ( ) W (E, a0 , b0 ) = W QEQT , Qa0 , Qb0
RI PT
(5)
that is, W(E, a0 , b0 ) must be an isotropic function with respect to its three
86
arguments simultaneously. Note that Eq. (4) stands for symmetry consider-
87
ations of the strain energy function while Eq. (5) stands for frame indepen-
88
dence so W(E, a0 , b0 ) must satisfy both of them. A material model satisfying
89
only Eq. (5) and not Eq. (4) has not necessarily orthotropic symmetries and
90
it would not be an acceptable formulation for an orthotropic material.
M AN U
SC
85
Another property that the function W(E, a0 , b0 ) must obey and which
92
seem to have never been suggested in the literature (to the best of the author’s
93
knowledge) is what we call material-symmetries congruency. This require-
94
ment states that for an orthotropic hyperelastic material, the model being
95
employed should tend to an isotropic formulation when slightly orthotropic,
96
nearly isotropic, materials are being characterized. That is, the tendency of
97
the orthotropic strain energy function when the limit of isotropic behavior is
98
considered should be such that
TE D
91
EP
( ) W (E, a0 , b0 ) ≈ W QEQT , a0 , b0
100
for any proper orthogonal tensor Q. Moreover, if an isotropic material is
AC C
99
(6)
characterized using an orthotropic model, then it is required that ) ( W (E, a0 , b0 ) = W QEQT , a0 , b0
6
(7)
ACCEPTED MANUSCRIPT
for any Q, as it effectively occurs for isotropic materials, i.e. W(E) =
102
W(QEQT ). Note the difference between Eq. (4), only forced to be ful-
103
filled for Qor ∈ Lor , and Eqs. (6) and (7), enforced for arbitrary rotations
104
Q ∈ L.
RI PT
101
We find in the literature some formulations satisfying the requirements
106
given in Eqs. (4) and (5), but which do not fulfill Eqs. (6) and (7) either
107
from a theoretical or from a, practical, numerical point of view. For example,
108
consider the incompressible uncoupled orthotropic model
SC
105
M AN U
W (E, a0 , b0 ) = ω11 (E11 )+ω22 (E22 )+ω33 (E33 )+2ω12 (E12 )+2ω23 (E23 )+2ω31 (E31 ) (8)
that we present in Ref. [2]. In Eq. (8), the arguments of the ω functions
110
are the components of the (deviatoric) Hencky strain tensor in the material
111
preferred basis X pr = {e1 , e2 , e3 } = {a0 , b0 , c0 }, that is Eij = ei · Eej , with
112
i, j = {1, 2, 3}. Note that these logarithmic strain components are deviatoric
113
(traceless) due to the incompressibility condition, and that the product of the
114
principal stretches is one because of the same reason. Thus, the invariance
115
relation given in Eq. (5) is automatically satisfied and the model is invariant
116
under simultaneous rotations of E, a0 and b0 . However, note that if any
117
of the shear terms ωij , i ̸= j, is not an even function of its corresponding
118
argument Eij , then the model does not satisfy the symmetry requirement of
119
Eq. (4) and the model would not be physically correct. Therefore, we only
120
consider even functions for the shear terms ωij as we properly address in
121
Ref. [2]. This fact shows that both symmetry and invariance requirements
122
must be satisfied for any orthotropic model and that Eq. (5) is not a suffi-
AC C
EP
TE D
109
7
ACCEPTED MANUSCRIPT
cient condition for a material model to reproduce an orthotropic behavior,
124
in contrast to what one could erroneously infer from, for example, Ref. [3]
125
(cf. page 274).
RI PT
123
It can be readily shown that if the orthotropic model of Eq. (8) is used
127
to characterize an assumed symmetric-Valanis-Landel-type incompressible
128
isotropic material from the corresponding set of six different experimental
129
data curves, then the computed model using the procedure detailed in Ref.
130
[2] reduces to
SC
126
M AN U
W (E, a0 , b0 ) = ω (E11 ) + ω (E22 ) + ω (E33 ) + 2ω (E12 ) + 2ω (E23 ) + 2ω (E31 ) (9)
with ω (E) being an even function of the strain component E. If as a result of
132
the model determination procedure the function ω contains quadratic terms
133
on its argument, it can be shown that the complete invariant E : E is
134
recovered in W multiplied by the proper coefficient. On the contrary, if as a
135
result of the model determination procedure the function ω contains fourth-
136
order terms on its argument, it can be shown that the complete invariant
137
E 4 : I is not recovered in W multiplied by the proper coefficient, so Eq. (4)
138
does not hold. Hence, although an arbitrary strain energy for an orthotropic
139
material can be approximated as closely as desired by a suitable polynomial of
140
the corresponding symmetry-preserving invariants (cf. Ref. [22], pages 211–
141
212), we remark that this “suitable” polynomial should contain the isotropic
142
formulation for the corresponding limit.
AC C
EP
TE D
131
143
Regarding the previous formulation, a possibility to modify the model to
144
obtain material-symmetries congruency simply consists in adding an isotropic
8
ACCEPTED MANUSCRIPT
145
contribution to the orthotropic one W (E, a0 , b0 ) = Wis (E) + Wor (E, a0 , b0 )
RI PT
(10)
146
where Wis (E) is the Valanis-Landel-type isotropic model expressed in terms
147
of principal logarithmic strains defined by Sussman and Bathe in Ref. [20],
148
i.e.
SC
Wis (E) = ω (E1 ) + ω (E2 ) + ω (E3 )
(11)
and Wor (E, a0 , b0 ) is given by Eq. (8). With respect to the modified model
150
provided in Eq. (10), the material-symmetries congruency requirements Eqs.
151
(6) and (7) will be effectively satisfied if both Wor (E, a0 , b0 ) tends to vanish
152
when slightly orthotropic, nearly isotropic, materials are characterized and
153
Wor (E, a0 , b0 ) = 0 in the limit of pure isotropic behavior, respectively. In
154
this case, due to the use of the spline-based methodology and the inversion
155
formula (cf. Refs. [19], [2], [20]), the equilibrium equations of the tests
156
are solved exactly for the isotropic model. This important fact implies that
157
Wis (E) exactly reproduces the prescribed isotropic behavior and that, as a
158
consequence, Wor (E, a0 , b0 ) vanishes for these type of materials.
159
3. Examples
160
3.1. Spline-based model according to Latorre and Mont´ans [19]
AC C
EP
TE D
M AN U
149
161
In Ref. [19] we present a spline-based transversely isotropic model which
162
is capable of reproducing the uniaxial experimental data of Diani et al. [23]
163
over transversely isotropic calendered rubber with a very good (exact in
164
practice) agreement. In that work, we also compare the results obtained with 9
ACCEPTED MANUSCRIPT
our model to those predicted by the model of Itskov and Aksel [1], showing a
166
better accuracy achieved by the spline model. Following similar reasonings to
167
those detailed above, it is readily shown that the spline-based model that we
168
presented in Ref. [19] does not have isotropic material-symmetry congruency.
169
However, as mentioned above, the addition of the isotropic contribution of
170
Eq. (11) circumvents this issue.
RI PT
165
In this example, and the next ones, we study two limit cases in which the
172
two curves presented by Diani et al coalesce to one of both curves. The aim
173
of these examples is to predict with transversely isotropic models two possi-
174
ble isotropic behaviors obtained as isotropic limits of a transversely isotropic
175
material and analyze the numerical material-symmetries congruency of the
176
models. Since the uniaxial tests are performed on an initially assumed trans-
177
versely isotropic material and along its “preferred” material directions, no
178
shear terms are needed and the extended spline-based transversely isotropic
179
model of Ref. [19] reduces to
TE D
M AN U
SC
171
W (E, a0 ) = Wis (E) + Wtr (E, a0 )
EP
= ω1# (E1 ) + ω1# (E2 ) + ω3# (E3 )
(12) (13)
where the simplifications ωn# (E) = ω (E) + ωn (E), n = 1, 3, are considered
181
for practical purposes. As in Ref. [19], the subscripts 1 and 2 are associated
182
to two orthogonal directions within the isotropic plane and the subscript 3
183
represents the assumed preferred direction of the material.
184
AC C
180
For the first case, we prescribe two discrete stress-stretch distributions
10
M AN U
SC
RI PT
ACCEPTED MANUSCRIPT
Figure 1: Nominal stresses P3 and P1 in terms of the stretches λ3 and λ1 , respectively, from both uniaxial tests (black solid dots) performed on a transversely isotropic material in a hypothetical limit of isotropic behavior and predictions obtained with the spline-based model. Original experimental data is obtained from the test in calendering direction from Diani et al [23].
which are identical to the original measured data points in the uniaxial test
186
performed along the calendering direction, i.e. direction 3, as it is shown in
187
Figure 1 with black solid dots. Using these “experimental” distributions as
188
initial data of the procedure detailed in Table 3 of Ref. [19] (where the func-
189
tions ωn are simply substituted by ωn# ), the solution functions are found to be
190
ω1# (E) = ω3# (E), with the transverse strain distribution E2 (E1 ) = −E1 /2
191
(ν12 = 1/2) computed also as a part of the solution. Clearly, the solution
192
converges to the corresponding isotropic solution. That is, considering the
193
relations ωn# (E) = ω (E) + ωn (E), n = 1, 3, we deduce
AC C
EP
TE D
185
11
ACCEPTED MANUSCRIPT
ωn (E) = 0
(14)
RI PT
ωn# (E) = ω (E) ;
These results indicate that the model given in Eqs. (12)–(13) has also
195
isotropic material-symmetry congruency from a numerical point of view, as
196
one would expect. The stress predictions in both directions provided by
197
the computed model are also shown in Figure 1. A strictly exact fit of the
198
experimental data is attained, as it should occur within the isotropic context
199
(cf. Refs. [20] and [19]).
200
3.2. Model of Itskov and Aksel [1]
M AN U
SC
194
Regarding the Itskov-Aksel work of Ref. [1], the orthotropic model pro-
202
vided therein in terms of generalized invariants, namely Eq. (58) of that
203
reference, includes the isotropic generalized Mooney-Rivlin model, Eq. (68)
204
of the same reference, as a particular case. Hence, the proposed strain energy
205
function has theoretical material-symmetries congruency. However, when cer-
206
tain slightly transversely isotropic, nearly isotropic, materials are character-
207
ized with the parameter identification procedure detailed in their work, the
208
tendency to the Mooney-Rivlin formulation may not be observed. Further-
209
more, when two identical stress-strain distributions are used to initialize the
210
parameter identification process, i.e. pure isotropic behavior is prescribed,
211
several specific transversely isotropic solutions are found to be more optimum
212
(in a least-squares sense) than the isotropic solution. This is due to the fact
213
that the optimization process employed to identify the material parameters
214
does not ensure that their limit values for an isotropic material are obtained.
AC C
EP
TE D
201
12
ACCEPTED MANUSCRIPT
This fact manifests that the Itskov-Aksel formulation, as given in their pa-
216
per, does not have material-symmetries congruency from a numerical point
217
of view.
218
219
RI PT
215
In order to show this unexpected result, the transversely isotropic incompressible model by Itskov and Aksel, namely Eq. (95) of Ref. [1]
[( ( )αr ] )βr 3 3 1 ∑ ∑ 1 (r) (r) W= µr wi λ2i −1 + wi λ−2 − 1 i 4 r=1 αr βr i=1 i=1 (r)
SC
3 1∑
where µr , αr , βr and wi
221
relations in Eq. (98) of the same Reference, i.e. (r)
are material parameters, along with the symmetry
(r)
w2 = w3 =
M AN U
220
(15)
) 1( (r) 1 − w1 , r = 1, 2, 3 2
(16)
are employed to characterize the isotropic materials of the preceding example
223
using the minimization procedure detailed in Ref. [1]. A standard least
224
squares curve-fitting algorithm (specifically, the Matlab function lsqcurvefit)
225
has been used to fit the two point distributions shown in Figure 1 using Eqs.
226
(100)–(102) of Ref. [1]. The optimization procedure has been initialized
227
with the set of material parameters given in Ref. [1], i.e. Eq. (104), and
228
the principal stretches in the calendering direction λ1 (j = 1, ..., m) for the
229
uniaxial test in the transversely isotropic direction 2 obtained as solution of
230
their parameter identification procedure (note that in this example we are
231
using the index numbering of Ref. [1], which is different from the one used in
232
the preceding example). All these material parameters and stretches initially
233
correspond to a transversely isotropic behavior. In this first attempt to assess
EP
TE D
222
AC C
(j)
13
ACCEPTED MANUSCRIPT
234
the numerical material-symmetries congruency of the model, the values of
235
the parameters which have been found as “best” solution using these initial
236
values are given in Table 1. The fact that w1 ̸= 1/3, r = 1, 2, 3 indicates
237
that the computed solution has not converged to the associated isotropic
238
solution (cf. Eqs. (67) and (98) in Ref. [1]). In Figure 2, it can be observed
239
that the stresses predicted by the model defined with the parameters given in
240
Table 1 reproduce the tendency of the “experimental” data, although little
241
differences between both curves (due to the lack of isotropy) are observable. (1)
w1 = 0.900 (2) w1 (3) w1
α1 = 4.27 α2 = 6.84 α3 = 1.0
β1 = 9.44 β2 = 7.78 β3 = 1.59
M AN U
µ1 = 3.84 × 10−3 MPa µ2 = 2.03 × 10−3 MPa µ3 = 5.44 MPa
SC
RI PT
(r)
= 0.0
= 0.323
Table 1: Computed material parameters for the transversely isotropic model of Ref. [1] corresponding to the first approach in Figure 2. Initial guess (associated to a transversely isotropic material) obtained from Ref. [1].
Clearly, the foregoing solution for the set of material parameters being
243
investigated corresponds to the nearest local minimum to the initial values
244
given in Eq. (104) of Ref. [1]. Obviously, a local solution of any optimization
245
procedure is strongly affected by the initial numerical values given to the
246
unknowns of the optimization problem. Hence, one may arguably think that
247
the isotropic limit solution could be obtained if another, more appropriate, set
248
of initial values nearer to the isotropic solution is taken. In a second attempt
249
to assess the numerical material-symmetries congruency of the model, we
250
have previously obtained the values of the nine material parameters µr , αr
251
and βr , r = 1, 2, 3, of the isotropic model (with w1 = 1/3, r = 1, 2, 3) which
252
best fit one of the two identical experimental distributions of Figure 1. For
AC C
EP
TE D
242
(r)
14
TE D
M AN U
SC
RI PT
ACCEPTED MANUSCRIPT
AC C
EP
Figure 2: Nominal stresses P3 and P1 in terms of the stretches λ3 and λ1 , respectively, from both uniaxial tests (black solid dots) performed on a transversely isotropic material in a hypothetical limit of isotropic behavior and predictions obtained with the model of Itskov and Aksel, Ref [1], with two different sets of material parameters. Original experimental data is obtained from the test in calendering direction from Diani et al [23]. Index numbering is consistent with Figure 1.
15
ACCEPTED MANUSCRIPT
253
this fist step, the experimental data points have been fitted by means of Eq.
254
(100) of Ref. [1] (with w1 = 1/3). The values of the parameters which have
255
been found as “best” solution for the isotropic model are given in Table 2.
256
These values provide a very good agreement of the isotropic model with the
257
experimental data. α1 = 4.495 α2 = 4.338 α3 = 1.003
β1 = 1.004 β2 = 9.017 β3 = 1.0
SC
µ1 = 2.30 × 10−1 MPa µ2 = 2.46 × 10−3 MPa µ3 = 5.39 MPa
RI PT
(r)
M AN U
Table 2: Material parameters of the isotropic model of Ref. [1] computed by means of the fitting of the experimental data (only one curve) of Figure 2.
258
Subsequently, as before, we have proceeded to fit the two point distribu-
259
tions of Figure 1 using Eqs. (100)–(102) of Ref. [1]. In this second step,
260
the optimization procedure has been initialized with the computed material
261
parameters given in Table 2 together with the values w1 = 1/3, r = 1, 2, 3,
262
and the principal stretches λ1
263
describe an isotropic model. Again, a solution corresponding to a trans-
264
versely isotropic description of the model, i.e. with w1 ̸= 1/3, r = 1, 2, 3, is
265
found to be the “best” solution for the considered initial values. In Figure
266
2, no noticeable differences are observable between the stress distributions
267
predicted by the model in both directions. However, this is just an apparent
268
isotropic response of the model in a specific situation (which has been inten-
269
tionally fitted). We want to emphasize that in other more generic situations
270
an isotropic mechanical response of the material will not be guaranteed if
271
the (transversely isotropic) model defined with the constants in Table 3 is
272
used. Finally, we have shown that this model, even though having theoretical
(r)
= (λ2 )−1/2 (j = 1, ..., m), which initially (j)
(r)
AC C
EP
TE D
(j)
16
ACCEPTED MANUSCRIPT
273
material-symmetries congruency, this congruency is not observed in practice
274
using the parameter identification procedure. (1)
w1 = 0.310 (2) w1 (3) w1
= 0.0 = 0.336
α1 = 4.354 α2 = 4.663 α3 = 1.0
β1 = 1.0 β2 = 8.523 β3 = 1.0
RI PT
µ1 = 2.16 × 10−1 MPa µ2 = 2.79 × 10−3 MPa µ3 = 5.39 MPa
SC
Table 3: Computed material parameters for the transversely isotropic model of Ref. [1] corresponding to the second approach in Figure 2. Initial guess (associated to an isotropic (r) (j) (j) material) given in Table 2 along with w1 = 1/3 (r = 1, 2, 3) and λ1 = (λ2 )−1/2 (j = 1, ..., m).
At this point, we note that, besides to the lack of numerical isotropic
276
material-symmetries congruency of this model, the optimization procedure
277
has resulted to be strongly dependent on the initial values given to the ma-
278
terial parameters. Hence, it is convenient that the material parameters of
279
a model have a clear physical meaning that facilitates the task of assigning
280
proper initial values. This fact is also encountered in isotropic materials [26].
281
Similar results are obtained if the two “experimental” discrete stress-
282
stretch distributions are assumed to be equal to the original measured data
283
points in the uniaxial test performed along the transverse direction, as it is
284
shown in Figure 3.
EP
TE D
M AN U
275
Again, the solution for the spline model is such that ω1# (E) = ω3# (E) =
286
ω (E) and E2 (E1 ) = −E1 /2. That is, the transversely isotropic contribution
287
Wtr (E, a0 ) in Eq. (12) vanishes and the model converges to the correspond-
288
ing isotropic solution. The exact predictions of the prescribed data points
289
provided by the computed isotropic limit model are shown in Figure 3.
AC C
285
290
In that Figure, we also represent the predictions given by the Itskov-
291
Aksel model with two different “best” solutions for the material parameters 17
ACCEPTED MANUSCRIPT
Experimental points
Prediction (Itskov-Aksel model - First approach)
Prediction (Spline model)
Prediction (Itskov-Aksel model - Second approach)
3
3
P (MPa)
P (MPa) 2
2
1
1
1
1.5
2
λ3
0
2.5
1
1.5
2
SC
0
RI PT
1
3
λ1
2.5
M AN U
Figure 3: Nominal stresses P3 and P1 in terms of the stretches λ3 and λ1 , respectively, from both uniaxial tests (black solid dots) performed on a transversely isotropic material in another hypothetical limit of isotropic behavior. Predictions obtained with the splinebased model and the Itskov-Aksel model of Ref. [1] with two different sets of material parameters. Original experimental data is obtained from the test in transverse direction from Diani et al [23]. Index numbering is consistent with previous Figures.
being investigated. The first set of them (first approach) corresponds to the
293
transversely isotropic solution given in Table 4. These material parameters
294
have been computed with the optimization procedure detailed in their work
295
using as initial values the results provided by these authors in Ref. [1],
296
i.e. Eq. (104) and the associated principal stretches λ1 (j = 1, ..., m) for
297
the test in the isotropic direction 2. A second set of transversely isotropic
298
material parameters (second approach) are provided in Table 5. In this case,
299
these material parameters have been obtained using as initial data the values
300
w1 = 1/3, r = 1, 2, 3, the stretches λ1 = (λ2 )−1/2 , j = 1, ..., m, and the
301
associated isotropic material parameters previously computed to fit only one
302
of the two stress-stretch distributions in Figure 3 using Eq. (100) of Ref. [1]
303
along with w1 = 1/3 (see Table 6). Again, we deduce from Tables 4 and
TE D
292
(r)
AC C
EP
(j)
(j)
(r)
18
(j)
ACCEPTED MANUSCRIPT
5 that numerical convergence to the values of the associated isotropic model
305
is not attained in any case, even though apparent isotropic predictions are
306
certainly observable in Figure 3. µ1 = 3.22 × 10−4 MPa µ2 = 1.17 × 10−3 MPa µ3 = 4.77 MPa
(1)
w1 = 0.0 (2) w1 (3) w1
= 0.136 = 0.330
α1 = 6.94 α2 = 7.18 α3 = 1.0
RI PT
304
β1 = 1.0 β2 = 11.1 β3 = 1.0
µ1 = 1.26 × 10−3 MPa µ2 = 1.82 × 10−3 MPa µ3 = 4.75 MPa
(1)
M AN U
SC
Table 4: Computed material parameters for the transversely isotropic model of Ref. [1] corresponding to the first approach in Figure 3. Initial guess (associated to a transversely isotropic material) obtained from Ref. [1].
w1 = 0.0 (2) w1 (3) w1
= 0.237 = 0.334
α1 = 6.35 α2 = 6.37 α3 = 1.0
β1 = 1.0 β2 = 12.4 β3 = 1.0
TE D
Table 5: Computed material parameters for the transversely isotropic model of Ref. [1] corresponding to the second approach in Figure 3. Initial guess (associated to an isotropic (j) (j) (r) material) given in Table 6 along with w1 = 1/3 (r = 1, 2, 3) and λ1 = (λ2 )−1/2 (j = 1, ..., m).
µ1 = 1.15 × 10−3 MPa µ2 = 3.80 × 10−3 MPa µ3 = 4.76 MPa
EP
α1 = 6.73 α2 = 6.72 α3 = 1.0
β1 = 1.01 β2 = 13.1 β3 = 1.0
307
AC C
Table 6: Material parameters of the isotropic model of Ref. [1] computed by means of the fitting of the experimental data (only one curve) of Figure 3.
3.3. Model of Holzapfel et al. [18]
308
Other types of strain energy functions widely encountered in the literature
309
are based on the notion of structural tensors and integrity basis, see Ref. 19
ACCEPTED MANUSCRIPT
[27]. For an incompressible transversely isotropic hyperelastic material, the
311
isochoric strain energy function is then formulated in terms of the following
312
four invariants I1 = C : I
I2 =
] 1[ (C : I)2 − C 2 : I 2
313
I5 = C 2 : A0
(17) (18)
SC
I4 = C : A0
RI PT
310
where C is the isochoric Cauchy-Green deformation tensor (i.e. det(C) =
315
I3 = 1) and A0 = a0 ⊗a0 is the structural tensor associated to the anisotropic
316
direction a0 . Hence, in the most general case W = W (I1 , I2 , I4 , I5 ), where
317
coupling terms between the four invariants may be present in the expression
318
of W. In biomechanics, it is convenient to split the strain energy function
319
into a part Wis associated with isotropic deformations and a part Wanis
320
associated with anisotropic deformations [28]. For the case at hand, we write
TE D
M AN U
314
W (I1 , I2 , I4 , I5 ) = Wis (I1 , I2 ) + Wtr (I4 , I5 )
(19)
Proceeding in that way, note that Eq. (12) has been retrieved, in this case
322
formulated in terms of C, i.e.
EP
321
(20)
AC C
W (C, a0 ) = Wis (C) + Wtr (C, a0 )
323
so material-symmetries congruency is guaranteed for these models from a
324
theoretical viewpoint. Eq. (19) is often simplified as (c.f. Ref. [29]) W (I1 , I4 ) = Wis (I1 ) + Wtr (I4 ) 20
(21)
ACCEPTED MANUSCRIPT
According to Refs. [30],[31] the additive, uncoupled, decomposition given in
326
Eq. (21) can be determined from a specific set of two experimental curves
327
as in Ref. [11]. That is, two independent behavior “curves” are used to
328
determine the constitutive functions Wis (I1 ) and Wtr (I4 ).
RI PT
325
In the following example we show that (reduced) models of the form of
330
Eq. (21) may also lack the isotropic material-symmetries congruency from
331
a numerical standpoint if the parameter fitting procedure is not properly
332
designed. As an example of Eq. (21) we take the transversely isotropic
333
counterpart of the well-known model of Holzapfel et al. [18] (see also Ref.
334
[32]), namely W (I1 , I4 ) =
M AN U
SC
329
( ) ) k1 ( c (I1 − 3) + exp k2 (I4 − 1)2 − 1 2 k2
(22)
This model is specifically intended for describing the mechanical behavior
336
of arterial walls, representing a simple starting point to model this type of
337
materials. However, one would expect that if two equal experimental data
338
points representations (obtained from a material tested in different directions,
339
as in the previous example) are used to determine the model constants c, k1
340
and k2 , then k1 should strictly vanish. This is, indeed, what we define as
341
numerical material-symmetries congruency.
EP
TE D
335
Following an analogous procedure to that of the previous example, we
343
have tried to fit the nominal stresses given by the model Eq. (22) to the
344
modified, purely isotropic, experimental data of Figure 1. For isochoric de-
345
formations λ1 λ2 λ3 = 1. Aside, for a uniaxial test performed in direction 1,
346
λ2 = λ3 if the plane 2 − 3 is isotropic, so λ3 = λ1
347
initial isotropic behavior for the initial guess values of the parameters, i.e.
AC C
342
−1/2
21
. We have assumed an
ACCEPTED MANUSCRIPT
348
for example c = 1 MPa and k1 = k2 = 0
(23) −1/2
along with the principal stretches in the “anisotropic” direction λ3 = λ1
350
for the uniaxial test performed over the isotropic direction (denoted herein
351
as axis 1). We note that the first partial derivative of W (I1 , I4 ) with respect
352
to I4 is given by
SC
( ) ∂W (I1 , I4 ) = 2k1 (I4 − 1) exp k2 (I4 − 1)2 ∂I4
RI PT
349
(24)
so the initial specification k2 = 0 gives no evaluation problems in the derived
354
stress expressions. In order to fit the stress predictions of the model to the
355
two point distributions P1 (λ1 ) and P3 (λ3 ) given in Figure 1, the Matlab
356
function lsqcurvefit has been employed. The objective function has been
357
defined as for the previous example. The solution for the material constants
358
that have been computed in this case are
TE D
M AN U
353
c = 1.659 MPa , k1 = 1.024 × 10−4 MPa and k2 = 0.211
361
EP
360
Moreover, as a result of the parameter identification procedure, the ob(j)
tained prediction for every “experimental” point (j) is also such that λ3 ̸= ( )−1/2 (j) λ1 for the uniaxial test performed over the “isotropic” direction.
AC C
359
(25)
362
From these results, it is apparent that the numerical procedure has not con-
363
verged to the corresponding isotropic solution included in the model given in
364
Eq. (22) because we have obtained Wtr (I4 ) ̸= 0. In Figure 4 it can be ob-
365
served the different responses predicted by the model over the two isotropic
22
ACCEPTED MANUSCRIPT
P (MPa) 4
Experimental points P3 (prediction)
RI PT
P1 (prediction)
3
SC
2
0 1
1.5
M AN U
1
2
λ
2.5
TE D
Figure 4: Nominal stresses P3 and P1 in terms of the stretches λ3 and λ1 , respectively, from both uniaxial tests (black solid dots) performed on a transversely isotropic material in a hypothetical limit of isotropic behavior and predictions obtained with the model of Holzapfel et al., Ref [18], initialized with isotropic material constants. Original experimental data is obtained from the test in calendering direction from Diani et al [23].
directions, respectively. We want to emphasize that the fact that k1 is rather
367
small, in relative terms (when compared to the isotropic material constant
368
c), causes that the isotropic behavior is dominant for moderate stretches in
369
both directions (so it is numerically congruent for moderately large strains).
370
However, for larger values of λ the anisotropic contribution becomes domi-
371
nant for the test in the “anisotropic” direction, which causes the response to
372
be different from the one in the transverse isotropic direction.
AC C
EP
366
373
Additionally, we show in Figure 5 the results obtained for the fitting of
374
this model to the two first Piola-Kirchhoff stress distributions of Figure 3.
23
ACCEPTED MANUSCRIPT
P (MPa) Experimental points
3
P3 (prediction)
RI PT
P1 (prediction) 2
0
1
1.5
2
SC
1
λ
2.5
M AN U
Figure 5: Nominal stresses P3 and P1 in terms of the stretches λ3 and λ1 , respectively, from both uniaxial tests (black solid dots) performed on a transversely isotropic material in another hypothetical limit of isotropic behavior. Predictions obtained with the model of Holzapfel et al., Ref [18], initialized with isotropic material constants. Original experimental data is obtained from the test in transverse direction from Diani et al [23].
Numerical convergence to the exact isotropic response of the model (with
376
k1 = 0) is not attained either, although better results are obtained for the
377
range of deformations under study.
TE D
375
Finally, we note that in order to prevent this undesired issue (specially
379
when slightly anisotropic, nearly isotropic materials are characterized with
380
anisotropic strain energy functions) it is possible to design a constrained
381
parameter fitting procedure so as to attain numerical material-symmetries
382
congruency for analytically congruent models.
383
4. Conclusions
AC C
EP
378
384
In this paper we have presented the material-symmetries congruency
385
property for anisotropic hyperelastic materials. This is a desirable material 24
ACCEPTED MANUSCRIPT
property for several reasons. For example, it guarantees that the isotropic be-
387
havior of the matrix is recovered when the volume of reinforcement vanishes,
388
as one should expect. It also guarantees for orthotropic materials that when
389
the volume of one of the fibers vanishes, the transversely isotropic behavior
390
is recovered.
RI PT
386
The material-symmetries congruency property is not only an analytical
392
issue, it is also a numerical one. The lack of numerical material-symmetries
393
congruency in analytically congruent materials is due to the design of the
394
parameter-fitting procedure inherent to the model from composite material
395
experiments. The use of material parameters with a clear physical inter-
396
pretation in the material models easies the development of procedures that
397
preserve the material-symmetries congruency not only from an analytical but
398
also from a numerical point of view.
TE D
Acknowledgements
M AN U
SC
391
Partial financial support for this work has been given by grant DPI201126635 from the Direcci´on General de Proyectos de Investigaci´on of the Min-
AC C
References
EP
isterio de Econom´ıa y Competitividad of Spain.
[1] M. Itskov, N. Aksel. A class of orthotropic and transversely isotropic hyperelastic constitutive models based on a polyconvex strain energy function. International Journal of Solids and Structures, 41, 3833–48, 2004. 25
ACCEPTED MANUSCRIPT
[2] M. Latorre, F.J. Mont´ans. What-You-Prescribe-Is-What-You-Get orthotropic hyperelasticity. Computational Mechanics, 53(6), 1279–1298,
RI PT
2014. [3] G.A. Holzapfel. Nonlinear Solid Mechanics. A Continuum Approach For Engineering. Wiley, Chichester, 2000.
SC
[4] G.A. Holzapfel, R.W. Ogden. Constitutive modelling of passive myocardium. A structurally-based framework for material characterization.
3445–3475, 2009.
M AN U
Philosophical Transactions of the Royal Society of London, Series A, 367,
[5] S. G¨oktepem S.N.S. Acharya, J. Wong, E. Kuhl. Computational modeling of passive myocardium. International Journal for Numerical Methods in Engineering, 27, 1–12, 2011.
TE D
[6] J. Schr¨oder, P. Neff, D. Balzani. A variational stable anisotropic hyperelasticity. International Journal of Solids and Structures, 42, 4352–4371, 2005.
EP
[7] M.P. Nash, P.J. Hunter. Computational mechanics of the heart. Journal of Elasticity, 61, 113–141, 2000.
AC C
[8] D. Tang, C. Yang, T. Geva, R. Rathod, H. Yamauchi, V. Gooty, A. Tang, M.H. Kural, K.L. Billiar, G. Gaudette, P.J. del Nido. A multiphysics modeling approach to develop right ventricle pulmonary valve replacement surgical procedures with a contracting band to improve ventricle ejection fraction. Computers & Structures, 122, 78–87, 2013.
26
ACCEPTED MANUSCRIPT
[9] C.O. Horgan, J.G. Murphy. Simple shearing of soft biological tissues. Proceedings of the Royal Society of London, Series A, 467, 760–777,
RI PT
2010. [10] K. Costa, J. Holmes, A. McCulloch. Modelling cardiac mechanical prop-
ciety of London, Series A, 359, 1233–1250, 2001.
SC
erties in three dimensions. Philosophical Transactions of the Royal So-
[11] T. Gasser, R.W. Ogden, G.A. Holzapfel. Hyperelastic modelling of arterial layers with distributed collagen fibre orientations. Journal of the
M AN U
Royal Society Interface, 3, 13–35, 2006.
[12] H. Schmid, P. O’Callagan, M.P. Nash, W. Lin, I.J. LeGrice, B.H. Smaill, A.A. Young,P.J. Hunter. Myocarcial material parameter estimation. A non-homogeneous finite element study from simple shear tests. Biome-
TE D
chanics and Modeling in Mechanobiology, 7, 161–173, 2008. [13] S. Rossi, T. Lassila, R. Ruiz-Baier, A. Sequeira, A. Quarteroni. Thermodynamically consistent orthotropic activation model capturing ventricular systolic wall thickning in cardiac electromechanics.
EP
European Journal of Mechanics, A/Solids, 2014, in Press. DOI: 10.1016/j.euromechsol.2013.10.009.
AC C
[14] F. Vogel, S. G¨oktepe, P. Steinmann, E.Kuhl. Modeling and simulation of viscous electro-active polymers. European Journal of Mechanics, A/Solids, 2014, In Press. DOI: 10.1016/j.euromechsol.2014.02.001. [15] J. Merodio, R.W. Ogden. Instabilities and loss of ellipticity in fiber-
27
ACCEPTED MANUSCRIPT
reinforced compressible non-linearly elastic solids under plane deformation. International Journal of Solids and Structures, 40, 4707–4727, 2003.
RI PT
[16] J. Merodio, R.W. Ogden. Material instabilities in fiber-reinforced non-
linearly elastic solids under plane deformation. Archives of Mechanics, 54, 525–552, 2002.
SC
[17] J.E. Bischoff, E.M. Arruda, K. Grosh. Finite element simulations of orthotropic hyperelasticity. Finite Elements in Analysis and Design, 38,
M AN U
983–998, 2002.
[18] G.A. Holzapfel, T. Gasser, R.W.Ogden. A new constitutive framework for arterial wall mechanics and a comparative study of material models. Journal of Elasticity, 61, 1–18, 2000.
[19] M. Latorre, F.J. Mont´ans. Extension of the Sussman–Bathe spline-based
TE D
hyperelastic model to incompressible transversely isotropic materials. Computers & Structures, 122, 13–26, 2013. [20] T. Sussman, K.J. Bathe. A Model of Incompressible Isotropic Hy-
EP
perelastic Material Behavior using Spline Interpolations of TensionCompression Test Data. Communications in Numerical Methods in En-
AC C
gineering, 25, 53–63, 2009.
[21] C. Truesdell, W. Noll. The non-linear field theories of mechanics. Handbuch der Physik, Vol. III/3, Springer, Berlin, New York, 2004. [22] R.W. Ogden. Nonlinear Elastic Deformations. Dover, Mineola, New York, 1997. 28
ACCEPTED MANUSCRIPT
[23] J. Diani, M. Brieu, J.M. Vacherand, A. Rezgui. Directional model for isotropic and anisotropic hyperelastic rubber-like materials. Mechanics
RI PT
of Materials, 36, 313–21, 2004. [24] L.L. Demer, F.C.P. Yin. Passive biaxial mechanical properties of isolated canine myocardium. Journal of Physiology, 339, 615–630, 1983.
[25] E.H. Twizell, R.W. Ogden. Non-linear optimization of the material con-
SC
stants in Ogden’s stress-deformation relation for incompressible isotropic elastic materials. J. Aust. Math. Soc. 24, 424–434, 1983.
M AN U
[26] R.W. Ogden, G. Saccomandi, I. Sgura. Fitting hyperelastic models to experimental data. Computational Mechanics, 34(6), 484–502, 2004. [27] A.J.M. Spencer. Constitutive theory for strongly anisotropic solids. In: A.J.M. Spencer (ed.), Continuum Theory of the Mechanics of FibreReinforced Composites, CISM Courses and Lectures No. 282, Inter-
1984.
TE D
national Centre for Mechanical Sciences, Springer-Verlag, Wien, 1–32,
[28] H.W. Weizs¨acker, G.A. Holzapfel, G.W. Desch, K. Pascale. Strain en-
EP
ergy density function for arteries from different topographical sites. Biomediz. Tech., 40, 139–141, 1995.
AC C
[29] J.D. Humphrey. Cardiovascular Solid Mechanics. Cells, Tissues, and Organs. Springer, New York, 2002. [30] J.D. Humphrey,F.C.P. Yin. On constitutive relations and finite deformations of passive cardiac tissue. Part I: A pseudo-strain energy function. Journal of Biomechanical Engineering, 109, 298–304, 1987. 29
ACCEPTED MANUSCRIPT
[31] G.A. Holzapfel, R.W. Ogden. On planar biaxial tests for anisotropic nonlinearly elastic solids. A continuum mechanical framework. Mathe-
RI PT
matics and Mechanics of Solids, 14, 474–489, 2009. [32] G.A. Holzapfel, R.W. Ogden. Constitutive modelling of arteries. Proc. R. Soc. A, 466, 1551–1597, 2010.
SC
[33] K.C.Valanis, R.F.Landel. The strain-energy function of a hyperelastic
material in terms of the extension ratios. Journal of Applied Physics,
AC C
EP
TE D
M AN U
38, 2997–3002, 1967.
30
ACCEPTED MANUSCRIPT Highlights: • •
AC C
EP
TE D
M AN U
SC
RI PT
•
The property of material-symmetries congruency for orthotropic and transversely isotropic hyperelastic materials is presented An anisotropic material model holding this property recovers an isotropic behavior when the experimental data used to calibrate the model corresponds to that of an isotropic material A model may have analytical material-symmetries congruency and still lack numerical one due to the procedure employed to obtain the material parameters from experimental data