Numerical calculation of early fireball radiation spectrum in strong explosion
-
摘要: 基于强爆炸火球光辐射的多群辐射流体力学方法, 采用算子分裂方法将方程组分裂为对流项和刚性源项, 其中源项部分根据方程形式, 进一步分裂为各群内的单独求解。数值计算表明:该方法克服了直接求解过程中辐射与流体耦合所带来的强不稳定性, 时间步长大幅提高, 给出的火球光辐射能谱特征与已有规律一致。可为定量分析光辐射能谱特征提供有效手段。Abstract: On the basis of multi-group radiation hydrodynamics method of fireball radiation in strong explosion, operator splitting method is used to split the equations into convection items and source items, which are split into radiation groups due to the equation formation and solved individually.Numerical calculations show that the method used here overcomes strong instability when solving the equations directly because of the coupling items between radiation and fluid.In the meantime, the time step in the calculation is increased obviously.Fireball radiation spectrum is obtained in fine accordance with the result in the literature.
-
在强爆炸过程中, 由于强爆炸所释放的巨大能量迅速加热周围冷空气, 从而形成高温高压的火球。火球在发展过程中, 不断向外辐射光和热, 称为光辐射(也称为热辐射)。光辐射是强爆炸的重要毁伤效应之一, 研究其能谱特征对于深入分析光辐射与物质相互作用、评估相应的毁伤效应具有重要意义[1]。
强爆炸火球光辐射是一个涵盖近紫外到可见光直至远红外的宽谱带, 同时伴随火球的发展变化过程呈现不同的特征。早期的研究是基于实验测量, 结合相应的理论分析, 给出了光辐射各谱段较为接近的经验公式[1-2]。随着相关理论及数值模拟研究的发展, 特别是对火球发展过程的深入分析, 较普遍认为采用灰体近似下的辐射流体力学方法可以较好的描述火球发展过程[3-7], 并给出了光辐射总能量的时空分布规律[8-11]。但是灰体近似在处理光辐射输运过程中, 假定辐射系数(空气吸收系数)与光子能量无关, 采用全谱范围内的平均计算, 因而无法给出光辐射能谱特征。
为能够给出光辐射能谱, 需采用多群方法, 即在不同能谱范围内求解光辐射输运过程, 分群越多, 计算得到光辐射能谱特征越精细[12]。显然, 这种处理方法的代价使得辐射方程的计算更为复杂, 需要辐射参数(分群吸收系数)更多, 与流体耦合后求解, 更容易由于方程的强非线性、强刚性[2](尤其是辐流耦合项)而出现计算不稳定(除非将时间步长限制非常小)。国外有关文献[4, 12]仅给出了部分结果, 但并未涉及方程求解中的具体处理过程。国内光辐射能谱方面的计算较少, 早期计算[13]中利用小步长(约10-15s)直接求解的方法, 给出了光辐射较短时间内的能量分布情况。
本文中基于辐射流体力学方法, 采用多群方法求解光辐射输运过程。方程求解过程中, 为克服方程的强非线性、辐流耦合项的强刚性可能产生的计算不稳定问题, 通过算子分裂方法, 将多群方程逐步分裂为对流项和源项单独求解[14-15], 其中源项部分的求解需要在各群内进一步分裂, 使之具有单群方程的形式。在这样的处理方法下, 多群辐流方程的求解比直接求解更稳定, 即使在大时间步长(约10-11s)下也能得到较好的结果, 有效提高了计算效率, 为大规模问题求解提供了条件。利用该计算方法, 计算得到了强爆炸早期光辐射能谱特征分布, 与经验关系给出的规律较为一致, 为定量分析光辐射能谱特征提供了有效手段。
1. 计算方程及参数
在局域热动力平衡(LTE)假定下, 描述强爆炸火球光辐射发展的多群辐射流体力学方程如下:
∂ρ∂t+∇⋅(ρv)=0 (1) ∂∂t(ρv)+∇p=−M(1) (2) ∂∂t(12ρv2+ρe)+∇⋅[(12ρv2+ρe+p)v]=−Q(0) (3) 式(1)~(3)分别为质量守恒、动量守恒和能量守恒方程。式中:ρ为空气密度,v为空气速度,p为空气压力,e为空气内能。M(1)为各群内光辐射输运的总动量, M(1)=∑gM(1)g=∑g1c2∂Fg∂t+∇⋅Pg;Q(0)为各群内光辐射输运的总能量, Q(0)=∑gQ(0)g=∑g∂Eg∂t+∇⋅Fg,上标“0”和“1”表示辐射输运方程的0阶矩和1阶矩方程, 下标g=1, 2, 3…即分群数目,Eg、Fg、Pg、Qg(0)、Mg(1)分别为第g群辐射能密度、辐射能流、辐射压力张量及光辐射输运能量和光辐射输运动量。具体计算如下:
∂Eg∂t+∇⋅Fg=cκg∫γgγg−1[4πcB(γ,T)−E(γ)]dγ+vcκg∫γgγg−1F(γ)dγ=cκg(Bg−Eg)+κgcv⋅Fg(4) (4) 1c2∂Fg∂t+∇⋅Pg=−1cκg∫γgγg−1F(γ)dγ+4πc2κgv∫γgγg−1B(γ,T)dγ+vcκg∫γgγg−1P(γ)dγ=κgc(−Fg+vBg+v⋅Pg) (5) 式中:c为光速, κg为第g群吸收系数, γ为光子频率。由于各群光子能量具有上限和下限, 因此采用γg、γg-1分别标记第g群光子频率的上限和下限。E (γ), F(γ)及P (γ)分别为辐射能密度、辐射能流及辐射压强张量函数, 与Eg、Fg、Pg关系为 Eg=∫γgγg−1E(γ)dγ,Fg=∫γgγg−1F(γ)dγ,Pg=∫γgγg−1P(γ)dγ; B(γ, T)为黑体辐射谱分布, h为普朗克常数, k为波尔兹曼常数, Bg为黑体辐射下第g群辐射能密度:
Bg=4πc∫γgγg−1B(γ,T)dγ=4πc∫γgγg−12hγ3c21ehγ/kT−1 dγ (6) 吸收系数采用文献[17-18]给出的21群Rosseland分群吸收系数κg。多群方法下, 吸收系数为空气温度、密度以及光子分群能量的函数。数值求解中, 通过双线性插值计算相应的温度、密度下, 空气的分群吸收系数。图 1给出空气密度为
, 分群吸收系数与光子能量(分群)和温度的关系。
对于各群内分群辐射压强张量Pg, 采用最大熵变Eddington因子近似[13, 16]进行计算, 空气的状态方程采用实际空气状态方程[9-11]。光辐射按照对应的光子能量pe分为21群[13, 17-18], 如表 1所示。
表 1 光子分群能量Table 1. Photon energy of each groupg pe/eV 1 0.01~0.5 2 0.5~1.0 3 1.0~1.8 4 1.8~2.1 5 2.1~2.5 6 2.5~3.1 7 3.1~4.0 8 4.0~7.0 9 7.0~10.0 10 10.0~20.0 11 20.0~40.0 12 40.0~70.0 13 40~70 14 100~200 15 200~400 16 400~1 000 17 1 000~2 000 18 2 000~5 000 19 5 000~10 000 20 10 000~20 000 21 20 000~80 000 - - - - - - 2. 数值求解
2.1 分裂方法过程
一维球对称下多群辐流方程组展开后也可写成如下矩阵形式:
∂f∂t+∂g∂r+Ψ+Φ=0 (7) 各项代表的物理量如下:
f=(ρρveE1F1c2⋮EgFgc2⋮),g=(ρvp+ρv2(e+p)vF1P1⋮FgPg⋮),Ψ=(2ρvr2ρv2r(e+p)vr2F1r3P1−E1r⋮3Pg−Egr⋮) Φ=(0κ1c(vP1+vB1−F1)+⋯+κgc(vPg+vBg−Fg)+⋯[cκ1(B1−E1)+κ1cvF1]+⋯+[cκg(Bg−Eg)+κgcvFg]+⋯cκ1(B1−E1)−κ1cvF1κ1c(F1−vP1−vB1)⋮cκg(Bg−Eg)−κgcvFgκgc(Fg−vPg−vBg)⋮) 式中:Ψ为对流项, Φ为辐射与流体耦合项, 即刚性源项部分。利用算子分裂方法, 将方程分2步求解:
{∂f(1)∂t+∂g∂r+Ψ=0∂f(2)∂t+Φ=0(f(2)|t=t0=f(1)∗) (8) 第1步求解对流项, 采用有限体积法, 构造五阶WENO格式, 数值通量的计算采用局部Lax-Friedrichs方法, 由于不含刚性源项, 因而时间步长可以提高到10-11s, 相比于直接法求解[13], 极大提高了计算效率。第2步求解源项, 其初始时刻(t=t0)的值为第1步对流项方程的解f(1)*。源项求解较为复杂, 首先根据方程特征, 写为:
Φ=φ1+φ2+⋯+φg (9) φg=(0κgc(vPg+vBg−Fg)cκg(Bg−Eg)+κgcvFgδ1g[cκg(Bg−Eg)−κgcvFg]δ1gκgc(Fg−vPg−vBg)⋮δngcκg[(Bg−Eg)−κgcvFg]δngκgc(Fg−vPg−vBg)⋮) (10) 其中
。在这种情况下, 源项方程的求解, 可以继续按照步骤进行(带“*”的参量均为上步方程的解):
{∂f(2)1∂t+φ1=0(f(2)1|t=t0=f(1)∗)∂f(2)2∂t+φ2=0(f(2)2|t=t0=f(2)∗1)⋮∂f(2)g∂t+φg=0(f(2)g|t=t0=f(2)∗g−1) (10) 分裂成为这样的步骤进行求解, 方程形式大为简化, 其过程本身也具有自身的物理意义:每步求解认为流体仅与该群光子进行能量和动量交换, 二者组成的体系中, 总能量和总动量守恒。这样的假设, 在方程分裂的数学处理过程中是严格满足的:因为其他群内辐射能密度与辐射能流随时间的微分都等于0。进一步的求解, 则转化为常微分方程组的求解, 可以通过多种方法实现快速、高精度求解。
2.2初始条件和边界条件
求解初始条件假定爆炸总能量集中于等压火球内, 辐射能与空气内能总合等于爆炸总能量, 计算边界条件采用对称边界。光辐射分群能量的初始分布通过Bg给出, 分群辐射能流为0。
在数值求解中, 计算区域的边界处, 假定物理量都处于未扰动状态, 即认为:流体速度以及辐射能流均为0, 而其他状态参量取初始值。
3. 计算结果与分析
取当量为1 kt, 高度在海平面(h=0 km), 求解相应的强爆炸火球光辐射输运过程, 空气初始状态参量如表 2所示。
表 2 不同高度空气初始状态参数Table 2. Air state parameters at different altitudeh/km p/kPa ρ/(kg·m-3) T/k cs/(m·s-1) 0 103.3 1.225 297.25 343.6 10 26.5 0.414 225.90 299.5 20 5.5 0.089 219.20 295.1 强爆炸火球发展过程中, 冲击波阵面和辐射波阵面是一个描述火球发展的重要参量。利用上述算子分裂方法计算得到的强爆炸火球阵面走时, 如图 2所示, 符合火球发展的物理过程, 与H.L.Brode[4]计算的结果符合较好, 说明该方法在处理过程中是稳定可靠的。
多群方法的运用, 能够给出火球光辐射在特定时刻向外辐射的光辐射能谱特征。图 3所示为t=0.023 s时刻火球光辐射分群(21群)能量分布, I为光辐射强度。从图 3中可以看出, 火球光辐射能量集中于第2~8群, 对应光子波长在0.2 ~2μm。
根据火球发展过程, 在光辐射强度第1个极大值后, 火球有效温度从约20 000 K降低至3 000 K, 在光辐射强度第2个极大值时, 火球有效温度略低于10 000 K。在这个温度范围内, 火球光辐射能谱大部分能量都集中与紫外(0.22~ 0.36μm)、可见(0.36~0.64μm)和红外(0.64~ 4.5μm)部分, 与图 3所示基本一致。
为具体分析火球光辐射能谱特征, 利用文献[1-2]给出的火球有效温度走时关系:
4πr2σT4e=Φe 式中:Te为火球有效温度, σ为斯特潘-玻尔兹曼常数, Φe为与火球辐射功率有关的函数表达式, r为火球半径, 各个参数可以通过已有规律计算给出。以t=0.01、0.02 s时刻为例, 利用上述理论方法计算得到该爆炸条件下火球有效温度约为3 500K和7 453 K, 据此给出0.01~2.5μm波长范围内光辐射强度随波长关系与对应有效温度下的黑体谱对比分析, 如图 4所示, 图中实线为本文中方法的计算结果, 虚线为由文献理论计算的黑体谱分布。
已有研究结果表明[1-2], 在火球光辐射强度第1个极大值后的整个发光阶段, 波长在(0.4~ 0.6μm)范围内, 可以近似看做黑体, 与图 4给出的结果在变化规律上基本一致, 数值大小上的差异, 与所采用的空气吸收系数等有关, 需要进一步工作进行细致改进。
4. 结论
(1) 基于强爆炸火球的多群辐射流体力学模型, 采用算子分裂方法对方程组进行数值求解。利用空气分群(21群)吸收系数, 计算给出了1 k T当量下火球光辐射能谱特征。数值计算结果验证了光辐射能谱主要集中在0.2~2μm(紫外、可见到红外波段), 与已有结果和经验规律符合的较为一致;
(2) 分裂求解的处理方法, 一方面克服了直接求解过程中辐射与流体耦合可能带来的强不稳定性, 另一方面扩大了时间步长, 提高了计算效率, 为类似方程的数值求解提供了一定借鉴。
-
表 1 光子分群能量
Table 1. Photon energy of each group
g pe/eV 1 0.01~0.5 2 0.5~1.0 3 1.0~1.8 4 1.8~2.1 5 2.1~2.5 6 2.5~3.1 7 3.1~4.0 8 4.0~7.0 9 7.0~10.0 10 10.0~20.0 11 20.0~40.0 12 40.0~70.0 13 40~70 14 100~200 15 200~400 16 400~1 000 17 1 000~2 000 18 2 000~5 000 19 5 000~10 000 20 10 000~20 000 21 20 000~80 000 - - - - - - 表 2 不同高度空气初始状态参数
Table 2. Air state parameters at different altitude
h/km p/kPa ρ/(kg·m-3) T/k cs/(m·s-1) 0 103.3 1.225 297.25 343.6 10 26.5 0.414 225.90 299.5 20 5.5 0.089 219.20 295.1 -
[1] 乔登江.核爆炸物理概论[M].北京: 国防工业出版社, 2003: 225-262. [2] 屠琴芬.核爆炸火球的辐射流体力学计算中的几个问题[R].西安: 西北核技术研究所, 1986. [3] Pomraning G C. The equations of radiation hydrodynamics[J]. Astrophysical Radiation Hydrodynamics, 1986(188): 45-69. [4] Brode H L. Fireball phenomenology[R]. AD0612197, 1964. [5] Crowley B K, Glenn H D, Marks R E. An analysis of marvel: A nuclear shock-tube experiment[J]. Journal of Geophysical Research, 1971, 76(14): 3356-3374. doi: 10.1029/JB076i014p03356 [6] Marrs R E, Moss W C, Whitlock B. Thermal radiation from nuclear detonations in urban environments[R]. UCRL-TR-231593, 2007. [7] Lowrie R B, Edwards J D. Radiative shock solutions with grey nonequilibrium diffusion[J]. Shock Waves, 2008, 18(2): 129-143. [8] 陈健华, 王心正, 谢龙生, 等.均匀大气中的强爆炸一维辐射流体力学数值解[J].爆炸与冲击, 1981, 1(2): 37-49.Chen Jian-hua, Wang Xin-zheng, Xie Long-sheng, et al. A one-dimensional radiation hydrodynamic numerical solution for a strong explosion in uniform atmosphere[J]. Explosion and Shock Waves, 1981, 1(2): 37-49. [9] 田宙, 乔登江, 郭永辉.不同当量强爆炸早期火球现象的数值模拟[J].爆炸与冲击, 2009, 29(4): 408-412.Tian Zhou, Qiao Deng-jiang, Guo Yong-hui. Numerical simulation on early fireball phenomenology of strong explosions for different yields[J]. Explosion and Shock Waves, 2009, 29(4): 408-412. [10] 田宙, 乔登江, 郭永辉.不同高度强爆炸早期火球数值研究[J].兵工学报, 2009, 30(8): 1078-1083.Tian Zhou, Qiao Deng-jiang, Guo Yong-hui. Numerical investigation of early fireball of strong explosion for different altitudes[J]. Acta Armamentarii, 2009, 30(8): 1078-1083. [11] 田宙, 乔登江, 郭永辉.强爆炸早期火球现象的一维数值研究[J].计算物理, 2010, 27(1): 9-14.Tian Zhou, Qiao Deng-jiang, Guo Yong-hui. A one dimensional numerical study on early fireball in strong explosion[J]. Chinese Journal of Computational Physics, 2010, 27(1): 9-14. [12] Symbalisty E M D, Zinn J, Whitaker R W. Radflo physics and algorithms[R]. LA-12988-MS, 1995. [13] 高银军, 田宙, 刘峰, 等.强爆炸早期火球光辐射能谱的分群计算[J].四川兵工学报, 2011, 32(3): 21-24.Gao Yin-jun, Tian Zhou, Liu Feng, et al. Calculation of energy spectrum of early fireball radiation in strong explosion with multi-group method[J]. Journal of Sichuan Ordnance, 2011, 32(3): 21-24. [14] 闫凯.二维辐射流体动力学方程组的数值求解[D].西安: 西北核技术研究所, 2011. [15] 闫凯.二维辐射流体动力学方程组的数值求解[C]//第六届全国青年计算物理学术会议.太原, 2011. [16] Minerbo G N. Maximum entropy eddington factors[J]. Journal of Quantitative Spectroscopy and Radiative Transfer, 1977, 20(6): 541-545. [17] 王文高, 张建泉.物质辐射不透明性Ⅰ[R].西安: 西北核技术研究所, 1978. [18] 王文高, 张建泉.物质辐射不透明性Ⅱ[R].西安: 西北核技术研究所, 1978. 期刊类型引用(5)
1. 王建国. 高空核爆炸磁流体动力学电磁脉冲. 强激光与粒子束. 2024(07): 18-24 . 百度学术
2. 王建国,刘利,牛胜利,左应红,高银军,朱金辉,张相华,李桠,李夏至. 高空核爆炸环境数值模拟. 现代应用物理. 2023(01): 3-14 . 百度学术
3. 李康,李守先,刘娜. 强爆炸火球辐射流体自适应网格高精度数值模拟. 计算物理. 2021(02): 146-152 . 百度学术
4. 高银军,高丽红,张相华,马壮,刘峰,彭国良,田宙. 强爆炸光辐射作用下材料的能量耦合特性. 中国光学. 2020(06): 1267-1275 . 百度学术
5. 高银军,田宙,闫凯,刘峰. 强爆炸光辐射脉冲辐照特征与爆炸当量的相关性. 爆炸与冲击. 2017(03): 549-553 . 本站查看
其他类型引用(0)
-