面向非牛顿流体仿真的边界处理方法

1、面向非牛顿流体仿真的边界处理方法摘 要：针对现有的光滑粒子流体动力学(SPH)非牛顿流体仿真方法的边界条件不合理问题，提出一种适用于非牛 顿流体仿真的边界处理方法.首先使用Cross模型实现牛顿流体和非牛顿流体的统一建模；然后通过对固体边界粒 子进行加权计算，为边界处的流体粒子设定边界力的法向量；最后采用Coulomb摩擦力边界条件对边界处流体粒子 的速度进行迭代更新，并将其融合到预测-校正算法框架中.实验结果表明，文中方法能够根据用户的需要调节流固 边界处的摩擦系数；与自由滑移边界下的仿真相比，该方法下非牛顿流体黏度随速度的降低而增大，并表现出固体 形态的物理特性.关键词：基于物理的动画；非

2、牛顿流体；光滑粒子流体动力学；预测-校正法；边界处理Boundary Handling for N on-N ewtonian Fluid SimulationAbstract: To solve the problem of unreasonable boundary conditions of the existing SPH non-Newtonian fluid simulation, a boundary handling method is proposed. First, a unified modeling of Newtonian and non-Newtonian flui

3、ds was implemented by Cross model. Then, by calculating the weighted sum of solid boundary particles, a normal vector of the boundary force was set for the fluid particles at the boundary. Finally, the velocity of fluid particles near the boundary was iteratively updated using Coulomb friction bound

4、ary conditions, which can be integrated into the predictive-corrective algorithm framework. Experiments show that the friction coefficient of the fluid-solid boundary can be adjusted according to the users experimental needs. Compared to simulation under free-slip boundary conditions, the non-Newton

5、ian fluid simulated by our method is able to increase in viscosity with the decrease of velocity and exhibit physical properties of solids.Key words: physically-based animation; non-Newtonian fluid; smoothed particle hydrodynamics; predictive-corrective method; boundary handling近年来，对非牛顿流体的仿真得到研究者的 关

6、注，在泥石流灾害仿真、影视特效、娱乐游戏等 领域具有广泛的应用价值.非牛顿流体的物理本 质非常复杂，它们具有独特的材质属性，同时具备 流体和固体的物理特征.以剪切变稀型非牛顿流 体为例，流体的黏度随速度的增加而降低，当流体 在较高的速度下流动时，黏度较低，其物理特性偏 向于牛顿流体，会出现飞溅、浪花等现象；当流体 受到较大的阻力而流速降低时，流体的黏度逐渐 增加，使其逐渐表现出固体的物理特性，流体粒子 会凝聚在一起，保持稳定.然而，目前已有的非牛顿流体仿真方法11中缺 失对流固边界的有效处理.本文实验证明，采用自 由滑移边界条件处理流固边界，会使得非牛顿流 体在与固体交互的过程中无法从固体边界

7、处获得 足够的阻力和反作用力，进而使得非牛顿流体黏 度随速度的降低而增大、表现出固体属性的物理特 征，无法通过仿真得到验证和体现.为了解决这一问题，本文对面向非牛顿流体 的边界处理进行研究.用光滑粒子流体动力学 （smoothed particle hydrodynamics，SPH）2处理非 牛顿流体边界问题存在许多困难：（1）流体边界重 构过程中边界法向量不精确.与网格法相比，由流 体粒子组成的流体边界难以找到精确的边界法向 量，而边界法向量的计算偏差会导致流体粒子到 真实边界的距离计算错误，进而影响SPH加权计 算中密度和受力的准确性.（2）边界摩擦力的计算. 采用恰当的摩擦力边界模型，

8、才能使非牛顿流体 在仿真过程中充分体现其物理特性.（3）将摩擦力 边界模型应用到已有的预测-校正算法中.为了使 非牛顿流体仿真能够在更大的时间步长下稳定运 行，已有的非牛顿流体仿真采用预测-校正方法. 为了维持原仿真算法的高效稳定性，需要将边界 处理方法结合到原有的预测-校正算法框架下.本文提出一种面向非牛顿流体的边界处理方 法.首先实现基于预测-校正的非牛顿流体仿真; 然后通过对固体边界粒子进行加权计算，对流体 粒子的边界法向量进行精确求解；最后采用 Coulomb摩擦力边界条件进行边界处的受力计算 和速度更新，并将边界处理方法融合到已有的预 测-校正算法框架中.1相关工作用拉格朗日粒子法进

9、行动画仿真一直都是计 算机图形学的研究热点A%在进行非牛顿流体和 高黏度牛顿流体的仿真过程中，研究人员提出很 多方法来处理边界问题.Gerszewski等通过找到 最接近相邻粒子在单个时间步长上的运动的仿射 变换，为每个粒子计算变形梯度，进而实现非牛顿 流体的本构模型.然而，当流体没有足够的邻居粒 子时梯度估计精确度下降，算法稳定性变差，这一 问题在缺失流体邻居粒子的流固边界格外明显. Paiva等81使用广义牛顿流体模型将非牛顿流体建 模为高黏度的固体形态，其流体特性只有在高温 时才能有所体现；该方法采用不滑移边界条件, 非牛顿流体在受力运动时无法体现其非牛顿物理 属性.Stomakhin等

10、91采用一种用户可控弹塑性本 构模型将雪建模为非牛顿流体，对于流固边界，当 流体和固体之间产生刚性碰撞时使用动摩擦力模 型对流体速度进行更新；当产生黏性碰撞时令边 界处流体速度为零，实现雪花的黏附效果.Zhou等 101使用基于物理的模型的隐式数值积分器对弹塑 性材料进行模拟，但当相邻粒子共面或共线时，积 分器中的变形矩阵的计算将受到影响，使得算法 稳定性变差，这一现象在流固交互界面尤其严重. Stomakhin等11使用材料点法模拟非牛顿流体的 融化和冻结，采用不滑移边界条件处理流固边界, 令与固体直接接触的流体都黏附在固体上，导致 黏度不稳定和额外的计算成本.Peer等121通过规 定涡量

11、扩散模拟高黏度流体的速度梯度场；采用 Akinci等131提出的针对牛顿流体的边界条件处理 流固边界，在该边界条件下，流体粒子只与固体边 界有受力计算，但不会黏附在固体边界上.近年来，研究人员提出了很多传统牛顿流体的 边界处理方法.例如，将固体边界离散为粒子，进行 边界处理时将固体粒子和流体粒子同等对待14-151. Muller等141设定惩罚力防止流体边界粒子穿透固 体边界，但惩罚力导致流体内部的压强变化增大， 需要很小的时间步长维持仿真稳定性.边界处理的 另一种方法是在流体边界处产生ghost粒子阳功， 它与相应的流体粒子具有相同的黏度、质量、密度 和压强力.这种方法在复杂边界处难以实施

12、，并且 仿真的数值稳定性受到边界粒子采样质量的约束. Band等181将移动最小二乘法结合到边界处理中, 通过拟合虚拟边界粒子到流体粒子的距离来计算 边界粒子的贡献，而不是评估真实边界粒子的密 度和压力.2非牛顿流体仿真计算2.1 SPH方法SPH方法是一种拉格朗日粒子法，其将流体在 三维空间上离散为一个个携带物理属性的粒子，每个粒子在X处的物理属性A(x)可以通过对其邻域 X-Xj内所有邻居粒子j进行权函数插值计算得到(A (x U 勺 awJ p 1其中，J W (x-Xyh) , W为光滑核函数，h为 核函数支持半径；mj为邻居粒子J的质量；pj为 邻居粒子J的密度.在SPH中，标量A的

13、梯度计算公式和矢量u 的散度计算公式少为土* *j 仔 + 3Wj(2)-m土* *j 仔 + 3Wj(2)P j 3 p j J2.2运动控制方程非牛顿流体属于连续介质，其运动控制方程 须遵循质量守恒原则和能量守恒原则.在拉格朗 日框架中，非牛顿流体的控制方程为-dP + V- v = 0，p dtdv1 1 = Vp H V - T + g .dt -pp其中，t为时间步长；v为速度；g为由外力产生 的加速度；p为压强；T为剪切应力张量；p为流 体密度；dp应控制为零，以保证流体的不可压缩 dt性；-1 Vp为由压强产生的流动加速度；-V-T pp为不同速度的流体粒子间的摩擦力产生的黏性加

14、 速度.使用通用牛顿模型1模拟剪切应力张量T = D P D可以实现牛顿流体和非牛顿流体的统一建模.其 中，D=1 (Vv + (Vv)T)为剪切速率；#为黏滞系数, 牛顿流体的Dp为常数.速度梯度场在SPH框 架下的表达公式回为_ mVv, =2里(vt - v j) 0VWy(5)j pj将式(5)代入式(2)，得到流动加速度的SPH表达式 土 耻=-1rmj p+p V,W,j将式(6)代入式(3)(4),得到黏性加速度的SPH表达式1V-弓1V-弓 对于非牛顿流体模型，目前存在很多种建模 方法，常用于非牛顿流体仿真的有幂律分布模型网 和Cross模型I.幂律分布模型可用于模拟假塑性 流

15、体和胀流型流体，适用于中等剪切速率的流体, 但当剪切速率较小或较大时模型的仿真效果较差24.Cross模型限制了流体的黏性下界和上界，它们使 流体在很小的剪切速率下表现牛顿流体的特性； 当剪切速率较大的时候，流体则表现出幂律分布 模型的特征.本文采用Cross模型对非牛顿流体进行建模#(| d|)=a +1+kd| (8)其中，K和为控制流体黏度的参数，根据所仿真 的非牛顿流体属性进行取值，通常为正数；四和 心分别为流体在低剪切率和高剪切率时的黏度上 界和黏度下界，0 W四W国，黏度的单位为m2/s . 当K = 0时本文采用Cross模型对非牛顿流体进行建模#(| d|)=a +1+kd|

16、(8)2.4非牛顿流体仿真算法传统的牛顿流体仿真算法通过对控制方程式 进行显式求解实现流体仿真，但由于非牛顿流体 的速度场和黏度互相限制，显式求解时时间步长 必须取很小的值来维持仿真算法稳定性.本文采 用预测-校正方法E对非牛顿流体进行建模.非牛 顿流体的预测-校正方法根据剪切应力和外力预测 流体粒子速度v：；然后求解压强项，为每个流体 粒子i设定一个独立的强度系数k，迭代地求解使 该流体粒子满足局部密度波动小于阈值们进而 维持整个流体系统的不可压缩性.根据剪切应力和外力得到流体粒子i的中间 速度v：后，代入式(1)，得到粒子i的中间密度p： = p,(t)+ mj(v； -v；)VWj(9)

17、J设定pi = k(pi -p0)，则由压强产生的流体速 度变化量1 7V ipkv, = Atp= -At-习m .VWjjmpi j(10)由压强校正而产生的密度变化值(10)Api = At mj (Avy Av j)NWjj将式(9)代入式(10)(11),流体粒子l校正后的密度Pi (t + At) = p*+ & At2%(11)(12)*其中0 %, %Pl2|mjVWj+|mjV(11)(12)*其中0 %, %Pl2将式(6)代入流动加速度，计算每个流体粒子 l在t + At时刻的速度v(t + At)和位置x(t + At),V (t + At) = v* At mj( +

18、 -)VWj(13)X (t + At) = x j (t) + At - vjt + At)(14)非牛顿流体仿真算法步骤如下：Step1.遍历所有流体粒子.对每个流体粒子l查 找邻居粒子N，设定时间步长At.Step2 .根据式计算每个流体粒子l的密度 P = mjWij，根据式(12)计算每个流体粒子i的系数% .Step3.根据式(8)计算流体黏度壮，根据式(7)计算 每个流体粒子i的黏性加速度-V. T.PStep4.根据每个流体粒子i的黏性加速度和外力 计算得到每个流体粒子i的中间速度v* .Step5.当流体平均密度和标准密度之间的偏差大 于或迭代次数小于3时，根据式(9)计算每

19、个流体粒子 i的中间密度p*,根据式(11)计算每个流体粒子i的强 度系数kt,对密度进行迭代求解.Step6.根据式(13)(14)计算每个流体粒子i在 t + At时刻的速度vjt + At)和位置xjt + At).3流固边界处理在流体仿真的过程中，由固体边界产生的边 界力Fb会使流体产生较大的密度波动.例如，高 速下落的流体撞到不可移动的固体边界时，流体 粒子在边界处迅速堆积.受到不可压缩条件约束, 通常需要使用较大的压强力和较小的时间步长维 持仿真的稳定性.此外，由于非牛顿流体具有爬 杆、黏附等物理特性，牛顿流体的自由滑移边界条 件不适用于非牛顿流体的仿真.本文提出一种针 对非牛顿流

20、体的边界处理方法，通过对固体边界 粒子进行加权，计算为每个流体粒子i设定边界力 的法向量n ,并采用Coulomb摩擦力边界条件对 边界处流体粒子的速度进行迭代更新.3.1固体边界处理使用快速泊松盘采样方法回将固体边界采样 为不可移动的、半径为0.5h的边界粒子，每个边界 粒子具有位置xb、法向量nb、质量mb、压强pb属 性.得到固体表面分布的边界粒子后，边界流体粒 子i的密度网为mi Wf + mbWib14 2 414 2 43pi Sdp Jbwdsy(15)Pi其中，PiFd为邻居流体粒子对密度的贡献, Pm 为边界固体粒子的密度贡献.图1所示 为流体边界粒子密度计算示意图.(15)

21、Pi.流体粒子 固体边界粒子图1流边界体粒子密度计算通常定义流体粒子i的法向量25ni = mbVWib -k其中，”,为边界加速度的方向，当边界为平面时, ”，应与平面垂直.然而，当固体边界采样不均匀 时会使”，存在计算偏差，进而影响边界处受力和 加速度的计算精度.为解决这一问题，本文采用加权平均方法为 流体粒子设定法向量对流体粒子和固体边界 进行碰撞检测，|x, -xb| W 0.5h时判定流体粒子穿 透固体边界，将流体粒子校正到距离它最近的4个 固体边界粒子的正中心，如图2所示.校正后俯视图流体粒子；固体边界粒子；固体边界粒子图2流体边界粒子碰撞校正图2中，绿色粒子和黄色粒子均为固体边界

22、粒 子.黄色粒子到流体粒子的距离知国=|(0.5h, 0.5h,0.5h)T ,绿色粒子到流体粒子的距离为 dgreen = | 1.5h, 0.5h, 0.5hJ|,由于黄色粒子距离流 体粒子较近，绿色粒子距离流体粒子较远，因此黄 色粒子比绿色粒子有更高的权重.此时，流体粒子 i的法向量叫=(衣四(烦瞄)+ 8NW(dgreen)nb(16)3.2摩擦力边界条件计算边界流体粒子所受压力和黏力时，根据牛 顿第三定律,Z Fiboundary + FbOundaryi=0，即所有 边界固体粒子给流体粒子i的力等于流体粒子i给 边界固体粒子的力.因此，边界处流体粒子i的流 动加速度af 1和黏性加

23、速度251计算公式分别为+ 4VW厂4 n(17)1 14 4 4 牛2 4P4 4 43p-JlUd忒boundaiya：=E %与 + 与卜 W+Z mb 竺(】8)1 4 444 2 Pf4 4 43 液 4 4p4 43aifluidboundary其中，此bounds为由边界固体粒子对流体粒子i 产生的流动加速度，a：n 为由边界固体粒子 对流体粒子i产生的黏性加速度.根据式(18)计算黏性加速度时，固体边界粒子 给流体粒子的黏力方向为流体速度的反方向，如 图3所示.然而，流体粒子所受边界摩擦力应沿边 界的切线方向与流体运动方向相反.因此，本文采 用式(18)计算黏性加速度会造成计算

24、偏差._ 根据速度得到的黏性加速度方向图3边界处黏性加速度计算方向示意图与剪切应力张量一样，摩擦力是由于粒子间 速度不同而造成的布朗运动所产生的，因此可以 用流体和固体之间的摩擦力计算261替代固体边界 对流体粒子产生的黏性加速度.文献27通过令 vi(t + At) = sv*(t + At)，实现摩擦力对边界流tan体粒子的影响，即计算中间速度v*(t + t)，将速度 的切向分量乘以摩擦系数g，以该值作为新的速度.本文对摩擦力边界的处理与文献27思路相 近.然而，文献27是针对牛顿流体进行的边界处 理.相比之下，非牛顿流体的剪切应力张量r与速 度的梯度相关，因此应在使用摩擦力边界条件更

25、新速度后用这个新速度作为Neumann边界，然后 通过预测-校正方法对流体的密度进行迭代更新， 最后求解出下一时间步长的流体粒子的速度和位置.本文采用Coulomb摩擦力模型28对边界处流 体粒子的速度进行更新/ s v,(t) -nJ)vi(t + At) = max 0,1 rvi(t)(19)vl(t) / 序临 I 其中，vl(t)tan为速度vt(t)沿边界切线方向的分量, s为固体和流体之间的摩擦系数.根据Coulomb 摩擦力模型，摩擦力公式包含静摩擦和动摩擦2种 情况.设由摩擦力产生的速度损耗为g|vi(t) -nJ , 当速度沿切向方向的分量小于速度损耗时为静摩 擦情况，流体

26、粒子i黏附在固体边界上不再移动, 其本质上是不滑移条件；当速度沿切向方向的分 量小于速度损耗时为动摩擦情况，流体粒子的切 向速度更新为被摩擦力损耗后的新速度.本文中， 固体和流体之间不存在弹力，因此流体粒子法向 速度更新为0.当s = 0时，为自由滑移边界.本文算法步骤如下：Step1.遍历所有流体粒子.对每个流体粒子i查 找邻居粒子N ,设定时间步长At.Step2.根据外力预测每个流体粒子i的速度v；和 位置X* .Step3.进行碰撞检测，根据式(16)计算边界处流体 粒子i的法向量nt.Step4.计算每个流体粒子i的密度p；和系数 .Step5.根据式(16)更新边界处流体粒子的速度

27、vf .Step6.根据式计算非边界流体粒子黏性加速度， 更新速度v； = v； + d.Step7.根据式(10)计算每个流体粒子i的中间密度 P；,根据式(12)计算每个流体粒子的强度系数k ,对 密度进行迭代求解.Step8.当流体平均密度和标准密度之间的偏差大 于或迭代次数小于3时，循环Step5Step7.Step9.将式(6)和式(17)代入粒子流动加速度，根 据式(13)(14)计算每个流体粒子i在t + At时刻的速度 vt(t + At)和位置 xi(t + At).在每个时间步长内，首先通过背景网格方法 为每个流体粒子寻找其邻居粒子；然后根据粒子 所受外力，预测粒子下一时间

28、步长时的速度和位 置；根据预测的位置进行流体和固体的碰撞检测， 对穿透固体边界的流体粒子进行位置校正，并根 据式(16)计算边界流体粒子的法向量n,.Step5Step7为考虑了摩擦力边界条件的预测-校正.首先计算每个流体粒子i此时的局部密度, 边界处的流体粒子根据式(15 )进行密度计算；边界 处流体粒子根据摩擦力边界条件更新速度，非边 界处的流体粒子通过计算黏性加速度更新速度. 以新的速度作为Neumann边界，通过式(10)计算 密度中间值，根据式(12)为每个流体粒子计算强度 系数kt ,并对每个粒子的局部密度进行循环迭代, 使平均局部密度逼近标准密度.Step9根据式(13)(14)

29、得到所有流体粒子在下 帧的位置和速度.4实验结果与讨论本文实验的硬件平台如下：Intel Xeon E5-2687W v4(3.0 GHz 12 Cores)处理器，80 GB 400 MHz ECC 内存，NVIDIA Quadro K4000 显卡. 仿真算法使用C+语言实现，使用OpenMP多线程 技术进行数据并行化.使用宽度为的。空间背景网 格对流体模拟中的邻居粒子进行哈希查找，光滑核 函数采用B-样条函数.采用各向异性核函数法对流 体粒子进行表面构造，并使用Marching Cubes算 法提取表面网格.使用开源物理引擎Bullet进行固 体边界模拟，使用Open Graphics

30、Library进行流 体粒子的实时显示.对每一帧的流体粒子进行表面 重构，生成流体的ply网格文件.使用Blender的光 线跟踪引擎Cycles进行渲染，将渲染好的每一帧图 片合成为视频，得到最终的流体模拟动画.为提高仿真效率，本文采用全局动态时间步 长进行仿真29,当流体在自由下落等密度波动 较小的情况下，时间步长可以取较大值；当流体 与固体发生激烈碰撞时，为维持仿真稳定性, 时间步长需取值较小.仿真动态时间步长受 Courant-Friedrich-Levy时间稳定性条件和黏性30 的共同限制，即(At W min 21其中，rmax = (At W min 21i值；安全系数21 1,

31、 22 1.首先按照第2节中的仿真算法实现牛顿流体 和剪切变稀型非牛顿流体的统一建模.在自由滑 移边界条件下1,牛顿流体和非牛顿流体的实验结 果如图4所示.a.牛顿流体enb.剪切变稀型非牛顿流体图a.牛顿流体enb.剪切变稀型非牛顿流体如第3.2节所述，当K = 0时，仿真模型为牛 顿流体.在接触到斜板固体边界后，牛顿流体很快 便失去了初始的兔子形状；当流体滑落至容器底 部时溅起水花，如图4a所示.相比之下，图4b中 的非牛顿流体则在触碰到斜板固体边界后短时间 内仍然可以维持一定的兔子形状，比牛顿流体用 了更长的时间才滑落至容器底部，并且在滑落至 容器底部时也不会产生大量水花.然而，在自由滑

32、 移边界条件下1,只能看出非牛顿流体比牛顿流体 具有更高的黏性，无法看出其黏度与速度梯度之间 的变化关系.因此，只有增加了摩擦力边界条件, 才能使仿真较好地体现出非牛顿流体的物理特性.图5本文方法非牛顿流体仿真效果本文实现了第3节所述流固边界处理方法.图 5 所示为本文提出的摩擦力边界条件下的非牛顿 流体仿真效果.在流体下落过程中，时间步长约为 7 ms;流体触碰到固体斜板边界后，时间步长约为 0.05 ms;仿真进行到第84帧时已经趋于稳定，流体停留在斜板上不再移动，流体形态也不再发生 改变.非牛顿流体由于受到斜板上的边界摩擦力, 速度降低，导致流体的黏度增加.流体黏度增加后 流体粒子间的速

33、度差异降低，因此流体会更倾向 于维持初始形态图5本文方法非牛顿流体仿真效果图6所示为图4和图5仿真实验的渲染结果图. 图6中，各组实验的动画时长均为4 s,动画帧数均 为96帧；兔子形状的流体由25 746个流体粒子组 成，固体斜板倾斜角度为兀/6.从图6a图6c可以 看出，在摩擦力边界条件下，非牛顿流体能够展现 出更明显的物理特性；在自由滑移边界条件下, 虽然能够看出牛顿流体和非牛顿流体的不同，但 由于没有摩擦力为非牛顿流体降速，非牛顿流体 在受重力下落，速度越来越快，黏度也随之降低, 在滑落至容器底部时已经呈现出低黏度牛顿流体 的物理特性；相比之下，在本文摩擦力边界条件 下，非牛顿流体受到斜板边界给的摩擦力阻碍，速 度降低、黏度增加，使得流体的非牛顿特性更加显 著，最终以趋近固体的状态停留在斜板上.从图 6d和图6e可以看出，在自由滑移边界条件下，非 牛顿流体沿