Modeling of viscoelastic deformation and rate-dependent fracture damage in rat bone

Abstract

Bone is a complex hierarchical structural material whose organ-level response is highly influenced by its constitutive behavior at the microstructural level, which can dictate the inelastic nonlinear deformation and fracture within the organ. In the current study, a combined experimental-computational approach was sought to first obtain the local constitutive properties. Later, a multiscale modeling framework utilizing a novel rate-dependent nonlinear viscoelastic cohesive zone (NVCZ) model was used to explore the fracture behavior at the microstructure of the bone and its influence on the global scale (organ-level) response. Toward this end, nanoindentation testing was conducted within the cross-section of a rat femur bone specimen. An inverse optimization process was used to identify the isotropic linear viscoelastic (LVE) properties of cortical bone by integrating the test results with a finite element model simulation of the nanoindentation testing. Model results using different numbers of spring-dashpot units in the generalized Maxwell model showed that four spring-dashpot units are sufficient to capture the LVE behavior, while solely LVE constitutive relation is limited to fully characterize the rat femur. The LVE constitutive properties were then used along with the rate-dependent NVCZ fracture within the representative volume element (RVE), which was two-way coupled to the global scale bone. A parametric study was conducted by varying the fracture properties of the NVCZ model. The model demonstrated the capability and features to represent inelastic deformation and nonlinear fracture that are linked between length scales. This further implies that the inelastic fracture model and the two-way coupled modeling can elucidate the complex multiscale deformation and fracture of bone, while model validation and further advancements with test results remain a follow-up study and are currently in progress.

Introduction

Most biological media have complex hierarchical structures that help or promote many essential features that drive the mechanical performance, and bone is no exception. Biological activities change the integrity of complex hierarchical structure in bone with time (i.e., age) and disease (e.g., osteoporosis) which determines its constitutive behavior and fracture resistance characteristics at multiple length scales (Sabet et al., 2016; Ural and Vashishth, 2014). Integrating and incorporating such complex features into numerical models for predicting the mechanical response of this tissue can be quite challenging. In a move towards a patient-specific analysis and diagnosis, numerical models generated from the imaging techniques play a crucial role and are indispensable tools (Alcântara et al., 2019; Podshivalov et al., 2014) in medical diagnostics. Critical to such models is incorporating the right constitutive behavior of its constituents at the characteristic length scale and to identify the right fracture properties and its damage process. Incorporating the constitutive behavior and fracture characteristics at different length scales into a multiscale modeling approach is an apparent direction in order to gain meaningful insight into the inelastic behavior of bone (Fish and Hu, 2017; Lloyd et al., 2015; Sabet et al., 2016; Ural, 2020; You et al., 2017). Such modeling approaches can help investigate the influence of diseases such as osteoporosis on bone stiffness, strength and fracture toughness.

In order to obtain the constitutive properties, several traditional approaches such as conducting compressive or tensile tests, pure shear tests, torsion tests, bending and fatigue tests have been conducted on bone specimens to characterize the global (or organ-level) constitutive properties. Using such traditional approaches, bone response at the organ-level has been characterized as isotropic elastic (Zysset et al., 1994), transversely isotropic (Xn and Guo, 2004), orthotropic elastic (Goldstein, 1987; Li et al., 2013), viscoelastic (Iyo et al., 2004; Yamashita et al., 2001), or elastic-plastic (Rincón-Kohli and Zysset, 2009). Nonetheless, complexity and difficulty increase when attempting to characterize the constitutive behavior at the lower length scales using similar traditional testing methods and approach (Groetsch et al., 2019). In recent years, the use of fractional calculus to model complex bodies such as biological materials has gained significant attention (Sumelka et al., 2020; Sun et al., 2018). Numerous investigations have demonstrated that fractional operators serve as effective tools in describing memory effects such as viscoelasticity (Bagley and Torvik, 1983; Shymanskyi and Sokolovskyy, 2021), plasticity (Sumelka, 2014), transport phenomena, multiscale phenomena in materials (Failla and Zingales, 2020), and modeling damage (Suzuki et al., 2021).

Of interest to the current study is the nanoindentation technique which has gained significant popularity over recent years and has been widely used for characterizing mechanical properties of various materials at the microscale (Oliver and Pharr, 2004). The nanoindentation experiment in itself is not a constitutive property deriving experiment in the sense that one cannot directly obtain the constitutive properties unlike traditional testing methods, which poses a serious restriction in obtaining the mechanical properties at the microscale. The displacement vs. load data from a nanoindentation experiment is a culmination of contact problem between the indenter and the microscale material response (Popov et al., 2019). The constitutive response of the microstructure is obtained through a solution to an initial boundary value problem (IBVP) based on contact theories. Analytical expressions have been derived for such IBVP for specific constitutive behavior (Sabet et al., 2016; Vandamme and Ulm, 2006; Yu and Blanchard, 1996). Deriving such analytical expressions for any given constitutive property of the microstructure can be quite tedious, and closed form solutions might not exist. Under such circumstances, characterizing the constitutive behavior of the microstructure based on an inverse optimization method using numerical modeling is an attractive approach which is explored in this study.

Apart from the constitutive behavior of the individual constituents, understanding the facture resistance characteristics and fracture mechanism of bone tissue is critical to the multiscale numerical analysis. Inelastic response of the bone tissue at the global-scale has been captured effectively (Johnson et al., 2010; Li et al., 2013; Mercer et al., 2006) and there is enough evidence to suggest that the inelastic behavior at the organ level is primarily due to the microstructural features of individual constituents and their corresponding fracture toughening and resistance characteristics (Burr et al., 1998; O'Brien et al., 2005; Ritchie et al., 2009; Taylor et al., 2007; Vashishth, 2007). Based on the above literature it is clear that the complex inelastic behavior of bone can be attributed mostly to its heterogenous hierarchical buildup at the microstructural level. Hence multiscale modeling approach that can capture such a complex response at the global-level by utilizing the microstructural features and the constitutive behavior of its heterogeneities along with the rate-dependent fracture characteristics is key to understand the complex bone behavior.

This study presents multiscale modeling of viscoelastic deformation and rate-dependent fracture in rat bone. The general multiscale modeling concept that we have used for different types of inelastic-heterogenous materials was applied to this study with several novel aspects included. Theoretically, we developed and implemented the rate-dependent fracture with Gaussian damage evolution in the two-way coupled multiscale framework. Experimentally, we used realistic rat bone test results (Liu et al., 2021) and their characterization to obtain local-scale properties (i.e., microstructure and nanoindentation results), and computationally we used the inverse optimization process to identify local viscoelastic properties of rat bone tissue and conducted the two-way coupled multiscale simulation with the rate-dependent nonlinear viscoelastic cohesive zone implemented.

More specifically, the viscoelastic local-scale bone properties were characterized using nanoindentation test and its finite element modeling incorporated with an inverse optimization method. Once the constitutive properties were identified through the inverse optimization approach, the local properties were incorporated within a two-way coupled multiscale modeling framework where a global scale bone was linked to its local scale heterogeneous representative volume element (RVE). Fracture response and damage growth within the microstructure and its influence on the global scale response was modeled through the nonlinear viscoelastic cohesive zone (NVCZ) which is incorporated with a novel Gaussian damage evolution. The efficacy of the Gaussian damage evolution criterion to simulate different fracture response and rate dependent behavior is shown through parametric case studies. The capability and features of the multiscale model that is sensitive to varying loads, geometry, and local material properties were demonstrated through parametric case studies.

Modeling framework

The two-way coupled multiscale modeling has been utilized to investigate different heterogenous materials (Kim et al., 2013; Rami et al., 2018; Souza and Allen, 2010) specifically to predict the global-scale response under the influence of the microstructural details. A similar approach was considered in this study. The multiscale modeling framework is based on the fundamental idea that heterogeneous materials such as bone, that have a statistically homogenous distribution of its heterogeneities, their global-scale behavior can be represented through their RVEs (Allen, 2001; Hashin, 1983; Hill, 1972; Hua Qin, 2008). In order to do so, one needs to identify an RVE of the microstructure that can essentially capture the overall global response, and the RVE size must be sufficiently large to describe the overall behavior of the global-scale object. The criterion for the selection of the RVE is described in detail by others (Gitman et al., 2007). Once the RVE is established the overall effective constitutive behavior of the global-scale object can be obtained through the volume averages of the field variables over the RVE under uniform boundary conditions.

Figure 1 shows a cross section of a rat femur (global object) investigated in this study along with its local RVE. The length scale associated with the RVE should be small compared to the global object under investigation and at the same time the RVE size should be large enough to include the heterogenous components of the global-scale object. It has to be pointed out here that although the rat femur utilized in the current study has significantly different microstructure compared to that of human bone, the analysis and methodology presented here can be easily extended to the human bone. Described next is the two-way coupled multiscale modeling framework that directly links the global response to the local scale RVE response. Theoretical details of the two-way coupled multiscale modeling framework can be found in previous studies applied to several different materials (Allen, 2001; Kim et al., 2013; Rami et al., 2018; Souza and Allen, 2010); however, key equations in a summarized form are presented in this paper for a complete discussion.

Figure 1.

Cross-section of the rat bone and its microstructure presenting different phases.

IBVP of global and local scale

In the multiscale framework, IBVP is solved separately at the global and local scale in succession. The global scale mechanical response is obtained by solving the IBVP that is a consequence of the boundary condition, balance laws (conservation of linear momentum and angular momentum), small strain-displacement relation and the stress-strain constitutive relation. Equations (1) to (4) show the appropriate relations describing each of the conditions stated above.

$${\sigma }_{ji,j}^G\left(x_m,t\right)=0,\mathrm{ }{x}_m\in V^G$$

(1)

$${\sigma }_{ij}^G\left(x_m,t\right)={\sigma }_{ji}^G\left(x_m,t\right),\mathrm{ }{x}_m\in V^G$$

(2)

$${\epsilon }_{ij}^G\left(x_m,t\right)=\frac{1}{2}\left(u_{i,j}^G(x_m,t)+u_{j,i}^G(x_m,t)\right),\mathrm{ }{x}_m\in V^G$$

(3)

$${\sigma }_{ij}^G\left(x_m,t\right)={\overline{\mathrm{\Omega }}}_{\tau =\mathrm{\infty }}^{\tau =t}\left\{{\epsilon }_{kl}^G\left(x_m,t\right)\right\}$$

(4)

where ${\sigma }_{ij}^G\left(x_m,t\right), {\epsilon }_{ij}^G\left(x_m,t\right)$ and $u^G(x_m,t)$ are the global scale Cauchy stress tensor, infinitesimal strain tensor and displacements, respectively at each material point ( $x_m$ ) and time (t), $V^G$ is the global scale volume, ${\overline{\mathrm{\Omega }}}_{\tau =\infty }^{\tau =t}$ is the functional that describes the homogenized constitutive relation between the stress tensor and strain tensor at the global scale, and the superscript G denotes global scale.

Similar to the global scale, one can describe the IBVP of local scale RVEs by utilizing the same set of conditions. It is natural to assume the local scale RVE to be governed under the same set of conditions as the global scale, although it is relatively smaller than the global body. Equations (5) to (8) show the local scale equations.

$${\sigma }_{ji,j}^L\left(x_m,t\right)=0,\mathrm{ }{x}_m\in V^L$$

(5)

$${\sigma }_{ij}^L\left(x_m,t\right)={\sigma }_{ji}^L\left(x_m,t\right),\mathrm{ }{x}_m\in V^L$$

(6)

$${\epsilon }_{ij}^L\left(x_m,t\right)=\frac{1}{2}\left(u_{i,j}^L(x_m,t)+u_{j,i}^L(x_m,t)\right),\mathrm{ }{x}_m\in V^L$$

(7)

$${\sigma }_{ij}^L\left(x_m,t\right)={\mathrm{\Omega }}_{\tau =\mathrm{\infty }}^{\tau =t}\left\{{\epsilon }_{ij}^L\left(x_m,t\right)\right\},\mathrm{ }{x}_m\in V^L$$

(8)

where ${\sigma }_{ij}^L\left(x_m,t\right), {\epsilon }_{ij}^L\left(x_m,t\right)$ and $u^L(x_m,t)$ are the local-scale stress tensor, strain tensor and displacements, respectively at the each material point ( $x_m$ ) and time (t), $V^L$ is the local-scale volume, ${\mathrm{\Omega }}_{\tau =\infty }^{\tau =t}$ is the functional that describes the constitutive relation between the Cauchy stress tensor and the strain tensor, and the superscript L denotes local scale.

Nonlinear viscoelastic cohesive zone (NVCZ) model

The deterioration in strength and loss of fracture toughness at the global scale is closely associated to local scale microstructural cracks that develop and propagate once a critical state is reached (Seref-Ferlengez et al., 2015; Zimmermann et al., 2015). Such loss in strength and fracture resistance characteristics can be effectively captured by using appropriate damage models within the local scale RVE. Damage in terms of crack initiation and propagation can be considered within individual phases or at the interface between the individual phases. In this study, bone damage within the interstitial matrix is considered by using the cohesive zone (CZ) model. CZ models have been used to model fracture of a number of different materials including bone (Alcântara et al., 2019; Ural and Mischinski, 2013; You et al., 2017) utilizing different damage evolution criteria. For the bone microstructure in the current study, a micromechanics-based nonlinear viscoelastic cohesive zone (NVCZ) model was utilized. Equations (9) and (10) describe the NVCZ model (Allen and Searcy, 2001) which represents a constitutive relation between cohesive zone traction and separation displacements within cohesive zone of a viscoelastic material:

$$T_i^L\left(t\right)=\frac{u_i^L\left(t\right)}{{\lambda }^L\left(t\right){\delta }_i^L}\left[1-{\alpha }^L\left(t\right)\right]\left[{\sigma }_i^L+{\int }_{t_o}^t{E}_{CZ}^L(t-\tau )\frac{\partial {\lambda }^L\left(t\right)}{\partial \tau }\text{d}\tau \right]on\mathrm{ }\partial V_{CZ}^L$$

(9)

$${\lambda }^L\left(t\right)=\sqrt{{\left(\frac{u_n^L}{{\delta }_n^L}\right)}^2+{\left(\frac{u_t^L}{{\delta }_t^L}\right)}^2}$$

(10)

where $T_i^L\left(t\right)$ is the CZ traction; $u_i^L\left(t\right)$ is the CZ opening displacement; ${\delta }_i^L$ is the CZ length scale parameter; the subscript $i$ is the intrinsic normal or tangential directions within the CZ; ${\alpha }^L\left(t\right)$ describes the internal damage evolution; ${\sigma }_i^L$ represents the critical stress state for the CZ to be activated; $E_{CZ}^L\left(t\right)\mathrm{ }$ is the linear viscoelastic relaxation modulus; ${\lambda }^L\left(t\right)$ is Euclidian norm of the CZ opening displacements. A rate-dependent evolution of the internal damage parameter, ${\alpha }^L\left(t\right)$ is considered in this study as shown in equations (11) and (12).

$$\left\{\begin{array}{l}\dot{\alpha }\left({\lambda }^L,{\dot{\lambda }}^L\right)=\beta \left({\lambda }^L\right) {\left[{\dot{\lambda }}^L\right]}^m, \text{when} {\dot{\lambda }}^L\left(t\right)>0\mathrm{ and\; \hspace{0.17em}\hspace{0.17em}}{\alpha }^L\left(t\right)<1\\ \dot{\alpha }\left({\lambda }^L,{\dot{\lambda }}^L\right)=0, \text{when} {\dot{\lambda }}^L\left(t\right)\le 0\mathrm{ or\; \hspace{0.17em}\hspace{0.17em}}{\alpha }^L\left(t\right)=1\end{array}\right.$$

(11)

$$\mathrm{\beta }\left({\lambda }^L\right)=\mathrm{A}\left[\exp \left(\frac{-{\left({\lambda }^L-\overline{\lambda }\right)}^2}{2{\Delta }^2}\right)\right]$$

(12)

where, $\dot{\alpha }\left({\lambda }^L,{\dot{\lambda }}^L\right)$ is the time derivative of the damage parameter which depends on the norm of the opening displacement $\left({\lambda }^L\right)$ and its time derivative $\left({\dot{\lambda }}^L\right)$ . As shown in equation (12), the factor $\mathrm{\beta }\left({\lambda }^L\right)$ has a Gaussian distribution with parameters A, $\overline{\lambda }$ and $\Delta$ . The advantage of such a damage evolution criterion is that rate-dependent inelastic behavior and fracture characteristics can be captured utilizing a single set of parameters for the damage evolution law. Equations (9) to (12) depicts the essential parts of the NVCZ model. It should be noted that the NVCZ elements are inserted adaptively within the interstitial matrix of the RVE when the initiation criterion (i.e., ${\sigma }_i^L$ : the critical stress state for the CZ to be activated) is met. Once adaptively inserted, their behavior is governed by equations (9) to (12) which essentially describes the traction evolution across the NVCZ element that transfers the load from one face to the other across the cohesive element. It can be noted in equation (9) that the traction response is time-rate and history-dependent owing to the viscoelastic constitutive response $E_{CZ}^L(t$ ) and rate-dependent damage evolution function. When the damage parameter ${\alpha }^L\left(t\right)$ becomes 1.0, then the traction is zero within the CZ element. This indicates a complete failure of the particular CZ element, which leads to a new surface (i.e., crack) created.

Homogenization for linking local to global scale

As described earlier, in the two-way coupled multiscale analysis, IBVP corresponding to both the global and local scales are solved in succession. The multiscaling process is basically accomplished by constructing a time-stepping algorithm. Briefly, the smaller length scale (i.e., local scale RVE) accounts for the material heterogeneity (due to the different phases placed) and various constitutive relations with fracture wherever the physics requires. RVE finite element models are constructed, and the results of this length scale are linked to the global scale (organ-level bone) finite element model via a homogenization theorem. Two-way coupling between the length scales is accomplished through a series of steps: 1) read inputs for different length scales; 2) obtain initial homogenized tangent tensor for the global scale problem; 3) solve the global scale IBVP at a given time-step; 4) apply a global scale solution to a local scale RVE via corresponding boundary conditions; 5) solve the local scale IBVP incorporated with its microscale fracture process at the time-step; 6) homogenize local scale RVE results affected by microscale fracture process; and 7) update homogenized local scale results to the global scale model at each integration point for the next time-step. This process occurs in a recursive manner until the end of the numerical procedure.

The global scale resultant solution in the form of the field variable, either the stress or the strain tensor is passed to the local scale to apply the uniform boundary conditions. Either uniform traction or linear displacement boundary conditions are applied to the surface of the local scale RVE as described in equations (13) and (14), respectively:

$$T_i={\sigma }_{ij}^o{n}_j$$

(13)

$$u_i={\epsilon }_{ij}^o{x}_j$$

(14)

where $T_i$ is the uniform traction; ${\sigma }_{ij}^o$ is the constant stress tensor; $\mathrm{ }{n}_j$ is the unit outward normal to the surface boundary of the RVE; $u_i$ is the resultant displacement; ${\epsilon }_{ij}^o$ is the constant strain tensor; x_j is the location at the boundary surface. In this study, linear displacement boundary conditions were applied at the surface of the local scale RVE.

Once the local scale IBVP is solved subjected to either of the conditions shown in equations (13) and (14), the global scale field variables can be obtained utilizing the mean-field homogenization theorem (39). The global scale stress tensor can be obtained using equation (15) and the global scale strain tensor can be obtained using equations (16) to (18):

$${\sigma }_{ij}^G={\overline{\sigma }}_{ij}^L=\frac{1}{V^L}{\int }_{\partial V_E^L}^{}{\sigma }_{ki}^L{n}_k^L{n}_j^L\text{d}S$$

(15)

$${\epsilon }_{ij}^G={\overline{\epsilon }}_{ij}^L={\overline{E}}_{ij}^L+{\overline{e}}_{ij}^L$$

(16)

$${\overline{E}}_{ij}^L\equiv \frac{1}{V^L}{\int }_{\partial V_E^L}^{}\frac{1}{2}\left(u_i^L{n}_j^L+u_j^L{n}_i^L\right)\text{d}S$$

(17)

$${\overline{e}}_{ij}^L\equiv \frac{1}{V^L}{\int }_{\partial V_I^L}^{}\frac{1}{2}\left(u_i^L{n}_j^L+u_j^L{n}_i^L\right)\text{d}S$$

(18)

where ${\overline{\sigma }}_{ij}^L$ and ${\overline{\epsilon }}_{ij}^L$ is the volume average of the stresses and strains over the RVE boundary, ${\overline{E}}_{ij}^L$ and ${\overline{e}}_{ij}^L$ are the volume average strains over the external $\left(\partial V_E^L\right)$ and internal $\left(\partial V_I^L\right)$ boundaries, respectively; $n_j^L$ is the outward normal to the internal or external boundary surface. Since the numerical problem of the multiscale analysis is conducted in incremental time steps, the incremental stress and strain at the global scale can be related through the stiffness tangent tensor which is obtained utilizing equations (15) to (18). A detailed description of the numerical scheme to obtain the tangent stiffness tensor is described elsewhere (Souza and Allen, 2010).

Rat bone specimen and test results

The bone specimen used in this study was a rat femur. A total of 40 specimens were scanned using a micro-CT (SkyScan 1172-D, Kontich, Belgium) to identify and quantify bone mineral density and other physical characteristics. Among the 40 specimens, one representative specimen was used to dissect a cross section for further investigation into its mechanical properties using nanoindentation experiment. The objective of the nanoindentation experiment here is to characterize the material constitutive behavior at the local scale. As depicted in Figure 1, the microstructure of the rat bone specimen can be assumed to be a two-phase system comprising of the interstitial matrix and voids (Haversian canals and lacunae). As such, the mechanical response from the nanoindentation testing can help capture the interstitial matrix properties.

The bone specimen for the indentation experiment was prepared in accordance with the sample preparation protocol established in an earlier study (Kim et al., 2015). The nanoindentation testing was performed using a Berkovich diamond tip indenter (Kim et al., 2010, 2013). To ensure a smooth surface, prior to the indentation experiment, a meticulous polishing process was carried out on the dissected cortical bone specimens. This process involved the use of silicon carbide abrasive papers with decreasing grit sizes (1200, 2400, and 4000 grits) and Al₂O₃ pastes (1 and 0.3 µm) on soft polishing cloths. Several indentations were made at the anterior and posterior locations as shown in Figure 2(a) to obtain statistically significant experimental data. The indentation test consisted of a loading stage, a hold stage and an unloading stage as shown in Figure 2(b). The test was performed in a displacement-controlled mode for the loading stage by increasing the load till a peak load of 5 mN, and then the indenter was held at the constant load for a period of 30 s (for the holding stage: creep) followed by the unloading stage (linearly decreasing load to zero). The typical load-displacement response from the nanoindentation testing is depicted in Figure 2(c). A total of 45 indentations were made, which resulted in 43 valid indentation results for further analysis after removing 2 indentations that presented erroneous data (e.g., not detectable or a very low modulus) due to probable indentations related to voids and/or debris. Figure 3 shows the 43 indentation results (displacement-time curve) with an average along with the error bars for the ± standard deviation. The average response curve was then utilized within the inverse optimization procedure described in the next section to obtain the LVE properties of the interstitial matrix phase.

Figure 2.

Nanoindentation testing of rat bone: (a) indentation points at the anterior-posterior regions; (b) applied load profile and (c) typical load-displacement results from the indentation.

Figure 3.

Nanoindentation displacement responses from 43 indents and the average of the results.

Inverse optimization of local scale constitutive properties

The average indentation response shown in Figure 3 was used to characterize the LVE constitutive behavior of the interstitial matrix. Direct estimation of constitutive behavior from the indentation results is restrictive due to the complex boundary conditions applied during the experiment and mostly due to the lack of analytical solutions for such experiment. Hence, an inverse optimization procedure that integrates nanoindentation testing procedure with its finite element model was used to identify the local-scale constitutive properties. One advantage of such an approach is that the constitutive properties at a local scale microstructure (component-level) can be obtained, and it can then be incorporated within a multiscale modeling framework by using the local scale properties of heterogeneities which are linked with global scale behavior. Another advantage is that any arbitrary loading patterns and boundary conditions can be used as the testing is simulated by the finite element modeling. To that end, a two-dimensional (2-D) axisymmetric model for the microstructure of the bone within a finite element framework was created to simulate the nanoindentation experiment. Isotropic LVE constitutive behavior of bone was attempted as the first trial to examine the method. To model LVE behavior, the generalized Maxwell mechanical analog represented by a Prony series was utilized. A nonlinear optimization protocol was developed and used to identify the isotropic LVE properties by minimizing the error between the nanoindentation test results and finite element model simulations.

Finite element modeling of nanoindentation testing

In this study, the bone and the indenter were modeled as 2-D axisymmetric using ABAQUS®. The bone was assigned linear viscoelastic properties while the Berkovich indenter was modeled as a rigid body with a conical shape and a half cone angle of 70.3°. The specific cone angle was selected as it applies the similar depth to contact area ratio between the Berkovich’s pyramidal and conical shape. The 2-D axisymmetric configuration, which is obviously limited to exactly representing the Berkovich’s pyramidal tip, can greatly reduce the computational effort, as it enables 2-D approximation, it was thus attempted in this study.

The isotropic LVE constitutive behavior of bone specimen is governed by the convolution integral given in equation (19):

$${\sigma }_{ij}\left(t\right)={\int }_0^{t}G\left(t-\tau \right)\frac{\partial {\epsilon }_{ij}}{\partial \tau }\text{d}\tau +{\int }_0^t\left[K\left(t-\tau \right)-\frac{2}{3}G\left(t-\tau \right)\right]{\delta }_{ij}\frac{\partial {\epsilon }_{kk}}{\partial \tau }\text{d}\tau$$

(19)

$$G\left(t\right)=G_{\infty }+\sum _{i=1}^{i=N}{G}_i{e}^{\left(\frac{-t}{{\rho }_i}\right)}$$

(20)

$$K\left(t\right)=K_{\infty }+\sum _{i=1}^{i=N}{K}_i{e}^{\left(\frac{-t}{{\rho }_i}\right)}$$

(21)

where G(t) and K(t) are the shear and bulk relaxation modulus, respectively, each of which can be represented by Prony series description shown in equations (20) and (21). $G_{\infty },G_i$ are spring constants in the shear relaxation modulus, $K_{\infty },K_i$ are spring constants in the bulk relaxation modulus, and ${\rho }_i$ are the relaxation times. A time independent Poisson’s ratio of 0.3 was assumed in the modeling, such that the shear and bulk modulus are linearly related. The nonlinear optimization protocol described earlier was utilized to obtain the Prony series parameters that represent LVE properties of bone interstitial matrix.

The mesh and the boundary conditions applied to the model are shown in Figure 4. As shown, the finite element mesh is composed of axisymmetric 4-node quadrilateral CAX4R elements and 3-node triangular CAX3 elements. A finer mesh assuring the numerical mesh convergence was used near the indenter tip for both the parts because of high deformations right below the tip and within regions where the indenter surface is likely to come into contact during the simulation. Horizontal displacements along the axis of symmetry and vertical displacements at the bottom of the specimen were constrained as shown in Figure 4(b). The linear viscoelastic properties of the bone are the key material properties that are inversely calibrated through the optimization process which is described next.

Figure 4.

Model simulation of nanoindentation testing: (a) model geometry and finite element mesh and (b) boundary conditions imposed.

Inverse optimization process

The inverse optimization process was developed in Matlab^TM as a sequence of function files that were executed in succession to control the model simulation in ABAQUS® and pass information to the optimization function. As described, the process is iterative and the Matlab^TM functions control: writing the input file for the finite element modeling, running the model and obtaining the simulation results and comparing with the experimental data. The simulation results include tip displacements and loads as a function of time. Displacements resulting from the model simulation are compared to the experimental data and the error between the two results is minimized using the optimization function described in equation (22).

$${\underset{\chi \in {\mathbb{R}}^N}{\min }{\mathrm{\Psi }}_u\left(\chi \right)=‖\mathit{u}-{\mathit{u}}_{\mathit{FEM}}\left(\chi \right)‖}_2$$

(22)

where, ${\mathrm{\Psi }}_u\left(\chi \right)$ is the objective function, $\chi$ is a vector of size $N$ material properties, ${\mathit{u}}_{\mathit{FEM}}\left(\chi \right)$ is tip displacement results from the finite element model simulation, $\mathit{u}$ is the experimental tip displacement results, and ${‖‖}_2$ is the Euclidian norm. Equation (22) is solved for the unknown constitutive property parameters using the Nelder-Mead algorithm implemented in the fminsearch function in Matlab^TM.

Results from the inverse optimization approach

Since the nanoindentation testing condition is based on the load-controlled mode, the same loading profile was applied to a reference node of the rigid indenter. The loading profile is divided into three stages: (i) monotonically increasing stage; (ii) holding stage; and (iii) linearly decreasing (unloading) stage. Since the number of Prony series terms (i.e., spring-dashpot (S-D) units) to be utilized for representing bone viscoelasticity is not known (equations (20) and (21)), a parametric study was conducted by considering three different numbers of S-D units (i.e., 4, 5 and 7). The finite element model simulation results in the form of displacement-time curves after optimization are presented in Figure 5, and the Mises stress contour at different loading stages resulting from the seven S-D units is depicted in Figure 5 as an exemplary illustration.

Figure 5.

Model simulation results of nanoindentation testing: (a) predicted tip displacements over time when different number of S-D units were used and (b) stress contours within the bone specimen at different stages of displacements.

The resulting viscoelastic model parameters are summarized in Table 1. From Figure 5, it can be generally noted that the inverse optimization method which integrates arbitrary testing with its finite element modeling can capture complex constitutive properties of inelastic materials such as bone. This implies that testing such as nanoindentation can be conducted to obtain small-scale (component-level or even smaller) response in heterogeneous microstructures, and the response is used to identify constitutive properties that are not often easy to obtain from conventional test methods. This can provide significant benefits such as component-level characterization and multiscale modeling-analysis to better understand the effects of small scale components on larger scale properties and performance. Figure 5 also indicates that four S-D units in the generalized Maxwell model are sufficient to capture the viscoelastic response of the bone interstitial matrix in this study and increasing the number of S-D units (i.e., Prony terms) did not significantly change the simulation results. Although some promising results were found, as shown in Figure 5, LVE constitutive relation is limited to fully characterize bone behavior, which was demonstrated by a much larger recovery in displacements observed from the finite element model simulation compared to nanoindentation test results. The discrepancy can be better captured using a constitutive relation such as viscoelastic-plastic that is expected to more accurately predict the larger irrecoverable displacements at the end of the unloading stage.

Table 1.

Prony series parameters resulting from the inverse optimization process.

No. of Terms	Four S-D Units		Five S-D Units		Seven S-D Units
	$G_{\infty }$ (MPa)		$G_{\infty }$ (MPa)		$G_{\infty }$ (MPa)
1	3362.89		3299.07		3303.46
	$G_i$ (MPa)	${\rho }_i$ (s)	$G_i$ (MPa)	${\rho }_i$ (s)	$G_i$ (MPa)	${\rho }_i$ (s)
2	3117.83	0.063	1104.12	0.0058	1290.35	0.053
3	2937.78	1.795	1054.28	0.0052	1244.70	0.114
4	1737.38	5.12	1309.68	0.6	1094.42	2.669
5	–	–	1633.42	10.5	686.53	12.831
6	–	–	–	–	188.77	48.129
7	–	–	–	–	0.070	668.276

Two-way coupled multiscale modeling with the NVCZ

The two-way coupled multiscale modeling was conducted to evaluate the NVCZ model and bone viscoelasticity for predicting rate-dependent behavior and local scale microcracking and its progression to global scale bone properties. Figure 6 shows the schematic of the global and local scale mesh of the rat femur obtained from a cross-sectional image.

Figure 6.

Geometry and finite element model for the two-way coupled multiscale simulation: (a) global scale femur cross-section and (b) local scale RVE and its microstructure.

The global scale image depicted in Figure 6 was taken right before the indentation experiment to ensure that the same geometry is used for the global scale finite element modeling. The same image was then used to construct the global scale mesh based on the information corresponding to the length associated with each pixel (0.8 µm per unit pixel) as depicted in Figure 6(a). The resolution of the femur cross-section was high enough to obtain a local scale microstructural details of the rat femur as shown in Figure 6(b) and construct a local scale mesh of the representative volume element (RVE). At the global scale, the bone model is regarded to have two parts comprising of the cortical bone and marrow, whereas the local scale RVE of the cortical bone was modeled with two microstructural phases: interstitial matrix and voids.

For the current study, the RVE size was determined using several microstructure images obtained from the global scale image of the femur cross-section. Since the rat bone microstructure comprises of only two phases (interstitial matrix and voids), a quantitative variation on the area fraction of voids was used to determine the RVE. Several square grids at different locations were used to obtain local images which were segmented to obtain the void percentage within each grid as depicted in Figure 7(a). Variation in the void % within each RVE was plotted against the RVE-size as shown in Figure 7(b). As shown in Figure 7(b) high variation in void % was observed when the RVE size was small, and this variation decreased with increase in RVE size. Based on the plot, an RVE size of 0.24 mm was selected that showed reduced fluctuations in void % and from the fact that an RVE size of 0.4 or 0.5 mm would increase the computational cost as several elements at the global-scale would be linked to RVE.

Figure 7.

(a) Different RVE sizes obtained at four different regions along the rat-femur cross section and (b) % void variation for different trial RVE sizes.

Figure 8 shows the boundary conditions applied to the global scale bone to simulate the multiscale behavior of rat bone influenced by local scale (i.e., RVE) microstructural properties and global scale loading conditions (such as loading rates in this study). In Figure 8, global scale elements around the loading point (as shown in red) were assigned for the two-way coupled multiscaling by linking with a local scale RVE whose microstructure is modeled as a two-phase system of the interstitial matrix and the haversian canal. Although the entire global scale elements can be linked to their corresponding local scale RVEs for multiscaling, this was not considered in this study, as it will increase the computational costs with a marginal level of modeling accuracy, and this study did not target model validation at this stage. Thus, only critical region in the global scale and a single type of RVE microstructure for the local scale was considered in this study, which can still demonstrate the capability and features of the modeling method.

Figure 8.

Schematic illustration of the multiscale bone modeling that connects global scale and local scale RVE with NVCZ for fracture.

Table 2 shows the constitutive properties that were considered for each phase at the global scale and at the local scale. As aforementioned, four S-D units in the generalized Maxwell model were sufficient to capture the viscoelastic response of the bone interstitial matrix, although it was limited to accurately characterize the irrecoverable behavior due to unloading. This might be better predicted by a viscoelastic-plastic constitutive relation which requires further investigation. Thus, for simplicity, the four S-D units were used to conduct the two-way coupled multiscale model simulation in this study.

Table 2.

Material properties used for the two-way coupled multiscale model simulation.

Global scale			Local scale
No. of terms	Four S-D units		No. of terms	Four S-D units
	$G_{\infty }$ (MPa)	$\nu$		$G_{\infty }$ (MPa)	$\nu$
1	3362.89	0.3	1	3362.89	0.3
	$G_i$ (MPa)	${\rho }_i$ (s)		$G_i$ (MPa)	${\rho }_i$ (s)
2	3117.83	0.063	2	3117.83	0.063
3	2937.78	1.795	3	2937.78	1.795
4	1737.38	5.12	4	1737.38	5.12
Bone marrow (Global)			Haversian canal
E [MPa]		$\nu$	E [MPa]		$\nu$
0.1		0.49	0.001		0.49

For the viscoelastic properties of global scale elements, it should be noted that the effective viscoelastic properties can be obtained by a homogenization process of the local scale RVE which includes both interstitial matrix and haversian canal; however, the 4-term Prony series used for the interstitial matrix were used for the global scale viscoelastic properties as the volume fraction of haversian canal in the local scale RVE was very low. This simplification did not induce any significant changes in the model simulation results. The bone marrow was assumed to be linear elastic with a small elastic modulus (0.1 MPa) and a Poisson’s ratio of 0.49.

Using the mechanical properties in Table 2, parametric analyses were then conducted by systematically changing the local scale fracture parameters of the NVCZ model. The interstitial matrix in the RVE allows insertion of NVCZ elements to simulate rate-dependent microcrack initiation and growth within the matrix. The NVCZ elements are adaptively inserted once the stress reaches a critical state in any locations within the matrix.

Two critical locations (i.e., RVE-1 and RVE-2) shown in Figure 8 were arbitrarily selected to illustrate local scale RVE behavior influenced by different NVCZ model parameters. RVE-1 was selected as it is under the maximum tensile stress, while RVE-2 experiences the highest compressive stress due to the load. Table 3 lists the cases attempted for the parametric study. As shown, for each case, only one parameter was varied keeping the others constant to investigate the influence of the varying parameter on both local and global scale behavior. All six cases in the table were investigated under the same loading rate (R1 = 0.0667 mm/s). It should also be noted that only Mode-I fracture (opening) was considered for simplicity by setting a high value for the tangential threshold ( ${\sigma }_t$ ) initiating Mode-II fracture.

Table 3.

NVCZ fracture parameters for the parametric analyses (loading rate of 0.0667 mm/s.

No.	Parameters	Mode-I fracture
		Case-1	Case-2	Case-3	Case-4	Case-5	Case-6
1	${\sigma }_t$ [MPa]	2E + 05	2E + 05	2E + 05	2E + 05	2E + 05	2E + 05
2	${\delta }_t$ [mm]	1	1	1	1	1	1
3	${\sigma }_n$ [MPa]	400; 300; 200; 150	200	200	200	200	200
4	${\delta }_n$ [mm]	1	1; 0.75; 0.5; 0.25	1	1	1	1
5	A	150	94.83	94.83; 125; 150; 200	150	150	150
6	m	1	1	1	1.0; 0.75; 0.5;0.25	1	1
7	$\overline{\lambda }$	0.005	0.005	0.005	0.005	0.03; 0.01; 0.008; 0.001	0.05
8	$\mathrm{\Delta }$	0.0025	0.005	0.005	0.005	0.01	0.0025; 0.0035; 0.005; 0.01

Note. The bold values are simply to highlight the specific parameter values that were changed for the model simulations while keeping the remaining parameters constant.

The parametric study was mainly conducted to assess the capability and features of the multiscale modeling utilized in this study to model complex inelastic-nonlinear fracture behavior of bone matter. Although the parametric analysis herein was not conducted with direct experimental data for model validation, it was based on parameters informed by studies (Arias-Moreno et al., 2020; Cory et al., 2010) or observations from experimental results in other studies (Crowninshield and Pope, 1974; Johnson et al., 2010; McElhaney, 1966). It is limited at this stage due to the lack of model validation; however, this parametric analysis can still demonstrate how the model can account for multiscale behavior that is sensitive to critical variables such as loads, geometry, and local material properties.

As an example, typical response at the global scale for Case-1 with ${\sigma }_n$  = 150 MPa is shown in Figure 9(a) in the form of a contour plot for the horizontal stress during different stages of loading. Figure 9(b) shows the corresponding cohesive zone development from RVE-1 as the loading continues. The global scale response is affected by the state of RVE in each time step, as the homogenized RVE effective properties are updated at each time step through the two-way coupled multiscaling. As damage within the local RVE microstructure progresses through the adaptive insertion of cohesive zone elements and their evolution, the homogenized RVE damage influences on the global scale response as depicted in Figure 9(a) where a significant softening behavior can be observed which results in reduced load carrying capacity.

Figure 9.

Example multiscale model simulation results (Case-1 with ${\sigma }_n$ = 150 MPa): (a) global scale stress response during different stages of loading and (b) corresponding RVEs with cohesive zone evolution at progressive loading times.

As described in equation (9), the parameter ${\sigma }_n$ depicts the critical state of stress at which the Mode-I fracture process is initiated in the interstitial matrix within the RVE. The range of the parameter ${\sigma }_n$ selected for Case-1 (i.e., 150 to 400 MPa) was based on the global scale experiments found in the literature (Arias-Moreno et al., 2020; Mustafy et al., 2019; Ritchie et al., 2008). Clearly, the parameter ${\sigma }_n$ influences the onset of damage accumulation in the global scale as observed in Figure 10(a), where significant deviation from the undamaged state occurs close to the critical stress state.

Figure 10.

Parametric analysis results for RVE-1 at a loading rate of 0.0667 mm/s: (a) the effects of NVCZ parameter ${\sigma }_n$ on the global stress-strain response and (b) the effects of ${\delta }_n$ on the global stress-strain response.

The influence of the parameter ${\delta }_n$ is shown in Figure 10(b). Usually, ${\delta }_n$ is a parameter that is associated with the critical crack opening displacement before complete failure, but one can see that the parameter also affects the rate at which damage within the RVE progresses. As shown in Figure 10(b), one can observe that decreasing the value of ${\delta }_n(<1)$ can increase the rate of damage evolution and change the fracture process from more compliant to a brittle fashion.

Figure 11 shows the influence of the four NVCZ damage parameters ( $A, m, \overline{\lambda }$ and $\mathrm{\Delta }$ ) that control the rate of crack evolution as described in equation (11). It appears that all damage parameters significantly influence the overall behavior of bone fracture and damage. Several different values were applied in order to investigate the sensitivity of the damage parameters for microscale crack growth and eventual macroscale failure of the bone. The microscale damage evolution rate is directly proportional to higher values of the microscale damage evolution coefficient, A. On the other hand, the microscale damage evolution rate increases as the microscale damage evolution exponent, m decreases, since ${\lambda }^L\left(t\right)$ in equation (10) is typically less than 1.0 before complete failure of the cohesive zone. The factor $\mathrm{\beta }\left(\lambda \right)$ shown in equation (12) is the Gaussian distribution function and it is related to two other parameters: $\overline{\lambda }$ and $\mathrm{\Delta }$ . The effects of the parameters are presented in Figure 11(c) and Figure 11(d), respectively.

Figure 11.

Parametric analysis results for RVE-1 at a loading rate of 0.0667 mm/s: (a) the effects of NVCZ parameter $\mathit{m}$ on the global stress-strain response; (b) the effects of $\mathit{A}$ on the global stress-strain response; (c) the effects of $\overline{\mathit{\lambda }}$ on the global stress-strain response and (d) the effects of $\mathbf{\Delta }$ on the global stress-strain response.

The capability of the NVCZ model to capture the rate-dependent fracture process within the local scale RVE was assessed by applying three different displacement loading rates on the global object. The three displacement rates applied along with the NVCZ model parameters are shown in Table 4. The resulting stress-strain plots are shown in Figure 12. It can be observed that as the loading rates increased, the stress within the RVE also increased.

Table 4.

Loading rates and NVCZ parameters used for case-7.

		Case-7
Name	Rate (mm/s)	NVCZ parameters
R1	0.0667	${\sigma }_t$ [MPa]	2E + 05
		${\delta }_t$ [mm]	1.00
R2	0.2	${\sigma }_n$ [MPa]	200
		${\delta }_n$ [mm]	1.00
		A	94.83
R3	5.0	m	1.00
		$\overline{\lambda }$	0.005
		$\mathrm{\Delta }$	0.005

Figure 12.

Parametric analysis results for RVE-1 at varying loading rates: (a) R1 = 0.0667 mm/s; (b) R2 = 0.2 mm/s and (c) R3 = 5 mm/s.

Summary and conclusions

This study presented multiscale modeling of viscoelastic deformation and rate-dependent fracture in cortical bone of rat femur. Toward that end, we characterized the viscoelastic properties of bone using the nanoindentation test and its finite element modeling incorporated with an inverse optimization method. Once the constitutive properties were identified through the inverse optimization approach, the local properties were incorporated within a two-way coupled multiscale modeling framework where a global scale bone was linked to its local scale heterogeneous RVE. Fracture and damage growth within the cortical bone microstructure and its influence on the global scale response were considered through an NVCZ model. A parametric study was conducted by varying the fracture properties of the NVCZ model. The following conclusions can be drawn.

It was observed that bone shows viscoelastic behavior at the tissue level, while solely LVE constitutive relation seems limited to fully characterize the bone response, which was demonstrated by a much larger recovery in displacements observed from the model simulation compared to test results. Constitutive models, such as viscoelastic-plastic models, that can better capture larger permanent deformation should be considered.

This study demonstrated that small-scale material property characterization (such as nanoindentation testing of bone microstructure) can be effectively used to provide core material properties to conduct multiscale modeling-analysis of heterogeneous bone that presents complex inelastic, rate-dependent constitutive behavior.

Rate-dependent damage growth within the bone microstructure can be incorporated by using a model such as the NVCZ model herein. A parametric study on the fracture parameters revealed that different fracture processes can be effectively captured through a set of parameters.

The model successfully demonstrated the capability of representing viscoelastic deformation and nonlinear fracture damage that are linked between length scales. This further implies that the viscoelastic fracture model and the two-way coupled modeling can elucidate the fundamental understanding of complex multiscale deformation and fracture of bone, along with potentially advancing the clinical practices related to bone fracture risk management and prevention.

At the current stage of modeling, several limitations of the model were inevitably allowed because of technical challenges within a limited work scope. Model validation is one of them due to the lack of direct experimental data. Other than cortical bone, cancellous bone, which presents very different microstructural characteristics compared to cortical bone tissue, can be modeled to encompass the entire bone tissue in the analysis process. Other species (e.g., bovine, ovine, even human) presenting more complex bone tissue microstructures can be modeled with the proper RVE and microstructural characteristics identified. The model validation and several advancements necessary to the challenges are currently in progress by the authors to improve modeling capability and validity. New results will be presented in future studies.