爆炸与冲击问题的大规模高精度计算(2)
图1 应用伪弧长方法物理量空间分布变化过程Fig.1 Physical quantities space distribution change process with pseudo arc-length method
对于多维非定常、无黏、气体导热、无气体导热的化学反应流动问题,考虑多维空间化学反应的混合气体欧拉方程组
式中,x是多维空间向量,F(w)为多维空间矩阵函数,s(w)为化学反应源项函数.
首先引入弧长参数ξ=ξ(x,t),其满足弧长表达式
进一步,可以定义监视函数
通过引入伪时间来得到多维空间中的网格移动速度的偏微分方程
对方程(1)在每个网格单元K积分得到并化简得
式中,为单元K的物理量均值,为单元K各外表面曲面.|K|为单元K体积,h(?,?,?)为数值流通量,计算中通常采用Lax--Friedrichs形式,为单元K外表面单位外法线向量,外表面内外侧函数值.
进而可对系统(4)和(5)采用龙格--库塔法耦合求解.耦合求解过程中,要保证网格物理量值插值重构的守恒性[15]
式中,wK和为重构前后的物理量值,AK和为重构前后的网格面积.
下面以一步化学反应过程为例给出计算实例
式中,R为化学反应质量分数,未反应区R为1,完全反应区R为0.指前因子为2566.4,活化能为50.一步化学反应状态方程为
式中,热释放率q为50,反应气体常数γ为1.2.以捷尔道维奇·纽曼--多灵(Zeldovich Neuman--Doring,ZND)模型解析解作为反应初值.
首先对二维爆轰波衍射问题进行研究.障碍物为反射边界,其余为无限边界.初始网格步长为?x=图2给出了当t=0.3时的模拟结果,结果表明网格可以随着爆轰波阵面而进行自适应变化.
图2 t=0.3时数值模拟结果Fig.2 Numerical results at t=0.3
三维前台阶问题:障碍物为反射边界,计算中采用并行算法编程,使用40万网格计算结果如图3.结果表明障碍物后方存在一个低密度区,并且伪弧长移动网格法能够实现三维空间中网格的自适应变化.
图3 t=0.178时,数值模拟结果Fig.3 Numerical results at t=0.178
2 基于附加龙格--库塔方法的高精度气相爆轰数值模拟
气相爆轰是流动和化学反应的强耦合过程,典型的一维爆轰波由前导激波和其后的化学反应区组成,而多维气相爆轰波波阵面则存在复杂的多波结构.近些年,爆轰波的基本理论并没有得到太大的发展,而由于观测手段的局限对爆轰波精细结构的实验研究也存在很大困难.而随着计算机技术的快速发展和计算能力的大幅提高,以基元反应模型为主的复杂模型成为气相爆轰数值模拟的重点,该模型求解的最大困难在于多组分反应欧拉方程的化学源项存在着很强的刚性,这种刚性主要体现在各个基元反应的时间尺度存在巨大的差异[16].如果使用显式格式对此类刚性方程组进行处理,格式精度所要求的时间步长将远小于由稳定性条件所决定的时间步长,这将大大增加计算时间,严重影响计算效率.由于显式方法在处理刚性问题中存在不足,隐式方法或者半隐式方法得到了很大的应用.考虑到多组分反应欧拉方程源项的刚性问题,采用附加龙格-库塔方法耦合方程中非刚性对流项和刚性反应源项.根据文献[17]中的工作,附加龙格--库塔方法能够求解下述形式的常微分方程[18-20]
式中,F(U)可以分为N项,对于反应欧拉方程,N=2,即对流项和反应源项.附加龙格--库塔方法的具体使用过程如下
式中,N项中的每一项被属于该项的s步龙格--库塔法积分.多组分反应欧拉方程可以写成Ut=Fns+Fs,其中Fns和Fs分别表示非刚性项和刚性项.附加龙格--库塔方法是一种隐--显式方法,由显式龙格--库塔方法和对角隐式龙格--库塔方法组合在一起使用,显式龙格--库塔方法积分非刚性项,对角隐式龙格--库塔方法积分刚性项.同时,非刚性对流项的空间离散通过5阶精度的加权本质无振荡(weighted essentially non-oscillatory,WENO)格式[21]完成.
应用5阶精度的加权本质无振荡格式和附加龙格--库塔方法对一维、二维气相爆轰波以及爆轰波在楔面和弯管内的传播过程进行了数值研究[19-20].气相爆轰的采用基元化学反应模型 (9组分 19反应[22]),并用70%的氩气稀释.初始压力和温度分别为6670Pa和298K.图4中给出了通过上述方法得到的二维爆轰波波阵面的精细结构,可以清楚地观察到前导激波,反应区,横波以及密度间断等结构.图5和图6分别是通过实验和数值模拟得到的胞格结构,两者都具有规则的胞格模式.
文章来源:《爆炸与冲击》 网址: http://www.bzycjzz.cn/qikandaodu/2021/0709/1298.html