从我们呼吸的空气到覆盖地球三分之二的海洋,流体在我们的身边随处可见,是我们所知道的一些最美丽和最令人印象深刻的现象的核心。从水的飞溅,到火焰和烟雾的旋转,流体已经成为计算机图形学的一个重要组成部分。这本书旨在涵盖模拟这些动画效果的基本知识。让我们来看看控制它们运动的基本方程。
动画中大多数有趣的流体流动都是由著名的incompressible Navier-Stokes方程控制的,这是一组贯穿整个流体的偏微分方程。方程通常被写为:
字母\(\vec{u}\)通常在流体力学中用来表示流体的速度。至于为什么不使用\(\vec{v}\)表示,这很去解释,但它符合另一个有用的约定,称为3D速度的三个分量\((u,v,w)\),就像位置x的三个分量通常被称为\((x,y,z)\)一样。
希腊语字母ρ代表流体的密度。水的密度大约是1000\(kg/m^3\),而对于一般海平面条件下的空气,这大约是1.3\(kg/m^3\),比例约为700:1。
值得强调的是,我坚持使用真实的单位(如米、公斤等):长期的经验表明,以公制单位隐式地保持求解器中的所有量是很值得的,而不是仅仅设置为任意值。当开始编写一个新的求解器时,只需为物理量(如密度)填充无单位值(如1),或者从表达式中完全删除它们,无论您是从a quick-and-dirty just-make-it-work的观点或更具有数学基础的无量纲化原理进行操作。但是,当模拟看起来不太正确、需要调整大小或以其他方式调整时,这通常会再次困扰您,而这些非物理参数中的哪一个需要调整无法明确得知。我们将在整本书中讨论这一点在算法设计中的影响。
字母\(p\)代表压力,即流体对任何东西施加的单位面积的力。
字母\(\vec{g}\)是重力加速度,通常取值为\((0,-9.81,0)m/s^2\)。在本书中,我们将约定y轴垂直向上,x轴和z轴是水平的。同时应该补充一点,在动画中,我们可以添加额外的控制加速(使流体以某种期望的方式运行)—— 我们将把所有这些归类为一个符号\(\vec{g}\)。更一般地说,人们将这些称为body forces,因为它们作用于整个流体,而不仅仅是表面。
字母\(ν\)在技术上称为运动粘度(kinematic viscosity)。它测量流体的粘度。像糖蜜等流体具有高粘度,而汞等流体具有低粘度:它用于衡量流体在流动时抵抗变形的程度(或者更直观地说,是搅拌的难度)。
第一个微分方程(1.1)实际上是三合一的向量方程,称为动量方程。这经典的牛顿方程\(\vec{F} = m\vec{a}\)的变相。 它告诉我们流体如何由于作用在其上的力而加速。 在进入第二个微分方程(1.2)之前,我们将尝试将其分解,这称为不可压缩条件。它告诉我们,流体是如何通过作用在它身上的力而加速的。在进入第二个微分方程(1.2)之前,我们将尝试对其进行分解,这被称为不可压缩性条件(incompressibility condition)。
让我们首先想象我们使用粒子系统模拟流体(在书的后面,我们将实际使用它作为一个实用的方法,但现在让我们把它作为一个思维实验)。每一个粒子都可能代表一小块流体。它的质量是\(m\),体积是\(V\),速度是\(\vec{u}\)。为了及时整合系统,我们需要找出作用在每个粒子上的力是什么:\(\vec{F} = m\vec{a}\)告诉我们粒子是如何加速的,我们从中得到它的运动。我们将以稍微奇怪的符号来写粒子的加速度(与上面的动量方程相关):
大D导数表示法被称为物质导数(material derivative 稍后会详细介绍)。牛顿定律现在表示为:
\[m\frac{D\vec{u}}{Dt} = \vec{F} \] 那么作用在粒子上的力是什么?最简单的当然是重力:\(m\vec{g}\)。然而,当我们考虑流体的其余部分如何施加力时,它就变得有趣了:即粒子如何与附近的其他粒子相互作用。
第一个流体力是压力。高压区推动低压区。请注意,我们真正关心的是粒子上的净外力:例如,如果压力在每个方向上都相等,那么净外力将为零,并且不会因压力而产生加速度。我们只有在不平衡时才能看到对流体粒子的影响,即粒子一侧的压力高于另一侧,导致力从高压指向低压。在附录中,我们展示了如何严格推导这一点,但现在让我们指出,测量粒子位置压力不平衡的最简单方法是仅采用压力的负梯度:\(-∇p\)。(回想一下微积分,梯度是“最陡峭的上升”方向,因此负梯度从高压区域指向低压区域。)我们需要将它整合到我们的流体体积上以获得压力。作为一个简单的近似值,我们将乘以体积\(V\)。但你可能会问自己,压力是什么?我们将跳过这一点,稍后我们讨论不可压缩性,但现在你可以认为它是保持流体体积恒定所需要的一切。
另一个流体力是粘度。粘性流体试图抵抗变形。稍后我们将更深入地推导出它,但现在让我们直观地将其表示为一种力,它试图使我们的粒子以附近粒子的平均速度移动,即,试图最小化相邻流体之间的速度差异。您可能还记得,在图像处理、数字几何处理、扩散或散热物理或许多其他领域中,衡量一个量与其周围平均值之间的距离的微分算子是拉普拉斯算子\(∇·∇\)。(现在是提一下附录中向量微积分的快速回顾的好时机,包括像拉普拉斯算子这样的微分算子。)一旦我们将它整合到blob的体积上,这将提供我们的粘性力。我们将使用动态粘度系数,用希腊字母\(µ\)表示(动态意味着我们从中获得了力;之前的运动粘度用于获得加速度)。我要在这里指出,在液体表面附近(在blob周围没有完整的邻域)以及具有可变粘度的流体,这个术语最终会稍微复杂一些;有关详细信息,请参阅第10章。
总体而言,一团液体的移动方式如下:
显然,当我们用有限的少量粒子来近似流体时,我们会导致错误。然后,当我们的粒子数量趋于无穷大并且每个blob的大小趋于零时,我们将达到极限。当然,这显然会产生另一种错误,因为真正的流体实际上是由(非常大的)有限数量的分子组成的!但是这个极限情况,我们称之为连续模型(continuum model),具有简洁的数学表达和独立于确切的blob数量的优点,并且已通过实验证明在大量场景中与现实非常接近。但是,采用连续极限情况确实给我们的粒子方程带来了问题,因为粒子的质量\(m\)和体积\(V\)必须为零,我们就没有留下任何意义的东西。我们可以通过首先将方程除以体积,然后取极限来解决这个问题。记住\(m/V\)只是流体密度\(ρ\),我们可以得到:
\[ρ\frac{D\vec{u}}{Dt} = ρ\vec{g}-∇p+µ∇·∇\vec{u} \]看着眼熟吗?我们将其除以密度并重新排列各项可以获得:
\[\frac{D\vec{u}}{Dt}+\frac{1}{ρ}∇p= \vec{g}+\frac{µ}{ρ}∇·∇\vec{u} \]为了进一步简化,我们将运动粘度定义为\(ν = µ/ρ\)可以得到:
\[\frac{D\vec{u}}{Dt}+\frac{1}{ρ}∇p= \vec{g}+ν∇·∇\vec{u} \]我们几乎回到了动量方程上!事实上,这种使用材料导数\(D/Dt\)的形式对我们在计算机图形学中更为重要,并将指导我们以数值方式求解方程。但我们仍然想了解材料导数是什么和如何将它与传统的动量方程联系起来。为此,我们需要了解拉格朗日观点和欧拉观点之间的区别。
当我们考虑运动的连续体(如流体或可变形固体)时,有两种方法可以跟踪这种运动:拉格朗日观点和欧拉观点。