一种保持强稳定性的变步长多步法时间离散算法

摘要

本发明提出了一种保持强稳定性的变步长多步法时间离散算法，其包括半离散Hermite本质无振荡有限体积算法的构造，初始步使用了三阶保强稳定的Runge-Kutta算法；本发明给出了两类时间离散的接口处理，特别是初始步的时间步长的选取。本发明的优点在于变时间步长和保强稳定性。通过实例展示了这个新算法在实例中的计算机模拟的高精度、高分辨率以及时间高效性等优势。这种技术还可以与多种空间算法组合，适用于可压缩流体模拟，最终能够推广到实际复杂流场的模拟等大规模科学计算问题。

说明书

技术领域

本发明涉及计算流体力学数值方法领域，特别是一种保持强稳定性的变步长多步法时间离散算法。

背景技术

在可压缩流体计算中，欲模拟的物理量通常快速变化，在狭窄区域内以大梯度变化，为准确计算机模拟流场的细节，许多高分辨技术被不断地提出和完善。

目前，求解的双曲守恒律的主流的高精度技术为间断Galerkin有限元方法，有限体积或者有限差分加权本质无振荡方法等技术，通常的这些技术在求解欧拉方程组时需要与保强稳定性Runge-Kutta方法，或者Lax-Wendroff方法相结合。最近，有人给出Lax-Wendroff方法与间断Galerkin有限元方法等高精度方法相结合的高精度算法设计，这些算法在计算时间效率上优于保强稳定性Runge-Kutta方法，然而算法实现比保强稳定性Runge-Kutta方法复杂。

又有人提出易于算法实现的总变差稳定多步法离散技术，与保强稳定性Runge-Kutta方法相比，总变差稳定多步法离散技术在计算效率上可以设计出比保强稳定性Runge-Kutta方法更高效率的方法，然而现有的保强稳定性的多步法时间离散技术基于等步长，不能用于实际问题计算。

发明内容

本发明的主要目的在于克服现有技术中的上述缺陷，提出一种应用于计算流体力学的实际问题的模拟的保持强稳定性的变步长多步法时间离散算法。

本发明采用如下技术方案：

一种保持强稳定性的变步长多步法时间离散算法，令现有的半离散有限体积Hermite加权本质无振荡算法表示为L(Q)，并得到其对应的半离散方程其特征在于：令t_n为时间节点，n为时间步，Δt_n＝t_n-t_n-1；给出对于所述半离散方程的以下离散算法

$>\left(\begin{array}{c}{Q}^{n+1}=\frac{(1+{\lambda }_2+{\lambda }_3-2{\lambda }_4){(1+{\lambda }_2+{\lambda }_3+{\lambda }_4)}^2}{{(1+{\lambda }_2+{\lambda }_3)}^3}{Q}^n+\frac{{(1+{\lambda }_2+{\lambda }_3+{\lambda }_4)}^2}{{(1+{\lambda }_2+{\lambda }_3)}^2}{h}_{4}L\left(Q^n\right)\\ +\frac{{\lambda }_4^2(3+3{\lambda }_2+3{\lambda }_3+2{\lambda }_4)}{{(1+{\lambda }_2+{\lambda }_3)}^3}{Q}^{n-3}+\frac{{\lambda }_4^2(1+{\lambda }_2+{\lambda }_3+{\lambda }_4)}{{(1+{\lambda }_2+{\lambda }_3)}^2}{h}_{1}L\left(Q^{n-3}\right)\end{array}\right);$>其中h₁＝Δt_n-2，h₂＝Δt_n-1，h₃＝Δt_n，h₄＝Δt_n+1；CFL数可取为1/3。

优选的，还包括另一离散算法

$>\left(\begin{array}{c}{Q}^{n+1}\text{=}\frac{(1+{\lambda }_2+{\lambda }_3+{\lambda }_4-2{\lambda }_5){(1+{\lambda }_2+{\lambda }_3+{\lambda }_4+{\lambda }_5)}^2}{{(1{+\lambda }_2+{\lambda }_3+{\lambda }_4)}^3}{Q}^n\\ +\frac{{(1+{\lambda }_2+{\lambda }_3+{\lambda }_4+{\lambda }_5)}^2}{{(1+{\lambda }_2+{\lambda }_3+{\lambda }_4)}^2}{h}_{5}L\left(Q^n\right)+\frac{{\lambda }_5^2(3+3{\lambda }_2+3{\lambda }_3+3{\lambda }_4+2{\lambda }_5)}{{(1+{\lambda }_2+{\lambda }_3+{\lambda }_4)}^3}{Q}^{n-4}\\ +\frac{{\lambda }_5^2(1+{\lambda }_2+{\lambda }_3+{\lambda }_4+{\lambda }_5)}{{(1+{\lambda }_2+{\lambda }_3+{\lambda }_4)}^2}{h}_{1}L\left(Q^{n-4}\right)\end{array}\right);$>其中

h₁＝Δt_n-3，h₂＝Δt_n-2，h₃＝Δt_n-1，h₄＝Δt_n，h₅＝Δt_n+1，

$>{\lambda }_2=\frac{h_2}{h_1},{\lambda }_3=\frac{h_3}{h_1},{\lambda }_4=\frac{h_4}{h_1},{\lambda }_5=\frac{h_5}{h_1};$>CFL数可取为0.5。

优选的，还包括另一离散算法

$>\left(\begin{array}{c}{Q}^{n+1}\text{=}\frac{(1+{\lambda }_2+{\lambda }_3+{\lambda }_4-2{\lambda }_5-2{\lambda }_6){(1+{\lambda }_2+{\lambda }_3+{\lambda }_4+{\lambda }_5+{\lambda }_6)}^2}{{(1{+\lambda }_2+{\lambda }_3+{\lambda }_4+{\lambda }_5)}^3}{Q}^n\\ +\frac{{(1+{\lambda }_2+{\lambda }_3+{\lambda }_4+{\lambda }_5+{\lambda }_6)}^2}{{(1+{\lambda }_2+{\lambda }_3+{\lambda }_4+{\lambda }_5)}^2}{h}_{6}L\left(Q^n\right)+\frac{{\lambda }_6^2(3+3{\lambda }_2+3{\lambda }_3+3{\lambda }_4+3{\lambda }_5+2{\lambda }_6)}{{(1+{\lambda }_2+{\lambda }_3+{\lambda }_4+{\lambda }_5)}^3}{Q}^{n-5}\\ +\frac{{\lambda }_6^2(1+{\lambda }_2+{\lambda }_3+{\lambda }_4+{\lambda }_5+{\lambda }_6)}{{(1+{\lambda }_2+{\lambda }_3+{\lambda }_4+{\lambda }_5)}^2}{h}_{1}L\left(Q^{n-5}\right)\end{array}\right);$>

其中h₁＝Δt_n-4，，h₂＝Δt_n-3，h₃＝Δt_n-2，h₄＝Δt_n-1，h₅＝Δt_n，h₆＝Δt_n+1；

$>{\lambda }_2=\frac{h_2}{h_1},{\lambda }_3=\frac{h_3}{h_1},{\lambda }_4=\frac{h_4}{h_1},{\lambda }_5=\frac{h_5}{h_1},{\lambda }_6=\frac{h_6}{h_1};$>CFL数可取为0.567。

优选的，所述离散算法的初始步使用经典保强稳定的Runge-Kutta时间离散技术，即：

$>\left(\begin{array}{c}{Q}^{\left(1\right)}=Q^n+\mathrm{\Delta tL}\left(Q^{\left(n\right)}\right)\\ Q^{\left(2\right)}=\frac{3}{4}{Q}^n+\frac{1}{4}(Q^{\left(1\right)}+\mathrm{\Delta tL}\left(Q^{\left(1\right)}\right))\\ Q^{n+1}=\frac{1}{3}{Q}^n+\frac{2}{3}(Q^{\left(2\right)}+\mathrm{\Delta tL}\left(Q^{\left(2\right)}\right))\end{array}\right)$>

初始步的时间步长可以选取正常计算的时间步长，其中n为整数。

由上述对本发明的描述可知，与现有技术相比，本发明具有如下有益效果：

本发明算法的时间步长可以根据发展方程的情况来自动调整时间步长；其中提供三种离散算法以供选择，而方法二和方法三的CFL数比方法一大，数值计算中CFL数越大，计算效率越高，因此这两种技术还具有更高的计算效率。另外，本发明所提出的算法可以与其它主流的守恒律的高精度技术结合。

附图说明

图1是30°斜坡示意图；

图2是本发明算法的整体流程示意图；

图3是本发明算法模拟双马赫问题的数值结果(网格2400*600等分的数值结果)；

图4是图3的局部放大图。

具体实施方式

以下通过具体实施方式对本发明作进一步的描述。

本发明对于现有的半离散有限体积Hermite加权本质无振荡算法，我们用记号L(Q)表示这部分算法，记号解释L是一种算法，而Q是算法L所要处理的变量。

一、我们叙述L的算法过程：

算法L用于模拟无粘绝热空气动力学，为了清楚描述所模拟的空气动力学物理过程，我们以二维气体动力学问题为例，我们知道问题必然满足三大物理基本定律，其数学模型为

U_t+F(U)_x+G(U)_y＝0；

其中

$>U=\left(\begin{array}{c}\rho \\ \mathrm{\rho u}\\ \mathrm{\rho v}\\ E\end{array}\right),F\left(U\right)\text{=}\left(\begin{array}{c}\mathrm{\rho u}\\ \mathrm{\rho uu}+P\\ \mathrm{\rho uv}\\ u(E+P)\end{array}\right),G\left(U\right)=\left(\begin{array}{c}\mathrm{\rho v}\\ \mathrm{\rho vu}\\ \mathrm{\rho vv}+P\\ v(E+P)\end{array}\right).$>

上述方程组中第一方程为质量守恒定律，第二个和第三个方程为动量守恒定律，第四个方程为能量守恒定律。相应的变量符号说明：单位体积的总能ρ，(u，v)，P分别为流体的密度，速度，压强；γ＝1.4为理想气体状态系数。

为了方便地描述模拟上述物理过程的算法，将其方程组用下面一般的符号描述双曲守恒律

符号说明：div(ab)＝a_x+b_y。

对上面的守恒律方程组作一个梯度grad(我们用在直角坐标下说明这个符号，$>\mathrm{grad}a=\left(\begin{array}{c}{a}_x\\ a_y\end{array}\right),$>得到

其中的含义为$>\left(\begin{array}{c}{a}_1\\ a_2\end{array}\right)\otimes \left(\begin{array}{c}{b}_1\\ b_2\end{array}\right)=\left(\begin{array}{cc}{a}_1{b}_1& a_1{b}_2\\ a_2{b}_1& a_2{b}_2\end{array}\right).$>

而Hermite有限体积算法的控制方程由上面两个方程(1)和(2)组成

Q_t+div H(Q)＝0，

其中Q＝(U，grad U)^T，

对上面的控制方程组在小区域Ω_j(物理表面记为)上积分，

$>\frac{1}{\left|{\Omega }_j\right|}{\int }_{{\Omega }_j}{Q}_t+\frac{1}{\left|{\Omega }_j\right|}{\int }_{{\Omega }_j}\mathrm{div}H\left(Q\right)=0;$>

再通过Gauss公式，把积分转化为带方向的积分可以得到半离散有限体积方法

$>\frac{d}{\mathrm{dt}}{\bar{Q}}_{{\Omega }_j}=-\frac{1}{\left|{\Omega }_j\right|}{\int }_{\partial {\Omega }_j}H\left(Q\right)·\mathrm{nds}---\left(3\right)$>

其中|Ω_j|是区域Ω_j的体积,ds可以理解为对某一物体表面细分后的面，对应的这个面的向外的方向为n。

另由S个组成，即那么对上面方程右端的积分就可以用在各个上的离散点数为q的高斯积分离散

$>{\int }_{\partial {\Omega }_j}H\left(Q\right)·\mathrm{nds}\approx \underset{s=1}{\overset{S}{\Sigma }}|\partial {\Omega }_{\mathrm{js}}|\underset{l=1}{\overset{q}{\Sigma }}{\omega }_{l}H\left(Q(G_{\mathrm{sl}},t)\right)·n;$>

其中：G_sl和ω_l是相应的高斯积分点和权重。(q点高斯积分离散指的是$>{\int }_{\partial {\Omega }_{\mathrm{js}}}H\left(Q\right)·\mathrm{nds}\approx |\partial {\Omega }_{\mathrm{js}}|\underset{l=1}{\overset{q}{\Sigma }}{\omega }_{l}H\left(Q(G_{\mathrm{sl}},t)\right)·n$>)

对于在某一高斯点的数值流通量H(Q(G_sl，t))·n可由近似或者精确Riemann解代替(这是由不同之间按照某一公式计算得到，这里算法不是本发明所要考虑的，不再详细叙述)。

到此为止，这部分关于算法L(L就是公式(3)中)叙述完了，这是因为上述过程中其实是基于的公式过程。因此这里我们记为简记为L(Q)。叙述完算法L，得到一个半离散方程

二、对于本发明提出的一种保持强稳定性的变步长多步法时间离散算法，给出了这类方法中的3种技术：

记号说明：t_n是时间节点，n为时间步，Δt_n＝t_n-t_n-1。

方法一、保强稳定性的变步长多步法时间离散技术：

$>\left(\begin{array}{c}{Q}^{n+1}=\frac{(1+{\lambda }_2+{\lambda }_3-2{\lambda }_4){(1+{\lambda }_2+{\lambda }_3+{\lambda }_4)}^2}{{(1+{\lambda }_2+{\lambda }_3)}^3}{Q}^n+\frac{{(1+{\lambda }_2+{\lambda }_3+{\lambda }_4)}^2}{{(1+{\lambda }_2+{\lambda }_3)}^2}{h}_{4}L\left(Q^n\right)\\ +\frac{{\lambda }_4^2(3+3{\lambda }_2+3{\lambda }_3+2{\lambda }_4)}{{(1+{\lambda }_2+{\lambda }_3)}^3}{Q}^{n-3}+\frac{{\lambda }_4^2(1+{\lambda }_2+{\lambda }_3+{\lambda }_4)}{{(1+{\lambda }_2+{\lambda }_3)}^2}{h}_{1}L\left(Q^{n-3}\right).\end{array}\right)---\left(4\right)$>

其中h₁＝Δt_n-2，h₂＝Δt_n-1，h₃＝Δt_n，h₄＝Δt_n+1和可以看到上面的时间步长h为变时间，如果Δt_n＝Δt,那么这一个方法变为等时间步长多步法

$>Q^{n+1}=\frac{16}{27}{Q}^n+\frac{16}{9}\mathrm{\Delta tL}\left(Q^n\right)+\frac{11}{27}{Q}^{n-3}+\frac{4}{9}L\left(Q^{n-3}\right).$>

这个方法的CFL数为1/3，因此本发明的多步法技术的CFL数也取为1/3。

CFL数是算法稳定相关的参数，CFL数越大，则计算时间效率越好。

方法二、保强稳定性的变步长多步法时间离散技术：

$>\left(\begin{array}{c}{Q}^{n+1}\text{=}\frac{(1+{\lambda }_2+{\lambda }_3+{\lambda }_4-2{\lambda }_5){(1+{\lambda }_2+{\lambda }_3+{\lambda }_4+{\lambda }_5)}^2}{{(1{+\lambda }_2+{\lambda }_3+{\lambda }_4)}^3}{Q}^n\\ +\frac{{(1+{\lambda }_2+{\lambda }_3+{\lambda }_4+{\lambda }_5)}^2}{{(1+{\lambda }_2+{\lambda }_3+{\lambda }_4)}^2}{h}_{5}L\left(Q^n\right)+\frac{{\lambda }_5^2(3+3{\lambda }_2+3{\lambda }_3+3{\lambda }_4+2{\lambda }_5)}{{(1+{\lambda }_2+{\lambda }_3+{\lambda }_4)}^3}{Q}^{n-4}\\ +\frac{{\lambda }_5^2(1+{\lambda }_2+{\lambda }_3+{\lambda }_4+{\lambda }_5)}{{(1+{\lambda }_2+{\lambda }_3+{\lambda }_4)}^2}{h}_{1}L\left(Q^{n-4}\right)\end{array}\right)---\left(5\right)$>

其中h₁＝Δt_n-3，h₂＝Δt_n-2，h₃＝Δt_n-1，h₄＝Δt_n，h₅＝Δt_n+1和这种技术的CFL数可取为1/2，同样选取相对应等步长方法的CFL数。

方法三、保强稳定性的变步长多步法时间离散技术：

$>\left(\begin{array}{c}{Q}^{n+1}\text{=}\frac{(1+{\lambda }_2+{\lambda }_3+{\lambda }_4-2{\lambda }_5-2{\lambda }_6){(1+{\lambda }_2+{\lambda }_3+{\lambda }_4+{\lambda }_5+{\lambda }_6)}^2}{{(1{+\lambda }_2+{\lambda }_3+{\lambda }_4+{\lambda }_5)}^3}{Q}^n\\ +\frac{{(1+{\lambda }_2+{\lambda }_3+{\lambda }_4+{\lambda }_5+{\lambda }_6)}^2}{{(1+{\lambda }_2+{\lambda }_3+{\lambda }_4+{\lambda }_5)}^2}{h}_{6}L\left(Q^n\right)+\frac{{\lambda }_6^2(3+3{\lambda }_2+3{\lambda }_3+3{\lambda }_4+3{\lambda }_5+2{\lambda }_6)}{{(1+{\lambda }_2+{\lambda }_3+{\lambda }_4+{\lambda }_5)}^3}{Q}^{n-5}\\ +\frac{{\lambda }_6^2(1+{\lambda }_2+{\lambda }_3+{\lambda }_4+{\lambda }_5+{\lambda }_6)}{{(1+{\lambda }_2+{\lambda }_3+{\lambda }_4+{\lambda }_5)}^2}{h}_{1}L\left(Q^{n-5}\right)\end{array}\right)---\left(6\right)$>

其中h₁＝Δt_n-4，h₂＝Δt_n-3，h₃＝Δt_n-2，h₄＝Δt_n-1，h₅＝Δt_n，h₆＝Δt_n+1和这种技术的CFL数可取为0.567，同样选取相对应等步长方法的CFL数。

如上面所述的方法一公式中得到Qⁿ⁺¹时需要Qⁿ和Q^n-3。我们计算从n＝0开始，因此无法直接使用多步法，需要使用以下的保强稳定性Runge-Kutta方法开得到Q⁰，Q¹，Q²，Q³。

$>\left(\begin{array}{c}{Q}^{\left(1\right)}=Q^n+\mathrm{\Delta tL}\left(Q^{\left(n\right)}\right)\\ Q^{\left(2\right)}=\frac{3}{4}{Q}^n+\frac{1}{4}(Q^{\left(1\right)}+\mathrm{\Delta tL}\left(Q^{\left(1\right)}\right))\\ Q^{n+1}=\frac{1}{3}{Q}^n+\frac{2}{3}(Q^{\left(2\right)}+\mathrm{\Delta tL}\left(Q^{\left(2\right)}\right))\end{array}\right)---\left(7\right)$>

初始步的时间步长可以选取正常计算的时间步长，符号说明本说明书中n作为上下标表示为整数。

以下使用本发明算法模拟计算流体力学中的双马赫问题：

1.问题描述：

这是二维气体动力学间断解问题，其数学模型为

U_t+F(U)_x+G(U)_y＝0，

其中

$>U=\left(\begin{array}{c}\rho \\ \mathrm{\rho u}\\ \mathrm{\rho v}\\ E\end{array}\right),F\left(U\right)\text{=}\left(\begin{array}{c}\mathrm{\rho u}\\ \mathrm{\rho uu}+P\\ \mathrm{\rho uv}\\ u(E+P)\end{array}\right),G\left(U\right)=\left(\begin{array}{c}\mathrm{\rho v}\\ \mathrm{\rho vu}\\ \mathrm{\rho vv}+P\\ v(E+P)\end{array}\right).$>

而单位体积的总能其中ρ，(u，v)，P分别为流体的密度，流体速度，压强，理想气体状态系数γ＝1.4。

初始时刻，Mach 10的正激波撞向一个30°斜坡，如图1。那么激波与斜坡之间成60°夹角。使用所发明的离散技术计算流场的演化，模拟出tprint＝0.2时刻流场的密度分布。

2.计算机模拟参数及边界条件：

斜坡平面作为x轴，建立平面直角坐标系，计算区域取为[0，4]×[0，1],反射壁面处于计算区域的内部，一个Mach数为10的斜强激波放置在处，并与x轴成60°夹角。在之前的底边壁面采用准确的激波波后条件，其他壁面采用反射边界条件。左侧波后区域的边界条件为入流边界条件，右侧为出流边界条件。上边界为精确边界条件即给出斜激波。

3.算法L过程：

如图2整体流程图的第一步，初值条件的赋值，时间层nnt＝0，tnum＝0：

1)对计算区域划分：

把计算区域[0，4]×[0，1]作在x轴和y轴上分别做2400等分和600等分，对划分单元编号，x轴和y轴的编号分别为i和j，则这些划分后的单元可用I_ij表示，单元关于坐标轴的长高分别为Δx，Δy。

2)赋初始条件，即赋计算量当y＞(x-1/6)tan(π/3)时，赋波后条件，否则赋波前条件。

每一步时间计算的边界条件：

3)赋边界条件:边界条件给定见上面2的描述。例如上边界条件，当时，赋波后条件，否则赋波前条件。

赋初始条件和边界条件后，求近似或者精确Riemann解,即说明中的公式(3)；到此完成了算法L(L(Q))。

4.多步法时间离散算法过程，这里叙述方法一的过程，方法二和方法三的过程与方法一一样，不一样之处在于他们的计算公式不同：

若采取上文中方法一，初始步数nnt为4，即需要保强稳定性Runge-Kutta方法(4)求出Q¹，Q²，Q³,如图2，当步数nnt<初始步数，保强稳定性Runge-Kutta方法执行初始步。nnt为时间计算步，tnum为计算时间，dt为计算步长度，tprint为最后的计算时间。

注意每执行完保强稳定性Runge-Kutta方法的中间步后要赋3)中的边界条件，时间更新tnum＝tnum+dt，其中dt为Δt。

执行完初始步后，进行方法一

$>\left(\begin{array}{c}{Q}^{n+1}=\frac{(1+{\lambda }_2+{\lambda }_3-2{\lambda }_4){(1+{\lambda }_2+{\lambda }_3+{\lambda }_4)}^2}{{(1+{\lambda }_2+{\lambda }_3)}^3}{Q}^n+\frac{{(1+{\lambda }_2+{\lambda }_3+{\lambda }_4)}^2}{{(1+{\lambda }_2+{\lambda }_3)}^2}{h}_{4}L\left(Q^n\right)\\ +\frac{{\lambda }_4^2(3+3{\lambda }_2+3{\lambda }_3+2{\lambda }_4)}{{(1+{\lambda }_2+{\lambda }_3)}^3}{Q}^{n-3}+\frac{{\lambda }_4^2(1+{\lambda }_2+{\lambda }_3+{\lambda }_4)}{{(1+{\lambda }_2+{\lambda }_3)}^2}{h}_{1}L\left(Q^{n-3}\right).\end{array}\right)$>

其中h₁＝Δt_n-2，h₂＝Δt_n-1，h₃＝Δt_n，h₄＝Δt_n+1和

反复通过上面公式计算从Qⁿ到Qⁿ⁺¹,直到达到指定时间tprint＝0.2，注意每执行完公式(7)就要赋3)中的边界条件,时间更新tnum＝tnum+dt，其中dt为Δt_n+1。

结果性能分析：

图3、图4为网格2400*600等分的数值结果，从图3可以看出发明的技术很好的捕捉了激波，从图4可以看出发明可以模拟出详细的涡结构。

本发明的优点在于变时间步长和保强稳定性。通过实例展示了这个新算法在实例中的计算机模拟的高精度、高分辨率以及时间高效性等优势。这种技术还可以与多种空间算法组合，适用于可压缩流体模拟，最终能够推广到实际复杂流场的模拟等大规模科学计算问题。

上述仅为本发明的具体实施方式，但本发明的设计构思并不局限于此，凡利用此构思对本发明进行非实质性的改动，均应属于侵犯本发明保护范围的行为。