- Email: [email protected]
The study of equivalent material parameters in a hyperelastic model
International Journal of Non-Linear Mechanics 89 (2017) 142–150
Contents lists available at ScienceDirect
International Journal of Non–Linear Mechanics journal homepage: www.elsevier.com/locate/nlm
The study of equivalent material parameters in a hyperelastic model c,⁎
c
L. Zanelli , A. Montanaro , E.L. Carniel a b c
a,b
, P.G. Pavan
a,b
, A.N. Natali
a,b
MARK
Department of Industrial Engineering, University of Padova, 35131 Padova, Italy Centre for Mechanics of Biological Materials, University of Padova, 35131 Padova, Italy Department of Mathematics, University of Padova, 35121 Padova, Italy
A R T I C L E I N F O
A B S T R A C T
Keywords: Hyperelasticity Constitutive parameters Polyconvexity Tissue modeling
We study a hyperelastic model of some biological soft tissues with emphasis on the problem of its matching with the material parameters acquired by experimental mechanical tests. First, we study the polyconvexity property of the hyperelastic model. Then, we explore the notion of equivalent sets of material parameters. We perform a numerical study of the regions of equivalent material parameters characterizing the curves predicted by the hyperelastic model that are close, within a prefixed tolerance, to those given by the experimental data. In the numerical study we use the quadratic variation and the Hausdorff distance. The study suggests that a qualitative knowledge of shape and volume of the regions of equivalent material parameters can provide both a criterion for the optimal match between the model with the experimental data and an indication on the reducibility of the number of parameters used in the model.
1. Introduction 1.1. An overview of soft tissue modeling Many biological soft tissues can be described by hyperelastic models, viscoelastic models, or even by more general models (see, e.g., [6,9,13,26,28,29] and their references). As a consequence, a lot of experiments can be numerically simulated and the related material parameters implemented in the model describing tissue properties can be adjusted in order to recover the experimental data. Nevertheless, there are several open issues about inverse approaches in hyperelasticity, as shown in the literature (see, e.g., [1] and references therein). We first point out that an inverse problem is usually related to a model which can be defined by different possible families of material parameters. Hence, once the experimental data and a related model have been fixed, one must adopt an optimal approach in order to localize the material parameters. The first general domain of parameters must provide physically reasonable material behavior of the model. It is well known that for hyperelastic models such a behavior is guaranteed by the polyconvexity condition of the strain energy function with respect to the deformation gradient (see, e.g., [7,20,25,2,3,26], and their references). The second more specific domain of parameters must match a known family of stress-strain curves given by the experimental data obtained for the specific tissue. With respect to this issue, various approaches can be given. Often a cost function and a suitable algorithm are employed to
⁎
evaluate and minimize the differences between experimental data and model simulations. For example, the particular simulated annealing algorithm, developed in [8], has been proved suitable for the complex behavior of the cost function. Optimization algorithms must be further enhanced because of the multimodal behavior of the cost function. Thus, coupled deterministic - stochastic algorithms have been developed (see for example [17,18] and the references therein). The approach allows both to minimize the discrepancy between experimental and model results and to explore the domain of admissible parameters. In another type of approach an easy fitting of polyconvex stored energies can be applied to soft tissues with no optimization procedure. This is shown for example in [3]. We stress that in both approaches the final target is to find at least one vector of constitutive parameters for one given tissue, thus avoiding a complete analysis of the model and its dependence from the parameters. 1.2. On the determination of material parameters In view of the different determination procedures of parameters outlined above, we stress that another related meaningful issue is that of uniqueness of the constitutive parameters for a given tissue. In fact, a given set of experiments on a soft tissue can provide only a certain number of limited information about the mechanical behavior of a material. The vector of the material parameters might be modified by increasing the set of the known experiments. This leads to the issue of
Corresponding author. E-mail address: [email protected] (L. Zanelli).
http://dx.doi.org/10.1016/j.ijnonlinmec.2016.12.014 Received 29 March 2016; Received in revised form 22 December 2016; Accepted 22 December 2016 Available online 26 December 2016 0020-7462/ © 2016 Elsevier Ltd. All rights reserved.
International Journal of Non-Linear Mechanics 89 (2017) 142–150
L. Zanelli et al.
dH (X , Y )≔max{d (X , Y ); d (Y , X )} d (X , Y )≔sup inf d (x, y ).
the non-uniqueness of the parameters inferred by the model. This is an open problem for many hyperelastic models of tissues, and it seems that there is not a general solution to it (see, e.g. the different approaches in [1,10,12,14,15,22,26,27] and the references therein). We also remark that increasing the number of experimental data could be useful to recover uniqueness of the material parameters, but this feature seems to be connected with the particular choice of the model. In some papers, the issue of non-uniqueness of the material parameters does not play any role because their goal is to represent the characteristic stress-strain curve for a given restricted set of experiments (see, e.g., Sect. 4.1 in [2]). For the hyperelastic model presented here, we consider the problems of the determination of the material parameters and also of their uniqueness (see Section 3.1). A fully incompressible hyperelastic model and an almost incompressible hyperelastic model are used to identify the constitutive parameters on different kind of tissues. This approach can be extended also to other hyperelastic models, viscoelastic models, and also to more general models.
x∈X y∈Y
It is used to compare the graphs of the curves relating the experimental data to the model. It guarantees a satisfactory analysis of the tissues. Now, let (ω) ⊂ 2 be the graph of the curve given by P11 (F, ω) for F = Fj and fixed ω, and let ⊂ 2 be the graph of the experimental curve with values yj; then the inequality ∼), ) ≤ d ( (ω ), ) d ( (ω (6) H
In what follows, we provide an overview of the main results of the paper with respect to the above discussed topics. In Section 3.1 we introduce the hyperelastic model with the strain energy function, related invariants, and constitutive parameters. In Section 3.2 we prove the property of polyconvexity of the strain energy function W = W (C, ω), where C = FT F is the right Cauchy-Green strain tensor, F is the deformation gradient and ω ∈ N is the vector whose components are the constitutive parameters (here N = 2, 4, 6 ). It turns out that W is twice continuously - differentiable with respect to C and continuous with respect to ω. It is well known that polyconvexity of the energy guarantees the physical behavior of the model. In our model such property holds when the set of constitutive parameters Ω ⊂ N has the form
2. Settings and preliminaries In this section we provide a summary of notations used in the paper together with some central definitions of hyperelasticity. We denote by Lin+ the set of all second-order tensors with positive determinant, to be identified with the family of 3×3 matrices n × n with positive determinant. The set Orth+ is the subset of Lin+ given by rigid rotation tensors R , namely such that RRT = I , where I is the unit tensor. The set Sym+ is the subset of Lin+ of symmetric tensors U , i.e. such that U = UT . Following the notations of continuum mechanics, here F ∈ Lin+ denotes the deformation gradient, C = FT F ∈ Sym+ denotes the right Cauchy-Green deformation tensor, the map F ∈ Lin+⟼W (FT F) ∈ is the strain energy function, and the first Piola-Kirchhoff stress tensor reads P = 2F∂CW . We refer the reader to the standard textbooks of continuum mechanics [13,23,24,28].
(1)
Thus within this domain we look for the regions of parameters that guarantee the match between the model (2)
P (F, ω) = 2F ∂CW (F, ω).
(first Piola-Kirchhoff tensor) and the experimental data. As it is known, in the fully incompressible case we have P = −p F−T + 2F ∂CW where p is a Lagrange multiplier. For a given family of deformation gradients {Fj : 1 ≤ j ≤ M}, let yj be the mean value of the experimental data set caused by the deformation Fj , which is taken on a given specimen. Then we look for a vector ω that minimizes the cost function
σ (ω, y )≔
1 M
M
∑
2−
j =1
yj P11 (Fj , ω)
−
P11 (Fj , ω) yj
Definition 2.1. A map F ∈ Lin+⟼W (F) ∈ is said to be convex if
2
W (τ F1 + (1 − τ ) F2) ≤ τW (F1) + (1 − τ ) W (F2) (3)
(7)
for every F1, F2 and τ with 0 ≤ τ ≤ 1. In the following, we recall some generalized convexity conditions. In particular, from [23] we recall
within a prefixed tolerance. This is done by a coupled stochastic-deterministic algorithm of the kind briefly discussed above and used for similar models in, e.g., [7,20] and in the references therein. In particular, a (stochastic) simulated annealing algorithm to find a first approximated local minimum. Then, we apply a (deterministic) simplex method to find a local minimum. After that, we need to study all the possible other material ∼ which are, in a fixed bounded region of Ω, equivalent parameters ω ∼ in such a way that to ω . We first select ω
∼, y ) ≤ σ (ω , y ). σ (ω
H
∼. As we will show in determines the second selection of the parameters ω Section 3.3, various two dimensional regions of parameters can be displayed by looking at different material parameters satisfying (4) and ∼ are minimizers of the cost (6). As explained above, the parameters ω function σ as well as the parameter ω . However, in order to determine ∼ we need first to determine at least a parameter the set of parameters ω ω and then we can write inequalities (4)–(6) and call a related routine ∼. This numerical study is that determines numerically the set of ω ® obtained through MATLAB (software house MathWorks). The related algorithm can be easily adjusted for different stress tensors and thus for different hyperelastic models. A qualitative analysis of such regions of parameters, for different experimental data, provides useful information about the mechanical properties of the tissue under study and also some indicators of the optimality of the chosen hyperelastic model. The novelty of our study is to show the possibility of a more complete analysis of the material parameters of a soft tissue (here in particular urethral tissue). This involves not only the determination of a (vector of) material parameter ω as done previously, but in addition the determination of the larger set of parameters providing the same minimization of cost function involving the model and the experiment.
1.3. Outline of the results
Ω = {(ω1, …, ωN ) = ω ∈ N ωi > 0}.
(5)
Definition 2.2. A map F ∈ Lin+⟼W (F) ∈ is said to be polyconvex if there exists a function P: 3×3 × 3×3 × → such that
W (F) = P (F, Adj[F], det[F])
(8)
∼ ∼ ∼ ∼ and (X͠ , Y , Z ) ∈ 19⟼P (X͠ , Y , Z ) ∈ is convex. The policonvexity property is usually used to select physically reasonable models (see, e.g., [26] and the references therein).
(4)
Definition 2.3. A twice differentiable function F ∈ Lin+⟼W (F) ∈ fulfills the Legendre-Hadamard condition if ∀ a, b ∈ 3, ∀ F ∈ Lin+ we have
In addition, as a second requirement we make use of the so-called Hausdorff distance dH between sets (see for example [4] and the references therein). Let us recall that the Hausdorff distance between two sets X , Y ⊂ 2 with Euclidean distance on 2 , is defined by
D2F W (F). (a ⊗ b, a ⊗ b ) ≥ 0. 143
(9)
International Journal of Non-Linear Mechanics 89 (2017) 142–150
L. Zanelli et al.
Because of the uniaxial configuration of the experimental tests, the constraints on the second and third diagonal terms P22 = P33 = 0 can be assumed. In the different settings for P as described above, these constrains provide the functional dependences λ22 = λ22 (λ11) and λ33 = λ33 (λ11). The numerical approach we used to get such two functions is obtained by the well known Newton-Raphson method (see for example [11] for the implementation by the use of MATLAB® ). With regard to the fully incompressible model, we recall that J = det(F) = 1 and
When the strain energy function W is twice continuously differentiable, the polyconvexity condition implies the weaker LegendreHadamard condition, which also guarantees a physically reasonable model (see, e.g., [2] and the references therein). 3. Main results 3.1. The hyperelastic model
P = −p F−T + 2F ∂CW
The strain energy function modeling the tissue that is considered here reads as
W (C) = Wm (C) + Wf (C, a 0 ⊗ a 0).
where p is a Lagrange multiplier determined by the boundary conditions. We stress that the aim of our paper is not to show that our energy functions in (12)–(13)–(14) are the optimal ones to modelize the urethral tissues under study with respect to other models used in the past. We recall that our primary target is to show that a complete analysis of the material parameters of a soft tissue (here in particular urethral tissue) should involve not only a single (vector of) material parameter ω but the larger set providing the same minimization of cost functions involving the model and the experiment. The analysis provided in this paper can be done with other energy functions and then one can compare the determination of the regions of material parameters with the same physical meaning.
(10)
This function W describes the behavior of a fibrous tissue with a family of fibers oriented along the local direction a 0 . (about strain energies for other kinds of tissue see, e.g., [5,7,16,20,21].) In our experiments, the tensile direction will be transversal or longitudinal. In particular, the first term in (10) reads
∼ Wm (C) = Um (J ) + W͠ m (I1)
(11)
with
Um (J ) =
Kv [(J − 1)2 + J −r + rJ − (r + 1)] 2 + r (r + 1)
C ∼ ∼ W͠ m (I1) = 1 {exp[α1 (I1 − 3)] − 1} α1
(12)
3.2. Polyconvexity
(13)
In this section, we provide the proof of the polyconvexity of the strain energy function W (C) as in (10) with respect to the deformation gradient F . For the sake of simplicity, we divide the computations into two parts: (A) Let us prove the polyconvexity of the map F⟼Wm (FT F) for Wm given by (11). As for the first term Um (J ), we recall that J = det(F). By recalling the definition of polyconvexity (see Definition 2.2) we now devote our attention to the one dimensional map J ⟼Um (J ) ∈ which a smooth function (i.e. C ∞) for J > 0 , and thus we try to prove its convexity property. The one dimensional setting ensures the convexity by the study of the second order derivative.
where J ≔(det C)1/2 = det F is the deformation Jacobian and ∼ I1 = tr(J −2/3C) is the first invariant of the isovolumetric part of C . The volumetric contribution Um of the strain energy function must be assumed to contribute in the almost incompressible case. The constitutive parameters Kv and r account for matrix compressibility, while C1 and α1 account for shear behavior. The second term in (10) reads as
Wf (C, a 0 ⊗ a 0) = Wf (I4 ) =
C4 {exp[α4 (I4 − 1)] − α4 (I4 − 1) − 1} (α4 )2 (14)
where I4 is the structural invariant defined by I4 = C: (a 0 ⊗ a 0). The parameters C4 and α4 account for the initial stiffness and increase of fibers stiffness by stretch, respectively. The decomposition (10) of the strain energy implies for the first Piola-Kirchhoff stress tensor
U ″ m (J ) =
Kv [2 + r (r + 1) J −r −2]. 2 + r (r + 1)
(21)
Thus, for r > 0 and Kv > 0 it follows
U ″m (J ) > 0,
(22)
∀ J > 0.
(15)
P = 2F ∂C W
This inequality ensures that J ⟼Um (J ) is a convex map on the domain (0, +∞). Hence, the composition of Um with J = det(F) implies that F⟼Um (det(F)) is polyconvex. ∼ As for the second term W͠ m (I1) in (11), simply notice that it is the composition (up to a constant) of the exponential (which is convex) and ∼ ∼ of the map I1 ⟼α1 (I1 − 3). More precisely,
the decomposition
P = Pm + Pf ,
(16)
where by (11)
Kv [2J (J − 1) − rJ−r + rJ ] F−T 2 + r (r + 1) 2∼ ∼ + C1 exp[α1 (I1 − 3)](2J −2/3F − I1 F−T ) 3
Pm = 2F∂CWm =
∼ I1 = tr(J −2/3F TF) = J −2/3tr(F TF) = (det(F))−2/3 ∥ F ∥2 .
C4 {exp[α4 (I4 − 1)] − 1} F (a 0 ⊗ a 0). α4
(17)
(18)
In the present paper we consider experimental data for uniaxial tensil tests along direction 1. A diagonal form of the deformation gradient can be assumed:
⎛ λ11 0 0 ⎞ ⎜ ⎟ F = ⎜ 0 λ22 0 ⎟ . ⎝ 0 0 λ33 ⎠
(23)
∼ Now we recall Lemma C.1 in [25], so that F⟼I1 (F) is polyconvex. By the sum of Um and W͠ m we can therefore conclude that F⟼Wm (FT F) is polyconvex. (B) Let us prove the convexity of the map F⟼Wf (FT F, a 0 ⊗ a 0) for Wf given by (14). This is a function of the structural invariant I4 = FT F: (a 0 ⊗ a 0) and in more details
and by (14)
Pf = 2F∂CWf = 2
(20)
I4≔tr((FT F)T (a 0 ⊗ a 0)) = tr((FT F)(a 0 ⊗ a 0))
FT F .
(24)
By recalling Lemma C.2 thanks to the symmetry of the matrix (point 1) in [25] we have that F⟼I4 (F) is convex. Now, it can be easily checked that I4 ⟼Wf (I4 ) ∈ is a smooth function and that the second order derivative reads
(19) 144
International Journal of Non-Linear Mechanics 89 (2017) 142–150
L. Zanelli et al.
Wf (I4 )″ =
C4 {[α4]2 exp[α4 (I4 − 1)]}. (α4 )2
(25)
For C4 > 0 , α4 > 0 we recover
Wf (I4 )″ = C4 exp[α4 (I4 − 1)] > 0,
∀ I4 ∈ ,
(26)
and hence also the convexity condition. Recalling Lemma B.9 in [25], we are now in the position to make the composition of the monotone increasing (when I4 ≥ 1) and convex function Wf with I4 as a function of F , and to get the convexity of F⟼Wf (FT F, a 0 ⊗ a 0). In the case 0 < I4 < 1, we first notice that if I4 is sufficiently closed to 1, then D2F Wf ≃ C4 D2F I4 which is semipositive definite, and whence Wf is convex. For the arbitrary case 0 < I4 < 1 we can always make a rescaling δ (I4 − 1) by a small δ > 0 so that D2F Wf becomes semipositive definite. Notice that such a transformation can be viewed as a rescaling of the parameter α4 > 0 , which remains positive; whence we can say that the hyperelastic model is not really changed and the match between the model and the experimental data involves such a rescaling of the parameter α4.
Fig. 1. Tensile test for distal urethral specimen comparison of experimental data and fully incompressible model results.
Remark 3.1. In view of the above computations, the constitutive parameters satisfying Kv > 0 , r > 0 , C1 > 0 , α1 > 0 , C4 > 0 , α4 > 0 , guarantee physical type responses for our hyperelastic model. As we will see in the next section, for different stress-strain curves, the involved region of equivalent material parameters will be bounded, namely there is not degeneracy of unbounded parameters.
ω = (C1, α1) = (1.09·10−3, 1.26). The two stress-strain curves coming from the experimental data (the mean values) and from the values of Pm,11 evaluated at ω are given in Fig. 1 Related to this case, Fig. 1 shows the two-dimensional region of equivalent material parameters with respect to (C1, α1): The above figure is determined by the check of the constraints (4)– (6); namely we look if they are fulfilled (or not) on the square net of points centered on ω = (C1, α1) given by ± 10 points with respect to the first and second coordinate, and length 0.5 between them. The same procedure is done for the other pairs of parameters and figures in the next part of the paper. We first notice that the region of parameters shown in Fig. 2 is connected, and no other disjoint regions appear. Physically, this means we have the uniqueness of parameters and that any rescaling of (C1, α1) changes the volume of the region. Moreover, we notice that the geometry of the region has a symmetry around a line. This fact suggests the introduction of a transformation given by anticlockwise rotation with center at (1.09·10−3, 1.26) and angle ϕ = 46°, whence from (C1, α1) 1, α1). We can say that α1 is varying slowly with respect to a new pair (C 1, since we can also introduce a further transformation given by an to C (arbitrary) rescaling of α1. Thus, it can be recognized that the model is more sensible with respect to one parameter with respect to another one. In this two-parameter model, such a behavior is expected and other arguments can be used to recover the same observation. We stress that when the parameters space has a dimension greater than two, then such an analysis can be more difficult to recover by a direct study of the
3.3. Domains of material parameters Our paper involves the analysis of experimental data from mechanical tests developed on urethral specimens (see [19]). Urethral tissues show different conformation depending on the specific location. Proximal urethra tissues are composed of an almost isotropic connective matrix reinforced by muscular fibers aligned along longitudinal direction. On the other hand, distal urethral tissues are composed by almost isotropic connective matrix. The mechanical characterization of urethral tissues took into account tensile tests performed on proximal and distal samples, according to both longitudinal and transversal directions. With specific regard to distal samples, results from tests along each different directions did not show statistical differences, enforcing the assumption of isotropic behavior. It follows that we provide six data sets given by experimental stress-strain curves. As we will describe in this section, for any of these experimental stress-strain curves, we provide the interpolation given by the hyperelastic model described in the previous section and a vector of material parameters. Later, we show that for any of these fixed vectors of material parameter we are able to localize all the equivalent material parameters within a (large) bounded region. For any of these figures representing regions of parameters we provide a qualitative analysis which leads to the possible reduction in the number of parameters, the determination of leading parameters, the study of the behavior of equivalent parameters with respect to data obtained by transversal or longitudinal tensils on specimens. Morever, the study of uniqueness or non-uniqueness of material parameters is given by the determination of non-connected regions of equivalent parameters; this approach is not affected by the rescaling of the parameters and whence is well posed, and allows to adjust for any type of material parameter the correct scale to represent its value. Let us begin by looking at the case of a distal urethra specimen by the fully incompressible model. The first Piola-Kirchhoff stress tensor we use here is given by (20) for the energy ∼ W = W͠ m (I1) (27) as given by (13) and with the material parameters (C1, α1). In this case, the minimization of the cost function (3) provides the parameters
Fig. 2. Graph of equivalent (C1, α1) given by computations over a net of 10×10 points with length 0.5·10−3 and 0.5 centered on (C1, α1) = (1.09·10−3, 1.2) .
145
International Journal of Non-Linear Mechanics 89 (2017) 142–150
L. Zanelli et al.
Fig. 3. Family of level sets of the parameters (C1, α1) for distal urethral specimen.
stress tensor and maybe also unexpected. The problem of the reduction of the number of parameters can thus be faced by using the above qualitative analysis of regions of equivalent parameters. Next, we provide a diagram involving a family of level sets (and neighborhoods of them) for the cost function dH in (6) involving (C1, α1) which are the parameters shown in Fig. 2. In this case, the other cost function σ in (4) is constant on the level sets of dH , since here the constraint dH is stronger than σ. In Fig. 3 we see the boundary of region shaped in the Fig. 2 and in addition the higher level sets of the cost function. This shows a unique local (and hence global) minimum of the cost function and implies that the region of parameters around this minimum is connected, namely no other disjoint regions of parameters appears. Similar work can be done for the other cases discussed in the paper. We now look at the proximal urethra specimen by the fully incompressible model, with energy function ∼ W = W͠ m (I1) + Wf (I4 ) (28)
Fig. 4. (a) - (b) Tensile test along transversal and longitudinal direction for proximal urethra specimen by the fully incompressible model.
see (13) and (14). The stress tensor of the fully incompressible model we use here is (20), and the material parameters are (C1, α1, C4, α4 ). The minimization of the cost function in (3) provides the parameters (C1, α1) = (1.12·10−3, 2.24) and (C4, α4 ) = (1.49·10−3, 5.16·10−1). The stress-strain curves are shown in Fig. 4 With respect to Fig. 4, we underline that the parameters space has dimension 4, whence the number of all the possible pairs of parameters is 6. In Fig. 5 we choose to display only the pairs (C1, α1) and (C4, α4 ). The two regions of equivalent parameters linked to Fig. 4 - (a) are given by Figs. 5–6. We first remark that in this case also (C4, α4 ) are varying in a square net of points centered on (1.49·10−3, 5.16·10−1) with length 0.5. Thus, Fig. 5 is the projection of the 4-dimensional set of all the equivalent parameters (C1, α1, C4, α4 ) into the subspace of the parameters (C1, α1). Notice that the region of points in the Fig. 5 is connected. Whence, also in this case we recover uniqueness of parameters up to a rescaling of the volume of this domain. We stress also that the region is strictly included with respect to the net of points where we made the computations. In fact, a possible nonconnected region which is far away can be localized by the same algorithm. We also observe that the above region can be qualitatively described by a sublevel set given by an inequality of type C1 ≤ f (α1) where f : + → + is some decreasing function, which can be easily interpolated thanks to the boundary of the domain in Fig. 5. In view of the above observations, we think that a possible general characterization of such sublevels in connection to some important features of the strain energy function is an interesting future target. These could be for example the choice of the invariants and the asymptotic behavior of the energy with respect to the parameters. We also stress that the determination of the above inequality from the analytical viewpoint without passing through the numerical study can
Fig. 5. Graph of equivalent (C1, α1) given by computations over a net of 10×10 points centered on (C1, α1) = (1.12·10−3, 2.24) with length 0.5·10−3 and 0,5.
be a rather non-trivial task. The region involving the other pair of parameters (C4, α4 ) is shown in Fig. 6. We remark that here we fix the parameters (C1, α1) and thus it is a different projection with respect to the previous one we used into the Fig. 5. Indeed, this is done because here we are looking at the projection in the space of (C4, α4 ) of a subset (i.e. with fixed (C1, α1)) of the 4dimensional domain of all the equivalent parameters which we are taken into account to recover Fig. 5. We see that the above domain in Fig. 6 is contained in the net of 10×10 square region of points, and this suggest that there not exist others domains to detect, whence the 146
International Journal of Non-Linear Mechanics 89 (2017) 142–150
L. Zanelli et al.
Fig. 8. Graph of equivalent (C4, α4 ) for fixed (C1, α1) given by computations over a net of 10×10 points centered on (C4, α4 ) = (1.49·10−3, 5.16·10−1) with length 0.5·10−3 and 0.5.
Fig. 6. Graph of equivalent (C4, α4 ) for fixed (C1, α1) given by computations over a net of 10×10 points centered on (C4, α4 ) = (1.49·10−3, 5.16·10−1) with length 0.5 × 10−3 and 0.5.
uniqueness of parameters up to rescalings of a single connected region. Obviously, different projections of the type in Fig. 6 can be given, by fixing any other equivalent pair of parameters (C1, α1) showed in Fig. 5, but analogous information are recovered. In Fig. 7 we provide the two regions of parameters linked to Fig. 4 (b). If we compare Figs. 7 and 5 we notice a very similar shape, but with a slight increase of the volume of the parameters. This fact suggests to consider the smallest region of parameters in Fig. 5 as the characterizing material parameters. By looking at Figs. 6 and 8 we notice a relevant change of the domain, which is in particular a subset of the previous one. Such a fact can be considered as an optimal feature of the model, namely that different experimental data on the same tissue allows to localize a subset of material parameters. This behavior allows to consider the smallest one as the characterizing material parameters. We underline the usefulness of this powerful approach which make use of transversal and longitudinal data to get characterizing material parameters. Let us now discuss the results on the distal urethral specimen for tensil in the longitudinal direction with the almost incompressible model. The energy function is
∼ W = Um (J ) + W͠ m (I1)
Fig. 9. Tensile test for distal urethral specimen, comparison of experimental data along longitudinal direction and almost incompressible model results.
parameters obtained by the minimization of the cost function are (K , r ) = (3.75, 2.77·10) and (C1, α1) = (1.09·10−3, 1.26). The stress-strain curves are given by Fig. 9. In what follows we consider all the 4-dimensional set of equivalent parameters (K , r , C1, α1) and we first display the projection into the 2dimensional space of the pair (K,r) Fig. 10. We stress that the region of the parameters is not completely
(29)
where the two terms on the righthand side are given by expressions (12) and (13). The first Piola-Kirchhoff stress tensor we use is Pm as defined in (17) where the involved parameters are (K,r) and (C1, α1). The
Fig. 7. Graph of equivalent (C1, α1) given by computations over a net of 10×10 points centered on (C1, α1) = (1.12·10−3, 2.24) with length 0.5·10−3 and 0.5.
Fig. 10. Graph of equivalent (K, r) given by computations over a net of 5×5 points centered on (K , r ) = (3.75, 2.77·10) with length 1.
147
International Journal of Non-Linear Mechanics 89 (2017) 142–150
L. Zanelli et al.
Fig. 11. Graph of equivalent (C1, α1) for fixed (K, r) given by computations over a net of 5×5 points centered on (C1, α1) = (1.09·10−3, 1.26) with length 10−3 and 1.
contained in the net of points we used for the comparison of the parameters. However, some simple arguments (which involve the asymptotics of Pm with respect to the parameters) show that the parameters must be bounded for a fixed stress-strain curve. About the other related couple of material parameters, we fix (K,r) and we provide the region of equivalent (C1, α1) in Fig. 11. In this case we have two non-connected components. However, a more complete description involves the projected parameters (C1, α1) for all the equivalent (K,r) given by the previous figure. This analysis provides, also in this case, a connected region of parameters. This local description is meaningful since it shows that not connected components can appear. The last case we discuss is about proximal urethral specimens in the almost incompressible setting for transversal and longitudinal tensils. ∼ W = Um (J ) + W͠ m (I1) + Wf (I4 ) (30)
Fig. 12. (a) - (b) Tensile test along transversal and longitudinal direction for proximal urethra specimen by the almost incompressible model.
where the three terms on the righthand side are given by (12), (13) and (14). The first Piola-Kirchhoff stress tensor here used is Pm + Pf as defined in (17)–(18) for parameters (K,r), (C1, α1) and (C4, α4 ). In particular, the parameters ω minimizing the cost function are (K , r ) = (6.45, 2.60) , (C1, α1) = (1.12·10−3, 2.24) and (C4, α4 ) = (1.49·10−3, 5.16·10−1). The stress-strain curves are given in Fig. 12. Since in this case the number of the parameter is 6, and whence larger than in the previous cases, we need to take into account the computational cost of the algorithm. Indeed, the comparison of 106 vector of parameters (instead of 102 or 104) implies a meaningful increase of the time for the computation. For this reason, here we are going to show only three pair of parameters when we have fixed the other four according to the vector of parameters ω . In Fig. 13 we provide the equivalent parameters of the transversal data sets shown into the above Fig. 12 - (a). In Fig. 13 we notice a large set of parameters which can be qualitatively described by inequalities of type K < g (r ) for an increasing function g: + → +. Also in this case, we can provide an interpolation of the boundary of this region in order to get an analytical description of such a function. We also underline that the above graph shows a slow dependence from r, and as a consequence this fact suggest to introduce a large rescaling of the energy W with respect to r. The second and the third pairs of equivalent parameters exhibit similar behavior when compared each other. The second pairs is given in Fig. 14, which shows that The behavior of α1 is rapidly decreasing when C1 is increasing. Thus, we can say that this hyperelastic model has a sensibility with respect to small variations of C1. Such a dependence with respect to its material parameter can be considered optimal, since it allows to directly localize small regions of parameters. In Fig. 14 we notice also a region of type
Fig. 13. Graph of equivalent (K, r) for fixed (C1, α1) and (C4, α4 ) given by computations over a net of 10×10 points centered on (K , r ) = (6.45, 2.60) with length 1.
C1 ≤ h1 (α1) for some decreasing function h1: + → + . The possible diverging behavior of α1 for very small values of C1 suggests to consider the region of parameters 0 < C1, min ≤ C1 for some small fixed material parameter C1, min . The same arguments used above here for Fig. 15 suggests the consideration of a lower bound of type 0 < C4, min ≤ C4 and inequality C4 ≤ h4 (α4 ) for some decreasing function h4 : + → + . With respect to the above experimental data in Fig. 12 - (b), the fist pair of material parameters are shown in Fig. 16. Notice that Fig. 16 is the same as Fig. 13 and hence the longitudinal (or transversal) tensile on this tissue does not change the determination of the material parameters (K,r). The second pair of parameters linked 148
International Journal of Non-Linear Mechanics 89 (2017) 142–150
L. Zanelli et al.
Fig. 14. Graph of equivalent (C1, α1) for fixed (K, r) and (C4, α4 ) given by computations over a net of 10×10 points centered on (C1, α1) = (1.12·10−3, 2.2) with length 10−3, 1.
Fig. 17. Graph of equivalent (C1, α1) for fixed (K, r) and (C4, α4 ) , given by computations over a net of 10×10 points centered on (C1, α1) = (1.12·10−3, 2.24) with length 10−3, 1.
Fig. 15. Graph of equivalent (C4, α4 ) for fixed (K, r) and (C1, α1) , given by computations over a net of 10×10 points centered on (C4, α4 ) = (1.49·10−3, 5.16·10−1) with length 10−3 and 1.
Fig. 18. Graph of equivalent (C4, α4 ) for fixed (K, r) and (C1, α1) , given by computations over a net of 10×10 points centered on (C4, α4 ) = (1.49·10−3, 5.16·10−1) with length 1.
In contrast to the previous two figures, here we stress that the Fig. 18 is different from Fig. 15, and in particular it is a small subset. This suggest to consider the points in Fig. 18 (and small domains around them which can be seen at a different scale) as the third couple of material parameters characterizing the tissue. 4. Conclusions In this paper we have provided a numerical study of the regions of equivalent material parameters for urethra type specimens. For any tissue of such kind, we displayed the curves coming from the hyperelastic model and we have shown the correspondence with the curves coming from experimental data up to some small values for the quadratic variation and the Hausdorff distance function. We stress that we used the Hausdorff distance function in order to recover a better match between model and experiments with respect to the one obtained only by the quadratic variation. We have shown that a qualitative study of the shape and the volume of such regions of parameters provides a criterion for the optimal match between the model and the experimental data, together with a check on the possible reduction of the number of parameters involved in the model. Indeed, for any of the figures involving the regions of parameters we provide a qualitative analysis which leads to the possible reduction of the number of parameters related to the determination of leading parameters, the study of the behavior of equivalent parameters with respect to data obtained by transversal or longitudinal tensils on specimens. Moreover, the study of uniqueness or non-uniqueness of material
Fig. 16. Graph of equivalent (K, r) for fixed (C1, α1) and (C4, α4 ) , given by computations over a net of 10×10 points centered on (K , r ) = (6.45, 2.60) with length 1.
to Fig. 12 - (b) is displayed in Fig. 17. As in the previous case, we stress that Figs. 14 and 17 are coinciding. We interpret this fact as a sort of stability result also for (C1, α1) which ensures an intrinsic character for this region of material parameters. We recall that a rescaling of this region is possible, whence a determination of a small volume just related to a different scale. Fig. 18 describes the third pair of material parameters. 149
International Journal of Non-Linear Mechanics 89 (2017) 142–150
L. Zanelli et al.
parameters is given by the determination of non-connected regions of equivalent parameters. Data from tensile tests only have been analyzed to identify the parameters of both fully incompressible and almost incompressible formulations. Physically, tensile tests do not provide information about the compressibility of the material and the proposed algorithm returns the known uniqueness of the parameters. Furthermore, a notion of stability of the hyperelastic model is given through the analysis of the behavior of the regions of equivalent parameters when the transversal and longitudinal tensils are considered on the same tissue. In order to complete the remarks on the paper, we underline that the approach here described can be used for any hyperelastic model, and that the features here we discussed can be used to select an optimal model for any given experimental data set. The proposed approach allows to evaluate relationships between model configuration and independent sets of experimental in order to obtain uniqueness of the parameters. Further studies can investigate the influence of more complex experimental situations, as biaxial or shear loading.
[9] [10]
[11]
[12]
[13] [14]
[15] [16]
[17] [18] [19]
References [1] S. Avril, S. Evans, K. Miller, Inverse problems and material identification in tissue biomechanics, Spec. Issue J. Mech. Behav. Biomed. Mater. 27 (November) (2013). [2] D. Balzani, J. Schröder, P. Neff, G.A. Holzapfel, A variational approach for materially stable anisotropic hyperelasticity, Int. J. Solids Struct. 42 (2005) 4352–4371. [3] D. Balzani, J. Schröder, P. Neff, G.A. Holzapfel, A Polyconvex framework for soft biological tissues. Adjustment to experimental data, Int. J. Solids Struct. 43 (2006) 6052–6070. [4] E. Belogay, C. Cabrelli, U. Molter, R. Shonkwiler, Calculating the hausdorff distance between curves, Inf. Process. Lett. 64 (1) (1997) 17–22. [5] J. Bonet, A.J. Burton, A simple orthotropic, transversely isotropic hyperelastic constitutive equation for large strain computation, Comput. Methods Appl. Mech. Eng. 162 (1–4) (1998) 161–164. [6] R. Cai, F. Holweck, Z. Feng, F. Peyraut, A new hyperelastic model for anisotropic hyperelastic materials with one fiber family, Int. J. Solids Struct. 84 (1) (2016) 1–16. [7] E.L. Carniel, V. Gramigna, C.G. Fontanella, C. Stefanini, A.N. Natali, Constitutive formulations for the mechanical investigation of colonic tissues, J. Biomed. Mater. Res A. 102 (5) (2014) 1243–1254. http://dx.doi.org/10.1002/jbm.a.34787 (Epub 2013 May 24). [8] A. Corana, N. Marchesi, C. Martini, S. Ridella, Minimizing multimodal functions of
[20]
[21]
[22] [23] [24] [25] [26] [27] [28] [29]
150
continuous variables with the Simulated Annealing algorithm, ACM Trans. Math. Softw. 13 (3) (1987) 262–280. S.C. Cowin, S.B. Doty, Tissue Mechanics, Springer Science and Business Media, 2007. L. Dubuis, S. Avril, J. Debayle, P. Badel, Identification of the material parameters of soft tissues in the compressed leg, Comput. Methods Biomech. Biomed. Eng. 15 (1) (2012) 3–11. W. Gander, M.J. Gander, F. Kwok, Scientific Computing - An Introduction using Maple and MATLAB, Texts in Computational Science and Engineering, 2014th Edition, Springer. M. Geerligs, L.C.A. van Breemen, G.W.M. Peters, P.A.J. Ackermans, F.P.T. Baaijens, C.W.J. Oomens, In vitro indentation to determine the mechanical properties of epidermis, J. Biomech. 44 (2011) 1176–1181. G.A. Holzapfel, Nonlinear Solid Mechanics, Wiley, New York, 2000. G.A. Holzapfel, T.C. Gasser, R.W. Ogden, A new constitutive framework for arterial wall mechanics and a comparative study of material models, J. Elast. 61 (2000) 1–48. J.D. Humphrey, Review Paper: continuum biomechanics of soft biological tissues, Proc. R. Soc. Lond. A 459 (2029) (2003) 3–46. M. Itskov, N. Askel, A class of orthotropic and transversely isotropic hyperelastic constitutive models based on a polyconvex strain energy function, Int. J. Solids Struct. 41 (2004) 3833–3848. A.N. Natali, E.L. Carniel, P. Pavan, Constitutive modelling of inelastic behavior of cortical bone, Med. Eng. Phys. 30 (2008) 905–912. A.N. Natali, C.G. Fontanella, E.L. Carniel, Constitutive formulation and analysis of heel pad tissues mechanics, Med. Eng. Phys. 32 (2010) 516–522. A.N. Natali, E.L. Carniel, A. Frigo, P.G. Pavan, S. Todros, P. Pachera, C.G. Fontanella, A. Rubini, L. Cavicchioli, Y. Avital, G.M. De Benedictis: Experimental Investigation of Urethral Tissues and Structures Biomechanics, Experimental Physiology, Volume 101, Issue 5, 1 May 2016, pp. 641–656. A.N. Natali, E.L. Carniel, P. Pavan, P. Dario, I. Izzo, Hyperelastic models for the analysis of soft tissue mechanics: definition of constitutive parameters, Biomed. Robot. Biomechatronics (2006). A.N. Natali, E.L. Carniel, H. Gregersen, Biomechanical behaviour of oesophageal tissues: material and structural configuration, experimental data and constitutive analysis, Med. Eng. Phys. 31 (2009) 1056–1062. S.A. Maas, B.J. Ellis, G.A. Ateshian, J.A. Weiss, FEBio: finite Elements for Biomechanics, J. Biomech. Eng. 134 (1) (2012) 011005. J.E. Marsden, J.R. Huges, Mathematical Foundations of Elasticity, Prentice-Hall, 1983. R.W. Ogden, Non-linear Elastic Deformations, Dover Publications, 1997. J. Schröder, P. Neff, Invariant formulation of hyperelastic transverse isotropy based on polyconvex free energy functions, Int. J. Solids Struct. 40 (2) (2003) 401–445. J. Schröder, P. Neff, Poly-, Quasi- and Rank-One Convexity in Applied Mechanics. CISM Courses and Lectures, vol. 516, SpringerWienNewYork, 2010. D.J. Steigmann, Frame-invariant polyconvex strain-energy functions for some anisotropic solids, Math. Mech. Solids 8 (2003) 497–506. C. Truesdell, W. Noll, The Non-Linear Field Theories of Mechanics, Third edition, Springer-Verlag, Berlin, 2004. J. Vincent, Structural Biomaterials, Third edition, Princeton University Press, 2012.
Contents lists available at ScienceDirect
International Journal of Non–Linear Mechanics journal homepage: www.elsevier.com/locate/nlm
The study of equivalent material parameters in a hyperelastic model c,⁎
c
L. Zanelli , A. Montanaro , E.L. Carniel a b c
a,b
, P.G. Pavan
a,b
, A.N. Natali
a,b
MARK
Department of Industrial Engineering, University of Padova, 35131 Padova, Italy Centre for Mechanics of Biological Materials, University of Padova, 35131 Padova, Italy Department of Mathematics, University of Padova, 35121 Padova, Italy
A R T I C L E I N F O
A B S T R A C T
Keywords: Hyperelasticity Constitutive parameters Polyconvexity Tissue modeling
We study a hyperelastic model of some biological soft tissues with emphasis on the problem of its matching with the material parameters acquired by experimental mechanical tests. First, we study the polyconvexity property of the hyperelastic model. Then, we explore the notion of equivalent sets of material parameters. We perform a numerical study of the regions of equivalent material parameters characterizing the curves predicted by the hyperelastic model that are close, within a prefixed tolerance, to those given by the experimental data. In the numerical study we use the quadratic variation and the Hausdorff distance. The study suggests that a qualitative knowledge of shape and volume of the regions of equivalent material parameters can provide both a criterion for the optimal match between the model with the experimental data and an indication on the reducibility of the number of parameters used in the model.
1. Introduction 1.1. An overview of soft tissue modeling Many biological soft tissues can be described by hyperelastic models, viscoelastic models, or even by more general models (see, e.g., [6,9,13,26,28,29] and their references). As a consequence, a lot of experiments can be numerically simulated and the related material parameters implemented in the model describing tissue properties can be adjusted in order to recover the experimental data. Nevertheless, there are several open issues about inverse approaches in hyperelasticity, as shown in the literature (see, e.g., [1] and references therein). We first point out that an inverse problem is usually related to a model which can be defined by different possible families of material parameters. Hence, once the experimental data and a related model have been fixed, one must adopt an optimal approach in order to localize the material parameters. The first general domain of parameters must provide physically reasonable material behavior of the model. It is well known that for hyperelastic models such a behavior is guaranteed by the polyconvexity condition of the strain energy function with respect to the deformation gradient (see, e.g., [7,20,25,2,3,26], and their references). The second more specific domain of parameters must match a known family of stress-strain curves given by the experimental data obtained for the specific tissue. With respect to this issue, various approaches can be given. Often a cost function and a suitable algorithm are employed to
⁎
evaluate and minimize the differences between experimental data and model simulations. For example, the particular simulated annealing algorithm, developed in [8], has been proved suitable for the complex behavior of the cost function. Optimization algorithms must be further enhanced because of the multimodal behavior of the cost function. Thus, coupled deterministic - stochastic algorithms have been developed (see for example [17,18] and the references therein). The approach allows both to minimize the discrepancy between experimental and model results and to explore the domain of admissible parameters. In another type of approach an easy fitting of polyconvex stored energies can be applied to soft tissues with no optimization procedure. This is shown for example in [3]. We stress that in both approaches the final target is to find at least one vector of constitutive parameters for one given tissue, thus avoiding a complete analysis of the model and its dependence from the parameters. 1.2. On the determination of material parameters In view of the different determination procedures of parameters outlined above, we stress that another related meaningful issue is that of uniqueness of the constitutive parameters for a given tissue. In fact, a given set of experiments on a soft tissue can provide only a certain number of limited information about the mechanical behavior of a material. The vector of the material parameters might be modified by increasing the set of the known experiments. This leads to the issue of
Corresponding author. E-mail address: [email protected] (L. Zanelli).
http://dx.doi.org/10.1016/j.ijnonlinmec.2016.12.014 Received 29 March 2016; Received in revised form 22 December 2016; Accepted 22 December 2016 Available online 26 December 2016 0020-7462/ © 2016 Elsevier Ltd. All rights reserved.
International Journal of Non-Linear Mechanics 89 (2017) 142–150
L. Zanelli et al.
dH (X , Y )≔max{d (X , Y ); d (Y , X )} d (X , Y )≔sup inf d (x, y ).
the non-uniqueness of the parameters inferred by the model. This is an open problem for many hyperelastic models of tissues, and it seems that there is not a general solution to it (see, e.g. the different approaches in [1,10,12,14,15,22,26,27] and the references therein). We also remark that increasing the number of experimental data could be useful to recover uniqueness of the material parameters, but this feature seems to be connected with the particular choice of the model. In some papers, the issue of non-uniqueness of the material parameters does not play any role because their goal is to represent the characteristic stress-strain curve for a given restricted set of experiments (see, e.g., Sect. 4.1 in [2]). For the hyperelastic model presented here, we consider the problems of the determination of the material parameters and also of their uniqueness (see Section 3.1). A fully incompressible hyperelastic model and an almost incompressible hyperelastic model are used to identify the constitutive parameters on different kind of tissues. This approach can be extended also to other hyperelastic models, viscoelastic models, and also to more general models.
x∈X y∈Y
It is used to compare the graphs of the curves relating the experimental data to the model. It guarantees a satisfactory analysis of the tissues. Now, let (ω) ⊂ 2 be the graph of the curve given by P11 (F, ω) for F = Fj and fixed ω, and let ⊂ 2 be the graph of the experimental curve with values yj; then the inequality ∼), ) ≤ d ( (ω ), ) d ( (ω (6) H
In what follows, we provide an overview of the main results of the paper with respect to the above discussed topics. In Section 3.1 we introduce the hyperelastic model with the strain energy function, related invariants, and constitutive parameters. In Section 3.2 we prove the property of polyconvexity of the strain energy function W = W (C, ω), where C = FT F is the right Cauchy-Green strain tensor, F is the deformation gradient and ω ∈ N is the vector whose components are the constitutive parameters (here N = 2, 4, 6 ). It turns out that W is twice continuously - differentiable with respect to C and continuous with respect to ω. It is well known that polyconvexity of the energy guarantees the physical behavior of the model. In our model such property holds when the set of constitutive parameters Ω ⊂ N has the form
2. Settings and preliminaries In this section we provide a summary of notations used in the paper together with some central definitions of hyperelasticity. We denote by Lin+ the set of all second-order tensors with positive determinant, to be identified with the family of 3×3 matrices n × n with positive determinant. The set Orth+ is the subset of Lin+ given by rigid rotation tensors R , namely such that RRT = I , where I is the unit tensor. The set Sym+ is the subset of Lin+ of symmetric tensors U , i.e. such that U = UT . Following the notations of continuum mechanics, here F ∈ Lin+ denotes the deformation gradient, C = FT F ∈ Sym+ denotes the right Cauchy-Green deformation tensor, the map F ∈ Lin+⟼W (FT F) ∈ is the strain energy function, and the first Piola-Kirchhoff stress tensor reads P = 2F∂CW . We refer the reader to the standard textbooks of continuum mechanics [13,23,24,28].
(1)
Thus within this domain we look for the regions of parameters that guarantee the match between the model (2)
P (F, ω) = 2F ∂CW (F, ω).
(first Piola-Kirchhoff tensor) and the experimental data. As it is known, in the fully incompressible case we have P = −p F−T + 2F ∂CW where p is a Lagrange multiplier. For a given family of deformation gradients {Fj : 1 ≤ j ≤ M}, let yj be the mean value of the experimental data set caused by the deformation Fj , which is taken on a given specimen. Then we look for a vector ω that minimizes the cost function
σ (ω, y )≔
1 M
M
∑
2−
j =1
yj P11 (Fj , ω)
−
P11 (Fj , ω) yj
Definition 2.1. A map F ∈ Lin+⟼W (F) ∈ is said to be convex if
2
W (τ F1 + (1 − τ ) F2) ≤ τW (F1) + (1 − τ ) W (F2) (3)
(7)
for every F1, F2 and τ with 0 ≤ τ ≤ 1. In the following, we recall some generalized convexity conditions. In particular, from [23] we recall
within a prefixed tolerance. This is done by a coupled stochastic-deterministic algorithm of the kind briefly discussed above and used for similar models in, e.g., [7,20] and in the references therein. In particular, a (stochastic) simulated annealing algorithm to find a first approximated local minimum. Then, we apply a (deterministic) simplex method to find a local minimum. After that, we need to study all the possible other material ∼ which are, in a fixed bounded region of Ω, equivalent parameters ω ∼ in such a way that to ω . We first select ω
∼, y ) ≤ σ (ω , y ). σ (ω
H
∼. As we will show in determines the second selection of the parameters ω Section 3.3, various two dimensional regions of parameters can be displayed by looking at different material parameters satisfying (4) and ∼ are minimizers of the cost (6). As explained above, the parameters ω function σ as well as the parameter ω . However, in order to determine ∼ we need first to determine at least a parameter the set of parameters ω ω and then we can write inequalities (4)–(6) and call a related routine ∼. This numerical study is that determines numerically the set of ω ® obtained through MATLAB (software house MathWorks). The related algorithm can be easily adjusted for different stress tensors and thus for different hyperelastic models. A qualitative analysis of such regions of parameters, for different experimental data, provides useful information about the mechanical properties of the tissue under study and also some indicators of the optimality of the chosen hyperelastic model. The novelty of our study is to show the possibility of a more complete analysis of the material parameters of a soft tissue (here in particular urethral tissue). This involves not only the determination of a (vector of) material parameter ω as done previously, but in addition the determination of the larger set of parameters providing the same minimization of cost function involving the model and the experiment.
1.3. Outline of the results
Ω = {(ω1, …, ωN ) = ω ∈ N ωi > 0}.
(5)
Definition 2.2. A map F ∈ Lin+⟼W (F) ∈ is said to be polyconvex if there exists a function P: 3×3 × 3×3 × → such that
W (F) = P (F, Adj[F], det[F])
(8)
∼ ∼ ∼ ∼ and (X͠ , Y , Z ) ∈ 19⟼P (X͠ , Y , Z ) ∈ is convex. The policonvexity property is usually used to select physically reasonable models (see, e.g., [26] and the references therein).
(4)
Definition 2.3. A twice differentiable function F ∈ Lin+⟼W (F) ∈ fulfills the Legendre-Hadamard condition if ∀ a, b ∈ 3, ∀ F ∈ Lin+ we have
In addition, as a second requirement we make use of the so-called Hausdorff distance dH between sets (see for example [4] and the references therein). Let us recall that the Hausdorff distance between two sets X , Y ⊂ 2 with Euclidean distance on 2 , is defined by
D2F W (F). (a ⊗ b, a ⊗ b ) ≥ 0. 143
(9)
International Journal of Non-Linear Mechanics 89 (2017) 142–150
L. Zanelli et al.
Because of the uniaxial configuration of the experimental tests, the constraints on the second and third diagonal terms P22 = P33 = 0 can be assumed. In the different settings for P as described above, these constrains provide the functional dependences λ22 = λ22 (λ11) and λ33 = λ33 (λ11). The numerical approach we used to get such two functions is obtained by the well known Newton-Raphson method (see for example [11] for the implementation by the use of MATLAB® ). With regard to the fully incompressible model, we recall that J = det(F) = 1 and
When the strain energy function W is twice continuously differentiable, the polyconvexity condition implies the weaker LegendreHadamard condition, which also guarantees a physically reasonable model (see, e.g., [2] and the references therein). 3. Main results 3.1. The hyperelastic model
P = −p F−T + 2F ∂CW
The strain energy function modeling the tissue that is considered here reads as
W (C) = Wm (C) + Wf (C, a 0 ⊗ a 0).
where p is a Lagrange multiplier determined by the boundary conditions. We stress that the aim of our paper is not to show that our energy functions in (12)–(13)–(14) are the optimal ones to modelize the urethral tissues under study with respect to other models used in the past. We recall that our primary target is to show that a complete analysis of the material parameters of a soft tissue (here in particular urethral tissue) should involve not only a single (vector of) material parameter ω but the larger set providing the same minimization of cost functions involving the model and the experiment. The analysis provided in this paper can be done with other energy functions and then one can compare the determination of the regions of material parameters with the same physical meaning.
(10)
This function W describes the behavior of a fibrous tissue with a family of fibers oriented along the local direction a 0 . (about strain energies for other kinds of tissue see, e.g., [5,7,16,20,21].) In our experiments, the tensile direction will be transversal or longitudinal. In particular, the first term in (10) reads
∼ Wm (C) = Um (J ) + W͠ m (I1)
(11)
with
Um (J ) =
Kv [(J − 1)2 + J −r + rJ − (r + 1)] 2 + r (r + 1)
C ∼ ∼ W͠ m (I1) = 1 {exp[α1 (I1 − 3)] − 1} α1
(12)
3.2. Polyconvexity
(13)
In this section, we provide the proof of the polyconvexity of the strain energy function W (C) as in (10) with respect to the deformation gradient F . For the sake of simplicity, we divide the computations into two parts: (A) Let us prove the polyconvexity of the map F⟼Wm (FT F) for Wm given by (11). As for the first term Um (J ), we recall that J = det(F). By recalling the definition of polyconvexity (see Definition 2.2) we now devote our attention to the one dimensional map J ⟼Um (J ) ∈ which a smooth function (i.e. C ∞) for J > 0 , and thus we try to prove its convexity property. The one dimensional setting ensures the convexity by the study of the second order derivative.
where J ≔(det C)1/2 = det F is the deformation Jacobian and ∼ I1 = tr(J −2/3C) is the first invariant of the isovolumetric part of C . The volumetric contribution Um of the strain energy function must be assumed to contribute in the almost incompressible case. The constitutive parameters Kv and r account for matrix compressibility, while C1 and α1 account for shear behavior. The second term in (10) reads as
Wf (C, a 0 ⊗ a 0) = Wf (I4 ) =
C4 {exp[α4 (I4 − 1)] − α4 (I4 − 1) − 1} (α4 )2 (14)
where I4 is the structural invariant defined by I4 = C: (a 0 ⊗ a 0). The parameters C4 and α4 account for the initial stiffness and increase of fibers stiffness by stretch, respectively. The decomposition (10) of the strain energy implies for the first Piola-Kirchhoff stress tensor
U ″ m (J ) =
Kv [2 + r (r + 1) J −r −2]. 2 + r (r + 1)
(21)
Thus, for r > 0 and Kv > 0 it follows
U ″m (J ) > 0,
(22)
∀ J > 0.
(15)
P = 2F ∂C W
This inequality ensures that J ⟼Um (J ) is a convex map on the domain (0, +∞). Hence, the composition of Um with J = det(F) implies that F⟼Um (det(F)) is polyconvex. ∼ As for the second term W͠ m (I1) in (11), simply notice that it is the composition (up to a constant) of the exponential (which is convex) and ∼ ∼ of the map I1 ⟼α1 (I1 − 3). More precisely,
the decomposition
P = Pm + Pf ,
(16)
where by (11)
Kv [2J (J − 1) − rJ−r + rJ ] F−T 2 + r (r + 1) 2∼ ∼ + C1 exp[α1 (I1 − 3)](2J −2/3F − I1 F−T ) 3
Pm = 2F∂CWm =
∼ I1 = tr(J −2/3F TF) = J −2/3tr(F TF) = (det(F))−2/3 ∥ F ∥2 .
C4 {exp[α4 (I4 − 1)] − 1} F (a 0 ⊗ a 0). α4
(17)
(18)
In the present paper we consider experimental data for uniaxial tensil tests along direction 1. A diagonal form of the deformation gradient can be assumed:
⎛ λ11 0 0 ⎞ ⎜ ⎟ F = ⎜ 0 λ22 0 ⎟ . ⎝ 0 0 λ33 ⎠
(23)
∼ Now we recall Lemma C.1 in [25], so that F⟼I1 (F) is polyconvex. By the sum of Um and W͠ m we can therefore conclude that F⟼Wm (FT F) is polyconvex. (B) Let us prove the convexity of the map F⟼Wf (FT F, a 0 ⊗ a 0) for Wf given by (14). This is a function of the structural invariant I4 = FT F: (a 0 ⊗ a 0) and in more details
and by (14)
Pf = 2F∂CWf = 2
(20)
I4≔tr((FT F)T (a 0 ⊗ a 0)) = tr((FT F)(a 0 ⊗ a 0))
FT F .
(24)
By recalling Lemma C.2 thanks to the symmetry of the matrix (point 1) in [25] we have that F⟼I4 (F) is convex. Now, it can be easily checked that I4 ⟼Wf (I4 ) ∈ is a smooth function and that the second order derivative reads
(19) 144
International Journal of Non-Linear Mechanics 89 (2017) 142–150
L. Zanelli et al.
Wf (I4 )″ =
C4 {[α4]2 exp[α4 (I4 − 1)]}. (α4 )2
(25)
For C4 > 0 , α4 > 0 we recover
Wf (I4 )″ = C4 exp[α4 (I4 − 1)] > 0,
∀ I4 ∈ ,
(26)
and hence also the convexity condition. Recalling Lemma B.9 in [25], we are now in the position to make the composition of the monotone increasing (when I4 ≥ 1) and convex function Wf with I4 as a function of F , and to get the convexity of F⟼Wf (FT F, a 0 ⊗ a 0). In the case 0 < I4 < 1, we first notice that if I4 is sufficiently closed to 1, then D2F Wf ≃ C4 D2F I4 which is semipositive definite, and whence Wf is convex. For the arbitrary case 0 < I4 < 1 we can always make a rescaling δ (I4 − 1) by a small δ > 0 so that D2F Wf becomes semipositive definite. Notice that such a transformation can be viewed as a rescaling of the parameter α4 > 0 , which remains positive; whence we can say that the hyperelastic model is not really changed and the match between the model and the experimental data involves such a rescaling of the parameter α4.
Fig. 1. Tensile test for distal urethral specimen comparison of experimental data and fully incompressible model results.
Remark 3.1. In view of the above computations, the constitutive parameters satisfying Kv > 0 , r > 0 , C1 > 0 , α1 > 0 , C4 > 0 , α4 > 0 , guarantee physical type responses for our hyperelastic model. As we will see in the next section, for different stress-strain curves, the involved region of equivalent material parameters will be bounded, namely there is not degeneracy of unbounded parameters.
ω = (C1, α1) = (1.09·10−3, 1.26). The two stress-strain curves coming from the experimental data (the mean values) and from the values of Pm,11 evaluated at ω are given in Fig. 1 Related to this case, Fig. 1 shows the two-dimensional region of equivalent material parameters with respect to (C1, α1): The above figure is determined by the check of the constraints (4)– (6); namely we look if they are fulfilled (or not) on the square net of points centered on ω = (C1, α1) given by ± 10 points with respect to the first and second coordinate, and length 0.5 between them. The same procedure is done for the other pairs of parameters and figures in the next part of the paper. We first notice that the region of parameters shown in Fig. 2 is connected, and no other disjoint regions appear. Physically, this means we have the uniqueness of parameters and that any rescaling of (C1, α1) changes the volume of the region. Moreover, we notice that the geometry of the region has a symmetry around a line. This fact suggests the introduction of a transformation given by anticlockwise rotation with center at (1.09·10−3, 1.26) and angle ϕ = 46°, whence from (C1, α1) 1, α1). We can say that α1 is varying slowly with respect to a new pair (C 1, since we can also introduce a further transformation given by an to C (arbitrary) rescaling of α1. Thus, it can be recognized that the model is more sensible with respect to one parameter with respect to another one. In this two-parameter model, such a behavior is expected and other arguments can be used to recover the same observation. We stress that when the parameters space has a dimension greater than two, then such an analysis can be more difficult to recover by a direct study of the
3.3. Domains of material parameters Our paper involves the analysis of experimental data from mechanical tests developed on urethral specimens (see [19]). Urethral tissues show different conformation depending on the specific location. Proximal urethra tissues are composed of an almost isotropic connective matrix reinforced by muscular fibers aligned along longitudinal direction. On the other hand, distal urethral tissues are composed by almost isotropic connective matrix. The mechanical characterization of urethral tissues took into account tensile tests performed on proximal and distal samples, according to both longitudinal and transversal directions. With specific regard to distal samples, results from tests along each different directions did not show statistical differences, enforcing the assumption of isotropic behavior. It follows that we provide six data sets given by experimental stress-strain curves. As we will describe in this section, for any of these experimental stress-strain curves, we provide the interpolation given by the hyperelastic model described in the previous section and a vector of material parameters. Later, we show that for any of these fixed vectors of material parameter we are able to localize all the equivalent material parameters within a (large) bounded region. For any of these figures representing regions of parameters we provide a qualitative analysis which leads to the possible reduction in the number of parameters, the determination of leading parameters, the study of the behavior of equivalent parameters with respect to data obtained by transversal or longitudinal tensils on specimens. Morever, the study of uniqueness or non-uniqueness of material parameters is given by the determination of non-connected regions of equivalent parameters; this approach is not affected by the rescaling of the parameters and whence is well posed, and allows to adjust for any type of material parameter the correct scale to represent its value. Let us begin by looking at the case of a distal urethra specimen by the fully incompressible model. The first Piola-Kirchhoff stress tensor we use here is given by (20) for the energy ∼ W = W͠ m (I1) (27) as given by (13) and with the material parameters (C1, α1). In this case, the minimization of the cost function (3) provides the parameters
Fig. 2. Graph of equivalent (C1, α1) given by computations over a net of 10×10 points with length 0.5·10−3 and 0.5 centered on (C1, α1) = (1.09·10−3, 1.2) .
145
International Journal of Non-Linear Mechanics 89 (2017) 142–150
L. Zanelli et al.
Fig. 3. Family of level sets of the parameters (C1, α1) for distal urethral specimen.
stress tensor and maybe also unexpected. The problem of the reduction of the number of parameters can thus be faced by using the above qualitative analysis of regions of equivalent parameters. Next, we provide a diagram involving a family of level sets (and neighborhoods of them) for the cost function dH in (6) involving (C1, α1) which are the parameters shown in Fig. 2. In this case, the other cost function σ in (4) is constant on the level sets of dH , since here the constraint dH is stronger than σ. In Fig. 3 we see the boundary of region shaped in the Fig. 2 and in addition the higher level sets of the cost function. This shows a unique local (and hence global) minimum of the cost function and implies that the region of parameters around this minimum is connected, namely no other disjoint regions of parameters appears. Similar work can be done for the other cases discussed in the paper. We now look at the proximal urethra specimen by the fully incompressible model, with energy function ∼ W = W͠ m (I1) + Wf (I4 ) (28)
Fig. 4. (a) - (b) Tensile test along transversal and longitudinal direction for proximal urethra specimen by the fully incompressible model.
see (13) and (14). The stress tensor of the fully incompressible model we use here is (20), and the material parameters are (C1, α1, C4, α4 ). The minimization of the cost function in (3) provides the parameters (C1, α1) = (1.12·10−3, 2.24) and (C4, α4 ) = (1.49·10−3, 5.16·10−1). The stress-strain curves are shown in Fig. 4 With respect to Fig. 4, we underline that the parameters space has dimension 4, whence the number of all the possible pairs of parameters is 6. In Fig. 5 we choose to display only the pairs (C1, α1) and (C4, α4 ). The two regions of equivalent parameters linked to Fig. 4 - (a) are given by Figs. 5–6. We first remark that in this case also (C4, α4 ) are varying in a square net of points centered on (1.49·10−3, 5.16·10−1) with length 0.5. Thus, Fig. 5 is the projection of the 4-dimensional set of all the equivalent parameters (C1, α1, C4, α4 ) into the subspace of the parameters (C1, α1). Notice that the region of points in the Fig. 5 is connected. Whence, also in this case we recover uniqueness of parameters up to a rescaling of the volume of this domain. We stress also that the region is strictly included with respect to the net of points where we made the computations. In fact, a possible nonconnected region which is far away can be localized by the same algorithm. We also observe that the above region can be qualitatively described by a sublevel set given by an inequality of type C1 ≤ f (α1) where f : + → + is some decreasing function, which can be easily interpolated thanks to the boundary of the domain in Fig. 5. In view of the above observations, we think that a possible general characterization of such sublevels in connection to some important features of the strain energy function is an interesting future target. These could be for example the choice of the invariants and the asymptotic behavior of the energy with respect to the parameters. We also stress that the determination of the above inequality from the analytical viewpoint without passing through the numerical study can
Fig. 5. Graph of equivalent (C1, α1) given by computations over a net of 10×10 points centered on (C1, α1) = (1.12·10−3, 2.24) with length 0.5·10−3 and 0,5.
be a rather non-trivial task. The region involving the other pair of parameters (C4, α4 ) is shown in Fig. 6. We remark that here we fix the parameters (C1, α1) and thus it is a different projection with respect to the previous one we used into the Fig. 5. Indeed, this is done because here we are looking at the projection in the space of (C4, α4 ) of a subset (i.e. with fixed (C1, α1)) of the 4dimensional domain of all the equivalent parameters which we are taken into account to recover Fig. 5. We see that the above domain in Fig. 6 is contained in the net of 10×10 square region of points, and this suggest that there not exist others domains to detect, whence the 146
International Journal of Non-Linear Mechanics 89 (2017) 142–150
L. Zanelli et al.
Fig. 8. Graph of equivalent (C4, α4 ) for fixed (C1, α1) given by computations over a net of 10×10 points centered on (C4, α4 ) = (1.49·10−3, 5.16·10−1) with length 0.5·10−3 and 0.5.
Fig. 6. Graph of equivalent (C4, α4 ) for fixed (C1, α1) given by computations over a net of 10×10 points centered on (C4, α4 ) = (1.49·10−3, 5.16·10−1) with length 0.5 × 10−3 and 0.5.
uniqueness of parameters up to rescalings of a single connected region. Obviously, different projections of the type in Fig. 6 can be given, by fixing any other equivalent pair of parameters (C1, α1) showed in Fig. 5, but analogous information are recovered. In Fig. 7 we provide the two regions of parameters linked to Fig. 4 (b). If we compare Figs. 7 and 5 we notice a very similar shape, but with a slight increase of the volume of the parameters. This fact suggests to consider the smallest region of parameters in Fig. 5 as the characterizing material parameters. By looking at Figs. 6 and 8 we notice a relevant change of the domain, which is in particular a subset of the previous one. Such a fact can be considered as an optimal feature of the model, namely that different experimental data on the same tissue allows to localize a subset of material parameters. This behavior allows to consider the smallest one as the characterizing material parameters. We underline the usefulness of this powerful approach which make use of transversal and longitudinal data to get characterizing material parameters. Let us now discuss the results on the distal urethral specimen for tensil in the longitudinal direction with the almost incompressible model. The energy function is
∼ W = Um (J ) + W͠ m (I1)
Fig. 9. Tensile test for distal urethral specimen, comparison of experimental data along longitudinal direction and almost incompressible model results.
parameters obtained by the minimization of the cost function are (K , r ) = (3.75, 2.77·10) and (C1, α1) = (1.09·10−3, 1.26). The stress-strain curves are given by Fig. 9. In what follows we consider all the 4-dimensional set of equivalent parameters (K , r , C1, α1) and we first display the projection into the 2dimensional space of the pair (K,r) Fig. 10. We stress that the region of the parameters is not completely
(29)
where the two terms on the righthand side are given by expressions (12) and (13). The first Piola-Kirchhoff stress tensor we use is Pm as defined in (17) where the involved parameters are (K,r) and (C1, α1). The
Fig. 7. Graph of equivalent (C1, α1) given by computations over a net of 10×10 points centered on (C1, α1) = (1.12·10−3, 2.24) with length 0.5·10−3 and 0.5.
Fig. 10. Graph of equivalent (K, r) given by computations over a net of 5×5 points centered on (K , r ) = (3.75, 2.77·10) with length 1.
147
International Journal of Non-Linear Mechanics 89 (2017) 142–150
L. Zanelli et al.
Fig. 11. Graph of equivalent (C1, α1) for fixed (K, r) given by computations over a net of 5×5 points centered on (C1, α1) = (1.09·10−3, 1.26) with length 10−3 and 1.
contained in the net of points we used for the comparison of the parameters. However, some simple arguments (which involve the asymptotics of Pm with respect to the parameters) show that the parameters must be bounded for a fixed stress-strain curve. About the other related couple of material parameters, we fix (K,r) and we provide the region of equivalent (C1, α1) in Fig. 11. In this case we have two non-connected components. However, a more complete description involves the projected parameters (C1, α1) for all the equivalent (K,r) given by the previous figure. This analysis provides, also in this case, a connected region of parameters. This local description is meaningful since it shows that not connected components can appear. The last case we discuss is about proximal urethral specimens in the almost incompressible setting for transversal and longitudinal tensils. ∼ W = Um (J ) + W͠ m (I1) + Wf (I4 ) (30)
Fig. 12. (a) - (b) Tensile test along transversal and longitudinal direction for proximal urethra specimen by the almost incompressible model.
where the three terms on the righthand side are given by (12), (13) and (14). The first Piola-Kirchhoff stress tensor here used is Pm + Pf as defined in (17)–(18) for parameters (K,r), (C1, α1) and (C4, α4 ). In particular, the parameters ω minimizing the cost function are (K , r ) = (6.45, 2.60) , (C1, α1) = (1.12·10−3, 2.24) and (C4, α4 ) = (1.49·10−3, 5.16·10−1). The stress-strain curves are given in Fig. 12. Since in this case the number of the parameter is 6, and whence larger than in the previous cases, we need to take into account the computational cost of the algorithm. Indeed, the comparison of 106 vector of parameters (instead of 102 or 104) implies a meaningful increase of the time for the computation. For this reason, here we are going to show only three pair of parameters when we have fixed the other four according to the vector of parameters ω . In Fig. 13 we provide the equivalent parameters of the transversal data sets shown into the above Fig. 12 - (a). In Fig. 13 we notice a large set of parameters which can be qualitatively described by inequalities of type K < g (r ) for an increasing function g: + → +. Also in this case, we can provide an interpolation of the boundary of this region in order to get an analytical description of such a function. We also underline that the above graph shows a slow dependence from r, and as a consequence this fact suggest to introduce a large rescaling of the energy W with respect to r. The second and the third pairs of equivalent parameters exhibit similar behavior when compared each other. The second pairs is given in Fig. 14, which shows that The behavior of α1 is rapidly decreasing when C1 is increasing. Thus, we can say that this hyperelastic model has a sensibility with respect to small variations of C1. Such a dependence with respect to its material parameter can be considered optimal, since it allows to directly localize small regions of parameters. In Fig. 14 we notice also a region of type
Fig. 13. Graph of equivalent (K, r) for fixed (C1, α1) and (C4, α4 ) given by computations over a net of 10×10 points centered on (K , r ) = (6.45, 2.60) with length 1.
C1 ≤ h1 (α1) for some decreasing function h1: + → + . The possible diverging behavior of α1 for very small values of C1 suggests to consider the region of parameters 0 < C1, min ≤ C1 for some small fixed material parameter C1, min . The same arguments used above here for Fig. 15 suggests the consideration of a lower bound of type 0 < C4, min ≤ C4 and inequality C4 ≤ h4 (α4 ) for some decreasing function h4 : + → + . With respect to the above experimental data in Fig. 12 - (b), the fist pair of material parameters are shown in Fig. 16. Notice that Fig. 16 is the same as Fig. 13 and hence the longitudinal (or transversal) tensile on this tissue does not change the determination of the material parameters (K,r). The second pair of parameters linked 148
International Journal of Non-Linear Mechanics 89 (2017) 142–150
L. Zanelli et al.
Fig. 14. Graph of equivalent (C1, α1) for fixed (K, r) and (C4, α4 ) given by computations over a net of 10×10 points centered on (C1, α1) = (1.12·10−3, 2.2) with length 10−3, 1.
Fig. 17. Graph of equivalent (C1, α1) for fixed (K, r) and (C4, α4 ) , given by computations over a net of 10×10 points centered on (C1, α1) = (1.12·10−3, 2.24) with length 10−3, 1.
Fig. 15. Graph of equivalent (C4, α4 ) for fixed (K, r) and (C1, α1) , given by computations over a net of 10×10 points centered on (C4, α4 ) = (1.49·10−3, 5.16·10−1) with length 10−3 and 1.
Fig. 18. Graph of equivalent (C4, α4 ) for fixed (K, r) and (C1, α1) , given by computations over a net of 10×10 points centered on (C4, α4 ) = (1.49·10−3, 5.16·10−1) with length 1.
In contrast to the previous two figures, here we stress that the Fig. 18 is different from Fig. 15, and in particular it is a small subset. This suggest to consider the points in Fig. 18 (and small domains around them which can be seen at a different scale) as the third couple of material parameters characterizing the tissue. 4. Conclusions In this paper we have provided a numerical study of the regions of equivalent material parameters for urethra type specimens. For any tissue of such kind, we displayed the curves coming from the hyperelastic model and we have shown the correspondence with the curves coming from experimental data up to some small values for the quadratic variation and the Hausdorff distance function. We stress that we used the Hausdorff distance function in order to recover a better match between model and experiments with respect to the one obtained only by the quadratic variation. We have shown that a qualitative study of the shape and the volume of such regions of parameters provides a criterion for the optimal match between the model and the experimental data, together with a check on the possible reduction of the number of parameters involved in the model. Indeed, for any of the figures involving the regions of parameters we provide a qualitative analysis which leads to the possible reduction of the number of parameters related to the determination of leading parameters, the study of the behavior of equivalent parameters with respect to data obtained by transversal or longitudinal tensils on specimens. Moreover, the study of uniqueness or non-uniqueness of material
Fig. 16. Graph of equivalent (K, r) for fixed (C1, α1) and (C4, α4 ) , given by computations over a net of 10×10 points centered on (K , r ) = (6.45, 2.60) with length 1.
to Fig. 12 - (b) is displayed in Fig. 17. As in the previous case, we stress that Figs. 14 and 17 are coinciding. We interpret this fact as a sort of stability result also for (C1, α1) which ensures an intrinsic character for this region of material parameters. We recall that a rescaling of this region is possible, whence a determination of a small volume just related to a different scale. Fig. 18 describes the third pair of material parameters. 149
International Journal of Non-Linear Mechanics 89 (2017) 142–150
L. Zanelli et al.
parameters is given by the determination of non-connected regions of equivalent parameters. Data from tensile tests only have been analyzed to identify the parameters of both fully incompressible and almost incompressible formulations. Physically, tensile tests do not provide information about the compressibility of the material and the proposed algorithm returns the known uniqueness of the parameters. Furthermore, a notion of stability of the hyperelastic model is given through the analysis of the behavior of the regions of equivalent parameters when the transversal and longitudinal tensils are considered on the same tissue. In order to complete the remarks on the paper, we underline that the approach here described can be used for any hyperelastic model, and that the features here we discussed can be used to select an optimal model for any given experimental data set. The proposed approach allows to evaluate relationships between model configuration and independent sets of experimental in order to obtain uniqueness of the parameters. Further studies can investigate the influence of more complex experimental situations, as biaxial or shear loading.
[9] [10]
[11]
[12]
[13] [14]
[15] [16]
[17] [18] [19]
References [1] S. Avril, S. Evans, K. Miller, Inverse problems and material identification in tissue biomechanics, Spec. Issue J. Mech. Behav. Biomed. Mater. 27 (November) (2013). [2] D. Balzani, J. Schröder, P. Neff, G.A. Holzapfel, A variational approach for materially stable anisotropic hyperelasticity, Int. J. Solids Struct. 42 (2005) 4352–4371. [3] D. Balzani, J. Schröder, P. Neff, G.A. Holzapfel, A Polyconvex framework for soft biological tissues. Adjustment to experimental data, Int. J. Solids Struct. 43 (2006) 6052–6070. [4] E. Belogay, C. Cabrelli, U. Molter, R. Shonkwiler, Calculating the hausdorff distance between curves, Inf. Process. Lett. 64 (1) (1997) 17–22. [5] J. Bonet, A.J. Burton, A simple orthotropic, transversely isotropic hyperelastic constitutive equation for large strain computation, Comput. Methods Appl. Mech. Eng. 162 (1–4) (1998) 161–164. [6] R. Cai, F. Holweck, Z. Feng, F. Peyraut, A new hyperelastic model for anisotropic hyperelastic materials with one fiber family, Int. J. Solids Struct. 84 (1) (2016) 1–16. [7] E.L. Carniel, V. Gramigna, C.G. Fontanella, C. Stefanini, A.N. Natali, Constitutive formulations for the mechanical investigation of colonic tissues, J. Biomed. Mater. Res A. 102 (5) (2014) 1243–1254. http://dx.doi.org/10.1002/jbm.a.34787 (Epub 2013 May 24). [8] A. Corana, N. Marchesi, C. Martini, S. Ridella, Minimizing multimodal functions of
[20]
[21]
[22] [23] [24] [25] [26] [27] [28] [29]
150
continuous variables with the Simulated Annealing algorithm, ACM Trans. Math. Softw. 13 (3) (1987) 262–280. S.C. Cowin, S.B. Doty, Tissue Mechanics, Springer Science and Business Media, 2007. L. Dubuis, S. Avril, J. Debayle, P. Badel, Identification of the material parameters of soft tissues in the compressed leg, Comput. Methods Biomech. Biomed. Eng. 15 (1) (2012) 3–11. W. Gander, M.J. Gander, F. Kwok, Scientific Computing - An Introduction using Maple and MATLAB, Texts in Computational Science and Engineering, 2014th Edition, Springer. M. Geerligs, L.C.A. van Breemen, G.W.M. Peters, P.A.J. Ackermans, F.P.T. Baaijens, C.W.J. Oomens, In vitro indentation to determine the mechanical properties of epidermis, J. Biomech. 44 (2011) 1176–1181. G.A. Holzapfel, Nonlinear Solid Mechanics, Wiley, New York, 2000. G.A. Holzapfel, T.C. Gasser, R.W. Ogden, A new constitutive framework for arterial wall mechanics and a comparative study of material models, J. Elast. 61 (2000) 1–48. J.D. Humphrey, Review Paper: continuum biomechanics of soft biological tissues, Proc. R. Soc. Lond. A 459 (2029) (2003) 3–46. M. Itskov, N. Askel, A class of orthotropic and transversely isotropic hyperelastic constitutive models based on a polyconvex strain energy function, Int. J. Solids Struct. 41 (2004) 3833–3848. A.N. Natali, E.L. Carniel, P. Pavan, Constitutive modelling of inelastic behavior of cortical bone, Med. Eng. Phys. 30 (2008) 905–912. A.N. Natali, C.G. Fontanella, E.L. Carniel, Constitutive formulation and analysis of heel pad tissues mechanics, Med. Eng. Phys. 32 (2010) 516–522. A.N. Natali, E.L. Carniel, A. Frigo, P.G. Pavan, S. Todros, P. Pachera, C.G. Fontanella, A. Rubini, L. Cavicchioli, Y. Avital, G.M. De Benedictis: Experimental Investigation of Urethral Tissues and Structures Biomechanics, Experimental Physiology, Volume 101, Issue 5, 1 May 2016, pp. 641–656. A.N. Natali, E.L. Carniel, P. Pavan, P. Dario, I. Izzo, Hyperelastic models for the analysis of soft tissue mechanics: definition of constitutive parameters, Biomed. Robot. Biomechatronics (2006). A.N. Natali, E.L. Carniel, H. Gregersen, Biomechanical behaviour of oesophageal tissues: material and structural configuration, experimental data and constitutive analysis, Med. Eng. Phys. 31 (2009) 1056–1062. S.A. Maas, B.J. Ellis, G.A. Ateshian, J.A. Weiss, FEBio: finite Elements for Biomechanics, J. Biomech. Eng. 134 (1) (2012) 011005. J.E. Marsden, J.R. Huges, Mathematical Foundations of Elasticity, Prentice-Hall, 1983. R.W. Ogden, Non-linear Elastic Deformations, Dover Publications, 1997. J. Schröder, P. Neff, Invariant formulation of hyperelastic transverse isotropy based on polyconvex free energy functions, Int. J. Solids Struct. 40 (2) (2003) 401–445. J. Schröder, P. Neff, Poly-, Quasi- and Rank-One Convexity in Applied Mechanics. CISM Courses and Lectures, vol. 516, SpringerWienNewYork, 2010. D.J. Steigmann, Frame-invariant polyconvex strain-energy functions for some anisotropic solids, Math. Mech. Solids 8 (2003) 497–506. C. Truesdell, W. Noll, The Non-Linear Field Theories of Mechanics, Third edition, Springer-Verlag, Berlin, 2004. J. Vincent, Structural Biomaterials, Third edition, Princeton University Press, 2012.