布料仿真中常用积分方法

来源:互联网 发布:3d建模软件 mac 编辑:程序博客网 时间:2024/03/29 14:16

1. 简介

布料仿真中,我们通常将布料剖分为三角形网格(或四边形网格),并用弹簧-质点模型构造动力学系统:质点即三角形的顶点,弹簧即三角形的边。质点在外力(如,重力)和内力(弹簧力)的作用下根据牛顿第二定律运动:

{x˙=vMv˙=Fexternal+Finternal

这就是我们需要求解的常微分方程,通常我们很难求其解析解,只能求其数值解。常用方法有:显式/隐式欧拉法,Symplectic Euler,Midpoint method,Leapfrog integration等。

2. 显式欧拉法(Explicit Euler)

形如,

{xn+1=xn+Δtvnvn+1=vn+M1Fn

称为显式积分,可以看出,位置和速度的积分用的都是显式方法。显式积分需控制积分步长Δt,以保证数值求解的稳定性。这里例证了显式条件稳定性。

3. 隐式欧拉法(Implicit Euler)

形如,

{xn+1=xn+Δtvn+1vn+1=vn+M1Fn+1

称为显式积分,可以看出,位置和速度的积分用的都是隐式方法。从理论上说,隐式欧拉法是无条件稳定的。即,积分步长可以任意大,数值求解都稳定(稳定不意味着精确)。然而,隐式方法需要求解非线性方程,通常要线性近似,因此实际中,隐式方法并不是无条件稳定的。即,仍要控制积分步长以保证数值求解的稳定。
这里例证了隐式欧拉法的无条件稳定性。
另注,隐式欧拉法并不是唯一的隐式方法,例如下面将要介绍的midpoint method也属于隐式方法。只是由于其简单及良好性能,隐式欧拉应用得更广泛。

4. Symplectic Euler

形如(1),

{xn+1=xn+Δtvnvn+1=vn+M1Fn+1

或(2),
{xn+1=xn+Δtvn+1vn+1=vn+M1Fn

称为Symplectic Euler。

  • 可以看出(1)位置的更新使用显式方法,速度的更新使用隐式方法;(2)位置的更新使用隐式方法,速度的更新使用显式方法。因此Symplectic Euler又称为semi-implicit Euler。
  • 由于(1)与隐式欧拉法的计算量相同,因此实际中很少使用(1),而直接使用隐式欧拉法;
  • (2)与显式欧拉法的计算量相同,但比显式欧拉法更稳定,可以使用更大的积分步长,因此实际中很少使用显式欧拉法,而使用(2)。从理论上分析,仿真谐振动系统(harmonic oscillation system),只要保证Δt<2ω=2m/k,形式(2)的Symplectic Euler就是稳定的。

5. 中点法(Midpoint method)

形如,

yn+1=yn+hf(tn+h2,yn+h2f(tn,yn))

称为显式中点法(Explicit Midpoint method),又称为改进的欧拉法(modified Euler method)。
形如,
yn+1=yn+hf(tn+h2,12(yn+yn+1))

称为隐式中点法(Implicit Midpoint method)
由上面两个式子可以看出,f总是计算在tntn+1中间时刻t=tn+h2的值f(tn+h2,y(tn+h2))
由两种方法估计f(tn+h2,y(tn+h2))的值:
1. 将y(tn+h2)tn处做一阶泰勒展开,得,
y(tn+h2)y(tn)+h2f(tn,y(tn))=yn+h2f(tn,yn)

2. 用y(tn)y(tn+h)的均值近似:
y(tn+h2)12(y(tn)+y(tn+h))=12(yn+yn+1)

  • 中点法的局部误差为O(h3),全局误差为O(h2)。因此比欧拉法的O(h2)局部误差,O(h)全局误差要精确。当然,隐式欧拉法的计算量更大。
  • 中点法实际上是Runge-Kutta method(以后会专门介绍)的一种特例。

以上几种方法又可以写成如下通用形式:

{1hΔx=(1β)vn+βvn+11hMΔv=(1α)Fn+αFn+1

- 当α=β=0时,为显式欧拉法
- 当α=β=1时,为隐式欧拉法
- 当α=0,β=1α=1,β=0时,为Symplectic Euler(Semi-implicit Euler)。
- 当α=β=12时,参考论文3 Cloth simulation using soft constraints 错误的将其认为是隐式中点法,而实际上中点法的核心是求tn+h2时刻的导数Fn+12,并不是tn时刻和tn+1时刻导数的均值12(Fn+Fn+1)。当α=β=12时,应称为 显式隐式混合法(mixed explicit-implicit),实际上除前三种情况外,都可以称为显式隐式混合法。


6. Leapfrog method

形如,

xn+1=xn+Δtvn+12vn+12=vn12+M1ΔtFn

称为Leapfrog method。此式又可改写成如下形式:
{xn+1=xn+Δtvn+12M1FnΔt2vn+1=vn+12M1(Fn+Fn+1)Δt

由此可知,Leapfrog method也是一种隐式方法。另外,Leapfrog method的局部误差为O(h3),全局误差为O(h2)。只要时间步长Δt1/ω,该方法就是稳定的。

7. 总结

  1. 一般来说显式法稳定性差,需要较小的积分步长,但能更好的保持能量(energy preserving);
  2. 隐式法稳定性好,可以使用更大的积分步长,但会增加能量耗散(energy dissipation)。
  3. 欧拉法(显式和隐式),Symplectic Euler的局部误差为O(h2),全局误差为O(h);中点法,Leapfrog method的局部误差为O(h3),全局误差为O(h2)
  4. 有时候我们会对各种力分别做显示和隐式处理,称为IMEX方法。
  5. 一般情况下,减小时间步长Δt,显式欧拉误差的下降速度要稍快与隐式欧拉。

8. 参考资料

  1. 维基百科 Leapfrog integration词条,Semi implicit Euler method词条,Midpoint method词条。
  2. Interactive simulation of elastic deformable materials
    Servin M, Lacoursiere C, Melin N. SIGRAD 2006.
  3. Frâncu M, Moldoveanu F. Cloth simulation using soft constraints[J]. 2015.
1 0
原创粉丝点击