单轴拉伸载荷下骨单位含径向微裂纹的密质骨力学模型

Mechanical Model of Cortical Bone with Radial Microcrack in Osteon Under Uniaxial Tensile Load

摘要

针对Haversian密质骨的力学行为，建立含有单条径向微裂纹的密质骨平面力学模型。基于线弹性断裂力学理论和位错方法，推导出骨单位应力场的数学表达式，并构建与应力强度因子相关的奇异积分方程组。数值结果表明，增大骨单位与间质骨间的剪切模量比是抑制裂纹扩展的关键，且当裂纹趋近黏合线时，该抑制效果会因几何屏蔽效应而进一步增强。研究结果可为生物医学工程中仿生材料的设计提供理论依据。

Abstract

This paper investigates the mechanical behavior of Haversian cortical bone by establishing a plane mechanical model of cortical bone containing a single radial microcrack. Based on linear elastic fracture mechanics and crystal dislocation theory， the mathematical expression of the stress field in the osteon is derived， and a set of singular integral equations related to the stress intensity factor is constructed. Through numerical examples， the influence of osteon size， shear modulus ratio， and the location and length of the microcrack on the stress intensity factor at the microcrack tip is systematically studied. The numerical results indicate that increasing the shear modulus ratio between the osteon and the interstitial bone is key to inhibiting crack propagation. Moreover， as the crack approaches the cement line， this inhibitory effect is further enhanced due to the geometric shielding effect. This study provides a theoretical basis for the design of biomimetic materials in biomedical engineering.

Graphical abstract

关键词

骨单位 / Haversian密质骨 / 微裂纹 / 位错力学理论 / 奇异积分方程

Key words

osteon / Haversian cortical bone / microcrack / dislocation theory / singular integral equation

骨骼作为人体和脊椎动物的重要支撑结构，在维持机体稳定性及运动功能中发挥着关键作用，其研究在生物医学和临床应用领域^［1-2］以及仿生材料与工程应用领域^［3-4］中一直备受关注。骨骼的力学性能与其微观结构^［5-6］密切相关，尤其在微观尺度下，骨的力学行为与损伤机理已成为骨生物力学研究的前沿课题。

自1960年Frost^［7］在活体骨中发现微裂纹并提出其与疲劳损伤密切相关以来， 学者们持续深入探索微裂纹对骨力学性能的影响。 Currey^［8］通过纳米压痕和原子力显微镜等技术进一步揭示， 微裂纹会显著降低骨的韧性与强度，从而增加骨折风险。 为系统研究微裂纹的生成与扩展机制，有限元分析等数值方法被广泛应用于模拟骨中微裂纹的力学行为、预测其扩展路径以及评估骨折风险^［9-10］。 同时，在密质骨微裂纹的理论研究方面取得了系列进展：Raeisi Najafi等^［11-12］建立了包含多条微裂纹的间质骨平面模型； Wang等^［13］采用奇异积分方程方法，分析了间质骨中微裂纹的长度、位置等因素对裂纹尖端应力强度因子的影响； Ji等^［14］通过相场模拟，研究了黏合线对微裂纹扩展的抑制作用。密质骨中微裂纹的研究多集中于间质骨区域，并已建立成熟的力学模型。与此相比，针对骨单位内部微裂纹的平面力学建模研究则尚未得到充分探索。骨裂纹研究中，宏观尺度关注整体骨结构的力学行为，通常采用四点弯曲、紧凑拉伸等实验测试方法观察断裂特性，侧重于工程应用；微观尺度聚焦密质骨微结构，采用理论模型解析单条微裂纹的扩展机制，侧重于剪切模量比、几何屏蔽效应等参数对裂纹尖端应力强度因子的影响。

基于此，论文旨在建立一种含径向微裂纹的骨单位密质骨力学模型，结合线弹性断裂力学理论与数值方法，系统探究微裂纹在骨组织中的力学行为与扩展规律，为深入理解骨损伤机理提供理论依据。

1 模型描述

密质骨主要由大量圆筒状的骨单位构成，这些骨单位之间填充着形状不规则的间质骨（图1）。文中将密质骨视为均匀弹性材料^［15］， 并将骨单位简化为均匀的圆形平面结构，假设其完全嵌入间质骨中， 且二者之间为焊接连接；同时忽略骨单位之间的相互作用^［11］。 模型中考虑一条沿径向分布于骨单位内部的微裂纹。由于微裂纹的尺度远小于间质骨的尺度，在理论分析中，可将间质骨视为无限大的平面，其平面示意图如图1（a）所示。 以骨单位中心为原点建立直角坐标系， 骨单位半径为R。模型在无穷远处承受沿y轴方向的单轴拉伸应力p的作用。 微裂纹长度为L， 其左端点为A，右端点为B。微裂纹右端点B到黏合线的距离为d。间质骨的剪切模量为G_{\mathrm{1}}， 弹性常数为{\kappa }_{\mathrm{1}}=\mathrm{3}-\mathrm{4}{v}_{\mathrm{1}}，骨单位的剪切模量为G_{\mathrm{2}}， 弹性常数为{\kappa }_{\mathrm{2}}=\mathrm{3}-\mathrm{4}{v}_{\mathrm{2}}， 其中v₁、v₂是泊松比。

2 问题求解

该模型的求解可分解为两个基本问题。 问题一如图1（b）所示：骨单位嵌入无限大间质骨中，在与图1（a）裂纹对应位置预设一条假想的虚拟裂纹。 模型在无穷远处承受y轴方向的单轴拉伸应力p。 问题二如图 1（c） 所示：模型构型与问题一相同，但虚拟裂纹变为实际裂纹，且无穷远处无载荷作用；取而代之的是， 在微裂纹面上施加一组与问题一中虚拟裂纹位置应力大小相等、方向相反的外载荷。 下文以应力函数的上标 “1” 和 “2” 分别标识问题一与问题二的解。

问题一中， 设复变量z=x+\mathrm{i}y, 根据Muskhelishvili 复变函数方法，应力分量可表示为

{\sigma }_{xx}+{\sigma }_{yy}=\mathrm{2}[{\phi }^{\text{'}}(z)+\overline{{\phi }^{\text{'}}(z)}],（1）

{\sigma }_{yy}-{\sigma }_{xx}+\mathrm{2}\mathrm{i}{\sigma }_{xy}=\mathrm{2}[\overline{z}{\phi }^{″}(z)+{\psi }^{\text{'}}(z)]\mathrm{。}（2）

其中：{\sigma }_{xx}\mathrm{、}{\sigma }_{xy}\mathrm{、}{\sigma }_{yy}为应力分量；\phi (z)和\psi (z)为待求的解析复应力函数。 在无穷远处， 当\sqrt[]{x^{\mathrm{2}}+y^{\mathrm{2}}}\to \mathrm{\infty }时， 边界条件为

{\sigma }_{xx}=\mathrm{0},{\sigma }_{yy}=p,{\sigma }_{xy}=\mathrm{0}\mathrm{。}（3）

利用平面线弹性断裂力学的幂级数解法^［16］，可得应力函数的表达式为

\phi (z)=\frac{p}{\mathrm{4}}\gamma z,\psi (z)=-\frac{p}{\mathrm{2}}\delta z\mathrm{。}（4）

其中：

\gamma =\frac{G_{\mathrm{2}}({\kappa }_{\mathrm{1}}+\mathrm{1})}{\mathrm{2}{G}_{\mathrm{2}}+G_{\mathrm{1}}({\kappa }_{\mathrm{1}}-\mathrm{1})},\delta =\frac{G_{\mathrm{2}}({\kappa }_{\mathrm{1}}+\mathrm{1})}{G_{\mathrm{1}}+G_{\mathrm{2}}{\kappa }_{\mathrm{1}}},

而G_{\mathrm{1}}\mathrm{、}{G}_{\mathrm{2}}分别为间质骨与骨单位的剪切模量，{\kappa }_{\mathrm{1}}\mathrm{、}{\kappa }_{\mathrm{2}}为对应的弹性常数。

结合式（1）~（2）与式（4），可得问题一的应力场表达式为

{\sigma }_{xx}^{\mathrm{1}}=p\frac{\gamma +\delta }{\mathrm{2}},{\sigma }_{yy}^{\mathrm{1}}=p\frac{\gamma -\delta }{\mathrm{2}},{\sigma }_{xy}^{\mathrm{1}}=\mathrm{0}\mathrm{。}（5）

根据以上公式，可计算出问题一中虚裂纹上任意一点的正应力{\sigma }_{\mathrm{n}}^{\mathrm{1}}(s)和切应力{\sigma }_{\mathrm{t}}^{\mathrm{1}}(s)。

问题二可采用位错方法进行求解。 该方法将裂纹模拟为连续分布的位错， 进而将裂纹问题转化为位错分布的积分方程求解。采用文献［17］中给出的边缘位错解作为格林函数，可建立问题二中应力扰动场的积分表达式。如图1（c）所示，考虑在骨单位微裂纹上一点s_{\mathrm{0}}\left(x_{\mathrm{0}},y_{\mathrm{0}}\right)处存在伯格斯矢量为b_x与b_y的两个刃型位错，则在骨单位内任意点P\left(x,y\right)处， 由该位错引起的应力分量{\sigma }_{xx}^{\mathrm{2}}\mathrm{、}{\sigma }_{xy}^{\mathrm{2}}\mathrm{、}{\sigma }_{yy}^{\mathrm{2}}可表示为

\frac{\mathrm{\pi }({\kappa }_{\mathrm{2}}+\mathrm{1})}{G_{\mathrm{2}}}{\sigma }_{xx}^{\mathrm{2}}(x,y,s_{\mathrm{0}})=h_{xxx}(x,y,s_{\mathrm{0}})b_x+h_{xxy}(x,y,s_{\mathrm{0}})b_y,（6）

\frac{\mathrm{\pi }({\kappa }_{\mathrm{2}}+\mathrm{1})}{G_{\mathrm{2}}}{\sigma }_{yy}^{\mathrm{2}}(x,y,s_{\mathrm{0}})=h_{yyx}(x,y,s_{\mathrm{0}})b_x+h_{yyy}(x,y,s_{\mathrm{0}})b_y,（7）

\frac{\mathrm{\pi }({\kappa }_{\mathrm{2}}+\mathrm{1})}{G_{\mathrm{2}}}{\sigma }_{xy}^{\mathrm{2}}(x,y,s_{\mathrm{0}})=h_{xyx}(x,y,s_{\mathrm{0}})b_x+h_{xyy}(x,y,s_{\mathrm{0}})b_y\mathrm{。}（8）

其中：系数函数表达式为

\begin{array}{l}{h}_{xxx}=-\mathrm{2}\left(\frac{y}{r_{\mathrm{1}}^{\mathrm{2}}}+\frac{\mathrm{2}{x}_{\mathrm{1}}^{\mathrm{2}}y}{r_{\mathrm{2}}^{\mathrm{4}}}\right)+\left(A+B\right)\left(\frac{y}{r_{\mathrm{2}}^{\mathrm{2}}}+\frac{\mathrm{2}{x}_{\mathrm{2}}^{\mathrm{2}}y}{r_{\mathrm{2}}^{\mathrm{4}}}\right)+\left(A-B\right)\left(\frac{\mathrm{2}{x}_{\mathrm{2}}^{\mathrm{2}}y}{r_{\mathrm{2}}^{\mathrm{4}}}\right)+\\ \frac{A\left(\mathrm{1}-{\beta }^{\mathrm{2}}\right)R}{{\beta }^{\mathrm{3}}}\left[\frac{\mathrm{4}{x}_{\mathrm{2}}y}{r_{\mathrm{2}}^{\mathrm{4}}}\right.-\frac{\mathrm{16}{x}_{\mathrm{2}}^{\mathrm{3}}y}{r_{\mathrm{2}}^{\mathrm{6}}}+\left.\frac{\left(\mathrm{1}-{\beta }^{\mathrm{2}}\right)R}{\beta }\left(\frac{\mathrm{2}y}{r_{\mathrm{2}}^{\mathrm{4}}}-\frac{\mathrm{8}{x}_{\mathrm{2}}^{\mathrm{2}}y}{r_{\mathrm{2}}^{\mathrm{6}}}\right)\right],\end{array}

\begin{array}{l}{h}_{xxy}=\mathrm{2}\left(\frac{\mathrm{2}{x}_{\mathrm{1}}^{\mathrm{3}}}{r_{\mathrm{1}}^{\mathrm{4}}}-\frac{x_{\mathrm{1}}}{r_{\mathrm{1}}^{\mathrm{2}}}\right)-\left(A+B\right)\left(\frac{\mathrm{2}{x}_{\mathrm{2}}^{\mathrm{3}}}{r_{\mathrm{2}}^{\mathrm{4}}}-\frac{x_{\mathrm{2}}}{r_{\mathrm{2}}^{\mathrm{2}}}\right)+\left(B-A\right)\left(\frac{\mathrm{2}{x}_{\mathrm{2}}^{\mathrm{3}}}{r_{\mathrm{2}}^{\mathrm{4}}}\right)+\mathrm{2}\overline{A}\left(\mathrm{1}+Q{\beta }^{\mathrm{2}}\right)\frac{\mathrm{1}}{\beta R}+\\ A\frac{\left(\mathrm{1}-{\beta }^{\mathrm{2}}\right)R}{{\beta }^{\mathrm{3}}}\left[\mathrm{2}{\beta }^{\mathrm{2}}\left(\frac{\mathrm{1}}{r_{\mathrm{2}}^{\mathrm{2}}}-\frac{\mathrm{2}{x}_{\mathrm{2}}^{\mathrm{2}}}{r_{\mathrm{2}}^{\mathrm{4}}}\right)\right.+\left(\frac{\mathrm{12}{x}_{\mathrm{2}}^{\mathrm{2}}}{r_{\mathrm{2}}^{\mathrm{4}}}\right.-\left.\left.\frac{\mathrm{16}{x}_{\mathrm{2}}^{\mathrm{4}}}{r_{\mathrm{2}}^{\mathrm{6}}}\right)+\frac{\left(\mathrm{1}-{\beta }^{\mathrm{2}}\right)R}{\beta }\left(\frac{\mathrm{6}{x}_{\mathrm{2}}}{r_{\mathrm{2}}^{\mathrm{4}}}-\frac{\mathrm{8}{x}_{\mathrm{2}}^{\mathrm{3}}}{r_{\mathrm{2}}^{\mathrm{6}}}\right)\right],\end{array}

\begin{array}{l}{h}_{xyx}=\mathrm{2}\left(\frac{\mathrm{2}{x}_{\mathrm{1}}^{\mathrm{3}}}{r_{\mathrm{1}}^{\mathrm{4}}}-\frac{x_{\mathrm{1}}}{r_{\mathrm{1}}^{\mathrm{2}}}\right)+\left(A+B\right)\left(\frac{x_{\mathrm{2}}}{r_{\mathrm{2}}^{\mathrm{2}}}-\frac{\mathrm{2}{x}_{\mathrm{2}}^{\mathrm{3}}}{r_{\mathrm{2}}^{\mathrm{4}}}\right)+\sqrt[]{\mathrm{2}}\left(B-A\right)\left(\frac{\mathrm{2}{x}_{\mathrm{2}}^{\mathrm{3}}}{r_{\mathrm{2}}^{\mathrm{4}}}-\frac{\mathrm{2}{x}_{\mathrm{2}}}{r_{\mathrm{2}}^{\mathrm{2}}}\right)+\\ \frac{A\left(\mathrm{1}-{\beta }^{\mathrm{2}}\right)R}{{\beta }^{\mathrm{3}}}\left[\frac{\mathrm{2}}{r_{\mathrm{2}}^{\mathrm{2}}}-\right.\left.\frac{\mathrm{16}{x}_{\mathrm{2}}^{\mathrm{2}}}{r_{\mathrm{2}}^{\mathrm{4}}}+\frac{\mathrm{16}{x}_{\mathrm{2}}^{\mathrm{4}}}{r_{\mathrm{2}}^{\mathrm{6}}}+\frac{\left(\mathrm{1}-{\beta }^{\mathrm{2}}\right)R}{\beta }\left(\frac{\mathrm{8}{x}_{\mathrm{2}}^{\mathrm{3}}}{r_{\mathrm{2}}^{\mathrm{6}}}-\frac{\mathrm{6}{x}_{\mathrm{2}}}{r_{\mathrm{2}}^{\mathrm{4}}}\right)\right],\end{array}

\begin{array}{l}{h}_{xyy}=\mathrm{2}\left(\frac{\mathrm{2}{x}_{\mathrm{1}}^{\mathrm{2}}y}{r_{\mathrm{1}}^{\mathrm{4}}}-\frac{y}{r_{\mathrm{1}}^{\mathrm{2}}}\right)+\left(A+B\right)\left(\frac{y}{r_{\mathrm{2}}^{\mathrm{2}}}-\frac{\mathrm{2}{x}_{\mathrm{2}}^{\mathrm{2}}y}{r_{\mathrm{2}}^{\mathrm{4}}}\right)+\left(B-A\right)\frac{\mathrm{2}{x}_{\mathrm{2}}^{\mathrm{2}}y}{r_{\mathrm{2}}^{\mathrm{4}}}+\\ \frac{A\left(\mathrm{1}-{\beta }^{\mathrm{2}}\right)R}{{\beta }^{\mathrm{3}}}\left[\left(\mathrm{8}-\mathrm{4}{\beta }^{\mathrm{2}}\right)\frac{x_{\mathrm{2}}y}{r_{\mathrm{2}}^{\mathrm{4}}}-\right.\left.\frac{\mathrm{16}{x}_{\mathrm{2}}^{\mathrm{3}}y}{r_{\mathrm{2}}^{\mathrm{6}}}+\frac{\left(\mathrm{1}-{\beta }^{\mathrm{2}}\right)R}{\beta }\left(\frac{\mathrm{2}y}{r_{\mathrm{2}}^{\mathrm{4}}}-\frac{\mathrm{8}{x}_{\mathrm{2}}^{\mathrm{2}}y}{r_{\mathrm{2}}^{\mathrm{6}}}\right)\right]\end{array}

\begin{array}{l}{h}_{yyx}=\mathrm{2}\left(\frac{\mathrm{2}{x}_{\mathrm{1}}^{\mathrm{2}}y}{r_{\mathrm{1}}^{\mathrm{4}}}-\frac{y}{r_{\mathrm{1}}^{\mathrm{2}}}\right)+\left(A+B\right)\left(\frac{y}{r_{\mathrm{2}}^{\mathrm{2}}}-\frac{\mathrm{2}{x}_{\mathrm{2}}^{\mathrm{2}}y}{r_{\mathrm{2}}^{\mathrm{4}}}\right)+\left(B-A\right)\left(\frac{\mathrm{2}{x}_{\mathrm{2}}^{\mathrm{2}}y}{r_{\mathrm{2}}^{\mathrm{4}}}-\frac{\mathrm{2}y}{r_{\mathrm{2}}^{\mathrm{2}}}\right)+\\ \frac{A\left(\mathrm{1}-{\beta }^{\mathrm{2}}\right)R}{{\beta }^{\mathrm{3}}}\left[\frac{\mathrm{16}{x}_{\mathrm{2}}^{\mathrm{3}}y}{r_{\mathrm{2}}^{\mathrm{6}}}-\right.\left.\frac{\mathrm{12}{x}_{\mathrm{2}}y}{r_{\mathrm{2}}^{\mathrm{4}}}+\frac{\left(\mathrm{1}-{\beta }^{\mathrm{2}}\right)R}{\beta }\left(\frac{\mathrm{8}{x}_{\mathrm{2}}^{\mathrm{2}}y}{r_{\mathrm{2}}^{\mathrm{6}}}-\frac{\mathrm{2}y}{r_{\mathrm{2}}^{\mathrm{4}}}\right)\right],\end{array}

\begin{array}{l}{h}_{yyy}=\mathrm{2}\left(\frac{\mathrm{8}{x}_{\mathrm{1}}}{r_{\mathrm{1}}^{\mathrm{2}}}-\frac{\mathrm{2}{x}_{\mathrm{1}}^{\mathrm{3}}}{r_{\mathrm{1}}^{\mathrm{4}}}\right)+\left(A+B\right)\left(\frac{\mathrm{2}{x}_{\mathrm{2}}^{\mathrm{3}}}{r_{\mathrm{2}}^{\mathrm{4}}}-\frac{\mathrm{8}{x}_{\mathrm{2}}}{r_{\mathrm{2}}^{\mathrm{2}}}\right)+\left(B-A\right)\left(\frac{\mathrm{2}{x}_{\mathrm{2}}}{r_{\mathrm{2}}}-\frac{\mathrm{2}{x}_{\mathrm{2}}^{\mathrm{3}}}{r_{\mathrm{2}}^{\mathrm{4}}}\right)+\frac{\mathrm{2}A\left(\mathrm{1}+Q{\beta }^{\mathrm{2}}\right)}{\beta R}+\\ \frac{A\left(\mathrm{1}-{\beta }^{\mathrm{2}}\right)R}{{\beta }^{\mathrm{3}}}\left[\mathrm{2}{\beta }^{\mathrm{2}}\left(\frac{\mathrm{2}{x}_{\mathrm{2}}}{r_{\mathrm{2}}^{\mathrm{4}}}-\frac{\mathrm{1}}{r_{\mathrm{2}}^{\mathrm{2}}}\right)+\frac{\mathrm{4}}{r_{\mathrm{2}}^{\mathrm{2}}}-\frac{\mathrm{20}{x}_{\mathrm{2}}^{\mathrm{2}}}{r_{\mathrm{2}}^{\mathrm{4}}}+\right.\left.\frac{\mathrm{16}{x}_{\mathrm{2}}^{\mathrm{4}}}{r_{\mathrm{2}}^{\mathrm{6}}}+\frac{\left(\mathrm{1}-{\beta }^{\mathrm{2}}\right)R}{\beta }\left(\frac{\mathrm{8}{x}_{\mathrm{2}}^{\mathrm{3}}}{r_{\mathrm{2}}^{\mathrm{6}}}-\frac{\mathrm{6}{x}_{\mathrm{2}}}{r_{\mathrm{2}}^{\mathrm{4}}}\right)\right]\mathrm{。}\end{array}

式中

A=\frac{\Gamma -\mathrm{1}}{\Gamma +{\kappa }_{\mathrm{2}}},B=\frac{\Gamma {\kappa }_{\mathrm{1}}-{\kappa }_{\mathrm{2}}}{\Gamma {\kappa }_{\mathrm{1}}+\mathrm{1}},Q=-\frac{\mathrm{2}A}{\mathrm{1}+A},\Gamma =\frac{G_{\mathrm{2}}}{G_{\mathrm{1}}},x_{\mathrm{1}}=x-x_{\mathrm{0}},x_{\mathrm{2}}=x-\frac{R^{\mathrm{2}}}{x_{\mathrm{0}}},\beta =\frac{x_{\mathrm{0}}}{R}\mathrm{。}

设{\sigma }_{\mathrm{n}}^{\mathrm{2}}和{\sigma }_{\mathrm{t}}^{\mathrm{2}}分别为微裂纹上的法向应力与切向应力。根据位错引起的应力场分布特征，法向应力与切向应力可表示为关于位错函数的代数形式。即

{\sigma }_{\mathrm{n}}^{\mathrm{2}}\left(s,s_{\mathrm{0}}\right)=h_{\mathrm{n}x}\left(s,s_{\mathrm{0}}\right)b_x+h_{\mathrm{n}y}\left(s,s_{\mathrm{0}}\right)b_y,（9）

{\sigma }_{\mathrm{t}}^{\mathrm{2}}\left(s,s_{\mathrm{0}}\right)=h_{\mathrm{t}x}\left(s,s_{\mathrm{0}}\right)b_x+h_{\mathrm{t}y}\left(s,s_{\mathrm{0}}\right)b_y。（10）

其中： h_{\mathrm{n}x}\mathrm{、}{h}_{\mathrm{n}y}\mathrm{、}{h}_{\mathrm{t}x}\mathrm{、}{h}_{\mathrm{t}y}是由h_{xxx}\mathrm{、}{h}_{xxy}\mathrm{、}{h}_{xyx}\mathrm{、}{h}_{xyy}\mathrm{、}{h}_{yyx}\mathrm{、}{h}_{yyy}等系数组合而成的解析函数。问题二的法向应力{\sigma }_{\mathrm{n}}^{\mathrm{2}}和切向应力{\sigma }_{\mathrm{t}}^{\mathrm{2}}亦可由应力分量{\sigma }_{xx}^{\mathrm{2}}\mathrm{、}{\sigma }_{xy}^{\mathrm{2}}\mathrm{、}{\sigma }_{yy}^{\mathrm{2}}表示，代入式（9）~（10），可求得系数函数h_{\mathrm{n}x}、h_{\mathrm{n}\mathrm{y}}、h_{\mathrm{t}x}、 h_{\mathrm{t}y}。

若假设位错沿微裂纹连续分布，即b_x(s_{\mathrm{0}})与b_y(s_{\mathrm{0}})在裂纹L上连续，则微裂纹面上的法向应力与切向应力可表示为

{\sigma }_{\mathrm{n}}^{\mathrm{2}}(s)={\int }_L\left[h_{\mathrm{n}x}(s,s_{\mathrm{0}})b_x(s_{\mathrm{0}})+h_{\mathrm{n}y}(s,s_{\mathrm{0}})b_y(s_{\mathrm{0}})\right]\mathrm{d}{s}_{\mathrm{0}},（11）

{\sigma }_{\mathrm{t}}^{\mathrm{2}}(s)={\int }_L\left[h_{\mathrm{t}x}(s,s_{\mathrm{0}})b_x(s_{\mathrm{0}})+h_{\mathrm{t}y}(s,s_{\mathrm{0}})b_y(s_{\mathrm{0}})\right]\mathrm{d}{s}_{\mathrm{0}}。（12）

为处理裂纹尖端的奇异性，将位错密度函数表示为

b_x\left(s_{\mathrm{0}}\right)=\frac{F_t\left(s_{\mathrm{0}}\right)}{\sqrt[]{\mathrm{1}-s_{\mathrm{0}}^{\mathrm{2}}}},b_y\left(s_{\mathrm{0}}\right)=\frac{F_n\left(s_{\mathrm{0}}\right)}{\sqrt[]{\mathrm{1}-s_{\mathrm{0}}^{\mathrm{2}}}},（13）

并满足单值性条件{\int }_L\frac{F_t\left(s_{\mathrm{0}}\right)}{\sqrt[]{\mathrm{1}-s_{\mathrm{0}}^{\mathrm{2}}}}\mathrm{d}{s}_{\mathrm{0}}=\mathrm{0},{\int }_L\frac{F_n\left(s_{\mathrm{0}}\right)}{\sqrt[]{\mathrm{1}-s_{\mathrm{0}}^{\mathrm{2}}}}\mathrm{d}{s}_{\mathrm{0}}=\mathrm{0}。

引入如下坐标变换以简化积分区间，即

x=x_{\mathrm{0}}+(\frac{L}{\mathrm{2}}\mathrm{c}\mathrm{o}\mathrm{s}{\alpha }_{\mathrm{0}})s,（14）

y=y_{\mathrm{0}}+(\frac{L}{\mathrm{2}}\mathrm{s}\mathrm{i}\mathrm{n}{\alpha }_{\mathrm{0}})s,（15）

\xi =x_{\mathrm{0}}+(\frac{L}{\mathrm{2}}\mathrm{c}\mathrm{o}\mathrm{s}\alpha )s_{\mathrm{0}},（16）

\eta =y_{\mathrm{0}}+(\frac{L}{\mathrm{2}}\mathrm{s}\mathrm{i}\mathrm{n}\alpha )s_{\mathrm{0}}。（17）

其中：L 为微裂纹长度。将式（13）~（17）代入式（11）~（12），可得积分方程为

{\sigma }_{\mathrm{n}}^{\mathrm{2}}\left(\mathrm{s}\right)={\int }_{-\mathrm{1}}^{\mathrm{1}}\left[h_{\mathrm{n}x}\frac{F_{\mathrm{t}}\left(s_{\mathrm{0}}\right)}{\sqrt[]{\mathrm{1}-s_{\mathrm{0}}^{\mathrm{2}}}}+h_{\mathrm{n}y}\frac{F_{\mathrm{n}}\left(s_{\mathrm{0}}\right)}{\sqrt[]{\mathrm{1}-s_{\mathrm{0}}^{\mathrm{2}}}}\right]\mathrm{d}{s}_{\mathrm{0}},{\sigma }_{\mathrm{t}}^{\mathrm{2}}\left(\mathrm{s}\right)={\int }_{-\mathrm{1}}^{\mathrm{1}}\left[h_{\mathrm{t}x}\frac{F_{\mathrm{t}}\left(s_{\mathrm{0}}\right)}{\sqrt[]{\mathrm{1}-s_{\mathrm{0}}^{\mathrm{2}}}}+h_{\mathrm{t}y}\frac{F_{\mathrm{n}}\left(s_{\mathrm{0}}\right)}{\sqrt[]{\mathrm{1}-s_{\mathrm{0}}^{\mathrm{2}}}}\right]\mathrm{d}{s}_{\mathrm{0}}。（18）

根据叠加原理可知，问题二与问题一在裂纹面上的应力满足

{\sigma }_{\mathrm{n}}^{\mathrm{2}}(s)=-{\sigma }_{\mathrm{n}}^{\mathrm{1}}(s),{\sigma }_{\mathrm{t}}^{\mathrm{2}}(s)=-{\sigma }_{\mathrm{t}}^{\mathrm{1}}(s)。（19）

据此，可求解未知函数F_{\mathrm{n}}、F_{\mathrm{t}}。

因此，微裂纹尖端的应力强度因子^［13］表述为

K_{\mathrm{Ⅰ}}(\pm \mathrm{1})=\pm \sqrt[]{\mathrm{\pi }l}\frac{\mathrm{2}{G}_{\mathrm{2}}}{{\kappa }_{\mathrm{1}}+\mathrm{1}}{F}_{\mathrm{n}}(\pm \mathrm{1}),K_{\mathrm{Ⅱ}}(\pm \mathrm{1})=\pm \sqrt[]{\mathrm{\pi }l}\frac{\mathrm{2}{G}_{\mathrm{2}}}{{\kappa }_{\mathrm{1}}+\mathrm{1}}{F}_{}(\pm \mathrm{1})。（20）

3 数值分析

采用高斯-切比雪夫积分方法离散积分方程组 （19），并应用广大超松弛迭代法求解^［18］，重点考察无穷远处沿y轴方向的拉伸载荷p=\mathrm{10}\mathrm{M}\mathrm{P}\mathrm{a}时，剪切模量比、裂纹长度、骨单位半径以及裂纹至黏合线的距离等因素对微裂纹尖端应力强度因子的影响。其中，考察裂纹长度、骨单位半径以及裂纹至黏合线的距离等因素时，计算中设定，间质骨的剪切模量和泊松比分别为G_{\mathrm{1}}=\mathrm{7.03}\mathrm{G}\mathrm{P}\mathrm{a}， v_{\mathrm{1}}=\mathrm{0.3}，骨单位的剪切模量和泊松比分别为 G_{\mathrm{2}}=\mathrm{8.08}\mathrm{G}\mathrm{P}\mathrm{a}，v_{\mathrm{2}}=\mathrm{0.3}。

为分析剪切模量比G_{\mathrm{2}}/G_{\mathrm{1}}对微裂纹端点应力强度因子的影响，设定骨单位半径R=\mathrm{150}\mathrm{\mu }\mathrm{m}， 微裂纹长度L=\mathrm{50}\mathrm{\mu }\mathrm{m}， 裂纹右端点B至黏合线的距离d=\mathrm{20}\mathrm{\mu }\mathrm{m}。 图2展示了微裂纹尖端标准化的SIF与剪切模量比的变化关系。结果表明，标准化应力强度因子随着G_{\mathrm{2}}/G_{\mathrm{1}}的增大呈严格单调递减趋势。 此结论表明，骨单位剪切模量高于间质骨剪切模量时，会降低裂纹尖端的应力强度因子，从而抑制微裂纹的扩展，这与文献［13］的结论相符。

为探究微裂纹尺寸效应，设定R=\mathrm{150}\mathrm{\mu }\mathrm{m}， d=\mathrm{20}\mathrm{\mu }\mathrm{m}。 图3展示了微裂纹尖端标准化的SIF与标准裂纹长度的变化关系。结果表明， 标准化应力强度因子随着标准化微裂纹长度的增加呈单调递增趋势。在裂纹较短时，其两端点的应力强度因子数值相近，但随着裂纹长度增加且端点B趋近黏合线，该端点处的应力强度因子增长速度减缓，表现出明显的屏蔽效应。将该数值结果与文献［14］对比可知，从骨单位萌生的内部微裂纹比间质骨中的外部微裂纹具有更强扩展倾向，更易引发骨单元结构破坏乃至宏观骨折。同时，与文献［19］中观察到的微裂纹在黏合线附近发生偏转或止裂的实验现象相吻合，从断裂力学角度证实了黏合线对微裂纹扩展具有抑制作用。

为分析标准骨单位半径R的变化对裂纹尖端A和B处应力强度因子的影响，设定 L=\mathrm{50}\mathrm{\mu }\mathrm{m}， d=\mathrm{20}\mathrm{\mu }\mathrm{m}。图4展示了微裂纹尖端标准化的SIF与标准化骨单位半径的变化关系。结果表明， 随着骨单位半径的增大，裂纹尖端的应力强度因子呈下降趋势。这表明，尺寸越大的骨单位对微裂纹扩展的抑制能力越强。

图5考察了标准化距离d/R对裂纹尖端A与B处应力强度因子的影响。计算中取R=\mathrm{150}\mathrm{\mu }\mathrm{m}， L=\mathrm{50}\mathrm{\mu }\mathrm{m}。结果显示， 当微裂纹靠近黏合线（即d较小）时， 尖端A的应力强度因子明显大于尖端B； 随着d增大， 两端应力强度因子逐渐趋近于相同数值。这表明，黏合线对裂纹扩展的抑制效应在靠近骨单位与间质骨交界处较强，并随着距离的增加而减弱，该结论与文献［14，20］的结果一致。

4 结论

文中建立了单轴拉伸载荷下骨单位含单条径向微裂纹的密质骨力学模型，并对推导出的奇异积分方程组进行了求解。结果表明，骨单位与间质骨的剪切模量比对裂纹尖端应力强度因子具有显著影响，刚度较高的骨单位有助于抑制微裂纹扩展，而刚度相对较低的骨单位则会促进裂纹定向扩展；当裂纹尖端靠近黏合线时，裂纹扩展表现出明显的几何屏蔽效应，通过偏转裂纹路径来阻碍其扩展。本研究从断裂力学角度揭示了密质骨的微观结构对微裂纹扩展的屏蔽机制，为理解密质骨在微观尺度下的损伤演化与断裂行为提供了理论依据，对骨生物力学研究、仿生材料设计以及骨质疏松等疾病的临床评估具有一定参考价值。未来可考虑引入更为复杂的微观密质骨结构，例如建立多骨单位平面模型，从而更真实地反映骨的复杂力学环境。

参考文献

[1] GAO Yuan， LI Bo， WANG Jianshan， et al. Fracture toughness analysis of helical fiber-reinforced biocomposites［J］. Journal of the Mechanics and Physics of Solids，2021，146：104206.DOI：10.1016/j.jmps.2020.104206 .

[2] MAGHAMI E， JOSEPHSON T O， MOORE J P， et al. Fracture behavior of human cortical bone： Role of advanced glycation end-products and microstructural features［J］. Journal of Biomechanics， 2021， 125： 110600. DOI：10.1016/j. jbiomech.2021.110600 .

[3] DAPAAH D， WILLETT T. A critical evaluation of cortical bone fracture toughness testing methods［J］. Journal of the Mechanical Behavior of Biomedical Materials，2022，134：105419.DOI：10.1016/j.jmbbm.2022.105419 .

[4] MAGHAMI E， NAJAFI A. Microstructural fatigue fracture behavior of glycated cortical bone［J］. Medical & Biological Engineering & Computing，2023，61（11）： 3021-3034.

[5] CARLSSON J， BRAESCH-ANDERSEN A， FERGUSON S J， et al. Fracture in porous bone analysed with a numerical phase-field dynamical model［J］. Journal of the Mechanical Behavior of Biomedical Materials，2023，139：105659.DOI：10.1016/j.jmbbm.2023.105659 .

[6] BLACK J D， TADROS B J. Bone structure：From cortical to calcium［J］. Orthopaedics and Trauma， 2020， 34（3）： 113-119.

[7] FROST H. Presence of microscopic cracks in vivo in bone［J］. Henry Ford Hospital Medical Bulletin， 1960， 8（1）：25-35.

[8] CURREY J D. Stress concentrations in bone［J］. Journal of Cell Science， 1962， S3-103（61）： 111-133.

[9] REILLY D T， BURSTEIN A H. The elastic and ultimate properties of compact bone tissue［J］. Journal of Biomechanics， 1975， 8（6）： 393-405.

[10] YIN Dagang， CHEN Bin， LIN Shiyun. Finite element analysis on multi-toughening mechanism of microstructure of osteon［J］. Journal of the Mechanical Behavior of Biomedical Materials， 2021， 117： 104408. DOI：10.1016/j.jmbbm. 2021.104408 .

[11] RAEISI NAJAFI A， ARSHI A R， ESLAMI M R， et al. Haversian cortical bone model with many radial microcracks： An elastic analytic solution［J］. Medical Engineering & Physics， 2007， 29（6）： 708-717.

[12] RAEISI NAJAFI A， ARSHI A R， SAFFAR K P， et al. A fiber-ceramic matrix composite material model for osteonal cortical bone fracture micromechanics： Solution of arbitrary microcracks interaction［J］. Journal of the Mechanical Behavior of Biomedical Materials， 2009， 2（3）： 217-223.

[13] WANG Xu， CHEN Yaogeng， DING Shenghu， et al. Cortical bone model with a microcrack under tensile loading［J］. Wuhan University Journal of Natural Sciences， 2022， 27（5）： 375-382.

[14] JI Chunhui， YANG Xiuyan， ZHANG Liang， et al. Microcrack behavior in bone： Stress field analysis at osteon cement line tips［J］. Proceedings of the Institution of Mechanical Engineers Part H， Journal of Engineering in Medicine， 2024， 238（8/9）： 909-921.

[15] ROBERTSON D M， ROBERTSON D， BARRETT C R. Fracture toughness， critical crack length and plastic zone size in bone［J］. Journal of Biomechanics， 1978， 11（8/9）： 359-364.

[16] MUSKHELISHVILI N I. Some basic problems of the mathematical theory of elasticity［M］.2nd ed. Groningen： P. Noordhoff， 1963.

[17] DUNDURS J， MURA T. Interaction between an edge dislocation and a circular inclusion［J］. Journal of the Mechanics and Physics of Solids， 1964， 12（3）： 177-189.

[18] XU Qiuyan. A kind of fast iterative methods with the application based on diagonal matrix splitting ［J］. Journal of Ningxia University （Natural Science Edition in Chinese and English）， 2026， 47（1）： 1-13.

[19] O’BRIEN F J， TAYLOR D， CLIVE LEE T. Bone as a composite material： The role of osteons as barriers to crack growth in compact bone［J］. International Journal of Fatigue， 2007， 29（6）： 1051-1056.

[20] GUO X E， LIANG L C， GOLDSTEIN S A. Micromechanics of osteonal cortical bone fracture［J］. Journal of Biomechanical Engineering， 1998， 120（1）： 112-117.

基金资助

宁夏自然科学基金资助项目(2023AAC03334)

国家自然科学基金资助项目(62566048)