计算力学学报

基于p4est的自适应间断伽辽金方法研究

PDF下载

周星

计算力学学报 | 研究简报 2025,42(5): 865-870

收起

计算力学学报 | 研究简报 2025, 42(5): 865-870

基于p4est的自适应间断伽辽金方法研究

全屏

周星

作者信息

南京航空航天大学　航空学院，南京　210016

周星^*（1991-），男，硕士（E-mail：1748278907@qq.com）.

Discontinuous Galerkin method on adaptive Cartesian grid based on p4set

Xing ZHOU

Affiliations

College of Aerospace, Nanjing University of Aeronautics and Astronautics, Nanjing 210016, China

出版时间: 2025-10-28 doi: 10.7511/jslx20240627002

文章导航

摘要

收起

间断伽辽金（DG）方法因其高精度、易并行等优良特性而得到广泛的应用，网格自适应（AMR）技术已得到广泛采用，与同等级全局细化相比，AMR可以大大降低计算量并提高计算效率。本文结合DG和AMR的优点，基于p4est开源库将一种新型混合限制器应用到自适应笛卡尔网格上的DG方法中，该限制器具有高精度、紧凑性和鲁棒性且其构造非常简单，并通过激波探测器来捕捉激波，本文比较了新型混合限制器和总变分有界（TVB）限制器的性能，发现前者的性能明显优于后者，通过一系列Euler方程和Navier-Stokes方程数值算例验证了本文方法的可行性和高效性，结果表明新型混合限制器在AMRDG方法中表现出良好的性能，具有较小的耗散和很好的激波捕捉能力，在保证高精度的前提下大大提高了计算效率。

关键词

DG / 自适应笛卡尔网格 / 新型混合限制器 / 激波捕捉

Abstract

收起

The discontinuous Galerkin (DG) method has been widely adopted due to its excellent properties such as high accuracy and ease of parallelization. The adaptive mesh refinement (AMR) technique has been widely adopted to improve computational efficiency with much less computational cost compared with uniform global refinement to the same level with AMR. This paper combines the advantages of DG and AMR, and a new hybrid limiter is applied to the DG method on adaptive Cartesian grid based on p4est, an open-source library. The limiter exhibits advantages of high precision, compactness, robustness, and ease of implementation. The shock wave is captured with a shock indictor and the performance of the new hybrid limiter is compared with that of the total variational bounded (TVB) limiter in this paper. The result shows that the performance of the former is significantly better than that of the latter. A series of numerical examples for Euler equations and Navier-Stokes equations are used to verify the feasibility and efficiency of the proposed method. The results show that the new hybrid limiter performs very well in the AMRDG method, it has lower dissipation and great shock capture ability, and the computational efficiency is greatly improved while the accuracy is guaranteed.

Key words

DG / adaptive Cartesian grid / new hybrid limiter / shock capturing

引用本文

周星. 基于p4est的自适应间断伽辽金方法研究. 计算力学学报, 2025 , 42 (5) : 865 -870 . DOI: 10.7511/jslx20240627002

Xing ZHOU. Discontinuous Galerkin method on adaptive Cartesian grid based on p4set[J]. Chinese Journal of Computational Mechanics, 2025 , 42 (5) : 865 -870 . DOI: 10.7511/jslx20240627002

正文

收起

1　引言

收起

间断有限元方法（DG）最早出现在Reed等^[1]求解稳态中子输运方程问题的论文中，提出了在中子输运框架下的第一个DG方法。随后Chavent等^[2]利用线性元离散空间、欧拉（Euler）向前法进行离散时间，将其推广到非线性守恒律的求解中，但其并不稳定，并且其时间步长很受限制。Cockburn等^[3]在空间离散化中使用DG并结合Runge-Kutta时间离散化方法，这一问题得到解决并且这些结果推广到了高阶精度和一维方程组。随后其建立能够轻松解决非线性时间依赖问题，如气体动力学的欧拉方程的框架^[4]。该框架中使用显式非线性稳定的高阶Runge-Kutta格式进行时间离散，使用黎曼求解器作为界面通量并利用总变分有界（TVB）的非线性限制器以实现对强激波的非振荡特性。DG自此得到广泛的应用，尤其在求解流场含有数值间断的问题中。DG方法针对强间断处的处理十分重要，目前仍是DG的难题，常用的方法包括限制器和人工黏性等。Cockburn等^[5]将斜率限制器和TVB限制器应用于二维Euler方程，成功抑制了间断处的数值震荡，但无法实现任意阶精度。Hoteit等^[6]提出了一种可以实现任意阶精度的新型坡度限制器，其思想是基于最小二乘法，通过对解的重构施加一些约束来定义局部最大原理区域，结果表明，其提出的坡度限制器具有较高的精度。Liu等^[7]为双曲守恒定律系统提出了一种本质上无振荡的间断有限元方法（OFDG），其通过引入人工黏性来控制伪振荡，结果表明了人工黏性在DG方法中的可行性和鲁棒性。Wei等^[8]提出了一种用于双曲守恒定律的DG方法中的一种混合限制器，本文称其为HYD型限制器，其结合了高阶DG近似和低阶有限解，保留了本质上的非振荡性质，实现了更好的多尺度结构的分辨率，其构造非常简单，且可以实现任意高阶精度，但其还未应用到网格自适应（AMR）和Navier-Stokes方程，本文将其应用到了AMR和Navier-Stokes方程上。

AMR已得到广泛应用，在精度相同的情况下与全局细化相比，AMR可以以更少的计算成本来提高计算效率。本文使用的p4est^[9]是一个大型开源的h型网格自适应库，其是基于树的网格自适应的一种，其基于MPI实现分布式大规模并行，同时拥有动态负载均衡能力，大大提高了计算效率，消除了开发人员考虑网格自适应复杂算法的麻烦。

本文结合DG方法和AMR技术的优点来提高求解性能。比较了TVB型限制器和HYD型限制器的性能，并实现了HYD型限制器在AMR和Navier-Stokes方程上的应用。

2　数值方法

收起

2.1　控制方程和DG方法

本文考虑Navier-Stokes方程，其控制方程为

(1)

式中U为守恒变量，F_c=（f_c，g_c），F_v=（f_v，g_v）分别为无黏通量和黏性通量，对于Euler方程F_v=0。变量的具体形式为^[10]

(2)

式中ρ为流体密度，u、v为x和y方向速度，p为压强，E=p/（γ－1）+0.5ρ（u²+v²），γ为气体比热比常数，除非特别说明取γ=1.4。

(3)

式中μ为动力黏度，δ_x和δ_y的具体形式为

(4)

式中Pr表示普朗特常数，e=E/ρ。

假设计算域为Ω，该区域剖分为单元K的集合，本文考虑矩形网格，即笛卡尔网格，对于复杂几何，结合浸入边界法便可以实现矩形网格下的复杂几何工况的计算^[11]。为简单起见，考虑标量U，而其他变量的离散过程是相似的。单元K上U的近似解

在由分段多项式组成的有限元空间

表示为

(5)

式中NP为自由度个数，NP=（k+1）（k+2）/2，

为自由度，

（x，y）为单元K上的局部基函数，本文选择勒让德多项式的张量积作为基函数，以P²为例需要的基函数为1，X，Y，X²－1/3，XY，Y²－1/3。其中X=2（x－x_c）/Δx，Y=2（y－y_c）/Δy，Δx、Δy分别为K单元x、y方向的网格尺寸，（x_c，y_c）为K单元的中心坐标，为了方便书写，下文将省略上标K和下标h。将式（1）乘以检验函数，并在K上进行积分，用高斯公式分部积分并用数值通量代替单元边界上的通量值得到DG的弱形式为^[10]

(6)

式中

和

分别为无黏和黏性数值通量，前者选用HLL格式^[12,13]求解，后者选用LDG方法^[14]。式（6）的积分项用高斯积分来逼近，具体可参见文献[5]。对于时间离散，本文采用强稳定性的三阶Runge-Kutta格式^[8]：

(7)

式中L（u）=∂u/∂t，Δt为时间步长，u为基函数系数。该DG全离散格式称为RKDG方法。

2.2　问题单元指示器

本文采用Wei等^[8]提出的Jump指示器来识别问题单元：

(8)

式中U为一所需标量，可以取熵、密度、压强等，本文选用的是密度，

和

为单元边界处内外单元高斯点处的U值，N_d为单元边界面高斯点的总数，单元K的Jump指示器为

(9)

式中h为K的外接圆半径。如果η_i>1，则单元K鉴定为问题单元，C_β为一常数，在一般情况下，设置C_β=1是足够有效的，α=1。

2.3　限制器

本文主要用到两种限制器，即TVB型限制器^[5]和HYD型限制器^[8]。下面将针对标量介绍其具体实现过程，对于方程组则需将变量投影到特征空间后进行限制，具体可参见文献[5]。

2.3.1　TVB型限制器

为了方便叙述，考虑第（i，j）个单元记为I_i_，_j，为了衡量数值解u是否出现震荡，计算出如下四个值：（u_i₊₁/2，_j－u_i_，_j，u_i，_j－u_i_－1/2，_j，u_i_，_j_+1/2－u_i_，_j，u_i，_j－u_i_，_j_－1/2），其中u_i_+1/2，_j为u在单元I_i_，_j右边的左侧极限在边上的积分平均值，其他三个符号类推，u_i_，_j为u的单元平均值，进行如下修正：（u_i_+1/2，_j－u_i_，_j）^new=m（u_i_+1/2，_j－u_i_，_j，u_i₊₁，_j－u_i_，_j，u_i_，_j－u_i_－1，_j），其余变量做类似修正。其中m为修正的minmod函数^[5]

(10)

式中M为常数，h为对应方向的网格尺寸。对于P²空间，若限制器发生作用，则令交叉项系数为0。当高于P²时，本文使用Hoteit等^[6]提出的基于最小二乘法实现任意阶精度的TVB型限制器。

2.3.2　HYD型限制器

首先应用Jump指示器识别问题单元，在任意问题单元Δ_i上，P^k的DG数值解

由低阶限制器限制为一个线性的低阶解

，为了减少这个低阶限制器引入的数值耗散，将这两个解在Δ_i上做如下混合

(11)

式中ω_i=ω（η_i）是高阶多项式的非线性权值，η_i为上文提到的Jump指示器，ω（η_i）计算方法为

(12)

可以看到其构造非常简单，下文的结果会展示其良好的性能。

2.4　自适应方法

2.4.1　AMR指示器

在存在激波的问题中，本文使用Person等^[15]提出的一种基于级数展开的激波探测器和上文提到的Jump指示器来识别需要细化或者粗化的单元。

在没有激波的算例中本文更关注流场中涡流的分布情况。

(13)

式中V为速度矢量，

，K_i为第i个单元的面积，N_c为总的单元数。若τ_ci>1.2σ_c，则该网格标记为需要细化的网格，若τ_ci<0.3σ_c，则该网格标记为需要粗化的网格。

2.4.2　新网格自由度

为了更好描述AMR过程，定义深度，初始网格深度为0，网格细化时会由父单元剖分成四个新的子单元，深度加1，粗化时则由四个子单元合并为一个新的父单元，深度减1，如图1所示。网格自适应后需要计算新单元内的自由度以支持后续的计算，为了保持近似精度和局部守恒性，采用L2投影方法得到了新生成的网格的自由度^[16]：

(14)

式中K^*为新的单元，ϕ为单元内的基函数。对于自适应后出现不对等边的情况，只需将式（7）中线积分项分段积分即可，其余无需改变。

3　数值结果

收起

采用数值模拟来验证本文提出的自适应笛卡尔网格RKDG方法的准确性和有效性。在所有的数值算例中，除了特别说明，均采用P²空间。

3.1　二维黎曼问题

二维黎曼问题包含了强激波、弱激波和接触不连续之间的相互作用^[17,18]，这可能包含了流动结构的各种解，包含了很多微小流场结构，常用于验证数值耗散和精度，其初始条件为（ρ，u，v，p）：

(15)

计算域为[0，1]×[0，1]，所有边界均采用出流边界，计算时间为t=0.8，采用400×400的均匀网格进行数值模拟。图2给出了t=0.8时采用四阶精度（P³），不同限制器计算得到的密度云图，与文献[17，18]的结果进行比较，结果十分吻合。可以看到HYD型限制器耗散明显小于TVB型限制器，有更多流场细节受到捕捉。

3.2　斜激波反射

考虑激波反射问题^[19]，其初始条件为（ρ，u，v，p）=（1.0，2.9，0，0.71429）。计算域为[0，4]×[0，1]，下边界和右边界分别为壁面和出流边界，左边界和上边界为狄利克雷边界，其值为

(16)

初始网格为120×30的均匀网格，并且最大允许加密深度为2，结果如图3所示，与宛一飞^[19]的结果相比本文的激波轮廓更锐利，图4的密度分布也说明了这一点，从网格分布可以看出激波得到准确捕捉。最终网格数为6027，远小于全局加密两次的57600。

3.3　双马赫反射

双马赫反射是二维Euler的一个经典算例^[5]，初始条件（ρ，u，v，p）为

(17)

计算域为[0，4]×[0，1]，左边界为入流边界，右边界为出流边界。对于下边界，对从x=0到x=1/6的部分施加精确的后激波条件，其余部分为壁面边界。对于上边界，对从x=0到

的部分施加后激波条件，其余部分采用前激波条件，计算时间为t=0.2，初始网格为240×60的均匀网格，最大允许加密深度为2，结果如图5所示。同Cockburn等^[5]480×120的结果对比，结果十分吻合，网格在激波面处精准加密。最终网格数为20643，这远小于全局加密两次的230400和Cockburn等^[5]的57600，极大提高了计算效率。

3.4　开尔文-亥姆霍兹不稳定问题

开尔文-亥姆霍兹不稳定性问题是流体动力学中最重要的不稳定性之一^[20]，其初始条件为（ρ，u，v，p）

(18)

式中w=0.1sin（4πx）（exp（－（y－0.25）²/0.0025）+exp（－（y－0.75）²/0.0025），γ=7/5。计算域为[0，1]×[0，1]，所有边界均为周期边界，计算时间为t=0.8，初始网格为80×80的均匀网格，最大允许加密深度为3，结果如图6所示。对比文献[20]，本文结果十分吻合，从网格分布均匀看出，本文准确捕捉到了密度梯度较大的位置。最终网格数为24535，远小于全局加密两次的409600，并且也远小于文献[20]的250000，极大提高了计算效率。

3.5　泰勒格林涡

上述验证的都是Euler方程，下面将验证Navier-Stokes方程，首先考虑二维等熵泰勒格林涡^[21]，其初始条件为（ρ，u，v，p）=（1，sinxcosy，－cosxsiny，1/（γMa²）－sin²y+cos²x），其中Ma为马赫数，取为0.05，雷诺数为10即粘性系数ν=0.1。计算域为[－π，π]×[－π，π]，所有边界均为周期边界，计算时间为t=4，初始网格为50×50的均匀网格，最大允许加密深度为1，结果如图7所示。与文献[21]对比，本文结果十分吻合，精确捕捉到了涡流区域并进行了加密，并且相比文献[21]400×400的均匀网格，本文最终网格数仅为3808，也远小于全局加密一次的10000。

4　结论

收起

本文基于P4est库将新型限制器HYD型限制器应用于AMR中，求解了二维Euler方程和Navier-Stokes方程。通过二维黎曼问题展示了HYD型限制器的高分辨率，并通过斜激波反射、双马赫反射、开尔文-亥姆霍兹不稳定问题、等熵泰勒格林涡等经典算例验证了HYD型限制器在AMR中的可行性和高效性，结果表明其能准确捕捉激波和涡流并对该位置进行网格加密。

参考文献

收起

文献

收起

参考文献引证文献

排序方式：

[1]

Reed

W H

, Hill

T R

. Triangular mesh methods for the neutron transport equation[J]. Los Alamos Report La, 1973.

[2]

Chavent

, Salzano

. A finite-element method for the 1-D water flooding problem with gravity[J]. Journal of Computational Physics, 1982, 45(3): 307-344.

[3]

Cockburn

, Hou

, Shu

C W

. The Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws. IV: The multidimensional case[J]. Mathematics of Computation, 1990, 54(190): 545-581.

[4]

Cockburn

, Shu

H C W

. The Runge-Kutta local projection P1-discontinuous-Galerkin finite element method for scalar conservation laws[A]. 1st National Fluid Dynamics Conference[C]. 1988.

[5]

Cockburn

, Shu

C W

. The Runge-Kutta discontinuous Galerkin method for conservation laws[J]. Journal of Computational Physics, 1998, 141(2): 199-224.

[6]

Hoteit

, Ackerer

, Mosé

, et al. New two-dimensional slope limiters for discontinuous Galerkin methods on arbitrary meshes[J]. International Journal for Numerical Methods in Engineering, 2004, 61(14): 2566-2593.

[7]

Liu

, Lu

J F

, Shu

C W

. An essentially oscillation-free discontinuous Galerkin method for hyperbolic systems[J]. SIAM Journal on Scientific Computing, 2022, 44(1): A230-A259.

[8]

Wei

, Xia

Y H

. An indicator-based hybrid limiter in discontinuous Galerkin methods for hyperbolic conservation laws[J]. Journal of Computational Physics, 2024, 498: 112676.

[9]

Burstedde

, Wilcox

L C

, Ghattas

. P4est: scalable algorithms for parallel adaptive mesh refinement on forests of octrees[J]. SIAM Journal on Scientific Computing, 2011, 33(3): 1103-1133.

[10]

Burbeau

, Sagaut

. Simulation of a viscous compressible flow past a circular cylinder with high-order discontinuous Galerkin methods[J]. Computers & Fluids, 2002, 31(8): 867-889.

[11]

任林娟, 周平章. 基于规则背景网格的数值积分[J]. 计算力学学报, 2025, 42(4): 615-622.

(REN

Linjuan

, ZHOU

Pingzhang

. Numerical integration based on regular background grid[J]. Chinese Journal of Computational Mechanics, 2025, 42(4): 615-622. (in Chinese))

[12]

胡立军, 朱可心, 黎烈龙. 一种求解可压缩流的低耗散中心格式[J]. 计算力学学报, 2025, 42(4): 572-580.

(HU

Lijun

, ZHU

Kexin

, LI

Lielong

. A low dissipative central scheme for compressible flow[J]. Chinese Journal of Computational Mechanics, 2025, 42(4): 572-580. (in Chinese))

[13]

Toro

E F

, Spruce

, Speares

. Restoration of the contact surface in the HLL-Riemann solver[J]. Shock Waves, 1994, 4(1): 25-34.

[14]

Cockburn

, Shu

C W

. The local discontinuous Galerkin method for time-dependent convection-diffusion systems[J]. SIAM Journal on Numerical Analysis, 1998, 35(6): 2440-2463.

[15]

Persson

P O

, Peraire

. Sub-cell shock capturing for discontinuous Galerkin methods[A]. 44^th AIAA Aerospace Sciences Meeting and Exhibit[C]. 2006.

[16]

Zhu

H Q

, Qiu

J X

. An h-adaptive RKDG method with troubled-cell indicator for two-dimensional hyperbolic conservation laws[J]. Advances in Computational Mathematics, 2013, 39(3): 445-463.

[17]

Zhang

, Huang

L T

, Jiang

Z H

, et al. A class of high-order improved fast weighted essentially non-oscillatory schemes for achieving optimal order at any critical points[J]. Physics of Fluids, 2022, 34(12): 126102.

[18]

Ning

J G

, Su

, Xu

X Z

. Improved fifth-order weighted essentially non-oscillatory scheme with low dissipation and high resolution for compressible flows[J]. Physics of Fluids, 2022, 34(5): 056105.

[19]

宛一飞. 双曲守恒律方程的混合WENO方法[D]. 中国科学技术大学, 2023.

(WAN

Yifei

. Hybrid WENO Methods for Hyperbolic Conservation Laws[D]. University of Science and Technology of China, 2023. (in Chinese))

[20]

Jagtap

A D

, Kumar

. Kinetic theory based multi-level adaptive finite difference WENO schemes for compressible Euler equations[J]. Wave Motion, 2020, 98: 102626.

[21]

张旭. 高阶隐式WCNS格式在计算流体力学中的应用[D]. 西南科技大学, 2022.

(ZHANG

. High-order Implicit WCNS Scheme Applied to Computational Fluid Dynamics[D]. Southwest University of Science and Technology, 2022. (in Chinese))

2025年第42卷第5期

PDF下载

引用本文

BibTeX

文章信息

doi: 10.7511/jslx20240627002

接收时间：2024-06-27
首发时间：2026-03-24
出版时间：2025-10-28

补充材料

相关文章

文章信息

作者

出版历史

收稿日期：2024-06-27
修回日期：2024-08-15

基金

作者信息

南京航空航天大学　航空学院，南京　210016

参考文献

分享链接

https://castjournals.cast.org.cn/joweb/jslxxb/CN/10.7511/jslx20240627002

分享至

全文二维码

扫描看全文

引用本文

BibTeX

本文的引用情况

2种不同金属材料的力学参数

科 Family	属数 Number of genus	种数 Number of species	占总种数比例 Percentage of total species (%)	属 Genus	种数 Number of species	占总种数比例 Percentage of total species (%)
鹅膏菌科Amanitaceae	2	11	5.26	鹅膏菌属 Amanita	10	4.78
小菇科 Mycenaceae	2	12	5.74	丝盖伞属 Inocybe	5	2.39
多孔菌科 Polyporaceae	8	14	6.70	蜡蘑属 Laccaria	5	2.39
红菇科 Russulaceae	3	23	11.00	小皮伞属 Marasmius	6	2.87
				小菇属 Mycena	11	5.26
				光柄菇属 Pluteus	5	2.39
				红菇属 Russula	17	8.13
				栓菌属 Trametes	5	2.39

关闭全屏

BibTeX
EndNote
RefWorks
TxT