[发明专利]一种软组织形变仿真方法

权利要求书

1.一种软组织形变仿真方法和技术，其特征在于包含以下步骤：

步骤1)：采集软组织的数据信息；

步骤2)：选择黏弹性模型，构建用于软组织形变仿真的方程组：

采用一个弹簧和一个黏壶并联的Kelvin黏弹性模型；

首先，构建三维格式的Kelvin黏弹性应力-应变本构方程：

依据Kelvin模型，各向同性材料的应力张量σ可分解成它的球形张量和偏斜张量部分，应变张量ε可分离为体积形变和等体积的形状畸变两部分：

σαβ=Sαβ+δαβσkk3---(1)]]>

ϵαβ=eαβ+δαβϵkk3---(2)]]>

其中，α，β＝x，y，z.σ_αβ为应力分量，ε_αβ为应变分量；δ_αβ为Kronecker符号，σ_kk＝σ_xx+σ_yy+σ_zz和ε_kk＝ε_xx+ε_yy+ε_zz分别为体积应力和体积应变；S_αβ和e_αβ分别为偏应力张量和偏应变张量的分量；

根据Kelvin模型，偏应力张量和偏应变张量之间、体积应力和体积应变之间的三维黏弹性本构关系可表示为：

Sαβ=E·eαβ+η·deαβdt,α,β=x,y,z]]>

σkk=E·ϵkk+η·dϵkkdt---(3)]]>

上式中，S_αβ和e_αβ分别为偏应力张量和偏应变张量的分量；σ_kk和ε_kk分别为体积应力和体积应变；E为材料的弹性模量，η为黏性系数；分别为偏应变分量和体积应变对时间的导数，即应变率，t为时间；

其次，构建应变-位移几何方程：

ϵxx=∂u∂x]]>ϵxy=∂u∂y+∂v∂x]]>

ϵyy=∂v∂y]]>ϵyz=∂v∂z+∂w∂y---(4)]]>

ϵzz=∂w∂z]]>ϵzx=∂w∂x+∂u∂z]]>

其中u，v，w为位移在三个坐标方向的分量；ε_xx，ε_yy，ε_zz为正应变分量；ε_xy，ε_yz，ε_zx为剪应变分量；

然后，构建用于加速度计算的动量方程，其形式如下：

dvdt=1ρ∂σ∂x---(5)]]>

其中，v为速度向量，t为时间，ρ为粒子的密度，x为坐标向量；

步骤3)：依据步骤1)采集的软组织的数据信息，利用各数据点的位置向量，构建一个没有网格连接的粒子模型，并初始化所述粒子模型中各粒子的位置、质量、速度、加速度、受作用力等信息，构建粒子模型的初始化状态；

步骤4)：定义空间栅格，采用链表搜索法搜索粒子的支持域，构建支持域内的光滑核函数；

核函数选取三次B-Spline光滑函数：

Wij=1πh3×1-32s2+34s30≤s<114(2-s)31≤s<20s≥2---(6)]]>

式中，W_ij为由邻近粒子p_j估计粒子p_i运动信息的光滑核函数，为粒子p_i与p_j之间的相对距离，r为粒子p_i与p_j之间的距离，h为光滑长度；

对核函数在各方向求导，得到核函数的一阶导数为：

∂Wij∂xβ=xiβ-xjβr·1πh4×-3s+94s20≤s<1-34(2-s)21≤s<20s≥2---(7)]]>

式中，上标β＝x，y，z表示坐标方向，s，r，h的含义同式(6)，和分别表示粒子p_i与p_j的位置坐标向量在各方向的分量；

步骤5)：应用光滑核函数W及其导数对参考粒子的支持域内所有粒子函数的加权平均近似的方法来构建步骤2)中各方程的SPH格式：

构建密度方程SPH格式

ρi=Σj=1NmjWij---(8)]]>

其中，ρ_i为粒子p_i的密度，m_j为p_i的支持域内邻近粒子p_j的质量，N为粒子p_i的支持域内粒子总数，W_ij的含义同式(6)；

将动量方程(5)转换为SPH格式

DvixDt=Σj=1N[mj(σixxρi2+σjxxρj2)∂Wij∂x+mj(σixyρi2+σjxyρj2)∂Wij∂y+mj(σixzρi2+σjxzρj2)∂Wij∂z]]]>

DviyDt=Σj=1N[mj(σiyyρi2+σjyyρj2)∂Wij∂y+mj(σiyxρi2+σjyxρj2)∂Wij∂x+mj(σiyzρi2+σjyzρj2)∂Wij∂z]---(9)]]>

DvizDt=Σj=1N[mj(σizzρi2+σjzzρj2)∂Wij∂z+mj(σizxρi2+σjzxρj2)∂Wij∂x+mj(σizyρi2+σjzyρj2)∂Wij∂y]]]>

其中，分别为粒子p_i的运动速度在各坐标方向的分量；为粒子p_i应力向量的正应力分量，为粒子p_i应力向量的剪应力分量；为粒子p_j应力向量的正应力分量，为粒子p_j应力向量的剪应力分量；ρ_i，ρ_j分别为粒子p_i和p_j的密度；m_j，W_ij，N的含义同公式(8)；

构建几何方程的SPH格式：

ϵixx=∂ui∂x=Σj=1Nmjρj·uj·∂Wij∂x]]>

ϵiyy=∂vi∂y=Σj=1Nmjρj·vj·∂Wij∂y]]>

ϵizz=∂wi∂z=Σj=1Nmjρj·wj·∂Wij∂z---(10)]]>

ϵixy=∂ui∂y+∂vi∂x=Σj=1Nmjρj·uj·∂Wij∂y+Σj=1Nmjρj·vj·∂Wij∂x]]>

ϵiyz=∂vi∂z+∂wi∂y=Σj=1Nmjρj·vj·∂Wij∂z+Σj=1Nmjρj·wj·∂Wij∂y]]>

ϵizx=∂wi∂x+∂ui∂z=Σj=1Nmjρj·wj·∂Wij∂x+Σj=1Nmjρj·uj·∂Wij∂z]]>

其中，为粒子p_i应变向量的正应变分量，为粒子p_i应变向量的剪应变分量；u_i，v_i，w_i为粒子p_i的运动位移在各坐标方向的分量，u_j，v_j，w_j为粒子p_j的运动位移在各坐标方向的分量；ρ_j为粒子p_j的密度，m_j，W_ij，N的含义同式(8)；

步骤6)：用显示积分法求解步骤5)的常微分方程，计算出粒子的密度、位置、速度等随时间的变化值；

步骤7)：循环执行上述步骤3)～步骤6)，计算出各粒子的状态；

步骤8)：将粒子模型的当前状态输出之屏幕，经过纹理和光照渲染后，得到软组织的动态形变过程。

2.根据权利要求1的软组织形变仿真方法，其特征在于步骤6)的求解过程如下：

1)对模型的某一质点p施加外力f_outp；

2)对模型中的各粒子p_i(i＝1，2，..，M)循环完成下述3)～7)步骤的运算处理；

3)对当前粒子p_i，以h为光滑半径，搜索支持域内的相邻粒子p_j(j＝1，2，..，N)；

4)使用公式(6)(7)，计算当前粒子p_i与支持域内各邻近粒子p_j(j＝1，2，..，N)间的光滑核函数W_ij及其导数；

5)使用方程(8)，采用密度求和法计算粒子的密度ρ_i；

6)计算粒子的加速度a_i，加速度的计算采用如下方法：

ai=-DviDt+σimi+Fimi]]>

其中，F_i表示粒子p_i所受的重力、外力等，σ_i表示粒子p_i所受的内力，即应力，m_i为粒子p_i的质量，使用公式(9)计算

7)计算粒子的速度和位移：

v_i(t)＝v_i(t)+dt·a_i

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

其中，v_i(t)为粒子p_i在t时刻的速度向量，x_i(t)为粒子p_i在t时刻的位置向量，a_i为粒子p_i的加速度，dt为时间增量；

8)计算各粒子与初始状态时的位移情况disp_i(t)＝x_i(t)-x_i(t₀)；

9)使用几何方程(4)和(10)，由上一步算出的位移计算各粒子的应变状态

10)记录各粒子当前的体积应变和形状畸变并使用方程(2)，由上一步计算出的应变计算各粒子的新的体积应变和形状畸变

11)计算并使用黏弹性应力-应变本构方程(3)，计算各粒子的体积应力和偏应力

12)使用方程(1)，计算出各粒子的应力状态