王之欢, 赵郑豪, 朱洪强

(南京邮电大学 理学院,江苏 南京 210023)

Runge-Kutta间断Galerkin(Runge-Kutta discontinuous Galerkin,RKDG)方法[1]是目前求解双曲守恒律方程的主流方法之一.它具有精度高,易于处理复杂计算区域和边界条件,并行效率高,易于进行hp自适应等优点.即使初值足够光滑,非线性双曲守恒律方程的解也可能含有间断,导致数值求解时容易出现非物理振荡,甚至导致求解失败.限制器是RKDG方法用来控制数值振荡的一个不可或缺的步骤,它通常由两部分组成:坏单元指示子和数值解重构.坏单元指示子用来探测坏单元,也就是间断附近的单元,然后通过数值解重构控制坏单元处的振荡.既简单又精准的坏单元指示子一直是RKDG方法的重要研究方向,因为坏单元探测过少会导致振荡控制失败,而探测过多会导致额外的数值计算成本.

目前已有大量的坏单元指示子的研究结果.Qiu和Shu[2]进行了一项较早的研究,他们详细地比较了7个常用的坏单元指示子,分别来自基于minmod函数的TVB限制器、BDF限制器、BSB限制器、保单调限制器、改进的保单调限制器、KXRCF激波探测器和子单元探测器.近年来,新的坏单元指示子研究工作不断出现.例如,Vuik和Ryan[3]对局部的指示变量值,利用箱型线进行异常值检测来探测坏单元.Fu和Shu[4]设计了一个做法简单、不含敏感参数的坏单元指示子.Zhu[5]等人将该工作推广至h自适应网格.Ray和Hesthaven[6-7]使用人工神经网络得到了一种新的坏单元指示子.

目前的坏单元指示子仍然面临一些挑战.首先,坏单元指示子的通用性还有待提高,目前没有一个坏单元指示子能在所有情况下都有优异表现.其次,大部分坏单元指示子针对结构网格设计,不能直接应用于非结构网格和其他特殊网格.再次,很多坏单元指示子含依赖于问题的敏感参数,实际应用中非常不便.

坏单元指示子中的指示变量值通常在连续区域和间断区域具有不同的量级.基于这一性质,Zhu等人[8]利用统计分析中常用的K均值聚类[9]算法,将局部区域中的坏单元和好单元分别聚成一类,然后通过简单的分类准则探测出坏单元.该坏单元指示子适用于各种指示变量[10]和网格,且不含对问题敏感的参数,在数值测试中表现良好.

K均值聚类虽然是较为简单且有效的聚类方法,但其算法是NP困难的,通常只能收敛到局部最优解.在K均值聚类坏单元指示子[8]中,用于聚类的指示变量数据是一维的,且最终目标是聚为两类,并不十分复杂,有望通过更简单的方法实现.因此,本文提出一种简单分类方法替代K均值聚类算法,实现将数据聚为两类的目标,从而在保持原方法优点的基础上,进一步简化算法,提高坏单元指示子的探测速度.

1 RKDG方法回顾

我们对双曲守恒律方程的初值问题简要地描述RKDG方法:

(1)

(2)

2 新型坏单元指示子

按照文献[8]的做法,新型坏单元指示子仍然由2个相互独立的部分组成.第一个部分是将所有的单元划分成J个模板Sj,j=1,2,…,J.这些模板必须满足2个条件:一是每个模板中至少包含一个好单元;二是模板尺寸越小越好.这使得任一模板要么仅包含好单元,要么既有好单元又有坏单元.第二个部分是检测出每个模板中的坏单元.下面给出详细的算法.

(3)

若满足

(4)

则模板Sj中无坏单元,检测结束,否则继续进行后面的操作.这一过滤条件是后面模板判断条件(6)的充分条件,可以用来加速判断过程.

(5)

(6)

则模板Sj中无坏单元,否则第Ⅰ类全为好单元,第Ⅱ类全为坏单元.将这一部分的做法总结为如下算法.

算法1给定坏单元指示变量I,对任意一个模板Sj(1≤j≤J)进行下列操作.

步骤3 若(4)式成立,则Sj中没有坏单元,算法结束.否则进行下列步骤:

●若(6)式成立,则Sj中没有坏单元;否则标记第Ⅰ类中的单元为好单元,第Ⅱ类中的单元为坏单元.

与K均值聚类坏单元指示子相比,这个新的坏单元指示子仅使用形式简单、计算速度极快的(5)式就实现了模板单元分类,避免了使用算法复杂、计算耗时的K均值聚类算法,从而大幅节约计算时间.

该坏单元指示子适用于各种指示变量,本文采用文献[10]中表现较好的单元边界跳量指示变量.假设目标单元为K0,它的s个直接邻居单元分别是K1,K2,…,Ks.pj(j=0,1,…,s)是单元Kj上的近似解,则单元边界跳量指示变量定义为

本文对一维问题采用均匀网格,二维问题采用均匀矩形网格.按照文献[8]的做法,取一维模板尺寸为5,二维模板尺寸为4×4.根据文献[10],单元边界跳量指示变量对应的参数为ε=0.001,β=2.0,对不同问题不敏感.对于本文新引入的简单分类准则(5)中的参数θ,我们使用一个算例进行反复试验得到其取值,然后将该取值用于其他算例.结果发现当θ=0.2时,本文所有算例均给出了满意的数值结果.

3 数值测试结果

我们以空气动力学中的欧拉方程组为模型方程,采用经典的数值算例测试新型坏单元指示子.方程的二维形式为

式中:ρ是密度,(u,v)是速度矢量,E是总能量,p是压力且满足p=(γ-1)[E-0.5ρ(u2+v2)],γ=1.4.

算例1Lax问题.考虑初始条件为

计算区域为[-5,5].我们取网格尺寸(单元长度)Δx=1/40,t=1.3时的密度数值解见图1,其中实线表示参考解,正方形代表新型坏单元指示子得到的数值解.图1也给出了新型指示子的坏单元历史图,并用小正方形表示.从图1可以看到,使用新型坏单元指示子的数值解没有出现振荡,且右边的激波能被很好地检测到.

a k=1时的密度解

b k=2时的密度解

c k=1时的坏单元历史图

d k=2时的坏单元历史图图1 Lax问题在t=1.3时的数值结果

算例2爆炸波问题.该测试问题涉及爆炸波的相互作用,其初始条件是

图2为t=0.038,Δx=1/800时的密度数值解和坏单元历史图.从图2可以观察到,新型坏单元指示子能精确捕捉到两个大的间断,且数值解无振荡.

a k=1时的密度解

b k=2时的密度解

c k=1时的坏单元历史图

d k=2时的坏单元历史图图2 爆炸波问题在t=0.038时的数值结果

算例3双马赫反射问题.计算区域为[0,4]×[0,1],底边1/6≤x≤4区域为固壁.初始时,一个与底边交于点(1/6,0),且与固壁形成60°夹角的激波以10马赫的速度向右移动撞击固壁.波前未受扰动的空气密度为1.4,压力为1.底边固壁区域使用反射边界条件,其余区域使用精确波后边界条件.顶部边界使用10马赫激波运动的精确边界条件,左、右边界分别使用入流和出流边界条件.

图3为t=0.2时刻,使用网格尺寸Δx=Δy=1/240得到的密度等高线图、坏单元图以及密度等高线的局部放大图.通过观察图3,我们可以发现探测到的坏单元区域与数值解间断区域保持一致,探测效果良好,且数值解无振荡.

a k=1时的密度等高线图

b k=2时的密度等高线图

c k=1时的坏单元图

d k=2时的坏单元图

e k=1时的密度等高线局部放大图

f k=2时的密度等高线局部放大图图3 双马赫反射问题在t=0.2时的数值结果

算例4前台阶问题.问题开始于长度为3,高度为1的风洞中,风洞底部的台阶高0.2,位于距离风洞左边端口0.6处并一直延伸到风洞右端口.风洞壁使用反射边界条件,左、右边界分别应用入流和出流边界条件.该问题运行到时间为4.0时,台阶角为奇点,会导致顺风方向的底部墙上出现错误的熵层.为避免此问题,我们采用了文献[4]中的处理方法.图4为t=4.0时刻使用网格尺寸Δx=Δy=1/160的密度等高线图以及坏单元图.可以观察到,坏单元探测准确,数值解无振荡.

a k=1时的密度等高线图

b k=2时的密度等高线图

c k=1时的坏单元图

d k=2时的坏单元图图4 前台阶问题在t=4.0时的数值结果

最后,我们将本文坏单元指示子与K均值聚类坏单元指示子进行定量比较.对二维算例的程序进行性能分析,记录坏单元指示子子程序的运行时间和调用次数,两者相除得到每次调用的平均时间,最终计算出本文坏单元指示子相对于K均值聚类坏单元指示子节省时间的百分比(表1).由表1可知,节省时间百分比对双马赫反射问题约为40%,对前台阶问题约为65%,计算时间节省效果显著.

表1 本文坏单元指示子与K均值聚类坏单元指示子运行时间比较

4 结语

本文使用简单分类方法改进了基于K均值聚类算法的坏单元指示子.新型坏单元指示子做法简单,易于数值实现,且不含依赖于问题的敏感参数.对一维和二维经典算例的数值试验表明,新的坏单元指示子能精准捕捉间断,有效控制数值解中的振荡,并显著节省计算时间.