基于坐标变换的强间断问题伪弧长算法_陈泽平.pdf

第 44 卷第 4 期兵 器 装 备 工 程 学 报2023 年 4月收稿日期:2022 05 23;修回日期:2022 06 13基金项目:国家自然科学基金项目(11822203)作者简介:陈泽平(1995)，男，硕士研究生，E-mail:chernzp_9637 outlook com。通信作者:马天宝(1981)，男，博士，教授，E-mail:madabal bit edu cn。doi:1011809/bqzbgcxb202304011基于坐标变换的强间断问题伪弧长算法陈泽平，王晨涛，李坤，马天宝(北京理工大学 爆炸科学与技术国家重点实验室，北京100081)摘要:针对爆炸冲击波数值模拟中经常要面临的强间断奇异性问题，在一维或多维空间中引入弧长参数，通过添加一个约束方程，致使原网格自适应移动变形而能够在强间断区域分布较多的网格点，以此提高大梯度物理量捕捉的分辨率。针对非均匀变形物理空间中格式构造困难的问题，提出了坐标变换，将物理空间中的控制方程映射至均匀正交的弧长计算空间。基于 MUSCL 格式编写了伪弧长算法相应的一维、二维程序，并对其开展人为解方法验证。将伪弧长算法(PALM)的数值结果同有限体积法进行比较，经过误差分析，验证了伪弧长算法在间断附近具有较高的分辨率。关键词:爆炸冲击波;数值模拟;伪弧长算法;坐标变换;精度验证本文引用格式:陈泽平，王晨涛，李坤，等 基于坐标变换的强间断问题伪弧长算法 J 兵器装备工程学报，2023，44(4):68 76Citation format:CHEN Zeping，WANG Chentao，LI Kun，et al esearch on pseudo arc-length algorithm for strong discon-tinuity based on coordinate transformation J Journal of Ordnance Equipment Engineering，2023，44(4):68 76中图分类号:O381文献标识码:A文章编号:2096 2304(2023)04 0068 09esearch on pseudo arc-length algorithm for strong discontinuitybased on coordinate transformationCHEN Zeping，WANG Chentao，LI Kun，MA Tianbao(State Key Laboratory of Explosive Science and Technology，Beijing Institute of Technology，Beijing 100081，China)Abstract:Aiming at strong discontinuity singularity often occurring in the numerical simulation ofexplosion shock wave，this paper introduces arc-length parameters in one-or multi-dimensional space Byadding a constraint equation，the original grid can move and deform adaptively to achieve distribution of agreat number of grid points in strong discontinuous areas so as to improve the resolution of large-gradientphysical quantity capture Concerning the difficulty in constructing a scheme in non-uniform deformationphysical space，a coordinate transformation is proposed to map the control equation in physical space intothe uniform orthogonal arc-length computational space In addition，the corresponding one-and two-dimensional programs of the pseudo arc-length algorithm are compiled based on MUSCL，and the method ofmanufactured solutions is developed to verify this algorithm The numerical results of the pseudo arc-lengthalgorithm are compared with the finite volume method Through error analysis，it is verified that the pseudoarc-length algorithm has high resolution near the discontinuityKey words:explosion shock wave;numerical simulation;pseudo arc-length algorithm;coordinatetransformation;accuracy verification0引言爆炸与冲击是在高温高压和相变等极端条件下，气液固多介质间强耦合作用的瞬态动力学问题，数值求解该问题模型对航空航天等国防工业及武器装备的研制开发均具有重要的基础应用价值1 2。计算爆炸力学很重要的一个问题在于精确捕捉爆轰波的波阵面及波阵面前后物理量的变化，但对于爆炸冲击波的双曲型守恒律方程，其解随着时间的演化往往会出现奇异性即存在强间断，例如激波、接触间断，或者稀疏波之类的弱间断。为了降低这种奇异性，提高间断分辨率，近年来发展了许多的理论和计算方法，像谱方法、奇异摄动理论、小波分析法3 以及某些高分辨率的数值格式等。在早期，数值计算捕捉激波通常采取一阶精度的差分格式，后来，Van Leer4 将一阶 Godunov 格式进行推广，率先提出了二阶精度 MUSCL(monotone upstream-centred schemes for conservation laws)，随后近 40 年，高精度、高分辨率数值格式蓬勃发展，出现诸如TVD、PPM、ENO、KDG、CE/SE、WENO、WENO-Z 及各类杂交格式等。高精度数值格式一般色散较强，容易产生非物理性振荡，通常需要额外的限制振荡的方法，比如人工粘性法，可人工粘性的添加有时会使数值解的耗散比较严重，间断的分辨率不够清晰，其也非从根本上解决振荡。此外，提升 Eulerian 法数值求解精度的另一个角度则是网格加密。若采用均匀网格计算求解，便需对整个计算域进行加密，使得计算资源极大地浪费，基于这个矛盾，网格自适应技术应运而生。网格自适应技术大致分三类，h-型(额外增加节点)、p-型(增加逼近多项式的阶数)和 r-型(移动网格节点来达到网格重新分布)，此外还有正在发展的结合性方法 h-p 方法及 h-r 方法等，这里着重介绍 r-型，也即移动网格方法。移动网格法发展至今，其所凭借的网格移动策略及网格泛函逐渐多样化，比如基于变量扩展的等势方法、调和映射、坐标变换的雅可比矩阵思想、基于所谓等分布和对齐条件的泛函等。最近几年，Luo 等5 又将 DG(discontinuousgalerkin)方法同移动网格偏微分方程法相组合，提出了一种针对双曲守恒律方程的拟拉格朗日移动网格间断 Galerkin法，从旧网格到新网格的物理变量并不需要插值;Lopez 等6 还提出一种基于 MMPDE 法的并行变分网格改进方案。本文中的伪弧长算法可归结为 r-型方法。由于伪弧长算法涉及网格的自适应移动，进而导致原始均匀正交的物理空间发生扭曲变形，这给格式的重构与插值带来了困难。为了避免直接在变形的非结构网格中重构数值格式(需要大量的模板，计算效率较低)，根据弧长映射关系，借助坐标变换将物理空间映射至均匀正交的弧长计算空间，然后在弧长计算空间中，基于维数分裂思想可以用较少的模板完成高精格式的重构，从而保证了伪弧长算法的计算效率。伪弧长算法能够将物理空间中的强间断问题转变成均分弧长空间中正常的弱奇异流体问题，而计算空间的转换更保证了原本数值格式的运用，在此过程中还能巧妙地避开虚假振荡。算法程序的可靠性一直是数值模拟领域亟待解决的重难点，诸如数理模型的简化、边界条件的设置近似以及迭代格式的选取都有可能对计算结果造成偏差7。人为解方法(method of manufactured solution，MMS)早期经典工作可现于Oberkampf、Trucano 及 oy 等。Blais 等8 提出了一种针对VANS 方程(the volume-averaged Navier-stokes equations)的人为解方法框架以验证其算法程序，并能同任何 CFD 技术相结合。Choudhary 等9 介绍了基于旋度而满足无散度约束的人为解方法，以验证两相不可压缩流控制方程。刘学哲等10 利用其构造的一类二维人为解模型针对性地验证辐射流体力学程序，该办法中动量、能量方程包含源项而质量方程无源项。本文中研究了二维情况下基于 MUSCL 的伪弧长算法模型建立过程，针对网格自适应移动造成的物理空间扭曲变形而不易构造高精度格式的难题，基于坐标变换的思想，将变形的物理空间映射至一套正交的弧长计算空间，从而在弧长计算空间中实现控制方程的求解与新网格守恒变量的插值过程，提高了计算效率。编写了该算法相应的一维、二维程序，在人为解方法基础上验证其精度，借助某些经典算例，将PALM 数值结果对比分析于有限体积法。1伪弧长算法本节伪弧长算法包含了 2 个部分。第一部分，借助有限体积法给出物理空间中控制方程的时空离散格式。如前描述，在考虑网格尺度影响的非均匀网格下，格式重构不易，尤其在二维及其以上的情况，因此，第二部分在网格自适应移动之后，采取坐标变换的策略将物理量映射至弧长计算空间中，并在该空间进行格式重构，求解完之后再逆变换映射回原物理空间。1 1离散格式考虑如下双曲守恒系统:wt+xF(w)+yG(w)=S(w)(1)式中:w 为质量、动量、能量组成的守恒变量;F(w)，G(w)，S(w)为 w 的函数。式(1)在网格单元 Ki，j上进行积分，得到:Ki，jwtd+Ki，jF(w)x+G(w)()yd=Ki，jS(w)d(2)式中:d 为面积微元。96陈泽平，等:基于坐标变换的强间断问题伪弧长算法设?wi，j为网格单元的物理量均值，利用 Green-Gauss 定理对式(2)进行变换，有:Ki，jt?wi，j+Ki，j(w)ni，jds=Ki，jS(?wi，j)(3)式中:Ki，j为网格 Ki，j的面积;Ki，j为网格的边界;ds 为网格边界长度微元;ni，j为边界Ki，j的单位外法向量;(w)=(F(w)，G(w)。数值通量借助局部 Lax-Friedrichs 格式，定义为:h(u，v，n)=12(u)n+(v)n c(v u)(4)式中:u 和 v 为守恒变量的内外重构值;n 为网格每一条边的单位外法向量;c=maxu，v(u)n。式(4)应符合守恒性跟相容性:h(u，v，n)=h(u，v，n)，h(u，u，n)=(u)n(5)式(3)写成半离散格式，有:Ki，jt?wi，j=4k=1kKi，jh(wint(k)i，j，wext(k)i，j，nki，j)ds+Ki，jS(?wi，j)=4k=1h(wint(k)i，j，wext(k)i，j，nki，j)kKi，j+Ki，jS(?wi，j)(6)式中:kKi，j为网格 Ki，j的第 k 条边，k=1，2，3，4;kKi，j为第 k 条边的长度;nki，j为边界kKi，j的单位外法向量;wint(k)i，j，wext(k)i，j为网格在边界的内外逼近值，如图 1 所示，且当wint(k)i，j=?wi，j，wext(k)i，j=?wi，j时，上式(6)退化为一阶格式。现对网格单元边界进行二阶重构，具体 MUSCL 形式11 如下:wint(1)i，j=?wi，j+12(?wi，j?wi，j+1，?wi，j1?wi，j)wext(1)i，j=?wi，j1+12(?wi，j1?wi，j，?wi，j2?wi，j1)wint(2)i，j=?wi，j+12(?wi，j?wi1，j，?wi+1，j?wi，j)wext(2)i，j=?wi+1，j+12(?wi+1，j wi，j，?wi+2，j?wi+1，j)wint(3)i，j=?wi，j+12(?wi，j?wi，j1，?wi，j+1?wi，j)wext(3)i，j=?wi，j+1+12(?wi，j+1?wi，j，?wi，j+2?wi，j+1)wint(4)i，j=?wi，j+12(?wi，j?wi+1，j，?wi1，j?wi，j)wext(4)i，j=?wi1，j+12(?wi1，j?wi，j，?wi2，j?wi1，j)(7)式中:为非线性限制器函数。在计算中 取4:(a，b)=(sign(a)+sign(b)|ab|a|+|b|+(8)式中:为小量，可取 =109。对于时间的离散，忽略掉网格索引 i，j，采取如下的三阶TVD-K 格式:?w(1)=?w(z)+tL(?w(z)?w(2)=34?w(z)+14?w(1)+14tL(?w(1)?w(z+1)=13?w(z)+23?w(2)+23tL(?w(2)L(w)=(4k=1kh(wint(k)，wext(k)，nk)ds+|K|S(w)|K|(9)式中:t 为时间步长;?w(z)为上步物理量，?w(z+1)为更新的物理量;L(?w(k)为微分算子，?w(1)、?w(2)为子步插值。图 1网格单元计算示意图Fig 1 Schematic diagram of grid cell calculation1 2网格自适应移动与计算空间的转换以 x=(x，y)和 =(，)分别代表物理空间坐标和弧长空间坐标，(xi 1，j 1，xi 1，j，xi，j，xi，j 1)为网格单元 Ki，j的4 个顶点，从计算空间 c到物理空间 p的一一对应坐标映射为(x，y)=(x(，)，y(，)。多维空间要同时照顾到网格的移动速度、方向和尺寸问题，且期望网格朝物理量梯度大之处进行移动，借助变分原理可知，网格的自适应移动将受泛函 E(x，y)最小值影响12。E(x，y)=12c(TxG1x+TyG2y)dd(10)式中:G1，G2为给定的正定对称矩阵;=(，)T，且有=/，=/，则式(10)对 应 的 Euler-Lagrange 方程13 为:(G1x)+(G1x)=0(G2y)+(G2y)=0(11)它对应于泛函临界点。定义伪弧长控制函数 =1+a1|w|2+a2|w|2，这里采用单一物理量(其可为，p，化学反应度等)及其梯度相结合的模型，其中 a1，a2为可调参数。借助低通滤波器进行光滑打磨，二维情况表达式为:i，j14i，j+18(i+1，j+i，j+1+i1，j+i，j1)+116(i+1，j+1+i+1，j1+i1，j+1+i1，j1)(12)07兵 器 装 备 工 程 学 报http:/bzxb cqut edu cn/具体光滑次数视情况而定，这里不再赘述。令 G=wI，如此，进一步简化式(11)，有:(x)=0(y)=0(13)对式(13)采用 Gauss-Seidel 迭代法进行计算，具体步骤如下:xz+1i，j=i+1/2，jxzi+1，j+i1/2，jxz+1i1，j+i，j+1/2xzi，j+1+i，j1/2xz+1i，j1i+1/2，j+i1/2，j+i，j+1/2+i，j1/2i1/2，j=12(w(z)i，j)+(w(z)i，j1)i+1/2，j=12(w(z)i+1，j)+(w(z)i+1，j1)i，j1/2=12(w(z)i，j)+(w(z)i1，j)i，j+1/2=12(w(z)i，j+1)+(w(z)i1，j)(14)式中:xz+1i，j为下一时刻网格的坐标，其中，N，N代表弧长空间，方向的网格总数，1iN，1jN，且对于边界(i=0，j=0，i=N，j=M)，网格保持不动。由式(14)知网格移动速度表达式:x=x(z+1)i，j x(z)i，j=x(z+1)i，j x(z)i，j t(15)网格移动完之后，需要得到物理量在新网格上的值，这一步称作物理量重映，本文中采用基于通量的重映方法。如图 2 二维新旧网格转化图所示。图 2二维新旧网格转换示意图Fig 2 Schematic diagram of old grid transformingto new grid in 2D设 Dk为经过一次迭代物理空间中网格变化导致边kKi，j所扫过的面积，则有:Dk=124n=1(xnyn+1 xn+1yn)(16)基于前人的工作14，考虑到新网格 xk+1 与旧网格 xk 之间通量守恒，得到如下守恒插值格式:K(z+1)i，j?w(z+1)i，j=K(z)i，j?w(z)i，j+4k=1Fk(wint(k)i，j，wext(k)i，j)(17)式中:K(z+1)i，j代表新网格的面积;w(z+1)i，j代表新网格的物理量;wint(k)i，j，wext(k)i，j的计算方法见式(7)，守恒插值算法的核心在于 wint(k)i，j，wext(k)i，j的精确构造，考虑到变形物理空间中不易构造高精度格式，因此将重构过程转换至均匀正交的弧长计算空间(详见图 2(b)中。为书写方便而不失一般性，此处略掉上下角标，Fk(，)，k=1，2，3，4，为?w 在 Dk上的积分，Fk(，)的近似求解方法有:Fk(m，n)=max Dk，0 n+min Dk，0 m(18)引入弧长坐标空间变换式(1)原始双曲守恒方程，得到:+F()+G()=S(w)F()=1JUUu+xpUv+yp(e+p)U tpG()=1JVVu+xp(e+p)V tp(19)式中:=t，U=t+ux+vy和 V=t+ux+vy为弧长空间的逆变速度，=w/J，J 为空间转换的 Jacobian 行列式，它的表达式为:1J=(x，y，t)(，)=xxxyyy001=xytxyt0011=xy xy(20)式中:x，x，y，y，x，y为度量系数，其中，x和 y为网格的移动速度，表达式见式(15)。Jacobian 矩阵 J 的计算可借助几何守恒条件来进行，若能给出 J 的初值，J 的推进便可采取式(21)来计算:Jt+(Jt)+(Jt)=0t=xy xyJ，t=yx xyJ(21)x的计算过程如式(22)，其余 3 个度量系数的计算相类似:17陈泽平，等:基于坐标变换的强间断问题伪弧长算法x=x+1/2 x1/2，=1，x+1/2=x+1+x2(22)这样通过式(19)求解完物理量之后，再进行逆变换 =w/J，便可得到原始空间物理量。综上所述，可以给出基于二阶 MUSCL 的伪弧长算法具体步骤:第 1 步:初始化，时间层 n=0 时，给出原始网格的坐标分布 x(0)i，j，其中 x(0)i，j=(x(0)，y(0)以及初始物理量 w(0)i，j。第 2 步:求解弧长监控函数，并根据式(12)进行光滑滤波打磨。第 3 步:根据式(14)迭代求解新位置网格坐标 x(z+1)i，j。迭代限制条件为:(x(z)i，j x(z+1)i，j)x(z)i，j，也即相邻时间步网格的距离差小于一个自定义的极小常数。第 4 步:新网格物理量的更新，根据式(17)更新新网格下的物理量，其中 wint(k)i，j，wext(k)i，j守恒变量的计算在均匀正交的弧长计算空间中完成。第 5 步:控制方程的时间步迭代更新，通过式(19)、式(21)及式(22)得到弧长计算坐标系中的控制方程，通过式(9)更新下一时刻物理量，求解完之后再根据逆变换将物理量映射回原物理空间。第 6 步:若求解达到终点时间，程序终止，否则转到第3 步。2精度验证出于方程跟问题的复杂性，精确解很难获取，因此，这里采用人为解方法对一维、二维程序进行验证。2 1一维情况第 1 章的理论基于二维情况，对一维而言需降维处理计算，这里不再赘述。针对一维情况构造一组人为解(算例中物理量为无量纲，下同)，计算域取 0，2，计算终止时间tend=2 0，初值条件为:(x，0)=1 0+0 2sin(x)u(x，0)=1 0，p(x，0)=1 0(23)计算采用周期性边界条件，则其精确解如下:(x，t)=1 0+0 2sin(x t)u(x，t)=1 0，p(x，t)=1 0(24)伪弧长控制函数取(a1，a2)=(10，20 25)，对于一维、二维情况范数误差公式及收敛阶计算式如下15:L1=Ni=1wexai wnumixL1=Ni=1Mj=1wexai，j wnumi，jxyO=ln(Ek1/Ek)ln(sk1/sk)(25)式中:L1误差范数反映的是数值解误差的均值;wexai为精确解;wnumi为数值解;sk代表网格单元长度;Ek代表网格步长为 sk时具体的 L1误差范数。验证结果如表 1 所示。可以看出伪弧长算法并没有因网格变形而损失太多的精度，随着计算网格数量的增多，其L1误差在逐渐减小，计算结果逐渐接近于精确解，并且，基于MUSCL 的伪弧长算法收敛阶也为二阶。将伪弧长算法的误差和收敛阶同有限体积法的进行比较，能看到相同网格数目下，PALM 的 L1误差一般比固定网格的要小。表 1有限体积法与伪弧长算法在不同网格数时的误差和精度(1D)Table 1 Error and accuracy of FVM and PALM(1D)网格数目fixed gridL1Opseudo arc-lengthL1O402717 4E-022705 8E-02807941 3E-031774 87919 4E-031772 61602235 6E-031828 72231 4E-031827 43206268 9E-041834 46203 5E-041846 822二维情况基于两步化学反应流模型，针对式(1)的双曲守恒系统，令:w=(，u，v，E，)TF(w)=(u，u2+p，uv，(E+p)u，u，u)TG(w)=(v，v2+p，uv，(E+p)v，v，v)TS(w)=(0，0，0，0，)(26)式中:，为化学反应变化速率，且有:=ddt=kexpE()T=ddt=0 0 kp22expE()T(1 )2expE+Q()T 0(27)式中:，分别代表未活化的物质占所有物质的比例和放热反应进行的程度;Q 为热释放率;为气体常数;k，k为反应率常数;E，E为激活能。状态方程:E=p 1+12(u2+v2)+Q(28)式中:为比热比。计算域取 0，2 0，2，计算终止时间 tend=2，给定初值条件:27兵 器 装 备 工 程 学 报http:/bzxb cqut edu cn/(x，0)=10+05sin(x+y)u(x，0)=v(x，0)=10，p(x，0)=10(x，0)=05+05sin(x+y)(x，0)=05+05sin(x+y)(29)计算采用周期性边界条件，其精确解如下:(x，t)=10+05sin(x+y t)u(x，t)=v(x，t)=10，p(x，t)=10(x，t)=05+05sin(x+y t)(x，0)=05+05sin(x+y t)(30)伪弧长控制函数取(a1，a2)=(6，16)。验证结果如表2 所示。二维情况下伪弧长算法比之固定网格方法同样能适当提高计算精度并减小 L1误差。须知，PALM 并未增加网格节点，它是根据物理量的分布状况致使网格自适应移动变形，从而在物理量梯度大的地方集聚了较多网格，减小了总体均值误差。随着网格单元数量增多，PALM 的L1误差在逐渐减小，且其收敛速度较快，收敛阶更接近于2。表 2有限体积法与伪弧长算法在不同网格数时的误差和精度(2D)Table 2 Error and accuracy of FVM and PALM(2D)网格数目fixed gridL1Opseudo arc-lengthL1O40 4080 80160 160320 3200324 780120 370034 130009 72143 201818 41812 00305 610109 340029 750007 821482 91877 91927 63数值算例下面基于伪弧长算法开展算例验证。31Sod 激波管问题初值条件:(，u，p)=(1，0，1)x 0(0125，0，01)x 0(31)计算域为5，5，计算终止时间 tmax=2 0，边界条件为自由输出边界条件，取 1 4，伪弧长控制函数为(a1，a2)=(17，200)，网格数 N=150。计算结果如图 3 所示，分别以Palm、Fixed 来代表伪弧长算法和固定网格下有限体积法的数值标识。该算例参照解(eference)在 5 000 个固定网格下获得。对比有限体积法，伪弧长算法能以较少的网格数获得更好的界面分辨率，在相同分辨率的要求下，PALM 显得更加高效;同时也发现，要想提高计算精度，相当数目的网格量是必须的。图3(b)给出了移动网格轨迹线图，表征伪弧长算法很好地实现了对奇异间断处的自适应捕捉，在弧长参数作用下，网格点可聚集在激波、接触间断处甚至是稀疏波的头部与尾部。图3Sod 激波管图Fig3 esults of Sod shock tube32一维爆炸波该双爆轰波碰撞问题初值条件如下:(，u，p)=(1，0，1 000)0 x 01(1，0，001)01 x 09(1，0，100)09 x 1(32)计算域取 0，1，计算终止时间为 tmax=0 038，置 CFL 系数为08，边界条件为反射边界，取 1 4。计算网格数 N=250，伪弧长控制函数以两套不同的参数进行比较，分别为:(2，20)和(2，200)，在图 4 中标识为 Palm-2-1 和 Palm-2-2。由于稀疏波的作用，计算过程中可能会出现负密度、负压力，从而导致程序终止，因此对该算法附加了保正性条件，详细理论见文献 16。此算例参照解在10 000 个固定网格下获得。图 4 将 3 种数值结果同参照解进行比较，说明在网格总数较少情况下伪弧长算法的数值结果显著优于有限体积法，而固定网格算法数值耗散较大，对极值点的捕捉能力很弱。在相同分辨率的要求下，PALM 比固定网格法更为高效。图5 的网格轨迹线刻画出冲击波在 t=0 027 时刻左右相撞，之后又沿 2 个方向继续传播，伪弧长算法对解变化剧烈区域之模拟效果更为逼真。另外能够验证，不同伪弧长监控函数模型对激波的捕捉能力存在差异，针对本算例第二套系数刻画的分辨率的确要优于第一套，其网格移动得更剧烈。37陈泽平，等:基于坐标变换的强间断问题伪弧长算法图4爆炸波问题数值结果比较Fig4 Numerical results of blast-wave problem图5不同参数下网格轨迹图Fig5 Grid trajectories under different parameters33二维 iemann 问题在计算域 0，1 0，1 内，二维 iemann 问题初始分布如图6 所示，其边界条件均为流入流出条件。图6二维 iemann 问题初始区域分布图Fig6 Initial region distribution of two-dimensionaliemann problem对于这样一个 iemann 问题 各有 2 个正负向滑移线的接触间断，初值条件为:(，u，v，p)=(1，075，05，1)ifx 05，y 05(3，075，05，1)ifx 05，y 05(1，075，05，1)ifx 05，y 05(2，075，05，1)ifx 05，y 05(33)计算终止时间为 tmax=03，伪弧长控制函数取(a1，a2)=(27，1 000)。数值结果如图7 所示。借助波传播的形态可以发现，采用 200 200 的固定网格算法，其模拟的精度较差，密度分布梯度线比较粗糙，而伪弧长算法能更锐利地追踪梯度界面。虽然网格自适应移动会导致额外的局部时间迭代，但移动网格的计算代价比之 h-型方法直接加密网格的代价要小，对于求解较大尺度问题，在保证精度的同时，还应减少 CPU 运行时间，合适的弧长参数下，伪弧长算法的模拟效率要优于固定网格方法。47兵 器 装 备 工 程 学 报http:/bzxb cqut edu cn/图7二维 iemann 问题在 t=03 时刻的密度云图和网格自适应分布图Fig7 Density cloud map and mesh adaptive distribution map of two-dimensional iemann problem at t=0334双马赫反射问题该算例是一个波系结构十分复杂的流场问题。为简化模型，将计算域设置在规则区域当中:0，4 0，1。其中底部(x 1/6 为固壁边界，x 1/6 为来流)一 Mach 数为 10 的斜强激波搁置于 x=1/6，y=0 处，与 x 轴成60角，其他壁面为反射边界条件，模型介绍这里不再赘述。初值条件:U=ULfory h(x，0)UotherwiseUL=(80，57159 7，33001 2，563544)TU=(14，00，00，25)Th(x，t)=3(x 1/6)20t(34)计算终止时间 tmax=0 2，伪弧长控制函数取(a1，a2)=(1，15)，该问题在320 80 网格数下计算，如图8 所示。图8双马赫反射问题在 t=02 时刻的密度云图和网格自适应分布图(包含局部细节)Fig8 Density cloud diagram and grid adaptive distributiondiagram(including local details)of dualMach reflection problem at t=02能够验证二者在流场的大尺度现象(例如跳变点、附壁射流、马赫杆及大致滑移面等)上模拟得相差无几，而图 8(b)能更真实地反映流场内波传播的强度与小尺度结构，诸如沿滑移线的小漩涡汇聚现象、第二个三波点等。从局部放大图可以看到，相对固定网格的计算，伪弧长算法对涡核附近以及x=25 激波相汇区域的计算更为细化。PALM 更好地模拟了高马赫数的传播与反射问题。4结论1)相对于固定网格方法，伪弧长算法在提高爆轰波强间断分辨率上彰显了很好的优越性，本文中基于坐标变换策略发展的弧长算法理论，让爆轰波的求解刻画有效地避开双曲守恒系统的奇异性问题，使得 MUSCL 格式在弧长计算空间中得到很好地运用，提高了数值求解的分辨率。2)对比分析不同伪弧长控制函数对计算精度的影响，验证了监控函数模型的选择将直接影响数值模拟的效果，弧长参数决定网格移动的合理性，网格移动得越剧烈，间断附近的分辨率就越高，但可能会牺牲光滑区域的分辨率且相应增加了迭代步计算。所以采取哪类网格细化准则，构造怎样的伪弧长控制函数(比如多物理量耦合决定、复杂多项式形式等)以及如何选择可调参数，从网格生成这方面来说是一个系统问题，特别在高维空间之时。参考文献:1 马天宝，任会兰，李健，等 爆炸与冲击问题的大规模高精度计算 J 力学学报，2016，48(3):599 608MA Tianbao，EN Huilan，LI Jian，et al Large scale highprecision computation for explosion and impact problems J Chinese Journal of Theoretical and Applied Mechan-ics，2016，48(3):599 60857陈泽平，等:基于坐标变换的强间断问题伪弧长算法 2 苏罗川，宜晨虹，刘文杰，等 轻质抗侵彻材料及结构研究现状 J 兵器装备工程学报，2018，39(1):157 167SU Luochuan，YI Chenhong，LIU Wenjie，et al Developmentof lightweight ballistic armor materials and structures J Journal of Ordnance and Equipment Engineering，2018，39(1):157 167 3 NASAB A K，KILIMAN A，BABOLIAN E，et al Waveletanalysis method for solving linear and nonlinear singularboundary value problems J Applied Mathematical Model-ling，2013，37(8):5876 5886 4 VAN LEE B Towards the ultimate conservative differencescheme V Asecond-order sequel to Godunovs method J Journal of Computational Physics，1979，32:101 136 5 LUO D M，HANG W Z，QIU J X A quasi-Lagrangian mov-ing mesh discontinuous Galerkin method for hyperbolic con-servation laws J Journal of Computational Physics，2019，396:544 578 6LOPEZ M，SHONTZ S M，HUANG W Z A parallel varia-tional mesh quality improvement method for tetrahedral me-shes based on the MMPDE methodJ Computer-AidedDesign，2022，148(prepublish)doi:10 1016/J CAD2022103242 7 张亮，艾邦成，陈智 非结构有限体积耗散格式精度分析 J 兵器装备工程学报，2018，39(11):184 193ZHANG Liang，AI Bangcheng，CHEN Zhi Accuracy analy-sis of unstructured finite volume discretization scheme fordiffusive flux J Journal of Ordnance and Equipment Engi-neering，2018，39(11):184 193 8 BLAIS B，BETAND F On the use of the method of man-ufactured solutions for the verification of CFD codes for thevolume-averaged Navier-Stokes equations J Computers Fluids，2015，114:121 129 9 CHOUDHAY A，OY C J，DIETIKE J F Code verifica-tion for multiphase flows using the method of manufacturedsolutionsJ International Journal of Multiphase Flow，2016，80:150 163 10 刘学哲，林忠，王瑞利，等