基于交错网格的中心五阶WENO格式的全流场模拟方法与流程

本发明涉及计算流体力学工程技术领域，具体而言涉及一种基于交错网格的中心五阶加权本质无振荡格式的全流场模拟方法。

背景技术：

在计算流体力学数值模拟中，高精度格式的构造及其应用一直是研究的重点内容。因为高精度格式能够对全流场进行精确模拟，能够很好的模拟出解的结构以及准确的捕捉激波位置。1983年，harten首次提出了tvd(totalvariationdiminishing)格式，并在此基础上与osher于1987年提出了eno(essentiallynon-oscillatory)高精度格式，eno格式的主要思想是在逐次扩展的模板中选出最光滑的模板来构造多项式求出单元边界处的值，进而达到高精度、高分辨率和在间断处基本无振荡的效果。但是，在方法的实现过程中，eno格式会造成计算结果的浪费，导致计算效率不高，因此，liu，osher和chan等提出了weno(weightedessentiallynon-oscillatory)格式，提高了计算结果的利用率并且使得r阶精度的eno格式提高到r+1阶精度。1996年，jiang和shu提出了一种新的weno格式，使得数值精度提高到2r-1阶。该经典weno格式的实现过程中，线性权依赖于母模板，且其求解过程相当复杂，因此，zhu和qiu于2016年改善了该weno格式的有限差分形式，在维持精度不减的情况下，随机选取大于零且总和为一的线性权，并于2017年提出了该weno格式的有限体积形式，简称为新形式的有限体积格式，这种新形式的有限体积格式依然在维持精度不减的情况下，随机选取大于零且总和为一的线性权。weno重构最初用于迎风框架，然后将它们应用到中心方案中，进而获得中心方案的优点。例如没有数值通量、稳定性相对较好、不需要应用通量分裂等。结合图2，与迎风框架相比，中心格式也是一类有效的双曲守恒律方案，且相对简单。第一个二阶中心格式是由nessyahu和tadmor于1990年提出的。之后，levy和puppo等人在2000年提出了二维情况下三阶的中心方案下的cweno格式，随后,qiu和shu于2001年提出了一维情况下五阶与九阶的中心方案下的cweno格式，以下简称为cweno-qs格式，之后，levy和puppo等人基于2000年提出了二维情况下三阶的中心方案下的cweno格式。他们又在2002年提出了二维情况下四阶的中心方案下的cweno格式，以下简称为cweno-lp格式。但是，在一维情况下，这些传统的中心方案的cweno格式对于部分算例效果并不理想，而且计算量较大；在二维情况下，这些传统的中心方案的cweno格式对于复杂流场的模拟也并不理想。

技术实现要素：

本发明目的在于提供一种基于在中心框架下的五阶基本加权无振荡格式的全流场模拟方法，在一维情况下能很好地克服上述困难，且在二维情况下能够很好的处理复杂流场的模拟。本发明首先给出一维情况下新型五阶基本加权无振荡格式的全流场模拟方法，与传统方法相比有着更好的精度，并且能够更好的捕捉间断，同时能够模拟部分传统算法不能够解决的问题；然后，给出二维情况下新型五阶基本加权无振荡格式的全流场模拟方法，与传统方法相比有着更好的精度，并且能够模拟传统算法所不能良好解决的双马赫反射问题与台阶问题，以下简称为cweno-zq格式。

为达成上述目的，结合图1，本发明提出一种基于交错网格的中心五阶weno格式的全流场模拟方法，适于在笛卡尔坐标下构造新形式的有限体积五阶weno格式以计算可压缩流场问题，所述全流场模拟方法包括：

s1：对双曲守恒律方程空间变量以及时间变量进行积分得到有限体积格式，采用新形式的有限体积五阶weno格式重构半点处的变量单元均值；

s2：对时间积分项使用gauss求积公式进行离散，得到通量在空间整点关于时间的高斯节点的组合式；

s3：基于单元均值采用新形式的有限体积五阶weno格式重构整点处的变量点值；

s4：基于步骤s3重构的点值采用新形式的有限体积五阶weno格式重构通量的导数在整点处的值；

s5：采用四阶的ncerk方法，基于空间变量的整点值与通量导数的整点值，通过迭代通量的导数在整点处的值进行五阶加权本质无振荡重构的过程，得到所有通量在空间整点关于时间的高斯节点；

组合半点处的单元均值变量与所有通量在空间整点关于时间的高斯节点，得到下一个时间层的半点处的单元均值；

依次迭代，得到计算区域内终止时刻流场的数值结果。

进一步的，步骤s1中，采用新形式的有限体积五阶weno格式重构半点处的变量单元均值的过程包括以下步骤：

s11：基于单元均值，选择母模板以及子模板，重构若干个不同精度的多项式；

s12：针对重构的若干个不同精度的多项式，取满足和为一的任意正数的线性权；

s13：计算光滑指示器，用于衡量重构多项式在目标单元上的光滑程度；

s14：计算非线性权；

s15：组合重构多项式以及非线性权，得到相应的半点处的变量单元均值的重构值。

进一步的，设一维双曲守恒律方程为：

所述全流场模拟方法包括以下步骤：

s1.1：对控制方程分别进行时间方向和空间方向的积分，得到:

其中，是时间层n+1的半点处的单元平均值，是时间层n的半点处的单元平均值，h是空间方向的步长，f(u(xi+1,t),t)是通量在变量于节点处的点值，是从时间层tⁿ到时间层tⁿ⁺¹的积分；

对半点处的单元均值变量进行五阶加权本质无振荡重构；

s1.2：采用下述公式对时间积分项使用gauss求积公式进行离散：

其中，δt是时间步长，αl与τl是高斯系数，得到通量于空间整点关于时间的高斯节点的组合式；

s1.3：基于单元均值使用新形式的有限体积五阶weno格式重构整点处的变量点值；

s1.4：基于步骤s1.3重构的点值使用新形式的有限体积五阶weno格式重构通量的导数在整点处的值；

s1.5：通过四阶的ncerk方法，基于空间变量的整点值与通量导数的整点值，通过迭代通量的导数在整点处的值进行五阶加权本质无振荡重构的过程，得到所有通量于空间整点关于时间的高斯节点，即得到：

其中，bj，cj，bj(τl)所得的值为中间系数，τl为高斯节点系数，δt是时间步长，g^(j)是基于通量的导数在整点处的重构值并通过ncerk方法迭代得到；

组合半点处的单元均值变量与所有通量于空间整点关于时间的高斯节点，得到下一个时间层的半点处的单元均值；

依次迭代，得到计算区域内终止时刻流场的数值结果。

进一步的，步骤s1.1中，所述对半点处的单元均值变量进行五阶加权本质无振荡重构的过程包括以下步骤：

s1.1.1：将目标单元iij以及周围共5个单元组成母模板，记为r1＝{ii-2,ii-1,ii,ii+1,ii+2}，且假设网格步长都为h，并记代表u在网格单元im的平均值(m＝i-2,…,i+2)；

s1.1.2：在母模板中选择2个子模板分别为：r2＝{ii-1,ii}，r3＝{ii,ii+1}；

s1.1.3：在母模板以及子模板上分别求出重构多项式pn(x,y),n＝1,2,3，使其满足:

n＝1,k＝i-2,i-1,i,i+1,i+2；n＝2,k＝i-1,i；n＝3,k＝i,i+1；

则有：

其中，表示u在网格单元ik上的平均值，xi表示网格节点处的值。

以上本发明的技术方案，与现有相比，其显著的有益效果在于：

(1)中心框架使用了交错网格的技术，相对于迎风框架而言不需要求解黎曼问题，进而不需要使用通量分裂。

(2)结构网格下的有限体积五阶新型weno格式可使用任意正线性权，相较于传统有限体积五阶weno格式而言算法结构更加简便，同时拥有更高的计算精度以及更优良的数值模拟效果。

(3)将新型五阶weno格式与中心框架结合，在一维情况下，与传统方法相比有着更好的数值精度，并且能够更好的捕捉间断，同时能够模拟部分低压低密问题，这是传统的方法所不能解决的；然后，给出二维情况下新型五阶基本加权无振荡格式的全流场模拟方法，与传统方法相比有着更好的精度，并且能够更好地数值模拟双马赫反射问题与台阶问题

应当理解，前述构思以及在下面更加详细地描述的额外构思的所有组合只要在这样的构思不相互矛盾的情况下都可以被视为本公开的发明主题的一部分。另外，所要求保护的主题的所有组合都被视为本公开的发明主题的一部分。

结合附图从下面的描述中可以更加全面地理解本发明教导的前述和其他方面、实施例和特征。本发明的其他附加方面例如示例性实施方式的特征和/或有益效果将在下面的描述中显见，或通过根据本发明教导的具体实施方式的实践中得知。

附图说明

附图不意在按比例绘制。在附图中，在各个图中示出的每个相同或近似相同的组成部分可以用相同的标号表示。为了清晰起见，在每个图中，并非每个组成部分均被标记。现在，将通过例子并参考附图来描述本发明的各个方面的实施例，其中：

图1是本发明的基于交错网格的中心五阶加权本质无振荡格式的全流场模拟方法的流程图。

图2是本发明的迎风框架和中心框架的对比图。

图3是本发明的一维母模板示意图。

图4是本发明的二维母模板示意图。

图5是本发明的t＝0.5/π时cweno-zq算法与cweno-qs算法各自在l1与l∞下的精度与误差表现示意图。

图6是本发明的t＝2时cweno-zq算法与cweno-qs算法各自在l1与l∞下的精度与误差表现示意图。

图7是本发明的t＝0.5/π时cweno-zq算法与cweno-lp算法各自在l1与l∞下的精度与误差表现示意图。

图8是本发明的t＝2时cweno-zq算法与cweno-lp算法各自在l1与l∞下的精度与误差表现示意图。

图9是本发明的算例五中的间断捕捉问题的测试结果示意图。

图10是本发明的算例六中的间断捕捉问题的测试结果示意图。

图11是本发明的算例七中的间断捕捉问题的测试结果示意图。

图12是本发明的算例八中的间断捕捉问题的测试结果示意图。

图13是本发明的算例九中的双马赫反射问题的测试结果示意图。

图14是本发明的算例十中的台阶问题的测试结果示意图。

具体实施方式

为了更了解本发明的技术内容，特举具体实施例并配合所附图式说明如下。

具体实施例一

下面结合两个具体的例子，详细阐述本发明提及的新形式的有限体积五阶weno格式、以及本发明所公开的全流场模拟方法分别在解决一维情况和二维情况下的具体过程。

一、一维情况的构造:

设一维euler方程为：

ut+f(u)x＝0(1)

其中，u＝(ρ,ρu,e)^t表示守恒变量,f(u)＝(ρu,ρu²+p,u(e+p))^t，f(u)表示通量，ut表示u对t求导，f(u)x表示f(u)对x求导，ρ、p、u、e分别表示流体密度、压强、水平方向速度和能量，t表示转置。假设网格单元步长是固定的，xi+1/2-xi-1/2＝δx＝h，网格中心为网格单元为ii＝[xi-1/2,xi+1/2]，其中下标i为坐标序号，xi处的值称为整点值，处的值称为半点值。对控制方程分别进行时间方向和空间方向的积分，得到:

其中，是时间层n+1的半点处的单元平均值，是时间层n的半点处的单元平均值，h是空间方向的步长，f(u(xi+1,t),t)是通量在变量于节点处的点值，是从时间层tⁿ到时间层tⁿ⁺¹的积分。

步骤1.1.对半点处的单元均值变量进行新型五阶加权本质无振荡重构,具体过程如下：

步骤1.1.1.如图3所示，将目标单元iij以及周围共5个单元组成母模板，记为r1＝{ii-2,ii-1,ii,ii+1,ii+2}，且假设网格步长都为h，并记代表u在网格单元im的平均值(m＝i-2,…,i+2)。

步骤1.1.2.在母模板中选择2个子模板分别为：r2＝{ii-1,ii}，r3＝{ii,ii+1}。

步骤1.1.3.在母模板以及子模板上分别求出重构多项式pn(x,y),n＝1,2,3，使其满足:

n＝1,k＝i-2,i-1,i,i+1,i+2；n＝2,k＝i-1,i；n＝3,k＝i,i+1；

则有：

其中，表示u在网格单元ik上的平均值，xi表示网格节点处的值。

步骤1.1.4.取满足和为一的任意正数的线性权，避免复杂的数值计算过程，我们取：

γ1＝10.0/12.0,γ2＝1.0/12.0,γ3＝1.0/12.0(7)

作为多项式p1(x)，p2(x)，p3(x)的线性权。

步骤1.1.5.计算光滑指示器，用于衡量重构多项式在目标单元上的光滑度。

计算光滑指示器βl，l＝1,2,3，用于衡量三角函数多项式pl(x)在区间[xi-1/2,xi+1/2]上的光滑度，计算公式为：

其中，h为空间步长，是多项式pj(x)的l次方导数，r是多项式pj(x,y)的最高次数。

步骤1.1.6.计算非线性权，计算公式为：

其中：

其中，ωn是非线性权，与τ是中间值，ε＝10^-10。

步骤1.1.7.求出半点处单元均值的表达式为：

其中，ω1,ω2,ω3是非线性权，γ1,γ2,γ3是线性权。

步骤1.2.对时间积分项利用gauss求积公式进行离散，即：

其中，δt是时间步长，αl与τl是高斯系数，

步骤1.3.对空间变量于整点值进行新型五阶加权本质无振荡重构，具体过程如下：

其得到的重构多项式以及非线性权与步骤1.1得到的插值多项式以及非线性权完全相同。求出半点处单元均值的表达式为：

其中，ω1,ω2,ω3是非线性权，pj(xi)是多项式pj(x)在整点xi处的值，j＝1,2,3，是变量在xi处的重构值，u(xi,tⁿ)是变量在xi处的精确值。

步骤1.4.对通量的导数于整点处进行新型五阶加权本质无振荡重构,其整体重构过程与步骤1.1完全相同，除了构造的插值多项式qn(x),n＝1,2,3，满足：

qn(xk)＝fk(13)

其中：

则有：

其中，fk分别表示f(u)在单元ik处网格中点的值，是多项式qj(x)在整点xi处的导数值，j＝1,2,3，是通量导数在xi处的重构值。

步骤1.5.基于空间变量的整点值与通量导数的整点值，通过迭代通量的导数于整点处的值进行新型五阶加权本质无振荡重构的过程，得到所有通量于空间整点关于时间的高斯节点,即：

其中，

且g⁽²⁾，g⁽³⁾，g⁽⁴⁾是分别基于进行步骤1.4重构得到。

其中：

b1(τl)＝2(1-4b1)τl³+3(3b1-1)τl²+τl,(18)

bj(τl)＝4(3cj-2)bjτl³+3(3-4cj)bjτl²,j＝2,3,4,(19)

步骤1.6.组合半点处的单元均值变量与所有通量于空间整点关于时间的高斯节点，得到下一个时间层的半点处的单元均值，从而依次得到稳定时的全流场数值模拟。

二、二维情况的构造

考虑二维euler方程：

ut+f(u)x+g(u)y＝0(21)

其中，t表示时间变量，x,y表示空间变量，u＝(ρ,ρu,ρv,e)表示守恒变量，f(u)＝(ρu,ρu²+p,ρuv,u(e+p))^t，g(u)＝(ρv,ρuv,ρv²+p,v(e+p))^t，f(u),g(u)表示通量，f(u)x表示f(u)对x求导，g(u)y表示g(u)对y求导，ρ，p，u，v，e分别表示流体密度、压强、水平方向速度、竖直方向速度以及能量，t表示转置。假设网格单元步长是固定的，xi+1/2-xi-1/2＝δx，yj+1/2-yj-1/2＝δy，网格中心为网格单元为iij＝[xi-1/2,xi+1/2]×[yj-1/2,yj+1/2]，其中下标i、j为坐标序号，(xi,yj)处的值称为整点值，处的值称为半点值。对控制方程分别进行时间方向和空间方向的积分，得到：

其中，是时间层n+1的半点处的单元平均值，是时间层n的半点处的单元平均值。

步骤2.1.对半点处的单元均值变量进行新型五阶加权本质无振荡重构,具体过程如下：

步骤2.1.1.如图4所示，将目标单元iij以及周围共25个单元组成母模板，记为r1，且假设网格步长都为h，并记代表u在网格单元im,n的平均值(m＝i-2,…,i+2；n＝j-2,…,j+2)。

步骤2.1.2.在母模板中选择5个子模板分别为：

r2＝{ii-1,j+1,ii,j+1,ii+1,j+1,ii-1,j,ii,j,ii+1,j,ii-1,j-1,ii,j-1,ii+1,j-1}；

r2-1＝{ii-1,j+1,ii,j+1,ii-1,j,ii,j}；

r2-2＝{ii,j+1,ii+1,j+1,ii,j,ii+1,j}；

r2-3＝{ii-1,j,ii,j,ii,j-1,ii-1,j-1}；

r2-4＝{ii,j,ii+1,j,ii,j-1,ii+1,j-1}。

步骤2.1.3.在母模板以及子模板上分别求出插值多项式pn(x,y),n＝5,3,(2_1),(2_2),(2_3),(2_4)，使其满足：

其中：

则有：

其中，表示u在ik,l的网格单元平均值，xi,yj表示网格节点处的值，δx,δy分别表示x方向和y方向的步长。

步骤2.1.4.取满足和为一的任意正数的线性权，避免复杂的数值计算过程，我们取：

其中γ3⁽³⁾,γ2-1⁽³⁾,γ2-2⁽³⁾,γ2-3⁽³⁾,γ2-4⁽³⁾是三阶多项式与二阶多项式结合时所需要的线性权，γ5⁽⁵⁾,γ3⁽⁵⁾,γ2-1⁽⁵⁾,γ2-2⁽⁵⁾,γ2-3⁽⁵⁾,γ2-4⁽⁵⁾是五阶多项式，三阶多项式和二阶多项式结合时所需要的线性权。

步骤2.1.5.计算光滑指示器，用于衡量重构多项式在目标单元上的光滑度。

计算光滑指示器βl，l＝5,3,(2_1),(2_2),(2_3),(2_4)，用于衡量三角函数多项式pl(x,y)在区间[xi-1/2,xi+1/2]×[yj-1/2,yj+1/2]上的光滑度，计算公式为：

其中，α为多重指标，d是求导算子，r是多项式pl(x,y)的最高次数。

步骤2.1.6.计算非线性权：

其中，ωn是非线性权，与τ是中间值，ε＝10^-10。

步骤2.1.7.求出半点处单元均值的表达式为：

步骤2.2.对时间积分项利用gauss求积公式进行离散，对空间积分进行外插处理，即：

其中，δt是时间步长，δx,δy分别是方向x空间步长与y方向空间步长，αl与τl是高斯系数，是外插系数，

步骤2.3.对空间变量于整点处进行新型五阶加权本质无振荡重构,具体过程如下：

其得到的重构多项式以及非线性权与步骤2.1得到的重构多项式以及非线性权完全相同。求出整点处单元点值的表达式为：

其中，ω5,ω3,ω2-1,ω2-2,ω2-3,ω2-4是非线性权，γl⁽⁵⁾，γll⁽³⁾是线性权，l＝5,3,2_1,2_2,2_3,2_4,ll＝3,2_1,2_2,2_3,2_4。

步骤2.4.分别对通量的导数于通量整点值f(u(xi,yj)),g(u(xi,yj))进行新型五阶加权本质无振荡重构,其整体重构过程与步骤2.1完全相同，除了构造的插值多项式qn(x,y),n＝5,3,(2-1),(2-2),(2-3),(2-4)，满足：

其中：

则有整点处通量导数的表达式：

其中，ω5,ω3,ω2-1,ω2-2,ω2-3,ω2-4是非线性权，γl⁽⁵⁾，γll⁽³⁾是线性权，l＝5,3,2_1,2_2,2_3,2_4,ll＝3,2_1,2_2,2_3,2_4。

步骤2.5.基于空间变量的整点值与通量导数的整点值，通过迭代通量的导数于整点处的值进行新型五阶加权本质无振荡重构的过程，得到所有通量于空间整点关于时间的高斯节点：

其中，g⁽²⁾，g⁽³⁾，g⁽⁴⁾是分别基于进行步骤2.4重构得到，高斯系数与步骤1.5相同。

步骤2.6.组合半点处的单元均值变量与所有通量于空间整点关于时间的高斯节点，得到下一个时间层的半点处的单元均值，从而依次得到稳定时的全流场数值模拟。

具体实施例二

下面通过十个算例对本发明所公开的方法做进一步说明。

算例一

关于精度问题的测试。我们考虑一维bugers方程：

初值取u(x,0)＝0.5+sin(πx)，使用周期边界条件，图5分别给出了t＝0.5/π时cweno-zq算法与cweno-qs算法各自在l1与l∞下的精度与误差表现。

算例二

关于精度问题的测试。我们考虑一维euler方程

ut+f(u)x＝0,(43)

其中，u＝(ρ,ρu,e)^t表示守恒变量,f(u)＝(ρu,ρu²+p,u(e+p))^t，f(u)表示通量，ut表示u对t求导，f(u)x表示f(u)对x求导，ρ，p，u，e分别表示流体密度、压强、水平方向速度和能量，t表示转置。初值取ρ(x,0)＝1+0.2sin(πx)，v(x,0)＝1，p(x,0)＝1，使用周期边界条件，图6分别给出了t＝2时cweno-zq算法与cweno-qs算法各自在l1与l∞下的精度与误差表现。

算例三

关于精度问题的测试。我们考虑二维bugers方程

初值取u(x,y,0)＝0.5+sin(π(x+y)/2)，使用周期边界条件，图7分别给出了t＝0.5/π时cweno-zq算法与cweno-lp算法各自在l1与l∞下的精度与误差表现。

算例四

关于精度问题的测试。我们考虑二维euler方程

ut+f(u)x+g(u)y＝0,(45)

其中，t表示时间变量，x,y表示空间变量，u＝(ρ,ρu,ρv,e)^t表示守恒变量，f(u)＝(ρu,ρu²+p,ρuv,u(e+p))^t，g(u)＝(ρv,ρuv,ρv²+p,v(e+p))^t，f(u),g(u)表示通量，f(u)x表示f(u)对x求导，g(u)y表示g(u)对y求导，ρ，p，u，v，e分别表示流体密度、压强、水平方向速度、竖直方向速度以及能量，t表示转置。初值取ρ(x,y,0)＝1+0.2sin(x+y)，u(x,y,0)＝1，v(x,y,0)＝1，p(x,y,0)＝1，使用周期边界条件，图8分别给出了t＝2时cweno-zq算法与cweno-lp算法各自在l1与l∞下的精度与误差表现。

算例五

关于间断捕捉问题的测试。我们一维euler方程(44)的lax问题,其中初值取：

图9从上到下分别是密度图和密度的局部放大图，从左到右分别是cweno-zq方法与cweno-qs。方法图中的实线是精确解，方格是数值解，均取200个网格点。

算例六

关于间断捕捉问题的测试。我们考虑一维euler方程(44)的shockdensitywaveinteraction问题,其中初值取：

图10从上到下分别是密度图和密度的局部放大图，从左到右分别是cweno-zq方法与cweno-qs。方法图中的实线是精确解，方格是数值解，均取400个网格点。

算例七

关于间断捕捉问题的测试。我们考虑一维euler方程(44)的twoblastwave问题,其中初值取：

使用反射边界条件，cweno-qs对于此算例并没有良好的表现，图11从左到右分别给出了t＝0.001时cweno-zq的密度,速度以及压力的分布图，图中的实线是精确解，方格是数值解，均取800个网格点。

算例八

关于间断捕捉问题的测试。我们考虑sedovblastwave问题,其中初值取：

使用外推边界条件，cweno-qs对于此算例并没有良好的表现，图12从左到右分别给出了t＝0.001时cweno-zq的密度、速度以及压力的分布图，图中的实线是精确解，方格是数值解，取800个网格点。

算例九

双马赫反射问题。该问题描述了一个与x轴成60°角的强激波射到反射墙上发生的变化，来流是马赫数为10的强激波。计算区域为[0,4]×[0,1]。区域底部从y＝0处开始为反射边界条件，其它的底部边界(从x＝0到那部分)为波前条件。cweno-lp算法对于此算例并没有良好的表现，图13从上到下分别给出了t＝0.2时cweno-zq算法在[0,3]×[0,1]区域的密度等值线图及其局部放大图，取900×300个网格点。

算例十

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

在本公开中参照附图来描述本发明的各方面，附图中示出了许多说明的实施例。本公开的实施例不必定义在包括本发明的所有方面。应当理解，上面介绍的多种构思和实施例，以及下面更加详细地描述的那些构思和实施方式可以以很多方式中任意一种来实施，这是因为本发明所公开的构思和实施例并不限于任何实施方式。另外，本发明公开的一些方面可以单独使用，或者与本发明公开的其他方面的任何适当组合来使用。

虽然本发明已以较佳实施例揭露如上，然其并非用以限定本发明。本发明所属技术领域中具有通常知识者，在不脱离本发明的精神和范围内，当可作各种的更动与润饰。因此，本发明的保护范围当视权利要求书所界定者为准。