流体模拟《Fluid Simulation for Computer Graphics》——从入门到放弃(五)综合与理解

这个系列前四篇看的时候啥都不会,很多很多都不懂,只是非常粗浅的就表面翻译进行了笨拙的理解,尤其是流体算法的部分,,what a mess,,,所以需要停下来做一个综合的整理和回顾,有些内容是直接摘录前几篇的说明,有些内容是增加了更系统更准确一点点的说明。

有错请尽管联系我,我理解错误什么的,,那可太正常了QAQ,,感激不尽:D


1 Navier-Stokes方程组

从不可压缩的Navier-Stokes方程组开始。

1.1 动量方程

方程组里的第一个方程也叫动量方程,其本质就是牛二定律F=ma

  • u表示速度场中的速度,在三维空间中的分解形式为(u,v,w)
  • ρ表示流体密度
  • p表示压力
  • g表示重力加速度
  • ν表示运动粘度(描述了流体在运动时抵抗变形的能力)

a=Du/Dt,根据牛二定律有m*Du/Dt=F

这个合力F包括内力/外力。

1.2 内力

内力是影响流体粒子与附近其他粒子相互作用的两种流体力,一种是压力,一种是由流体粘度引起的力。

首先是压力,因为外力施加在流体上的时候,外力并没有马上传导到整个流体,而是在受力区附近产生高压,然后高压区推向低压区(所以要取负压力梯度:−∇p),这个值乘以体积V(这里研究的都是不可压缩流体,所以V为常量),就是别的粒子施加给它的压力和−V∇p

第二种内力由流体粘度引起,粘性流体试图抵抗变形,这种内力试图使这个粒子以周围粒子的平均速度运动,即邻近粒子之间的速度差异最小化。这种粘度造成了动量的扩散,也可以称为数值粘性/人工粘度/假扩散。比如下图这样的随着移动边界模糊。

这个粘度的力是这么计算的,用拉普拉斯算子∇·∇衡量一个量与周围平均数的差,µ是一个粘滞系数,然后乘以体积V,就得到了这一项:Vµ∇·∇u

1.3 外力

外力就是均匀的作用于整个流体的力,所以主要就是考虑重力,重力似乎没啥好说的:mg

1.4 进一步推导

把上面这些作用力都加起来,就可以得到由牛二而来的这个等式:m*Du/Dt=mg−V∇p+Vµ∇·∇u

Du/Dt是一个全微分项,根据链式法则,可以得到:,用Du/Dt替换F=ma里的加速度a,就可以得到动量方程。

1.5 不可压缩

真实的流体,甚至是水这样的液体体积也会改变,但通常变化不大。这样的研究对象称为可压缩流,研究起来很复杂,而宏观上的影响又不大,所以就应用层面来说没必要考虑。

所以在流体模拟中一个前提条件是流体体积变化率为0(不可压缩),对应着Navier-Stokes方程组里的第二个式子:∇·u=0

具体推导过程如下:

对于流体中一个确定的区域Ω,边界是∂Ω,不可压缩流体,所以这个区域Ω体积变化率为0,即:)。由高斯公式有:,也就是对于任意Ω都不可压缩,所以整个流体的体积变化率也是0。然后就可以得到∇·u=0

1.6 Material derivative

最后,再补上一个叫物质导数的东西,根据链式法则得到下式。

其中q表示通过流体传输的模拟量之一(比如温度),t表示时间。∂q/∂t就表示了q在固定点变化的速度。第二项u⋅∇q纠正多少变化是由于流体流过的差异。

2 Simulation Viewpoints(模拟视角)

跟踪运动的流体或可变形固体(连续体)的运动有两种途径:欧拉视角拉格朗日视角

2.1 拉格朗日视角

把流体看成许多点,每个点都标记为一个单独的粒子,有各自的位置x,速度u

2.2 欧拉视角

在空间上保持固定点,关注这些空间中的固定点,查看这些点所处的流体(速度/温度等)如何随时间变化。

在这种视角下,更容易在固定的欧拉网格上对一些流体量的空间导数进行分析及对这些导数进行数值逼近(而拉格朗日视角下是一堆位置任意的点集,相对获取空间导数难很多)。

2.3 两种视角的关系

这两种视角之间还是可以联系在一起的,联系上面的1.6节的式子。

1.6节表明从拉格朗日视角出发探究每个粒子q值的变化情况,链式法则推导可知每个粒子对应的物理量q在流体运动过程中只是被传递,而不会发生改变。其中等式右边的∂q/∂t就表示了q在固定点变化的速度,也就是在欧拉视角下,这个q是个不变量。

3 流体算法

3.1 拉格朗日方法

拉格朗日方法显然是在拉格朗日视角下的方法,把流体模拟成很多离散的粒子,只要能描述任意一个粒子的运动状况(位置随时间变化的规律等),那么就可以获得整个流体的运动状况。

还是从不可压缩的Navier-Stokes方程组开始:

接下来就是要求解这个方程,引入SPH方法(光滑粒子流体动力学方法)。

首先引入一个光滑核的概念,每个粒子属性都会扩散到它周围,影响周围的粒子,越远那么影响就越小,这种随着距离衰减的函数就是光滑核函数

然后考虑流体中某点r(点和粒子区分,某点不一定有粒子),以这个点为圆心,光滑核半径h为半径的圆中有其他粒子r0,r1,…,rj,那么点r处就有一个关于属性A的累加公式:,其中Aj/mj/ρj代表其他粒子的属性A/质量/密度,W是光滑核函数(其实所有粒子质量一样,都是m)。

再考虑流体中一个位置为ri的点,由牛二推导的式子可得:),其中u(ri)/p(ri)/ρ(ri)分别代表此处的速度/压力/密度。

接下来就是每一项分开求解的过程了,比如说要求密度,那么把关于属性A的累加公式里的A换成密度ρ,然后取一个光滑核函数,经过推导得到密度计算公式。(此处太过复杂,文末给出的那个“SPH算法系列”中大神博主已给出了非常详细的推导)

最后,虽然把流体想象成很多粒子的想法很合理,但这是一种不规则离散,在实际中很难计算(2.2节中做了说明),需要对无数多的粒子进行跟踪;而核函数的选取也很重要;并且【本系列第四篇第8节】里的用于表面追踪的水平集方法在构造光滑液体表面的问题上大大优于粒子提取出的表面;但是这种方法不需要网格划分。如果是不规则的边界条件或者多种流体混合作用,就还是用拉格朗日方法。

更常用的是接下来的欧拉方法。

3.2 欧拉方法

欧拉方法就是欧拉视角下的方法。欧拉方法用网格法对Navier-Stokes方程组进行离散化,研究空间中的固定点,用有限差分求解。

还是从不可压缩的Navier-Stokes方程组开始:

对于任意一个确定的位置x和确定的时刻t,用矢量场v(x,y,z,t)描述流体的速度,标量场p(x,t,z,t)描述流体的压强。

网格法有两种,最流行的是MAC法(标记网格法),在【本系列第二篇文章4.2节】中描述了。看下图就很直观:压力定义在格子中心,r方向速度分量定义在格子左右边界中点,y方向速度分量定义在格子上下边界中点,三维同理。

还有一种方法是把所有数据都标记在网格节点上。

然后回到欧拉方法,有限差分求解三维Navier-Stokes方程组,半拉格朗日方法求解动量方程的对流项【本系列第三篇文章的第5节】,隐式方法求解粘度和压强。具体步骤是:初始为u0(x),先增加合力得u1(x),然后求解对流项得u2(x),再求解粘度得u3(x),最后投影得u4(x)

欧拉方法和用于表面追踪的水平集方法结合,可以达到很好的效果。


最后放上所有厉害资料集合,非常特别宇宙无敌感谢!!!

打赏
  • 版权声明: 本博客所有文章除特别声明外,著作权归作者所有。转载请注明出处!
  • Copyrights © 2018-2020 LeFlacon

奶茶一杯 快乐起飞

支付宝
微信