1. 拉格朗日视角 vs 欧拉视角
拉格朗日视角:其中的元素跟随一起移动。
欧拉视角 :其中的元素固定在对应的地方,网格的点对应的速度。
2. 弹簧质点系统
胡克定律
f i j = − k ( ∥ x i − x j ∥ 2 l i j ) x i j ^ f i = ∑ j j ≠ i f i j \begin{aligned} \bold{f}_{ij} &= -k(\Vert \bold{x}_i - \bold{x}_j \Vert_2 l_{ij})\hat{\bold{x}_{ij}} \\ \bold{f}_i &= \sum^{j \neq i}_j {\bold{f}_{ij}} \end{aligned} fijfi=−k(∥xi−xj∥2lij)xij^=j∑j=ifij
牛顿第二定律
∂ v i ∂ t = 1 m i f i ∂ x i ∂ t = v i \begin{aligned} \cfrac{\partial \bold{v}_i}{\partial t} &= \cfrac{1}{m_i}\bold{f}_i \\ \cfrac{\partial \bold{x}_i}{\partial t} &= \bold{v}_i \end{aligned} ∂t∂vi∂t∂xi=mi1fi=vi
2.1 时间积分
前欧拉(显式)积分器 :根据现有状态推测以后的状态
v t + 1 = v t + Δ t f t m x t + 1 = x t + Δ t v t \begin{aligned} \bold{v}_{t+1} &= \bold{v}_t + \Delta t \cfrac{\bold{f}_t}{m} \\ \bold{x}_{t+1} &= \bold{x}_t + \Delta t \bold{v}_t \\ \end{aligned} vt+1xt+1=vt+Δtmft=xt+Δtvt
半隐式欧拉积分器 :根据现有状态推测以后的状态,速度根据推测出的速度计算
v t + 1 = v t + Δ t f t m x t + 1 = x t + Δ t v t + 1 \begin{aligned} \bold{v}_{t+1} &= \bold{v}_t + \Delta t \cfrac{\bold{f}_t}{m} \\ \bold{x}_{t+1} &= \bold{x}_t + \Delta t \bold{v}_{t+1} \\ \end{aligned} vt+1xt+1=vt+Δtmft=xt+Δtvt+1
2.2 代码
- 计算速度: v t + 1 = v t + Δ t f t m \bold{v}_{t+1} = \bold{v}_t + \Delta t \cfrac{\bold{f}_t}{m} vt+1=vt+Δtmft
- 与地面碰撞(防止计算位移之后陷入地下)
- 计算位移: x t + 1 = x t + Δ t v t + 1 \bold{x}_{t+1} = \bold{x}_t + \Delta t \bold{v}_{t+1} xt+1=xt+Δtvt+1
@ti.kernel
def substep():
# Compute force and new velocity
n = num_particles[None]
for i in range(n):
v[i] *= ti.exp(-dt * damping[None]) # 阻尼
total_force = ti.Vector(gravity) * particle_mass
for j in range(n):
# 两个粒子有弹簧时
if rest_length[i, j] != 0:
x_ij = x[i] - x[j]
total_force += -spring_stiffness[None] * (x_ij.norm() - rest_length[i, j]) * x_ij.normalized()
v[i] += dt * total_force / particle_mass
# Collide with ground
for i in range(n):
if x[i].y < bottom_y:
x[i].y = bottom_y
v[i].y = 0
# Compute new position
for i in range(num_particles[None]):
x[i] += v[i] * dt
NOTE :
- 其中
Δ
t
\Delta t
Δt 需要满足一下情况:
Δ t = c m k ( c ∼ 1 ) \Delta t = c\sqrt{\cfrac{m}{k}} (c \sim 1) Δt=ckm(c∼1)
2.3 质点系统
隐式时间积分 :
x t + 1 = x t + Δ t v t + 1 v t + 1 = v t + Δ t M − 1 f ( x t + 1 ) \begin{aligned} \bold{x}_{t+1} &= \bold{x}_t + \Delta t \bold{v}_{t+1} \\ \bold{v}_{t+1} &= \bold{v}_t + \Delta t \bold{M}^{-1}\bold{f}(\bold{x}_{t+1}) \\ \end{aligned} xt+1vt+1=xt+Δtvt+1=vt+ΔtM−1f(xt+1)
NOTE :其中 $ \bold{M} $ 是质量矩阵。
消除 x t + 1 \bold{x}_{t+1} xt+1 :
v t + 1 = v t + Δ t M − 1 f ( x t + Δ t v t + 1 ) \bold{v}_{t+1} = \bold{v}_t + \Delta t \bold{M}^{-1}\bold{f}(\bold{x}_t + \Delta t \bold{v}_{t+1}) \\ vt+1=vt+ΔtM−1f(xt+Δtvt+1)
线性化(牛顿法的一步):
v t + 1 = v t + Δ t M − 1 [ f ( x t ) + ∂ f ∂ x ( x t ) Δ t v t + 1 ] \bold{v}_{t+1} = \bold{v}_t + \Delta t \bold{M}^{-1} \left[ {\bold{f}(\bold{x}_t) + \cfrac{\partial \bold{f}}{\partial \bold{x}}(\bold{x}_t)\Delta t\bold{v}_{t+1} } \right] vt+1=vt+ΔtM−1[f(xt)+∂x∂f(xt)Δtvt+1]
3. 形变
形变映射 ϕ \phi ϕ :静止材料位置 ↦ \mapsto ↦ 形变材料位置
x d e f o r m e d = ϕ ( x r e s t ) \bold{x_{deformed}} = \phi (\bold{x_{rest}}) xdeformed=ϕ(xrest)
形变梯度 F \bold{F} F :
F ≔ ∂ x d e f o r m e d ∂ x r e s t \bold{F} \coloneqq \cfrac{\partial \bold{x_{deformed}}}{\partial \bold{x_{rest}}} F:=∂xrest∂xdeformed
Note :
- 平移过程中形变梯度保持不变: ϕ 1 = ϕ ( x r e s t ) \phi_1 = \phi(\bold{x_{rest}}) ϕ1=ϕ(xrest) 和 ϕ 1 = ϕ ( x r e s t ) + c \phi_1 = \phi(\bold{x_{rest}}) + c ϕ1=ϕ(xrest)+c 有相同的形变梯度
变形/静止体积比 J J J :
J = det ( F ) J = \det(\bold{F}) J=det(F)
Note :
- 三维空间中,矩阵行列式的性质,即体积变化倍数
4. 弹性体
4.1 超弹性体(Hyperelasticity)
超弹性材料:应力 – 应变关系 由 应变能密度函数 定义 :
Ψ = Ψ ( F ) \Psi = \Psi(\bold{F}) Ψ=Ψ(F)
直观理解: Ψ \Psi Ψ 是惩罚形变的势函数。
应力 :材料的内部弹性力。
应变 :现在用 形变梯度 F \bold{F} F 代替。
Note :
- Ψ \Psi Ψ 是应变能量密度函数
- ϕ \phi ϕ 是形变映射
4.2 应力张量
应力:代表材料微元施加在其周围为材料微元的内力。
4.2.1 三种常用的应力张量
- Piola - Kirchhoff 应力张量 (PK1):
P ( F ) = ∂ Ψ ( F ) ∂ F \bold{P(F)} = \cfrac{\partial \Psi(\bold{F})}{\partial \bold{F}} P(F)=∂F∂Ψ(F)
(计算简单,是在静止空间计算,需要变换到形变空间)
- 基尔霍夫应力(Kirchhoff stress) :
τ
\tau
τ
- 柯西应力张量(Cauchy stress tensor) :
σ
\sigma
σ
(对称,因为角动量守恒)
4.2.2 关系式
- τ = J σ = P F T \tau = J\sigma = \bold{PF}^T τ=Jσ=PFT
- P = J σ F − T \bold{P} = J\sigma\bold{F}^{-T} P=JσF−T
Note:
- F − T \bold{F}^{-T} F−T :补偿材料变形
- 用 F − T \bold{F}^{-T} F−T 替代 F − 1 \bold{F}^{-1} F−1 :因为变换的是法向量 n \bold{n} n 而不是 x \bold{x} x
4.2.3 牵引力
- t = σ T n \bold{t} = \sigma^T\bold{n} t=σTn
直观来说,将法向量乘以应力张量即可获得材料向周围微元施加的力。
4.3 弹性模量(各向同性材料)
-
杨氏模量 :应力张量与应变张量的比值
E = σ ϵ E = \cfrac{\sigma}{\epsilon} E=ϵσ -
体积模量 :压强关于体积的变化,常用于液体
K = − V d P d V K = - V \cfrac{dP}{dV} K=−VdVdP -
泊松比 :
ν ∈ [ 0.0 , 0.5 ) \nu \in [0.0,0.5) ν∈[0.0,0.5)
Note :
- 辅助项具有负泊松比;
- 泊松比为 0 时,拉长物体时,截面积不会发生变化;
- 泊松比较大时,物体会尽量保持体积不变,例如在拉长物体时,物体会变细。
拉梅常数(Lamé parameters) :
- Lamé 第一参数:
λ
\lambda
λ
表示材料的压缩性,等价与体弹性模量或者杨氏模量
- Lamé 第二参数:
μ
\mu
μ
表示材料的剪切模量,用 G 表示
换算公式 :
-
K
=
E
3
(
1
−
2
ν
)
K = \cfrac{E}{3(1 - 2\nu)}
K=3(1−2ν)E (常用于可压缩液体)
-
λ
=
E
ν
(
1
+
ν
)
(
1
−
2
ν
)
\lambda = \cfrac{E\nu}{(1 + \nu)(1 - 2\nu)}
λ=(1+ν)(1−2ν)Eν
- μ = E 2 ( 1 + ν ) \mu = \cfrac{E}{2(1 + \nu)} μ=2(1+ν)E
广义胡克定律 :均匀和各向同性的材料
σ = 2 μ ϵ + λ t r ( ϵ ) \sigma = 2 \mu \epsilon + \lambda tr(\epsilon) σ=2μϵ+λtr(ϵ)
5. 常见的超弹性模型
- 经典的 MPM 方法将每个粒子看做材料的一个微元,材料的本构模型会有一个能量密度函数 Ψ \Psi Ψ ;
- 对能量密度函数 Ψ \Psi Ψ 关于整个模型求积分,得到整个材料的势能;
- 势能对材料点的形变梯度进行求导: P ( F ) = ∂ Ψ ∂ F P(\bold{F}) = \cfrac{\partial \Psi}{\partial \bold{F}} P(F)=∂F∂Ψ ;
- 物理意义上来说,势能对位置求导的结果就是力, P ( F ) P(\bold{F}) P(F) 可以看做材料点的受力。
6.1 Neo-Hookean 模型
适用于各向同性材料
-
Ψ
(
F
)
=
μ
2
∑
i
[
(
F
T
F
)
i
i
−
1
]
−
μ
log
(
J
)
+
λ
2
log
2
(
J
)
\Psi(\bold{F}) = \cfrac{\mu}{2} \sum_i [(\bold{F}^T\bold{F})_{ii} - 1] - \mu \log(J) + \cfrac{\lambda}{2} \log^2(J)
Ψ(F)=2μ∑i[(FTF)ii−1]−μlog(J)+2λlog2(J)
- P ( F ) = ∂ Ψ ∂ F = μ ( F − F T ) + λ log ( J ) F − T P(\bold{F}) = \cfrac{\partial \Psi}{\partial \bold{F}} = \mu(\bold{F} - \bold{F}^T) + \lambda\log(J)\bold{F}^{-T} P(F)=∂F∂Ψ=μ(F−FT)+λlog(J)F−T
因为 Neo-Hookean 模型容易造成能量流失,这时可以考虑 Corotated 模型。
FEM 中应用 Neo-Hookean 模型的示例代码 :
dim = 2
E, nu = 1000, 0.3
la = E * nu / ((1 + nu) * (1 - 2 * nu))
mu = E / (2 * (1 + nu))
element_V = 0.01
x = ti.Vector(dim, dt=real, shape=n_nodes, needs_grad=True)
v = ti.Vector(dim, dt=real, shape=n_nodes)
total_energy = ti.var(dt=real, shape=(), needs_grad=True)
# 计算势能
@ti.kernel
def compute_total_energy():
for i in range(n_elements):
......# get F
# NeoHookean
I1 = (F @ F.transpose()).trace()
# 防止 J 的值过小引起错误
J = max(0.2, F.determinant())
element_energy_density = 0.5 * mu * (I1 - dim) - mu * ti.log(J) + 0.5 * la * ti.log(J)**2
total_energy[None] += element_energy_density * element_V
# 渲染循环
while True:
for s in range(30):
# 调用 taichi 的自动微分器
# 定义的损失函数:compute_total_energy()
# 函数值:total_energy
# 变量值:x
# 微分结果(i 点上的力) :x.grad[i]
# f = - \partial (total_energy) / \partial x
with ti.Tape(total_energy):
compute_total_energy()
......
Note :
- 在 FEM 中的势能(后续有详细说明):
U ( e ) = ∫ e Ψ ( F ( x ) ) d x = V e Ψ ( F e ) U(e) = \int_e\bold{\Psi(F(x))\mathrm{d}x} = V_e\bold{\Psi(F_e)} U(e)=∫eΨ(F(x))dx=VeΨ(Fe)
6.2 (Fixed)Corotated 模型
-
Ψ
(
F
)
=
μ
∑
i
(
σ
i
−
1
)
2
+
λ
2
(
J
−
1
)
2
\Psi(\bold{F}) = \mu \sum_i(\sigma_i - 1)^2 + \cfrac{\lambda}{2}(J - 1)^2
Ψ(F)=μ∑i(σi−1)2+2λ(J−1)2
- P ( F ) = ∂ Ψ ∂ F = 2 μ ( F − R ) + λ ( J − 1 ) J F − T P(\bold{F}) = \cfrac{\partial \Psi}{\partial \bold{F}} = 2\mu(\bold{F} - \bold{R}) + \lambda(J - 1)J\bold{F}^{-T} P(F)=∂F∂Ψ=2μ(F−R)+λ(J−1)JF−T
Note :
- σ i \sigma_i σi 是 F \bold{F} F 的奇异值。
6.3 MPM 教程
《The Material Point Method for Simulating Continuum Materials》
6. 有限元基础
有限元法 :建立离散模型的 Galerkin 离散格式。
6.1 线性四面体(三角形)有限元法
线性四面体有限元(用于弹性)假设 形变映射 p h i phi phi 是一个仿射变换,因此 形变梯度 F \bold{F} F 在单个四面体单元内是恒定的,对单个元素来说:
x d e f o r m e d = F x r e s t + p \bold{x_{deformed}} = \bold{Fx_{rest}} + \bold{p} xdeformed=Fxrest+p
对于每个元素 e e e ,对能量密度函数求体积的积分,可以计算其弹性势能:
U ( e ) = ∫ e Ψ ( F ( x ) ) d x = V e Ψ ( F e ) U(e) = \int_e\bold{\Psi(F(x))\mathrm{d}x} = V_e\bold{\Psi(F_e)} U(e)=∫eΨ(F(x))dx=VeΨ(Fe)
- 其中 形变梯度 F \bold{F} F 在元素 e e e 上是个常数 F e \bold{F_e} Fe ,即 Ψ ( F e ) \bold{\Psi(F_e)} Ψ(Fe) 也为常数,因此可以直接得到积分结果。
6.2 计算形变梯度
x d e f o r m e d = F x r e s t + p \bold{x_{deformed}} = \bold{Fx_{rest}} + \bold{p} xdeformed=Fxrest+p
在 2D 三角形元素(三维空间中是四面体元素)中,设:
- 静止时的顶点位置: a r e s t \bold{a_{rest}} arest , b r e s t \bold{b_{rest}} brest , c r e s t \bold{c_{rest}} crest ;
- 变形后的顶点位置: a d e f o r m e d \bold{a_{deformed}} adeformed , b d e f o r m e d \bold{b_{deformed}} bdeformed , c d e f o r m e d \bold{c_{deformed}} cdeformed 。
因为在线性三角形元素中 F \bold{F} F 是常数,则有:
a d e f o r m e d = F a r e s t + p b d e f o r m e d = F b r e s t + p c d e f o r m e d = F c r e s t + p \begin{aligned} \bold{a_{deformed}} &= \bold{Fa_{rest}} + \bold{p} \\ \bold{b_{deformed}} &= \bold{Fb_{rest}} + \bold{p} \\ \bold{c_{deformed}} &= \bold{Fc_{rest}} + \bold{p} \end{aligned} adeformedbdeformedcdeformed=Farest+p=Fbrest+p=Fcrest+p
消除 p \bold{p} p :
( a d e f o r m e d − c d e f o r m e d ) = F ( a r e s t − c r e s t ) ( b d e f o r m e d − c d e f o r m e d ) = F ( b r e s t − c r e s t ) \begin{aligned} (\bold{a_{deformed}} - \bold{c_{deformed}}) &= \bold{F}(\bold{a_{rest}} - \bold{c_{rest}}) \\ (\bold{b_{deformed}} - \bold{c_{deformed}}) &= \bold{F}(\bold{b_{rest}} - \bold{c_{rest}}) \end{aligned} (adeformed−cdeformed)(bdeformed−cdeformed)=F(arest−crest)=F(brest−crest)
现在 F 2 × 2 \bold{F}_{2\times 2} F2×2 有四个线性约束:
B = [ a r e s t − c r e s t ∣ b r e s t − c r e s t ] − 1 D = [ a d e f o r m e d − c d e f o r m e d ∣ b d e f o r m e d − c d e f o r m e d ] F = D B \begin{aligned} \bold{B} &= [\bold{a_{rest}} - \bold{c_{rest}} | \bold{b_{rest}} - \bold{c_{rest}}]^{-1} \\ \bold{D} &= [\bold{a_{deformed}} - \bold{c_{deformed}} | \bold{b_{deformed}} - \bold{c_{deformed}}] \\ \bold{F} &= \bold{D}\bold{B} \end{aligned} BDF=[arest−crest∣brest−crest]−1=[adeformed−cdeformed∣bdeformed−cdeformed]=DB
其中 B \bold{B} B 在整个物理过程中是常数。因此可进行预计算。
B = ti.Matrix(dim, dim, dt=real, shape=n_elements)
vertices = ti.var(dt=ti.i32, shape=(n_elements, 3))
@ti.func
def compute_D(i):
a = vertices[i, 0]
b = vertices[i, 1]
c = vertices[i, 2]
return ti.Matrix.cols([x[b] - x[a], x[c] - x[a]])
@ti.kernel
def compute_B():
for i in range(n_elements):
B[i] = compute_D(i).inverse()
@ti.kernel
def compute_total_energy():
for i in range(n_elements):
D = compute_D(i)
F = D @ B[i]
......
6.3 显式时间积分
v t + 1 , i = v t , i + Δ t f t , i m i x t + 1 , i = x t , i + Δ t v t + 1 , i \begin{aligned} \bold{v}_{t+1,i} &= \bold{v}_{t,i} + \Delta t\cfrac{\bold{f}_{t,i}}{m_i} \\ \bold{x}_{t+1,i} &= \bold{x}_{t,i} + \Delta t\bold{v}_{t+1,i} \end{aligned} vt+1,ixt+1,i=vt,i+Δtmift,i=xt,i+Δtvt+1,i
- v t , i \bold{v}_{t,i} vt,i 和 x t , i \bold{x}_{t,i} xt,i 存储在顶点中。
弹性势能对位置求导,结果的相反值即为顶点的受力:
f t , i = − ∂ U ∂ x i = − ∑ e ∂ U ( e ) ∂ x i = − ∑ e V e ∂ Ψ ( F e ) ∂ F e ∂ F e ∂ x i = − ∑ e V e P ( F e ) ∂ F e ∂ x i \begin{aligned} \bold{f}_{t,i} &= - \cfrac{\partial U}{\partial \bold{x}_i} \\ &= -\sum_e \cfrac{\partial U(e)}{\partial \bold{x}_i} \\ &= -\sum_e V_e \cfrac{\partial\Psi(\bold{F}_e)}{\partial\bold{F}_e} \cfrac{\partial\bold{F}_e}{\partial{\bold{x}_i}} \\ & = -\sum_e V_e \bold{P}(\bold{F}_e)\cfrac{\partial\bold{F}_e}{\partial{\bold{x}_i}} \end{aligned} ft,i=−∂xi∂U=−e∑∂xi∂U(e)=−e∑Ve∂Fe∂Ψ(Fe)∂xi∂Fe=−e∑VeP(Fe)∂xi∂Fe
@ti.kernel
def integrate():
for p in x:
...... # 碰撞检测
# 显式时间积分
v[p] = (v[p] + ((-x.grad[p] / node_mass) + ti.Vector([0, -10])) * dt) * math.exp(dt * -6)
x[p] += dt * v[p]
- 其中:
v t + 1 , i = ( v t , i + f t , i + m i g m i Δ t ) e − 6 Δ t \bold{v}_{t+1,i} = \left(\bold{v}_{t,i} + \cfrac{\bold{f}_{t,i}+m_i g}{m_i} \Delta t \right) e^{-6 \Delta t} vt+1,i=(vt,i+mift,i+migΔt)e−6Δt
6.4 隐式时间积分
[ I − Δ t 2 M − 1 ∂ f ∂ x ( x t ) ] v t + 1 = v t + Δ t M − 1 f ( x t ) \left[ \bold{I} - \Delta t^2 \bold{M}^{-1} \cfrac{\partial \bold{f}}{\partial \bold{x}}(\bold{x}_t) \right] \bold{v}_{t+1} = \bold{v}_t + \Delta t \bold{M}^{-1} \bold{f} (\bold{x}_t) [I−Δt2M−1∂x∂f(xt)]vt+1=vt+ΔtM−1f(xt)
力的微分计算:
∂ f ∂ x = − ∂ 2 Ψ ∂ x 2 \cfrac{\partial \bold{f}}{\partial \bold{x}} = - \cfrac{\partial^2 \bold{\Psi}}{\partial \bold{x}^2} ∂x∂f=−∂x2∂2Ψ