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

结论: Riemann 问题的基本解

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

在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

线性对流方程的黎曼问题,特征线左边恒为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(具体推导见下面第四部分)

非线性例子:Burgers方程,可以发现在最中间处产生了激波!
注意时空平面(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}

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

非线性例子:Burgers方程,的稀疏波例子

注意时空平面(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

作者:知乎用户(登录查看详情)

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