基于PBD和体积约束的薄壳弹性碰撞变形模拟

基于PBD和体积约束的薄壳弹性碰撞变形模拟
张耘齐;孙红岩;孙晓鹏
【摘要】A technique for simulating deformation of elastic thin shell materials was presented,such as the balloon and the inflatable toys.The framework of position-based dynamics was used in this simulation and the details of the local deformation were controlled directly by volume constraint which improved the constraint solving.The geometric model and the elastic physical model of the thin shell were constructed.The stretch constraint was constructed and the volume constraint was constructed combining the internal ball filling method.Based on the position dynamics method and the number of iterations,the deformation intensity and the local details of deformation were controlled.Experimental results show that the proposed method is pared with the explicit Euler method and the implicit Euler method,it has higher stability,higher computation efficiency and lower computation.%针对薄壳弹性碰撞变形问题,提出基于位置动力学的体积约束法.构造薄壳的几何模型和弹性约束物理模型,构造拉伸约束,基于内部球填充方法构造体积约束,基于位置动力学方法求解,以迭代次数控制薄壳模型在碰撞时的弹性强度和局部细节变形.实验结果表明,该算法在弹性变形的局部细节控制方面具有效性,采用位置动力学方法实现薄壳弹性变形,与显式欧拉法等相比,其具有较高的稳定,与隐式欧拉法相比,其计算效率较高,计算开销较低.
【期刊名称】《计算机工程与设计》
【年(卷),期】2018(039)004
【总页数】5页(P1190-1194)
【关键词】薄壳;弹性变形;体积约束;刚度;位置动力学
【作者】张耘齐;孙红岩;孙晓鹏
【作者单位】辽宁师范大学计算机与信息技术学院计算机系统研究所,辽宁大连116029;辽宁师范大学计算机与信息技术学院计算机系统研究所,辽宁大连116029;辽宁师范大学计算机与信息技术学院计算机系统研究所,辽宁大连116029;北京邮电大学智能通信软件与多媒体北京市重点实验室,北京 100876
【正文语种】中文
【中图分类】TP391.9
0 引言
为了解决变形模拟中不符合物理规律的模拟效果问题，基于物理的物体运动模拟逐渐成为了热门研究领域，基于物理的模拟有刚体模拟、弹性物体模拟和流体模拟等，本文主要研究的是弹性薄壳模拟。

目前很多模拟物体运动变形的方法都具有较高的计算精度，但计算开销大，缺乏实时交互性，随着电脑游戏与虚拟现实技术的迅速发展，人们对物理模拟的实时交互性提出了更高的要求，实时模拟应用更注重保持视觉合理前提下的稳定性、健壮性和计算效率。

与以内外力计算更新速度、加速度的数值积分，从而更新运动位置的传统方法不同，基于位置动力学方法(position-based dynamics，PBD)可忽略内力直接计算点的位置、运动过程可控性强，算法较为简单、运动帧更新较快，适用于游戏等实时交互性较强的应用环境[1]。

本文
针对薄壳弹性碰撞变形问题[2]，结合内部球填充方法提出基于位置动力学的体积
约束法，模拟复杂薄壳模型在碰撞时的精确变形，以迭代次数控制薄壳的局部变形
细节，与显式欧拉法相比更加稳定，与隐式欧拉法相比计算效率更高，计算开销更低。

1 相关工作
Bender等基于有限元思想[3]，提出一种基于连续介质力学的运动控制方法，以应变能作为约束函数，隐式积分和小时间步长控制模拟的健壮性，该方法优势在于此方法模拟效果稳定，可以处理复杂的物理模拟，不足之处在于计算开销比较大，模拟速度较慢实时性差[4]。

Müller等通过Green-St Venant应变张量控制应变离散化方向，结合了计算开销较大的FEM方法与PBD方法的优点，但在体积约束刚度过大(接近1时)等极端形变下会产生抖动[5]。

Tournier等提出带有二阶导数的运动约束，将计算效率较低的非线性形变改造为线性形变，不足之处在于健壮性较差，几何刚度只包含一个较精确的一阶算子，不能完全保证形变过程的稳定性[6]。

基于PBD方法的形变方程求解能够提高模拟精度。

Liu等基于局部-全局的方法求解弹簧质点模型变形，以局部方法处理非线性项，以全局方法预分解拉伸变换，其优势是能在短时间内获得精确牛顿解，运用于实时应用，不足之处在于碰撞力较大时缺乏细节[7]。

Wang等提出切比雪夫加速方法，基于前迭代结果加速收敛且算法可以并行化，不足在于以雅克比迭代求解方程组时产生发散，难以加速非线性形变[8]。

Narain等基于交替方向乘子算法处理隐式时间积分和非线性本构模型的刚性约束，其优势在于可并行、健壮性较高，不足在于当求解物体质量非常小的方程时，可能无法收敛[9]。

Lin等对碰撞检测方法进行了分类，系统介绍多种几何模型的表面碰撞检测算法以及碰撞检测算法的核心问题[10]。

Teschner等针对虚拟手术应用背景，详细分析介绍了层次包围盒、空间划分、图像空间技术和随机方法等碰撞检测方法以及基于GPU加速图像空间技术[11]。

2 约束求解
本文以单层三角网格模型G={V，E}表示薄壳的几何结构，其中V为待变形物体的N个顶点集合，设每个顶点的质量为mi，位置为xi，当前速度为vi，记为质点，其中i=1,2,…,N。

设薄壳模型G的质点运动约束为Cj，j=1,2,…,M。

约束Cj具有如下属性：质点下标集合index和质点数量nj，表达式Cj，定义约束强度大小的刚度参数kj∈[0,1]。

若约束表达式有等式约束Cj(x1,…,xnj)=0和不等式约束Cj(x1,…,xnj)≥0两种形式[1]。

对于给定时间步长Δt，薄壳G在碰撞后的变形模拟处理流程如下：
(1)初始化所有的质点的位置xi、速度vi、质量mi；
(2)以显式积分法[12]预测所有质点的下一位置pi；
(3)对xi到pi的运动路径进行碰撞检测，并进行碰撞处理[10]；
(4)将pi代入约束，以非线性高斯-赛德尔迭代法[1]求解约束，修正预测质点位置pi；
(5)更新质点速度(pi-xi)/Δt→vi，更新质点位置pi→xi，跳至(3)。

为修正xi的预测点pi，使其满足所有约束，需要求解由约束组成的方程组，方程
组由M个方程和3N个坐标分量组成，如果M>3N(M<3N)，方程组超定(欠定)，所求方程一般为非线性。

假设约束只作用在单个质点pi上，需要求解方程组为
其中，是1×N维向量，由所有参数(即P的N个)的导数构成，和都是常量。

若M=3N且为等式约束，方程组可以用任意线性方法求解。

若M≠3N，方程组矩阵无逆，无法直接求解方程组。

本文基于位置动力学方法，使用非线性高斯-赛德
尔法[1]，将单质点约束拓展到多质点约束C(p1,…,pn)上
其中
式中：wi=1/mi，代入式(1)即可求出各约束下相关质点的位置。

3 形变约束
本文薄壳碰撞变形约束的类型包括拉伸约束和体积约束两种。

3.1 拉伸约束
在拉伸约束中，距离约束质点迭代移动方向如图1所示，对于给定距离d，令其距离约束函数为C(p1,p2)=(p1-p2)2-d2，令约束函数梯度为▽p1C(p1,p2)=(p1-
p2)/|p1-p2|以及和▽p2C(p1,p2)=-(p1-p2)/|p1-p2|(文献[1])，代入式(1)，得到
图1 距离约束
3.2 体积约束
本文在考虑封闭薄壳碰撞变形的过程中加入了体积约束，以提高薄壳形变过程模拟的真实性。

以封闭球壳下落撞击地面弹起为例，球壳撞击地面后发生形变，球壳内气体体积瞬间压缩，为了恢复原有体积，球壳在垂直于地面的方向上产生向上的弹力，将球壳弹起并恢复原有体积。

对于以三角网格表示的球壳表面，连接球心与球壳表面任意三角形的3个顶点，从而得到以封闭四面体(如图2所示)。

对于复杂模型，本文基于文献[13]方法实现球体逼近填充，进而实现复杂封闭壳的四面体体积约束(如图3所示)。

图2 球体表面三角形与球心组成一个四面体
本文基于文献[14]的方法，将单个四面体体积的约束函数定义为
(p4-p1)-V0
图3 以内部球填充薄壳模型
式中：p1，p2，p3和p4是四面体4个顶点，V0是四面体静止体积。

4 实验结果及分析
本文实验的软硬件环境为************************，16 GB RAM，64位Windows 7操作系统。

所用的GPU为NVIDIA Quadro 2000，核心数192，显存1 GB GDDR5。

编程平台为Microsoft Visual Studio 2008，CUDA6.5以及OpenGL。

考虑计算效率，本文基于文献[10]给出的包围球法处理各种碰撞。

图4为未加入体积约束的球壳碰撞弹起形变关键帧，图5为加入体积约束的球壳碰撞弹起形变关键帧。

对比图4和图5的视觉效果可知，封闭薄壳形变过程在加入体积约束和未加入体积约束两种情况下，存在显著的差异。

仅从视觉观察即可发现，在加体积约束后球壳碰撞变形的模拟效果更加真实，更符合客观的物理效果。

正常情况下，封闭球壳内具有一定体积的气体，球壳落地碰撞变形过程中，加入体积约束的球壳内保持一定的体积，而不加体积约束的球壳在碰撞变形壳内的体积减少，导致在变形球壳与地面接触面的边缘部分看起来更加“锐利”，不符合客观自然情况下，球壳碰撞变形时接触面边缘部分较为“圆润”的实际效果。

(由于截图问题，球体略显大小不一)。

图4 未加入体积约束的球壳碰撞弹起形变关键帧
图5 加入体积约束的球壳碰撞弹起形变关键帧
本文方法的迭代求解次数控制薄壳变形细节，同时能够表现球壳的软硬差异。

本文举例了3种不同模型在不同迭代次数下的模拟效果。

如图6所示，cloth模型在迭代次数较少时顶端变形拉伸较大，褶皱细节较多，在迭代次数较多时，顶端变形拉伸较小，褶皱较少，显得更加“僵硬”。

sphere模型在迭代次数较少时(如图7所示迭代3次)由于自身重力无法保持原有形状，迭代次数较多时可以很好地保持原有的形状。

图6 cloth模型在不同约束迭代次数的模拟效果
图7 sphere模型在不同约束迭代次数的模拟效果
duck模型在迭代次数较少时(如图8所示迭代50次)薄壳较软、变形较大、细微变形较多；随着迭代次数的增加模型刚度也越来越大，在迭代次数较多时(如图9所示迭代150次)，球壳碰撞时所产生的细微形变较少。

本文以球壳碰撞变形时的质点加速度度量球壳上每个质点的变形程度，从而得到约束迭代50次落地碰撞变形关键帧b、d和f和约束迭代150次落地碰撞变形关键帧b、d和f的能量分布(如图10所示)。

显然，在迭代次数较多时，薄壳表面能量分布更为集中，表明随着迭代次数的增加，薄壳整体性较强、变形区域较小，局部细微变形减少。

图8 duck模型约束迭代50次落地碰撞变形关键帧
图9 duck模型约束迭代150次落地碰撞变形关键帧
图10 duck模型迭代50次和150次时关键帧b、d和f的能量分布对比
3种模型的复杂程度不同所需要的迭代求解次数也不相同，模型顶点越多、受力程度越大、变形越复杂所需的迭代次数也越多，反之，则所需迭代次数越少。

迭代次数太多不仅会使物体太“僵硬”从而损失变形的细节更增大了计算开销，迭代次数太少会使物体太“柔软”，不能很好地保持原有形状，所以在模拟过程中要选择合适的迭代次数，除了模型复杂度之外还需要根据模拟的对象调整迭代次数。

根据3种不同模型，合适的迭代次数范围见表1。

表1 3种模型的合适迭代次数模型名称顶点数面片数迭代次数范围
cloth2344081-3sphere642128030-50duck21084212150-200
5 结束语
本文基于位置动力学和体积约束，针对薄壳碰撞变形问题给出了高精度的细节变形模拟，相关实验结果表明了本文算法的有效性，相比质点弹簧等方法，本文方法计算量较小，稳定性较好。

本文简单地以包围球法方法实现碰撞检测，碰撞检测精度
有待提高；另外，以球体填充逼近复杂薄壳以构造四面体、实现体积约束的过程中，球体填充方法紧凑型较差，对最终形变效果的逼真性有所影响。

参考文献：
[1]Bender J,Müller M,Otaduy MA,et al.A survey on position-
based simulation methods in computer graphics[J].Computer Graphics For um,2015,33(6):228-251.
[2]Terzopoulos D,Platt J,Barr A,et al.Elastically deformable models[C]//Conf erence on Computer Graphics and Interactive Techniques.NY:ACM,1987:20 5-214.
[3]O′Brien JF,Hodgins JK.Graphical modeling and animation of brittle fractu re[C]//Conference on Computer Graphics and Interactive Techniques,1999: 137-146.
[4]Bender J,Dan K,Charrier P,et al.Position-
based simulation of continuous materials[J].Computers & Graphics,2014,44 :1-10.
[5]Müller M,Chentanez N,Kim TY,et al.Strain based dyna-
mics[C]//Proceedings of the ACM SIGGRAPH/Eurograp-
hics Symposium on Computer Animation.Switzerland:Eurographics Associa tion,2014:149-157.
[6]Tournier M,Nesme M,Gilles B.Stable constrained dynamics[J].Acm Transa ctions on Graphics,2015,34(4):1-10.
[7]Liu T,Bargteil AW,O′Brien JF,et al.Fast simulation of mass-
spring systems[J].Acm Transactions on Graphics,2013,32(6):214.
[8]Wang H.A chebyshev semi-
iterative approach for accelerating projective and position-
based dynamics[J].Acm Transactions on Graphics,2015,34(6):246.
[9]Narain R,Overby M,Brown GE.ADMM⊇ projective dynamics:Fast simulati on of general constitutive models[C]//Proceedings of the ACM SIGGRAPH/ Eurographics Symposium on Computer Animation.Switzerland:Eurographic s Association,2016:21-28.
[10]Lin MC,Gottschalk S.Collision detection between geometric[J].Visual Co mputer,1998,11(10):542-561.
[11]Teschner M,Kimmerle S,Heidelberger B,et al.Collision detection for def ormable objects[J].Computers and Grap-hics,2006,30(3):432-438.
[12]Breen DE,House DH,Wozny MJ.Predicting the drape of woven cloth usi ng interacting particles[C]//Proceedings of the 21st Annual Conference on Computer Rraphics and Interactive Techniques.NY:ACM,1994:365-372. [13]Weller R,Zachmann G.Inner sphere trees and their application to collisi on detection[C]//Virtual Realities-Dagstuhl Seminar.DBLP,2008:181-201.
[14]Diziol R,Bender J,Bayer D.Robust real-
time deformation of incompressible surface meshes[C]//Acm Siggraph/Eur ographics Symposium on Computer Animation,2011:237-246.。