一种通过高阶WENO格式降阶对可压缩流动问题进行数值模拟的方法与流程

本发明涉及可压缩流动的数值模拟技术领域，尤其是一种通过高阶weno格式降阶对可压缩流动问题进行数值模拟的方法。

背景技术：

在可压缩流动的数值模拟中，可能同时存在间断和其他连续的小尺度结构，高阶精度weno(weightedessentiallynon-oscillatory，加权本质无振荡)格式既可以捕捉间断，又能在光滑区域达到高阶精度，常用来模拟此类问题。通常更高阶的格式能在流场连续区域获得更精细的流动结构，并且能更加锐利地捕捉激波。

但是，对于强非线性间断问题，采用高阶weno格式更容易出现非物理震荡，甚至可能计算失败。原因之一是高阶格式的耗散更低，因而可能无法有效地抑制振荡。其二是它们所用的较宽模板可能包含多于一个的间断。那么，所有子模板可能都不光滑，因此无法得到本质无振荡的解。其三是高阶插值或重构中的runge现象。偏离中心的宽模板上的数值通量可能导致较大的数值误差，进而影响稳定性。对于欧拉方程，由于不同特征方向的相互作用，情况会变得更糟。非常高阶格式的解中可能会出现严重的数值振荡，并且计算更可能会发散。对于非常高阶weno格式，现有方法采用递归降阶来解决这个问题，该方法需要密度和压强在单元边界处的重构值。然而，高阶weno差分格式的重构通常是在通量上进行的。因此，该方法不能直接应用到守恒型差分格式中，而需要构造其它方法。

技术实现要素：

本发明所要解决的技术问题是：针对上述存在的问题，提供一种通过高阶weno格式降阶对可压缩流动问题进行数值模拟的方法。

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

一种通过高阶weno格式降阶对可压缩流动问题进行数值模拟的方法，包括如下步骤：

s1，计算某一位置的数值通量时，对于所给的最高精度weno格式，根据一个准则判断是否需要降阶；

s2，如果需要降阶，则对于降阶后的次高精度weno格式，根据步骤s1中的同一准则判断是否需要降阶；

s3，循环执行步骤s2，直到不需要降阶或者到达最低阶精度的weno格式；然后采用相应精度的weno格式计算该位置的数值通量，以进一步计算得到该位置对含激波的可压缩流动问题的数值模拟值。

进一步的，记所求数值通量的位置为xj+1/2，weno格式用到的模板为s＝{xj-m+1,xj-m+2,…,xj+m}，则步骤s1与步骤s2中的所述准则为同时满足如下两个要求：

(1)当模板点数2m≥4时，压强的光滑因子需满足：

β^(2m)≤4β^(2m-2)；

光滑因子β^(2m)的计算需用到模板s上所有2m个点；

光滑因子β^(2m-2)的计算则用到模板s去掉两个端点后的所有2m-2个点；

(2)压强p和速度u需满足：

i∈(j-m+1,…,j+m-1),l∈(j-m+1,…,j+m-2)；

其中，γ为气体比热比；h及δt分别为当前方向空间网格尺度和时间步长；d为所求问题的空间维数。

进一步的，所述光滑因子β^(2m)或光滑因子β^(2m-2)的计算公式如下：

其中，β⁽ⁿ⁾表示光滑因子β^(2m)或光滑因子β^(2m-2)；

xc为模板的中点，即xj+1/2；

σ,δ的定义如下：

其中，h为当前方向空间网格尺度。

进一步的，所述光滑因子β^(2m)或光滑因子β^(2m-2)的计算公式如下：

β⁽ⁿ⁾＝(bn)²+|an·cn|

其中，β⁽ⁿ⁾表示光滑因子β^(2m)或光滑因子β^(2m-2)；

xc为模板的中点，即xj+1/2。

进一步的，an,bn,cn均可写为如下形式：

其中，k为光滑因子所用模板最左端点的下标，z为an,bn或cn；系数αn,l的取值如下：

进一步的，n＞3，且n为整数。

综上所述，由于采用了上述技术方案，本发明的有益效果是：

1、本发明可以在连续区域采用非常高阶的weno格式，从而获得更精细的流场结构；

2、本发明在强激波附近会启用相对低阶的weno格式，因此具有良好的激波稳定性，即激波附近出现的非物理振荡的可能性更小，以及计算中更不易由于出现负压而导致计算失败；

3、本发明可用于守恒型差分格式，并且由于预先确定所用格式，因此计算量较小。

附图说明

为了更清楚地说明本发明实施例的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，应当理解，以下附图仅示出了本发明的某些实施例，因此不应被看作是对范围的限定，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他相关的附图。

图1为本发明通过高阶weno格式降阶对可压缩流动问题进行数值模拟的方法的流程框图。

图2为本发明算例一中网格点数为n＝200，计算终止时间t＝1.8时的几种格式计算结果的对比图。

图3为本发明算例二中网格点数为n＝400，计算终止时间t＝0.038时的几种格式计算结果的对比图。

图4a为本发明算例三中计算区域是[0,1]×[0,1]，网格点数为n＝400，计算终止时间t＝0.8时，weno5-m格式的计算结果的密度云图。

图4b为本发明算例三中计算区域是[0,1]×[0,1]，网格点数为n＝400，计算终止时间t＝0.8时，porweno9-s格式的计算结果的密度云图。

图4c为本发明算例三中计算区域是[0,1]×[0,1]，网格点数为n＝400，计算终止时间t＝0.8时，porweno17-s格式的计算结果的密度云图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明，即所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。通常在此处附图中描述和示出的本发明实施例的组件可以以各种不同的配置来布置和设计。因此，以下对在附图中提供的本发明的实施例的详细描述并非旨在限制要求保护的本发明的范围，而是仅仅表示本发明的选定实施例。基于本发明的实施例，本领域技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例，都属于本发明保护的范围。

说明：对于欧拉方程的守恒型差分方法求解，通常的数值求解步骤如下：

步骤1，网格划分，根据计算要求和计算条件生成网格；

步骤2，赋初值，根据初值条件在网格点上给出节点初值；

步骤3，设置时间步长，根据问题和所使用方法结合cfl条件得到；

步骤4，开始循环如下时间推进过程，直到问题所要求的终止时间：

(1)对于内点，采用空间格式进行数值模拟计算；

(2)根据边界条件，采用合适的数值边界处理；

(3)内点和边界的数值计算空间离散完成后，采用合适的时间格式进行推进。

本发明主要涉及步骤4中第(1)分步骤采用空间格式进行数值模拟计算的方法，由此本发明公开了一种通过高阶weno格式降阶对可压缩流动问题进行数值模拟的方法，包括如下步骤：

s1，计算某一位置的数值通量时，对于所给的最高精度weno格式，根据一个准则判断是否需要降阶；

s2，如果需要降阶，则对于降阶后的次高精度weno格式，根据步骤s1中的同一准则判断是否需要降阶；

s3，循环执行步骤s2，直到不需要降阶或者到达最低阶精度的weno格式；然后采用相应精度的weno格式计算该位置的数值通量，以进一步计算得到该位置对含激波的可压缩流动问题的数值模拟值。

以一维欧拉方程为例对本发明的主要过程进行说明。一维欧拉方程的控制方程为：

ut+f(u)x＝0；

u＝(ρ,ρu,e)^t，f＝(ρu,ρu²+p,u(e+p))^t；

其中，ρ,u,p,e分别为密度、速度、压强与总能；状态方程为：

其中，γ为气体比热比。一维欧拉方程的半离散守恒型差分形式为：

其中，为数值通量，可由各种重构格式(例如本发明涉及的weno格式)获得。

记所求数值通量的位置为xj+1/2，weno格式用到的模板为s＝{xj-m+1,xj-m+2,…,xj+m}，则步骤s1与步骤s2中的所述准则为同时满足如下两个要求：

(1)当模板点数2m≥4时，压强的光滑因子需满足：

β^(2m)≤4β^(2m-2)；

光滑因子β^(2m)的计算需用到模板s上所有2m个点；

光滑因子β^(2m-2)的计算则用到模板s去掉两个端点后的所有2m-2个点；

该要求(1)基于模板上压强的光滑程度，能有效抑制非物理振荡的出现。

(2)压强p和速度u需满足：

i∈(j-m+1,…,j+m-1),l∈(j-m+1,…,j+m-2)；

其中，γ为气体比热比；h及δt分别为当前方向空间网格尺度和时间步长；d为所求问题的空间维数。该要求(2)基于压强的控制方程，能有效防止计算过程中负压的出现。

设采用的一组weno格式包括：九阶的weno9-s、七阶的weno7-s和五阶的weno5-m格式。则本发明的的主要过程为：在计算每一个数值通量时，先取m＝5，根据准则判断是否同时满足上述两个要求，若满足则采用weno9-s格式；否则再令m＝4，再次根据准则判断是否同时满足上述两个要求，若满足则采用weno7-s，否则采用weno5-m格式。

所述光滑因子β^(2m)或光滑因子β^(2m-2)的计算公式可以采用以下两种：

第一种：

其中，β⁽ⁿ⁾表示光滑因子β^(2m)或光滑因子β^(2m-2)；

xc为模板的中点，即xj+1/2；

σ,δ的定义如下：

其中，h为当前方向空间网格尺度。

第二种：

β⁽ⁿ⁾＝(bn)²+|an·cn|

其中，β⁽ⁿ⁾表示光滑因子β^(2m)或光滑因子β^(2m-2)；

xc为模板的中点，即xj+1/2。

对于第二种，进一步地，an,bn,cn均可写为如下形式：

其中，k为光滑因子所用模板最左端点的下标，z为an,bn或cn；系数αn,l的取值如下：

实际上，该表中n的取值可以是n＞3，且n为整数，可以推论相应的参数值，不应仅以该表限定本发明。

以下通过具体算例对本发明的特征和性能作进一步的详细描述。计算中时间格式采用3阶tvd龙格库塔方法，空间格式所用的weno格式组合为weno(2r-1)-s、weno(2r-3)-s、…、weno7-s、weno5-m，为了便于说明，将本发明通过高阶weno格式降阶构造的格式以最高阶为(2r-1)阶格式记为porweno(2r-1)-s格式。作为对比，同时还给出了weno5-m的结果。计算中还采用了局部lax–friedrichs分裂方法和特征投影技术，另外cfl数统一取0.25。

算例一：激波与熵波相互作用问题

该问题的控制方程为一维欧拉方程，一个mach数为3的激波与熵波相互干扰，诱发高频振动。图2中给出了网格点数为n＝200，计算终止时间t＝1.8时的几种格式计算结果的对比图，由于该问题没有解析解，这里采用n＝8000时经典五阶weno格式的计算结果作为参考精确解，由图2中实线“exact”表示；w5表示weno5-m的计算结果，w(2r-1)表示porweno(2r-1)-s的结果，(2r-1)≥9。从图2中可以看出，本发明构造的porweno(2r-1)-s格式在高频振荡区获得了较好的结果。

算例二：双爆轰波问题

该问题的控制方程为一维欧拉方程，其初值条件为：

此问题对格式的稳定性要求较高，直接采用非常高阶的weno格式会出现计算失败。在采用了本发明的降阶方法后，这些格式都得到了稳定可靠的结果。图2中给出了网格点数为n＝400，计算终止时间t＝0.038时的几种格式计算结果的对比图，两端采用反射边界条件。由于该问题没有解析解，这里采用n＝8000时经典五阶weno格式的计算结果作为参考精确解，由图3中实线“exact”表示；w5表示weno5-m的计算结果，w(2r-1)表示porweno(2r-1)-s的结果，(2r-1)≥9。尽管不同格式的数值结果非常接近，可以在x＝0.78附近看到它们的差距。从这里的放大图中，可以看到本发明构造的porweno(2r-1)-s格式获得了比weno5-m格式更高的峰值。

算例三：二维黎曼问题

该问题的控制方程为二维欧拉方程，其初值条件为

图4a、4b、4c分别给出了计算区域是[0,1]×[0,1]，网格点数为n＝400，计算终止时间t＝0.8时，weno5-m格式、porweno9-s和porweno17-s的计算结果的密度云图。porweno9-s和porweno17-s比weno5-m格式能获得更细致的流场结构。

通过上述内容可知，本发明具有的有益效果如下：

1、本发明可以在连续区域采用非常高阶的weno格式，从而获得更精细的流场结构；

2、本发明在强激波附近会启用相对低阶的weno格式，因此具有良好的激波稳定性，即激波附近出现的非物理振荡的可能性更小，以及计算中更不易由于出现负压而导致计算失败；

3、本发明可用于守恒型差分格式，并且由于预先确定所用格式，因此计算量较小。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。