Numerical simulation of influence of target dimension and laser spot size on laser impulse
-
摘要: 采用二维雷诺平均N-S方程,数值模拟研究了大气条件下短脉冲激光与固体靶相互作用所产生等离子体的动力学过程。采用k-ε两方程模型用于湍流的数值模拟,分别利用ROE格式和二阶中心格式对对流通量和粘性通量进行离散处理;用高斯-赛德尔隐式格式对方程进行时间推进求解。数值模拟给出了激光引发靶蒸气等离子体侧向膨胀、稀疏等二维流体动力学过程的物理图像,讨论了靶与光斑尺寸对脉冲激光冲量的影响。结果表明,不同宽度固体靶受到的激光冲量有很大差异,固体靶宽度越大,受到的激光冲量也越大。Abstract: The plasma flow field of solid target in atmosphere under impulse laser irradiation was computed. The flow is governed by the 2-D Reynolds averaged Navier-Stokes (NS) equations. The k-ε turbulence model was used for turbulence simulations. To split the viscosity flux and the convective flux of the NS equations, the second order central scheme and the ROE scheme were adopted respectively. With the implicit Gauss-Seidel scheme, the code was advanced in time. The 2-D hydrokinetics process of laser plasma was presented in this paper by numerical simulation. The influence of target dimension and laser spot size on laser impulse was discussed. The results indicated that laser impulse increases generally with target dimension.
-
Key words:
- numerical simulation /
- laser impulse /
- target /
- laser spot size
-
蒙特卡罗(简称“蒙卡”)方法在核反应堆堆芯设计与分析中起着越来越重要的作用。鉴于此方法能够对堆芯几何结构进行精确描述并且能够使用连续能量截面的优势,基于该方法的程序经常被用于三维全堆输运计算并给出参考结果。蒙特卡罗算法的优点亦可以用到现有的两步法堆芯计算中(燃料组件计算+堆芯计算),而少群参数的制作是两步法计算中的关键。目前,针对蒙卡方法计算截面参数的研究仍然比较受限,大部分研究还是集中在适用于扩散程序的少群均匀化参数上。
在组件计算时,边界条件经常设为全反射或者周期性边界条件。热中子的平均自由程在1~2 cm之间,基本可满足全反射边界条件的假设。相反,快中子的平均自由程远超过单个燃料组件的尺寸。相邻组件的联系在快中子之间更为紧密,故而全反射边界条件对其并不适用。因此,在考虑快中子堆组件计算时需要考虑泄漏效应。考虑泄漏效应的少群常数制作可以通过在计算中嵌入泄漏修正模型来实现,此方法期望通过理论模拟在即使不清楚反应堆堆芯真实运行条件下也可以构造出组件临界状态。基于蒙卡方法的组件泄漏模型大致分为三种。方法一来自于确定论中常用的基模修正法,该方法首次用于韩国KAERI的蒙卡程序McCARD中[1],芬兰VTT实验室开发的蒙卡燃耗程序SERPENT[2]也使用此方法做泄漏修正。第二种泄漏修正方法来自于韩国KAIST实验室,这种方法的主要思想是通过改变边界反照率来搜索组件临界状态[3],此方法已经在MCNP5中实现并且与MONTEBURNS程序耦合做燃耗计算。第三种泄漏模型由日本京都大学的Yamamoto教授提出,通过给中子权重赋予复数形式,新颖地将中子泄漏项与其传播方向结合起来[4]。与前两种方法相比,第三种方法能够给出每个中子的泄漏概率,更合理地应用于蒙卡程序的模拟过程中。尽管此方法占据了理论上的优势,但它的几何适用范围受限,为了消除裂变源项中的虚数部分,该方法只适用于对称的几何结构上。
本文介绍了法国蒙卡程序TRIPOLI4中嵌入的泄漏修正模型。此泄漏模型的开发主要用于少群截面参数的制作[5]。目前,TRIPOLI4程序中选择的泄漏修正模型为基模修正方法。另外,基于蒙卡方法制得的少群截面参数的验证也是本文的一个重点内容。本文设计了一套自然泄漏的数值计算模型用以验证TRIPOLI4程序中少群截面参数的制作功能。此外,该程序产生的截面参数也将与已经相对成熟的快堆组件程序ECCO[6-7]和另一套蒙卡程序SERPENT进行对比。最后,这些少群截面参数还将被用于堆芯计算中并分析比较其计算结果。
1. 模型与算法
主要叙述基模修正模型及其在程序中使用的具体算法,包括详细的针对基模方程的物理近似,以及对此模型在程序中的流程嵌入和代码实现。
1.1 B1均匀化方程
输运方程在均匀无穷介质结构下的解即通常所称的基态模式。组件计算一般在全反射或者周期性边界条件下进行。因而,泄漏模型必须要加进燃料组件计算中才能强制有效增殖系数Keff=1。基模的解通常由两部分组成,第一部分是基态通解φ(r, E, Ω)对于任意位置的组件通解保持不变;第二部分f(r)被用于堆芯中的通量幅值形状,并且,此形状函数f(r)是拉普拉斯方程(1)的解。
∇2f(r)+B2f(r)=0 (1) 式中:B2为堆芯的几何曲率。该方程对应的通解为f(r)=αeiB·r。对于每个独立组件,基态通解已经不再随空间变化。因此,组件内的角通量密度可写为
φ(r,E,Ω)=φ(E,Ω)eiB⋅r (2) 式中:B代表曲率向量,其值满足B·B=B2。
把式(2)变量分离后的通量密度表达式代入连续能量的Boltzmann方程中,并采用线性各向异性近似假设,即可得到均匀化的B1基模方程
[iΩ⋅B+Σt(E)]φ(r,E)=∫∞014π{Σs0(E′→E)φ(E′)+3Σs1(E′→E)Ω⋅J(E)}dE′+14πKeff∫∞0χ(E′→E)νΣf(E′)φ(E)dE′ (3) 式中:Σt,Σs0,Σs1分别代表总截面、零阶散射截面、一阶散射截面;ν为中子裂变增殖系数;χ为中子裂变谱;J (E)代表基态中子流,表示为
J(E)=∫4πΩφ(E,Ω)d2Ω (4) 为了消除曲率和中子流的方向性,假设J和B具有相同的方向,故而有J(E)·B =J(E)B。
通过对方程(3)进行角度积分以及乘以1iΩ⋅B+∑t(E)后再进行角度积分,可分别得到下面两组与标通量密度相关的方程
iBJ(E)+Σt(E)φ(E)=∫∞0Σs0(E′→E)φ(E′)dE′+1Keff∫∞0χ(E′→E)νΣf(E′)φ(E′)dE′ (5) iJ(E)B=1γ[B,Σt(E)]Σt(E){13φ(E)+∫∞0Σs1(E′→E)iJ(E′)B dE′} (6) 式中:γ[B, Σt(E)]为新引入的变量,其定义如下
γ[B,Σt(E)]={13B∑t(E)arctanB∑t(E)1−B∑t(E)arctanB∑t(E),B2<0B23∑t(E)ln∑t(E)+Im(B)∑t(E)−Im(B)2Im(B)−∑t(E)ln∑t(E)+Im(B)∑t(E)−Im(B),B2>01+415(B∑t(E))2−12175(B∑t(B))4+⋯,B2>0 (7) 需要特别指出的是,当曲率接近为零时,iJ(E)/B应该一直保持为实数并且有限。实际上,根据所在介质的不同,J(E)可以为虚数或为实数。当堆芯处于超临界状态时,J(E)为虚数并且B为正数;反之,当堆芯处于次临界状态时,J(E)为实数而且B为虚数。
于是,通过菲克定律,可以得到扩散系数为
d(B,E)=1BiJ(E)φ(E) (8) 将式(8)分别代入方程(5)与方程(6)中,可得到如下变形后的B1方程组
[Σt(E)+d(B,E)B2]φ(E)=∫∞0Σs0φ(E′)dE′+1Keff∫∞0χ(E′→E)νΣf(E′)φ(E′)dE′ (9) d(B,E)=1γ[B,Σt(E)]Σt(E){13+∫∞0Σs1(E′→E)d(B,E′)φ(E′)φ(E)dE′} (10) 有了上述连续能量的B1方程组之后,可以很容易通过并群得到其少群形式下的基态方程组
Σt,gφg+dgB2φg=∑Gg′=1∑g′→gs0φ′g+1Keff ∑Gg′=1χ(g′→g)νΣf,g′φ′g (11) dg=1∑t,gγg{13+∑Gg′=1∑g′→gs1φ′gφg d′g} (12) 式中:φg为g群内通量密度,φg=∫E∈gφ(E)dE;dg为g群内扩散系数,dg=1φg∫E∈gd(B,E)φ(E)dE;γg为系数,可简化为近似表达形式γg≈γ(B, Σt, g);Σi, g为一维少群截面,例如总截面、吸收截面、裂变截面等,∑i,g=∫E∈g∑i(E)φ(E)dE∫E∈gφ(E)dE;∑g′→gs0为多群散射截面的零阶项,∑g′→gs0=(Σt,g′−Σa,g′)∗P(g′→g),其中,P(g′→g)为群间散射概率;∑g′→gs1为多群散射截面的一阶项,∑g′→gs1=∑g′→gs0ˉμg′→g,其中,ˉμg′→g代表从g′到g群的平均散射角余弦。
值得指出的是,少群截面参数是由连续能量的TRIPOLI4模拟计算提供,而方程组(11)~(12)中未知量为dg,φg和γg,其具体求解方法将在下一节详细介绍。
1.2 B1算法流程
本节给出TRIPOLI4程序中求解B1均匀化方程的计算流程,其求解过程由内迭代和外迭代两部分交替完成。内迭代过程即在固定B2值下求出收敛后的通量密度以及Keff值;外迭代过程即通过方程(12)变化B2值来使得Keff=1。
B2n+1=B2n+(1Kexp eff −1Kneff )¯ν∑nfˉdn (13) 式中:n为内迭代小标数;Keffexp为用户设定的Keff=1值;νΣfn为第n次内迭代后,由通量密度权重得到的平均裂变产生截面;dn为第n次内迭代后,由通量密度权重得到的扩散系数。
图 1给出了完整的迭代流程。图中左边部分代表连续能量截面的蒙卡模拟过程,右边部分代表求解B1均匀化方程的过程。左右两部分通过两条线连接。其一为收敛后的少群截面参数,将作为已知量被用于求解B1均匀化方程中;其二为由求解B1方程后得到的泄漏截面dgB2,此泄漏系数将被代入蒙卡的连续截面模拟过程中。由于该泄漏修正模型涉及到基理完全不同的蒙卡程序和确定论程序的耦合,因此在收敛准则上需要注意匹配。此处主要涉及到三种收敛准则。第一,用户根据需求可以给出蒙卡模拟过程中的总粒子数(批数量×每批粒子数)以确保每个能群下计算得到的截面参数都是收敛值;第二,求解B1均匀化方程时对Keff的收敛标准也控制得比较严格,一般在|Keff-1| < 1×10-5;第三,最终的蒙卡模拟得到的Keff值与用户设定的期望值Keffexp可保持在如下收敛范围内:|Keff-Keffexp| < 100×10-5+3σ。公式中σ为蒙卡模拟时给出的自身不确定性。
2. 验证
介绍了TRIPOLI4程序中泄漏修正模块的验证工作。首先,通过蒙卡连续能量模拟得到的一个一维临界几何结构被设计出来用于给出各项参考量,例如临界通量密度谱、反应率以及少群截面参数等。与TRIPOLI4程序中泄漏模块进行对比的有蒙卡程序SERPENT以及已经较成熟的用于快堆截面制作的确定论程序ECCO。本验证工作选择了快堆组件Zero Power Plutonium Reactor Double Column Fuel(ZPPRDCF)的材料成分作为一维临界几何模型来展开计算以及分析对比工作。此组件模型的Kinf高达1.67,也就意味着此临界模型的中子泄漏率达40%。这种高泄漏结构也就对泄漏模型的处理要求更加严格。另外,为了使各种程序的计算结果有可比性,各软件统一使用JEFF3.1.1作为输入截面库。下面将分别给出计算结果参数的对比分析。
2.1 临界几何曲率B2
现今包含泄漏模型的组件计算程序给出几何曲率已经不再困难,而如何找到一套可靠的参考值却是值得推敲的问题。其实,该问题可以转移到如何确定方程(1)中空间分布函数f(r)。多数情况下,f(r)被认为是矢量通量密度分布φ(r)。分布函数f(r)实际上是两个不同通量密度的内积函数
f(r)=∫∞0φ+(E)φ(r,E)dE (14) 式中:φ+ (E)为基态模式下的共轭通量密度;φ(r, E)为堆芯计算下的直接通量密度。此两者通量密度皆由TRIPOLI4连续蒙卡模拟计算给出。因而,在已知f(r)轴向分布的情况下,可以拟合出几何曲率B2。
表 1给出了各种不同程序计算得到的几何曲率值以及他们的比较,其中,由通量密度内积拟合出的曲率值作为参考结果。可看出,TRIPOLI4的计算结果与参考值最为接近,相差1.09%,而其他两套程序的计算结果与参考值相比偏差稍大。此结果可归功于TRIPOLI4的泄漏修正模块中通过连续截面的蒙卡模拟过程与B1方程求解的交叉迭代直至达到真正临界状态,而这样的计算模型在其他两套程序中是不具备的。
表 1 临界曲率比较结果Table 1. Comparison of critical buckling valuesB2/(10-3cm-2) relative diff/% numerical curve 2.854 - TRIPOLI4 2.823 -1.09 SERPENT 2.807 -1.65 ECCO 2.747 -3.75 2.2 归一化临界通量密度谱
图 2给出了多种包含泄漏修正模型程序的归一化临界通量密度计算结果。其中的实线(Numerical Exp)代表由连续能量截面的TRIPOLI4模拟给出的参考结果。图 2(a)给出了各程序计算得到的归一化通量密度值,其中归一化原则为∑Gg=1φg=1,能群总数G可取33。可以看出,这三套程序中的泄漏修正模型都能够给出与参考值趋势一致的结论。图中灰色区域为重要能群区,故而将此区域数据单独描在图 2(b)中。此外,三套程序与参考值的相对误差也在图 2(c)中给出。可以看出,在高能区部分,TRIPOLI4,SERPENT以及ECCO程序中的泄漏修正模型均低估了临界通量密度水平。相反,在较低能量区域,三套程序中的泄漏修正模型却高估了通量密度水平。
2.3 泄漏率
本节给出了各程序中泄漏修正模型模拟得到的泄漏率(dgB2φg)及其与参考值的对比。参考值依然来自于TRIPOLI4连续能量截面模拟得到的中子泄漏率。图 3中描述了各软件计算得到的泄漏率与参考值的偏差。可以看到,这三套程序模拟的泄漏率与参考值相比在每个能群内偏差都控制在500×10-5以内。此外,这三套程序中都存在高能群区与低能群区的互补现象,这一现象主要由基态模型的近似造成。在所使用的基态模型中,认为中子谱不随组件在堆芯中的位置变化。实际上,在一个真实堆芯中,堆芯外缘的谱要比堆芯中部的谱稍硬,因而在堆芯外缘处,快中子泄漏的机率偏小,这也就解释了图 1中在2 MeV处中子通量密度估算值偏高的现象。
3. 结论
本文介绍了法国蒙卡程序TRIPOLI4中嵌入的泄漏修正模型以及利用此模型制作适用于快堆的少群截面参数的工作。此泄漏模型利用通过蒙卡产生的少群参数作为B1均匀化方程的已知量从而求得相应的少群泄漏参数。该泄漏模型的优点在于TRIPOLI4程序能够利用以上解得的泄漏系数在连续蒙卡模拟过程中直至迭代成功到真正临界状态。本文也给出了关于此泄漏模型的基础验证工作。通过对比几何曲率、临界中子谱、泄漏率等参数,可以看到TRIPOLI4程序中的泄漏修正模型能够给出参考值相符合的结果,甚至比同类程序(SERPENT,ECCO)的计算精度更高。
-
[1] 孙承纬. 激光辐照效应[M]. 北京: 国防工业出版社, 2002.Sun Chengwei. Laser irradiation effects. Beijing: National Defense industry Press, 2002 [2] Phipps C R, Turner T P, Harrison R F, et al. Impulse coupling to targets in vacuum by KrF, HF, and CO2, single-pulse lasers[J]. Journal of Applied Physics, 1988, 64(3): 1083-1096. doi: 10.1063/1.341867 [3] Nehls M, Edwards D, Gray P. Ablative laser propulsion using multi-layered material systems[C]//Proc of 33rd AIAA Plasma Dynamics and Lasers Conference. 2002: 20-23. [4] Bass M, Nassar M A, Swimm R T. Impulse coupling to aluminum resulting from ND: glass laser irradiation induced material removal[J]. J Appl Phys, 1987, 61(3): 1137-1144. [5] Pakhmov A V, Lin J, Kenneth A H. Effect of air pressure on propulsion with TEA CO2[C]//Proc of SPIE. 2004, 5448: 1017-1027. [6] 章玉珠, 王广安, 沈中华, 等. 等离子体屏蔽和稀疏波对冲量耦合系数的影响[J]. 强激光与粒子束, 2007, 19(4): 585-588. http://www.hplpb.com.cn/article/id/3146Zhang Yuzhu, Wang Guang'an, Shen Zhonghua, et al. Influences of plasma shielding and rarefaction wave on impulse coupling coefficient. High Power Laser and Particle Beams, 2007, 19(4): 585-588 http://www.hplpb.com.cn/article/id/3146 [7] Phipps C, Luke J, Funk D, et al. Laser impulse coupling at 130 fs[J]. Applied Surface Science, 2006, 252(13): 4838-4844. doi: 10.1016/j.apsusc.2005.07.079 [8] 英敏菊, 王潇潇, 程伟, 等. 激光烧蚀靶材蒸发波模型的气体动力学与冲量耦合系数计算[J]. 北京师范大学学报(自然科学版), 2015, 51(1): 40-44. https://www.cnki.com.cn/Article/CJFDTOTAL-BSDZ201501011.htmYing Minju, Wang Xiaoxiao, Cheng Wei, et al. Calculation of gas dynamic parameters and impulse-coupling coefficient for laser ablation of a target. Journal of Beijing Normal University(Natural Science), 2015, 51(1): 40-44 https://www.cnki.com.cn/Article/CJFDTOTAL-BSDZ201501011.htm [9] 张黎, 贺佳, 谭福利. 激光加热金属板流固耦合数值模拟[J]. 强激光与粒子束, 2011, 23(4): 866-869. http://www.hplpb.com.cn/article/id/5128Zhang Li, He Jia, Tan Fuli. Numerical simulation of metal plates under laser irradiation based on fluid-solid coupling. High Power Laser and Particle Beams, 2011, 23(4): 866-869 http://www.hplpb.com.cn/article/id/5128 期刊类型引用(2)
1. 吴宏春,杨红义,郑友琦,曹良志,杜夏楠,杨勇,刘一哲,胡赟. 快中子反应堆堆芯物理分析方法的研究现状与发展建议. 原子能科学技术. 2024(03): 513-527+507 . 百度学术
2. 李金洲,张滕飞,何东豪,潘清泉,刘晓晶. MORPHY程序在铅冷快堆中的应用. 强激光与粒子束. 2024(06): 142-153 . 本站查看
其他类型引用(0)
-