CN108694290A - 一种基于八叉树网格的有限元模型的软组织变形方法 - Google Patents

Abstract

本发明提供一种基于八叉树网格的有限元模型的软组织变形方法。本发明方法，包括：绘制软组织的三维模型，基于AABB法构建多个均匀的六面体网格，在六面体网格基础上基于八叉树的网格生成算法生成八叉树网格，对六面体网格进行有限元方法建模，求解软组织的变形过程，将相邻单元节点上的单元刚度矩阵组装成离散域的总刚度矩阵，在动力平衡情况下通过时间积分对各矩阵的数值求解动力平衡微分方程，得到随时间变化的节点位移，通过对节点位移的渲染，显示软组织受力变形场景。本发明逼真的模拟虚拟手术中任意形状软组织表皮受拉力变形的过程，有很高的实时性，减少了计算量，解决了以往有限元网格数量复杂，不能实时仿真软组织变形过程的问题。

Description

一种基于八叉树网格的有限元模型的软组织变形方法

技术领域

本发明涉及虚拟手术仿真领域，尤其涉及一种基于八叉树网格的有限元模型的软组织变形方法。

背景技术

传统临床医学通常使用橡胶人体模型和人类尸体、小白鼠、青蛙等生物体作为手术训练对象。橡胶人体模型结构简单、功能单一，仿真度存在误差；人类尸体的材料活性与活体组织存在较大差异，复用性低，且随着我国社会的发展，尸体的获得越来越难；动物与人体的生理结构不同，作为临床手术训练对象并不准确，同时存在着迫害动物等诸多道德问题。

随着计算机硬件处理性能的提升，虚拟手术仿真系统得到了广泛研究。而软组织变形是虚拟手术中的重要组成部分。目前常用的软组织变形方法从物理计算模型角度可以分为三类：第一类，基于有限元方法的仿真学方法，由于包含软组织的材料特性，能够逼真模拟软组织的形变，使得其结果精度高，但有限元的形变计算是基于大量的网格单元，因此计算代价很高。第二类，基于位置的动力学方法，虽然计算速度快、不存在稳定性问题，但其产生的变形效果仅仅是貌似符合物理定律，而非真正符合物理定律。第三类，无网格方法，这种方法适用于大形变的情况，但由于采样点比较密集，计算负担比较重。

发明内容

根据上述提出的技术问题，而提供一种在不降低仿真逼真度的同时减少计算量的基于八叉树网格的有限元模型的软组织变形方法。本发明采用的技术手段如下：

一种基于八叉树网格的有限元模型的软组织变形方法,包括如下步骤：

S1、绘制软组织的三维模型，基于轴对齐边界包围盒法，即AABB法，构建多个均匀的六面体，并过滤掉中心不在软组织几何模型空间的六面体；

S2、基于BON八叉树原理，创建BON八叉树，对各六面体编号，依次判断三维模型顶点是否在六面体内部，若是，则不属于边界元；若否，则表明三维模型顶点在六面体边界上，将该顶点插入对应的最底层六面体的八叉树结构中；

S3、判断三维模型顶点组成的三角网格是否共面，若是，则删除两个三角形的公共点，将原本的两个三角形转换为大三角形，将两个相邻六面体合并为中型六面体，提取合并后的中型六面体的边界点，直到所有相邻三角形不共面、不能进行收缩为止，进而得到最终的六面体网格；

S4、对最终的六面体网格进行有限元方法建模，求解软组织的变形过程，具体包括：通过应力平衡方程的积分形式与材料的本构关系相结合，建立各个有限元的节点的位移和载荷之间的表达式，从而得到单元刚度矩阵；

S5、将相邻单元节点上的单元刚度矩阵、单元阻尼矩阵、单元质量矩阵组装成离散域的总刚度矩阵、总阻尼矩阵、总质量矩阵，在动力平衡情况下对各矩阵的数值求解动力平衡微分方程，得到随时间变化的节点位移；

S6、通过对节点位移的渲染，不断更新软组织受力变形场景。

进一步地，步骤S1中，通过如下方法绘制软组织的三维模型：

S11、通过模型加载库读取模型文件的顶点坐标、法线坐标、纹理坐标和面片的相关索引信息以及模型材质库信息；

S12、通过节点之间拓扑关系以及材质信息，利用OpenGL图形库绘制出软组织的三维模型。

进一步地，所述基于AABB法生成均匀的六面体网格，具体包括如下步骤：

S13、模型任意点的三维坐标分别为x、y、z，AABB包围盒内的最小点坐标为p_min＝(x_min，y_min，z_min)，最大点坐标为p_max＝(x_max，y_max，z_max)，AABB包围盒内点满足以下条件：

x_min≤x≤x_max

y_min≤y≤y_max

z_min≤z≤z_max，

中心点O＝(x_o，y_o，z_o)为两个顶点的中点，表示包围盒的质点，其满足如下关系：

x_o＝(x_min+x_max)*0.5

y_o＝(y_min+y_max)*0.5

z_o＝(z_min+z_max)*0.5，

以点p_min、点p_max建立第一个最大的正六面体包围盒，设定剖分深度N，随后以每个边的中点为界，依次剖分N次，得到2的幂次方数量的六面体，通过改变N的值可以得到不同分辨率的模型；

S14、通过射线投影法，检查S13中生成的正六面体中心点是否在物体几何空间的外部，如果正六面体的中心点在物体几何空间的外部，那么把该中心点删除。

进一步地，步骤S2中，所述基于BON八叉树原理，创建BON八叉树具体为：

以openGL的右手坐标系对每个六面体进行编号，编号顺序采用先下后上，先左后右的顺序，编号方式采用0-7的二进制位000-111进行分配。

进一步地，所述步骤S3中，通过如下方法判断三维模型顶点组成的三角网格是否共面：如果相邻六面体中的两个三角形至少有一个边相邻，且法向量相同，则两个三角形共面。

进一步地，步骤S4中，所述通过应力平衡方程的积分形式与材料的本构关系相结合具体包括如下步骤：

S401、读取此时模型的位置信息；

S41、构建柯西应变张量，具体为：以X表示体中初始状态下的一点，在变形后，该点映射到了新的点x，体中每个点都会映射到新的点，这个映射表示为：

该点的位移为：

u＝x-X

则变形梯度为：

变形前后长度的改变为：

||dx||²-||dX||²＝dX^T(F^TF-I)dX

其中，F^F表示变形梯度的转置，I表示单位矩阵，即未变形状态，因此Green-Lagrange应变张量为：

在未变形配置F＝I处，对Green-Lagrange应变张量进行泰勒展开：

则柯西应变张量，即无限小应变张量为

S42、构建旋转应变，具体为：共旋弹性力将线性材料应力-应变的简单性与非线性特征相结合，确保旋转不变性，

利用极坐标分解F＝RS，构造∈_c＝S-I替代中的ε，得到能量密度，即：

其中，μ∈_c即为旋转应变，∑是F的奇异值，F＝U∑V^T，tr²(∈_c)表示矩阵∈_c的迹的平方，tr²(∑-I)表示矩阵∑-I的迹的平方，

μ，λ表示拉梅常数，

k表示杨氏模量，与拉伸阻力有关；v表示泊松比，与不可压缩性有关，对于共旋弹性力的1阶Piola-Kirchhoff应力张量P为：

P(F)＝2μ(F-R)+λtr(R^TF-I)R；

S43、软组织模型离散化，具体为：

假设变形体不受内部任何摩擦力且变形体质量只分布在网格节点上，

在连续的软组织模型中，由能量守恒定律得：

其中，E_total表示软组织整体的总能量，E(x)表示软组织整体的应变能，K(v)表示软组织整体的动能。N表示节点数量，m_i表示每个节点的质量，表示每个节点的速度，

对E_total求偏导，并由第二牛顿定律得：

或

其中，表示单个网格节点相关的弹性力，f表示所有网格节点相关的弹性力合力，

离散后，与N_i个节点相邻的独立元素Ω_e，离散后的每一个节点总作用力为

其中，e∈N_i，E^e(x)表示离散后每个节点的能量；

S44、建立材料的本构关系，

设在x轴方向上，形变映射每一个顶点满足则有

得到形变矩阵D_s和参考形状矩阵D_m之间的关系为：

D_s＝FD_m

得到未发生形变时的四面体体积为：

当W≠0，且D_m非奇异时，

由知，四面体的四个顶点的弹力为

通过四面体的弹力计算推导出出六面体各顶点的弹力。

进一步地，所述步骤S4中，所述建立各个有限元的节点的位移和载荷之间的表达式，从而得到单元刚度矩阵具体包括如下包括如下步骤：

S45、通过伽辽金法得到单元刚度矩阵，具体为：

取试函数

其中，N_i为取自完备函数集的线性无关基函数，a_i为待定系数，即广义坐标，

设增量为R，取权函数W_i＝N_i，在边界上，

由加权余量积分得，

由虚功原理，

δW_int＝δW_ext

可得单元刚度矩阵；

S46、求解单元Rayleigh阻尼矩阵，具体为：假设结构的阻尼矩阵是质量矩阵和刚度矩阵的组合，

[C]＝a₀[M]+a₁[K]

其中，a₀、a₁是两个比例系数，

使用振型正交化条件，

C_n＝a₀M_n+a₁K_n

引入模态阻尼比，

C_n＝2ω_nM_nξ_n

得任意第n阶振型等效阻尼比得表达式为

假设结构各阶振型得阻尼比是相同的，即ξ_i＝ξ_j＝ξ，则待定系数

由C_n＝a₀M_n+a₁K_n可得所需要的Rayleigh阻尼矩阵表达式。

进一步地，所述步骤S5具体包括如下步骤：

S51、将单元刚度矩阵、单元阻尼矩阵、单元质量矩阵组装成离散域的总刚度矩阵[K]、总阻尼矩阵[C]、总质量矩阵[M]；

S52、采用Newmark全隐式积分，求解动力平衡微分方程，具体包括：

S521、给定初始值{u}₀、

S522、选择积分步长Δt、参数β、γ，并计算积分常数

α₆＝Δt(1-β)，α₇＝βΔt

S523、形成有效刚度矩阵

S524、计算t+Δt时刻的有效荷载：

S525、求解t+Δt时刻的位移：

S526、计算t+Δt时刻的速度和加速度

S527、基于步骤S526位移的求解，利用PCG方法，任选初始向量x₀，计算r₀＝b-Ax₀，解方程组Ms₀＝r₀，令p₀＝s₀，对于k＝0,1,2,3…，计算：

其中，M为预处理矩阵；

S528、基于步骤S527的PCG方法中的预处理矩阵M，采用矩阵分裂法中的SSOR分裂求解，具体公式如下：

将A分裂其中D＝diag(α₁₁，α₂₂，Λ，α_nn)，C_L为严格下三角矩阵，,由A相应元素取负号以后构成，

取

其中，ω使参数，0≤ω≤2，默认为1，

预处理矩阵M为：

S529、利用预处理共轭梯度算法，求解线性方程组，利用Richardson法则判断步骤S527中PCG方法是否收敛，若是，则迭代结束，得到下一时刻的位移，如否，则继续迭代，具体公式如下：

如果对于所有的x，

且|x₀-x^*|＜δ

那么

|x_k+1-x^*|≤γ|x_k-x^*|

或者

lim_k→∞|x_k+1-x^*|＝lim_k→∞γ^k|x₀-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包围盒内的最小点坐标为p_min＝(x_min，y_min，z_min)，最大点坐标为p_max＝(x_max，y_max，z_max)，AABB包围盒内点满足以下条件：

x_min≤x≤x_max

y_min≤y≤y_max

z_min≤z≤z_max，

中心点O＝(x_o，y_o，z_o)为两个顶点的中点，代表了包围盒的质点，其满足如下关系：

x_o＝(x_min+x_max)*0.5

y_o＝(y_min+y_max)*0.5

z₀＝(z_min+z_max)*0.5，

以点p_min、点p_max建立第一个最大的正六面体包围盒，设定剖分深度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||²-||dX||²＝dX^T(F^TF-I)dX

其中，F^T表示变形梯度的转置，I表示单位矩阵，即未变形状态，因此Green-Lagrange应变张量为：

但是是变形的二次非线性函数，增加了构造本构模型的复杂性，并且导致节点力离散化。因此在未变形配置F＝I处，对Green-Lagrange应变张量进行泰勒展开：

则柯西应变张量，即无限小应变张量为

S42、构建旋转应变，具体为：共旋弹性力将线性材料应力-应变的简单性与非线性特征相结合，确保旋转不变性，

利用极坐标分解F＝RS，构造∈_c＝S-I替代中的ε，得到能量密度，即：

其中，μ∈_c即为旋转应变，∑是F的奇异值，F＝U∑V^T，tr²(∈_c)表示矩阵∈_c的迹的平方，tr²(∑-I)表示矩阵∑-I的迹的平方，

μ，λ表示拉梅常数，

k表示杨氏模量，与拉伸阻力有关；v表示泊松比，与不可压缩性有关，对于共旋弹性力的1阶Piola-Kirchhoff应力张量P为：

P(F)＝2μ(F-R)+λtr(R^TF-I)R；

S43、软组织模型离散化，具体为：

假设变形体不受内部任何摩擦力且变形体质量只分布在网格节点上，

在连续的软组织模型中，由能量守恒定律得：

其中，E_total表示软组织整体的总能量，E(x)表示软组织整体的应变能，K(v)表示软组织整体的动能。N表示节点数量，m_i表示每个节点的质量，表示每个节点的速度，

对E_total求偏导，并由第二牛顿定律得：

或

其中，表示单个网格节点相关的弹性力，f表示所有网格节点相关的弹性力合力，

离散后，与N_i个节点相邻的独立元素Ω_e，离散后的每一个节点总作用力为

其中，e∈N_i，E^e(x)表示离散后每个节点的能量；

S44、建立材料的本构关系，

设在x轴方向上，形变映射每一个顶点满足则有

得到形变矩阵D_s和参考形状矩阵D_m之间的关系为：

D_s＝FD_m

得到未发生形变时的四面体体积为：

当W≠0，且D_m非奇异时，

由知，四面体的四个顶点的弹力为

通过四面体的弹力计算推导出出六面体各顶点的弹力。

所述步骤S4中，所述建立各个有限元的节点的位移和载荷之间的表达式，从而得到单元刚度矩阵具体包括如下包括如下步骤：

S45、通过伽辽金法得到单元刚度矩阵，具体为：

取试函数

其中，N_i为取自完备函数集的线性无关基函数，a_i为待定系数，即广义坐标，

设增量为R，取权函数W_i＝N_i，在边界上，

由加权余量积分得，

由虚功原理，

δW_int＝δW_ext

可得单元刚度矩阵；

S36、求解单元Rayleigh阻尼矩阵，具体为：假设结构的阻尼矩阵是质量矩阵和刚度矩阵的组合，

[C]＝a₀[M]+a₁[K]

其中，a₀、a₁是两个比例系数，

使用振型正交化条件，

C_n＝a₀M_n+a₁K_n

引入模态阻尼比，

C_n＝2ω_nM_nξ_n

得任意第n阶振型等效阻尼比得表达式为

假设结构各阶振型得阻尼比是相同的，即ξ_i＝ξ_j＝ξ，则待定系数

由C_n＝a₀M_n+a₁K_n可得所需要的Rayleigh阻尼矩阵表达式。

所述步骤S5具体包括如下步骤：

S51、将单元刚度矩阵、单元阻尼矩阵、单元质量矩阵组装成离散域的总刚度矩阵[K]、总阻尼矩阵[C]、总质量矩阵[M]；

S52、采用Newmark全隐式积分，求解动力平衡微分方程，具体包括：

S521、给定初始值{u}₀、

S522、选择积分步长Δt、参数β、γ，并计算积分常数

α₆＝Δt(1-β)，α₇＝βΔt

S523、形成有效刚度矩阵

S524、计算t+Δt时刻的有效荷载：

S525、求解t+Δt时刻的位移：

S526、计算t+Δt时刻的速度和加速度

S527、基于步骤S526位移的求解，利用PCG方法，任选初始向量x₀，计算r₀＝b-Ax₀，解方程组Ms₀＝r₀，令p₀＝s₀，对于k＝0,1,2,3…，计算：

其中，M为预处理矩阵；

S528、基于步骤S527的PCG方法中的预处理矩阵M，采用矩阵分裂法中的SSOR分裂求解，具体公式如下：

将A分裂其中D＝diag(α₁₁，α₂₂，Λ，α_nn)，C_L为严格下三角矩阵，,由A相应元素取负号以后构成，

取

其中，ω使参数，0≤ω≤2，默认为1，

预处理矩阵M为：

S529、利用预处理共轭梯度算法，求解线性方程组，利用Richardson法则判断步骤S527中PCG方法是否收敛，若是，则迭代结束，得到下一时刻的位移，如否，则继续迭代，具体公式如下：

如果对于所有的x，

且|x₀-x^*|＜δ

那么

|x_k+1-x^*|≤γ|x_k-x^*|

或者

lim_k→∞|x_k+1-x^*|＝lim_k→∞γ^k|x₀-x^*|＝0

达到线性收敛。

最后应说明的是：以上各实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述各实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分或者全部技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的范围。

Claims (8)

1.一种基于八叉树网格的有限元模型的软组织变形方法，其特征在于，包括如下步骤：

S1、绘制软组织的三维模型，基于轴对齐边界包围盒法，即AABB法，构建多个均匀的六面体，并过滤掉中心不在软组织几何模型空间的六面体；

S2、基于BON八叉树原理，创建BON八叉树，对各六面体编号，依次判断三维模型顶点是否在六面体内部，若是，则不属于边界元；若否，则表明三维模型顶点在六面体边界上，将该顶点插入对应的最底层六面体的八叉树结构中；

S3、判断三维模型顶点组成的三角网格是否共面，若是，则删除两个三角形的公共点，将原本的两个三角形转换为大三角形，将两个相邻六面体合并为中型六面体，提取合并后的中型六面体的边界点，直到所有相邻三角形不共面、不能进行收缩为止，进而得到最终的六面体网格；

S4、对最终的六面体网格进行有限元方法建模，求解软组织的变形过程，具体包括：通过应力平衡方程的积分形式与材料的本构关系相结合，建立各个有限元的节点的位移和载荷之间的表达式，从而得到单元刚度矩阵；

S5、将相邻单元节点上的单元刚度矩阵、单元阻尼矩阵、单元质量矩阵组装成离散域的总刚度矩阵、总阻尼矩阵、总质量矩阵，在动力平衡情况下对各矩阵的数值求解动力平衡微分方程，得到随时间变化的节点位移；

S6、通过对节点位移的渲染，不断更新软组织受力变形场景。

2.根据权利要求1所述的方法，其特征在于，所述步骤S1中，通过如下方法绘制软组织的三维模型：

S11、通过模型加载库读取模型文件的顶点坐标、法线坐标、纹理坐标和面片的相关索引信息以及模型材质库信息；

S12、通过节点之间拓扑关系以及材质信息，利用OpenGL图形库绘制出软组织的三维模型。

3.根据权利要求2所述的方法，所述步骤S1中，基于AABB法生成均匀的六面体网格，具体包括如下步骤：

S13、模型任意点的三维坐标分别为x、y、z，AABB包围盒内的最小点坐标为p_min＝(x_min，y_min，z_min)，最大点坐标为p_max＝(x_max，y_max，z_max)，AABB包围盒内点满足以下条件：

x_min≤x≤x_max

y_min≤y≤y_max

z_min≤z≤z_max，

中心点O＝(x_o，y_o，z_o)为两个顶点的中点，表示包围盒的质点，其满足如下关系：

x_o＝(x_min+x_max)*0.5

y_o＝(y_min+y_max)*0.5

z_o＝(z_min+z_max)*0.5，

以点p_min、点p_max建立第一个最大的正六面体包围盒，设定剖分深度N，随后以每个边的中点为界，依次剖分N次，得到2的幂次方数量的六面体，通过改变N的值可以得到不同分辨率的模型；

S14、通过射线投影法，检查S13中生成的正六面体中心点是否在物体几何空间的外部，如果正六面体的中心点在物体几何空间的外部，那么把该中心点删除。

4.根据权利要求3所述的方法，其特征在于，所述步骤S2中，所述基于BON八叉树原理，创建BON八叉树具体为：

以openGL的右手坐标系对每个六面体进行编号，编号顺序采用先下后上，先左后右的顺序，编号方式采用0-7的二进制位000-111进行分配。

5.根据权利要求4所述的方法，其特征在于，所述步骤S3中，通过如下方法判断三维模型顶点组成的三角网格是否共面：如果相邻六面体中的两个三角形至少有一个边相邻，且法向量相同，则两个三角形共面。

6.根据权利要求5所述的方法，其特征在于，步骤S4中，所述通过应力平衡方程的积分形式与材料的本构关系相结合具体包括如下步骤：

S401、读取此时模型的位置信息；

S41、构建柯西应变张量，具体为：以X表示体中初始状态下的一点，在变形后，该点映射到了新的点x，体中每个点都会映射到新的点，这个映射表示为：

该点的位移为：

u＝x-X

则变形梯度为：

变形前后长度的改变为：

||dx||²-||dX||²＝dX^T(F^TF-I)dX

其中，F^T表示变形梯度的转置，I表示单位矩阵，即未变形状态，Green-Lagrange应变张量为：

在未变形配置F＝I处，对Green-Lagrange应变张量进行泰勒展开：

则柯西应变张量，即无限小应变张量为

S42、构建旋转应变，具体为：共旋弹性力将线性材料应力-应变的简单性与非线性特征相结合，确保旋转不变性，

利用极坐标分解F＝RS，构造∈_c＝S-I替代中的ε，得到能量密度，即：

其中，μ∈_c即为旋转应变，∑是F的奇异值，F＝U∑V^T，tr²(∈_c)表示矩阵∈_c的迹的平方，tr²(∑-I)表示矩阵∑-I的迹的平方，

μ，λ表示拉梅常数，

k表示杨氏模量，与拉伸阻力有关；v表示泊松比，与不可压缩性有关，对于共旋弹性力的1阶Piola-Kirchhoff应力张量P为：

P(F)＝2μ(F-R)+λtr(R^TF-I)R；

S43、软组织模型离散化，具体为：

假设变形体不受内部任何摩擦力且变形体质量只分布在网格节点上，

在连续的软组织模型中，遵循能量守恒定律

其中，E_total表示软组织整体的总能量，E(x)表示软组织整体的应变能，K(v)表示软组织整体的动能，N表示节点数量，m_i表示每个节点的质量，表示每个节点的速度，

对E_total求偏导，并由第二牛顿定律得：

或

其中，表示单个网格节点相关的弹性力，f表示所有网格节点相关的弹性力合力，

离散后，与N_i个节点相邻的独立元素Ω_e，离散后的每一个节点总作用力为：

其中，e∈N_i，E^e(x)表示离散后每个节点的能量；

S44、建立材料的本构关系，

设在x轴方向上，形变映射每一个顶点满足则有

得到形变矩阵D_s和参考形状矩阵D_m之间的关系为：

D_s＝FD_m

得到未发生形变时的四面体体积为：

当W≠0，且D_m非奇异时，

由知，四面体的四个顶点的弹力为

通过四面体的弹力计算推导出六面体各顶点的弹力。

7.根据权利要求6所述的方法，其特征在于，所述步骤S4中，所述建立各个有限元的节点的位移和载荷之间的表达式，从而得到单元刚度矩阵具体包括如下包括如下步骤：

S45、通过伽辽金法得到单元刚度矩阵，具体为：

取试函数

其中，N_i为取自完备函数集的线性无关基函数，a_i为待定系数，即广义坐标，

设增量为R，取权函数W_i＝N_i，在边界上，

由加权余量积分得，

由虚功原理，得到单元刚度矩阵：

δW_int＝δW_ext；

S46、求解单元Rayleigh阻尼矩阵，具体为：假设结构的阻尼矩阵是质量矩阵和刚度矩阵的组合，

[C]＝a₀[M]+a₁[K]

其中，a₀、a₁是两个比例系数，

使用振型正交化条件，

C_n＝a₀M_n+a₁K_n

引入模态阻尼比，

C_n＝2ω_nM_nξ_n

得任意第n阶振型等效阻尼比得表达式为

假设结构各阶振型得阻尼比是相同的，即ξ_i＝ξ_j＝ξ，则待定系数

由C_n＝a₀M_n+a₁K_n可得所需要的Rayleigh阻尼矩阵表达式。

8.根据权利要求7所述的方法，其特征在于，所述步骤S5具体包括如下步骤：

S51、将单元刚度矩阵、单元阻尼矩阵、单元质量矩阵组装成离散域的总刚度矩阵[K]、总阻尼矩阵[C]、总质量矩阵[M]；

S52、采用Newmark全隐式积分，求解动力平衡微分方程，具体包括：

S521、给定初始值{u}₀、

S522、选择积分步长Δt、参数β、γ，并计算积分常数

α₆＝Δt(1-β)，α₇＝βΔt

S523、形成有效刚度矩阵

S524、计算t+Δt时刻的有效荷载：

S525、求解t+Δt时刻的位移：

S526、计算t+Δt时刻的速度和加速度

S527、基于步骤S426位移的求解，利用PCG方法，任选初始向量x₀，计算r₀＝b-Ax₀，解方程组Ms₀＝r₀，令p₀＝s₀，对于k＝0,1,2,3…，计算：

其中，M为预处理矩阵；

S528、基于步骤S527的PCG方法中的预处理矩阵M，采用矩阵分裂法中的SSOR分裂求解，具体公式如下：

将A分裂其中D＝diag(α₁₁，α₂₂，Λ，α_nn)，C_L为严格下三角矩阵，,由A相应元素取负号以后构成，

取

其中，ω使参数，0≤ω≤2，默认为1，

预处理矩阵M为：

S529、利用预处理共轭梯度算法，求解线性方程组，利用Richardson法则判断步骤S527中PCG方法是否收敛，若是，则迭代结束，得到下一时刻的位移，如否，则继续迭代，具体公式如下：

如果对于所有的x，

且|x₀-x^*|＜δ

那么

|x_k+1-x^*|≤γ|x_k-x^*|

或者

lim_k→∞|x_k+1-x^*|＝lim_k→∞γ^k|x₀-x^*|＝0

达到线性收敛。