技术领域technical field
本发明涉及计算流体力学工程技术领域，具体而言涉及一种基于交错网格的中心五阶加权本质无振荡格式的全流场模拟方法。The invention relates to the technical field of computational fluid dynamics engineering, in particular to a full flow field simulation method based on a staggered grid-based fifth-order weighted intrinsically non-oscillating format.
背景技术Background technique
在计算流体力学数值模拟中，高精度格式的构造及其应用一直是研究的重点内容。因为高精度格式能够对全流场进行精确模拟，能够很好的模拟出解的结构以及准确的捕捉激波位置。1983年，Harten首次提出了TVD(Total Variation Diminishing)格式，并在此基础上与Osher于1987年提出了ENO(Essentially Non-Oscillatory)高精度格式，ENO格式的主要思想是在逐次扩展的模板中选出最光滑的模板来构造多项式求出单元边界处的值，进而达到高精度、高分辨率和在间断处基本无振荡的效果。但是，在方法的实现过程中，ENO格式会造成计算结果的浪费，导致计算效率不高，因此，Liu，Osher和Chan等提出了WENO(Weighted Essentially Non-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格式对于复杂流场的模拟也并不理想。In the numerical simulation of computational fluid dynamics, the construction and application of high-precision formats have always been the focus of research. Because the high-precision format can accurately simulate the full flow field, the structure of the solution can be well simulated and the shock position can be accurately captured. In 1983, Harten first proposed the TVD (Total Variation Diminishing) format, and on this basis and Osher in 1987 proposed the ENO (Essentially Non-Oscillatory) high-precision format. The main idea of the ENO format is in the successively expanded template. The smoothest template is selected to construct a polynomial to obtain the value at the element boundary, thereby achieving high precision, high resolution and basically no oscillation at discontinuities. However, in the implementation process of the method, the ENO format will cause waste of calculation results, resulting in low calculation efficiency. Therefore, Liu, Osher, and Chan et al. proposed the WENO (Weighted Essentially Non-oscillatory) format to improve the utilization of the calculation results. rate and make the r-order precision ENO format improved to r+1-order precision. In 1996, Jiang and Shu proposed a new WENO format, which improved the numerical accuracy to order 2r-1. In the implementation process of the classic WENO format, the linear weights depend on the parent template, and the solution process is quite complicated. Therefore, Zhu and Qiu improved the finite difference form of the WENO format in 2016. Under the condition that the accuracy is not reduced, Randomly select linear weights greater than zero and sum to one, and proposed a finite volume form of the WENO format in 2017, which is referred to as a new form of finite volume format. This new form of finite volume format is still maintaining accuracy. case, randomly select linear weights greater than zero and sum to one. WENO reconstructions were initially used for upwind frames, and then they were applied to the central scheme to gain the advantages of the central scheme. For example, no numerical flux, relatively good stability, no need to apply flux splitting, etc. Combined with Fig. 2, compared with the upwind frame, the central scheme is also a kind of effective hyperbolic conservation law scheme, and it is relatively simple. The first second-order centroid scheme was proposed by Nessyahu and Tadmor in 1990. After that, in 2000, Levy and Puppo et al. proposed the CWENO scheme under the central scheme of the third order in the two-dimensional case. Subsequently, Qiu and Shu proposed the central scheme of the fifth order and the ninth order in the one-dimensional case in 2001. CWENO format, hereinafter referred to as CWENO-QS format, after that, based on 2000, Levy and Puppo et al. proposed the CWENO format under the center scheme of the third order in the two-dimensional case. In 2002, they proposed the CWENO format under the center scheme of the fourth order in the two-dimensional case, hereinafter referred to as the CWENO-LP format. However, in the one-dimensional case, the CWENO format of these traditional central schemes is not ideal for some cases, and the amount of calculation is large; Simulation isn't ideal either.
发明内容SUMMARY OF THE INVENTION
本发明目的在于提供一种基于在中心框架下的五阶基本加权无振荡格式的全流场模拟方法，在一维情况下能很好地克服上述困难，且在二维情况下能够很好的处理复杂流场的模拟。本发明首先给出一维情况下新型五阶基本加权无振荡格式的全流场模拟方法，与传统方法相比有着更好的精度，并且能够更好的捕捉间断，同时能够模拟部分传统算法不能够解决的问题；然后，给出二维情况下新型五阶基本加权无振荡格式的全流场模拟方法，与传统方法相比有着更好的精度，并且能够模拟传统算法所不能良好解决的双马赫反射问题与台阶问题，以下简称为CWENO-ZQ格式。The purpose of the present invention is to provide a full flow field simulation method based on the fifth-order basic weighted non-oscillation scheme under the central frame, which can well overcome the above difficulties in the one-dimensional case, and can be very good in the two-dimensional case. Handle simulations of complex flow fields. The invention first provides a new fifth-order basic weighted non-oscillation full flow field simulation method in one-dimensional case, which has better accuracy than traditional methods, and can better capture discontinuities. The problem that can be solved; then, a new fifth-order basic weighted non-oscillation full flow field simulation method in two-dimensional case is given, which has better accuracy than the traditional method, and can simulate the two-dimensional problem that cannot be solved well by the traditional algorithm. Mach reflection problem and step problem, hereinafter referred to as CWENO-ZQ format.
为达成上述目的，结合图1，本发明提出一种基于交错网格的中心五阶WENO格式的全流场模拟方法，适于在笛卡尔坐标下构造新形式的有限体积五阶WENO格式以计算可压缩流场问题，所述全流场模拟方法包括：In order to achieve the above purpose, with reference to Fig. 1, the present invention proposes a full flow field simulation method based on a staggered grid center fifth-order WENO format, which is suitable for constructing a new form of finite volume fifth-order WENO format under Cartesian coordinates to calculate For the compressible flow field problem, the full flow field simulation method includes:
S1：对双曲守恒律方程空间变量以及时间变量进行积分得到有限体积格式，采用新形式的有限体积五阶WENO格式重构半点处的变量单元均值；S1: Integrate the space variables and time variables of the hyperbolic conservation law equation to obtain the finite volume format, and use the new form of the fifth-order WENO format of the finite volume to reconstruct the variable element mean at the half point;
S2：对时间积分项使用Gauss求积公式进行离散，得到通量在空间整点关于时间的高斯节点的组合式；S2: Use the Gauss quadrature formula to discretize the time integral term to obtain the combination of the Gaussian nodes of the flux at the whole point of space with respect to time;
S3：基于单元均值采用新形式的有限体积五阶WENO格式重构整点处的变量点值；S3: Reconstruct the variable point value at the whole point based on the element mean using a new form of the finite volume fifth-order WENO format;
S4：基于步骤S3重构的点值采用新形式的有限体积五阶WENO格式重构通量的导数在整点处的值；S4: based on the point value reconstructed in step S3, the new form of the finite volume fifth-order WENO format is used to reconstruct the value of the derivative of the flux at the whole point;
S5：采用四阶的NCERK方法，基于空间变量的整点值与通量导数的整点值，通过迭代通量的导数在整点处的值进行五阶加权本质无振荡重构的过程，得到所有通量在空间整点关于时间的高斯节点；S5: Using the fourth-order NCERK method, based on the integral point value of the spatial variable and the integral point value of the flux derivative, the fifth-order weighted intrinsically non-oscillating reconstruction process is performed by iterating the value of the derivative of the flux at the integral point to obtain All fluxes are Gaussian nodes with respect to time at the whole point in space;
组合半点处的单元均值变量与所有通量在空间整点关于时间的高斯节点，得到下一个时间层的半点处的单元均值；Combine the element mean variable at the half point with the Gaussian nodes of all fluxes at the whole point in space with respect to time to obtain the element mean at the half point of the next time layer;
依次迭代，得到计算区域内终止时刻流场的数值结果。Iterate sequentially to obtain the numerical results of the flow field at the termination time in the calculation area.
进一步的，步骤S1中，采用新形式的有限体积五阶WENO格式重构半点处的变量单元均值的过程包括以下步骤：Further, in step S1, the process of reconstructing the mean value of the variable element at the half point by adopting the new finite volume fifth-order WENO format includes the following steps:
S11：基于单元均值，选择母模板以及子模板，重构若干个不同精度的多项式；S11: Based on the unit mean, select the mother template and the sub-template, and reconstruct several polynomials with different precisions;
S12：针对重构的若干个不同精度的多项式，取满足和为一的任意正数的线性权；S12: For several reconstructed polynomials with different precisions, take the linear weight of any positive number satisfying the sum of one;
S13：计算光滑指示器，用于衡量重构多项式在目标单元上的光滑程度；S13: Calculate the smoothness indicator, which is used to measure the smoothness of the reconstructed polynomial on the target unit;
S14：计算非线性权；S14: Calculate the nonlinear weight;
S15：组合重构多项式以及非线性权，得到相应的半点处的变量单元均值的重构值。S15: Combine the reconstruction polynomial and the nonlinear weight to obtain the reconstruction value of the variable unit mean value at the corresponding half point.
进一步的，设一维双曲守恒律方程为：Further, let the one-dimensional hyperbolic conservation law equation be:
所述全流场模拟方法包括以下步骤：The full flow field simulation method includes the following steps:
S1.1：对控制方程分别进行时间方向和空间方向的积分，得到:S1.1: Integrate the control equations in the time direction and the space direction respectively, and get:
其中，是时间层n+1的半点处的单元平均值，是时间层n的半点处的单元平均值，h是空间方向的步长，f(u(xi+1,t),t)是通量在变量于节点处的点值，是从时间层tn到时间层tn+1的积分；in, is the cell mean at the half-point of time layer n+1, is the cell mean value at the half point of the time layer n, h is the step size in the spatial direction, f(u(xi+1 ,t),t) is the point value of the flux at the node of the variable, is the integral from time horizon tn to time horizon tn+1 ;
对半点处的单元均值变量进行五阶加权本质无振荡重构；unit mean variable at half point Perform fifth-order weighted intrinsically oscillation-free reconstruction;
S1.2：采用下述公式对时间积分项使用Gauss求积公式进行离散：S1.2: Use the following formula to discretize the time integral term using the Gauss quadrature formula:
其中，Δt是时间步长，αl与τl是高斯系数，得到通量于空间整点关于时间的高斯节点的组合式；Among them, Δt is the time step, αl and τl are the Gaussian coefficients, and the combined formula of the Gaussian node of the flux at the whole point of space with respect to time is obtained;
S1.3：基于单元均值使用新形式的有限体积五阶WENO格式重构整点处的变量点值；S1.3: Reconstruct the variable point value at the whole point using a new form of finite volume fifth-order WENO scheme based on the element mean;
S1.4：基于步骤S1.3重构的点值使用新形式的有限体积五阶WENO格式重构通量的导数在整点处的值；S1.4: Reconstruct the value of the derivative of the flux at the whole point based on the point value reconstructed in step S1.3 using the new form of the finite volume fifth-order WENO scheme;
S1.5：通过四阶的NCERK方法，基于空间变量的整点值与通量导数的整点值，通过迭代通量的导数在整点处的值进行五阶加权本质无振荡重构的过程，得到所有通量于空间整点关于时间的高斯节点，即得到：S1.5: Through the fourth-order NCERK method, based on the integer value of the spatial variable and the integer value of the flux derivative, the fifth-order weighted intrinsically non-oscillating reconstruction process is performed by iterating the value of the derivative of the flux at the integer point. , get all the Gaussian nodes of the flux at the whole point of space with respect to time, that is, get:
其中，bj，cj，bj(τl)所得的值为中间系数，τl为高斯节点系数，Δt是时间步长，g(j)是基于通量的导数在整点处的重构值并通过NCERK方法迭代得到；where bj , cj , bj (τl ) are the intermediate coefficients, τl is the Gaussian node coefficient, Δt is the time step, and g(j) is the weight of the flux-based derivative at the whole point The value is constructed and obtained iteratively through the NCERK method;
组合半点处的单元均值变量与所有通量于空间整点关于时间的高斯节点，得到下一个时间层的半点处的单元均值；Combining the element mean variable at the half point with all the Gaussian nodes of the flux at the whole point in space with respect to time, the element mean at the half point of the next time layer is obtained;
依次迭代，得到计算区域内终止时刻流场的数值结果。Iterate sequentially to obtain the numerical results of the flow field at the termination time in the calculation area.
进一步的，步骤S1.1中，所述对半点处的单元均值变量进行五阶加权本质无振荡重构的过程包括以下步骤：Further, in step S1.1, the unit mean variable at the half point The process of performing a fifth-order weighted intrinsically oscillation-free reconstruction consists of the following steps:
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.1: The target unit Iij and 5 surrounding units form a master template, denoted as R1 ={Ii-2 ,Ii-1 ,Ii ,Ii+1 ,Ii+2 }, and Assuming that the grid steps are all h, and record Represents the average value of U in the grid unit Im (m =i-2,...,i+2);
S1.1.2：在母模板中选择2个子模板分别为：R2＝{Ii-1,Ii}，R3＝{Ii,Ii+1}；S1.1.2: Select two sub-templates in the master template: R2 ={Ii-1 ,Ii }, R3 ={Ii ,Ii+1 };
S1.1.3：在母模板以及子模板上分别求出重构多项式pn(x,y),n＝1,2,3，使其满足:S1.1.3: Find the reconstruction polynomial pn (x, y), n=1, 2, 3 on the mother template and the child template respectively, so that it satisfies:
n＝1,k＝i-2,i-1,i,i+1,i+2；n＝2,k＝i-1,i；n＝3,k＝i,i+1；n=1, k=i-2, i-1, i, i+1, i+2; n=2, k=i-1, i; n=3, k=i, i+1;
则有：Then there are:
其中，表示U在网格单元Ik上的平均值，xi表示网格节点处的值。in, represents the average value of U over grid element Ik , andxi represents the value at grid nodes.
以上本发明的技术方案，与现有相比，其显著的有益效果在于：The above technical scheme of the present invention, compared with the existing ones, has the following significant beneficial effects:
(1)中心框架使用了交错网格的技术，相对于迎风框架而言不需要求解黎曼问题，进而不需要使用通量分裂。(1) The central frame uses a staggered grid technique, which does not require solving the Riemann problem compared to the upwind frame, and thus does not need to use flux splitting.
(2)结构网格下的有限体积五阶新型WENO格式可使用任意正线性权，相较于传统有限体积五阶WENO格式而言算法结构更加简便，同时拥有更高的计算精度以及更优良的数值模拟效果。(2) The new fifth-order WENO scheme of finite volume under structured grid can use any positive linear weight. Compared with the traditional fifth-order WENO scheme of finite volume, the algorithm structure is simpler, and it has higher calculation accuracy and better numerical value. Simulation effect.
(3)将新型五阶WENO格式与中心框架结合，在一维情况下，与传统方法相比有着更好的数值精度，并且能够更好的捕捉间断，同时能够模拟部分低压低密问题，这是传统的方法所不能解决的；然后，给出二维情况下新型五阶基本加权无振荡格式的全流场模拟方法，与传统方法相比有着更好的精度，并且能够更好地数值模拟双马赫反射问题与台阶问题(3) Combining the new fifth-order WENO scheme with the central frame, in the one-dimensional case, it has better numerical accuracy than the traditional method, and can better capture discontinuities, and can simulate some low-voltage and low-density problems. can not be solved by traditional methods; then, a new fifth-order basic weighted non-oscillation full flow field simulation method in two-dimensional case is given, which has better accuracy and better numerical simulation than traditional methods. Double Mach reflection problem and step problem
应当理解，前述构思以及在下面更加详细地描述的额外构思的所有组合只要在这样的构思不相互矛盾的情况下都可以被视为本公开的发明主题的一部分。另外，所要求保护的主题的所有组合都被视为本公开的发明主题的一部分。It is to be understood that all combinations of the foregoing concepts, as well as additional concepts described in greater detail below, are considered to be part of the inventive subject matter of the present disclosure to the extent that such concepts are not contradictory. Additionally, all combinations of the claimed subject matter are considered to be part of the inventive subject matter of this disclosure.
结合附图从下面的描述中可以更加全面地理解本发明教导的前述和其他方面、实施例和特征。本发明的其他附加方面例如示例性实施方式的特征和/或有益效果将在下面的描述中显见，或通过根据本发明教导的具体实施方式的实践中得知。The foregoing and other aspects, embodiments and features of the present teachings can be more fully understood from the following description when taken in conjunction with the accompanying drawings. Other additional aspects of the invention, such as features and/or benefits of the exemplary embodiments, will be apparent from the description below, or learned by practice of specific embodiments in accordance with the teachings of this invention.
附图说明Description of drawings
附图不意在按比例绘制。在附图中，在各个图中示出的每个相同或近似相同的组成部分可以用相同的标号表示。为了清晰起见，在每个图中，并非每个组成部分均被标记。现在，将通过例子并参考附图来描述本发明的各个方面的实施例，其中：The drawings are not intended to be drawn to scale. In the drawings, each identical or nearly identical component that is illustrated in various figures may be represented by the same reference numeral. For clarity, not every component is labeled in every figure. Embodiments of various aspects of the present invention will now be described by way of example and with reference to the accompanying drawings, wherein:
图1是本发明的基于交错网格的中心五阶加权本质无振荡格式的全流场模拟方法的流程图。FIG. 1 is a flow chart of a method for simulating a full flow field based on a staggered grid-based fifth-order weighted essentially non-oscillating format.
图2是本发明的迎风框架和中心框架的对比图。Figure 2 is a comparison diagram of the windward frame and the center frame of the present invention.
图3是本发明的一维母模板示意图。3 is a schematic diagram of a one-dimensional mother template of the present invention.
图4是本发明的二维母模板示意图。4 is a schematic diagram of a two-dimensional master template of the present invention.
图5是本发明的T＝0.5/π时CWENO-ZQ算法与CWENO-QS算法各自在L1与L∞下的精度与误差表现示意图。5 is a schematic diagram showing the performance of the precision and error of the CWENO-ZQ algorithm and the CWENO-QS algorithm under L1 and L∞ respectively when T=0.5/π according to the present invention.
图6是本发明的t＝2时CWENO-ZQ算法与CWENO-QS算法各自在L1与L∞下的精度与误差表现示意图。6 is a schematic diagram showing the performance of the precision and error of the CWENO-ZQ algorithm and the CWENO-QS algorithm under L1 and L∞ respectively when t=2 of the present invention.
图7是本发明的T＝0.5/π时CWENO-ZQ算法与CWENO-LP算法各自在L1与L∞下的精度与误差表现示意图。7 is a schematic diagram showing the performance of the precision and error of the CWENO-ZQ algorithm and the CWENO-LP algorithm under L1 and L∞ respectively when T=0.5/π of the present invention.
图8是本发明的t＝2时CWENO-ZQ算法与CWENO-LP算法各自在L1与L∞下的精度与误差表现示意图。FIG. 8 is a schematic diagram showing the performance of the accuracy and error of the CWENO-ZQ algorithm and the CWENO-LP algorithm under L1 and L∞ respectively when t=2 of the present invention.
图9是本发明的算例五中的间断捕捉问题的测试结果示意图。FIG. 9 is a schematic diagram of the test result of the discontinuous capture problem in the fifth calculation example of the present invention.
图10是本发明的算例六中的间断捕捉问题的测试结果示意图。FIG. 10 is a schematic diagram of the test result of the discontinuous capture problem in the sixth calculation example of the present invention.
图11是本发明的算例七中的间断捕捉问题的测试结果示意图。FIG. 11 is a schematic diagram of the test result of the discontinuous capture problem in the seventh calculation example of the present invention.
图12是本发明的算例八中的间断捕捉问题的测试结果示意图。FIG. 12 is a schematic diagram of the test result of the discontinuous capture problem in the eighth calculation example of the present invention.
图13是本发明的算例九中的双马赫反射问题的测试结果示意图。FIG. 13 is a schematic diagram of the test result of the double Mach reflection problem in the ninth calculation example of the present invention.
图14是本发明的算例十中的台阶问题的测试结果示意图。FIG. 14 is a schematic diagram of the test result of the step problem in the tenth calculation example of the present invention.
具体实施方式Detailed ways
为了更了解本发明的技术内容，特举具体实施例并配合所附图式说明如下。In order to better understand the technical content of the present invention, specific embodiments are given and described below in conjunction with the accompanying drawings.
具体实施例一Specific embodiment one
下面结合两个具体的例子，详细阐述本发明提及的新形式的有限体积五阶WENO格式、以及本发明所公开的全流场模拟方法分别在解决一维情况和二维情况下的具体过程。The following is a detailed description of the new form of the finite volume fifth-order WENO scheme mentioned in the present invention and the full flow field simulation method disclosed in the present invention in conjunction with two specific examples. .
一、一维情况的构造:First, the construction of the one-dimensional case:
设一维Euler方程为：Let the one-dimensional Euler equation be:
Ut+f(U)x＝0 (1)Ut +f(U)x =0 (1)
其中，U＝(ρ,ρu,E)T表示守恒变量,f(U)＝(ρu,ρu2+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处的值称为整点值，处的值称为半点值。对控制方程分别进行时间方向和空间方向的积分，得到:Among them, U=(ρ, ρu, E)T represents the conserved variable, f(U)=(ρu, ρu2 +p, u(E+p))T , f(U) represents the flux, Ut represents U Deriving from t, f(U)x represents the derivation of f(U) against x, ρ, p, u, and E represent fluid density, pressure, velocity and energy in the horizontal direction, respectively, and T represents transposition. Assuming that the grid cell step size is fixed, xi+1/2 -xi-1/2 =Δx=h, the grid center is The grid unit is Ii =[xi-1/2 ,xi+1/2 ], where the subscript i is the coordinate number, and the value at xi is called the integer value, The value at is called the half point value. Integrate the governing equations in the time direction and the space direction, respectively, to get:
其中，是时间层n+1的半点处的单元平均值，是时间层n的半点处的单元平均值，h是空间方向的步长，f(u(xi+1,t),t)是通量在变量于节点处的点值，是从时间层tn到时间层tn+1的积分。in, is the cell mean at the half-point of time layer n+1, is the cell mean value at the half point of the time layer n, h is the step size in the spatial direction, f(u(xi+1 ,t),t) is the point value of the flux at the node of the variable, is the integral from time horizontn to time horizon tn+1 .
步骤1.1.对半点处的单元均值变量进行新型五阶加权本质无振荡重构,具体过程如下：Step 1.1. Unit mean variable at half point The new fifth-order weighted intrinsically oscillation-free reconstruction is carried out, and the specific process is as follows:
步骤1.1.1.如图3所示，将目标单元Iij以及周围共5个单元组成母模板，记为R1＝{Ii-2,Ii-1,Ii,Ii+1,Ii+2}，且假设网格步长都为h，并记代表U在网格单元Im的平均值(m＝i-2,…,i+2)。Step 1.1.1. As shown in Figure 3, the target unit Iij and 5 surrounding units form a master template, denoted as R1 ={Ii-2 ,Ii-1 ,Ii ,Ii+1 , Ii+2 }, and assuming that the grid steps are all h, and record Represents the average value of U in grid cell Im (m =i-2, . . . , i+2).
步骤1.1.2.在母模板中选择2个子模板分别为：R2＝{Ii-1,Ii}，R3＝{Ii,Ii+1}。Step 1.1.2. Select two sub-templates in the master template, respectively: R2 ={Ii-1 ,Ii }, R3 ={Ii ,Ii+1 }.
步骤1.1.3.在母模板以及子模板上分别求出重构多项式pn(x,y),n＝1,2,3，使其满足:Step 1.1.3. Find the reconstruction polynomial pn (x, y), n=1, 2, 3 on the mother template and the child template respectively, so that it satisfies:
n＝1,k＝i-2,i-1,i,i+1,i+2；n＝2,k＝i-1,i；n＝3,k＝i,i+1；n=1, k=i-2, i-1, i, i+1, i+2; n=2, k=i-1, i; n=3, k=i, i+1;
则有：Then there are:
其中，表示U在网格单元Ik上的平均值，xi表示网格节点处的值。in, represents the average value of U over grid element Ik , andxi represents the value at grid nodes.
步骤1.1.4.取满足和为一的任意正数的线性权，避免复杂的数值计算过程，我们取：Step 1.1.4. Take the linear weight of any positive number satisfying the sum of one to avoid complicated numerical calculation process, we take:
γ1＝10.0/12.0,γ2＝1.0/12.0,γ3＝1.0/12.0 (7)γ1 =10.0/12.0, γ2 =1.0/12.0, γ3 =1.0/12.0 (7)
作为多项式p1(x)，p2(x)，p3(x)的线性权。As the linear weights of the polynomials p1 (x), p2 (x), p3 (x).
步骤1.1.5.计算光滑指示器，用于衡量重构多项式在目标单元上的光滑度。Step 1.1.5. Calculate the smoothness indicator, which measures the smoothness of the reconstructed polynomial on the target element.
计算光滑指示器βl，l＝1,2,3，用于衡量三角函数多项式pl(x)在区间[xi-1/2,xi+1/2]上的光滑度，计算公式为：Calculate the smoothness indicator βl , l=1,2,3, which is used to measure the smoothness of the trigonometric function polynomial pl (x) on the interval [xi-1/2 ,xi+1/2 ], the calculation formula for:
其中，h为空间步长，是多项式pj(x)的l次方导数，r是多项式pj(x,y)的最高次数。where h is the space step, is the l-th derivative of the polynomial pj (x), and r is the highest degree of the polynomial pj (x,y).
步骤1.1.6.计算非线性权，计算公式为：Step 1.1.6. Calculate the nonlinear weight, the calculation formula is:
其中： in:
其中，ωn是非线性权，与τ是中间值，ε＝10-10。where ωn is the nonlinear weight, With τ being an intermediate value, ε=10−10 .
步骤1.1.7.求出半点处单元均值的表达式为：Step 1.1.7. The expression for finding the unit mean at half point is:
其中，ω1,ω2,ω3是非线性权，γ1,γ2,γ3是线性权。Among them, ω1 , ω2 , and ω3 are nonlinear weights, and γ1 , γ2 , and γ3 are linear weights.
步骤1.2.对时间积分项利用Gauss求积公式进行离散，即：Step 1.2. Use the Gauss quadrature formula to discretize the time integral term, namely:
其中，Δt是时间步长，αl与τl是高斯系数，where Δt is the time step, αl and τl are the Gaussian coefficients,
步骤1.3.对空间变量于整点值进行新型五阶加权本质无振荡重构，具体过程如下：Step 1.3. Perform a new fifth-order weighted non-oscillation reconstruction of the space variable at the whole point value. The specific process is as follows:
其得到的重构多项式以及非线性权与步骤1.1得到的插值多项式以及非线性权完全相同。求出半点处单元均值的表达式为：The obtained reconstruction polynomial and nonlinear weight are exactly the same as the interpolation polynomial and nonlinear weight obtained in step 1.1. The expression to find the unit mean at half point is:
其中，ω1,ω2,ω3是非线性权，pj(xi)是多项式pj(x)在整点xi处的值，j＝1,2,3，是变量在xi处的重构值，u(xi,tn)是变量在xi处的精确值。Among them, ω1 , ω2 , ω3 are nonlinear weights, pj (xi ) is the value of the polynomial pj (x) at the whole point xi , j=1, 2, 3, is the reconstructed value of the variable at xi and u(xi ,tn ) is the exact value of the variable at xi .
步骤1.4.对通量的导数于整点处进行新型五阶加权本质无振荡重构,其整体重构过程与步骤1.1完全相同，除了构造的插值多项式qn(x),n＝1,2,3，满足：Step 1.4. Perform a new fifth-order weighted intrinsically non-oscillating reconstruction of the derivative of the flux at the whole point. The overall reconstruction process is exactly the same as that of Step 1.1, except for the constructed interpolation polynomial qn (x), n = 1, 2 , 3, satisfy:
qn(xk)＝fk (13)qn (xk )=fk (13)
其中：in:
则有：Then there are:
其中，fk分别表示f(U)在单元Ik处网格中点的值，是多项式qj(x)在整点xi处的导数值，j＝1,2,3，是通量导数在xi处的重构值。Among them, fk represents the value of f(U) at the grid midpoint of unit Ik , respectively, is the derivative value of the polynomial qj (x) at the integer point xi , j=1,2,3, is the reconstructed value of the flux derivative at xi .
步骤1.5.基于空间变量的整点值与通量导数的整点值，通过迭代通量的导数于整点处的值进行新型五阶加权本质无振荡重构的过程，得到所有通量于空间整点关于时间的高斯节点,即：Step 1.5. Based on the integral point value of the space variable and the integral point value of the flux derivative, perform a new fifth-order weighted intrinsically non-oscillating reconstruction process by iterating the value of the derivative of the flux at the integral point to obtain all the fluxes in space. The whole point is a Gaussian node about time, namely:
其中，in,
且g(2)，g(3)，g(4)是分别基于进行步骤1.4重构得到。And g(2) , g(3) , g(4) are based on Perform step 1.4 to reconstruct.
其中：in:
b1(τl)＝2(1-4b1)τl3+3(3b1-1)τl2+τl, (18)b1 (τl )=2(1-4b1 )τl3 +3(3b1 -1)τl2 +τl , (18)
bj(τl)＝4(3cj-2)bjτl3+3(3-4cj)bjτl2,j＝2,3,4, (19)bj (τl )=4(3cj -2)bj τl3 +3(3-4cj )bj τl2 ,j=2,3,4,(19)
步骤1.6.组合半点处的单元均值变量与所有通量于空间整点关于时间的高斯节点，得到下一个时间层的半点处的单元均值，从而依次得到稳定时的全流场数值模拟。Step 1.6. Combine the element mean variable at the half point and all the Gauss nodes of the flux at the whole point of space with respect to time to obtain the element mean at the half point of the next time layer, so as to obtain the numerical simulation of the full flow field when it is stable in turn.
二、二维情况的构造2. The structure of the two-dimensional case
考虑二维Euler方程：Consider the two-dimensional Euler equation:
Ut+f(U)x+g(U)y＝0 (21)Ut +f(U)x +g(U)y =0 (21)
其中，t表示时间变量，x,y表示空间变量，U＝(ρ,ρu,ρv,E)表示守恒变量，f(U)＝(ρu,ρu2+p,ρuv,u(E+p))T，g(U)＝(ρv,ρuv,ρv2+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)处的值称为整点值，处的值称为半点值。对控制方程分别进行时间方向和空间方向的积分，得到：Among them, t represents time variables, x, y represent spatial variables, U=(ρ, ρu, ρv, E) represents conserved variables, f(U)=(ρu, ρu2 +p, ρuv, u(E+p) )T , g(U)=(ρv,ρuv,ρv2 +p,v(E+p))T , f(U),g(U) represents flux, f(U)x represents f(U) Derivative to x, g(U)y represents the derivative of g(U) to y, ρ, p, u, v, E represent fluid density, pressure, horizontal velocity, vertical velocity and energy, respectively, T represents rotation set. Assuming that the grid cell step size is fixed, xi+1/2 -xi-1/2 =Δx, yj+1/2 -yj-1/2 =Δy, and the grid center is The grid unit is Iij =[xi-1/2 ,xi+1/2 ]×[yj-1/2 ,yj+1/2 ], where the subscripts i and j are the coordinate numbers, ( The value at xi , yj ) is called the integer value, The value at is called the half point value. Integrate the governing equations in the time direction and the space direction, respectively, to get:
其中，是时间层n+1的半点处的单元平均值，是时间层n的半点处的单元平均值。in, is the cell mean at the half-point of time layer n+1, is the cell mean at half-point in time layer n.
步骤2.1.对半点处的单元均值变量进行新型五阶加权本质无振荡重构,具体过程如下：Step 2.1. Unit mean variable at half point The new fifth-order weighted intrinsically oscillation-free reconstruction is carried out, and the specific process is as follows:
步骤2.1.1.如图4所示，将目标单元Iij以及周围共25个单元组成母模板，记为R1，且假设网格步长都为h，并记代表U在网格单元Im,n的平均值(m＝i-2,…,i+2；n＝j-2,…,j+2)。Step 2.1.1. As shown in Figure 4, the target unit Iij and a total of 25 surrounding units form a master template, denoted as R1 , and the grid step size is assumed to be h, and denoted as Represents the average value of U in grid cellsIm,n (m=i-2,...,i+2; n=j-2,...,j+2).
步骤2.1.2.在母模板中选择5个子模板分别为：Step 2.1.2. Select 5 sub-templates in the parent template:
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 ={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-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-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-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}。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)，使其满足：Step 2.1.3. Calculate the interpolation polynomial pn (x, y), n=5, 3, (2_1), (2_2), (2_3), (2_4) on the mother template and the child template respectively, so that it satisfies :
其中：in:
则有：Then there are:
其中，表示U在Ik,l的网格单元平均值，xi,yj表示网格节点处的值，Δx,Δy分别表示x方向和y方向的步长。in, Represents the grid cell average value of U in Ik,l , xi , yj represent the values at the grid nodes, Δx, Δy represent the step size in the x and y directions, respectively.
步骤2.1.4.取满足和为一的任意正数的线性权，避免复杂的数值计算过程，我们取：Step 2.1.4. Take the linear weight of any positive number satisfying the sum of one to avoid complicated numerical calculation process, we take:
其中γ3(3),γ2-1(3),γ2-2(3),γ2-3(3),γ2-4(3)是三阶多项式与二阶多项式结合时所需要的线性权，γ5(5),γ3(5),γ2-1(5),γ2-2(5),γ2-3(5),γ2-4(5)是五阶多项式，三阶多项式和二阶多项式结合时所需要的线性权。Among them γ3(3) , γ2-1(3) , γ2-2(3) , γ2-3(3) , γ2-4(3) are the third-order polynomial and second-order polynomial required when combining The linear weights of γ5(5) , γ3(5) , γ2-1(5) , γ2-2(5) , γ2-3(5) , γ2-4(5) are the fifth order Polynomial, linear weights required when combining third-order polynomials and second-order polynomials.
步骤2.1.5.计算光滑指示器，用于衡量重构多项式在目标单元上的光滑度。Step 2.1.5. Calculate the smoothness indicator, which measures the smoothness of the reconstructed polynomial on the target element.
计算光滑指示器β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]上的光滑度，计算公式为：Calculate the smoothness indicator βl , l=5,3,(2_1),(2_2),(2_3),(2_4), to measure the trigonometric function polynomial pl (x,y) in the interval [xi-1/ 2 ,xi+1/2 ]×[yj-1/2 ,yj+1/2 ]The smoothness is calculated as:
其中，α为多重指标，D是求导算子，r是多项式pl(x,y)的最高次数。Among them, α is the multiple index, D is the derivation operator, and r is the highest degree of the polynomial pl (x, y).
步骤2.1.6.计算非线性权：Step 2.1.6. Calculate nonlinear weights:
其中，ωn是非线性权，与τ是中间值，ε＝10-10。where ωn is the nonlinear weight, With τ being an intermediate value, ε=10−10 .
步骤2.1.7.求出半点处单元均值的表达式为：Step 2.1.7. The expression for finding the unit mean at half point is:
步骤2.2.对时间积分项利用Gauss求积公式进行离散，对空间积分进行外插处理，即：Step 2.2. Use the Gauss quadrature formula to discretize the time integral term, and extrapolate the spatial integral, namely:
其中，Δt是时间步长，Δx,Δy分别是方向x空间步长与y方向空间步长，αl与τl是高斯系数，是外插系数，Among them, Δt is the time step, Δx, Δy are the x space step and the y direction space step, respectively, αl and τl are the Gaussian coefficients, is the extrapolation coefficient,
步骤2.3.对空间变量于整点处进行新型五阶加权本质无振荡重构,具体过程如下：Step 2.3. Perform a new fifth-order weighted essential non-oscillation reconstruction of the spatial variable at the whole point. The specific process is as follows:
其得到的重构多项式以及非线性权与步骤2.1得到的重构多项式以及非线性权完全相同。求出整点处单元点值的表达式为：The obtained reconstruction polynomial and nonlinear weight are exactly the same as the reconstruction polynomial and nonlinear weight obtained in step 2.1. The expression to find the value of the unit point at the whole point is:
其中，ω5,ω3,ω2-1,ω2-2,ω2-3,ω2-4是非线性权，γl(5)，γll(3)是线性权，l＝5,3,2_1,2_2,2_3,2_4,ll＝3,2_1,2_2,2_3,2_4。Among them, ω5 , ω3 , ω2-1 , ω2-2 , ω2-3 , ω2-4 are nonlinear weights, γl(5) ,γ 1(3) are linear weights, 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)，满足：Step 2.4. Perform a novel fifth-order weighted intrinsically non-oscillating reconstruction of the flux derivatives at the integral point values f(u(xi , yj )), g(u(xi , yj )) respectively, which The overall reconstruction process is exactly the same as step 2.1, except that the constructed interpolation polynomial qn (x, y), n=5, 3, (2-1), (2-2), (2-3), (2- 4), satisfy:
其中：in:
则有整点处通量导数的表达式：Then there is an expression for the flux derivative at the whole point:
其中，ω5,ω3,ω2-1,ω2-2,ω2-3,ω2-4是非线性权，γl(5)，γll(3)是线性权，l＝5,3,2_1,2_2,2_3,2_4,ll＝3,2_1,2_2,2_3,2_4。Among them, ω5 , ω3 , ω2-1 , ω2-2 , ω2-3 , ω2-4 are nonlinear weights, γl(5) ,γ 1(3) are linear weights, l=5, 3,2_1,2_2,2_3,2_4,ll=3,2_1,2_2,2_3,2_4.
步骤2.5.基于空间变量的整点值与通量导数的整点值，通过迭代通量的导数于整点处的值进行新型五阶加权本质无振荡重构的过程，得到所有通量于空间整点关于时间的高斯节点：Step 2.5. Based on the integral point value of the space variable and the integral point value of the flux derivative, perform a new fifth-order weighted intrinsically non-oscillating reconstruction process by iterating the value of the derivative of the flux at the integral point to obtain all the fluxes in the space. The whole point of the Gaussian node about time:
其中，g(2)，g(3)，g(4)是分别基于进行步骤2.4重构得到，高斯系数与步骤1.5相同。in, g(2) , g(3) , g(4) are based on After reconstruction in step 2.4, the Gaussian coefficient is the same as that in step 1.5.
步骤2.6.组合半点处的单元均值变量与所有通量于空间整点关于时间的高斯节点，得到下一个时间层的半点处的单元均值，从而依次得到稳定时的全流场数值模拟。Step 2.6. Combine the element mean variable at the half point with the Gaussian nodes of all fluxes at the whole point in space with respect to time to obtain the element mean at the half point of the next time layer, thereby obtaining the numerical simulation of the full flow field when stable.
具体实施例二Specific embodiment two
下面通过十个算例对本发明所公开的方法做进一步说明。The method disclosed in the present invention is further described below through ten calculation examples.
算例一Example 1
关于精度问题的测试。我们考虑一维Bugers方程：A test on precision issues. We consider the one-dimensional Bugers equation:
初值取u(x,0)＝0.5+sin(πx)，使用周期边界条件，图5分别给出了T＝0.5/π时CWENO-ZQ算法与CWENO-QS算法各自在L1与L∞下的精度与误差表现。The initial value is taken as u(x,0)=0.5+sin(πx), and the periodic boundary conditions are used. Figure 5 shows that the CWENO-ZQ algorithm and the CWENO-QS algorithm are at L1 and L∞ respectively when T=0.5/π. accuracy and error performance.
算例二Example 2
关于精度问题的测试。我们考虑一维Euler方程A test on precision issues. We consider the one-dimensional Euler equation
Ut+f(U)x＝0, (43)Ut +f(U)x =0, (43)
其中，U＝(ρ,ρu,E)T表示守恒变量,f(U)＝(ρu,ρu2+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∞下的精度与误差表现。Among them, U=(ρ, ρu, E)T represents the conserved variable, f(U)=(ρu, ρu2 +p, u(E+p))T , f(U) represents the flux, Ut represents U Derived from t, f(U)x represents the derivation of f(U) against x, ρ, p, u, E represent fluid density, pressure, horizontal velocity and energy, respectively, and T represents transposition. The initial value is ρ(x,0)=1+0.2sin(πx), v(x,0)=1, p(x,0)=1, using periodic boundary conditions, Figure 6 gives t=2 respectively The accuracy and error performance of CWENO-ZQ algorithm and CWENO-QS algorithm under L1 and L∞ respectively.
算例三Example 3
关于精度问题的测试。我们考虑二维Bugers方程A test on precision issues. We consider the two-dimensional Bugers equation
初值取u(x,y,0)＝0.5+sin(π(x+y)/2)，使用周期边界条件，图7分别给出了T＝0.5/π时CWENO-ZQ算法与CWENO-LP算法各自在L1与L∞下的精度与误差表现。The initial value is u(x,y,0)=0.5+sin(π(x+y)/2), and periodic boundary conditions are used. Figure 7 shows the CWENO-ZQ algorithm and CWENO- Accuracy and error performance of LP algorithms under L1 and L∞ respectively.
算例四Example 4
关于精度问题的测试。我们考虑二维Euler方程A test on precision issues. We consider the two-dimensional Euler equation
Ut+f(U)x+G(U)y＝0, (45)Ut +f(U)x +G(U)y =0, (45)
其中，t表示时间变量，x,y表示空间变量，U＝(ρ,ρu,ρv,E)T表示守恒变量，f(U)＝(ρu,ρu2+p,ρuv,u(E+p))T，g(U)＝(ρv,ρuv,ρv2+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∞下的精度与误差表现。Among them, t represents the time variable, x, y represents the space variable, U=(ρ, ρu, ρv, E)T represents the conserved variable, f(U)=(ρu, ρu2 +p, ρuv, u(E+p ))T , g(U)=(ρv,ρuv,ρv2 +p,v(E+p))T , f(U),g(U) represents flux, f(U)x represents f(U ) is derived from x, g(U)y means g(U) is derived from y, ρ, p, u, v, E represent fluid density, pressure, horizontal velocity, vertical velocity and energy, respectively, T represents Transpose. The initial value is ρ(x,y,0)=1+0.2sin(x+y), u(x,y,0)=1, v(x,y,0)=1, p(x,y, 0)=1, using periodic boundary conditions, Figure 8 shows the accuracy and error performance of the CWENO-ZQ algorithm and the CWENO-LP algorithm under L1 and L∞ respectively when t=2.
算例五Example 5
关于间断捕捉问题的测试。我们一维Euler方程(44)的Lax问题,其中初值取：A test for the discontinuous capture problem. The Lax problem of our one-dimensional Euler equation (44), where the initial value is taken as:
图9从上到下分别是密度图和密度的局部放大图，从左到右分别是CWENO-ZQ方法与CWENO-QS。方法图中的实线是精确解，方格是数值解，均取200个网格点。Figure 9 shows the density map and the partial enlarged view of the density from top to bottom, and the CWENO-ZQ method and CWENO-QS from left to right, respectively. The solid line in the method diagram is the exact solution, the square is the numerical solution, and 200 grid points are taken.
算例六Example 6
关于间断捕捉问题的测试。我们考虑一维Euler方程(44)的shock density waveinteraction问题,其中初值取：A test for the discontinuous capture problem. We consider the shock density waveinteraction problem of the one-dimensional Euler equation (44), where the initial value is:
图10从上到下分别是密度图和密度的局部放大图，从左到右分别是CWENO-ZQ方法与CWENO-QS。方法图中的实线是精确解，方格是数值解，均取400个网格点。Figure 10 shows the density map and the partial enlarged view of the density from top to bottom, and the CWENO-ZQ method and CWENO-QS from left to right. The solid line in the method diagram is the exact solution, the square is the numerical solution, and 400 grid points are taken.
算例七Example 7
关于间断捕捉问题的测试。我们考虑一维Euler方程(44)的two blast wave问题,其中初值取：A test for the discontinuous capture problem. We consider the two blast wave problem of the one-dimensional Euler equation (44), where the initial value is:
使用反射边界条件，CWENO-QS对于此算例并没有良好的表现，图11从左到右分别给出了t＝0.001时CWENO-ZQ的密度,速度以及压力的分布图，图中的实线是精确解，方格是数值解，均取800个网格点。Using the reflection boundary condition, CWENO-QS does not perform well for this example. Figure 11 shows the density, velocity and pressure distributions of CWENO-ZQ at t=0.001 from left to right. The solid line in the figure is the exact solution, and the squares are numerical solutions, all taking 800 grid points.
算例八Example 8
关于间断捕捉问题的测试。我们考虑Sedov blast wave问题,其中初值取：A test for the discontinuous capture problem. We consider the Sedov blast wave problem, where the initial value is:
使用外推边界条件，CWENO-QS对于此算例并没有良好的表现，图12从左到右分别给出了t＝0.001时CWENO-ZQ的密度、速度以及压力的分布图，图中的实线是精确解，方格是数值解，取800个网格点。Using the extrapolated boundary conditions, CWENO-QS does not perform well for this example. Figure 12 shows the density, velocity, and pressure distributions of CWENO-ZQ at t=0.001 from left to right. Lines are exact solutions, squares are numerical solutions, taking 800 grid points.
算例九Example 9
双马赫反射问题。该问题描述了一个与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个网格点。Double Mach reflection problem. The problem describes the changes that occur when a strong shock wave hits a reflecting wall at an angle of 60° to the x-axis, and the incoming flow is a strong shock wave with a Mach number of 10. The calculation area is [0,4]×[0,1]. the bottom of the area from The reflection boundary condition starts at y=0, and the other bottom boundaries (from x=0 to that part) is the wavefront condition. The CWENO-LP algorithm does not perform well for this example. Figure 13 shows the density contours of the CWENO-ZQ algorithm in the [0,3]×[0,1] region at t=0.2 from top to bottom. Figure and its partial enlargement, take 900×300 grid points.
算例十Example ten
台阶问题。该问题是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个网格点。Step problem. This problem is a classical example proposed by Emery in 1968 to test the nonlinear hyperbolic conservation law scheme. The initial data is that the horizontal incoming Mach number is 3, the density is 1.4, the horizontal velocity is 3, the vertical velocity is 0, the pressure is 1, the pipeline area is [0,3]×[0,1], and the distance is 0.6 from the left boundary. There is a step with a height of 0.2, and the step extends to the end of the pipe. The upper and lower boundaries are the reflection boundaries, the left boundary is the incoming flow boundary, and the right boundary is the outflow boundary. The CWENO-LP algorithm does not perform well for this example. Figure 14 shows the density contour map of the CWENO-ZQ algorithm at t=4, taking 600×200 grid points.
在本公开中参照附图来描述本发明的各方面，附图中示出了许多说明的实施例。本公开的实施例不必定义在包括本发明的所有方面。应当理解，上面介绍的多种构思和实施例，以及下面更加详细地描述的那些构思和实施方式可以以很多方式中任意一种来实施，这是因为本发明所公开的构思和实施例并不限于任何实施方式。另外，本发明公开的一些方面可以单独使用，或者与本发明公开的其他方面的任何适当组合来使用。Aspects of the invention are described in this disclosure with reference to the accompanying drawings, in which a number of illustrative embodiments are shown. Embodiments of the present disclosure are not necessarily defined to include all aspects of the invention. It should be understood that the various concepts and embodiments described above, as well as those described in greater detail below, can be implemented in any of a variety of ways, as the concepts and embodiments disclosed herein do not limited to any implementation. Additionally, some aspects of the present disclosure may be used alone or in any suitable combination with other aspects of the present disclosure.
虽然本发明已以较佳实施例揭露如上，然其并非用以限定本发明。本发明所属技术领域中具有通常知识者，在不脱离本发明的精神和范围内，当可作各种的更动与润饰。因此，本发明的保护范围当视权利要求书所界定者为准。Although the present invention has been disclosed above with preferred embodiments, it is not intended to limit the present invention. Those skilled in the art to which the present invention pertains can make various changes and modifications without departing from the spirit and scope of the present invention. Therefore, the protection scope of the present invention should be determined according to the claims.
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910710350.XACN110457806A (en) | 2019-08-02 | 2019-08-02 | Full flow field simulation method based on staggered grid with central fifth-order WENO scheme |
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910710350.XACN110457806A (en) | 2019-08-02 | 2019-08-02 | Full flow field simulation method based on staggered grid with central fifth-order WENO scheme |
| Publication Number | Publication Date |
|---|---|
| CN110457806Atrue CN110457806A (en) | 2019-11-15 |
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201910710350.XAPendingCN110457806A (en) | 2019-08-02 | 2019-08-02 | Full flow field simulation method based on staggered grid with central fifth-order WENO scheme |
| Country | Link |
|---|---|
| CN (1) | CN110457806A (en) |
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111159956A (en)* | 2019-12-10 | 2020-05-15 | 北京航空航天大学 | A feature-based method for capturing discontinuity in flow field |
| CN111177965A (en)* | 2019-12-25 | 2020-05-19 | 南京航空航天大学 | A fixed-point fast scan method based on multi-resolution WENO format for solving steady-state problems |
| CN111563314A (en)* | 2020-03-23 | 2020-08-21 | 空气动力学国家重点实验室 | Construction method of seven-point WENO format |
| CN111783343A (en)* | 2020-07-08 | 2020-10-16 | 中国人民解放军国防科技大学 | Method of selecting direction template for unstructured grid and solving method of flow field |
| CN112100835A (en)* | 2020-09-06 | 2020-12-18 | 西北工业大学 | High-efficiency high-precision numerical simulation method suitable for complex flow |
| CN114638173A (en)* | 2022-01-25 | 2022-06-17 | 中国空气动力研究与发展中心计算空气动力研究所 | High-order nonlinear shock wave capturing space dispersion method |
| CN116702571A (en)* | 2023-08-07 | 2023-09-05 | 中国空气动力研究与发展中心计算空气动力研究所 | Numerical simulation method and device based on multiple smoothness measurement factors |
| CN117922840A (en)* | 2024-01-22 | 2024-04-26 | 南京航空航天大学 | Wing performance test method based on ALW-MR-WENO algorithm |
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US7516054B2 (en)* | 2003-07-31 | 2009-04-07 | Canon Kabushiki Kaisha | Method and apparatus for numerical calculation related to fluids, program for implementing the method, and storage medium storing the program |
| CN105760602A (en)* | 2015-12-30 | 2016-07-13 | 南京航空航天大学 | Total flow field numerical simulation method for finite volume weighted essentially non-oscillatory scheme |
| CN110069854A (en)* | 2019-04-22 | 2019-07-30 | 南京航空航天大学 | Multiple resolution TWENO format is to the analogy method that can press flow field problems |
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US7516054B2 (en)* | 2003-07-31 | 2009-04-07 | Canon Kabushiki Kaisha | Method and apparatus for numerical calculation related to fluids, program for implementing the method, and storage medium storing the program |
| CN105760602A (en)* | 2015-12-30 | 2016-07-13 | 南京航空航天大学 | Total flow field numerical simulation method for finite volume weighted essentially non-oscillatory scheme |
| CN110069854A (en)* | 2019-04-22 | 2019-07-30 | 南京航空航天大学 | Multiple resolution TWENO format is to the analogy method that can press flow field problems |
| Title |
|---|
| JIANXIAN QIU等: "On the construction,comparison,and local characteristic decomposition for high-order central WENO schemes", 《JOURNAL OF COMPUTATIONAL PHYSICS》* |
| 吕梦迪: "求解双曲守恒律方程的五阶CWENO型熵稳定格式研究", 《中国优秀硕士学位论文全文数据库基础科学辑》* |
| 徐振礼等: "双曲守恒律方程的加权本质无振荡格式新进展", 《力学进展》* |
| 郭彦等: "一维Eular方程的特征有限体积格式", 《应用数学和力学》* |
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111159956B (en)* | 2019-12-10 | 2021-10-26 | 北京航空航天大学 | Feature-based flow field discontinuity capturing method |
| CN111159956A (en)* | 2019-12-10 | 2020-05-15 | 北京航空航天大学 | A feature-based method for capturing discontinuity in flow field |
| CN111177965A (en)* | 2019-12-25 | 2020-05-19 | 南京航空航天大学 | A fixed-point fast scan method based on multi-resolution WENO format for solving steady-state problems |
| CN111177965B (en)* | 2019-12-25 | 2022-06-17 | 南京航空航天大学 | A fixed-point fast scan method based on multi-resolution WENO format for solving steady-state problems |
| CN111563314B (en)* | 2020-03-23 | 2023-09-12 | 空气动力学国家重点实验室 | Seven-point WENO format construction method |
| CN111563314A (en)* | 2020-03-23 | 2020-08-21 | 空气动力学国家重点实验室 | Construction method of seven-point WENO format |
| CN111783343A (en)* | 2020-07-08 | 2020-10-16 | 中国人民解放军国防科技大学 | Method of selecting direction template for unstructured grid and solving method of flow field |
| CN111783343B (en)* | 2020-07-08 | 2024-12-24 | 中国人民解放军国防科技大学 | Unstructured grid orientation template selection method and flow field solution method |
| CN112100835A (en)* | 2020-09-06 | 2020-12-18 | 西北工业大学 | High-efficiency high-precision numerical simulation method suitable for complex flow |
| CN112100835B (en)* | 2020-09-06 | 2022-06-14 | 西北工业大学 | A high-efficiency and high-precision numerical simulation method for airfoil flow around an airfoil suitable for complex flow |
| CN114638173A (en)* | 2022-01-25 | 2022-06-17 | 中国空气动力研究与发展中心计算空气动力研究所 | High-order nonlinear shock wave capturing space dispersion method |
| CN114638173B (en)* | 2022-01-25 | 2023-06-02 | 中国空气动力研究与发展中心计算空气动力研究所 | Space discrete method for capturing high-order nonlinear shock waves |
| CN116702571A (en)* | 2023-08-07 | 2023-09-05 | 中国空气动力研究与发展中心计算空气动力研究所 | Numerical simulation method and device based on multiple smoothness measurement factors |
| CN116702571B (en)* | 2023-08-07 | 2023-10-20 | 中国空气动力研究与发展中心计算空气动力研究所 | Numerical simulation method and device based on multiple smoothness measurement factors |
| CN117922840A (en)* | 2024-01-22 | 2024-04-26 | 南京航空航天大学 | Wing performance test method based on ALW-MR-WENO algorithm |
| CN117922840B (en)* | 2024-01-22 | 2025-07-22 | 南京航空航天大学 | Wing performance test method based on ALW-MR-WENO algorithm |
| Publication | Publication Date | Title |
|---|---|---|
| CN110457806A (en) | Full flow field simulation method based on staggered grid with central fifth-order WENO scheme | |
| CN107220399A (en) | Weight the whole flow field analogy method of non-oscillatory scheme substantially based on Hermite interpolation | |
| CN108763683B (en) | New WENO format construction method under trigonometric function framework | |
| Restelli et al. | A conservative discontinuous Galerkin semi-implicit formulation for the Navier–Stokes equations in nonhydrostatic mesoscale modeling | |
| Troch et al. | An active wave generating–absorbing boundary condition for VOF type numerical model | |
| CN112100835B (en) | A high-efficiency and high-precision numerical simulation method for airfoil flow around an airfoil suitable for complex flow | |
| CN106767780B (en) | The extension ellipsoid set-membership filtering method approached based on Chebyshev polynomial interopolations | |
| CN109446471B (en) | A data transfer method for fluid-structure interaction interface considering load uncertainty | |
| CN112861445A (en) | Simulation method of grid-free numerical model | |
| CN105760602A (en) | Total flow field numerical simulation method for finite volume weighted essentially non-oscillatory scheme | |
| CN108280273A (en) | A kind of limited bulk Flow Field Numerical Calculation method based under non equidistance grid analysis | |
| Kurganov et al. | Adaptive moving mesh central-upwind schemes for hyperbolic system of PDEs: Applications to compressible Euler equations and granular hydrodynamics | |
| CN112307684B (en) | A fixed-point fast scanning method based on multi-resolution WENO format combined with ILW boundary processing | |
| Sonawane et al. | High resolution incompressible flow computations over unstructured mesh using SDWLS gradients | |
| Thai-Quang et al. | High-order alternating direction implicit method based on compact integrated-RBF approximations for unsteady/steady convection-diffusion equations | |
| CN111177965A (en) | A fixed-point fast scan method based on multi-resolution WENO format for solving steady-state problems | |
| CN114611423A (en) | Three-dimensional multiphase compressible fluid-solid coupling rapid calculation method | |
| Ahmed et al. | Finite Volume Method for a Time‐Dependent Convection‐Diffusion‐Reaction Equation with Small Parameters | |
| CN114492238B (en) | Depth coupling method for river and lake one-dimensional and plane two-dimensional hydrodynamic model | |
| Mandal et al. | High resolution schemes for genuinely two-dimensional HLLE Riemann solver | |
| CN105046324A (en) | Height anomaly fitting interpolation calculation method based on mobile neural network | |
| Zhan et al. | Meshfree lattice Boltzmann flux solver for compressible inviscid flows | |
| Zhang et al. | A new r‐ratio formulation for TVD schemes for vertex‐centered FVM on an unstructured mesh | |
| Hashemabadi et al. | Efficient Gridless Method Using Constrained Weights Optimization for Two-Dimensional Unsteady Inviscid Flows at Low Angles of Attack | |
| Wang et al. | Hermite radial basis collocation method for unsaturated soil water movement equation |
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| RJ01 | Rejection of invention patent application after publication | Application publication date:20191115 | |
| RJ01 | Rejection of invention patent application after publication |