A stochastic-structurally based three dimensional finite-strain damage model for fibrous soft tissue
José F. Rodríguez, , Fernando Cacho, José A. Bea and Manuel Doblaré
Group of Structural Mechanics and Materials Modeling, Aragón Institute of Engineering Research (I3A), Torres Quevedo Building, María de Luna 3, Zaragoza 50018, Spain
Received 4 May 2005; revised 14 October 2005; accepted 17 October 2005. Available online 1 December 2005.
Abstract
A fully three dimensional finite-strain damage model for fibrous soft tissue is developed. The model assumes uncoupled contributions for the matrix and collagen fibers, and uncoupled bulk and deviatoric response over any range of deformations. A simple isotropic damage mechanism within the framework of continuum damage mechanics has been used to describe the softening behavior under deformation for the matrix. On the other hand, statistical aspects related to the length distribution of the reinforcing fibers lead to a damage model for the reinforcing material. As a result, a general theoretical framework for constitutive modeling of biological soft tissue with continuum damage is obtained. A theoretical example consisting of a biaxial test of a soft tissue reinforced with two families of collagen fibers has been considered to demonstrate the capabilities of the proposed model and to study the sensitivity to changes in the statistical parameters associated with the reinforcing material. Also, a preliminary numerical example is included to demonstrate the model on a inhomogeneous boundary value problem. Results show that the model is able to capture the typical stress–strain behavior observed in fibrous soft tissue and seems to confirm the soundness of the proposed formulation.
Keywords: Soft tissues; Continuous damage; Reinforced material; Statistical fiber distribution
Article Outline
1. Introduction
2. Basic results in continuum mechanics for hyperelastic media
3. Constitutive modeling of soft tissues
4. Elastic damage model for the matrix
5. Elastic damage model for the fibers
5.1. Bundle of fibers
5.2. Oriented fiber bundles
6. Numerical examples
6.1. Biaxial test: analytical example
6.2. Torsion–extension test: inhomogeneous boundary value problem
7. Discussion and conclusions
Acknowledgements
Appendix A. First derivatives for
Appendix B. Probability density function for
References
1. Introduction
Constitutive modeling of soft tissue has been an area of extensive research in the last few years (Fung, 1993, Holzapfel et al., 2000 and Humphrey, 2002). Accurate mechanical models of soft tissue coupled with appropriate numerical approaches can potentially aid in areas such as tissue engineering, in the study of atherosclerosis, heart functioning, and the simulation of surgical interventions or accident trauma. When modeling the mechanical behavior of soft tissue, particular difficulties arise not only from the highly nonlinear and anisotropic stress–strain response of the biological material, but also from the large variability in mechanical properties exhibited by these material (Sacks, 2000) due to its complex structure and composition. This variability can lead to wrong estimated values for the parameters in the constitutive equations as pointed by Chew et al. (1986) and Sacks (2003). As is usual in continuum mechanics for these type of material, the description of their constitutive behavior relies on the identification of an appropriate strain energy density function (SEF) from which stress–strain relations and local elasticity tensors can be derived. In the past, most SEFs proposed in soft tissue mechanics were purely phenomenological (Guccione and McCulloch, 1991, Fung, 1993 and Delfino et al., 1997). Even though these SEFs were successful for particular applications and for describing material properties, they were limited, in most cases, to a given range of physiological loading, depended on a large number of parameters, or were overly simplified and therefore only capable of describing partial aspects of the tissue's mechanical response. In general, these SEFs were unable to elucidate the underlying mechanics of soft tissue. In this regard, structural constitutive models attempt to integrate information of tissue composition and structure into the definition of the SEF in order to avoid ambiguities in material characterization, and to offer a better insight into the physical meaning of the material constants associated with these tissues. However, the degree of structure present in a SEF will depend on the level of detail used for its definition. A number of structural SEFs have been developed for a variety of soft tissues. Lanir (1983) developed the first structural theory for planar soft tissue, Farquhar et al. (1990) for cartilage, Belkoff and Haut (1991) for mature skin, Mijailovich et al. (1993) for lung, Weiss et al. (1996) for ligaments, Holzapfel et al. (2000) for arteries, and Billiar and Sacks (2000) for heart valves.
The large variability in structure and composition exhibited by biological soft tissue makes it necessary to include, to some extent, this information in the definition of the SEF. In this regard, Lanir (1983) developed a stochastic structural constitutive model for soft tissues. In this model, the total strain energy was assumed to be the sum of the individual fiber strain energies (all referred to the global tissue coordinates). The orientation of the fibers was modeled using a statistical distribution, whose parameters were estimated numerically by fitting mechanical tests data. Zulliger et al. (2004) have recently proposed a model which accounts for the composition of the tissue and the waviness of the collagen fibers in the definition of the strain energy function. This structural aspect of the soft tissue is modeled using a log-logistic probability density function with distribution parameters estimated numerically as in the case of Lanir (1983). In this model the fiber orientation is taken to be deterministic.
All models mentioned above have been successfully used to analyze the mechanical response of various soft tissues, as well as having offered valid alternatives to elucidate the role played by different structural aspects of the tissue on its macromechanical behavior. However, most of them are limited to analyzing the mechanical response in the elastic range and not through the failure region. Kwan and Woo, 1989 and Hurschler et al., 1997 and Liao and Belkoff (1999) have proposed models to describe damage in soft tissue;however, these models have been restricted to damage of the collagenous component only. Additionally, some of these models contain multiple parameters, and it is not clear if these parameters can be determined uniquely.
In this paper we propose a structural SEF to model damage in fibrous soft tissue where damage for the matrix and collagen fibers are accounted for separately. A simple isotropic damage model within the framework of continuum damage mechanics describes the softening behavior under deformation for the matrix. On the other hand, statistical aspects related to the length distribution of collagen fibers introduced in the formulation define the damage model for the reinforcing material. As a result, a general theoretical framework for constitutive modeling of biological soft tissue with continuum damage is obtained. This paper is organized as follows: the next section gives basic results from continuum mechanics necessary to develop a constitutive theory for hyperelastic materials. Sections 3–5 present the SEF developed in this work giving particular expressions for the uncoupled deviatoric and dilatational responses. Section 6 provides examples to demonstrate the capabilities of the developed SEF. Finally, Section 7 presents some discussion and concluding remarks.
2. Basic results in continuum mechanics for hyperelastic media
Let be a continuum body defined as a set of points in a certain assumed reference configuration. It will be also assumed that there exists a one-to-one mapping continuously differentiable (as well as its inverse χ-1) which puts into correspondence with some region , the deformed configuration, in the Euclidean space. This one-to-one mapping χ transforms a material point to a position in the deformed configuration.
The deformation gradient F is defined as
(1)
with J(X)=det(F)>0 the local volume ratio. It is sometimes useful to consider the multiplicative split of F
(2)
into dilatational and distortional (isochoric) parts, where I is the second-order identity tensor. Note that . From this, it is now possible to define the right and left Cauchy–Green deformation tensors, C and b, respectively, and their corresponding isochoric counterparts and
(3)
it being straightforward to show that
(4)
where is the fourth-order identity tensor, where and are non-standard dyadic product operators (Rüter and Stein, 2000) defined by
(5)
For isothermal and reversible processes, there exists a SEF, Ψ, from which the hyperelastic constitutive equations
(6)
for the second Piola–Kirchhoff stress tensor, and
(7)
for the elasticity tensor are obtained (Eringen, 1989).
A weighted push-forward operation of the second Piola–Kirchhoff stress tensor, S yields the Cauchy stress tensor σ,
(8)σ=J-1FSFT.
A similar weighted push-forward operation of the elasticity tensor leads to its spatial counterpart, ,
(9)
For further details and alternative formulations, the interested reader is referred to Holzapfel (2000).
3. Constitutive modeling of soft tissues
Fibrous soft tissues such as ligaments, cardiac muscle, arteries, and veins are materials composed primarily of connective tissue proteins, elastin and collagen, and smooth cells. These dense connective tissue consist mainly of fibered collagenous tissuegrouped in one or several families, embedded in a highly compliant solid matrix (i.e., a ground substance made of proteoglycans, water, collagen and glycoproteins) (Fung, 1993) as shown in Fig. 1. It is the alignment of collagen fibers along preferred directions which gives the typical anisotropic behavior to these material, with the solid matrix being responsible for their incompressible response.
Full-size image (18K)
Fig. 1. Fibrous structure of the soft tissue.
Based on this observation, most fibrous soft tissues are assumed to be a continuum fiber reinforced and sometimes layered material (Streeter and Bassett, 1966, Rhodin, 1979, Fung, 1993, Weiss et al., 1996 and Holzapfel et al., 2000). Additionally, since most soft tissues have a negligible volume change in the physiological range of deformation (Carew et al., 1968 and Costa et al., 1996), constitutive modeling is usually based on the definition of a quasi-incompressible hyperelastic fibered material. In this kind of constitutive model, the SEF is usually expressed with its deviatoric and volumetric contributions decoupled, hence
(10)
where U and are purely volumetric and isochoric contributions to the SEF, respectively, and n0 and m0 are unit vectors along preferred directions of the material as shown in Fig. 1. Directions n0 and m0 are the elements through which the tissue's structure is introduced into the model. Further, the isochoric component is additively split into a part associated with isotropic deformations and a part associated with anisotropic deformations (Holzapfel et al., 2000). The isotropic part is related to the mechanical response of the matrix, or the non-collagenous material. The anisotropic part, on the other hand, is responsible for the tissue's stiffness to stretch in certain directions, and is entirely due to the collagen fibers. Hence, the strain energy potential Ψ is written as
(11)
where n0n0 and m0m0 are symmetric structural tensors. Spencer (1980) showed that the irreducible integrity bases for the three symmetric second-order tensors , n0n0, and m0m0 are given in terms of the eight invariants
(12)
Invariants and are directly related to the deformation of the matrix, while and are the square of the stretches in the directions of n0 and m0, respectively, so they represent a direct stretch measure of the collagen fibers. and also arise from the anisotropy introduced by the two reinforcing families of fibers and are related to fiber shear deformation. These invariants are usually not included in the formulations directly due to their strong correlation with and which leads to an ill posed parameter estimation problem for most experiments. Finally, is related to the interaction between both families of fibers. Therefore, the strain energy function can be written as
(13)
Two basic assumptions are considered in the development that follows. First, the cross-sectional area of the fiber is small as compared to the fiber length, and second, all components of the soft tissue act in parallel (no significant interaction between fiber families exists). The first assumption implies that the effect of invariants and are captured by invariants and , while the second leads to neglecting the functional dependence with and allows the additive split of into
(14)
These observations have been pointed out by several authors when formulating structural SEFs for soft tissues (Holzapfel et al., 2000). Therefore, the basic SEF proposed here has the form
(15)
4. Elastic damage model for the matrix
In this section, the basic structure of the SEF for the matrix will be developed. It will be assumed that material damage associated with the matrix in soft tissues is related to maximum distortional energy and is independent of hydrostatic pressure (Simo, 1987). Therefore, following the model proposed by Simo (1987), let be a damage-dependent deviatoric free energy density function of the form
(16)
where D[0,1] is a damage parameter. From the Clausius–Duhem inequality, it follows that
(17)
(18)
where DEV[·]=[·]-1/3([·]:C)C-1. Inequality (17) tells that damage is a dissipative process and that is the thermodynamic variable conjugate to D. Therefore, the evolution of the damage parameter is characterized by an irreversible equation of evolution in terms of (Simo, 1987). In order to establish the law of evolution, first define an equivalent strain Ξs
(19)
where is the deviatoric strain tensor at time , and E0 a constant with units of stress and of the order of . Now, let be the maximum value of Ξs over the past history up to current time
(20)
With Ξs and at hand, we can now define a damage surface in the strain space
(21)
Hence, a non-increasing damage criterion in the strain space is defined by the condition . Denote , the normal to the damage surface in the strain space. Therefore, the following situations can occur: either
(22)
where stands for an admissible variation of the deviatoric strain tensor. So, according to (22), and borrowing terminology typically from classical plasticity, it will be spoken of unloading, neutral loading or loading from a damaged state, depending on the sign of . Finally, the rate of the damage variable D is specified by the irreversible rate equation
(23)
where is a given function which characterizes the damage process in the material and is specified from available experimental data. For the present model, we have adopted the following form for :
(24)
with α and β regarded as experimental parameters. It is straightforward to prove that (24) leads to
(25)
This particular law for the damage variable allows control of damage initiation as well as the rate at which damage develops with strain as shown in Fig. 2, where D has been plotted for different values of α and β. This characteristic of (25) is in good agreement with tensile experiments performed along the transverse direction on the interosseous ligament of the forearm (Stabile et al., 2004).
Full-size image (41K)
Fig. 2. Damage function for different values of parameters α and β.
5. Elastic damage model for the fibers
Most constitutive models developed for fibrous soft tissue assume all structural parameters related to the fiber arrangement to be deterministic. However, histological studies performed in a number of soft tissues (Sacks et al., 1994, Canham et al., 1997, Hsu et al., 1998 and Dingemans et al., 2000) have shown that collagen fibers appear to be wavy and distributed about preferential directions (Lanir, 1983), thus as the load is applied, more and more collagen fibers start to bear load. However, the degree of straightening of each fiber will also depend upon its orientation relative to the loading. In addition, it should be pointed out that the waviness, as well as the fiber preferred grouping directions are highly random parameters in real soft tissues. Hence, their distributions should be accounted for while developing an accurate SEF for these materials. In this paper we only focus on treating the wavy nature of collagen fibers letting the orientation of the fibers be a deterministic parameter. Future developments in this area will introduce this statistical aspect into the model.
5.1. Bundle of fibers
First, the fibrous part will be considered to be composed of a bundle of fibers of different length grouped such that their average orientation are perfectly aligned, as shown in Fig. 3.
Full-size image (16K)
Fig. 3. A bundle of aligned collagen fibers.
Each fiber in the bundle is assumed to behave following a worm-like chain model, with the SEF defined as
(26)
with
(27)
where εs is the deformation of the bundle at time , is the deformation at the fully extended length of the fibril, is a parameter related to collagen fibril failure, and C1 is a constant related to the density of collagen fibril and the persistence length of the fibril. Usually, this constant is given as C1=Nkθ/A, where k is the Boltzmann's constant, N the fibrils density, θ the temperature and A the persistence length of a fibril chain. This model was introduced by Kratky and Porod (1949) and has been used in modeling the stretching of type II collagen (Sun et al., 2004) and in molecular biology (Bustamante et al., 2003), as well as for modeling soft tissue growth (Garikipati et al., 2004).
Eqs. (26) and (27) consider all fibers to have the same length and be perfectly aligned along their average directions. Different lengths of collagen fibers are introduced in the model by considering as a random variable. The probability density function associated with is taken to be a Beta probability density function (Larson, 1982)
(28)
where γ and η are related to the mean, , and variance, , of the distribution as
(29)
We have chosen this particular probability density function since it appears as a posterior distribution for the probability of failure of collagen fibers in the tissue as shown in Appendix B.
Let now be the maximum strain over the past history up to time . The damage of the fiber bundle increases whenever since new fibers fail (those with ). For this damage criterion, the SEF of the fiber bundle can be identified with the addition of the means of individual fiber SEFs or, equivalently, with the mean of the whole SEF of the bundle, and written as
(30)
Note that definition (30) only accounts for those fibers which are in tension and have not failed; therefore it naturally introduces a damage mechanism in the formulation, with the amount of damage controlled by the probability density function (28). Reversal of the order of integration in (30) and using expressions (27) and (28) lead to the following expressions for and
(31)
(32)
where Beta(x,a,b) is the incomplete Beta function, F1(a,b,c,d,x,y) is the Appell hypergeometric function (Abramowitz and Stegun, 1970), and
(33)
5.2. Oriented fiber bundles
Eqs. (31)–(32) define the strain energy for fiber bundles when stretched along their longitudinal axis. In a general fibrous soft tissue, for a family of fibers aligned along a preferred direction n0, the fiber deformation at time s, is given as
(34)
Therefore, the anisotropic component of the strain energy function associated with this family of fibers is given by
(35)
Let the maximum strain over the past history up to time t be
(36)
Therefore, the non-increasing damage criterion in the strain space is defined by the condition that at any time of the loading process
(37)
with defining damaged surface in the strain space. Proceeding as in the case for the matrix, let , be the normal to the damage surface in the strain space. Similarly, the following situations can occur: either
(38)
where is again an arbitrary admissible variation of the deviatoric strain tensor.
The second Piola–Kirchhoff stress tensor at time t, is then
(39)
The Cauchy stress tensor is found by a weighted push-forward operation (8) of the previous expression. Equivalent expressions for the second family of fibers are obtained by substituting by .
6. Numerical examples
The theory developed in the previous sections has been implemented in a numerical formulation in order to explore the capabilities of the model as well as to perform a sensitivity analysis for the different model parameters.1 Two examples are presented in this section. First, a homogeneous boundary value problem consisting on a biaxial test of planar tissue with two families of fibers arranged symmetrically about the longitudinal axis of the specimen. The second example corresponds to an inhomogeneous boundary value problem. For this case, a fully three dimensional finite element simulation of a nonuniform cross-section specimen subjected to tension and torsion is considered.
6.1. Biaxial test: analytical example
Fig. 4 shows the fiber reinforced planar tissue under homogeneous deformation. The specimen has two families of fibers symmetrically distributed about its longitudinal axis.
Full-size image (18K)
Fig. 4. Biaxial test of planar tissue.
This particular example has been chosen for its versatility since it allows full testing of the constitutive model for the matrix (θ=90) and fibers (θ=0) independently, as well as their individual contributions for different fiber orientation, giving ample room for numerical testing. It is important to point out that for θ=0, the problem reduces to a planar tissue with one family of fibers only. The resulting mathematical formulation has been implemented and solved in Matlab 6.0 (The MathWorks, Inc.).
The material will be treated as incompressible , with the matrix modeled using the strain energy proposed by Ballyk et al. (1998)
(40)
where c1m and c2m are material constants. The Cauchy stress tensor is now
(41)
where n and m are the fiber directions in the current configuration, and we have used the fact that the orientation of the fibers is symmetric about the longitudinal axis (i.e., ) and J=1.
For the simulations, it will be taken λ2=1.0 which leads to , and σ=diag(σ11,σ22,0), with
(42)
where is given in Appendix A.
Table 1 summarizes the value for the different model parameters used in the computations. c1m, c2m, α and β have been obtained by fitting stress–strain data of transverse interosseous ligament specimens from Stabile et al. (2004). Similarly, stress–strain data from longitudinal interosseous ligament specimens from Stabile et al. (2004) have been used to fit values for C1, , and κ. The values for , and have been assumed based on experimental studies conducted by Sun et al. (2004) on collagen and the findings of Haut (1983) on tensile experiments conducted on rat tail tendons.
Table 1. Model parameters used for the computations
Parameter Symbol Value Units
Matrix material constant c1m 86.5 kPa
Matrix material constant c2m 69.5 kPa
Matrix damage constant α 0.49 —
Matrix damage constant β 5.78 —
Fibers material constant C1 19.47 MPa
Fiber bundle failure parameter 1.08 —
Fiber bundle maximum strain mean 0.125 —
Fiber bundle maximum strain standard deviation 0.045 —
Upper limit for the fiber bundle maximum strain κ 0.1863 —
Fig. 5 shows the fiber and matrix contribution to stress σ11 for the case when the fibers are along the longitudinal axis. As expected, the fibrous part dominates the mechanical response of the tissue in this case, since the contribution of the matrix is four orders of magnitude less than the contribution of the fibers.
Full-size image (40K)
Fig. 5. Fiber and matrix contribution to σ11 for θ=0. The specimen is stretched to λ1=1.15 and then stretched back to λ1=1.0.
Fig. 6 shows longitudinal and transverse stresses, σ11 and σ22, respectively, for different fiber orientations. These figures show how as the orientation angle increases, the amount of damage induced in the fibers decreases since the collagen contributes less to load bearing. Note, however, that the material response is similar in both directions, with basically a change in the magnitude of the stress. In this regard, it is important to point out the lack of symmetry at 45 shown in the figure. It is because of the fact that it is a biaxial test with λ2=1, and not an equi-biaxial test where symmetry would therefore be obtained for this angle.
Full-size image (43K)
Fig. 6. σ11 and σ22 for different fiber orientation angles. In all cases, the specimen is stretched up to λ1=1.15 and back to λ1=1.0.
Fig. 7 shows the response of the tissue (θ=10) under cyclic straining with increasing amplitude. The specimen is stretched up to λ1=1.15 and back to λ1=1.0 in the first cycle. In the second cycle, the specimen is stretched to λ1=1.165 and back to λ1=1.0. The results from this simulation clearly show the Mullins effect as well as a significant loss of stiffness as damage develops in the tissue. Note also that the stress–strain path followed by the material after restraining replicates the previous unloading path, a typical behavior for materials with no fading memory.
Full-size image (37K)
Fig. 7. σ11 for cyclic straining with increasing amplitude (θ=10). In the first cycle the specimen is stretched to λ1=1.15 and then stretched back to λ1=1.0. In the second cycle, the specimen is stretched up to λ1=1.165 and back to λ1=1.0.
To gain a better insight into the effect of the different model parameters on the predicted material response, a sensitivity analysis has been carried out. The study comprises the effect of the mean, , and standard deviation, , of the distribution , the fiber bundle failure parameter, , and the upper limit for the fiber bundle maximum strain, κ.
Fig. 8 depicts the sensitivity of σ11 to changes on the mean and standard deviation of the distribution for θ=10.
Full-size image (47K)
Fig. 8. σ11 for different values of and (θ=10).
The results show that the deformation at which the maximum stress occurs, as well as the maximum stress itself, increases as increases. However, the material stiffness exhibits a relatively small change. On the other hand, larger values of , lead to larger values of the deformation at which maximum stress occurs while the maximum stress in the stress–strain curve decreases. Also, the stiffness of the material (slope of the stress–strain curve) seems to decrease as increases. However, for relatively small strains (λ1<1.05) material stiffness does not appear to change significantly. This particular effect is expected since larger values of introduce more dispersion in the contribution to load bearing of each fibril within the bundle for a particular strain.
Fig. 9 shows the sensitivity of σ11 to changes in the parameter for θ=10.
Full-size image (33K)
Fig. 9. σ11 for different values of (θ=10).
A remarkable influence of this parameter on the material response is obtained. A considerable increment of the material stiffness, the maximum stress reached in the stress–strain curve, as well as the deformation at which this maximum stress occurs are observed as increases. This result is expected since as each fiber bundle bears load for a large amount of deformation, increasing the overall stiffness and strength of the tissue.
Fig. 10 shows σ11 for different values of κ for θ=10. The figure shows little influence of this parameter over the stiffness of the tissue but not over the strength. Also, this parameter appears to have considerable influence over the rate at which the tissue degrades once failure starts to occur (note the abrupt loss of stiffness for κ=0.176 after the maximum stress has been surpassed).
Full-size image (37K)
Fig. 10. σ11 for different values of κ (θ=10).
6.2. Torsion–extension test: inhomogeneous boundary value problem
In this example we consider a numerical simulation of a cylindrical specimen with nonuniform cross-section with fibers running longitudinally. The specimen is simultaneously subjected to an extension of 4% and an end-to-end rotation of 90. Fig. 11 shows a detail of the mesh (6534 linear tetrahedral elements). The theory discussed in the previous section has been implemented within the nonlinear finite element program ABAQUS. The details will be presented in a forthcoming communication that will focus upon computational aspects and further numerical examples.
Full-size image (21K)
Fig. 11. Finite element mesh of the specimen for the tensile–torsion test.
The data for the simulations correspond to Table 1. However, for the particular example considered in this paper, only damage in the fibers is considered (α=0.0). In what follows, damage of the fibrous part of the material is quantified as
(43)
This expression corresponds to the cumulative probability for the Beta probability density function (28) which relates to the amount of fiber failure within the material for a given loading condition.
Fig. 12 shows the maximum principal stress and distribution of damage in the specimen at 100% of the test strain. The figure shows a maximum stress concentration in the mid section of the specimen with damage growing from the outer surface to the interior of the specimen as expected.
Full-size image (58K)
Fig. 12. (a) Contours of principal tensile stress for the torsion–extension test for the 100% of the test strain. (b) Damage distribution in the fibers for the torsion–extension test (100% of the test strain).
Fig. 13 shows the maximum shear stress–strain and the damage evolution curves for the torsion–extension test for an element located in the mid-section of the specimen over its outer surface.
Full-size image (34K)
Fig. 13. Maximum shear stress–strain and damage evolution curves for the tensional–torsional test.
The figure shows the lost of stiffness experienced by the material after damage starts to develop within the material. The results also show the rapid growth of damage after an apparent threshold value of the strain is attained. This behavior is due to the statistical nature of the damage mechanism introduced by the model. It is important to point out that this threshold value is controlled by the mean, , and standard deviation, , of the probability density function (28).
7. Discussion and conclusions
We have developed a fully three dimensional finite-strain constitutive model for fibrous soft tissue accounting for damage in both the matrix and the fibers. Uncoupling between volumetric and deviatoric anisotropic response is obtained as a result of the multiplicative split of the deformation gradient into volumetric and deviatoric parts. A simple isotropic damage model within the framework of continuum damage mechanics is introduced in the model to incorporate damage of the matrix. On the other hand, damage of the fibrous part is incorporated through the statistical distribution of the deformation at the fully extended length of collagen fiber bundles, a structural characteristic of soft tissue. In both cases, damage is characterized by the maximum value previously attained by the strain energy of the undamaged material.
The model contains nine parameters, four of which are related to the structure of the tissue and five to the mechanical properties of the constituents. As indicated by the model analysis on a biaxial test, the stress–strain behavior observed in recent experiments conducted in fibrous soft tissue, appears to be captured by the developed framework. The model is well suited for predicting abrupt and prolonged failure regions of highly aligned collagenous tissue, a characteristic which has not been obtained in previous damage models for fibrous soft tissue (Liao and Belkoff, 1999). However, neither fiber–fiber interaction nor fiber–matrix interaction have been considered in the development of the present framework.
At the fiber level, this model differs from those proposed by Lanir (1983), Hurschler et al. (1997), Liao and Belkoff (1999) in which the fiber bundle is assumed to be linear elastic and wavy with constant rupture energy (e.g., all fibers fail when they reach a given limit strain). In our model the fiber bundle response is nonlinear with each fiber bundle having different strain energy to fracture (a fact observed on tests conducted on collagen, Sun et al., 2004). The latter aspect is introduced in the model by considering the rupture strain of the bundle as a random variable following a Beta probability density function. On the other hand, the fact that the Beta distribution is bounded ensures that the probability of a fiber bundle failing at 0 strain is 0, as well as for fiber bundles strained above the distribution limit κ, with damage occurring continuously as the bundle is stretched, a unique aspect of the present framework.
The sensitivity analysis conducted on the model shows a remarkable influence of the structural parameters on the predicted strength and stiffness of the tissue. Therefore, these parameters should be determined from tests performed on fiber bundles and from histological studies conducted on individual collagen fibers within a given bundle. This implies that model calibration requires working at two scale levels. At a mesoscopic scale where structural parameters are identified from individual fiber bundles, and at a macroscopic scale where material parameters are determined from uniaxial or biaxial tests conducted on the tissue. Also, since the number of model parameters to be identified at the macroscopic level is reduced, the likelihood that these parameters may be determined uniquely increases.
The numerical example presented demonstrates the stochastic-structurally based damage model presented in this paper. The example consists of a cylindrical specimen with nonuniform cross-section subjected to simultaneous torsion and extension. The results, while only preliminary, illustrate the model and seems to confirm the soundness of the present formulation.
Acknowledgements
Financial support for this research has been provided by the Spanish Ministry of Science and Technology through the research project DPI2004-07410-C03-01 and the Sixth Framework Programme through the project FP6-2002-SME-1-513226. We thank the reviewers for their helpful comments.
References
Abramowitz and Stegun, 1970 M. Abramowitz and I.A. Stegun, Handbook of Mathematical Functions, Dover, New York (1970).
Ballyk et al., 1998 P.D. Ballyk, C. Walsh, J. Butany and M. Ojha, Compliance mismatch may promote graft-artery hyperplasia by altering suture-line stresses, J. Biomech. 31 (1998), pp. 229–237. View Record in Scopus | Cited By in Scopus (84)
Belkoff and Haut, 1991 S.M. Belkoff and R.C. Haut, A structural model used to evaluate the changing microstructure of maturing rat skin, J. Biomech. 24 (1991) (8), pp. 711–720. Abstract | PDF (1623 K) | View Record in Scopus | Cited By in Scopus (32)
Billiar and Sacks, 2000 K.L. Billiar and M.S. Sacks, Biaxial mechanical properties of fresh and glutaraldehyde treated porcine aortic valve cusps: part II—a structurally guided constitutive model, ASME J. Biomech. Eng. 122 (2000) (4), pp. 327–335. Full Text via CrossRef
Bustamante et al., 2003 C. Bustamante, Z. Bryant and S.B. Smith, Ten years of tension: single-molecule DNA mechanics, Nature 421 (2003), pp. 423–427. Full Text via CrossRef | View Record in Scopus | Cited By in Scopus (438)
Canham et al., 1997 P.B. Canham, H.M. Finlay and D.R. Boughner, Contrasting structure of the saphenous vein and internal mammary artery used as coronary bypass vessels, Cardiovasc. Res. 34 (1997), pp. 557–567. View Record in Scopus | Cited By in Scopus (29)
Carew et al., 1968 T.E. Carew, R.N. Vaishnav and D.J. Patel, Compressibility of the arterial wall, Circ. Res. 23 (1968), pp. 61–68. View Record in Scopus | Cited By in Scopus (95)
Chew et al., 1986 P.H. Chew, F.C. Yin and S.L. Zeger, Biaxial stress–strain properties of canine pericardium, J. Mol. Cell. Cardiol. 18 (1986) (6), pp. 567–578. Abstract | PDF (846 K) | View Record in Scopus | Cited By in Scopus (28)
Costa et al., 1996 K.D. Costa, P.J. Hunter, L.K. Waldman, J.M. Guccione and A.D. McCulloc, A three-dimensional finite element method for large elastic deformations of ventricular myocardium: part II—prolate-spherical coordinates, J. Biomech. Eng. 118 (1996), pp. 464–472. Full Text via CrossRef | View Record in Scopus | Cited By in Scopus (84)
Delfino et al., 1997 A. Delfino, N. Stergiopulos, J.E. Moore and J.J. Meisters, Residual strain effects on the stress field in a thick wall finite element model of the human carotid bifurcation, J. Biomech. 30 (1997), pp. 777–786. Abstract | Article | PDF (4603 K) | View Record in Scopus | Cited By in Scopus (130)
Dingemans et al., 2000 K.P. Dingemans, P. Teeling, J.H. Lagendijk and A.E. Becker, Extracellular matrix of the human aortic media: an ultrastructural histochemical and immunohistochemical study of the adult aortic media, Anat. Rec. 258 (2000) (1), pp. 1–14. Full Text via CrossRef | View Record in Scopus | Cited By in Scopus (75)
Eringen, 1989 A.C. Eringen, Mechanics of Continua, Krieger, Florida (1989).
Farquhar et al., 1990 T. Farquhar, P.R. Dawson and P.A. Torzilli, A microstructural model for the anisotropic drained stiffness of articular cartilage, ASME J. Biomech. Eng. 112 (1990) (4), pp. 414–425. Full Text via CrossRef | View Record in Scopus | Cited By in Scopus (41)
Fung, 1993 Y.C. Fung, Biomechanics: Mechanical Properties of Living Tissue, Springer, New York (1993).
Garikipati et al., 2004 K. Garikipati, E.M. Arruda, K. Grosh, H. Narayanan and S. Calve, A continuum treatment of growth in biological tissue: the coupling of mass transport and mechanics, J. Mech. Phys. Solids 52 (2004) (7), pp. 1595–1625. Article | PDF (445 K) | View Record in Scopus | Cited By in Scopus (59)
Guccione and McCulloch, 1991 J.M. Guccione and A.D. McCulloch, Finite element modeling of ventricular mechanics. In: L. Glass, P. Hunter and A. McCulloch, Editors, Theory of Heart: Biomechanics, Biophysic, and Nonlinear Dynamics of Cardiac Function, Springer, New York (1991).
Haut, 1983 R.C. Haut, Age-dependent influence of strain rate on the tensile failure of rat-tail tendon, J. Biomech. Eng. 105 (1983), pp. 296–299. Full Text via CrossRef | View Record in Scopus | Cited By in Scopus (57)
Holzapfel, 2000 G.A. Holzapfel, Nonlinear Solid Mechanics. A Continuum Approach for Engineering, Wiley, Chichester (2000).
Holzapfel et al., 2000 G.A. Holzapfel, T.C. Gasser and R.W. Ogden, A new constitutive framework for arterial wall mechanics and a comparative study of material models, J. Elasticity 61 (2000), pp. 1–48. Full Text via CrossRef | View Record in Scopus | Cited By in Scopus (490)
Hsu et al., 1998 E.W. Hsu, A.L. Muzikant, S.A. Matulevicius, R.C. Penland and C.S. Henriquez, Magnetic resonance myocardial fiber-orientation mapping with direct histological correlation, Am. J. Physiol. 274 (1998) (Heart Circ. Physiol. 43), pp. H1627–H1634. View Record in Scopus | Cited By in Scopus (133)
Humphrey, 2002 J.D. Humphrey, Cardiovascular Solid Mechanics, Springer, New York (2002).
Hurschler et al., 1997 C. Hurschler, B. Loitz-Ramage and R.A. Vanderby, A structurally based stress–stretch relationship for tendon and ligament, J. Biomech. Eng. 119 (1997), pp. 392–399. Full Text via CrossRef | View Record in Scopus | Cited By in Scopus (63)
Kratky and Porod, 1949 O. Kratky and G. Porod, Rötgenuntersushung gelöster fagenmoleküle, Rec. Trav. Chim. Pays-Bas. 68 (1949), pp. 1106–1123.
Kwan and Woo, 1989 M.K. Kwan and S.L. Woo, A structural model to describe the nonlinear stress–strain behavior for parallel-fibered collagenous tissues, J. Biomech. Eng. 111 (1989), pp. 361–363. Full Text via CrossRef | View Record in Scopus | Cited By in Scopus (28)
Lanir, 1983 Y. Lanir, Constitutive equations for fibrous connective tissues, J. Biomech. 16 (1983), pp. 1–12. Abstract | PDF (1211 K) | View Record in Scopus | Cited By in Scopus (158)
Larson, 1982 H.J. Larson, Introduction to Probability Theory and Statistical Inference, Wiley, New York (1982).
Liao and Belkoff, 1999 H. Liao and S.M. Belkoff, A failure model for ligaments, J. Biomech. 32 (1999), pp. 183–188. Abstract | PDF (145 K) | View Record in Scopus | Cited By in Scopus (35)
Mijailovich et al., 1993 S.M. Mijailovich, D. Stamenovic and J.J. Fredberg, Toward a kinematic theory of connective tissue micromechanics, J. Appl. Physiol. 74 (1993) (2), pp. 665–681. View Record in Scopus | Cited By in Scopus (52)
Rhodin, 1979 Rhodin, J.A.G., 1979. Architecture of the vessel wall. In: Sparks Jr., H.V., Bohr, D.F., Somlyo, A.D., Geiger, S.R. (Eds.), Handbook of Physiology, The Cardiovascular System, vol. 2. American Physiological Society, Bethesda, Maryland, pp. 1–31.
Rüter and Stein, 2000 Rüter, M., Stein, E., 2000. Analysis, finite element computation and error estimation in transversely isotropic nearly incompressible finite elasticity. Comput. Methods Appl. Mech. Eng. 519–541.
Sacks, 2000 M.S. Sacks, Biaxial mechanical evaluation of planar biological materials, J. Elasticity 61 (2000), pp. 199–246. Full Text via CrossRef | View Record in Scopus | Cited By in Scopus (105)
Sacks, 2003 M.S. Sacks, Incorporation of experimentally-derived fiber orientation into a structural constitutive model for planar collagenous tissues, ASME J. Biomech. Eng. 125 (2003), pp. 280–287. Full Text via CrossRef | View Record in Scopus | Cited By in Scopus (86)
Sacks et al., 1994 M.S. Sacks, C.J. Chuong and R. More, Collagen fiber architecture of bovine pericardium, ASAIO J. 40 (1994), pp. M632–M637. Full Text via CrossRef | View Record in Scopus | Cited By in Scopus (49)
Simo, 1987 J.C. Simo, On a fully three-dimensional finite-strain viscoelastic damage model: formulation and computational aspects, Comput. Methods Appl. Mech. Eng. 60 (1987), pp. 153–173. Abstract | PDF (1672 K) | View Record in Scopus | Cited By in Scopus (231)
Spencer, 1980 A.J.M. Spencer, Continuum Mechanics, Longman Scientific & Technical, Essex (1980).
Stabile et al., 2004 K.J. Stabile, J. Pfaeffle, J.A. Weiss, K. Fischer and M.M. Tomaino, Bi-directional mechanical properties of the human forearm interosseus ligament, J. Orthop. Res. 22 (2004), pp. 607–612. Article | PDF (418 K) | Full Text via CrossRef | View Record in Scopus | Cited By in Scopus (14)
Streeter and Bassett, 1966 D.D.J. Streeter and D.L. Bassett, An engineering analysis of myocardial fiber orientation in pig's left ventricle in systole, Anat. Rec. 155 (1966), pp. 503–511. Full Text via CrossRef
Sun et al., 2004 Y.L. Sun, Z.P. Luo, A. Fertala and K.N. An, Stretching type II collagen with optical tweezers, J. Biomech. 37 (2004) (11), pp. 1665–1669. Article | PDF (243 K) | View Record in Scopus | Cited By in Scopus (45)
Weiss et al., 1996 J. Weiss, N.M. Bradley and S. Govindjee, Finite element implementation of incompressible, transversely isotropic hyperelasticity, Comput. Methods Appl. Mech. Eng. 135 (1996), pp. 107–128. Abstract | Article | PDF (1705 K) | View Record in Scopus | Cited By in Scopus (200)
Zulliger et al., 2004 M.A. Zulliger, P. Fridez, K. Hayashi and N. Stergiopulos, A strain energy function for arteries accounting for wall composition and structure, J. Biomech. 37 (2004) (7), pp. 989–1000. Abstract | PDF (547 K) | View Record in Scopus | Cited By in Scopus (81)
Appendix A. First derivatives for
According to the definition of the strain energy function (15), its derivatives with respect to invariants and are given by
(A.1)
and
(A.2)
Similarly, for the second derivative, the expressions for the derivative of and with respect to and required in (A.1) and (A.2) are
(A.3)
(A.4)
where i={4,6} and is either or . Also note that, the integral form of the Appell hypergeometric function has been used for convenience in the implementation.
Appendix B. Probability density function for
The functional form of the probability density function for the parameter will be calculated using bayesian statistics. When using these techniques, the unknown parameter θ is treated as a random variable with a known prior probability density function. The prior probability is related to the degree of believe regarding the unknown parameter before observing a sample of a random variable X which probability depends on the unknown parameter θ. Therefore, a prior distribution with small variance implies a low degree of uncertainty about its value. After definition of the prior distribution, the sample values are observed and used to calculate the posterior distribution for the unknown parameter using Bayes theorem as
(B.1)
where fΘ|X(θ|x) is the posterior distribution for the unknown parameter θ, fX,Θ(x,θ) is the joint probability for the sample and the unknown parameter, and fX(x) is the marginal probability of the sample. Therefore, the posterior distribution is just the conditional probability for Θ given the sample.
Following the previous procedure, let q be the probability of failure of an individual fiber within the tissue (our unknown parameter), and suppose that individual fiber failure occurs independently. In order to estimate q, we first consider a uniform prior probability density function for q in (0,1)
(B.2)
Note that this assumption on fQ(q) implies total ignorance about the parameter, or that all possible values of q have the same likelihood.
If we consider n fibers, with n large, and denote by X the total number of failures (our random variable), then X follows a binomial distribution with parameters n and q
(B.3)
The joint probability density function for X and Q is
(B.4)
the marginal for X is
(B.5)
for x=0,1,…,n. Therefore, all values of X have the same likelihood when averaging over all possible values of q, all of which have the same likelihood. Hence, the posterior distribution for q is
(B.6)
a Beta probability density function with parameters x+1 and n-x+1. To relate this result to it suffices to realize that q is given by
(B.7)
where is the strain at which the fiber fails, and F(·) its cumulative probability. Therefore, the probability density function for is related to a Beta probability density function.
Corresponding author. Tel.:+34 9767600005230; fax:+34 976762578.
