imex型半隐式法在求解材料中的应用.docxVIP

imex型半隐式法在求解材料中的应用 1 爆击过程数值模拟 爆炸波和普通爆炸波之间的区别在于，首次爆炸波经过处理后，将激发未反应物，刺激高速反应，形成高温高压的爆炸物，并释放大量能量。爆炸波可以衰减传播的原因是爆炸物产生量的释放。可见, 爆轰过程是一个非常复杂的过程, 它不仅是一个流体动力学过程, 而且也是一个复杂的化学反应动力学过程, 两者互相影响、互相耦合。20世纪80年代以来, 随着计算流体力学的迅猛发展, 以及爆轰实验的高花费, 使得利用数值模拟方法对爆轰过程进行研究成为一种主流。作为气相爆轰数值研究所用的控制方程一般采用带反应源项的欧拉方程进行描述, 但是由于源项中化学反应之间以及反应和流动之间的特征尺度相差悬殊, 从而导致了数值计算的刚性问题。在处理由此引起的刚性问题时, 如果采用显式方法对控制方程进行积分, 则会严重影响计算效率;而如果采用全隐式方法, 则会占用大量的计算时间和内存。因此, 在实际应用中, 附加的半隐式方法 (The Additive Semi-Implicit Method) 得到了广泛应用, 它的原理是把ODEs分为刚性部分和非刚性部分, 刚性部分需用隐式格式计算, 而非刚性部分则用显式方法处理即可。Kennedy设计了一种高精度、L稳定的ARK2方法, 应用于对流扩散反应方程的计算中, 该方法最高有5阶精度。 基于上述分析, 本研究采用ARK2方法并结合的加权本质无震荡 (WENO) 格式对反应欧拉方程进行数值模拟。在此基础上, 对气相爆轰问题进行数值模拟, 着重研究一维、二维爆轰波的传播、波阵面复杂结构与胞格结构。 2 方程和数值法的控制 2.1 密度对欧氧基苯甲酰胺f的影响 在二维笛卡尔坐标系下, 忽略热传导, 黏性和扩散效应的守恒形式多组分反应欧拉方程为 式中:U为守恒量;F、G为通量;ρ为密度;u、v为x、y方向的速度;E为总能, E=ρh-p+ρu2/2, p为压力, h为比焓, 并且满足h=e+p/ρ (e为比内能) ;N为组分数, Yi为第i种组分的质量分数, 并且;S (U, x, y, t) 为欧拉方程中的反应源项;ωi为多组分反应中第i组分的质量生成率。 2.2 刚性问题的处理 反应欧拉方程中由于化学反应项源项之间, 化学反应和对流项之间时间尺度差别比较大, 因此在数值计算中通常会出现刚性问题, 本研究采用Additive Runge-Kutta (ARK) 算法处理非刚性对流项和刚性反应源项 式中:代数常微分方程的右端项L (U) 为非刚性项, S (U) 为刚性项, 对前者显式积分, 后者隐式处理。对流项采用5阶精度的加权本质无震荡 (WENO) 格式进行离散。 3 测试计算示例 3.1 不同网格尺寸数值模拟 为了验证本研究数值格式的可行性, 采用标准算例冲击波遇到刚壁的马赫反射问题, 选用4种不同网格尺寸进行数值模拟, 研究网格分辨率对模拟结果的影响。网格尺寸x分别为, 取相同的初始条件。计算结果如图1所示。 对比采用不同网格尺寸计算得到的马赫反射等值线分布可以发现, 本研究使用的格式能够很好地描述双马赫反射的三波结构以及马赫杆后的复杂流场。 3.2 小体积、小接触参数的控制策略 主要研究化学当量比为H2和O2、并用70%的Ar稀释的混合气体。采用一维ZND定常解作为初始条件 (见图2) , 同时在波阵面前面设置密度扰动, 使一维ZND爆轰能够迅速地形成二维胞格爆轰。边界条件采用无催化、无滑移绝热固壁条件, 对称面取物理量法向梯度为零的镜面条件。 在计算过程中, 为了降低误差, 需要对控制方程中的各物理量进行无量纲化, 无量纲化参数取特征长度L0, 未反应区声速c0、温度T0、密度ρ0, 压力p0。无量纲方法为 式中:L0=0.012m, T0=298K, p0=6.670kPa。 4 胞格结构的三维运动 图3显示了计算得到的爆轰波阵面三波结构。典型的三波结构由马赫杆、入射激波和横波组成。从图3和图4中可以看出, 三波点附近的压力、密度和温度在胞格结构中都是最高的。入射激波要比马赫杆平整得多, 然而强度却低于马赫杆, 其后的密度梯度和温度梯度也比马赫杆后小很多, 这种情况说明, 入射激波后的化学反应比马赫杆后的延迟, 相应的反应阵面滞后。同时根据马赫杆阵面的压力分布 (见图5) 可以发现, 马赫杆中心的压力最低, 在两侧方向的压力逐渐升高并在三波点附近达到峰值, 相对而言, 入射激波阵面的压力分布要均匀很多。 对于二维爆轰波来说, 胞格结构的运动模式为平面二维模式。图6为6个时刻的压力云图, 可以清晰地表明胞格内部复杂波系的运动规律。从图6中可以清晰地看出前导激波由间隔排列的马赫杆和入射激波构成, 与同一个马赫杆相交的相邻的两个横波相互远离, 与同一个入射激波相交