CN109344500A - 一种非弹性质量弹簧模型及改进欧拉算法的软组织变形模拟方法 - Google Patents

Abstract

本发明涉及一种非弹性质量弹簧模型及改进欧拉算法的软组织变形模拟方法，属于软组织形变模型建模技术领域。先针对需要模拟形变的软组织，给出质点弹簧模型和一个与之相似的有限元模型，在相同受力情况下，将质点弹簧模型与有限元模型进行比较，通过有限元法确定质点弹簧模型初始参数，再对处于临界应变后变形过程中的软组织调整质点弹簧模型参数，然后用改进的欧拉法求解质点的速度、位移矢量及所受合力。本发明通过对质点弹簧模型中的弹簧初始参数的设置及软组织形变过程中弹簧参数的变化的分析，提高了质点弹簧模型的精度；通过改进的欧拉法，用显式欧拉法对v_i进行迭代求解，用隐式欧拉法对x_i进行求解，减少了计算量，保证了精度和稳定性。

Description

一种非弹性质量弹簧模型及改进欧拉算法的软组织变形模拟 方法

技术领域

本发明涉及一种非弹性质量弹簧模型及改进欧拉算法的软组织变形模拟方法，属于软组织形变模型建模技术领域。

背景技术

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

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

现有质点弹簧模型几乎都设置相同的参数，而软组织通常具有粘弹性、非线性、各向异性等复杂特征，为所有弹簧模型设置相同的参数缺乏理论基础。另外，软组织变形过程中弹簧的不变特性在极端情况下并不成立。本发明针对现有技术的不足，着力解决质点弹簧模型中弹簧参数的设置，提出了一种利用有限元法确定质点弹簧参数的改进方法，同时在解决质点弹簧模型的动力学方程中采用了改进的欧拉算法。

发明内容

本发明要解决的技术问题是针对现有技术质点弹簧模型参数设置不合理、精度不高或计算复杂等不足，提出一种非弹性质量弹簧模型及改进欧拉算法的软组织变形模拟方法，解决质点弹簧模型中弹簧参数的设置，在保证精度和稳定性的同时，简化计算、提高效率。

为解决上述技术问题，本发明提供一种非弹性质量弹簧模型及改进欧拉算法的软组织变形模拟方法，先针对需要模拟形变的软组织，给出质点弹簧模型和一个与之相似的有限元模型，在相同受力情况下，将质点弹簧模型与有限元模型进行比较，通过有限元法确定质点弹簧模型初始参数，再对处于临界应变后变形过程中的软组织调整质点弹簧模型参数，然后用改进的欧拉法求解质点的速度、位移矢量及所受合力。

所述非弹性质量弹簧模型及改进欧拉算法的软组织变形模拟方法的具体步骤如下：

(1)针对需要模拟形变的软组织，给出Kelvin质点弹簧模型和一个与之相似的有限元模型，分别将Kelvin质点弹簧模型与有限元模型的两个节点链接单元标记为1，2，位置坐标为在软组织弹簧的节点1和节点2上施加外力F₁和F₂，在相同受力情况下，将Kelvin质点弹簧模型与有限元模型进行比较，通过有限元法，按下式确定质点弹簧模型初始参数：

其中，k_s是软组织弹性系数，E是软组织的杨氏模量，A是连接单元的横截面积，是质点弹簧模型与有限元模型节点1和节点2之间的原始距离；

(2)对于处于临界应变后变形过程中的软组织弹簧，取P₀(ε₀，σ₀)和P₁(ε₁，σ₁)两个点，P₀为比例线与指数曲线的连接点，P₁为指数曲线上的另一点，按下式调整质点弹簧模型参数：

其中，σ是应力，σ₀是临界应力，σ₁是实际应力，ε是应变，ε₀是临界应变，ε₀≈0.25，ε₁是实际应变；当应变低于ε₀时，肌肉组织的应力与应变成正比，肌肉组织的应力与应变呈比例关系(比例线)；当应变高于ε₀时，肌肉组织的应力与应变呈指数函数关系(指数曲线)；

(3)按以下步骤求解动力学方程

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

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

(3.3)按下式改进的欧拉法求解质点i的速度与位移矢量：

其中，Δt是预先定义的迭代时间步，k代表不同的空间点，m_i是质点i的质量，x_i表示质点i的位置矢量，v_i是质点i的速度；

先用显式欧拉法求解v_i：

再用隐式欧拉法求解x_i：

(3.4)按下式计算质点i的合力f_ext：

其中，f_ext表示质点i所受合力，包括摩擦力、重力、用户施加的力等，m_i是质点i的质量，x_i表示质点i的位置矢量，x_j表示质点j的位置矢量，|x_j-x_i|是两个质点间位置矢量之差， |x_j-x_i|⁰是发生形变前弹簧的原长，k_s为弹簧的弹性系数，k_d为弹簧的阻尼系数，v_i是质点i 的速度，v_j是质点j的速度；

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

(3.6)更新并绘制当前模型；

(3.7)跳至步骤(3.3)，重新计算质点i的位置矢量x_i、速度v_i、合力f_ext。

本发明的原理是：

(一)弹簧初始参数的设置

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

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

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

为了达到与有限元模型相同的变形量，需要找出质点弹簧模型与有限元模型的关系，即弹簧参数(包括阻尼系数k_d和弹性系数k_s)与软组织力学参数(包括杨氏模量和泊松比)之间的关系，以及有限元单元的几何特征(包括节点的坐标、面积、厚度等)。如图2所示，给出Kelvin模型与有限元的链接单元，将它们进行比较可知，Kelvin模型与有限元的链接单元都有两个节点，分别标记为1，2。位置坐标为F₁和F₂是施加在弹簧的节点1和节点2上的外力，包括重力以及与节点相连的其它弹簧的弹簧力。根据公式(1)可得到弹簧的阻尼力，如公式(2)所示：

其中，分别表示节点1，节点2的速度；是从节点1到节点2的矢量。弹簧的弹力由公式(3)给出：

其中，代表节点1和节点2之间的原始距离。

对于有限元的链接单元的两个节点的动力学方程，可以由下面的公式(4)和公式(5) 表示：

其中，A是连接单元的横截面积，σ是应力，ε是应变，并且E是软组织的杨氏模量。为了得到相同的变形，将公式(1)与公式(4)进行比较，得到公式(6)：

再将公式(2)，(3)，(5)插入公式(6)，得到公式(7)：

消除相同的项，得到公式(8)：

方程(8)显示了弹簧参数与软组织力学参数之间的关系，以及有限元单元的几何特征。此外，由于阻尼系数k_d通常远小于弹性系数k_s，因此可以忽略不计，即将阻尼系数k_d设为零，最终得到弹性系数k_s如公式(9)所示：

(二)变形过程中弹簧参数的变化

在大多数质点弹簧模型中，弹簧的弹性特征在软组织变形过程中保持不变。但是，根据材料的真实力学性能可知，这是不正确的。如图3所示的是肌肉组织的应力-应变曲线。从图中可以看出，当应变低于ε₀(ε₀≈0.25)时，肌肉组织的应力与应变成正比，当应变高于ε₀时，应力与应变的函数关系呈指数型。假设应变和应力之间的函数关系如公式(10)所示：

在图3中，取两点P₀(ε₀，σ₀)，P₁(ε₁，σ₁)，其中P₀为比例线与指数曲线的连接点，P₁为指数曲线上的另一点，则可以得到公式(11)：

将公式(10)，(11)插入公式(9)中，则弹性系数k_s如公式(12)所示：

(三)动力学方程及求解

图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是预先定义的迭代时间步，k代表不同的空间点。该算法的具体步骤如下：

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

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

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

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

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

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

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

本发明一方面通过对质点弹簧模型中的弹簧初始参数的设置以及软组织形变过程中弹簧参数的变化的分析，提高了质点弹簧模型的精度；另一方面改进的欧拉法用显式欧拉法对v_i进行迭代求解，并直接使用v_i的计算结果，用隐式欧拉法对x_i进行求解。与单纯的隐式欧拉法相比，减少了计算量，同时保证了隐式欧拉法的精度和稳定性。

附图说明

图1是本发明Kelvin模型示意图。

图2是本发明Kelvin质点弹簧模型的链接单元示意图。

图3是本发明与质点弹簧模型对应的有限元模型的链接单元示意图。

图4是本发明肌肉组织的应力-应变曲线图。

图5是本发明实施例基于四面体拓扑结构的质点-弹簧体模型图。

具体实施方式

下面结合附图对本发明的具体实施方式作进一步详尽描述，实施例中未注明的技术或产品，均为现有技术或可以通过购买获得的常规产品。

实施例1：如图1-5所示，本非弹性质量弹簧模型及改进欧拉算法的软组织变形模拟方法，先针对需要模拟形变的软组织，给出质点弹簧模型和一个与之相似的有限元模型，在相同受力情况下，将质点弹簧模型与有限元模型进行比较，通过有限元法确定质点弹簧模型初始参数，再对处于临界应变后变形过程中的软组织调整质点弹簧模型参数，然后用改进的欧拉法求解质点的速度、位移矢量及所受合力。具体步骤如下：

(1)针对需要模拟形变的软组织，给出Kelvin质点弹簧模型和一个与之相似的有限元模型，分别将Kelvin质点弹簧模型与有限元模型的两个节点链接单元标记为1，2，位置坐标为在软组织弹簧的节点1和节点2上施加外力F₁和F₂，在相同受力情况下，将Kelvin质点弹簧模型与有限元模型进行比较，通过有限元法，按下式确定质点弹簧模型初始参数：

其中，k_s是软组织弹性系数，E是软组织的杨氏模量，A是连接单元的横截面积，是质点弹簧模型与有限元模型节点1和节点2之间的原始距离；

(2)对于处于临界应变后变形过程中的软组织弹簧，取P₀(ε₀，σ₀)和P₁(ε₁，σ₁)两个点，P₀为比例线与指数曲线的连接点，P₁为指数曲线上的另一点，按下式调整质点弹簧模型参数：

其中，σ是应力，σ₀是临界应力，σ₁是实际应力，ε是应变，ε₀是临界应变，ε₀≈0.25，ε₁是实际应变；当应变低于ε₀时，肌肉组织的应力与应变成正比，肌肉组织的应力与应变呈比例关系(比例线)；当应变高于ε₀时，肌肉组织的应力与应变呈指数函数关系(指数曲线)；

(3)按以下步骤求解动力学方程

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

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

(3.3)按下式改进的欧拉法求解质点i的速度与位移矢量：

其中，Δt是预先定义的迭代时间步，k代表不同的空间点，m_i是质点i的质量，x_i表示质点i的位置矢量，v_i是质点i的速度；

先用显式欧拉法求解v_i：

再用隐式欧拉法求解x_i：

(3.4)按下式计算质点i的合力f_ext：

其中，f_ext表示质点i所受合力，包括摩擦力、重力、用户施加的力等，m_i是质点i的质量，x_i表示质点i的位置矢量，x_j表示质点j的位置矢量，|x_j-x_i|是两个质点间位置矢量之差， |x_j-x_i|⁰是发生形变前弹簧的原长，k_s为弹簧的弹性系数，k_d为弹簧的阻尼系数，v_i是质点i 的速度，v_j是质点j的速度；

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

(3.6)更新并绘制当前模型；

(3.7)跳至步骤(3.3)，重新计算质点i的位置矢量x_i、速度v_i、合力f_ext。

上面结合附图对本发明的技术内容作了说明，但本发明的保护范围并不限于所述内容，在本领域的普通技术人员所具备的知识范围内，还可以在不脱离本发明宗旨的前提下对本发明的技术内容做出各种变化，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (2)

1.一种非弹性质量弹簧模型及改进欧拉算法的软组织变形模拟方法，其特征在于：先针对需要模拟形变的软组织，给出质点弹簧模型和一个与之相似的有限元模型，在相同受力情况下，将质点弹簧模型与有限元模型进行比较，通过有限元法确定质点弹簧模型初始参数，再对处于临界应变后变形过程中的软组织调整质点弹簧模型参数，然后用改进的欧拉法求解质点的速度、位移矢量及所受合力。

2.根据权利要求1所述的非弹性质量弹簧模型及改进欧拉算法的软组织变形模拟方法，其特征在于：所述软组织变形模拟方法的具体步骤如下：

(1)针对需要模拟形变的软组织，给出Kelvin质点弹簧模型和一个与之相似的有限元模型，分别将Kelvin质点弹簧模型与有限元模型的两个节点链接单元标记为1，2，位置坐标为在软组织弹簧的节点1和节点2上施加外力F₁和F₂，在相同受力情况下，将Kelvin质点弹簧模型与有限元模型进行比较，通过有限元法，按下式确定质点弹簧模型初始参数：

其中，k_s是软组织弹性系数，E是软组织的杨氏模量，A是连接单元的横截面积，是质点弹簧模型与有限元模型节点1和节点2之间的原始距离；

(2)对于处于临界应变后变形过程中的软组织弹簧，取P₀(ε₀，σ₀)和P₁(ε₁，σ₁)两个点，P₀为比例线与指数曲线的连接点，P₁为指数曲线上的另一点，按下式调整质点弹簧模型参数：

其中，σ是应力，σ₀是临界应力，σ₁是实际应力，ε是应变，ε₀是临界应变，ε₀≈0.25，ε₁是实际应变；当应变低于ε₀时，肌肉组织的应力与应变成正比，肌肉组织的应力与应变呈比例关系(比例线)；当应变高于ε₀时，肌肉组织的应力与应变呈指数函数关系(指数曲线)；

(3)按以下步骤求解动力学方程

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

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

(3.3)按下式改进的欧拉法求解质点i的速度与位移矢量：

其中，Δt是预先定义的迭代时间步，k代表不同的空间点，m_i是质点i的质量，x_i表示质点i的位置矢量，v_i是质点i的速度；

先用显式欧拉法求解v_i：

再用隐式欧拉法求解x_i：

(3.4)按下式计算质点i的合力f_ext：

其中，f_ext表示质点i所受合力，m_i是质点i的质量，x_i表示质点i的位置矢量，x_j表示质点j的位置矢量，|x_j-x_i|是两个质点间位置矢量之差，|x_j-x_i|⁰是发生形变前弹簧的原长，k_s为弹簧的弹性系数，k_d为弹簧的阻尼系数，v_i是质点i的速度，v_j是质点j的速度；

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

(3.6)更新并绘制当前模型；

(3.7)跳至步骤(3.3)，重新计算质点i的位置矢量x_i、速度v_i、合力f_ext。