A preliminary investigation of the dynamic viscoelastic relaxation of bovine cortical bone

A new experimental approach is proposed to characterize the dynamic viscoelastic relaxation behaviour of cortical bone. Theoretical models are presented to show that a linear viscoelastic material, when allowed to relax between two long elastic bars, will produce stress, strain and strain rate histories that contain characteristic features. Furthermore, typical experimental results are presented to show that these characteristic features are observed during split Hopkinson bar tests on bovine cortical bone using a Cone-in-Tube striker. The interpretation of this behaviour in the context of a standard linear viscoelastic model is discussed.


Introduction
The design of protective structures to prevent injury during impact loading events, such as sporting accidents and vehicle collisions, requires a detailed understanding of the mechanisms of dynamic bone deformation and fracture.Traditional approaches to quantify responses of the human body to impacts have typically included full scale experimental testing involving the use of a post mortem human subject (PMHS) [14].However, such tests are expensive and time consuming which has lead to the increasing use of numerical modelling [6], for which detailed material models of bone are required.
At the macroscopic level, bone is classified into cortical (hard/compact) bone and cancellous/trabecular (soft/spongy) bone, of which the former is the primary skeletal load-bearing component.At a microscopic level, bone is a complex composite material composed of mineral crystals (hydroxyapatite) embedded in an organic matrix (primarily collagen fibres).This combination of organic and mineral components gives rise to the viscoelastic quasi-brittle behaviour of bone, i.e. bone displays highly rate sensitive viscoelastic behaviour but fails via progressive brittle fracture.Additional viscous contributions may arise during dynamic loading events due to the fluid component contained within the bone cavities.
Given its viscoelastic nature, the development of a material model for cortical bone requires knowledge of its mechanical properties over a wide range of strain rates, from quasi-static (QS) to high strain rate (HSR), i.e. 10 −3 s −1 − 10 3 s −1 .Johnson et al. [8] reported an extensive survey of the rate dependence of the apparent Young's modulus of cortical bone.They note that most of the work has focused on strain rates in the range of 10 −3 s −1 − 0.1 s −1 , which are relevant to activities such as walking and running [2].By contrast, high strain rate bone characterization tests, using drop testing [10] or split Hopkinson bar (SHB) [1,3] techniques, are typically restricted to strain rates above 300 s −1 , which correspond to severe impact loading conditions.Consequently, the data showed that the intermediate strain rate (ISR) range (i.e. 1 s −1 − 100 s −1 ) is sparsely populated.Furthermore, despite the relatively large quantity of HSR data available, there are some concerns regarding the quality of the older data since it is only within the last decade that techniques for conducting dynamic compression tests on bone at nearconstant strain rates have been developed [1,3].
The issue of the lack of data in the ISR regime has recently been addressed by Cloete et al. [4], who presented results for the dynamic compression of bovine cortical bone.A novel Cone-in-Tube (CiT) striker design was introduced that allows constant strain rate compression tests to be conducted in the upper ISR regime (i.e.50 s −1 − 100 s −1 ) using a conventional split Hopkinson bar arrangement, as shown in Fig. 1, while a Wedge-Bar apparatus was used in the lower ISR regime (i.e. 10 s −1 − 50 s −1 ).The results show a rapid transition in the viscoelastic strain rate sensitivity of bovine cortical bone under ISR conditions, as shown in Fig. 2.
Several viscoelastic models for cortical bone have been proposed based on compression data such as that shown in Fig. 2 [3,8].However, for a more complete description of bone behaviour, knowledge of the viscoelastic relaxation response is also required.This can be accomplished using a stress relaxation tests or a cyclic test using a dynamic mechanical analyser (DMA) in which a small EPJ Web of Conferences  specimen is placed under a cyclic load [13].An analysis of the displacement response provides values for the storage and loss moduli, which are viewed as indicative of the elastic and viscous components of the material.However, these techniques are typically limited to low strain rates and are therefore not suitable for the HSR or ISR regimes [12].Furthermore, Saraf et al. [11] note that techniques involving oscillatory loading may examine different mechanisms than those developed by the single HSR loading pulse typically generated during an impact.
In this paper, a new experimental approach is proposed to investigate the dynamic viscoelastic relaxation of cortical bone in the upper ISR regime.This work is a continuation of that reported by Cloete et al. [4] and makes use of the same split Hopkinson setup with a CiT striker.It will be shown that a dynamically loaded bone specimen, when allowed to relax between two Hopkinson bars, produces stress, strain and strain rate histories that contain characteristic features which can be interpreted in the context of a standard linear viscoelastic model.Finally, a few preliminary experimental results are presented to illustrate the proposed approach.

Specimen preparation
Whole fresh femurs of slaughter-age cattle were sourced from a local butcher and stored in a frozen condition.During specimen preparation, the bones were thawed and cortical specimens were extracted from the mid diaphyses of two femurs.The specimens were machined to nominal dimensions of 4 mm diameter and 5 mm length and individually inspected, weighed, measured and refrozen.The results reported in this paper are for specimens originally orientated in the longitudinal direction relative to the bone axis.As mentioned above, this paper is a continuation of previous work reported by Cloete et al. [4] and will focus on the subset of the data in the upper ISR regime.In total, five specimens were tested at a strain rate of approximately 130 s −1 of which two did not fail and can be used to investigate dynamic viscoelastic relaxation.

Experimental method
Dynamic bone compression and relaxation tests were performed using the classic SHB technique [7] with a CiT [4], as shown in Fig. 1.As the name suggests, a CiT striker consists of a slender cone made from 6063-T6 aluminium coaxially mounted inside a tube of the same material, where the large end of the cone is threaded to the tube.The CiT striker had a length of 870 mm, with inner cone end diameters of 9 mm and 25 mm.The outer tube has a standard section and was chosen to have approximately the same cross-sectional area as the large end of the inner cone.The Hopkinson bars had a length of 3 m and a diameter of 19.05 mm and were machined from maraging steel with a density of 7950 kg m −3 and a longitudinal wave propagation speed of 4575 ms −1 .Diametrically opposed foil strain gauges (2 mm, 120 ) were bonded to the midpoint of the Hopkinson bars.The signals were amplified using custom-built amplifiers (gain of 1000, bandwidth of 100 kHz) based on the Burr Brown INA110 chip and digitized using an ADLINK PCI-9812 card (10 MHz, 12 bit).
The CiT striker design is intended to produce an incidence pulse that allows the apparent Young's modulus of a viscoelastic specimen to be measured at small strains and near-constant strain rates in the upper ISR range.This requires a shaped incidence pulse that has a steeply rising slope with the same gradient as that of the specimen response.Furthermore, the end of the pulse is sharply truncated to allow rapid load removal, which has several benefits.Firstly, it allows the test duration to be maximized for a given SHB length, i.e. the strain rate can be minimized for a given final strain.Secondly, it would allow for the recovery of an intact specimen for microstructural analysis.Thirdly, it allows for the possibility of conducting dynamic viscoelastic relaxation tests.In other words, the cone-in-tube striker design provides a means to control both the shape and duration of the incidence wave such that a bone specimen can be loaded at a constant strain rate up to a point just before fracture occurs and then fully relax without any further induced loading.This has not been possible with previously published pulse shaping concepts which have only been able to meet the steep slope requirement [1,3].The CiT striker is launched in a conventional manner using a gas gun.However, only the inner cone impacts upon the incident bar, whereas the outer tube impacts upon a large reaction mass mounted coaxially around the incident bar, as shown in Fig. 1.The resulting stress waves in the cone and tube simultaneously reach the free end of the striker where the cone and tube are connected.In isolation, the stress wave produced in a slender cone with a large ratio between its end diameters is insufficient to reverse the velocity of the large end, resulting in further loading of the incident bar.However, in the CiT striker, the reflected stress waves of the cone and tube interact and, if the tube is large enough, result in the velocity of the large end of the conical striker being reversed.This velocity reversal induces a tensile stress wave in the cone that, upon arrival at the impact face, causes the conical striker to separate from the incident bar, thus abruptly terminating the incidence pulse.In this way, a steeply rising and sharply truncated incidence pulse, with a large ratio between the initial and final stress wave magnitudes, can be created, as shown in Fig. 3.
A minor limitation of the CiT striker technique is that the required cone angle is not known a priori, and, in general, a set of inner cones with a range of cone angles would be required to cover possible variation in specimen responses.The cone angle and diameter ratios used for the CiT striker was selected based on previous experimental results and a simple spreadsheet-based one-dimensional mass-spring code and did not require any iterations [3].

Experimental results
In addition to the shaped incidence pulse, Fig. 3 shows typical reflected and transmitted signals for the case where a cortical bone specimen does not fracture.Note that after the passage of the reflected wave the stress signals in both Hopkinson bars show the same decay behaviour.This is due to the specimen exerting equal but opposite forces on the bar faces during the unloading phase of the test and shows that the specimen is in quasi-equilibrium.A similar degree of specimen equilibrium was obtained during the loading phase as shown in Fig. 4, where the stress histories at the specimen -bar interfaces have been obtained using standard SHB data processing formulae, such as those given by Gray [7].
The processed results corresponding to the raw data in Fig. 3 are given in Fig. 5, which shows the resulting stress and strain rate curves as a function of strain.The stress curve has an initial gradient of 8.5 GPa up to a strain of 0.01, followed by a loading gradient of 15 GPa.Similar initial gradients of varying duration were observed in three of the other four tests but was absent in the fourth.By contrast, there was less variation in the main portion of the loading curves where all five experiments resulted in gradients of 15.2 ± 0.7 GPa.This implies that the latter part of the loading curve represents the actual mechanical properties of the bone specimen, while the initial gradient represent transient "settling" behaviour as the specimen experiences slight adjustments until it is in full contact with the Hopkinson bar faces.This transient behaviour is also evident in the strain rate curve which rises slightly during the initial loading and then remains near-constant at approximately 140 s −1 .
During unloading the bone specimen displays significant hysteresis.After a peak stress of 290 MPa, the strain rate changes rapidly to a value of −100 s −1 , while the stress initially falls with a gradient of 21 GPa.There is a noticeable linear relationship between the strain and strain rate during the unloading phase, which will be referred to as the relaxation gradient and discussed in Sect. 5. Due to the limited test duration of a SHB apparatus, the complete specimen unloading behaviour could not be captured.In lieu of this, a simple extrapolation of the relaxation gradient would appear to indicate a residual strain of approximately 0.01, although the "settling" behaviour during the loading phase suggests a value closer to 0.005.Furthermore, preliminary measurement of the recovered specimen immediately after the test indicated that there was some residual deformation, though not as much as indicated by the SHB data.However, further measurement of the specimen several days later showed that the specimen had relaxed completely to its original length.The implication is that the deformation of cortical bone is a complex combination of several competing deformation mechanisms with some acting within less than a millisecond while others act over hours or even days.

Theoretical model
While there is clear stress relaxation behaviour shown in Figs. 3 to 5, it does not conform to the strict definition of a viscoelastic relaxation test for which, per definition, the strain is kept constant [5].Hence, none of the standard equations for assessing viscoelastic relaxation can be applied, which possibly explains the scarcity of data in the literature regarding the dynamic viscoelastic relaxation of bone.In lieu of this, it is proposed that the SHB test results can be evaluated using a standard linear viscoelastic model, which consists of a Kelvin-Voigt element (E K V and µ K V ) in parallel with a Maxwell element (E M and µ M ), as shown in Fig. 6.
During the unloading phase, the specimen applies equal but opposite forces to the Hopkinson bars, which generate stress waves, as shown in Fig. 3. Hence, based on the well known stress wave formulae [9], the particle velocities in the Hopkinson bars will be proportional to the stress applied by the specimen.In other words, the Hopkinson bars will behave like perfect linear dampers and their effect can be modelled by adding an equivalent damper µ S H B , as shown in Fig. 6.Note, however, that the equivalent SHB damper does not represents a component of the material model for the specimen and is, therefore, not included when calculating the specimen stress.
The magnitude of the equivalent damper coefficient can be obtained by expressing the relationship between stress and particle velocity in a Hopkinson bar in terms of the specimen stress and strain, which gives, where the subscripts s and b denote respectively the specimen and Hopkinson bar stress σ , cross-sectional area A, density ρ, wave speed c, strain rate ε and length L. Note that only half the specimen length is considered since both of the Hopkinson bars will be moving with the same speed v b but in opposite directions.
Hence the equivalent damper coefficient can be expressed as, From the data given in Sects. 2 and 3, the value of the equivalent damper coefficient is 2.062 MPa s −1 , which is used throughout the work presented in this paper.
Given the above, an explanation for the linear relationship between the strain and the strain rate during the unloading phase is proposed.Consider an unloading scenario where the strain rate is sufficiently high to cause the Maxwell damper µ M to "lock-up", i.e. its deformation relative to the Kelvin-Voigt element is negligible such that the Maxwell spring E M is effectively operating in parallel to the Kelvin-Voigt spring E K V .Since there is no external force, equilibrium dictates that the spring and damper forces must be equal and opposite, ( Hence, the ratio between the specimen strain and the strain rate is constant and has the form, Assuming that the loading gradient of 15 GPa in Fig. 5 is approximately equal to a combination of the Kelvin-Voigt and Maxwell springs, while neglecting the Kelvin-Voigt damper as small in comparison to equivalent SHB damper, gives an estimate for the relaxation gradient of −7274 s −1 .
Given the simplicity of the model, this result compares well with the measured value of −8300 s −1 .

Theoretical model results
A discrete version of the standard linear viscoelastic model described above was implemented in a spreadsheet and numerically integrated to replicate the observed behaviour of bovine cortical bone in dynamic relaxation tests.During the loading phase the strain rate was specified and the equivalent SHB damper was not considered, whereas during the unloading phase the equivalent SHB damper was included and the strain rate was determined via the interaction of the spring and damper components.Initial hypothetical test cases were simulated to illustrate the effect that the various parameters have on the overall response.Figure 7 shows the result for a simulation with parameter values words, the effect of the Maxwell element is removed while the effect of the Kelvin-Voigt damper is emphasized.The Kelvin-Voigt damper causes positive and negative stress jumps to occur at the beginning of the loading and unloading phases respectively.The apparent unloading gradient is slightly smaller than that of the loading phase, which is due to the decreasing strain rate during the unloading phase.The relaxation gradient is linear, as expected following the argument presented in Sect.5, and the residual strain tends to zero at the end of the simulation.
Figure 8 shows the result for a simulation with parameter values E K V = 0 MPa, E M = 17500 MPa, µ K V = 0 MPa s −1 and µ M = 5 MPa s −1 .In contrast to the previous test case, the effect of the Maxwell element has been emphasized while the effect of the Kelvin-Voigt element has been removed.Due to the absence of the Kelvin-Voigt damper, no stress jumps occur at the beginning of the loading or unloading phases.Furthermore, unlike the previous case, the apparent unloading gradient is significantly larger than that of the loading phase, due to the residual strain remaining at the end of the simulation.Note that the relaxation gradient is again linear, and, during the loading phase, where the strain rate is constant, the gradient varies slightly whereas during the unloading phase, where the strain decreases, the gradient appears to be constant.Figure 10 shows a similar simulation replicating the experimental result of a different specimen from the test series reported by Cloete et al. [4].The specimen was nominally identical to that used to obtain the results in Fig. 9 with parameter values E K V = 10500 MPa, E M = 7000 MPa, µ K V = 0.1 MPa s −1 and µ M = 0.75 MPa s −1 .
In both Figs. 9 and 10 the origin of the numerical results has been shifted to compensate for the "settling" behaviour in the experimental curve.

Discussion of results
The numerical results in both Figs. 9 and 10 show good correlation with the experimental results and have consistent parameter values.
In both tests there is a relatively small drop in the stress after the point of maximum stress, which indicates that the contribution of Kelvin-Voigt µ K V damper is small at strain rates up to 100 s −1 .This result is consistent with the data in Fig. 2 and similar results reported in the literature [3,8] where the apparent modulus begins to increase exponentially above the ISR regime.
Significant hysteresis is observed in both tests, which is reflected by the large contribution made by the Maxwell 03004-p.5 EPJ Web of Conferences element in both simulations.This result is also consistent with the data in Fig. 2 where there is a sudden transition in the apparent modulus in the ISR regime.
The numerical results begin to deviate appreciably from experimental results toward to end of the unloading phase.In particular, the simulated relaxation curves do not remain linear below a strain rate of 50 s −1 .Furthermore, the numerical simulations do not predict a significant residual strain such as that observed in the bone specimens immediately after the SHB tests.Based on the correlation in Fig. 2, which required the use of a nonlinear Maxwell element [4], it is likely that the model presented in this paper requires further development including additional elements with nonlinear properties.

Conclusions and recommendations
An experimental approach to characterize the dynamic viscoelastic relaxation behaviour of cortical bone was proposed.A standard SHB incorporating a CiT striker was used to perform dynamic relaxation tests on bovine cortical bone samples.The use of a CiT striker allowed the bone sample to be compressed at a near constant strain rate during the loading phase.During the unloading phase, the specimen relaxed between the Hopkinson bars without any further interaction with stress waves resulting from the loading phase.The technique displayed good specimen equilibrium and provided consistent results.
The Hopkinson bars were shown to act like linear viscous dampers during the unloading phase, which provided an explanation for the observed nearlinear relaxation curves.A standard linear viscoelastic material model was presented and used to replicate the experimental results.The model was able to capture several important features of the stress and strain rate curves, although detailed features toward the end of the unloading phase were not captured and further development of the model is required.Nevertheless, the preliminary results are encouraging and further use of CiT striker concept, or similar techniques, to conduct dynamic relaxation test on bovine cortical bone is recommended.

Figure 1 .
Figure 1.Layout of the split Hopkinson bar setup using a CiT striker.

Figure 2 .
Figure 2.The strain rate sensitivity of the apparent Young's modulus of bovine cortical bone specimen as reported by Cloete et al.[4].

Figure 3 .
Figure 3. Raw SHB data for a dynamic relaxation test on a bovine bone specimen.

Figure 4 .Figure 5 .
Figure 4. Specimen equilibrium during a dynamic relaxation test on a bovine bone specimen.

Figure 6 .
Figure 6.Schematic diagram of the standard linear viscoelastic model used in this work.The damper element µ S H B represents the role of the Hopkinson bars during unloading and is deactivated during the loading phase.

Figure 9
Figure9shows the result of a simulation replicating the experimental result discussed in Sect. 4 with parameter values, E K V = 9800 MPa, E M = 7500 MPa, µ K V = 0.03 MPa s −1 and µ M = 1 MPa s −1 .Figure10shows a similar simulation replicating the experimental result of a different specimen from the test series reported by Cloete et al.[4].The specimen was nominally identical to that used to obtain the results in Fig.9with parameter values E K V = 10500 MPa, E M = 7000 MPa, µ K V = 0.1 MPa s −1 and µ M = 0.75 MPa s −1 .In both Figs. 9 and 10 the origin of the numerical results has been shifted to compensate for the "settling" behaviour in the experimental curve.