CN108491619B - 基于物理与非物理混合的复杂场景流固耦合高效模拟方法 - Google Patents
基于物理与非物理混合的复杂场景流固耦合高效模拟方法 Download PDFInfo
- Publication number
- CN108491619B CN108491619B CN201810226034.0A CN201810226034A CN108491619B CN 108491619 B CN108491619 B CN 108491619B CN 201810226034 A CN201810226034 A CN 201810226034A CN 108491619 B CN108491619 B CN 108491619B
- Authority
- CN
- China
- Prior art keywords
- particle
- rigid body
- fluid
- particles
- dfsph
- 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 112
- 239000007787 solid Substances 0.000 title claims abstract description 103
- 238000004088 simulation Methods 0.000 title claims abstract description 67
- 238000010168 coupling process Methods 0.000 title claims abstract description 50
- 230000008878 coupling Effects 0.000 title claims abstract description 37
- 238000005859 coupling reaction Methods 0.000 title claims abstract description 37
- 239000002245 particle Substances 0.000 claims abstract description 309
- 239000012530 fluid Substances 0.000 claims abstract description 108
- 239000012634 fragment Substances 0.000 claims description 38
- 230000000694 effects Effects 0.000 claims description 18
- 239000013598 vector Substances 0.000 claims description 17
- 230000009471 action Effects 0.000 claims description 12
- 230000008569 process Effects 0.000 claims description 11
- 230000006870 function Effects 0.000 claims description 10
- 230000008859 change Effects 0.000 claims description 9
- 239000000203 mixture Substances 0.000 claims description 8
- 239000000463 material Substances 0.000 claims description 6
- 239000000243 solution Substances 0.000 claims description 6
- 238000005070 sampling Methods 0.000 claims description 5
- 239000011259 mixed solution Substances 0.000 claims description 4
- 239000000126 substance Substances 0.000 claims description 4
- 230000008447 perception Effects 0.000 abstract description 3
- 230000011218 segmentation Effects 0.000 abstract description 3
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 11
- 230000003993 interaction Effects 0.000 description 8
- 238000010586 diagram Methods 0.000 description 6
- 230000008901 benefit Effects 0.000 description 4
- 238000004364 calculation method Methods 0.000 description 3
- 230000002457 bidirectional effect Effects 0.000 description 2
- 230000002146 bilateral effect Effects 0.000 description 1
- 238000005094 computer simulation Methods 0.000 description 1
- 230000005484 gravity Effects 0.000 description 1
- 230000003116 impacting effect Effects 0.000 description 1
- 238000013507 mapping Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000035515 penetration Effects 0.000 description 1
- 230000000704 physical effect Effects 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
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- General Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明公开了一种基于物理与非物理混合的复杂场景流固耦合高效模拟方法。其步骤为：1)不可压缩流体模拟，结合了欧拉方法和拉格朗日方法，提出了一种基于隐式粒子法的散度为零的光滑粒子流体动力学方法。2)流固耦合中的动力学问题求解，根据对象属性，将模拟对象的运动问题细分为三个子问题，分别使用不同方法对这三个子问题进行求解，实现一种多维度的分合计算框架。3)将断裂力学中的应变能密度概念和Voronoi空间分割相结合，实现了一种物理感知的破碎方法。本发明解决了在复杂流固耦合场景模拟中难以解决的模拟效率问题。对比以往流固耦合方法，本方法在相同系统资源下能够模拟细节更加丰富的流体场景，并且能够满足模拟流体冲击下固体发生破碎的需求。
Description
技术领域
本发明涉及基于粒子的流体模拟和流固双向耦合模拟领域，尤其涉及一种基于物理与非物理混合的复杂场景流固耦合高效模拟方法。
背景技术
图形学领域，流固耦合方法的使用往往伴随着流体模拟。从计算对象上可以分为固体对流体单向耦合(one-way solid-to-fluid coupling)，流体对固体单向耦合(one-way fluid-to-solid coupling)以及双向耦合(two-way coupling)。在固体对流体单向耦合方法中，固体按照设定的状态运动，流体不会对固体的运动路径产生影响。相对应的，流体对固体单向耦合方法中，通过对接触面上流体压力和扭转力的积分确定流体对固体的运动影响。然而这两种方法只能简单地作为固体或流体的边界条件，无法真实地模拟固体与流体交互。而双向耦合方法由于交互过程在物理上十分复杂，不存在解析解，不同的方法之间存在各种差异，适用范围小。
下面先介绍已有的流固双边耦合模拟方法：
1)任意拉格朗日-欧拉(Arbitrary Lagrangian-Eulerian)方法
任意拉格朗日-欧拉方法结合了传统的拉格朗日和欧拉方法，允许边界网格或面片顶点的任意运动，从而有效跟踪物质边界的运动，在内部网格上，根据定义参数求解网格，在独立于物质实体同时使得网格不发生严重的畸变。
2)内嵌边界(Immersed Boundary)方法
IB方法，通过在固体点所在的位置计算质量力来模拟流体域固体在边界的交互。质量力产生的运动约束被用来维持接触面上速度和应力的连续性。
以上的方法都不能满足复杂流固耦合模拟的需求，其原因在于模拟复杂场景时，场景内包含了大量固体与流体。
发明内容
本发明的目的在于解决现有现有图形学领域流固耦合模拟受到计算效率限制仅局限在较小的场景规模，而且无法在耦合过程中模拟固体破碎等效果的问题，提出了一种基于物理与非物理混合的复杂场景下流固耦合高效模拟方法。
本发明具体采用的技术方案如下：
基于物理与非物理混合的复杂场景流固耦合高效模拟方法，包括以下步骤：
1)将流体离散成2种粒子：DFSPH粒子和FLIP粒子；开始时，通过插值邻域FLIP粒子上个时间片最终的速度得到当前时间片DFSPH粒子的初始速度；然后，通过散度为零的光滑粒子流体动力学(Divergence-free SPH，简称DFSPH)方法精确求解粘性不可压缩流体方程，得到当前时间片DFSPH粒子最终的速度和位置；最后，通过插值邻域DFSPH粒子最终的速度求解得到当前时间片FLIP粒子最终的速度和位置；该步骤可以通过FLIP粒子丰富DFSPH方法模拟的流体细节的同时确保了模拟的效率。
2)将刚体质量平均分配到刚体表面的刚体粒子集合上；将刚体粒子集合标记为特殊的SPH粒子(赋予与原集合中其他粒子不同的标记)加入到1)中的DFSPH粒子集合中；通过累加刚体粒子集合的受力和转矩矢量求解得到刚体在当前时间片最终的位置、速度和旋转；通过直行树方法更新流体粒子位置和速度，防止穿透刚体表面；
3)根据时间步长，循环执行步骤1)和2)，得到每个时间片的刚体和流体粒子位置和速度。
上述方案中，若待模拟的场景中包含在预定时间破碎的固体，还可以采用一优选方式，即：当到达预定时间时，在执行步骤2)之前先执行步骤A)，步骤A)具体为：
计算该时刻固体表面的应变能分布；通过质心Voronoi方法最小化应变能分布计算得到所有可能的中心点集合P；最后输出最小化的中心点集合N作为破碎碎片的分布；将生成的碎片作为步骤2)中的刚体。
基于上述两个方案，步骤1)和2)可具体采用如下方式实现：
所述的步骤1)可以具体采用如下方法：
将流体离散成DFSPH粒子和FLIP粒子后，依次执行步骤(1)～(15)：
(1)对每个DFSPH粒子iD，搜索粒子iD所处位置半径为2h内的所有邻域FLIP粒子jF，以及粒子iD所处位置半径为h内除粒子iD之外的所有DFSPH粒子jD；
(2)由邻域的FLIP粒子上个时间片最终的速度vjF插值得到DFSPH粒子iD的初始速度viD：
式中WiDjF是一种经典的高斯形式核函数，它定义为：
式中：xiD、xjF分别是DFSPH粒子iD、FLIP粒子jF的坐标；d是空间维数，h是光滑核半径；
(3)根据邻域DFSPH粒子jD上个时间片的密度ρjD，通过核函数插值求解该DFSPH粒子iD当前时间片的密度ρiD：
式中：mjD是DFSPH粒子jD的质量；
(4)计算DFSPH粒子iD除了压力外的合力
(5)根据合力
计算中间速度
式中：Δt为模拟过程中设定的时间步长；
(6)计算系数αiD：
式中：
表示梯度；
(7)根据密度不变算子矫正粒子iD的中间速度
得到v′iD：
式中
ρ0设定的流体密度；
(8)根据矫正后的速度v′iD，更新DFSPH粒子iD的位置x′iD：
x′iD＝xiD+v′iDΔt
(9)更新粒子位置后，对所有DFSPH粒子iD，重新搜索邻域DFSPH粒子集合jD,计算粒子密度ρiD,计算常量系数αiD；
(10)根据散度为零算子，得到该时刻粒子iD的最终速度：
式中
为ρiD的物质导数；
(11)对所有FLIP粒子iF,搜索半径2h内的所有DFSPH粒子jD；
(12)重新加权插值得到FLIP粒子的速度
(13)根据FLIP方法的速度插补公式，得到
式中：viF为当前时间片FLIP粒子iF的初始速度；
(14)然后，通过正则化参数
将
和
混合求解得到当前时间片FLIP粒子的速度v″iF：
(15)根据求解得到的FLIP粒子的速度v″iF和DFSPH粒子的速度v″iD分别更新对应的位置：
x″iD＝x′iD+v″iDΔt
x″iF＝xiF+v″iFΔt。
所述的步骤2)可具体采用如下步骤实现：
(16)使用泊松盘方法对刚体表面进行采样，得到刚体粒子l；
(17)将刚体粒子l标记为特殊的DFSPH粒子，使用所述步骤(3)～(10)求解得到流体作用下的刚体粒子速度v″l；
(18)根据表面的刚体粒子速度集合U求解对应刚体的运动变化，求解过程如18.a)～18.e)：
18.a)对刚体粒子i，根据冲量定理，它的冲量Ji是动量的变化：
Ji＝FiΔt＝mi(v″l-vl)
式中：vl为刚体粒子l在上个时间片最终的速度；
18.b)计算耦合过程中刚体粒子i的受力Fi：
18.c)根据定义计算得到刚体粒子i转矩矢量τi：
τi＝(xi-xctr)×Fi
xi为刚体粒子的坐标，xctr是刚体质量中心的坐标；
18.d)将属于该刚体上刚体粒子受力和转矩矢量累加，得到施加在该刚体上的合力F以及扭力τ：
18.e)得到刚体所受合力F和扭力τ后，进行一次时间积分迭代计算合力与扭力对刚体产生的作用效果，刚体的线速度u″和角速度ω定义为：
ω＝I-1L
式中：uold为上个时间片刚体的线速度；M为刚体总质量，I是惯性张量，L是角动量；在每个时间步长里，L根据公式：
L＝Lold+τ
式中：Lold为上个时间片刚体的角动量；
惯性张量I定义为：
其中，
Ixy＝Iyx＝∑imixiyi，Ixz＝Izx＝∑imixizi，Iyz＝Izy＝∑imiyizi。xi,yi,zi分别对应刚体粒子i空间位置的X轴、Y轴、Z轴的坐标；
(19)使用直行树方法，根据刚体位置按照19.a)～19.d)约束流体粒子的运动：
19.a)对于流体粒子p，运动速度为vp，以及包含水平集数据的目标刚体，在流体粒子p所处位置x，假设局部水平集值为
那么其向外的法向量
目标刚体的局部速度为vs；
19.b)在当前的水平集值
基础上，预测在时间步长Δτ后粒子p的水平集值：
19.c)通过将公式中的
设为0，得到粒子p与目标刚体表面刚好发生接触的时间Δτ*：
19.d)计算得到Δτ*后，根据不同固体表面的弹性系数α对流体粒子的位置和速度进行补偿，避免流体粒子穿透表面渗入刚体内部的问题出现，从而实现在刚体约束下的流体粒子运动模拟。
所述的步骤A)可采用如下具体步骤：
A1)对于体积为Ω的固体，第k个中心点为pk，应变能密度函数W(x)，固体内的任意点x∈Ω，计算固体形变能量分布ED，k：
式中δ为预设的补偿项参数；
为基于距离的形变能量分布：
式中：dist(x,pk)2是x到中心点pk距离的平方；
A2)根据应变能量场来分割固体空间生产碎片，划分原则为碎片中心点位置使得破碎前的能量分布ED,k的和ED最小，对于区域中心点集合
需要满足：
P*＝argminED
A3)最后输出N＝min|P*|，N即是对原刚性固体满足ED<ε的空间划分，每个划分区域表示一个碎片，ε为预设的破碎阈值，当固体应变能量场大于ε发生破碎。
本发明的优点在于：
流固耦合问题涵盖了计算流体力学、数学、计算机等多门学科，往往包含大量的物理性质和数学运算，需要进行合理准确的抽象和建模，并通过计算机高效准确地模拟仿真。图形学领域，传统的流固耦合方法受到计算效率的限制，对于该类问题仅局限在较小的场景规模。
本发明结合了欧拉方法和拉格朗日方法，提出了一种基于隐式粒子法的散度为零SPH方法。第一步，使用DFSPH方法求解得出SPH粒子速度场分布；第二步，使用加权插值方法将SPH粒子速度场映射到FLIP隐式粒子上；最后，同时更新SPH粒子和FLIP粒子的位置。一方面通过使用粗粒度的SPH方法，确保流体模拟的不可压缩性，另一方使用细粒度的FLIP隐式粒子丰富流体模拟的细节效果。使参与模拟的粒子数达到百万的量级，从而使得复杂场景的同效模拟成为可能。
为进一步提高流固耦合效率，本发明提出了一种流体与固体在复杂场景下运动问题的解决方法。处理流体与固体交互：第一步，对刚体表面采样得到刚体粒子，将流体与固体接触面上的刚体粒子作为特殊的流体粒子参与到SPH速度场求解过程中，从而更新刚体粒子的速度场和流体粒子的速度场。第二步，通过累加刚体粒子的受力和扭矩矢量还原刚体的运动。通过直接运动树方法，解决流体粒子在固体边界出现的穿透问题。处理固体与固体的交互时，使用基于位置的动力学方法，求解碰撞约束，从而解决大量固体碎片的交互问题，实现了一种多维度的复杂场景下流固耦合分合计算框架。
为模拟流体冲击下的固体破碎效果，本发明提出了一种物理感知的高效破碎方法。将断裂力学中的应变能密度概念和Voronoi空间分割相结合，实现了一种物理感知的破碎方法。使得模拟结果在视觉上真实可知，更加符合人们对物理学上破碎效果的预期，从而高效实现了流体冲击下的固体破碎效果。
最后，本发明将上述方法整合，形成了一个复杂场景的流固耦合高效模拟框架，成功解决了传统流固耦合方法难以模拟固体在流体作用下发生破碎等问题。本发明中最后对高速水流冲破墙壁及洪水冲毁堤坝等场景进行了模拟，与传统的DFSPH方法相比，本发明在百万量级粒子参与的情况下，在保持了一定的场景真实感同时，大幅提高了模拟效率。
附图说明
图1是传统DFSPH方法与本发明步骤1)方法对比图；
图2是本发明步骤1)流体模拟结果图
图3是本发明步骤2)粒子维度下流固交互示意图；
图4是本发明结合步骤2)流固耦合场景模拟结果图；
图5传统破碎方法与本发明步骤3)破碎效果对比图；
图6是本发明实际应用：水流冲破墙壁模拟效果图；
图7是本发明实际应用：水流冲击大坝模拟效果图,a)～g)为逐步的结果展示。
具体实施方式
基于物理与非物理混合的复杂场景下流固耦合高效模拟方法包括以下三部分内容：
一、结合了欧拉方法和拉格朗日方法，提出了一种基于隐式粒子法(FLIP)的散度为零的光滑粒子流体动力学(DFSPH)方法，实现了高效的流体模拟方法。
二、流固耦合中的动力学问题求解，根据对象的属性，将模拟对象的运动问题细分为三个子问题，分别使用不同的方法对这三个子问题进行求解，实现了一种多维度的分合计算框架；
三、将断裂力学中的应变能密度概念和Voronoi空间分割相结合，实现了一种物理感知的破碎方法。
下面结合四个实施例来说明本发明中三部分内容所产生的技术效果。
实施例1
本实施例主要用于说明步骤1)流体模拟的优点。
本发明中流固耦合高效模拟时，首先将流体离散成2种粒子：DFSPH粒子和FLIP粒子；开始模拟时，通过插值邻域FLIP粒子上个时间片最终的速度(初次模拟时，上个时间片最终的速度可选取预设的初始速度)得到当前时间片DFSPH粒子的初始速度；然后，通过散度为零的光滑粒子流体动力学(Divergence-free SPH，简称DFSPH)方法精确求解粘性不可压缩流体方程，得到当前时间片DFSPH粒子最终的速度和位置；最后，通过插值邻域DFSPH粒子最终的速度求解得到当前时间片FLIP粒子最终的速度和位置。
本实施例中，将流体离散成DFSPH粒子和FLIP粒子后，依次执行具体步骤(1)～(15)：
(1)对每个DFSPH粒子iD，搜索粒子iD所处位置半径为2h内的所有邻域FLIP粒子jF，以及粒子iD所处位置半径为h内除粒子iD之外的所有DFSPH粒子jD；
(2)由邻域的FLIP粒子上个时间片最终的速度vjF插值得到DFSPH粒子iD的初始速度viD：
式中WiDjF是一种经典的高斯形式核函数，它定义为：
式中：xiD、xjF分别是DFSPH粒子iD、FLIP粒子jF的坐标；d是空间维数，h是光滑核半径；
(3)根据邻域DFSPH粒子jD上个时间片的密度ρjD，通过核函数插值求解该DFSPH粒子iD当前时间片的密度ρiD：
式中：mjD是DFSPH粒子jD的质量；
(4)计算DFSPH粒子iD除了压力外的合力
(5)根据合力
计算中间速度
式中：Δt为模拟过程中设定的时间步长；
(6)计算系数αiD：
式中：
表示梯度；
(7)根据密度不变算子矫正粒子iD的中间速度
得到v′iD：
式中
ρ0设定的流体密度；
(8)根据矫正后的速度v′iD，更新DFSPH粒子iD的位置x′iD：
x′iD＝xiD+v′iDΔt
(9)更新粒子位置后，对所有DFSPH粒子iD，重新搜索邻域DFSPH粒子集合jD,计算粒子密度ρiD,计算常量系数αiD；
(10)根据散度为零算子，得到该时刻粒子iD的最终速度：
式中
为ρiD的物质导数；
(11)对所有FLIP粒子iF,搜索半径2h内的所有DFSPH粒子jD；
(12)重新加权插值得到FLIP粒子的速度
(13)根据FLIP方法的速度插补公式，得到
式中：viF为当前时间片FLIP粒子iF的初始速度；
(14)然后，通过正则化参数
将
和
混合求解得到当前时间片FLIP粒子的速度v″iF：
(15)根据求解得到的FLIP粒子的速度v″iF和DFSPH粒子的速度v″iD分别更新对应的位置：
x″iD＝x′iD+v″iDΔt
x″iF＝xiF+v″iFΔt。
图1是传统DFSPH方法与本发明步骤1)方法的结果对比图；图2是本发明步骤1)流体模拟结果图。图1在(a)的模拟中，使用了约20万SPH粒子和约80万FLIP粒子，通过本发明方法1)求解流体粒子运动；在(b)的模拟中，只使用了约20万SPH粒子，通过传统的DFSPH方法求解流体粒子运动；在(c)的模拟中，使用了约100万SPH粒子，通过传统的DFSPH方法求解流体粒子运动，可以发现(a)与(b)，都使用了20万SPH粒子，而(a)使用本文的FLIP-DFSPH方法通过80万FLIP粒子丰富了流体的细节，最终得到的效果远好于(b)。可以发现，(a)的最终效果与使用了100万粒子的DFSPH模拟结果(c)在视觉上相差不多。而在运行环境，CPU：i7-4790K、GPU：NVidia GTX 960、内存8GB下，(a)平均4.2s秒得到一帧，而传统方法(c)平均27.8秒得到一帧足以证明本发明在模拟较多流体粒子时效率高于传统的DFSPH方法，并且本发明方法可以通过设置FLIP和DFSPH粒子的数目，更加灵活地在计算精度和模拟效果之间调节，更适合应用于模拟细节丰富的复杂流固耦合场景。
实施例2
本实施例主要用于说明步骤2)流固耦合中的动力学问题求解的优点。
本发明中，流固耦合时，先将刚体质量平均分配到刚体表面的刚体粒子集合上；将刚体粒子集合标记为特殊的SPH粒子加入到1)中的DFSPH粒子集合中；通过累加刚体粒子集合的受力和转矩矢量求解得到刚体在当前时间片最终的位置、速度和旋转；通过直行树方法更新流体粒子位置和速度，防止穿透刚体表面。该步骤的具体实现过程为：
(16)使用泊松盘方法对刚体表面进行采样，得到刚体粒子l；
(17)将刚体粒子l标记为特殊的DFSPH粒子，使用所述步骤(3)～(10)相同的方法求解得到流体作用下的刚体粒子速度v″l；
(18)根据表面的刚体粒子速度集合U求解对应刚体的运动变化，求解过程如18.a)～18.e)：
18.a)对刚体粒子i，根据冲量定理，它的冲量Ji是动量的变化：
Ji＝FiΔt＝mi(v″l-vl)
式中：vl为刚体粒子l在上个时间片最终的速度；
18.b)计算耦合过程中刚体粒子i的受力Fi：
18.c)根据定义计算得到刚体粒子i转矩矢量τi：
τi＝(xi-xctr)×Fi
xi为刚体粒子的坐标，xctr是刚体质量中心的坐标；
18.d)将属于该刚体上刚体粒子受力和转矩矢量累加，得到施加在该刚体上的合力F以及扭力τ：
18.e)得到刚体所受合力F和扭力τ后，进行一次时间积分迭代计算合力与扭力对刚体产生的作用效果，刚体的线速度u″和角速度ω定义为：
ω＝I-1L
式中：uold为上个时间片刚体的线速度；M为刚体总质量，I是惯性张量，L是角动量；在每个时间步长里，L根据公式：
L＝Lold+τ
式中：Lold为上个时间片刚体的角动量；
惯性张量I定义为：
其中，
Ixy＝Iyx＝∑imixiyi，Ixz＝Izx＝∑imixizi，Iyz＝Izy＝∑imiyizi。xi,yi,zi分别对应刚体粒子i空间位置的X轴、Y轴、Z轴的坐标；
(19)使用直行树方法，根据刚体位置按照19.a)～19.d)约束流体粒子的运动：
19.a)对于流体粒子p，运动速度为vp，以及包含水平集数据的目标刚体，在流体粒子p所处位置x，假设局部水平集值为
那么其向外的法向量
目标刚体的局部速度为vs；
19.b)在当前的水平集值
基础上，预测在时间步长Δτ后粒子p的水平集值：
19.c)通过将公式中的
设为0，得到粒子p与目标刚体表面刚好发生接触的时间Δτ*：
19.d)计算得到Δτ*后，根据不同固体表面的弹性系数α对流体粒子的位置和速度进行补偿，避免流体粒子穿透表面渗入刚体内部的问题出现，从而实现在刚体约束下的流体粒子运动模拟。
图3是执行上述步骤2)后粒子维度下流固交互示意图；图4给出了本发明方法的应用实例，是本发明结合步骤2)流固耦合场景模拟结果图。如图4中的(a)～(d)所示，随着时间的迭代，刚体在流体作用下发生了位移和旋转，成功模拟了刚体在水体的漂浮、沉没场景。同时可以观察到流体在刚体的作用下也发生了波动起伏。从而证明本发明的方法能够成功地应用于复杂流固耦合场景中流体和刚体的运动模拟。
实施例3
本实施例主要用于说明步骤A)所能产生的技术效果。该步骤仅在模拟场景中具有需要破碎的固体(例如被水冲击后破碎的大坝等)时方才执行，若场景中不存在此类固体时，无需执行。
模拟场景中包含的固体需在预定时间破碎，因此当到达预定时间时，在执行步骤2)之前先执行步骤A)，步骤A)具体为：
计算该时刻固体表面的应变能分布；通过质心Voronoi方法最小化应变能分布计算得到所有可能的中心点集合P；最后输出最小化的中心点集合N作为破碎碎片的分布；将生成的碎片作为步骤2)中的刚体。该步骤的具体实现方式如下：
A1)对于体积为Ω的固体，第k个中心点为pk，应变能密度函数W(x)，固体内的任意点x∈Ω，计算固体形变能量分布ED,k：
式中δ为预设的补偿项参数；
为基于距离的形变能量分布：
式中：dist(x,pk)2是x到中心点pk距离的平方；
A2)根据应变能量场来分割固体空间生产碎片，划分原则为碎片中心点位置使得破碎前的能量分布ED,k的和ED最小，对于区域中心点集合
需要满足：
P*＝argminED
A3)最后输出N＝min|P*|，N即是对原刚性固体满足ED<ε的空间划分，每个划分区域表示一个碎片，ε为预设的破碎阈值，当固体应变能量场大于ε发生破碎。
图5比较了上述步骤3)的破碎方法和传统随机Voronoi破碎方法的结果，其中：(a)为本发明破碎方法，(b)为传统随机Voronoi破碎方法。可以看出本发明步骤3)生成的破碎碎片与碰撞点相关，而传统随机Voronoi破碎方法生成的碎片与真实物理碰撞没有关系。所以本发明步骤3)的破碎模拟方法相比于传统的随机Voronoi破碎方法在视觉上更加真实，能够与物体实际受力相结合。
实施例4
本实施例以水库中洪水水流作为流体，水库大坝作为固体，其主要实现方式为：为通过本发明步骤1)求解流体粒子的运动，通过本发明步骤2)计算场景内流体粒子和刚体粒子因为相互作用而产生的运动变化，根据本发明步骤A)模拟固体在流体粒子作用下发生破碎，破碎后的碎片被作为刚体加入系统中参与到本发明的步骤2)模拟中。
本实施例中，基于物理与非物理混合的复杂场景流固耦合高效模拟方法，包括以下步骤：
1)将流体离散成2种粒子：DFSPH粒子和FLIP粒子；开始时，通过插值邻域FLIP粒子上个时间片最终的速度得到当前时间片DFSPH粒子的初始速度；然后，通过散度为零的光滑粒子流体动力学(Divergence-free SPH，简称DFSPH)方法精确求解粘性不可压缩流体方程，得到当前时间片DFSPH粒子最终的速度和位置；最后，通过插值邻域DFSPH粒子最终的速度求解得到当前时间片FLIP粒子最终的速度和位置；该步骤可以通过FLIP粒子丰富DFSPH方法模拟的流体细节的同时确保了模拟的效率。
再具体来说，在将流体离散成DFSPH粒子和FLIP粒子后，依次执行步骤(1)～(15)：
(1)对每个DFSPH粒子iD，搜索粒子iD所处位置半径为2h内的所有邻域FLIP粒子jF，以及粒子iD所处位置半径为h内除粒子iD之外的所有DFSPH粒子jD；
(2)由邻域的FLIP粒子上个时间片最终的速度vjF插值得到DFSPH粒子iD的初始速度viD：
式中WiDjF是高斯形式核函数，它定义为：
式中：xiD、xjF分别是DFSPH粒子iD、FLIP粒子jF的坐标；d是空间维数，h是光滑核半径；
(3)根据邻域DFSPH粒子jD上个时间片的密度ρjD，通过核函数插值求解该DFSPH粒子iD当前时间片的密度ρiD：
式中：mjD是DFSPH粒子jD的质量；
(4)计算DFSPH粒子iD除了压力外的合力
本实施例中只考虑重力；
(5)根据合力
计算中间速度
式中：Δt为模拟过程中设定的时间步长；
(6)计算系数αiD：
式中：
表示梯度；
(7)根据密度不变算子矫正粒子iD的中间速度
得到v′iD：
式中
ρ0设定的流体密度；
(8)根据矫正后的速度v′iD，更新DFSPH粒子iD的位置x′iD：
x′iD＝xiD+v′iDΔt
(9)更新粒子位置后，对所有DFSPH粒子iD，重新搜索邻域DFSPH粒子集合jD,计算粒子密度ρiD,计算常量系数αiD；
(10)根据散度为零算子，得到该时刻粒子iD的最终速度：
式中
为ρiD的物质导数；
(11)对所有FLIP粒子iF,搜索半径2h内的所有DFSPH粒子jD；
(12)重新加权插值得到FLIP粒子的速度
(13)根据FLIP方法的速度插补公式，得到
式中：viF为当前时间片FLIP粒子iF的初始速度；
(14)然后，通过正则化参数
将
和
混合求解得到当前时间片FLIP粒子的速度v″iF：
(15)根据求解得到的FLIP粒子的速度v″iF和DFSPH粒子的速度v″iD分别更新对应的位置：
x″iD＝x′iD+v″iDΔt
x″iF＝xiF+v″iFΔt。
需注意的是，在初次执行步骤1)时，由于不存在上个时间片的各参数，因此需预设初始参数。
经过上述步骤(1)～(15)，即完成了当前时间片的流体模拟。由于模拟场景中包含大坝，而大坝在水流冲击一定时间后会破碎，因此当到达预定破碎时间时，需在执行步骤2)之前先执行步骤A)：计算该时刻大坝固体表面的应变能分布；通过质心Voronoi方法最小化应变能分布计算得到所有可能的中心点集合P；最后输出最小化的中心点集合N作为破碎碎片的分布；将生成的碎片作为步骤2)中的刚体。
本实施例中，该步骤的碎片的具体生成方法如下：
A1)对于体积为Ω的固体，第k个中心点为pk，应变能密度函数W(x)，固体内的任意点x∈Ω，计算固体形变能量分布ED,k：
式中δ为预设的补偿项参数；
为基于距离的形变能量分布：
式中：dist(x,pk)2是x到中心点pk距离的平方；
A2)根据应变能量场来分割固体空间生产碎片，划分原则为碎片中心点位置使得破碎前的能量分布ED,k的和ED最小，对于区域中心点集合
需要满足：
P*＝argminED
A3)最后输出N＝min|P*|，N即是对原刚性固体满足ED<ε的空间划分，每个划分区域表示一个碎片，ε为预设的破碎阈值，当固体应变能量场大于ε发生破碎。
由此，即可得到大坝破碎后的碎片情况，当大坝破碎后其会与洪水水流产生耦合，因此需要进行流固耦合。
2)将前一步中刚体(即大坝碎片)质量平均分配到刚体表面的刚体粒子集合上；将刚体粒子集合标记为特殊的SPH粒子加入到1)中的DFSPH粒子集合中；通过累加刚体粒子集合的受力和转矩矢量求解得到刚体在当前时间片最终的位置、速度和旋转；通过直行树方法更新流体粒子位置和速度，防止穿透刚体表面。其具体实现步骤为：
(16)使用泊松盘方法对刚体表面进行采样，得到刚体粒子l；
(17)将刚体粒子l标记为特殊的DFSPH粒子，使用所述步骤(3)～(10)求解得到流体作用下的刚体粒子速度v″l；
(18)根据表面的刚体粒子速度集合U求解对应刚体的运动变化，求解过程如18.a)～18.e)：
18.a)对刚体粒子i，根据冲量定理，它的冲量Ji是动量的变化：
Ji＝FiΔt＝mi(v″l-l)
式中：vl为刚体粒子l在上个时间片最终的速度；
18.b)计算耦合过程中刚体粒子i的受力Fi：
18.c)根据定义计算得到刚体粒子i转矩矢量τi：
τi＝(xi-xctr)×Fi
xi为刚体粒子的坐标，xctr是刚体质量中心的坐标；
18.d)将属于该刚体上刚体粒子受力和转矩矢量累加，得到施加在该刚体上的合力F以及扭力τ：
18.e)得到刚体所受合力F和扭力τ后，进行一次时间积分迭代计算合力与扭力对刚体产生的作用效果，刚体的线速度u″和角速度ω定义为：
ω＝I-1L
式中：uold为上个时间片刚体的线速度；M为刚体总质量，I是惯性张量，L是角动量；在每个时间步长里，L根据公式：
L＝Lold+τ
式中：Lold为上个时间片刚体的角动量；
惯性张量I定义为：
其中，
Ixy＝Iyx＝∑imixiyi，Ixz＝Izx＝∑imixizi，Iyz＝Izy＝∑imiyizi。xi,yi,zi分别对应刚体粒子i空间位置的X轴、Y轴、Z轴的坐标；
(19)使用直行树方法，根据刚体位置按照19.a)～19.d)约束流体粒子的运动：
19.a)对于流体粒子p，运动速度为vp，以及包含水平集数据的目标刚体，在流体粒子p所处位置x，假设局部水平集值为
那么其向外的法向量
目标刚体的局部速度为vs；
19.b)在当前的水平集值
基础上，预测在时间步长Δτ后粒子p的水平集值：
19.c)通过将公式中的
设为0，得到粒子p与目标刚体表面刚好发生接触的时间Δτ*：
19.d)计算得到Δτ*后，根据不同固体表面的弹性系数α对流体粒子的位置和速度进行补偿，避免流体粒子穿透表面渗入刚体内部的问题出现，从而实现在刚体约束下的流体粒子运动模拟。
经过步骤1)和2)，即可得到模拟过程中第一个时间片的刚体和流体粒子位置和速度。
3)根据模拟中设置的时间步长，即可循环迭代执行步骤1)和2)，得到每个时间片的刚体和流体粒子位置和速度。
通过以上步骤，图6和图7给出了本发明实际应用模拟效果图：图6中步骤(a)～(d)模拟了水流冲击墙壁，墙壁破碎后碎片和流体的运动情况；图7中步骤(a)～(g)模拟了洪水冲击大坝，大坝崩溃后，洪水与大坝碎片的运动情况。对比以往流固耦合方法，本方法在相同系统资源下能够模拟细节更加丰富的流体场景，并且能够满足模拟流体冲击下固体发生破碎的需求。
以上列举的仅是本发明的具体实施例。显然，本发明不限于以上实施例，还可以有许多变形。本领域的普通技术人员能从本发明公开的内容直接导出或联想到的所有变形，均应认为是本发明的保护范围。
Claims (4)
1.一种基于物理与非物理混合的复杂场景流固耦合高效模拟方法，其特征在于，包括以下步骤：
1)将流体离散成2种粒子：DFSPH粒子和FLIP粒子；开始时，通过插值邻域FLIP粒子上个时间片最终的速度得到当前时间片DFSPH粒子的初始速度；然后，通过散度为零的光滑粒子流体动力学方法精确求解粘性不可压缩流体方程，得到当前时间片DFSPH粒子最终的速度和位置；最后，通过插值邻域DFSPH粒子最终的速度求解得到当前时间片FLIP粒子最终的速度和位置；
2)将刚体质量平均分配到刚体表面的刚体粒子集合上；将刚体粒子集合标记为特殊的SPH粒子加入到1)中的DFSPH粒子集合中；通过累加刚体粒子集合的受力和转矩矢量求解得到刚体在当前时间片最终的位置、速度和旋转；通过直行树方法更新流体粒子位置和速度，防止穿透刚体表面；
3)根据时间步长，循环执行步骤1)和2)，得到每个时间片的刚体和流体粒子位置和速度；
所述的步骤1)具体为：
将流体离散成DFSPH粒子和FLIP粒子后，依次执行步骤(1)～(15)：
(1)对每个DFSPH粒子iD，搜索粒子iD所处位置半径为2h内的所有邻域FLIP粒子jF，以及粒子iD所处位置半径为h内除粒子iD之外的所有DFSPH粒子jD；
(2)由邻域的FLIP粒子上个时间片最终的速度vjF插值得到DFSPH粒子iD的初始速度viD：
式中WiDjF是高斯形式核函数，它定义为：
式中：xiD、xjF分别是DFSPH粒子iD、FLIP粒子jF的坐标；d是空间维数，h是光滑核半径；
(3)根据邻域DFSPH粒子jD上个时间片的密度ρjD，通过核函数插值求解该DFSPH粒子iD当前时间片的密度ρiD：
式中：mjD是DFSPH粒子jD的质量；
(4)计算DFSPH粒子iD除了压力外的合力
(5)根据合力
计算中间速度
式中：Δt为模拟过程中设定的时间步长；
(6)计算常量系数αiD：
式中：
表示梯度；
(7)根据密度不变算子矫正粒子iD的中间速度
得到v′iD：
式中
ρ0设定的流体密度；
(8)根据矫正后的速度v′iD，更新DFSPH粒子iD的位置x′iD：
x′iD＝xiD+v′iDΔt
(9)更新粒子位置后，对所有DFSPH粒子iD，重新搜索邻域DFSPH粒子集合jD,计算粒子密度ρiD,计算常量系数αiD；
(10)根据散度为零算子，得到当前时刻粒子iD的最终速度：
式中
为ρiD的物质导数；
(11)对所有FLIP粒子iF,搜索半径2h内的所有DFSPH粒子jD；
(12)重新加权插值得到FLIP粒子的速度
(13)根据FLIP方法的速度插补公式，得到
式中：viF为当前时间片FLIP粒子iF的初始速度；
(14)然后，通过正则化参数
将
和
混合求解得到当前时间片FLIP粒子的速度v″iF：
(15)根据求解得到的FLIP粒子的速度v″iF和DFSPH粒子的速度v″iD分别更新对应的位置：
x″iD＝x′iD+v″iDΔt
x″iF＝xiF+v″iFΔt。
2.如权利要求1所述的基于物理与非物理混合的复杂场景流固耦合高效模拟方法，其特征在于，模拟场景中包含在预定时间破碎的固体，当到达预定时间时，在执行步骤2)之前先执行步骤A)，步骤A)具体为：
计算当前时刻固体表面的应变能分布；通过质心Voronoi方法最小化应变能分布计算得到所有可能的中心点集合P；最后输出最小化的中心点集合N作为破碎碎片的分布；将生成的碎片作为步骤2)中的刚体。
3.根据权利要求1所述的基于物理与非物理混合的复杂场景流固耦合高效模拟方法，其特征在于所述的步骤2)具体为：
(16)使用泊松盘方法对刚体表面进行采样，得到刚体粒子l；
(17)将刚体粒子l标记为特殊的DFSPH粒子，使用所述步骤(3)～(10)求解得到流体作用下的刚体粒子速度v″l；
(18)根据表面的刚体粒子速度集合U求解对应刚体的运动变化，求解过程如18.a)～18.e)：
18.a)对刚体粒子i，根据冲量定理，它的冲量Ji是动量的变化：
Ji＝FiΔt＝mi(v″l-vl)
式中：vl为刚体粒子l在上个时间片最终的速度；
18.b)计算耦合过程中刚体粒子i的受力Fi：
18.c)根据定义计算得到刚体粒子i转矩矢量τi：
τi＝(xi-xctr)×Fi
xi为刚体粒子的坐标，xctr是刚体质量中心的坐标；
18.d)将属于该刚体上刚体粒子受力和转矩矢量累加，得到施加在该刚体上的合力F以及扭力τ：
18.e)得到刚体所受合力F和扭力τ后，进行一次时间积分迭代计算合力与扭力对刚体产生的作用效果，刚体的线速度u″和角速度ω定义为：
ω＝I-1L
式中：uold为上个时间片刚体的线速度；M为刚体总质量，I是惯性张量，L是角动量；在每个时间步长里，L根据公式：
L＝Lold+τ
式中：Lold为上个时间片刚体的角动量；
惯性张量I定义为：
其中，
Ixy＝Iyx＝∑imixiyi，Ixz＝Izx＝∑imixizi，Iyz＝Izy＝∑imiyizi，xi,yi,zi分别对应刚体粒子i空间位置的X轴、Y轴、Z轴的坐标；
(19)使用直行树方法，根据刚体位置按照19.a)～19.d)约束流体粒子的运动：
19.a)对于流体粒子p，运动速度为vp，以及包含水平集数据的目标刚体，在流体粒子p所处位置x，假设局部水平集值为
那么其向外的法向量
目标刚体的局部速度为vs；
19.b)在当前的水平集值
基础上，预测在时间步长Δτ后粒子p的水平集值：
19.c)通过将公式中的
设为0，得到粒子p与目标刚体表面刚好发生接触的时间Δτ*：
19.d)计算得到Δτ*后，根据不同固体表面的弹性系数α对流体粒子的位置和速度进行补偿，避免流体粒子穿透表面渗入刚体内部的问题出现，从而实现在刚体约束下的流体粒子运动模拟。
4.根据权利要求3所述的一种基于物理与非物理混合的复杂场景流固耦合高效模拟方法，其特征在于所述的步骤A)具体为：
A1)对于体积为Ω的固体，第k个中心点为pk，应变能密度函数W(x)，固体内的任意点x∈Ω，计算固体形变能量分布ED,k：
式中δ为预设的补偿项参数；
为基于距离的形变能量分布：
式中：dist(x,pk)2是x到中心点pk距离的平方；
A2)根据应变能量场来分割固体空间生产碎片，划分原则为碎片中心点位置使得破碎前的能量分布ED,k的和ED最小，对于区域中心点集合
需要满足：
P*＝argmin ED
A3)最后输出N＝min|P*|，N即是对原刚性固体满足ED<ε的空间划分，每个划分区域表示一个碎片，ε为预设的破碎阈值，当固体应变能量场大于ε发生破碎。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810226034.0A CN108491619B (zh) | 2018-03-19 | 2018-03-19 | 基于物理与非物理混合的复杂场景流固耦合高效模拟方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810226034.0A CN108491619B (zh) | 2018-03-19 | 2018-03-19 | 基于物理与非物理混合的复杂场景流固耦合高效模拟方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN108491619A CN108491619A (zh) | 2018-09-04 |
| CN108491619B true CN108491619B (zh) | 2020-05-26 |
Family
ID=63318432
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201810226034.0A Expired - Fee Related CN108491619B (zh) | 2018-03-19 | 2018-03-19 | 基于物理与非物理混合的复杂场景流固耦合高效模拟方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN108491619B (zh) |
Families Citing this family (9)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109902376B (zh) * | 2019-02-25 | 2021-01-15 | 北京理工大学 | 一种基于连续介质力学的流固耦合高精度数值模拟方法 |
| CN109992830B (zh) * | 2019-02-26 | 2020-11-13 | 浙江大学 | 基于物质点方法的山体滑坡灾害场景模拟方法 |
| CN111241742B (zh) * | 2019-12-27 | 2021-11-19 | 西安交通大学 | 一种多相流计算方法 |
| CN111695282A (zh) * | 2020-06-08 | 2020-09-22 | 河海大学 | 一种基于流固耦合模拟的液舱晃荡预测和控制分析方法 |
| CN111814279B (zh) * | 2020-09-14 | 2020-12-11 | 四川轻化工大学 | 一种基于sph的齿轮齿条动态啮合及传动过程分析方法 |
| CN112307668A (zh) * | 2020-11-02 | 2021-02-02 | 浙江工业大学 | 一种基于粒子的粘液类效果模拟方法 |
| CN113192567B (zh) * | 2021-04-30 | 2022-12-09 | 西安交通大学 | 一种核反应堆板状燃料熔化流固耦合无网格分析方法 |
| CN113283066B (zh) * | 2021-05-14 | 2022-09-09 | 北京大学 | 带表面张力的固液强耦合模拟方法、装置、设备及介质 |
| CN116595849B (zh) * | 2023-05-19 | 2024-01-19 | 长安大学 | 一种金属结构冲击损伤破坏模型的构建方法及装置 |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104143027A (zh) * | 2014-08-01 | 2014-11-12 | 北京理工大学 | 一种基于sph算法的流体热运动仿真系统 |
| CN104318598A (zh) * | 2014-10-17 | 2015-01-28 | 中国科学技术大学 | 一种三维流固单向耦合的实现方法及系统 |
| CN107203652A (zh) * | 2017-04-01 | 2017-09-26 | 浙江科技学院（浙江中德科技促进中心） | 地震液化中地下结构上浮离心机试验的精细化模拟方法 |
-
2018
- 2018-03-19 CN CN201810226034.0A patent/CN108491619B/zh not_active Expired - Fee Related
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104143027A (zh) * | 2014-08-01 | 2014-11-12 | 北京理工大学 | 一种基于sph算法的流体热运动仿真系统 |
| CN104318598A (zh) * | 2014-10-17 | 2015-01-28 | 中国科学技术大学 | 一种三维流固单向耦合的实现方法及系统 |
| CN107203652A (zh) * | 2017-04-01 | 2017-09-26 | 浙江科技学院（浙江中德科技促进中心） | 地震液化中地下结构上浮离心机试验的精细化模拟方法 |
Non-Patent Citations (3)
| Title |
|---|
| Divergence-Free SPH for Incompressible and Viscous Fluids;Jan Bender,etc;《IEEE TRANSCATION ON VISUALIZAITON AND COMPUTER GRAPHICS》;20160608;第23卷(第3期);第1193页至1206页 * |
| GPU中的流体场景实时模拟算法;陈曦等;《计算机辅助设计与图像学学报》;20100330;第22卷(第3期);第396页至405页 * |
| Physical Solid Fracture Simulation Based on Random Voronoi Tessallation;Shaoxiong Zhang,etc;《Advances in Computer Science Research》;20161130;第52卷;第209页至213页 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN108491619A (zh) | 2018-09-04 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN108491619B (zh) | 基于物理与非物理混合的复杂场景流固耦合高效模拟方法 | |
| Jiang et al. | The material point method for simulating continuum materials | |
| Chentanez et al. | Simultaneous coupling of fluids and deformable bodies | |
| Boyd et al. | MultiFLIP for energetic two-phase fluid simulation | |
| Yokoi | A practical numerical framework for free surface flows based on CLSVOF method, multi-moment methods and density-scaled CSF model: Numerical simulations of droplet splashing | |
| US7505883B2 (en) | Computer simulation of body dynamics including a solver that solves in linear time for a set of constraints | |
| Gerszewski et al. | A point-based method for animating elastoplastic solids | |
| Li et al. | A fully 3D simulation of fluid-structure interaction with dynamic wetting and contact angle hysteresis | |
| CN110992456B (zh) | 一种基于位置动力学的雪崩模拟方法 | |
| Mihalef et al. | Simulation of two‐phase flow with sub‐scale droplet and bubble effects | |
| CN106650064A (zh) | 一种基于粒子模型的凝结现象仿真方法 | |
| Barreiro et al. | Conformation constraints for efficient viscoelastic fluid simulation | |
| Li et al. | A finite difference method with subsampling for immersed boundary simulations of the capsule dynamics with viscoelastic membranes | |
| US8392154B2 (en) | System and method for real-time cloth simulation | |
| Ma et al. | Learning neural constitutive laws from motion observations for generalizable pde dynamics | |
| Chang et al. | A particle-based method for viscoelastic fluids animation | |
| Zhang et al. | Simulation system for collisions and two-way coupling of non-Newtonian fluids and solids | |
| CN105183965B (zh) | 用于预测雾化过程的大涡模拟方法 | |
| Gao et al. | Accelerating liquid simulation with an improved data‐driven method | |
| Jones et al. | Physically-based droplet interaction | |
| US20180089878A1 (en) | Inertial damping for enhanced simulation of elastic bodies | |
| CN115393545B (zh) | 基于可形变网格的碰撞处理方法、系统、设备及介质 | |
| Si et al. | Thin-feature-aware transport-velocity formulation for SPH-based liquid animation | |
| Lazo et al. | Real-time physical engine for floating objects with two-way fluid-structure coupling | |
| Shao et al. | Unified particle-based coupling algorithm for stable soft tissue-blood interaction |
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 |
Granted publication date: 20200526 Termination date: 20210319 |
|
| CF01 | Termination of patent right due to non-payment of annual fee |