- Email: [email protected]
Viscoelastic characterization of woven Dacron for aortic grafts by using direction-dependent quasi-linear viscoelasticity
Journal of the Mechanical Behavior of Biomedical Materials 82 (2018) 282–290
Contents lists available at ScienceDirect
Journal of the Mechanical Behavior of Biomedical Materials journal homepage: www.elsevier.com/locate/jmbbm
Viscoelastic characterization of woven Dacron for aortic grafts by using direction-dependent quasi-linear viscoelasticity
T
⁎
Marco Amabili , Prabakaran Balasubramanian, Ivan Breslavsky, Giovanni Ferrari, Eleonora Tubaldi Department of Mechanical Engineering, McGill University, 817 Sherbrooke Street West, Montreal, Canada H3A 0C3
A R T I C LE I N FO
A B S T R A C T
Keywords: Viscoelasticity Quasi-linear Relaxation Dacron Graft Dynamic modulus
In case of direction-dependent viscoelasticity, a simplified formulation of the three-dimensional quasi-linear viscoelasticity has been obtained manipulating the original Fung equation. The experimental characterization of the static hyperelastic behaviour, the relaxation, the dynamic modulus and the loss factor of woven Dacron from a commercial aortic prosthesis has been performed. An 11% difference of the reduced relaxation (after infinite time) between axial and circumferential directions has been observed for the woven Dacron. A very large increase in stiffness is obtained in case of harmonic loading with respect to the static loading. These findings are particularly relevant for dynamic modelling of currently used aortic grafts.
1. Introduction The quasi-linear viscoelasticity (QLV) was introduced in biomechanics by Fung (1973, 1993), even if it was known before for rubber materials. The theory of linear viscoelasticity applies for small deformations. In case of finite deformations, like the ones reached in arteries due to the pressure and flow pulsation, it is necessary to use nonlinear (hyperelastic) stress-strain relationships. In order to address this issue, Fung introduced the QLV. The superposition principle is assumed to hold true in the QLV. The superposition is a characteristic of linear systems; in the QLV it is applied to nonlinear elasticity and named nonlinear superposition. This is the reason for labelling the theory as quasi-linear: it extends the linear viscoelasticity to hyperelastic materials. The QLV presents the simplification that the reduced relaxation function is independent of the strain level. The QLV theory is still very popular, as proved by its wide use in biomechanics literature. Here just a short review is reported. Dortmans et al. (1994) obtained the exact formulation for the reduced creep function in the case of continuous spectrum of relaxation within the framework of the QLV. Puso and Weiss (1998) formulated a finite element implementation of anisotropic QLV using a discrete spectrum approximation and a single relaxation function. An implementation of three-dimensional QLV has been developed by Bischoff (2006) and Giles et al. (2007), by using a single relaxation function for stresses in different directions. Giles et al. (2007) also obtained
⁎
interesting results by using a different nonlinear viscoelastic model based on the study conducted by Holzapfel et al. (2002). Sarver et al. (2003) derived stress normalization methods for QLV modelling of soft tissue. Gimbel et al. (2004) studied the effect of overshoot on estimated QLV parameters. Abramovitch and Woo (2004) developed an improved technique for fitting experimental relaxation experiments to the QLV relaxation function for continuous spectrum. Troyer et al. (2012) introduced a correction method for stress relaxation experiments, applicable also to QLV. Even in recent experiments on relaxation, data are fitted by using the reduced relaxation function obtained by Fung with the QLV (Castile et al., 2016). Babaei et al. (2015) and Babaei et al. (2017) introduced a discrete spectral analysis for determining the QLV properties of biological materials. In the present study, a simple formulation of the three-dimensional QLV is obtained manipulating the original Fung (1993) equation. It represents a relevant simplification in case of direction-dependent viscoelasticity. Experiments on woven Dacron used for commercial aortic prostheses (Hemashield Platinum by Maquet) have been carried out to characterize the static behaviour, the relaxation, the dynamic modulus and the loss factor. No previous studies on relaxation of woven Dacron (PET: polyethylene terephthalate), except the one by Lee and Wilson (1986) at low strain rates, are known to the authors, while the static mechanical properties have been well investigated (e.g. Hasegawa and Azuma, 1979; Lee and Wilson, 1986; How, 1992; Yeoman et al., 2010; Bustos et al., 2016). To the authors’ knowledge, no data are available in
Correspondence to: Department of Mechanical Engineering, McGill University, Macdonald Engineering Building, 817 Sherbrooke Street West, Montreal, PQ, Canada H3A 0C3. E-mail address: [email protected] (M. Amabili).
https://doi.org/10.1016/j.jmbbm.2018.03.038 Received 13 February 2018; Received in revised form 9 March 2018; Accepted 29 March 2018 Available online 30 March 2018 1751-6161/ © 2018 Elsevier Ltd. All rights reserved.
Journal of the Mechanical Behavior of Biomedical Materials 82 (2018) 282–290
M. Amabili et al.
the literature for the loss factor and dynamic modulus of woven Dacron.
σ͠ kl =
2. Materials and methods
The three-dimensional QLV was introduced by Fung in Section 7.13 of his monograph (1993) by using the equation (here a different notation is used) 3
3
∑ ∑ ⎧⎨σ͠ kl [E (0 + )] Gijkl (t ) + ∫0 k=1 l=1
t
Gijkl (t − τ )
⎩
∂σ͠ kl [E (τ )] ⎫ dτ , ⎬ ∂τ ⎭
−1
⎧ ⎡ ⎛ t ⎞ ⎛ t ⎞⎤⎫⎡ ⎛ τ2ij ⎞ ⎤ Gij (t ) = 1 + aij ⎢g ⎜ ⎟ − g ⎜ ⎟ ⎥ ⎢1 + aij ln ⎜ ⎟ ⎥ ⎨ τ τ2ij τ1ij ⎬ ⎝ ⎠⎦⎭⎣ ⎝ 1ij ⎠ ⎦ ⎣ ⎝ ⎠ ⎩
(1) where σij are the components of the 2nd Piola-Kirchhoff stress tensor (named simply Kirchhoff by Fung) σ(t ) , E is the Green's strain tensor, σ͠ kl are the components of the instantaneous 2nd Piola-Kirchhoff stress tensor σ͠ corresponding to the static strain tensor E, and Gijkl are the components of the 4th-order tensor of the reduced relaxation functions. The symbol 0+ indicates the limit for t going to zero from positive values of t. Many of the components of Gijkl are not independent. If the derivatives of both Gijkl and σ͠ kl are continuous functions, Eq. (1) can be integrated by parts 3
σij (t ) =
3
∑ ∑ ⎧Gijkl (t = 0) σ͠ kl [E (t )] + ∫0
t
∑ ∑ ⎧Gijkl (t = 0) σ͠ kl [E (t )] + ∫0
t
k=1 l=1 3 3
=
k=1 l=1
⎨ ⎩
⎨ ⎩
σ͠ kl [E (t − τ )] σ͠ kl [E (τ )]
∂Gijkl (τ ) ∂τ
∂Gijkl (t − τ ) ∂ (t − τ )
g (t ) =
σij (t ) = σ͠ ij [E (t )] +
∫0
t
Gij (t − τ )
σ͠ ij [E (τ )]
∂σ͠ ij [E (τ )]
∂Gij (t − τ ) ∂ (t − τ )
∂τ
dτ .
dτ ,
dτ ⎫ ⎬ ⎭
∂W [E (t )] , ∂E (t )
e−x dx , for t ≥ 0. x
(9)
. (10)
The reduced relaxation function G(t) given in Eq. (8) is plotted in Fig. 1 versus time for τ1= 0.01 s, τ2 = 100 s and three different values of the coefficient a: 0.05, 0.1, 0.25.
dτ ⎫ , ⎬ ⎭
2.2. Experimental procedure 2.2.1. Material Strips of dimension about 45 × 10.5 mm, taken in axial and circumferential directions, were cut from an aortic woven Dacron graft (Maquet 175428 P Hemashield Platinum Woven Double Velour, nominal diameter 32 mm), as shown in Fig. 2(a-c). The width of the axial strip used in the present study is 10.59 mm, with a standard deviation of 0.134 mm, while the width of the circumferential strip is 10.50 mm, with a standard deviation of 0.178 mm. The thickness of the strips was measured in ten different points by using a Mitutoyo digital micrometer of diameter a quarter of an inch under a 0.1 MPa tensile stress. The resulting thickness is 0.328 mm with a standard deviation of 0.023 mm. The thickness, if measured with a smaller tip, is variable with the position being the material a fabric. The cross-section of the test specimen is assumed to be rectangular of constant thickness and the material is conveniently modelled as a continuum, since it is practical in numerical codes. However, the space between the fibers is actually occupied by air. This is not taken into account in the present material model, which is derived at macro-scale.
(4)
(5)
Eqs. (4) and (5), with respect to Eq. (2), do not contain summations anymore, so each stress has its own reduced relaxation function, which is much simpler and more intuitive. In addition, the 4th-order tensor of the reduced relaxation functions has been substituted with a 2nd-order tensor G(t) of components Gij(t). This simplification is a consequence of considering the instantaneous elastic stress σ͠ instead of the strain in the QLV theory. Therefore, the same simplification cannot be achieved for the three-dimensional linear viscoelasticity. The number of independent reduced relaxation functions in the three-dimensional QLV theory is 6, even in the most general anisotropic case. Even more important, there is no cross-link between relaxation in different directions. The reduced relaxation functions Gij (t ) can be identified by three uniaxial and three shear relaxation tests. The instantaneous elastic stresses σ͠ kl [E (t )] due to the Green's strain tensor E(t ) can be obtained from the hyperelastic model of the material
σ͠ =
∞
−1
where δ indicates the Kronecker delta. Eq. (3) shows that many of the functions Gijkl(t) are zero at t = 0; Gijkl(t) decrease with time and are positive-valued functions, so they must be zero at any time if their initial value is zero. Therefore, the relationship Gijkl(t) = 0 holds true if i ≠ k and j ≠ l . This property can be used to rewrite Eqs. (1) and (2) in a simplified way t
∫t
⎡ ⎛ τ2ij ⎞ ⎤ Gij (∞) = ⎢1 + aij ln ⎜ ⎟ ⎥ τ ⎝ 1ij ⎠ ⎦ ⎣
(3)
∫0
(8)
For t → ∞, the function g(t) tends to zero and Gij(t) takes the limit value
where the last passage has been obtained by using the variable change ξ = t − τ and then changing back the variable symbol ξ to τ in the resulting expression. Eqs. (1) and (2) are equivalent. Eq. (2) must be valid also for t = 0, for which σij (t = 0) = σ͠ ij [E (0)]; this gives
σij (t ) = σ͠ ij [E (0 + )] Gij (t ) +
,
where g(t) is the exponential integral function defined by
(2)
Gijkl (t = 0) = δik δjl,
(7)
In Eqs. (6) and (7) W is the strain energy density function of the material obtained from one of the many hyperelastic models available in the literature for anisotropic materials. The following expression, obtained by Fung (1993) using the continuous spectrum of relaxation, can be applied for each one of the six independent functions Gij(t), so that they differ from each other only for the three material coefficients aij, τ1ij and τ2ij
2.1. Three-dimensional QLV theory
σij (t ) =
∂W [ε11 (t ), ε12 (t ), ε13 (t ), ...] . ∂εkl (t )
1.0
Gt
0.8
a = 0.05
0.6
a = 0.1
0.4
a = 0.25
0.2 0.0 0.001
0.010
0.100
1
10
100
1000
Time s Fig. 1. Reduced relaxation function G(t) versus time for τ1= 0.01 s, τ2 = 100 s; curves for three different values of a are plotted: continuous line, a = 0.05; dashed line, a = 0.1; dot-dashed line, a = 0.25.
(6)
which can be rewritten in terms of components 283
Journal of the Mechanical Behavior of Biomedical Materials 82 (2018) 282–290
M. Amabili et al.
Fig. 2. Photos of the traction experiments. (a) Woven Dacron aortic prosthesis of 32 mm diameter with visible axial corrugation; (b) particular of the Dacron fabric; (c) axial and circumferential strips and one cut prosthesis; (d) a circumferential strip during traction test.
Test System with MTEST Quattro controller and Interface WMC10 (50 N) load cell at room temperature, as shown in Fig. 2(d). Both tensile and relaxation tests were carried out on axial and circumferential strips after 10 preconditioning cycles at the same strain rate and reaching the same peak stress of the actual tests. Several strips were tested from the same graft. A similar outcome was obtained from different samples; results are reported for one axial (initial length between the jaws 30.17 mm, with standard deviation of 0.05 mm) and one circumferential strip (initial length between the jaws 30.31 mm, with standard deviation of 0.05 mm). During the tensile tests, a strain rate of 0.025 mm/s was used. The distance between the jaws and the tensile forces were recorded directly by the test system with accuracy of ± 1 µm and ± 0.01 N, respectively. The relaxation test was performed with a fast-loading linear tensile ramp, obtained for constant strain rate about 0.825 mm/s corresponding to the maximum rate possible for the test system, followed by a hold at constant strain for 10,000 s.
Piezo load cell Shaker
Laser spot Strip holder Fixed surface
Fig. 3. Experimental setup for the dynamic stiffness measurement of a circumferential strip.
The graft and the corresponding strips are corrugated in axial direction. Since these corrugations tend to permanently decrease when the prosthesis is implanted or loaded for a long enough time, the graft was left submerged in water with a weight of 420 g for 4 weeks before any testing. The peak-to-peak average amplitude of the residual corrugations is 0.52 mm and the wavelength is 2.32 mm in the axial direction, obtained from image processing. Since the material (and not the global properties of the corrugated tube) is of interest in the present study, the axial strip was modestly preloaded during all tests to remove the residual corrugation.
2.2.3. Dynamic stiffness test The experimental setup for the dynamic stiffness measurement is shown in Fig. 3. A circumferential/axial strip was glued by cyanoacrylate to two L-shaped strip holders made of aluminum to minimize inertia. It was verified that they replace the jaws with great accuracy for the applied loads and the glue did not introduce any significant change in the stiffness. A strip holder was connected to a fixed surface through a B&K 8203 piezoelectric load cell. The other strip holder was fixed to the armature of a B&K 4824 electrodynamic exciter (shaker), connected to a B&K 2732 power amplifier and a LMS/Siemens Test-Lab system, and its displacement was measured by a Polytec OFV-505 laser Doppler vibrometer with displacement decoder. The initial free length between
2.2.2. Uniaxial tensile and relaxation tests Uniaxial tensile tests were performed on an eXpert 4000 Admet Micro 284
Journal of the Mechanical Behavior of Biomedical Materials 82 (2018) 282–290
M. Amabili et al.
the strip holders was 30.05 mm, with a standard deviation of 0.15 mm, for the axial strip and 30.22 mm, with a standard deviation of 0.24 mm, for the circumferential strip. Initially, a preconditioning of the strip was performed by increasing and decreasing the DC voltage applied to the shaker for 10 cycles. After that, a constant DC voltage was supplied to ensure the desired stretch. Both strips were tested for three different stretch levels. After allowing 300 s for relaxation, a small AC voltage was added to generate a harmonic excitation of amplitude about 0.1 N, giving a dynamic elongation of the order of 10 µm, at various frequencies. The frequency range of harmonic excitation was from 1 to 60 Hz with 0.1–0.5 Hz adaptive resolution. At every frequency, the initial 50 cycles were discarded to eliminate the transient vibration. The dynamic elongation of the strip Δl and the harmonic force Δf applied to the strip were measured for all the excitation frequencies for 10 cycles. From the measured dynamic force, elongation and the phase between them, the storage modulus and loss factor of the strips are calculated. The amplitude of the complex modulus E* is given by
Δf / A , Δl / L
E* =
(11)
where, L is the free length of the strip and A its cross-section area. The storage modulus E ′ and the loss modulus E ″ are given by
E ′ = E * cos ϕ,
E ″ = E * sin ϕ,
(12a,b)
where ϕ is the measured phase angle between the harmonic force and the dynamic elongation. The loss factor η is
η = tan ϕ.
(13) Fig. 4. Uniaxial tensile tests on strips from a woven Dacron aortic graft; 2nd Piola-Kirchhoff stresses, Green strains. The blue dotted curves are those taken to fit the material model. (a) Circumferential strip, 9 repeated tests, mean values and standard deviations; (b) axial strip, 9 repeated tests, mean values and standard deviations. (For the interpretation of references to colour in this figure legend, teh reader is referred to the web version of the article.)
2.3. Constitutive modelling 2.3.1. Hyperelastic law The following anisotropic Fung-type strain energy density function is applied to the woven Dacron 2
W=
∑ Ci (exp(Qi) − 1),
2 Qi = 4μi (εxx cos4 βi + 2εxx εθθ cos2 βi sin2 βi
Eq. (16) is minimized by a genetic algorithm (Mitchell, 1996), which can be briefly described as follows. We consider a population of AT candidate solutions (called “individuals”), each of them represented by a “genetic code”, i.e. a vector of 6 parameters {C1, μ1 , β1, C2, μ 2 , β2} . At the initial stage, a large enough pool of AT individuals with random genetic codes is created. Then, iterations of the algorithm start. Initially, for every individual “fitting”, represented by Eq. (16), is evaluated. Afterwards, all individuals are arranged according to their fitting. The AS ( AS < < AT ) most fitted individuals, i.e. those with the smallest corresponding value of fs obtained by Eq. (16), “survive” and they are included in the list of candidates for the next iteration. Then, AT − AS “offspring” individuals are added to the pool. The offspring individuals are generated in two different ways: (i) from one single parent, randomly varying the parameters of every of the AS surviving individuals, and (ii) from two parents. In the latter case, the offspring genetic code contains parameters with random values comprised between those of the parents. The new pool of individuals contains AS survived and AT − AS offspring individuals, so the total number of candidate solutions in the pool remains AT . With this new pool, a new iteration starts. Iterations continue until the solution with “small enough value” of the error function (16) is found. Generally, the genetic algorithm is a good optimization method of complicated multidimensional problems with complex shape of the error function. The sensitivity of the algorithm to the initial guesses is low due to the very large number AT of candidate solutions in the initial pool.
i=1
+
2 εθθ
sin4 βi ),
(14)
where Ci > 0 are stress-like material parameters, μi are non-dimensional material parameters, and βi are angles with the axial direction. Eq. (14) is a formulation for a material with two families of different fibers. The woven Dacron is made of circumferential and axial fibers that are not completely independent due to weaving. The convexity of W is verified for Ci, μi > 0 (Holzapfel, 2006). The 2nd Piola-Kirchhoff principal stresses take the following expressions
σ͠ xx =
σ͠ θθ =
2
∂W = ∂εxx
∑ Ci exp(Qi)8μi (εxx cos4 βi + εθθ cos2 βi sin2 βi),
∂W = ∂εθθ
∑ Ci exp(Qi)8μi (εxx cos2 βi sin2 βi + εθθ sin4 βi).
i=1
(15a)
2 i=1
(15b)
In case of uniaxial extension in axial direction x, Eq. (15a) is used to fit the experimental data and (15b) is used to obtain εθθ as a function of εxx considering σ͠ θθ = 0 . Each one of Eq. (15a,15b) fits the corresponding experimental data set, in order to obtain the unknown material parameters Ci, μi ,βi . The material parameters that guarantee the best fit of the experimental data are obtained by minimizing the following nonlinear stress based function N
fs =
⎧ ∂W
∑ ⎨ ⎜⎛ ∂ε
n=1
⎩⎝
xx n
2
2
∂W (n) ⎞ − σxx +⎛⎜ ⎟ ∂ ⎠(axial) ⎝ εθθ
n
⎫ (n) ⎞ − σθθ , ⎟ ⎬ ⎠(circ ) ⎭
(16)
2.3.2. QLV for relaxation In order to fit relaxation experiments with the reduced relaxation function given in Eq. (8), the following procedure is applied
where n indicates the n-th experimental data and the subscripts “axial” and “circ” indicate which strip is considered (axial or circumferential). 285
Journal of the Mechanical Behavior of Biomedical Materials 82 (2018) 282–290
M. Amabili et al.
Fig. 5. Fitting of the hyperelastic material model to the experimental tensile test data. •, experimental data; —, material model; 2nd Piola-Kirchhoff stresses, Green strains. (a) Circumferential strip; (b) axial strip.
Table 1 Material parameters for the constitutive law of the woven Dacron obtained from tensile test. C1 (kPa)
μ1
β1 (rad)
C2 (kPa)
μ2
β2 (rad)
2.455
1216.51
− 1.233
25.891
49.28
0
Fig. 7. Fitting of the QLV model to relaxation experiments. •, experimental data; — , constitutive material model; 2nd Piola-Kirchhoff stresses versus time. (a) Circumferential strip, t0θθ = 0.55 s, γθθ = 0.030 s−1; (b) axial strip, t0xx = 1.535 s, γxx = 0.030 s−1. Table 2 Reduced relaxation parameters of the woven Dacron according to the QLV theory. Fig. 6. Relaxation tests on strips from a woven Dacron aortic graft; 2nd PiolaKirchhoff stresses versus time. (a) Circumferential strip, 8 repeated tests, mean values and standard deviations; (b) axial strip, 8 repeated tests, mean values and standard deviations.
τ1 (s)
τ2 (s)
a
Circumferential Axial
0.049 0.049
1545.4 2670
0.038 0.023
where γxx is the constant strain rate. The time derivative of (15a) in the interval 0 < t < t0xx is
(Abramowitch and Woo, 2004). The instantaneous hyperelastic stresses in axial and circumferential directions are given by Eq. (15a,15b). In the time interval 0 < t < t0xx , i.e. during the ramp for the axial strip, the strain versus time is given by
εxx (τ ) = γxx τ ,
Direction
(17) 286
Journal of the Mechanical Behavior of Biomedical Materials 82 (2018) 282–290
M. Amabili et al.
λ = 1.012 λ = 1.010 λ = 1.008
(a)
λ = 1.042 λ = 1.039 λ = 1.037
Fig. 8. Storage modulus versus loading frequency at three different initial prestretches λ; the error bar at 30 Hz quantify the estimated measurement error. (a) Circumferential strip; (b) axial strip.
(b) Fig. 10. Storage modulus versus loading frequency at three different initial prestretches λ; ○, experimental data; —, generalized Maxwell material model with N = 6. (a) Circumferential strip; (b) axial strip.
λ = 1.008 λ = 1.010 λ = 1.012
(a)
λ = 1.037 λ = 1.039 λ = 1.042
Fig. 9. Loss factor versus loading frequency at three different initial pre-stretches λ; the error bar at 30 Hz quantify the estimated measurement error. (a) Circumferential strip; (b) axial strip.
(b) Fig. 11. Loss factor versus loading frequency at three different initial prestretches λ; ○, experimental data; —, generalized Maxwell material model with N = 6. (a) Circumferential strip; (b) axial strip.
287
Journal of the Mechanical Behavior of Biomedical Materials 82 (2018) 282–290
M. Amabili et al.
Table 3 Material parameters for the generalized Maxwell model with N = 6. Axial strip
λ = 1.008
Circumferential strip
λ = 1.037
λ = 1.039
λ = 1.042
λ = 1.008
λ = 1.01
β1∞
0.273
0.543
0.56
1.619
1.582
1.673
β2∞
0.183
0.19
0.17
0.9
0.772
0.64
λ = 1.010
λ = 1.012
λ = 1.012
β3∞
0.207
0.2
0.184
0.957
0.929
0.768
β4∞
0.011
0.004
0.006
0.033
0.04
0.053
β5∞
0.028
0.122
0.078
0.313
0.437
0.613
β6∞
0.202
0.321
0.474
1.422
1.477
2.162
τ1 τ2 τ3 τ4 τ5 τ6
0.001 0.016 0.136 7.058 59.803 1478.05
0.0006 0.015 0.122 6.923 91.364 1115.07
0.0006 0.014 0.113 7.07 90.309 1632.99
0.002 0.016 0.124 7.896 79.65 849.29
0.002 0.014 0.125 7.287 87.364 1418.6
0.001 0.012 0.105 0.334 104.68 1094.53
(a)
λ = 1.042
λ = 1.039 λ = 1.012
λ = 1.042
λ = 1.010 λ = 1.008
(b) Fig. 13. Loss factor versus loading frequency at three different initial prestretches λ; ○, experimental data; —, generalized Maxwell material model with N = 6 and the same parameters βα and τα for all the pre-stretches. (a) Circumferential strip; (b) axial strip.
(a)
Table 4 Material parameters for the generalized Maxwell model with N = 6 and the same parameters βα and τα for all the pre-stretches.
λ = 1.042 λ = 1.039 λ = 1.037
(b) Fig. 12. Storage modulus versus loading frequency at three different initial prestretches λ; ○, experimental data; —, generalized Maxwell material model with N = 6 and the same parameters βα and τα for all the pre-stretches. (a) Circumferential strip; (b) axial strip.
∂σ͠ xx = ∂τ
Axial strip
Circumferential strip
β1∞ β2∞ β3∞ β4∞ β5∞ β6∞
0.352
1.507
0.172
0.816
0.188
0.817
0.01
0.114
0.07
0.921
0.319
1.108
τ1 τ2 τ3 τ4 τ5 τ6
0.001 0.015 0.119 9.102 66.113 1100.32
0.002 0.019 0.136 5.214 67.991 1171.948
2
∑ 64μi2 Ci γxx exp(Qi) ⎡⎢ (γxx τ cos4 βi + εθθ (γxx τ )cos2 βi sin2 βi)2 i=1
σxx (t ) =
⎣
+ (γxx τ cos4 βi + εθθ (γxx τ )cos2 βi sin2 βi)(γxx τ cos2 βi sin2 βi
1 ′ (γxx τ )cos2 βi sin2 βi ) ⎤, (cos4 βi + 2εθθ ⎥ 8μi ⎦
∂σ͠ xx / ∂τ
t
1 + axx ln
( ) τ 2xx τ 1xx
⎧ ⎛t − τ ⎞ ⎛ t − τ ⎞⎤⎫ 1 + axx ⎡ ⎢g ⎜ τ2 ⎟ − g ⎜ τ1 ⎟ ⎥ ⎬ dτ , ⎨ xx ⎝ ⎠ ⎝ xx ⎠ ⎦ ⎭ ⎣ ⎩
0 < t < t0xx .
′ (γxx τ ) + εθθ (γxx τ )sin4 βi ) εθθ +
∫0
(19)
In the time interval t > t0xx , i.e. at constant strain εxx (τ ) = γxx t0xx , the stress is given by (18)
σxx (t ) =
where the prime represents the derivative with respect to the argument. Eq. (4) can be rewritten as
∫0
t 0xx
∂σ͠ xx / ∂τ 1 + axx ln
t > t0xx ,
( ) τ 2xx τ 1xx
⎧ ⎛t − τ ⎞ ⎛ t − τ ⎞⎤⎫ 1 + axx ⎡ ⎢g ⎜ τ2 ⎟ − g ⎜ τ1 ⎟ ⎥ ⎬ dτ , ⎨ xx ⎠ ⎝ xx ⎠ ⎦ ⎭ ⎣ ⎝ ⎩ (20)
since after t = t0xx the derivative of Eq. (15a) with respect to time is 288
Journal of the Mechanical Behavior of Biomedical Materials 82 (2018) 282–290
M. Amabili et al.
3.2. Uniaxial relaxation test
zero. Experiments give the values of the measured stresses Ri at discrete times ti.. Eqs. (19) or (20), where the choice between the two expressions is made according to the value of ti, is evaluated at those specific times, giving σi = σ (t ) t = ti . Then, the following least square objective function fxx, function of the reduced relaxation parameters, is built
The results of eight repeated relaxation tests are presented in Fig. 6(a,b) for the first few seconds (to highlight the ramp) with their mean values and standard deviations. The linear ramp is at constant strain rate γxx = γθθ = 0.030 s−1, and results are almost within the standard deviation of those obtained in Section 3.1 (mostly of the ramp curve is within the standard deviation range of the tensile test data). The different strain rate does not change the results significantly, so the material parameters in Table 1 are used to fit the relaxation experiments. The experimental results and the stress obtained by the constitutive material model versus time are compared in Fig. 7(a,b). Differences between the fitted material model and experiments in the time history are very small, except at the very beginning of the relaxation phase since the test machine cannot stop the ramp instantaneously, but it has a short deceleration phase. The peak value of stress in Fig. 7(a,b) corresponds to the maximum stress compatible with normal physiological conditions in thoracic aortic replacements. The 2nd Piola-Kirchhoff hoop stress σθ in a thin circular cylindrical shell under transmural pressure p and small cylindrical deformation is given by
N
fxx (axx , τ1xx , τ2xx ) =
∑ [Ri − σi]2 ,
(21)
i=1
where N is the number of experimental stresses measured. The three reduced relaxation parameters axx , τ1xx , τ2xx have been identified by minimizing Eq. (21) by using a genetic algorithm. A similar objective function is built for the circumferential stress.
2.3.3. Viscoelasticity for harmonic oscillation In this section, a continuum viscoelastic formulation with internal variables is applied (Holzapfel et al., 2002). The formulation is valid for finite strains and small perturbations away from the static equilibrium configuration. Therefore, it has hypotheses similar to the quasi-linear viscoelasticity previously introduced. It is here restricted to incompressible anisotropic materials. A generalized Maxwell material model is introduced in each of the two orthogonal in-plane directions with N being the number of parallel spring/dashpots arms. It is possible to derive the storage (dynamic) modulus E ′ of the material model as
⎡ E ′ (ω, λ ) = E∞ (λ ) ⎢1 + ⎣
N
(ω τ )2
∑ βα∞ 1 + (ωα τ
α)
α=1
⎤
2⎥
⎦
σθ =
N
∑α = 1 βα∞ N
ω τα 1 + (ω τα )2
1 + ∑α = 1 βα∞
(ω τα )2 1 + (ω τα )2
(24)
where the vessel radius R = 14 mm and thickness h = 0.327 mm are assumed. The stress value σθ = 0.685 MPa is estimated considering the peak value of the physiological pressure p = 120 mmHg (converted in Pa); this stress is close to the peak value in Fig. 7(a,b). The coefficients of the two reduced relaxation functions are given in Table 2. Inserting these values in Eq. (10), the reduced relaxations for t → ∞ are obtained; they are 0.718 and 0.797 in circumferential and axial direction, respectively.
, (22)
where E∞ is the static modulus, ω is the oscillation frequency, βα∞ > 0 are nondimensional stiffness coefficients and τα are time constants. Both E ′ and E∞ are functions of the initial static stretch λ of the material; in general, the same is true for βα∞ and τα . However, their variation with λ may be modest to the point that constant values can be assumed in some applications. The integer N must be chosen large enough to fit well the experimental data. The loss factor, which coincides with the loss tangent in case of linear viscoelasticity, is given by
η (ω) =
pR , h
3.3. Storage modulus and loss factor The storage modulus in circumferential and axial directions is plotted in Fig. 8(a,b) versus the loading frequency. Results for zero frequency, taken from the data in Section 3.1, are reported for comparison. The estimated measurement error at 30 Hz is indicated. It is evident that a very large increase in stiffness is obtained in case of harmonic loading with respect to the static loading. The stiffness moderately increases with the loading frequency between 1 and 60 Hz, which is the range of interest in physiological applications (the heart rate can vary from 0.66 Hz to 3.33 Hz, but the pressure and flow pulse waves introduce higher harmonics with frequency several times higher than the heart rate). Three curves are presented for each direction at different initial stretches since the material is hyperelastic; the modulus increases with the stretch, as expected. The corresponding loss factors are reported in Fig. 9(a,b); they are relatively flat between 1 and 60 Hz. The loss factor of the circumferential strip is larger (about double) than the axial strip. No significant difference in the loss factor is observed at different stretch for the axial strip, while a decrease with increased stretch is observed in circumferential direction. The storage modulus and loss factor are necessary for dynamic models of aortic prostheses (Tubaldi et al., 2018; Jayendiran et al., 2018). The experimental results and the storage moduli obtained by the generalized Maxwell material model with N = 6 versus frequency are compared in Fig. 10(a,b), while the corresponding loss factors are presented in Fig. 11(a,b). The material model fits with good accuracy the experiments for both the storage modulus and loss factors for the three pre-stretch levels investigated. In particular, the circumferential strip is better approximated than the axial strip. The corresponding material parameters βα and τα are given in Table 3. For each direction, the parameters for the three pre-stretch levels are somehow different; this is due to the simultaneous optimization of βα and τα for the best fitting. However, it is possible to use the same values of βα and τα to fit
. (23)
A minimization procedure based on a genetic algorithm is used also in this case to identify the material parameters βα and τα , while E∞ (λ ) is computed from Eq. (15a,15b). Eqs. (22) and (23) are then specified for x and θ directions. 3. Results and discussion 3.1. Uniaxial tensile test The results of nine repeated uniaxial tensile tests are presented in Fig. 4(a,b) and show very different stiffness values in circumferential and axial directions and relatively small standard deviation, indicating good repeatability after preconditioning. The experimental data from the blue dotted curves in Fig. 4(a,b) and the hyperelastic material model are presented in Fig. 5(a,b) for Dacron in circumferential and axial directions, respectively. The identified hyperelastic material parameters are given in Table 1. A good agreement between the hyperelastic law and the experimental data is observed in Fig. 5(a,b) for the full range of values presented, which covers the full range of physiological σ-ε values for an aortic implant. Results agree to those by Lee and Wilson (1986) for the circumferential direction, while in axial direction the corrugation has been removed in the present study by applying a small initial pre-stress, so the comparison is not possible. 289
Journal of the Mechanical Behavior of Biomedical Materials 82 (2018) 282–290
M. Amabili et al.
the three stretch levels. This suggests the possibility of building a material model with ratio between dynamic and static modulus and loss factor independent of the initial pre-stretch, as introduced in (Holzapfel et al., 2002), with great simplification. The results of this simplified material model are presented in Fig. 12(a,b) for the storage moduli and in Fig. 13(a,b) for the loss factors. In this case, a single loss factor curve is obtained for all the pre-stretches. The material parameters identified in this case are given in Table 4. Clearly, the material model does not compare to the experimental results as well as the one with parameters βα and τα identified at each pre-stretch level, but it is of much simpler implementation in numerical codes.
biological materials. J. R. Soc. Interface 12, 20150707. Babaei, B., Velasquez-Mao, A.J., Thomopoulos, S., Elson, E.L., Abramowitch, S.D., Genin, G.M., 2017. Discrete quasi-linear viscoelastic damping analysis of connective tissues, and the biomechanics of stretching. J. Mech. Behav. Biomed. Mater. 69, 193–202. Bischoff, J.E., 2006. Reduced parameter formulation for incorporating fiber level viscoelasticity into tissue level biomechanical models. Annu. Biomed. Eng. 34, 1164–1172. Bustos, C.A., García-Herrera, C.M., Celentano, D.J., 2016. Mechanical characterisation of Dacron graft: experiments and numerical simulation. J. Biomech. 49, 13–18. Castile, R.M., Skelly, N.W., Babaei, B., Borphy, R.H., Lake, S.P., 2016. Microstructural properties and mechanics vary between bundles of the human anterior cruciate ligament during stress-relaxation. J. Biomech. 49, 87–93. Dortmans, L.J.M.G., van de Ven, A.A.F., Sauren, A.A.H.J., 1994. A note on the reduced creep function corresponding to the quasi-linear visco-elastic model proposed by Fung. J. Biomech. Eng. 116, 373–375. Fung, Y.C., 1973. Biorheology of soft tissues. Biorheology 10, 139–155. Fung, Y.C., 1993. Biomechanics: Mechanical Properties of Living Tissues, 2nd ed. Springer, New York, USA. Giles, J.M., Black, A.E., Bischoff, J.E., 2007. Anomalous rate dependence of the preconditioned response of soft tissue during load controlled deformation. J. Biomech. 40, 777–785. Gimbel, J.A., Sarver, J.J., Soslowsky, L.J., 2004. The effect of overshooting the target strain on estimating viscoelastic properties from stress relaxation experiments. J. Biomech. Eng. 126, 844–848. Hasegawa, M., Azuma, T., 1979. Mechanical properties of synthetic arterial grafts. J. Biomech. 12, 509–517. Holzapfel, G.A., 2006. Determination of material models for arterial walls from uniaxial extension tests and histological structure. J. Theor. Biol. 238, 290–302. Holzapfel, G.A., Gasser, T.C., Stadler, M., 2002. A structural model for the viscoelastic behavior of arterial walls: continuum formulation and finite element analysis. Eur. J. Mech. A/Solids 21, 441–463. How, T.V., 1992. Mechanical properties of arteries and arterial grafts. In: Hastings, G.W. (Ed.), Cardiovascular Materials. Springer, Berlin, Germany, pp. 1–35. Jayendiran, R., Nour, B.M., Ruimi, A., 2018. Dacron graft as replacement to dissected aorta: a three-dimensional fluid-structure-interaction analysis. J. Mech. Behav. Biomed. Mater. 78, 329–341. Lee, J.M., Wilson, G.J., 1986. Anisotropic tensile viscoelastic properties of vascular graft materials tested at low strain rates. Biomaterials 7, 423–431. Mitchell, M., 1996. An Introduction to Genetic Algorithms. MIT Press, Cambridge, MA, USA. Puso, M.A., Weiss, J.A., 1998. Finite element implementation of anisotropic quasi-linear viscoelasticity using a discrete spectrum approximation. J. Biomech. Eng. 120, 62–70. Sarver, J.J., Robinson, P.S., Elliott, D.M., 2003. Methods for quasi-linear viscoelastic modeling of soft tissue: application to incremental stress–relaxation experiments. J. Biomech. Eng. 125, 754–758. Troyer, K.L., Estep, D.J., Puttlitz, C.M., 2012. Viscoelastic effects during loading play an integral role in soft tissue mechanics. Acta Biomater. 8, 234–243. Tubaldi, E., Païdoussis, M.P., Amabili, M., 2018. Nonlinear dynamics of Dacron aortic prostheses conveying pulsatile flow. J. Biomech. Eng. 240 art. 061004 (12 pages). Yeoman, M.S., Reddy, B.D., Bowles, H.C., Bezuidenhout, D., Zilla, P., Franz, T., 2010. A constitutive model for the warp-weft coupled non-linear behavior of knitted biomedical textiles. Biomaterials 31, 8484–8493.
4. Conclusions Eqs. (4) and (5) have probably been obtained in the present study for the first time; they allow a relevant simplification in case of direction-dependent viscoelasticity. A 11% difference of the reduced relaxations for t → ∞ between axial and circumferential directions has been observed for the woven Dacron; differences in the loss factors are much larger. This indicates that investigation of direction-dependent viscoelasticity is worth the effort, and the present simplified formulation of the three-dimensional QLV facilitates the applications. The experimental data for storage moduli and loss factors versus loading frequency and the corresponding material parameters obtained by the generalized Maxwell model are particularly relevant for dynamic models of aortic grafts. The dynamic modulus is much larger than the static modulus in both axial and circumferential directions, which is very significant in the evaluation of the radial deformation of aortic grafts under pulsatile blood pressure. Acknowledgments The authors acknowledge the financial support of the NSERC Discovery Grant (372493-13), the Canada Research Chair program and the Qatar grant NPRP 7–032-2–016. References Abramovitch, S.D., Woo, S.L.-Y., 2004. An improved method to analyze the stress relaxation of ligaments following a finite ramp time based on the quasi-linear viscoelastic theory. J. Biomech. Eng. 126, 92–97. Babaei, B., Abramowitch, S.D., Elson, E.L., Thomopoulos, S., Genin, G.M., 2015. A discrete spectral analysis for determining quasi-linear viscoelastic properties of
290
Contents lists available at ScienceDirect
Journal of the Mechanical Behavior of Biomedical Materials journal homepage: www.elsevier.com/locate/jmbbm
Viscoelastic characterization of woven Dacron for aortic grafts by using direction-dependent quasi-linear viscoelasticity
T
⁎
Marco Amabili , Prabakaran Balasubramanian, Ivan Breslavsky, Giovanni Ferrari, Eleonora Tubaldi Department of Mechanical Engineering, McGill University, 817 Sherbrooke Street West, Montreal, Canada H3A 0C3
A R T I C LE I N FO
A B S T R A C T
Keywords: Viscoelasticity Quasi-linear Relaxation Dacron Graft Dynamic modulus
In case of direction-dependent viscoelasticity, a simplified formulation of the three-dimensional quasi-linear viscoelasticity has been obtained manipulating the original Fung equation. The experimental characterization of the static hyperelastic behaviour, the relaxation, the dynamic modulus and the loss factor of woven Dacron from a commercial aortic prosthesis has been performed. An 11% difference of the reduced relaxation (after infinite time) between axial and circumferential directions has been observed for the woven Dacron. A very large increase in stiffness is obtained in case of harmonic loading with respect to the static loading. These findings are particularly relevant for dynamic modelling of currently used aortic grafts.
1. Introduction The quasi-linear viscoelasticity (QLV) was introduced in biomechanics by Fung (1973, 1993), even if it was known before for rubber materials. The theory of linear viscoelasticity applies for small deformations. In case of finite deformations, like the ones reached in arteries due to the pressure and flow pulsation, it is necessary to use nonlinear (hyperelastic) stress-strain relationships. In order to address this issue, Fung introduced the QLV. The superposition principle is assumed to hold true in the QLV. The superposition is a characteristic of linear systems; in the QLV it is applied to nonlinear elasticity and named nonlinear superposition. This is the reason for labelling the theory as quasi-linear: it extends the linear viscoelasticity to hyperelastic materials. The QLV presents the simplification that the reduced relaxation function is independent of the strain level. The QLV theory is still very popular, as proved by its wide use in biomechanics literature. Here just a short review is reported. Dortmans et al. (1994) obtained the exact formulation for the reduced creep function in the case of continuous spectrum of relaxation within the framework of the QLV. Puso and Weiss (1998) formulated a finite element implementation of anisotropic QLV using a discrete spectrum approximation and a single relaxation function. An implementation of three-dimensional QLV has been developed by Bischoff (2006) and Giles et al. (2007), by using a single relaxation function for stresses in different directions. Giles et al. (2007) also obtained
⁎
interesting results by using a different nonlinear viscoelastic model based on the study conducted by Holzapfel et al. (2002). Sarver et al. (2003) derived stress normalization methods for QLV modelling of soft tissue. Gimbel et al. (2004) studied the effect of overshoot on estimated QLV parameters. Abramovitch and Woo (2004) developed an improved technique for fitting experimental relaxation experiments to the QLV relaxation function for continuous spectrum. Troyer et al. (2012) introduced a correction method for stress relaxation experiments, applicable also to QLV. Even in recent experiments on relaxation, data are fitted by using the reduced relaxation function obtained by Fung with the QLV (Castile et al., 2016). Babaei et al. (2015) and Babaei et al. (2017) introduced a discrete spectral analysis for determining the QLV properties of biological materials. In the present study, a simple formulation of the three-dimensional QLV is obtained manipulating the original Fung (1993) equation. It represents a relevant simplification in case of direction-dependent viscoelasticity. Experiments on woven Dacron used for commercial aortic prostheses (Hemashield Platinum by Maquet) have been carried out to characterize the static behaviour, the relaxation, the dynamic modulus and the loss factor. No previous studies on relaxation of woven Dacron (PET: polyethylene terephthalate), except the one by Lee and Wilson (1986) at low strain rates, are known to the authors, while the static mechanical properties have been well investigated (e.g. Hasegawa and Azuma, 1979; Lee and Wilson, 1986; How, 1992; Yeoman et al., 2010; Bustos et al., 2016). To the authors’ knowledge, no data are available in
Correspondence to: Department of Mechanical Engineering, McGill University, Macdonald Engineering Building, 817 Sherbrooke Street West, Montreal, PQ, Canada H3A 0C3. E-mail address: [email protected] (M. Amabili).
https://doi.org/10.1016/j.jmbbm.2018.03.038 Received 13 February 2018; Received in revised form 9 March 2018; Accepted 29 March 2018 Available online 30 March 2018 1751-6161/ © 2018 Elsevier Ltd. All rights reserved.
Journal of the Mechanical Behavior of Biomedical Materials 82 (2018) 282–290
M. Amabili et al.
the literature for the loss factor and dynamic modulus of woven Dacron.
σ͠ kl =
2. Materials and methods
The three-dimensional QLV was introduced by Fung in Section 7.13 of his monograph (1993) by using the equation (here a different notation is used) 3
3
∑ ∑ ⎧⎨σ͠ kl [E (0 + )] Gijkl (t ) + ∫0 k=1 l=1
t
Gijkl (t − τ )
⎩
∂σ͠ kl [E (τ )] ⎫ dτ , ⎬ ∂τ ⎭
−1
⎧ ⎡ ⎛ t ⎞ ⎛ t ⎞⎤⎫⎡ ⎛ τ2ij ⎞ ⎤ Gij (t ) = 1 + aij ⎢g ⎜ ⎟ − g ⎜ ⎟ ⎥ ⎢1 + aij ln ⎜ ⎟ ⎥ ⎨ τ τ2ij τ1ij ⎬ ⎝ ⎠⎦⎭⎣ ⎝ 1ij ⎠ ⎦ ⎣ ⎝ ⎠ ⎩
(1) where σij are the components of the 2nd Piola-Kirchhoff stress tensor (named simply Kirchhoff by Fung) σ(t ) , E is the Green's strain tensor, σ͠ kl are the components of the instantaneous 2nd Piola-Kirchhoff stress tensor σ͠ corresponding to the static strain tensor E, and Gijkl are the components of the 4th-order tensor of the reduced relaxation functions. The symbol 0+ indicates the limit for t going to zero from positive values of t. Many of the components of Gijkl are not independent. If the derivatives of both Gijkl and σ͠ kl are continuous functions, Eq. (1) can be integrated by parts 3
σij (t ) =
3
∑ ∑ ⎧Gijkl (t = 0) σ͠ kl [E (t )] + ∫0
t
∑ ∑ ⎧Gijkl (t = 0) σ͠ kl [E (t )] + ∫0
t
k=1 l=1 3 3
=
k=1 l=1
⎨ ⎩
⎨ ⎩
σ͠ kl [E (t − τ )] σ͠ kl [E (τ )]
∂Gijkl (τ ) ∂τ
∂Gijkl (t − τ ) ∂ (t − τ )
g (t ) =
σij (t ) = σ͠ ij [E (t )] +
∫0
t
Gij (t − τ )
σ͠ ij [E (τ )]
∂σ͠ ij [E (τ )]
∂Gij (t − τ ) ∂ (t − τ )
∂τ
dτ .
dτ ,
dτ ⎫ ⎬ ⎭
∂W [E (t )] , ∂E (t )
e−x dx , for t ≥ 0. x
(9)
. (10)
The reduced relaxation function G(t) given in Eq. (8) is plotted in Fig. 1 versus time for τ1= 0.01 s, τ2 = 100 s and three different values of the coefficient a: 0.05, 0.1, 0.25.
dτ ⎫ , ⎬ ⎭
2.2. Experimental procedure 2.2.1. Material Strips of dimension about 45 × 10.5 mm, taken in axial and circumferential directions, were cut from an aortic woven Dacron graft (Maquet 175428 P Hemashield Platinum Woven Double Velour, nominal diameter 32 mm), as shown in Fig. 2(a-c). The width of the axial strip used in the present study is 10.59 mm, with a standard deviation of 0.134 mm, while the width of the circumferential strip is 10.50 mm, with a standard deviation of 0.178 mm. The thickness of the strips was measured in ten different points by using a Mitutoyo digital micrometer of diameter a quarter of an inch under a 0.1 MPa tensile stress. The resulting thickness is 0.328 mm with a standard deviation of 0.023 mm. The thickness, if measured with a smaller tip, is variable with the position being the material a fabric. The cross-section of the test specimen is assumed to be rectangular of constant thickness and the material is conveniently modelled as a continuum, since it is practical in numerical codes. However, the space between the fibers is actually occupied by air. This is not taken into account in the present material model, which is derived at macro-scale.
(4)
(5)
Eqs. (4) and (5), with respect to Eq. (2), do not contain summations anymore, so each stress has its own reduced relaxation function, which is much simpler and more intuitive. In addition, the 4th-order tensor of the reduced relaxation functions has been substituted with a 2nd-order tensor G(t) of components Gij(t). This simplification is a consequence of considering the instantaneous elastic stress σ͠ instead of the strain in the QLV theory. Therefore, the same simplification cannot be achieved for the three-dimensional linear viscoelasticity. The number of independent reduced relaxation functions in the three-dimensional QLV theory is 6, even in the most general anisotropic case. Even more important, there is no cross-link between relaxation in different directions. The reduced relaxation functions Gij (t ) can be identified by three uniaxial and three shear relaxation tests. The instantaneous elastic stresses σ͠ kl [E (t )] due to the Green's strain tensor E(t ) can be obtained from the hyperelastic model of the material
σ͠ =
∞
−1
where δ indicates the Kronecker delta. Eq. (3) shows that many of the functions Gijkl(t) are zero at t = 0; Gijkl(t) decrease with time and are positive-valued functions, so they must be zero at any time if their initial value is zero. Therefore, the relationship Gijkl(t) = 0 holds true if i ≠ k and j ≠ l . This property can be used to rewrite Eqs. (1) and (2) in a simplified way t
∫t
⎡ ⎛ τ2ij ⎞ ⎤ Gij (∞) = ⎢1 + aij ln ⎜ ⎟ ⎥ τ ⎝ 1ij ⎠ ⎦ ⎣
(3)
∫0
(8)
For t → ∞, the function g(t) tends to zero and Gij(t) takes the limit value
where the last passage has been obtained by using the variable change ξ = t − τ and then changing back the variable symbol ξ to τ in the resulting expression. Eqs. (1) and (2) are equivalent. Eq. (2) must be valid also for t = 0, for which σij (t = 0) = σ͠ ij [E (0)]; this gives
σij (t ) = σ͠ ij [E (0 + )] Gij (t ) +
,
where g(t) is the exponential integral function defined by
(2)
Gijkl (t = 0) = δik δjl,
(7)
In Eqs. (6) and (7) W is the strain energy density function of the material obtained from one of the many hyperelastic models available in the literature for anisotropic materials. The following expression, obtained by Fung (1993) using the continuous spectrum of relaxation, can be applied for each one of the six independent functions Gij(t), so that they differ from each other only for the three material coefficients aij, τ1ij and τ2ij
2.1. Three-dimensional QLV theory
σij (t ) =
∂W [ε11 (t ), ε12 (t ), ε13 (t ), ...] . ∂εkl (t )
1.0
Gt
0.8
a = 0.05
0.6
a = 0.1
0.4
a = 0.25
0.2 0.0 0.001
0.010
0.100
1
10
100
1000
Time s Fig. 1. Reduced relaxation function G(t) versus time for τ1= 0.01 s, τ2 = 100 s; curves for three different values of a are plotted: continuous line, a = 0.05; dashed line, a = 0.1; dot-dashed line, a = 0.25.
(6)
which can be rewritten in terms of components 283
Journal of the Mechanical Behavior of Biomedical Materials 82 (2018) 282–290
M. Amabili et al.
Fig. 2. Photos of the traction experiments. (a) Woven Dacron aortic prosthesis of 32 mm diameter with visible axial corrugation; (b) particular of the Dacron fabric; (c) axial and circumferential strips and one cut prosthesis; (d) a circumferential strip during traction test.
Test System with MTEST Quattro controller and Interface WMC10 (50 N) load cell at room temperature, as shown in Fig. 2(d). Both tensile and relaxation tests were carried out on axial and circumferential strips after 10 preconditioning cycles at the same strain rate and reaching the same peak stress of the actual tests. Several strips were tested from the same graft. A similar outcome was obtained from different samples; results are reported for one axial (initial length between the jaws 30.17 mm, with standard deviation of 0.05 mm) and one circumferential strip (initial length between the jaws 30.31 mm, with standard deviation of 0.05 mm). During the tensile tests, a strain rate of 0.025 mm/s was used. The distance between the jaws and the tensile forces were recorded directly by the test system with accuracy of ± 1 µm and ± 0.01 N, respectively. The relaxation test was performed with a fast-loading linear tensile ramp, obtained for constant strain rate about 0.825 mm/s corresponding to the maximum rate possible for the test system, followed by a hold at constant strain for 10,000 s.
Piezo load cell Shaker
Laser spot Strip holder Fixed surface
Fig. 3. Experimental setup for the dynamic stiffness measurement of a circumferential strip.
The graft and the corresponding strips are corrugated in axial direction. Since these corrugations tend to permanently decrease when the prosthesis is implanted or loaded for a long enough time, the graft was left submerged in water with a weight of 420 g for 4 weeks before any testing. The peak-to-peak average amplitude of the residual corrugations is 0.52 mm and the wavelength is 2.32 mm in the axial direction, obtained from image processing. Since the material (and not the global properties of the corrugated tube) is of interest in the present study, the axial strip was modestly preloaded during all tests to remove the residual corrugation.
2.2.3. Dynamic stiffness test The experimental setup for the dynamic stiffness measurement is shown in Fig. 3. A circumferential/axial strip was glued by cyanoacrylate to two L-shaped strip holders made of aluminum to minimize inertia. It was verified that they replace the jaws with great accuracy for the applied loads and the glue did not introduce any significant change in the stiffness. A strip holder was connected to a fixed surface through a B&K 8203 piezoelectric load cell. The other strip holder was fixed to the armature of a B&K 4824 electrodynamic exciter (shaker), connected to a B&K 2732 power amplifier and a LMS/Siemens Test-Lab system, and its displacement was measured by a Polytec OFV-505 laser Doppler vibrometer with displacement decoder. The initial free length between
2.2.2. Uniaxial tensile and relaxation tests Uniaxial tensile tests were performed on an eXpert 4000 Admet Micro 284
Journal of the Mechanical Behavior of Biomedical Materials 82 (2018) 282–290
M. Amabili et al.
the strip holders was 30.05 mm, with a standard deviation of 0.15 mm, for the axial strip and 30.22 mm, with a standard deviation of 0.24 mm, for the circumferential strip. Initially, a preconditioning of the strip was performed by increasing and decreasing the DC voltage applied to the shaker for 10 cycles. After that, a constant DC voltage was supplied to ensure the desired stretch. Both strips were tested for three different stretch levels. After allowing 300 s for relaxation, a small AC voltage was added to generate a harmonic excitation of amplitude about 0.1 N, giving a dynamic elongation of the order of 10 µm, at various frequencies. The frequency range of harmonic excitation was from 1 to 60 Hz with 0.1–0.5 Hz adaptive resolution. At every frequency, the initial 50 cycles were discarded to eliminate the transient vibration. The dynamic elongation of the strip Δl and the harmonic force Δf applied to the strip were measured for all the excitation frequencies for 10 cycles. From the measured dynamic force, elongation and the phase between them, the storage modulus and loss factor of the strips are calculated. The amplitude of the complex modulus E* is given by
Δf / A , Δl / L
E* =
(11)
where, L is the free length of the strip and A its cross-section area. The storage modulus E ′ and the loss modulus E ″ are given by
E ′ = E * cos ϕ,
E ″ = E * sin ϕ,
(12a,b)
where ϕ is the measured phase angle between the harmonic force and the dynamic elongation. The loss factor η is
η = tan ϕ.
(13) Fig. 4. Uniaxial tensile tests on strips from a woven Dacron aortic graft; 2nd Piola-Kirchhoff stresses, Green strains. The blue dotted curves are those taken to fit the material model. (a) Circumferential strip, 9 repeated tests, mean values and standard deviations; (b) axial strip, 9 repeated tests, mean values and standard deviations. (For the interpretation of references to colour in this figure legend, teh reader is referred to the web version of the article.)
2.3. Constitutive modelling 2.3.1. Hyperelastic law The following anisotropic Fung-type strain energy density function is applied to the woven Dacron 2
W=
∑ Ci (exp(Qi) − 1),
2 Qi = 4μi (εxx cos4 βi + 2εxx εθθ cos2 βi sin2 βi
Eq. (16) is minimized by a genetic algorithm (Mitchell, 1996), which can be briefly described as follows. We consider a population of AT candidate solutions (called “individuals”), each of them represented by a “genetic code”, i.e. a vector of 6 parameters {C1, μ1 , β1, C2, μ 2 , β2} . At the initial stage, a large enough pool of AT individuals with random genetic codes is created. Then, iterations of the algorithm start. Initially, for every individual “fitting”, represented by Eq. (16), is evaluated. Afterwards, all individuals are arranged according to their fitting. The AS ( AS < < AT ) most fitted individuals, i.e. those with the smallest corresponding value of fs obtained by Eq. (16), “survive” and they are included in the list of candidates for the next iteration. Then, AT − AS “offspring” individuals are added to the pool. The offspring individuals are generated in two different ways: (i) from one single parent, randomly varying the parameters of every of the AS surviving individuals, and (ii) from two parents. In the latter case, the offspring genetic code contains parameters with random values comprised between those of the parents. The new pool of individuals contains AS survived and AT − AS offspring individuals, so the total number of candidate solutions in the pool remains AT . With this new pool, a new iteration starts. Iterations continue until the solution with “small enough value” of the error function (16) is found. Generally, the genetic algorithm is a good optimization method of complicated multidimensional problems with complex shape of the error function. The sensitivity of the algorithm to the initial guesses is low due to the very large number AT of candidate solutions in the initial pool.
i=1
+
2 εθθ
sin4 βi ),
(14)
where Ci > 0 are stress-like material parameters, μi are non-dimensional material parameters, and βi are angles with the axial direction. Eq. (14) is a formulation for a material with two families of different fibers. The woven Dacron is made of circumferential and axial fibers that are not completely independent due to weaving. The convexity of W is verified for Ci, μi > 0 (Holzapfel, 2006). The 2nd Piola-Kirchhoff principal stresses take the following expressions
σ͠ xx =
σ͠ θθ =
2
∂W = ∂εxx
∑ Ci exp(Qi)8μi (εxx cos4 βi + εθθ cos2 βi sin2 βi),
∂W = ∂εθθ
∑ Ci exp(Qi)8μi (εxx cos2 βi sin2 βi + εθθ sin4 βi).
i=1
(15a)
2 i=1
(15b)
In case of uniaxial extension in axial direction x, Eq. (15a) is used to fit the experimental data and (15b) is used to obtain εθθ as a function of εxx considering σ͠ θθ = 0 . Each one of Eq. (15a,15b) fits the corresponding experimental data set, in order to obtain the unknown material parameters Ci, μi ,βi . The material parameters that guarantee the best fit of the experimental data are obtained by minimizing the following nonlinear stress based function N
fs =
⎧ ∂W
∑ ⎨ ⎜⎛ ∂ε
n=1
⎩⎝
xx n
2
2
∂W (n) ⎞ − σxx +⎛⎜ ⎟ ∂ ⎠(axial) ⎝ εθθ
n
⎫ (n) ⎞ − σθθ , ⎟ ⎬ ⎠(circ ) ⎭
(16)
2.3.2. QLV for relaxation In order to fit relaxation experiments with the reduced relaxation function given in Eq. (8), the following procedure is applied
where n indicates the n-th experimental data and the subscripts “axial” and “circ” indicate which strip is considered (axial or circumferential). 285
Journal of the Mechanical Behavior of Biomedical Materials 82 (2018) 282–290
M. Amabili et al.
Fig. 5. Fitting of the hyperelastic material model to the experimental tensile test data. •, experimental data; —, material model; 2nd Piola-Kirchhoff stresses, Green strains. (a) Circumferential strip; (b) axial strip.
Table 1 Material parameters for the constitutive law of the woven Dacron obtained from tensile test. C1 (kPa)
μ1
β1 (rad)
C2 (kPa)
μ2
β2 (rad)
2.455
1216.51
− 1.233
25.891
49.28
0
Fig. 7. Fitting of the QLV model to relaxation experiments. •, experimental data; — , constitutive material model; 2nd Piola-Kirchhoff stresses versus time. (a) Circumferential strip, t0θθ = 0.55 s, γθθ = 0.030 s−1; (b) axial strip, t0xx = 1.535 s, γxx = 0.030 s−1. Table 2 Reduced relaxation parameters of the woven Dacron according to the QLV theory. Fig. 6. Relaxation tests on strips from a woven Dacron aortic graft; 2nd PiolaKirchhoff stresses versus time. (a) Circumferential strip, 8 repeated tests, mean values and standard deviations; (b) axial strip, 8 repeated tests, mean values and standard deviations.
τ1 (s)
τ2 (s)
a
Circumferential Axial
0.049 0.049
1545.4 2670
0.038 0.023
where γxx is the constant strain rate. The time derivative of (15a) in the interval 0 < t < t0xx is
(Abramowitch and Woo, 2004). The instantaneous hyperelastic stresses in axial and circumferential directions are given by Eq. (15a,15b). In the time interval 0 < t < t0xx , i.e. during the ramp for the axial strip, the strain versus time is given by
εxx (τ ) = γxx τ ,
Direction
(17) 286
Journal of the Mechanical Behavior of Biomedical Materials 82 (2018) 282–290
M. Amabili et al.
λ = 1.012 λ = 1.010 λ = 1.008
(a)
λ = 1.042 λ = 1.039 λ = 1.037
Fig. 8. Storage modulus versus loading frequency at three different initial prestretches λ; the error bar at 30 Hz quantify the estimated measurement error. (a) Circumferential strip; (b) axial strip.
(b) Fig. 10. Storage modulus versus loading frequency at three different initial prestretches λ; ○, experimental data; —, generalized Maxwell material model with N = 6. (a) Circumferential strip; (b) axial strip.
λ = 1.008 λ = 1.010 λ = 1.012
(a)
λ = 1.037 λ = 1.039 λ = 1.042
Fig. 9. Loss factor versus loading frequency at three different initial pre-stretches λ; the error bar at 30 Hz quantify the estimated measurement error. (a) Circumferential strip; (b) axial strip.
(b) Fig. 11. Loss factor versus loading frequency at three different initial prestretches λ; ○, experimental data; —, generalized Maxwell material model with N = 6. (a) Circumferential strip; (b) axial strip.
287
Journal of the Mechanical Behavior of Biomedical Materials 82 (2018) 282–290
M. Amabili et al.
Table 3 Material parameters for the generalized Maxwell model with N = 6. Axial strip
λ = 1.008
Circumferential strip
λ = 1.037
λ = 1.039
λ = 1.042
λ = 1.008
λ = 1.01
β1∞
0.273
0.543
0.56
1.619
1.582
1.673
β2∞
0.183
0.19
0.17
0.9
0.772
0.64
λ = 1.010
λ = 1.012
λ = 1.012
β3∞
0.207
0.2
0.184
0.957
0.929
0.768
β4∞
0.011
0.004
0.006
0.033
0.04
0.053
β5∞
0.028
0.122
0.078
0.313
0.437
0.613
β6∞
0.202
0.321
0.474
1.422
1.477
2.162
τ1 τ2 τ3 τ4 τ5 τ6
0.001 0.016 0.136 7.058 59.803 1478.05
0.0006 0.015 0.122 6.923 91.364 1115.07
0.0006 0.014 0.113 7.07 90.309 1632.99
0.002 0.016 0.124 7.896 79.65 849.29
0.002 0.014 0.125 7.287 87.364 1418.6
0.001 0.012 0.105 0.334 104.68 1094.53
(a)
λ = 1.042
λ = 1.039 λ = 1.012
λ = 1.042
λ = 1.010 λ = 1.008
(b) Fig. 13. Loss factor versus loading frequency at three different initial prestretches λ; ○, experimental data; —, generalized Maxwell material model with N = 6 and the same parameters βα and τα for all the pre-stretches. (a) Circumferential strip; (b) axial strip.
(a)
Table 4 Material parameters for the generalized Maxwell model with N = 6 and the same parameters βα and τα for all the pre-stretches.
λ = 1.042 λ = 1.039 λ = 1.037
(b) Fig. 12. Storage modulus versus loading frequency at three different initial prestretches λ; ○, experimental data; —, generalized Maxwell material model with N = 6 and the same parameters βα and τα for all the pre-stretches. (a) Circumferential strip; (b) axial strip.
∂σ͠ xx = ∂τ
Axial strip
Circumferential strip
β1∞ β2∞ β3∞ β4∞ β5∞ β6∞
0.352
1.507
0.172
0.816
0.188
0.817
0.01
0.114
0.07
0.921
0.319
1.108
τ1 τ2 τ3 τ4 τ5 τ6
0.001 0.015 0.119 9.102 66.113 1100.32
0.002 0.019 0.136 5.214 67.991 1171.948
2
∑ 64μi2 Ci γxx exp(Qi) ⎡⎢ (γxx τ cos4 βi + εθθ (γxx τ )cos2 βi sin2 βi)2 i=1
σxx (t ) =
⎣
+ (γxx τ cos4 βi + εθθ (γxx τ )cos2 βi sin2 βi)(γxx τ cos2 βi sin2 βi
1 ′ (γxx τ )cos2 βi sin2 βi ) ⎤, (cos4 βi + 2εθθ ⎥ 8μi ⎦
∂σ͠ xx / ∂τ
t
1 + axx ln
( ) τ 2xx τ 1xx
⎧ ⎛t − τ ⎞ ⎛ t − τ ⎞⎤⎫ 1 + axx ⎡ ⎢g ⎜ τ2 ⎟ − g ⎜ τ1 ⎟ ⎥ ⎬ dτ , ⎨ xx ⎝ ⎠ ⎝ xx ⎠ ⎦ ⎭ ⎣ ⎩
0 < t < t0xx .
′ (γxx τ ) + εθθ (γxx τ )sin4 βi ) εθθ +
∫0
(19)
In the time interval t > t0xx , i.e. at constant strain εxx (τ ) = γxx t0xx , the stress is given by (18)
σxx (t ) =
where the prime represents the derivative with respect to the argument. Eq. (4) can be rewritten as
∫0
t 0xx
∂σ͠ xx / ∂τ 1 + axx ln
t > t0xx ,
( ) τ 2xx τ 1xx
⎧ ⎛t − τ ⎞ ⎛ t − τ ⎞⎤⎫ 1 + axx ⎡ ⎢g ⎜ τ2 ⎟ − g ⎜ τ1 ⎟ ⎥ ⎬ dτ , ⎨ xx ⎠ ⎝ xx ⎠ ⎦ ⎭ ⎣ ⎝ ⎩ (20)
since after t = t0xx the derivative of Eq. (15a) with respect to time is 288
Journal of the Mechanical Behavior of Biomedical Materials 82 (2018) 282–290
M. Amabili et al.
3.2. Uniaxial relaxation test
zero. Experiments give the values of the measured stresses Ri at discrete times ti.. Eqs. (19) or (20), where the choice between the two expressions is made according to the value of ti, is evaluated at those specific times, giving σi = σ (t ) t = ti . Then, the following least square objective function fxx, function of the reduced relaxation parameters, is built
The results of eight repeated relaxation tests are presented in Fig. 6(a,b) for the first few seconds (to highlight the ramp) with their mean values and standard deviations. The linear ramp is at constant strain rate γxx = γθθ = 0.030 s−1, and results are almost within the standard deviation of those obtained in Section 3.1 (mostly of the ramp curve is within the standard deviation range of the tensile test data). The different strain rate does not change the results significantly, so the material parameters in Table 1 are used to fit the relaxation experiments. The experimental results and the stress obtained by the constitutive material model versus time are compared in Fig. 7(a,b). Differences between the fitted material model and experiments in the time history are very small, except at the very beginning of the relaxation phase since the test machine cannot stop the ramp instantaneously, but it has a short deceleration phase. The peak value of stress in Fig. 7(a,b) corresponds to the maximum stress compatible with normal physiological conditions in thoracic aortic replacements. The 2nd Piola-Kirchhoff hoop stress σθ in a thin circular cylindrical shell under transmural pressure p and small cylindrical deformation is given by
N
fxx (axx , τ1xx , τ2xx ) =
∑ [Ri − σi]2 ,
(21)
i=1
where N is the number of experimental stresses measured. The three reduced relaxation parameters axx , τ1xx , τ2xx have been identified by minimizing Eq. (21) by using a genetic algorithm. A similar objective function is built for the circumferential stress.
2.3.3. Viscoelasticity for harmonic oscillation In this section, a continuum viscoelastic formulation with internal variables is applied (Holzapfel et al., 2002). The formulation is valid for finite strains and small perturbations away from the static equilibrium configuration. Therefore, it has hypotheses similar to the quasi-linear viscoelasticity previously introduced. It is here restricted to incompressible anisotropic materials. A generalized Maxwell material model is introduced in each of the two orthogonal in-plane directions with N being the number of parallel spring/dashpots arms. It is possible to derive the storage (dynamic) modulus E ′ of the material model as
⎡ E ′ (ω, λ ) = E∞ (λ ) ⎢1 + ⎣
N
(ω τ )2
∑ βα∞ 1 + (ωα τ
α)
α=1
⎤
2⎥
⎦
σθ =
N
∑α = 1 βα∞ N
ω τα 1 + (ω τα )2
1 + ∑α = 1 βα∞
(ω τα )2 1 + (ω τα )2
(24)
where the vessel radius R = 14 mm and thickness h = 0.327 mm are assumed. The stress value σθ = 0.685 MPa is estimated considering the peak value of the physiological pressure p = 120 mmHg (converted in Pa); this stress is close to the peak value in Fig. 7(a,b). The coefficients of the two reduced relaxation functions are given in Table 2. Inserting these values in Eq. (10), the reduced relaxations for t → ∞ are obtained; they are 0.718 and 0.797 in circumferential and axial direction, respectively.
, (22)
where E∞ is the static modulus, ω is the oscillation frequency, βα∞ > 0 are nondimensional stiffness coefficients and τα are time constants. Both E ′ and E∞ are functions of the initial static stretch λ of the material; in general, the same is true for βα∞ and τα . However, their variation with λ may be modest to the point that constant values can be assumed in some applications. The integer N must be chosen large enough to fit well the experimental data. The loss factor, which coincides with the loss tangent in case of linear viscoelasticity, is given by
η (ω) =
pR , h
3.3. Storage modulus and loss factor The storage modulus in circumferential and axial directions is plotted in Fig. 8(a,b) versus the loading frequency. Results for zero frequency, taken from the data in Section 3.1, are reported for comparison. The estimated measurement error at 30 Hz is indicated. It is evident that a very large increase in stiffness is obtained in case of harmonic loading with respect to the static loading. The stiffness moderately increases with the loading frequency between 1 and 60 Hz, which is the range of interest in physiological applications (the heart rate can vary from 0.66 Hz to 3.33 Hz, but the pressure and flow pulse waves introduce higher harmonics with frequency several times higher than the heart rate). Three curves are presented for each direction at different initial stretches since the material is hyperelastic; the modulus increases with the stretch, as expected. The corresponding loss factors are reported in Fig. 9(a,b); they are relatively flat between 1 and 60 Hz. The loss factor of the circumferential strip is larger (about double) than the axial strip. No significant difference in the loss factor is observed at different stretch for the axial strip, while a decrease with increased stretch is observed in circumferential direction. The storage modulus and loss factor are necessary for dynamic models of aortic prostheses (Tubaldi et al., 2018; Jayendiran et al., 2018). The experimental results and the storage moduli obtained by the generalized Maxwell material model with N = 6 versus frequency are compared in Fig. 10(a,b), while the corresponding loss factors are presented in Fig. 11(a,b). The material model fits with good accuracy the experiments for both the storage modulus and loss factors for the three pre-stretch levels investigated. In particular, the circumferential strip is better approximated than the axial strip. The corresponding material parameters βα and τα are given in Table 3. For each direction, the parameters for the three pre-stretch levels are somehow different; this is due to the simultaneous optimization of βα and τα for the best fitting. However, it is possible to use the same values of βα and τα to fit
. (23)
A minimization procedure based on a genetic algorithm is used also in this case to identify the material parameters βα and τα , while E∞ (λ ) is computed from Eq. (15a,15b). Eqs. (22) and (23) are then specified for x and θ directions. 3. Results and discussion 3.1. Uniaxial tensile test The results of nine repeated uniaxial tensile tests are presented in Fig. 4(a,b) and show very different stiffness values in circumferential and axial directions and relatively small standard deviation, indicating good repeatability after preconditioning. The experimental data from the blue dotted curves in Fig. 4(a,b) and the hyperelastic material model are presented in Fig. 5(a,b) for Dacron in circumferential and axial directions, respectively. The identified hyperelastic material parameters are given in Table 1. A good agreement between the hyperelastic law and the experimental data is observed in Fig. 5(a,b) for the full range of values presented, which covers the full range of physiological σ-ε values for an aortic implant. Results agree to those by Lee and Wilson (1986) for the circumferential direction, while in axial direction the corrugation has been removed in the present study by applying a small initial pre-stress, so the comparison is not possible. 289
Journal of the Mechanical Behavior of Biomedical Materials 82 (2018) 282–290
M. Amabili et al.
the three stretch levels. This suggests the possibility of building a material model with ratio between dynamic and static modulus and loss factor independent of the initial pre-stretch, as introduced in (Holzapfel et al., 2002), with great simplification. The results of this simplified material model are presented in Fig. 12(a,b) for the storage moduli and in Fig. 13(a,b) for the loss factors. In this case, a single loss factor curve is obtained for all the pre-stretches. The material parameters identified in this case are given in Table 4. Clearly, the material model does not compare to the experimental results as well as the one with parameters βα and τα identified at each pre-stretch level, but it is of much simpler implementation in numerical codes.
biological materials. J. R. Soc. Interface 12, 20150707. Babaei, B., Velasquez-Mao, A.J., Thomopoulos, S., Elson, E.L., Abramowitch, S.D., Genin, G.M., 2017. Discrete quasi-linear viscoelastic damping analysis of connective tissues, and the biomechanics of stretching. J. Mech. Behav. Biomed. Mater. 69, 193–202. Bischoff, J.E., 2006. Reduced parameter formulation for incorporating fiber level viscoelasticity into tissue level biomechanical models. Annu. Biomed. Eng. 34, 1164–1172. Bustos, C.A., García-Herrera, C.M., Celentano, D.J., 2016. Mechanical characterisation of Dacron graft: experiments and numerical simulation. J. Biomech. 49, 13–18. Castile, R.M., Skelly, N.W., Babaei, B., Borphy, R.H., Lake, S.P., 2016. Microstructural properties and mechanics vary between bundles of the human anterior cruciate ligament during stress-relaxation. J. Biomech. 49, 87–93. Dortmans, L.J.M.G., van de Ven, A.A.F., Sauren, A.A.H.J., 1994. A note on the reduced creep function corresponding to the quasi-linear visco-elastic model proposed by Fung. J. Biomech. Eng. 116, 373–375. Fung, Y.C., 1973. Biorheology of soft tissues. Biorheology 10, 139–155. Fung, Y.C., 1993. Biomechanics: Mechanical Properties of Living Tissues, 2nd ed. Springer, New York, USA. Giles, J.M., Black, A.E., Bischoff, J.E., 2007. Anomalous rate dependence of the preconditioned response of soft tissue during load controlled deformation. J. Biomech. 40, 777–785. Gimbel, J.A., Sarver, J.J., Soslowsky, L.J., 2004. The effect of overshooting the target strain on estimating viscoelastic properties from stress relaxation experiments. J. Biomech. Eng. 126, 844–848. Hasegawa, M., Azuma, T., 1979. Mechanical properties of synthetic arterial grafts. J. Biomech. 12, 509–517. Holzapfel, G.A., 2006. Determination of material models for arterial walls from uniaxial extension tests and histological structure. J. Theor. Biol. 238, 290–302. Holzapfel, G.A., Gasser, T.C., Stadler, M., 2002. A structural model for the viscoelastic behavior of arterial walls: continuum formulation and finite element analysis. Eur. J. Mech. A/Solids 21, 441–463. How, T.V., 1992. Mechanical properties of arteries and arterial grafts. In: Hastings, G.W. (Ed.), Cardiovascular Materials. Springer, Berlin, Germany, pp. 1–35. Jayendiran, R., Nour, B.M., Ruimi, A., 2018. Dacron graft as replacement to dissected aorta: a three-dimensional fluid-structure-interaction analysis. J. Mech. Behav. Biomed. Mater. 78, 329–341. Lee, J.M., Wilson, G.J., 1986. Anisotropic tensile viscoelastic properties of vascular graft materials tested at low strain rates. Biomaterials 7, 423–431. Mitchell, M., 1996. An Introduction to Genetic Algorithms. MIT Press, Cambridge, MA, USA. Puso, M.A., Weiss, J.A., 1998. Finite element implementation of anisotropic quasi-linear viscoelasticity using a discrete spectrum approximation. J. Biomech. Eng. 120, 62–70. Sarver, J.J., Robinson, P.S., Elliott, D.M., 2003. Methods for quasi-linear viscoelastic modeling of soft tissue: application to incremental stress–relaxation experiments. J. Biomech. Eng. 125, 754–758. Troyer, K.L., Estep, D.J., Puttlitz, C.M., 2012. Viscoelastic effects during loading play an integral role in soft tissue mechanics. Acta Biomater. 8, 234–243. Tubaldi, E., Païdoussis, M.P., Amabili, M., 2018. Nonlinear dynamics of Dacron aortic prostheses conveying pulsatile flow. J. Biomech. Eng. 240 art. 061004 (12 pages). Yeoman, M.S., Reddy, B.D., Bowles, H.C., Bezuidenhout, D., Zilla, P., Franz, T., 2010. A constitutive model for the warp-weft coupled non-linear behavior of knitted biomedical textiles. Biomaterials 31, 8484–8493.
4. Conclusions Eqs. (4) and (5) have probably been obtained in the present study for the first time; they allow a relevant simplification in case of direction-dependent viscoelasticity. A 11% difference of the reduced relaxations for t → ∞ between axial and circumferential directions has been observed for the woven Dacron; differences in the loss factors are much larger. This indicates that investigation of direction-dependent viscoelasticity is worth the effort, and the present simplified formulation of the three-dimensional QLV facilitates the applications. The experimental data for storage moduli and loss factors versus loading frequency and the corresponding material parameters obtained by the generalized Maxwell model are particularly relevant for dynamic models of aortic grafts. The dynamic modulus is much larger than the static modulus in both axial and circumferential directions, which is very significant in the evaluation of the radial deformation of aortic grafts under pulsatile blood pressure. Acknowledgments The authors acknowledge the financial support of the NSERC Discovery Grant (372493-13), the Canada Research Chair program and the Qatar grant NPRP 7–032-2–016. References Abramovitch, S.D., Woo, S.L.-Y., 2004. An improved method to analyze the stress relaxation of ligaments following a finite ramp time based on the quasi-linear viscoelastic theory. J. Biomech. Eng. 126, 92–97. Babaei, B., Abramowitch, S.D., Elson, E.L., Thomopoulos, S., Genin, G.M., 2015. A discrete spectral analysis for determining quasi-linear viscoelastic properties of
290