Sage Journals: Discover world-class research

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.

Keywords

Computational modeling rat bone viscoelasticity fracture cohesive zone

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.

σ_{j i, j}^{G} (x_{m}, t) = 0, x_{m} \in V^{G}

(1)

σ_{i j}^{G} (x_{m}, t) = σ_{j i}^{G} (x_{m}, t), x_{m} \in V^{G}

(2)

ε_{i j}^{G} (x_{m}, t) = \frac{1}{2} (u_{i, j}^{G} (x_{m}, t) + u_{j, i}^{G} (x_{m}, t)), x_{m} \in V^{G}

(3)

σ_{i j}^{G} (x_{m}, t) = {\bar{Ω}}_{τ = \infty}^{τ = t} \{ε_{k l}^{G} (x_{m}, t)\}

(4)

where

σ_{i j}^{G} (x_{m}, t), ε_{i j}^{G} (x_{m}, t)

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,

{\bar{Ω}}_{τ = \infty}^{τ = 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.

σ_{j i, j}^{L} (x_{m}, t) = 0, x_{m} \in V^{L}

(5)

σ_{i j}^{L} (x_{m}, t) = σ_{j i}^{L} (x_{m}, t), x_{m} \in V^{L}

(6)

ε_{i j}^{L} (x_{m}, t) = \frac{1}{2} (u_{i, j}^{L} (x_{m}, t) + u_{j, i}^{L} (x_{m}, t)), x_{m} \in V^{L}

(7)

σ_{i j}^{L} (x_{m}, t) = Ω_{τ = \infty}^{τ = t} \{ε_{i j}^{L} (x_{m}, t)\}, x_{m} \in V^{L}

(8)

where

σ_{i j}^{L} (x_{m}, t), ε_{i j}^{L} (x_{m}, t)

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,

Ω_{τ = \infty}^{τ = 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} (t) = \frac{u_{i}^{L} (t)}{λ^{L} (t) δ_{i}^{L}} [1 - α^{L} (t)] [σ_{i}^{L} + \int_{t_{o}}^{t} E_{C Z}^{L} (t - τ) \frac{\partial λ^{L} (t)}{\partial τ} d τ] o n \partial V_{C Z}^{L}

(9)

λ^{L} (t) = \sqrt{{(\frac{u_{n}^{L}}{δ_{n}^{L}})}^{2} + {(\frac{u_{t}^{L}}{δ_{t}^{L}})}^{2}}

(10)

where

T_{i}^{L} (t)

is the CZ traction;

u_{i}^{L} (t)

is the CZ opening displacement;

δ_{i}^{L}

is the CZ length scale parameter; the subscript

i

is the intrinsic normal or tangential directions within the CZ;

α^{L} (t)

describes the internal damage evolution;

σ_{i}^{L}

represents the critical stress state for the CZ to be activated;

E_{C Z}^{L} (t)

is the linear viscoelastic relaxation modulus;

λ^{L} (t)

is Euclidian norm of the CZ opening displacements. A rate-dependent evolution of the internal damage parameter,

α^{L} (t)

is considered in this study as shown in equations (11) and (12).

\{\begin{array}{l} \dot{α} (λ^{L}, {\dot{λ}}^{L}) = β (λ^{L}) {[{\dot{λ}}^{L}]}^{m}, when {\dot{λ}}^{L} (t) > 0 and α^{L} (t) < 1 \\ \dot{α} (λ^{L}, {\dot{λ}}^{L}) = 0, when {\dot{λ}}^{L} (t) \leq 0 or α^{L} (t) = 1 \end{array}

(11)

β (λ^{L}) = A [\exp (\frac{- {(λ^{L} - \bar{λ})}^{2}}{2 Δ^{2}})]

(12)

where,

\dot{α} (λ^{L}, {\dot{λ}}^{L})

is the time derivative of the damage parameter which depends on the norm of the opening displacement

(λ^{L})

and its time derivative

({\dot{λ}}^{L})

. As shown in equation (12), the factor

β (λ^{L})

has a Gaussian distribution with parameters A,

\bar{λ}

and

Δ

. 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.,

σ_{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_{C Z}^{L} (t

) and rate-dependent damage evolution function. When the damage parameter

α^{L} (t)

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} = σ_{i j}^{o} n_{j}

(13)

u_{i} = ε_{i j}^{o} x_{j}

(14)

where

T_{i}

is the uniform traction;

σ_{i j}^{o}

is the constant stress tensor;

n_{j}

is the unit outward normal to the surface boundary of the RVE;

u_{i}

is the resultant displacement;

ε_{i j}^{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):

σ_{i j}^{G} = {\bar{σ}}_{i j}^{L} = \frac{1}{V^{L}} \int_{\partial V_{E}^{L}}^{} σ_{k i}^{L} n_{k}^{L} n_{j}^{L} d S

(15)

ε_{i j}^{G} = {\bar{ε}}_{i j}^{L} = {\bar{E}}_{i j}^{L} + {\bar{e}}_{i j}^{L}

(16)

{\bar{E}}_{i j}^{L} \equiv \frac{1}{V^{L}} \int_{\partial V_{E}^{L}}^{} \frac{1}{2} (u_{i}^{L} n_{j}^{L} + u_{j}^{L} n_{i}^{L}) d S

(17)

{\bar{e}}_{i j}^{L} \equiv \frac{1}{V^{L}} \int_{\partial V_{I}^{L}}^{} \frac{1}{2} (u_{i}^{L} n_{j}^{L} + u_{j}^{L} n_{i}^{L}) d S

(18)

where

{\bar{σ}}_{i j}^{L}

and

{\bar{ε}}_{i j}^{L}

is the volume average of the stresses and strains over the RVE boundary,

{\bar{E}}_{i j}^{L}

and

{\bar{e}}_{i j}^{L}

are the volume average strains over the external

(\partial V_{E}^{L})

and internal

(\partial V_{I}^{L})

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):

σ_{i j} (t) = \int_{0}^{t} G (t - τ) \frac{\partial ε_{i j}}{\partial τ} d τ + \int_{0}^{t} [K (t - τ) - \frac{2}{3} G (t - τ)] δ_{i j} \frac{\partial ε_{k k}}{\partial τ} d τ

(19)

G (t) = G_{\infty} + \sum_{i = 1}^{i = N} G_{i} e^{(\frac{- t}{ρ_{i}})}

(20)

K (t) = K_{\infty} + \sum_{i = 1}^{i = N} K_{i} e^{(\frac{- t}{ρ_{i}})}

(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

ρ_{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).

{\min_{χ \in R^{N}} Ψ_{u} (χ) = ‖u - u_{FEM} (χ)‖}_{2}

(22)

where,

Ψ_{u} (χ)

is the objective function,

χ

is a vector of size

N

material properties,

u_{FEM} (χ)

is tip displacement results from the finite element model simulation,

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
No. of Terms	$G_{\infty}$ (MPa)		$G_{\infty}$ (MPa)		$G_{\infty}$ (MPa)
1	3362.89		3299.07		3303.46
	$G_{i}$ (MPa)	$ρ_{i}$ (s)	$G_{i}$ (MPa)	$ρ_{i}$ (s)	$G_{i}$ (MPa)	$ρ_{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)	$ν$		$G_{\infty}$ (MPa)	$ν$
1	3362.89	0.3	1	3362.89	0.3
	$G_{i}$ (MPa)	$ρ_{i}$ (s)		$G_{i}$ (MPa)	$ρ_{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]		$ν$	E [MPa]		$ν$
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 ( $σ_{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
No.	Parameters	Case-1	Case-2	Case-3	Case-4	Case-5	Case-6
1	$σ_{t}$ [MPa]	2E + 05	2E + 05	2E + 05	2E + 05	2E + 05	2E + 05
2	$δ_{t}$ [mm]	1	1	1	1	1	1
3	$σ_{n}$ [MPa]	400; 300; 200; 150	200	200	200	200	200
4	$δ_{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	$\bar{λ}$	0.005	0.005	0.005	0.005	0.03; 0.01; 0.008; 0.001	0.05
8	$Δ$	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 $σ_{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 $σ_{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 $σ_{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 $σ_{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 $σ_{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 $σ_{n}$ on the global stress-strain response and (b) the effects of $δ_{n}$ on the global stress-strain response.

The influence of the parameter $δ_{n}$ is shown in Figure 10(b). Usually, $δ_{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 $δ_{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, \bar{λ}$ and $Δ$ ) 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 $λ^{L} (t)$ in equation (10) is typically less than 1.0 before complete failure of the cohesive zone. The factor $β (λ)$ shown in equation (12) is the Gaussian distribution function and it is related to two other parameters: $\bar{λ}$ and $Δ$ . 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 $m$ on the global stress-strain response; (b) the effects of $A$ on the global stress-strain response; (c) the effects of $\bar{λ}$ on the global stress-strain response and (d) the effects of $Δ$ 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	$σ_{t}$ [MPa]	2E + 05
R1	0.0667	$δ_{t}$ [mm]	1.00
R2	0.2	$σ_{n}$ [MPa]	200
		$δ_{n}$ [mm]	1.00
		A	94.83
R3	5.0	m	1.00
		$\bar{λ}$	0.005
		$Δ$	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.

Footnotes

Acknowledgements

The authors thank to Dr. Gabriel Nsengiyumva for helping the algorithm to conduct the inverse optimization process for this study.

Data availability statement

Some or all data, models, or code generated or used during the study are available from the corresponding author upon reasonable request.

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: Texas A&M Engineering Experiment Station (TEES) Faculty Research Fund.

ORCID iD

Yong-Rak Kim

References

Alcântara

Assis

Prada

, et al. (2019) Patient-specific bone multiscale modelling, fracture simulation and risk analysis – A survey. Materials 13(1): 106.

Allen

(2001) Homogenization principles and their application to continuum damage mechanics. Composites Science and Technology 61(15): 2223–2230.

Allen

Searcy

(2001) A micromechanical model for a viscoelastic cohesive zone. International Journal of Fracture 107(2): 159–176.

Arias-Moreno

Ito

van Rietbergen

(2020) Accuracy of beam theory for estimating bone tissue modulus and yield stress from 3-point bending tests on rat femora. Journal of Biomechanics 101: 109654.

Bagley

Torvik

(1983) A theoretical basis for the application of fractional calculus to viscoelasticity. Journal of Rheology 27(3): 201–210.

Burr

Turner

Naick

, et al. (1998) Does microdamage accumulation affect the mechanical properties of bone? Journal of Biomechanics 31(4): 337–345.

Cory

Nazarian

Entezari

, et al. (2010) Compressive axial mechanical properties of rat bone as functions of bone volume fraction, apparent density and micro-CT based mineral density. Journal of Biomechanics 43(5): 953–960.

Crowninshield

Pope

(1974) The response of compact bone in tension at various strain rates. Annals of Biomedical Engineering 2(2): 217–225.

Failla

Zingales

(2020) Advanced materials modelling via fractional calculus: Challenges and perspectives. Philosophical Transactions. Series A, Mathematical, Physical, and Engineering Sciences 378(2172): 20200050.

10.

Fish

(2017) Multiscale modeling of femur fracture. International Journal for Numerical Methods in Engineering 111(1): 3–25.

11.

Gitman

Askes

Sluys

(2007) Representative volume: existence and size determination. Engineering Fracture Mechanics 74(16): 2518–2534.

12.

Goldstein

(1987) The mechanical properties of trabecular bone: Dependence on anatomic location and function. Journal of Biomechanics 20(11–12): 1055–1061.

13.

Groetsch

Gourrier

Schwiedrzik

, et al. (2019) Compressive behaviour of uniaxially aligned individual mineralised collagen fibres at the micro-and nanoscale. Acta Biomaterialia 89: 313–329.

14.

Hashin

(1983) Analysis of composite materials – A survey. Journal of Applied Mechanics 50(3): 481–505.

15.

Hill

(1972) On constitutive macro-variables for heterogeneous solids at finite strain. Proceedings of the Royal Society of London A Mathematical and Physical Sciences 326: 131–147.

16.

Hua Qin

(2008) Macro-Micro Theory on Multifield Coupling Behavior of Heterogeneous Materials. Beijing: Higher Education Press.

17.

Iyo

Maki

Sasaki

, et al. (2004) Anisotropic viscoelastic properties of cortical bone. Journal of Biomechanics 37(9): 1433–1437.

18.

Johnson

Socrate

Boyce

(2010) A viscoelastic, viscoplastic model of cortical bone valid at low and high strain rates. Acta Biomaterialia 6(10): 4073–4080.

19.

Kim

D-G

Huja

Lee

, et al. (2010) Relationships of viscosity with contact hardness and modulus of bone matrix measured by nanoindentation. Journal of Biomechanical Engineering 132(2): 024502.

20.

Kim

D-G

Huja

Navalgund

, et al. (2013) Effect of estrogen deficiency on regional variation of a viscoelastic tissue property of bone. Journal of Biomechanics 46(1): 110–115.

21.

Kim

D-G

Jeong

Y-H

Kosel

, et al. (2015) Regional variation of bone tissue properties at the human mandibular condyle. Bone 77: 98–106.

22.

Kim

Y-R

Souza

Teixeira

JESL

(2013) A two-way coupled multiscale model for predicting damage-associated performance of asphaltic roadways. Computational Mechanics 51(2): 187–201.

23.

Demirci

Silberschmidt

(2013) Variability and anisotropy of mechanical behavior of cortical bone in tension and compression. Journal of the Mechanical Behavior of Biomedical Materials 21: 109–120.

24.

Liu

Kim

, et al. (2021) Multiscale characterization of ovariectomized rat femur. Journal of Biomechanics 122: 110462.

25.

Lloyd

Wang

Donnelly

(2015) Multiscale contribution of bone tissue material property heterogeneity to trabecular bone mechanical behavior. Journal of Biomechanical Engineering 137(1): 0108011–0108018.

26.

McElhaney

(1966) Dynamic response of bone and muscle tissue. Journal of Applied Physiology 21(4): 1231–1236.

27.

Mercer

Wang

, et al. (2006) Mechanisms governing the inelastic deformation of cortical bone and application to trabecular bone. Acta Biomaterialia 2(1): 59–68.

28.

Mustafy

Londono

Moldovan

, et al. (2019) High impact exercise improves bone microstructure and strength in growing rats. Scientific Reports 9(1): 13128–14.

29.

O'Brien

F J

Taylor

Clive Lee

(2005) The effect of bone microstructure on the initiation and growth of microcracks. Journal of Orthopaedic Research: Official Publication of the Orthopaedic Research Society 23(2): 475–480.

30.

Oliver

Pharr

(2004) Measurement of hardness and elastic modulus by instrumented indentation: Advances in understanding and refinements to methodology. Journal of Materials Research 19(1): 3–20.

31.

Podshivalov

Fischer

Bar-Yoseph

(2014) On the road to personalized medicine: Multiscale computational modeling of bone tissue. Archives of Computational Methods in Engineering 21(4): 399–479.

32.

Popov

Heß

Willert

(2019) Handbook of Contact Mechanics. Berlin: Springer.

33.

Rami

K Z

Kim

Y-R

Khedmati

, et al. (2018) Two-way linked multiscale method integrated with nanomechanical tests and cohesive zone fracture to model highly heterogeneous binding materials. Journal of Engineering Mechanics 144(10): 04018095.

34.

Rincón-Kohli

Zysset

(2009) Multi-axial mechanical properties of human trabecular bone. Biomechanics and Modeling in Mechanobiology 8(3): 195–208.

35.

Ritchie

Buehler

Hansma

(2009). Plasticity and toughness in bone. Physics Today 62(6): 41–47.

36.

Ritchie

Koester

Ionova

, et al. (2008) Measurement of the toughness of bone: A tutorial with special reference to small animal studies. Bone 43(5): 798–812.

37.

Sabet

Raeisi Najafi

Hamed

, et al. (2016) Modelling of bone fracture and strength at different length scales: A review. Interface Focus 6(1): 20150055.

38.

Seref-Ferlengez

Kennedy

Schaffler

(2015) Bone microdamage, remodeling and bone fragility: How much damage is too much damage? BoneKEy Reports 4: 644.

39.

Shymanskyi

Sokolovskyy

(2021) Finite element calculation of the linear elasticity problem for biomaterials with fractal structure. The Open Bioinformatics Journal 14(1): 114–122.

40.

Souza

Allen

(2010) Multiscale modeling of impact on heterogeneous viscoelastic solids containing evolving microcracks. International Journal for Numerical Methods in Engineering 82(4): 464–504.

41.

Sumelka

(2014) Fractional viscoplasticity. Mechanics Research Communications 56: 31–36.

42.

Sumelka

Łuczak

Gajewski

, et al. (2020) Modelling of AAA in the framework of time-fractional damage hyperelasticity. International Journal of Solids and Structures 206: 30–42.

43.

Sun

Zhang

Baleanu

, et al. (2018) A new collection of real world applications of fractional calculus in science and engineering. Communications in Nonlinear Science and Numerical Simulation 64: 213–231.

44.

Suzuki

Zhou

D’Elia

, et al. (2021) A thermodynamically consistent fractional visco-elasto-plastic model with memory-dependent damage for anomalous materials. Computer Methods in Applied Mechanics and Engineering 373: 113494.

45.

Taylor

Hazenberg

Lee

(2007) Living with cracks: damage and repair in human bone. Nature Materials 6(4): 263–268.

46.

Ural

(2020) Advanced modeling methods – Applications to bone fracture mechanics. Current Osteoporosis Reports 18(5): 568–576.

47.

Ural

Mischinski

(2013) Multiscale modeling of bone fracture using cohesive finite elements. Engineering Fracture Mechanics 103: 141–152.

48.

Ural

Vashishth

(2014) Hierarchical perspective of bone toughness–from molecules to fracture. International Materials Reviews 59(5): 245–263.

49.

Vandamme

Ulm

F-J

(2006) Viscoelastic solutions for conical indentation. International Journal of Solids and Structures 43(10): 3142–3165.

50.

Vashishth

(2007) Hierarchy of bone microdamage at multiple length scales. International Journal of Fatigue 29(6): 1024–1033.

51.

Guo

(2004) The dependence of transversely isotropic elasticity of human femoral cortical bone on porosity. Journal of Biomechanics 37: 1281–1287.

52.

Yamashita

Furman

Rawls

, et al. (2001) The use of dynamic mechanical analysis to assess the viscoelastic properties of human cortical bone. Journal of Biomedical Materials Research 58(1): 47–53.

53.

You

Kim

Park

(2017) Two-way coupled multiscale model for predicting mechanical behavior of bone subjected to viscoelastic deformation and fracture damage. Journal of Engineering Materials and Technology 139(2): 021016.

54.

Blanchard

(1996) An elastic-plastic indentation model and its solutions. Journal of Materials Research 11(9): 2358–2367.

55.

Zimmermann

Busse

Ritchie

(2015) The fracture mechanics of human bone: Influence of disease and treatment. BoneKEy Reports 4: 743.

56.

Zysset

Sonny

Hayes

(1994) Morphology-mechanical property relations in trabecular bone of the osteoarthritic proximal tibia. The Journal of Arthroplasty 9(2): 203–216.

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

Abstract

Keywords

Introduction

Modeling framework

IBVP of global and local scale

Nonlinear viscoelastic cohesive zone (NVCZ) model

Homogenization for linking local to global scale

Rat bone specimen and test results

Inverse optimization of local scale constitutive properties

Finite element modeling of nanoindentation testing

Inverse optimization process

Results from the inverse optimization approach

Two-way coupled multiscale modeling with the NVCZ

Summary and conclusions

Footnotes

Acknowledgements

Data availability statement

Declaration of conflicting interests

Funding

ORCID iD

References