一种基于有限元的超声理疗中软组织影响下的骨组织内部应力分布数值模拟方法

摘要

一种基于有限元的超声理疗中软组织影响下的骨组织内部应力分布数值模拟方法，本发明涉及骨组织内部应力分布数值模拟方法。本发明是要解决无法进行准确测量和得到超声引起骨组织的准确应力分布数据的问题而提出的一种基于有限元的超声理疗中软组织影响下的骨组织内部应力分布数值模拟方法，该方法是通过一、建立有限元模型；二、定义骨组织机械特性；三、进行耦合处理；四、设定理疗超声参数和作用方式；五、得到骨组织模型单元内部应力-应变分布；六、绘制应力-应变云图；七、量化测量路径应力分布结果等步骤实现的；本发明应用于生物力学研究领域。

说明书

技术领域

本发明涉及基于有限元的超声理疗中软组织影响下的骨组织内部应力分布数值模拟方法。 

背景技术

自从1880年P&JCurie发现压电效应以来，超声被广泛的应用于医学的各个领域。与CT、MRI等其他医学技术相比，超声技术操作简便、灵活、无放射性且具有较高的经济性。在骨骼超声检测方面，除了透射式超声和背散射式超声能够用于骨质健康状况检测以外，超声还能够作为一种有效的理疗工具来促进骨折愈合、骨组织新陈代谢。具体的实施方式是把骨骼组织置于具有一定能量的功率超声照射中。 

骨组织超声理疗主要是利用了超声声流引起的三种生物效应：热效应、空化效应和机械效应。热效应主要是由于粘滞性骨组织吸收了声流中的能量引起的。William指出超声能够引起组织0.86℃/min(1W/cm2,1-MHz)的温升。空化效应和机械效应是一种来自超声波传播的物理性应力效应，这种应力效应能够促进细胞再生、细胞新陈代谢和新骨形成。但是随着功率和照射时间的增强，理疗超声能够对骨组织产生破坏作用，这种破坏作用与超声引起的骨组织内部应力分布紧密相关，而基于当前的研究成果和实验条件，很难利用现有设备对生物组织中的超声能流分布进行在体测量，从而无法得到超声引起骨组织的准确的内部应力分布数据。一些研究方法中采用替代性生物材料来间接研究超声照射下的骨组织应力-应变及热能分布，但是这种材料在材料特性及形状方面跟真实骨组织差别很大，不能准确反映真实骨组织中的应力分布。 

发明内容

本发明的目的是为了解决利用现有实验条件和实验方法无法对理疗超声照射下的骨组织应力-应变分布进行准确测量和无法得到超声引起骨组织的准确应力分布数据的问题而提出的一种基于有限元的超声理疗中软组织影响下的骨组织内部应力分布数值模拟方法。 

上述的发明目的是通过以下技术方案实现的： 

步骤一、建立骨组织和骨组织外覆软组织有限元模型； 

步骤二、定义骨组织机械特性和骨组织外覆软组织机械特性； 

步骤三、对骨组织和骨组织外覆软组织之间界面进行耦合处理； 

步骤四、设定骨组织和骨组织外覆软组织应力分布数值模拟中的理疗超声参数和理疗 超声作用方式； 

步骤五、在骨组织和外覆软组织内部对线弹性波动方程进行有限元求解，根据骨组织和外覆软组织单元各节点所得应力-应变计算结果，利用Lagrangian插值函数得到骨组织和外覆软组织模型单元内部应力-应变分布； 

步骤六、输出利用Lagrangian插值函数得到的骨组织单元内部应力-应变分布值，绘制骨组织应力-应变云图； 

步骤七、在不同骨组织层中定义应力测量路径，量化测量路径应力分布结果；即完成了一种基于有限元的超声理疗中软组织影响下的骨组织内部应力分布数值模拟方法。 

发明效果： 

本发明提出的一种基于有限元的超声理疗中软组织影响下的骨组织内部应力分布数值模拟方法能够对骨组织内部超声诱导应力分布进行准确求解，为今后的骨骼超声理疗提供相关数据基础，从而为超声理疗医师提供指导，由于超声理疗与超声引起的骨组织及外敷软组织内部应力分布紧密相关，因此本发明在骨组织和骨组织外覆软组织内部对线弹性波动方程（骨组织介质中的应力传播方程）进行有限元求解，根据骨组织和外覆软组织单元各节点所得应力-应变计算结果，利用Lagrangian插值函数得到骨组织和外覆软组织模型得到的单元内部应力-应变分布值（冯米斯应力），绘制骨组织应力-应变云图；根据骨组织应力-应变云图量化不同骨组织层中应力分布结果，如图8、9所示。对超声理疗下的不同骨组织层的应力-应变分布进行准确的描述和量化如图11、12、13和14。这种应力-应变分布可进一步用来对超声能流分布及其所引起的生物效应进行预测。从而使得本发明为生物医学工程应用提供了一种有效地辅助研究方法，并且可以作为一种对骨组织内部声流分布的虚拟测量方法。 

附图说明

图1是具体实施方式一中提出一种基于有限元的超声理疗中软组织影响下的骨组织内部应力分布数值模拟方法流程图； 

图2具体实施方式二数值模拟所采用骨组织及骨组织外覆软组织尺寸示意图，其中1为骨组织和骨组织外覆的软组织的耦合界面、2为骨组织外覆的软组织、3为骨组织、Rb为骨组织的中空圆柱结构的内半径、Rmb为骨组织的中空圆柱结构的外半径、Rm为骨组织外覆软组织的中空圆柱结构外半径； 

图3是具体实施方式二提出的实际模型建立为中空圆柱结构的一半； 

图4是具体实施方式二提出的骨组织和骨组织外覆软组织的中空圆柱结构对称面上 施加对称约束示意图； 

图5是具体实施方式二提出的骨组织及骨组织外覆软组织离散化模型示意图； 

图6是具体实施方式五提出的超声应力波在骨组织及骨组织外覆软组织中的传播过程示意图； 

图7是具体实施方式四提出的骨组织及其外覆软组织界面节点耦合示意图，其中Ux为骨组织模型轴向方向的自由度值，Uy为骨组织模型外表面法线方向的自由度值，Uz为骨组织模型外表面切线方向自由度值；U’x为骨组织外覆软组织模型轴向方向的自由度值，U’y为骨组织外覆软组织模型内表面法线方向的自由度值，U’z为骨组织外覆软组织模型内表面切线方向自由度值； 

图8是具体实施方式七提出的骨组织及其外覆软组织应变场分布示意图，其中图8正下方渐变图例代表骨组织及其外覆软组织内部位应变分布值范围为-0.796E-16～-0.720E-10m； 

图9是具体实施方式七提出的骨组织及其外覆软组织应力场分布示意图，其中图9正下方渐变图例代表骨组织及其外覆软组织内部位应力分布值范围为0.240E-5～3646.01Pa； 

图10是具体实施方式八中提出的定义的不同骨组织层中的应力计算路径示意图，1为超声波发生器，2为骨组织外层，3为骨组织中间层； 

图11是具体实施方式八中提出的以轴向距离为横坐标，应力幅值为纵坐标的计算路径A的骨组织内顶端外层轴向应力分布量化结果曲线图； 

图12是具体实施方式八中提出的以轴向距离为横坐标，应力幅值为纵坐标的计算路径B的骨组织内顶端外层轴向应力分布量化结果曲线图； 

图13是具体实施方式八中提出的提出的以轴向距离为横坐标，应力幅值为纵坐标的计算路径C的骨组织内顶端外层轴向应力分布量化结果曲线图； 

图14是具体实施方式八中提出的以轴向距离为横坐标，应力幅值为纵坐标的计算路径D的骨组织内顶端外层轴向应力分布量化结果曲线图。 

具体实施方式

具体实施方式一：本实施方式的一种基于有限元的超声理疗中软组织影响下的骨组织内部应力分布数值模拟方法具体是按照以下步骤制备的： 

步骤一、建立骨组织和骨组织外覆软组织有限元模型； 

步骤二、定义骨组织机械特性和骨组织外覆软组织机械特性； 

步骤三、对骨组织和骨组织外覆软组织之间界面进行耦合处理； 

步骤四、设定骨组织和骨组织外覆软组织应力分布数值模拟中的理疗超声参数和理疗超声作用方式； 

步骤五、在骨组织和外覆软组织内部对线弹性波动方程（骨组织介质中的应力传播方程）进行有限元求解，根据骨组织和外覆软组织单元各节点所得应力-应变计算结果，利用Lagrangian插值函数得到骨组织和外覆软组织模型单元内部应力-应变分布； 

步骤六、输出利用Lagrangian插值函数得到的骨组织单元内部应力-应变分布值，绘制骨组织应力-应变云图； 

步骤七、在不同骨组织层中定义应力测量路径，量化测量路径应力分布结果；即完成了一种基于有限元的超声理疗中软组织影响下的骨组织内部应力分布数值模拟方法；如图1。 

本实施方式效果： 

本实施方式提出的一种基于有限元的超声理疗中软组织影响下的骨组织内部应力分布数值模拟方法能够对骨组织内部超声诱导应力分布进行准确求解，为今后的骨骼超声理疗提供相关数据基础，从而为超声理疗医师提供指导，由于超声理疗与超声引起的骨组织及外敷软组织内部应力分布紧密相关，因此本实施方式在骨组织和骨组织外覆软组织内部对线弹性波动方程（骨组织介质中的应力传播方程）进行有限元求解，根据骨组织和外覆软组织单元各节点所得应力-应变计算结果，利用Lagrangian插值函数得到骨组织和外覆软组织模型得到的单元内部应力-应变分布值（冯米斯应力），绘制骨组织应力-应变云图；根据骨组织应力-应变云图量化不同骨组织层中应力分布结果，如图8、9所示。对超声理疗下的不同骨组织层的应力-应变分布进行准确的描述和量化如图11、12、13和14。这种应力-应变分布可进一步用来对超声能流分布及其所引起的生物效应进行预测。从而使得本实施方式为生物医学工程应用提供了一种有效地辅助研究方法，并且可以作为一种对骨组织内部声流分布的虚拟测量方法。 

具体实施方式二：本实施方式与具体实施方式一不同的是：步骤一中建立骨组织和骨组织外覆软组织有限元模型的具体过程如图2为： 

(1)建立骨组织和骨组织外覆软组织有限元模型，其中骨组织和骨组织外覆软组织模型分别建立为中空圆柱结构（骨组织及其外覆软组织均等效为中空圆柱结构）；由于声传播介质的形状及软组织的存在对超声能流分布具有很大的影响，二者必须予以考虑；在对骨组织进行建模时需要尽可能接近真实骨骼生理结构。 

(2)将骨组织的中空圆柱结构的内半径和外半径分别设为4mm和8mm，将骨组织外覆软组织的中空圆柱结构内半径和外半径分别设为8mm和16mm；所用于进行数值模拟分析的骨组织样本长度取为50mm； 

(3)由于骨组织及外覆软组织为对称结构，为了减少计算量，在进行有限元计算时实际模型建立为中空圆柱结构的一半，如图3并在骨组织和骨组织外覆软组织的中空圆柱结构对称面上施加对称约束如图4； 

(4)数值模拟前，需要对建立的有限元模型进行离散化如图5，对建立的骨组织和骨组织外覆软组织的中空圆柱模型进行离散化，离散化过程采用显示欧拉向前积分法，将离散化单元采用3-D8节点结构实体单元，每个单元节点具有三个方向上的位移自由度Ux，Uy，Uz，为了保证计算精度并尽可能地减少计算时间，设定模型离散化单元尺寸时需满足： 

$\frac{1}{10}{\lambda }_{\mathrm{be}}\le L\le \frac{1}{20}{\lambda }_{\mathrm{be}}---\left(1\right)$

式中，L为离散化单元尺寸，λ_be为骨组织中超声波长。其它步骤及参数与具体实施方式一相同。 

具体实施方式三：本实施方式与具体实施方式一或二不同的是：步骤二中定义骨组织机械特性和骨组织外覆软组织机械特性为： 

进行模拟计算时需要对骨组织及其外覆软组织材料特性进行准确定义，骨组织主要由胶原蛋白和无机矿物质基质组成，其中胶原蛋白是一种聚合物，具有J型应力-应变曲线关系，骨组织应力在超过应变阈值时会随应变迅速增加；软组织是一种高弹体，宏观上应力-应变关系呈现非线性，但在应变<3%-5%时，可等效为线性关系。将骨组织和骨组织外覆软组织置于理疗超声辐射中时，骨组织和骨组织外覆软组织应变范围<1%（骨组织和骨组织外覆软组织质点位移范围很小<1%），本发明实施例中骨组织及其外覆软组织在小应变有限元计算时近似为各向同性的线性弹性体，（在微小的应变尺度下，骨组织机械特性和骨组织外覆软组织机械特性具有线性关系），骨组织机械特性和骨组织外覆软组织应力-应变关系近似满足胡克定律，应力计算误差可以忽略。其它步骤及参数与具体实施方式一或二相同。 

具体实施方式四：本实施方式与具体实施方式一至三之一不同的是：对骨组织和骨组织外覆软组织之间界面进行耦合处理主要过程如图7为： 

(1)设定Ux为骨组织模型轴向方向的自由度值，Uy为骨组织模型外表面法线方向的 自由度值，Uz为骨组织模型外表面切线方向自由度值； 

(2)设定U’x为骨组织外覆软组织模型轴向方向的自由度值，U’y为骨组织外覆软组织模型内表面法线方向的自由度值，U’z为骨组织外覆软组织模型内表面切线方向自由度值； 

(3)将骨组织和骨组织外覆软组织之间进行界面耦合，使骨组织节点Bn与骨组织外覆软组织节点Sn在模型表面法线方向上具有相同的自由度值Uy=U’y，对骨组织与骨组织外覆软组织的接触界面进行单元离散化时，骨组织和骨组织外覆软组织的界面两侧的节点数量和节点位置相对应，也就是骨组织节点Bn与骨组织外覆软组织节点Sn相对应，骨组织节点Bn自由度Ux、Uz与软组织节点Sn自由度U’x、U’z不相关。其它步骤及参数与具体实施方式一至三之一相同。 

具体实施方式五：本实施方式与具体实施方式一至四之一不同的是：步骤四中设定骨组织和骨组织外覆软组织应力分布数值模拟中的理疗超声参数和理疗超声作用方式为：数值模拟中理疗超声采用2-MHz谐响应作用模式，激励峰-峰值为2-μm，其中理疗超声激励点作用于软组织表面上靠近端面的5*5个节点上如图6。其它步骤及参数与具体实施方式一至四之一相同。 

具体实施方式六：本实施方式与具体实施方式一至五之一不同的是：步骤五中在骨组织和外覆软组织内部对线弹性波动方程（骨组织介质中的应力传播方程）进行有限元求解，根据骨组织和外覆软组织单元各节点所得应力-应变计算结果，利用Lagrangian插值函数得到骨组织和外覆软组织模型单元内部应力-应变分布的过程为： 

(1)数值模拟是基于对无外力作用下的骨组织和外覆软组织内部线弹性波动方程（骨组织介质中的应力传播方程）进行求解，线弹性波动方程表示为： 

$(\lambda +2\mu )▿▿·u+\mu {▿}^{2}u=\rho \frac{{\partial }^{2}u}{{\partial }^{2}t}---\left(2\right)$

式中，λ和μ为Lame常量，u为骨组织或外覆软组织内部应力场向量，ρ为骨组织密度； 

为了利用声速对位移场进行求解引入骨组织或外覆软组织中超声纵波速度VL和剪切波速度VS分别为： 

$\left(\begin{array}{cc}{V}_L{\left(\frac{\lambda +2\mu }{\rho }\right)}^{\frac{1}{2}}& V_S={\left(\frac{\mu }{\rho }\right)}^{\frac{1}{2}}\end{array}\right)---\left(3\right)$

公式（3）带入方程（4），可将对线弹性方程的求解改为求解以下方程： 

${V_L}^2▿▿·u+{V_S}^2{▿}^{2}u\frac{{\partial }^{2}u}{{\partial }^{2}t}=0---\left(4\right)$

u为骨组织或外覆软组织内部应力场向量，t为时间； 

进一步引入能量泛函求解应力-应变分布: 

$\Pi (u,t)=\frac{1}{2}\int {\left[ϵ\right]}^T\left[\sigma \right]\mathrm{dv}-{\int }_v{u}^t\mathrm{fdv}---\left(5\right)$

式中，σ为应力，ε为应变，f为力密度， 

(2)根据线弹性方程和能量泛函，对应力σ和应变ε进行求解，求出第一主应力、第二主应力和第三主应力σ1，σ2，σ3；根据骨组织和外覆软组织各节点应力-应变计算求解结果，利用Lagrangian插值函数对单元内部应力-应变分布进行计算获得骨组织和外覆的软组织单元内部应力分布值。其它步骤及参数与具体实施方式一至五之一相同。 

具体实施方式七：本实施方式与具体实施方式一至六之一不同的是步骤六输出利用Lagrangian插值函数得到的骨组织单元内部应力（如图9）-应变（如图8）分布值，绘制骨组织应力-应变云图过程为： 

根据线弹性波动方程（骨组织介质中的应力传播方程）求得骨组织和外覆软组织的第一主应力、第二主应力和第三主应力σ1，σ2，σ3，利用第四强度理论求得骨组织和外覆软组织的冯米斯应力（当量应力）σ_e为： 

${\sigma }_e={(\frac{1}{2}·[{(\sigma 1-\sigma 2)}^2+{(\sigma 2-\sigma 3)}^2+{(\sigma 3-\sigma 1)}^2\left]\right)}^{\frac{1}{2}}$

根据骨组织单元内部应力σ_e-应变分布值绘制骨组织应力-应变云图。其它步骤及参数与具体实施方式一至六之一相同。 

具体实施方式八：本实施方式与具体实施方式一至七之一不同的是：步骤七在不同骨组织层中定义应力测量路径，量化测量路径应力分布结果具体过程如图10： 

(1)定义不同骨组织层中的应力计算路径：在骨组织端面上沿波束传播方向上选取A、C两点；在骨组织端面上垂直波速方向选取B、D两点；其中C、D为骨组织中间层的点，A、B为骨组织外层上的点； 

(2)过骨组织外层和骨组织中间层选取的A、B、C和D四个点沿骨轴向定义四条路径；根据步骤六冯米斯应力计算结果σ_e及定义的四条路径，计算出在不同骨层定义的四条路径上的应力量化值；如图11、12、13、14。其它步骤及参数与具体实施方式一至七之一相同。 

实施例一 

本实施例一种基于有限元的超声理疗中软组织影响下的骨组织内部应力分布数值模拟方法，具体是按照以下步骤制备的： 

步骤一、建立骨组织和骨组织外覆软组织有限元模型； 

(1)建立骨组织和骨组织外覆软组织有限元模型，其中骨组织和骨组织外覆软组织模型分别建立为中空圆柱结构模型的建立采用CYL4命令实现； 

(2)将骨组织的中空圆柱结构的内半径和外半径分别设为4mm和8mm，将骨组织外覆软组织的中空圆柱结构内半径和外半径分别设为8mm和16mm；所用于进行数值模拟分析的骨组织样本长度取为50mm； 

(3)在进行有限元计算时实际模型建立为中空圆柱结构的一半，如图3并在骨组织和骨组织外覆软组织的中空圆柱结构对称面上施加对称约束如图4其建立方式分别为： 

CYL4,0,0,4E-3,0,8E-3,-180,5E-2 

CYL4,0,0,8E-3,0,16E-3,-180,5E-2； 

并在骨组织和骨组织外覆软组织的中空圆柱结构对称面上施加对称约束方式为： 

CSYS,0 

NSEL,S,LOC,Y,0 

DSYM,SYMM,Y,O 

CSYS,1；； 

(4)对建立的骨组织和骨组织外覆软组织的中空圆柱模型进行离散化，离散化过程采用显示欧拉向前积分法，将离散化单元采用3-D8节点结构实体单元SOLID185，每个单元节点具有三个方向上的位移自由度Ux,Uy,Uz，设定模型离散化单元尺寸时需满足如图5： 

$\frac{1}{10}{\lambda }_{\mathrm{be}}\le L\le \frac{1}{20}{\lambda }_{\mathrm{be}}---\left(1\right)$

式中，L为离散化单元尺寸，λ_be为骨组织中超声波长，离散化后共获得节点数为1,304,704，单元数为1,200,207； 

模型离散化具体实现方式为： 

周向离散化： 

LSEL,S,LOC,Z,0 

LSEL,R,LOC,X,8E-3 

LESIZE,ALL,,,CNUM,,1,,,0 

径向离散化： 

LSEL,S,LOC,X,(4E-3+8E-3)/2 

LSEL,A,LOC,X,(8E-3+16E-3)/2 

LESIZE,ALL,,,TNUM,,1,,,0 

轴向离散化： 

LSEL,S,LOC,X,8E-3 

LSEL,R,LOC,Z,5E-2/2 

LESIZE,ALL,LLEN,,,,1,,,0 

步骤二、定义骨组织机械特性和骨组织外覆软组织机械特性为： 

由于骨组织主要由胶原蛋白和无机矿物质基质组成，其中胶原蛋白是一种聚合物，聚合物的应力-应变曲线为J型曲线，骨组织应力在超过应变阈值时会随应变迅速增加；骨组织外覆软组织是一种高弹体，宏观上应力-应变关系呈现非线性，但在骨组织和骨组织外覆软组织应变<3%-5%时，可等效为线性关系；将骨组织和骨组织外覆软组织置于理疗超声辐射中时，骨组织和骨组织外覆软组织应变范围<1%（骨组织和骨组织外覆软组织质点位移范围很小<1%），可将骨组织机械特性和骨组织外覆软组织机械特性在小应变有限元计算时近似为各向同性的线性弹性体，其应力-应变关系近似满足胡克定律，应力计算误差可以忽略，其中骨组织及外覆软组织材料属性如下表所示： 

表1 

骨组织机械特性定义方式为： 

MP,DENS,1,1200 

MP,EX,1,10E9 

MP,NUXY,1,0.4 

外覆软组织机械特性定义为： 

MP,DENS,2,1600 

MP,EX,2,8.2E6 

MP,NUXY,2,0.47 

步骤三、对骨组织和骨组织外覆软组织之间界面进行耦合处理主要过程为： 

(1)设定Ux为骨组织模型轴向方向的自由度值，Uy为骨组织模型外表面法线方向的自由度值，Uz为骨组织模型外表面切线方向自由度值； 

(2)设定U’x为骨组织外覆软组织模型轴向方向的自由度值，U’y为骨组织外覆软组织模型内表面法线方向的自由度值，U’z为骨组织外覆软组织模型内表面切线方向自由度值； 

(3)将骨组织和骨组织外覆软组织之间进行界面耦合，使骨组织节点Bn与骨组织外覆软组织节点Sn在模型表面法线方向上具有相同的自由度值Uy=U’y，对骨组织与骨组织外覆软组织的接触界面进行单元离散化时，骨组织和骨组织外覆软组织的界面两侧的节点数量和节点位置相对应，也就是骨组织节点Bn与骨组织外覆软组织节点Sn相对应，在骨组织端切线方向Uz和骨组织端面轴向方向Ux上则骨组织节点Bn与软组织节点Sn自由度不相关； 

步骤四、设定骨组织和骨组织外覆软组织应力分布数值模拟中的理疗超声参数和理疗超声作用方式为：数值模拟中理疗超声采用2-MHz谐响应作用模式，激励峰-峰值为2-μm，其中理疗超声激励点作用于软组织表面上靠近端面的5*5个节点上如图6； 

其中激励峰峰值定义如下： 

ALLSEL 

CMSEL,S,SENSORNODE,NODE!添加振动点 

NROTAT,ALL 

D,ALL,UX,2e-6 

谐响应模式定义方式如下： 

ANTYPE,HARMIC 

HARF,2e6,2e6 

步骤五、在骨组织和外覆软组织内部对线弹性波动方程（骨组织介质中的应力传播方程）进行有限元求解，根据骨组织和外覆软组织单元各节点所得应力-应变计算结果，利用Lagrangian插值函数得到骨组织和外覆软组织模型单元内部应力-应变分布为： 

(1)数值模拟是基于对无外力作用下的线性骨组织内部对应力传播方程（骨组织介质中的线弹性波动方程）进行有限元求解，线弹性波动方程表示为： 

$(\lambda +2\mu )▿▿·u+\mu {▿}^{2}u=\rho \frac{{\partial }^{2}u}{{\partial }^{2}t}---\left(2\right)$

式中，λ和μ为Lame常量，u为骨组织内部应力场向量，ρ为骨组织密度； 

引入骨组织中超声纵波速度V_L和剪切波速度V_S分别为： 

$\left(\begin{array}{cc}{V}_L{\left(\frac{\lambda +2\mu }{\rho }\right)}^{\frac{1}{2}}& V_S={\left(\frac{\mu }{\rho }\right)}^{\frac{1}{2}}\end{array}\right)---\left(3\right)$

公式带入方程，可将对线弹性方程的求解改为求解以下方程： 

${V_L}^2▿▿·u+{V_S}^2{▿}^{2}u\frac{{\partial }^{2}u}{{\partial }^{2}t}=0---\left(4\right)$

u为骨组织内部应力场向量，t为时间； 

进一步引入能量泛函求解应力-应变分布: 

$\Pi (u,t)=\frac{1}{2}\int {\left[ϵ\right]}^T\left[\sigma \right]\mathrm{dv}-{\int }_v{u}^t\mathrm{fdv}---\left(5\right)$

式中，σ为应力，ε为应变，f为力密度， 

(2)根据线弹性方程和能量泛函，对应力σ和应变ε进行求解，求出第一主应力、第二主应力和第三主应力σ1，σ2，σ3；根据骨组织与外覆的软组织各节点应力-应变计算求解结果，利用Lagrangian插值函数对单元内部应力-应变分布进行计算获得单元内部应力分布值； 

步骤六、输出利用Lagrangian插值函数计算得到的骨组织单元内部应力-应变分布值，绘制骨组织应力（如图9）-应变（如图8）云图过程为： 

根据骨组织和软组织中应力传播方程求得的第一主应力、第二主应力和第三主应力σ1，σ2，σ3，利用第四强度理论求得骨组织和软组织的冯米斯应力（当量应力）σ_e为： 

${\sigma }_e={(\frac{1}{2}·[{(\sigma 1-\sigma 2)}^2+{(\sigma 2-\sigma 3)}^2+{(\sigma 3-\sigma 1)}^2\left]\right)}^{\frac{1}{2}}$

采用等效冯米斯应力(vonMisesStress)绘制应力云图如图9； 

步骤七、在不同骨组织层中定义应力测量路径，从骨组织应力-应变云图中提取测量路径，量化测量路径应力分布结果过程为： 

(1)定义不同骨组织层中的应力计算路径：在骨组织端面上沿波数传播方向上竖直选取A、C两点；在骨组织端面上垂直波数传播方向并过圆心选取B、D两点；其中C、D为骨组织中间层的点，A、B为骨组织外层上的点； 

(2)将骨组织外层和骨组织中间层选取的A、B、C和D四个点的不同位置沿骨轴向定 义四条路径，其中A、C为顶端位置，B、D为侧面位置； 

(3)分别对顶端为A、B、C和D四条路径上的应力分布进行计算：根据步骤六冯米斯应力计算结果σ_e及定义的四条路径，计算出在不同骨层定义的四条路径上的应力量化值；根据四条路径上的应力量化值来计算骨组织内部四条路径的位置，以及沿每条径线上的不同位置的应力分布结果，并将该路径的应力分布结果提取出来如图11、图12、图13和图14所示。