宁建国 原新鹏 马天宝李 健

(北京理工大学爆炸科学与技术国家重点实验室,北京100081)

-生物、工程及交叉力学

计算动力学中的伪弧长方法研究1)

宁建国 原新鹏 马天宝2)李 健

(北京理工大学爆炸科学与技术国家重点实验室,北京100081)

动力学问题通常采用微分方程来描绘,但由于工程实际问题的复杂性,微分方程模型常伴随着解的不连续性、刚性或激波间断奇异性特点,传统方法很难求解,奇异性问题是计算动力学难点,同时也是国内外学者研究的热点.伪弧长数值算法是针对计算动力学中的奇异性问题所提出的,其基本思想为通过在解曲线上引入伪弧长参数,并增加一个约束方程,在伪弧长参数作用下,使得原始离散单元发生扭曲形变,从而达到消除或减弱奇异性的目的.本文首先介绍伪弧长方法求解定常对流--扩散方程的奇异性问题,并提出针对双曲守恒定律的局部伪弧长算法,其思想在于首先通过间断解的梯度变换来确定强间断所处位置,进而通过局部网格点重构以及数值修正来达到强间断处奇异性消除与降低的目的.针对高维问题,提出全局伪弧长方法,通过对整个计算区域内的网格点进行重构,使得所有网格点向奇异间断点处移动,从而降低间断点的影响域,达到降低奇异性的目的.重点讨论了三维全局伪弧长算法问题的计算难点,即三维空间网格扭曲大变形导致的数值算法不收敛，并提出在算法设计过程中采用分块重构与整体计算相结合的策略，实现了三维空间中的伪弧长数值算法,最后通过数值实验来验证伪弧长算法对于奇异性问题的有效性.

计算动力学,伪弧长,自适应,爆炸冲击波,三维

引言

计算动力学是利用科学计算方法来研究动力学现象与特性的科学[1].但由于动力学问题较复杂,许多实验的开展面临着周期长、代价大等困难,通过理论分析又难以获得大多数现实工程问题的解析解,数值方法就成为一种有效的替代手段.特别是在第二次世界大战以后,电子计算机的出现以及随后数十年中大型计算机的迅速普及和数值计算方法的迅猛发展[2],为计算科学的发展提供了物质与理论基础.因为动力学问题常常伴随着时间和空间的演化,所以计算动力学问题通常采用基于连续介质力学问题的微分方程模型来描述.而且由于工程实际问题的复杂性,动力系统微分方程模型常伴随着解的不连续性、刚性或激波间断等奇异性特点.奇异性问题是计算力学中的难点,同时也是国内外学者研究的热点[35].1979年,Riks[6]首次应用弧长法求解了带有参数的非线性方程中的极值问题.Boor证明了利用参数变换求解微分方程的可行性,并且提出高阶逼近函数的方法[7].之后Crisfiel[8]通过变分法将非线性几何问题的本构方程变换为非线性代数方程,进而引入弧长参数,建立了弧长算法的基本理论.此外Chan[9]在1984年通过引入伪弧长控制参数将其引申为伪弧长方法,并将其应用于求解非线性方程的奇异性问题中.在此基础之上,忻孝康等[10]、陈国谦等[11]将伪弧长方法推广，用于求解常微分方程动力系统的奇异性问题.受此启发,武际可等[12]将常微分方程动力系统的伪弧长算法推广,用来求解微分动力系统中的双曲偏微分方程的Burgers奇异性问题,总结出伪弧长算法的基本思想.通过在解曲线上引入伪弧长参数,并增加一个约束方程,在伪弧长参数作用下,使得原始离散单元发生扭曲形变,达到消除或减弱奇异性的目的.此时双曲偏微分方程的伪弧长算法的框架思想已经基本形成.与此同时,出现了一种以网格自适应为目的,通过网格点的重新分布与数值解的重新插值来求解奇异双曲问题的数值算法——移动网格方法.该方法在1983年由Harten等[13]第一次提出.Ren等[14]对该方法进行了改进,提出利用弧长参数来控制网格点的自适应移动.在此基础之上,Huang等[15]根据均分原则[16],提出了结合伪弧长控制参数的移动网格控制方程.此后,大量学者利用移动网格控制方程与双曲控制方程构成的耦合系统来研究移动网格方法[15,1718].对比分析伪弧长方法与移动网格方法的特点,可以发现,两种数值方法具有一定联系,移动网格法就是一种具有伪弧长性质的双曲方程数值算法,属于计算动力学中的伪弧长数值算法的一部分.对于双曲偏微分方程来说,虽然伪弧长方法与移动网格方法的出发点不一样,但是移动网格方法具有伪弧长数值算法的特点,符合伪弧长算法的定义,可以被称作是伪弧长移动网格算法.此外，由于移动网格方法的出发点是网格的移动重构,缺乏网格的整体性与直观性,特别是在三维问题的计算过程中,认为六面体网格在移动后,某个面发生畸变,变形为非六面体,从而导致计算失败.当前对三维问题都是采用工程的简化计算方法来完成的,所获得的计算结果不可避免地会失去原有工程问题的物理本质.再者大多数工程实际问题都是三维问题，所以针对三维问题的大规模高效算法的研究是迫切与必要的.

近年来,宁建国等将伪弧长算法用于求解爆炸冲击波问题,先后发展出局部伪弧长算法与全局伪弧长移动网格算法[1924],通过伪弧长变换来捕捉爆炸冲击波阵面,建立了爆炸冲击波问题的伪弧长算法的基础理论体系.先引入一般形式的Runge-Kutta法,得到双曲守恒控制方程与网格映射控制方程所构成的限制微分动力系统[2526],进一步利用牛顿迭代法,将非线性问题转化为类线性问题来处理,继而对类线性问题的系数矩阵进行正则化归约分析,建立了伪弧长数值算法的稳定性理论,得到了伪弧长算法的收敛性结论[24].针对反应气体动力学模型源项的刚性特点以及反应区域的冲击波与稀疏波的复杂作用导致计算过程密度、压力等物理量出现负值的问题[2728],提出了伪弧长数值算法的保正性理论[23].

本文在研究计算动力学问题中的三类伪弧长算法的基础之上,从伪弧长算法降低奇异性的特点出发,提出了三维空间中六面体网格单元分割与整体结合的思想,实现了三维空间中爆炸冲击问题的伪弧长法的数值模拟.进而针对二维、三维爆轰波对于板型障碍物的爆炸冲击问题进行数值模拟,实现了网格对于爆轰波阵面的捕捉,对障碍物后方的多个标定点处的压力大小进行分析,得到障碍物后方的安全位置.

1 定常对流--扩散问题中的伪弧长方法

定常对流--扩散问题广泛存在于化学工程、传热传质、大气和水质污染等领域中[2930].其模型通常采用常微分方程描绘,在数值计算中常常会出现伪振荡、非物理的负浓度解、激波层或边界层被拉宽等奇异性现象.采用伪弧长算法求解此类问题,可以有效地避免与降低常微分方程的奇异性问题.定常对流--扩散方程的一维无量纲形式为

其中，u为对流速度,q为源汇项,c为污染物的浓度,Pe称为Peclet数,定义为

式中，U为特征对流速度,L为特征长度,α为扩散系数.

在大Peclet数情况下,即Pe≫1,可令

方程(1)可以写成奇异摄动型的二阶常微分方程式

其中，y=c,f,g,h均可是x和y(c)的函数,α,β为常数.引入解曲线的弧长s,则有

令v=dy/dx,则式(4)可转化为如下一阶常微分方程组

假设x=a端的弧长为零,x=b端的弧长为Smax,于是式(5)化为如下的边界条件

因此,奇异摄动两点边值问题式(4)和式(5)就转化为一阶常微分方程的边值问题式(7)和式(8).对于方程组(7)通常采用Rosenbrock格式求解[31].对于边值问题式(8)可采用打靶法进行求解[32].

2 非定常对流问题的伪弧长方法

非定常对流问题通常采用双曲型偏微分方程来刻画[33].这类问题通常伴随着激波间断性而导致解出现强间断奇异性[34].数值求解这类奇异间断性问题的核心在于对间断处的精确捕捉与计算.对于这类问题的求解,伪弧长方法包括局部伪弧长方法和全局伪弧长方法.局部伪弧长方法为伪弧长方法的早期形式,其核心在于首先通过间断解的梯度变换来确定强间断所处位置,进而通过局部网格点重构以及数值修正来达到强间断处奇异性消除与降低的目的,如图1(a)所示.而全局伪弧长方法(伪弧长移动网格法)则是通过对整个计算区域内的网格点进行重构,使得所有网格点向奇异间断点处进行移动,从而降低间断点的影响域,进而达到降低奇异性的目的,如图1(b)所示.

2.1 局部伪弧长方法

考虑一维双曲守恒系统

图1 双曲微分方程伪弧长方法示意图Fig.1 Schematic diagram for the hyperbolic di ff erential equations pseudo arc-length method

利用雅克比矩阵A(U)=∇UF,守恒型方程可以转换到如下形式

引入弧长参量ξ,其满足关系式

其中u是U的分量形式,进而由上式可得

进而联立控制方程(10)与伪弧长限制方程(13)可以得到

中,采用下式计算

其中ε是一个小的正数.进而对式(14)进行空间离散,得

在计算空间中可采用均匀网格,ui=u(ξi,t),Δξ=ξi+1-ξi为计算空间网格步长,物理空间网格步长为Δxi=x(ξi,τ+ Δτ)-x(ξi,τ).进一步可对式 (16)利用Runge-Kutta[3536]进行时间离散求解.

为提高上述过程计算效率,令参数

其中,Δu=u(xi+1,t)-u(xi,t),Δx为物理空间网格步长,在求解过程中,对于每个离散单元,检验Φ值,当Φ≤Φ0(其中Φ0为一个小的正数),采用引入弧长变量ξ的方程,而其他的单元仍采用原有的方程.

2.2 全局伪弧长方法

由于在高维空间中,间断点的跟踪与控制方程在间断处的变形修正都是十分困难的,因此相比于局部伪弧长算法,全局伪弧长方法更适合处理二维、三维问题.在二维空间中,网格变形如图2所示,在计算过程中需要保证每个单元网格面积大于零,避免网格面积为零或为负值的出现,相比于二维问题,三维问题要复杂很多,由于在三维空间中六面体单元畸变以后,其外表面各点会出现不共面的情况,如图3所示,这样会导致计算过程中的物理量重构不守恒,以及扭曲变形后网格单元的控制方程演化失败.为解决这一难题,本文采用分块重构与整体计算相结合的策略,即在物理量重构阶段,将六面体单元每个空间外表面拆分成两个面(如图4过程),则每个六面体在X,Y,Z每个方向上被拆分成两个三棱柱分别进行处理(如图5所示).在网格移动与控制方程的演化计算阶段再回归整体单元区域进行计算.分块重构过程可以避免网格的重构过程中由于多个点不在一个平面导致物理量重构不守恒问题.整体移动与计算可以避免多个块体分别计算导致计算量过大以及多个块体的同时移动导致网格错位以及交叉网格的出现.

图2 二维空间网格变形示意图Fig.2 Two-dimensional spatial grid deformation

图3 三维空间网格变形示意图Fig.3 Three-dimensional grid deformation

图4 三维曲面计算示意图Fig.4 Computational diagram of three-dimensional curved surface

对于三维非定常、无黏、无热传导反应流体动力学问题,考虑三维空间中一步反应流体欧拉方程组

其中,x为三维物理空间向量；三维空间矩阵函数

物理量 w=(ρ,ρu,ρv,ρw,E,ρR)T；化学反应源项函数为s(w)=(0,0,0,0,0,ω)T.对于一步Arrhenius化学反应模型ω

一步化学反应状态方程为

引入弧长参数ξ=ξ(x,t),其满足弧长表达式

进一步,我们可以定义监视函数

通过引入伪时间来得到多维空间中的网格移动速度的偏微分方程

对上式可以采用Gauss-Seidel循环迭代求解

其中

对方程(18)在每个网格单元上K积分得到并化简得

其中

进而可对系统式 (22)和式 (23)采用 Runge-Kutta[3536]耦合求解.耦合求解过程中,要保证网格物理量值插值重构的守恒性

其中,wK,K为重构前后的物理量值,AK,K为重构前后的网格体积.令Dl为式(24)一次循环求解后,外表面所扫过的体积.所以得到对于单个网格的守恒插值策略

其中,Dl为外表面扫过区域产生的体积改变量,这里以计算D1为例,D1为图6所示的三棱柱块体.D1可以分解为三个四面体进行求解,因此可得

通过四面体体积公式求解出每个四面体有向体积VOABC，即当ABC为逆时针排列时,其值为正,顺时针排列时,其值为负.

图6 单个网格面扫过的体积Fig.6 Swept volume of the typical mesh surface

除此之外,对于化学反应流问题,在计算过程中要应用保正性策略,保证计算过程中的压力、密度、化学反应质量分数等为正.笔者在文献[23]已经讨论了二维空间中保正性问题,建立了全局伪弧长保正性算法的体系结构.全局伪弧长算法的保正性理论包括两个方面：(1)控制方程离散的保正性;(2)物理量自适应重构的保正性.其实现步骤也包括两个方面：(1)通过添加时间步长和网格移动步长的约束条件实现网格均值的保正性;(2)通过补充算法实现网格的每个点物理量值的保正性.

3 数值算例

算例1 首先通过具有解析解的一维Sod’s问题来说明全局伪弧长算法能够降低激波间断奇异性.对于Euler方程(18)的一维初值问题

计算域取[0,1],最终计算时间为T=0.2.计算过程中取监视函数为

表1 中对比给出伪弧长法与固定网格(fi mesh)方法的计算效果.从表 1中可以看出,利用伪弧长(pseudo arc-length)法50个网格点就可以获得比固定网格方法100个网格点更好的计算效果,而且计算时间并没有增加.利用100个网格点的不同的弧长参数的伪弧长法可以获得比固定网格方法150,200,300个更好的计算效果.特别是L2与L∞误差范数显著降低,说明伪弧长算法对于激波间断奇异性问题有很好的计算效果.

表1 伪弧长法与固定网格方法计算效果对比Table 1 Computational comparison for pseudo arc-length method and fi ed grid method

算例2 考虑二维板型障碍物爆轰衍射问题.计算域如图7所示,红色为障碍物区域,其空间位置为[1.6,1.9]×[0,1.0].初始反应区选定为X方向的小于1的区域.计算过程中取未反应区R为1,完全反应区R为0.指前因子˜K为2566.4,活化能˜T为50.一步化学反应状态方程为

其中,热释放率q为50,反应气体常数γ为1.2.以Zeldovich Neuman-Doring(ZND)模型[3738]解析解作为反应初值.初始计算域如图7所示,障碍物为刚性反射边界,地面为刚性反射边界,计算域前后及上方均为无限边界,初始网格步长为Δx=Δy=图8给出了当t=0.35时的模拟结果.结果表明,网格随着爆轰波阵面而进行自适应变化,实现了网格对爆炸冲击波阵面的捕捉.选定板型障碍物前后12个位置点,根据位置点与障碍物的位置关系分成4组，考察其各点处的压力变化,12点的位置如图7中所示,其压力变化如图9所示.选定一般研究认为的压力安全线3.0为参考压力线.从图9中可以看出,在障碍物前方区域压力较高,为危险区域.在障碍物后方区域中,位置较高点的爆轰波压力较大,位置较低的点压力相对较小,且越靠近障碍物的后方区域,其安全性越高.在本文所选取的12个点中,安全性最高的是e(2.0,0.5)位置的点.f(2.0,0.1)位置由于地面反射波的原因,稍有危险.所以在爆炸发生时,最安全的区域是障碍物后方的中部位置.

图7 二维爆轰模拟初态示意图Fig.7 Initial state of two-dimensional detonation

图8 二维板型障碍物绕射问题Fig.8 Two-dimensional di ff raction problem for plate-shape obstacle

图9 二维板型障碍物绕射各位置点压力变化Fig.9 Pressure history of typical locations in two-dimensional plate-shape obstacle di ff raction problem

图9 二维板型障碍物绕射各位置点压力变化(续)Fig.9 Pressure history of typical locations in two-dimensional plate-shape obstacle di ff raction problem(continued)

算例3 下面分析考察三维空间中的板型障碍物绕射问题.计算域如图10(a)所示,红色为障碍物区域,其空间位置为[1.6,1.9]×[0,1.0]×[0,1.0].初始反应区选定为X方向的小于1的区域,其余为未反应区.计算过程中取未反应区R为1,完全反应区R为0.计算参数与算例2相同,并且同样以Zeldovich Neuman-Doring(ZND)模型解析解作为反应初值.障碍物为刚性反射边界,地面为刚性反射边界,计算域前后左右及上方均为无限边界.计算中采用100万网格.图11中给出的是三维空间中压力与密度云图,计算结果表明障碍物后方存在一个低密度区.图12中给出障碍物处切片的网格分布图,特别是图12(b)中给出图12(a)中A部分的网格变形前后的对比图,黑色为网格畸变前的,红色表示网格畸变后的.图13中给出非障碍物处切片的网格分布图,可以明显看出三维空间中网格的移动变形.计算结果表明,伪弧长算法能够实现三维空间中网格的自适应变化.图14中,给出图11中点C处的密度计算效果对比,其中参照解是采用2000万网格的计算效果,从图中可以看出,伪弧长算法的计算效果与参照解的计算效果基本一致,且明显优于固定网格算法的100万网格的计算效果.

图10 三维爆轰模拟初态示意图Fig.10 Initial state diagram of three-dimensional detonation

图11 三维板型障碍物绕射问题Fig.11 Three-dimensional di ff raction problem for plate-shape obstacle

图12 三维板型障碍物绕射问题障碍物处切片Fig.12 Slices at obstacle for three-dimensional plate-shape obstacle di ff raction problem

图13 三维板型障碍物绕射问题非障碍物处切片Fig.13 Slices at non-obstacle for three-dimensional plate-shape obstacle di ff raction problem

图14 点C处的密度计算效果对比Fig.14 Comparison of densities at point C

进而选取障碍物前后X方向的a,b,c,d四个横切面,如图10(a)所示,继而在每个横切面上选取5个位置点,各点的选取如图10(b)所示.根据数值模拟结果测定各个位置点处的压力变化,如图15所示.图15(a)中给出的是障碍物前方a位置处各点的压力变化曲线.从图中可以看出,板前区域的压力普遍较大,均远高于压力安全线3.0.图15(b)中给出的是障碍物后方临近障碍物b位置处的若干点的压力变化情况.从图中可以看出位置b4处的压力是处于安全线以下的.图15(c)和图15(d)给出的是障碍物后方c位置和d位置各点的压力变化曲线.从图中可以看出,其各点均处于危险区域中,但是d位置处的最大压力峰值相对滞后.图16中给出了更为直观的各个位置点处的最大压力分布图,从图中可以明显看出第一组的峰值压力高于其他各组,最小的峰值压力在第二组中.

图15 各位置点处压力变化Fig.15 Pressure history of typical locations

图16 各位置点处最大压力值Fig.16 Maximum pressure at typical locations

4 总结与讨论

本文针对于计算动力学问题,分析讨论了定常对流--扩散问题的伪弧长法以及求解非定常对流问题的伪弧长法.其中非定常对流问题的伪弧长法包括局部伪弧长法和全局伪弧长法.针对三维爆炸冲击奇异性问题的计算难点,提出了爆炸冲击问题的三维全局伪弧长算法,与前人研究的移动网格方法不同,由于移动网格方法的出发点是网格的移动重构,缺乏网格的整体性与直观性,认为六面体网格移动后,某个面发生畸变,变形为非六面体,从而导致计算失败.而伪弧长算法的计算目标不是网格移动,而是减少奇异性.算例表明,本文提出的六面体网格单元分割与整体结合的计算方法成功实现了三维空间中的爆炸与冲击问题的数值模拟.针对二维,三维爆轰波对于板型障碍物的爆炸冲击问题,本文得到了以下结论：

(1)伪弧长算法可以有效实现网格对于爆轰波阵面的捕捉,从而提高计算精度,其计算精度要明显优于固定网格方法.

(2)障碍物前方区域均为危险区域.

(3)障碍物后方临近障碍物的区域会存在某些安全区域.

(4)在障碍物后方临近障碍物的区域中,远离障碍物边缘与地面的地方为最安全区域.在二维空间中,最安全的区域是障碍物后方的中部位置.在三维空间中,最安全的区域是障碍物后方离地面与边缘最远的位置.

1马天宝,任会兰,李健等.爆炸与冲击问题的大规模高精度计算.力学学报,2016,48(3)：599-608(Ma Tianbao,Ren Huilan,Li Jian,et al.Large scale high precision computation for explosion and impactproblems.ChineseJournalofTheoreticalandAppliedMechanics,2016,48(3)：599-608(in Chinese))

2杨露斯,黎炼.论计算机发展史及展望.信息与电脑,2010,6(1)：188-188(Yang Lusi,Li Lian.Theory of computer development and prospects.China Computer&Communication,2010,6(1)：188-188(in Chinese))

3 Luo ST,Tran HV,Yu Y.Some inverse problems in periodic homogenization of hamilton-jacobi equations.Archive for Rational Mechanics and Analysis,2016,221(3)：1585-1617

4 Deleuze Y,Chiang CY,Thiriet M,et al. Numerical study of plume patterns in a chemotaxis–di ff usion–convection coupling system.Computers&Fluids,2016,126(1)：58-70

5 Tian DS.Multiple positive periodic solutions for second-order differential equations with a singularity.Acta Applicandae Mathematicae,2016,144(1)：1-10

6 Riks E.An incremental approach to the solution of snapping and buckling problems.International Journal of Solids&Structures,1979,15(7)：529-551

7 Boor CD.On local spline approximation by moments.Journal of Mathematics and Mechanics,1968,17(1)：729-735

8 Crisfiel MA.An arc-length method including line searches and accelerations.International Journal for Numerical Methods in Engineering,1983,19(1)：1268-1289

9 Chan TF.Newton-like pseudo-arclength methods for computing simple turning points.Siam Journal on Scientifi&Statistical Computing,1984,5(1)：135-148

10忻孝康,唐登海.一维定常对流–扩散方程的伪弧长参数法.应用力学学报,1988,5(1)：45-54(Xin Xiaokang,Tang Denghai.An arc-length parameter technique for steady di ff usion-convection equation.Chinese Journal of Applied Mechanics,1988,5(1)：45-54(in Chinese))

11陈国谦,高智.对流扩散方程的迎风变换及相应有限差分方法.力学学报,1991,23(4)：418-425(Chen Guoqian,Gao Zhi.Transformation of the convective di ff usion equation with corresponding finit di ff erence method.Chinese Journal of Theoretical and Applied Mechanics,1991,23(4)：418-425(in Chinese))

12武际可,许为厚,丁红丽.求解微分方程初值问题的一种弧长法.应用数学和力学,1999,20(8)：875-880(Wu Jike,Xu Weihou,Ding Hongli.Arc-lengthmethodfordi ff erentialequations.AppliedMathematics and Mechanics,1999,20(8)：875-880(in Chinese))

13 Harten A,Hyman JM.Self adjusting grid methods for onedimensional hyperbolic conservation laws.Journal of Computational Physics,1983,50(2)：235-269

14 Ren YH,Russell RD.Moving mesh techniques based upon equidistribution,and their stability.Siam Journal on Scientifi&Statistical Computing,1992,13(6)：1265-1286

15 Huang WZ,Ren YH,Russell RD.Moving mesh partial di ff erential equations(MMPDEs)based upon the equidistribution principle.Siam Journal on Numerical Analysis,1994,31(3)：709-730

16 White AB.On selection of equidistributing meshes for two-point boundary-value problems.Siam Journal on Numerical Analysis,1979,16(3)：472-502

17 Stockie JM,Mackenzie JA,Russell RD.A moving mesh method for one-dimensional hyperbolic conservation laws.Siam Journal on Scientifi Computing,2001,22(5)：1791-1813

18 Tang HZ,Tang T.Adaptive mesh methods for one-and twodimensional hyperbolic conservation laws.Siam Journal on Numerical Analysis,2003,41(2)：487-515

19 Ning JG,Wang X,Ma TB,et al.Numerical simulation of shock wave interaction with a deformable particle based on the pseudo arclength method.Science China Technological Sciences,2015,58(5)：848-857

20王星,马天宝,宁建国.双曲偏微分方程的局部伪弧长方法研究.计算力学学报,2014,31(3)：384-389(Wang Xing,Ma Tianbao,Ning Jianguo.Local pseudo arc-length method for hyperbolic partial di ff erential equation.Chinese Journal of Computational Mechanics,2014,31(3)：384-389(in Chinese))

21 Wang X,Ma TB,Ning JG.A pseudo arc-length method for numerical simulation of shock waves.Chinese Physics Letter,2014,31(3)：030201-030201

22 Wang X,Ma TB,Ren HL,et al.A local pseudo arc-length method for hyperbolic conservation laws.Acta Mechanica Sinica,2014,30(6)：956-965

23 Ning JG,Yuan XP,Ma TB,et al.Positivity-preserving moving mesh scheme for two-step reaction model in two dimensions.Computers&Fluids,2015,123(1)：72-86

24 Yuan XP,Ning JG,Ma TB,et al.Stability of newton tvd runge–kutta scheme for one-dimensional Euler equations with adaptive mesh.Applied Mathematics&Computation,2016,282(1)：1-16

25 Sun LP,Cong YH,Kuang JX.Asymptotic behavior of nonlinear delay di ff erential–algebraic equations and implicit Euler methods.Applied Mathematics&Computation,2014,228(1)：395-403

26 Harwood SM,Kai H,Barton PI.Efficient solution of ordinary differential equations with a parametric lexicographic linear program embedded.Numerische Mathematik,2016,133(4)：623-653

27 Perthame B,Shu CW.On positivity preserving finitvolume schemes for Euler equations.Numerische Mathematik,1996,73(1)：119-130

28 Wang C,Zhang X,Shu CW,et al.Robust high order discontinuous galerkin schemes for two-dimensional gaseous detonations.Journal of Computational Physics,2012,231(2)：653-665

29魏涛,许明田,汪引.求解对流扩散问题的积分方程法.化工学报,2015,66(10)：3888-3894(Wei Tao,Xu Mingtian,Wang Yin.Integral equation approach to convection-di ff usion problems.Ciesc Journal,2015,66(10)：3888-3894(in Chinese))

30 Zhou XF,Guo J,Li F.Stability,accuracy and numerical di ff usion analysis of nodal expansion method for steady convection di ff usion equation.Nuclear Engineering&Design,2015,295(1)：567-575

31 Sandu A.Rosenbrock methods with an explicit firs stage.International Journal of Computer Mathematics,2015,93(6)：1-18

32冯帆,王自发,唐晓.一个基于打靶法的大气污染源反演自适应算法.大气科学,2016,40(4)：719-729(Feng Fan,Wang Zifa,Tang Xiao.Development of an adaptive algorithm based on the shooting method and its application in the problem of estimating air pollutant emissions.Chinese Journal of Atmospheric Sciences,2016,40(4)：719-729(in Chinese))

33 Chen Z,Huang H,Yan J.Third order maximum-principle-satisfying direct discontinuous galerkin methods for time dependent convection di ff usion equations on unstructured triangular meshes.Journal of Computational Physics,2015,308(1)：198-217

34刘朝阳,王振国,孙明波,等.一种求解激波问题的中心差分--WENO混合方法研究.推进技术,2016,37(8)：1522-1528(Liu Chaoyang,Wang Zhenguo,Sun Mingbo,et al.Investigation of a hybrid central di ff erence-WENO method for shock wave problems.Journal of Propulsion Technology,2016,37(8)：1522-1528(in Chinese))

35 Shu CW,Osher S.Efficient implementation of essentially nonoscillatory shock-capturing schemes. Journal of Computational Physics,1989,83(1)：32-78

36 Gotlieb S,Shu CW.Total variation diminishing Runge-Kutta schemes.Mathematics of Computation,1996,67(221)：73-85

37 Doring W.On detonation processes in gases.Annals of Physics,1943,43(9)：421-436

38 Fickett W,Wood WW.Flow calculations for pulsating onedimensional detonations.The Physics of Fluids,1966,9(5)：903-916

PSEUDO ARC-LENGTH NUMERICAL ALGORITHM FOR COMPUTATIONAL DYNAMICS1)

Ning Jianguo Yuan Xinpeng Ma Tianbao2)Li Jian
(State Key Laboratory of Explosion Science and Technology，Beijing Institute of Technology，Beijing 100081，China)

Di ff erential equations are often used to describe the dynamic problems.Classical approaches are always hard to solve it in engineering practice due to its characteristics of strong discontinuity,rigidity and shock singularity,among which singularity problem is one of the most difficult and hot issues among scholars.Pseudo arc-length numerical algorithm is proposed for singularity problems in computational dynamics,whose basic idea is to introduce a pseudo arc-length parameter in the original governing equations so that a constraint equation is added.Under the action of a pseudo arc-length parameter,the original discrete elements are distorted to achieve the goal of eliminating or weakening singularity.Firstly,this paper gave an introduction about the pseudo arc-length method for solving the singularity problem in steady di ff usion-convection equations.Then the local pseudo arc-length algorithm is proposed for hyperbolic conservation laws.This algorithm has two steps.The firs step is to determine the location of the strong discontinuity and the second one aims to eliminate or reduce the singularity by reconstructing the local mesh.Secondly,a global pseudoarc-length algorithm is put forward for high dimensional problems,which can reconstruct the mesh in whole area.Since the reconstructed mesh can adaptively catch the singularity points,the singularity is greatly reduced.Thirdly,the threedimensional global pseudo arc-length algorithm and its computational difficulties in numerical algorithm convergence caused by large grid distortion in three-dimensional space are presented.Then the combination of block reconstruction and integral calculation strategy is adopted in the algorithm design process to achieve the pseudo arc-length numerical algorithm in three-dimensional space.Finally,numerical examples were employed to verify the validity of the pseudo arc-length algorithm.

computational dynamics,pseudo arc-length,adaptive,explosion shock wave,three-dimensional

O302,O175

：A

10.6052/0459-1879-16-385

2016–12–19 收稿，2017–03–14 录用,2017–03–14 网络版发表.

1)国家自然科学基金资助项目(11390363,11532012).

2)马天宝,副教授,主要研究方向：计算爆炸力学.E-mail：madabal@bit.edu.cn

宁建国,原新鹏,马天宝,李健.计算动力学中的伪弧长方法研究.力学学报,2017,49(3)：703-715

Ning Jianguo,Yuan Xinpeng,Ma Tianbao,Li Jian.Pseudo arc-length numerical algorithm for computational dynamics.Chinese Journal of Theoretical and Applied Mechanics,2017,49(3)：703-715