CN108694290B - 一种基于八叉树网格的有限元模型的软组织变形方法 - Google Patents
一种基于八叉树网格的有限元模型的软组织变形方法 Download PDFInfo
- Publication number
- CN108694290B CN108694290B CN201810566190.1A CN201810566190A CN108694290B CN 108694290 B CN108694290 B CN 108694290B CN 201810566190 A CN201810566190 A CN 201810566190A CN 108694290 B CN108694290 B CN 108694290B
- Authority
- CN
- China
- Prior art keywords
- matrix
- soft tissue
- hexahedron
- point
- deformation
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Expired - Fee Related
Links
- 238000000034 method Methods 0.000 title claims abstract description 63
- 210000004872 soft tissue Anatomy 0.000 title claims abstract description 54
- 239000011159 matrix material Substances 0.000 claims abstract description 102
- 238000006073 displacement reaction Methods 0.000 claims abstract description 28
- 230000010354 integration Effects 0.000 claims abstract description 10
- 230000008569 process Effects 0.000 claims abstract description 8
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 6
- 238000009877 rendering Methods 0.000 claims abstract description 5
- 238000013016 damping Methods 0.000 claims description 33
- 239000000463 material Substances 0.000 claims description 21
- 230000014509 gene expression Effects 0.000 claims description 16
- 238000007781 pre-processing Methods 0.000 claims description 12
- 239000006185 dispersion Substances 0.000 claims description 9
- 238000004364 calculation method Methods 0.000 claims description 8
- 238000013334 tissue model Methods 0.000 claims description 6
- 239000013598 vector Substances 0.000 claims description 6
- 230000008859 change Effects 0.000 claims description 4
- 230000001133 acceleration Effects 0.000 claims description 3
- 238000001914 filtration Methods 0.000 claims description 3
- 238000013507 mapping Methods 0.000 claims description 3
- 239000002245 particle Substances 0.000 claims description 3
- 239000004576 sand Substances 0.000 claims description 3
- 238000012360 testing method Methods 0.000 claims description 3
- 238000006467 substitution reaction Methods 0.000 claims description 2
- 210000002615 epidermis Anatomy 0.000 abstract description 2
- 238000010586 diagram Methods 0.000 description 8
- 238000004088 simulation Methods 0.000 description 6
- 230000008602 contraction Effects 0.000 description 4
- 241001465754 Metazoa Species 0.000 description 2
- 230000007812 deficiency Effects 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 238000012549 training Methods 0.000 description 2
- 241000269350 Anura Species 0.000 description 1
- 241000699670 Mus sp. Species 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 239000003814 drug Substances 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 238000012886 linear function Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000012805 post-processing Methods 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 238000005070 sampling Methods 0.000 description 1
- 210000001519 tissue Anatomy 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
- G06T17/20—Finite element generation, e.g. wire-frame surface description, tesselation
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Geometry (AREA)
- General Physics & Mathematics (AREA)
- Computer Graphics (AREA)
- Software Systems (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- General Engineering & Computer Science (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明提供一种基于八叉树网格的有限元模型的软组织变形方法。本发明方法,包括:绘制软组织的三维模型,基于AABB法构建多个均匀的六面体网格,在六面体网格基础上基于八叉树的网格生成算法生成八叉树网格,对六面体网格进行有限元方法建模,求解软组织的变形过程,将相邻单元节点上的单元刚度矩阵组装成离散域的总刚度矩阵,在动力平衡情况下通过时间积分对各矩阵的数值求解动力平衡微分方程,得到随时间变化的节点位移,通过对节点位移的渲染,显示软组织受力变形场景。本发明逼真的模拟虚拟手术中任意形状软组织表皮受拉力变形的过程,有很高的实时性,减少了计算量,解决了以往有限元网格数量复杂,不能实时仿真软组织变形过程的问题。
Description
技术领域
本发明涉及虚拟手术仿真领域,尤其涉及一种基于八叉树网格的有限元模型的软组织变形方法。
背景技术
传统临床医学通常使用橡胶人体模型和人类尸体、小白鼠、青蛙等生物体作为手术训练对象。橡胶人体模型结构简单、功能单一,仿真度存在误差;人类尸体的材料活性与活体组织存在较大差异,复用性低,且随着我国社会的发展,尸体的获得越来越难;动物与人体的生理结构不同,作为临床手术训练对象并不准确,同时存在着迫害动物等诸多道德问题。
随着计算机硬件处理性能的提升,虚拟手术仿真系统得到了广泛研究。而软组织变形是虚拟手术中的重要组成部分。目前常用的软组织变形方法从物理计算模型角度可以分为三类:第一类,基于有限元方法的仿真学方法,由于包含软组织的材料特性,能够逼真模拟软组织的形变,使得其结果精度高,但有限元的形变计算是基于大量的网格单元,因此计算代价很高。第二类,基于位置的动力学方法,虽然计算速度快、不存在稳定性问题,但其产生的变形效果仅仅是貌似符合物理定律,而非真正符合物理定律。第三类,无网格方法,这种方法适用于大形变的情况,但由于采样点比较密集,计算负担比较重。
发明内容
根据上述提出的技术问题,而提供一种在不降低仿真逼真度的同时减少计算量的基于八叉树网格的有限元模型的软组织变形方法。本发明采用的技术手段如下:
一种基于八叉树网格的有限元模型的软组织变形方法,包括如下步骤:
S1、绘制软组织的三维模型,基于轴对齐边界包围盒法,即AABB法,构建多个均匀的六面体,并过滤掉中心不在软组织几何模型空间的六面体;
S2、基于BON八叉树原理,创建BON八叉树,对各六面体编号,依次判断三维模型顶点是否在六面体内部,若是,则不属于边界元;若否,则表明三维模型顶点在六面体边界上,将该顶点插入对应的最底层六面体的八叉树结构中;
S3、判断三维模型顶点组成的三角网格是否共面,若是,则删除两个三角形的公共点,将原本的两个三角形转换为大三角形,将两个相邻六面体合并为中型六面体,提取合并后的中型六面体的边界点,直到所有相邻三角形不共面、不能进行收缩为止,进而得到最终的六面体网格;
S4、对最终的六面体网格进行有限元方法建模,求解软组织的变形过程,具体包括:通过应力平衡方程的积分形式与材料的本构关系相结合,建立各个有限元的节点的位移和载荷之间的表达式,从而得到单元刚度矩阵;
S5、将相邻单元节点上的单元刚度矩阵、单元阻尼矩阵、单元质量矩阵组装成离散域的总刚度矩阵、总阻尼矩阵、总质量矩阵,在动力平衡情况下对各矩阵的数值求解动力平衡微分方程,得到随时间变化的节点位移;
S6、通过对节点位移的渲染,不断更新软组织受力变形场景。
进一步地,步骤S1中,通过如下方法绘制软组织的三维模型:
S11、通过模型加载库读取模型文件的顶点坐标、法线坐标、纹理坐标和面片的相关索引信息以及模型材质库信息;
S12、通过节点之间拓扑关系以及材质信息,利用OpenGL图形库绘制出软组织的三维模型。
进一步地,所述基于AABB法生成均匀的六面体网格,具体包括如下步骤:
S13、模型任意点的三维坐标分别为x、y、z,AABB包围盒内的最小点坐标为pmin=(xmin,ymin,zmin),最大点坐标为pmax=(xmax,ymax,zmax),AABB包围盒内点满足以下条件:
xmin≤x≤xmax
ymin≤y≤ymax
zmin≤z≤zmax,
中心点O=(xo,yo,zo)为两个顶点的中点,表示包围盒的质点,其满足如下关系:
xo=(xmin+xmax)*0.5
yo=(ymin+ymax)*0.5
zo=(zmin+zmax)*0.5,
以点pmin、点pmax建立第一个最大的正六面体包围盒,设定剖分深度N,随后以每个边的中点为界,依次剖分N次,得到2的幂次方数量的六面体,通过改变N的值可以得到不同分辨率的模型;
S14、通过射线投影法,检查S13中生成的正六面体中心点是否在物体几何空间的外部,如果正六面体的中心点在物体几何空间的外部,那么把该中心点删除。
进一步地,步骤S2中,所述基于BON八叉树原理,创建BON八叉树具体为:
以openGL的右手坐标系对每个六面体进行编号,编号顺序采用先下后上,先左后右的顺序,编号方式采用0-7的二进制位000-111进行分配。
进一步地,所述步骤S3中,通过如下方法判断三维模型顶点组成的三角网格是否共面:如果相邻六面体中的两个三角形至少有一个边相邻,且法向量相同,则两个三角形共面。
进一步地,步骤S4中,所述通过应力平衡方程的积分形式与材料的本构关系相结合具体包括如下步骤:
S401、读取此时模型的位置信息;
S41、构建柯西应变张量,具体为:以X表示体中初始状态下的一点,在变形后,该点映射到了新的点x,体中每个点都会映射到新的点,这个映射表示为:
该点的位移为:
u=x-X
则变形梯度为:
变形前后长度的改变为:
||dx||2-||dX||2=dXT(FTF-I)dX
其中,FF表示变形梯度的转置,I表示单位矩阵,即未变形状态,因此Green-Lagrange应变张量为:
在未变形配置F=I处,对Green-Lagrange应变张量进行泰勒展开:
则柯西应变张量,即无限小应变张量为
S42、构建旋转应变,具体为:共旋弹性力将线性材料应力-应变的简单性与非线性特征相结合,确保旋转不变性,
其中,μ∈c即为旋转应变,∑是F的奇异值,F=U∑VT,tr2(∈c)表示矩阵∈c的迹的平方,tr2(∑-I)表示矩阵∑-I的迹的平方,
μ,λ表示拉梅常数,
k表示杨氏模量,与拉伸阻力有关;v表示泊松比,与不可压缩性有关,对于共旋弹性力的1阶Piola-Kirchhoff应力张量P为:
P(F)=2μ(F-R)+λtr(RTF-I)R;
S43、软组织模型离散化,具体为:
假设变形体不受内部任何摩擦力且变形体质量只分布在网格节点上,
在连续的软组织模型中,由能量守恒定律得:
对Etotal求偏导,并由第二牛顿定律得:
离散后,与Ni个节点相邻的独立元素Ωe,离散后的每一个节点总作用力为
其中,e∈Ni,Ee(x)表示离散后每个节点的能量;
S44、建立材料的本构关系,
得到形变矩阵Ds和参考形状矩阵Dm之间的关系为:
Ds=FDm
当W≠0,且Dm非奇异时,
进一步地,所述步骤S4中,所述建立各个有限元的节点的位移和载荷之间的表达式,从而得到单元刚度矩阵具体包括如下包括如下步骤:
S45、通过伽辽金法得到单元刚度矩阵,具体为:
取试函数
其中,Ni为取自完备函数集的线性无关基函数,ai为待定系数,即广义坐标,
由加权余量积分得,
由虚功原理,
δWint=δWext
可得单元刚度矩阵;
S46、求解单元Rayleigh阻尼矩阵,具体为:假设结构的阻尼矩阵是质量矩阵和刚度矩阵的组合,
[C]=a0[M]+a1[K]
其中,a0、a1是两个比例系数,
使用振型正交化条件,
Cn=a0Mn+a1Kn
引入模态阻尼比,
Cn=2ωnMnξn
得任意第n阶振型等效阻尼比得表达式为
假设结构各阶振型得阻尼比是相同的,即ξi=ξj=ξ,则待定系数
由Cn=a0Mn+a1Kn可得所需要的Rayleigh阻尼矩阵表达式。
进一步地,所述步骤S5具体包括如下步骤:
S51、将单元刚度矩阵、单元阻尼矩阵、单元质量矩阵组装成离散域的总刚度矩阵[K]、总阻尼矩阵[C]、总质量矩阵[M];
S52、采用Newmark全隐式积分,求解动力平衡微分方程,具体包括:
S522、选择积分步长Δt、参数β、γ,并计算积分常数
S523、形成有效刚度矩阵
S524、计算t+Δt时刻的有效荷载:
S525、求解t+Δt时刻的位移:
S527、基于步骤S526位移的求解,利用PCG方法,任选初始向量x0,计算r0=b-Ax0,解方程组Ms0=r0,令p0=s0,对于k=0,1,2,3…,计算:
其中,M为预处理矩阵;
S528、基于步骤S527的PCG方法中的预处理矩阵M,采用矩阵分裂法中的SSOR分裂求解,具体公式如下:
取
其中,ω使参数,0≤ω≤2,默认为1,
预处理矩阵M为:
S529、利用预处理共轭梯度算法,求解线性方程组,利用Richardson法则判断步骤S527中PCG方法是否收敛,若是,则迭代结束,得到下一时刻的位移,如否,则继续迭代,具体公式如下:
如果对于所有的x,
那么
|xk+1-x*|≤γ|xk-x*|
或者
limk→∞|xk+1-x*|=limk→∞γk|x0-x*|=0
达到线性收敛。
本发明通过把仿真区域离散成有限数量的六面体生成改进的八叉树网格,通过应力平衡方程的积分形式与材料的本构关系结合,建立各个有限元的节点的位移和载荷之间的表达式,从而得到单元刚度矩阵,通过联立方程组,在动力平衡情况下通过时间积分求解动力平衡微分方程即可得到节点位移随时间的变化,最后通过后处理渲染阶段显示软组织受力变形场景,逼真的模拟虚拟手术中任意形状软组织表皮受拉力变形的过程,并具有很高的实时性,大大减少了计算量,有效解决了以往有限元网格数量复杂,计算量大,不能实时仿真软组织变形过程的问题。
基于上述理由本发明可在虚拟手术仿真广泛推广。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案,下面将对实施例或现有技术描述中所需要使用的附图做以简单地介绍,显而易见地,下面描述中的附图是本发明的一些实施例,对于本领域普通技术人员来讲,在不付出创造性劳动性的前提下,还可以根据这些附图获得其他的附图。
图1为本发明主体流程流程图示意图。
图2为本发明具体流程流程图示意图。
图3为本发明基于AABB包围盒划分均匀六面体的示意图。
图4为本发明实施例基于八叉树网格的有限元模型的软组织变形方法的八叉树网格生成算法示意图,(a)为相邻的八个体元中的三角片示意图,(b)为体元合并以后,删除在六面体单元中共面的三角形面片,得到八叉树网格示意图。
图5为本发明实施例小立方体收缩成大立方体示意图,图(a)为收缩前,图(b)为收缩后。
具体实施方式
为使本发明实施例的目的、技术方案和优点更加清楚,下面将结合本发明实施例中的附图,对本发明实施例中的技术方案进行清楚、完整地描述,显然,所描述的实施例是本发明一部分实施例,而不是全部的实施例。基于本发明中的实施例,本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例,都属于本发明保护的范围。
如图1、图2所示,一种基于八叉树网格的有限元模型的软组织变形方法,包括如下步骤:
S1、绘制软组织的三维模型,基于轴对齐边界包围盒法,即AABB法,构建多个均匀的六面体,并过滤掉中心不在软组织几何模型空间的六面体;
S2、基于BON八叉树原理,创建BON八叉树,对各六面体编号,依次判断三维模型顶点是否在六面体内部,若是,则不属于边界元;若否,则表明三维模型顶点在六面体边界上,将该顶点插入对应的最底层六面体的八叉树结构中;
由AABB包围盒剖分可知,剖分的深度为N,以openGL的右手坐标系对每个子立方体进行编号,编号顺序采用先下后上,先左后右的顺序。编号方式采用0-7的二进制位000-111进行分配。因为S14中剔除了中心点在物体几何空间的外部的正六面体,所以有的节点并不是满8个子节点。
S3、判断三维模型顶点组成的三角网格是否共面,若是,则删除两个三角形的公共点,将原本的两个三角形转换为大三角形,将两个相邻六面体合并为中型六面体,提取合并后的中型六面体的边界点,直到所有相邻三角形不共面、不能进行收缩为止,进而得到最终的六面体网格;
如图4(a)所示,ΔAED、ΔBEF、ΔCDF与ΔEFD均共面,则删除点D、E、F,同时删除它们之间的公共面,得到如图4(b)所示的收缩后的六面体,如图5(a)所示,中层表示编号示意图中的大立方体,下层表示大立方体中的八个小立方体,分别从0-7依次储存,因为ΔAED、ΔBEF、ΔCDF与ΔEFD均共面,执行收缩,删除了0-7的小立方体,所以中层的父节点吸收了8个子节点。
S4、对最终的六面体网格进行有限元方法建模,求解软组织的变形过程,具体包括:通过应力平衡方程的积分形式与材料的本构关系相结合,建立各个有限元的节点的位移和载荷之间的表达式,从而得到单元刚度矩阵;
S5、将相邻单元节点上的单元刚度矩阵、单元阻尼矩阵、单元质量矩阵组装成离散域的总刚度矩阵、总阻尼矩阵、总质量矩阵,在动力平衡情况下对各矩阵的数值求解动力平衡微分方程,得到随时间变化的节点位移;
S6、通过对节点位移的渲染,不断更新软组织受力变形场景。
步骤S1中,通过如下方法绘制软组织的三维模型:
如图3所示,S11、利用Assimp模型加载库从模型文件中读取顶点坐标,法线坐标、纹理坐标、面片的相关索引信息以及模型材质库信息。它将整个模型加载进一个场景(Scene)对象,它会包含导入的模型/场景中的所有数据。Assimp会将场景载入为一系列的节点(Node),每个节点包含了场景对象中所储存数据的索引,每个节点都可以有任意数量的子节点;
S12、通过节点之间拓扑关系以及材质信息,利用OpenGL图形库绘制出软组织的三维模型。
S13、模型任意点的三维坐标分别为x、y、z,AABB包围盒内的最小点坐标为pmin=(xmin,ymin,zmin),最大点坐标为pmax=(xmax,ymax,zmax),AABB包围盒内点满足以下条件:
xmin≤x≤xmax
ymin≤y≤ymax
zmin≤z≤zmax,
中心点O=(xo,yo,zo)为两个顶点的中点,代表了包围盒的质点,其满足如下关系:
xo=(xmin+xmax)*0.5
yo=(ymin+ymax)*0.5
z0=(zmin+zmax)*0.5,
以点pmin、点pmax建立第一个最大的正六面体包围盒,设定剖分深度N,随后以每个边的中点为界,依次剖分N次,得到2的幂次方数量的六面体,通过改变N的值可以得到不同分辨率的模型;
S14、通过射线投影法,检查S13中生成的正六面体中心点是否在物体几何空间的外部,如果正六面体的中心点在物体几何空间的外部,那么把该中心点删除。
步骤S2中,所述基于八叉树的网格生成算法生成八叉树网格,具体包括以下步骤:
如图4所示,S21、基于BON八叉树原理,创建BON八叉树具体为:
以openGL的右手坐标系对每个六面体进行编号,编号顺序采用先下后上,先左后右的顺序,编号方式采用0-7的二进制位000-111进行分配。
进一步地,所述步骤S3中,通过如下方法判断三维模型顶点组成的三角网格是否共面:如果相邻六面体中的两个三角形至少有一个边相邻,且法向量相同,则两个三角形共面。
步骤S4中,所述通过应力平衡方程的积分形式与材料的本构关系相结合具体包括如下步骤:
S401、读取此时模型的位置信息;
S41、构建柯西应变张量,具体为:以X表示体中初始状态下的一点,在变形后,该点映射到了新的点x,体中每个点都会映射到新的点,这个映射表示为:
该点的位移为:
u=x-X
则变形梯度为:
变形前后长度的改变为:
||dx||2-||dX||2=dXT(FTF-I)dX
其中,FT表示变形梯度的转置,I表示单位矩阵,即未变形状态,因此Green-Lagrange应变张量为:
则柯西应变张量,即无限小应变张量为
S42、构建旋转应变,具体为:共旋弹性力将线性材料应力-应变的简单性与非线性特征相结合,确保旋转不变性,
其中,μ∈c即为旋转应变,∑是F的奇异值,F=U∑VT,tr2(∈c)表示矩阵∈c的迹的平方,tr2(∑-I)表示矩阵∑-I的迹的平方,
μ,λ表示拉梅常数,
k表示杨氏模量,与拉伸阻力有关;v表示泊松比,与不可压缩性有关,对于共旋弹性力的1阶Piola-Kirchhoff应力张量P为:
P(F)=2μ(F-R)+λtr(RTF-I)R;
S43、软组织模型离散化,具体为:
假设变形体不受内部任何摩擦力且变形体质量只分布在网格节点上,
在连续的软组织模型中,由能量守恒定律得:
对Etotal求偏导,并由第二牛顿定律得:
离散后,与Ni个节点相邻的独立元素Ωe,离散后的每一个节点总作用力为
其中,e∈Ni,Ee(x)表示离散后每个节点的能量;
S44、建立材料的本构关系,
得到形变矩阵Ds和参考形状矩阵Dm之间的关系为:
Ds=FDm
当W≠0,且Dm非奇异时,
所述步骤S4中,所述建立各个有限元的节点的位移和载荷之间的表达式,从而得到单元刚度矩阵具体包括如下包括如下步骤:
S45、通过伽辽金法得到单元刚度矩阵,具体为:
取试函数
其中,Ni为取自完备函数集的线性无关基函数,ai为待定系数,即广义坐标,
由加权余量积分得,
由虚功原理,
δWint=δWext
可得单元刚度矩阵;
S36、求解单元Rayleigh阻尼矩阵,具体为:假设结构的阻尼矩阵是质量矩阵和刚度矩阵的组合,
[C]=a0[M]+a1[K]
其中,a0、a1是两个比例系数,
使用振型正交化条件,
Cn=a0Mn+a1Kn
引入模态阻尼比,
Cn=2ωnMnξn
得任意第n阶振型等效阻尼比得表达式为
假设结构各阶振型得阻尼比是相同的,即ξi=ξj=ξ,则待定系数
由Cn=a0Mn+a1Kn可得所需要的Rayleigh阻尼矩阵表达式。
所述步骤S5具体包括如下步骤:
S51、将单元刚度矩阵、单元阻尼矩阵、单元质量矩阵组装成离散域的总刚度矩阵[K]、总阻尼矩阵[C]、总质量矩阵[M];
S52、采用Newmark全隐式积分,求解动力平衡微分方程,具体包括:
S522、选择积分步长Δt、参数β、γ,并计算积分常数
S523、形成有效刚度矩阵
S524、计算t+Δt时刻的有效荷载:
S525、求解t+Δt时刻的位移:
S527、基于步骤S526位移的求解,利用PCG方法,任选初始向量x0,计算r0=b-Ax0,解方程组Ms0=r0,令p0=s0,对于k=0,1,2,3…,计算:
其中,M为预处理矩阵;
S528、基于步骤S527的PCG方法中的预处理矩阵M,采用矩阵分裂法中的SSOR分裂求解,具体公式如下:
取
其中,ω使参数,0≤ω≤2,默认为1,
预处理矩阵M为:
S529、利用预处理共轭梯度算法,求解线性方程组,利用Richardson法则判断步骤S527中PCG方法是否收敛,若是,则迭代结束,得到下一时刻的位移,如否,则继续迭代,具体公式如下:
如果对于所有的x,
那么
|xk+1-x*|≤γ|xk-x*|
或者
limk→∞|xk+1-x*|=limk→∞γk|x0-x*|=0
达到线性收敛。
最后应说明的是:以上各实施例仅用以说明本发明的技术方案,而非对其限制;尽管参照前述各实施例对本发明进行了详细的说明,本领域的普通技术人员应当理解:其依然可以对前述各实施例所记载的技术方案进行修改,或者对其中部分或者全部技术特征进行等同替换;而这些修改或者替换,并不使相应技术方案的本质脱离本发明各实施例技术方案的范围。
Claims (3)
1.一种基于八叉树网格的有限元模型的软组织变形方法,其特征在于,包括如下步骤:
S1、绘制软组织的三维模型,基于轴对齐边界包围盒法,即AABB法,构建多个均匀的六面体,并过滤掉中心不在软组织几何模型空间的六面体;
S2、基于BON八叉树原理,创建BON八叉树,对各六面体编号,依次判断三维模型顶点是否在六面体内部,若是,则不属于边界元;若否,则表明三维模型顶点在六面体边界上,将该顶点插入对应的最底层六面体的八叉树结构中;
S3、判断三维模型顶点组成的三角网格是否共面,若是,则删除两个三角形的公共点,将原本的两个三角形转换为大三角形,将两个相邻六面体合并为中型六面体,提取合并后的中型六面体的边界点,直到所有相邻三角形不共面、不能进行收缩为止,进而得到最终的六面体网格;
S4、对最终的六面体网格进行有限元方法建模,求解软组织的变形过程,具体包括:通过应力平衡方程的积分形式与材料的本构关系相结合,建立各个有限元的节点的位移和载荷之间的表达式,从而得到单元刚度矩阵;
S5、将相邻单元节点上的单元刚度矩阵、单元阻尼矩阵、单元质量矩阵组装成离散域的总刚度矩阵、总阻尼矩阵、总质量矩阵,在动力平衡情况下对各矩阵的数值求解动力平衡微分方程,得到随时间变化的节点位移;
S6、通过对节点位移的渲染,不断更新软组织受力变形场景;
所述步骤S1中,通过如下方法绘制软组织的三维模型:
S11、通过模型加载库读取模型文件的顶点坐标、法线坐标、纹理坐标和面片的相关索引信息以及模型材质库信息;
S12、通过节点之间拓扑关系以及材质信息,利用OpenGL图形库绘制出软组织的三维模型;
所述步骤S1中,基于AABB法生成均匀的六面体网格,具体包括如下步骤:
S13、模型任意点的三维坐标分别为p=(a、b、c),AABB包围盒内的最小点坐标为pmin=(amin,bmin,cmin),最大点坐标为pmax=(amax,bmax,cmax),AABB包围盒内点满足以下条件:
amin≤a≤amax
bmin≤b≤bmax
cmin≤c≤cmax,
中心点O=(ao,bo,co)为两个顶点的中点,表示包围盒的质点,其满足如下关系:
ao=(amin+amax)*0.5
bo=(bmin+bmax)*0.5
co=(cmin+cmax)*0.5,
以点pmin、点pmax建立第一个最大的正六面体包围盒,设定剖分深度Mh,随后以每个边的中点为界,依次剖分Mh次,得到2的幂次方数量的六面体,通过改变Mh的值可以得到不同分辨率的模型;
S14、通过射线投影法,检查S13中生成的正六面体中心点是否在物体几何空间的外部,如果正六面体的中心点在物体几何空间的外部,那么把该中心点删除;
所述步骤S2中,所述基于BON八叉树原理,创建BON八叉树具体为:
以openGL的右手坐标系对每个六面体进行编号,编号顺序采用先下后上,先左后右的顺序,编号方式采用0-7的二进制位000-111进行分配;
所述步骤S3中,通过如下方法判断三维模型顶点组成的三角网格是否共面:如果相邻六面体中的两个三角形至少有一个边相邻,且法向量相同,则两个三角形共面;
步骤S4中,所述通过应力平衡方程的积分形式与材料的本构关系相结合具体包括如下步骤:
S401、读取此时模型的位置信息;
S41、构建柯西应变张量,具体为:以X表示六面体中初始状态下的一点的坐标,其中X的坐标为X=(X1,X2,X3,X4)在变形后,该点映射到了新的点x新,即此时该点的位置坐标为x新,其中x新的坐标为x新=(x1,x2,x3,x4)六面体中每个点都会映射到新的点,这个映射表示为:
该点的位移为:
u=x新-X
则变形梯度为:
变形前后长度的改变为:
||dx新||2-||dX||2=dXT(FTF-I)dX
其中,FT表示变形梯度的转置,I表示单位矩阵,即未变形状态,Green-Lagrange应变张量为:
在未变形配置F=I处,对Green-Lagrange应变张量进行泰勒展开:
则柯西应变张量,即无限小应变张量为
S42、构建旋转应变,具体为:共旋弹性力将线性材料应力-应变的简单性与非线性特征相结合,确保旋转不变性,
其中,μ∈c即为旋转应变,∑是F的奇异值,F=U∑VT,tr2(∈c)表示矩阵∈c的迹的平方,tr2(∑-I)表示矩阵∑-I的迹的平方,
μ,λ表示拉梅常数,
k表示杨氏模量,与拉伸阻力有关;v表示泊松比,与不可压缩性有关,对于共旋弹性力的1阶Piola-Kirchhoff应力张量P为:
P(F)=2μ(F-R0)+λtr(RTF-I)R0;
S43、软组织模型离散化,具体为:
假设变形体不受内部任何摩擦力且变形体质量只分布在网格节点上,
在连续的软组织模型中,遵循能量守恒定律
其中,Etotal表示软组织整体的总能量,E(x新)表示软组织整体的应变能,K(v新)表示软组织整体的动能,其中v新的坐标为v新=(v1,v2,v3,v4),N表示节点数量,mi表示每个节点的质量,表示每个节点的速度,
对Etotal求偏导,并由第二牛顿定律得:
离散后,与Ni个节点相邻的独立元素Ωe,离散后的每一个节点总作用力为:
其中,e∈Ni,Ee(x新)表示离散后每个节点的能量;
S44、建立材料的本构关系,
得到形变矩阵Ds和参考形状矩阵Dm之间的关系为:
Ds=FDm
当W≠0,且Dm非奇异时,
2.根据权利要求1所述的方法,其特征在于,所述步骤S4中,所述建立各个有限元的节点的位移和载荷之间的表达式,从而得到单元刚度矩阵具体包括如下包括如下步骤:
S45、通过伽辽金法得到单元刚度矩阵,具体为:
取试函数
其中,Ni为取自完备函数集的线性无关基函数,ai为待定系数,即广义坐标,
由加权余量积分得,
由虚功原理,得到单元刚度矩阵:
δWint=δWext;
S46、求解单元Rayleigh阻尼矩阵,具体为:假设结构的总阻尼矩阵是总质量矩阵和总刚度矩阵的组合,
[C]=a0[M]+a1[K]
其中,a0、a1是两个比例系数,
使用振型正交化条件,
Cn=a0Mn+a1Kn
引入模态阻尼比,
Cn=2ωnMnξn
得任意第n阶振型等效阻尼比得表达式为
假设结构各阶振型得阻尼比是相同的,即ξi=ξj=ξ,则待定系数
由Cn=a0Mn+a1Kn可得所需要的Rayleigh阻尼矩阵表达式。
3.根据权利要求2所述的方法,其特征在于,所述步骤S5具体包括如下步骤:
S51、将单元刚度矩阵、单元阻尼矩阵、单元质量矩阵组装成离散域的总刚度矩阵[K]、总阻尼矩阵[C]、总质量矩阵[M];
S52、采用Newmark全隐式积分,求解动力平衡微分方程,具体包括:
S522、选择积分步长Δt、参数β、γ,并计算积分常数
S523、形成有效刚度矩阵
S524、计算t+Δt时刻的有效荷载:
S525、求解t+Δt时刻的位移:
S527、基于步骤S526位移的求解,利用PCG方法,任选初始向量x0,计算r0=b-Ax0,解方程组Ms0=r0,令p0=s0,对于k=0,1,2,3…,计算:
其中,M为预处理矩阵;
S528、基于步骤S527的PCG方法中的预处理矩阵M,采用矩阵分裂法中的SSOR分裂求解,具体公式如下:
取
其中,ω是参数,0≤ω≤2,默认为1,
预处理矩阵M为:
S529、利用预处理共轭梯度算法,求解线性方程组,利用Richardson法则判断步骤S527中PCG方法是否收敛,若是,则迭代结束,得到下一时刻的位移,如否,则继续迭代,具体公式如下:
如果对于所有的x,
那么
|xk+1-x*|≤γb|xk-x*|
或者
达到线性收敛。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201810566190.1A CN108694290B (zh) | 2018-06-05 | 2018-06-05 | 一种基于八叉树网格的有限元模型的软组织变形方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201810566190.1A CN108694290B (zh) | 2018-06-05 | 2018-06-05 | 一种基于八叉树网格的有限元模型的软组织变形方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN108694290A CN108694290A (zh) | 2018-10-23 |
CN108694290B true CN108694290B (zh) | 2021-12-07 |
Family
ID=63848146
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201810566190.1A Expired - Fee Related CN108694290B (zh) | 2018-06-05 | 2018-06-05 | 一种基于八叉树网格的有限元模型的软组织变形方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN108694290B (zh) |
Families Citing this family (14)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN109408977B (zh) * | 2018-10-30 | 2022-07-22 | 河海大学 | 一种基于距离势函数可变形三维凸多面体块体离散单元法 |
CN109992856B (zh) * | 2019-03-19 | 2023-04-18 | 东北大学 | 一种适用于不同涂层厚度的硬涂层叶盘有限元建模方法 |
CN110729051B (zh) * | 2019-10-10 | 2022-11-22 | 中国科学院深圳先进技术研究院 | 一种介入手术中的导丝力学分析方法、系统及电子设备 |
CN110889893B (zh) * | 2019-10-25 | 2021-10-08 | 中国科学院计算技术研究所 | 表达几何细节和复杂拓扑的三维模型表示方法和系统 |
CN110826153B (zh) * | 2019-12-04 | 2022-07-29 | 中国直升机设计研究所 | 应用于直升机水中稳定性计算的水作用力模拟、实现方法 |
CN111814382B (zh) * | 2020-07-23 | 2023-09-22 | 中国工程物理研究院总体工程研究所 | 非平面波在多胞材料中传播的波阵面识别方法 |
CN112613211B (zh) * | 2020-12-22 | 2022-10-21 | 郑州大学 | 平面结构中任意三角形单元的变形分解方法 |
CN114969992A (zh) * | 2021-02-26 | 2022-08-30 | 湖南大学 | 基于节点积分的物体非线性大变形问题仿真方法 |
CN113689568B (zh) * | 2021-08-03 | 2023-05-23 | 南昌威爱信息科技有限公司 | 一种基于云端渲染的三维效果图高精度建模方法 |
CN113610875A (zh) * | 2021-09-03 | 2021-11-05 | 成都天地罔极科技有限公司 | 软组织填充效果的仿真方法及系统 |
CN113870434A (zh) * | 2021-09-23 | 2021-12-31 | 上海石指健康科技有限公司 | 一种基于有限元的生物组织仿真方法、装置和电子设备 |
CN114330034B (zh) * | 2022-03-09 | 2022-05-31 | 中国空气动力研究与发展中心计算空气动力研究所 | 一种预测可压-不可压复合材料弹性行为的计算方法 |
CN115017637B (zh) * | 2022-05-10 | 2023-03-28 | 西北工业大学 | 一种航天用张拉整体模块构件展开过程的动力学特性分析方法 |
CN116127611B (zh) * | 2023-04-13 | 2023-06-20 | 中国人民解放军国防科技大学 | 一种水下航行器动态仿真方法 |
Family Cites Families (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
WO2011055740A1 (ja) * | 2009-11-05 | 2011-05-12 | 国立大学法人名古屋工業大学 | 軟組織弾性率分布計測方法および軟組織弾性率分布計測装置 |
CN103699714B (zh) * | 2013-12-01 | 2016-08-31 | 北京航空航天大学 | 一种基于有限元和无网格耦合的柔性物体实时切割仿真方法 |
CN105302974B (zh) * | 2015-11-06 | 2018-06-08 | 北京航空航天大学 | 一种基于有限元和时变模态分析的柔性物体实时切割仿真方法 |
-
2018
- 2018-06-05 CN CN201810566190.1A patent/CN108694290B/zh not_active Expired - Fee Related
Also Published As
Publication number | Publication date |
---|---|
CN108694290A (zh) | 2018-10-23 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN108694290B (zh) | 一种基于八叉树网格的有限元模型的软组织变形方法 | |
Lloyd et al. | Identification of spring parameters for deformable object simulation | |
Teran et al. | Finite volume methods for the simulation of skeletal muscle | |
Kiendl et al. | Isogeometric shape optimization of shells using semi-analytical sensitivity analysis and sensitivity weighting | |
Hammer et al. | Mass-spring model for simulation of heart valve tissue mechanical behavior | |
Zou et al. | A new deformation model of biological tissue for surgery simulation | |
CN105654543B (zh) | 面向激光点云数据的阔叶树真实叶片建模与形变方法 | |
CN110289104B (zh) | 软组织按压和形变恢复的模拟方法 | |
Lu et al. | Cylindrical element: Isogeometric model of continuum rod | |
Qin et al. | A novel modeling framework for multilayered soft tissue deformation in virtual orthopedic surgery | |
CN111341449B (zh) | 一种虚拟血管介入手术训练的模拟方法 | |
CN113409443B (zh) | 一种基于位置约束和非线性弹簧的软组织建模方法 | |
CN108710735A (zh) | 一种实时交互的无网格软组织形变模拟方法 | |
Zhang et al. | An optimized mass-spring model with shape restoration ability based on volume conservation | |
CN103699716A (zh) | 一种个性化三维医学图像驱动的器官虚拟显示方法 | |
CN104794742B (zh) | 一种基于有限元方法的气球膨胀动画模拟方法 | |
Baudet et al. | Integrating tensile parameters in mass-spring system for deformable object simulation | |
Xu et al. | An improved realistic mass-spring model for surgery simulation | |
KR101350732B1 (ko) | 변형체의 실시간 시뮬레이션을 위한 다해상도 무요소법 | |
JP2007293540A (ja) | 連続弾性体変形シミュレーション方法、プログラム、および記録媒体 | |
Ge et al. | Blending isogeometric and Lagrangian elements in three-dimensional analysis | |
Zhong et al. | An autowave based methodology for deformable object simulation | |
Santhanam et al. | Physiologically-based modeling and visualization of deformable lungs | |
Yan et al. | Soft tissue deformation simulation in virtual surgery using nonlinear finite element method | |
Huang et al. | Virtual surgery deformable modelling employing gpu based computation |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
PB01 | Publication | ||
PB01 | Publication | ||
SE01 | Entry into force of request for substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant | ||
CF01 | Termination of patent right due to non-payment of annual fee | ||
CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20211207 |