技术领域
[0001] 本
发明涉及计算
流体力学数值方法领域，特别是一种保持强稳定性的变步长多步法时间离散方法。
背景技术
[0002] 在可压缩流体计算中，欲模拟的物理量通常快速变化，在狭窄区域内以大梯度变化，为准确计算机模拟流场的细节，许多高分辨技术被不断地提出和完善。
[0003] 目前，求解的双曲守恒律的主流的高
精度技术为间断Galerkin有限元方法，有限体积或者有限差分加权本质无振荡方法等技术，通常的这些技术在求解欧拉方程组时需要与保强稳定性Runge-Kutta方法，或者Lax-Wendroff方法相结合。最近，有人给出Lax-Wendroff方法与间断Galerkin有限元方法等高精度方法相结合的高精度
算法设计，这些算法在计算时间效率上优于保强稳定性Runge-Kutta方法，然而算法实现比保强稳定性Runge-Kutta方法复杂。
[0004] 又有人提出易于算法实现的总变差稳定多步法离散技术，与保强稳定性Runge-Kutta方法相比，总变差稳定多步法离散技术在计算效率上可以设计出比保强稳定性Runge-Kutta方法更高效率的方法，然而现有的保强稳定性的多步法时间离散技术基于等步长，不能用于实际问题计算。
发明内容
[0005] 本发明的主要目的在于克服
现有技术中的上述
缺陷，提出一种应用于计算
流体力学的实际问题的模拟的保持强稳定性的变步长多步法时间离散方法。
[0006] 本发明采用如下技术方案：
[0007] 一种保持强稳定性的变步长多步法时间离散方法，令现有的半离散有限体积Hermite加权本质无振荡算法表示为L(Q)，并得到其对应的半离散方程 其特征在于：令tn为时间
节点，n为时间步，Δtn＝tn-tn-1；给出对于所述半离散方程的以下离散算法
[0008]其中h1＝Δtn-2，h2＝Δtn-1，h3＝Δtn，h4＝Δtn+1； CFL数可取
为1/3。
[0009] 优选的，还包括另一离散算法
[0010]其中h1＝Δtn-3，h2＝Δtn-2，h3＝Δtn-1，h4＝Δtn，h5＝Δtn+1，
[0011] CFL数可取为0.5。
[0012] 优选的，还包括另一离散算法
[0013]
[0014] 其中h1＝Δtn-4，h2＝Δtn-3，h3＝Δtn-2，h4＝Δtn-1，h5＝Δtn，h6＝Δtn+1；
[0015] CFL数可取为0.567。
[0016] 优选的，所述离散算法的初始步使用经典保强稳定的Runge-Kutta时间离散技术，即：
[0017] Q(1)＝Qn+ΔtL(Qn)
[0018]
[0019]
[0020] 初始步的时间步长可以选取正常计算的时间步长，其中n为整数。
[0021] 由上述对本发明的描述可知，与现有技术相比，本发明具有如下有益效果：
[0022] 本发明方法的时间步长可以根据发展方程的情况来自动调整时间步长；其中提供三种离散算法以供选择，而方法二和方法三的CFL数比方法一大，数值计算中CFL数越大，计算效率越高，因此这两种技术还具有更高的计算效率。另外，本发明所提出的算法可以与其它主流的守恒律的高精度技术结合。
附图说明
[0023] 图1是30°斜坡示意图；
[0024] 图2是本发明方法的整体流程示意图；
[0025] 图3是本发明方法模拟双
马赫问题的数值结果(网格2400*600等分的数值结果)；
[0026] 图4是图3的局部放大图。
具体实施方式
[0027] 以下通过具体实施方式对本发明作进一步的描述。
[0028] 本发明对于现有的半离散有限体积Hermite加权本质无振荡算法，我们用记号L(Q)表示这部分算法，记号解释L是一种算法，而Q是算法L所要处理的变量。
[0029] 一、我们叙述L的算法过程：
[0030] 算法L用于模拟无粘绝热空
气动力学，为了清楚描述所模拟的
空气动力学物理过程，我们以二维
气体动力学问题为例，我们知道问题必然满足三大物理基本定律，其数学模型为
[0031] Ut+F(U)x+G(U)y＝0；
[0032] 其中
[0033]
[0034] 上述方程组中第一方程为
质量守恒定律，第二个和第三个方程为动量守恒定律的流通量，第四个方程为
能量守恒定律。相应的变量符号说明：单位体积的总能ρ，(u，v)，P分别为流体的
密度，速度，压强；γ＝1.4为理想气体状态系数，Q为U和U的梯度组成的量，U表示为密度、动量和能量构成的向量。
[0035] 为了方便地描述模拟上述物理过程的算法，将其方程组用下面一般的符号描述双曲守恒律
[0036]
[0037] 符号说明：
[0038] 对上面的守恒律方程组作一个梯度grad(我们用在直
角坐标下说明这个符号，得到
[0039]
[0040] 其中 的含义为
[0041] 而Hermi te有限体积算法的控制方程由上面两个方程(1)和(2)组成[0042] Qt+div H(Q)＝0，
[0043] 其中Q＝(U，grad U)T，
[0044] 对上面的控制方程组在小区域Ωj(物理表面记为 )上积分，
[0045]
[0046] 再通过Gauss公式，把积分 转化为带方向的积分可以得到半离散有限体积方法
[0047]
[0048] 其中|Ωj|是区域Ωj的体积, ds可以理解为对某一物体表面细分后的面，对应的这个面的向外的方向为n。
[0049] 另 由S个 组成，即 那么对上面方程右端的积分就可以用在各个 上的离散点数为q的高斯积分离散
[0050]
[0051] 其中：Gsl和ωl是相应的高斯积分点和权重。(q点高斯积分离散指的是[0052] 对于在某一高斯点的数值流通量H(Q(Gsl，t))·n可由近似或者精确Riemann解代替(这是由不同 之间按照某一公式计算得到，这里算法不是本发明所要考虑的，不再详细叙述)。
[0053] 到此为止，这部分关于算法L(L就是公式(3)中 )叙述完了，这是因为上述过程中其实是基于 的公式过程。因此这里我们记为 简记为L(Q)。
叙述完算法L，得到一个半离散方程
[0054] 二、对于本发明提出的一种保持强稳定性的变步长多步法时间离散方法，给出了这类方法中的3种技术：
[0055] 记号说明：tn是时间节点，n为时间步，Δtn＝tn-tn-1。
[0056] 方法一、保强稳定性的变步长多步法时间离散技术：
[0057]
[0058] 其中h1＝Δtn-2，h2＝Δtn-1，h3＝Δtn，h4＝Δtn+1和 可以看到上面的时间步长h为变时间，如果Δtn=Δt,那么这一个方法变为等时间步长多步法[0059]
[0060] 这个方法的CFL数为1/3，因此本发明的多步法技术的CFL数也取为1/3。
[0061] CFL数是算法稳定相关的参数，CFL数越大，则计算时间效率越好。
[0062] 方法二、保强稳定性的变步长多步法时间离散技术：
[0063]
[0064] 其中h1＝Δtn-3，h2＝Δtn-2，h3＝Δtn-1，h4＝Δtn，h5＝Δtn+1[0065] 和 这种技术的CFL数可取为1/2，同样选取相对应等步长方法的CFL数。
[0066] 方法三、保强稳定性的变步长多步法时间离散技术：
[0067]
[0068] 其中h1＝Δtn-4，h2＝Δtn-3，h3＝Δtn-2，h4＝Δtn-1，h5＝Δtn，h6＝Δtn+1和这种技术的CFL数可取为0.567，同样选取相对应等步长方法的CFL数。
[0069] 如上面所述的方法一公式中得到Qn+1时需要Qn和Qn-3。我们计算从n＝0开始，因此无法直接使用多步法，需要使用以下的保强稳定性Runge-Kutta方法开得到Q0，Q1，Q2，Q3。
[0070]
[0071] 初始步的时间步长可以选取正常计算的时间步长，符号说明本
说明书中n作为上下标表示为整数。
[0072] 以下使用本发明方法模拟计算流体力学中的双马赫问题：
[0073] 1.问题描述：
[0074] 这是二维气体动力学间断解问题，其数学模型为
[0075] Ut+F(U)x+G(U)y＝0，
[0076] 其中
[0077]
[0078] 而单位体积的总能 其中ρ，(u，v)，P分别为流体的密度，流体速度，压强，理想气体状态系数γ＝1.4。
[0079] 初始时刻，Mach 10的正激波撞向一个30°斜坡，如图1。那么激波与斜坡之间成60°夹角。使用所发明的离散技术计算流场的演化，模拟出tprint＝0.2时刻流场的密度分布。
[0080] 2.计算机模拟参数及边界条件：
[0081] 斜坡平面作为x轴，建立平面直角
坐标系，计算区域取为[0，4]×[0，1],反射壁面处于计算区域的内部，一个Mach数为10的斜强激波放置在 y＝0处，并与x轴成60°夹角。在 之前的底边壁面采用准确的激波波后条件，其他壁面采用反射边界条件。左侧波后区域的边界条件为入流边界条件，右侧为出流边界条件。上边界为精确边界条件即给出斜激波。
[0082] 3.算法L过程：
[0083] 如图2整体
流程图的第一步，初值条件的赋值，时间层nnt＝0，tnum＝0：
[0084] 1)对计算区域划分：
[0085] 把计算区域[0，4]×[0，1]作在x轴和y轴上分别做2400等分和600等分，对划分单元编号，x轴和y轴的编号分别为i和j，则这些划分后的单元可用Iij表示，单元关于坐标轴的长高分别为Δx，Δy。
[0086] 2)赋初始条件，即赋计算量 当y＞(x-1/6)tan(π/3)时，赋波后条件，否则赋波前条件。
[0087] 每一步时间计算的边界条件：
[0088] 3)赋边界条件:边界条件给定见上面2的描述。例如上边界条件，当时，赋波后条件，否则赋波前条件。
[0089] 赋初始条件和边界条件后，求近似或者精确Riemann解,即说明中的公式(3)；到此完成了算法L(L(Q))。
[0090] 4.多步法时间离散算法过程，这里叙述方法一的过程，方法二和方法三的过程与方法一一样，不一样之处在于他们的计算公式不同：
[0091] 若采取上文中方法一，初始步数nnt为4，即需要保强稳定性Runge-Kutta方法(4)求出Q1，Q2，Q3,如图2，当步数nnt<初始步数，保强稳定性Runge-Kutta方法执行初始步。nnt为时间计算步，tnum为计算时间，dt为计算步长度，tprint为最后的计算时间。
[0092] 注意每执行完保强稳定性Runge-Kutta方法的中间步后要赋3)中的边界条件，时间更新tnum＝tnum+dt，其中dt为Δt。
[0093] 执行完初始步后，进行方法一
[0094]
[0095] 其中h1＝Δtn-2，h2＝Δtn-1，h3＝Δtn，h4＝Δtn+1和
[0096] 反复通过上面公式计算从Qn到Qn+1,直到达到
指定时间tprint＝0.2，注意每执行完公式(7)就要赋3)中的边界条件,时间更新tnum＝tnum+dt，其中dt为Δtn+1。
[0097] 结果性能分析：
[0098] 图3、图4为网格2400*600等分的数值结果，从图3可以看出发明的技术很好的捕捉了激波，从图4可以看出发明可以模拟出详细的涡结构。
[0099] 本发明的优点在于变时间步长和保强稳定性。通过实例展示了这个新算法在实例中的计算机模拟的高精度、高
分辨率以及时间高效性等优势。这种技术还可以与多种空间算法组合，适用于可压缩流
体模拟，最终能够推广到实际复杂流场的模拟等大规模科学计算问题。
[0100] 上述仅为本发明的具体实施方式，但本发明的设计构思并不局限于此，凡利用此构思对本发明进行非实质性的改动，均应属于侵犯本发明保护范围的行为。
