1.一种质点弹簧模型软组织变形模拟方法,其特征在于,所述方法包括:S1:利用有限元法确定质点弹簧模型的参数,包括质点弹簧模型初始参数设置、软组织变形过程中质点弹簧模型弹簧参数的确定以及软组织各项异性特征的参数化;
S2:采用改进的Verlet算法计算质点弹簧模型的质点位置和速度;
S3:采用轴向包围盒AABB和方向包围盒OBB相结合的混合包围盒算法进行碰撞检测,以判断虚拟对象与包围盒之间是否发生碰撞,所述虚拟对象包括非刚体对象和刚体对象,其中,非刚体对象采用轴向包围盒AABB,刚体对象采用方向包围盒OBB;
所述采用改进的Verlet算法计算质点弹簧模型的质点位置和速度的过程包括以下步骤:S20:推导质点弹簧模型中质点受力、速度、加速度的动力学表达方程式,推导过程包括以下步骤:S201:定义质点弹簧模型中每一质点所受的力包括内力和外力,所述内力包括弹簧形变产生的弹力fs和质点运动产生的阻尼力fd,所述外力包括作用在质点上的强迫力fe;
S202:采用如下公式表示质点所受合力F:
F=fe‑fs‑fd;
S203:按照胡克定律,采用如下公式表示质点i所受弹力:其中,kij表示质点i和j之间弹簧的弹性系数,lij=Xi‑Xj表示质点i和j之间矢量距离,表示质点i和j之间弹簧的原始长度,lij/||lij||为单位向量;
S204:采用如下公式表示质点i所受的阻尼力:其中,ηij表示质点i和j之间的阻尼系数, 表示质点i和j之间的速度差值;
S205:采用牛顿第二定律,得到质点i的动力学方程为:miai=fe‑fs‑fd
其中,mi、ai为质点i的质量和加速度;
S206:结合步骤S203中的fs表达式、步骤S204中的fd表达式,得到质点i受力、速度、加速度的表达式:S21:定义x0、v0、h、F0、a0分别代表质点弹簧模型中质点的初始位置、速度、步长时间、所受合力和加速度,设定质点受力变化的经验阈值为μ;
S22:采用2阶Taylor级数法对改进的Verlet算法进行首次迭代计算,得到质点首次迭代计算后的位移、速度分别x1、v1,其中:2
x1=x0+h·v0+a0·h/2
v1=v0+a0·h;
S23:定义正整数k,令k=1;
S24:采用步骤S206中所述的受力公式计算并更新质点的受力Fk值;
S25:令ΔF=Fk‑Fk‑1,将ΔF的值与预设质点受力变化的经验阈值μ进行比较,若ΔF>μ,跳转至步骤S26,否则,跳转至步骤27;
S26:采用Beeman算法完成质点所受合外力改变较大情况下,质点位置和速度的迭代计算,得到位置和速度的值:2
xk+1=xk+hvk+h(4ak‑ak‑1)/6vk+1=vk+h(2ak+1+5ak‑ak‑1)/6其中,xk+1、xk分别表示质点在下一时刻和当前时刻的空间位置,vk+1、vk分别表示质点在下一时刻和当前时刻的速度,ak+1、ak、ak‑1分别表示质点在下一时刻、当前时刻和前一时刻的加速度,h表示步长时间,跳转至步骤28;
S27:采用常规Verlet算法完成质点所受合外力改变较小情况下,质点位置和速度的迭代计算,得到位置和速度的值:2
xk+1=2xk‑xk‑1+ak·h
vk+1=(xk+1‑xk‑1)/(2h)
其中,xk+1、xk、xk‑1分别表示质点在下一时刻、当前时刻和前一时刻的空间位置,vk+1表示质点在下一时刻的速度,ak表示质点在当前时刻的加速度,h表示步长时间,跳转至步骤28;
S28:处理弹簧变形约束,更新质点的位置、速度值,并令k加1;
S29:重复步骤S24至步骤S28,直至接收到模拟结束命令。
2.根据权利要求1所述的基于改进Verlet算法的质点弹簧模型软组织变形模拟方法,其特征在于,步骤S1中,所述质点弹簧模型初始参数设置的过程包括以下步骤:S101:选取Kelvin模型作为质点弹簧模型,Kelvin模型由弹性元件和阻尼元件组成,采用如下拉格朗日运动方程定义Kelvin模型的受力与质量、位移的关系:其中,m表示质点的质量,X表示质点的位移,Fe表示施加在每个质点上的力,kd表示阻尼系数,ks表示弹性系数;
略去阻尼系数相关的一项,得到如下简化公式:
S102:定义一相同受力的有限元模型,该有限元模型与质点弹簧Kelvin模型的变形量相同,计算有限元模型与质点弹簧Kelvin模型的参数关系,有限元模型与质点弹簧Kelvin模型的参数关系包括弹簧参数与软组织力学参数之间的关系、有限元模型的单元几何特征。
3.根据权利要求2所述的一种质点弹簧模型软组织变形模拟方法,其特征在于,步骤S102中,所述计算有限元模型与质点弹簧Kelvin模型的参数关系的过程包括以下步骤:S102a:定义Kelvin模型与有限元模型链接单元上的两个节点分别为节点1、节点2,F1、F2为施加在节点1、节点2上的外力,所述外力包括重力、与节点相连的其它弹簧的弹簧力,定义L和ΔL表示弹簧的初始长度和长度变化;
S102b:采用如下公式表示质点弹簧Kelvin模型在仅有弹力情况下的受力F1:其中,m表示质点的质量,X表示质点的位移,ks表示弹性系数,ΔL表示弹簧的长度变化;
采用如下公式表示有限元模型的链接单元中节点1的动力学方程为:其中,A表示有限元模型的链接单元横截面积,σ表示应力,σ=Eε,E表示杨氏模量,ε表示应变,ε=ΔL/L;
S102c:步骤S102b中的质点弹簧Kelvin模型与有限元模型的节点1受力F1相同,结合二个表达式得到:ksΔL=Aσ
再结合应力、应变表达式,得到质点弹簧模型的初始弹性参数ks为:
4.根据权利要求1所述的一种质点弹簧模型软组织变形模拟方法,其特征在于,步骤S1中,所述软组织变形过程中质点弹簧模型弹簧参数的确定过程包括以下步骤:S111:采用如下公式表示软组织在变形过程中遵守的应力‑应变曲线关系:bε
σ=ae
其中,a和b表示系数,采用最小二乘拟合算法估计得到,σ表示应力,ε表示应变;
S112:将步骤S111中得到的应力‑应变关系,代入杨氏模量定义公式,可得:其中,E表示杨氏模量,σ表示应力,ε表示应变,a和b表示系数;
S113:将杨氏模量表达式代入步骤S102c中质点弹簧模型的初始弹性参数公式,得到质点弹簧模型在软组织变形过程中的弹性参数ks'为:其中,ks'表示质点弹簧模型在软组织变形过程中的弹性系数,A表示有限元模型的链接单元横截面积,L表示弹簧的初始长度,ε表示应变,a和b表示系数。
5.根据权利要求1所述的一种质点弹簧模型软组织变形模拟方法,其特征在于,步骤S1中,所述软组织各项异性特征的参数化的过程包括以下步骤:S121:根据软组织的各项异性特征,定义质点弹簧模型在轴向和径向的弹性系数分别为 采用如下公式表示:α β
其中,A表示有限元模型的链接单元在轴向的横截面积、A 表示有限元模型的链接单元α β α在径向的横截面积,L表示弹簧在轴向的初始长度、L表示弹簧在径向的初始长度,ε表示β α β α β轴向的应变、ε表示径向的应变,a、a、b、b表示系数;
α β
S122:定义质点弹簧模型在轴向和径向的受力分别为F、F,采用如下公式表示:α β
其中,ΔL和ΔL分别表示ΔL的轴向变化分量和径向变化分量;
S123:定义质点弹簧模型的弹簧轴向与软组织轴向的夹角为θ,采用如下公式表示弹簧上的总作用力F,用以实现软组织各项异性特征的参数化:α β
F=Fcosθ+Fsinθ。
6.根据权利要求1所述的一种质点弹簧模型软组织变形模拟方法,其特征在于,所述采用轴向包围盒AABB和方向包围盒OBB相结合的混合包围盒算法进行碰撞检测,以判断虚拟对象与包围盒之间是否发生碰撞,所述虚拟对象包括非刚体对象和刚体对象,其中,非刚体对象采用轴向包围盒AABB,刚体对象采用方向包围盒OBB的过程包括以下步骤:S31:创建虚拟对象模型,构造树状层次结构的轴向包围盒AABB树、方向包围盒OBB树;
S32:设置虚拟对象位置,判断虚拟对象是刚体对象或者非刚体对象;
S33:根据非刚体对象采用轴向包围盒AABB,刚体对象采用方向包围盒OBB的原则对虚拟对象进行碰撞检测,判断轴向包围盒AABB、方向包围盒OBB是否碰撞,若发生碰撞,遍历轴向包围盒AABB、方向包围盒OBB判断是否相交,进而得出虚拟对象的碰撞信息,跳转至步骤S34;
若未发生碰撞,直接跳转至步骤S34;
S34:重复步骤S32至S33,直至收到结束指令。