1 前言
作为力学出生的人，对自己力学专业的一些方程有特别的情感。因此这篇文章简单记录一下流体力学的几个方程，以及其在线性分布假设情况下的一种离散形式。
不管是无人驾驶当中的多传感器融合问题，还是力学模型仿真问题，它本身可以看成一个优化问题。既然是优化问题，就需要确定出优化变量以及变量之间的约束方程。整个模型可以看成一个大的优化模型，而变量之间通过约束相互耦合，进而所有变量通过约束相互影响。力学中控制方程的离散形式，决定了约束方程建立的方式。
在离散公式得出前，先约定如下变量含义：
\(\rho\): 密度，单位 \(kg/m^3\).
\(p\): 压力，单位\(Pa\).
\(e\): 单位质量内能，单位\(m^2/s^2\).
\(c\): 摩尔浓度，单位 \(mol/m^3\).
\(v_x\): \(x\)方向流速，单位\(m/s\).
\(v_y\): \(y\)方向流速，单位\(m/s\).
\(v_z\): \(z\)方向流速，单位\(m/s\).
\(v\)：速度，其值满足\(v^2=v_x^2+v_y^2+v_z^2\).
\(g\): 重力加速度，单位 \(m/s^2\).
\(f_x\): \(x\)方向体积力加速度，单位\(m/s^2\).
\(f_y\): \(y\)方向体积力加速度，单位\(m/s^2\).
\(f_z\): \(z\)方向体积力加速度，单位\(m/s^2\).
\(T\)：温度，单位°C.
\(k\)：导热系数(thermal conductivity).
\(D\)：溶质在液体里的扩散系数，单位为\(m^2/s\).
\(\mu\)：运动粘度系数(molecular viscosity coefficient)，单位为\(Pa \cdot s\).
\(\gamma\)：体积粘度系数(bulk viscosity coefficient), 其值满足\(\gamma = -\dfrac{2}{3}\mu\).
2 线性分布假设
对于平面，可以将一个多边形区域，使用三角形单元进行网格划分。不失一般性，设三角形单元三个节点坐标为\((x_0,y_0)\)、\((x_1,y_1)\)、\((x_2,y_2)\)，设某物理量(这里以压力为例)在三个节点上的值分别为\(p_0、p_1、p_2\)。现在假设物理量在三角形单元内为线性分布，即对三角形单元内任意一点\((x,y)\)，其对应的物理量\(p\)满足如下形式
(2-1)
$$
p=a_p x + b_b y + c_p
$$
公式里\(a_p、b_p、c_p\)为变量\(p\)在单元里的线性分布系数。这里设置如下矩阵\(M\)
(2-2)
$$
M = \begin{bmatrix}
x_0 & y_0 & 1
\\
x_1 & y_1 & 1
\\
x_2 & y_2 & 1
\end{bmatrix}
$$
则系数\((a_p,b_p,c_p)\)可以用如下形式表示
(2-3)
$$
\begin{bmatrix}
a_p
\\
b_p
\\
c_p
\end{bmatrix} = M^{-1}
\begin{bmatrix}
p_0
\\
p_1
\\
p_2
\end{bmatrix}
$$
这里，设三角形单元形心坐标为
(2-4)
$$
\begin{cases}
x_c=\dfrac{x_0+x_1+x_2}{3}
\\
y_c=\dfrac{y_0+y_1+y_2}{3}
\end{cases}
$$
然后，设置如下中间变量
(2-5)
$$
\mathbf{m}^T =\begin{bmatrix}m_0 & m_1 & m_2 \end{bmatrix}= \begin{bmatrix}
x_c & y_c & 1
\end{bmatrix}M^{-1}
$$
则形心上的物理量\(p_c\)可表示为
(2-6)
$$
\begin{aligned}
p_c &= a_p x_c+b_py_c+c_p
\\
&=\begin{bmatrix}
x_c & y_c & 1
\end{bmatrix}M^{-1}
\begin{bmatrix}
p_0
\\
p_1
\\
p_2
\end{bmatrix}
\\
&=m_0p_0+m_1p_1+m_2p_2
\end{aligned}
$$
很明显，(2-6)显示了形心上的物理量与三个节点上物理量的关系。
进一步，结合(2-1)到(2-6)，很明显如下偏导成立
(2-7)
$$
\begin{cases}
\dfrac{\partial p_c}{\partial x} = a_p
\\
\\
\dfrac{\partial p_c}{\partial y} = b_p
\\
\\
\dfrac{\partial p_c}{\partial p_0} = m_0
\\
\\
\dfrac{\partial p_c}{\partial p_1} = m_1
\\
\\
\dfrac{\partial p_c}{\partial p_2} = m_2
\\
\\
\dfrac{\partial p_c}{\partial t} = m_0\dfrac{\partial p_0}{\partial t} + m_1\dfrac{\partial p_1}{\partial t} + m_2\dfrac{\partial p_2}{\partial t}
\end{cases}
$$
根据线性分布假设，结合(2-7)可以得出\(m_0+m_1+m_2=1\)。
另外，在计算中，可能涉及到系数\(a_p、b_p、c_p\)对变量\(p_0、p_1、p_2\)求导的情况，因此根据(2-3)可以得到
(2-8)
$$
\begin{bmatrix}
\dfrac{\partial a_p}{\partial p_0} & \dfrac{\partial a_p}{\partial p_1} & \dfrac{\partial a_p}{\partial p_2}
\\
\\
\dfrac{\partial b_p}{\partial p_0} & \dfrac{\partial b_p}{\partial p_1} & \dfrac{\partial b_p}{\partial p_2}
\\
\\
\dfrac{\partial c_p}{\partial p_0} & \dfrac{\partial c_p}{\partial p_1} & \dfrac{\partial c_p}{\partial p_2}
\end{bmatrix}=M^{-1}
$$
同理，对于密度、速度物理量在三角形单元内均做线性分布假设，则有上面类似的结论。这里假设密度的线性分布系数为\((a_\rho,b_\rho,c_\rho)\)，\(x\)方向流速线性分布系数为\((a_x,b_x,c_x)\)，\(y\)方向流速线性分布系数为\((a_y,b_y,c_y)\).
3 控制方程
这里的控制方程，是指物质点上物理量需要满足的约束方程。为了简化计算，一般会在空间中选取多个控制点（比如常见的单元上的节点），让控制点满足某些约束，而控制点之间通过某种方式进行关联起来，进而可以模拟一个较大区域的物理变化情况。流体力学中，一般涉及到质量守恒、动量守恒、能量守恒控制方程，有时也会使用扩散方程，或者其它形式的控制方程。
3.1 质量守恒
质量守恒方程也叫连续性方程，这个控制方程描述了任意微元流体在流动过程中质量保持不变的性质。这里使用文献[1]里的守恒型质量控制方程来描述
(3.1-1)
$$
\dfrac{\partial \rho}{\partial t} + \dfrac{\partial(\rho v_x)}{\partial x} + \dfrac{\partial(\rho v_y)}{\partial y} + \dfrac{\partial(\rho v_z)}{\partial z} = 0
$$
对于二维情况，可以得到如下变量表达式
(3.1-2)
$$
\begin{aligned}
0&=\dfrac{\partial \rho}{\partial t}+\dfrac{\partial(\rho v_x)}{\partial x}+\dfrac{\partial(\rho v_y)}{\partial y}
\\
&=\dfrac{\partial \rho}{\partial t}+v_x\dfrac{\partial \rho}{\partial x}+v_y\dfrac{\partial \rho}{\partial y}+\rho(\dfrac{\partial v_x}{\partial x}+\dfrac{\partial v_y}{\partial y})
\end{aligned}
$$
因此，对二维模型使用线性分布假设，则可以得到如下的质量守恒方程
(3.1-3)
$$
\dfrac{\partial \rho}{\partial t} + v_x a_\rho + v_y b_\rho + \rho(a_x+b_y)=0
$$
3.2 动量守恒
动量守恒方程也叫运动方程，流体的动量守恒方程描述了流体总动量保持不变的特性。这里使用文献[1]里的守恒型动量控制方程来描述
(3.2-1)
$$
\begin{cases}
\dfrac{\partial (\rho v_x)}{\partial t} + \dfrac{\partial(\rho v_x v_x)}{\partial x} + \dfrac{\partial(\rho v_x v_y)}{\partial y} + \dfrac{\partial(\rho v_x v_z)}{\partial z} = -\dfrac{\partial p}{\partial x} + \dfrac{\partial \tau_{xx}}{\partial x} + \dfrac{\partial \tau_{yx}}{\partial y} + \dfrac{\partial \tau_{zx}}{\partial z} + \rho f_x
\\
\\
\dfrac{\partial (\rho v_y)}{\partial t} + \dfrac{\partial(\rho v_y v_x)}{\partial x} + \dfrac{\partial(\rho v_y v_y)}{\partial y} + \dfrac{\partial(\rho v_y v_z)}{\partial z} = -\dfrac{\partial p}{\partial y} + \dfrac{\partial \tau_{xy}}{\partial x} + \dfrac{\partial \tau_{yy}}{\partial y} + \dfrac{\partial \tau_{zy}}{\partial z} + \rho f_y
\\
\\
\dfrac{\partial (\rho v_z)}{\partial t} + \dfrac{\partial(\rho v_z v_x)}{\partial x} + \dfrac{\partial(\rho v_z v_y)}{\partial y} + \dfrac{\partial(\rho v_z v_z)}{\partial z} = -\dfrac{\partial p}{\partial z} + \dfrac{\partial \tau_{xz}}{\partial x} + \dfrac{\partial \tau_{yz}}{\partial y} + \dfrac{\partial \tau_{zz}}{\partial z} + \rho f_z
\end{cases}
$$
上面公式中\(\tau_{xx}、\tau_{yy}、\tau_{zz}、\tau_{xy}、\tau_{yx}、\tau_{xz}、\tau_{zx}、\tau_{yz}、\tau_{zy}\)取如下值
(3.2-2)
$$
\begin{cases}
\tau_{xx} = 2\mu \dfrac{\partial v_x}{\partial x} + \gamma(\dfrac{\partial v_x}{\partial x} + \dfrac{\partial v_y}{\partial y} + \dfrac{\partial v_z}{\partial z})
\\
\\
\tau_{yy} = 2\mu \dfrac{\partial v_y}{\partial y} + \gamma(\dfrac{\partial v_x}{\partial x} + \dfrac{\partial v_y}{\partial y} + \dfrac{\partial v_z}{\partial z})
\\
\\
\tau_{zz} = 2\mu \dfrac{\partial v_z}{\partial z} + \gamma(\dfrac{\partial v_x}{\partial x} + \dfrac{\partial v_y}{\partial y} + \dfrac{\partial v_z}{\partial z})
\\
\\
\tau_{xy} = \tau_{yx}=\mu(\dfrac{\partial v_x}{\partial y} + \dfrac{\partial v_y}{\partial x})
\\
\\
\tau_{yz} = \tau_{zy}=\mu(\dfrac{\partial v_y}{\partial z} + \dfrac{\partial v_z}{\partial y})
\\
\\
\tau_{zx} = \tau_{xz}=\mu(\dfrac{\partial v_z}{\partial x} + \dfrac{\partial v_x}{\partial z})
\end{cases}
$$
对于2维情况，且y取竖直方向，则可简化为如下公式
(3.2-3)
$$
\begin{cases}
v_x\dfrac{\partial \rho}{\partial t}+\rho\dfrac{\partial v_x}{\partial t}
+v_x^2\dfrac{\partial \rho}{\partial x}+2\rho v_x\dfrac{\partial v_x}{\partial x}+v_xv_y\dfrac{\partial\rho}{\partial y}+\rho(v_y\dfrac{\partial v_x}{\partial y}+ v_x\dfrac{\partial v_y}{\partial y}) = -\dfrac{\partial p}{\partial x} + \dfrac{\partial \tau_{xx}}{\partial x} + \dfrac{\partial \tau_{yx}}{\partial y}
\\
\\
v_y\dfrac{\partial \rho}{\partial t}+\rho\dfrac{\partial v_y}{\partial t}+
v_y^2\dfrac{\partial \rho}{\partial y}+2\rho v_y\dfrac{\partial v_y}{\partial y}+v_xv_y\dfrac{\partial \rho}{\partial x}+\rho(v_y\dfrac{\partial v_x}{\partial x}+v_x\dfrac{\partial v_y}{\partial x})= -\dfrac{\partial p}{\partial y} + \dfrac{\partial \tau_{xy}}{\partial x} + \dfrac{\partial \tau_{yy}}{\partial y} + \rho g
\end{cases}
$$
这里使用线性假设分布，因为速度满足线性分布假设，则可以得到
(3.2-4)
$$
\dfrac{\partial \tau_{xx}}{\partial x}=0,
\dfrac{\partial \tau_{yy}}{\partial y}=0,
\dfrac{\partial \tau_{xy}}{\partial x}=0,
\dfrac{\partial \tau_{yx}}{\partial y}=0
$$
因此，对二维模型使用线性分布假设，结合(3.2-3)、(3.2-4)可以得到如下动量守恒方程
(3.2-5)
$$
\begin{cases}
v_x\dfrac{\partial \rho}{\partial t}+\rho\dfrac{\partial v_x}{\partial t}+v_x^2 a_\rho+2\rho v_x a_x + v_xv_yb_\rho+\rho(v_yb_x+v_xb_y)+a_p = 0
\\
\\
v_y\dfrac{\partial \rho}{\partial t}+\rho\dfrac{\partial v_y}{\partial t}+v_y^2b_\rho+2\rho v_yb_y+v_xv_ya_\rho+\rho(v_ya_x+v_xa_y)+b_p-\rho g = 0
\end{cases}
$$
3.3 能量守恒
能量守恒也叫能量方程，是指在考虑流速、密度、温度、内能变化时, 微元流体总能量不变的特性。对于能量方程，这里采用[1]里的守恒型能量控制方程来描述
(3.3-1)
$$
\begin{aligned}
&\dfrac{\partial}{\partial t}(\rho (e + v^2/2)) + \dfrac{\partial}{\partial x}(\rho (e + v_x v^2/2)) + \dfrac{\partial}{\partial y}(\rho (e + v_y v^2/2))+ \dfrac{\partial}{\partial z}(\rho (e + v_z v^2/2))\\=&\rho q-\dfrac{\partial (pv_x)}{\partial x} -\dfrac{\partial (pv_y)}{\partial y} -\dfrac{\partial (pv_z)}{\partial z}
\\
&+ \rho(f_xv_x + f_yv_y+f_zv_z)
\\
&+\dfrac{\partial}{\partial x}(k\dfrac{\partial T}{\partial x}) + \dfrac{\partial}{\partial y}(k\dfrac{\partial T}{\partial y})+ \dfrac{\partial}{\partial z}(k\dfrac{\partial T}{\partial z})
\\
&+\dfrac{\partial (v_x \tau_{xx})}{\partial x} +\dfrac{\partial (v_x \tau_{yx})}{\partial y} +\dfrac{\partial (v_x \tau_{zx})}{\partial z}
\\
&+\dfrac{\partial (v_y \tau_{xy})}{\partial x} +\dfrac{\partial (v_y \tau_{yy})}{\partial y} +\dfrac{\partial (v_y \tau_{zy})}{\partial z}
\\
&+\dfrac{\partial (v_z \tau_{xz})}{\partial x} +\dfrac{\partial (v_z \tau_{yz})}{\partial y} +\dfrac{\partial (v_z \tau_{zz})}{\partial z}
\end{aligned}
$$
上面公式里\(q\)为单位质量体积热能变化率，其在 \(x,y,z\)三个方向的分量分别为
(3.3-2)
$$
\begin{cases}
q_x = -k\dfrac{\partial T}{\partial x}
\\
\\
q_y = -k\dfrac{\partial T}{\partial y}
\\
\\
q_z = -k\dfrac{\partial T}{\partial z}
\end{cases}
$$
对于常温的二维模型，(3.3-1)可变为如下形式
(3.3-3)
$$
\begin{aligned}
&\dfrac{\partial}{\partial t}(\rho (e + v^2/2)) + \dfrac{\partial}{\partial x}(\rho (e + v_x v^2/2)) + \dfrac{\partial}{\partial y}(\rho (e + v_y v^2/2))
\\=&-\dfrac{\partial (pv_x)}{\partial x} -\dfrac{\partial (pv_y)}{\partial y}
+\dfrac{\partial (v_x \tau_{xx})}{\partial x} +\dfrac{\partial (v_x \tau_{yx})}{\partial y}
+\dfrac{\partial (v_y \tau_{xy})}{\partial x} +\dfrac{\partial (v_y \tau_{yy})}{\partial y} + \rho g v_y
\\
=&\dfrac{\partial (v_x \tau_{xx}+v_y\tau_{xy} – pv_x)}{\partial x}+\dfrac{\partial (v_y \tau_{yy}+v_x\tau_{yx} – pv_y)}{\partial y}+\rho g v_y
\end{aligned}
$$
因此，对于二维模型，使用线性分布假设，由(3.3-3)可以得到能量守恒方程
(3.3-4)
$$
\begin{aligned}
&\dfrac{\partial}{\partial t}(\rho (e + v^2/2)) + \rho (\dfrac{\partial e}{\partial x} + \dfrac{\partial e}{\partial y})\\=
&-(p(a_x+b_y)+v_xp_a+v_yp_b)
\\
&-a_\rho(e+v_xv^2/2) – b_\rho(e+v_yv^2/2) – \dfrac{\rho}{2}(3a_xv_x^2+a_xv_y^2+2a_yv_xv_y + 3b_y^2v_y^2+b_yv_x^2+2v_xv_yb_x)
\\
&+a_x\tau_{xx}+b_x\tau_{xy}+a_y\tau{xy}+b_y\tau_{yy}+\rho g v_y
\end{aligned}
$$
3.4 扩散方程
扩散方程表征了流动系统的质量传递规律，这里将文献[2]里的数学模型进行速度梯度项扩展来描述
(3.4-1)
$$
\dfrac{\partial c}{\partial t} + c(\dfrac{\partial v_x}{\partial x}+ \dfrac{\partial v_y}{\partial y} + \dfrac{\partial v_z}{\partial z})+(v_x\dfrac{\partial c}{\partial x}+ v_y\dfrac{\partial c}{\partial y} + v_z\dfrac{\partial c}{\partial z})=D(\dfrac{\partial^2 c}{\partial x^2}+\dfrac{\partial^2 c}{\partial y^2}+\dfrac{\partial^2 c}{\partial z^2})
$$
而对于二维情况，可以得到如下模型
(3.4-2)
$$
\dfrac{\partial c}{\partial t} + c(\dfrac{\partial v_x}{\partial x}+ \dfrac{\partial v_y}{\partial y})+(v_x\dfrac{\partial c}{\partial x}+ v_y\dfrac{\partial c}{\partial y})=D(\dfrac{\partial^2 c}{\partial x^2}+\dfrac{\partial^2 c}{\partial y^2})
$$
在方程里，由于涉及到摩尔浓度，而浓度一般和密度有很好的相关性，因此这里使用盐水浓度与密度关系[3]里的二阶模型将浓度与密度进行关联。即浓度与密度满足如下方程
(3.4-3)
$$
c = b_0 + b_1\rho + b_2\rho^2
$$
这里取如下中间变量
(3.4-4)
$$
\eta=\dfrac{\partial c}{\partial \rho} = b_1+2b_2\rho
$$
进而可以使用密度替换掉浓度变量，这里将(3.4-3)与(3.4-4)代入(3.4-2)有
(3.4-5)
$$
\eta \dfrac{\partial \rho}{\partial t}+c(\dfrac{\partial v_x}{\partial x}+\dfrac{\partial v_y}{\partial y})+\eta(v_x\dfrac{\partial \rho}{\partial x}+v_y\dfrac{\partial \rho}{\partial y})=D(\partial(\eta\dfrac{\rho}{\partial x})/\partial x+\partial(\eta\dfrac{\rho}{\partial y})/\partial y)
$$
因为假设密度满足线性分布，因此密度对位置的二阶导为0，因此上式可以变为
(3.4-6)
$$
\begin{aligned}
\eta \dfrac{\partial \rho}{\partial t}+c(\dfrac{\partial v_x}{\partial x}+\dfrac{\partial v_y}{\partial y})+\eta(v_x\dfrac{\partial \rho}{\partial x}+v_y\dfrac{\partial \rho}{\partial y})&=D(\partial(\eta\dfrac{\rho}{\partial x})/\partial x+\partial(\eta\dfrac{\rho}{\partial y})/\partial y)
\\
&=D(\dfrac{\partial \eta}{\partial x}\dfrac{\partial \rho}{\partial x}+\dfrac{\partial \eta}{\partial y}\dfrac{\partial \rho}{\partial y})
\\
&=2b_2D((\dfrac{\partial \rho}{\partial x})^2+(\dfrac{\partial \rho}{\partial y})^2)
\end{aligned}
$$
为了让方程误差量级尽量与(3.1-3)一致，这里根据(3.4-6)建立如下约束方程
(3.4-7)
$$
\dfrac{\partial \rho}{\partial t}+\dfrac{c}{\eta}(\dfrac{\partial v_x}{\partial x}+\dfrac{\partial v_y}{\partial y})+(v_x\dfrac{\partial \rho}{\partial x}+v_y\dfrac{\partial \rho}{\partial y})
-\dfrac{2b_2D}{\eta} ((\dfrac{\partial \rho}{\partial x})^2+(\dfrac{\partial \rho}{\partial y})^2)=0
$$
最终代入线性分布系数可以得到扩散方程
(3.4-8)
$$
\dfrac{\partial \rho}{\partial t}+\dfrac{c}{\eta}(a_x+b_y)+v_xa_\rho+v_yb_\rho-\dfrac{2b_2D}{\eta}(a_\rho^2+b_\rho^2)=0
$$
4 方程离散
上面的公式(3.1-3)、(3.2-5)、(3.3-4)、(3.4-8)为需要离散的控制方程。可以发现，这几个公式，因为线性假设的关系，已经没有对位置偏导项的存在。这时仅需考虑对时间项的偏导处理，而时间偏导可以用差分代替。
由于(3.1-3)反映了质量守恒规律，而(3.4-8)反映了质量传递规律，因此在实际使用中按情况可使用其中一个方程。而公式(3.3-4)为能量守恒方程，在考虑温度与能量情况下，需要对此方程进行离散，否则就不需要离散。
接下来，本文将提供一种显示离散方式，其中使用前(后)差分格式代替求导。在提供离散公式前，这里约定：\(t_1=t_0+\delta t\)，其中\(\delta t\)表示时间步长，变量带上划线表示\(t_1\)时刻的值，变量不带上划线表示\(t_0\)时刻的值。例如，在时刻\(t_0\)时，压力、密度、\(x\)方向速度、\(y\)方向速度变量为\(p、\rho、v_x、v_y\)；在\(t_1\)时，压力、密度、\(x\)方向速度、\(y\)方向速度变量为\(\bar{p}、\bar{\rho}、\bar{v}_x、\bar{v}_y\)。
4.1 密度求解
如果不考虑浓度梯度引起的质量传递，则在\(t_0\)时刻，物理变量需要满足控制方程(3.1-3)，即
(4.1-1)
$$
\begin{aligned}
0 &= \dfrac{\partial \rho}{\partial t} + v_xa_\rho+v_yb_\rho+\rho(a_x+b_y)
\\
&=\dfrac{\bar{\rho}-\rho}{\delta t} + v_xa_\rho+v_yb_\rho+\rho(a_x+b_y)
\end{aligned}
$$
由上面公式可以得到
(4.1-2)
$$
\bar{\rho} = \rho – \delta t(v_xa_\rho+v_yb_\rho+\rho(a_x+b_y))
$$
如果考虑浓度梯度引起的质量传递，则在\(t_0\)时刻，物理变量需要满足控制方程(3.4-8)，即
(4.1-3)
$$
\begin{aligned}
0 &= \dfrac{\partial \rho}{\partial t} + \dfrac{c}{\eta}(a_x+b_y)+v_xa_\rho+v_yb_\rho-\dfrac{2b_2D}{\eta}(a_\rho^2+b_\rho^2)
\\
&= \dfrac{\bar{\rho} – \rho}{\delta t} + \dfrac{c}{\eta}(a_x+b_y)+v_xa_\rho+v_yb_\rho-\dfrac{2b_2D}{\eta}(a_\rho^2+b_\rho^2)
\end{aligned}
$$
由上面公式可以得到
(4.1-4)
$$
\bar{\rho}=\rho + \delta t(\dfrac{c}{\eta}(a_x+b_y)+v_xa_\rho+v_yb_\rho-\dfrac{2b_2D}{\eta}(a_\rho^2+b_\rho^2))
$$
4.2 速度求解
在\(t_0\)时刻，物理变量需要满足控制方程(3.2-5)，即
(4.2-1)
$$
\begin{cases}
\begin{aligned}
0&=v_x\dfrac{\partial \rho}{\partial t}+\rho\dfrac{\partial v_x}{\partial t}+v_x^2 a_\rho+2\rho v_x a_x + v_xv_yb_\rho+\rho(v_yb_x+v_xb_y)+a_p
\\
&=v_x \dfrac{\bar{\rho} – \rho}{\delta t} + \rho \dfrac{\bar{v}_x – v_x}{\delta t}+v_x^2 a_\rho+2\rho v_x a_x + v_xv_yb_\rho+\rho(v_yb_x+v_xb_y)+a_p
\\
\\
0&=v_y\dfrac{\partial \rho}{\partial t}+\rho\dfrac{\partial v_y}{\partial t}+v_y^2b_\rho+2\rho v_yb_y+v_xv_ya_\rho+\rho(v_ya_x+v_xa_y)+b_p-\rho g
\\
&=v_y \dfrac{\bar{\rho} – \rho}{\delta t} + \rho \dfrac{\bar{v}_y-v_y}{\delta t}+v_y^2b_\rho+2\rho v_yb_y+v_xv_ya_\rho+\rho(v_ya_x+v_xa_y)+b_p-\rho g
\end{aligned}
\end{cases}
$$
由上面公式可以得到
(4.2-2)
$$
\begin{cases}
\begin{aligned}
\bar{v}_x &= v_x – \dfrac{1}{\rho}(v_x(\bar{\rho} – \rho) + \delta t (v_x^2 a_\rho+2\rho v_x a_x + v_xv_yb_\rho+\rho(v_yb_x+v_xb_y)+a_p))
\\
\bar{v}_y &= v_y – \dfrac{1}{\rho}(v_y(\bar{\rho} – \rho) + \delta t(v_y^2b_\rho+2\rho v_yb_y+v_xv_ya_\rho+\rho(v_ya_x+v_xa_y)+b_p-\rho g))
\end{aligned}
\end{cases}
$$
4.3 压力梯度求解
在\(t_1\)时刻，物理变量需要满足控制方程(3.2-5)，即
(4.3-1)
$$
\begin{cases}
\begin{aligned}
0&=\bar{v}_x\dfrac{\partial \bar{\rho}}{\partial t}+\bar{\rho}\dfrac{\partial \bar{v}_x}{\partial t}+\bar{v}_x^2 \bar{a}_\rho+2\bar{\rho} \bar{v}_x \bar{a}_x + \bar{v}_x\bar{v}_y\bar{b}_\rho+\bar{\rho}(\bar{v}_y\bar{b}_x+\bar{v}_x\bar{b}_y)+\bar{a}_p
\\
&=\bar{v}_x\dfrac{\bar{\rho} – \rho}{\delta t}+\bar{\rho}\dfrac{\bar{v}_x – v_x}{\delta t}+\bar{v}_x^2 \bar{a}_\rho+2\bar{\rho} \bar{v}_x \bar{a}_x + \bar{v}_x\bar{v}_y\bar{b}_\rho+\bar{\rho}(\bar{v}_y\bar{b}_x+\bar{v}_x\bar{b}_y)+\bar{a}_p
\\
\\
0&=\bar{v}_y\dfrac{\partial \bar{\rho}}{\partial t}+\bar{\rho}\dfrac{\partial \bar{v}_y}{\partial t}+\bar{v}_y^2\bar{b}_\rho+2\bar{\rho} \bar{v}_y\bar{b}_y+\bar{v}_x\bar{v}_y\bar{a}_\rho+\bar{\rho}(\bar{v}_y\bar{a}_x+\bar{v}_x\bar{a}_y)+\bar{b}_p-\bar{\rho} g
\\
&=\bar{v}_y\dfrac{\bar{\rho} – \rho}{\delta t}+\bar{\rho}\dfrac{\bar{v}_y – v_y}{\delta t}+\bar{v}_y^2\bar{b}_\rho+2\bar{\rho} \bar{v}_y\bar{b}_y+\bar{v}_x\bar{v}_y\bar{a}_\rho+\bar{\rho}(\bar{v}_y\bar{a}_x+\bar{v}_x\bar{a}_y)+\bar{b}_p-\bar{\rho} g
\end{aligned}
\end{cases}
$$
由上面公式可以得到
(4.3-2)
$$
\begin{cases}
\begin{aligned}
\bar{a}_p &=-(\bar{v}_x\dfrac{\bar{\rho} – \rho}{\delta t}+\bar{\rho}\dfrac{\bar{v}_x – v_x}{\delta t}+\bar{v}_x^2 \bar{a}_\rho+2\bar{\rho} \bar{v}_x \bar{a}_x + \bar{v}_x\bar{v}_y\bar{b}_\rho+\bar{\rho}(\bar{v}_y\bar{b}_x+\bar{v}_x\bar{b}_y))
\\
\\
\bar{b}_p &= \bar{\rho} g – (\bar{v}_y\dfrac{\bar{\rho} – \rho}{\delta t}+\bar{\rho}\dfrac{\bar{v}_y – v_y}{\delta t}+\bar{v}_y^2\bar{b}_\rho+2\bar{\rho} \bar{v}_y\bar{b}_y+\bar{v}_x\bar{v}_y\bar{a}_\rho+\bar{\rho}(\bar{v}_y\bar{a}_x+\bar{v}_x\bar{a}_y))
\end{aligned}
\end{cases}
$$
5 过程求解
5.1 控制方程
这里使用三角形单元对平面进行单元划分，所有待求解物理量均挂在单元节点上，然后单元内的变量符合线性分布假设。
这里的边界条件有速度边界与密度(浓度)边界条件。因此求解就变成了对每个迭代步\(\delta t\)求解。
通过观察4里的离散方程: （1）对于质量方程这里考虑浓度，因此使用公式(4.1-4)来对密度进行更新，当然，如果不考虑浓度，也可以使用公式(4.1-2）来对密度进行更新; （2）速度的更新使用公式(4.2-2)；(3) 而对于压力则不直接更新，这里只更新压力梯度，因为在公式里，只要知道每个时刻压力梯度，则密度、速度就能完全关联起来，因此对于压力梯度的更新使用公式(4.3-2)。
因此，最终需要求解的变量为：密度、速度、压力梯度。
5.2 常规求解
上面的离散公式均是针对单元内某个点的离散公式，也就是说单元内的点位上的物理属性满足上面的控制方程。不失一般性，这里设单元形心点位上的物理量满足上述控制方程。而形心上的物理量，可以直接使用类似(2-6)的公式来表示，因此最终将所有方程组合在一起，将求解一个大型稀疏线性方程组。
5.3 加速递推求解
每次迭代求解时，如果均需要求解一个大型稀疏线性方程组，消耗的时间与计算成本较高。那是否可以规避求解一个大型稀疏线性方程组，而直接求解这些控制方程呢？答案是可以的。
通过观察(4.1-4)可以发现，当知道\(t_0\)时刻的物理量，那通过(4.1-4)则可以得到\(t_1\)时刻每个单元中心点的密度\(\bar{\rho}\)。注意，这时只得到\(t_1\)时刻单元中心点的密度，而单元当前密度分布系数以及各个节点的密度值并不知道。然而，单元密度有线性分布的假设，在\(t_1\)时刻，只要知道单元上任意两个节点的密度，那结合单元中心点的密度值，就可以得到单元密度的线性分布系数，进而第三个节点的密度值可以通过分布系数计算出来，也就是说，对每个单元，只要确定了单元相邻的任意两个节点上的密度，那这个单元的密度分布系数以及相连的三个节点密度就可以确定下来。只要有一个单元的两个节点有密度边界约束，那这个单元的所有节点密度就可以确定下来，而和这个单元有共同边的单元，因为共同边上的节点密度已经计算出，因此共同边所在的单元也能确定出所有节点密度，进而通过单元拓扑关系，所有连通域的单元密度均能确定下来。
当所有单元在\(t_1\)时刻的密度以及密度线性分布系数确定下来后，可以通过(4.2-2)确定出每个单元\(t_1\)时刻形心点的速度值。和密度节点更新类似，使用速度边界，然后根据拓扑关系，可以更新所有单元在\(t_1\)时刻的速度以及速度线性分布系数。
而对于\(t_1\)时刻的压力梯度，则对每个单元直接使用(4.3-2)进行更新。由于更新的是单元内压力梯度值，因此不需要再通过拓扑关系进行更新压力。
最终由\(t_0\)时刻迭代到\(t_1\)时刻的密度、速度、压力梯度值，按如下步骤计算：
step1：使用(5.3-1)更新中心密度，接着通过密度(浓度)边界条件启动计算边界节点密度，然后通过拓扑关系更新所有节点密度。
(5.3-1)
$$
\bar{\rho}=\rho + \delta t(\dfrac{c}{\eta}(a_x+b_y)+v_xa_\rho+v_yb_\rho-\dfrac{2b_2D}{\eta}(a_\rho^2+b_\rho^2))
$$
step2：使用(5.3-2)更新中心速度，接着通过速度边界条件启动计算边界节点速度，然后通过拓扑关系更新所有节点速度。
(5.3-2)
$$
\begin{cases}
\begin{aligned}
\bar{v}_x &= v_x – \dfrac{1}{\rho}(v_x(\bar{\rho} – \rho) + \delta t (v_x^2 a_\rho+2\rho v_x a_x + v_xv_yb_\rho+\rho(v_yb_x+v_xb_y)+a_p))
\\
\bar{v}_y &= v_y – \dfrac{1}{\rho}(v_y(\bar{\rho} – \rho) + \delta t(v_y^2b_\rho+2\rho v_yb_y+v_xv_ya_\rho+\rho(v_ya_x+v_xa_y)+b_p-\rho g))
\end{aligned}
\end{cases}
$$
step3：使用(5.3-3)更新压力梯度。
(5.3-3)
$$
\begin{cases}
\begin{aligned}
\bar{a}_p &=-(\bar{v}_x\dfrac{\bar{\rho} – \rho}{\delta t}+\bar{\rho}\dfrac{\bar{v}_x – v_x}{\delta t}+\bar{v}_x^2 \bar{a}_\rho+2\bar{\rho} \bar{v}_x \bar{a}_x + \bar{v}_x\bar{v}_y\bar{b}_\rho+\bar{\rho}(\bar{v}_y\bar{b}_x+\bar{v}_x\bar{b}_y))
\\
\\
\bar{b}_p &= \bar{\rho} g – (\bar{v}_y\dfrac{\bar{\rho} – \rho}{\delta t}+\bar{\rho}\dfrac{\bar{v}_y – v_y}{\delta t}+\bar{v}_y^2\bar{b}_\rho+2\bar{\rho} \bar{v}_y\bar{b}_y+\bar{v}_x\bar{v}_y\bar{a}_\rho+\bar{\rho}(\bar{v}_y\bar{a}_x+\bar{v}_x\bar{a}_y))
\end{aligned}
\end{cases}
$$
注意：
(1) 这里提到的拓扑关系求解，需要使用外部边界条件进行启动(密度速度边界)，另外需要模型在一个连通域内，并且这种求解方式：对于5.2中建立的大型非线性稀疏方程组，如果这个线性方程组是满秩的，则两种方法求到的结果是一致的。
(2) 5.3的方法不需要求解大型稀疏线性方程组，它本身是一种递推求解，因此也适合多线程加速求解，其求解效率上要比直接求解线性方程组效率高。
6 后话
力学仿真，控制方程虽然是核心，但是非规则模型、网格划分、单元预处理、网格单元自动增删、以及特殊边界处理、以及能量沙漏(固体力学)、边界反射(固体力学)、求解计算等问题，可能会被大部分研究者给忽视，但是对一个完整的仿真系统而言，有时一个细小处理过程的差异会引导着结果朝不可预期的方向发展。
参考
[1] https://users.metu.edu.tr/csert/me582/ME582%20Ch%2001.pdf
[2] https://cushman.host.dartmouth.edu/courses/engs43/Diffusion-Advection.pdf
[3] 盐水浓度与密度关系