Computational Investigation of the Mechanical Behavior of a Bone-Inspired Nanocomposite Material

Abstract

Bioinspired nanocomposites aim to mimic the structure of natural materials. These materials exhibit excellent mechanical properties such as high strength, toughness, and stiffness. Using modeling and simulation, we can gain insight into the underlying mechanisms that control the properties of these materials, study the impact of various parameters on their performance, and design new materials with high performance. This study investigates a bone-inspired nanocomposite that consists of two subunits: Subunit-A (Mineralized Collagen Fibril) and Subunit-B (Extrafibrillar Matrix). Subunit-B provides the composite with stiffness before yielding. After yielding, Subunit-A stretches to accommodate the deformation up to the final failure. The adhesive material in the interface plays an important role in this nanocomposite’s failure. The composite’s toughness is enhanced by multiple mechanisms: diffuse damage in Subunit-B, strain relaxation around crack tips through horizontal interface delamination between the subunits, and the crack bridging role of Subunit-A. This study provides insight into the mechanical behavior of bone-inspired nanocomposites under tensile loading conditions, highlighting the importance of the adhesive phase in optimizing the material performance in various applications.

1. Introduction

In nature, there are many examples of tough, strong, lightweight, and multifunctional structures, which can be an inspiration for the design of high performance materials, such as nacre and bone [1,2,3,4,5]. Bone is a complex biological material that serves as the primary structural support for the human body [6]. It is a composite material consisting of mineralized collagen fibrils and non-collagenous proteins, and its mechanical properties are governed by its hierarchical structure spanning from the macro- down to nano-scale [7,8,9,10]. This hierarchical structure is essential for the unique combination of stiffness, strength, and toughness that bones exhibit. Understanding the hierarchical structure of bone is critical to the development of biomimetic materials that can mimic the mechanical properties of bone. By mimicking the structure of bone, it is possible to design new materials that can be used in a wide range of applications in industry [11,12,13,14,15,16,17]. The hierarchical structure of bone is shown in Figure 1. At the macroscopic scale, bones have a well-defined shape and structure, with a compact portion called cortical bone and spongy portion called cancellous bone [18,19,20]. At the mesoscale, the basic building blocks of cortical bone are osteons, also known as Haversian systems, which are the basic functional unit of compact bone tissue and consist of concentric layers of bone tissue arranged around a central canal. These concentric layers are known as lamellae. At the nanoscale, each lamella is composed of Mineralized Collagen Fibril (MCF) and Extrafibrillar Matrix (EFM) [21,22,23,24]. The MCF is a complex structure composed of collagen and hard hydroxyapatite (HA) mineral platelets. The EFM consists of hard HA crystals bonded through a thin layer of Non-Collagenous Proteins (NCPs), which can be considered as a tough bioadhesive [25]. The layer of NCPs consists of Osteocalcin (OC), Osteopontin (OPN), H${}_2$O, and Proteoglycans (PGs). The MCF and EFM make up the building blocks of bone tissue, providing it with the unique combination of strength, toughness, and resilience required for its many functions in the body. Understanding the composition and structure of lamellae is crucial to study the mechanics and properties of bone.

Bioinspired nanocomposites are a class of materials that aim to mimic the structure and properties of natural materials [29,30,31]. These materials have been found to exhibit excellent mechanical properties, including high strength, toughness, and stiffness. Due to their unique properties and potential applications, biomimetic materials have attracted extensive attention in recent years. Computational modeling and simulation techniques are critical tools in the design and development of these materials. The interfaces in these materials play a significant role in determining their material properties. The cohesive zone model (CZM) is widely used to describe and analyze the damage mode of material interface during the failure process [32,33,34,35,36,37]. The CZM has proven to be highly effective in capturing the intricate interactions between different material phases and predicting the initiation and propagation of cracks, debonding, and other failure mechanisms. Using modeling and simulation, we can gain insight into the underlying mechanisms that control the properties of bioinspired nanocomposites, study the impact of various parameters on the performance, and design new materials with higher performance, leading to new and exciting applications in many fields.

2. Computational Model

2.1. Geometric Model

In this study, we develop a two-dimensional (2D) computational model of bone-inspired nanocomposite material. The structural model of the lamella-like nanocomposite consists of two different subunits (Figure 2), Subunit-A (MCF) and Subunit-B (EFM), which are aligned in the horizontal direction and connected by a thin layer of glue-like interface material. Subunit-A is a soft matrix that contains stiff plates with a size of $40\mathrm{n}\mathrm{m}\times 5\mathrm{n}\mathrm{m}$, while Subunit-B is composed of 200 randomly distributed stiff grains via Voronoi tessellation and connected by a glue-like interface material. According to the structural description of bone tissue from previous studies [21,22,23,24], the MCF in the bone structure is modeled as a periodic arrangement. The Subunit-A mimics the MCF in bone, which exhibits a periodic arrangement. Subunit-B resembles the EFM of bone tissue, which is modeled as a random aggregate of polycrystals of a stiff material (similar to hydroxyapatite grains) bounded through a thin adhesive substance. The height of Subunit-A is 150 nm and the height of Subunit-B is 140 nm, while both subunits have the same length of 1750 nm. The total dimension of the composite structure is $1750\mathrm{n}\mathrm{m}\times 1620\mathrm{n}\mathrm{m}$.

The structural model of the proposed nanocomposite is shown in Figure 2. The composite material is studied under uniaxial tensile loading in the X-axis direction. The right and left edges of the nanocomposite were subjected to horizontal tensile displacement boundary conditions, while being constrained for vertical movement. Throughout the loading step, the total kinetic energy of the model was monitored and compared with respect to the total strain energy of the model. The ratio was less than 5%, ensuring an acceptable state of a quasi-static loading condition. The model verification and a convergence study on nanocomposite size and mesh size were conducted in our previous studies in composite materials with a similar structure [16,17,24,38].

2.2. Material Property

The bone-inspired nanocomposite material consists of a soft matrix, stiff plates, stiff grains, and adhesive. The volume fraction of each component can be found in

Table 1. Volume fraction of different constituents of the nanocomposite [16].

Components	Volume Fraction (%)
Soft matrix	49.5
Stiff platelets	6.1
Stiff grains	38.6
Adhesive	5.8

. The material parameters used for the modeling matrix, platelets, and grains are listed in

Table 2. Material parameters used for modeling matrix, platelets, and grains in base model [40,41,42].

Material Property	Soft Matrix	Stiff Platelets	Stiff Grains
Young’s modulus (E) (GPa)	3	100	100
Poission ratio	0.4	0.3	0.3
Failure strain ${\epsilon }_f$ (%)	15	-	-

and used to accurately simulate the mechanical properties and behavior of the material under mechanical loading conditions. Polymethyl methacrylate (PMMA) is by far the most frequently used bone substitute material for vertebroplasty [39]. The soft matrix properties were set to the following: Young’s modulus $E=3$ GPa, Poission ratio $\nu =0.4$, and failure strain ${\epsilon }_f=15\%$ [40]. According to the material properties of hard HA mineral platelet crystals reported in the literature [41,42], the material parameters of stiff plates and stiff grains were set to the following: Young’s modulus $E=100$ GPa and Poission ratio $\nu =0.3$. By tuning the relevant material parameters, the computational model is able to replicate the mechanical properties and behavior of the bone-inspired nanocomposite materials. The detailed characterization of the material properties and the accurate representation of the different components in the model can lead to a better understanding of the performance of the bone-inspired nanocomposite material and how it can be optimized for specific applications.

2.3. Cohesive Zone Model for Material Interface Modeling

In this study, a finite element (FE) model of the bone-inspired nanocomposite was created with 142,000 elements. Those elements are comprised of 126,000 linear triangle elements, 5000 linear quadrilateral elements, and 11,000 user-defined linear cohesive elements. The geometry and FE model were created by ABAQUS software. A generalized cohesive zone model was implemented in a VUEL subroutine (user-defined element) in the ABAQUS package. A central difference dynamic explicit solution scheme was chosen over a Newton–Raphson implicit scheme mainly because of its capability to ensure a converged solution, especially when material softening occurs in the system (in our case, at cohesive elements). The increment time step was selected as ${\delta }_t=2\ast {10}^{-12}$ s, and mass scaling (adjusting the density of elements) was activated during the step to ensure that the minimum stable time increment of all elements was greater than the specified ${\delta }_t$. The average element size was found to be 7 nm from mesh sensitivity analysis. The finite element mesh could capture the complex microstructure of the nanocomposite material, including its different components and their interactions. To provide a clear understanding of the mesh structure, a zoomed-in view of the finite element mesh is presented in Figure 3. Using a large number of elements and a detailed mesh structure in the finite element model allowed us to make sure of the high-level accuracy of the simulation results. This in turn provided a more comprehensive understanding of the mechanical behavior of the bone-inspired nanocomposite material.

To describe the adhesive material deformation behavior, a cohesive zone model (CZM) was proposed to mimic various interface behaviors. Figure 4 shows the traction–separation relations in normal and tangential directions, which were used to govern the adhesive material behaviors.

The cohesive zone model is widely used to study the composite material failure process [38,43,44,45,46,47,48,49,50,51]. The adhesive tractions in the surface normal $T_n$ and tangential direction $T_t$ take the following form:

$$T_n=\left\{\begin{array}{cc}{\sigma }_c\left(\frac{{\Delta }_n-{\delta }_0}{{\delta }_{dn}-{\delta }_0}\right)exp\left(1\right)exp\left(-\frac{{\Delta }_n-{\delta }_0}{{\delta }_{dn}-{\delta }_0}\right)\hfill & {\Delta }_n\le {\delta }_{dn}\\ {\sigma }_c{\left(\frac{{\delta }_{fn}-{\Delta }_n}{{\delta }_{fn}-{\delta }_n}\right)}^{p_n}\hfill & {\delta }_{dn}<{\Delta }_n<{\delta }_{fn}\\ 0\hfill & {\Delta }_n\ge {\delta }_{fn}\end{array}\right.$$

(1)

$$T_t=\left\{\begin{array}{cc}{\tau }_c\sqrt{exp\left(1\right)}\left(\frac{{\Delta }_t}{{\delta }_{dt}}\right)exp\left(-\frac{{\Delta }_t^2}{2{\delta }_{dt}^2}\right)\hfill & \hfill 0\le \left|{\Delta }_t\right|\le {\delta }_{dt}\\ {\tau }_c\frac{{\Delta }_t}{\left|{\Delta }_t\right|}{\left(\frac{{\delta }_{ft}-\left|{\Delta }_t\right|}{{\delta }_{ft}-{\delta }_{dt}}\right)}^{p_t}\hfill & \hfill {\delta }_{dt}<\left|{\Delta }_t\right|<{\delta }_{ft}\\ 0\hfill & \hfill {\Delta }_n\ge {\delta }_{ft}\end{array}\right.$$

(2)

where ${\sigma }_c$ and ${\tau }_c$ are the cohesive strength in the normal and tangential directions, respectively, and ${\delta }_{dn}$ and ${\delta }_{dt}$ are critical separation. ${\delta }_{fn}$ and ${\delta }_{ft}$ are normal and tangential final failure separation. ${\delta }_0$ is the equilibrium position and $exp\left(1\right)=2.71828$. The shape parameters ($p_n,p_t$) are introduced to describe different adhesive material behaviors (e.g., brittle, ductile, elastic–perfectly plastic, rubber-like behavior, etc). The variables ${\Delta }_n$ and ${\Delta }_t$ represent normal and tangential separation, respectively. The main advantage of the model is that one can independently control the interface properties (e.g., stiffness, strength, toughness). The stiffness ($K_n$, $K_t$) of the adhesive at the equilibrium position is defined as

$$K_n={\left.\frac{d{T}_n}{d{\Delta }_n}\right|}_{{\Delta }_n={\delta }_0}=\frac{{\sigma }_c{e}^{q_n}}{{\delta }_{dn}-{\delta }_0}$$

(3)

$$K_t={\left.\frac{d{T}_t}{d{\Delta }_t}\right|}_{{\Delta }_t=0}=\frac{{\tau }_c{e}^{\frac{q_t}{2}}}{{\delta }_{dt}}$$

(4)

The adhesive toughness (${⌀}_n$,${⌀}_t$) can be determined from the integrated area under the traction–separation curve. They can be defined as

$${⌀}_n={\int }_{{\delta }_0}^{{\delta }_{\mathrm{fn}}}\mathbf{T}\mathrm{d}\delta$$

(5)

$${⌀}_t={\int }_0^{{\delta }_{\mathrm{ft}}}\mathbf{T}\mathrm{d}\delta$$

(6)

The adhesive material behaviors were governed by the cohesive zone model. The cohesive zone parameters can be obtained from experimental measurements (e.g., double cantilever beam), molecular dynamics, or other numerical simulations [52]. The cohesive zone parameters used in our model were estimated based on the experimental data or numerical simulation results reported in the literature [41,53,54,55]. Cohesive zone parameters used for modeling adhesive material in the computational model are listed in

Table 3. Cohesive zone parameters used for modeling adhesive material in the computational model [41,53,54,55].

${\mathit{\sigma }}_{\mathit{c}}\left(\mathbf{MPa}\right)$	${\mathit{K}}_{\mathit{n}}\left(\mathbf{MPa}\right)$	${\mathit{⌀}}_{\mathit{n}}\left(\mathbf{J}/\mathbf{m}{}^2\right)$	${\mathit{p}}_{\mathit{n}}$	${\mathit{\tau }}_{\mathit{c}}\left(\mathbf{MPa}\right)$	${\mathit{K}}_{\mathit{t}}\left(\mathbf{MPa}\right)$	${\mathit{⌀}}_{\mathit{t}}\left(\mathbf{J}/\mathbf{m}{}^2\right)$	${\mathit{p}}_{\mathit{t}}$
200	3	600	1	200	3	600	1

.

3. Results and Discussion

The average stress versus applied strain for the bone-inspired nanocomposite is plotted in Figure 5. The initial stiffness of the composite $E_{initial}$ is 27.5 GPa and yield stress ${\sigma }_y$ was found to be 100 MPa. After yielding, the bulk stress proceeded with a gradual increase with a much lower slope and reached the ultimate strength, which may indicate a strain hardening process. The post-yield stiffness $E_{post-yield}$ is 3.13 GPa, and the ultimate stress ${\sigma }_{uts}$ is 232 MPa. The failure strain ${\epsilon }_f$ is 4.7%.

Before the yielding point, both subunits resist the applied strain in a cooperative manner, but Subunit-B is the primary load carrying unit and is stressed more than Subunit-A. Figure 6a also shows that Subunit-Bs have the highest tensile stress regions of the nanocomposite in the pre-yield stage. In the post-yield stage of deformation, Subunit-A begins to increasingly carry a greater portion of the load and contributes to the bulk mechanical response of the composite.

To investigate the mechanical performance of the nanocomposite, diffuse damage and crack evolution are two important factors. The nanocomposite damage showed in the figure is based on ductile damage modelling (DUCTCRT). FE simulations captured the distribution and evolution of damage in the model (Figure 7). The red color represents adhesive regions that have reached the critical cohesive strength (${\delta }_{dn}<{\Delta }_n<{\delta }_{fn}$). Before yielding, the averaged stress in Subunit-B was found to be higher than that in Subunit-A. Once the stress in Subunit-B reached the critical strength of the adhesive phase, yielding ensued in the bulk stress–strain curve. The adhesive continued softening and experienced failure until the bulk stress reached its ultimate strength. Thereafter, stress quickly dropped to zero.

The crack propagation pattern in the post-yield stage is displayed in Figure 8. When ${\epsilon }_a$ reached 3.5%, many interface regions in Subunit-B were under softening. In this stage, Subunit-B was the major load-bearing component, which is the reason why it was particularly vulnerable to damage. When ${\epsilon }_f$ reached 4%, damage began in Subunit-B and created many cracks to dissipate energy. As the stress continued to increase in the soft matrix, it eventually reached its critical strength, and the cracks began to propagate into the soft matrix, leading to the formation of the final failure crack. The final failure process is observed at a range of applied strains between 4.2% and 4.5%.

Through the analysis of the fracture pattern in the nanocomposite, we found that the cracks were initiated at the adhesive perpendicular to the loading direction when the model was loaded in tension (Figure 7 and Figure 8). When the damage at the adhesive began to be initiated, the stress was released within the damaged cohesive zone region. As the loading kept increasing, cracks nucleated and coalesced and finally traveled vertically, passing through the nanocomposite (Figure 8). The interface sliding (Mode II) occurred along the grain boundaries in Subunit-B. It seems that there is a tendency for crack initiation, nucleation and coalescence mainly in mode I.

The stress–strain curve for the bulk nanocomposite and its two individual subunits A and B is plotted in Figure 9. Upon comparison of the three stress–strain curves, it is observed that at the initial stage of loading, Subunit-B carried the majority of the load and provided a high stiffness in bulk. Once the average stress in Subunit-B reached its critical strength, yielding occurred in the bulk stress–strain curve. During the post-yield deformation, the bulk stress continued to increase with a much lower slope. Once diffuse damage and cracks happened in Subunit-B, Subunit-A was observed to carry a load.

This analysis shows the shift in load-bearing responsibilities in the composite. Subunit-B made a significant contribution in the pre-yield loading regime. It is clear that high stresses were present in Subunit-B. As the loading progressed, Subunit-B started to experience damage, and cracks grew. The soft matrix in Subunit-A began softening, and the stiff platelets in Subunit-A then experienced high stresses. This shift in load-bearing responsibilities highlights the importance of understanding the behavior of individual subunits in a composite material and their contributions to the overall behavior of the composite.

To investigate the stiff plate effect, Figure 10 shows the strain contour at maximum applied strain. It is notable that the high strain concentration (${\epsilon }_{xx}\approx 24\%$) around the crack tips eventually initiated the failure of the soft matrix. An inhomogeneous stress distribution in Subunit-A created a pattern of low and high strain bands. Specifically, the low strain band had stiff plates, which experienced high level stress compared to the high strain band. The low strain band effectively absorbed energy and slowed down the material failure.

To study the cohesive zone parameter effect on the composite, we designed a case to study the material response under different adhesive strengths. Figure 11 shows different strain–stress curves under different adhesive strengths. The snapshots were captured at the final failure stage. L.1 shows the composite deformation under ${\sigma }_c={\tau }_c=100$ Mpa, L.2 shows the composite deformation under ${\sigma }_c={\tau }_c=200$ Mpa, and L.3 shows the composite deformation under ${\sigma }_c={\tau }_c=300$ Mpa. The parametric study revealed that the behavior of the composite material was strongly influenced by the strength of the interface. Specifically, the results indicate that a weaker interface in the composite led to more crack formation, increasing the composite’s ultimate stress and toughness. One might expect a stronger interface to result in a stronger composite material; however, a strong interface can also cause stress concentrations, which can lead to the premature failure of the material. In contrast, a weaker interface allows for more crack formation and redistribution of stress, which can increase the overall toughness of the bone-inspired nanocomposite.

It is important to note that the behavior of bone-inspired composite materials is complex and can be influenced by a variety of factors, including composition, processing conditions, and external loads. Nevertheless, the findings from this study provide valuable insights into the role of adhesive parameters on the behavior of composite materials. These insights can inform the design and optimization of composite structures for various applications.

4. Conclusions

This bone-inspired nanocomposite is made up of two subunits. Before yielding, Subunit-B, which comprises 86% of the total stiff phase content, provides the composite with stiffness. After yielding, Subunit-A stretches to accommodate the composite deformation up to 5%. Stiff platelets in Subunit-A also contribute to the post-yield stiffness of the composite. The strength of the adhesive material plays a significant role in the yield stress of the composite. Additionally, the failure strain of the soft matrix is a prime contributor to the ultimate stress of the composite. The toughness of the composite material is enhanced by several mechanisms:

• Energy dissipation through diffuse damage in Subunit-B,
• Strain relaxation around crack tips through horizontal interface delamination between the subunits, and
• The crack bridging role of Subunit-A.

The failure mechanism of the composite material is a multi-stage process. The final failure begins with damage nucleation in Subunit-A from the tips of the stiff platelets. This leads to the formation of a crack in Subunit-A. As the crack grows, it can coalesce with other cracks in the subunits and eventually form a crack across the entire composite. This can ultimately result in the failure of the material. Overall, this study provides insight into the mechanical behavior of the bone-inspired nanocomposite material under tensile loading conditions. The results highlight the importance of the different components and their interactions in determining the overall mechanical performance of the material. Understanding the stiffness, strengthening, toughening, and failure mechanisms of composite materials is crucial for their design and optimization for various applications.