CN109271696B - 基于mpm的血液凝固模拟方法及系统 - Google Patents
基于mpm的血液凝固模拟方法及系统 Download PDFInfo
- Publication number
- CN109271696B CN109271696B CN201811045362.7A CN201811045362A CN109271696B CN 109271696 B CN109271696 B CN 109271696B CN 201811045362 A CN201811045362 A CN 201811045362A CN 109271696 B CN109271696 B CN 109271696B
- Authority
- CN
- China
- Prior art keywords
- speed
- particle
- mesh point
- blood
- nonce
- 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
- 238000000034 method Methods 0.000 title claims abstract description 74
- 230000023555 blood coagulation Effects 0.000 title claims abstract description 54
- 239000002245 particle Substances 0.000 claims abstract description 143
- 210000004369 blood Anatomy 0.000 claims abstract description 46
- 239000008280 blood Substances 0.000 claims abstract description 46
- 238000004088 simulation Methods 0.000 claims abstract description 39
- 239000012530 fluid Substances 0.000 claims description 52
- 239000007787 solid Substances 0.000 claims description 21
- 238000004364 calculation method Methods 0.000 claims description 10
- 238000007711 solidification Methods 0.000 abstract description 16
- 230000008023 solidification Effects 0.000 abstract description 16
- 230000008569 process Effects 0.000 abstract description 12
- 239000007788 liquid Substances 0.000 description 16
- 230000006872 improvement Effects 0.000 description 11
- 239000000126 substance Substances 0.000 description 11
- 239000012071 phase Substances 0.000 description 8
- 238000010586 diagram Methods 0.000 description 7
- 230000008859 change Effects 0.000 description 6
- 230000001133 acceleration Effects 0.000 description 5
- 230000008901 benefit Effects 0.000 description 5
- 230000009466 transformation Effects 0.000 description 5
- 230000005484 gravity Effects 0.000 description 4
- 210000004027 cell Anatomy 0.000 description 3
- 238000006243 chemical reaction Methods 0.000 description 3
- 238000013461 design Methods 0.000 description 3
- 238000005516 engineering process Methods 0.000 description 3
- 239000007791 liquid phase Substances 0.000 description 3
- 239000000463 material Substances 0.000 description 3
- 230000007704 transition Effects 0.000 description 3
- 102000009123 Fibrin Human genes 0.000 description 2
- 108010073385 Fibrin Proteins 0.000 description 2
- BWGVNKXGVNDBDI-UHFFFAOYSA-N Fibrin monomer Chemical compound CNC(=O)CNC(=O)CN BWGVNKXGVNDBDI-UHFFFAOYSA-N 0.000 description 2
- 230000015271 coagulation Effects 0.000 description 2
- 238000005345 coagulation Methods 0.000 description 2
- 229950003499 fibrin Drugs 0.000 description 2
- 239000011159 matrix material Substances 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 238000012545 processing Methods 0.000 description 2
- 102000008946 Fibrinogen Human genes 0.000 description 1
- 108010049003 Fibrinogen Proteins 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- 238000006073 displacement reaction Methods 0.000 description 1
- 229940012952 fibrinogen Drugs 0.000 description 1
- 230000003993 interaction Effects 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 239000000155 melt Substances 0.000 description 1
- 238000012544 monitoring process Methods 0.000 description 1
- 238000005293 physical law Methods 0.000 description 1
- 210000002381 plasma Anatomy 0.000 description 1
- 230000009024 positive feedback mechanism Effects 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 230000004044 response Effects 0.000 description 1
- 239000007921 spray Substances 0.000 description 1
- 238000010257 thawing Methods 0.000 description 1
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 1
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
- G06T13/00—Animation
- G06T13/20—3D [Three Dimensional] animation
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2210/00—Indexing scheme for image generation or computer graphics
- G06T2210/24—Fluid dynamics
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- General Engineering & Computer Science (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明公开了一种基于MPM的血液凝固模拟方法，包括：S1，初始化通用数据；S2，将血液离散化成粒子和网格；S3，根据物质点法模拟血液凝固。本发明还公开了一种基于MPM的血液凝固模拟方法系统。采用本发明，可基于MPM框架，根据凝固的特点结合连续介质力学以及流体力学N‑S方程，并且根据血液凝固速度设计拉梅系数函数，来精确地模拟血液凝固过程。
Description
技术领域
本发明涉及计算机图形学技术领域，尤其涉及一种基于MPM的血液凝固模拟方法及一种基于MPM的血液凝固模拟系统。
背景技术
由于现实中难以使用大量的真实血液来进行电影特效的创作以及游戏动画中血液的模拟，所以血液凝固的模拟通常被电影特效工业以及游戏动画行业所使用，以计算机模拟的血液对电影场景中的血液进行替代，其中主要涉及的技术包括流体模拟技术以及相变模拟技术。
流体模拟是指使用计算机对现实世界中流体进行模拟，以在计算机中对现实世界中难以实现的大规模流体、特殊效果以及数值监测进行数值模拟或视觉呈现。它通常将流体遵循的物理定律使用数学公式表示，并针对模拟进行具体的公式简化或者公式变形，并通过计算机编程实现。流体模拟的难点在于如何快速而又不失真实感地使用计算机实现，以用于实时模拟。
在电影工业，越来越多的流体被应用到电影特效的制作中，从烟雾缭绕的场景到巨浪拍下，从漂亮的水花溅起到熔岩的缓慢流下，由于娱乐产业上的应用，流体模拟渐渐成为计算机图形学重要的研究方向。流体模拟都需要从纳维斯托克斯方程(Navier-Stokesequation，简称为N-S方程)开始建模，它是对自然界流体微元进行受力分析所得到的偏微分方程。
目前流体模拟的方法主要分为三种：欧拉方法(Euler)，拉格朗日方法(Lagrange)以及两种方法的混合方法。
如图1所示，在欧拉方法中，N-S方程已经构建好了微分方程框架，需要进行流体模拟的有限空间被离散化为有限数量在空间上固定的网格点，那么整个模拟过程流体总是会流过网格点，那么在某一个时间点每个网格点就可以存储流经该点的流体的属性值，包括速度、温度、压力梯度、密度等。由于N-S方程已经构建好了微分方程框架，只需通过求微分方程数值解的通用方法，把N-S方程离散化为差分方程，再使用显式、隐式或半隐式的方法结合有限差分法来进行多时间步迭代计算，每一个时间步迭代更新所有网格点的属性值，即每一个时间步生成一帧流体画面。另一方面，欧拉方法可以很容易地和流体不可压缩条件结合，而整体框架不需要更改。
如图2所示，拉格朗日方法则是从另一个角度去解决N-S方程，流体被看成是一个粒子系统，流体是由许多细小的粒子组成的，流体的流动即是每一个粒子的流动，每一个流体粒子都有其固有的属性值，包括速度、位置等，那么每当流体流动时，所有粒子的速度和位置随之变化。与欧拉方法对比，欧拉方法是通过空间中固定的网格点来表示此位置流体的性质，网格点不会随着流体流动而变换位置，拉格朗日方法则追踪每一个流体粒子来模拟流体。
混合方法结合了欧拉方法便于处理边界条件和微分的优点和拉格朗日方法数值耗散小的优点，同时使用网格和粒子，并通过网格和粒子的相互转换来模拟流体。该方法在网格点计算微分，这是在拉格朗日方法中难以精确计算的，然后把在网格点得到的值映射回粒子上，再对粒子进行速度的积分，即可得到位移。由此便可以形成一帧画面，循环计算既可以形成流体动画。
相变模拟与流体模拟类似，是通过计算机对现实世界中的相变进行数值模拟或视觉实现。目前这个方向的模拟主要有流体转化成固体的模拟，例如水的凝固、岩浆的凝固等，以及固体转化为流体的模拟，例如冰的融化、固体油受热融化等，通常情况下固液的互相转换可以由同一个模型进行描述，这样固体转化为液体以及逆向转换只需要对模型参数进行逆向变化即可实现。
现有的相变模拟技术多数是基于流体模拟中的粒子方法实现的，最典型的的粒子模拟方法是SPH(Smoothed Particle Hydrodynamics，光滑粒子流体动力学)方法，相变模拟在粒子模拟方法的基础上对N-S方程中的粘度系数进行重新建模，结合具体相变模拟情况，来对粘度系数进行数值上的限制以及建立函数关联，以此实现固液之间的转换。
现有技术主要是对原来模拟流体的SPH方法的推广和改进，使得SPH方法可以用来模拟固液之间的相变。Afonso Paiva等人在使用SPH模拟固液相变方面提出了一种方法，使得SPH方法可以模拟固体以及液体以及固液之间的相变。他的方法是基于粘塑性流体和场景交互对SPH框架进行推广，其中使用通用的牛顿流体模型框架，使得SPH框架可以模拟粘塑性物质，并且使用了模拟非牛顿流体的粘度模型来构造粘度函数，这个模型只需要使用一个参数jump number来调整粘塑性模型，易于操控，不需要考虑更多的技术参数，包括刚度、可压缩性、塑性、粘性、内聚力等。即在牛顿流体模型框架上加入了非牛顿流体粘度函数模型，以此来模拟非牛顿流体以及粘塑性物质，并以SPH框架进行实现
如图3所示，现有SPH模拟固液相变方法的具体流程如下：
(1)根据现实中的物理量，初始化粒子的位置、速度、质量、体积等基本信息，为运行程序做准备；
(2)对所有粒子更新N-S方程的不同项，即每个粒子的速度梯度、变形张量率、密度微分、压力、粘度；
(3)由上一步所得到的数值，根据N-S方程更新所有粒子的加速度，并加入人工粘度以防止数值上的不稳定；
(4)对所有粒子使用蛙跳格式更新粒子速度和粒子密度，使用XSPH校正粒子速度以防止粒子相互渗透，根据速度更新粒子经过最后碰撞后的位置；
(5)根据CFL条件更新时间步长，并把上一步的时间步长加入计算所耗时间中；
(6)回到第(1)步，直到计算耗时达到设定的时间上限为止。
但是，该方法使用了高粘度非牛顿流体模型来模拟粘塑性固体物质，并没有使用真实反应物质特性的连续介质力学模型，问题在于高粘度非牛顿流体模型只能用于模拟流体，不能真实地反映出粘塑性固体的物质特性，只能以一种近似的方式模拟粘塑性固体，并且不能模拟粘弹性固体物质。同时，该方法中，粒子的速度、速度梯度、压力梯度、密度等粒子的基本量是通过邻域粒子的值来计算的，所以每次更新一个粒子的值时，需要搜索这个粒子对应的邻域粒子，才能对这个粒子的相应值进行更新，因此，每次更新所有粒子的特定值时，都需要重新搜索，这样时间复杂度非常高，需要花费额外的搜索时间。另外，该方法中粘塑性仅仅是通过改变N-S流体方程中的粘性系数，并且加上一个粘性与jump number的关联函数来实现粘塑性效果，并不能针对血液凝固进行凝固过程中从液态到中间态再到固态的精确模拟。
发明内容
本发明所要解决的技术问题在于，提供一种基于MPM的血液凝固模拟方法及系统，可基于MPM框架，根据凝固的特点结合连续介质力学以及流体力学N-S方程，并且根据血液凝固速度设计拉梅系数函数，来精确地模拟血液凝固过程。
为了解决上述技术问题，本发明提供了一种基于MPM的血液凝固模拟方法，包括：
S1，初始化通用数据；
S2，将血液离散化成粒子和网格；
S3，根据物质点法模拟血液凝固。
作为上述方案的改进，所述步骤S3包括：
S31，初始化粒子数据；
S32，将粒子的数据映射到网格点上，并计算粒子的密度和体积；
S33，进行网格点上的计算；
S34，更新粒子变形梯度；
S35，更新粒子的速度；
S36，处理粒子的碰撞；
S37，更新粒子的位置。
作为上述方案的改进，所述步骤S33包括：
S331，计算网格点的所受的力；
S332，根据所受的力更新网格点的速度临时值；
S333，处理网格点的碰撞，并更新网格点的速度临时值；
S334，更新速度临时值为网格点的最终速度。
作为上述方案的改进，所述步骤S332包括：当血液为流体时，采用N-S方程更新网格点的速度临时值；当血液为固体时，采用连续介质力学更新网格点的速度临时值。
作为上述方案的改进，所述步骤S334包括：根据血液状态更新速度临时值为网格点的最终速度。
相应地，本发明还提供了一种基于MPM的血液凝固模拟系统，包括：初始化模块，用于初始化通用数据；离散化模块，用于将血液离散化成粒子和网格；模拟模块，用于根据物质点法模拟血液凝固。
作为上述方案的改进，所述模拟模块包括：初始化单元，用于初始化粒子数据；映射单元，用于将粒子的数据映射到网格点上，并计算粒子的密度和体积；计算单元，用于进行网格点上的计算；梯度更新单元，用于更新粒子变形梯度；速度更新单元，用于更新粒子的速度；碰撞单元，用于处理粒子的碰撞；位置更新单元，用于更新粒子的位置。
作为上述方案的改进，所述计算单元包括：受力计算子单元，用于计算网格点的所受的力；临时值计算子单元，用于根据所受的力更新网格点的速度临时值；碰撞计算子单元，用于处理网格点的碰撞，并更新网格点的速度临时值；速度计算子单元，用于更新速度临时值为网格点的最终速度。
作为上述方案的改进，所述临时值计算子单元包括：第一临时值计算子单元，用于当血液为流体时，采用N-S方程更新网格点的速度临时值；第二临时值计算子单元，用于当血液为固体时，采用连续介质力学更新网格点的速度临时值。
作为上述方案的改进，所述速度计算子单元包括：第一速度计算子单元，用于当血液为流体时，更新速度临时值为网格点的最终速度；第二速度计算子单元，用于当血液为固体时，更新速度临时值为网格点的最终速度。
本发明可基于MPM框架，根据凝固的特点结合连续介质力学以及流体力学N-S方程，并且根据血液凝固速度设计拉梅系数函数，来精确地模拟血液凝固过程。具体地，本发明具有以下有益效果：
(1)结合能精确模拟血液凝固纤维蛋白的连续介质力学模型和模拟血液的流体力学模型来模拟血液凝固的整个过程，其中，使用流体力学模型来模拟血液的流体状态，使用连续介质力学模型来模拟血液的凝胶体状态，并且针对血液凝固的非线性特性，提出了流体力学模型向连续介质力学转换的非线性模型，使得凝血过程中的固态、液态可以精确地模拟出来；
(2)使用MPM来模拟血液，MPM使用了网格作为计算的对象，由粒子把物理量映射到网格上，在网格上完成计算后再把速度映射回粒子上，粒子再根据返回的速度进行位置的更新，避免SPH方法的高时间复杂度的邻域搜索。
(3)针对血液凝固的速度，设计凝固随时间变化的函数，即连续介质力学中的拉梅系数随时间变化的函数，使得本发明可以模拟血液凝固中液态转化成固态的中间态，而不是直接由液态转化成固态。
附图说明
图1是二维情况下欧拉方法的示意图；
图2是二维情况下拉格朗日方法的示意图；
图3是SPH模拟固液相变方法的流程图；
图4是本发明基于MPM的血液凝固模拟方法的流程图；
图5是图4中根据物质点法模拟血液凝固的流程图；
图6是图5中进行网格点上的计算的流程图；
图7是本发明基于MPM的血液凝固模拟系统的结构示意图；
图8是图7中模拟模块的结构示意图；
图9是图8中计算单元的结构示意图；
图10是图9中临时值计算子单元的结构示意图；
图11是图9中速度计算子单元的结构示意图。
具体实施方式
为使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明作进一步地详细描述。仅此声明，本发明在文中出现或即将出现的上、下、左、右、前、后、内、外等方位用词，仅以本发明的附图为基准，其并不是对本发明的具体限定。
参见图4，图4显示了本发明基于MPM的血液凝固模拟方法的流程图，其包括：
S1，初始化通用数据；
初始化凝固开始时间点为tcogBegin，凝固结束时间点为tcogEnd，单个网格的边长h，粒子之间的间距dparticle，拉梅系数拉梅系数迭代时间步长Δt。
S2，将血液离散化成粒子和网格；
根据血液所占空间，把该空间离散化为整齐排列的粒子，并离散化为有限个三维整齐排列的网格。
例如，总空间为一个边长为1的正方体，单个网格边长为0.1，粒子之间的间距为0.01，那么总空间离散化为10×10×10个网格，100×100×100个粒子，其中网格的位置随时间变化是固定的，粒子的位置随时间变化是流动的。
S3，根据物质点法(MPM，Material Point Method)模拟血液凝固。
如图5所示，所述步骤S3包括：
S31，初始化粒子数据；
初始化粒子质量mp，粒子速度vp，变形张量Fp＝I,FEp＝I,FPp＝I，p表示粒子的编号。
S32，将粒子的数据映射到网格点上，计算网格点的质量和速度并计算粒子的密度ρp和体积Vp；
根据公式(1)计算第n步网格点的质量
其中，mp为编号为p的粒子的质量，为第i个网格点和第p个粒子之间的权重，
其中，xp,yp,zp为粒子的三维坐标，i,j,k为网格的三维顺序编号，h为正方体网格的大小，N(x)为三维b样条函数，
根据公式(4)计算网格点的速度
其中，是第n步粒子的速度，mp、mi、的意义与公式(1)中一致。
若迭代步为第一步，根据公式(5)计算第一步粒子的密度
其中，h的意义与公式(1)(2)中的意义一致。
若迭代步为第一步，根据公式(6)计算粒子的体积
其中mp、的意义与公式(5)中的意义一致。
S33，进行网格点上的计算；
如图6所示，所述步骤S33包括：
S331，计算网格点的所受的力
根据公式(7)计算网格i所受的力
其中，为粒子p的体积，为粒子p的能量密度，Fp为变形张量，FEp为变形张量的弹性部分，FPp为变形张量的塑性部分，他们之间满足关系Fp＝FEpFPp，为权重的梯度值，为新的变形梯度的弹性部分，
其中，为网格点的虚拟位置，为权重的梯度，为旧的变形梯度的弹性部分。
Ψ(FE,FP)为能量密度函数，
其中，μ和λ为拉梅系数，RE为FE的极分解得到的正交矩阵，JE为FE的行列式，
其中，tcogBegin为凝固开始的时间，tcogEnd为凝固结束的时间，λ0为血液凝固前的拉梅系数，μ0、λ1为血液凝固后的拉梅系数。
S332，根据所受的力更新网格点的速度临时值vi *；
具体地，所述步骤S332包括：
当血液为流体时，采用N-S方程更新网格点的速度临时值v*；
若t<tcogEnd，那么采用N-S方程更新速度，根据公式(12)计算网格点的速度临时值v*，
其中，vn为上一迭代步的速度，Δt为每一步迭代的时间步长，ρn为上一迭代步的密度，σμ为对应的柯西应力张量，g为重力加速度。
根据公式(13)计算Ψμ对应的柯西应力张量σμ，
其中，Ψμ为能量密度函数的塑性惩罚部分，J为变形张量F的行列式，JE为变形张量弹性部分FE的行列式，FE为变形张量的弹性部分。
当血液为固体时，采用连续介质力学更新网格点的速度临时值
若t≥tcogEnd，则采用连续介质力学更新速度，根据公式(14)计算网格点的速度临时值
其中，为上一迭代步的速度，Δt为每一步迭代的时间步长，mi为网格点的质量，为网格点所受的力，g为重力加速度。
S333，处理网格点的碰撞，并更新网格点的速度临时值
若t<tcogEnd，则跳过这步，否则使用水平集方法判断物质距离碰撞物体的距离φ，当距离为负数φ<0时，根据公式(15)计算相对速度在法线方向上的分速度vn，
vn＝vrel·n (15)
其中，vrel＝v-vco为物质与碰撞物体的相对速度，vco为碰撞物体的速度，n为法线方向。若vn大于0，则不发生碰撞，若vn小于0，那么根据公式(16)计算更新的相对速度v′rel，
v′rel＝vt+μfricvnvt/||vt||
(16)
其中，vt＝v-vn为物质在切线方向上的速度，μfric为动摩擦系数。
根据公式(17)计算碰撞后的速度v′，
v′＝v′rel+vco (17)
其中，vco为碰撞物体的速度。
S334，更新速度临时值为网格点的最终速度vi。
具体地，所述步骤S334包括：根据血液状态更新速度临时值为网格点的最终速度。
若t<tcogEnd，根据公式(18)计算下一迭代步网格点速度vn+1，
其中，Δt为迭代时间步长，ρn为密度，JP为变形张量FP的行列式，λ为拉梅系数，JE为变形张量FE的行列式。
若t≥tcogEnd，根据公式(19)计算下一迭代步网格点速度vn+1，
其中，为网格点的速度临时值。
S34，更新粒子变形梯度Fp，FEp，FPp；
根据公式(20)计算下一个迭代步的粒子变形张量
其中，Δt为迭代时间步长，为速度梯度，分别为变形张量的弹性部分和塑性部分。
根据公式(21)计算下一个迭代步的粒子变形张量的弹性部分
其中，Δt为迭代时间步长，为速度梯度，分别为上一个迭代步变形张量的弹性部分。
根据公式(22)计算下一个迭代步的粒子变形张量的塑性部分
其中，为新的变形张量的弹性部分，为新的粒子变形张量。
S35，更新粒子的速度vp；
根据公式(23)计算所有粒子下一个迭代步的速度，
其中，为权重值，为网格点下一个迭代步的速度，分别为网格点和粒子上一个迭代步的速度。
S36，处理粒子的碰撞；
与处理网格点碰撞后速度的步骤一样，只需把网格点的速度替换成粒子的速度这里不再详细说明。
S37，更新粒子的位置xp。
根据公式(24)计算粒子在速度作用下的新位置
其中，为粒子上一步的位置，Δt为迭代时间步长，为粒子下一迭代步的速度。
由上可知，本发明的步骤S331模拟血液凝固的中间态，根据血液凝固过程血浆中的纤维蛋白原转化为纤维蛋白是一个正反馈机制的特点，设计了连续介质力学模型中的拉梅系数函数，使得凝固的速度随时间变化越来越快，符合血液凝固的特点；本发明的步骤S332、S333及S334根据血液凝固的时间特点，灵活结合连续介质力学模型和流体力学N-S模型，在血液为流体时使用流体力学N-S模型模拟血液，在血液为固体时使用连续介质力学模型模拟凝血纤维蛋白，并且结合MPM框架避免了拉格朗日方法计算微分不准确的问题，提高了血液凝固模拟的真实感。
参见图7，图7显示了本发明基于MPM的血液凝固模拟系统100的具体结构，其包括初始化模块1、离散化模块2及模拟模块3，具体地：
初始化模块1，用于初始化通用数据；所示通用数据包括：开始时间点tcogBegin，凝固结束时间点tcogEnd，单个网格的边长h，粒子之间的间距dparticle，拉梅系数拉梅系数迭代时间步长Δt，但不以此为限制。
离散化模块2，用于将血液离散化成粒子和网格；即根据血液所占空间，把该空间离散化为整齐排列的粒子，并离散化为有限个三维整齐排列的网格。
模拟模块3，用于根据物质点法(MPM，Material Point Method)模拟血液凝固。
如图8所示，所述模拟模块3包括初始化单元31、映射单元32、计算单元33、梯度更新单元34、速度更新单元35、碰撞单元36及位置更新单元37，具体地：
初始化单元31，用于初始化粒子数据；所示粒子数据包括：粒子质量mp，粒子速度vp，变形张量Fp＝I,FEp＝I,FPp＝I，粒子的编号p，但不以此为限制。
映射单元32，用于将粒子的数据映射到网格点上，并计算粒子的密度和体积；
具体地，映射单元32根据公式(1)计算第n步网格点的质量
其中，mp为编号为p的粒子的质量，为第i个网格点和第p个粒子之间的权重，
其中，xp,yp,zp为粒子的三维坐标，i,j,k为网格的三维顺序编号，h为正方体网格的大小，N(x)为三维b样条函数，
根据公式(4)计算网格点的速度
其中，是第n步粒子的速度，mp、mi、的意义与公式(1)中一致。
若迭代步为第一步，根据公式(5)计算第一步粒子的密度
其中，h的意义与公式(1)(2)中的意义一致。
若迭代步为第一步，根据公式(6)计算粒子的体积
其中mp、的意义与公式(5)中的意义一致。
计算单元33，用于进行网格点上的计算；
如图9所示，所述计算单元33包括受力计算子单元331、临时值计算子单元332、碰撞计算子单元333及速度计算子单元334，具体地：
受力计算子单元331，用于计算网格点的所受的力；
受力计算子单元331根据公式(7)计算网格i所受的力
其中，为粒子p的体积，为粒子p的能量密度，Fp为变形张量，FEp为变形张量的弹性部分，FPp为变形张量的塑性部分，他们之间满足关系Fp＝FEpFPp，为权重的梯度值，为新的变形梯度的弹性部分，
其中，为网格点的虚拟位置，为权重的梯度，为旧的变形梯度的弹性部分。
Ψ(FE,FP)为能量密度函数，
其中，μ和λ为拉梅系数，RE为FE的极分解得到的正交矩阵，JE为FE的行列式，
其中，tcogBegin为凝固开始的时间，tcogEnd为凝固结束的时间，λ0为血液凝固前的拉梅系数，μ0、λ1为血液凝固后的拉梅系数。
临时值计算子单元332，用于根据所受的力更新网格点的速度临时值；
如图10所示，所述临时值计算子单元332包括第一临时值计算子单元3321及第二临时值计算子单元3322，其中：
第一临时值计算子单元3321，用于当血液为流体时，采用N-S方程更新网格点的速度临时值；
若t<tcogEnd，那么采用N-S方程更新速度，根据公式(12)计算网格点的速度临时值v*，
其中，vn为上一迭代步的速度，Δt为每一步迭代的时间步长，ρn为上一迭代步的密度，σμ为对应的柯西应力张量，g为重力加速度。
根据公式(13)计算Ψμ对应的柯西应力张量σμ，
其中，Ψμ为能量密度函数的塑性惩罚部分，J为变形张量F的行列式，JE为变形张量弹性部分FE的行列式，FE为变形张量的弹性部分。
第二临时值计算子单元3322，用于当血液为固体时，采用连续介质力学更新网格点的速度临时值。
若t≥tcogEnd，则采用连续介质力学更新速度，根据公式(14)计算网格点的速度临时值
其中，为上一迭代步的速度，Δt为每一步迭代的时间步长，mi为网格点的质量，为网格点所受的力，g为重力加速度。
碰撞计算子单元333，用于处理网格点的碰撞，并更新网格点的速度临时值；
若t<tcogEnd，则跳过这步，否则使用水平集方法判断物质距离碰撞物体的距离φ，当距离为负数φ<0时，根据公式(15)计算相对速度在法线方向上的分速度vn，
vn＝vrel·n (15)
其中，vrel＝v-vco为物质与碰撞物体的相对速度，vco为碰撞物体的速度，n为法线方向。若vn大于0，则不发生碰撞，若vn小于0，那么根据公式(16)计算更新的相对速度v′rel，
v′rel＝vt+μfricvnvt/||vt||
(16)
其中，vt＝v-vn为物质在切线方向上的速度，μfric为动摩擦系数。
根据公式(17)计算碰撞后的速度v′，
v′＝v′rel+vco (17)
其中，vco为碰撞物体的速度。
速度计算子单元334，用于更新速度临时值为网格点的最终速度。
如图11所示，所述速度计算子单元334包括：
第一速度计算子单元3341，用于当血液为流体时，更新速度临时值为网格点的最终速度；
若t<tcogEnd，根据公式(18)计算下一迭代步网格点速度vn+1，
其中，Δt为迭代时间步长，ρn为密度，JP为变形张量FP的行列式，λ为拉梅系数，JE为变形张量FE的行列式。
第二速度计算子单元3342，用于当血液为固体时，更新速度临时值为网格点的最终速度。
若t≥tcogEnd，根据公式(19)计算下一迭代步网格点速度vn+1，
其中，为网格点的速度临时值。
梯度更新单元34，用于更新粒子变形梯度；
根据公式(20)计算下一个迭代步的粒子变形张量
其中，Δt为迭代时间步长，为速度梯度，分别为变形张量的弹性部分和塑性部分。
根据公式(21)计算下一个迭代步的粒子变形张量的弹性部分
其中，Δt为迭代时间步长，为速度梯度，分别为上一个迭代步变形张量的弹性部分。
根据公式(22)计算下一个迭代步的粒子变形张量的塑性部分
其中，为新的变形张量的弹性部分，为新的粒子变形张量。
速度更新单元35，用于更新粒子的速度；
根据公式(23)计算所有粒子下一个迭代步的速度，
其中，为权重值，为网格点下一个迭代步的速度，分别为网格点和粒子上一个迭代步的速度。
碰撞单元36，用于处理粒子的碰撞；
与处理网格点碰撞后速度的步骤一样，只需把网格点的速度替换成粒子的速度这里不再详细说明。
位置更新单元37，用于更新粒子的位置。
根据公式(24)计算粒子在速度作用下的新位置
其中，为粒子上一步的位置，Δt为迭代时间步长，为粒子下一迭代步的速度。
由上可知，本发明可基于MPM框架，根据凝固的特点结合连续介质力学以及流体力学N-S方程，并且根据血液凝固速度设计拉梅系数函数，来精确地模拟血液凝固过程。具体地，本发明具有以下有益效果：
(1)结合能精确模拟血液凝固纤维蛋白的连续介质力学模型和模拟血液的流体力学模型来模拟血液凝固的整个过程，其中，使用流体力学模型来模拟血液的流体状态，使用连续介质力学模型来模拟血液的凝胶体状态，并且针对血液凝固的非线性特性，提出了流体力学模型向连续介质力学转换的非线性模型，使得凝血过程中的固态、液态可以精确地模拟出来；
(2)使用MPM来模拟血液，MPM使用了网格作为计算的对象，由粒子把物理量映射到网格上，在网格上完成计算后再把速度映射回粒子上，粒子再根据返回的速度进行位置的更新，避免SPH方法的高时间复杂度的邻域搜索。
(3)针对血液凝固的速度，设计凝固随时间变化的函数，即连续介质力学中的拉梅系数随时间变化的函数，使得本发明可以模拟血液凝固中液态转化成固态的中间态，而不是直接由液态转化成固态。
以上所述是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也视为本发明的保护范围。
Claims (4)
1.一种基于MPM的血液凝固模拟方法，其特征在于，包括：
S1，初始化通用数据；
S2，将血液离散化成粒子和网格；
S3，根据物质点法模拟血液凝固；
所述步骤S3包括：
S31，初始化粒子数据；
S32，将粒子的数据映射到网格点上，并计算粒子的密度和体积；
S33，进行网格点上的计算；其中，所述步骤S33包括：S331，计算网格点的所受的力；S332，根据所受的力更新网格点的速度临时值，具体地，当血液为流体时，采用N-S方程更新网格点的速度临时值；当血液为固体时，采用连续介质力学更新网格点的速度临时值；S333，处理网格点的碰撞，并更新网格点的速度临时值；S334，更新速度临时值为网格点的最终速度；
S34，更新粒子变形梯度；
S35，更新粒子的速度；
S36，处理粒子的碰撞；
S37，更新粒子的位置。
2.如权利要求1所述的基于MPM的血液凝固模拟方法，其特征在于，所述步骤S334包括：根据血液状态更新速度临时值为网格点的最终速度。
3.一种基于MPM的血液凝固模拟系统，其特征在于，包括：
初始化模块，用于初始化通用数据；
离散化模块，用于将血液离散化成粒子和网格；
模拟模块，用于根据物质点法模拟血液凝固；
所述模拟模块包括：
初始化单元，用于初始化粒子数据；
映射单元，用于将粒子的数据映射到网格点上，并计算粒子的密度和体积；
计算单元，用于进行网格点上的计算；其中，所述计算单元包括：受力计算子单元，用于计算网格点的所受的力；临时值计算子单元，用于根据所受的力更新网格点的速度临时值，具体地，当血液为流体时，采用N-S方程更新网格点的速度临时值；当血液为固体时，采用连续介质力学更新网格点的速度临时值；碰撞计算子单元，用于处理网格点的碰撞，并更新网格点的速度临时值；速度计算子单元，用于更新速度临时值为网格点的最终速度；
梯度更新单元，用于更新粒子变形梯度；
速度更新单元，用于更新粒子的速度；
碰撞单元，用于处理粒子的碰撞；
位置更新单元，用于更新粒子的位置。
4.如权利要求3所述的基于MPM的血液凝固模拟系统，其特征在于，所述速度计算子单元包括：
第一速度计算子单元，用于当血液为流体时，更新速度临时值为网格点的最终速度；
第二速度计算子单元，用于当血液为固体时，更新速度临时值为网格点的最终速度。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201811045362.7A CN109271696B (zh) | 2018-09-07 | 2018-09-07 | 基于mpm的血液凝固模拟方法及系统 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201811045362.7A CN109271696B (zh) | 2018-09-07 | 2018-09-07 | 基于mpm的血液凝固模拟方法及系统 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN109271696A CN109271696A (zh) | 2019-01-25 |
| CN109271696B true CN109271696B (zh) | 2019-07-23 |
Family
ID=65188096
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201811045362.7A Active CN109271696B (zh) | 2018-09-07 | 2018-09-07 | 基于mpm的血液凝固模拟方法及系统 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN109271696B (zh) |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN114841047A (zh) * | 2022-04-25 | 2022-08-02 | 北京航空航天大学 | 基于物质点模型的血液凝固现象仿真方法 |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103699715A (zh) * | 2013-12-01 | 2014-04-02 | 北京航空航天大学 | 一种基于光滑粒子流体动力学和非线性有限元的流固耦合方法 |
| CN104318598A (zh) * | 2014-10-17 | 2015-01-28 | 中国科学技术大学 | 一种三维流固单向耦合的实现方法及系统 |
| CN107633123A (zh) * | 2017-09-13 | 2018-01-26 | 浙江工业大学 | 一种用于光滑粒子流体动力学模拟出血及处理加速的方法 |
| CN108461149A (zh) * | 2018-01-30 | 2018-08-28 | 中山大学 | 一种基于pbf的血液模拟方法 |
Family Cites Families (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| KR101267627B1 (ko) * | 2011-06-13 | 2013-05-27 | 한국과학기술원 | 멀티 레벨 소용돌이를 위한 ｓｐｈ 유체 시뮬레이션 방법, 시스템 및 이를 위한 기록 매체 |
| CN105956262B (zh) * | 2016-04-28 | 2019-08-09 | 清华大学 | 基于sph方法的多组分固体和流体模拟方法及系统 |
-
2018
- 2018-09-07 CN CN201811045362.7A patent/CN109271696B/zh active Active
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103699715A (zh) * | 2013-12-01 | 2014-04-02 | 北京航空航天大学 | 一种基于光滑粒子流体动力学和非线性有限元的流固耦合方法 |
| CN104318598A (zh) * | 2014-10-17 | 2015-01-28 | 中国科学技术大学 | 一种三维流固单向耦合的实现方法及系统 |
| CN107633123A (zh) * | 2017-09-13 | 2018-01-26 | 浙江工业大学 | 一种用于光滑粒子流体动力学模拟出血及处理加速的方法 |
| CN108461149A (zh) * | 2018-01-30 | 2018-08-28 | 中山大学 | 一种基于pbf的血液模拟方法 |
Non-Patent Citations (1)
| Title |
|---|
| 基于物质点法的水流运动三维数值仿真;王新宇 等;《科学技术与工程》;20180131;第293-297页 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN109271696A (zh) | 2019-01-25 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Solenthaler et al. | A unified particle model for fluid–solid interactions | |
| Chen et al. | Dynamics-aware numerical coarsening for fabrication design | |
| CN105069826B (zh) | 弹性物体变形运动的建模方法 | |
| Müller et al. | Detailed rigid body simulation with extended position based dynamics | |
| Xu et al. | Integrating viscoelastic mass spring dampers into position-based dynamics to simulate soft tissue deformation in real time | |
| Hoshyari et al. | Vibration-minimizing motion retargeting for robotic characters | |
| JP7431029B2 (ja) | 大規模環境用のマルチインスタンス型シミュレーション | |
| WO2002003173A2 (en) | Equal order method for fluid flow simulation | |
| JP2021072123A (ja) | スカラ輸送についてのガリレイ不変を強制する格子ボルツマンベースのスカラ輸送を使用して物理的プロセスをシミュレートするためのコンピュータシステム | |
| Civit‐Flores et al. | Robust treatment of degenerate elements in interactive corotational fem simulations | |
| Li et al. | Kinetic-based multiphase flow simulation | |
| CN109271696B (zh) | 基于mpm的血液凝固模拟方法及系统 | |
| Masterjohn et al. | Velocity level approximation of pressure field contact patches | |
| Lee et al. | Volumetric Object Modeling Using Internal Shape Preserving Constraint in Unity 3D. | |
| Ferrentino et al. | Finite element analysis-based soft robotic modeling: Simulating a soft actuator in sofa | |
| CN107526863A (zh) | 用于不可压缩瞬态和稳态navier‑stokes方程的最佳压力投影方法 | |
| Li et al. | Soft articulated characters in projective dynamics | |
| CN111898297B (zh) | 荷载和环境条件下跳台滑雪运动仿真模拟方法 | |
| Huang et al. | A survey on fast simulation of elastic objects | |
| Hamamoto et al. | Free-flight analysis of dragonfly hovering by fluid–structure interaction analysis based on an arbitrary Lagrangian–Eulerian method | |
| Volino et al. | Stop-and-go cloth draping | |
| CN108536954A (zh) | 一种基于交点间断伽辽金的高精度格子波尔兹曼方法 | |
| Alexiades et al. | Tin melting: Effect of grid size and scheme on the numerical solution. | |
| JP6832967B2 (ja) | リアルタイムアプリケーションにおける複数の接続されたボディのシミュレーション | |
| JP4006316B2 (ja) | 樹脂流動解析方法及びその装置 |
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 |