Next Article in Journal
Dynamics of Classical Solutions of a Two-Stage Structured Population Model with Nonlocal Dispersal
Previous Article in Journal
Improved Properties of Positive Solutions of Higher Order Differential Equations and Their Applications in Oscillation Theory
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Anisotropic Hyperelastic Material Characterization: Stability Criterion and Inverse Calibration with Evolutionary Strategies

by
Claudio Canales
1,
Claudio García-Herrera
1,*,
Eugenio Rivera
1,
Demetrio Macías
2 and
Diego Celentano
3
1
Departamento de Ingeniería Mecánica, Universidad de Santiago de Chile, Santiago 9170020, Chile
2
Laboratory Light, Nanomaterials & Nanotechnologies—L2n, University of Technology of Troyes & CNRS EMR 7004, 12 Rue Marie Curie, CS 42060, 10004 Troyes, France
3
Departamento de Ingeniería Mecánica y Metalúrgica, Pontificia Universidad Católica de Chile, Av. Vicuña Mackenna 4860, Santiago 7820436, Chile
*
Author to whom correspondence should be addressed.
Mathematics 2023, 11(4), 922; https://doi.org/10.3390/math11040922
Submission received: 27 December 2022 / Revised: 1 February 2023 / Accepted: 7 February 2023 / Published: 11 February 2023
(This article belongs to the Section Computational and Applied Mathematics)

Abstract

:
In this work, we propose a reliable and stable procedure to characterize anisotropic hyperelastic materials. For this purpose, a metaheuristic optimization method known as evolutionary strategies is used. The advantage of this technique with respect to traditional methods used for non-linear optimization, such as the Levenberg–Marquardt Method, is that this metaheuristic algorithm is oriented to the global optimization of a problem, is independent of gradients and allows to solve problems with constraints. These features are essential when characterizing hyperelastic materials that have non-linearities and are conditioned to regions of stability. To characterize the mechanical behavior of the arteries analyzed in this work, the anisotropic hyperelastic models of Holzapfel–Gasser–Ogden and Gasser–Holzapfel–Ogden are used. An important point of the analysis is that these models may present a non-physical behavior: this drawback is overcome by defining a new criterion of stabilization in conjunction with the evolutionary strategies. Finally, the finite element simulations are used in conjunction with the evolutionary strategies to characterize experimental data of the artery pressurization test, ensuring that the parameters obtained are stable and representative of the material response.

1. Introduction

Hyperelastic materials are used in multiple areas of interest, such as industrial applications [1,2,3], soft robotics [4,5], or soft tissue modeling [6]. They describe the behavior of materials that exhibit large deformations and usually a non-linear stress–strain relationship. It is possible to classify these materials into two main groups: isotropic and anisotropic. The latter is generally composed of an isotropic matrix reinforced with fibers that induce anisotropy [7]. Soft tissues are an excellent example of this type of material, such as arteries [8,9], veins [10], corneas [11], vaginal tissue [12], muscles [13] or skin among others [14]. Among soft tissues, arteries are one of the best-documented materials and also present multiple constitutive models such as Holzapfel–Gasser–Ogden (HGO) [15] and Gasser–Holzapfel–Ogden Model (GHO) [16].
There is significant research interest in adequately characterizing these materials to increase the accuracy of computational models. The characterization problem can be seen as an optimization one, which determines the constitutive model’s parameters that minimize the discrepancy concerning the experimental data. A wide range of experimental techniques are available to characterize hyperelastic materials, such as uniaxial tests in different directions [17], biaxial [18], bulge [19], pressurization [20,21], and AFM nanoindentation [22]. In each of these cases, the stress is modeled as a function of the material parameters. There are different alternatives for this, such as describing the deformation gradient as homogeneous [23], identifying the displacement field by optical means [24], or numerical simulations of the mechanical tests [25]. Describing the deformation gradient as homogeneous allows for analytical equations, but it is insufficient to describe the tests adequately in many cases. With optical methods, the strain gradient can be obtained with experimental measurements, but the experimental setup becomes complex and requires preprocessing and postprocessing the obtained displacement fields. The numerical simulations describe the mathematical model in detail but require extensive computational resources and therefore time-consuming. Regardless of the method used to model the mechanical tests, it is necessary to define an objective function to quantify the discrepancy between the mathematical model and the experimental values. The material parameters that minimize the objective function will be those that characterize the mechanical behavior of the material. In addition, it is necessary that these parameters correctly model different modes of deformation (e.g., uniaxial, equiaxial, pressurization), capture anisotropy if it exists, and meet certain (specific) stability criteria (e.g., Drucker [26] in isotropic materials).
Several authors have defined the objective function as the sum of the quadratic errors of the discrepancy between the mathematical model and the experimental data [27]. In order to minimize this function, gradient-dependent algorithms such as non-linear least-squares [28], Levenberg Marquardt [29] have been used to minimize the function. However, these algorithms only guarantee finding local optima, and the function must be smooth. Furthermore, minimizing the objective function suffers from several difficulties, such as the impossibility of obtaining an analytical derivative (e.g., numerical simulations), the presence of valley regions (derivatives are almost zero), stability constraints, and the existence of multiple local optima. To deal with these problems, some authors have used metaheuristic optimization algorithms to characterize hyperelastic materials to deal with these problems. Fernández et al. [30] have used genetic algorithms to characterize the three-parameter Mooney Rivlin model considering only a uniaxial test and assuming that the strain gradient is homogeneous. Another recent publication of López-Campos [31], uses a genetic algorithm to characterize tendons as a hyperelastic material with damage. Ramzanpour et al. [32] have proposed a constrained particle swarm optimization algorithm (C-PSO) adapted to consider the elastic compatibility equations and the Drucker stability criterion. This method characterizes isotropic hyperelastic materials with a viscoelastic model.
In the case of anisotropic materials, it is necessary to consider multiple mechanical tests to reveal the mechanical behavior. Due to the simplicity of uniaxial testing, it is possible to model the deformation gradient as homogeneous [33]. In this model type, one generally has a highly non-linear objective function. However, in the case of other tests, analytical models are restricted to ideal conditions assumptions such as the thin wall condition [20]. In addition to the difficulties of mathematical modeling, a study revealed a nonphysical response when the energy density function is defined as a sum of a deviatoric and a volumetric component [34]. Furthermore, the relevance of transverse deformations has been studied, where under a uniaxial tensile condition, the HGO and GHO models can present work increasing or decreasing transverse elongations. Some hyperelastic models allow prescribing the behavior of transverse elongations as What You Prescribe Is What You Get (WYPIWYG) [35]. However, it is possible to prescribe the behavior of transverse stretches in formulating the objective function.
The novel aspects of this work focus on exploring the use of evolution strategies to characterize anisotropic hyperelastic materials. In addition, a stabilization criterion is established to deal with unexpected physical responses, such as transverse elongations growing under a uniaxial tensile condition. This instability is studied analytically, and it is shown that it originates when the anisotropic stiffness grows disproportionately with respect to the isotropic stiffness. The instability equations show that the HGO model will always be unstable for large deformations and the GHO model will present stable parameters belonging to a domain defined by the stability constraint. Finally, it is presented how arteries can be characterized using multiple mechanical tests to capture anisotropy, numerical simulations to describe inhomogeneous deformation gradients, evolutionary strategies as a global optimization algorithm, and the stability criterion proposed in this work. This procedure can be extended to any anisotropic material. This study formulates the stability criterion and studies in depth the use of evolutionary strategies in the characterization of hyperelastic materials, previously used in an experimental study by Rivera et al. [21].

2. Materials and Methods

2.1. Anisotropic Hyperelastic Modeling

Multiple types of materials such as elastomers or soft tissue can be modeled as hyperelastic because they exhibit elastic behavior and large deformations [36]. The deformation of these materials can be represented locally by the deformation gradient, denoted by F , and the energy density function can describe the elastic response of the material W, which generally, depends on the deformation gradient F . For W to comply with the principle of objectivity, the dependence of F is through the right Cauchy–Green tensor C = F T F . Therefore, W = W ( C ) and Cauchy’s stress tensor σ for an incompressible material is given by [37]:
σ = p I + F W C F T
where p is a Lagrange multiplier associated with the constraint of incompressibility. In order to describe the mechanical behavior of fiber-reinforced materials, it is necessary to define the orientation of the fibers, and generally, one or two families of fibers are used. The unit vectors M and M are used to describe the orientation of the fiber families in the initial configuration. The vectors of the fibers in the deformed configuration are represented by the vectors m and m , where m = FM and m = FM which are not a unitary vector. According to the work of Spencer [38], a strain energy function of a compressible elastic material with two preferred unit directions M and M can be expressed as W ( C , M M , M M ) . As stated by Shariff [39], this strain energy function W can be expressed in terms of a set of invariants I i i = 1 , . . . , 10 . However, multiple authors of hyperelastic models such as Lin and Yin [40] or Gasser et al. [16] neglect most of the invariants to reduce the dependency of W and to facilitate the correlation of the experimental data [41].
In this work, the GHO model is used to characterize anisotropic hyperelastic material (arteries), and it is defined by:
W ( I 1 , I 4 , I 6 ) = μ 2 ( I 1 3 ) + k 1 2 k 2 i = 4 , 6 ( exp ( k 2 ( I 1 κ + ( 1 3 κ ) I i 1 ) 2 ) 1 ]
This model has three material parameters ( μ , k 1 , k 2 ) and a parameter that defines how dispersed the fibers are in the space, which is κ and belongs between [ 0 , 1 / 3 ] . It is possible to observe that only the invariants I 1 = t r ( C ) , I 4 = M · CM and I 6 = M · CM are considered. It is important to emphasize that the material will be considered as incompressible I 3 = d e t ( C ) = 1 . If the value of κ = 0 then the GHO model is equivalent to the well-known Holzapfel—Gasser–Ogden model (HGO) [42].

2.2. Material Characterization via the Tensile Test

This section presents the analytical equations of the uniaxial tensile test for a density energy function W that depends on I 1 , I 4 , I 6 [43]. Furthermore, an objective function is provided that minimizes the differences between the analytical model and the experimental data.
In order to model the uniaxial tensile test, it is assumed that the anisotropic material is incompressible and has a homogeneous deformation gradient:
F = λ 1 0 0 0 λ 2 0 0 0 λ 3
where λ i ( i = 1 , 2 , 3 ) are the principal stretches of the test sample. Figure 1 schematizes the uniaxial test. It is important to remember that as the material is incompressible and λ 3 = ( λ 1 λ 2 ) 1 . For this kind of constitutive models, it is not possible to obtain an explicit expression of σ 11 (stress associated with the main direction of a tensile test). According to the procedure described by [20], the equations that model the uniaxial tensile test analytically are:
σ 11 = 2 ( λ 1 2 λ 1 2 λ 2 2 ) ( W 1 ) + 2 m 1 2 W 4 + 2 m 1 2 W 6
σ 22 = 2 ( λ 2 2 λ 1 2 λ 2 2 ) ( W 1 ) + 2 m 2 2 W 4 + 2 m 2 2 W 6
where W i = W I i ( i = 1 , 4 , 6 ) and the values of m i and m i ( i = 1 , 2 , 3 ) are the components of the two fiber vectors in the deformed configuration. It is possible to observe that σ 11 depends on λ 1 and λ 2 so it is necessary to use Equation (5) to find the value of λ 2 , which is done by solving this nonlinear equation through the Newton–Raphson method.
Arteries are modeled recurrently as an anisotropic hyperelastic material, and they are used in this work to exemplify the characterization methodology. However, all the described procedures can be extended to characterize any anisotropic hyperelastic material with one or two preferential directions. For example, arteries generally possess two families of collagen fibers that extend in a helical shape, so locally the arterial wall has a plane of anisotropy normal to the radial direction, which is conformed by the fibers. Families of collagen fibers share mechanical properties, and they are symmetrical with the longitudinal axis, which forms an angle Θ , in the material configuration. The symmetry is shown in Figure 2.
Since the families of fibers are symmetrical, it is possible to deduce for the uniaxial traction case that the invariant I 4 = I 6 . Therefore, the energy density function can be written as W ( I 1 , I 4 ) . The Cauchy stress Equations (4) and (5) are expressed as:
σ 11 = 2 ( λ 1 2 λ 1 2 λ 2 2 ) ( W 1 ) + 4 m 1 2 W 4
σ 22 ( λ 1 , λ 2 , x ) = 0 = 2 ( λ 2 2 λ 1 2 λ 2 2 ) ( W 1 ) + 4 m 2 2 W 4
With these equations, it is possible to characterize the material employing an implicit model of the uniaxial test. Thus, the tensile test samples of an artery are obtained from the anisotropic plane and are loaded as illustrated in Figure 3. The objective function is defined from the standardized quadratic errors of the mechanical tests.
J ( y , y ^ ) = 1 n i = 1 n ( y i y i ^ ) 2 | m a x ( y ) m i n ( y ) |
where n is the number of experimental points, y are the experimental values, and y ^ are the values predicted by the model. With this metric, it is possible to define the following objective function:
min x A f ( x ) = J ( σ c i r , σ ^ c i r ) + J ( σ l o n , σ ^ l o n )
where σ c i r and σ l o n are the uniaxial stresses from the circumferential and longitudinal tensile test, respectively. The analytical equation of the uniaxial tensile test is Equation( 4) and is important to remember that this equation is conditioned to the direction of traction. The objective function carries out the key task of linking the mathematical model with the mechanical response of the material.

2.3. Transversal Isotropy Stabilization Criterion

Most anisotropic hyperelastic materials are proposed based on a decomposition of the density energy function between the isotropy and the anisotropy [6]. One of the main challenges for the characterization of anisotropic materials is to correctly determine the proportion of the strain energy from the matrix (isotropy) or the fibers (anisotropy). If the energy of the preferential directions is unbalanced concerning the isotropic energy, then a nonphysical behavior may occur. Furthermore, it is possible to obtain a similar strain energy function with different proportions of anisotropy in some conditions, leading to a similar stress behavior but with different deformation modes.
To illustrate the non-physical behavior we present two sets of parameters x = ( μ , κ , k 1 , k 2 , θ ) of the GHO constitutive model. The first set of parameters is x 1 = ( 10.4 [ kPa ] , 0.27 , 30.6 [ kPa ] , 3.94 , 61.3 ) and the second one is x 2 = ( 0.73 [ kPa ] , 0 , 13.9 [ kPa ] , 2.13 , 48.1 ) . If we plot the longitudinal and the circumferential behavior of these parameters, we can see in Figure 4a that they have similar stress responses. However, if we see the transversal stretch λ 3 in Figure 4, it is possible to observe that the second set of parameters has a non-physical behavior because it tends to grow under traction.

2.3.1. Uniaxial Transversal Stretch

The vast majority of the anisotropic hyperelastic models propose an additive decomposition between the isotropic and anisotropic energy strain functions. When the experimental data are fitted with the analytical Equations (4) and (5), it is only possible to fit the stress response, but it is not verified whether the balance between the isotropic and anisotropic energy density functions are right.
One way to quantify the anisotropy of the material is through the uniaxial tensile test (see Figure 1). The main stretch λ 1 and the force exerted on the specimen F are measured; however, the transversal stretches λ 2 , or λ 3 are not. If one of these stretches is quantified, then the balance between the isotropy and anisotropy of the strain energy function could be determined. We will explain next how to quantify the anisotropy with the transverse stretches. First, the fibers will be oriented according:
M = c o s ( Θ 1 ) E ^ 1 + s i n ( Θ 1 ) E ^ 2 ; M = c o s ( Θ 2 ) E ^ 1 + s i n ( Θ 2 ) E ^ 2
The variables Θ 1 and Θ 2 are the direction of the fibers in the anisotropy plane in the material configuration. Employing the deformation Gradient (3) the deformed fiber vectors are obtained.
FM = m = λ 1 c o s ( θ 1 ) e ^ 1 + λ 2 s i n ( θ 1 ) e ^ 2 ; FM = m = λ 1 c o s ( θ 2 ) e ^ 1 + λ 2 s i n ( θ 2 ) e ^ 2
The variables θ 1 and θ 2 stand for the fiber orientations in the spatial configuration and are illustrated in Figure 1. If we consider that the strain energy function depends on W ( I 1 , I 4 , I 6 ) and we replace Equation (11) in (7), we can obtain that:
0 = ( λ 2 2 λ 1 2 λ 2 2 ) W 1 + [ λ 2 s i n ( θ 1 ) ] 2 W 4 + [ λ 2 s i n ( θ 2 ) ] 2 W 6
Grouping λ 2 from the previous equation, the following equation is obtained:
λ 2 = [ W 1 + s i n 2 ( θ 1 ) W 4 + s i n 2 ( θ 2 ) W 6 W 1 ] 1 / 4 λ 1 1 / 2 = Ψ 1 / 4 λ 1 1 / 2
In order to simplify the expression of Equation (13) the variable Ψ was introduced. Nonetheless, Ψ hides a significant physical meaning which reveals the anisotropy balance through the rigidity ratio between the anisotropic plane and the transversal isotropy. In other words, Ψ allows quantifying the degree of anisotropy of the material. The value of Ψ is always 1 , and if Ψ = 1 , it means that the material remains isotropic.
Because the material is considered incompressible, it must be fulfilled that λ 3 = 1 λ 1 λ 2 . Therefore, Equation (13) can be written in terms of λ 3 as:
λ 3 = Ψ 1 / 4 λ 1 1 / 2
Thus, characterizing a material based on the transverse stretches makes it possible to know the balance between the isotropic and anisotropic energies, as long as the angle of the fibers is known.

2.3.2. Stability Criterion

We will call an incompressible and transversely isotropic material stable under uniaxial stress if the stretch aligned with the transversal isotropy direction decrease. If we consider that the deformation gradient is homogeneous for the case of uniaxial traction as in Equation (3), it is conceivable the definition of a transverse stability criterion:
d λ 3 d λ 1 0 Transversal stability criterion
The strain energy function W ( I 1 , I 4 , I 6 ) can be expressed with the principal stretches.
W ( I 1 , I 4 , I 6 ) = W ( λ 1 , λ 2 , λ 3 )
Since the material is incompressible and λ 3 = 1 / ( λ 1 λ 2 ) , then the energy density function can be written as:
W ( I 1 , I 4 , I 6 ) = W ( λ 1 , λ 2 )
The variables θ 1 and θ 2 determine the direction of the fibers in the anisotropy plane. The deformed fiber vector comes from Equation (11) and the invariants ( I 1 , I 4 , I 6 ) are calculated:
I 1 = λ 1 2 + λ 2 2 + ( λ 1 λ 2 ) 2 I 4 = λ 1 2 c o s 2 ( θ 1 ) + λ 2 2 s i n 2 ( θ 1 ) I 6 = λ 1 2 c o s 2 ( θ 2 ) + λ 2 2 s i n 2 ( θ 2 )
With these equations proposed, it is possible to correlate the constitutive model with the principal stretches. From Equation (14), the transverse elongation λ 3 can be written as:
λ 3 = Ψ 1 / 4 λ 1 1 / 2 = [ W 1 + s i n 2 ( θ 1 ) W 4 + s i n 2 ( θ 2 ) W 6 W 1 ] λ 1 1 / 2
Accordingly, if we replace Equation (19) in Equation (20), we get that:
d λ 3 d λ 1 = 1 4 Ψ 3 / 4 λ 1 1 / 2 d Ψ d λ 1 1 2 Ψ 1 / 4 λ 1 3 / 2 0
From the policonvexivity conditions of Walton and Wilber [44], we know that W k > 0 for k = ( 1 , 2 , 4 ) and hence Ψ 1 . We also know that λ 1 > 0 and the inequity (15) can be expressed as:
d Ψ d λ 1 2 Ψ λ 1
It is noteworthy to mention that if Ψ is constant for λ 1 > 1 , then the stability criterion is satisfied (remember that Ψ 1 ). It is relevant to note that Equation (21) states that transverse stability can be controlled through the growth rate of the material’s anisotropic energy. An alternative way to present Equation (21) is to express the derivative of Ψ in terms of λ 1 and λ 2 :
Ψ 1 [ λ 1 2 λ 2 2 ] + Ψ 4 [ λ 1 2 c o s 2 ( θ 1 ) λ 2 2 s i n 2 ( θ 1 ) ] + Ψ 6 [ λ 1 2 c o s 2 ( θ 2 ) λ 2 2 s i n 2 ( θ 2 ) ] Ψ .
Finally, an inequity is obtained that makes it possible to determine whether the material is transversally stable as a function of the invariants.

2.3.3. Transverse Stability of HGO and GHO Constitutive Models

This section shows that the GHO model [16] is transversely stable for an elongation λ 1 that tends to infinity and that the HGO model is not. GHO model is given by:
W ( I 1 , I 4 , I 6 ) = μ 2 ( I 1 3 ) + k 1 2 k 2 i = 4 , 6 ( exp ( k 2 ( I 1 κ + ( 1 3 κ ) I i 1 ) 2 ) 1 ]
In order to proof the stability of the model, we assume that there is a critical angle θ , which coincides with both preferred directions. Therefore, it is only necessary to consider the invariant I 4 . Moreover, it is possible to observe that Ψ depends on I 1 and I 4 :
Ψ ( I 1 , I 4 ) = W 1 ( I 1 , I 4 ) + s i n 2 ( θ ) W 4 ( I 1 , I 4 ) W 1 ( I 1 , I 4 )
By evaluating the partial derivatives of the energy function in Equation (24), then Ψ is:
Ψ ( I 1 , I 4 ) = 1 + k 1 E 4 ( 1 3 κ ) s i n 2 ( θ ) μ / 2 + e k 2 E 4 2 + k 1 κ E 4
and the value of E 4 is defined by:
E 4 = I 1 κ + ( 1 3 κ ) I 4 1
If λ 1 , then I 1 , I 4 and E 4 . Therefore, the value of Ψ when λ 1 tends to infinity is:
lim λ 1 Ψ = 1 + ( 1 3 κ ) s i n 2 ( θ ) κ ; κ 0
The value of κ must be between ( 0 , 1 / 3 ] . It is possible to observe that when κ = 0 , GHO model is identical to HGO model and exhibits transversal instability. However, if κ belongs to the interval of ( 0.1 / 3 ] , then Ψ converges to a positive real value when λ 1 . As the value of Ψ converges to a positive constant, then it is obtained that:
lim λ 1 d Ψ d λ 1 = 0 κ ( 0 , 1 / 3 ]
Consequently, with the results of (27) and (28), the inequity is evaluated:
lim λ 1 d Ψ d λ 1 lim λ 1 2 Ψ λ 1 0 0
As the inequity (29) is fulfilled, then the GHO model is transversally stable. In other words, there is a large enough λ x , that for all λ 1 > λ x the GHO model will be stable.
Remember that when the value of κ = 0 , then the model of GHO is equivalent to the model of HGO [42]. In this case, it is shown that the model of Holzapfel–Gasser–Ogden is transversally unstable.

2.3.4. Penalty and Objective Function for Transverse Stability

One of the easiest ways to obtain a set of transversely stable parameters is to consider in the fitting procedure the transverse stretches ( λ 2 or λ 3 ). However, measuring experimentally transverse stretches can be complex due to the experimental setup.
A numerical alternative is presented to define the feasible decision space, which is transversal stable in the range of the experimental data. Based on the inequity (21), a modified stability criterion is defined:
d Ψ d λ 1 2 Ψ λ 1 h ( λ 1 )
where h ( λ 1 ) is a function ∈ [0,1] and controls the degree of anisotropy of the material in the domain of the stable parameters, with this stability Criterion (30) a penalty term can be defined, and the following objective function can be set:
min x A n f ( x ) = J ( σ c i r , σ ^ c i r ) + J ( σ l o n , σ ^ l o n ) + g ( x , λ i l o n , λ i c i r )
where g represents the penalty function and is defined as:
g ( x , λ i l o n , λ i c i r ) = C , Si d Ψ d λ 1 > 2 Ψ λ 1 h ( λ 1 ) λ 1 λ i c i r , λ i l o n 0 , Otherwise
The penalty constant C is an immense value, penalizing the objective function when the transversal stability is not met. With the objective Function (31) it is possible to find a set of parameters that adjust the uniaxial traction data and satisfy the transversal stability condition.
The function h ( λ 1 ) can be simplified as a constant. A recommended value for the GHO model is h = 0.7 . If h = 1 , the set of parameters obtained may be critical, and it is advisable to search with a sub-critical constraint. On the other hand, if h = 0 , it will obtain a set of parameters that exhibit isotropic behavior. To characterize the materials throughout this work, h = 0.7 is employed.

2.4. Inverse Finite Element Characterization

The analytical model of the uniaxial tensile test, which considers a homogeneous deformation gradient, gives a response similar to the one obtained with numerical models (see [24]). However, the analytical model of the pressurization test is too simplistic and is limited to the thin-walled case [20]. In reality, the strain gradients are not homogeneous, and there is a non-linear stress distribution. Therefore, finite element simulations of the pressurization tests are used to calibrate the parameters based on numerical simulations. To illustrate the inverse finite element characterization procedure, experimental data from the mechanical behavior of lambs’ arteries will be used in conjunction with the GHO model. The experimental data comes from the work of Rivera et al. [45] and we will characterize only the CN group. It is important to remember that this work focuses on the characterization process and not on bio-mechanics. Therefore, the numerical simulation is presented in detail, followed by the definition of the objective function.

2.4.1. Pressurization Test and Numerical Simulation

The pressurization test is used to replicate In vivo blood vessel loading conditions. In this mechanical test, a sample of an artery is subjected to uniaxial elongation in a traction machine followed by the application of internal pressure employing a fluid (calcium-free saline solution) that travels through the interior of the blood vessel and generates a radial load. The ends of the cylindrical sample of the thoracic aorta are fixed to a metal nozzle. The experimental set-up is similar to the one described by Guinea et al. [46]. Figure 5 shows the experimental set-up used in the work of Rivera et al.
The aorta’s internal pressure and external diameter are recorded during the entire mechanical test using video; these data are processed to obtain the internal pressure vs. diametrical stretch curves. The circumferential (or diametrical) elongation is defined as D / D 0 , where D and D 0 denote the deformed diameter and the undeformed diameter, respectively. Ten loading cycles are performed up to a pressure of 170 [mmHg] to precondition the samples. The last cycle is used to perform the test analysis.
The numerical simulation to be used consists of an axisymmetric 2D simulation using finite elements because it idealizes the geometry of the aorta as a perfect cylinder and simplifies the problem with symmetry. The simulation is done through two steps:
The first step consists of subjecting the aorta to uniaxial traction until a physiological lengthening is achieved. In Rivera’s work, they determine that the value of λ z = 1.2 , this lengthening remains constant during the whole pressurization process. The second step corresponds to the process of pressurization (inflation), where internal pressure is applied to the arterial wall that covers the physiological range up to hypertension; that is, the pressure applied to the arterial wall varies from 0 [mmHg] to 170 [mmHg]. Figure 6a represents the pressurization test.
The shape of the artery is simplified as a rectangular surface by making a longitudinal cut to the cylindrical tube. Each dimension of the artery is obtained through the average of the measurements of the samples available. The dimensions of the internal radius, thickness, and length are (CN: R 0 = 3.71 mm, t = 2.06 mm, L = 14 mm )
The mesh used has 2626 nodes and 2500 quadrilateral elements (25 × 100). The mesh presents a refinement at the base of the artery, where it is fixed (see Figure 6b). With this numerical simulation, the circumferential elongation can be determined as a function of the parameters of the GHO model x and the internal pressure P as
λ θ ( x , P ) = 1 + u r ( x , P ) D 0 ,
where u r is the radial displacement of the node located at the outer diameter at the edge of the longitudinal symmetry (red dot Figure 6a), with this expression, the circumferential stretch of the objective function to be presented in the following section is evaluated.

2.4.2. Inverse Calibration Procedure

The procedure to characterize the material is straightforward; an objective function is defined, combining the uniaxial analytical test with the numerical pressurization test. This objective function is as follows.
min x A f ( x ) = J ( σ c i r , σ ^ c i r ) + J ( σ l o n , σ ^ l o n ) + 2 J ( λ θ p r e s u , λ ^ θ p r e s u ) + g 1 ( x , λ i l o n , λ i c i r )
It is possible to observe that each curve presents a standardized error, and that the pressurization test is weighted by a factor of two so that it has the same weight as the two tensile tests. This objective function is numerically calculated following the block diagram presented in Figure 7. Initially, the quadratic errors of the longitudinal and circumferential tensile tests are obtained, and it is determined whether the set of parameters exhibits transverse stability in the experimental range. If this restriction is not satisfied, the objective function is penalized and takes a substantial value. If the parameters are not stable, there is no need to perform the simulation and save computation time. Similarly, if the pressurization simulation does not converge for a given set of parameters, the value of the target function is penalized by an immense value. If the set of parameters of the GHO model is stable in the experimental range and the pressurization simulation converges, then the objective function is calculated.
The objective function determines the optimization problem; it can be calculated following the block diagram in Figure 7 and it is essential to recognize that the simulation is evaluated at the end of the block to save computation time. The objective function is subject to stability regions, multiple parameters, and discontinuities associated with the convergence of the simulation and, therefore, it is necessary to use some optimizer capable of solving such a problem. This is where the evolutionary strategies come in.

2.5. Evolutionary Strategies

Evolution Strategies (ES) belong to the class meta-heuristic optimization methods and have proven successful for the solution optimization problems in continuous spaces. Furthermore, ES allows to solve problems in which, for example, there is not explicit functional expression or derivatives that can be calculated. ES imitates the biological principle of evolution and serves as an approach to machine learning and global optimization methods. These strategies are based on three mechanisms of the Darwinian evolutionary process: selection, recombination, and mutation.
The optimization problem presented in Equation (34) is non-linear and is subject to multiple constraints and a discontinuous space. ES has proven to be a successful solution to this type of inverse problem [47,48], in which it is not necessary to satisfy continuity and calculate derivatives. For these reasons, they are used in this work. For a more detailed depiction of how this algorithm works, we refer to the interested reader to references [49,50]. Nevertheless, a detailed block diagram illustrating the link between the physics, the model, and the optimizer is presented in Figure 8.
As shown in Figure 8, the optimization process begins with the creation of a set of potential solutions to the problem studied. Since ES mimics the evolution of a population through time, each element of the initial set is modified, in the variation stage, through the application of the recombination and mutation operators. The former serves to exchange the genetic information between the elements of the initial set. The latter, on the other hand, introduces some random variation that helps the ES to explore the search space. The throughput of the variation stage is a secondary set of modified solutions that will be evaluated, through some predefined metric that relates the optimization algorithm with the physical problem, and only those elements whose fitness values may lead the ES to the optimal solution will be retained to constitute the initial population for the next iteration of the optimization loop which, as shown in Figure 8, will end once the termination criterion has been satisfied. In the present work, we denote x for each element of the initial set. Furthermore, we make use of the so-called Elitist Evolution Strategy, which selects from both the initial and the modified sets only those elements with the best fitness values. It is noteworthy to mention that there is another selection scheme known as Non-Elitist Evolution Strategy. On the contrary to the Elitist Evolution Strategy, it only selects the best elements from the modified set. The results of extensive numerical experiments showed that the Non-Elitist Strategy was not the most suitable option for the present problem.

3. Results

3.1. Assessment of Evolution Strategies

The characterization of an anisotropic material is often more complex than that of an isotropic material. The main difference between these two problems is that the anisotropy of the material must be captured, and achieved with more complex hyperelastic models and more mechanical tests. Although ES has been applied for several optimization problems, it is unknown how efficient and robust they are when used to characterize the constitutive models of hyperelastic anisotropic materials.
A benchmark test is performed to assess the suitability of ES to characterize an anisotropic material. This test consists of finding the characteristic parameters of the material through the simulation of experimental data of an artery under traction. In other words, the GHO model is used with a set of known parameters to generate synthetic experimental data curves. Generally, when the experimental tests are carried out, longitudinal and circumferential tensile tests are performed.
The equations that model the mechanical behavior of the material are given by Equations (4) and (5). The constitutive model of GHO [16] is used to generate the simulated experimental curves, consisting of the parameters x = ( μ , κ , k 1 , k 2 , θ ) . The objective is to generate the curves with a set of known parameters and verify the possibility of retrieving the parameters by making use of the ES and the objective function. Therefore, the solution to the optimization problem (Equation (9)) has to match the set of known parameters.
The set of known parameters is presented in Table 1. The search domain and the parameters of the evolutionary strategies are also detailed. The solution obtained from the evolutionary algorithm involves stochastic processes, and it is possible that in some runs of the optimization algorithm convergence is achieved to the global optimum and in others runs not. Therefore, 100 realizations of the optimization problem are made, and a statistical analysis of 100 solutions are obtained.
The curves generated with the known parameters and the adjustment obtained with the evolutionary strategies are illustrated in Figure 9. It can be seen that the set of parameters obtained correctly adjusts the two curves. In this case, the average of the 100 realizations coincides with the problem’s solution and presents a minimal standard deviation (see Table 2). Therefore, in this case, all the solutions obtained are consistent with the global optimum, showing the resolution capabilities of the ES.
From Table 2, the elitist strategy exhibits better results than the non-elitist strategy. In this case, both variants find the global optimum, but the elitist strategy converges more consistently.
One way to illustrate that the set of parameters obtained is the right one is to compare the isocontours of the energy density function W ( λ 1 , λ 2 , λ 3 ) of the exact solution and the solution obtained through the adjustment. Since the material is incompressible, one principal stretch could be determined by the others, and the energy density function can be expressed as W ( λ 1 , λ 2 ) . In Figure 10, it can be seen that the isocontours are superimposed and that the energy density function is practically the same.
It is important to notice that the green curve and the cyan curve illustrate how the stretches behave in the uniaxial tensile test. The path of this curve is not constant as in the isotropic case, and it is subject to the solution of a non-linear equation. Another equivalent way to determine the second stretch is with the following equation:
W λ 2 = 0
This equation can be obtained through the constitutive Model (23), and it can be seen that it is satisfied in Figure 10. In other words, it is verified that the gradient is zero at the intersection of the energy isocontours. In the circumferential curve, it is also satisfied that the gradient W λ 1 = 0 , the only difference is that the principal stretch is λ 2 and the second principal stretch is λ 1 .
This section has shown that ES can consistently characterize an anisotropic hyperelastic material using the GHO model [16]. It is noteworthy to mention that the success of the ES depends on the search parameters and the complexity of the problem to be solved. Nevertheless, it is essential to recall that adjusting uniaxial stress data does not guarantee a coherent physical response.

3.2. Characterization of Experimental Data with Stabilization Criterion

In this section, evolutionary strategies are used with the stabilization criterion to obtain the constitutive parameters of the GHO model with the experimental data published in the work of García et al. [51]. In García’s work, experiments and modeling of the passive mechanical response of the descending thoracic aorta of humans are carried out. For this purpose, uniaxial traction tests were performed on healthy samples of the arteries of newborns, young people, and adults. Subsequently, data from the uniaxial tensile test were used by García to calibrate HGO model (2000) with Levenberg Marquardt optimization algorithm.
There are three groups of experimental data; Group A comes from newborns, Group B is conformed by young people, Group C is adult arteries. Each of these groups represents a set of parameters. Each group tensile test represents the average behavior of that group of several individuals. However, we only characterize Group C which is the most unstable and will not detail the individuals or the dispersion of the data, and will only use the average curve.
Therefore, with the data published by Garcia of the uniaxial traction test, we proceed to calibrate the parameters of GHO model x 1 = ( μ , κ , k 1 , k 2 , θ ) with the evolutionary strategies and the stabilization criterion with the objective Function (34).
The reported parameters of Garcia were obtained using the Levenberg Marquardt method and the HGO model. In Table 3, these parameters are written in terms of the GHO model, considering that κ = 0 . Furthermore, the parameters obtained are presented with the evolutionary strategies, stabilization term and the GHO model. The value of the objective function determines the degree of fitness of the set of parameters obtained. With the proposed methodology the objective function is 32.7 times lower than the parameters retrieved by García and these parameters satisfy the stability condition. The confidence interval of the evolutionary strategies is reasonably narrow and therefore the solution obtained is consistent.
The experimental fit of the uniaxial traction and the traversal stretches are presented in Figure 11. It is clear that the set of parameters obtained adjusts better the data and satisfy the stability condition. In contrast, the results presented by Garcia are inherently unstable for large stretches since it was shown that the Holzapfel–Gasser–Ogden model will always be unstable for elongations tending to infinity.
Finally, using the ES, with the GHO model and the proposed penalty function, it is possible to obtain parameters that are transversally stable and an excellent fit of the experimental data.

3.3. Inverse Finite Element Characterization

In this procedure, arterial tissue is characterized using experimental data from three mechanical tests: two uniaxial tension and one of pressurization. The experimental data is fitted in two different ways, with and without the penalty function for transverse stability. This procedures characterize the material considering the FEM simulation of the pressurization test and the analytical equation for the tensile test. This task is achieved by minimizing Equation (34) with the evolutionary strategies.
Some of the difficulties of calculating the objective function with finite element simulations is the associated computational cost. In this case, the aim is to reduce as much as possible the number of evaluations of the objective function. The search parameters of the evolutionary strategy are presented in Table 4.
The numerical implementation of ES is compiled with GNU FORTRAN and is parallelized with the OPENMP library. In order for the simulations to reduce the possibility of having problems related to the convexity of the energy density function, the minimum value of μ is 10 [kPa]. This value guarantees an isotropic base component, which helps to stabilize the pressurization finite element simulations.
In order to study the influence of the stabilizing term, the arterial tissue is characterized using the inverse methodology with and without the stabilizer. The results of the adjustments are shown in Figure 12, where it is possible to appreciate that the parameters obtained adequately adjust each of the mechanical tests and that they belong to the confidence interval. The standardized quadratic errors of the mechanical tests are presented in Table 5 and in each of the curves presents an error that is less or equal to 1 × 10 4 . When the model does not consider the stability constraint, the obtained parameters correctly capture the mechanical tests but with a transversely unstable set of parameters since the transversal stretch λ 3 increases. However, when considering the stabilization term, the transversal stretch λ 3 decreases and a stable solution is obtained that accurately adjusts each mechanical test. It is essential to notice that a nonlinear constraint defines the decision space of this function. Therefore evolutionary strategies are an appropriate tool to find the solution in that space.
There is a fundamental relationship between material stability and the balance between isotropy and anisotropy. A parameter that controls the degree of anisotropy of the material in the GHO model is κ , if this is equal to zero, then all fibers are aligned in the preferential directions, and if κ = 1 / 3 then the fiber arrangement is uniform and therefore the behavior of the material is isotropic. From Table 6, it can be seen that the anisotropic component is more significant in the unstable set of parameters than in the stable case, simply by looking at the value of κ . Furthermore, the constant μ is proportional to the material’s isotropic energy, which is higher when the material is stable. However, it is essential to notice that the angle obtained is insensitive to the stability constraint.
Finally, it is demonstrated that this inverse methodology can characterize soft tissue that is constrained by multiple mechanical tests and a stability domain with FEM.

4. Discussion

This work has focused on the characterization of anisotropic hyperelastic materials and on stabilizing the mechanical response associated with the constitutive parameters, considering the HGO and GHO models as examples. The analysis performed can be summarized in three parts: (i) analysis of the stability problem; (ii) application of ES in anisotropic hyperelasticity; (iii) inverse FEM characterization procedure. These points are proposed to find the material constitutive parameters that are stable, model the deformation gradients as inhomogeneous, and that the optimization algorithm is oriented to global optimization.
From the first point is derived the definition of a new stability criterion that allows obtaining a physically expected mechanical behavior. This criterion establishes that if the material is being subjected to uniaxial tension, the elongation normal to the plane of the fibers must decrease.
The contribution of this criterion is important since the instability of a material occurs due to the excessive growth of the anisotropic stiffness with respect to the transverse isotropic stiffness. This imbalance imposes that the transverse elongation has to grow to satisfy static equilibrium and quasi-incompressibility.
In this study, it shows that the HGO model is unconditionally unstable for large deformations, since, the energy density function grows exponentially for anisotropic invariants and linearly for isotropic ones. On the other hand, the GHO model presents a stability domain where it is stable for large deformations because the energy density function grows exponentially for isotropic and anisotropic invariants, keeping a certain balance.
In the second point, an ES were applied to obtain constitutive parameters of an anisotropic hyperelastic material. To evaluate the suitability of evolutionary strategies a benchmark problem with a known solution was defined with synthetic experimental data of arteries in the circumferential and longitudinal direction. The results show that the application of an ES successfully found the global optimum in each of the one hundred realizations of the optimization problem for both elitist and non-elitist selection. Therefore, the effectiveness of this algorithm to characterize anisotropic hyperelastic materials, where the nonlinearity is large, with multiple constraints and lack of an analytical expression of the gradients of the objective function, was verified.
From the third point, robust characterization of anisotropic hyperelastic materials was obtained. The stability criterion is included in the characterization through a penalty function that defines the domain of the stable parameters. Evolutionary strategies play a fundamental role here, because with this algorithm it is only necessary to sum this penalty term to the objective function. After characterizing experimental data of arteries found in the literature, it is demonstrated that this technique provides a superior fit to the one reported using Levenberg–Marquardt and, in addition, it presents a transversely stable mechanical response. Due to the nature of this optimization algorithm, it cannot be guaranteed that the global optimum of the problem will be found, but a stable and good candidate set of parameters will be found.
Finally, this work shows that the correct characterization of material is conditioned by the experimental data, the optimization algorithm, the stability of the mathematical model and simplifications made to the deformation kinematic. The inverse characterization procedure with finite elements, the stabilization term and the evolutionary strategies; proposes a solution for each of these points, since it has been shown to effectively characterize the mechanical behavior of complex materials, such as arteries, considering non-homogeneous deformation gradients, multiple mechanical tests, stability constraints and an optimization algorithm oriented to the global optimization. The major drawback of this procedure is the computational cost since it is necessary to evaluate the objective function several times. However, model order reduction techniques could be a potential solution to this inconvenience.

Author Contributions

C.C., C.G.-H., E.R., D.M. and D.C. developed the concept and design of the work. C.C. and C.G.-H. developed the stability criterion. C.C., C.G.-H. and E.R. performed postprocessing of data and numerical simulation. C.C., C.G.-H. and E.R. analyzed experimental and computational results. C.C., C.G.-H., E.R., D.M. and D.C. provided suggestions and editing assistance. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the Chilean National Agency of Research and Development (ANID) through project FONDECYT No. 1220956 and the Faculty of Engineering (FING) USACH.

Data Availability Statement

Not applicable.

Acknowledgments

The authors are grateful for the support of the Chilean National Agency of Research and Development (ANID) through project FONDECYT No. 1220956 and the Faculty of Engineering (FING) USACH.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Hernandez, J.A.; Al-Qadi, I.L. Tire–pavement interaction modelling: Hyperelastic tire and elastic pavement. Road Mater. Pavement Des. 2017, 18, 1067–1083. [Google Scholar] [CrossRef]
  2. Abubakar, I.; Myler, P.; Zhou, E. Constitutive modelling of elastomeric seal material under compressive loading. Model. Numer. Simul. Mater. Sci. 2016, 6, 28–40. [Google Scholar] [CrossRef]
  3. Ucar, H.; Basdogan, I. Dynamic characterization and modeling of rubber shock absorbers: A comprehensive case study. J. Low Freq. Noise Vib. Act. Control 2018, 37, 509–518. [Google Scholar] [CrossRef]
  4. Goury, O.; Duriez, C. Fast, generic, and reliable control and simulation of soft robots using model order reduction. IEEE Trans. Robot. 2018, 34, 1565–1576. [Google Scholar] [CrossRef]
  5. Park, Y.L.; Majidi, C.; Kramer, R.; Bérard, P.; Wood, R.J. Hyperelastic pressure sensing with a liquid-embedded elastomer. J. Micromech. Microeng. 2010, 20, 125029. [Google Scholar] [CrossRef]
  6. Chagnon, G.; Rebouah, M.; Favier, D. Hyperelastic energy densities for soft biological tissues: A review. J. Elast. 2015, 120, 129–160. [Google Scholar] [CrossRef]
  7. Nolan, D.R.; Gower, A.L.; Destrade, M.; Ogden, R.W.; McGarry, J. A robust anisotropic hyperelastic formulation for the modelling of soft tissue. J. Mech. Behav. Biomed. Mater. 2014, 39, 48–60. [Google Scholar] [CrossRef]
  8. Avril, S.; Badel, P.; Duprey, A. Anisotropic and hyperelastic identification of in vitro human arteries from full-field optical measurements. J. Biomech. 2010, 43, 2978–2985. [Google Scholar] [CrossRef]
  9. García-Herrera, C.M.; Celentano, D.J.; Cruchaga, M.A. Bending and pressurisation test of the human aortic arch: Experiments, modelling and simulation of a patient-specific case. Comput. Methods Biomech. Biomed. Eng. 2013, 16, 830–839. [Google Scholar] [CrossRef]
  10. Veselỳ, J.; Hornỳ, L.; Chlup, H.; Adámek, T.; Krajíček, M.; Žitnỳ, R. Constitutive modeling of human saphenous veins at overloading pressures. J. Mech. Behav. Biomed. Mater. 2015, 45, 101–108. [Google Scholar] [CrossRef]
  11. Whitford, C.; Movchan, N.V.; Studer, H.; Elsheikh, A. A viscoelastic anisotropic hyperelastic constitutive model of the human cornea. Biomech. Model. Mechanobiol. 2018, 17, 19–29. [Google Scholar] [CrossRef] [PubMed]
  12. Martins, P.; Peña, E.; Calvo, B.; Doblaré, M.; Mascarenhas, T.; Natal Jorge, R.; Ferreira, A. Prediction of nonlinear elastic behaviour of vaginal tissue: Experimental results and model formulation. Comput. Methods Biomech. Biomed. Eng. 2010, 13, 327–337. [Google Scholar] [CrossRef] [PubMed]
  13. Nemavhola, F.; Ngwangwa, H.M.; Pandelani, T. An investigation of uniaxial mechanical properties of excised sheep heart muscle fibre–fitting of different hyperelastic constitutive models. Preprints 2021, 2021080566. [Google Scholar]
  14. Annaidh, A.N.; Bruyère, K.; Destrade, M.; Gilchrist, M.D.; Otténio, M. Characterization of the anisotropic mechanical properties of excised human skin. J. Mech. Behav. Biomed. Mater. 2012, 5, 139–148. [Google Scholar] [CrossRef] [PubMed]
  15. Holzapfel, A.G. Nonlinear Solid Mechanics II; Wiley: Hoboken, NJ, USA, 2000. [Google Scholar]
  16. Gasser, T.C.; Ogden, R.W.; Holzapfel, G.A. Hyperelastic modelling of arterial layers with distributed collagen fibre orientations. J. R. Soc. Interface 2006, 3, 15–35. [Google Scholar] [CrossRef]
  17. Dodson, R.B.; Martin, J.T.; Hunter, K.S.; Ferguson, V.L. Determination of hyperelastic properties for umbilical artery in preeclampsia from uniaxial extension tests. Eur. J. Obstet. Gynecol. Reprod. Biol. 2013, 169, 207–212. [Google Scholar] [CrossRef]
  18. Labus, K.M.; Puttlitz, C.M. An anisotropic hyperelastic constitutive model of brain white matter in biaxial tension and structural–mechanical relationships. J. Mech. Behav. Biomed. Mater. 2016, 62, 195–208. [Google Scholar] [CrossRef]
  19. Tonge, T.K.; Atlan, L.S.; Voo, L.M.; Nguyen, T.D. Full-field bulge test for planar anisotropic tissues: Part I–Experimental methods applied to human skin tissue. Acta Biomater. 2013, 9, 5913–5925. [Google Scholar] [CrossRef]
  20. Holzapfel, G.A.; Ogden, R.W. Modelling the layer-specific three-dimensional residual stresses in arteries, with an application to the human aorta. J. R. Soc. Interface 2010, 7, 787–799. [Google Scholar] [CrossRef]
  21. Rivera, E.; Canales, C.; Alarcón, M.P.; García-Herrera, C.; Macías, D.; Celentano, D.J.; Herrera, E.A. Biomechanical characterization of the passive response of the thoracic aorta in chronic hypoxic newborn lambs using an evolutionary strategy. Sci. Rep. 2021, 11, 13875. [Google Scholar] [CrossRef]
  22. Ding, Y.; Xu, G.K.; Wang, G.F. On the determination of elastic moduli of cells by AFM based indentation. Sci. Rep. 2017, 7, 45575. [Google Scholar] [CrossRef] [PubMed]
  23. Chanda, A.; Chatterjee, S.; Gupta, V. Soft composite based hyperelastic model for anisotropic tissue characterization. J. Compos. Mater. 2020, 54, 4525–4534. [Google Scholar] [CrossRef]
  24. Sasso, M.; Palmieri, G.; Chiappini, G.; Amodio, D. Characterization of hyperelastic rubber-like materials by biaxial and uniaxial stretching tests based on optical methods. Polym. Test. 2008, 27, 995–1004. [Google Scholar] [CrossRef]
  25. Íñiguez-Macedo, S.; Lostado-Lorza, R.; Escribano-García, R.; Martínez-Calvo, M.Á. Finite element model updating combined with multi-response optimization for hyper-elastic materials characterization. Materials 2019, 12, 1019. [Google Scholar] [CrossRef]
  26. Lapeer, R.; Gasson, P.; Karri, V. Simulating plastic surgery: From human skin tensile tests, through hyperelastic finite element models to real-time haptics. Prog. Biophys. Mol. Biol. 2010, 103, 208–216. [Google Scholar] [CrossRef]
  27. Gendy, A.; Saleeb, A. Nonlinear material parameter estimation for characterizing hyper elastic large strain models. Comput. Mech. 2000, 25, 66–77. [Google Scholar] [CrossRef]
  28. Ogden, R.W.; Saccomandi, G.; Sgura, I. Fitting hyperelastic models to experimental data. Comput. Mech. 2004, 34, 484–502. [Google Scholar] [CrossRef]
  29. Martins, P.; Natal Jorge, R.; Ferreira, A. A comparative study of several material models for prediction of hyperelastic properties: Application to silicone-rubber and soft tissues. Strain 2006, 42, 135–147. [Google Scholar] [CrossRef]
  30. Fernández, J.; López-Campos, J.; Segade, A.; Vilán, J. A genetic algorithm for the characterization of hyperelastic materials. Appl. Math. Comput. 2018, 329, 239–250. [Google Scholar] [CrossRef]
  31. López-Campos, J.; Ferreira, J.; Segade, A.; Fernández, J.; Natal, R. Characterization of hyperelastic and damage behavior of tendons. Comput. Methods Biomech. Biomed. Eng. 2020, 23, 213–223. [Google Scholar] [CrossRef]
  32. Ramzanpour, M.; Hosseini-Farid, M.; Ziejewski, M.; Karami, G. A constrained particle swarm optimization algorithm for hyperelastic and visco-hyperelastic characterization of soft biological tissues. Int. J. Comput. Methods Eng. Sci. Mech. 2020, 21, 169–184. [Google Scholar] [CrossRef]
  33. Mansouri, M.R.; Darijani, H.; Baghani, M. On the correlation of FEM and experiments for hyperelastic elastomers. Exp. Mech. 2016, 57, 195–206. [Google Scholar] [CrossRef]
  34. Helfenstein, J.; Jabareen, M.; Mazza, E.; Govindjee, S. On non-physical response in models for fiber-reinforced hyperelastic materials. Int. J. Solids Struct. 2010, 27, 2056–2061. [Google Scholar] [CrossRef]
  35. Latorre, M.; Montáns, F.J. What-you-prescribe-is-what-you-get orthotropic hyperelasticity. Comput. Mech. 2014, 53, 1279–1298. [Google Scholar] [CrossRef]
  36. Marckmann, G.; Verron, E. Comparison of hyperelastic models for rubber-like materials. Rubber Chem. Technol. 2006, 79, 835–858. [Google Scholar] [CrossRef]
  37. Bonet, J.; Wood, R.D. Nonlinear Continuum Mechanics for Finite Element Analysis; Cambridge University Press: Cambridge, UK, 2008. [Google Scholar]
  38. Spencer, A.J.M. Constitutive theory for strongly anisotropic solids. In Continuum Theory of the Mechanics of Fibre-Reinforced Composites; Springer: Nottingham, UK, 1984; pp. 1–32. [Google Scholar]
  39. Shariff, M.; Bustamante, R. On the independence of strain invariants of two preferred direction nonlinear elasticity. Int. J. Eng. Sci. 2015, 97, 18–25. [Google Scholar] [CrossRef]
  40. Lin, D.; Yin, F. A multiaxial constitutive law for mammalian left ventricular myocardium in steady-state barium contracture or tetanus. J. Biomech Eng. 1998, 120, 504–517. [Google Scholar] [CrossRef]
  41. Holzapfel, G.A.; Ogden, R.W. Constitutive modelling of arteries. Proc. R. Soc. A Math. Phys. Eng. Sci. 2010, 466, 1551–1597. [Google Scholar] [CrossRef]
  42. Holzapfel, G.A.; Gasser, T.C.; Ogden, R.W. A new constitutive framework for arterial wall mechanics and a comparative study of material models. J. Elast. Phys. Sci. Solids 2000, 61, 1–48. [Google Scholar]
  43. Ogden, R.W. Nonlinear elasticity, anisotropy, material stability and residual stresses in soft tissue. In Biomechanics of Soft Tissue in Cardiovascular Systems; Springer: Glasgow, UK, 2003; pp. 65–108. [Google Scholar]
  44. Walton, J.R.; Wilber, J.P. Sufficient conditions for strong ellipticity for a class of anisotropic materials. Int. J. Non-Linear Mech. 2003, 38, 441–455. [Google Scholar] [CrossRef]
  45. Rivera, E.; García-Herrera, C.; González-Candia, A.; Celentano, D.J.; Herrera, E. Effects of melatonin on the passive mechanical response of arteries in chronic hypoxic newborn lambs. J. Mech. Behav. Biomed. Mater. 2020, 112, 104013. [Google Scholar] [CrossRef]
  46. Atienza, J.M.; Guinea, G.V.; Rojo, F.J.; Burgos, R.J.; García-Montero, C.; Goicolea, F.J.; Aragoncillo, P.; Elices, M. The Influence of Pressure and Temperature on the Behavior of the Human Aorta and Carotid Arteries. Rev. Esp. Cardiol. 2007, 3, 259–267. [Google Scholar] [CrossRef]
  47. Yang, J.M.; Chen, Y.P.; Horng, J.T.; Kao, C.Y. Applying family competition to evolution strategies for constrained optimization. In Proceedings of the International Conference on Evolutionary Programming, Indianapolis, IN, USA, 13–16 April 1997; pp. 201–211. [Google Scholar]
  48. Mezura-Montes, E.; Coello, C.A.C. An empirical study about the usefulness of evolution strategies to solve constrained optimization problems. Int. J. Gen. Syst. 2008, 37, 443–473. [Google Scholar] [CrossRef]
  49. Macías, D.; Vial, A.; Barchiesi, D. Application of evolution strategies for the solution of an inverse problem in near-field optics. JOSA A 2004, 21, 1465–1471. [Google Scholar] [CrossRef] [PubMed]
  50. Beyer, H.G. The Theory of Evolution Strategies; Springer Science & Business Media: Dortmund, Germany, 2001. [Google Scholar]
  51. García-Herrera, C.M.; Celentano, D.J.; Cruchaga, M.A.; Rojo, F.J.; Atienza, J.M.; V, G.G.; Goicolea, J.M. Mechanical characterisation of the human thoracic descending aorta: Experiments and modelling. Comput. Methods Biomech. Biomed. Eng. 2012, 15, 185–193. [Google Scholar] [CrossRef]
Figure 1. Uniaxial tension test of anisotropic two fibers family in spatial configuration. No symmetry is implied.
Figure 1. Uniaxial tension test of anisotropic two fibers family in spatial configuration. No symmetry is implied.
Mathematics 11 00922 g001
Figure 2. Longitudinal symmetry of the fiber family in the material configuration.
Figure 2. Longitudinal symmetry of the fiber family in the material configuration.
Mathematics 11 00922 g002
Figure 3. Graphical representation of the test tubes obtained in an artery.
Figure 3. Graphical representation of the test tubes obtained in an artery.
Mathematics 11 00922 g003
Figure 4. Material responses obtained with unstable and stable model parameters. (a) Tensile test; (b) Transversal Stretch λ 3 .
Figure 4. Material responses obtained with unstable and stable model parameters. (a) Tensile test; (b) Transversal Stretch λ 3 .
Mathematics 11 00922 g004
Figure 5. Experimental assembly. (a) Unstretched thoracic aorta (b) Pressurized thoracic aorta with a pre-stretch.
Figure 5. Experimental assembly. (a) Unstretched thoracic aorta (b) Pressurized thoracic aorta with a pre-stretch.
Mathematics 11 00922 g005
Figure 6. (a) 2D representation of the boundary conditions for the simulation of the pressurization test. (b) Structured mesh of the finite element simulation with 2500 quadrilateral elements.
Figure 6. (a) 2D representation of the boundary conditions for the simulation of the pressurization test. (b) Structured mesh of the finite element simulation with 2500 quadrilateral elements.
Mathematics 11 00922 g006
Figure 7. Block diagram illustrating how the objective function is computed.
Figure 7. Block diagram illustrating how the objective function is computed.
Mathematics 11 00922 g007
Figure 8. Evolutionary strategy block diagram.
Figure 8. Evolutionary strategy block diagram.
Mathematics 11 00922 g008
Figure 9. Adjustment of a known solution of an anisotropic material with evolutionary strategies.
Figure 9. Adjustment of a known solution of an anisotropic material with evolutionary strategies.
Mathematics 11 00922 g009
Figure 10. Exact solution and the solution obtained with ES elitist.
Figure 10. Exact solution and the solution obtained with ES elitist.
Mathematics 11 00922 g010
Figure 11. Adjustments of the experimental data of unixial traction with the GHO model and the transverse stretch λ 3 corresponding to each of the adjustments. (a) Tensile test; (b) Transversal Stretch λ 3 .
Figure 11. Adjustments of the experimental data of unixial traction with the GHO model and the transverse stretch λ 3 corresponding to each of the adjustments. (a) Tensile test; (b) Transversal Stretch λ 3 .
Mathematics 11 00922 g011
Figure 12. CN group adjustments. On the left are the results of the non-stabilized characterization and on the right the characterization with the stability term. (a) Pressurization; (b) Pressurization; (c) Tensile Test; (d) Tensile Test; (e) Transversal Stretch; (f) Transversal Stretch.
Figure 12. CN group adjustments. On the left are the results of the non-stabilized characterization and on the right the characterization with the stability term. (a) Pressurization; (b) Pressurization; (c) Tensile Test; (d) Tensile Test; (e) Transversal Stretch; (f) Transversal Stretch.
Mathematics 11 00922 g012aMathematics 11 00922 g012b
Table 1. Search parameters of ES and known solution.
Table 1. Search parameters of ES and known solution.
ParametersKnown SolutionLower LimitUpper Limit
μ [ kPa ] 22.589 0.00 1.00 × 10 5
k 1 [ kPa ] 224.217 0.00 1.00 × 10 6
k 2 2.464 0.00 1.00 × 10 1
κ 0.285 0.00 0.3333
θ ° 39.624 0.00 90
Generations150Population1140
Table 2. Results of the adjustment of a known solution using the GHO model.
Table 2. Results of the adjustment of a known solution using the GHO model.
ParametersES ElitistES Non-Elitist
Average95% ICAverage95% IC
μ [kPa] 2.2260 × 10 1 8.3100 × 10 2 2.1930 × 10 1 8.9184 × 10 2
k 1 [kPa] 2.1303 × 10 2 2.3413 × 10 0 1.9990 × 10 2 1.2902 × 10 0
k 2 2.3705 × 10 0 1.9879 × 10 2 2.2650 × 10 0 1.4794 × 10 2
κ 2.7391 × 10 1 2.3596 × 10 3 2.6000 × 10 1 1.8826 × 10 3
θ ° 4.0400 × 10 1 1.7648 × 10 1 4.8742 × 10 1 7.5583 × 10 3
Objective Function 5.06 × 10 5 7.06 × 10 5
Table 3. Parameters of the fittings obtained from the experimental data of García et al. with evolutionary strategies, stabilization criterion and GHO model.
Table 3. Parameters of the fittings obtained from the experimental data of García et al. with evolutionary strategies, stabilization criterion and GHO model.
Parameters μ [kPa] κ k 1 [kPa] k 2 θ °Objective
Function
García et al. 24.655 0 45.055 5.3279 42.19 48.29
Es
Elitist
Av. 5.996 0.2985 307.27 1.258 72.07 1.473
IC 95% 8.20 × 10 2 5.94 × 10 2 1.37 × 10 0 5.80 × 10 3 3.32 × 10 1
Table 4. ES parameters to characterize a GHO model with finite elements.
Table 4. ES parameters to characterize a GHO model with finite elements.
ParametersLower BoundaryUpper Boundary
μ [kPa]10.0 1.0 × 10 3
k 1 [kPa]0.00 1.0 × 10 4
k 2 0.006.00
κ 0.001/3
θ °0.0090.0
Generations100
Population200
Table 5. Standardized Quadratic Errors of the CN Group.
Table 5. Standardized Quadratic Errors of the CN Group.
Standarized Quadratic ErrorsCN
Without StabilizationWith Stabilization
Presurization 1.28 × 10 5 8.28 × 10 5
Tensile test longitudinal 1.00 × 10 4 2.35 × 10 5
Tensile test circumferential 2.60 × 10 4 1.10 × 10 4
Table 6. GHO model parameters obtained with and without stabilization term.
Table 6. GHO model parameters obtained with and without stabilization term.
Parameters μ [kPa] k 1 [kPa] k 2 κ θ °Computing
Time [h]
Without
Stabilization
10.01520.3150.17086 5.1799 × 10 3 45.8312.4
With
Stabilization
18.66824.9050.29120.136746.4412.0
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Canales, C.; García-Herrera, C.; Rivera, E.; Macías, D.; Celentano, D. Anisotropic Hyperelastic Material Characterization: Stability Criterion and Inverse Calibration with Evolutionary Strategies. Mathematics 2023, 11, 922. https://doi.org/10.3390/math11040922

AMA Style

Canales C, García-Herrera C, Rivera E, Macías D, Celentano D. Anisotropic Hyperelastic Material Characterization: Stability Criterion and Inverse Calibration with Evolutionary Strategies. Mathematics. 2023; 11(4):922. https://doi.org/10.3390/math11040922

Chicago/Turabian Style

Canales, Claudio, Claudio García-Herrera, Eugenio Rivera, Demetrio Macías, and Diego Celentano. 2023. "Anisotropic Hyperelastic Material Characterization: Stability Criterion and Inverse Calibration with Evolutionary Strategies" Mathematics 11, no. 4: 922. https://doi.org/10.3390/math11040922

APA Style

Canales, C., García-Herrera, C., Rivera, E., Macías, D., & Celentano, D. (2023). Anisotropic Hyperelastic Material Characterization: Stability Criterion and Inverse Calibration with Evolutionary Strategies. Mathematics, 11(4), 922. https://doi.org/10.3390/math11040922

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop