计算流体力学——从原理到代码（三）：非线性方程的黎曼问题与初等波（超长预警）

结论: Riemann 问题的基本解

线性常数方程组，每个波是一个速度 $s_i = \lambda_i$ 行进的间断（ $\lambda_i$ 为线性退化场）
对非线性系统，波可以是间断（激波，接触间断）和光滑过度（稀疏波）

在Riemann问题的解中，出现的波的类型，关键依赖封闭条件，（例如对Euler方程，我们只考虑状态方程，使得从出现的波仅是激波+接触间断+稀疏波）

激波：两个常数状态 $U_L, U_R$ 通过单个跳跃间断连续，且满足：

(1）RH条件 $s_i[U] = [F]$

(2) 熵条件 $\lambda_i(U_L) > s_i >\lambda_i(U_R)$

接触间断

(1) RH条件 $[F] = s_i[U]$

(2) 广义Riemann不变量，即 $\frac { d u _ { 1} } { R _ { 1} ^ { ( i ) } } = \cdots = \frac { d u _ { m } } { R _ { m } ^ { ( i ) } }$ ，且有m-1个独立的首次积分

(3) 并行特征条件 $\lambda_i(U_L) = \lambda_i(U_R) = s_i$

稀疏波

(1) ) 广义Riemann不变量 $\frac { d u _ { 1} } { R _ { 1} ^ { ( i ) } } = \cdots = \frac { d u _ { m } } { R _ { m } ^ { ( i ) } }$

(2) 特征发散 $\lambda_i(U_L) < \lambda_i(U_R)$

零、Review for Beginner（读过之前文章的，可以跳至一）

回顾，第一篇文章

徐政豪：计算流体力学——从原理到代码（一）：对流方程的格式与实现

里提到的计算流体力学——从原理到代码（一）：对流方程的格式与实现里提到的线性对流方程： $\partial _ { t } u + A \partial _ { x } u = 0$ 其中，A是一个常系数（矩阵）。

我们这里，考虑u只有一维的情况，那么A就是常数。

所有问题都要考虑自变量因变量，所以我们考虑，不同时间t上，给定空间的不同点x上，方程值u的变化。

发现，所有的u沿着x-t平面上的射线保持不变，称之为特征线： $u ( x ,t ) = u _ { 0} ( x - A t )$

此外

（默认已知内容，不懂的查资料，很多）从流体力学三大守恒（质量、动量、能量）推出基本方程： $\frac { \partial } { \partial t } \left( \begin{array} { c } { \rho } \\ { \rho u } \\ { E } \end{array} \right) + \frac { \partial } { \partial x} \left( \begin{array} { c } { \rho u } \\ { \rho u ^ { 2} + p } \\ { u ( E + p ) } \end{array} \right) = 0$

这是一个双曲拟线性守恒方程，换句话说：

单位体积微元，每时间单位体积里（质量、动量、能量）的变化 $i.e. \frac { \partial } { \partial t } \vec{u}$ ，等于体积微元与外界交换的流的变化 $i.e. \frac { \partial } { \partial t } \vec{u} + \frac { \partial } { \partial t } \vec{f}_{lux}(\vec{u}) = \vec{0}$ 。

第二篇文章里

徐政豪：计算流体力学——从原理到代码（二）：从线性对流方程到黎曼问题

我们开始接触黎曼问题，一般的流体基本方程的Riemann问题的解的形式，即当初值是有限常数状态， $\frac { \partial u } { \partial t } + \frac { \partial f ( u ) } { \partial x } = 0 ,u ( x ,0) = \left\{ \begin{array} { l l } { u _ { L } ,} & { x < 0} \\ { u _ { R } ,} & { x > 0} \end{array} \right.$ ，这就叫做Riemann问题。

为什么研究分片常数初值的PDE的解这么重要？还是那回事，我们使用离散的算法毕竟原先的连续数学问题，那么，计算机计算的时候，相当每次间断处都在求解一个分段函数初值的PDE，我们怎么相信计算机解的是“好”的，就要好好地研究CFD的Riemann问题？

对于线性对流方程，只有一个变量u的黎曼问题 $\left\{\begin{array} { l } { u = u _ { L } ,( x < 0) } \\ { u = u _ { R } ,( x \geq 0) } \end{array} \right.$ , $x=0$ 处解不连续，作特征线，特征线上解断了，如图，特征线左右值恒为ul，ur

接着，我们通过双曲性，推得了，比如上一节的最后，我们得到了一般性的流体方程Riemann问题的解： $u ( x ,t ) = \sum _ { i = 1} ^ { m } w _ { i } ^ { ( 0) } \left( x - \lambda _ { i } t \right) R ^ { ( i ) }, w ^ { ( 0) } ( x ) = R ^ { - 1} u _ { 0} ( x )$

一、走向非线性

前面我们一直考虑的是线性系统： $\partial _ { t } u + A \partial _ { x } u = 0$ 其中，A是一个常系数（矩阵）。

现在考虑守恒律的拟线性形式：

$\partial _ { t } u + f ^ { \prime } ( u ) \partial _ { x } u = 0\quad \text{ with } u ( x ,t = 0) = u _ { 0} ( x )$

它的特征线（即沿着该线解为常数）满足

$X ^ { \prime } ( t ) = f ^ { \prime } [ u ( X ( t ) ,t ) ]$

(证明由定义推一下就可以)

因为 $u(x,t)$ 在 $X(t)$ 曲线上是常数，所以当然可以积分，得到从 $X = X _ { 0} ,t = 0$ 出发的特征线为： $X ( t ) = X _ { 0} + f ^ { \prime } \left[ u _ { 0} \left( X _ { 0} \right) \right] t$ ,

如果 $f'(u)$ 是一个标量，那么只要解还是光滑的，特征线一定都是射线。

二、经典例子：Burgers 方程

PDE中大名鼎鼎的 Burgers 方程：

$\\\partial _ { t } u + \partial _ { x } \frac{u^2}{2} = 0$

特征线 $X ( t ) = X _ { 0} + u _ { 0} \left( X _ { 0} \right) t$ ，但是呢，如果 $\partial _ { x } u _ { 0} > 0$ ，特征线就会相交，但是物理量不可能相互穿过而无影响（保持光滑），另一方面，如果 $\partial _ { x } u _ { 0} < 0$ ，特征线会发散。我们具体来看

(一) 激波

于是，我们特别的，考虑他的Riemann问题 $\left \{\begin{array} { l } { u = u _ { L } = - u _ { 0} / 2,( x < 0) } \\ { u = u _ { R } = u _ { 0} / 2,( x \geq 0) } \end{array} \right.$ ，特征线相交与射线： $x=0,t>0$ ，在这里形成一个激波。

守恒律告诉我们： $\left( u _ { L } - u _ { R } \right) \Delta x = \left[ f \left( u _ { L } \right) - f \left( u _ { R } \right) \right] \Delta t$

于是，我们得到激波的波速 $s = \frac { f \left( u _ { r } \right) - f \left( u _ { l } \right) } { u _ { r } - u _ { l } }$

这个结果可以进一步而被推广到多元的方程组，我们称为Rankine-Hugoniot jump conditions（具体推导见下面第四部分）

注意时空平面（x-t）上，u值相同的射线斜率等于u值，这意味着他们之间会相互接触，波会接触，产生了激波——这就是非线性的独特之处

（二）稀疏波

另一方面，如果 $\partial _ { x } u _ { 0} < 0$ ，对于黎曼问题 $\left \{\begin{array} { l } { u = u _ { L } = - u _ { 0} / 2,( x < 0) } \\ { u = u _ { R } = u _ { 0} / 2,( x \geq 0) } \end{array} \right.$

（和前面不同之处在于，x=0左右谁高谁低换了），会在x=0这条射线上，没有任何特征线通过，只有特征线发散出去，意味着没有物理量会通过这里，只会出去。这就是稀疏波，

此时，问题的一个弱解为： $u ( x ,t ) = \left\{ \begin{array} { l l } { u _ { L } ,} & { x < u _ { L } t } \\ { x / t ,} & { u _ { L } t \leq x \leq u _ { R } t } \\ { u _ { R } ,} & { x > u _ { R } t } \end{array} \right.$

但是，可以验证，下面这个也是一个弱解： $u ( x ,t ) = \left\{ \begin{array} { l l } { u _ { L } ,} & { x < s _ { m } t } \\ { u _ { m } ,} & { s _ { m } t \leq x \leq u _ { m } t } \\ { \frac { x } { t } ,} & { u _ { m } t \leq x \leq u _ { R } t } \\ { u _ { R } ,} & { x > u _ { R } t } \end{array} \right.$ ，其中 $u _ { m } \in \left[ u _ { L } ,u _ { R } \right]$ 是任意的， $s_m = \frac { u _ { L } + u _ { m } } { 2}$

因此，守恒律的弱解是不唯一的，需要加上额外的条件，以便从众多弱解中，选取出物理意义的解，这个条件被称为熵条件，或者叫做可容许条件

注意时空平面（x-t）上，u值相同的射线斜率等于u值，这意味着他们之间会相互分散，产生了稀疏波——这就是非线性的独特之处

三、弱解与Rankine-Hugoniot条件

徐政豪：计算流体力学——从原理到代码（四）：有限体积法与Godunov Scheme 初步

里，我简单的提到了弱解的不唯一性，我们回顾一下：

之前，我们一直研究的是流体守恒方程的微分形式: $\\ \frac{\partial }{\partial { t }} u + A(u) \nabla_x u = 0$

上一节，我们讨论了非线性的流体问题，例如，Burgers equation:

$\frac{\partial}{\partial t} u + \frac{ \partial u^2/2}{\partial x} =0 \\\Leftrightarrow \frac{\partial}{\partial t} u + u\frac{\partial}{\partial x} u = 0$

在 $x=0$ 处的Riemann问题会产生解的非光滑性。而非光滑的地方，拟线性PDE则失效了。

但是这样的解，可以用积分形式来表示——因为积分形式的方程是流体守恒的更基本的描述（流体基本方程的推导过程，大家可以参考z站上的其他答案）

拟线性双曲守恒律的积分形式：在区域 $\Omega = \left[ x _ { L } ,x _ { R } \right] \times \left[ t _ { 1} ,t _ { 2} \right]$ 上

对微分形式重积分，得 $\int _ { t _ { 1} } ^ { t _ { 2} } \left[ \frac { d } { d t } \int _ { x _ { L } } ^ { x _ { R } } U ( x ,t ) d x \right] d t = \int _ { x _ { 1} } ^ { x _ { 2} } \left[ u \left( x ,t _ { 2} \right) - u \left( x ,t _ { 1} \right) \right] d x = \int _ { t _ { 1} } ^ { t _ { 2} } F \left[ u \left( x _ { 1} ,t \right) \right] - F \left[ u \left( x _ { 2} ,t \right) \right] d t\\ \text{ for any } t _ { 1} ,t _ { 2} ,x _ { 1} ,x _ { 2}$

（更严格来说，不仅仅要对区域 $\Omega = \left[ x _ { L } ,x _ { R } \right] \times \left[ t _ { 1} ,t _ { 2} \right]$ 成立，而是要对任意控制体都成立，大家明白意思就好）

由散度定理，控制体上量的变化，等于时间微元内边缘的流的变化，即

$\\\int _ { \partial \Omega } U d x - F ( U ) d t = 0$

这样的方程的解 $U(x,t)$ 被称为弱解。

此外，解的爆破我们介绍一下

在特征线发生相交之前，单值解都可以用特征线的方法给出，如前所述的 $u ( x ,t ) = u _ { 0} ( x - a ( u ) t )$

当特征线开始相交时，我们就说解爆破（Blow Up）此时解的偏导数 $\partial _x U$ 变为无穷大，爆破发生在由 $x_0$ 发散出来的特征线 $x = x_0 + a(u_0 (x_0)) t$ ，其中 x0满足： $a _ { x } \left( x _ { 0} \right) = a ^ { \prime } u _ { 0} ^ { \prime } < 0$ ，且 $| a _ { x } \left( x _ { 0} \right) |$ 是最大值的点

下面我们来建立Rankine-Hugoniot间断跳跃条件：

间断线 $x = \xi ( t )$ 上, 令 $\Gamma _ \delta$ 是下面四条曲线构成的闭回路，即在间断线的任意的一点 $\left( \xi \left( t _ { p } \right) ,t _ { p } \right)$ 附近取小邻域： $t = t _ { p } - \delta ,t = t _ { p } + \delta ,x = \xi ( t ) - \varepsilon ,x = \xi ( t ) + \varepsilon ,\delta > 0,\varepsilon > 0$

由守恒律的积分形式： $\int _ { \Gamma _ { \delta } } u d x - f d t = 0$

我们将 $\Gamma _\delta$ 拆成四条曲线，分别计算积分，得到： $\int _ { t _ { p } - \delta } ^ { t _ { p } + \delta } ( u d \xi ( t ) - f d t ) _ { x = \xi ( t )} + \varepsilon + \int _ { \xi ( t ) + \varepsilon } ^ { \xi ( t ) - \varepsilon } ( u d x - f d t ) _ { t = t _ { p } + \delta } \\ + \int _ { t _ { p } + \delta } ^ { t _ { p } - \delta } ( u d \xi ( t ) - f d t ) _ { x = \xi ( t ) - \varepsilon } + \int _ { \xi ( t ) - \varepsilon } ^ { \xi ( t ) + \varepsilon } ( u d x - f d t ) _ { t = t _ { p } - \delta } = 0$

注意到

$\int _ { \xi ( t ) \pm \varepsilon } ^ { \xi ( t ) \mp \varepsilon } ( u d x - f d t ) _ { t = t _ { p } \pm \delta }$ 中 $f d t$ 因为t是固定的，所以这一项是0

那么，上面的式子可以简化为

$\begin{array} { l } { \int _ { t _ { p } - \delta } ^ { t _ { p } - \delta } \left( u ( \xi ( t ) + \varepsilon ,t ) \frac { d \xi } { d t } - f ( u ( \xi ( t ) + \varepsilon ,t ) ) \right) d t } \\ { + \int _ { t _ { p } - \delta } \left( u ( \xi ( t ) - \varepsilon ,t ) \frac { d \xi } { d t } - f ( u ( \xi ( t ) - \varepsilon ,t ) ) \right) d t } \\ { + \int _ { \xi ( t ) - \varepsilon } ^ { \xi ( t ) + \varepsilon } u \left( x ,t _ { p } - \delta \right) d x + \int _ { \xi ( t ) + \varepsilon } ^ { \xi ( t ) - \varepsilon } u \left( x ,t _ { p } + \delta \right) d x = 0} \end{array}$

数分里的常见思想，现在令 $\epsilon \rightarrow 0$

所以

$\int _ { t _ { p } - \delta } ^ { t _ { p } + \delta } \left( u ^ { + } \frac { d \xi } { d t } - f ^ { + } \right) d t - \int _ { t _ { p } - \delta } ^ { t _ { p } + \delta } \left( u ^ { - } \frac { d \xi } { d t } - f ^ { - } \right) d t = 0$

由于 $t_p , \delta$ 是任意的，所以在间断线 $x = \xi ( t )$ 上， $\left( u ^ { + } - u ^ { - } \right) \frac { d \xi } { d t } = f ^ { + } - f ^ { - } ,or ~~\vec{ s} [ u ] = [ f ]$

其中 $s = \frac { d \xi } { d t }$ 是间断速度， $[ u ] : = u^{+} - u^-=u ( x + 0,t ) - u ( x - 0,t )$

这样我们就得到了非线性守恒律激波和接触间断处解满足的条件，即Rankine-HugonIot条件

给定 $U _ { t } + F ( U ) _ { x } = 0$ 和一个对应 $\lambda_i$ 的特征场的速度为 $\frac { d x } { d t } = s _i$ 的间断或波，则跨越间断线的RH条件：

$s _ { i } [ U ] = [ F ( U ) ]$

举个例子，根据这个关系可以找出线性常系数双曲方程的Riemann问题的解

比如说：

线性常系数双曲方程的Riemann问题的解：中间部分的密度和速度记为rho_* , u_*

对于 $U = \left( \begin{array} { c } { u _ { 1} } \\ { u _ { 2} } \end{array} \right) = \left( \begin{array} { c } { \rho } \\ { u } \end{array} \right) ,\quad A = \left( \begin{array} { c c } { 0} & { \rho _ { 0} } \\ { a ^ { 2} / \rho_ 0} & { 0} \end{array} \right)$ ，特征值 $\lambda _ { 1} = - a ,\lambda _ { 2} = a$

跨越速度为 $s _ { 1} = \lambda _ { 1}$ 的1-波的RH条件为： $\left( \begin{array} { c c } { 0} & { \rho _ { 0} } \\ { \frac { a ^ { 2} } { \rho _ { 0} } } & { 0} \end{array} \right) \left( \begin{array} { l } { \rho _ { * } - \rho _ { L } } \\ { u _ { * } - u _ { L } } \end{array} \right) = - a \left( \begin{array} { c } { \rho_ * - \rho _ { L } } \\ { u _ { * } - u _ { L } } \end{array} \right)$

由此： $u _ { * } = u _ { L } - \left( \rho * - \rho _ { R } \right) \frac { a } { \rho _ { 0} }$

同理2-波处，由RH condition： $u _ { * } = u _ { R } + \left( \rho * - \rho _ { R } \right) \frac { a } { \rho _ { 0} }$

联立二元方程组，就得到了Riemann问题的中间状态的：

$\left.\begin{array} { l } { \rho _ { * } = \frac { 1} { 2} \left( \rho _ { L } + \rho _ { R } \right) - \frac { 1} { 2} \left( u _ { R } - U _ { L } \right) \frac { \rho _ { 0} } { a } } \\ { u _ { * } = \frac { 1} { 2} \left( u _ { L } + u _ { R } \right) - \frac { 1} { 2} \left( \rho _ { R } - \rho _ { L } \right) \frac { a } { \rho _ { 0} } } \end{array} \right.$

码子码latex太累了……

没讲完的有熵条件，Riemann不变量

感觉写得太细了，大家没什么人想看……

大家了解一下什么是初等波，初等波的一些结论，以及一些定理和术语就好

这一块的东西太PDE了

就这样吧，如果想弄清细节的话，大家可以看Toro的书

来源：知乎 www.zhihu.com

作者：知乎用户（登录查看详情）

【知乎日报】千万用户的选择，做朋友圈里的新鲜事分享大牛。
点击下载