CN109492250A - 一种非线性各向异性质量弹簧模型的软组织变形模拟方法 - Google Patents

Abstract

本发明的基于卷非线性各向异性质量弹簧模型的软组织变形模拟方法，包括如下步骤；步骤1）构建质点弹簧模型，通过有限元法确定质点弹簧模型的初始弹性参数；步骤2）根据肌肉组织的应力、应变关系设置变形过程中弹簧参数；步骤3）通过构建弹簧的力学方程将软组织的各向异性特征进行参数化的描述；步骤4）构建质点弹簧模型的动力学方程并求解。有益效果：本发明着力解决质点弹簧模型中弹簧参数的设置，并将各向异性参数化，同时在解决质点弹簧模型的动力学方程中采用了改进的欧拉算法。

Description

一种非线性各向异性质量弹簧模型的软组织变形模拟方法

技术领域

本发明属于软组织形变模型建模技术领域，特别涉及了一种非线性各向异性质量弹簧模型的软组织变形模拟方法。

背景技术

随着科学技术的不断发展，虚拟手术模拟系统可以极大地降低医学研究的难度和医疗实践的成本，因此成为近年来的热门话题，虚拟手术模拟系统也得到了前所未有的发展。模拟软组织变形是虚拟手术模拟系统中尚未解决的问题，一方面在于医学领域中真实的软组织在形变的过程中会表现出复杂的生物力学特性，另一方面则是软组织模型在虚拟手术训练中需要很高的计算效率和物理精度，这给软组织建模带来了巨大困难。

因此，近年来关于软组织形变的模拟成为医学虚拟训练中研究的焦点。常见的虚拟手术软组织形变的模型有质点弹簧模型和有限元模型。其中质点弹簧模型是最早提出的物理模型，也是较为经典的一种形变模型，由于其计算简单，建模速度较快而被广泛的应用于虚拟手术软组织仿真研究中。而有限单元法是求解复杂微分方程近似解的有效工具，是工程设计和分析的可靠工具。因其建模精度较高，稳定性较好，模型参数易于调节的优势，逐渐被应用到软组织建模中。综合比较两种方法的优缺点，开发一个计算成本低，模拟效果好的模型至关重要。

对于大多数现有的质点弹簧模型，每个弹簧的设置参数都是一样的。而软组织通常具有粘弹性、非线性、各向异性等复杂特征。因此，为所有弹簧设置相同的参数缺乏理论基础。另外，软组织变形过程中弹簧的不变特性在极端情况下是有问题的。

发明内容

本发明的目的是克服上述背景技术的不足，着力解决质点弹簧模型中弹簧参数的设置，提出了一种利用有限元法确定质点弹簧参数的改进方法，并将各向异性参数化，同时在解决质点弹簧模型的动力学方程中采用了改进的欧拉算法。为了实现上述技术目的，本发明的技术方案为：

所述非线性各向异性质量弹簧模型的软组织变形模拟方法，包括如下步骤；

步骤1)构建质点弹簧模型，通过有限元法确定质点弹簧模型的初始弹性参数；

步骤2)根据肌肉组织的应力、应变关系设置变形过程中弹簧参数；

步骤3)通过构建弹簧的力学方程将软组织的各向异性特征进行参数化的描述；

步骤4)构建质点弹簧模型的动力学方程并求解。

所述非线性各向异性质量弹簧模型的软组织变形模拟方法的进一步设计在于，所述步骤1)通过对比弹簧参数与软组织力学参数之间的关系得到初始弹性参数；所述弹簧参数包括阻尼系数k_d和弹性系数k_s，软组织力学参数包括杨氏模量和泊松比。

所述非线性各向异性质量弹簧模型的软组织变形模拟方法的进一步设计在于，所述步骤1)具体包括如下步骤:

步骤1-1)根据式(1)构建基于Kelvin模型的动力学方程，

其中m表示质量，X是位移，F_e是施加在每个质点上的力，k_d和k_s分别表示阻尼系数和弹性系数；

步骤1-2)设定Kelvin模型与有限元的链接单元都含有设于模型两端的两个节点，分别标记为第一节点与第二节点，设定F₁和F₂分别对应为施加在弹簧的第一节点和第二节点上的外力，令L和ΔL表示弹簧的原始长度和长度变化；

根据质点弹簧模型得到F₁的表达式如式(2)：

结合式(2)对于限元链接单元得到F₁的表达式如式(3)：

式(3)中，A是有限元链接单元的横截面积，σ是应力，通过式(5)计算得出：

σ＝Eε (5)

式(5)中，E是杨氏模量，ε是应变，且ε＝ΔL/L；

根据式(1)与式(3)可得弹簧的初始弹性参数如式(5)

所述非线性各向异性质量弹簧模型的软组织变形模拟方法的进一步设计在于，所述步骤2)具体包括如下步骤：

步骤2-1)设定应变和应力之间的函数关系如式(6)：

σ＝ae^bε (6)

其中，a和b是系数，通过最小二乘拟合算法进行估算；

步骤2-2)根据杨氏模量和式(6)，通过式(7)计算E：

进而根据式(8)得到变形过程中的弹簧参数：

所述非线性各向异性质量弹簧模型的软组织变形模拟方法的进一步设计在于，所述步骤3)包括如下步骤：

步骤3-1)根据式(8)计算弹簧在轴向和径向的弹性系数和

步骤3-2)根据式(9)计算施加在弹簧上的力的轴向分量F^α和径向分量F^β：

其中，ΔL^α和ΔL^β分别表示ΔL的轴向分量和径向分量；

步骤3-3)根据式(10)计算弹簧上的作用力F，结合式(9)形成软组织的各向异性特征的表达：

F＝F^αcosθ+F^βsinθ (10)

其中，θ弹表示簧的轴向与软组织的轴向之间的夹角。

所述非线性各向异性质量弹簧模型的软组织变形模拟方法的进一步设计在于，所述步骤4)中构建质点弹簧模型的动力学方程的过程中设定质点弹簧模型基于四面体拓扑结构，并且由n个质点组成的，在任意时刻，模型的状态都由每个质点的位置x_i、速度v_i确定，任意质点i满足式(11)，

其中，m_i是质点i的质量，x_i表示质点i的位置矢量，f_in表示与质点i相连的所有弹簧施加于质点i的内力，由弹性力f_s和阻尼力f_d组成，f_s、f_d分别根据式(12)、式(13)所示，

f_d＝k_d(v_j-v_i) (13)

结合式(12)、式(13)得到式(14)形成质点弹簧模型的动力学方程：

其中，|x_j-x_i|表示表示两个质点间位置矢量之差，|x_j-x_i|⁰是发生形变前弹簧的原长，k_s为弹簧的弹性系数，k_d为弹簧的阻尼系数。

所述非线性各向异性质量弹簧模型的软组织变形模拟方法的进一步设计在于，所述步骤4)中对所述质点弹簧模型的动力学方程基于改进的欧拉法求解对所述动力学方程的速度与位移矢量进行求解，如式(15)，

其中，t为初始时间点，Δt为预先定义的迭代时间步长，k表示不同的空间点；所述改进的欧拉法包括如下步骤：

步骤4-1)采用链表形式存放软组织数据，构建质点弹簧体模型，同时设置模型的初始状态为静态，并求出各个质点间弹簧的初始长度L₀；

步骤4-2)用外接力反馈设备为软组织施加作用力，并返回发生碰撞的质点的序列号i；

步骤4-3)结合式(14)计算质点i的合力，并通过显式欧拉法v_i求解，如式(16)：

步骤4-4)根据步骤4-3)的计算结果，根据式(17)用隐式欧拉法求解x_i的值：

x_i(t+Δt)＝x_i(t)+Δt·v_i(t+Δt) (17)

步骤4-5)计算质点i邻域内其他质点的受力和运动状态；

步骤4-6)更新并绘制当前模型；

步骤4-7)重新计算质点i的合力，跳至步骤4-3)。

本发明的有益效果为：

本发明的非线性各向异性质量弹簧模型的软组织变形模拟方法，一方面通过对质点弹簧模型中的弹簧初始参数的设置,软组织形变过程中弹簧参数的变化以及各向异性参数化的分析，提高了质点弹簧模型的精度；另一方面，用显示欧垃积分对v_i进行推进，而采用隐式欧拉积分对x_i进行推进。本发明的技术方案与传统显示欧垃积分相比，既没有增加求解计算量，又因为在计算x_i(t+Δt)时使用了未来时刻物理量v_i(t+Δt)，所以它具有较显示欧拉积分法更好的稳定性。

附图说明

图1为Kelvin模型的示意图。

图2为Kelvin模型与有限元的链接单元的示意图。

图3肌肉组织的应力应变曲线示意图。

图4为基于四面体网络拓扑结构的质点弹簧模型示意图。

具体实施方式

下面对本申请作进一步详细的说明。

本实施例的非线性各向异性质量弹簧模型的软组织变形模拟方法包括如下步骤：

步骤1)构建质点弹簧模型，通过有限元法确定质点弹簧模型的初始弹性参数；

质点弹簧模型有许多不同的类型，最常用的一种是Kelvin模型，这种模型全部由弹性元件和阻尼元件组成，如图1所示。Kelvin模型的演化符合动力学基本定律的控制，即符合拉格朗日运动方程，如公式(1)所示：

其中m表示质量，X是位移，F_e是施加在每个质点上的力，k_d和k_s分别表示阻尼系数和弹性系数。由于阻尼系数通常远小于弹性系数，所以方程(1)可以简化如下：

和质点弹簧模型相比，有限元模型通常在软组织变形模拟中具有更高的精度。因此，本实施例在相同受力情况下，利用一个与质点弹簧模型相似的有限元模型，即可以通过有限元法确定质点弹簧模型参数，从而提高质点弹簧模型的精度。

为了达到与有限元模型相同的变形量，需要找出质点弹簧模型与有限元模型的关系，即弹簧参数(包括阻尼系数k_d和弹性系数k_s)与软组织力学参数(包括杨氏模量和泊松比)之间的关系，以及有限元单元的几何特征(包括节点的坐标、面积、厚度等)。如图2所示，给出Kelvin模型与有限元的链接单元。经比较，Kelvin模型与有限元的链接单元都有两个节点，分别标记为1,2。F₁和F₂是施加在弹簧的节点1和节点2上的外力，包括重力以及与节点相连的其它弹簧的弹簧力。令L和ΔL表示弹簧的原始长度和长度变化。一方面，对于质点弹簧来说，只考虑弹力的情况下，可以计算F₁：

另一方面，对于有限元链接单元，节点1的动力学方程可写为：

其中，A是有限元链接单元的横截面积，σ是应力，可以通过计算得出：

σ＝Eε (5)

其中，E是杨氏模量，ε是应变，且ε＝ΔL/L。

比较式(2)和式(4)，可以得到式(6)：

k_sX＝Aσ (6)

然后，弹簧的初始弹性参数可以表示为：

步骤2)根据肌肉组织的应力、应变关系设置变形过程中弹簧参数；

在大多数质点弹簧模型中，弹簧的弹性特征在软组织变形过程中保持不变。但是，根据材料的真实力学性能可知，弹性特征在软组织变形过程中保持不变是存在缺陷的。如图3，假设应变和应力之间的函数关系如公式(8)所示：

σ＝ae^bε (8)

其中，a和b是系数。可以使用最小二乘拟合算法来估计它们。图3中的曲线是该函数的表达式。

根据杨氏模量和方程(8)的定义，E可以通过下式计算，得到公式(9)：

然后，公式(7)可以进一步写成：

步骤3)通过构建弹簧的力学方程将软组织的各向异性特征进行参数化的描述。

各向异性是软组织的另一个特性，意味着杨氏模量E在不同方向上是不同的。本实施例只考虑正交方向，即轴向和径向。所以，k_s在这两个方向上也是不同的。假设和表示弹簧在轴向和径向的弹性系数。它们可以根据方程(10)计算。然后，施加在弹簧上的力的轴向分量和径向分量(由F^α和F^β表示)可写为：

其中，ΔL^α和ΔL^β分别表示ΔL的轴向分量和径向分量。假设弹簧的轴向与软组织的轴向之间的夹角为θ，则弹簧上的作用力(用F表示)可以通过下式计算：

F＝F^αcosθ+F^βsinθ (12)

通过方程(11)的分解和方程(12)的组成，软组织的各向异性特征被明确地表示出来。

步骤4)构建质点弹簧模型的动力学方程并求解。

本实施例的质点弹簧体模型基于四面体拓扑结构，参见图4。设整个模型是由n个质点组成的，在任意时刻，模型的状态都由每个质点的位置x_i、速度v_i确定。任意质点i满足公式(13)：

其中，m_i是质点i的质量，x_i表示质点i的位置矢量，f_in表示与质点i相连的所有弹簧施加于质点i的内力，由弹性力f_s和阻尼力f_d组成，如公式(14)和公式(15)所示，f_ext表示质点i所受外力，包括摩擦力、重力、用户施加的力等。

f_d＝k_d(v_j-v_i) (15)

将公式(14)和公式(15)代入公式(13)整理得公式(16)：

其中，|x_j-x_i|表示表示两个质点间位置矢量之差，|x_j-x_i|⁰是发生形变前弹簧的原长，k_s为弹簧的弹性系数，k_d为弹簧的阻尼系数。

目前，求解二阶微分方程的主要方法有欧拉法、龙格—库塔法。其中欧拉法分为显示欧拉法和隐式欧拉法。显式欧拉法计算量小，但是一般只一阶收敛，精度不高，为实现逼真的模拟效果，时间步长需要设置很小，从而导致整个形变过程延长，并且由于存在进退性冲击波，所以不适合用于质点弹簧模型。隐式欧拉法精度高，稳定性好，但是计算量却较大。龙格—库塔法也具有较高的精度和稳定性，但是，在一步计算中需要计算四次函数的值，因此计算量大，实时性差。因此提出一种改进的欧拉法来求解速度与位移矢量，如公式(17)所示：

其中，t是初始时间点，Δt是预先定义的迭代时间步长，k代表不同的空间点。

改进的欧拉法的具体步骤如下：

1)采用链表形式存放软组织数据，构建质点弹簧体模型，同时设置模型的初始状态为静态，并求出各个质点间弹簧的初始长度L₀；

2)用外接力反馈设备为软组织施加作用力，并返回发生碰撞的质点的序列号i；

3)利用公式(16)计算质点i的合力，并用显式欧拉法v_i求解：

4)利用步骤3)的计算结果，用隐式欧拉法求解x_i的值：

x_i(t+Δt)＝x_i(t)+Δt·v_i(t+Δt)

5)计算质点i邻域内其他质点的受力和运动状态；

6)更新并绘制当前模型；

7)重新计算质点i的合力，跳至步骤3)。

本实施例的非线性各向异性质量弹簧模型的软组织变形模拟方法，一方面通过对质点弹簧模型中的弹簧初始参数的设置,软组织形变过程中弹簧参数的变化以及各向异性参数化的分析，提高了质点弹簧模型的精度；另一方面，用显示欧垃积分对v_i进行推进，而采用隐式欧拉积分对x_i进行推进。本发明的技术方案与传统显示欧垃积分相比，既没有增加求解计算量，又因为在计算x_i(t+Δt)时使用了未来时刻物理量v_i(t+Δt)，所以它具有较显示欧拉积分法更好的稳定性。

以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应该以权利要求的保护范围为准。

Claims (7)

1.一种非线性各向异性质量弹簧模型的软组织变形模拟方法，其特征在于包括如下步骤；

步骤1)构建质点弹簧模型，通过有限元法确定质点弹簧模型的初始弹性参数；

步骤2)根据肌肉组织的应力、应变关系设置变形过程中弹簧参数；

步骤3)通过构建弹簧的力学方程将软组织的各向异性特征进行参数化的描述；

步骤4)构建质点弹簧模型的动力学方程并求解。

2.根据权利要求1所述的非线性各向异性质量弹簧模型的软组织变形模拟方法，其特征在于所述步骤1)通过对比弹簧参数与软组织力学参数之间的关系得到初始弹性参数；所述弹簧参数包括阻尼系数k_d和弹性系数k_s，软组织力学参数包括杨氏模量和泊松比。

3.根据权利要求2所述的非线性各向异性质量弹簧模型的软组织变形模拟方法，其特征在于所述步骤1)具体包括如下步骤:

步骤1-1)根据式(1)构建基于Kelvin模型的动力学方程，

其中m表示质量，X是位移，F_e是施加在每个质点上的力，k_d和k_s分别表示阻尼系数和弹性系数；

步骤1-2)设定Kelvin模型与有限元的链接单元都含有设于模型两端的两个节点，分别标记为第一节点与第二节点，设定F₁和F₂分别对应为施加在弹簧的第一节点和第二节点上的外力，令L和ΔL表示弹簧的原始长度和长度变化；

根据质点弹簧模型得到F₁的表达式如式(2)：

结合式(2)对于限元链接单元得到F₁的表达式如式(3)：

式(3)中，A是有限元链接单元的横截面积，σ是应力，通过式(5)计算得出：

σ＝Eε (5)

式(5)中，E是杨氏模量，ε是应变，且ε＝ΔL/L；

根据式(1)与式(3)可得弹簧的初始弹性参数如式(5)

4.根据权利要求3所述的非线性各向异性质量弹簧模型的软组织变形模拟方法，其特征在于所述步骤2)具体包括如下步骤：

步骤2-1)设定应变和应力之间的函数关系如式(6)：

σ＝ae^bε (6)

其中，a和b是系数，通过最小二乘拟合算法进行估算；

步骤2-2)根据杨氏模量和式(6)，通过式(7)计算E：

进而根据式(8)得到变形过程中的弹簧参数：

5.根据权利要求4所述的非线性各向异性质量弹簧模型的软组织变形模拟方法，其特征在于所述步骤3)包括如下步骤：

步骤3-1)根据式(8)计算弹簧在轴向和径向的弹性系数和

步骤3-2)根据式(9)计算施加在弹簧上的力的轴向分量F^α和径向分量F^β：

其中，ΔL^α和ΔL^β分别表示ΔL的轴向分量和径向分量；

步骤3-3)根据式(10)计算弹簧上的作用力F，结合式(9)形成软组织的各向异性特征的表达：

F＝F^αcosθ+F^βsinθ (10)

其中，θ弹表示簧的轴向与软组织的轴向之间的夹角。

6.根据权利要求5所述的非线性各向异性质量弹簧模型的软组织变形模拟方法，其特征在于所述步骤4)中构建质点弹簧模型的动力学方程的过程中设定质点弹簧模型基于四面体拓扑结构，并且由n个质点组成的，在任意时刻，模型的状态都由每个质点的位置x_i、速度v_i确定，任意质点i满足式(11)，

其中，m_i是质点i的质量，x_i表示质点i的位置矢量，f_in表示与质点i相连的所有弹簧施加于质点i的内力，由弹性力f_s和阻尼力f_d组成，f_s、f_d分别根据式(12)、式(13)所示，

f_d＝k_d(v_j-v_i) (13)

结合式(12)、式(13)得到式(14)形成质点弹簧模型的动力学方程：

其中，|x_j-x_i|表示表示两个质点间位置矢量之差，|x_j-x_i|⁰是发生形变前弹簧的原长，k_s为弹簧的弹性系数，k_d为弹簧的阻尼系数。

7.根据权利要求6所述的非线性各向异性质量弹簧模型的软组织变形模拟方法，其特征在于所述步骤4)中对所述质点弹簧模型的动力学方程基于改进的欧拉法求解对所述动力学方程的速度与位移矢量进行求解，如式(15)，

其中，t为初始时间点，Δt为预先定义的迭代时间步长，k表示不同的空间点；所述改进的欧拉法包括如下步骤：

步骤4-1)采用链表形式存放软组织数据，构建质点弹簧体模型，同时设置模型的初始状态为静态，并求出各个质点间弹簧的初始长度L₀；

步骤4-2)用外接力反馈设备为软组织施加作用力，并返回发生碰撞的质点的序列号i；

步骤4-3)结合式(14)计算质点i的合力，并通过显式欧拉法v_i求解，如式(16)：

步骤4-4)根据步骤4-3)的计算结果，根据式(17)用隐式欧拉法求解x_i的值：

x_i(t+Δt)＝x_i(t)+Δt·v_i(t+Δt) (17)

步骤4-5)计算质点i邻域内其他质点的受力和运动状态；

步骤4-6)更新并绘制当前模型；

步骤4-7)重新计算质点i的合力，跳至步骤4-3)。