技术领域Technical Field
本发明属于计算流体力学工程技术领域，具体涉及一种多重分辨率WENO格式结合ILW边界处理的定点快速扫描方法。The invention belongs to the technical field of computational fluid dynamics engineering, and in particular relates to a fixed-point fast scanning method combining a multi-resolution WENO format with ILW boundary processing.
背景技术Background Art
双曲守恒律方程的稳态问题是流体力学领域中经常遇到的数学问题，在工程应用中也占据重要的位置。因此，构造解决此类问题的高鲁棒性和高精度数值模拟方法变得很重要，很多学者也在这一领域攻坚克难。由于实际的问题当中计算区域多种多样，导致其物理量残差很难收敛到机器零。因此设计一种可以使残差下降到机器零的物面边界的高精度算法显得格外重要。在计算大型定常问题时，虽然大型超级计算机的广泛应用可以使计算对时间没那么敏感，但在算法设计中，程序运行的效率依然非常重要。The steady-state problem of hyperbolic conservation law equations is a mathematical problem often encountered in the field of fluid mechanics, and it also occupies an important position in engineering applications. Therefore, it is very important to construct a highly robust and high-precision numerical simulation method to solve such problems, and many scholars are also working hard in this field. Due to the variety of calculation areas in actual problems, it is difficult for the residuals of physical quantities to converge to machine zero. Therefore, it is particularly important to design a high-precision algorithm that can reduce the residual to the boundary of the object surface to machine zero. When calculating large-scale steady-state problems, although the widespread use of large supercomputers can make calculations less sensitive to time, the efficiency of program operation is still very important in algorithm design.
为捕捉流场内的激波，很多激波捕捉格式被设计出来。Godunov最先提出了一阶精度的激波捕捉格式。一阶精度的数值方法能抑制非物理振荡但会把大梯度处过度抹平，而强间断对物理现象的描述十分重要，因此需构造高阶精度的数值格式来更精确捕捉强间断。Harten于1983年首次提出了全变差减少(TVD)格式，之后Osher接着提出了本质无振荡格式(ENO)格式。ENO格式的主要思想是在逐次扩展的模板中选用最光滑的模板构造多项式求出单元半点的值，进而在光滑区域达到高阶精度，同时在间断附近实现基本无振荡的效果。但通过格式的构造过程可以看到，ENO格式是选用所有候选模板中的最优模板，其他的模板全部浪费掉，且数值精度越高浪费得越多，严重影响了计算效率。为了提高模板更有效的使用，Liu，Osher和Chan等于1994年提出了加权本质无振荡(WENO)格式，提高了一阶的计算精度。1996年，Jiang和Shu进一步改善了WENO格式，使得数值精度能够提高到2r-1阶，并设计出了光滑因子和非线性权的构造框架。WENO格式的主要思想是通过低阶重构通量的线性凸组合获得高阶近似。但该经典WENO格式的实现过程中，线性权的计算复杂，且在很多定常问题中残差无法下降到机器零。因此，2018年Zhu和Shu提出了多重分辨WENO格式，通过设计了不等距模板，让线性权可以任取为和为1的正数，且在光滑区域格式的数值精度保持最优，可使许多经典的定常问题算例残差可以下降到接近机器零。In order to capture shock waves in the flow field, many shock wave capture formats have been designed. Godunov first proposed a shock wave capture format with first-order accuracy. The numerical method with first-order accuracy can suppress non-physical oscillations but will over-smooth the large gradients. Strong discontinuities are very important for describing physical phenomena, so it is necessary to construct a numerical format with high-order accuracy to capture strong discontinuities more accurately. Harten first proposed the total variation reduction (TVD) format in 1983, and then Osher proposed the essentially non-oscillatory format (ENO) format. The main idea of the ENO format is to select the smoothest template in the successively expanded templates to construct a polynomial to obtain the value of the unit half point, thereby achieving high-order accuracy in the smooth area and achieving a basically non-oscillatory effect near the discontinuity. However, through the construction process of the format, it can be seen that the ENO format selects the best template among all candidate templates, and all other templates are wasted. The higher the numerical accuracy, the more waste, which seriously affects the computational efficiency. In order to improve the more effective use of templates, Liu, Osher and Chan et al. proposed the weighted essential non-oscillatory (WENO) format in 1994, which improved the first-order calculation accuracy. In 1996, Jiang and Shu further improved the WENO format, so that the numerical accuracy can be improved to the 2r-1 order, and designed a construction framework for smooth factors and nonlinear weights. The main idea of the WENO format is to obtain a high-order approximation through a linear convex combination of low-order reconstruction fluxes. However, in the implementation of this classic WENO format, the calculation of linear weights is complicated, and the residuals cannot be reduced to machine zero in many steady-state problems. Therefore, in 2018, Zhu and Shu proposed a multi-resolution WENO format. By designing unequal-distance templates, the linear weights can be arbitrarily taken as positive numbers with a sum of 1, and the numerical accuracy of the format in the smooth region is kept optimal, which can reduce the residuals of many classic steady-state problem examples to close to machine zero.
但是对于复杂的计算区域的数值模拟，有限差分方法最大的困难就是物面边界的处理。目前非贴体网格的主要方法有浸入边界法、嵌入边界法、SharpInterface方法、虚拟单元法。但这些方法都有各自的局限性，而且有一个共同的缺点，数值精度达不到高阶。于是2010年Sirui Tan提出了ILW边界处理方法，通过泰勒展开构造虚拟流场点上的值，耦合高阶WENO外推成功模拟了静止和运动物体的流动问题，数值结果显示对于静止问题可以达到五阶精度，对于运动问题可以达到三阶精度。之后ILW边界处理方法又被推广到了求带源项的欧拉方程，使得求解带粘性的NS方程也成为可能。后来又被改进为了简化的ILW边界处理方法，即高阶导数可直接由多项式外推得到，不用ILW过程依然不影响精度。However, for numerical simulations of complex computational areas, the biggest difficulty of the finite difference method is the processing of object surface boundaries. At present, the main methods for non-body-fitting grids are immersed boundary method, embedded boundary method, SharpInterface method, and virtual unit method. However, these methods have their own limitations and a common disadvantage, that is, the numerical accuracy cannot reach a high order. Therefore, in 2010, Sirui Tan proposed the ILW boundary processing method, which constructs the values on the virtual flow field points through Taylor expansion, and successfully simulates the flow problems of stationary and moving objects by coupling high-order WENO extrapolation. The numerical results show that the fifth-order accuracy can be achieved for stationary problems and the third-order accuracy can be achieved for moving problems. After that, the ILW boundary processing method was extended to solve the Euler equations with source terms, making it possible to solve the NS equations with viscosity. Later, it was improved to a simplified ILW boundary processing method, that is, high-order derivatives can be directly obtained by polynomial extrapolation, and the accuracy is not affected without the ILW process.
对于欧拉方程经典的三阶TVD Runge-Kutta时间离散，迭代次数较多和迭代的CPU时间较长，迭代效率有些低。为了提高迭代效率，于是在时间离散上提出了很多新的离散方法，比如快速行进算法和快速扫描算法。快速行进算法是让解严格按照递增或递减的顺序更新点值，这又会增加排序的时间复杂度。为了继续加快迭代效率就提出了快速扫描算法。与快速行进法相比，快速扫描法是一种并行的算法。由迎风性扫描覆盖了每一个特征方向，这样就省去了排序的时间，再结合Gauss-Seidel迭代，可以很大程度加快迭代的收敛。快速扫描法最开始使被用来求解静态的Hamilton-Jacobi方程。2016年Wu和Zhang将该快速扫描算法应用到求解双曲守恒律方程，也可以明显地加快格式的迭代速度。2020年Li和Zhu又将该算法与多重分辨WENO结合求解简单计算区域的定常流问题，达到了快速的完全收敛。For the classic third-order TVD Runge-Kutta time discretization of the Euler equation, the number of iterations is large and the CPU time of the iteration is long, and the iteration efficiency is somewhat low. In order to improve the iteration efficiency, many new discretization methods have been proposed in time discretization, such as the fast marching algorithm and the fast scanning algorithm. The fast marching algorithm allows the solution to update the point value strictly in increasing or decreasing order, which will increase the time complexity of sorting. In order to further speed up the iteration efficiency, the fast scanning algorithm is proposed. Compared with the fast marching method, the fast scanning method is a parallel algorithm. The upwind scanning covers every characteristic direction, which saves the sorting time. Combined with the Gauss-Seidel iteration, it can greatly speed up the convergence of the iteration. The fast scanning method was first used to solve the static Hamilton-Jacobi equation. In 2016, Wu and Zhang applied the fast scanning algorithm to solve the hyperbolic conservation law equation, which can also significantly speed up the iteration speed of the format. In 2020, Li and Zhu combined the algorithm with multi-resolution WENO to solve steady flow problems in simple computational areas, achieving rapid complete convergence.
发明内容Summary of the invention
本发明所要解决的技术问题是针对上述现有技术的不足，提供一种多重分辨率WENO格式结合ILW边界处理的定点快速扫描方法。The technical problem to be solved by the present invention is to provide a fixed-point fast scanning method combining a multi-resolution WENO format with ILW boundary processing in view of the above-mentioned deficiencies in the prior art.
为实现上述技术目的，本发明采取的技术方案为：In order to achieve the above technical objectives, the technical solution adopted by the present invention is:
一种多重分辨率WENO格式结合ILW边界处理的定点快速扫描方法，其特征在于，所述方法用于针对多种复杂计算区域的可压定常流场问题进行高精度数值模拟，包括：A fixed-point fast scanning method combining a multi-resolution WENO format with ILW boundary processing, characterized in that the method is used to perform high-precision numerical simulations on compressible constant flow field problems in various complex calculation areas, including:
步骤1.在笛卡尔坐标系下，把定常双曲守恒律问题转化为依赖时间的双曲守恒律问题，用新型的五阶ILW边界处理方法处理物面边界；Step 1. In the Cartesian coordinate system, the steady hyperbolic conservation law problem is transformed into a time-dependent hyperbolic conservation law problem, and the object surface boundary is processed using a novel fifth-order ILW boundary processing method;
步骤2.将双曲守恒律方程空间部分用有限差分多重分辨加权基本无振荡格式进行离散；Step 2. Discretize the spatial part of the hyperbolic conservation law equation using a finite-difference multi-resolution weighted essentially non-oscillatory scheme;
步骤3.对控制方程中的时间部分用三阶龙格库塔方法和定点快速扫描法离散成全离散的有限差分格式；Step 3. Discretize the time part of the control equation into a fully discrete finite difference format using the third-order Runge-Kutta method and the fixed-point fast scanning method;
步骤4.根据时空全离散方法得到下一时间层每一个点的近似值，依次迭代，得到计算区域内守恒变量的残差在趋于稳定时的数值模拟结果。Step 4. Obtain the approximate value of each point in the next time layer according to the full space-time discretization method, iterate in sequence, and obtain the numerical simulation results when the residual of the conserved variables in the calculation area tends to be stable.
为优化上述技术方案，采取的具体措施还包括：To optimize the above technical solutions, the specific measures taken also include:
上述的步骤1中，对于一维双曲守恒律方程：In step 1 above, for the one-dimensional hyperbolic conservation law equation:
其半离散格式的形式为：Its semi-discrete format is:
其中，U＝(ρ,ρu,E)T表示守恒变量，f(U)＝(ρu,ρu2+p,u(E+p))T表示通量，Ut表示U对t求导，f(U)x表示f(U)对x求导，ρ，u，p，E分别表示流体密度，速度，压强，能量，T表示转置，U0表示初始状态值，L(U)表示-fx(U)的空间离散形式；Among them, U＝(ρ,ρu,E)T represents the conserved variable, f(U)＝(ρu,ρu2 +p,u(E+p))T represents the flux, Ut represents the derivative of U with respect to t, f(U)x represents the derivative of f(U) with respect to x, ρ, u, p, E represent the fluid density, velocity, pressure, and energy respectively, T represents the transpose, U0 represents the initial state value, and L(U) represents the spatial discrete form of -fx (U);
把空间离散成统一长度的网格单元单元长度单元中心为其中i为坐标序号，a+h/2＝x0＜x1＜...＜xN＝b－h/2。Discretize space into grid cells of uniform length Unit length The unit center is Where i is the coordinate number, a+h/2=x0 <x1 <...<xN =b-h/2.
上述的步骤1中，用新型的五阶ILW边界处理方法处理物面边界，物面边界外的虚拟点用ILW边界处理,虚拟点由物面边界点经泰勒展开求出，包括：In the above step 1, the new fifth-order ILW boundary processing method is used to process the object surface boundary. The virtual points outside the object surface boundary are processed by ILW boundary processing. The virtual points are obtained by Taylor expansion of the object surface boundary points, including:
步骤1.1虚拟点用泰勒展开公式求出，U的第m个分量在点xN+1处的值为：Step 1.1 The virtual point is found using Taylor expansion formula. The value of the mth component of U at point xN+1 is:
然后通过已知的边界条件确定该边界是出流还是入流，如果U1(b,t)＝g1(t),U2(b,t)＝g2(t)，那么第一个和第二个分量为入流边界；Then, the known boundary conditions are used to determine whether the boundary is an outflow or an inflow. If U1 (b, t) = g1 (t), U2 (b, t) = g2 (t), then the first and second components are inflow boundaries;
入流边界用时间导数转换成空间导数求，出流边界用多项式外推求得；The inflow boundary is obtained by converting the time derivative into the space derivative, and the outflow boundary is obtained by polynomial extrapolation;
记L(UN)为的左特征向量；Let L(UN ) be The left eigenvector of ;
记(V3)j＝l3(UN)Uj，j＝N,N-1,...,N-4.用这五个点通过新的WENO外推求出k＝0,1,2Let (V3 )j = l3 (UN )Uj , j = N, N-1, ..., N-4. Use these five points to extrapolate the new WENO to obtain k＝0,1,2
3,4.用g1,g2和通过原方程求出边界点的各阶导数；3,4. Useg1 ,g2 and Calculate the derivatives of each order of the boundary point through the original equation;
步骤1.2通过原方程时间和空间的转换得到方程：Step 1.2: The equation is obtained by converting the original equation in time and space:
通过(4)这个线性方程求出边界值的一阶导数，其余各阶导数也同样的方式求出；The first-order derivative of the boundary value is obtained through the linear equation (4), and the other derivatives are obtained in the same way;
然后代入(3)式求出物面边界附近虚拟点的值。Then substitute into equation (3) to find the value of the virtual point near the object surface boundary.
上述的Vm*(k)的具体计算步骤如下：The specific calculation steps of the above Vm*(k) are as follows:
步骤a：选取三个模板T3＝[IN-4,IN-3,IN-2,IN-1,IN]，T2＝[IN-2,IN-1,IN]，T1＝[IN]。在每个模板上分别重构代数多项式q1(x)、q2(x)和q3(x)，使得其在单元边界有一阶，三阶，五阶精度；Step a: Select three templates T3 = [IN-4 ,IN-3 ,IN-2 ,IN-1 ,IN ], T2 = [IN-2 ,IN-1 ,IN ], T1 = [IN ]. Reconstruct the algebraic polynomials q1 (x), q2 (x) and q3 (x) on each template respectively, so that they have first-order, third-order and fifth-order accuracy at the unit boundary;
步骤b：计算光滑指示器βl，用于衡量重构多项式pl(x)在目标单元上的光滑度，计算公式为：Step b: Calculate the smoothness indicator βl , which is used to measure the smoothness of the reconstructed polynomial pl (x) on the target unit. The calculation formula is:
其中l＝2,3表示对应模板序号，表示多项式pl(x)对x的α阶导数，r＝2，但是β1＝0；Where l=2,3 represents the corresponding template number, represents the α-order derivative of the polynomial pl (x) with respect to x, r = 2, but β1 = 0;
步骤c：通过线性权γl和光滑指示器βl计算非线性权ωl，其计算公式为：Step c: Calculate the nonlinear weight ωl through the linear weight γl and the smooth indicator βl , and the calculation formula is:
其中l＝1,2,3表示对应模板序号，τ为计算过程中的过渡值，βl为光滑指示器，ε＝10-6防止分母为零；Where l = 1, 2, 3 represents the corresponding template number, τ is the transition value in the calculation process, βl is the smoothness indicator, and ε = 10-6 prevents the denominator from being zero;
步骤d：通过非线性权和多项式函数p,可以得到V函数的各阶导数：Step d: Through the nonlinear weight and polynomial function p, the derivatives of the V function can be obtained:
上述的步骤a具体过程如下：The specific process of step a above is as follows:
在三个模板T1、T2和T3上分别构造代数多项式q1(x)，q2(x)和q3(x)，使其满足:Construct algebraic polynomials q1 (x), q2 (x) and q3 (x) on three templates T1 , T2 and T3 respectively, so that they satisfy:
取线性权为:γ12＝1/11,γ22＝10/11，γ13＝1/111，γ23＝10/111，γ33＝100/111.The linear weights are: γ12 = 1/11, γ22 = 10/11, γ13 = 1/111, γ23 = 10/111, γ33 = 100/111.
重新构造出p1(x)，p2(x)和p3(x)，满足：Reconstruct p1 (x), p2 (x) and p3 (x) to satisfy:
p1(x)＝q1(x)， (8)p1 (x) = q1 (x), (8)
上述的步骤1中，对于二维双曲守恒律方程：In step 1 above, for the two-dimensional hyperbolic conservation law equation:
其半离散格式的形式为Its semi-discrete format is
其中，U＝(ρ,ρu,ρv,E)T表示守恒变量，F(U)，G(U)是通量，F(U)＝(ρu,ρu2+p,ρuv,u(E+p))TG(U)＝(ρu,ρuv,ρv2+p,v(E+p))T，Ut表示U对t求导，f(U)x表示f(U)对x求导，G(U)y表示f(U)对y求导，ρ，u，v,p,E分别表示流体密度，速度，压强，能量，T表示转置，U0表示初始状态值，L(U)表示-fx(U)-gy(U)的空间离散形式，其空间离散用一维多分辨率WENO离散的方法进行逐维离散；Among them, U＝(ρ,ρu,ρv,E)T represents the conserved variable, F(U), G(U) are fluxes, F(U)＝(ρu,ρu2 +p,ρuv,u(E+p))T G(U)＝(ρu,ρuv,ρv2 +p,v(E+p))T , Ut represents the derivative of U with respect to t, f(U)x represents the derivative of f(U) with respect to x, G(U)y represents the derivative of f(U) with respect to y, ρ, u, v, p, E represent the fluid density, velocity, pressure, and energy respectively, T represents the transpose, U0 represents the initial state value, L(U) represents the spatial discretization form of -fx (U)-gy (U), and its spatial discretization is discretized dimension by dimension using the one-dimensional multi-resolution WENO discretization method;
计算过程中需要用到物面边界外的虚拟点的值，根据所求的虚拟点的不同变换方程，对于点U(i,j)，记该点为P点；The value of the virtual point outside the object surface boundary is needed in the calculation process. According to the different transformation equations of the virtual point, for point U(i,j), the point is recorded as point P;
做P点与附近物面边界的外法线，外法线与物面边界交点记为P0，外法线与x轴之间的角度记为θ，进行下面坐标变换：Draw the external normal line between point P and the nearby object surface boundary. The intersection point of the external normal line and the object surface boundary is recorded as P0 . The angle between the external normal line and the x-axis is recorded as θ. Perform the following coordinate transformation:
原方程可变换成：The original equation can be transformed into:
Ut+F(U)x+G(U)y＝0.Ut +F(U)x +G(U)y = 0.
其中：in:
上述的步骤1中，用新型的五阶ILW边界处理方法处理物面边界，物面边界外的虚拟点用ILW边界处理,虚拟点由物面边界点经泰勒展开求出，包括：In the above step 1, the new fifth-order ILW boundary processing method is used to process the object surface boundary. The virtual points outside the object surface boundary are processed by ILW boundary processing. The virtual points are obtained by Taylor expansion of the object surface boundary points, including:
步骤1.1点P的值可用泰勒展开公式求出，U的第m个分量在点P处的值为：Step 1.1 The value of point P can be found using the Taylor expansion formula. The value of the mth component of U at point P is:
其中Δ为点P0到点P1的距离，然后确定边界出流还是入流，如果U2在P0点的值已知为g2(t)；入流边界用时间导数转换成空间导数求，出流边界用多项式外推求得；记P0点在计算区域内靠的最近的点为P1，按这个最近的点将流场内空间模板上的点都进行特征分解，然后记的左特征向量；记(Vm)r,s＝lmUr,s，(r,s)∈Θ.Θ是构造新的WENO外推所需要用到的模板，求出用g2和可以通过原方程求出边界外虚拟点的各阶导数；Where Δ is the distance from point P0 to point P1. Then determine whether the boundary is outflow or inflow. If the value of U2 at point P0 is known as g2 (t), the inflow boundary is converted from the time derivative to the space derivative, and the outflow boundary is obtained by polynomial extrapolation. The nearest point to point P0 in the calculation area is P1. According to this nearest point, all points on the space template in the flow field are characteristically decomposed, and then recorded as The left eigenvector of ; Let (Vm )r,s = lm Ur,s , (r,s)∈Θ.Θ is the template needed to construct a new WENO extrapolation. Use g2 and The derivatives of various orders of virtual points outside the boundary can be obtained through the original equation;
步骤1.2通过原方程时间和空间的转换得到方程：Step 1.2: The equation is obtained by converting the original equation in time and space:
通过(29)这个线性方程求出边界值的一阶导数，其余各阶导数以同样的方式求出；The first-order derivative of the boundary value is obtained through the linear equation (29), and the remaining derivatives are obtained in the same way;
然后代入(28)式就可以求出物面边界附近虚拟点的值。Then, by substituting into equation (28), we can find the value of the virtual point near the boundary of the object surface.
上述的步骤2中，一维情形时将方程写成守恒形式：In step 2 above, the equation is written in conservation form in one dimension:
其中，和分别表示通量f(U)在目标单元Ii的边界和处的五阶近似的数值通量，Ui(t)表示U在网格单元Ii内点xi处的值U(xi,t)。in, and They represent the flux f(U) at the boundary of target unit Ii and The numerical flux of the fifth-order approximation at , Ui (t) represents the value U(xi ,t) of U at point xi in grid cell Ii .
上述的步骤2中，求通量f(U)在目标单元Ii的边界和处的五阶近似值具体步骤如下：In step 2 above, find the flux f(U) at the boundary of the target unit Ii and The fifth-order approximation at The specific steps are as follows:
步骤2.1用Lax-Friedrichs分裂把通量分裂为其中f+(U)在点处的重构过程及定义为Step 2.1 Use Lax-Friedrichs splitting to split the flux into in f+ (U) at point The reconstruction process and definition of
步骤2.2将目标单元Ii构造大模板T3＝[Ii-2,Ii-1,Ii,Ii+1,Ii+2]，小模板T1＝[Ii]和小模板T2＝[Ii-1,Ii,Ii+1]，其中Ii为对应序号的网格单元；Step 2.2 constructs the target unit Ii into a large template T3 = [Ii-2 , Ii-1 , Ii , Ii+1 , Ii+2 ], a small template T1 = [Ii ] and a small template T2 = [Ii-1 , Ii , Ii+1 ], where Ii is the grid unit of the corresponding sequence number;
步骤2.3在每个模板上分别重构代数多项式q1(x)、q2(x)和q3(x)，使得其在单元边界有一阶、三阶和五阶精度；Step 2.3 reconstructs the algebraic polynomials q1 (x), q2 (x) and q3 (x) on each template respectively, so that they have first-order, third-order and fifth-order accuracy at the unit boundary;
步骤2.4任意取线性权为：γ12，γ22，γ13，γ23，γ33；Step 2.4: arbitrarily select linear weights as: γ12 , γ22 , γ13 , γ23 , γ33 ;
将多项式q1(x)、q2(x)和q3(x)重组为新的多项式p1(x)，p2(x)和p3(x)，满足Reorganize the polynomials q1 (x), q2 (x) and q3 (x) into new polynomials p1 (x), p2 (x) and p3 (x) that satisfy
步骤2.5计算光滑指示器βl，用于衡量重构多项式pl(x)在目标单元上的光滑度，计算公式为：Step 2.5: Calculate the smoothness indicator βl , which is used to measure the smoothness of the reconstruction polynomial pl (x) on the target unit. The calculation formula is:
其中l＝2,3表示对应模板序号，表示多项式pl(x)对x的α阶导数，r＝2；Where l=2,3 represents the corresponding template number, represents the α-order derivative of the polynomial pl (x) with respect to x, r = 2;
其中：γ1,1＝1-γ0,1，in: γ1,1 ＝1-γ0,1 ，
步骤2.6通过线性权γl和光滑指示器βl计算非线性权ωl，其计算公式为：Step 2.6 Calculate the nonlinear weight ωl through the linear weight γl and the smooth indicator βl , and the calculation formula is:
其中l＝1,2,3表示对应模板序号，τ为计算过程中的过渡值，βl为光滑指示器，ε＝10-6防止分母为零；Where l = 1, 2, 3 represents the corresponding template number, τ is the transition value in the calculation process, βl is the smoothness indicator, and ε = 10-6 prevents the denominator from being zero;
步骤2.7求出数值通量分裂f+(U)在点处的多重分辨WENO重构值：Step 2.7 Find the numerical flux splitting f+ (U) at the point The multi-resolution WENO reconstruction value at:
从而得到空间离散的常微分方程。This results in a spatially discretized ordinary differential equation.
上述的步骤3为：Step 3 above is:
使用三阶TVD Runge-Kutta时间离散公式和定点快速扫描公式将半离散有限差分格式离散成时空全离散有限差分格式，所述时空全离散有限差分格式为关于时间层的迭代公式；Using a third-order TVD Runge-Kutta time discretization formula and a fixed-point fast scanning formula, the semi-discrete finite difference format is discretized into a time-space fully discrete finite difference format, wherein the time-space fully discrete finite difference format is an iterative formula with respect to the time layer;
所述步骤3中，时间离散的三阶TVD Runge-Kutta时间离散公式：In step 3, the time-discrete third-order TVD Runge-Kutta time-discrete formula is:
其中，U(1),U(2)为中间过渡值，Δt为时间步长，上标n表示第n时间层，L(Un),L(U(1)),L(U(2))为高阶空间离散形式的近似值；Among them, U(1) , U(2) are intermediate transition values, Δt is the time step, the superscript n represents the nth time layer, and L(Un ), L(U(1) ), L(U(2) ) are approximate values of the high-order spatial discrete form;
时间离散的定点快速扫描法格式为：The format of the time-discrete fixed-point fast scanning method is:
其中*代表有n+1层的新值就用n+1层的新值，没有新值就用第n层的旧值，扫描的顺序为：The * indicates that if there is a new value in the n+1 layer, the new value in the n+1 layer is used, and if there is no new value, the old value in the nth layer is used. The scanning order is:
i＝1,…,N和i＝N,…,1交替扫描。i＝1,…,N and i＝N,…,1 are scanned alternately.
本发明具有以下有益效果：The present invention has the following beneficial effects:
相比于经典的ILW边界处理方法，新的WENO(Multi-resolution WeightedEssentially Non-Oscillatory，分多分辨加权基本无振荡格式)外推改进了线性权与网格尺度相关的问题，线性权可以任意取值。而且节省了使用的空间模板数量，减少了非线性权中参数的选取。使用起来更加方便简单，也有更广泛的适用性。Compared with the classic ILW boundary processing method, the new WENO (Multi-resolution Weighted Essentially Non-Oscillatory) extrapolation improves the problem of linear weights and grid scales. The linear weights can be arbitrarily selected. It also saves the number of spatial templates used and reduces the selection of parameters in the nonlinear weights. It is more convenient and simple to use and has a wider range of applicability.
相比于经典WENO格式，该多重分辨WENO格式通过采用一系列不等距中心模版的信息使定常问题的残差下降得更快且其值能接近于机器零。该格式能精确捕捉激波，且在解的光滑区域能保持最优数值精度。能任意选择线性权的取值，在减少计算量的同时不降低格式在解的光滑区域的数值精度。Compared with the classic WENO format, the multi-resolution WENO format uses a series of unequally spaced central template information to make the residual of the steady-state problem decrease faster and its value can be close to the machine zero. The format can accurately capture shock waves and maintain the optimal numerical accuracy in the smooth region of the solution. The value of the linear weight can be arbitrarily selected, which reduces the amount of calculation while not reducing the numerical accuracy of the format in the smooth region of the solution.
相比于经典的三阶Runge-Kutta时间离散，快速扫描法可以取较大的CFL数，可极大减少格式的迭代次数和节省大量的CPU时间，并且该方法易于推广到高维情形。Compared with the classic third-order Runge-Kutta time discretization, the fast scanning method can take a larger CFL number, which can greatly reduce the number of iterations of the format and save a lot of CPU time, and this method can be easily extended to high-dimensional situations.
新型ILW(Inverse Lax-Wendroff)边界处理方法与原ILW边界处理方法相比，能在保持物面边界附近计算格式高阶精度不降低的优点，同时能有效减少空间模板的数量和减少需要人为调整的参数，在求解定常问题时物理量的残差可以快速收敛。此新型ILW边界处理方法在结合了定点快速扫描算法和多分辨加权基本无振荡格式后不但能保持原有高精度数值计算格式在激波捕捉时的各种优点，还可以明显加快物理量残差收敛到机器零的速度。最后，本发明采用上述方法对多个经典的二维定常可压流场问题进行了高精度数值模拟，充分验证了本发明的有效性和可靠性。Compared with the original ILW boundary processing method, the new ILW (Inverse Lax-Wendroff) boundary processing method can maintain the advantage of not reducing the high-order accuracy of the calculation format near the object surface boundary, and at the same time can effectively reduce the number of spatial templates and reduce the parameters that need to be manually adjusted. When solving steady-state problems, the residuals of physical quantities can converge quickly. This new ILW boundary processing method, after combining the fixed-point fast scanning algorithm and the multi-resolution weighted substantially non-oscillation format, can not only maintain the various advantages of the original high-precision numerical calculation format in shock wave capture, but also significantly speed up the convergence of the physical quantity residuals to machine zero. Finally, the present invention uses the above method to perform high-precision numerical simulations on multiple classic two-dimensional steady compressible flow field problems, which fully verifies the effectiveness and reliability of the present invention.
附图说明BRIEF DESCRIPTION OF THE DRAWINGS
图1为实施例一数值结果，其中图1a和图1c是实施例一中采用多重分辨WENO空间离散，龙格库塔时间离散，得到的CFL＝0.6时的压力等值线和残差的下降曲线；图1b和图1d是实施例一中采用多重分辨WENO空间离散，定点快速扫描法空间离散，得到的CFL＝0.6时的压力等值线和残差的下降曲线。Figure 1 is the numerical results of Example 1, wherein Figures 1a and 1c are the pressure contour lines and residual drop curves when CFL=0.6 obtained by using multi-resolution WENO spatial discretization and Runge-Kutta time discretization in Example 1; Figures 1b and 1d are the pressure contour lines and residual drop curves when CFL=0.6 obtained by using multi-resolution WENO spatial discretization and fixed-point fast scanning method spatial discretization in Example 1.
图2为实施例二数值结果，其中图2a和图2c是实施例二中采用多重分辨WENO空间离散，龙格库塔时间离散，得到的CFL＝0.5时的压力等值线和残差的下降曲线；图2b和图2d是实施例二中采用多重分辨WENO空间离散，定点快速扫描法空间离散，得到的CFL＝0.5时的压力等值线和残差的下降曲线。Figure 2 is the numerical results of Example 2, where Figures 2a and 2c are the pressure contour lines and residual drop curves when CFL=0.5 obtained by using multi-resolution WENO spatial discretization and Runge-Kutta time discretization in Example 2; Figures 2b and 2d are the pressure contour lines and residual drop curves when CFL=0.5 obtained by using multi-resolution WENO spatial discretization and fixed-point rapid scanning method spatial discretization in Example 2.
图3为实施例三数值结果，其中图3a和图3c是实施例三中采用多重分辨WENO空间离散，龙格库塔时间离散，得到的CFL＝1.0时的压力等值线和残差的下降曲线；图3b和图3d是实施例三中采用多重分辨WENO空间离散，定点快速扫描法空间离散，得到的CFL＝1.0时的压力等值线和残差的下降曲线。Figure 3 is the numerical results of Example 3, where Figures 3a and 3c are the pressure contour lines and residual drop curves when CFL=1.0 obtained by using multi-resolution WENO spatial discretization and Runge-Kutta time discretization in Example 3; Figures 3b and 3d are the pressure contour lines and residual drop curves when CFL=1.0 obtained by using multi-resolution WENO spatial discretization and fixed-point fast scanning method spatial discretization in Example 3.
图4为实施例四数值结果，其中图4a和图4c是实施例四中采用多重分辨WENO空间离散，龙格库塔时间离散，得到的CFL＝0.9时的压力等值线和残差的下降曲线；图4b和图4d是实施例四中采用多重分辨WENO空间离散，定点快速扫描法空间离散，得到的CFL＝1.0时的压力等值线和残差的下降曲线。Figure 4 is the numerical results of Example 4, where Figures 4a and 4c are the pressure contour lines and residual drop curves when CFL=0.9 obtained by using multi-resolution WENO spatial discretization and Runge-Kutta time discretization in Example 4; Figures 4b and 4d are the pressure contour lines and residual drop curves when CFL=1.0 obtained by using multi-resolution WENO spatial discretization and fixed-point fast scanning method spatial discretization in Example 4.
图5为实施例五数值结果，其中图5a和图5c是实施例五中采用多重分辨WENO空间离散，龙格库塔时间离散，得到的CFL＝1.7时的压力等值线和残差的下降曲线；图5b和图5d是实施例五中采用多重分辨WENO空间离散，定点快速扫描法空间离散，得到的CFL＝1.3时的压力等值线和残差的下降曲线。Figure 5 is the numerical results of Example 5, where Figures 5a and 5c are the pressure contour lines and residual drop curves when CFL=1.7 obtained by using multi-resolution WENO spatial discretization and Runge-Kutta time discretization in Example 5; Figures 5b and 5d are the pressure contour lines and residual drop curves when CFL=1.3 obtained by using multi-resolution WENO spatial discretization and fixed-point fast scanning method spatial discretization in Example 5.
图6为实施例六数值结果，其中图6a和图6c是实施例五中采用多重分辨WENO空间离散，龙格库塔时间离散，得到的CFL＝1.6时的压力等值线和残差的下降曲线；图6b和图6d是实施例五中采用多重分辨WENO空间离散，定点快速扫描法空间离散，得到的CFL＝1.2时的压力等值线和残差的下降曲线。Figure 6 is the numerical results of Example 6, where Figures 6a and 6c are the pressure contour lines and residual drop curves when CFL=1.6 obtained by using multi-resolution WENO spatial discretization and Runge-Kutta time discretization in Example 5; Figures 6b and 6d are the pressure contour lines and residual drop curves when CFL=1.2 obtained by using multi-resolution WENO spatial discretization and fixed-point rapid scanning method spatial discretization in Example 5.
图7为实施例七数值结果，其中图7a和图7c是实施例五中采用多重分辨WENO空间离散，龙格库塔时间离散，得到的CFL＝0.7时的压力等值线和残差的下降曲线；图7b和图7d是实施例五中采用多重分辨WENO空间离散，定点快速扫描法空间离散，得到的CFL＝0.7时的压力等值线和残差的下降曲线。Figure 7 is the numerical results of Example 7, where Figures 7a and 7c are the pressure contour lines and residual drop curves when CFL=0.7 obtained by using multi-resolution WENO spatial discretization and Runge-Kutta time discretization in Example 5; Figures 7b and 7d are the pressure contour lines and residual drop curves when CFL=0.7 obtained by using multi-resolution WENO spatial discretization and fixed-point fast scanning method spatial discretization in Example 5.
图8是本发明方法流程图。FIG8 is a flow chart of the method of the present invention.
具体实施方式DETAILED DESCRIPTION
以下结合附图对本发明的实施例作进一步详细描述。The embodiments of the present invention are further described in detail below in conjunction with the accompanying drawings.
参见图8，本发明的一种多重分辨率WENO格式结合ILW边界处理的定点快速扫描方法，用于针对多种复杂计算区域的可压定常流场问题进行高精度数值模拟，包括：Referring to FIG8 , a fixed-point fast scanning method of the present invention using a multi-resolution WENO format combined with ILW boundary processing is used to perform high-precision numerical simulations on compressible constant flow field problems in various complex calculation areas, including:
步骤1.在笛卡尔坐标系下，把定常双曲守恒律问题转化为依赖时间的双曲守恒律问题，用新型的五阶ILW边界处理方法处理物面边界；Step 1. In the Cartesian coordinate system, the steady hyperbolic conservation law problem is transformed into a time-dependent hyperbolic conservation law problem, and the object surface boundary is processed using a novel fifth-order ILW boundary processing method;
先考虑一维双曲守恒律方程：First consider the one-dimensional hyperbolic conservation law equation:
其半离散格式的形式为Its semi-discrete format is
其中，U＝(ρ,ρu,E)T表示守恒变量，f(U)＝(ρu,ρu2+p,u(E+p))T表示通量，Ut表示U对t求导，f(U)x表示f(U)对x求导，ρ，u，p，E分别表示流体密度，速度，压强，能量，T表示转置，U0表示初始状态值，L(U)表示-fx(U)的空间离散形式。Among them, U＝(ρ,ρu,E)T represents the conserved variable, f(U)＝(ρu,ρu2 +p,u(E+p))T represents the flux, Ut represents the derivative of U with respect to t, f(U)x represents the derivative of f(U) with respect to x, ρ, u, p, E represent fluid density, velocity, pressure, and energy respectively, T represents transpose, U0 represents the initial state value, and L(U) represents the spatial discrete form of -fx (U).
把空间离散成统一长度的网格单元单元长度单元中心为其中i为坐标序号，a+h/2＝x0＜x1＜...＜xN＝b－h/2。Discretize space into grid cells of uniform length Unit length The unit center is Where i is the coordinate number, a+h/2=x0 <x1 <...<xN =b-h/2.
用新型的五阶ILW边界处理方法处理物面边界，物面边界外的虚拟点用ILW边界处理,虚拟点由物面边界点经泰勒展开求出，包括：The new fifth-order ILW boundary processing method is used to process the object surface boundary. The virtual points outside the object surface boundary are processed by ILW boundary processing. The virtual points are obtained by Taylor expansion of the object surface boundary points, including:
步骤1.1这些虚拟点用泰勒展开公式求出，U的第m个分量在点xN+1处的值为：Step 1.1 These virtual points are solved using the Taylor expansion formula. The value of the mth component of U at point xN+1 is:
然后通过已知的边界条件确定该边界是出流还是入流，如果U1(b,t)＝g1(t),U2(b,t)＝g2(t)。那么第一个和第二个分量为入流边界。Then, the boundary is determined to be outflow or inflow by the known boundary conditions. If U1 (b, t) = g1 (t), U2 (b, t) = g2 (t), then the first and second components are inflow boundaries.
入流边界用时间导数转换成空间导数求，出流边界用多项式外推求得。记L(UN)为的左特征向量。记(V3)j＝l3(UN)Uj，j＝N,N-1,...,N-4.用这五个点通过新的WENO外推求出k＝0,1,23,4.用g1,g2和可以通过原方程求出边界点的各阶导数。The inflow boundary is obtained by converting the time derivative into the space derivative, and the outflow boundary is obtained by polynomial extrapolation. Let L(UN ) be The left eigenvector of . Let (V3 )j = l3 (UN )Uj , j = N, N-1, ..., N-4. Use these five points to extrapolate the new WENO to obtain k＝0,1,23,4. Useg1 ,g2 and The derivatives of each order at the boundary point can be obtained through the original equation.
步骤1.2通过原方程时间和空间的转换可以得到方程：Step 1.2 The equation can be obtained by converting the original equation in time and space:
通过(4)这个线性方程求出边界值的一阶导数，其余各阶导数也可以同样的方式求出。然后代入(3)式就可以求出物面边界附近虚拟点的值。The first-order derivative of the boundary value can be obtained through the linear equation (4), and the other derivatives can be obtained in the same way. Then, substituting it into equation (3), the value of the virtual point near the boundary of the object surface can be obtained.
求的具体步骤如下：beg The specific steps are as follows:
步骤a：选取三个模板T3＝[IN-4,IN-3,IN-2,IN-1,IN]，T2＝[IN-2,IN-1,IN]，T1＝[IN]。在每个模板上分别重构代数多项式q1(x)、q2(x)和q3(x)，使得他们在单元边界有一阶，三阶，五阶精度。Step a: Select three templates T3 = [IN-4 , IN-3 ,IN-2 , IN-1 ,IN ], T2 = [IN-2 ,IN-1 ,IN ], T1 = [IN ]. Reconstruct the algebraic polynomials q1 (x), q2 (x) and q3 (x) on each template so that they have first-order, third-order and fifth-order accuracy on the unit boundary.
其具体过程如下：The specific process is as follows:
在三个模板T1、T2和T3上分别构造代数多项式q1(x)，q2(x)和q3(x)，使其满足:Construct algebraic polynomials q1 (x), q2 (x) and q3 (x) on three templates T1 , T2 and T3 respectively, so that they satisfy:
任意取线性权为:γ12＝1/11,γ22＝10/11，γ13＝1/111，γ23＝10/111，γ33＝100/111.重新构造出p1(x)，p2(x)和p3(x)，满足：Arbitrarily take the linear weights as: γ12 = 1/11, γ22 = 10/11, γ13 = 1/111, γ23 = 10/111, γ33 = 100/111. Reconstruct p1 (x), p2 (x) and p3 (x) to satisfy:
p1(x)＝q1(x)， (8)p1 (x) = q1 (x), (8)
步骤b：计算光滑指示器βl，用于衡量重构多项式pl(x)在目标单元上的光滑度，计算公式为：Step b: Calculate the smoothness indicator βl , which is used to measure the smoothness of the reconstructed polynomial pl (x) on the target unit. The calculation formula is:
其中l＝2,3表示对应模板序号，表示多项式pl(x)对x的α阶导数，r＝2，但是β1＝0。Where l=2,3 represents the corresponding template number, represents the α-order derivative of the polynomial pl (x) with respect to x, r=2, but β1 =0.
步骤c：通过线性权γl和光滑指示器βl计算非线性权ωl，其计算公式为：Step c: Calculate the nonlinear weight ωl through the linear weight γl and the smooth indicator βl , and the calculation formula is:
其中l＝1,2,3表示对应模板序号，τ为计算过程中的过渡值，βl为光滑指示器，ε＝10-6防止分母为零。Where l = 1, 2, 3 represents the corresponding template number, τ is the transition value in the calculation process, βl is the smoothness indicator, and ε=10-6 prevents the denominator from being zero.
步骤d：通过非线性权和多项式函数p,可以得到V函数的各阶导数。Step d: Through the nonlinear weights and polynomial function p, the derivatives of the V function can be obtained.
得到的各阶导数，代入式(3)中可以得到所需的所有虚拟点的点值。Substituting the obtained derivatives of various orders into formula (3) can obtain the point values of all the required virtual points.
再考虑二维双曲守恒律方程：Consider the two-dimensional hyperbolic conservation law equation:
其半离散格式的形式为Its semi-discrete format is
其中，U＝(ρ,ρu,ρv,E)T表示守恒变量，F(U)，G(U)是通量，F(U)＝(ρu,ρu2+p,ρuv,u(E+p))TG(U)＝(ρu,ρuv,ρv2+p,v(E+p))T，Ut表示U对t求导，f(U)x表示f(U)对x求导，G(U)y表示f(U)对y求导，ρ，u，v,p,E分别表示流体密度，速度，压强，能量，T表示转置，U0表示初始状态值，L(U)表示-fx(U)-gy(U)的空间离散形式。其空间离散用一维多分辨率WENO离散的方法进行逐维离散。Among them, U＝(ρ,ρu,ρv,E)T represents the conserved variables, F(U), G(U) are fluxes, F(U)＝(ρu,ρu2 +p,ρuv,u(E+p))T G(U)＝(ρu,ρuv,ρv2 +p,v(E+p))T , Ut represents the derivative of U with respect to t, f(U)x represents the derivative of f(U) with respect to x, G(U)y represents the derivative of f(U) with respect to y, ρ, u, v, p, E represent the fluid density, velocity, pressure, and energy respectively, T represents the transposition, U0 represents the initial state value, and L(U) represents the spatial discretization form of -fx (U)-gy (U). Its spatial discretization is discretized dimension by dimension using the one-dimensional multi-resolution WENO discretization method.
计算过程中需要用到物面边界外的虚拟点的值，根据所求的虚拟点的不同变换方程。以求点U(i,j)为例，记该点为P点。做P点与附近物面边界的外法线，外法线与物面边界交点记为P0，外法线与x轴之间的角度记为θ。进行下面坐标变换：The calculation process requires the value of the virtual point outside the object surface boundary, and the equation is transformed according to the different virtual points. Take point U(i,j) as an example, and record this point as point P. Draw the external normal between point P and the nearby object surface boundary, the intersection of the external normal and the object surface boundary is recorded as P0 , and the angle between the external normal and the x-axis is recorded as θ. Perform the following coordinate transformation:
原方程可变换成：The original equation can be transformed into:
Ut+F(U)x+G(U)y＝0.Ut +F(U)x +G(U)y = 0.
其中：in:
步骤1中，用新型的五阶ILW边界处理方法处理物面边界，物面边界外的虚拟点用ILW边界处理,虚拟点由物面边界点经泰勒展开求出，包括：In step 1, the object surface boundary is processed using a novel fifth-order ILW boundary processing method. Virtual points outside the object surface boundary are processed using the ILW boundary processing. Virtual points are obtained by Taylor expansion of the object surface boundary points, including:
步骤1.1点P的值可用泰勒展开公式求出，U的第m个分量在点P处的值为：Step 1.1 The value of point P can be found using the Taylor expansion formula. The value of the mth component of U at point P is:
其中Δ为点P0到点P1的距离，然后确定边界出流还是入流，如果U2在P0点的值已知为g2(t)。入流边界用时间导数转换成空间导数求，出流边界用多项式外推求得。记P0点在计算区域内靠的最近的点为P1，按这个最近的点将流场内空间模板上的点都进行特征分解，然后记l(U)为的左特征向量。记(Vm)r,s＝lmUr,s，(r,s)∈Θ.Θ是构造新的WENO外推所需要用到的模板，求出k＝0,1,2,3,4.用g2和可以通过原方程求出边界外虚拟点的各阶导数。Where Δ is the distance from point P0 to point P1. Then determine whether the boundary is outflow or inflow, if the value of U2 at point P0 is known as g2 (t). The inflow boundary is converted from the time derivative to the space derivative, and the outflow boundary is obtained by polynomial extrapolation. The nearest point of point P0 in the calculation area is denoted as P1. According to this nearest point, all points on the space template in the flow field are characteristically decomposed, and then l(U) is denoted as The left eigenvector of . Let (Vm )r,s = lm Ur,s , (r,s)∈Θ.Θ is the template needed to construct a new WENO extrapolation. k＝0,1,2,3,4. Use g2 and The derivatives of various orders of virtual points outside the boundary can be obtained through the original equation.
步骤1.2通过原方程时间和空间的转换可以得到方程：Step 1.2 The equation can be obtained by converting the original equation in time and space:
通过(29)这个线性方程求出边界值的一阶导数，其余各阶导数也可以同样的方式求出。然后代入(28)式就可以求出物面边界附近虚拟点的值。The first-order derivative of the boundary value can be obtained through the linear equation (29), and the other derivatives can be obtained in the same way. Then, substituting it into equation (28), the value of the virtual point near the boundary of the object surface can be obtained.
求Vm*(k)的具体步骤如下：The specific steps to calculate Vm*(k) are as follows:
步骤a：选取三个模板T1，T2，T3分别包括1个点，9个点，25个点。模板T1的点就选择点P1，模板T2的点选择计算区域内P0附近的九个点，每层三个点，模板T3的点选择计算区域内P0附近的二十五个点，每层五个点。在每个模板上像一维一样分别重构出二维代数多项式q1(x,y),q2(x,y)和q3(x,y)。Step a: Select three templates T1 , T2 , and T3 , including 1 point, 9 points, and 25 points respectively. The point of template T1 is point P1 , the point of template T2 is nine points near P0 in the calculation area, three points in each layer, and the point of template T3 is twenty-five points near P0 in the calculation area, five points in each layer. Reconstruct the two-dimensional algebraic polynomials q1 (x, y), q2 (x, y), and q3 (x, y) on each template as if they were one-dimensional.
任意取线性权为：γ12＝1/11，γ22＝10/11，γ13＝1/111，γ23＝10/111，γ33＝100/111。重新构造出p1(x)，p2(x)和p3(x)，满足：Arbitrarily choose linear weights: γ12 = 1/11, γ22 = 10/11, γ13 = 1/111, γ23 = 10/111, γ33 = 100/111. Reconstruct p1 (x), p2 (x) and p3 (x) to satisfy:
p1(x,y)＝q1(x,y), (17)p1 (x, y) = q1 (x, y), (17)
步骤b：计算光滑指示器βl，用于衡量重构多项式pl(x,y)在目标单元上的光滑度，计算公式为：Step b: Calculate the smoothness indicator βl , which is used to measure the smoothness of the reconstructed polynomial pl (x, y) on the target unit. The calculation formula is:
其中l＝2,3表示对应模板序号，r＝2,β1＝0。Wherein l=2, 3 represents the corresponding template number, r=2, β1 =0.
步骤c：通过线性权γl和光滑指示器βl计算非线性权ωl，其计算公式为：Step c: Calculate the nonlinear weight ωl through the linear weight γl and the smooth indicator βl , and the calculation formula is:
其中l＝1,2,3表示对应模板序号，τ为计算过程中的过渡值，βl为光滑指示器，ε＝10-6防止分母为零。Where l = 1, 2, 3 represents the corresponding template number, τ is the transition value in the calculation process, βl is the smoothness indicator, and ε=10-6 prevents the denominator from being zero.
步骤d：通过非线性权和多项式函数p,得到V函数的各阶导数。Step d: Obtain the derivatives of each order of the V function through the nonlinear weight and the polynomial function p.
代入的式子(15)、(16)，可以得到所需的所有虚拟点的点值。Substituting into equations (15) and (16), we can obtain the point values of all the required virtual points.
步骤2.将双曲守恒律方程空间部分用有限差分多重分辨加权基本无振荡格式进行离散；包括：Step 2. Discretize the spatial part of the hyperbolic conservation law equation using a finite difference multi-resolution weighted essentially non-oscillatory format; including:
一维情形时将方程写成守恒形式In one dimension, the equation is written in conservation form
其中，和分别表示通量f(U)在目标单元Ii的边界和处的五阶近似的数值通量，Ui(t)表示U在网格单元Ii内点xi处的值U(xi,t)。in, and They represent the flux f(U) at the boundary of target unit Ii and The numerical flux of the fifth-order approximation at , Ui (t) represents the value U(xi ,t) of U at point xi in grid cell Ii .
求通量f(U)在目标单元Ii的边界和处的五阶近似值和具体步骤如下：Find the flux f(U) at the boundary of the target unit Ii and The fifth-order approximation at and The specific steps are as follows:
步骤2.1用Lax-Friedrichs分裂把通量分裂为其中为了简单起见，本发明只描述f+(U)在点处的重构过程并将其定义为Step 2.1 Use Lax-Friedrichs splitting to split the flux into in For simplicity, the present invention only describes f+ (U) at point The reconstruction process at and defines it as
步骤2.2将目标单元Ii构造大模板T3＝[Ii-2,Ii-1,Ii,Ii+1,Ii+2]，小模板T1＝[Ii]和小模板T2＝[Ii-1,Ii,Ii+1]，其中Ii为对应序号的网格单元。Step 2.2 constructs the target unitIi into a large templateT3 = [Ii-2 , Ii-1 ,Ii , Ii+1 , Ii+2 ], a small templateT1 = [Ii ] and a small templateT2 = [Ii-1 ,Ii ,Ii+1 ], whereIi is the grid unit of the corresponding sequence.
步骤2.3在每个模板上分别重构代数多项式q1(x)、q2(x)和q3(x)，使得他们在单元边界有一阶、三阶和五阶精度。Step 2.3 reconstructs the algebraic polynomials q1 (x), q2 (x) and q3 (x) on each template respectively, so that they have first-order, third-order and fifth-order accuracy at the cell boundary.
步骤2.4任意取线性权为：γ12＝1/11，γ22＝10/11，γ13＝1/111，γ23＝10/111，γ33＝100/111。将多项式q1(x)、q2(x)和q3(x)重组为新的多项式p1(x)，p2(x)和p3(x)，满足p1(x)＝q1(x),Step 2.4 Arbitrarily select linear weights: γ12 = 1/11, γ22 = 10/11, γ13 = 1/111, γ23 = 10/111, γ33 = 100/111. Reorganize the polynomials q1 (x), q2 (x) and q3 (x) into new polynomials p1 (x), p2 (x) and p3 (x), satisfying p1 (x) = q1 (x),
步骤2.5计算光滑指示器βl，用于衡量重构多项式pl(x)在目标单元上的光滑度，计算公式为：Step 2.5: Calculate the smoothness indicator βl , which is used to measure the smoothness of the reconstruction polynomial pl (x) on the target unit. The calculation formula is:
其中l＝2,3表示对应模板序号，表示多项式pl(x)对x的α阶导数，r＝2。但是β1比较特殊，要独立定义：Where l=2,3 represents the corresponding template number, represents the α-order derivative of the polynomial pl (x) with respect to x, r = 2. However, β1 is special and needs to be defined independently:
其中：γ1,1＝1-γ0,1，in: γ1,1 ＝1-γ0,1 ，
步骤2.6通过线性权γl和光滑指示器βl计算非线性权ωl，其计算公式为：Step 2.6 Calculate the nonlinear weight ωl through the linear weight γl and the smooth indicator βl , and the calculation formula is:
其中l＝1,2,3表示对应模板序号，τ为计算过程中的过渡值，βl为光滑指示器，ε＝10-6防止分母为零。Where l = 1, 2, 3 represents the corresponding template number, τ is the transition value in the calculation process, βl is the smoothness indicator, and ε=10-6 prevents the denominator from being zero.
步骤2.7求出数值通量分裂f+(U)在点处的多重分辨WENO重构值：Step 2.7 Find the numerical flux splitting f+ (U) at the point The multi-resolution WENO reconstruction value at:
从而得到空间离散的常微分方程。对于二维和高维情况，逐维离散即可。Thus, we can obtain the spatially discretized ordinary differential equation. For two-dimensional and high-dimensional cases, we can just discretize dimension by dimension.
步骤3.对控制方程中的时间部分用三阶龙格库塔方法和定点快速扫描法离散成全离散的有限差分格式；Step 3. Discretize the time part of the control equation into a fully discrete finite difference format using the third-order Runge-Kutta method and the fixed-point fast scanning method;
所述步骤3为：The step 3 is:
使用三阶TVD Runge-Kutta时间离散公式和定点快速扫描公式将半离散有限差分格式离散成时空全离散有限差分格式，所述时空全离散有限差分格式为关于时间层的迭代公式。The semi-discrete finite difference format is discretized into a space-time fully discrete finite difference format using a third-order TVD Runge-Kutta time discretization formula and a fixed-point fast scanning formula. The space-time fully discrete finite difference format is an iterative formula with respect to the time layer.
所述步骤3，时间离散的三阶TVD Runge-Kutta时间离散公式：In step 3, the time-discrete third-order TVD Runge-Kutta time-discrete formula is:
其中，U(1),U(2)为中间过渡值，Δt为时间步长，上标n表示第n时间层，L(Un),L(U(1)),L(U(2))为高阶空间离散形式的近似值；Among them, U(1) , U(2) are intermediate transition values, Δt is the time step, the superscript n represents the nth time layer, and L(Un ), L(U(1) ), L(U(2) ) are approximate values of the high-order spatial discrete form;
时间离散的定点快速扫描法格式为：The format of the time-discrete fixed-point fast scanning method is:
其中*代表有n+1层的新值就用n+1层的新值，没有新值就用第n层的旧值，扫描的顺序为：i＝1,…,N和i＝N,…,1交替扫描。Among them, * represents that if there are new values in the n+1 layer, the new values in the n+1 layer will be used, and if there are no new values, the old values in the nth layer will be used. The scanning order is: i=1,…,N and i=N,…,1 are scanned alternately.
步骤4.根据时空全离散方法得到下一时间层每一个点的近似值，依次迭代，得到计算区域内守恒变量的残差在趋于稳定时的数值模拟结果。Step 4. Obtain the approximate value of each point in the next time layer according to the full space-time discretization method, iterate in sequence, and obtain the numerical simulation results when the residual of the conserved variables in the calculation area tends to be stable.
所述步骤4为：The step 4 is:
将所述时空全离散有限差分格式初始状态值设为精确解或者近似解，通过迭代公式求出下一时间层的近似值，依次迭代得到残差稳定时计算区域内守恒量的数值模拟值。The initial state value of the space-time fully discrete finite difference format is set to an exact solution or an approximate solution, the approximate value of the next time layer is calculated by an iterative formula, and the numerical simulation value of the conserved quantity in the calculation area is obtained by iterating in sequence when the residual is stable.
这样就得到时空全离散有限差分格式，上述的时空全离散有限差分格式为关于时间层的迭代公式，初始状态值已知，通过迭代公式求出下一时间层的近似值，依次得到残差稳定时计算区域内的数值模拟值。残差ResA的定义如下：In this way, the space-time fully discrete finite difference format is obtained. The above space-time fully discrete finite difference format is an iterative formula for the time layer. The initial state value is known. The approximate value of the next time layer is obtained through the iterative formula, and the numerical simulation value in the calculation area when the residual is stable is obtained in turn. The definition of the residual ResA is as follows:
其中：对于我们的例子欧拉方程，其中in: For our example Euler equation, in
下面给出几个算例作为本发明所公开方法的具体实施算例。Several examples are given below as specific implementation examples of the method disclosed in the present invention.
实施例一.经典的正规激波反射算例。下边界为反射边界，左边界和上边界为狄利克雷边界，右边界为超音速出口。计算区域为[0,4]×[0,1].其数值结果如图1所示。Example 1. Classical regular shock wave reflection calculation example. The lower boundary is the reflection boundary, the left boundary and the upper boundary are Dirichlet boundaries, and the right boundary is the supersonic exit. The calculation area is [0,4]×[0,1]. The numerical results are shown in Figure 1.
实施例二.超音速圆柱绕流算例。计算区域为[-0.5,0]×[-0.5,0.5],圆柱半径0.05，圆心在原点，入流速度为3马赫，其数值结果如图2所示。Example 2. Supersonic flow around a cylinder. The calculation area is [-0.5, 0] × [-0.5, 0.5], the cylinder radius is 0.05, the center is at the origin, and the inflow velocity is Mach 3. The numerical results are shown in FIG2 .
实施例三.超音速NACA0012翼型绕流算例。计算区域为[-1.5,1.5]×[-1.5,1.5]，入流速度为2马赫，冲击角度为1度。其数值结果如图3所示。Example 3. Supersonic NACA0012 airfoil flow calculation example. The calculation area is [-1.5, 1.5] × [-1.5, 1.5], the inflow velocity is Mach 2, and the impact angle is 1 degree. The numerical results are shown in Figure 3.
实施例四.超音速NACA0012翼型绕流算例。计算区域为[-1.0,2.0]×[-1.5,1.5]，入流速度为3马赫，冲击角度为5度。其数值结果如图4所示。Example 4. Supersonic NACA0012 airfoil flow calculation example. The calculation area is [-1.0, 2.0] × [-1.5, 1.5], the inflow velocity is Mach 3, and the impact angle is 5 degrees. The numerical results are shown in Figure 4.
实施例五.亚音速NACA0012翼型绕流算例。计算区域为[-1.0,2.0]×[-1.5,1.5]，入流速度为0.8马赫，冲击角度为1度。其数值结果如图5所示。Example 5. Flow around a subsonic NACA0012 airfoil. The calculation area is [-1.0, 2.0] × [-1.5, 1.5], the inflow velocity is 0.8 Mach, and the impact angle is 1 degree. The numerical results are shown in FIG5 .
实施例六.亚音速NACA0012翼型绕流算例。计算区域为[-1.0,2.0]×[-1.5,1.5]，入流速度为0.85马赫，冲击角度为1.25度。其数值结果如图6所示。Example 6. Flow around a subsonic NACA0012 airfoil. The calculation area is [-1.0, 2.0] × [-1.5, 1.5], the inflow velocity is 0.85 Mach, and the impact angle is 1.25 degrees. The numerical results are shown in FIG6 .
实施例七.亚音速圆柱绕流算例。计算区域为[-5,5]×[-5,5]。圆柱半径为0.5，圆心在原点。入流速度为0.38马赫。其数值结果如图7所示。Example 7. Calculation example of subsonic flow around a cylinder. The calculation area is [-5,5]×[-5,5]. The radius of the cylinder is 0.5, and the center is at the origin. The inflow velocity is 0.38 Mach. The numerical results are shown in Figure 7.
综上所述，本发明定点扫描算法推广到了计算复杂几何区域。ILW边界方法是高阶精度的边界处理方法，利用建立不等距模板的新型的WENO外推使得边界能够达到高阶精度且能有效避免振荡。快速扫描算法扫过物面边界时不需要再进行特殊处理，在网格点上有新点值就用新点值计算，没有新点值就用旧点值计算。通过在每个方向按坐标的拓扑结构顺序，上升和下降的顺序各扫描计算一次，可以极大加快迭代方法的迭代收敛速度。In summary, the fixed-point scanning algorithm of the present invention is extended to the calculation of complex geometric areas. The ILW boundary method is a high-order precision boundary processing method. The new WENO extrapolation by establishing unequal-distance templates enables the boundary to achieve high-order precision and effectively avoid oscillation. The fast scanning algorithm does not require special processing when scanning the boundary of the object surface. If there is a new point value on the grid point, the new point value is used for calculation, and if there is no new point value, the old point value is used for calculation. By scanning and calculating once in each direction according to the topological structure order of the coordinates, in the ascending and descending order, the iterative convergence speed of the iterative method can be greatly accelerated.
本发明提出了一种新型ILW边界处理方法，并提出了一种将它与有限差分多重分辨加权基本无振荡格式结合的快速扫描法，能针对多种复杂计算区域的可压定常流场问题进行高精度数值模拟。本发明给出了该算法的具体的构造过程。相比于经典的ILW物面边界处理方法，通过设计不等距模板，利用多分辨率WENO的思想，使得线性权可以任取和为1的数。该方法不仅能保持物面边界附近计算格式高阶精度不降低的优点，同时能有效减少空间模板的使用数量和减少需要人为调整的参数，线性权也不需要求解可以任取，在求解定常问题时物理量的残差可以快速收敛，实施更加简单且易于实现。对于空间离散，多分辨率WENO离散相比于经典WENO空间离散通过设计不等距模板使得间断处精度可以降的更低，进而使得残差能够下降到机器零。快速扫描算法相比于三阶TVD Runge-Kutta时间离散格式，可以明显提高迭代效率，节约一半左右的迭代CPU时间。该快速扫描法是通过在所有方向上依次扫描，覆盖所有的迎风方向，再结合Gauss-Seidel迭代的特点，有新值用新值没有新值用旧值，可以提高欧拉向前时间离散格式的迭代效率，在保证格式收敛的情况下CFL数可以取较大值。Runge-Kutta时间离散格式由于多计算了两个虚拟时间层，相对计算效率比较低下。此格式的构造简单，数值精度较高，易于推广到多维情形。The present invention proposes a novel ILW boundary processing method, and proposes a fast scanning method combining it with a finite difference multi-resolution weighted basic non-oscillation format, which can perform high-precision numerical simulations for compressible steady-state flow field problems in various complex calculation areas. The present invention provides a specific construction process of the algorithm. Compared with the classical ILW object surface boundary processing method, by designing unequal-distance templates and utilizing the idea of multi-resolution WENO, the linear weight can be arbitrarily selected with a sum of 1. This method can not only maintain the advantage of high-order accuracy of the calculation format near the object surface boundary, but also effectively reduce the number of spatial templates used and the parameters that need to be manually adjusted. The linear weight does not need to be solved and can be arbitrarily selected. When solving steady-state problems, the residual of the physical quantity can converge quickly, and the implementation is simpler and easier to implement. For spatial discretization, multi-resolution WENO discretization can reduce the accuracy of the discontinuity to a lower level by designing unequal-distance templates compared to the classical WENO spatial discretization, thereby reducing the residual to machine zero. Compared with the third-order TVD Runge-Kutta time discretization format, the fast scanning algorithm can significantly improve the iteration efficiency and save about half of the iteration CPU time. The fast scanning method scans in all directions in turn to cover all windward directions. Combined with the characteristics of Gauss-Seidel iteration, the new value is used if there is a new value, and the old value is used if there is no new value. This can improve the iteration efficiency of the Euler forward time discretization format, and the CFL number can take a larger value while ensuring the convergence of the format. The Runge-Kutta time discretization format has a relatively low computational efficiency due to the calculation of two more virtual time layers. This format is simple in structure, has high numerical accuracy, and is easy to generalize to multi-dimensional situations.
以上仅是本发明的优选实施方式，本发明的保护范围并不仅局限于上述实施例，凡属于本发明思路下的技术方案均属于本发明的保护范围。应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理前提下的若干改进和润饰，应视为本发明的保护范围。The above are only preferred embodiments of the present invention. The protection scope of the present invention is not limited to the above embodiments. All technical solutions under the concept of the present invention belong to the protection scope of the present invention. It should be pointed out that for ordinary technicians in this technical field, some improvements and modifications without departing from the principle of the present invention should be regarded as the protection scope of the present invention.
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202011106710.4ACN112307684B (en) | 2020-10-16 | 2020-10-16 | A fixed-point fast scanning method based on multi-resolution WENO format combined with ILW boundary processing |
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202011106710.4ACN112307684B (en) | 2020-10-16 | 2020-10-16 | A fixed-point fast scanning method based on multi-resolution WENO format combined with ILW boundary processing |
| Publication Number | Publication Date |
|---|---|
| CN112307684A CN112307684A (en) | 2021-02-02 |
| CN112307684Btrue CN112307684B (en) | 2024-10-01 |
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202011106710.4AActiveCN112307684B (en) | 2020-10-16 | 2020-10-16 | A fixed-point fast scanning method based on multi-resolution WENO format combined with ILW boundary processing |
| Country | Link |
|---|---|
| CN (1) | CN112307684B (en) |
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN114707254B (en)* | 2022-06-01 | 2022-08-26 | 中国空气动力研究与发展中心计算空气动力研究所 | Two-dimensional boundary layer grid generation method and system based on template construction method |
| CN115329696B (en)* | 2022-10-13 | 2022-12-20 | 中国空气动力研究与发展中心计算空气动力研究所 | Conservation type fixed wall boundary numerical simulation method and device based on non-body-attached grid |
| CN119150754B (en)* | 2024-11-18 | 2025-01-24 | 西南科技大学 | Numerical simulation method and device based on high-precision immersion boundary method |
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| 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 |
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| 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 |
|---|---|---|---|---|
| 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 |
| Title |
|---|
| Inverse Lax–Wendroff Procedure for Numerical Boundary Conditions of Hyperbolic Equations: Survey and New Developments;Sirui Tan;Fields Institute Communications;第1-20页* |
| Publication number | Publication date |
|---|---|
| CN112307684A (en) | 2021-02-02 |
| Publication | Publication Date | Title |
|---|---|---|
| CN112307684B (en) | A fixed-point fast scanning method based on multi-resolution WENO format combined with ILW boundary processing | |
| CN109902376B (en) | A high-precision numerical simulation method of fluid-structure interaction based on continuum mechanics | |
| CN114757070B (en) | A new WENO format construction method for numerical simulation based on trigonometric function framework | |
| Forrester et al. | Optimization using surrogate models and partially converged computational fluid dynamics simulations | |
| CN112100835B (en) | A high-efficiency and high-precision numerical simulation method for airfoil flow around an airfoil suitable for complex flow | |
| CN107220399A (en) | Weight the whole flow field analogy method of non-oscillatory scheme substantially based on Hermite interpolation | |
| Daxini et al. | Parametric shape optimization techniques based on Meshless methods: A review | |
| CN110457806A (en) | Full flow field simulation method based on staggered grid with central fifth-order WENO scheme | |
| Ye et al. | Numerical investigation on the aerothermoelastic deformation of the hypersonic wing | |
| CN108280273A (en) | A kind of limited bulk Flow Field Numerical Calculation method based under non equidistance grid analysis | |
| CN113591020A (en) | Nonlinear system state estimation method based on axial symmetry box space filtering | |
| CN111177965B (en) | A fixed-point fast scan method based on multi-resolution WENO format for solving steady-state problems | |
| CN117922840A (en) | Wing performance test method based on ALW-MR-WENO algorithm | |
| Hazra et al. | Simultaneous pseudo-timestepping for aerodynamic shape optimization problems with state constraints | |
| Chen et al. | Algorithms of isogeometric analysis for MIST-based structural topology optimization in MATLAB | |
| Sheng et al. | A strategy to implement high-order WENO schemes on unstructured grids | |
| US20240233294A1 (en) | Topological preserving deformation method for 3d model based on multiple volumetric harmonic field | |
| Boschitsch et al. | High accuracy computation of fluid-structure interaction in transonic cascades | |
| Mittal et al. | Mixed-Order Meshes through Rp-Adaptivity for Surface Fitting to Implicit Geometries | |
| CN120429962B (en) | Method for propelling disturbance domain of pneumatic simulation of blunt aircraft by overlapping grid shock wave assembly | |
| Felici et al. | Analysis of vortex induced vibration of a thermowell by high fidelity FSI numerical analysis based on RBF structural modes embedding | |
| Aristova et al. | Characteristic scheme for the solution of the transport equation on an unstructured grid with barycentric interpolation | |
| Li et al. | The p-weighted limiter for the discontinuous Galerkin method in solving compressible flows on tetrahedral grids | |
| Ren et al. | On High-Resolution Entropy-Consistent Flux with Slope Limiter for Hyperbolic Conservation Laws | |
| Li et al. | High order compact generalized finite difference methods for solving inviscid compressible flows |
| 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 | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |