Fluid Simulation for Computer Graphics - 第一章(The Equations of Fluids)学习

  从我们呼吸的空气到覆盖地球三分之二的海洋,流体在我们的身边随处可见,是我们所知道的一些最美丽和最令人印象深刻的现象的核心。从水的飞溅,到火焰和烟雾的旋转,流体已经成为计算机图形学的一个重要组成部分。这本书旨在涵盖模拟这些动画效果的基本知识。让我们来看看控制它们运动的基本方程。
  动画中大多数有趣的流体流动都是由著名的incompressible Navier-Stokes方程控制的,这是一组贯穿整个流体的偏微分方程。方程通常被写为:

图片名称
  乍一看,这些问题可能会显得很复杂。我们很快就会把它们分解成易于理解的部分(在附录B中提供了一个更严格的解释),但首先让我们开始定义每个符号的含义。

1.1 符号定义

  字母u通常在流体力学中用来表示流体的速度。至于为什么不使用v表示,这很去解释,但它符合另一个有用的约定,称为3D速度的三个分量(u,v,w),就像位置x的三个分量通常被称为(x,y,z)一样。
  希腊语字母ρ代表流体的密度。水的密度大约是1000kg/m3,而对于一般海平面条件下的空气,这大约是1.3kg/m3,比例约为700:1。
  值得强调的是,我坚持使用真实的单位(如米、公斤等):长期的经验表明,以公制单位隐式地保持求解器中的所有量是很值得的,而不是仅仅设置为任意值。当开始编写一个新的求解器时,只需为物理量(如密度)填充无单位值(如1),或者从表达式中完全删除它们,无论您是从a quick-and-dirty just-make-it-work的观点或更具有数学基础的无量纲化原理进行操作。但是,当模拟看起来不太正确、需要调整大小或以其他方式调整时,这通常会再次困扰您,而这些非物理参数中的哪一个需要调整无法明确得知。我们将在整本书中讨论这一点在算法设计中的影响。
  字母p代表压力,即流体对任何东西施加的单位面积的力。
  字母g重力加速度,通常取值为(0,9.81,0)m/s2。在本书中,我们将约定y轴垂直向上,x轴和z轴是水平的。同时应该补充一点,在动画中,我们可以添加额外的控制加速(使流体以某种期望的方式运行)—— 我们将把所有这些归类为一个符号g。更一般地说,人们将这些称为body forces,因为它们作用于整个流体,而不仅仅是表面。
  字母ν在技术上称为运动粘度(kinematic viscosity)。它测量流体的粘度。像糖蜜等流体具有高粘度,而汞等流体具有低粘度:它用于衡量流体在流动时抵抗变形的程度(或者更直观地说,是搅拌的难度)。

1.2 动量方程

  第一个微分方程(1.1)实际上是三合一的向量方程,称为动量方程。这经典的牛顿方程F=ma的变相。 它告诉我们流体如何由于作用在其上的力而加速。 在进入第二个微分方程(1.2)之前,我们将尝试将其分解,这称为不可压缩条件。它告诉我们,流体是如何通过作用在它身上的力而加速的。在进入第二个微分方程(1.2)之前,我们将尝试对其进行分解,这被称为不可压缩性条件(incompressibility condition)。
  让我们首先想象我们使用粒子系统模拟流体(在书的后面,我们将实际使用它作为一个实用的方法,但现在让我们把它作为一个思维实验)。每一个粒子都可能代表一小块流体。它的质量是m,体积是V,速度是u。为了及时整合系统,我们需要找出作用在每个粒子上的力是什么:F=ma告诉我们粒子是如何加速的,我们从中得到它的运动。我们将以稍微奇怪的符号来写粒子的加速度(与上面的动量方程相关):

aDuDt

  大D导数表示法被称为物质导数(material derivative 稍后会详细介绍)。牛顿定律现在表示为:

mDuDt=F

  那么作用在粒子上的力是什么?最简单的当然是重力:mg。然而,当我们考虑流体的其余部分如何施加力时,它就变得有趣了:即粒子如何与附近的其他粒子相互作用。
  第一个流体力是压力。高压区推动低压区。请注意,我们真正关心的是粒子上的净外力:例如,如果压力在每个方向上都相等,那么净外力将为零,并且不会因压力而产生加速度。我们只有在不平衡时才能看到对流体粒子的影响,即粒子一侧的压力高于另一侧,导致力从高压指向低压。在附录中,我们展示了如何严格推导这一点,但现在让我们指出,测量粒子位置压力不平衡的最简单方法是仅采用压力的负梯度:p。(回想一下微积分,梯度是“最陡峭的上升”方向,因此负梯度从高压区域指向低压区域。)我们需要将它整合到我们的流体体积上以获得压力。作为一个简单的近似值,我们将乘以体积V。但你可能会问自己,压力是什么?我们将跳过这一点,稍后我们讨论不可压缩性,但现在你可以认为它是保持流体体积恒定所需要的一切。
  另一个流体力是粘度。粘性流体试图抵抗变形。稍后我们将更深入地推导出它,但现在让我们直观地将其表示为一种力,它试图使我们的粒子以附近粒子的平均速度移动,即,试图最小化相邻流体之间的速度差异。您可能还记得,在图像处理、数字几何处理、扩散或散热物理或许多其他领域中,衡量一个量与其周围平均值之间的距离的微分算子是拉普拉斯算子·。(现在是提一下附录中向量微积分的快速回顾的好时机,包括像拉普拉斯算子这样的微分算子。)一旦我们将它整合到blob的体积上,这将提供我们的粘性力。我们将使用动态粘度系数,用希腊字母µ表示(动态意味着我们从中获得了力;之前的运动粘度用于获得加速度)。我要在这里指出,在液体表面附近(在blob周围没有完整的邻域)以及具有可变粘度的流体,这个术语最终会稍微复杂一些;有关详细信息,请参阅第10章
  总体而言,一团液体的移动方式如下:

mDuDt=mgVp+Vµ·u

  显然,当我们用有限的少量粒子来近似流体时,我们会导致错误。然后,当我们的粒子数量趋于无穷大并且每个blob的大小趋于零时,我们将达到极限。当然,这显然会产生另一种错误,因为真正的流体实际上是由(非常大的)有限数量的分子组成的!但是这个极限情况,我们称之为连续模型(continuum model),具有简洁的数学表达和独立于确切的blob数量的优点,并且已通过实验证明在大量场景中与现实非常接近。但是,采用连续极限情况确实给我们的粒子方程带来了问题,因为粒子的质量m和体积V必须为零,我们就没有留下任何意义的东西。我们可以通过首先将方程除以体积,然后取极限来解决这个问题。记住m/V只是流体密度ρ,我们可以得到:

ρDuDt=ρgp+µ·u

  看着眼熟吗?我们将其除以密度并重新排列各项可以获得:

DuDt+1ρp=g+µρ·u

  为了进一步简化,我们将运动粘度定义为ν=µ/ρ可以得到:

DuDt+1ρp=g+ν·u

  我们几乎回到了动量方程上!事实上,这种使用材料导数D/Dt的形式对我们在计算机图形学中更为重要,并将指导我们以数值方式求解方程。但我们仍然想了解材料导数是什么和如何将它与传统的动量方程联系起来。为此,我们需要了解拉格朗日观点欧拉观点之间的区别。

1.3 拉格朗日与欧拉观点

  当我们考虑运动的连续体(如流体或可变形固体)时,有两种方法可以跟踪这种运动:拉格朗日观点和欧拉观点。
  拉格朗日方法,以你最熟悉的法国数学家拉格朗日的名字命名。它就像对待粒子系统一样对待连续体。流体或固体中的每个点都被标记为一个单独的粒子,位置为x和速度为u。你甚至可以把每个粒子看作是液体中的一个分子。这里没什么特别的!固体几乎总是以拉格朗日的方式模拟,其中一组离散的粒子通常连接在一个网格中。
  欧拉方法,以瑞士数学家欧拉的名字命名,它采用了一种通常用于流体的不同策略。我们不是跟踪每个粒子,而是观察空间中的固定点,看看在这这些点上流体量的测量,如密度、速度、温度等,是如何随时间变化的。流体可能流过这些点,造成了一种变化:例如,当温暖的流体经过冷的流体时,空间中固定点的温度就会降低——即使流体中任何单个粒子的温度不会变化!此外,流体变量可能在流体中发生变化,从而导致另一种可能在固定点测量的变化:例如,在空间中固定点测量的温度会随着各地的流体冷却而降低。
  考虑这两种观点的一种方法是做一个天气报告。在拉格朗日的观点中,你在一个随风漂流的气球中,测量在你身边流动的空气的压力、温度和湿度等等。从欧拉的观点来看,你被困在地上,测量流过的空气的压力、温度和湿度等。这两种测量都可以创建一个条件如何变化的图,但这些图可以是完全不同的,因为它们正在以根本不同的方式测量变化率。
  数值上,拉格朗日观点对应于一个粒子系统,有或没有连接粒子的网格,欧拉观点对应于使用一个固定的网格,即使流体流过它,它在空间中也不会改变。
  似乎欧拉方法过于复杂:为什么不坚持使用拉格朗日粒子系统呢?事实上,确实有一些方案,例如涡旋方法和平滑粒子流体动力学(SPH)。然而,即使这些都依赖于欧拉派生的流体力方程,在本书中,我们将在很大程度上坚持使用欧拉方法,原因如下:

  • 使用欧拉观点中的压力梯度和粘度项等空间导数进行分析工作更容易。
  • 在固定的欧拉网格上对这些空间导数进行数值近似比在任意移动的粒子云上要容易得多。

  连接这两种观点的关键是物质导数。我们将从拉格朗日描述开始:有位置为x和速度为u的粒子。让我们看一个我们称之为q 的通用量:每个粒子都有一个q值。(q可能是密度、速度、温度或许多其他东西。)特别是,函数q(t,x)告诉我们恰好在位置x的粒子在时间tq值:这是一个欧拉变量,因为它是空间的函数,而不是粒子的函数。那么对于位置由x(t)作为时间函数给出的粒子,q的变化速度有多快,即拉格朗日问题?只需取总导数(a.k.a.链式法则):

ddtq(t,x(t))=qt+q·dxdt=qt+q·uDqDt

  这就是物质导数!
  让我们回顾一下物质导数中的两个术语。第一个是q/t,这就是q在空间中固定点的变化速度,这是欧拉测量。第二项q·u正在纠正这种变化有多少仅仅是由于流过的流体的差异(例如,温度变化是因为热空气被冷空气取代,而不是因为任何分子的温度正在改变)。
  为了完整起见,让我们写出物质导数,以及所有的偏导数:

DqDt=qt+uqx+vqy+wqz

  显然,在2D中,我们可以去掉wz项。
请注意,我一直在谈论数量、分子或粒子如何随着速度场u移动。这被称为平流(有时是对流传输;它们的含义相同)。平流方程只是一个使用物质导数的方程,最简单的是将其设置为零:

DqDt=0i.e.,qt+u·q=0

  在拉格朗日观点中这只是意味着数量在移动,并没有改变。

1.3.1 一个例子

  为了解决这个问题,让我们通过一个一维的例子来讲解。我们将使用T代替q来表示温度。我们会说,在某一时刻,温度曲线是:

T(x)=10x;

  也就是说,它在原点结冰,随着我们进一步向右时变得更暖,在x=10处的温度为100时。现在假设有一股速度为c的稳定风在吹,即,流体速度在任意点都是c

u=c.

  我们假设每个空气粒子的温度没有变化——它们只是在移动。因此,用拉格朗日观点衡量事物的物质导数表示变化为零:

DTDt=0.

  如果我们扩展它,我们可以得到:

Tt+u·T=0Tt+10·c=0Tt=10·c

  也就是说,在空间的固定点,温度以-10c的速率变化。如果风已经停止,c=0,没有任何改变。如果风以速度c=1方向向右吹动,在一个固定点的温度将以−10的速率下降。如果风速以c=−2的速度向左加速,则在一个固定点上的温度将以20的速度上升。所以即使拉格朗日导数是零,在这种情况下,欧拉导数可以是任何东西,这取决于流动的速度和运动方向。

1.3.2 Advecting Vector Quantities

  一个常见的混淆点是,当物质导数应用于矢量量时意味着什么,比如RGB的颜色,或者最令人困惑的是,速度场u本身。简单的答案是:分别对待每个组件。让我们写出颜色向量 C=(R,G,B)的物质导数:

DCDt=[DR/DtDG/DtDB/Dt]=[R/t+u·RG/t+u·GB/t+u·B]=Ct+u·C.

  因此,尽管符号u·C可能没有严格意义上的意义(向量的梯度是矩阵吗?向量与矩阵的点积是多少?)但不难弄清楚我们是否把向量分割成标量分量。
  让我们对速度本身做同样的事情,除了u出现在两个地方之外,它实际上并没有什么不同,作为流体移动的速度场和平流的流体量。人们有时称此为自我平流,以强调速度出现在两个不同的角色中,但公式的工作原理与颜色完全相同。所以仅仅通过复制和粘贴,这是速度u=(uvw)的平流:

DrDt=[Du/DtDv/DtDw/Dt]=[u/t+u·uv/t+u·vw/t+u·w]=ut+u·u,

  或者,如果您想深入了解偏导数的具体细节,

DuDt=[ut+uux+vuy+wuzvt+uvx+vvy+wvzwt+uwx+vwy+wwz].

  如果您想更进一步,将矩阵数量平流化,也没有什么不同:只需分别处理每个组件。

1.4 不可压缩性

  真正的流体,甚至像水这样的液体,它们的体积都会被改变。事实上,这正是声波的本质:流体体积的扰动,进而导致密度和压力的扰动。你可能曾经被教导过,液体和气体之间的区别在于气体会改变它们的体积,而液体不会,但事实并非如此:否则你将无法在水下听到声音!
  然而,关键是流体通常不会改变它们的体积。 即使有一个非常强大的泵,也几乎不可能改变水的体积。即使空气也不会改变它的体积,除非你把它放在泵里,或者处理非常极端的情况,比如音爆和冲击波。对流体在这些情况下如何表现的研究通常称为可压缩流动。模拟起来既复杂又昂贵,而且除了声学之外,它并没有明显地进入日常生活——甚至声波在体积中也是如此微小的扰动,对流体在宏观水平上的移动方式影响如此之小(水的晃,烟雾滚滚等),它们几乎与动画无关。
  这意味着在动画中,我们通常可以将所有流体(包括液体和气体)视为不可压缩的,这意味着它们的体积不会改变。这在数学上意味着什么?附录B中有更严格的解释,但我们现在可以草拟一个简短的论点。
  选择任意一块液体来观察一段时间。 我们称这个体积为Ω和它的边界面为Ω(这些是传统符号)。我们可以通过积分边界周围速度的法向分量来测量这块流体的体积变化的速度:

ddtvolume(Ω)=Ωu·n^

  如果曲面切向移动,则不会影响体积;如果它在向外/向内移动,体积会按比例增加/减少。对于不可压缩的流体,体积最好保持不变,即这个变化率必须为零:

Ωu·n^=0.

  现在我们可以使用散度定理将其变为体积积分。基本上,这是微积分基本定理的一个多维版本:如果你积分一个函数的导数,你会得到在积分范围内评估的原始函数(如果你需要复习你的向量微积分,请参阅附录A进行回顾)。在这种情况下,我们得到:

Ω·u=0.

现在,神奇的部分来了:这个方程应该适用于任何Ω的选择,任何流体区域。唯一独立于积分区域积分为零的连续函数是零本身。因此被积函数在任何地方都必须为零:

·u=0.

  这是不可压缩条件,不可压缩 Navier-Stokes 方程的另一部分。
  出于显而易见的原因,满足不可压缩条件的矢量场被称为无散度。模拟不可压缩流体的棘手部分之一是确保速度场保持无发散。这就是压力的来源。
  考虑压力的一种方法是,它正是保持速度无散度所需的力。如果您熟悉约束动力学,您可以将不可压缩条件视为一个约束,而将压力场视为满足该约束所需的拉格朗日乘子,该约束遵循零虚功原则。如果你不是,别担心。让我们确切地推导出压力必须是多少。
  压力只出现在动量方程中,我们想以某种方式将其与速度的散度联系起来。因此,让我们取动量方程两边的散度:

·ut+·(u·u)+·1ρp=·(g+v·u).(1.3)

  这与我们的数值模拟并不完全相关,但值得一看,因为当我们离散化时,我们将经历几乎完全相同的步骤,从查看体积变化的速度到压力方程。

1.5 Dropping Viscosity

  在某些情况下,粘性力非常重要:例如,模拟蜂蜜或非常小规模的流体流动。但在我们希望制作动画的大多数其他情况下,粘度只起次要作用,因此我们经常放弃它:方程越简单越好。事实上,大多数模拟流体的数值方法不可避免地会引入误差,这些误差可以在物理上重新解释为粘度(稍后会详细介绍)——所以即使我们在方程中降低粘度,我们仍然会得到看起来像它的东西。事实上,计算流体动力学的一大挑战是尽可能避免这种虚假的粘性误差。因此,对于本书的其余部分,除了第10章侧重于高粘度或者粘度变化的流体外,我们将假设粘度已经被放弃。
  没有粘性的Navier-Stokes方程称为Euler方程,而这种没有粘性的理想流体称为无粘性流体。为了清楚地说明已删除的内容,这里是不可压缩的欧拉方程,使用物质导数来强调简单性:

DuDt+1ρp=g,·u=0.

  我们将主要使用这些方程。不要忘记它们是进一步的近似,并不是说水和空气实际上是理想的无粘性流体——只是它们的粘度对数值模拟的贡献通常与模拟中的其他误差相形见绌,因此不值得建模.

1.6 Boundary Conditions

  数值模拟流体的大部分“乐趣”在于使边界条件正确。到目前为止,我们只讨论了流体内部发生了什么:那么在边界处发生了什么呢?
  现在,让我们只关注两个边界条件,固体壁自由表面
  固体壁边界是流体与固体接触的地方。用速度来表达这一点最简单:流体最好不要流入或流出固体,因此速度的法向分量必须为零:

u·n^=0.

  如果固体本身也在移动,事情会稍微复杂一些。一般来说,我们需要流体速度的法向分量来匹配固体速度的法向分量,这样相对速度的法向分量为零:

u·n^=usolid·n^.

  在这两个方程中,n^当然是实体边界的法线。这有时被称为不粘条件,因为我们只限制速度的法向分量,允许流体沿切线方向自由滑过。这是要记住的重要一点:流体的切向速度可能与固体的切向速度完全没有关系。
  这就是速度的作用:固体壁上的压力又是什么情况呢?我们再次回到压力是“使流体不可压缩所需要的一切”的想法。我们将对此进行补充,“并强制实施实体壁边界条件。” 动量方程中的p/ρ项甚至适用于边界,因此对于在实心壁上控制u·n^的压力,显然这是关于 p·n^,也称为压力的正常导数:∂ p/∂^n。我们将等到对边界条件进行数值处理后再进行更精确的处理。
  这就是非粘性流体的固体壁边界的全部内容。如果我们有粘性,情况就会变得更加复杂。在这种情况下,固体的粘性通常会影响流体速度的切向分量,迫使其匹配。这称为无滑移边界条件,简单地说就是:

u=0,

  或者如果固体在移动,

u=usolid.

  对于具有非常低但非零粘度的流体,在观察固体旁边的流动的微观细节时,这实际上比不粘条件更准确,但通常的粘度的影响和粘性条件只出现在固体周围一个非常薄的边界层,并且在流体的其他地方,我们就好像有一个不粘边界。边界层是一个困难的研究领域,但由于它们通常比我们在动画模拟中可以解决的要薄得多,所以我们不用担心它们并且只要在除了高粘度的情况下使用不粘条件。同样,在我们进入数值实现之前,我们将避免讨论确切的细节。
  如果您在从固体上落下一滴液体的情况下考虑这些边界条件,您可能会感到困惑:如果固体表面的速度法向分量为零,它如何真正与固体分离?在某种程度上,答案是连续统模型并没有真正正确地处理分离:这是一个分子层面的过程,不能被抽象化。通常,我们要么从模拟方法中的数值误差中获得合理的分离,要么如果看起来不合理,我们就将其破解。有人在图形中引入了一种更有原则的“宏观”分离模型,该模型具有不等式条件,允许流体离开墙壁但不进入墙壁:

u·n^usolid·n^.

  这种情况的后果,无论是物理的还是算法的,仍然很不清楚,所以我们不会在本书中进一步跟进,尽管它是一个主要的研究方向。
  顺便说一句,有时固体壁实际上是流体可以通过的通风口或排水口:在这种情况下,我们显然希望u·n^不同于壁速度。它应该是流体在该点被泵入或泵出模拟的速度。
  我们感兴趣的另一个边界条件是自由表面。这是我们停止对流体建模的地方。例如,如果我们模拟水溅到周围,那么没有与实心墙接触的水面是自由表面。在这种情况下,确实存在另一种流体,空气,但我们可能也不想要模拟空气的麻烦。而且由于空气比水轻700倍,因此无论如何它都不会对水产生那么大的影响(除了一些明显的例外,比如气泡!)。因此,我们简化了模型,将空气表示为具有恒定大气压的区域。实际上,由于只有压力差(在不可压缩流动中),我们可以将气压设置为任意常数:最方便的是零。 因此,自由表面是p=0的表面,我们不以任何特定方式控制速度。
  出现自由表面的另一种情况是我们试图模拟属于更大域一部分的流体:例如,模拟露天烟雾。我们显然无法模拟地球的整个大气层,所以我们只需要制作一个网格来覆盖我们期望“有趣”的区域。(我在这里先声明一下,要模拟烟雾,您还需要模拟附近的无烟空气,而不仅仅是烟雾区域本身——但是,我们可以避免模拟远离所有动作的空气。) 流体继续流过模拟区域的边界,但我们没有跟踪它;我们允许流体随意进出该区域,因此很自然地认为这是一个自由表面,p=0,即使实际上没有可见表面。
  关于自由表面的最后一点说明:对于较小规模的液体,表面张力可能非常重要。在潜在的分子水平上,表面张力的存在是因为不同类型的分子之间的吸引强度不同。例如,水分子对其他水分子的吸引力比对空气分子的吸引力更大:因此,分离水和空气的表面的水分子试图移动以尽可能地被水包围。 从几何的角度来看,物理化学家将其建模为试图最小化表面积的力,或者等效地,试图降低表面的平均曲率。您可以将第一个想法(最小化表面积)解释为不断试图缩小表面的张力,因此名称为表面张力;使用平均曲率使用第二种方法可能更方便一些。(稍后,在第 8 章中,我们将讨论如何实际测量平均曲率以及它的确切含义。)简而言之,模型是两种流体之间实际上存在压力跳跃,与平均曲率成正比 :

[p]=γκ.

  符号[p]表示压力的跳跃,即在水侧测量的压力与在空气侧测量的压力差,γ是您可以查找的表面张力系数(对于室温下的水和空气它约为γ0.073N/m),κ 是平均曲率,以 1/m 为单位测量。对于具有表面张力的自由表面,这意味着水表面的压力是气压(我们假设为零)加上压力跳跃:

p=γκ.(1.4)

  自由表面确实存在一个主要问题:气泡会立即坍塌(内部没有压力来阻止它们失去体积)。虽然空气比水轻得多,因此通常可能无法将大量动量传递给水,但它仍然具有不可压缩性约束:水中的气泡在很大程度上保持其体积。用自由表面对气泡建模将使气泡崩溃并消失。 要处理这种情况,您需要基于将气泡粒子添加到自由表面流中的技巧,或者更一般地模拟空气和水(称为两相流,因为涉及两种相或两种类型的流体)。

posted @   律四  阅读(311)  评论(0编辑  收藏  举报
相关博文:
阅读排行:
· 清华大学推出第四讲使用 DeepSeek + DeepResearch 让科研像聊天一样简单!
· 推荐几款开源且免费的 .NET MAUI 组件库
· 实操Deepseek接入个人知识库
· 易语言 —— 开山篇
· Trae初体验
点击右上角即可分享
微信分享提示