UnrealEngine 4/5 学习笔记:Niagara 流体模拟(2)
对流算法将对流数值计算的算法记为 q^{n+1}=advect(\vec u,\Delta t,q^{n}) ,其中
\vec u:在 MAC 网格上的离散化的速度场,必须保证是无散度的,否则模拟结果会出现一些奇怪的失真现象;
\Delta t :时间步长;
q^{n} :当前的物理场(如流体的密度、速度等);
q^{n+1} :经过对流后得到的新物理场;
半拉格朗日对流算法(Semi-Lagrangian Advection)
公式(7)写成偏微分的形式为 \frac{\partial q}{\partial t}+u\frac{\partial q}{\partial x}=0 ,如果分别采用前向欧拉差分法计算对时间的偏导和中心差分法计算对空间的偏导,可以得到
\frac{q_{i}^{n+1}-q_{i}^{n}}{\Delta t}+u_{i}^{n}\frac{q_{i+1}^{n}-q_{i-1}^{n}}{2\Delta x}=0 \\ 转换为 q_{i}^{n+1} 为计算目标的显式公式为
q_{i}^{n+1}=q_{i}^{n}-\Delta tu_{i}^{n}\frac{q_{i+1}^{n}-q_{i-1}^{n}}{2\Delta x} \\ 前向欧拉法是无条件不稳定的空间离散方法,时间上会出现累计误差发散,空间上会出现零空间问题,导致模拟效果中出现许多奇怪的高频摆动和震荡。研究者们提出了一个新的方法——半拉格朗日法。假设我们的目标是求解网格点 \vec{x}_{G} 的在第 n+1 个时间步时关于物理量 q 的新值,记为 q^{n+1}_G 。在拉格朗日的视角下,我们可以寻找在第 n+1 时间步之前,是空间中的哪一个点上的流体粒子在速度场 \vec u 的作用下“流向”了 \vec{x}_{G} ,我们记这个粒子在第 n 个时间步时的网格位置为 \vec x_P ,则第 n+1 个时间步时 \vec x_G 的 q^{n+1}_G 即为第 n 个时间步时 \vec x_P 的 q^{n}_P 。
半拉格朗日对流法的第一步就是要找出 \vec{x}_{P} ,为此我们根据 \vec{x}_{G} 做反向的追踪,粒子位置对时间的导数就是速度场: \frac{d\vec{x}}{dt}=\vec{u}(\vec{x}) 。经过一个时间步长 \Delta{t} 之后,粒子由 \vec{x}_{P} 移动到 \vec{x}_{G} ,为了得到 \vec{x}_{P} ,最简单的方法就是采用前向欧拉法进行倒推: \vec{x}_{P}=\vec{x}_{G}-\Delta{t}\vec{u}(\vec{x}_{G}) ,前向欧拉法只有一阶的精度,为了在不改变\Delta t 的情况下提高精度,采用二阶的龙格库塔法(Runge-Kutta method)如下所示:
\vec{x}_{mid}=\vec{x}_{G}-\frac{1}{2}\Delta{t}\vec{u}(\vec{x}_{G}) \\ \vec{x}_{P}=\vec{x}_{G}-\Delta{t}\vec{u}(\vec{x}_{mid}) \\
三维模拟通常采用三线性插值,而二维的则采用双线性插值:
q_{G}^{n+1}=interpolate(q_{n},\vec x_{p}) \\
边界情况
如果 \vec{x}_{P} 处于边界之外,一种方法是将这个外部流入的物理量作为返回值即可,还有一种方法是根据边界上的最近点外推出所求得物理量。
时间步长大小
\Delta{t}\leq{C}\frac{\Delta{x}}{\left| \vec{u} \right|} \\ 其中, C 是一个小的整数常量。
数值耗散
q=w_{-1}(t)q_{j-1}+w_{0}(t)q_{j}+w_{1}(t)q_{j+1}+w_{2}(t)q_{j+2} \\
页:
[1]