关于 Euler 方程组的一些常用的基本知识。

控制方程与状态方程

一维的可压缩流体满足 Euler 方程组，主要由三个控制方程组成

\[ \begin{pmatrix} \rho \\ \rho u \\ E \end{pmatrix}_t + \begin{pmatrix} \rho u \\ \rho u^2 + p \\ u(E + p) \end{pmatrix}_x = 0 \]

分别对应质量守恒、动量守恒和能量守恒。涉及的记号含义为：

  • 密度（density） \(\rho=\rho(x,t)\)；
  • 速度（velocity）\(u=u(x,t)\)；
  • 压强（pressure）\(p=p(x,t)\)；
  • 单位体积的总能量（total energy）\(E=E(x,t)\)。

这里还缺少一个方程与控制方程组一起组成封闭的方程组。

可压缩流体的单位体积总能量\(E\)包括动能和内能两部分 \[ E = \frac12 \rho u^2 + \rho e \] 其中\(e\)是比内能（specific internal energy），由流体的压强和密度决定：\(e = e(\rho,p)\)。

对于理想气体，比内能和总能量的具体形式为 \[ \begin{aligned} e &= e(\rho,p) = \frac{p}{(\gamma-1)\rho}\\ E &= \frac12 \rho u^2 + \frac{p}{\gamma-1} \end{aligned} \] 这里 \(\gamma = c_p/c_v\) 是绝热指数（等压热容和等容热容的比值），通常取为 \(\gamma = \frac75 = 1.4\)。

因此添加下面的状态方程，就可以得到封闭的方程组 \[ E = \frac12 \rho u^2 + \frac{p}{\gamma-1} \]

一些物理记号

为了下文推导的方便，我们还需要一些具有物理背景的概念与记号。

在前面我们引入的几个物理量：密度、速度、压强、能量和比内能，除了速度允许用正负值来代表方向之外，其它的物理量都是必须大于零的，有时等于零在物理上也是合理的，但是相对于一种极端情形，我们不作考虑。

在可压缩流体中的声速（sound speed）为 \[ c = \sqrt{\frac{\gamma\, p}{\rho}} \] 这里根式内部的值在物理上一定是非负的，可以进行开方。但是这在数值计算中是不一定的，我们可能需要用保正限制器来专门抑制数值振荡，确保开方操作的合法性。

与总能量/比内能类似，在物理上还有总焓和比焓的概念，通常用\(H\)代表总焓，用\(h\)代表比焓，具体表达式为 \[ \begin{aligned} h &= h(\rho,p) = e + \frac{p}{\rho} \\ H &= \frac12 u^2 + h = \frac{E+p}{\rho} \end{aligned} \]

一个更重要的概念是熵（entropy），用记号\(s\)表示 \[ s = c_v \ln(\frac{p}{\rho^\gamma}) + C_0 \] 我们考虑的理想气体的（无黏）流动实际上是等熵流，这意味着 \[ \frac{D s}{Dt} = s_t + u s_x = 0 \] 这个方程并不独立，它可以通过前几个控制方程和状态方程推出，也可以用来替换能量守恒方程。 等熵流表明：沿着流体粒子的运动路径，我们有如下的守恒量（等熵定律） \[ p = C \rho^\gamma \] 这里的 \(C = C(s_0)\)由初始点的熵\(s_0\)决定，只要流动保持平滑，沿着路径的熵就是恒定的。反例：流体粒子在穿过激波时不满足上述规律。

守恒与非守恒形式

一维 Euler 方程组是典型的守恒律方程组，可以整理为如下的守恒形式

\[ \mathbf{U}_t + \mathbf{F}(\mathbf{U})_x = 0 \]

其中守恒变量\(U=(\xi_1,\xi_2,\xi_3)^T\)，对应的变量称为原始变量。 守恒变量与原始变量间的转换如下： \[ \xi_1 = \rho, \xi_2 = \rho u, \xi_3 = E \]

记 Jacobi 矩阵\(\mathbf{A}(\mathbf{U}) = \nabla \mathbf{F}(\mathbf{U})\)（注意这里的求导是对守恒分量进行的），我们可以将方程改写为非守恒形式

\[ \mathbf{U}_t + \mathbf{A}(\mathbf{U}) \mathbf{U}_x = 0 \]

Euler 方程组的 Jacobi 矩阵具体为

\[ \mathbf{A}(\mathbf{U}) = \begin{pmatrix} 0 & 1 & 0 \\ \frac12(\gamma-3)u^2 & (3-\gamma)u & \gamma - 1 \\ \frac12(\gamma-1)u^3 - H u & H - (\gamma-1)u^2 & \gamma u \end{pmatrix} \]

Euler 方程组是标准的双曲模型，Jacobi 矩阵可以进行相似对角化，我们将特征值记作\(\lambda_1 \le \lambda_2 \le \lambda_3\)，组成的对角阵记作

\[ \Lambda(\mathbf{U})= \begin{pmatrix} \lambda_1 & & \\ & \lambda_2 & \\ & & \lambda_3 \end{pmatrix} = \begin{pmatrix} u-c & & \\ & u & \\ & & u+c \end{pmatrix} \]

Euler 方程组的（右）特征向量和左特征向量具体为

\[ \begin{aligned} \mathbf{R}(\mathbf{U}) &= \begin{pmatrix} \mathbf{r}_1(\mathbf{U}) & \mathbf{r}_2(\mathbf{U}) & \mathbf{r}_3(\mathbf{U}) \end{pmatrix} = \begin{pmatrix} 1 & 1 & 1 \\ u-c & u & u+c \\ H - u c & \frac12 u^2 & H + u c \end{pmatrix} \\ \mathbf{L}(\mathbf{U}) &= \begin{pmatrix} \mathbf{l}_1(\mathbf{U}) \\ \mathbf{l}_2(\mathbf{U}) \\ \mathbf{l}_3(\mathbf{U}) \end{pmatrix} = \frac{1}{c^2} \begin{pmatrix} \frac12 u c + \frac14(\gamma-1)u^2 & -\frac12(\gamma-1)u -\frac12 c & \frac12(\gamma-1) \\ c^2-\frac12(\gamma-1)u^2 & (\gamma-1)u & 1-\gamma \\ -\frac12 u c + \frac14(\gamma-1)u^2 & -\frac12(\gamma-1)u + \frac12 c & \frac12(\gamma-1) \end{pmatrix} \end{aligned} \]

注：（右）特征向量即满足 \(\mathbf{A} \mathbf{r}_i = \lambda_i \mathbf{r}_i\) 的列向量，左特征向量即满足 \(\mathbf{l}_i \mathbf{A} = \lambda_i \mathbf{l}_i\) 的行向量，并且易知 \(\mathbf{L} \mathbf{R} = \mathbf{I}\)。

最终，可以得到 Jacobi 矩阵的对角化分解

\[ \mathbf{A}(\mathbf{U}) = \mathbf{R}(\mathbf{U}) \mathbf{\Lambda}(\mathbf{U}) \mathbf{L}(\mathbf{U}) \]

守恒律方程组可以基于特征分解进行变形（注意到\(\mathbf{R}^{-1}=\mathbf{L}\)）

\[ \begin{aligned} \mathbf{U}_t + \mathbf{A} \mathbf{U}_x &= 0\\ \mathbf{U}_t + \mathbf{R} \mathbf{\Lambda} \mathbf{L} \mathbf{U}_x &= 0\\ \mathbf{L} \mathbf{U}_t + \mathbf{\Lambda} \mathbf{L} \mathbf{U}_x &= 0\\ \end{aligned} \]

基于局部冻结系数的思想，在局部近似地将\(\mathbf{R},\mathbf{\Lambda},\mathbf{L}\)视为常矩阵，那么我们就可以引入特征变量\(\mathbf{V} = \mathbf{L} \mathbf{U}\)，特征变量满足解耦的、线性化的方程组

\[ \mathbf{V}_t + \mathbf{\Lambda} \mathbf{V}_x = 0 \]

每一个特征分量分别满足一个线性方程

\[ v^{(i)}_t + \lambda_i v^{(i)}_x = 0,\, (i=1,2,3) \]

对于 Euler 方程组，可以理解为在\(x\)位置，存在三个波分别以速度\(u-c,u,u+c\)进行传播，其中\(u\)是流体自身的速度，而\(c>0\)是流体中的声速。 这里进行一个简单的分类讨论：

  • 流速始终保持\(|u|<c\)，这是通常的情景，此时\(u-c<0\)视作向左传播，\(u+c>0\)始终向右传播；
  • 流速始终保持\(|u|>c\)，超声速情景，此时\(u-c,u,u+c\)的传播方向都和\(u\)的方向相同，三个波同向传播；
  • 其它情形，即流速可能跨声速，此时的情景更加复杂，不做讨论。

边界条件

考虑一维有界区域\((x_l,x_r)\)的 Euler 方程组，我们需要给予适定的边界条件：

  • 从方程组的角度，与边界处的三个特征值正负有关；
  • 从物理含义角度，与流体速度和声速的关系有关。

以左边界\(x=x_l\)为例，我们根据\(x=x_l\)处的特征值正负进行分类讨论，并且可以按照特征分量对应的波方程进行理解：

  1. \(0<u-c \le u \le u+c\)，三个波都是向右传播，需要三个边界条件；（超声速入口）
  2. \(u-c \le 0 < u \le u+c\)，一个波向左传播，两个波向右传播，需要两个边界条件；（亚声速入口）
  3. \(u-c \le u \le 0 < u+c\)，两个波向左传播，一个波向右传播，需要一个边界条件；（亚声速出口）
  4. \(u-c \le u \le u+c \le 0\)，三个波都是向左传播，不需要边界条件。（超声速出口）

对于右边界\(x=x_r\)，类似可得：

  1. \(0<u-c \le u \le u+c\)，三个波都是向右传播，不需要边界条件；（超声速出口）
  2. \(u-c \le 0 < u \le u+c\)，一个波向左传播，两个波向右传播，需要一个边界条件；（亚声速出口）
  3. \(u-c \le u \le 0 < u+c\)，两个波向左传播，一个波向右传播，需要两个边界条件；（亚声速入口）
  4. \(u-c \le u \le u+c \le 0\)，三个波都是向左传播，需要三个边界条件。（超声速入口）

通过局部特征分解，我们可以确定需要的边界条件个数。

为了简化问题，我们假设流体的边界在任意时刻始终只处于一个分类情景中，通常给定守恒分量的精确值作为边界条件：

  • 如果\(x=x_l\)需要一个边界条件，给定第一个守恒分量的值 \[ U^{(1)}(x_l,t) = \rho = g^{(1)}(t) \]
  • 如果\(x=x_l\)需要两个边界条件，给定前两个守恒分量的值 \[ \begin{aligned} U^{(1)}(x_l,t) &= \rho = g^{(1)}(t)\\ U^{(2)}(x_l,t) &= \rho u = g^{(2)}(t) \end{aligned} \]
  • 如果\(x=x_l\)需要三个边界条件，给定所有守恒分量的值 \[ \begin{aligned} U^{(1)}(x_l,t) &= \rho = g^{(1)}(t)\\ U^{(2)}(x_l,t) &= \rho u = g^{(2)}(t)\\ U^{(3)}(x_l,t) &= E = g^{(3)}(t)\\ \end{aligned} \]

还有一些更复杂的情景，给出的边界条件可能不是关于守恒变量的，而是关于压强或者温度的约束，略。

精度测试算例

这里列出几个常见的验证数值格式的数值精度的经典算例。 由于 Euler 方程组是非线性方程组，一般情形下的解析解很难获取，除了使用右端加源项的处理，还可以考虑某些特殊的退化情形，此时我们可以获取解析解或精确解。

例一（线性方程）

计算区域 \(\Omega = [x_l,x_r]\)，使用周期性边界条件， 初值和精确值为

\[ \left\{ \begin{aligned} \rho(x,t) &= \alpha + \beta \sin(\omega(x-u_0 t) + \phi)\\ u(x,t) &= u_0 \\ p(x,t) &= p_0 \\ \end{aligned} \right. \]

其中的参数 \(\alpha,\beta,\omega,\rho,u_0,p_0\)可以自由选取，要求 \(\omega >0\) 并且满足 \(\omega(x_r-x_l)=2\pi\)。

由于物理上\(\rho \ge 0\)，因此必须保持\(\alpha > |\beta|\)。

声速

\[ \begin{aligned} c = \sqrt{\gamma\frac{p}{\rho}} \in [ \sqrt{\frac{\gamma p_0}{\alpha + |\beta|}}, \sqrt{\frac{\gamma p_0}{\alpha - |\beta|}} ] \end{aligned} \]

可以与常量\(u=u_0\)比较，得到三个特征值的正负，通常可以选取

\[ 0 < u_0 < \sqrt{\frac{\gamma p_0}{\alpha + |\beta|}} \]

使得下式恒成立

\[ u-c < 0 < u < u+c \]

对于周期边界，我们则不需要注意边界条件以及特征值的正负分类。

将此算例的初值或精确解代入 Euler 方程组可得

\[ \begin{aligned} \rho_t + u_0 \rho_x &= 0 \\ u_0 (\rho_t + u_0 \rho_x) &= 0\\ \frac12 u_0^2 (\rho_t + u_0 \rho_x) &= 0 \end{aligned} \]

相当于全部退化为了同一个线性方程，因此我们可以获取解析解。

例二（Burgers 方程）

为了凑出一个退化为 Burgers 方程的特例，我们取常数 \(\gamma=3\)，假设满足如下关系

\[ \left\{ \begin{aligned} u(x,t) &= \sqrt{3} \rho(x,t) \\ p(x,t) &= \rho(x,t)^3\\ \end{aligned} \right. \]

代入可得

\[ \left\{ \begin{aligned} \rho_t + (\sqrt{3}\rho^2)_x &= 0\\ \sqrt{3} (\rho^2)_t + 4(\rho^3)_x &= 0\\ 2(\rho^3)_t + 3\sqrt{3}(\rho^4)_x &= 0 \end{aligned} \right. \]

整理易知，三个式子都等价于

\[ \rho_t + 2\sqrt{3} \rho \rho_x = 0 \]

这个式子就是 Burgers 方程的简单变形：取\(\rho = \frac{1}{2\sqrt{3}}\mu\)代换，得到

\[ \mu_t + \mu \mu_x = 0 \]

因此\(\mu(x,t)\)就是 Burgers 方程的解（可以基于初值和特征线方法求得）。

我们可以基于 Burgers 方程的求解给出一组欧拉方程组的精确解：取计算区域\(\Omega=[x_l,x_r]\)，取\(\mu\)满足三角函数初值的 Burgers 方程

\[ \left\{ \begin{aligned} \mu_t + \left(\frac{\mu^2}2\right)_x &= 0\\ \mu(x,0) &= \alpha + \beta \sin(\omega x+\phi) \end{aligned} \right. \]

周期性边界条件，系数要求\(\omega>0\)，满足\(\omega(x_r-x_l)=2\pi\)。 由此衍生出欧拉方程组的初值和精确解

\[ \left\{ \begin{aligned} \rho(x,t) &= \frac{1}{2\sqrt{3}} \mu(x,t)\\ u(x,t) &= \sqrt{3} \rho(x,t) = \frac12 \mu(x,t)\\ p(x,t) &= \rho(x,t)^3 = \frac{1}{24\sqrt{3}} \mu(x,t)^3\\ \end{aligned} \right. \]

相应的最早出现间断的时刻\(T_b\)也根据 Burgers 方程来判断，使用\(T_b\)之前的时刻进行精度测试。

由于物理上要求\(\rho \ge 0\)，因此 Burgers 方程的\(\mu\)并不能是任意的，需要保持 \(\mu(x,t)>0\)。

关于特征值的处理更加麻烦，声速

\[ \begin{aligned} c^2 = \frac{\gamma p}{\rho} = \frac{3\rho^3}{\rho} = 3 \rho^2 \end{aligned} \]

因此恰有 \(u=c\)，处处有

\[ 0 = u-c < u < u+c \]

一个特征值恒为零，处理上可能有麻烦？验证一下

\[ \begin{aligned} u = \sqrt{3} \rho, p = \rho^3, E = 2 \rho^3 \\ c = \sqrt{3} \rho, H = 3 \rho^2 \end{aligned} \]

此时

\[ \mathbf{A}(\mathbf{U}) = \begin{pmatrix} 0 & 1 & 0 \\ \frac12(\gamma-3)u^2 & (3-\gamma)u & \gamma - 1 \\ \frac12(\gamma-1)u^3 - H u & H - (\gamma-1)u^2 & \gamma u \end{pmatrix} = \begin{pmatrix} 0 & 1 & 0 \\ 0 & 0 & 2 \\ 0 & -3\rho^2 & 3\sqrt{3}\rho \end{pmatrix} \]

三个特征值为

\[ \Lambda(\mathbf{U})= \begin{pmatrix} \lambda_1 & & \\ & \lambda_2 & \\ & & \lambda_3 \end{pmatrix} = \begin{pmatrix} 0 & & \\ & \sqrt{3}\rho & \\ & & 2\sqrt{3}\rho \end{pmatrix} \]

两个矩阵

\[ \mathbf{R}(\mathbf{U}) = \begin{pmatrix} 1 & 1 & 1 \\ 0 & \sqrt{3}\rho & 2\sqrt{3}\rho \\ 0 & \frac32 \rho^2 & 6\rho^2 \end{pmatrix} , \mathbf{L}(\mathbf{U}) = \frac{1}{3\rho^2} \begin{pmatrix} 3\rho^2 & -\frac{3\sqrt{3}}2\rho & 1 \\ 0 & 2\sqrt{3}\rho & -2 \\ 0 & -\frac{\sqrt{3}}2\rho & 1 \end{pmatrix} \]