- Email: [email protected]
Influence of viscoelasticity of protein on the toughness of bone
J O U R N A L O F T H E M E C H A N I C A L B E H AV I O R O F B I O M E D I C A L M AT E R I A L S
3 (2010) 260–267
available at www.sciencedirect.com
journal homepage: www.elsevier.com/locate/jmbbm
Research paper
Influence of viscoelasticity of protein on the toughness of bone S. Anup a,∗ , S.M. Sivakumar a , G.K. Suraishkumar b a Department of Applied Mechanics, Indian Institute of Technology Madras, India b Department of Bio-technology, Indian Institute of Technology Madras, India
A R T I C L E
I N F O
A B S T R A C T
Article history:
Bone is an ultrafine composite of protein (collagen) and mineral (hydroxyapatite). An anal-
Received 23 May 2009
ysis to determine the influence of the viscoelasticity of protein on the toughness of bone at
Received in revised form
the ultrafine scale is conducted by developing a discrete lattice model appropriate for the
22 October 2009
ultrafine scale called the incremental continuous damage random fuse model (ICDRFM).
Accepted 24 October 2009
Collagen viscoelasticity at ultrafine scale is shown to contribute significantly to the tough-
Published online 4 November 2009
ness of bone. The results obtained are important in the design of biomimetic ultratough artificial composites. c 2009 Elsevier Ltd. All rights reserved.
1.
Introduction
Bio-composites such as bone, nacre and dentin possess excellent mechanical properties such as toughness (fracture energy), when compared to their constituents. Several reasons have been given by researchers as to why they exhibit such excellent properties. One of the explanations put forward is the presence of a hierarchical structure (Ji and Gao, 2004) with dimensions ranging from millimeters (macroscale) to the order of nm (ultrafine scale). The presence of sacrificial bonds in protein, leading to a saw tooth stress–strain curve has also been shown in the literature to be a mechanism which enhances the strength and toughness of bio-composites (Ji and Gao, 2004; Smith et al., 1999; Thompson et al., 2001). These studies point to the importance to be given to the investigations on the reasons for toughness of bio-composites. The information obtained is useful in the development of artificial composites as well as in selecting appropriate drugs to cure
diseases (in the case of bone) (Nalla et al., 2006; Anup et al., 2008). Among the various bio-composites, bone has been extensively studied due to its clinical importance. A more extensive bibliography on these studies has been given by Nalla et al. (2004). The ultrafine structure of bone is expected to play an important role in its overall mechanical behavior (Gao et al., 2003; Ji and Gao, 2004). Therefore, the scope of the study is limited to ultrafine structure level (of the order of 10−7 m to 10−8 m), where a network of collagen molecules inter-spaced with mineral crystals is considered (Fratzl et al., 2004; Ji and Gao, 2004). Platelets are made up of hydroxyapatite. Collagen molecules have a triple helix structure and each of these is made up of polypeptides. Collagen molecules contain sacrificial bonds which break instead of the back bone, and this contributes to toughness of bone (Thompson et al., 2001). The idealized arrangement assumed at the ultrafine level consisting of a staggered structure of mineral platelets embedded
∗ Corresponding address: Indian Institute of Space Science and Technology, Department of Aerospace Engineering, Thiruvananthapuram 695022, India. Tel.: +91 9847395059. E-mail addresses: [email protected], [email protected] (S. Anup). c 2009 Elsevier Ltd. All rights reserved. 1751-6161/$ - see front matter doi:10.1016/j.jmbbm.2009.10.007
J O U R N A L O F T H E M E C H A N I C A L B E H AV I O R O F B I O M E D I C A L M AT E R I A L S
a
2h_1
3 (2010) 260–267
261
b
PLATELET
PLATELET TENSION ELEMENTS
MATRIX
2l1
MATRIX SHEAR ELEMENTS
q
2l2
MATRIX TENSION ELEMENTS
h3
Fig. 1 – (a) The idealized arrangement at the ultrafine scale. (b) A part of the random fuse lattice used for simulation.
within a protein matrix is shown in Fig. 1 (a). The present study focuses on ultrafine scale, though the toughness of bone can be influenced by mechanisms, operating at other scales. However, the processes taking pace at ultrafine scale have a profound influence on the way bone behaves at the macroscale. Bio-composites like bone exhibit significant viscoelasticity (Sasaki et al., 1993; Iyo et al., 2004; Ji and Gao, 2004). Viscoelasticity can dissipate energy in materials during deformation and fracture (Rahulkumar et al., 2000), thereby providing higher toughness. Among the constituents of bone, collagen is found to be viscoelastic from experiments while mineral (hydroxyapatite) is found to be nearly linear elastic (Ji and Gao, 2004). One of the important ways in which the effect of viscoelasticity can be observed is the variation in the stress–strain response with strain rate (Wineman and Rajagopal, 2000). When failure of the material occurs, failure energy dissipation (toughness) varies with strain rate (Behiri and Bonfield, 1980) due to viscoelasticity. The viscoelastic energy dissipation in bone have been attributed to various mechanisms including viscous like cement line motion in osteon interfaces, thermoelastic damping, damping in lamallae and fluid flow in bone (Garner et al., 2000). Scant open literature is available on the relative contributions of these mechanisms in viscoelastic energy dissipation characteristics (toughness) in bone during failure. Collagen, being viscoelastic is a probable cause of such toughness in bone. Therefore, investigations into the role played by collagen in the viscoelastic failure dissipation in bone at the ultrafine scale is considered be an important research direction. Discrete lattice models have been proposed for simplified numerical as well as analytical studies of failure of materials with disorder (Sahimi, 2003; Batrouni, 1998; Herrmann and Roux, 1990). In discrete lattice models, the medium investigated is divided into a discrete set of elastic bonds. These elastic bonds fail at random threshold values (Nukala and Simunovic, 2005). The main features that characterize discrete lattice models are disorder (randomness), stress–strain characteristics and a rule by which each of the bonds (elements) of the lattice are broken. The stress–strain and
breaking response characteristic of each element corresponds to the response at the microstructural level in the material (Nukala and Simunovic, 2005). Discrete lattice models are expected to capture the essential physics of the problem including mechanisms that influence the toughness such as crack deviation, branching and stress concentration effects at crack tips (Nukala and Simunovic, 2005). A scalar form of the discrete lattice model, called the random fuse model can be used to simulate failure, where a voltage is applied to a network of resistors (de Arcangelis et al., 1985; Kahng et al., 1988). A resistor network is a scalar analog of the elastic medium. The analysis is simpler in such a model. However, the important characteristics of the problem are not lost due to this simplification (Nukala and Simunovic, 2005). A modified version of the random fuse model, the continuous damage random threshold fuse model (CDRFM), has also been used to simulate the failure in materials if damage evolution takes place (Zapperi et al., 1997). CDRFM is extremely useful for comparing the toughness of bio-composites such as nacre as demonstrated by Nukala and Simunovic (2005). However, almost all analysis involving continuous damage random fuse models are used for materials with linear constitutive laws only as could be seen from literature (Nukala and Simunovic, 2005; Hansen, 2005). In this context, a computationally inexpensive and easily implementable, but useful method would be to simulate failure under quasistatic loads employing a random fuse model, while incorporating viscoelastic nature of the constituents. Therefore, an incremental continuous damage random fuse model (ICDRFM) is developed in the present work to include viscoelastic constitutive law and is used to investigate the influence of collagen viscoelasticity on the toughness of bone at the ultrafine scale. Therefore, the objectives of the present study are to • Develop an incremental continuous damage random fuse model to simulate the failure of viscoelastic materials such as bone. • Investigate the influence of collagen viscoelasticity at the ultrafine scale on the toughness of bone.
262 2.
J O U R N A L O F T H E M E C H A N I C A L B E H AV I O R O F B I O M E D I C A L M AT E R I A L S
Formulation of the incremental CDRFM
The various dimensions used in the study are as given in Table 1 (Kotha et al., 2000; Rho et al., 1998). The dimensions are normalized with respect to the thickness of the platelets (2h_1). A unit normalized width, w is assumed in the out of plane direction. The mineral platelets are assumed to carry tensile stresses only and matrix is assumed to transfer shear stresses between platelets (Ji and Gao, 2004; Kotha et al., 2000). The shear stresses are assumed to be a constant, while the tensile stresses are assumed to change along the length of the platelets which are consistent with observations made by other researchers (Ji and Gao, 2004; Kotha et al., 2000). The random fuse models can take care of the redistribution of stresses due to the presence of platelets automatically (Nukala and Simunovic, 2005). The idealized staggered arrangement of mineral platelets embedded in a protein matrix is modeled as a square lattice network of L × L size. A part of the lattice network used for simulation is given in Fig. 1(b). A current–voltage analog of the model is employed (Nukala and Simunovic, 2005). This is done since the stress–strain relation for an elastic medium is similar to the current voltage relationship of a resistor network. Moreover, in the case of an elastic medium, the total force in any direction at any node is equal to zero (Force equilibrium). Similarly, in the case of an electrical network, sum of the electric current at any node is equal to zero (Kirchhoff’s Law). The matrix shear elements (shear fuses) connect the longer sides of mineral platelets while matrix tension elements (tension fuses) connect the shorter sides. Mineral platelets are modeled as mineral tension elements (tension fuses). The bond breaking thresholds (strengths) for both of these are prescribed to be constant. The conductance of an element is taken to be its stiffness. The platelet is assumed to be linear elastic and brittle, represented by tension elements. The conductance of a platelet tension element is made zero once the current passing through it exceeds the critical threshold value. The protein (collagen) in bone has been shown experimentally to contain sacrificial bonds that protect the collagen backbone and break one after the other in response to stress (Thompson et al., 2001). Matrix is assumed to be a material with evolving damage due to breaking of sacrificial bonds and is represented by shear elements with conductance decreasing in discrete steps (Nukala and Simunovic, 2005). Each time the current (stress) in the matrix shear fuse exceeds the critical threshold value, it is considered to have failed, the conductance is decreased by multiplying by a factor a where, 0 < a < 1, and this is termed matrix damage. The conductance of a matrix element is calculated to be ai Cm , where i is the number of times the matrix element has damaged. Once the maximum shear strain is reached, the conductance is made zero to simulate complete failure of the matrix and is termed as matrix tearing. Thus matrix failure consists of both matrix damage and matrix tearing. This stress–strain behavior of the matrix and the lattice network used in the present study is similar to the CDRFM model that has been employed by Nukala and Simunovic (2005) for the simulation of failure of nacre. In traditional CDRFM, numerically, a unit voltage difference is applied and current flowing through all the elements is found out. For each element, the ratio of the current
3 (2010) 260–267
flowing through each of the elements to the threshold value of each element is found out. The element having the maximum value of this ratio, rmax is considered to have failed (Delaplace et al., 1996). It is assumed that the time taken for current relaxation is much higher than the time taken for bond breakage. Hence, only one bond is broken at a time. Since the equations are linear, the external voltage and current are multiplied by 1/rmax , so that the current in the broken element is equal to its threshold value. This means that the externally applied voltage (displacement) is either increased or decreased so the stress in one of the bonds is just enough to break it. Once the breakage of a bond occurs, the stiffness matrix is altered and the equations are solved for the new system with the changed stiffness matrix. However, in the case of incremental CDRFM, which is developed in the present study, an incremental voltage difference is applied between the bus bars and the Kirchhoff’s equations are solved to find out the current (stress) in each of the elements of the lattice. If the stresses in none of the matrix elements exceeds the critical value, the voltage difference is again increased incrementally. This process is repeated until the current (stress) in at least one of the elements exceeds the critical threshold value. If the number of elements having a current above the threshold value is more than one, then the element having highest value of current is assumed to have failed. In that case, the stiffness of the particular failed element is reduced by a factor and the process is continued. This means that the voltage is always increased in incremental steps and is similar to a displacement controlled test, with a constant displacement rate. This process of breaking of bonds one by one goes on until the lattice separates completely into two parts. Every time, the equations are solved, the current passing through the system (Ic ) and voltage (Vc ) applied are stored. By plotting Ic /L vs Vc /L, the stress–strain graph is obtained. The toughness of the composite (bone) is found out from the simulation, in the sense of work of fracture, which is the area under the stress–strain curve multiplied by the length of the specimen. In the present work, the matrix material alone is assumed to fail since this mode of failure would provide maximum overall toughness (Anup et al., 2008). In addition to viscoelasticity, matrix is assumed to have evolving damage. This coupled behavior is approximated by the following steps. Initially, the matrix element is assumed to be elastic with evolving damage and the number of damages required to reach the maximum strain is found out, as given in literature (Nukala and Simunovic, 2005). In the next case, matrix is considered to be viscoelastic. The number of damages that a matrix element undergoes for the viscoelastic case is assumed to be the same as that found for the linear elastic case. This assumption is based on physical grounds, since the number of sacrificial bonds remains the same whether the material is considered to be elastic or viscoelastic. Therefore, when material property of the matrix is changed, the number of sacrificial bonds is kept a constant. Matrix is also assumed to be linear viscoelastic and hence is modeled as a standard solid with a single relaxation time. The response of viscoelastic materials subjected to load has been studied by making use of the finite element method. Zocher et al. (1997) proposed a methodology to simulate this viscoelastic response, which has been used to model the
J O U R N A L O F T H E M E C H A N I C A L B E H AV I O R O F B I O M E D I C A L M AT E R I A L S
263
3 (2010) 260–267
Table 1 – Dimensions at ultrafine scale: lengths in nm, area in nm2 (Kotha et al., 2000; Rho et al., 1998). Platelet area (2h1 )
Dimension Value Normalized value Number of elements
2.5 1 –
Gap (2l1 ) 15 6 6
response of viscoelastic materials as well as composites (Patel et al., 2007). This methodology has been modified by Dai and Sadd (2004) for simulating damage in viscoelastic materials. This modified form (Dai and Sadd, 2004), is adapted in the present work to simulate the viscoelastic response of the matrix elements, and is given in detail elsewhere (Anup et al., in preparation). The basics of the methodology is given below. For viscoelastic materials, the shear stress τ(t) for current time t is expressed as a hereditary integral (Wineman and Rajagopal, 2000; Zocher et al., 1997; Schreurs, 2007) as shown in the following equation. Z t τ(t) = G(t − s)γ(s)ds, ˙ (1) 0
where G(t) is the relaxation function of a matrix element at current time t and γ˙ is the strain rate. Putting G(t) = G(t)/h3 , and displacement rate u˙ = γh ˙ 3 for a matrix element, the Eq. (1) becomes Z t ˙ τ(t) = G(t − s)u(s)ds, (2) 0
where G(t) is the relaxation function of the shear stiffness (conductance). Here, G(t) is expressed as a Prony series as given by G(t) = Gα + Gd e−t/TR ,
(3)
where Gα is the shear stiffness corresponding to the long term elastic modulus, Gd is the term corresponding to the difference between short term and long term shear stiffness and TR is the relaxation time. The total time interval [0, t] is discretized into intervals (time steps). It is assumed that the intervals are of equal length. The hereditary integral is split as an integral over [0, tn ] and an integral over the current or last increment [tn , tn+1 = t] (Schreurs, 2007). The subscript for t indicates the number of the time step. Substituting Eq. (3) into Eq. (2) and simplifying, we get Z tn Z t t−τ − t−τ ˙ ˙ τ(t) = Gα u(t) + Gd e TR u(s)ds + Gd e TR u(s)ds. (4) tn
0
Therefore, the stress is obtained as − ∆t T
τ(t) = Gα u(t) + e where Pd =
TR ∆t
R
τd (tn ) + Gd Pd ∆u, !
1−e
− ∆t T R
.
(5)
Here, u is the total displacement. The stress increment ∆τ is obtained by subtraction as given below
Platelet length (2l2 )
Matrix thickness (h3 )
Matrix area (q)
50 20 20
1.35 0.54 1
2.5 1 –
relation, − ∆t T
τd (tn ) = Gdc Pd ∆u + τd (tn−1 )e
R
(8)
where Gdc is the current value of Gd , which decreases with the number of times the matrix element fails (damages). For maintaining equilibrium, the following procedure is adopted. The global current (force) is obtained through an appropriate assembly of element contributions as given by the following equation. [Ks ] {∆u} = {χ}
(9)
where [Ks ] is the global tangent stiffness (conductance) matrix and {χ} is the global force. The Eq. (9) is solved to obtain the nodal voltages. The displacement in each element is the voltage difference across an element. The procedure employed for the simulation is outlined below. 1. Form the tangent stiffness matrix Ks employing the value of Ts (See Eq. (7)). For iteration = 0. 2. Find out the displacements (voltage), ∆u1 using Eq. (9). For iteration = n + 1. 3. ∆un+1 , τ(tn ) and τd (tn ) are obtained from previous time step. 4. Calculate ∆τ employing Eq. (6). 5. Calculate τ(t) = τ(tn ) + ∆τ. 6. Check if any matrix element exceeds the threshold value; then if the number of total failures (damages) has been reached, stress in that element is made zero else reduce the stress in the failed element by multiplying by a factor a. Update Gd and Gα for the failed element. Change the tangent stiffness matrix incorporating change in stiffness of the failed element. 7. Calculate τd (tn+1 ) employing Eq. (8). 8. Find out the incremental displacements ∆un+2 using Eq. (9) 9. Go to the next iteration. Periodic boundary conditions are imposed in the horizontal direction to simulate an infinite system. This is a standard approach of RVE to obtain material continuity. An incrementally increasing voltage difference is applied between the top and the bottom bus bars. For solving the Kirchhoff’s equations obtained as a result, the Matlab package (MathWorks, 2000) is employed in this work.
∆τ = τ(tn+1 ) − τ(tn ), ∆τ = [Gα + Gd Pd ]∆u − (1 − e
− ∆t T R
)τd (tn ).
(6)
Tangent stiffness (conductance) Ts is given as Ts = [Gα + Gd Pd ].
(7)
The total stress is obtained by adding up the stress increments in each time step. Here τd (tn ) is given by the recursive
3.
Results and discussion
Experiments show that the maximum tensile strain that collagen can undergo is about 10% (Wang et al., 2002). In this work, without much loss of generality, the maximum shear strain, εmax that matrix can undergo is also taken to be
J O U R N A L O F T H E M E C H A N I C A L B E H AV I O R O F B I O M E D I C A L M AT E R I A L S
fore stiffness becomes zero, nmax =
log(ε0 /εmax ) ≈ 33, where log(a)
ε0 = h3 Ym /G (Nukala and Simunovic, 2005). The value of a is chosen in the present work to be 0.9 as an example, though other values of a does not seem to affect the results (Nukala and Simunovic, 2005). The matrix tension elements are assumed to carry no stresses and are given zero conductance values. The basic value of parameters relating to viscoelastic behavior of the matrix was obtained from the stress relaxation experiments conducted (Atkinson et al., 1999) on human patellar tendon. The values obtained from tendon could be assumed to be that of collagen since patellar tendon consists mostly of collagen. A linear viscoelastic standard solid (a Prony series with a single relaxation time) is used in the present work to fit the data and the parameters were found out. The values got from the fitted curve are, relaxation time, TR = 20, and the ratio of short term modulus to long term modulus (enhancement factor), λ = 4. The short term modulus is taken to be the shear modulus of collagen available in literature (3.3 GPa) (Cusack and Miller, 1979). This gives the values of long term modulus, Gα = 0.825 GPa. There is no data available for shear strength of collagen, and there is a wide variation in the tensile strength obtained by various researchers, a typical value being 20 MPa (Ji and Gao, 2004; Kotha and Guzelsu, 2007). Hence, this value is taken as the threshold matrix shear strength in this study. Young’s modulus of the platelet is given a value of 100 GPa (Ji and Gao, 2004). A system size of L = 16 is used for analysis. The stress–strain graphs obtained using ICDRFM and the linear method are compared as shown in Figs. 2 and 3. For a viscoelastic material, when relaxation time is very small, i.e. when t → 0, the second term in Eq. (3) becomes equal to Gd , and G(t) = Gα + Gd which is the long term modulus of the material. Therefore, the predictions of the viscoelastic model with t → 0 should agree with the predictions of a model with linear elastic stress–strain law having a modulus equal to the long term modulus of the viscoelastic material. Similarly, for large relaxation time, predictions of the viscoelastic model should agree with the predictions of a model with linear elastic stress–strain law having a modulus equal to the short term modulus of the viscoelastic material. ICDRFM developed in this study is compared with the predictions of a model with linear elastic stress–strain law and is found to agree for two extreme cases of small and large relaxation times. The fact that fracture energy (area under the stress–strain curves) obtained for these two cases for the viscoelastic model agrees with that obtained for corresponding model with linear elastic modulus is also consistent with the fact that viscoelastic energy dissipation is found to tend to zero for very high and small relaxation times (Wineman and Rajagopal, 2000). In Fig. 2, the stress–strain graph obtained by the linear method (shear modulus, Gm = 3.3 GPa) and that obtained by ICDRFM (material being viscoelastic, TR = 30 000 s) are compared and seen to be in good agreement. In Fig. 3, the stress–strain graph obtained by the linear method (shear
3 (2010) 260–267
0.1
0.08 STRESS (GPa)
10% (Nukala and Simunovic, 2005). Each time the matrix damages, stiffness is decreased to 90% of its previous value, to simulate the breaking of sacrificial bonds. This is accomplished by putting value of a = 0.9. This gives the number of damages that each matrix shear element undergoes be-
0.06
0.04
DIRECT METHOD (Nukala et al., 2005)
0.02
INCREMENTAL METHOD
0
0
1
2
3
4 5 STRAIN
6
7
8 × 10–3
Fig. 2 – Comparison of stress–strain graphs for linear elastic case and viscoelastic case with very large relaxation time. The toughness (area under the curve) are in agreement for the two cases.
0.1
0.08 STRESS (GPa)
264
DIRECT METHOD (Nukala et al., 2005)
0.06
INCREMENTAL METHOD 0.04
0.02
0
0
0.005
0.01
0.015 STRAIN
0.02
0.025
Fig. 3 – Comparison of stress–strain graphs for linear elastic case and viscoelastic case with very small relaxation time. The toughness (area under the curve) are in agreement for the two cases.
modulus, Gm = 0.825 GPa) and that obtained by ICDRFM (material being viscoelastic, TR = 1/10 s) are compared and seen to be in good agreement. Fig. 4 shows the variation of fracture energy dissipation with relaxation time. From Fig. 4, it could be seen that the fracture energy dissipation (area under the stress–strain curve) increases with decrease in relaxation time, reaches a maximum and then decreases. This trend is consistent with the case during deformation of other viscoelastic materials (Wineman and Rajagopal, 2000). Therefore, these results indicate that ICDRFM developed in this work is a valid method for simulating failure, and for comparison of toughness (Fracture energy) values. Fig. 5 shows the effect of strain rate on fracture energy dissipation. The viscous component of fracture energy dissipation is almost zero below a strain rate of 6.25 × 10−6 s−1 . Increase in strain rate increases the fracture energy dissipation until a strain rate of 6.25 × 10−4 s−1 . Thereafter, with increase in strain rate, fracture energy dissipation (toughness) decreases considerably. Therefore, the toughness of bone at ultrafine scale shows an increasing–decreasing trend with increase in strain rate.
J O U R N A L O F T H E M E C H A N I C A L B E H AV I O R O F B I O M E D I C A L M AT E R I A L S
3
0.1
2.5
0.08
2 0.06
W / W0
STRESS (GPa)
265
3 (2010) 260–267
0.04
1.5 1
0.02
0
0
EXPERIMENT (Behiri et al.,1980)
0.5
0.005 0.01
0.015 0.02 0.025 0.03 STRAIN
0.035 0.04
Fig. 4 – Stress–strain curves for different values of relaxation time. Note that area under the stress–strain curve increases, reaches a maximum value and then decreases.
0
SIMULATION 0
0.2
0.4 0.6 STRAIN RATE ( s–1)
0.8
1 × 10–3
Fig. 6 – The variation of work of fracture with strain rate is shown in the figure. The work of fracture increases with strain rate, reaches a maximum and then decreases. Especially note the strain rate at which maximum work of fracture occurs, which is close together for the two curves.
0.1
STRESS (GPa)
0.08
0.06
STRAIN RATE 6.25 ×10–6 s–1 STRAIN RATE 6.25 ×10–5 s–1
0.04
STRAIN RATE 3.12 ×10–4 s–1 STRAIN RATE 6.25 ×10–4 s–1
0.02
0
0
0.005 0.01
0.015 0.02 0.025 0.03 STRAIN
0.035 0.04
Fig. 5 – Stress–strain curves for different values of strain rate. Note that area under the stress–strain curve increases, reaches a maximum value and then decreases.
The variation in work of fracture (fracture energy dissipation) in the case of bovine bone for different cross head speeds was given by Behiri and Bonfield (1980). Taking the distance between loading platens in their work to be approximately 20 mm, the strain rates can be calculated to vary from 8.5 × 10−6 s−1 to 8.5 × 10−4 s−1 . Increase in strain rate to a value of 8.5 × 10−4 s−1 , produced a considerable decrease in the work of fracture, according to their experimental results (Behiri and Bonfield, 1980). The results of the present study are compared with that of Behiri and Bonfield (1980) in Fig. 6. Here, W is the work of fracture obtained as the area under the stress–strain curve and W0 is the value of the work of fracture for the smallest strain rate (6.25×10−6 s−1 ). The simulated results of the present study and experimental results obtained from literature (Behiri and Bonfield, 1980) follow the same trend, with a good match over the strain rate at which the peak energy dissipation is obtained. Natali and Meroi (1989) conducted experiments on bone under a larger range of strain rates (from 10−5 s−1 to 1 s−1 ). The maximum values of work dissipated in their study is seen to be in the 10−3 s−1 to 10−2 s−1 range which is the same strain rate range over which maximum viscoelastic dissipation is observed in the simulation in the present
study. Melnis and Knets (1982) also have data relating to the influence of strain rate on the specific deformation energy over a range of moisture contents. For smaller moisture contents (2.5%), a maximum value of the deformation energy is reached at a strain rate of 10−3 s−1 which is similar to that obtained in the simulation in the present study. These results suggest that collagen viscoelasticity at ultrafine scale contributes significantly to the energy dissipation during fracture in bone. The viscoelastic dissipation in bone is attributed to various mechanisms such as viscous like cement line motion in osteon interfaces, thermoelastic damping, damping in lamallae and fluid flow in bone (Garner et al., 2000). The results obtained from the present study is significant since it bring forth evidence in support of collagen viscoelasticity at the ultrafine scale amongst the various mechanisms suggested (Garner et al., 2000), as the important contributor to viscoelastic fracture energy dissipation (toughness) in bone. This finding is significant since artificial composites could be modeled based on the structure and properties of bone at ultrafine scale (staggered arrangement of the mineral platelets embedded in a viscoelastic matrix) so that maximum energy dissipation and hence toughness is achieved. The model employed in the present work is an over simplified version of the processes taking place at the ultrafine structure in bone. The viscoelastic properties of collagen is modeled using a standard linear solid rather than nonlinear viscoelastic material models that are more appropriate for other such biological materials. Nevertheless, this simple model gives insights into the way in which viscoelasticity affects fracture energy dissipation. The requirement of matrix mode of failure so as to provide maximum toughness has been described earlier (Anup et al., 2008). There seems to be an additional requirement for matrix mode of failure in bone at the ultrafine scale so as to make use of the viscoelasticity in the fracture energy dissipation. The structural arrangement at the ultrafine scale, while helping to maintain a stiffness close to that of the mineral (Ji and Gao, 2004), aids in higher fracture energy dissipation, by permitting matrix mode of failure alone. The model takes into account the
266
J O U R N A L O F T H E M E C H A N I C A L B E H AV I O R O F B I O M E D I C A L M AT E R I A L S
mechanisms acting at the ultrafine scale only. However, the effect of mechanisms acting on other scales is seen to be not appreciable given the fact that the macroscopic experimental and simulated curves follow the same trend. Moreover, many researchers have put forward the idea that the mechanisms at the ultrafine scale predominantly affect the stress–strain response at the macroscale (Kotha and Guzelsu, 2002). Since ultrafine scale properties are being considered in this study, use of information obtained from nanoindentation experiments at various loading rates could have been more appropriate. The data related to nanoindentation on bone is extensive. However, the studies have been mostly conducted with respect to how the properties vary with respect to different tissues, composition and site rather than loading rate (Ebenstein and Puritt, 2006). However, there is scant open literature on the effect on loading rate on the toughness of bone. It should be noted that apart from the modes of failure that has been assumed in the present study (platelet breakage and matrix failure), there could be more failure modes like debonding at the interface between the matrix and the mineral platelet. The platelets and matrix were assumed to be isotropic, though in reality they may possess anisotropic properties. Research on the effect of other modes of failure and anisotropy on the toughness of bone is being pursued currently.
4.
Conclusions
In this work, an incremental continuous damage random fuse model is developed for use at the ultrafine scale in bone and is employed to examine the way in which viscoelasticity of protein influences the toughness of bone at ultrafine scale. The method developed is significant in the fact that the viscoelastic properties of the matrix could also be included in the model, while retaining a simple and computationally inexpensive numerical analysis. The toughness of bone at ultrafine scale shows an increasing–decreasing trend with increase in strain rate. There is a very good correlation between experimental results obtained from the literature and the simulation results of the present study, especially on the strain rate at which maximum toughness is achieved. Therefore, the results obtained from the present study bring forth evidence in support of the collagen viscoelasticity as one of the important contributors to the toughness (fracture energy) in bone. The results are significant in the design of biomimetic composites.
Acknowledgement SA gratefully acknowledges the financial support under the QIP scheme of the AICTE, New Delhi, India. REFERENCES
Ji, B., Gao, H., 2004. Mechanical properties of nanostructure of biological materials. J. Mech. Phys. Solids 52 (9), 1963–1990.
3 (2010) 260–267
Smith, B., Schaffer, T., Viani, M., Thompson, J., Frederick, N., Kindt, J., Belcher, , Stucky, G., Morse, D., Hansma, P., 1999. Molecular mechanistic origin of the toughness of natural adhesives, fibres and composites. Nature 399 (6738), 761–763. Thompson, J., Kindt, J., Drake, B., Hansma, H., Morse, D., Hansma, P., 2001. Bone indentation recovery time correlates with bond reforming time. Nature 414 (6865), 773–776. Nalla, R.K., Kruzic, J.J., Kinney, J.H., Balooch, M., Ager III, J.W., Ritchie, R.O., 2006. Role of microstructure in the aging-related deterioration of the toughness of human cortical bone. Mat. Sci. Eng. C-Bio. S 26 (8), 1251–1260. Anup, S., Sivakumar, S.M., Suraishkumar, G.K., 2008. Influence of relative strength of constituents on the overall strength and toughness of bone. J. Mech. Med. Biol. 8 (4), 527–539. Nalla, R.K., Kruzic, J.J., Ritchie, R.O., 2004. On the origin of the toughness of mineralized tissue: Microcracking or crack bridging?. Bone 34 (5), 790–798. Gao, H., Ji, B., Jager, I.L., Arzt, E., Fratzl, P., 2003. From the Cover: Materials become insensitive to flaws at nanoscale: Lessons from nature. PP Natl. Acad. Sci. USA 100 (10), 5597–5600. Fratzl, P., Gupta, H.S., Paschalis, E.P., Roschger, P., 2004. Structure and mechanical quality of the collagen-mineral nanocomposite in bone. J. Mater Chem. 14 (14), 2115–2123. Sasaki, N., Nakayama, Y., Yoshimawa, M., Enyo, A., 1993. Stress relaxation function of bone and bone collagen. J. Biomech. 26 (12), 1369–1376. Iyo, T., Maki, Y., Sasaki, N., Nakata, M., 2004. Anisotropic viscoelastic properties of cortical bone. J. Biomech. 37 (9), 1433–1437. Rahulkumar, P., Jagota, A., Bennison, S., Saigal, S., 2000. Cohesive element modeling of viscoelastic fracture: Application to peel testing of polymers. Int. J. Solids. Struct. 37 (13), 1873–1897. Wineman, A.S., Rajagopal, K.R., 2000. Mechanical Response of Polymers: An Introduction. Cambridge university press. Behiri, J.C., Bonfield, W., 1980. Crack velocity dependence of longitudinal fracture in bone. J. Mater Sci. 15 (7), 1841–1849. Garner, E., Lakes, R., Lee, T., Swan, C., Brand, R., 2000. Viscoelastic dissipation in compact bone: Implications for stress-induced fluid flow in bone. J. Biomech. Eng-T ASME 122 (2), 166–172. Sahimi, M., 2003. Heterogeneous Materials II, Nonlinear and Breakdown Properties and Atomistic Modeling. Springer, New York. Batrouni, G.G., 1998. A. Hansen, Fracture in three-dimensional fuse networks. Phys. Rev. Lett. 80 (2), 325–328. Herrmann, H.J., Roux, S. (Eds.), 1990. Statistical Models for the Fracture of Disordered Materials. Elsevier, Amsterdam. Nukala, P.K.V.V., Simunovic, S., 2005. A continuous damage random thresholds model for simulating the fracture behavior of nacre. Biomaterials 26 (30), 6087–6098. Nukala, P.K.V.V., Simunovic, S., 2005. Statistical physics models for nacre fracture simulation. Phys. Rev. E 72 (4), 041919. de Arcangelis, L., Redner, S., Herrmann, H., 1985. A random fuse model for breaking processes. J. Phys. Lett-Paris 46 (13), 585–590. Kahng, B., Batrouni, G., Redner, S., de Arcangelis, L., Herrmann, H., 1988. Electrical breakdown in a fuse network with random, continuously distributed breaking strengths. Phys. Rev. B 37 (13), 7625–7637. Zapperi, S., Vespignani, A., Stanley, H., 1997. Plasticity and avalanche behaviour in microfracturing phenomena. Nature 388 (6643), 658–660. Hansen, A., 2005. Physics and fracture. Comput. Sci. Eng. 44 (1), 90–95. Kotha, S., Kotha, S., Guzelsu, N., 2000. A shear-lag model to account for interaction effects between inclusions in composites reinforced with rectangular platelets. Compos. Sci. Technol. 60 (11), 2147–2158.
J O U R N A L O F T H E M E C H A N I C A L B E H AV I O R O F B I O M E D I C A L M AT E R I A L S
Rho, J.-Y., Kuhn-Spearing, L., Zioupos, P., 1998. Mechanical properties and the hierarchical structure of bone. Med. Eng. Phys. 20 (2), 92–102. Delaplace, A., Pijaudier-Cabot, G., Roux, S., 1996. Progressive damage in discrete models and consequences on continuum modelling. J. Mech. Phys. Solids 44 (1), 99–136. Zocher, M., Groves, S., Allen, D., 1997. Three-dimensional finite element formulation for thermoviscoelastic orthotropic media. Int. J. Numer. Meth. Eng. 40 (12), 2267–2288. Patel, R.K., Bhattacharya, B., Basu, S., 2007. A finite element based investigation on obtaining high material damping over a large frequency range in viscoelastic composites. J. Sound. Vib. 303 (3–5), 753–766. Dai, Q., Sadd, M.H., 2004. Micromechanical modelling of viscoelastic behaviour of asphalt materials. In: 17th ASCE Eng. Mech. Conf., pp. 17–24. Anup, S., Sivakumar, S.M., Suraishkumar, G.K., Method for simulating fracture of non-linear materials by random fuse model, (in preparation). MathWorks, http://www.mathworks.com/(2000). Wang, X., Li, X., Bank, R., Agrawal, C., 2002. Effects of collagen unwinding and cleavage on the mechanical integrity of the collagen network in bone. Calcified Tissue Int. 71 (2), 186–192.
3 (2010) 260–267
267
Atkinson, T., Ewers, B., Haut, R., 1999. The tensile and stress relaxation responses of human patellar tendon varies with specimen cross-sectional area. J. Biomech. 32 (9), 907–914. Cusack, S., Miller, A., 1979. Determination of the elastic constants of collagen by brillouin light scattering. J. Mol. Biol. 135 (1), 29–51. Kotha, S., Guzelsu, N., 2007. Tensile behavior of cortical bone: Dependence of organic matrix material properties on bone mineral content. J. Biomech. 40 (1), 36–45. Schreurs, P.J., 2007. Material models 2007. URL http://www.mate. tue.nl/piet/edu/mcm/htm2/Pcmm3.html. Natali, A.N., Meroi, E.A., 1989. A review of the biomechanical properties of bone as a material. J. Biomed. Eng. 11 (4), 266–276. Melnis, A.E., Knets, I.V., 1982. Effect of the rate of deformation on the mechanical properties of compact bone tissue. Mech. Comp. Mat. 18, 358–362. Kotha, S.P., Guzelsu, N., 2002. Modeling the tensile mechanical behavior of bone along the longitudinal direction. J. Theoret. Bio. 219, 269–279. Ebenstein, D.M., Puritt, L.M., 2006. Nanoindentation of Biological Materials. Nanotoday 1 (3), 26–33.
3 (2010) 260–267
available at www.sciencedirect.com
journal homepage: www.elsevier.com/locate/jmbbm
Research paper
Influence of viscoelasticity of protein on the toughness of bone S. Anup a,∗ , S.M. Sivakumar a , G.K. Suraishkumar b a Department of Applied Mechanics, Indian Institute of Technology Madras, India b Department of Bio-technology, Indian Institute of Technology Madras, India
A R T I C L E
I N F O
A B S T R A C T
Article history:
Bone is an ultrafine composite of protein (collagen) and mineral (hydroxyapatite). An anal-
Received 23 May 2009
ysis to determine the influence of the viscoelasticity of protein on the toughness of bone at
Received in revised form
the ultrafine scale is conducted by developing a discrete lattice model appropriate for the
22 October 2009
ultrafine scale called the incremental continuous damage random fuse model (ICDRFM).
Accepted 24 October 2009
Collagen viscoelasticity at ultrafine scale is shown to contribute significantly to the tough-
Published online 4 November 2009
ness of bone. The results obtained are important in the design of biomimetic ultratough artificial composites. c 2009 Elsevier Ltd. All rights reserved.
1.
Introduction
Bio-composites such as bone, nacre and dentin possess excellent mechanical properties such as toughness (fracture energy), when compared to their constituents. Several reasons have been given by researchers as to why they exhibit such excellent properties. One of the explanations put forward is the presence of a hierarchical structure (Ji and Gao, 2004) with dimensions ranging from millimeters (macroscale) to the order of nm (ultrafine scale). The presence of sacrificial bonds in protein, leading to a saw tooth stress–strain curve has also been shown in the literature to be a mechanism which enhances the strength and toughness of bio-composites (Ji and Gao, 2004; Smith et al., 1999; Thompson et al., 2001). These studies point to the importance to be given to the investigations on the reasons for toughness of bio-composites. The information obtained is useful in the development of artificial composites as well as in selecting appropriate drugs to cure
diseases (in the case of bone) (Nalla et al., 2006; Anup et al., 2008). Among the various bio-composites, bone has been extensively studied due to its clinical importance. A more extensive bibliography on these studies has been given by Nalla et al. (2004). The ultrafine structure of bone is expected to play an important role in its overall mechanical behavior (Gao et al., 2003; Ji and Gao, 2004). Therefore, the scope of the study is limited to ultrafine structure level (of the order of 10−7 m to 10−8 m), where a network of collagen molecules inter-spaced with mineral crystals is considered (Fratzl et al., 2004; Ji and Gao, 2004). Platelets are made up of hydroxyapatite. Collagen molecules have a triple helix structure and each of these is made up of polypeptides. Collagen molecules contain sacrificial bonds which break instead of the back bone, and this contributes to toughness of bone (Thompson et al., 2001). The idealized arrangement assumed at the ultrafine level consisting of a staggered structure of mineral platelets embedded
∗ Corresponding address: Indian Institute of Space Science and Technology, Department of Aerospace Engineering, Thiruvananthapuram 695022, India. Tel.: +91 9847395059. E-mail addresses: [email protected], [email protected] (S. Anup). c 2009 Elsevier Ltd. All rights reserved. 1751-6161/$ - see front matter doi:10.1016/j.jmbbm.2009.10.007
J O U R N A L O F T H E M E C H A N I C A L B E H AV I O R O F B I O M E D I C A L M AT E R I A L S
a
2h_1
3 (2010) 260–267
261
b
PLATELET
PLATELET TENSION ELEMENTS
MATRIX
2l1
MATRIX SHEAR ELEMENTS
q
2l2
MATRIX TENSION ELEMENTS
h3
Fig. 1 – (a) The idealized arrangement at the ultrafine scale. (b) A part of the random fuse lattice used for simulation.
within a protein matrix is shown in Fig. 1 (a). The present study focuses on ultrafine scale, though the toughness of bone can be influenced by mechanisms, operating at other scales. However, the processes taking pace at ultrafine scale have a profound influence on the way bone behaves at the macroscale. Bio-composites like bone exhibit significant viscoelasticity (Sasaki et al., 1993; Iyo et al., 2004; Ji and Gao, 2004). Viscoelasticity can dissipate energy in materials during deformation and fracture (Rahulkumar et al., 2000), thereby providing higher toughness. Among the constituents of bone, collagen is found to be viscoelastic from experiments while mineral (hydroxyapatite) is found to be nearly linear elastic (Ji and Gao, 2004). One of the important ways in which the effect of viscoelasticity can be observed is the variation in the stress–strain response with strain rate (Wineman and Rajagopal, 2000). When failure of the material occurs, failure energy dissipation (toughness) varies with strain rate (Behiri and Bonfield, 1980) due to viscoelasticity. The viscoelastic energy dissipation in bone have been attributed to various mechanisms including viscous like cement line motion in osteon interfaces, thermoelastic damping, damping in lamallae and fluid flow in bone (Garner et al., 2000). Scant open literature is available on the relative contributions of these mechanisms in viscoelastic energy dissipation characteristics (toughness) in bone during failure. Collagen, being viscoelastic is a probable cause of such toughness in bone. Therefore, investigations into the role played by collagen in the viscoelastic failure dissipation in bone at the ultrafine scale is considered be an important research direction. Discrete lattice models have been proposed for simplified numerical as well as analytical studies of failure of materials with disorder (Sahimi, 2003; Batrouni, 1998; Herrmann and Roux, 1990). In discrete lattice models, the medium investigated is divided into a discrete set of elastic bonds. These elastic bonds fail at random threshold values (Nukala and Simunovic, 2005). The main features that characterize discrete lattice models are disorder (randomness), stress–strain characteristics and a rule by which each of the bonds (elements) of the lattice are broken. The stress–strain and
breaking response characteristic of each element corresponds to the response at the microstructural level in the material (Nukala and Simunovic, 2005). Discrete lattice models are expected to capture the essential physics of the problem including mechanisms that influence the toughness such as crack deviation, branching and stress concentration effects at crack tips (Nukala and Simunovic, 2005). A scalar form of the discrete lattice model, called the random fuse model can be used to simulate failure, where a voltage is applied to a network of resistors (de Arcangelis et al., 1985; Kahng et al., 1988). A resistor network is a scalar analog of the elastic medium. The analysis is simpler in such a model. However, the important characteristics of the problem are not lost due to this simplification (Nukala and Simunovic, 2005). A modified version of the random fuse model, the continuous damage random threshold fuse model (CDRFM), has also been used to simulate the failure in materials if damage evolution takes place (Zapperi et al., 1997). CDRFM is extremely useful for comparing the toughness of bio-composites such as nacre as demonstrated by Nukala and Simunovic (2005). However, almost all analysis involving continuous damage random fuse models are used for materials with linear constitutive laws only as could be seen from literature (Nukala and Simunovic, 2005; Hansen, 2005). In this context, a computationally inexpensive and easily implementable, but useful method would be to simulate failure under quasistatic loads employing a random fuse model, while incorporating viscoelastic nature of the constituents. Therefore, an incremental continuous damage random fuse model (ICDRFM) is developed in the present work to include viscoelastic constitutive law and is used to investigate the influence of collagen viscoelasticity on the toughness of bone at the ultrafine scale. Therefore, the objectives of the present study are to • Develop an incremental continuous damage random fuse model to simulate the failure of viscoelastic materials such as bone. • Investigate the influence of collagen viscoelasticity at the ultrafine scale on the toughness of bone.
262 2.
J O U R N A L O F T H E M E C H A N I C A L B E H AV I O R O F B I O M E D I C A L M AT E R I A L S
Formulation of the incremental CDRFM
The various dimensions used in the study are as given in Table 1 (Kotha et al., 2000; Rho et al., 1998). The dimensions are normalized with respect to the thickness of the platelets (2h_1). A unit normalized width, w is assumed in the out of plane direction. The mineral platelets are assumed to carry tensile stresses only and matrix is assumed to transfer shear stresses between platelets (Ji and Gao, 2004; Kotha et al., 2000). The shear stresses are assumed to be a constant, while the tensile stresses are assumed to change along the length of the platelets which are consistent with observations made by other researchers (Ji and Gao, 2004; Kotha et al., 2000). The random fuse models can take care of the redistribution of stresses due to the presence of platelets automatically (Nukala and Simunovic, 2005). The idealized staggered arrangement of mineral platelets embedded in a protein matrix is modeled as a square lattice network of L × L size. A part of the lattice network used for simulation is given in Fig. 1(b). A current–voltage analog of the model is employed (Nukala and Simunovic, 2005). This is done since the stress–strain relation for an elastic medium is similar to the current voltage relationship of a resistor network. Moreover, in the case of an elastic medium, the total force in any direction at any node is equal to zero (Force equilibrium). Similarly, in the case of an electrical network, sum of the electric current at any node is equal to zero (Kirchhoff’s Law). The matrix shear elements (shear fuses) connect the longer sides of mineral platelets while matrix tension elements (tension fuses) connect the shorter sides. Mineral platelets are modeled as mineral tension elements (tension fuses). The bond breaking thresholds (strengths) for both of these are prescribed to be constant. The conductance of an element is taken to be its stiffness. The platelet is assumed to be linear elastic and brittle, represented by tension elements. The conductance of a platelet tension element is made zero once the current passing through it exceeds the critical threshold value. The protein (collagen) in bone has been shown experimentally to contain sacrificial bonds that protect the collagen backbone and break one after the other in response to stress (Thompson et al., 2001). Matrix is assumed to be a material with evolving damage due to breaking of sacrificial bonds and is represented by shear elements with conductance decreasing in discrete steps (Nukala and Simunovic, 2005). Each time the current (stress) in the matrix shear fuse exceeds the critical threshold value, it is considered to have failed, the conductance is decreased by multiplying by a factor a where, 0 < a < 1, and this is termed matrix damage. The conductance of a matrix element is calculated to be ai Cm , where i is the number of times the matrix element has damaged. Once the maximum shear strain is reached, the conductance is made zero to simulate complete failure of the matrix and is termed as matrix tearing. Thus matrix failure consists of both matrix damage and matrix tearing. This stress–strain behavior of the matrix and the lattice network used in the present study is similar to the CDRFM model that has been employed by Nukala and Simunovic (2005) for the simulation of failure of nacre. In traditional CDRFM, numerically, a unit voltage difference is applied and current flowing through all the elements is found out. For each element, the ratio of the current
3 (2010) 260–267
flowing through each of the elements to the threshold value of each element is found out. The element having the maximum value of this ratio, rmax is considered to have failed (Delaplace et al., 1996). It is assumed that the time taken for current relaxation is much higher than the time taken for bond breakage. Hence, only one bond is broken at a time. Since the equations are linear, the external voltage and current are multiplied by 1/rmax , so that the current in the broken element is equal to its threshold value. This means that the externally applied voltage (displacement) is either increased or decreased so the stress in one of the bonds is just enough to break it. Once the breakage of a bond occurs, the stiffness matrix is altered and the equations are solved for the new system with the changed stiffness matrix. However, in the case of incremental CDRFM, which is developed in the present study, an incremental voltage difference is applied between the bus bars and the Kirchhoff’s equations are solved to find out the current (stress) in each of the elements of the lattice. If the stresses in none of the matrix elements exceeds the critical value, the voltage difference is again increased incrementally. This process is repeated until the current (stress) in at least one of the elements exceeds the critical threshold value. If the number of elements having a current above the threshold value is more than one, then the element having highest value of current is assumed to have failed. In that case, the stiffness of the particular failed element is reduced by a factor and the process is continued. This means that the voltage is always increased in incremental steps and is similar to a displacement controlled test, with a constant displacement rate. This process of breaking of bonds one by one goes on until the lattice separates completely into two parts. Every time, the equations are solved, the current passing through the system (Ic ) and voltage (Vc ) applied are stored. By plotting Ic /L vs Vc /L, the stress–strain graph is obtained. The toughness of the composite (bone) is found out from the simulation, in the sense of work of fracture, which is the area under the stress–strain curve multiplied by the length of the specimen. In the present work, the matrix material alone is assumed to fail since this mode of failure would provide maximum overall toughness (Anup et al., 2008). In addition to viscoelasticity, matrix is assumed to have evolving damage. This coupled behavior is approximated by the following steps. Initially, the matrix element is assumed to be elastic with evolving damage and the number of damages required to reach the maximum strain is found out, as given in literature (Nukala and Simunovic, 2005). In the next case, matrix is considered to be viscoelastic. The number of damages that a matrix element undergoes for the viscoelastic case is assumed to be the same as that found for the linear elastic case. This assumption is based on physical grounds, since the number of sacrificial bonds remains the same whether the material is considered to be elastic or viscoelastic. Therefore, when material property of the matrix is changed, the number of sacrificial bonds is kept a constant. Matrix is also assumed to be linear viscoelastic and hence is modeled as a standard solid with a single relaxation time. The response of viscoelastic materials subjected to load has been studied by making use of the finite element method. Zocher et al. (1997) proposed a methodology to simulate this viscoelastic response, which has been used to model the
J O U R N A L O F T H E M E C H A N I C A L B E H AV I O R O F B I O M E D I C A L M AT E R I A L S
263
3 (2010) 260–267
Table 1 – Dimensions at ultrafine scale: lengths in nm, area in nm2 (Kotha et al., 2000; Rho et al., 1998). Platelet area (2h1 )
Dimension Value Normalized value Number of elements
2.5 1 –
Gap (2l1 ) 15 6 6
response of viscoelastic materials as well as composites (Patel et al., 2007). This methodology has been modified by Dai and Sadd (2004) for simulating damage in viscoelastic materials. This modified form (Dai and Sadd, 2004), is adapted in the present work to simulate the viscoelastic response of the matrix elements, and is given in detail elsewhere (Anup et al., in preparation). The basics of the methodology is given below. For viscoelastic materials, the shear stress τ(t) for current time t is expressed as a hereditary integral (Wineman and Rajagopal, 2000; Zocher et al., 1997; Schreurs, 2007) as shown in the following equation. Z t τ(t) = G(t − s)γ(s)ds, ˙ (1) 0
where G(t) is the relaxation function of a matrix element at current time t and γ˙ is the strain rate. Putting G(t) = G(t)/h3 , and displacement rate u˙ = γh ˙ 3 for a matrix element, the Eq. (1) becomes Z t ˙ τ(t) = G(t − s)u(s)ds, (2) 0
where G(t) is the relaxation function of the shear stiffness (conductance). Here, G(t) is expressed as a Prony series as given by G(t) = Gα + Gd e−t/TR ,
(3)
where Gα is the shear stiffness corresponding to the long term elastic modulus, Gd is the term corresponding to the difference between short term and long term shear stiffness and TR is the relaxation time. The total time interval [0, t] is discretized into intervals (time steps). It is assumed that the intervals are of equal length. The hereditary integral is split as an integral over [0, tn ] and an integral over the current or last increment [tn , tn+1 = t] (Schreurs, 2007). The subscript for t indicates the number of the time step. Substituting Eq. (3) into Eq. (2) and simplifying, we get Z tn Z t t−τ − t−τ ˙ ˙ τ(t) = Gα u(t) + Gd e TR u(s)ds + Gd e TR u(s)ds. (4) tn
0
Therefore, the stress is obtained as − ∆t T
τ(t) = Gα u(t) + e where Pd =
TR ∆t
R
τd (tn ) + Gd Pd ∆u, !
1−e
− ∆t T R
.
(5)
Here, u is the total displacement. The stress increment ∆τ is obtained by subtraction as given below
Platelet length (2l2 )
Matrix thickness (h3 )
Matrix area (q)
50 20 20
1.35 0.54 1
2.5 1 –
relation, − ∆t T
τd (tn ) = Gdc Pd ∆u + τd (tn−1 )e
R
(8)
where Gdc is the current value of Gd , which decreases with the number of times the matrix element fails (damages). For maintaining equilibrium, the following procedure is adopted. The global current (force) is obtained through an appropriate assembly of element contributions as given by the following equation. [Ks ] {∆u} = {χ}
(9)
where [Ks ] is the global tangent stiffness (conductance) matrix and {χ} is the global force. The Eq. (9) is solved to obtain the nodal voltages. The displacement in each element is the voltage difference across an element. The procedure employed for the simulation is outlined below. 1. Form the tangent stiffness matrix Ks employing the value of Ts (See Eq. (7)). For iteration = 0. 2. Find out the displacements (voltage), ∆u1 using Eq. (9). For iteration = n + 1. 3. ∆un+1 , τ(tn ) and τd (tn ) are obtained from previous time step. 4. Calculate ∆τ employing Eq. (6). 5. Calculate τ(t) = τ(tn ) + ∆τ. 6. Check if any matrix element exceeds the threshold value; then if the number of total failures (damages) has been reached, stress in that element is made zero else reduce the stress in the failed element by multiplying by a factor a. Update Gd and Gα for the failed element. Change the tangent stiffness matrix incorporating change in stiffness of the failed element. 7. Calculate τd (tn+1 ) employing Eq. (8). 8. Find out the incremental displacements ∆un+2 using Eq. (9) 9. Go to the next iteration. Periodic boundary conditions are imposed in the horizontal direction to simulate an infinite system. This is a standard approach of RVE to obtain material continuity. An incrementally increasing voltage difference is applied between the top and the bottom bus bars. For solving the Kirchhoff’s equations obtained as a result, the Matlab package (MathWorks, 2000) is employed in this work.
∆τ = τ(tn+1 ) − τ(tn ), ∆τ = [Gα + Gd Pd ]∆u − (1 − e
− ∆t T R
)τd (tn ).
(6)
Tangent stiffness (conductance) Ts is given as Ts = [Gα + Gd Pd ].
(7)
The total stress is obtained by adding up the stress increments in each time step. Here τd (tn ) is given by the recursive
3.
Results and discussion
Experiments show that the maximum tensile strain that collagen can undergo is about 10% (Wang et al., 2002). In this work, without much loss of generality, the maximum shear strain, εmax that matrix can undergo is also taken to be
J O U R N A L O F T H E M E C H A N I C A L B E H AV I O R O F B I O M E D I C A L M AT E R I A L S
fore stiffness becomes zero, nmax =
log(ε0 /εmax ) ≈ 33, where log(a)
ε0 = h3 Ym /G (Nukala and Simunovic, 2005). The value of a is chosen in the present work to be 0.9 as an example, though other values of a does not seem to affect the results (Nukala and Simunovic, 2005). The matrix tension elements are assumed to carry no stresses and are given zero conductance values. The basic value of parameters relating to viscoelastic behavior of the matrix was obtained from the stress relaxation experiments conducted (Atkinson et al., 1999) on human patellar tendon. The values obtained from tendon could be assumed to be that of collagen since patellar tendon consists mostly of collagen. A linear viscoelastic standard solid (a Prony series with a single relaxation time) is used in the present work to fit the data and the parameters were found out. The values got from the fitted curve are, relaxation time, TR = 20, and the ratio of short term modulus to long term modulus (enhancement factor), λ = 4. The short term modulus is taken to be the shear modulus of collagen available in literature (3.3 GPa) (Cusack and Miller, 1979). This gives the values of long term modulus, Gα = 0.825 GPa. There is no data available for shear strength of collagen, and there is a wide variation in the tensile strength obtained by various researchers, a typical value being 20 MPa (Ji and Gao, 2004; Kotha and Guzelsu, 2007). Hence, this value is taken as the threshold matrix shear strength in this study. Young’s modulus of the platelet is given a value of 100 GPa (Ji and Gao, 2004). A system size of L = 16 is used for analysis. The stress–strain graphs obtained using ICDRFM and the linear method are compared as shown in Figs. 2 and 3. For a viscoelastic material, when relaxation time is very small, i.e. when t → 0, the second term in Eq. (3) becomes equal to Gd , and G(t) = Gα + Gd which is the long term modulus of the material. Therefore, the predictions of the viscoelastic model with t → 0 should agree with the predictions of a model with linear elastic stress–strain law having a modulus equal to the long term modulus of the viscoelastic material. Similarly, for large relaxation time, predictions of the viscoelastic model should agree with the predictions of a model with linear elastic stress–strain law having a modulus equal to the short term modulus of the viscoelastic material. ICDRFM developed in this study is compared with the predictions of a model with linear elastic stress–strain law and is found to agree for two extreme cases of small and large relaxation times. The fact that fracture energy (area under the stress–strain curves) obtained for these two cases for the viscoelastic model agrees with that obtained for corresponding model with linear elastic modulus is also consistent with the fact that viscoelastic energy dissipation is found to tend to zero for very high and small relaxation times (Wineman and Rajagopal, 2000). In Fig. 2, the stress–strain graph obtained by the linear method (shear modulus, Gm = 3.3 GPa) and that obtained by ICDRFM (material being viscoelastic, TR = 30 000 s) are compared and seen to be in good agreement. In Fig. 3, the stress–strain graph obtained by the linear method (shear
3 (2010) 260–267
0.1
0.08 STRESS (GPa)
10% (Nukala and Simunovic, 2005). Each time the matrix damages, stiffness is decreased to 90% of its previous value, to simulate the breaking of sacrificial bonds. This is accomplished by putting value of a = 0.9. This gives the number of damages that each matrix shear element undergoes be-
0.06
0.04
DIRECT METHOD (Nukala et al., 2005)
0.02
INCREMENTAL METHOD
0
0
1
2
3
4 5 STRAIN
6
7
8 × 10–3
Fig. 2 – Comparison of stress–strain graphs for linear elastic case and viscoelastic case with very large relaxation time. The toughness (area under the curve) are in agreement for the two cases.
0.1
0.08 STRESS (GPa)
264
DIRECT METHOD (Nukala et al., 2005)
0.06
INCREMENTAL METHOD 0.04
0.02
0
0
0.005
0.01
0.015 STRAIN
0.02
0.025
Fig. 3 – Comparison of stress–strain graphs for linear elastic case and viscoelastic case with very small relaxation time. The toughness (area under the curve) are in agreement for the two cases.
modulus, Gm = 0.825 GPa) and that obtained by ICDRFM (material being viscoelastic, TR = 1/10 s) are compared and seen to be in good agreement. Fig. 4 shows the variation of fracture energy dissipation with relaxation time. From Fig. 4, it could be seen that the fracture energy dissipation (area under the stress–strain curve) increases with decrease in relaxation time, reaches a maximum and then decreases. This trend is consistent with the case during deformation of other viscoelastic materials (Wineman and Rajagopal, 2000). Therefore, these results indicate that ICDRFM developed in this work is a valid method for simulating failure, and for comparison of toughness (Fracture energy) values. Fig. 5 shows the effect of strain rate on fracture energy dissipation. The viscous component of fracture energy dissipation is almost zero below a strain rate of 6.25 × 10−6 s−1 . Increase in strain rate increases the fracture energy dissipation until a strain rate of 6.25 × 10−4 s−1 . Thereafter, with increase in strain rate, fracture energy dissipation (toughness) decreases considerably. Therefore, the toughness of bone at ultrafine scale shows an increasing–decreasing trend with increase in strain rate.
J O U R N A L O F T H E M E C H A N I C A L B E H AV I O R O F B I O M E D I C A L M AT E R I A L S
3
0.1
2.5
0.08
2 0.06
W / W0
STRESS (GPa)
265
3 (2010) 260–267
0.04
1.5 1
0.02
0
0
EXPERIMENT (Behiri et al.,1980)
0.5
0.005 0.01
0.015 0.02 0.025 0.03 STRAIN
0.035 0.04
Fig. 4 – Stress–strain curves for different values of relaxation time. Note that area under the stress–strain curve increases, reaches a maximum value and then decreases.
0
SIMULATION 0
0.2
0.4 0.6 STRAIN RATE ( s–1)
0.8
1 × 10–3
Fig. 6 – The variation of work of fracture with strain rate is shown in the figure. The work of fracture increases with strain rate, reaches a maximum and then decreases. Especially note the strain rate at which maximum work of fracture occurs, which is close together for the two curves.
0.1
STRESS (GPa)
0.08
0.06
STRAIN RATE 6.25 ×10–6 s–1 STRAIN RATE 6.25 ×10–5 s–1
0.04
STRAIN RATE 3.12 ×10–4 s–1 STRAIN RATE 6.25 ×10–4 s–1
0.02
0
0
0.005 0.01
0.015 0.02 0.025 0.03 STRAIN
0.035 0.04
Fig. 5 – Stress–strain curves for different values of strain rate. Note that area under the stress–strain curve increases, reaches a maximum value and then decreases.
The variation in work of fracture (fracture energy dissipation) in the case of bovine bone for different cross head speeds was given by Behiri and Bonfield (1980). Taking the distance between loading platens in their work to be approximately 20 mm, the strain rates can be calculated to vary from 8.5 × 10−6 s−1 to 8.5 × 10−4 s−1 . Increase in strain rate to a value of 8.5 × 10−4 s−1 , produced a considerable decrease in the work of fracture, according to their experimental results (Behiri and Bonfield, 1980). The results of the present study are compared with that of Behiri and Bonfield (1980) in Fig. 6. Here, W is the work of fracture obtained as the area under the stress–strain curve and W0 is the value of the work of fracture for the smallest strain rate (6.25×10−6 s−1 ). The simulated results of the present study and experimental results obtained from literature (Behiri and Bonfield, 1980) follow the same trend, with a good match over the strain rate at which the peak energy dissipation is obtained. Natali and Meroi (1989) conducted experiments on bone under a larger range of strain rates (from 10−5 s−1 to 1 s−1 ). The maximum values of work dissipated in their study is seen to be in the 10−3 s−1 to 10−2 s−1 range which is the same strain rate range over which maximum viscoelastic dissipation is observed in the simulation in the present
study. Melnis and Knets (1982) also have data relating to the influence of strain rate on the specific deformation energy over a range of moisture contents. For smaller moisture contents (2.5%), a maximum value of the deformation energy is reached at a strain rate of 10−3 s−1 which is similar to that obtained in the simulation in the present study. These results suggest that collagen viscoelasticity at ultrafine scale contributes significantly to the energy dissipation during fracture in bone. The viscoelastic dissipation in bone is attributed to various mechanisms such as viscous like cement line motion in osteon interfaces, thermoelastic damping, damping in lamallae and fluid flow in bone (Garner et al., 2000). The results obtained from the present study is significant since it bring forth evidence in support of collagen viscoelasticity at the ultrafine scale amongst the various mechanisms suggested (Garner et al., 2000), as the important contributor to viscoelastic fracture energy dissipation (toughness) in bone. This finding is significant since artificial composites could be modeled based on the structure and properties of bone at ultrafine scale (staggered arrangement of the mineral platelets embedded in a viscoelastic matrix) so that maximum energy dissipation and hence toughness is achieved. The model employed in the present work is an over simplified version of the processes taking place at the ultrafine structure in bone. The viscoelastic properties of collagen is modeled using a standard linear solid rather than nonlinear viscoelastic material models that are more appropriate for other such biological materials. Nevertheless, this simple model gives insights into the way in which viscoelasticity affects fracture energy dissipation. The requirement of matrix mode of failure so as to provide maximum toughness has been described earlier (Anup et al., 2008). There seems to be an additional requirement for matrix mode of failure in bone at the ultrafine scale so as to make use of the viscoelasticity in the fracture energy dissipation. The structural arrangement at the ultrafine scale, while helping to maintain a stiffness close to that of the mineral (Ji and Gao, 2004), aids in higher fracture energy dissipation, by permitting matrix mode of failure alone. The model takes into account the
266
J O U R N A L O F T H E M E C H A N I C A L B E H AV I O R O F B I O M E D I C A L M AT E R I A L S
mechanisms acting at the ultrafine scale only. However, the effect of mechanisms acting on other scales is seen to be not appreciable given the fact that the macroscopic experimental and simulated curves follow the same trend. Moreover, many researchers have put forward the idea that the mechanisms at the ultrafine scale predominantly affect the stress–strain response at the macroscale (Kotha and Guzelsu, 2002). Since ultrafine scale properties are being considered in this study, use of information obtained from nanoindentation experiments at various loading rates could have been more appropriate. The data related to nanoindentation on bone is extensive. However, the studies have been mostly conducted with respect to how the properties vary with respect to different tissues, composition and site rather than loading rate (Ebenstein and Puritt, 2006). However, there is scant open literature on the effect on loading rate on the toughness of bone. It should be noted that apart from the modes of failure that has been assumed in the present study (platelet breakage and matrix failure), there could be more failure modes like debonding at the interface between the matrix and the mineral platelet. The platelets and matrix were assumed to be isotropic, though in reality they may possess anisotropic properties. Research on the effect of other modes of failure and anisotropy on the toughness of bone is being pursued currently.
4.
Conclusions
In this work, an incremental continuous damage random fuse model is developed for use at the ultrafine scale in bone and is employed to examine the way in which viscoelasticity of protein influences the toughness of bone at ultrafine scale. The method developed is significant in the fact that the viscoelastic properties of the matrix could also be included in the model, while retaining a simple and computationally inexpensive numerical analysis. The toughness of bone at ultrafine scale shows an increasing–decreasing trend with increase in strain rate. There is a very good correlation between experimental results obtained from the literature and the simulation results of the present study, especially on the strain rate at which maximum toughness is achieved. Therefore, the results obtained from the present study bring forth evidence in support of the collagen viscoelasticity as one of the important contributors to the toughness (fracture energy) in bone. The results are significant in the design of biomimetic composites.
Acknowledgement SA gratefully acknowledges the financial support under the QIP scheme of the AICTE, New Delhi, India. REFERENCES
Ji, B., Gao, H., 2004. Mechanical properties of nanostructure of biological materials. J. Mech. Phys. Solids 52 (9), 1963–1990.
3 (2010) 260–267
Smith, B., Schaffer, T., Viani, M., Thompson, J., Frederick, N., Kindt, J., Belcher, , Stucky, G., Morse, D., Hansma, P., 1999. Molecular mechanistic origin of the toughness of natural adhesives, fibres and composites. Nature 399 (6738), 761–763. Thompson, J., Kindt, J., Drake, B., Hansma, H., Morse, D., Hansma, P., 2001. Bone indentation recovery time correlates with bond reforming time. Nature 414 (6865), 773–776. Nalla, R.K., Kruzic, J.J., Kinney, J.H., Balooch, M., Ager III, J.W., Ritchie, R.O., 2006. Role of microstructure in the aging-related deterioration of the toughness of human cortical bone. Mat. Sci. Eng. C-Bio. S 26 (8), 1251–1260. Anup, S., Sivakumar, S.M., Suraishkumar, G.K., 2008. Influence of relative strength of constituents on the overall strength and toughness of bone. J. Mech. Med. Biol. 8 (4), 527–539. Nalla, R.K., Kruzic, J.J., Ritchie, R.O., 2004. On the origin of the toughness of mineralized tissue: Microcracking or crack bridging?. Bone 34 (5), 790–798. Gao, H., Ji, B., Jager, I.L., Arzt, E., Fratzl, P., 2003. From the Cover: Materials become insensitive to flaws at nanoscale: Lessons from nature. PP Natl. Acad. Sci. USA 100 (10), 5597–5600. Fratzl, P., Gupta, H.S., Paschalis, E.P., Roschger, P., 2004. Structure and mechanical quality of the collagen-mineral nanocomposite in bone. J. Mater Chem. 14 (14), 2115–2123. Sasaki, N., Nakayama, Y., Yoshimawa, M., Enyo, A., 1993. Stress relaxation function of bone and bone collagen. J. Biomech. 26 (12), 1369–1376. Iyo, T., Maki, Y., Sasaki, N., Nakata, M., 2004. Anisotropic viscoelastic properties of cortical bone. J. Biomech. 37 (9), 1433–1437. Rahulkumar, P., Jagota, A., Bennison, S., Saigal, S., 2000. Cohesive element modeling of viscoelastic fracture: Application to peel testing of polymers. Int. J. Solids. Struct. 37 (13), 1873–1897. Wineman, A.S., Rajagopal, K.R., 2000. Mechanical Response of Polymers: An Introduction. Cambridge university press. Behiri, J.C., Bonfield, W., 1980. Crack velocity dependence of longitudinal fracture in bone. J. Mater Sci. 15 (7), 1841–1849. Garner, E., Lakes, R., Lee, T., Swan, C., Brand, R., 2000. Viscoelastic dissipation in compact bone: Implications for stress-induced fluid flow in bone. J. Biomech. Eng-T ASME 122 (2), 166–172. Sahimi, M., 2003. Heterogeneous Materials II, Nonlinear and Breakdown Properties and Atomistic Modeling. Springer, New York. Batrouni, G.G., 1998. A. Hansen, Fracture in three-dimensional fuse networks. Phys. Rev. Lett. 80 (2), 325–328. Herrmann, H.J., Roux, S. (Eds.), 1990. Statistical Models for the Fracture of Disordered Materials. Elsevier, Amsterdam. Nukala, P.K.V.V., Simunovic, S., 2005. A continuous damage random thresholds model for simulating the fracture behavior of nacre. Biomaterials 26 (30), 6087–6098. Nukala, P.K.V.V., Simunovic, S., 2005. Statistical physics models for nacre fracture simulation. Phys. Rev. E 72 (4), 041919. de Arcangelis, L., Redner, S., Herrmann, H., 1985. A random fuse model for breaking processes. J. Phys. Lett-Paris 46 (13), 585–590. Kahng, B., Batrouni, G., Redner, S., de Arcangelis, L., Herrmann, H., 1988. Electrical breakdown in a fuse network with random, continuously distributed breaking strengths. Phys. Rev. B 37 (13), 7625–7637. Zapperi, S., Vespignani, A., Stanley, H., 1997. Plasticity and avalanche behaviour in microfracturing phenomena. Nature 388 (6643), 658–660. Hansen, A., 2005. Physics and fracture. Comput. Sci. Eng. 44 (1), 90–95. Kotha, S., Kotha, S., Guzelsu, N., 2000. A shear-lag model to account for interaction effects between inclusions in composites reinforced with rectangular platelets. Compos. Sci. Technol. 60 (11), 2147–2158.
J O U R N A L O F T H E M E C H A N I C A L B E H AV I O R O F B I O M E D I C A L M AT E R I A L S
Rho, J.-Y., Kuhn-Spearing, L., Zioupos, P., 1998. Mechanical properties and the hierarchical structure of bone. Med. Eng. Phys. 20 (2), 92–102. Delaplace, A., Pijaudier-Cabot, G., Roux, S., 1996. Progressive damage in discrete models and consequences on continuum modelling. J. Mech. Phys. Solids 44 (1), 99–136. Zocher, M., Groves, S., Allen, D., 1997. Three-dimensional finite element formulation for thermoviscoelastic orthotropic media. Int. J. Numer. Meth. Eng. 40 (12), 2267–2288. Patel, R.K., Bhattacharya, B., Basu, S., 2007. A finite element based investigation on obtaining high material damping over a large frequency range in viscoelastic composites. J. Sound. Vib. 303 (3–5), 753–766. Dai, Q., Sadd, M.H., 2004. Micromechanical modelling of viscoelastic behaviour of asphalt materials. In: 17th ASCE Eng. Mech. Conf., pp. 17–24. Anup, S., Sivakumar, S.M., Suraishkumar, G.K., Method for simulating fracture of non-linear materials by random fuse model, (in preparation). MathWorks, http://www.mathworks.com/(2000). Wang, X., Li, X., Bank, R., Agrawal, C., 2002. Effects of collagen unwinding and cleavage on the mechanical integrity of the collagen network in bone. Calcified Tissue Int. 71 (2), 186–192.
3 (2010) 260–267
267
Atkinson, T., Ewers, B., Haut, R., 1999. The tensile and stress relaxation responses of human patellar tendon varies with specimen cross-sectional area. J. Biomech. 32 (9), 907–914. Cusack, S., Miller, A., 1979. Determination of the elastic constants of collagen by brillouin light scattering. J. Mol. Biol. 135 (1), 29–51. Kotha, S., Guzelsu, N., 2007. Tensile behavior of cortical bone: Dependence of organic matrix material properties on bone mineral content. J. Biomech. 40 (1), 36–45. Schreurs, P.J., 2007. Material models 2007. URL http://www.mate. tue.nl/piet/edu/mcm/htm2/Pcmm3.html. Natali, A.N., Meroi, E.A., 1989. A review of the biomechanical properties of bone as a material. J. Biomed. Eng. 11 (4), 266–276. Melnis, A.E., Knets, I.V., 1982. Effect of the rate of deformation on the mechanical properties of compact bone tissue. Mech. Comp. Mat. 18, 358–362. Kotha, S.P., Guzelsu, N., 2002. Modeling the tensile mechanical behavior of bone along the longitudinal direction. J. Theoret. Bio. 219, 269–279. Ebenstein, D.M., Puritt, L.M., 2006. Nanoindentation of Biological Materials. Nanotoday 1 (3), 26–33.