一种三角函数框架下新WENO格式构造方法与流程

本发明属于计算流体力学工程技术领域，具体涉及了一种三角函数框架下新weno格式构造方法。

背景技术

在工程应用中，流场问题常常出现，比如气体动力系统和浅水建模等。因此，制定解决这类问题的鲁棒的、精确的、高效的数值模拟方法至关重要，也吸引了很多研究者的兴趣。1959年，godunov为解流场问题提出了一阶精度的数值模拟格式。一阶精度的数值模拟方法在捕捉激波时不会出现非物理的数值振荡但会过度抹平强间断，而往往强间断对问题的后续研究有着重要意义，因此需引进高精度数值计算格式模拟强间断类问题。

为了提高格式的精度，模拟解的结构以及准确的捕捉激波位置，harten于1983年首次提出了tvd(totalvariationdiminishing)格式，并在此基础上与osher于1987年提出了eno(essentiallynon-oscillatory)高精度格式。eno格式的主要思想是在逐次扩展的模板中选用最光滑的模板构造多项式求出单元边界处的值，进而在光滑区域达到高阶精度、高分辨率，同时在间断附近实现基本无振荡的效果。但是，在方法的实现过程中，eno格式最终只选用所有候选模板中最优的模板，造成计算结果的浪费，且构造的精度越高浪费的越多，导致计算效率不高。因此，liu，osher和chan等于1994年提出了weno(weightedessentiailynon-oscillatory)格式，提高了计算结果的利用率并且使得r阶精度的eno格式提高到r+1阶精度。1996年，jiang和shu进一步改善了weno格式，使得数值精度能够提高到2r-1阶，并设计出新光滑因子和非线性权的构造框架。weno格式的主要思想是通过低阶重构通量的线性凸组合获得高阶近似。但该经典weno格式的实现过程中，线性权依赖于母模板，且其求解过程相当复杂，因此，2016年zhu和qiu改善了该weno格式，在维持精度不减的情况下，随机选取大于零且总和为一的线性权。这些格式已被成功地用到很多应用领域，特别是包含激波和复杂解结构的问题，比如模拟可压缩湍流系统和空气声学系统等。

波类和高频振荡类问题在工程应用中常常出现。然而，对模拟此类问题的更适合的三角函数多项式插值weno格式的研究较少。虽然baron于1976年研究三角函数插值呈现了neville类方法，muhlbach提出了牛顿三角函数插值，但是这些成果不能直接应用于eno类插值格式。为此，christofi于1996年提出了能直接用于eno格式中的三角函数重构方法。zhu和qiu于2010年提出了用三角函数多项式重构weno格式的方法，但计算复杂，不易实现。

技术实现要素：

本发明针对现有技术中的不足，提供一种三角函数框架下新weno格式构造方法，能针对各种可压流场问题，进行高精度数值模拟。

为实现上述目的，本发明采用以下技术方案：

一种三角函数框架下新weno格式构造方法，在笛卡尔坐标系下，利用tweno格式对可压流场问题进行数值模拟，其特征在于，包括以下步骤：

步骤一、把双曲守恒律方程离散为空间半离散的有限差分格式，采用tweno格式重构通量的近似值；

步骤二、对控制方程中的时间导数使用三阶tvdrunge-kutta离散公式将半离散有限差分格式离散成时空全离散有限差分格式；

步骤三、根据时空全离散有限差分格式得到下一时间层上的近似值，依次迭代，得到计算区域内终止时刻流场的数值模拟值。

为优化上述技术方案，采取的具体措施还包括：

所述步骤一中，双曲守恒律方程为：

其半离散格式的形式为：

其中，u＝(ρ，ρu，e)^t表示守恒变量，f(u)＝(ρu，ρu²+p，u(e+p))^t表示通量，ut表示u对t求导，f(u)x表示f(u)对x求导，t表示时间变量，x表示空间变量，ρ、u、p、e分别表示流体密度、速度、压强、能量，t表示转置，u0表示初始状态值，l(u)表示-fx(u)的空间离散形式；

把空间离散成统一长度的网格单元单元长度单元中心为其中i为坐标序号，有：

其中，和分别表示通量f(u)在目标网格单元ii的边界和处的五阶近似的数值通量，ui(t)表示u在目标网格单元ii内点xi处的值u(xi，t)。

所述步骤一中，求通量f(u)在目标网格单元ii的边界和处的五阶近似值和具体步骤如下：

步骤1、采用lax-friedrichs分裂把通量分裂为其中，

步骤2、将目标网格单元ii以及其周围共五个网格单元组成一个大模板t1＝[ii-2，ii-1，ii，ii+1，ii+2]，从大模板中选择两个包含两个单元的小模板t2＝[ii-1，ii]和t3＝[ii，ii+1]；

步骤3、在t1、t2、t3每个模板上分别重构三角函数多项式p1(x)、p2(x)和p3(x)，使得：

p2(x)，p3(x)∈span{1，sin(x-xi)}；

步骤4、任意取三组线性权：

γ1＝0.98，γ2＝0.01，γ3＝0.01；

γ1＝1/3，γ2＝1/3，γ3＝1/3；

γ1＝0.01，γ2＝0.495，γ3＝0.495；

步骤5、计算光滑指示器βl，用于衡量重构多项式pl(x)在目标单元上的光滑度，计算公式为：

其中，l＝1，2，3表示对应模板序号，表示多项式pl(x)对x的α阶导数，r1＝4，r2＝1，r3＝1；

步骤6、通过线性权γl和光滑指示器βl计算非线性权ωl，其计算公式为：

其中，τ为计算过程中的过渡值，ε＝10^-6；

步骤7、求出数值通量分裂f⁺(u)在点处的tweno重构值：

同理，求出数值通量分裂f^-(u)在点处的tweno重构值、数值通量分裂f⁺(u)在点处的tweno重构值、数值通量分裂f^-(u)在点处的tweno重构值；

将计算结果代入含有时间导数项的半离散有限差分格式，得到关于时间导数的常微分方程。

所述步骤3中，具体步骤如下：

步骤3.1、在三个模板t1、t2和t3上分别构造三角函数多项式p1(x)、p2(x)和p3(x)，使其满足：

步骤3.2、得到每个模板上的三角函数插值多项式p1(x)、p2(x)和p3(x)，如下：

其中，ii-2、ii-1、ii、ii+1、ii+2分别表示第i-2、i-1、i、i+1、i+2个网格单元，分别表示f⁺(u)在点xi-2、xi-1、xi、xi+1、xi+2的值。

所述步骤二中，利用三阶tvdrunge-kutta离散公式：

得到时空全离散有限差分格式，其中，u⁽¹⁾，u⁽²⁾为中间过渡值，δt为时间步长，上标n表示第n时间层，l(u^m)，l(u⁽¹⁾)，l(u⁽²⁾)为-fx(u)的高阶空间离散形式的近似值。

所述步骤三中，时空全离散有限差分格式为关于时间层的迭代公式，初始状态值已知，通过迭代公式求出下一时间层的近似值，依次得到终止时刻计算区域内的数值模拟值。

本发明的有益效果是：相比于weno格式，该tweno格式通过把三角函数多项式而不是代数多项式作为有限差分tweno格式的构建模块，模拟了波类和高频振荡类可压流场问题，同时在光滑区域能够达到高阶精度；相比于已有三角函数多项式重构格式，该tweno格式得到的全局l¹截断误差与l^∞截断误差更小，同时也避免了在强激波和接触间断处产生非物理振荡，该新五阶tweno格式中的相关线性权不再需要通过复杂的计算得到而是被设为和为一的任意正数，因此该新tweno格式具有更简便更易拓展到高维空间的优势。

附图说明

图1a-1c是实施例一中的台阶问题，利用本发明的有限差分tweno格式得到的密度等值线图，分别采用线性权①、②、③。

图2a-2c是实施例二中的双马赫问题，利用本发明的有限差分tweno格式得到的密度等值线图，分别采用线性权①、②、③。

图3a-3c是实施例三中的激波和涡流相互干扰问题，利用本发明的有限差分tweno格式得到的t＝0.35时的压强等值线图，分别采用线性权①、②、③。

图4a-4c是实施例三中的激波和涡流相互干扰问题，利用本发明的有限差分tweno格式得到的t＝0.6时的压强等值线图，分别采用线性权①、②、③。

图5a-5c是实施例三中的激波和涡流相互干扰问题，利用本发明的有限差分tweno格式得到的t＝0.8时的压强等值线图，分别采用线性权①、②、③。

图6a-6c是实施例四中初值条件为(19)的二维euler黎曼问题，利用本发明的有限差分tweno格式得到的t＝0.25时的密度等值线图，分别采用线性权①、②、③。

图7a-7c是实施例四中初值条件为(20)的二维euler黎曼问题，利用本发明的有限差分tweno格式得到的t＝0.25时的密度等值线图，分别采用线性权①、②、③。

图8a-8c是实施例四中初值条件为(21)的二维euler黎曼问题，利用本发明的有限差分tweno格式得到的t＝0.3时的密度等值线图，分别采用线性权①、②、③。

图9a-9c是实施例四中初值条件为(22)的二维euler黎曼问题，利用本发明的有限差分tweno格式得到的t＝0.2时的密度等值线图，分别采用线性权①、②、③。

图10a-10c是实施例四中初值条件为(23)的二维euler黎曼问题，利用本发明的有限差分tweno格式得到的t＝0.3时的密度等值线图，分别采用线性权①、②、③。

具体实施方式

现在结合附图对本发明作进一步详细的说明。

本发明给出了笛卡尔网格下解可压流场问题的新型五阶有限差分tweno高精度数值计算格式的构造过程，相比于经典wen0格式，该tweno格式通过把重构的三角函数多项式而不是重构的代数多项式作为有限差分weno格式的构建模块，解决了波类和高频振荡类的可压流场问题的数值模拟，且在光滑区域能够达到高阶精度近似，捕捉到尖锐和无振荡激波的转换。新的tweno格式采用的线性权不再需要通过繁冗的数值计算得到，可设为满足和为一的任意正数，此格式方法简单、精度高，易于推广到多维空间中。该方法在笛卡尔坐标系下，利用tweno格式对可压流场问题进行数值模拟，具体步骤如下：

一、把双曲守恒律方程离散为空间半离散的有限差分格式，采用tweno格式重构通量的近似值。

考虑一维双曲守恒律方程：

其半离散格式的形式为：

其中，u＝(ρ，ρu，e)^t表示守恒变量，f(u)＝(ρu，ρu²+p，u(e+p))^t表示通量，ut表示u对t求导，f(u)x表示f(u)对x求导，t表示时间变量，x表示空间变量，ρ、u、p、e分别表示流体密度、速度、压强、能量，t表示转置，u0表示初始状态值，l(u)表示-fx(u)的空间离散形式。

把空间离散成统一长度的网格单元单元长度单元中心为其中i为坐标序号，有：

其中，和分别表示通量f(u)在目标单元ii的边界和处的五阶近似的数值通量，ui(t)表示u在网格单元ii内点xi处的值u(xi，t)。

求通量f(u)在目标单元ii的边界和处的五阶近似值和具体步骤如下：

步骤1、用最简单的lax-friedrichs分裂把通量分裂为其中为了简单起见，本发明只描述f⁺(u)在点处的重构过程并将其定义为

步骤2、将目标单元ii以及其周围共5个网格单元组成一个大模板t1＝[ii-2，ii-1，ii，ii+1，ii+2]，从大模板中选择两个包含两个单元的小模板t2＝[ii-1，ii]和t3＝[ii，ii+1]，其中ii为对应序号的网格单元。

步骤3、在每个模板上分别重构三角函数多项式p1(x)、p2(x)和p3(x)，使得：

p2(x)，p3(x)∈span{1，sin(x-xi)}。

其具体过程如下：

步骤3.1、在三个模板t1、t2和t3上分别构造三角函数多项式p1(x)、p2(x)和p3(x)，使其满足：

步骤3.2、得到每个模板上的三角函数插值多项式p1(x)、p2(x)和p3(x)，如下：

其中，ii-2、ii-1、ii、ii+1、ii+2分别表示第i-2、i-1、i、i+1、i+2个单元，分别表示f⁺(u)在点xi-2、xi-1、xi、xi+1、xi+2的值，h为网格步长。

步骤4、任意取三组线性权：

①γ1＝0.98，γ2＝0.01，γ3＝0.01；

②γ1＝1/3，γ2＝1/3，γ3＝1/3；

③γ1＝0.01，γ2＝0.495，γ3＝0.495。

步骤5、计算光滑指示器βl，用于衡量重构多项式pl(x)在目标单元上的光滑度，计算公式为：

其中，l＝1，2，3表示对应模板序号，表示多项式pl(x)对x的α阶导数，r1＝4，r2＝1，r3＝1。

步骤6、通过线性权γl和光滑指示器βl计算非线性权ωl，其计算公式为：

其中，l＝1，2，3表示对应模板序号，τ为计算过程中的过渡值，βl为光滑指示器，ε＝10^-6防止分母为零。

步骤7、求出数值通量分裂f⁺(u)在点处的twen0重构值：

其次，将计算结果代入含有时间导数项的半离散有限差分格式，得到关于时间导数的常微分方程。

二、对控制方程中的时间导数使用三阶tvdrunge-kutta离散公式将半离散有限差分格式离散成时空全离散有限差分格式。

利用三阶tvdrunge-kutta离散公式：

得到时空全离散有限差分格式，其中，u⁽¹⁾，u⁽²⁾为中间过渡值，δt为时间步长，上标n表示第n时间层，l(u^m)，l(u⁽¹⁾)，l(u⁽²⁾)为-fx(u)的高阶空间离散形式的近似值。

三、根据时空全离散有限差分格式得到下一时间层上的近似值，依次迭代，得到计算区域内终止时刻流场的数值模拟值。

时空全离散有限差分格式为关于时间层的迭代公式，初始状态值已知，通过迭代公式求出下一时间层的近似值，依次得到终止时刻计算区域内的数值模拟值。对于二维问题，逐维用上面的重构过程。

下面给出几个算例作为本发明所公开方法的具体实施例。

实施例一、台阶问题。该问题是emery于1968年提出的一个用于检验非线性双曲型守恒律格式的经典算例。初始数据为水平来流马赫数为3，密度为1.4，水平速度为3，竖直速度为0，压强为1，管道区域为[0，3]×[0，1]，在距离左边界0.6处有一高度为0.2的台阶，且台阶延伸到管道的尽头。上下边界为反射边界，左边界为来流边界，右边界为出流边界。图1a-1c给出了t＝4时的密度等值线图。

实施例二、双马赫反射问题。该问题描述了一个与x轴成60°角的强激波射到反射墙上发生的变化，来流是马赫数为10的强激波。计算区域为[0，4]×[0，1]。区域底部从y＝0处开始为反射边界条件，其它的底部边界(从x＝0到那部分)为波前条件。图2a-2c给出了t＝0.2时在[0，3]×[0，1]区域的密度等值线图。

实施例三、激波和涡流相互干扰问题。马赫数为1.1的激波位于x＝0.5处且垂直于x轴。激波初始状态是小涡流位于该激波的左边且其中心位于(xc，yc)＝(0.25，0.5)处。涡流可看成平均流的速度，温度和熵的扰动，表示为：

其中，τ＝r/rc，ε＝0.3，rc＝0.05，α＝0.204，γ＝1.4，计算区域为[0，2]×[0，1]。图3a-3c给出了t＝0.35时在[0，1]×[0，1]区域的压强等值线图。图4a-4c给出了t＝0.6时在[0.4，1.45]×[0，1]区域的压强等值线图。图5a-5c给出了t＝0.8时在[0，2]×[0，1]区域的压强等值线图。

实施例四、二维euler黎曼问题。计算区域为[0，1]×[0，1]，初值条件分别设为：

图6a-6c给出了初值条件为(18)时二维euler黎曼问题在t＝0.25时刻的密度等值线图。图7a-7c给出了初值条件为(19)时二维euler黎曼问题在t＝0.25时刻的密度等值线图。图8a-8c给出了初值条件为(20)时二维euler黎曼问题在t＝0.3时刻的密度等值线图。图9a-9c给出了初值条件为(21)时二维euler黎曼问题在t＝0.2时刻的密度等值线图。图10a-10c给出了初值条件为(22)时二维euler黎曼问题在t＝0.3时刻的密度等值线图。从图中可以看出基于三角函数多项式空间的有限差分tweno格式对本发明任意取得的线性权都可以很好的捕捉到黎曼问题的大部分流动特性。

以上仅是本发明的优选实施方式，本发明的保护范围并不仅局限于上述实施例，凡属于本发明思路下的技术方案均属于本发明的保护范围。应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理前提下的若干改进和润饰，应视为本发明的保护范围。