Numerical computation of shock wave using wavelet methods
-
摘要: 基于自适应小波配点法和人工黏性技术,构造出一种简单稳定的冲击波数值计算方法。采用小波阈值滤波,生成适应流场分布的多尺度自适应网格,并利用密度场最细尺度的小波系数构造幂函数形式的冲击波定位函数,用以判断冲击波位置。联合人工黏性与冲击波定位函数,自动根据流场梯度严格控制人工黏性的大小和分布。对强/弱冲击波管问题进行计算,结果表明,该方法能够准确捕捉冲击波和有效抑制数值振荡,并且使用简单、分辨率高、计算量小。Abstract: A simple and stable wavelet method, which is based on adaptive wavelet collocation methods and artificial viscosity techniques, was proposed to compute shock waves. Dynamic multiscale grids generated by wavelet threshold filtering adaptive to the flow field were used. The shock waves can be checked out by the shock locator functions with power formula, which are constructed through using the magnitudes of the wavelet coefficients on the finest level in the density fields. Then, the artificial viscous terms including viscosity and shock locator functions strictly control the magnitudes and distributions of the artificial viscosity according to the gradients in the flow field. A strong and a weak shock tubes were computed, which shows that the method can accurately capture shock fronts and effectively restrain numerical oscillations. By the way, it is easy to manipulate, high of resolution and small of computational costs.
-
传统的冲击波捕捉方法为了获得高分辨率、高精度解,需要非常细密的网格,以便能捕捉清晰的冲击波结构。但是,精度和分辨率越高,计算量则越大,特别是高维问题,除非使用非规则(irregular)网格,否则因计算量巨大的原因将很难实现[1]。小波数值方法是基于多分辨分析(multi resolution analysis, MRA)发展的新方法,由于小波函数具有紧支撑特性,因此能够对流场数据进行压缩,生成捕捉流场局部结构的适应性网格,适合描述局部流动特征显著的问题[2]。
小波数值方法主要有两类,即小波-迦辽金(wavelet-Galerkin)[3]和小波配点(wavelet collocation)法[4-6]。小波-迦辽金法不适合处理非线性算子和任意边界条件问题,而自适应小波配点法在这两方面均具有优势,特别是二代小波,在真实物理域中进行变换,可以方便处理任意边界条件,因此发展迅速[7-8]。
Harten[9]最早将小波用于压缩算法,从而减少高精度计算格式所需的网格数;Bürger等[10]、王昱[11]和孙阳等[12]采用小波压缩联合WENO或TVD等格式,各自实现了传统高精度格式的高分辨率计算,表明自适应多分辨率格式为传统的冲击波捕捉格式节约了计算量。
由于小波数值方法基于嵌套的动态自适应网格,而在动态网格上使用迎风格式或其他高精度格式非常复杂[1]。赵勇等[13-14]提出了双重小波收缩法和第二黏性法,避免了在动态网格上使用迎风格式。Regele等[1, 15-17]、Kassoy等[18]和Schneider等[19]根据传统的构造冲击波定位函数和人工黏性的方法,利用最细层小波系数的绝对值和守恒变量构造冲击波定位函数,然后利用冲击波定位函数控制人工黏性的大小和分布,成功求解了冲击波、爆轰波等问题。构造的冲击波定位函数能够准确判断冲击波位置,但由于需要对所有守恒量进行小波分解和计算范数,因此计算稍显复杂,不利于解决复杂的工程问题。
本文给出一种使用简单、计算稳定的冲击波定位函数和控制人工黏性的方法;详细介绍利用小波多尺度分解生成自适应网格、构造冲击波定位函数、控制人工黏性的原理和过程,并对强、弱冲击波管问题进行计算验证,分析计算冲击波的能力和特点。
1. 自适应小波配点法
由多分辨分析理论可知,任意函数f(x)可以多尺度分解为如下形式:
fVJ(x)≈∑k∈Zsj0,kφj0,k(x)+jmax∑j=j0∑k∈Zdj,kψj,k(x) (1) 式中:φj0和sj0分别为尺度函数和尺度函数系数,ψj和dj分别为小波函数和小波系数。
尺度空间V和小波空间W互为正交补空间,即
Vj+1=Vj+Wj ,因此,当小波函数的尺度因子j取至jmax时,尺度空间VJ的尺度因子J=jmax+1。为了获得式(1)中的系数,需要对函数f(x)进行采样,获得(2J+1)个配点的值,然后进行离散小波变换。小波配点法以小波分解为基础,通过式(1)计算偏微分方程中的空间导数项,从而将偏微分方程化为关于时间的常微分方程;然后利用Euler法、R-K法等常微分方程的数值方法步进求解。
冲击波等具有局部特性的问题,小波系数的绝对值仅在变化急剧的局部区域较大,而在连续区域很小。根据小波系数相对于某个阈值ε(ε>0)的大小,将式(1)分为两部分:
f(x)=f⩾(x)+f<(x) (2) 其中
f⩾(x)=∑k∈Zsj0,kφj0,k(x)+J−1∑j=j0∑|dj,k|⩾εdj,kψj,k(x) (3) f<(x)=J−1∑j=j0∑|dj,k|<εdj,kψj,k(x) (4) 对于正则方程,舍去式(4)中的小波,误差的上限满足[7]
‖f(x)−f⩾(x)‖≤C1ε‖f(x)‖ (5) 因此能够保障计算精度。由于小波和配点一一对应,因此舍去小波的同时删除了对应位置的配点,从而生成自适应网格。
生成自适应网格的过程需要注意两个细节:(1) 由于局部流动结构发生变化和移动,因此,根据当前流场生成的网格应具有预测下一时刻流场特性的能力,所以需要适当保留局部特征附近的网格点[7];(2) 为了能够在自适应网格上进行小波变换,需要保留小波变换所需的上一层的网格点。
为了方便处理边界和减少离散小波变换的计算量,本文采用一种二代小波—提升插值小波进行小波变换,预测和更新分别利用两侧各两个点,具体可参考文献[7]。
2. 冲击波捕捉方法
利用小波对流场进行多尺度分解,对于最细尺度,系数最大的小波对应冲击波位置,并且系数越大表明梯度越大。因此通过最细尺度的小波系数构造冲击波定位函数。
2.1 冲击波定位函数
一维守恒方程:
∂U∂t+∂F∂x=∂∂x(υ(Φ)∂U∂x) (6) 式中:U为守恒变量,F为通量,υ为人工黏性,Φ为冲击波定位函数(或通量限制器),等号右侧为人工黏性项。最细尺度小波系数最大的区域能反映冲击波位置,位置误差为最小网格单元。而对于j<jmax,冲击波定位函数不能为判断冲击波位置提供准确信息。因此选择最细尺度jmax对应的小波系数构造冲击波定位函数:[1]
Φk=min(|djmaxk|u,1) (7) 式中:
djmaxk 为最细尺度的小波系数,并利用对应的守恒量的范数对djmaxk,max 归一化。Euler方程组包含多个守恒变量,理论上需要对每个变量进行小波变换,并分别计算冲击波定位函数,并选择最大值作为最终的定位函数值[1]。这种方法适用范围广,但计算复杂,且增加了计算量,并不利于在工程应用中推广。考虑到冲击波问题和一般Riemann问题的共性,选择守恒变量密度ρ进行小波变换、生成适应性网格并构造冲击波定位函数。这有利于简化计算、节约计算量。但对于强冲击波问题,由于人工黏性不足,并不能有效抑制数值振荡。由于式(7)给出的定位函数为不大于1的正数,因此,采用幂函数形式的定位函数能够控制黏性分布的宽度:Φk=min[(|djmaxk|ρ)α,1] (8) 当α<1时,冲击波波阵面两侧的定位函数值增大,数值控制的冲击波区域变宽;当α>1时,波阵面两侧定位函数值减小,区域变窄。
式(7)和式(8)均未给出尺度j<jmax对应配点的Φ值。由于j<jmax和jmax对应的配点满足二分性,因此j<jmax对应配点的Φ值可通过j=jmax对应的Φ值插值给定。对于如图1(a)所示的冲击波,参照阈值ε删除部分配点后,保留的小波系数分布如图1(b)所示,箭头代表小波系数,箭头尾部纵坐标为尺度因子,箭头长短代表小波系数的大小。从图1(b)中可以看出,在冲击波波阵面附近小波系数较大,因此保留了大量配点,而在远离波阵面的区域仅保留了少量的配点。按式(8)计算冲击波定位函数,并取α=1/3,如图1(c)所示,可见集中在波阵面附近的7个点处Φ值较大,远离波阵面处,Φ迅速趋于0。Φ值较大的这7个点构成数值计算的冲击波区,人工黏性将在该区域内起作用。
2.2 计算格式
通量的空间导数选择简单的中心差分格式,二阶差分格式可以计算冲击波,并作为基本计算格式,但其精度较低。当采用四阶或更高阶格式时,可以使接触间断的精度更高。时间积分均采用一阶步进格式。
人工黏性项的离散采用中心差分格式,因此式(6)的二阶计算格式为[20]
U(n+1)−U(n)Δt=−F(n)i+1−F(n)i−12Δx+υi+1/2Ui+1−UiΔx−υi−1/2Ui−Ui−1ΔxΔx (9) 式中:Δx为空间步长,由尺度因子J决定,U为守恒变量。根据非线性稳定性条件确定的最小黏性,式(9)标准的显示格式为
U(n+1)−U(n)Δt=−F(n)i+1−F(n)i−12Δx+12(|ai+1/2|Ui+1−UiΔx−|ai−1/2|Ui−Ui−1Δx) (10) 式中:a为对流速度。对于Euler方程组,为了计算简便,a取Jacobian矩阵的最大特征值c+|u|,其中c为当地声速,c=(γp/ρ)1/2,γ为比热比,p为压力,ρ为密度;u为速度[20]。这种方法的数值稳定性好,但对接触间断有一定抹平[1]。比较式(9)和(10)的人工黏性项,可得人工黏性表达式为
υ=12(c+|u|)Δx (11) 为了保证人工黏性仅在大梯度区起作用,因此将冲击波定位函数Φ与人工黏性υ联合,由Φ控制υ的大小分布,形成对整个计算域通用的人工黏性项:
υ=12Φ(c+|u|)Δx (12) 半点i±1/2的值取相邻整点i和i±1的平均值:
υi±1/2=12(υi+υi±1) (13) 3. 计算结果和分析
为了验证计算冲击波的能力和特点,选择可得精确解的一维冲击波管问题进行计算,控制方程为一维Euler方程,状态方程采用理想气体状态方程,气体比热比γ=1.4。
3.1 弱冲击波问题
冲击波管初始条件设为
{ρ,u,p}={{1,0,1}x∈[−0.5,0]{0.125,0,0.1}x∈(0,0.5] (14) 选取不同的尺度因子进行计算,当J<9时,冲击波抹平严重,J分别取9、10、11或更大值时可以得到较好的数值解。以N0表示有效配点数量,N表示使用的配点数量,具体数值如表1所示。
表 1 弱冲击波管计算参数Table 1. Computational parameters for weak shock tubeJ N α ε/10−4 9 513 1.5 1.0 10 1 025 1.5 1.0 11 2 049 1.5 1.0 当t=0.24时,气体密度分布分别如图2所示。图2表明,分辨率越高,计算的冲击波越接近Riemann解。在相同分辨率条件下,接触间断处的抹平较为明显,精度较低,采用四阶或更高阶差分格式可以在一定程度上提高接触间断处的精度。
基础网格点数量N0=2J+1,J增加1,分辨率和基础网格数量增加一倍,但所用配点数量并未相应增加一倍。J=9(N0=513)时,计算中使用202个配点,当J=10(N0=1 025)时,增加了91个配点,而J=11(N0=2 049)时又增加了81个配点,压缩比(N0/N)分别为2.54、3.50和5.51,说明分辨率越高,数据压缩程度和相对计算效率越高。
图3为不同分辨率条件下小波函数对应的配点分布,尺度函数空间V2的5个配点{−0.5, −0.25, 0, 0.25, 0.5}均不删除,因此未在图中展示。从图3可以看出,随着分辨率的提高,增加的配点集中在冲击波波阵面、接触间断、稀疏波波头和波尾4个梯度较大的位置。
在不同分辨率条件下,冲击波定位函数显示的冲击波位置如图4所示。可以看出,冲击波位置随时间的变化与Riemann解一致。由二代小波分解的特点可知,最细尺度的小波系数判别冲击波位置,产生的最大位置误差为一个网格宽度。因此,分辨率越高,位置误差越小。
3.2 强冲击波问题
初始条件设为
{ρ,u,p}={{1,0,1},x∈[−0.5,0]{0.01,0,0.01},x∈(0,0.5] (15) 当J<10,冲击波抹平严重,J≥10时,可以得到较好的结果,因此分别取J=10、11和12进行计算。对于强冲击波问题,需令参数α<1,使得冲击波区的宽度适当增加,从而能够有效抑制振荡。表2给出了J=10, 11, 12时使用的计算参数,图5对应给出了t=0.05时刻的密度分布。J=10时,计算使用284个配点;当J=11时,分辨率提升一倍,使用的配点增加了58个;而J=12(N0=4 097)时,又增加了60个配点,压缩比(N0/N)分别为3.61、6.00和10.20,分辨率越高,压缩效率越高。
表 2 强冲击波管计算参数Table 2. Computational parameters for strong shock tubeJ N α ε/10−5 10 1 025 0.3 1.0 11 2 049 0.3 1.0 12 4 097 0.3 1.0 3.3 计算效率
为了说明自适应小波配点法在计算效率上的优势,比较几种传统的冲击波捕捉格式在相同分辨率和CFL条件下的运算时间。分别选取基于Steger-Warming矢通量分裂的一阶迎风(Up-wind)格式、三阶ENO格式和五阶WENO格式对3.1节和3.2节中的冲击波管问题进行计算。四种格式的计算时间记为tAWCM、tUp-wind、tENO和tWENO,自适应小波配点法相对于迎风格式、ENO格式和WENO格式的计算时间分别记为τA_U、τA_E和τA_W,定义为
τA_U=tAWCMtUp-wind (15) τA_E=tAWCMtENO (16) τA_W=tAWCMtWENO (17) 尺度因子选择为9~14,对应6种分辨率,计算网格数范围为512~16 384。自适应小波配点法的相对计算时间如图6所示。若CFL条件与自适应小波配点法相同,采用WENO格式时冲击波后出现数值振荡,需要通过减小时间步长消除,取
ΔtWENO=0.75Δtawcm 时,基本无数值振荡,但导致更长的计算时间。相对于一阶迎风格式,自适应小波配点法节省40%~70%的计算时间;相对于ENO和WENO高精度格式,能节约超过70%的计算时间。分辨率提升一倍时,相对计算效率更高。
需要指出的是,本文的自适应配点法没有考虑流动特性,采用了简单的中心差分格式,对于局部间断问题,不利于获得高精度解。提升分辨率后,计算结果的质量低于ENO和WENO格式,特别是在接触间断附近,如图7所示。因此,后续可对考虑流动特性并利用小波构造高精度格式进行研究。
4. 结 论
基于自适应小波配点法和人工黏性技术构造了冲击波数值计算格式;利用小波系数构建了指数形式的冲击波定位函数,并用以控制人工黏性的分布;对强/弱冲击波管问题进行计算,得出如下结论:
(1)利用小波系数构建的指数形式的冲击波定位函数能够准确捕捉冲击波位置,最大位置误差为一个网格宽度;
(2)冲击波定位函数的指数可以控制冲击波区域的宽度,指数越小,冲击波区域越宽,计算越稳定;
(3)自适应配点法利用小波阈值滤波删除大量网格点,比传统方法计算高效,且分辨率越高,相对计算效率越高。
-
表 1 弱冲击波管计算参数
Table 1. Computational parameters for weak shock tube
J N α ε/10−4 9 513 1.5 1.0 10 1 025 1.5 1.0 11 2 049 1.5 1.0 表 2 强冲击波管计算参数
Table 2. Computational parameters for strong shock tube
J N α ε/10−5 10 1 025 0.3 1.0 11 2 049 0.3 1.0 12 4 097 0.3 1.0 -
[1] REGELE J D, VASILYEV O V. An adaptive wavelet-collocation method for shock computations [J]. International Journal of Computational Fluid Dynamics, 2009, 23(7): 503–518. DOI: 10.1080/10618560903117105. [2] CAI W, WANG J Z. Adaptive multiresolution collocation methods for initial-boundary value problems of nonlinear PDEs [J]. SIAM Journal on Numerical Analysis, 1996, 33(3): 937–970. DOI: 10.1137/0733047. [3] 唐玲艳. 双曲型方程数值解的小波方法研究[D]. 长沙: 国防科学技术大学, 2007: 43−84. [4] 李慧敏. 小波分析在双曲型守恒律方程数值解中的应用研究[D]. 南京: 南京航空航天大学, 2002: 6−28.LI H M. Application of wavelet analysis to solving hyperbolic conservation laws [D]. Nanjing: Nanjing University of Aeronautics and Astronautics, 2002: 6−28. [5] COHEN A, KABER S M, MÜLLER S, et al. Fully adaptive multiresolution finite volume schemes for conservation laws [J]. Mathematics of Computation, 2003, 72(241): 183–226. DOI: 10.1090/s0025-5718-01-01391-6. [6] 王一博, 唐建候, 杨慧珠. 自适应小波方法求解波传问题 [J]. 石油地球物理勘探, 2005, 40(6): 628–631. DOI: 10.3321/j.issn:1000-7210.2005.06.004.WANG Y B, TANG J H, YANG H Z. Using adaptive wavelet method to resolve issue of wave propagation [J]. Oil Geophysical Prospecting, 2005, 40(6): 628–631. DOI: 10.3321/j.issn:1000-7210.2005.06.004. [7] VASILYEV O V, BOWMAN C. Second-generation wavelet collocation method for the solution of partial differential equations [J]. Journal of Computational Physics, 2000, 165(2): 660–693. DOI: 10.1006/jcph.2000.6638. [8] VASILYEV O V. Solving multi-dimensional evolution problems with localized structures using second generation wavelets [J]. International Journal of Computational Fluid Dynamics, 2003, 17(2): 151–168. DOI: 10.1080/1061856021000011152. [9] HARTEN A. Adaptive multiresolution schemes for shock computations [J]. Journal of Computational Physics, 1994, 115(2): 319–338. DOI: 10.1006/jcph.1994.1199. [10] BÜRGER R, KOZAKEVICIUS A. Adaptive multiresolution WENO schemes for multi-species kinematic flow models [J]. Journal of Computational Physics, 2007, 224(2): 1190–1222. DOI: 10.1016/j.jcp.2006.11.010. [11] 王昱. 偏微分方程的小波求解法及其在燃烧计算中的初步应用[D]. 长沙: 国防科学技术大学, 2008: 81−128. [12] 孙阳, 吴勃英, 冯国泰. 利用多小波自适应格式求解流体力学方程 [J]. 力学学报, 2008, 40(6): 744–751. DOI: 10.3321/j.issn:0459-1879.2008.06.004.SUN Y, WU B Y, FENG G T. The AUSMPW scheme based on adaptive algorithm of interpolating multiwavelets applied to solve Euler equations [J]. Chinese Journal of Theoretical and Applied Mechanics, 2008, 40(6): 744–751. DOI: 10.3321/j.issn:0459-1879.2008.06.004. [13] 赵勇. 小波数值方法在船舶流体力学中的若干应用[D]. 大连: 大连理工大学, 2011: 29−56.ZHAO Y. Wavelet numerical method with some applications to marine hydrodynamics [D]. Dalian: Dalian University of Technology, 2011: 29−56. [14] 赵勇, 宗智, 王天霖. 一种抑制冲击波计算中数值震荡现象的双重小波数值方法 [J]. 应用数学和力学, 2014, 35(6): 620–629. DOI: 10.3879/j.issn.1000-0887.2014.06.004.ZHAO Y, ZONG Z, WANG T L. A dual wavelet shrinkage procedure for suppressing numerical oscillation in shock wave calculation [J]. Applied Mathematics and Mechanics, 2014, 35(6): 620–629. DOI: 10.3879/j.issn.1000-0887.2014.06.004. [15] REGELE J D, KASSOY D R, VASILYEV O V. Numerical modeling of acoustic timescale detonation initiation: AIAA 2008- 1037 [R]. USA: NASA, 2008. DOI: 10.1080/13647830.2011.647090. [16] REGELE J, RABINOVITCH J, COLONIUS T, et al. Numerical modeling and analysis of early shock wave interactions with a dense particle cloud: AIAA 2012-3161 [R]. USA: NASA, 2012. DOI: 10.2514/6.2012-3161. [17] REGELE J D. Purely gasdynamic multidimentioanal indirect detonation initiation using localized acoustic timescale power deposition: AIAA 2013-1172 [R]. USA: NASA, 2013. DOI: 10.2514/6.2013-1172. [18] KASSOY D, REGELE J, VASILYEV O. Detonation initiation on the microsecond time scale: one and two dimensional results obtained from adaptive wavelet-collocation numerical methods: AIAA 2007-986 [R]. USA: NASA, 2007. DOI: 10.1080/13647830.2014.971058. [19] SCHNEIDER K, VASILYEV O V. Wavelet methods in computational fluid dynamics [J]. Annual Review of Fluid Mechanics, 2010, 42(1): 473–503. DOI: 10.1146/annurev-fluid-121108-145637. [20] LANEY C B. Computational gasdynamics [M]. Cambridge: Cambridge University Press, 1998: 249−299. DOI: 10.1017/cbo9780511605604 期刊类型引用(2)
1. 许志宇,谭永华,李小明. 基于人工黏性的二维激波小波多尺度数值计算. 推进技术. 2021(10): 2229-2236 . 百度学术
2. 许志宇,李小明,谭永华,李永锋,胡攀,董万峰. 热量损失和气相凝结对电爆阀建压影响研究. 火工品. 2020(01): 1-5 . 百度学术
其他类型引用(0)
-