含新型间断探测器的混合WCNS格式在间断无粘可压流的应用

张 昊,邓小刚,2*

(1. 国防科技大学 空天科学学院, 湖南 长沙 410073; 2. 军事科学院 系统工程研究院, 北京 100082)

超声速无粘可压缩流动由双曲守恒律描述,在各种实际问题中得到了广泛应用,相应方程也属于计算流体力学(computational fluid dynamics,CFD)重点求解的模型问题。它的一个主要特征是即使初始条件光滑,解也可能发展出间断[1]。具体来说,当超声速来流流过飞行器外形或受到各种扰动影响后,多样化的流场结构随之产生,例如激波、接触间断、涡与湍流,以及它们的相互作用等。激波等间断会使得流动变量值发生突变,能否对其进行精确还原关系着数值模拟方法是否可靠。在求解过程中,格式不仅需要最小化色散耗散误差达到高分辨率,来精细还原流场中的小尺度涡,还需要在间断附近保持足够的数值耗散,从而抑制非物理振荡,确保计算的稳定性[2]。在这些近乎矛盾的要求之下,如何构造高保真的高阶精度数值格式,同时具备良好的激波捕捉能力,来准确精细地模拟复杂流动特征,成为CFD领域长期的重大挑战。

为了能实现以上要求,高精度格式应运而生。由于内含的数值耗散更低,在同样的网格规模下,高阶格式相比于二阶格式在计算精度上有显著提升[3]。特别地,高阶有限差分方法的数学形式直接且构造方式简单,因而受到了广泛关注。此外对于多维问题,有限差分方法由于可以维数分解,相比于其他方法也更加高效适用[4]。当前发展最多的高阶激波捕捉格式之一就是加权紧致非线性格式(weighted compact nonlinear scheme, WCNS),它由Deng等[5]提出,然后经过一系列改进,是有效结合求解流场光滑与间断区域的方法。根据解的局部光滑程度,WCNS格式可以通过非线性权动态调整子模板在插值中的权重,实现无振荡地捕捉间断。

然而,这种非线性机制也降低了WCNS在光滑区的谱特性。对于无粘可压缩流的小尺度结构,格式的数值耗散可能超过物理耗散[6],使其在光滑区的色散耗散误差经常大于其应恢复到的线性格式的水平[7]。而且在计算含间断的双曲守恒律系统时,需要对全场使用局部特征分解,将问题变换到特征空间求解,以解耦非线性系统从而减少扰动。这虽然可以避免间断周围解的振荡,但也大大增加了计算时间,格式的效率和鲁棒性都有待进一步提高。通常有两种方法可以改善非线性加权格式的色散耗散特性:增强非线性机制、构造混合格式[8]。这里仅对后一类进行讨论。

混合格式由于其易实现性,正成为越来越多学者的研究对象。通过将间断与光滑区二分,它可以更好地利用激波捕捉格式及低耗散高保真算法各自的优势特性,从而相比于完全采用特征分解的单一方法具有更高的分辨率、更优的计算效率[9]。Deng等[10]通过边界变差递减方法将可调节耗散算法用于尺度分辨模拟,构造了激波捕捉的迎风-中心型混合格式。Shen等[11]提出了一类求解双曲守恒律的加权紧致中心格式,能实现任意一致的时空高阶,并借由限制器完成间断捕捉。Wan等[12]提出了一种新的混合策略,将模板分为光滑、不光滑或过渡区域,各自应用相适应的格式来模拟Euler方程的定常解。

不过,在分别处理光滑区和间断之前,必须能够正确分辨它们。间断探测器正是实现此类识别的一种有效手段。目前已经发展出许多种类的问题单元探测器,其性质对整个格式的分辨率和鲁棒性都会产生很大影响,成为发展混合格式的重要挑战[8]。Jiang等[13]以模板上非线性权的最小值与最大值之比作为探测器构造依据,实现了较好的间断捕捉能力。类似地,Tang[14]提出了使用子模板光滑度量的最小值与最大值之比的探测器,进而自适应地调节非线性权中的指数。Ruan等[15]设计了单调保持-加权群速度控制间断传感器,分别使用光滑和间断探测器,将流场分为光滑、大间断、小间断三类区域分别求解。Takagi等[16]基于无参数间断探测标准,将目标本质无振荡格式的非线性权值组合作为探测器,发展了新的激波捕捉框架。Xue等[17]为混合格式发展了一种简化的多层感知光滑探测器,并在光滑与不光滑之外引入高频区,提高了格式的谱分辨率。

本文基于结构网格的高精度有限差分方法,总结在间断无粘可压缩流动计算中的间断探测与混合格式应用问题;给出混合格式的通用两步构造框架,并归纳基于导数组合、光滑度量方法这两种常见构造类型的间断探测器,进而设计出一种新的问题单元探测器,实现对光滑与间断波形的识别,并将其用于发展混合WCNS格式以完成对间断单元的非线性处理;针对波形识别和一维、二维Euler方程算例,从多个角度展示计算结果,并与其他探测器进行对比,对各自的优势与不足进行更深入的分析。

1 控制方程与WCNS离散

将研究限制在五阶WCNS格式,对于其他阶精度可以用同样的方法实现。为了描述算法,简化无粘流动模拟,考虑一维标量双曲守恒律

ut+f(u)x=0

(1)

式中,u=u(x,t)表示流动变量,f(u)是u的通量,x∈[a,b]和t∈[0,∞)分别是空间与时间坐标。在均匀网格单元上,使用有限差分方法进行空间离散。xj=jh表示点j的坐标,其中j=0,1,…,N,网格间距h=(b-a)/N。然后可以得到守恒律方程的半离散形式

(2)

式中,uj(t)和F′j分别代表变量u(xj,t)和通量导数f′j=∂f/∂x|x=xj的数值近似。F′j使用六阶显式差分计算。

(3)

式中,dm(m=0,1,2)被称为最优权或线性权,d0=1/16,d1=5/8,d2=5/16。式(3)的展开式也对应着五阶显式迎风线性插值

(4)

用非线性权ωm代替最优权dm,就得到了非线性插值格式

(5)

通过设计非线性权,在光滑区令ωm接近dm,使目标格式的精度接近线性格式;而在间断区则令ωm趋于0,从而避免跨间断插值,抑制数值振荡。在ωm的选择上,使用耗散较低的Z型权[18]进行计算,从而更加凸显探测器对格式分辨率的提高作用,其形式为

(6)

其中:下标m=0,1,2;q=1为耗散控制参数;全局光滑度量τ=|β0-β2|;小量ε=10-40用来避免分母为零;子模板光滑度量βm的具体形式亦可参考文献[18]。

矢量扩展式(1),有模拟无粘流的Euler方程

(7)

式中,守恒变量U的表达式为

U=(ρ,ρu,E)T

(8)

通量F的表达式为

F=(ρu,ρu2+p,(E+p)u)T

(9)

ρ、u、p分别代表密度、速度、压力;总能E=e+ρu2/2,其中e为内能。为了使方程封闭,再引入理想气体状态方程

p=(γ-1)ρE

(10)

式中,γ为比热比,对于非高温空气一般取γ=1.4。

当使用非均匀网格时,需要从物理坐标系x向计算坐标系ξ做网格变换,该变换关系为

ξ=ξ(x)

(11)

坐标变换后的Euler方程为

(12)

(13)

式(13)中所用到网格导数的表达式为

(14)

其中,1/J为Jacobian矩阵的行列式,即网格导数雅可比

(15)

2 混合WCNS算法与间断探测器

2.1 混合WCNS算法

步骤1：在每个时间步开始时,在整个计算域应用间断探测器来识别不光滑单元。将在半节点xj+1/2处的探测结果记为σj+1/2,对于光滑流动其值为1,而在间断处则等于0。显然,该指标对混合格式的效率和准确性起着重要作用,这将在后文进一步研究。

步骤2：要构造混合WCNS格式,单元j边界上的插值公式需变为

(16)

至此阐明了混合WCNS格式模拟无粘流动的一般框架,可以发现其重点在于能够自动准确分辨间断,从而自适应地采用非线性插值求解,并在光滑节点使用更高效、更高分辨率的线性格式。此外将引入一些经典的间断探测器作为对比,各探测器的表达式将在2.2节一并给出。

需要注意,插值过程在求解方程组时按变量分量逐个进行,而在多维问题中是逐维展开[8-9]。但将探测器应用到多维时,会存在多维耦合效应,给出以下几点原因做定性分析:第一,混合格式所用到的WCNS本身在推广到多维时就考虑了几何守恒律,需要探测器及混合格式与其相匹配;第二,实际流动结构是多维的,流动变量在各方向同时变化,只在一维层面上识别间断可能忽略这种耦合效应,影响结果的准确性;第三,根据文献[19]对依赖域的论述,多维情况下所用到的其他方向节点,虽然处于上游节点的影响域内,但实际上超出了目标点按流动特征所形成的依赖域。此外,在每个空间节点上做多维耦合探测,相比于解耦探测并不会显著增加间断单元识别的计算量,具体细节这里不做讨论。

边界条件方面,由于主要关注内点空间格式构造,且采用笛卡尔网格计算,所以在边界使用虚拟点方法。在超声速入口、出口边界,分别使用自由流条件和内点外推。对于无粘固壁边界,虚拟点的密度、压力等标量与对应内点相同,速度矢量的切向、法向分量则都与对应内点大小相等、方向相反,以保证壁面处速度为0。至于对称边界,将固壁边界中的速度切向分量改为与对应内点方向相同即可。

2.2 新型间断探测器

在介绍新型间断探测器之前,先回顾几种常用的探测器,新方法的提出与它们内含的思想一脉相承。此外每种情况下,在探测前都先对变量做归一化处理,这有助于在光滑区中当解的绝对值很大但相对变化小时消除误判[8]。

Harten[20]基于流场梯度的变化,考虑三点模板上一二阶导数的比值,提出

(17)

(18)

其中:φj+1/2用于反映界面xj+1/2附近的流场光滑程度;阈值φc=0.3;ε为防止分母为0的小量,这里取ε=10-3。

Li等[8]改进了Harten指示器,通过将模板扩展到五点,并重新组合导数,改善了探测器比值项对间断特异性捕捉的能力,大大降低了在极值点附近的误判。具体公式如下:

ψj+1/2=min(ψj,ψj+1)

(19)

(20)

(21)

其中,

(22)

ξ=10-2。从而有间断探测器

(23)

式中,ψc为阈值常数,若ψj+1/2≥ψc则认为流场光滑,否则包含间断。通常ψc越大,混合格式的鲁棒性越好;ψc越小,混合格式的分辨率越高。这里取ψc=0.4。

除了使用相邻点导数组合,还发展出了基于非线性加权格式光滑因子的间断探测器。Fu[9]发现TENO格式通过改进非线性机制,能够在谱空间中有效区分光滑区和间断,并由此提出一种新的基于六点模板S6={xj-2,xj-1,xj,xj+1,xj+2,xj+3}的问题单元探测器。具体形式如下:

(24)

式中,q=6,ε=10-4,βk按照光滑度量公式在四个三点子模板上计算得到。之后对γk进行归一化

(25)

得到探测结果

(26)

(27)

及全场一阶导平均值

(28)

构造组合参数,来反映解的变化情况。

(29)

其中小量ε=10-4用来避免分母为0。进而得到探测结果

(30)

阈值C=3.0。

不难看出,这里新构造的探测器相比Li、Fu探测器,最核心的不同是重点关注子模板一阶导数,并引入一阶导全场平均值作为比较基准,从而敏锐地捕捉间断信息,其效果将在下一节给出。而且在确定阈值C的过程中发现,这种ηk形式对不同流场结构的识别结果差异明显,容易通过设定阈值进行区分,实现对间断的准确识别。

3 数值结果

对上一节提到的间断探测器进行测试,算例使用的所有物理量均为无量纲量。首先用波形来检验探测器的识别能力和适用性,然后将其用于混合格式,求解模拟无粘流的Euler方程对比结果。关于识别变量的选择,波形测试中直接使用波形函数f,Euler方程中使用密度ρ。在所有Euler方程算例中,应用混合Roe-Rusanov通量格式,并采用三阶强稳定性保持Runge-Kutta方法[21]进行时间推进,CFL数均取0.6。

3.1 一维波形测试

首先通过不同波数和幅值的正弦波形来检验各探测器的分辨能力。将计算域x∈[-1,1]用400个均匀网格离散,波形函数如下:

(31)

不同探测器的正弦波形测试结果如图1所示。其中正弦波形用红色实线标出,蓝色点(坐标在右轴)表示探测结果,且调换了σj+1/2的0、1值位置,此时被灰色阴影所覆盖的区域即为识别到的间断,下同。从图中可以看出,四种探测器都识别出了x=0附近的间断,其中新型探测器与Li、Fu探测器的结果相似,除了较高波数时,没有像Harten探测器在变化的幅值和频率下产生误判。

(a) Harten探测器(a) Harten detector

(b) Li探测器(b) Li detector

(c) Fu探测器(c) Fu detector

(d) 新型探测器(d) New detector

然后使用Sod问题[22]结果,进一步考察探测器对不同类型流动结构的识别能力,此时密度波形中从左至右依次分布着膨胀波、接触间断、激波。待识别波形为五阶原始WCNS在2 000个单元上的计算结果,也用2 000个均匀网格离散计算域x∈[-5,5]。

不同探测器的Sod波形测试结果如图2所示。四种探测器都识别出了x=3.5附近的间断,其中Harten探测器的阴影区相比后三者的更细。

(a) Harten探测器(a) Harten detector

(b) Li探测器(b) Li detector

(c) Fu探测器(c) Fu detector

(d) 新型探测器(d) New detector

此外,只有新型探测器判别出了x=2周围的接触间断,表明其在间断识别方面具有更好的鲁棒性。同时,新型探测器也对少部分膨胀波做了判别。

3.2 一维Euler方程测试

下面通过混合WCNS格式计算进一步考察探测器的性能。首先选取Sod激波管问题[22],可以较全面地考察探测器对非定常稀疏波、接触间断、激波等常见流动结构的识别能力。初始条件为

(32)

用N=200个均匀网格离散,推进到无量纲时间t=2.0。图3给出了密度计算结果,其中精确解为五阶WCNS在2 000个单元上的结果,右上方框内为间断区域的放大。几种探测器都能准确识别各种光滑与间断结构,实现混合格式效果。

图3 Sod问题密度计算结果Fig.3 Density results of Sod problem

为了进一步测试混合格式的分辨率与鲁棒性,下个算例使用经典的Shu-Osher问题[23],即一道马赫数为3的激波与带密度扰动的正弦熵波发生非定常相互作用,其也广泛应用于测试格式在捕捉激波的同时对小尺度结构的分辨能力。初始条件为

(33)

使用N=200个均匀网格离散,时间推进到t=1.8。图4给出了流场的密度分布,其中精确解由五阶WCNS在2 000个单元上计算得到,左下方框内为相互作用区的放大。通过对比可知各格式都能实现较好模拟,在间断位置均无过冲。新型探测器具有优于Harten与Fu探测器的分辨率,稍弱于Li探测器的结果,体现出其分辨率与鲁棒性兼顾的识别状态,适合嵌入混合格式。

图4 Shu-Osher问题密度计算结果Fig.4 Density results of Shu-Osher problem

接下来由表1给出应用各探测器的混合格式计算所花费的CPU时间、占原WCNS-Z格式计算时间的比例,以及在最后时刻所判别到的间断单元在全流场所占的百分比。数据结果体现出新型探测器的识别范围适中,且对应混合格式的计算效率较高。

表1 Shu-Osher问题计算时间及比例

3.3 二维Euler方程测试

首先计算双马赫反射问题[24]。该算例描述了在二维无粘流动中,马赫数10的激波以60°角入射到静止平板并继续传播,能考核格式精确捕捉强激波和回流区小尺度结构等的能力。初始条件为

(34)

式中,v表示y方向速度。计算域(x,y)∈[0,4]×[0,1],用400×100个均匀网格离散,并使用消息传递接口(message passing interface,MPI)给定4个核进行并行计算,推进到t=1.8。通过沿滑移线产生的涡结构的丰富程度,可以判断格式的耗散大小。图5给出了四种探测器在回流区附近得到的密度等值线。所有结果都捕捉到了三波点附近的马赫杆、反射激波等,而且Fu和新型探测器识别到了更丰富的滑移线涡结构。

(a) Harten探测器(a) Harten detector

(b) Li探测器(b) Li detector

(c) Fu探测器(c) Fu detector

(d) 新型探测器(d) New detector

在最后的时间步中,四种格式识别到的间断分布如图6所示。可以发现,四种探测器都能识别到马赫反射特有的入射激波、反射激波和马赫杆等结构。其中Li和新型探测器可以准确分辨出滑移线,而且特别判断到了回流区内靠近壁面的卷起涡,Harten的则有较多无效识别。

(a) Harten探测器(a) Harten detector

(b) Li探测器(b) Li detector

(c) Fu探测器(c) Fu detector

(d) 新型探测器(d) New detector

再由表2给出计算时间与比例。新型探测器的间断单元数是除Li探测器之外最少的,与图6结果相对应,而且其计算时间最短,体现出它在识别与求解两阶段的高效。此外,本算例中几种混合格式的计算时间占比都在60%左右。

表2 双马赫反射计算时间及比例

另一个算例是Rayleigh-Taylor(R-T)不稳定[25]问题。最初通过上下放置两种不同密度的流体形成位于y=0.5的接触间断,然后撤掉中间隔板,在重力影响下密度更大的流体会流入另一侧,R-T不稳定就在交界面处产生,并出现小尺度涡结构。由于算例无粘,所以流动细节的多少就可以作为比较格式数值耗散的依据。初始条件为

(35)

其中,c表示声速,计算域为[0,0.25]×[0,1]。在128×512个网格上使用4核MPI并行计算到t=1.95。应用四种探测器的混合格式得到的密度等值线见图7。可以看出,四种混合格式都能得到蘑菇帽状尖峰、剪切层等主要流场结构,且都出现了尖峰两侧下方及沿剪切层的卷起涡,体现出较低的耗散水平,表明探测器对于小尺度结构具有较好的分辨率特性。尤其后三者的结果更相似,且对小尺度结构的分辨更精细,具有更高的可信度。

(a) Harten探测器(a) Harten detector (b) Li探测器 (b) Li detector (c) Fu探测器(c) Fu detector (d) 新型探测器(d) New detector

然后结合图8间断单元分布做进一步对比。Harten探测器对剪切层识别过少,在尖峰下方又大量判定为间断,导致其难以做出准确模拟;后三者在剪切层附近的识别相似且合理。在两侧卷起涡处的问题单元,Li和新型探测器的都较少,且能有效区分涡的边缘和内部,而Fu探测器则将此区域全部作为间断处理,这也再次体现了新型探测器对于混合格式分辨率与鲁棒性的平衡。

(a) Harten探测器(a) Harten detector (b) Li探测器(b) Li detector (c) Fu探测器(c) Fu detector (d) 新型探测器(d) New detector

计算时间与比例由表3展示。由于Harten探测器的结果不够准确,暂不将其列入比较。余下探测器中,Li探测器的间断单元比例最低,然后是新型探测器,Fu探测器的最多,与图8结果相对应;而且新型探测器的计算时间也介于Li和Fu之间,占WCNS-Z格式的时间比例都在60%左右。

表3 R-T不稳定计算时间及比例

为了更多地考察四种探测器的特性,使用对应的混合格式计算二维黎曼问题[26]中的构型12。给定初值

(36)

计算终止时间为t=0.25。这里重点关注间断单元,其在最后时刻的分布见图9。出于鲁棒性的考虑,新型探测器的间断单元占比更高,对剪切层结构的包覆也更完整;与之相对应,Li探测器的间断单元数最少,沿剪切层的包络线最细,甚至在间断面上还出现了一些小缺口;Harten探测器的结果则跟Fu探测器的较为相似,单元比例和分布介于前两者之间。

(a) Harten探测器(a) Harten detector

(b) Li探测器(b) Li detector

(c) Fu探测器(c) Fu detector

(d) 新型探测器(d) New detector

4 结论

本文研究了在同时包含间断与光滑区域的超声速无粘可压流情况下,对激波、接触间断等结构的识别与混合求解问题;提出了一种基于模板导数组合的新型间断探测器,能够自主识别光滑与间断区域,并以此为依据构造了混合WCNS格式,分别采用耗散水平相适应的线性与非线性格式求解,可用于模拟双曲守恒律问题,实现准确高效的激波捕捉。

通过一系列波形识别与一维、二维Euler方程算例,对比了不同探测器的间断识别能力,以及它们对应混合格式的特性。结果表明新型探测器具有不弱于常用间断探测器的识别能力,能兼顾计算的分辨率与鲁棒性,而且在识别和求解阶段的效率较高,更是优于在全场使用特征插值的原始WCNS格式,是混合格式设计中值得推广的方法。后续研究将针对间断探测器更深层机理,并采用更全面的算例展开。