流体模拟渲染基础
- 前言
- 矢量微积分
- Naiver-Stokes偏微分方程组
- N-S方程的分步求解
- 对流算法
前言
本文主要参考文献《FLUID SIMULATION SIGGRAPH 2007 Course Notes》,结合我的理解单纯地讲述一下流体渲染的一些基础知识,本人水平有限,如有错误,欢迎指出。本文只是单纯针对流体模拟领域,可能一些地方不太严谨,但是对于虚拟模拟来说是可行的。即便如此,本文涉及到大量的数学方法。
一、矢量微积分
高等数学中太多数讨论的是一维的微积分,而矢量微积分则是一维微积分的高维扩展。矢量微积分的三个基础算子:梯度(符号为 ∇ ∇ ∇),散度(符号为 ∇ ⋅ ∇\cdot ∇⋅),旋度(符号为 ∇ × ∇\times ∇×),在此基础上流体力学中经常用到的还有拉普拉斯算子。
1、梯度(Gradient)
梯度实际上就是矢量的空间偏导数,且结果依然是一个矢量, 2 2 2维的梯度如下:
(1.1) ∇ f ( x , y ) = ( ∂ f ∂ x , ∂ f ∂ y ) ∇f(x,y)=(\frac{\partial f}{\partial x},\frac{\partial f}{\partial y}) \tag {1.1} ∇f(x,y)=(∂x∂f,∂y∂f)(1.1)
依此类推, 3 3 3维的梯度有如下形式:
(1.2) ∇ f ( x , y , z ) = ( ∂ f ∂ x , ∂ f ∂ y , ∂ f ∂ z ) ∇f(x,y,z)=(\frac{\partial f}{\partial x},\frac{\partial f}{\partial y},\frac{\partial f}{\partial z}) \tag {1.2} ∇f(x,y,z)=(∂x∂f,∂y∂f,∂z∂f)(1.2)
有时也会采用如下形式来表示梯度:
(1.3) ∇ f = ∂ f ∂ x ⃗ ∇f=\frac{\partial f}{\partial \vec x} \tag {1.3} ∇f=∂x∂f(1.3)
梯度通常用来近似计算函数值(实际上就是一维形式的推广):
(1.4) f ( x ⃗ + Δ x ⃗ ) ≈ f ( x ⃗ ) + ∇ f ( x ⃗ ) ⋅ Δ x ⃗ f(\vec x+\Delta \vec x)\approx f(\vec x)+∇f(\vec x)\cdot \Delta \vec x \tag {1.4} f(x+Δx)≈f(x)+∇f(x)⋅Δx(1.4)
同样的,多个函数的梯度就构成了一个矩阵:
(1.5) ∇ F ⃗ = ∇ ( f , g , h ) = ( ∂ f ∂ x ∂ f ∂ y ∂ f ∂ z ∂ g ∂ x ∂ g ∂ y ∂ g ∂ z ∂ h ∂ x ∂ h ∂ y ∂ h ∂ z ) = ( ∇ f ∇ g ∇ h ) ∇\vec F=∇(f,g,h)=\left( \begin{matrix} \frac{\partial f}{\partial x} & \frac{\partial f}{\partial y} & \frac{\partial f}{\partial z} \\ \frac{\partial g}{\partial x} & \frac{\partial g}{\partial y} & \frac{\partial g}{\partial z} \\ \frac{\partial h}{\partial x} & \frac{\partial h}{\partial y} & \frac{\partial h}{\partial z} \\ \end{matrix} \right) =\left( \begin{matrix}∇f\\ ∇g\\ ∇h\\ \end{matrix} \right) \tag {1.5} ∇F=∇(f,g,h)=⎝⎜⎛∂x∂f∂x∂g∂x∂h∂y∂f∂y∂g∂y∂h∂z∂f∂z∂g∂z∂h⎠⎟⎞=⎝⎛∇f∇g∇h⎠⎞(1.5)
2、散度(Divergence)
散度算子仅仅应用于向量场,它衡量在某一点出相应的矢量聚集或者发散程度,测量方向为径向,结果为标量。 2 2 2维、 3 3 3维形式的散度算子如下所示:
∇ ⋅ u ⃗ = ∇ ⋅ ( u , v ) = ∂ u ∂ x + ∂ v ∂ y ∇\cdot \vec u=∇\cdot (u,v)=\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y} ∇⋅u=∇⋅(u,v)=∂x∂u+∂y∂v
(1.6) ∇ ⋅ u ⃗ = ∇ ⋅ ( u , v , w ) = ∂ u ∂ x + ∂ v ∂ y + ∂ w ∂ z ∇\cdot \vec u=∇\cdot (u,v,w)=\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}+\frac{\partial w}{\partial z} \tag {1.6} ∇⋅u=∇⋅(u,v,w)=∂x∂u+∂y∂v+∂z∂w(1.6)
输入是矢量,而输出为标量。类比梯度,散度符号 ∇ ⋅ u ⃗ ∇\cdot \vec u ∇⋅u可以理解为梯度 ∇ ∇ ∇与矢量 u ⃗ \vec u u的点乘:
(1.7) ∇ ⋅ u ⃗ = ( ∂ ∂ x , ∂ ∂ y , ∂ ∂ z ) ⋅ ( u , v , w ) = ∂ ∂ x u + ∂ ∂ y v + ∂ ∂ z w ∇\cdot \vec u=(\frac{\partial}{\partial x},\frac{\partial}{\partial y},\frac{\partial}{\partial z})\cdot (u,v,w)=\frac{\partial}{\partial x}u+\frac{\partial}{\partial y}v+\frac{\partial}{\partial z}w \tag {1.7} ∇⋅u=(∂x∂,∂y∂,∂z∂)⋅(u,v,w)=∂x∂u+∂y∂v+∂z∂w(1.7)
若矢量场散度为 0 0 0,则称该矢量场无散度。
3、旋度(Curl)
旋度衡量围绕某一点的旋转速度,测量方向为切向,三维形式的旋度是一个向量:
(1.8) ∇ × u ⃗ = ∇ × ( u , v , w ) = ( ∂ w ∂ y − ∂ v ∂ z , ∂ u ∂ z − ∂ w ∂ x , ∂ v ∂ x − ∂ u ∂ y ) ∇\times \vec u=∇\times (u,v,w) =(\frac{\partial w}{\partial y}-\frac{\partial v}{\partial z}, \frac{\partial u}{\partial z}-\frac{\partial w}{\partial x}, \frac{\partial v}{\partial x}-\frac{\partial u}{\partial y}) \tag {1.8} ∇×u=∇×(u,v,w)=(∂y∂w−∂z∂v,∂z∂u−∂x∂w,∂x∂v−∂y∂u)(1.8)
倒推到 2 2 2维,我们取上式中的 w = 0 w=0 w=0,即矢量场为 ( u , v , 0 ) (u,v,0) (u,v,0), 2 2 2维向量场的旋度是一个标量:
(1.9) ∇ × u ⃗ = ∇ × ( u , v ) = ∂ v ∂ x − ∂ u ∂ y ∇\times \vec u=∇\times (u,v)=\frac{\partial v}{\partial x}-\frac{\partial u}{\partial y} \tag {1.9} ∇×u=∇×(u,v)=∂x∂v−∂y∂u(1.9)
同样地,旋度符号 ∇ × u ⃗ ∇\times \vec u ∇×u我们可以理解为梯度 ∇ ∇ ∇与矢量场 u ⃗ \vec u u的叉乘:
(1.10) ∇ × u ⃗ = ( ∂ ∂ x , ∂ ∂ y , ∂ ∂ z ) × ( u , v , w ) ∇\times \vec u= (\frac{\partial }{\partial x}, \frac{\partial }{\partial y}, \frac{\partial }{\partial z})\times(u,v,w) \tag {1.10} ∇×u=(∂x∂,∂y∂,∂z∂)×(u,v,w)(1.10)
若矢量场旋度为 0 0 0,则称该矢量场无旋度。
4、拉普拉斯算子(Laplacian)
拉普拉斯算子定义为梯度的散度,符号表示为 ∇ ⋅ ∇ ∇\cdot∇ ∇⋅∇,显然 ∇ ⋅ ∇\cdot ∇⋅是散度,而后面的 ∇ ∇ ∇则为梯度,故拉普拉斯算子即梯度的散度,是一个二阶微分算子。 2 2 2维、 3 3 3维形式分别如下:
∇ ⋅ ∇ f = ∂ 2 f ∂ x 2 + ∂ 2 f ∂ y 2 ∇\cdot∇f=\frac{\partial^2f}{\partial x^2}+\frac{\partial^2f}{\partial y^2} ∇⋅∇f=∂x2∂2f+∂y2∂2f
(1.11) ∇ ⋅ ∇ f = ∂ 2 f ∂ x 2 + ∂ 2 f ∂ y 2 + ∂ 2 f ∂ z 2 ∇\cdot∇f=\frac{\partial^2f}{\partial x^2}+\frac{\partial^2f}{\partial y^2}+\frac{\partial^2f}{\partial z^2} \tag {1.11} ∇⋅∇f=∂x2∂2f+∂y2∂2f+∂z2∂2f(1.11)
简言之,拉普拉斯算子定义如下:
(1.12) ∇ ⋅ ∇ f = Σ i = 1 n ∂ 2 f ∂ x i 2 ∇\cdot∇f=\Sigma_{i=1}^n\frac{\partial^2f}{\partial x_i^2} \tag {1.12} ∇⋅∇f=Σi=1n∂xi2∂2f(1.12)
偏微分方程 ∇ ⋅ ∇ f = 0 ∇\cdot ∇f=0 ∇⋅∇f=0被称为拉普拉斯方程;而如果右边为某个非 0 0 0常数,即 ∇ ⋅ ∇ f = q ∇\cdot ∇f=q ∇⋅∇f=q,我们称之为泊松方程。更一般地,如果梯度再乘上一个标量 a a a(如 1 / ρ 1/\rho 1/ρ),即 ∇ ⋅ ( a ∇ f ) = q ∇\cdot (a∇f)=q ∇⋅(a∇f)=q,我们依旧称之为泊松问题。
二、 N a i v e r − S t o k e s Naiver-Stokes Naiver−Stokes偏微分方程组
流体模拟器的构建主要围绕著名的不可压缩 N a v i e r − S t o k e s Navier-Stokes Navier−Stokes方程展开,它是一个流体力学领域的偏微分方程,方程形式如下:
(2.1) ∂ u ⃗ ∂ t + u ⃗ ⋅ ∇ u ⃗ + 1 ρ ∇ p = g ⃗ + ν ∇ ⋅ ∇ u ⃗ \frac{\partial \vec u}{\partial t}+\vec u\cdot ∇\vec u+\frac1\rho∇p=\vec g+\nu∇\cdot∇\vec u \tag {2.1} ∂t∂u+u⋅∇u+ρ1∇p=g+ν∇⋅∇u(2.1)
(2.2) ∇ ⋅ u ⃗ = 0 ∇\cdot\vec u=0 \tag {2.2} ∇⋅u=0(2.2)
这个方程组看起非常地复杂,接下来我们就把它剖析成一个个比较简单的部分来理解。
1、符号标记
我们有必要定义一些物理量的符号用以标记:
符号 u ⃗ \vec u u在流体力学中通常表示为流体的速度矢量,记 3 3 3维的速度矢量 u ⃗ = ( u , v , w ) \vec u=(u,v,w) u=(u,v,w);
希腊字符 ρ \rho ρ是流体的密度,对于水,该值大约为 1000 k g / m 3 1000kg/m^3 1000kg/m3,而空气则大约为 1.3 k g / m 3 1.3kg/m^3 1.3kg/m3;
字符 p p p代表压力,流体对任何物体施加的单位面积力;
字符 g ⃗ \vec g g则是我们熟悉的重力加速度,通常取 ( 0 , − 9.81 , 0 ) m / s 2 (0,-9.81,0)m/s^2 (0,−9.81,0)m/s2。我们约定 y y y轴向上,而 x x x轴和 z z z轴在水平面上。另外补充一点,我们把其他的一些类似的力都累加到 g ⃗ \vec g g上,也就是我们统一用 g ⃗ \vec g g表示所有类似力之和,这类力我们称之为体积力(称之为体积力是因为它们的力是作用到整个流体而不只是流体的表面);
希腊字符 ν \nu ν是流体的运动粘度,它衡量流体的黏滞性。糖蜜之类的流体有非常高的粘度,而像酒精之类的流体有很低的粘度;
其它一些矢量微积分的符号算子前面已经提到过,不再赘述。
2、动量方程
偏微分方程 ( 2.1 ) (2.1) (2.