基于OpenFOAM的三维内弹道两相流程序开发及其应用

肖辉鹏,王 浩,陶如意

(1.南京理工大学 能源与动力工程学院,江苏 南京 210094;2.中国船舶重工集团公司750试验场,云南 昆明 650051)

近年来,随着计算机硬件的快速发展与多相流模型及其数值算法的日臻成熟,计算内弹道学为探究武器发射系统燃烧室内流场变化规律提供了重要研究手段。计算内弹道从最初的零维集总参数法发展到双一维两相流计算内弹道,甚至到现在的二维或者三维计算内弹道两相流技术,在这个过程中内弹道两相流数值计算取得了诸多的研究成果[1-4]。

目前,在应用内弹道两相流算法针对不同的武器发射系统进行内弹道仿真分析时,计算内弹道两相流算法往往是独立编程实现的。独立编程实现内弹道两相流程序存在如下问题:①仿真周期长,独立编写的内弹道两相流程序都是针对每一个问题专门编写一个程序,因此针对工程问题仿真周期比较长;②代码复用率低,代码一般采用面向过程的编程思想,而且针对不同结构燃烧室开发不同的仿真程序;③不具备复杂计算域的求解能力,在代码开发前期往往需要将模型简化为一维或者二维模型,因此所开发的程序不能对同结构燃烧室进行三维内弹道仿真。

为了克服上述问题,本文拟在开源计算流体软件OpenFOAM框架上开发三维内弹道两相流求解器,减少内弹道仿真周期,提高代码复用率以及该程序复杂计算域求解的能力。

1 数学物理模型

内弹道两相流基本假设详见文献[1],采用双欧拉描述的三维内弹道两相流数学物理模型如下所示。

1)气相质量方程。

(1)

2)固相质量方程。

(2)

3)气相动量方程。

(3)

式中:p为当地气相压力,Fdrag为相间阻力。

4)固相动量方程。

(4)

式中:Rp为颗粒间的应力。

5)气相能量方程。

(5)

式中:hg为火药燃气比总焓,即单位质量火药气体的总焓;Wg为火药气体对火药颗粒所做的功;Qg为相间热传递热量。

相关辅助方程,限于篇幅本文不再列出,详细可见文献[5]。

2 压力修正算法

压力修正算法是以迭代的方式求解非线性耦合的动量方程与压力方程的数值方法,到目前为止,已经将压力修正算法成功地应用于多相流全声速流动问题。本文的三维内弹道两相流求解器主体是在开源计算流体软件OpenFOAM中的twoPhaseEuler求解器框架上实现的。twoPhaseEuler求解器是采用两相流压力修正算法实现的,该求解算法主要以学者Rushe和Weller所做的工作为理论依据[6-7]。

2.1 控制方程的离散

控制方程的离散包括质量方程的离散、动量方程的离散、能量方程的离散以及压力方程的离散。通过质量方程离散可以得到相分数方程,通过动量方程离散可以得到速度预测方程。上述方程的离散及辅助方程具体细节参考文献[5]所做的工作。本文着重讨论火药着火模型数值计算方法和火药颗粒特征尺寸数值计算方法,该数值模型是实现内弹道两相流求解器的关键。

2.1.1 火药颗粒特征尺寸数值计算方法

内弹道两相流数值问题,常常需要计算与火药颗粒特征尺寸有关的物理量,例如:计算雷诺数、放热系数时都需要预先得到火药颗粒的等效直径;计算火药生成量需要预先得到火药颗粒表面积。

由于火药颗粒燃烧采用等面燃烧的假设,因此只要计算得到火药已燃相对厚度,其他的火药特征参数就可以导出。本文采用火药已燃相对厚度作为待求解的火药特征参数,其他火药特征参数可由火药已燃相对厚度导出。假设空间任意一点的火药已燃相对厚度为位于该点控制体内的所有火药颗粒已燃厚度的平均值,为了能在欧拉坐标系下捕捉火药已燃相对厚度,根据文献[8]计算火药特征等效直径思想,建立已燃相对厚度的输运方程:

(6)

2.1.2 火药着火模型数值计算方法

火药表面温度着火模型是采用拉格朗日坐标系描述的,为了使火药着火模型在内弹道双欧拉两相假设下实现,将火药着火方程转化为欧拉描述:

(7)

式中:ap为火药导温系数,λp为火药导热系数,δt为时间增量。对于单相不可压缩流体而言,速度v的散度为0,即·v=0,将其带入方程(7)中可得:

(8)

火药表面温度方程采用全显式求解,根据线性叠加原理,火药表面温度升高为由输运引起的温度升高与源项引起的温度升高之和,火药表面温度方程可写成如下形式:

(9)

(10)

2.2 内弹道两相流求解器的实现

内弹道两相流求解的实现需在开源计算流体OpenFOAM框架中加入与火药相关的物理模型,及对两相流twoPhaseEuler求解器的修改,主要包括在相应的位置对源项、相间阻力项、颗粒间应力项的处理,以及火药表面温度计算与火药颗粒特征参数计算。

2.2.1 相关物理模型的实现

火药气体模型:在OpenFOAM中没有火药气体,火药气体模型主要提供火药气体可压缩性的计算以及气体密度的计算。

火药颗粒间阻力模型:火药颗粒阻力模型主要提供系数Ad的计算,因为相间阻力项的离散采用半隐藏离散,因此需要计算系数Ad,该系数定义为

(11)

式中:Cf为火药颗粒床摩擦系数。

火药颗粒应力模型:该模型提供火药颗粒之间应力的计算,火药颗粒间应力采用显式离散。

2.2.2 内弹道两相流程序求解流程

内弹道两相流程序计算流程主要求解过程如下:

①开始PIMPLE压力修正算法循环。

②通过固相的相分数方程求解固相相分数,用得到的固相相分数φp更新火药气体的相分数(空隙率),气体相分数通过关系式φg=1-φp得到。

③更新与相间阻力有关的系数Ad以及火药颗粒间的应力Rp。

④用前一次迭代计算得到的速度构建动量方程。

⑤求解气相能量方程,用计算得到的火药气体能量更新火药气体的温度场、密度场。

⑥压力校正循环。

a. 分别用固相通量预测方程与气相通量预测方程计算预测通量;

b. 构造并求解压力方程;

c. 修正固相通量方程与气相通量方程;

d. 重构固相速度与气相速度。

⑦计算与火药颗粒有关的参数。

a. 求解火药表面温度输运方程,计算火药表面温度;

b. 求解火药已燃相对厚度输运方程,计算火药已燃相对厚度;

c. 更新火药颗粒表面积,更新火药体积、火药等效直径、火药燃烧线速度;

d. 求解火药生成的质量。

图1给出了内弹道两相流程序计算流程图。图中,n(Correctors),n(OutCorrectors)为OpenFOAM计算手册中规定的标准参数,分别表示压力修正循环需要迭代的次数和PIMPLE循环需要迭代的次数。

图1 内弹道两相流程序计算流程

3 算例验证

本文从三方面对基于OpenFOAM实现的三维内弹道两相流程序进行了验证:一是压力波捕捉能力验证,二是全局守恒性验证,三是内弹道试验验证。

3.1 压力波捕捉能力验证

内弹道两相流流动问题常常伴随压力波的传播,精确捕捉压力波传播过程中波形与位置是实现内弹道两相流程序的重点[4]。为了捕捉压力波的传播,内弹道两相流程序需具备高分辨压力波捕捉能力。为了验证基于OpenFOAM实现的内弹道两相流程序的压力捕捉能力,本文拟以“经典算例”一维激波管问题为例,将一维激波管问题的精确解与内弹道两相流程序得到的数值解进行对比。

激波管的结构示意图如图2所示,激波管总长为10 m,左端高压区域长度为5 m,右端低压区域长度为5 m,低压区域与高压区域的初始条件如表1所示。

图2 激波管示意图

表1 激波管初始条件

区域p/kPaT/Kρ/(kg·m-3)高压区域1003001.158 623 0低压区域103000.115 862 3

图3分别给出了t=7 ms时,一维激波管问题OpenFOAM计算得到的数值解与解析解的压力p和速度vx曲线的对比图。从图中可以得到,精确解与OpenFOAM数值解无论是在激波管的膨胀区域还是间断区域都符合较好,说明基于OpenFOAM的内弹道两相流程序具有足够的精度,对压力波捕捉具有高分辨率。

图3 t=7 ms时激波管问题数值解与精确解的对比

3.2 全局守恒性验证

分别采用基于OpenFOAM的内弹道两相流程序与经典内弹道集总参数法仿真密闭空间内火药颗粒被引燃的问题,验证密闭空间总质量守恒性以及密闭空间的总能量守恒性,最后对内弹道两相流求解器与集总参数法计算得到的内弹道物理量曲线进行对比。

在10 cm×10 cm×10 cm密闭空间内装填46.12 g的2#小粒黑火药。采用基于OpenFOAM实现的三维内弹道两相流程序模拟上述装填条件下火药在密闭空间内燃烧过程,整体计算域网格示意图如图4(a)所示。网格总数为8 000,每个正方形单元网格边长为0.5 mm,2#小粒黑火药均匀分布在计算域中心的4 cm×4 cm×4 cm正方形区域内,如图4(b)所示。

图4 计算域网格

表2给出了2种方法计算得到的物理量,表中,te表示火药燃烧结束时间,me表示火药燃烧结束时燃气总质量,Ee表示火药燃烧结束时燃气总能量。密闭空间中生成的火药燃气总质量结果分别为46.33 g与46.12 g,误差为0.45%,可以认为OpenFOAM实现的三维内弹道两相流程序满足全局质量守恒性。密闭空间生成火药燃气总能量值分别为47 725 J与47 710 J,误差为0.03%,可以认为火药燃烧释放的能量等于火药燃气增加的总能,这验证了OpenFOAM实现的三维内弹道两项流程序满足全局能量守恒性。

表2 计算结果的比较

图5为OpenFOAM内弹道两相流与集总参数法计算得到的已燃相对厚度曲线与压力曲线对比图,其中OpenFOAM计算得到的已燃相对厚度与压力是取所有单元的平均值。由图可知,2种方法计算得到的对应物理量值几乎近似相等,其原因如下:集总参数法将火药的三维密闭燃烧问题转换为零维问题,所有物理量取其平均值,不考虑压力波的传播;而OpenFOAM内弹道两相流方法,虽然考虑了压力波在三维密闭空间的传播,但是火药全部被点燃的初始条件、装药条件以及密闭空间形状,共同决定了空间压力扰动能迅速传到三维空间每个计算单元所在的位置。

图5 OpenFOAM内弹道两相流与集总参数计算结果对比曲线

3.3 试验验证

图6为中心炸管式抛撒系统结构示意图。本节拟采用基于OpenFOAM的内弹道两相流程序模拟抛撒机构定容阶段内弹道过程,将测压孔处试验测得的压力曲线与数值计算得到的压力曲线进行对比,分析抛撒机构的点传火过程。

图6 中心炸管式抛撒系统结构示意图

3.3.1 仿真模型的建立

本文将全部火药燃气可流动区域作为求解计算域。为了减少计算量,建立抛撒系统燃气可流动区域的轴对称模型,计算域网格如图7所示。

图7 求解计算域网格示意图

3.3.2 仿真与试验结果的对比

图8给出了测压孔位置测试结果与计算结果的压力-时间对比曲线。

图8 测压孔位置压力-时间对比曲线

由图8可见,在0～2 ms时测压孔位置的压力曲线一直基本保持水平。在中心炸管式抛撒系统进行点传火过程中,当点火具达到破膜压力后,点火具内火药燃气从点火基座的小孔喷射到传火管中,传火管中高温高压火药燃气一部分从传火孔中流出,进入燃烧室并引燃抛撒药,另一部分在传火管内沿着轴向流动;在t=2 ms时,从I区压力-时间曲线局部放大图中可以看出,测压孔处的压力开始上升,说明点传火过程已经从中心管顶端位置进行到中心管底端位置;在t=5 ms时,中心管内的抛撒药进入全部燃烧阶段,抛撒药迅速燃烧,中心管内的压力呈指数增长的趋势迅速增加,直至测压孔压力达到37.6 MPa时中心管炸裂。

通过比较测压孔位置的测试压力曲线和计算压力曲线,可以看出计算结果与试验结果符合较好,说明基于OpenFOAM的内弹道两相流程序能够准确地模拟中心炸管系统实际点传火过程,计算程序能为类似结构燃气发生器优化提供良好的理论依据。

3.3.3 点传火过程特性分析

图9给出了不同时刻中心炸管式抛撒系统燃烧室内火药着火程度云图,图中,T*为无量纲温度。云图中主要有3种颜色分布,红色区域表示该位置的火药颗粒已经达到着火温度,对于已经达到着火温度(大于600 K)的火药颗粒在云图中用T*=1表示;灰白色区域表示火药颗粒表面温度为常温(300 K),在云图取T*=0.5;蓝色区域没有火药颗粒的存在,在编程处理时,只有存在火药颗粒的区域火药颗粒表面温度才有值,对于没有火药颗粒存在的区域火药表面温度的值取T*=0。

图9 不同时刻的火药着火程度云图

图9(a)～9(c)说明点火基座小孔破膜后,高温火药燃气与火药颗粒从点火基座上的小孔流入传火管中,一部分高温点火药燃气与已燃火药颗粒从传火孔流入中心管,另一部分继续沿着传火管轴向流动。由于火药燃气到达每个传火孔的时间具有延迟性,因此靠近首孔位置的抛撒药首先被引燃,然后依次将靠近下一个传火孔处的抛撒药引燃,说明点火过程具有延迟性,延迟性的大小取决于传火管内点火药燃气的流动速度。

由图9(d)～9(f)可见,首孔位置附近区域的抛撒药已经全部达到着火温度,此处的抛撒药燃烧产生高温、高压燃气,高压的燃气将抛撒药向四周推散,抛撒药床与中心管内壁的空间属于“自由空间”,已燃的抛撒药与抛撒药燃气主要向自由空间流动,而部分抛撒药与抛撒药燃气向未燃烧的抛撒药床区域推进。由于中心管管壁附近存在自由空间,流入此处的已燃抛撒药与高温燃气快速地向中心管底部推进,从而加速引燃靠近管壁附近的抛撒药床。因此,在设计中心炸管式抛撒系统装药结构时,抛撒药与管壁之间要预留一定自由空间,以利于提高抛撒药引燃的速度。

从图9(g)～9(h)可以看出,末端传火孔处的抛撒药已燃区域要大于前面几个传火孔处的抛撒药已燃区域,其原因是,在传火管内高速流动燃气到达末端后,撞击壁面后经过反射,使末端传火孔处的压力以及能量迅速增大,导致末端传火孔处的抛撒药燃烧区域要大于前面几个传火孔附近抛撒药已燃区域。在2.69～3.5 ms时间段内,在末尾11个传火孔处的抛撒药床已燃区域是从上至下扩大的,说明该处抛撒药床主要是被管壁附近的抛撒药燃气所引燃,弥补了传火能量在此处能量的不足。

4 结论

本文给出了三维内弹道两相流的数值方法实现过程,并基于OpenFOAM开发了三维内弹道两相流求解器,最后进行了求解器的正确性验证,得出以下结论:

①该求解器可以针对不同结构燃烧室内弹道两相流问题进行仿真分析,减少了对工程问题仿真分析周期,借助OpenFOAM的前处理能力,所开发的内弹道两相流求解器可以分析复杂的三维计算域火药燃烧流动问题。

②内弹道两相流求解器具有足够的精度,能够精确捕捉压力波传播过程中的波形与位置,无需加入人工黏性项,对压力波捕捉具有高分辨率。

③内弹道两相流求解器具有全局守恒性:火药燃气减少的总质量等于火药燃气增加的总质量,火药燃烧产生的总能量等于全部内能和动能的总和。

④基于OpenFOAM的内弹道两相流求解器能够准确地模拟中心炸管式抛撒系统实际点传火过程,计算程序能为类似结构的抛撒装置结构优化提供良好的理论依据。