【地质】连续性方程在地球动力学流体模拟中的意义

1. 介绍
地球的结构和其他类地行星类似，是层状的，他们可以通过物理，化学及流体学的特性来进行区分。从大的尺度来讲，地球拥有一个富含硅质的地壳，一层较软，熔融状的地幔，一层液体的外核和一个固体的内核。一直以来的科学研究表明，地壳和地幔，可以被看作是连续的地质介质，或者被称为不包含孔洞及间隙的地质材料。科学家们通常采用场变量（field variables）来描述连续材料的物理特性，比如
1. 标量（Scalars）：压力，温度，质量等
2. 矢量（Vectors）：速度，热通量等
3. 张量（Tensors）：应力，应变速率等

图片1：场变量可以用两种方式来描述，如图一所示 (a)完全连续的场变量（b）离散-连续的场变量（线性插值的连续性表现的精度取决于节点的数量）（Gerya，2009）

所以如果想要计算场变量里面的数值，那么就需要用到线性或者非线性的插值算法来计算节点之间的差值。地质介质的连续性通常意味着物质的移动的错位并不产生孔洞和间隙，这一连续的物质行为很近似于水的的形变，或者说是流体的形态变化。我们可以通过这一特性，推断出地壳和地幔的形变是等同于长期的缓速蔓延的流体动力学来进行演化的。

图片2：模拟刚性的方块在粘性介质中沉降的形变，对应地壳俯冲进入地幔的变化（例如，日本下方的俯冲带）。关键参数为密度，密度差距过大时几乎没有形变发生，密度差距较小时（10倍以内）会出现较大的形变。（模型代码：I2VIS）

关于如何将对于介质连续性的定性理解转化为定量数学形式，这一转化以连续性方程的形式广泛的用于地球动力学的模型建立。连续性方程可以被用来描述连续介质位移期间的质量守恒。作为自然界普遍存在的基本定律之一，质量守恒定律指出对于任何物质和能量全部转移的系统来说，系统的质量必须随着时间的推移保持不变，因为系统质量不能改变，不能增加或消除。因此，质量随着时间的推移而保持不变。质量守恒定律的连续性方程（以及许多其他依赖于时间的守恒方程）可以通过欧拉或者格拉朗日对于流体运动不同的观察方式推导而出。

所以在这里必须要先了解流体力学中到欧拉瑞恩和拉格朗日对于流体运动的不同描述方法。欧拉为瑞士数学家，物理学家。1707年4月15日生于瑞士巴塞尔，1735年因劳累导致右目失明， 1766年双目全盲后仍继续从事科学研究，1783年9月13日卒于俄国彼得堡。在1765年出版的《刚体运动理论》一书中。欧拉提出了质点动力学微分方程可以应用于液体。他使用两种方法来描述流体的运动，即分别根据空间固定点(1755年)和根据确定流体质点(1759年)描述流体速度场。这两种方法通常分别称为“欧拉描述方法”和“拉格朗日描述法”。欧拉奠定了理想流体(假设流体不可压缩，且其粘性可忽略)的运动理论基础，给出反映质量守恒的连续性方程(1752年)和反映动量变化规律的流体动力学方程(1755年)。

拉格朗日是欧拉的继任者，1736年1月25日生于意大利都灵，1813年4月10日卒于巴黎，他最先提出速度势和流函数的概念来奠定了流体无旋运动理论的基础。他在1788年出版的《分析力学》一书中提出了了“拉格朗日描述法”，即着眼于流体质点，描述每个流体质点自始至终的运动过程。

2. 欧拉瑞恩描述法
首先，欧拉瑞恩法侧重于“场”，把流体的性质（质量密度、速度等）定义为空间位置+时间的函数（王力乐，2014）。观测者位于空间中的某个不动点（固定点）来对流体的状态进行观测，用来记录每一时刻流过这个点的流体的速度，比如说t1时刻质点1流过这个空间点，我们就记录他的速度v1，t2时刻质点2流过这个点，我们记录速度v2（耿添，2014）。在关于物质守恒的问题中，我们主要关注的物理参数包括质量，速度，还有时间。所以欧拉连续方程式可以表示如下：

(1) $\begin{equation*} \frac{\partial \rho}{\partial t}+\operatorname{div}(\rho \vec{v})=0. \end{equation*}$

式中， $\frac{\partial}{\partial t}$ 为欧拉时间导数； $\rho$ 为局域密度，用来表示每个单位空间内有多少的物质质量 $\left(\mathrm{kg} /\mathrm{m}^{3}\right)$ ; $\vec{v}$ 表示局部速度；最后 $\operatorname{div}()$ 则是散度算子。可以用以下公式来定义，

(2) $\begin{equation*} \operatorname{div}(\vec{v})=\frac{\partial v_{x}}{\partial x}+\frac{\partial v_{y}}{\partial y}+\frac{\partial v_{z}}{\partial z} \end{equation*}$

本文在这里只列出三维的表达式，即包括 $x$ , $y$ , $z$ 两个笛卡尔坐标（Cartesian coordinates）和 $v_{x}$ , $v_{y}$ 两个平行于坐标轴的速度矢量 $\vec{v}$ 。简单的来讲，当散度为正数的时候，矢量场是发散的，即流体离开观察点向周围四散而去。当散度为负数的时候，矢量场是收敛的，即流体从四面八方涌向观察点。最后当散度等于0的时候，流体就这么自然的流动着，既没有发散也没有收敛（图片3）。

图片3：二维速度场环绕质点的图例（a）发散，（b）收敛，和（c）自然。箭头指示方向。

从质量守恒的角度来看，我们想象空间是一个小屋子，里面有若干物体，每个物体都有质量，当 $\operatorname{div}(\vec{v})>0$ （发散）的时候，就意味着这些物体不断地离开小屋子，屋子内的空间越来越空旷，这说明屋子里的物质密度随着时间下降了（例如： $\frac{\partial \rho}{\partial t}<0$ ）。所以 $\rho \vec{v}$ 在这里没有别的意思，它仅仅表明了局部物质密度的变化（比如水流的多少）。在空间中可以被分解为三个坐标上的量，公式如下：

(3) $\begin{equation*} \rho \vec{v}=\left(\rho v_{x}, \rho v_{y}, \rho v_{z}\right) \end{equation*}$

空间内密度的变化公式如下

(4) $\begin{equation*} \rho_{a}=\frac{m_{a}}{\Delta x \Delta y \Delta z} \end{equation*}$

通过分析流体在单位时间内通过固定的小房子的量（在场内的固定点观察流体密度的变化）就可以推导出欧拉连续方程，即质量 $m$ 是随着时间变化，空间维度 $\Delta x, \Delta y, \Delta z.$ 保持不变。最终当 $\Delta t, \Delta x, \Delta y, \Delta z.$ 全部趋近于0的时候，欧拉连续方程可以写作：

(5) $\begin{equation*} \frac{\partial \rho}{\partial t}+\frac{\partial\left(\rho v_{x}\right)}{\partial x}+\frac{\partial\left(\rho v_{y}\right)} {\partial y}+\frac{\partial\left(\rho v_{z}\right)}{\partial z}=0 \end{equation*}$

3. 拉格朗日描述法
拉格朗日表述侧重于“质点”，其关注的是质点的初始位置和当前时刻。用这个时刻坐标来标记质点，记录这个质点每时每刻所在的位置，简单的说就是站在质点上观察流体。图片4中对于这两种观察方法进行了模拟和视觉化处理。

图片4：绿色为质点，箭头表示速度大小，黄色线条为粒子路径，白色箭头为矢量场。（1）欧拉观察：通过质子在空间内的移动观察到了黄色的轨迹。（2）拉格朗日观察：质子移动的速度相对于周围的场一直在发生变化，如绿色箭头所示。使用YAFFA进行模拟。使用Paraview进行后期处理（Tommi Mikkola，2015）。

所以拉格朗日连续方程式可以表示如下：

(6) $\begin{equation*} \frac{D \rho}{D t}+\rho \operatorname{div}(\vec{v})=0 \end{equation*}$

式中 $\frac{D}{D t}$ 为拉格朗日时间导数，其他的参数在前文中以及提到过，这里就不再详细描述（详见公式1）。

格拉朗日描述法有着不一样的视角，如以下公式所示：

(7) $\begin{equation*} \rho_{a}=\frac{m}{\Delta x_{a} \Delta y_{a} \Delta z_{a}} \end{equation*}$

此处 $a$ 为变量，因为我们在质点上观察流体，所以质量 $m$ 是恒定的。而此处空间维度 $\Delta x, \Delta y, \Delta z.$ 会随着时间以及流体的运动而发生变化。最终拉格朗日连续方程如下所示。

(8) $\begin{equation*} \frac{D \rho}{D t}+\rho \frac{\partial v_{x}}{\partial x}+\rho \frac{\partial v_{y}}{\partial y}+\rho \frac{\partial v_{z}}{\partial z}=0 \end{equation*}$

4. 欧拉描述和格拉朗日描述的对比-对流方程
因为这两种观察法实际上是观测的同一种流体运动，只是站在不同的角度之上，所以理论上来讲这两种坐标是可以互相转化的。这里欧拉的坐标基于时间，拉格朗日的坐标注重位置的变化，基于位移等于速度乘以时间，这里提出了对流方程。
这里对流方程如下所示：

(9) $\begin{equation*} \frac{D \rho}{D t}=\frac{\partial \rho}{\partial t}+\vec{v} \operatorname{grad}(\rho) \end{equation*}$

公式中， $\vec{v} \operatorname{grad}(\rho)$ 对应着欧拉连续方程内固定点的物质密度改变。拉格朗日连续方程并不依赖于物质密度改变。这里举一个情景，流体以一个稳定的速度 $v_{x}$ 在小房间里通过。在这个时间点 $t_{0}$ 欧拉站在屋子里的一个固定点观察在流体，拉格朗日则是选择在流体中观察。随着流体的移动，每个时间点都对应着一个新的观察结果（密度 $/rho$ )
欧拉的观察结果：

(10) $\begin{equation*} \frac{\partial \rho}{\partial t} \approx \frac{\Delta \rho_{t}}{\Delta t}=\frac{\rho_{\mathrm{B}}-\rho_{\mathrm{A}}}{t_{1}-t_{0}} \end{equation*}$

拉格朗日的观察结果：

(11) $\begin{equation*} \frac{\partial \rho}{\partial x} \approx \frac{\Delta \rho_{x}}{\Delta x}=\frac{\rho_{\mathrm{A}}-\rho_{\mathrm{B}}}{x_{1}-x_{0}} \end{equation*}$

式中， $\Delta \rho_{t}$ 对应的是密度随时间的改变（欧拉）， $\Delta \rho_{x}$ 则是对应密度在不同位置的改变（拉格朗日）。当时间的区间 $\Delta t=t_{1}-t_{0}$ 很小时，我们可以得出流体在点A和B之间的位移 $\Delta x=x_{1}-x_{0} )$ 同样很小。此时

$\begin{array}{c}{\Delta x=v_{x} \Delta t} \\ {\Delta \rho_{x}=-\Delta \rho_{t}}\end{array}$

我们可以得到一维的表达式：

(12) $\begin{equation*} \frac{\Delta \rho_{t}}{\Delta t}=-v_{x} \frac{\Delta \rho_{x}}{\Delta x} \end{equation*}$

对流公式可以用来描述物质在非均匀连续介质中的移动的特点，比如密度，温度，成分等的改变等。

5. 特殊情况
对于地球的地壳和地幔，还有一种可能的情况是材料的密度不会随着时间来发生改变，即介质不可以被压缩
这个时候，如果压力和温度的变化不大，并且没有相位转变来改变空间容积，那么表达式可以写作：

(13) $\begin{equation*} \frac{D \rho}{D t}=\frac{\partial \rho}{\partial t}+\vec{v} \operatorname{grad}(\rho)=0 \end{equation*}$

此时欧拉和拉格朗日不可压缩连续方程是一致的，即：
/[\operatorname{div}(\vec{v})=0/] 这种特殊情况对应了一些地球动力学实例，比如俯冲板块之上的位于地幔楔里的拐角流（Turcotte and Schubert, 2002）和俯冲隧道中的物质混杂环流(Cloos, 1982).

小结：连续性方程是理解流体运动中的基础，，其中两种形式，分别对应了欧拉和拉格朗日两种对于流体的观察方式，这两者之间是可以互相转化的，本质上是还是在利用质量守恒定律。这两种描述最大的区别则是在于对流方程，借助对流方程，地质学家可以对地球动力学流体模拟中物质分布变化的问题展开研究。最后，本文仅代表作者的观点和浅显的理解，依旧有很多不足以及误解的部分，欢迎有兴趣的朋友来一起参与讨论和一起学习进步

6.参考文献：
Gerya, T. (2009). Introduction to Numerical Geodynamic Modelling (Vol. 9780521887540, pp. 1–345). Cambridge: Cambridge University Press. https://doi.org/10.1017/CBO9780511809101

Turcotte, D. L., Schubert, G. (2002) Geodynamics. Cambridge University Press.

Cloos, M. (1982) Flow melanges – numerical modeling and geologic constraints on their
origin in the Franciscan subduction complex, California. Geological Society of
America Bulletin, 93, 330–45.

耿添.（2014）怎样理解流体力学中的拉格朗日描述和欧拉描述. 知乎. 链接：https://www.zhihu.com/question/26129680/answer/32231957

王力乐, （2014）怎样理解流体力学中的拉格朗日描述和欧拉描述. 知乎. 链接：https://www.zhihu.com/question/26129680/answer/32457812

Bowen2019年5月20日

1,488 1 minute read

【地质】连续性方程在地球动力学流体模拟中的意义

Bowen

【美食】蜂蜜柚子茶