CN105302974B - 一种基于有限元和时变模态分析的柔性物体实时切割仿真方法 - Google Patents
一种基于有限元和时变模态分析的柔性物体实时切割仿真方法 Download PDFInfo
- Publication number
- CN105302974B CN105302974B CN201510751613.3A CN201510751613A CN105302974B CN 105302974 B CN105302974 B CN 105302974B CN 201510751613 A CN201510751613 A CN 201510751613A CN 105302974 B CN105302974 B CN 105302974B
- Authority
- CN
- China
- Prior art keywords
- simulation
- finite element
- time
- subspace
- cutting
- 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.)
- Active
Links
Landscapes
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明涉及一种基于有限元和时变模态分析的柔性物体实时切割仿真方法，主要针对虚拟手术实时交互式软组织形变与切割仿真算法进行研究并实现，该方法包括以下步骤：对模型的仿真域进行体素化和有限元物理建模、基于先验模态分析的子空间构建、切割操作响应、子空间基向量的显著性分析和筛选、子空间时变有效性估计。本方法对仿真算法进行并行化设计，借助GPU强大的计算能力，达到实时交互式效率。
Description
技术领域
本发明属于物理仿真技术领域，具体涉及基于有限元和时变模态分析的柔性物体实时切割仿真方法。
背景技术
在计算机图形学领域，形变模型最早由Terzopolous等人引入。随着计算机图形技术的不断发展，特别是在最近的十年，由于生物力学仿真的外科手术规划与辅助决策在现代医学中的作用日益突出，高保真和高效率的物理变形/切割模拟已经受到了越来越多的关注。
针对弹性物体的形变与任意切割操作所涉及的物理以及弹性力学建模，有限元方法已经成为一个强有力的工具。Dick等人提出了一种利用在GPU上求解的几何多重网格计算框架来求解使用六面体有限元建模的线性弹性物体。
切割形变物体的四面体划分方法有Bielser等人引入。为了减少病态元的数量，Nienhuys和van der Stappen提出了沿元素表面切割。通过将划分限制为几种细分模式，切割所造成的额外模拟元素的数量可以减少。此方法的多分辨率解法由Ganovelli等人提出。Dick等人提出了自适应正六面体的近似方法被用来建模切割操作。
模态化简(Model Reduction)方法，也称为维度模态化简或者维度化简，是一种对使用微分方程进行运动仿真系统的化简技术。对于几何非线性实时形变物体(使用二次格林-拉格朗日应变张量)，Barbic和James提出了一种使用模型推导技术自动选择一组低维基向量方法。An等人提出了一种支持任意非线性材质的模态化简模型。
对于离线有限元形变仿真应用，迫切地需要一种对于初始设置没有依赖的模态化简方法来解决实际应用问题。Ryckelynck等人提出了一种基于先验知识的模态化简方法(APHR)用来求解非线性有限元仿真问题。Kim和James通过将在线模态化简技术根据系统状态自适应地替换离线仿真，降低了计算开销。
现有技术中，通常采用模态化简(Model Reduction)方法来解决大规模网格模型仿真问题，此类方法需要经过复杂的预处理(Preprocess)操作来得到满足仿真精度的模态化简模型，而在仿真过程中模态化简模型不会发生变化。显然，对于仿真过程中由于交互操作所产生的大形变以及切割操作所产生的网格拓扑结构发生的变化，上述方法受限于其预处理(Preprocess)操作从而不能做出准确及时的响应。
发明内容
本发明要解决的技术问题是：克服现有方法对大规模网格实时计算效能低，不能对仿真过程中所发生的网格更新操作进行动态自适应调整和更新的不足，提供一种基于有限元和时变模态分析的柔性物体实时切割仿真方法，并利用GPU硬件的计算能力，提高了计算与绘制效率，该方法主要利用有限元方法对仿真域进行建模，同时对有限元仿真输出进行基于先验模态分析的子空间动态构建。在形变和切割的物理仿真过程中，通过将两种仿真方法进行交替达到了对大尺度仿真模型的化简，并且通过对算法的并行优化使得物理仿真达到交互式效率。
本发明解决上述技术问题的技术方案为：一种基于有限元和时变模态分析的柔性物体实时切割仿真方法，包括如下步骤：
步骤(1)、将仿真对象的材质域依据仿真精度要求进行层次细节体素化并生成基于八叉树六面体网格；
步骤(2)、根据六面体网格进行有限元建模；
步骤(3)、在仿真过程中，对全物理域的有限元仿真输出进行基于先验模态分析的子空间动态构建；
步骤(4)、当子空间的仿真结果与有限元仿真结果的误差低于阀值时，使用子空间仿真方法替换有限元仿真；
步骤(5)、在仿真过程中，计算每个仿真顶点的局部旋转，对线性系统仿真误差进行累计修正；
步骤(6)、在仿真过程中，根据子空间时变有效性对子空间的基向量进行显著性分析和更新；
步骤(7)、当切割操作发生时，对切割所涉及到的有限元区域进行网格细分并通过克隆网格体现拓扑变化，并在新生成的有限元区域上重新进行子空间动态构建。
进一步的，步骤(1)包括：读取OBJ模型数据，根据材质边界曲面函数对模型建立域边界划分投票函数(voting function)；建立模型整体的AABB包围盒作为基于八叉树的六面体网格根节点，如果当前六面体空间内包含模型元素：顶点、三角面片和材质边界曲面，且没有达到最高层次精度，则对该六面体进行剖分并对新生成的8个子六面体重复上述操作；当所有六面体网格均符合上述条件时，停止操作，并对八叉树的叶子节点上的六面体元素使用其中心点作为输入参数调用边界划分投票函数确定其是否所属的仿真域；该步骤(1)最终生成了一个反应模型几何和材质特征的层次化八叉树六面体网格。
进一步的，步骤(2)包括：为了提高计算效率，采用线性基，并生成六面体有限元形函数以及形函数的三维空间导数，并以此构建形函数矩阵和应变矩阵；根据当前域所属的材质属性构建材质矩阵，所述的材质属性为杨氏模量和泊松比，依据有限元公式组织局部刚度矩阵、质量矩阵和施力向量。
进一步的，步骤(3)包括：利用先验模态分析方法不需要在初始时刻设置形变模型的基向量的优点，针对大尺度形变和切割操作所引发的形变模型的变化，模型通过对新产生的仿真结果进行合成，更新形变模型的基向量，使其能够准确地体现参照模型所发生的变化。
进一步的，步骤(4)包括：将子空间仿真结果投影到参照空间，并对子空间的仿真结果进行误差分析；当子空间的仿真结果与有限元仿真结果的误差低于阀值时，使用子空间仿真方法替换有限元仿真，从而加速形变模型的仿真效率。
进一步的，步骤(5)包括：针对所采用的线性弹性有限元仿真模型(线性柯西应变张量，Cauchy Strain Tensor)只适用于小范围形变的不足，利用局部应变旋转技术，将仿真元素形变过程中所包含的位移与旋转分离，支持大范围形变情况。
进一步的，步骤(6)包括：在有限元仿真和子空间仿真交替计算的过程中，子空间仿真会产生累计误差，特别是在动态仿真过程中，累计误差会被记录在系统的速度(Velocity)和加速度(Acceleration)中，通过建立子空间仿真的力残差和最近一次有限元仿真输出的对应关系模型，可以估计当前子空间仿真结果误差的变化趋势。
进一步的，步骤(7)包括：将切割工具建模为有三角面片组成的曲面，当切割操作发生时，通过六面体网格与三角面片做碰撞检测，能够得到被切割的网格，对于粗层次的网格首先要对其进行网格细分以达到仿真精度，对于被切割的六面体进行克隆操作以满足切割所造成的拓扑变化，最后对受影响的区域进行步骤(2)和步骤(3)，更新仿真方程。
本发明的原理在于：
基于有限元和时变模态分析的柔性物体实时切割仿真方法，包括如下步骤：
1)仿真模型体素化和物理建模
步骤(1)、读取OBJ模型数据，根据材质边界曲面函数对模型建立域边界划分投票函数(voting function)。对于一个网格元素而言，当它只满足一个投票函数时，它被认为是一个内部元素(Interior Element)，当它满足多于一个投票函数时，它被认为是边界元素(Boundary Element),否则被认为是外部元素(Exterior Element)不参与运算。建立模型整体的AABB包围盒作为基于八叉树的六面体网格根节点，如果当前六面体空间内包含模型元素(顶点、三角面片、材质边界曲面)且没有达到最高层次精度，则对该六面体进行剖分并对新生成的8个子六面体重复上述操作，在剖分的过程中，为了防止出现悬挂点(HangingNode),需要保证每一个八叉树叶子节点与其邻接节点的层次精度至多相差一个级别，否则也要进行剖分。当所有六面体网格均符合上述条件时，停止操作，并对八叉树的叶子节点上的六面体元素使用其中心点作为输入参数调用投票函数(voting function)确定其是否属于仿真域。该步骤最终生成了一个反应模型表面特征的层次化八叉树六面体网格。
步骤(2)、为了提高计算效率，本发明采用线性基(p＝[1,x,y,z])，并生成六面体有限元形函数以及形函数的三维空间导数，并以此构建形函数(Shape Function)矩阵和应变张量(Strain Tensor)。根据仿真域所属的材质属性(杨氏模量和泊松比)构建材质矩阵。依据有限元公式组织局部刚度矩阵、质量矩阵和施力向量。由于步骤(1)所生成的六面体网格形状是相同的只是尺寸不同，本发明只需要计算上述矩阵一次，通过乘以不同的放缩因子得到不同尺寸的对应矩阵。
2)基于先验模态分析的子空间动态构建
步骤(3)、利用先验模态分析方法不需要在初始时刻设置形变模型的基向量的优点，针对大尺度形变和切割操作所引发的形变模型的变化，模型通过对新产生的仿真结果进行合成，更新形变模型的基向量，并依据仿真结果的误差，动态获得更新基向量的规模，使其能够准确地体现参照模型所发生的变化。
步骤(4)、将子空间仿真结果投影到参照空间，并对子空间的仿真结果进行误差分析(Error Analysis)；当子空间的仿真结果与有限元仿真结果的误差低于阀值(Threshold)时，使用子空间仿真方法替换有限元仿真，从而加速形变模型的仿真效率。对于步骤(2)所得到的局部矩阵根据全局自由度发布到全局(Global Assembling)，由于有限元方法具有紧支特性，得到的全局矩阵具有稀疏特性，本发明开发了符合GPU存储特点的大型稀疏矩阵的自由存储数据结构(见图5)，利用预条件共轭梯度方法(PreconditionedConjugate Gradient Method)求解。对于步骤(3)所得到的化简形变子空间(ReducedDeformation Subspace)可以得到一个规模较小的对称稠密线性系统(Dense SymmetricLinear System),使用QR分解技术直接求解形变子空间的模态振幅，并投影变换为参照空间位移。为了进行动态仿真，采用NewMark隐式积分方法，利用所得到的位移场(Displacement Field)计算采样点的速度场(Velocity Field)和加速度场(AccelerationField)。
步骤(5)、针对所采用的线性弹性有限元仿真模型(线性柯西应变张量，CauchyStrain Tensor)只适用于小范围形变(即形变幅度不大的形变，如人体呼吸，相对于大范围形变，如人体跳跃)的不足，利用局部应变旋转(Modal Warping)技术，将仿真元素形变过程中所包含的位移与旋转分离，支持大范围形变情况。
步骤(6)、在有限元仿真和子空间仿真交替计算的过程中，子空间仿真会产生累计误差，特别是在动态仿真过程中，累计误差会被记录在系统的速度(Velocity)和加速度(Acceleration)中，通过建立子空间仿真的力残差和最近一次有限元仿真输出的对应关系模型，可以估计当前子空间仿真结果误差的变化趋势。
3)切割操作响应
步骤(7)、本发明将切割工具建模为有三角面片组成的曲面，当切割操作发生时，通过六面体网格与三角面片做碰撞检测，可以得到被切割的网格，对于粗层次的网格首先要对其进行网格细分以达到仿真精度，对于被切割的六面体进行克隆操作以满足切割所造成的拓扑变化，最后对受影响的区域进行步骤(2)和步骤(3)，更新全局方程组。对于切割造成的模型表面的几何变化采用Delaunay三角化进行更新。
本发明与现有技术相比的有点在于：
本发明能够对大尺度复杂几何属性柔性物体的变形和切割进行实时物理仿真。本发明主要有两点贡献：第一，提出了一种不受预计算限制的基于先验模态分析的子空间构建方法，该方法可以在大尺度形变和切割操作所带来的模态几何和拓扑结构发生变化的情况下，有效地对形变模型进行快速分析并产生相应地修正基向量。第二，建立一种子空间仿真与参照空间仿真之间力残差与最近一次有限元仿真输出的对应关系分析模型，从而能够对随后若干仿真步的误差范围进行预测，避免子空间仿真结果与参照模型仿真结果的误差过大，从而产生子空间重构等问题。
附图说明
图1为本发明整体流程图；
图2为模型初始化流程图；其中，(a)初始模型；(b)六面体网格模型；(c)基于八叉树体素化网格2D示意图；
图3为先验模态分析方法执行流程示意图；
图4为通过先验模态分析技术得到的仿真输出序列的振动频率，其中，(a),(b)(c)，(d)分别是通过先验模态分析技术对不同仿真输出(图中三维模型)得到的振动频率分布结果；
图5为无结构的大型稀疏矩阵内存布局示意图；
图6为模型切割响应示意图；其中，(a)初始八叉树网格模型；(b)切割发生后的六面体网格模型；(c)三角表面网格切割效果图；
图7为仿真过程执行流水线示意图；
图8为基于有限元和时变模态分析切割仿真输出：其中(a)为系统的初始模型；(b)为系统对图8中(a)的形变仿真结果；(c)为切割操作发生后系统对图8中(a)所做的涉及几何、物理的更新操作；(d)为系统使用形变子空间对图8中(c)进行仿真的结果。
具体实施方式
图1给出了基于有限元和时变模态分析的柔性物体实时切割仿真过程的总体处理流程，下面结合附图及具体实施方式进一步说明本发明。
本发明提供一种基于有限元和时变模态分析的柔性物体实时切割仿真方法，主要步骤介绍如下：
步骤(1)根据输入的OBJ模型(见图2中的(a))和模型边界方程建立模型边界投票函数(voting function，v(x,y,z))。其中投票函数以空间点坐标为输入参数，当空间坐标包含在所表示的模型仿真域内返回true，否则返回false。建立模型整体的AABB包围盒作为基于八叉树(Octree)的六面体网格根节点，如果当前六面体空间内包含模型元素(顶点、三角面片、材质边界曲面)且没有达到最高层次精度，则对该六面体进行剖分并对新生成的8个子六面体重复上述操作，在剖分的过程中，为了防止出现悬挂点(Hanging Node),需要保证每一个八叉树叶子节点与其邻接节点的层次精度至多相差一个级别，否则也要进行剖分。当所有六面体网格均符合上述条件时，停止操作，并对八叉树的叶子节点上的六面体元素使用其中心点作为输入参数调用前述的投票函数v(x,y,z)确定是否作为有效的仿真元素。该步骤最终生成了一个反应模型表面的层次化八叉树六面体网格(见图2中的(b))。图2中的(c)表示了层次化八叉树六面体网格的二维示意图。
步骤(2)本发明采用线性基(p＝[1,x,y,z])，并生成六面体有限元形函数(公式1)以及形函数的三维空间导数(公式2)，其中ζ,η,μ为局部坐标，f为施加的外力，并以此构建形函数矩阵和应变矩阵。
根据当前域所属的材质属性(杨氏模量和泊松比)构建材质矩阵。依据有限元公式组织局部刚度矩阵(公式3)、质量矩阵(公式4)和施力向量(公式5)，其中公式(3)中的C是材料矩阵，公式(4)中的ρ是仿真对象的密度值。
由于所生成的六面体网格形状是相同的只是尺寸不同，本发明只需要计算上述矩阵一次，通过乘以不同的放缩因子(Scale Facor)得到不同尺寸的对应矩阵。
步骤(3)在进行物理仿真过程中，针对物理仿真的输出进行先验模态分析，生成当前仿真模型的形变子空间，从而实现模型的化简。先验模态分析方法需要及时计算和更新代表当前有限元模型的基向量矩阵U(Basis Matrix)。在仿真过程中，基向量矩阵U的每一列代表了形变子空间的基向量，这些空间基向量是由从仿真输出结果中提取的经验特征向量(Empirical Eigenvectors)构成的。空间基向量的个数(即形变子空间的维度)会根据当前系统的状态进行动态地改变，目的是为了更好地反应当前的物理模型(见图3)。一般地，当前子空间处于构建过程中，空间基向量个数会根据仿真输出的个数进行递增，见算法1；当形变子空间仿真结果的误差大于误差阀值或者物理模型发生了切割操作导致已有的子空间不能表示物理模型时，形变子空间需要进行自适应更新，此时会根据形变空间基向量对物理模型的显著性(Dominant Basis)进行选择，抛弃无效的基向量，并重新启动子空间构建过程，见算法2。
对于形变子空间的构造，本发明采用如下策略：根据已有的子空间特征向量U(n)和新生成的仿真输出结果u，将u对U(n)中的特征向量进行正交操作，得到正交向量uortho，计算uortho的2范数与u的2范数的比值，以确定u对增广当前子空间的贡献是否超过阈值，如果超过阈值，将uortho作为新的特征向量追加到U(n)，得到U(n+1)；没有超过阈值时，不对U(n)进行改变。
对于形变子空间的更新，本发明采用如下策略：根据已有的子空间特征向量U(n)，对U(n)中的每一列c计算其势能式中的K为当前系统的全局刚度矩阵，并对列c根据其势能值进行降序排列，将子空间特征向量个数减半，保留势能最大的特征向量，生成新的U(n+1)。
步骤(4)在进行形变子空间仿真过程中，需要将子空间仿真结果投影到物理模型空间(利用子空间基向量矩阵)，见公式6，其中qsub代表子空间仿真输出，usub代表qsub投影到物理空间的结果，K为当前系统的全局刚度矩阵，DoFs是当前系统的全局自由度；同时，对子空间的仿真结果进行误差分析，当子空间的仿真结果与有限元仿真结果的误差τerror低于阀值τthreshold时，使用子空间仿真方法替换有限元仿真，从而加速形变模型的仿真效率。
usub＝Uqsub (6)
DoFs是当前系统的全局自由度，即全局刚度矩阵K的行数。
步骤(5)针对所采用的线性弹性有限元仿真模型(线性柯西应变张量，CauchyStrain Tensor)只适用于小范围形变的不足，利用局部应变旋转技术，将仿真元素形变过程中所包含的位移与旋转分离，支持大范围形变情况。根据所得到的位移场(DisplacementField)，从中提取仿真点在局部坐标系中的旋转向量，并将距离场中的旋转部分从位移场中分离，见公式(8-12)，其中：k表示当前时间步，代表分离旋转算子，是一个块对角矩阵(Block Diagonal Matrix)，其每一个子矩阵块是一个3×3矩阵(见公式9)，代表每一个顶点的旋转分离算子。公式(10)列出了顶点旋转分离算子的计算公式，其中I表示单位矩阵，wi代表根据顶点i的位移旋度(Curl)得到的局部旋转向量，其计算公式见公式(13)。公式(11)是公式(10)中的展开形式。公式(12)是对公式(10)中向量的斜对称算子×的展开公式。公式(13)是对于有限元素e的由位移引起的旋度计算公式，式中Ne表示当前元素e的形函数矩阵，表示当前元素e的形函数的梯度矩阵，ue表示元素的位移。
步骤(6)针对有限元仿真和子空间仿真交替计算的过程中，子空间仿真会产生累计误差，特别是在动态仿真过程中，累计误差会被记录在系统的速度(Velocity)和加速度(Acceleration)中，通过建立子空间仿真的力残差和最近一次有限元仿真输出的对应关系模型，可以估计当前子空间仿真结果误差的变化趋势。通过公式14，可以得到当前子空间的累计误差增长趋势在未来的step仿真步之内，可以满足误差阀值，式中τthreshold表示系统允许的最大误差阀值，uortho表示最近一次有限元仿真输出向当前子空间基向量的投影残差。在未来的step仿真步之内，系统可以不进行误差分析，直接采用形变子空间仿真输出，如见图7所示。
步骤(7)将切割工具建模为有三角面片组成的曲面，当切割操作发生时，通过六面体网格与三角面片做碰撞检测，可以得到被切割的网格，对于粗层次的网格首先要对其进行网格细分以达到仿真精度，对于被切割的六面体进行克隆操作以满足切割所造成的拓扑变化(见图6中的a-b)，最后对受影响的区域进行步骤(2)和步骤(3)，更新全局方程组。对于切割造成的模型表面的几何变化采用Delaunay三角化进行更新(见图6中的c)。
本发明的实现使用的软件平台为Microsoft visual studio 2010与OpenGL，使用了CUDA来加速并行算法的计算效率。硬件平台为3.4GHz Inter(R)Core(TM)i7-2600CPU、4GB内存以及NVIDIA GeForce GTX570 GPU。方法效果图如图8所示，图8中的(a)显示了系统的初始模型，图8中的(b)显示了系统对图8中的(a)的形变仿真结果，图8中的(c)显示了切割操作发生后系统对图8中的(a)所做的涉及几何、物理的更新操作，此时由于切割操作导致已经生成的形变子空间无法表示切割后的有限元系统，系统选择对图8中的(c)使用有限元系统进行仿真，图8(d)显示了在对新生成的有限元系统进行仿真的同时，根据新生成的形变数据所构建的形变子空间已经满足系统对误差的要求后，系统选择使用形变子空间方法来对图8中的(d)进行仿真。
图5中，左侧的“稀疏矩阵图”描述了稀疏矩阵的逻辑结构(实际并不存在)，中间的“列”和“值”描述了本发明所采用的无结构矩阵的实际数据分布方式，右侧是对这种数据分布方式的示例，四个“局部矩阵”表示四个共享同一顶点的有限元局部矩阵，在左侧的逻辑结构中，四个彩色部分数据是累加为一个数值存储在逻辑结构中，而无结构矩阵并不对其之进行累加，而是将其分别存储在值列表中，从而对四个局部矩阵进行解耦。
Claims (8)
1.一种基于有限元和时变模态分析的柔性物体实时切割仿真方法，其特征在于：包括如下步骤：
步骤(1)将仿真对象的材质域依据仿真精度要求进行层次细节体素化，并生成基于八叉树的六面体网格；
步骤(2)对步骤(1)中的六面体网格采用线性基进行描述，并利用线性柯西应变张量进行有限元建模，构建有限元系统；由于六面体网格的规则性以及线性基和线性柯西应变张量在仿真过程中的不变性，所得到的线性弹性有限元系统能够进行快速组装和仿真；
步骤(3)在仿真过程中，对步骤(2)所构建的有限元系统利用GPU所提供的并行计算能力，并使用支持并行计算的共轭梯度算法进行仿真计算；同时采用Newmark隐式积分技术使得系统支持无条件稳定特性；对有限元系统的仿真输出进行基于先验模态分析的时变子空间动态构建；
步骤(4)当由步骤(3)所构建的时变子空间的仿真结果与步骤(3)的有限元仿真结果的误差值低于阀值时，采用所构建的子空间仿真结果替换步骤(3)的有限元仿真；
步骤(5)根据步骤(3)、(4)输出的系统位移，计算每个仿真顶点的局部旋转，并根据得到的局部旋转对步骤(3)、(4)输出的系统位移进行修正，得到修正后的系统位移；
步骤(6)当没有发生切割操作时，根据步骤(5)输出的系统位移，对当前子空间的仿真精度进行测量，得出子空间基向量的显著性权重；并根据仿真进度的变化计算子空间时变有效性区间；
步骤(7)当切割操作发生时，首先对切割操作所涉及的区域进行识别，通过网格与切割工具进行碰撞检测，对得到的切割所涉及的网格区域进行网格细分以满足系统的仿真精度并通过对所涉及网格进行克隆并沿着切割的轨迹进行网格分裂以体现网格模型的拓扑变化，并在新生成的网格区域上重新建立有限元仿真系统，同时在随后的时间里对新生成的有限元仿真系统进行形变子空间的动态构建。
2.根据权利要求1所述的基于有限元和时变模态分析的柔性物体实时切割仿真方法，其特征在于：所述步骤(1)将仿真对象的材质域依据仿真精度要求进行层次细节体素化并生成基于八叉树的六面体网格的具体实现为：
(11)读取OBJ模型数据，根据仿真对象材质属性的边界曲面函数对模型建立域边界划分投票函数(voting function)，通过该投票函数区分仿真对象的材质域，从而得到该投票函数区所拥有的材质属性；
(12)建立模型整体的AABB包围盒作为基于八叉树的六面体网格根节点，如果当前六面体空间内包含模型元素：顶点、三角面片和材质边界曲面，且没有达到最高层次精度，则对该六面体进行剖分并对新生成的8个子六面体重复上述操作；当所有六面体网格均符合上述条件时，停止操作，所述上述条件为所有声称的六面体包含模型元素：顶点、三角面片和材质边界曲面，且达到最高层次精度；并对八叉树的叶子节点上的六面体元素使用其中心点作为输入参数调用边界划分投票函数确定当前六面体元素是否所属的仿真域；最终生成一个反应模型几何和材质特征的层次化八叉树的六面体网格。
3.根据权利要求1所述的基于有限元和时变模态分析的柔性物体实时切割仿真方法，其特征在于：所述步骤(2)包括：对步骤(1)中生成的六面体网格中的每一个六面体元素采用线性基加以表示，生成六面体元素的形函数以及形函数的三维空间导数，并以此构建形函数矩阵和应变矩阵；根据当前域所属的材质属性构建材质矩阵，所述的材质属性为杨氏模量和泊松比；依据有限元公式组织局部刚度矩阵、局部质量矩阵和局部施力向量；根据六面体元素间的连接关系，将每一个六面体元素的局部刚度矩阵、局部质量矩阵和局部施力向量组装成全局的刚度矩阵、质量矩阵和施力向量，从而得到线性弹性有限元系统。
4.根据权利要求1所述的基于有限元和时变模态分析的柔性物体实时切割仿真方法，其特征在于：所述步骤(3)包括：在有限元系统的仿真过程中，采用支持并行计算的共轭梯度算法对有限元系统方程进行并行求解，采用Newmark隐式积分技术使得仿真系统支持无条件稳定特性；在子空间构建过程中，针对大尺度形变和切割操作所引发的有限元模型的更新，采用先验模态分析技术对子空间进行随时间变化的动态构建与更新，通过将新产生的仿真结果与当前形变子空间合成，动态更新形变子空间的基向量，使其能够准确地反映有限元模型所发生的变化。
5.根据权利要求1所述的基于有限元和时变模态分析的柔性物体实时切割仿真方法，其特征在于：所述步骤(4)包括：当由步骤(3)所构建的时变子空间的仿真结果与步骤(3)的有限元仿真结果的误差值低于阀值时，采用所构建的子空间仿真结果替换步骤(3)的有限元仿真结果，提高系统的计算效能，子空间仿真的计算采用基于Householder变换的QR分解方法直接求解，同时采用Newmark隐式积分技术使得系统支持无条件稳定特性。
6.根据权利要求1所述的基于有限元和时变模态分析的柔性物体实时切割仿真方法，其特征在于：步骤(5)还包括：由于线性弹性系统不能表示旋转运动，导致步骤(3)、(4)的仿真输出中包含了旋转位移，通过计算每个仿真顶点的局部旋转，对仿真输出中的旋转运动加以分离，并将修正后的结果作为系统仿真的位移。
7.根据权利要求1所述的基于有限元和时变模态分析的柔性物体实时切割仿真方法，其特征在于：步骤(6)包括：步骤(5)输出的系统位移是由系统根据对子空间仿真误差的测量结果，在步骤(3)的有限元仿真输出和步骤(4)的子空间仿真输出之间选择得到的；对子空间仿真误差的测量是通过建立子空间仿真的力残差和最近一次有限元仿真输出的对应关系模型，得到子空间时变有效性区间；当仿真时间处于时变有效性区间内时，子空间的仿真精度满足要求，无需对子空间进行更新操作；当仿真时间超出时变有效性区间时，对子空间基向量的显著性进行分析，抛弃无效的基向量，并调用步骤(3)对子空间进行更新，重新计算新的子空间的时变有效性区间。
8.根据权利要求1所述的基于有限元和时变模态分析的柔性物体实时切割仿真方法，其特征在于：步骤(7)包括：当切割操作发生时，将切割工具建模为由三角面片组成的曲面，通过六面体网格与三角面片作碰撞检测，能够对被切割的网格区域进行识别；在被切割的网格区域内进行切口生成的过程中，当被切割的网格处在所对应的八叉树结构中粗层次时，首先要对其进行步骤(1)中的网格细分以达到最小仿真精度，以便对切口部分进行满足仿真精度的表示；对于被切割的六面体进行克隆操做，分裂克隆的新六面体和被克隆的原六面体，使它们分别置于表示切割工具的三角曲面两侧，以满足切割所造成的拓扑变化；最后将进行细分操作和克隆操作的区域更新到有限元系统中，使得有限元仿真系统得到更新；对新生成的有限元仿真系统调用步骤(3)进行仿真计算。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201510751613.3A CN105302974B (zh) | 2015-11-06 | 2015-11-06 | 一种基于有限元和时变模态分析的柔性物体实时切割仿真方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201510751613.3A CN105302974B (zh) | 2015-11-06 | 2015-11-06 | 一种基于有限元和时变模态分析的柔性物体实时切割仿真方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN105302974A CN105302974A (zh) | 2016-02-03 |
| CN105302974B true CN105302974B (zh) | 2018-06-08 |
Family
ID=55200242
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201510751613.3A Active CN105302974B (zh) | 2015-11-06 | 2015-11-06 | 一种基于有限元和时变模态分析的柔性物体实时切割仿真方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN105302974B (zh) |
Families Citing this family (15)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2015035178A2 (en) * | 2013-09-06 | 2015-03-12 | Brigham And Women's Hospital, Inc. | System and method for a tissue resection margin measurement device |
| CN105929943A (zh) * | 2016-04-15 | 2016-09-07 | 郭清锁 | 一种用于计算机虚拟装配的隐式交互系统 |
| CN106227922B (zh) * | 2016-07-14 | 2019-08-20 | 燕山大学 | 在Laplace-Beltrami形状空间基于样例的弹性材料的实时仿真方法 |
| CN106202783B (zh) * | 2016-07-20 | 2021-07-09 | 吉林大学 | 一种驾驶座椅设计方法 |
| CN106991220A (zh) * | 2017-03-23 | 2017-07-28 | 西华大学 | 基于ansys的微型车载转运作业平台的模型构建方法 |
| CN107229804B (zh) * | 2017-06-26 | 2020-09-01 | 中国航发湖南动力机械研究所 | 直升机主减速器润滑系统动态仿真分析方法及装置 |
| CN108089457B (zh) * | 2017-11-29 | 2020-08-14 | 北京航空航天大学 | 一种基于在线有限元仿真的过程质量控制方法 |
| CN108694290B (zh) * | 2018-06-05 | 2021-12-07 | 东北大学 | 一种基于八叉树网格的有限元模型的软组织变形方法 |
| CN108763841B (zh) * | 2018-07-24 | 2022-04-05 | 北京航空航天大学青岛研究院 | 一种基于对偶边界元和应变能优化分析的弹性断裂仿真方法 |
| CN109033641B (zh) * | 2018-07-29 | 2023-05-23 | 南京信息工程大学 | 一种基于硅胶愈合模型的虚拟切割算法 |
| CN110309861B (zh) * | 2019-06-10 | 2021-05-25 | 浙江大学 | 一种基于生成对抗网络的多模态人类活动识别方法 |
| CN110941494A (zh) * | 2019-12-02 | 2020-03-31 | 哈尔滨工程大学 | 一种面向深度学习的gpu并行计算的数据处理方法 |
| CN112487610B (zh) * | 2020-11-09 | 2021-10-08 | 河北工业大学 | 具有复杂几何特征分析对象的形变确定方法及系统 |
| CN113204051B (zh) * | 2021-06-10 | 2022-04-15 | 成都理工大学 | 一种基于变分模态分解的低秩张量地震数据去噪方法 |
| CN114462265A (zh) * | 2021-12-21 | 2022-05-10 | 深圳先进技术研究院 | 一种用于可变形物体剪切断裂的仿真方法及材料仿真方法 |
Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US6236738B1 (en) * | 1998-04-09 | 2001-05-22 | Board Of Trustees Of The Leland Stanford Junior University | Spatiotemporal finite element method for motion analysis with velocity data |
| CN103699714A (zh) * | 2013-12-01 | 2014-04-02 | 北京航空航天大学 | 一种基于有限元和无网格耦合的柔性物体实时切割仿真方法 |
-
2015
- 2015-11-06 CN CN201510751613.3A patent/CN105302974B/zh active Active
Patent Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US6236738B1 (en) * | 1998-04-09 | 2001-05-22 | Board Of Trustees Of The Leland Stanford Junior University | Spatiotemporal finite element method for motion analysis with velocity data |
| CN103699714A (zh) * | 2013-12-01 | 2014-04-02 | 北京航空航天大学 | 一种基于有限元和无网格耦合的柔性物体实时切割仿真方法 |
Non-Patent Citations (2)
| Title |
|---|
| Real-time physical deformation and cutting of heterogeneous objects via hybrid coupling of meshless approach and finite element method;Chen Yang等;《COMPUTER ANIMATION AND VIRTUAL WORLDS》;20141231;全文 * |
| 基于八叉树的三维网格模型体素化方法;吴晓军等;《工程图学学报》;20051231;全文 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN105302974A (zh) | 2016-02-03 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN105302974B (zh) | 一种基于有限元和时变模态分析的柔性物体实时切割仿真方法 | |
| Nesme et al. | Preserving topology and elasticity for embedded deformable models | |
| Teran et al. | Adaptive physics based tetrahedral mesh generation using level sets | |
| CN103699714B (zh) | 一种基于有限元和无网格耦合的柔性物体实时切割仿真方法 | |
| CN101719284B (zh) | 一种基于层次模型的虚拟人皮肤物理变形方法 | |
| Iben et al. | Generating surface crack patterns | |
| CN104268943A (zh) | 一种基于欧拉-拉格朗日耦合方法的流体仿真方法 | |
| CN101398942A (zh) | 三维试衣仿真系统 | |
| Guo et al. | Meshless thin-shell simulation based on global conformal parameterization | |
| CN106934192B (zh) | 一种参数优化的浅水方程模型水体建模方法 | |
| CN105354879A (zh) | 基于质点弹簧结构的通用服装三维模型仿真方法及系统 | |
| CN101833785A (zh) | 一种具有物理真实感的可控动态形状插值方法 | |
| CN103268634A (zh) | 一种外存模型基于顶点聚类的快速并行自适应简化方法 | |
| CN103426196B (zh) | 一种流体环境下的关节动画建模方法 | |
| CN111341449A (zh) | 一种虚拟血管介入手术训练的模拟方法 | |
| KR101178443B1 (ko) | 의상 시뮬레이션을 위한 물리법칙에 기반한 멀티그리드 방법 및 그 프로그램이 기록된 컴퓨터가 판독가능한 기록매체 | |
| Merrell et al. | Constraint-based model synthesis | |
| CN109118561B (zh) | 一种基于位置的层次化动态模拟方法 | |
| Lee et al. | A Fast and Compact Solver for the Shallow Water Equations. | |
| CN101655832B (zh) | 一种基于标量场梯度的物理变形方法 | |
| Aldrich et al. | Collision-Driven Volumetric Deformation on the GPU. | |
| Oh et al. | Practical simulation of hierarchical brittle fracture | |
| Guo et al. | Meshless methods for physics-based modeling and simulation of deformable models | |
| Li et al. | An object-oriented system for dynamics-based 3D cloth simulation | |
| McDonnell et al. | Dynamic subdivision-based solid modeling |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| C06 | Publication | ||
| PB01 | Publication | ||
| C10 | Entry into substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |