针对双曲守恒型方程, 一些计算工作者提出了大量的数值计算格式, 如Godunov方法^[1], TVD^[2], ENO^[3], WENO^[4], DG^[5]等.所有的这些格式都是通过适当的引入数值粘性以达到计算的稳定性, 但是数值粘性会使间断被磨损了.对于非线性间断, 激波, 由于两边有压力差存在, 其磨损不是很严重, 特别是一些高精度高分辨率格式在网格足够细的情况下磨损变得非常小.然而, 对于线性间断, 一维的切向间断和二维的滑移线, 由于其两边没有压力差, 其磨损非常严重, 即使是用当今流行的格式在非常细的网格下进行计算, 其磨损还是非常大不能被人们所接受.

为克服线性间断的磨损, 人们发展了很多的方法和技巧, 例如Yang的人工压力方法^[6], Harten的subcell技巧^[7]等. Ultra-bee格式是由Roe针对线性方程提出^[8], 后由Després进一步研究和发掘^[9], 这是一个一阶的TVD格式, 其TVD限制器对应Sweby著名的TVD的外边框, 该格式一个吸引人的特点是它能精确地分辨线性间断.

Bouchut曾在文[10]中对Ultra-bee格式做了一个熵修正, 然而相当复杂, 很难将其推广到Euler方程组.最近几年来, 茅德康等针对发展偏微分方程的数值模拟提出了Entropy格式, 见文[11, 12].

李红霞, 王志刚和茅德康在文[13]中将Entropy格式和Ultra-bee格式结合起来针对线性传输方程建立了一个所谓的Entropy-Ultra-bee格式, 该格式中的数值解在各网格中由两个数值实体来刻画, 一个是数值解的网格平均$u_j^n$, 另外一个是数值熵的网格平均$U_j^n$.

本文主要利用Godunov格式、Entropy格式、Ultra-bee格式和Entropy-Ultra-bee格式计算线性传输方程.数值初值包含光滑部分和间断的方型波, 数值实验表明Entropy-Ultra-bee格式有比较高的分辨率, 而且没有非物理振荡产生.

标量守恒方程的一般形式为

$u_t+f(u)_x=0,$

(2.1)

这里$f(u)$是关于$u$的连续函数.为了方便起见, 我们假设$f''(u)\geq 0,$即关于$u$非凹.众所周知, 方程 (2.1) 的解会发生间断, 此时解函数$u$是按弱的意义下满足方程 (2.1) 的.但是方程 (2.1) 的解不唯一.为了保证解的唯一性, 弱解还必须满足如下的熵条件

$U(u)_t+F(u)_x\leq 0,$

(2.2)

其中$U(u)$是$u$的任一凸函数, $U''(u)>0$, $U(u)$称为$u$的熵函数, 而$F(u)=\displaystyle\int f'(u)U'(u)dx$成为熵流函数.此时方程 (2.2) 也是在弱的意义下满足的, 详细描述见文[14].注意, 当$u$是光滑时, 方程 (2.1) 按照经典意义满足, 此时方程 (2.2) 的等号成立, 且也是按照经典意义下满足.

当$f(u)=au$, 其中$a$为一固定常数时, 方程 (2.1) 退化为

$u_t+au_x=0.$

(2.3)

方程 (2.3) 是一个线性传输方程, 此时方程有如下形式的通解

$u(x,t)=\varphi(x-at).$

(2.4)

同样方程 (2.1) 容许间断解, 且解按照弱的意义下满足方程.此时熵条件 (2.2) 中等号成立, 也即

$U(u)_t+aU(u)_x= 0,$

(2.5)

对弱解也成立.

台阶重构的Godunov格式是对标准的Godunov格式的一个推广.如同标准的Godunov格式一样, 数值解${u_j^n}$为对$t_n$时刻精确解网格平均的近似, 即

$u_j^n\simeq\frac{1}{h}\int_{x_{j-\frac{1}{2}}}^{x_{j+\frac{1}{2}}} u(x,t_n)dx.$

(2.6)

如同标准的Godunov格式一样, 台阶重构的Godunov格式也按重构, 发展和网格平均三步来进行.

i) 重构.和标准的Godunov格式不同的是, 每一个网格内重构解是包含两片常数的阶梯函数,

$R(x;\cdot)=u_j^n+\left\{ \begin{array}{ll} -d_j^n, {x_{j-1/2}< x\leq x_j}, \\ +d_j^n, {x_j< x\leq x_{j+1/2}}, \end{array} \right.$

(2.7)

这里$d_j^n$称为重构半台阶, 我们称这样的重构为台阶重构.不难看出上述重构自动满足

$\frac{1}{h}\int_{x_{j-1/2}}^{x_{j+1/2}}R(x;\cdot)dx=u_j^n.$

(2.8)

ii) 发展.以重构函数$R(x;\cdot)$作为$t_n$层的初值, 求解以下初值问题

$\left\{\begin{array}{lc} v_{t}+f(v)_{x}=0 , -\infty < x < \infty, t_{n}< t \leq t_{n+1}, \\ v(x,t_n)=R(x;\cdot), -\infty<x<\infty. \end{array} \right.$

(2.9)

定解问题 (2.9) 的精确解为$v(x,t)$.

iii) 网格平均.在$t=t_{n+1}$时的数值解计算为

$u_{j}^{n+1}=\displaystyle\frac{1}{h}\displaystyle\int_{x_{j- \frac{1}{2}}}^{x_{j+\frac{1}{2}}} v(x,t_{n+1})dx.\quad$

(2.10)

通常情况下我们采用以下的通量形式来计算数值解

$u_{j}^{n+1}=u_{j}^{n}-\lambda(\hat{f}_{j+\frac{1}{2}}^{n}- \hat{f}_{j-\frac{1}{2}}^{n}),$

(2.11)

其中数值流函数为

${\hat f}^n_{j+1/2}=f(u^\ast(u^n_j+d_j^n,u^n_{j+1}-d_{j+1}^n)).$

(2.12)

这里, $u^\ast(v,w)$表示分别以$v$和$w$为左右状态的Riemann问题的解在$x=0$处的值.特别是当$f(u)=au$, 方程 (2.1) 退化为方程 (2.3), 此时

${\hat f}^n_{j+1/2}=\left \{ \begin{array}{ll} a(u^n_{j+1}-d^n_{j+1}), a < 0, \\ a(u^n_j+d^n_j),a \geq 0. \end{array} \right.$

(2.13)

以下介绍的所有格式都是台阶重构的Godunov格式.

标准的Godunov格式是一个半台阶$d_j^n=0$的台阶重构格式, 此时, 重构函数和$u^n$有关, 即

$R(x;u^n)=u_j^n.$

(2.14)

Entropy格式是针对线性传输方程 (2.3) 的.该格式也是一个台阶重构的Godunov格式.该格式引入了一个数值熵, 它是对解的熵的网格平均的逼近, 即

$U_{j}^{n}\simeq\frac{1}{h}\int_{x_{j-\frac{1}{2}}}^{x_{j+\frac{1}{2}}}U(u(x,t_{n}))dx.$

(2.15)

因此该格式有两个数值实体$u_j^n$和$U_j^n$.此时台阶重构函数与$u^n$和$U^n$都有关, 因此方程 (2.7) 为

$R(x;u^n,U^n)=u_j^n+\left\{ \begin{array}{ll} -d_j^{n,e}, {x_{j-1/2}< x\leq x_j}, \\ +d_j^{n,e}, {x_j< x\leq x_{j+1/2}}, \end{array} \right.$

(2.16)

这里$d_j^{n,e}$是熵的半台阶.为了计算$d_j^{n,e}$, 要求$ \frac{1}{h}\displaystyle\int_{x_{j-1/2}}^{x_{j+1/2}}U(R(x;u^n,U^n))dx=U_j^n, $即重构解的熵的网格平均等于该网格的数值熵.将式 (2.16) 代入式 (2.17) 不难得到

$\frac{1}{2}[U(u^n_j-d_j^{n,e})+U(u^n_j+d_j^{n,e})]=U^n_j.$

(2.17)

Ultra-bee格式的重构函数只和$u^n$有关, 此时方程 (2.7) 式成为

$R(x;u^n)=u_j^n+\left\{ \begin{array}{ll} -d_j^{n,ub}, {x_{j-1/2}< x\leq x_j}, \\ +d_j^{n,ub}, {x_j< x\leq x_{j+1/2}}, \end{array} \right.$

(2.18)

这里$d_j^{n,ub}$是重构的Ultra-bee半台阶.就可得$a>0$时的半台阶$d_j^{n,ub}$的计算为

$d_j^{n,ub}=\left\{ \begin{array}{ll} 0, & (u_j^n-u_{j-1}^n)(u_{j+1}^n-u_j^n)<0, \\ u_{j+1}^n-u_j^n , & |u_j^n-u_{j-1}^n|\geq\frac{\nu}{1-\nu}|u_{j+1}^n-u_j^n|, \\ \frac{1-\nu}{\nu}(u_j^n-u_{j-1}^n), & |u_j^n-u_{j-1}^n|<\frac{\nu}{1-\nu}|u_{j+1}^n-u_j^n|. \end{array} \right.$

(2.19)

当$a\leq 0$, 可得Ultra-bee台阶计算为

$d_j^{n,ub}=\left\{ \begin{array}{ll} 0, & (u_j^n-u_{j-1}^n)(u_{j+1}^n-u_j^n)<0, \\ u_j^n-u_{j-1}^n , & |u_{j+1}^n-u_j^n|\geq\frac{-\nu}{1+\nu}|u_j^n-u_{j-1}^n|, \\ \frac{1+\nu}{-\nu}(u_{j+1}^n-u_j^n), & |u_{j+1}^n-u_j^n|<\frac{-\nu}{1+\nu}|u_j^n-u_{j-1}^n|. \end{array} \right.$

(2.20)

台阶重构和$u^n$和$U^n$都有关$ R(x;u^n,U^n)=u_j^n+\left\{ \begin{array}{ll} -d_j^n, & {x_{j-1/2}< x\leq x_j}, \\ +d_j^n, & {x_j< x\leq x_{j+1/2}}, \end{array} \right. $这里$d_j^n$是重构半台阶.

为了结合Entropy格式和Ultra-bee格式的优点, 我们将格式的半台阶计算为

$d_j^n=\mbox{sgn}(u_{j+1}^n-u_{j-1}^n)\mbox{min} \{|d_j^{n,e}|,|d_j^{n,ub}|\},$

(2.21)

其中熵半台阶$d_j^{n,e}$的计算如2.3中所述, Ultra-bee的半台阶$d_j^{n,ub}$的计算如2.4中所述.

算例3.1 考虑以下初值问题

$\begin{array}{l} {u_t} + {u_x} = 0,0 \le x \le 1,\\ u(x,0) = {e^{ - 200{{(x - 0.3)}^2}}} + \left\{ {\begin{array}{*{20}{c}} 1,&{0.6 < x < 0.8,}\\ 0,&{其它.} \end{array}} \right. \end{array}$

本算例取之于文[14], 问题的精确解含一光滑波和一个方波 (两个间断).取$h=0.01$, 100个网格, CFL取为0.45. Entropy格式中取$U(u)=u^2$.圆圈表示数值解, 实线表示精确解.图 3.1-3.2显示的是在$t=1$(一个周期) 时的数值解.从图 3.1-3.2可以看出Godunov格式将解的光滑部分和间断都磨损得非常厉害, Entropy格式改善了解的光滑部分的精度, 但是在间断处出现了非物理振荡, Ultra-bee格式在间断处有比较高的分辨率, 但是在光滑区出现了梯形状结构, Entropy-Ultra-bee格式在光滑区和间断附近都有较高的分辨率, 而且没有出现非物理振荡.

Entropy-Ultra-bee格式结合了Entropy格式和Ultra-bee格式的优点, 格式严格地保证了熵守恒, 在保证格式稳定性的前提下最大程度的克服了磨损, 网格加密可以维持其过渡点不增加, 提高了整个计算区域的分辨率, 而且没有出现非物理振荡.