【www.zhangdahai.com--方案模板格式】
武从海 李虎 刘旭亮 罗勇 张树海
(中国空气动力研究与发展中心空气动力学国家重点实验室,四川绵阳 621000)
(中国空气动力研究与发展中心计算空气研究所,四川绵阳 621000)
可压缩流动中存在激波、接触间断等结构.为了捕捉这些流动结构,激波捕捉格式应运而生.早期的激波捕捉格式以总变差减小(total variation diminishing,TVD)格式[1]、无振荡、无自由参数耗散(non-oscillatory and non-free parameter dissipation,NND)格式[2]为主.TVD 方法是当前工业界最常用的激波捕捉方法.经典TVD 方法在极值点附近降为一阶精度,整体最多只能达到二阶精度,不利于获得更细致的流场结构信息.湍流的大涡模拟、直接数值模拟以及气动声学问题的数值模拟中通常需要捕捉这些细微结构,这要求格式具有良好的分辨率性质.如果问题中存在间断,则可采用具有良好分辨率性质的高精度激波捕捉格式进行模拟.
高精度激波捕捉格式最受关注的一个类别是加权本质无振荡 (weighted essentially non-oscillatory,WENO)格式.WENO 格式最初由Liu 等[3]提出,后由Jiang 等[4]改进并达到了模板上的最优精度.后续研究者们针对WENO 格式做了大量的改进研究,如为了获得更好流场分辨率的WENO-M[5],WENOZ[6-7],WENO-NS[8],WENO-CU[9]等,为了获得更好收敛性或稳定性的ZSWENO[10],ESWENO[11]等.加权紧致非线性格式(weighted compact nonlinear scheme,WCNS)[12-13]是另一类借鉴了WENO 思想的激波捕捉方法,张树海[14]对两类格式进行了比较.而关于非结构网格WENO 格式的构造可参看文献[15-16].
WENO 方法采用多个候选模板,每个候选模板上均有一个逼近值,对这些逼近值进行加权得到最终的逼近值.候选模板上函数的光滑程度直接影响这些逼近值的权重.含间断候选模板的权重几乎为0,因此可以达到捕捉间断的效果.在连续区域中,这些权值需尽量接近于最优权值,从而减小数值离散误差.实际计算中,变量在连续区域小尺度的波动使得WENO 格式可能明显偏离线性格式,给数值模拟带来不利影响.
Wu 等[17]提出了光滑因子单频波保常设计准则,即光滑因子满足对单频波为常数.这种光滑因子在一般极值点附近的变化幅度小,从而使WENO 格式更接近线性格式.按照这一要求,Wu 等[17-18]构造了一类光滑因子,并设计了相应的WENO-S 格式.由于光滑因子满足单频波保常准则,因此这类格式的近似色散关系[19]与线性基底格式一致.同时,该类光滑因子在小尺度波动区域变化小,间断附近变化大,因而WENO-S 格式能更有效地区分间断与小尺度波动,在数值模拟中能获得良好的结果.
WENO-S 格式的光滑因子比经典形式的光滑因子[4]更简洁,所需浮点运算少,因而具有较高的计算效率.由于其光滑因子仅与子模板上几个函数值相关,与子模板在大模板中的位置无关.对于线性对流方程,计算相邻数值通量时部分光滑因子相同,无需重复计算.因此,可以通过消除冗余的光滑因子计算来提高计算效率.
本文以7 阶WENO-S 格式为例介绍其构造方法,包括光滑因子的设计准则,然后对格式特点进行分析,接下来讨论计算效率,针对线性对流方程给出基于消除冗余光滑因子计算的算法,然后讨论该方法的使用条件,并将其推广到其他问题的数值模拟中.最后通过数值实验证明该算法的有效性.
1.1 光滑因子单频波保常设计准则
常见WENO 格式对于小尺度波动通常表现出过多的耗散.其原因可归结于极值点附近其权值偏离线性权值较大.如最常用的Jiang-Shu 型r点子模板光滑因子[4]
式中,q(x) 为该子模板上的重构函数,xj为大模板的中点,h为空间步长.在远离临界点的区域,可由泰勒展开得到临界点表示一阶导数为0 的点,极值点则是一类常见的临界点.在临界点附近,则有显然,该类光滑因子在极值点附近下降明显,这使得相应的模板权重变化较大,对数值模拟产生不利影响.
针对上述情况,需要减小连续区域光滑因子的变化幅度.为此,Wu 等[17]提出了光滑因子单频波保常准则.依此准则设计的光滑因子对于单频波为常数,相应WENO 格式非线性权保持不变,退化为线性权.而对于其他一般连续波形,非线性权的变化也明显减小.该文认为,光滑因子应该对单频波为常数的原因如下:
(1)单频波可视为一个单频振荡的函数,每一部分都具有相同的频率和最大振幅;
(2)如图1 所示,几何意义下圆周处处同样光滑,而单频波可作为表示圆周横坐标(或者纵坐标)关于角度的函数,可视为函数意义下处处同样光滑.
图1 圆周和正弦函数的光滑程度示意图Fig.1 Schematic diagram of the smoothness of a circle and a sine function
Wu 等[17-18]给出了一系列满足该准则的光滑因子.以4 点等距模板为例,对应的光滑因子为
1.2 7 阶WENO-S 格式
考虑一维双曲守恒律方程
对于等距网格xj=jh,其半离散守恒型差分格式为
下文中公式均假设特征速度 ∂f/∂u为正.对于特征速度为负的情况,可通过对称性获得相应的结果.若无法确定正负,则需要采用通量分裂方法,如Lax-Friedrichs 分裂,时间格式可以采用高阶Runge-Kutta 方法,相关处理可参见文献[4]及其参考文献.
2.1 收敛精度
7 阶WENO-S 格式在连续区域通常为7 阶精度.但是非线性格式在临界点附近可能发生降阶,这里临界点指一阶导数为0 的点.临界点的阶数用ncp表示,如ncp=1 表示f′=0,f′′≠0,ncp=2表示f′=0,f′′=0,f′′′≠0.根据Wu 等[17]的分析,WENO7-S 在二阶及以下临界点可保持7 阶精度,而在3 阶及以上临界点可能会发生降阶.
对于不小于3 阶的临界点,如果要保持7 阶收敛精度,可在7 阶WENO-S 格式非线性权的计算中采用与空间步长相关的 ε值[20].
2.2 极值点附近的模拟
极值点附近波形的数值模拟是激波捕捉方法的难点.正是由于极值点的存在,经典TVD 格式只能达到二阶精度第1.1 节中光滑因子单频波保常设计准则的出发点之一也是降低极值点附近光滑因子的变化幅度[17].
极值点一般位于连续区域,为获得更好的模拟结果,通常需要WENO 格式的非线性权系数尽量接近于线性权系数.绝大多数极值点属于一阶临界点.由Taylor 展开可知,经典光滑因子 β(JS)[4]在一阶临界点附近和一般连续区域分别为本文中简略起见,光滑因子在不涉及到具体点时省去下标.而7 阶WENO-S 的光滑因子 β(S)则保持为即光滑因子变化较小.那么 β(S)不仅对单频波为常数,而且在一般波形的极值点附近变化也较小.因此,7 阶WENO-S 在极值点附近更接近线性格式,从而数值模拟误差通常更小.
2.3 间断稳定性
为了消除或减小捕捉间断时的振荡现象,要求WENO格式中含间断子模板的权系数尽量接近于0.光滑因子 β(S)在间断区域为O(1),连续区域为而经典光滑因子 β(JS)在两类区域分别为O(1)和那么光滑因子 β(S)在两类区域之间的差距更大.若权系数公式其他部分相同,则7 阶WENO-S格式中含间断子模板的权系数更小.这使得该格式具有良好的间断稳定性.
2.4 分辨率性质
Lele[21]阐述了分辨率性质对差分格式的意义,并采用Fourier 方法分析了紧致差分格式的分辨率特性.但是,Fourier 方法不适用于非线性格式.针对此问题,文献中最常用的方法是近似色散关系(approximate dispersion relation,ADR)分析[18].该方法首先使用数值方法对不同无量纲波数 κ的单频波进行短时间的推进,然后采用离散傅里叶变换结合初值计算得到有效波数.有效波数的实部代表数值解的相速度,通常也用来分析色散性能,而虚部代表耗散率.Li 等[22]和毛枚良等[23]中将ADR 方法进行了简化,去掉了时间推进过程.
由于WENO-S 格式采用了对单频波为常数的光滑因子,该类格式在计算单频波时退化为线性格式,因而其近似色散关系与线性基底格式一致,优于绝大多数同阶WENO 格式.文献[17]中图2 给出了几种7 点WENO 格式及其基底格式的近似色散关系,其中7 阶WENO-S 与线性基底格式完全重合.
Zhao 等[24]和Li 等[25-26]中在ADR 方法的基础上考虑了非线性响应,即单频波在非线性格式下产生的其他波数的杂波.Cravero 等[27]中借用热力学中的“温度”来表示非线性格式的这种性质.由于WENO-S 格式对于单频波退化为线性格式而不产生杂波,因此在这类分析方法中,该类格式不存在非线性响应,具有温度0 的特征.
7 阶WENO-S 格式具有良好的计算效率.其光滑因子 β(S)比经典光滑因子 β(JS)更简洁,所需计算量更小.但是该格式的计算效率还需考虑参数 τ的计算.Wu 等[17]对比了几种WENO 格式单个数值通量的计算量,并通过数值实验说明了其效率优势.但是,对于一些特殊问题,还可进一步提升该格式的计算效率.
3.1 线性对流方程的效率提升方法
采用经典的WENO 格式(3 阶WENO 除外)时,各子模板光滑因子的计算公式形式不同,在相邻数值通量的计算中不会出现相同的光滑因子.而WENO-S 的光滑因子在各子模板上的计算公式除下标不同外形式一致.那么对于线性对流方程,计算相邻数值通量时,部分光滑因子相同.因而可以通过消除WENO-S 的冗余光滑因子计算来提高计算效率.
具体计算中,可将所有光滑因子先行计算并存储.另外,考虑到大模板上的光滑因子 τ(S)的计算,同时需要存储的还有=−fj+3fj+1−3fj+2+fj+3组成的序列.那么,当网格点数较多时,对于7 阶WENO-S格式,子模板光滑因子的计算量约为原来的1/4.
考虑一维线性对流方程
假设对流速度a>0.设计算网格点为 {x1,x2,···,xN}.采用虚拟网格点的方法处理边界条件,对于7 阶WENO 格式,左端需设置4 个虚拟点,右端需设置3 个虚拟点.计算点和虚拟点上的通量组成的序列为 {f−3,f−2,f−1,···,fN+3}.对于a<0 的情况,可通过对称性得到.
简便起见,这里用WENO7-JS 表示经典的7 阶WENO 格式,而WENO7-Z 在WENO7-JS 的基础上采用了Z 型权系数[28],其中幂因子取1,WENO7-S表示7 阶WENO-S 格式.
表1 几种WENO 格式计算N+1 个数值通量所需浮点运算量Table 1 The number of floating point operations required to calculate N+1 numerical fluxes for several WENO schemes
表1 中可以看到,采用该快速算法可以有效降低WENO-S 格式的计算量,N很大时数值通量计算量降低 5/13≈ 38.46%.实际程序运行时,由于还有其他部分的运算,其整体计算效率提高的程度应小于该值.在后文的数值算例测试中,其计算时间约节省20%.
3.2 其他适用情况
上述效率提升方法也适用于其他一些问题.其关键在于一条网格线上用于WENO 重构或插值的量可以表示为一个序列,进行网格线上不同位置的WENO 重构或插值时,选取序列中相应连续几个位置的值,然后进行WENO 重构或插值.
除7 阶W E N O-S 格式外,还有一些其他WENO 格式也可采用这种方法提升计算效率.其要求是光滑因子在各子模板上的计算公式除下标不同外形式一致.如更高阶的WENO-S 格式[18]、WENO-XS格式[31]、WENO-LOC 格式[32]和FWENO 格式[30]等.
在稍复杂流动问题的数值模拟中,一些额外的处理过程可能会破坏这一条件.
(1) 通量分裂
采用守恒型差分方法求解非线性方程时,考虑到迎风性,通常采用通量分裂的方法.如果要采用这种效率提升方法,则必须采用全局通量分裂.而局部特征分裂将破坏上述条件,从而无法采用该方法提升效率.
对于有限体积法和WCNS 型[12]的差分方法,可采用Godunov 类方法计算通量.该类方法中进行非线性插值或重构的量不是通量,而是原始变量或守恒变量,避免了上述全局通量分裂的要求,可采用该加速方法计算.
(2) 特征投影
对于双曲型偏微分方程组,由于信号沿特征方向传播的特点,数值求解采用特征投影处理更符合方程的特点,在含间断问题中结果更准确.对于非线性问题,投影时一般采用局部特征矩阵.
不采用特征投影而直接对各个方程进行数值离散求解的方法称为组元型方法.组元型方法一般仅用于较弱激波的捕捉,对于强激波,该处理可能会产生较强的波后非物理振荡,对计算结果产生明显不利影响.但是局部特征投影破坏了该效率提升方法的前提.那么,这种效率提升方法只能结合组元型求解方法使用.
某些特定问题中部分方程具有明确的特征方向,即使剩下的方程组成的方程组采用特征投影求解,这些方程也可以采用该方法减少计算量.如多组分问题中的组分方程[33]、多介质流问题中表示界面的方程[34]等.
本节先进行极值点附近权系数偏离情况测试,然后采用多个算例对7 阶WENO-S 的性能进行测试验证,并且同时测试所提效率提升方法的有效性.所选算例包括一维对流问题、球面波传播问题、二维旋转问题、二维小扰动传播问题及一维和二维无黏流动问题.空间离散采用WENO7-JS、WENO7-Z和WENO7-S 格式,时间推进采用3 阶TVD Runge-Kutta 方法[4].在进行时间效率对比时,WENO7-JS和WENO-7S 在格式的后面加0 或1 以分别代表经典计算方法和快速计算方法.
4.1 极值点附近权系数偏离程度测试
WENO 格式的权系数在极值点附近常常明显偏离线性权系数,这可能导致数值模拟误差的增大.本算例比较极值点附近几种WENO 格式的权系数偏离程度.下面用表示某点通量的权系数偏离值,统计极值点附近若干通量点的最大偏离值以及平均偏离值.所选取的3 个函数为:(1)f1(x)=,(2)f2(x)=ex−x−1,(3)f3(x)=sin4(πx/4).极值点均取x=0.该点为f1(x)和f2(x)的一阶临界点,f3(x)的3 阶临界点.
计算中取该极值点作为网格点之一,统计该点左右各4 个通量点的情况.网格间距均取h=0.1,其结果放在表2 中.可以看到表中WENO7-S 的权系数偏离值在极值点相对更小.其中一阶临界点附近小一个量级甚至更多,这与第2.2 小节的结论相符.
表2 权系数偏离值的平均值与最大值Table 2 The average values and maximum values of the deviations of weights
对于f3(x)=sin4(πx/4) 的3 阶临界点x=0,几种格式的权系数偏离值均较大,这表明WENO 格式在这种临界点附近与线性基底格式有明显差异,最终导致它们可能在此问题中发生降阶.注意当网格间距进一步减小时,WENO7-JS 由于更大的ε值,其权系数偏离值可能会小于其他两个格式,使得其更容易在密网格下恢复7 阶精度.
4.2 一维对流问题
考虑一组包括高斯波、方波、尖三角波和椭圆波在内的复杂组合包以速度1 向右传播的问题[4],其控制方程为
其中的常数a=0.5,z=−0.7,δ=0.005,α=10,β=.两端采用周期边界条件.
这里考虑波形的远距传输,计算终止时间取T=2000,网格点数取400.为减小时间离散误差的影响,CFL 数取0.01.计算结果显示在图2 中.对于这种远距传输问题,WENO7-JS 耗散过大,而WENO7-Z 则在间断附近出现了明显的数值振荡,但是在三角波峰值附近表现最好.WENO7-S 格式对于其他几种波形均获得了优良的结果.WENO7-JS,WENO7-Z 和WENO7-S 计算结果的L1误差依次为2.51 × 10−1,4.67 × 10−2和4.79 × 10−2.此问题中WENO7-Z 的误差最小,这可能是因为其三角波附近的优良表现和更小的间断抹平程度.
图2 组合波的远距传输问题,网格点数N=400,计算终止时间T=2000Fig.2 Long-distance advection of combined waves with the number of grid points N=400 and the end time T=2000
表3 给出了几种格式的计算时间,结果表明WENO7-S0 比WENO7-JS1 略快,后续算例也大多如此,而表1 中的浮点计算统计则刚好相反.WENO7-JS0 和WENO7-Z 之间也有类似表现.在实际计算中,运行时间与计算平台、浮点计算单元使用率和缓存命中率等均有较大关系[35].本文计算均采用Fortran 程序在CPU 为AMD®Ryzen®R5 4500 U 的计算机上运行.WENO7-S0 的光滑因子具有相同的形式,编译器可能将其处理为向量计算,这会提升其计算速度.WENO7-Z 相比WENO7-JS0 在光滑因子部分的计算时间减少了,而τ及权系数的计算时间增大了.可能在程序运行中,增大的这一部分时间大于减小的时间.
表3 组合波远距传输问题的计算时间Table 3 Computational time of long-distance advection of combined waves
当消除WENO7-S 的冗余光滑因子计算后,其效率得到进一步提高,该问题中WENO7-S1 比WENO7-S0 提高24.64%.
4.3 球面波传播问题
球面波传播问题是一个计算气动声学的标准测试算例[36],其控制方程为
计算区域为[5,450],初始条件为u=0,在r=5处的边界条件为u=sin(ωt).
取 ω=π/4,计算终止时间为T=400,空间步长为 dr=1.为减小时间推进误差的影响,时间步长取dt=0.1.计算结果显示在图3 中.图中几种格式计算结果的波幅有明显差距,WENO7-JS 耗散最大,而WENO7-S 由于其良好的分辨率性质,获得了最优的计算结果.
图3 球面波传播问题,计算终止时间为T=400,空间步长为dr=1.时间步长为dt=0.1Fig.3 The spherical wave propagation problem with the end time T=400,the grid length dr=1,and the time step dt=0.1
由于该问题计算时间过短,计算时间差距过小.因此,表3 给出了时间步长 dt取0.01 时的计算时间统计.结果表明,WENO7-S 计算效率较高,而消除冗余计算后,效率进一步提升21.82%.
表4 球面波传播问题的计算时间,dt=0.01Table 4 Computational time of spherical wave propagation problem with dt=0.01
4.4 二维旋转问题
这个问题来自于Zalesak[37]的文章.该问题中,一个割开的圆柱绕着附近的轴旋转.该运动的控制方程为
这里u=−Ω(y−y0),v=Ω(x−x0),其中 Ω(=2π/360)是旋转角速度.这里计算区域为 (x,y)∈[0,10]×[0,10],(x0,y0)=(5,5)是旋转轴的坐标.初始时刻,割开圆柱的中心坐标为 (xc,yc)=(5,7.5),圆柱半径为1.5,两条割线的横坐标分别为4.75 和5.25,割开截止位置的纵坐标为8.5;割开圆柱上的 ϕ值为3.0,其他区域则是1.0.根据这个控制方程,旋转过程中割开圆柱应该保持原始的形状.本问题将该割开圆柱旋转5 个周期.
该问题中两个分量速度既可能为正又可能为负,通常同时包含正负速度时需要进行通量分裂.但是该问题每条网格线上特征速度可确定为正或负,无需通量分裂处理.
采用200 × 200 均匀网格进行计算,网格点坐标为xi=(i−0.5)·∆x,yj=(j−0.5)·∆y,其中Δx和Δy分别为两个方向的网格间距.图4 给出了计算结果,其中高度代表 ϕ值.图中可以看到,WENO7-JS 对圆柱边缘处的抹平最为严重,顶部区域有塌陷现象.WENO7-Z 的边缘更锐利,但是附近有轻微的过冲现象.WENO7-S 则避免了过度的抹平和边缘处过冲现象.
为了更清晰地显示计算结果的差异,截取直线y=y150=7.475和x=x100=4.975上的结果进行对比.图5(a)为两条直线的位置示意图,图5(b)和图5(c)给出了两条线上的结果.从图5(b) 中可以看到WENO7-JS在割开部分表现出明显的耗散,WENO7-Z 和WENO7-S 则出现 了下冲,其 中WENO7-Z 的下冲程度相对较大,这也与图5(c)中区间[6,8.5]的结果相符.图5(b)中WENO7-Z 在x=6.5附近出现了微弱的上冲.图5(c)中区间[8.5,9]处WENO7-JS 的峰值明显小于精确值3,这也对应图4(b)中的顶部塌陷现象.
图4 二维旋转问题计算结果与精确解Fig.4 The computational results and the exact solution of the twodimensional rotation problem
图5 二维旋转问题两条线上的结果对比图Fig.5 Comparison of results on two lines of the two-dimensional rotation problem
采用
来表示上冲/下冲幅度.WENO7-JS,WENO7-Z和WENO7-S 的 δ值分别为1.158 × 10−4,0.156和5.944 × 10−2.这说明WENO7-Z 更容易出现上冲/下冲现象.而3 种格式的L1误差依次为2.298,1.774 和1.959.说明尽管WENO7-Z 有较明显的上冲/下冲,但是其误差相对最小.从图5(b)和图5(c)中也可看出,WENO7-Z 的间断抹平程度最小,而间断附近贡献了大部分误差量,这使得其误差小于WENO7-S.
为了考察不同网格下的计算情况,采用200 ×200,400 × 400 和800 × 800 这3 套均匀网格对该问题进行模拟,其时间步长依次设为0.1,0.05 和0.025.为减少模拟时间,这里仅将圆柱旋转一个周期.计算结果的L1误差和上冲/下冲幅度由表5 给出.显然WENO7-Z 的L1误差依然最小.而由于间断的存在,3 种格式的L1误差收敛精度均约为1 阶.对于上冲/下冲幅度 δ,200 × 200 和400 × 400 网格时WENO7-JS 最小,而800 × 800 时WENO7-S 最小.随着网格的加密,WENO7-JS 和WENO7-Z 没有表现出明显的变化趋势,而WENO7-S 的 δ值迅速减小,体现了其良好的间断稳定性.
表5 几种格式不同网格下计算结果的L1 误差和上冲/下冲幅度Table 5 The L1 error and up/down overshooting amplitude of the computational results with different schemes and grids
表6 给出了200 × 200 网格旋转5 个周期的计算时间对比.结果表明WENO7-S 计算效率较高,新的计算方法进一步提升了其计算效率,该问题中提升23.48%.
表6 二维旋转问题的计算时间Table 6 Computational time for the two-dimensional rotation problem
4.5 二维小扰动传播问题
控制方程采用二维线化Euler 方程
Mx=0.5,My=0是x和y方向的平均流马赫数.U′表示密度、速度及压强的扰动量.计算区域为(x,y)∈[−100,100]×[−100,100].初值为
该问题为计算气动声学的一个标准算例[36,38].
根据迎风格式的要求,采用了每条网格线上的全局通量分裂.x方向计算时正负通量分别为方向计算时正负通量分别为
计算采用的网格为200 × 200,终止时间为T=30.图6 给出了水平中心线y=0上的密度分布图.可以看到几种WENO 格式均取得了与精确解非常接近的结果.相比而言,WENO7-S 在峰值附近取得了最佳的结果.表7 给出了几种格式的计算时间.本算例中WENO7-S 计算效率最高,且还可受益于本文计算效率提升方法,该问题中提升17.84%.
表7 二维小扰动传播问题的计算时间Table 7 Computational time of two-dimensional small disturbance propagation problem
图6 二维小扰动传播问题计算结果,网格200 × 200,终止时间T=30Fig.6 Computational results of two-dimensional small disturbance propagation problem,with grid 200 × 200 and end time T=30
4.6 一维无黏流动问题
控制方程为一维Euler 方程
式中 ρ,u,p,E表示密度、速度、压强和内能.气体状态方程为
其中比热比 γ=1.4.
这里计算Shu-Osher 问题,该问题中马赫数3 的激波与正弦波相互干扰,诱发高频振动,其初值条件为
计算终止时间T=1.8,这里取网格点为200.
分别采用组元型和特征型方法结合几种几种WENO 格式进行计算.采用5 阶WENO-JS 在4000 个均匀网格点的计算结果作为参考解,图7(a)为参考解的密度分布图.图7(b)对比了几种WENO 格式的计算结果.图例中几种WENO 格式省略了名称前缀WENO7,并用后缀Comp 表示组元型方法,Char 表示特征型方法.
图7 Shu-Osher 问题计算结果的密度对比Fig.7 The density distributions of computational results of Shu-Osher problem
从图7(b)中可以看出,在高频震荡区域,特征型的结果更好,这在WENO7-JS 格式的结果中体现得更为明显.在图中高频震荡左侧,WENO7-Z 和WENO7-S 的组元型结果相对更接近参考解.几种WENO 格式中,WENO7-S 的结果最好,其中特征型的结果比组元型的结果略好.
该问题计算时间较短,为此将网格数设为2000,此时特征型WENO7-JS 的计算时间为5.317 s.对于组元型计算方法,表8 给出了几种方法的计算时间.WENO7-JS 特征型的计算时间约为组元型的1.6 倍.组元型方法中WENO7-S 格式计算时间最短,其中本文加速方法提升了其25.26%的计算速度.
表8 Shu-Osher 问题计算时间,N=2000Table 8 Computational time of Shu-Osher problem,N=2000
该问题中激波马赫数为3,并非弱激波.此时组元型和特征型的结果之间的差距可能相对较明显,如WENO7-JS 格式.但是用WENO7-S 格式计算时两种方法的差距较小,可考虑采用组元型方法并结合本文加速方法.
4.7 二维无黏流动问题
控制方程为守恒形式的二维Euler 方程
式 中 ρ,u,v,p,E表示密 度、x向速度、y向速度、压强和内能.气体状态方程为
其中比热比 γ=1.4.
这里采用二维激波旋涡相互作用问题作为测试算例[4],该问题描述了一个静止的激波和旋涡的相互作用.计算区域为 [ 0,2]×[0,1].一个马赫数为1.1 的静止激波放置在x=0.5 处且与x轴垂直.左状态为 (ρ,u,v,p)=右状态由Rankine-Hugoniot 关系给出.一个中心在 (xc,yc)=(0.25,0.5)的小旋涡放在激波左侧.这个旋涡可被视为速度 (u,v)、熵温度T=p/ρ的小扰动,即
由于特征投影方法采用了流场当地的投影矩阵,WENO7-S1 的效率提升方法无法使用.这里采用组元型方法,该方法计算量比特征投影方法小,但是不宜用于强激波问题.该问题来流马赫数为1.1,说明激波较弱,计算结果对比也表明两种方法差距很小.另外,特征型WENO 格式中特征投影及逆投影所需计算量很大.本算例采用WENO7-JS1 结合特征型方法时计算时间为41.893 s,约为该格式结合组元型方法计算时间13.195 s 的3 倍.这高于一维问题中的1.6 倍.原因是特征投影矩阵在二维问题中是四维,其乘法计算量显著高于一维问题中的三维矩阵.可以预计,在三维问题中,特征型算法的计算时间相比组元型多3 倍以上.
采用组元型方法计算时,x方向正负通量分别为其 中λx为 该x向网格 线上所有矩阵 ∂F/∂U特征值绝对值的最大值.y方向计算时正负通量分别为其中λy为该y向网格线上所有矩阵 ∂G/∂U特征值绝对值的最大值.
图8 中给出了计算结果的压强云图,其等值线范围为1.19～ 1.37,共90 条.其中图8(a)为特征型WENO7-JS 的结果,可以看到在较弱的激波下,组元型和特征型结果相差较小.组元型方法的计算结果中,WENO7-JS 和WENO7-S 较为相近,而WENO7-Z在激波附近聚集了更多等值线.这可能是由于WENO7-Z 在间断附近耗散更小,激波波后非物理振荡更为明显,而这种振荡会导致相应值的等值线多次出现.
图8 激波旋涡相互作用问题,压强等值线范围1.19-1.37,共90 条,网格200 × 100,终止时间T=0.6Fig.8 Shock vortex interaction problem.90 pressure isolines ranging from 1.19 to 1.37.Component-wise seventh-order WENO schemes with grid 200 × 100 and end time T=0.6
表9 给出了几种格式的计算时间.结果表明WENO7-S 格式效率相对较高,WENO7-S1 的效率提升方法进一步减少了21.66%的计算时间.
表9 激波旋涡相互作用问题的计算时间Table 9 Computational time of shock vortex interaction problem
WENO-S 格式具有许多优良的性质.特别是其光滑因子对单频波为常数,使得其是一类满足近似色散关系与线性基底格式一致的激波捕捉格式.另外,该格式能有效区分间断与小尺度波动,在连续区域误差小,间断附近也能保证良好的稳定性,在数值模拟中能获得良好的结果.
由于WENO-S 格式的光滑因子在各子模板上的计算公式除下标不同外形式一致,在计算线性对流方程的相邻数值通量时,部分光滑因子相同.本文基于7 阶WENO-S 格式,介绍了如何通过避免冗余光滑因子计算提高格式计算效率.除WENO-S 格式外,还有一些其他WENO 格式也可采用这种方法提升计算效率.
这种方法可推广至多种问题,其应用条件是一条网格线上用于WENO 重构或插值的量可以表示为一个序列.因此在非线性方程中使用通量分裂时只能采用全局通量分裂.非线性方程的另一类处理方法是对原始变量或守恒变量直接进行WENO 重构或插值,然后采用Godunov 类方法求得数值通量.这时也可以采用避免冗余光滑因子计算的算法.
对于可压缩流动问题,局部特征投影求解会破坏这种消除冗余计算的方法.因此,在可压缩流动模拟中,该方法只能结合组元型方法求解.某些流动问题中部分方程的求解不受该限制,仍可采用该方法减小计算量,如多组分问题中的组分方程、多介质流问题中表示界面的方程等.
数值实验结果表明7 阶WENO-S 格式同时具有良好的小尺度结构分辨率和激波捕捉能力.在计算效率方面,本文所提方法能有效减少7 阶WENOS 格式约20%的计算时间.
猜你喜欢激波通量数值冬小麦田N2O通量研究农业灾害研究(2022年1期)2022-05-07数值大小比较“招招鲜”中学生数理化·高一版(2021年11期)2021-09-05一种基于聚类分析的二维激波模式识别算法航空学报(2020年8期)2020-09-10基于HIFiRE-2超燃发动机内流道的激波边界层干扰分析航空发动机(2020年3期)2020-07-24斜激波入射V形钝前缘溢流口激波干扰研究实验流体力学(2018年3期)2018-10-10适于可压缩多尺度流动的紧致型激波捕捉格式北京航空航天大学学报(2017年8期)2017-12-20基于Fluent的GTAW数值模拟焊接(2016年2期)2016-02-27基于MATLAB在流体力学中的数值分析哈尔滨师范大学自然科学学报(2015年6期)2015-04-23春、夏季长江口及邻近海域溶解甲烷的分布与释放通量中国海洋大学学报（自然科学版）(2014年12期)2014-02-28保护地土壤N2O排放通量特征研究植物营养与肥料学报(2011年4期)2011-10-26本文来源：http://www.zhangdahai.com/qiyewenhua/fanganmobangeshi/2023/0915/654715.html
