摘要： 针对激光辐照热效应模拟中传统的有限元法对复杂几何结构和多次散射光处理能力较差的问题，在三维多物理场流固耦合问题并行计算程序CFS中实现射线追踪有限元模拟方法.与传统有限元法相比，该方法在有限元网格划分的基础上，采用射线追踪技术将各网格单元所吸收的辐照激光能量加入方程求解系统.对2A12硬铝材料腔体结构件激光辐照热效应的计算表明，该方法提高程序对计算模型几何结构复杂性的适应能力，解决多次散射光存在时传统方法难以准确描述辐照激光能量在求解区域中分布的问题，从而有效扩展CFS程序的使用范围.
　　关键词： 射线追踪； 激光辐照； 热效应； 多次散射光； 有限元； CFS
　　中图分类号： TN244；TB115.1 文献标志码： B
　　Hybrid simulation method for thermal effects of
　　laser irradiation
　　ZHANG Xianghua1, SHU Qingbang1, DU Taijiao1,2, WANG Yuheng1,2, CHEN Zhihua1,2
　　（ 1. Northwest Institute of Nuclear Technology, Xi’an 710024, China;2. State Key Laboratory of Laser Interaction with Matter, Xi’an 710024, China）
　　Abstract： As to the problem that the handling capacity on complex geometrical structure and repetitious scattering rays is poor in the simulation on thermal effects of laser irradiation using traditional finite element method, ray tracing vs finite element simulation method is implemented in the parallel computing program CFS which is for solving fluidsolid coupling problem in 3D multiphysical fields. Compared with the traditional finite element method, based on the finite element meshing and ray tracing technology, the system is solved by adding laser irradiation energy absorbed by elements in equations. The calculation of thermal effects of laser irradiation on cavity structure which is made of hard Aluminum 2Al2 shows that the method improves the capacity of computational model program in adopting for geometrical structure complexity, and also can solve the problem that it is difficult to describe the laser irradiation energy distribution accurately in the solving area using traditional method while repetitious scattering rays exist. Then the use of program CFS can be extended efficiently.
　　Key words： ray tracing; laser irradiation; thermal effect; repetitious scattering ray; finite element; CFS
　　0 引 言
　　激光辐照热效应是激光与物质相互作用的重要研究内容之一，多年来一直受到国内外众多学者的关注.ZHAO等［1］研究激光辐照下HgCdTe探测器的热效应；王玉恒等［2］和刘峰等［3］研究激光辐照下结构件的热响应；刘全喜等［4］采用有限元法研究光伏型光电探测器的激光加热效应；唐芳等［5］用数值模拟研究激光辐照下的皮肤热效应；孙承伟［6］归纳和总结国内外学者的研究成果.上述研究工作主要基于几何结构简单的计算模型，未涉及入射激光在模型中多次散射的问题.
　　射线追踪有限元模拟方法综合射线追踪技术和有限元法的优点，可以有效模拟入射能量在复杂结构中的传输与分布.近年来，关于射线追踪有限元模拟方法的研究逐渐增多，GRUHLKE［7］利用该方法成功实现NPSAT1卫星的热分析；FIALA等［8］采用该方法在MATLAB中开发光学照明系统的数值模拟程序；ZIEMELIS［9］介绍该方法在电磁散射问题中的应用；ZHANG等［10］采用该方法研究面泵浦激光二极管内三维应力场和温度场的分布.虽然射线追踪有限元模拟方法已经在多个领域取得成功应用，但是在激光辐照热效应研究方面尚未见到相关的文献报道.
　　本文在三维多物理场流固耦合问题并行计算程序CFS中实现激光辐照热效应的射线追踪有限元模拟方法，并以2A12硬铝材料腔体结构件为例，采用该方法计算其在激光辐照下的温度场分布，分析该方法与传统有限元法计算结果的差异，并探讨存在多次散射光时射线追踪有限元模拟方法的优越性. 　　1 基本理论
　　射线追踪有限元模拟方法是一种将射线追踪技术与有限元法联合起来对物理问题进行求解的数值模拟方法.在激光辐照热效应问题求解中，在激光光源位置随机产生大量沿激光传输方向的射线；追踪这些射线，获取射线在计算模型有限元网格中的沉积情况，从而得到每个有限元网格所吸收的激光能量，进而实现物理问题的求解.
　　1.1 射线追踪技术简介
　　射线追踪技术实质上是一种近似模拟方法，适用于空间特征尺寸远大于电磁辐射波长条件下的电磁辐射的散射、传输和透射等问题.在激光辐照热效 应问题求解中，如果入射激光仅被计算模型表面吸收，那么可按图1所示的流程进行射线追踪计算.
　　由图1可知，在射线追踪计算流程中存在激光光源的随机抽样、射线与边界单元表面第一个交点的计算和交点处新射线的产生等3个关键步骤.激光光源的随机抽样要求随机抽样的位置应能准确反映激光的强度分布信息，如对于均匀激光，应产生服从均匀分布的高质量随机数.［11］在计算射线与边界单元表面第一个交点时，主要困难是计算区域可能存在多个边界单元，若直接求解，则会面临大量的空间直线与空间直线(或空间平面)相交问题的计算，导致计算速度缓慢，因此，在计算射线与边界单元表面第一个交点时引入加速算法.以三维问题为例，本文通过八叉树技术将边界单元进行分组，使得需要求解的空间直线与空间平面相交问题的数目由N降至log8N，大大提高计算效率.至于交点处新射线的产生，则要保证其方向服从一定的物理规律，本文假定计算模型具有灰体性质，交点处新产生射线的方向遵守Lambert定律.
　　1.2 热传导问题的有限元求解
　　对于一般的三维热传导问题，瞬态温度场的场变量 T(x,y,z,t)在直角坐标系中满足式(1)所示的偏微分控制方程［1213］ ρcTt－xkxTx－ykyTy－zkzTz－ρQ=0 在计算域Ω内
　　(1)
　　式(1)的边界条件有3类，分别为
　　T= T -在Г1边界上
　　kxTxnx+kyTyny+kzTznz=q在Г2边界上
　　kxTxnx+kyTyny+kzTznz=h(Ta－T) 在Г3边界上
　　(2) 式中：ρ为材料密度，kg/m3；c为材料比热容，J/(kg·K)；t为时间，s；kx，ky和kz分别为材料沿x，y和z方向的热导率，W/(m·K)；Q=Q(x,y,z,t)为体热源密度，W/kg；nx，ny和nz为边界法线的方向余弦；T= T - (Γ,t)为在Г1边界上的给定温度，K；q=q(Г,t)为在Г2边界上的给定热流密度，W/m2；h为表面传热系数，W/(m2·K)；对于Г3边界，Ta=Ta(Г,t)在自然对流条件下为外界环境温度，在强迫对流条件下为边界层的绝热壁温度.
　　此外，求解瞬态温度场还应给定初始条件 T(x,y,z,0)=T0(x,y,z)
　　(3) 由变分原理和Galerkin近似，经推导可以得到式(1)的积分形式为： ∫Ω δTρcTt+ δTx kxTx+δTykyTy+δTzkzTz－δTρQdΩ－
　　∫Γ2δTqdΓ－∫Γ3δTh(Ta－T)dΓ=0
　　(4)
　　将求解域Ω离散为有限个单元体，单元内的温度T可近似用单元节点温度Ti插值得到，即 T=nei=1Ni(x,y,z)Ti(t)
　　(5) 式中：ne为单元节点数目；Ni(x,y,z,t)为插值函数.将式(5)代入式(4)，并考虑δTi的任意性，即可得到用于确定n个节点温度Ti的有限元求解方程 C T · + KT= P
　　(6) 式中： C为热容矩阵； K为热导矩阵； P为温度载荷列阵； T为节点温度列阵.
　　2 数值算例
　　假设一个腔体结构由六面体一端挖去一个小六面体形成，其示意见图2，外围六面体几何尺寸为0.20 m×0.20 m×0.30 m，挖去部分几何尺寸为0.10 m×0.10 m×0.25 m.假定腔体由2A12硬铝材料构成，其激光吸收因子取为0.25，不同温度下该材料的比热容、热导率和发射系数等选取文献［14］和［15］中的给定值. 在腔体模型中建立原点位于六面体中心、坐标轴方向如图2所示的坐标系.假设功率密度为1.5 kW/cm2的均匀激光沿(45.0°，90.0°，-45.0°)方向辐照腔体，并假设激光光斑形状为矩形，分别采用传统有限元法和射线追踪有限元模拟方法计算激光辐照下腔体温度场分布情况.
　　采用传统有限元法时，首先用几何方法确定腔体的激光辐照区域，然后将入射激光能量等效为热流边界施加到该区域；采用射线追踪有限元模拟方法时，通过射线追踪方法从光源出发追踪大量射线，得到入射激光能量在腔体内的准确分布.20.0 s时腔体三维温度场分布见图3.
　　由图3可知，与传统有限元法相比，射线追踪有限元模拟方法计算所得腔体最高温度大约增加58.6 ℃，且具有更大范围的高温区域.这主要是因为2A12硬铝材料对激光吸收系数较小(计算时取0.25)，当入射激光首次到达腔体表面时，只有小部分能量被直接吸收，大部分能量以散射光的形式在腔体内重新分布或逸出腔体，使腔体接收激光能量的实际面积和吸收的实际激光能量均大于传统有限元法的描述.
　　为进一步研究多次散射光的影响，选取光斑中心(0.05，0.0，0.1)及其yOz平面对称点(-0.05，0.0，0.1)绘制节点温度时间曲线，见图4.
　　由图4可知，光斑中心温度的射线追踪有限元模拟方法计算结果明显大于传统有限元法计算结果，说明腔体内部存在较强的多次散射光，光斑中心所吸收的能量是入射激光能量与多次散射光能量之和.采用射线追踪有限元模拟方法计算时，光斑中心yOz平面对称点温度随时间的增加而上升，20.0 s内累积温升约为173.7 ℃；采用传统有限元法计算时，该点温度在20.0 s计算时间内保持初始温度(22.7 ℃)不变，说明20.0 s时该点尚未受到热传导影响，节点温升由多次散射光引起，再次证明腔体内存在较强的多次散射光. 　　图3和4表明，当激光辐照如图2所示的2A12硬铝材料腔体结构件时，由于腔体内存在较强的多次散射光，使得传统有限元法无法准确描述入射激光能量在腔体内的分布与沉积，导致计算偏差较大.因此，当由吸收系数较小的材料构成的腔体结构件被激光辐照时，腔体内会产生较强的多次散射光，腔体实际吸收激光能量的表面积和吸收的激光能量均难以直接描述，传统有限元法已不再适用，而射线追踪有限元模拟方法对复杂入射边界具有很好的适应能力，可方便地求解该类问题.
　　3 结束语
　　与传统有限元法相比，射线追踪有限元模拟方法对计算模型几何复杂性具有更强的适应能力，可用于求解入射激光在计算模型中多次散射的问题，适用范围更广泛.
　　对2A12硬铝材料腔体结构激光辐照热效应的计算表明：当材料激光吸收系数较小、腔内多次散射光较强时，射线追踪有限元模拟方法比传统有限元法在准确描述计算区域内入射激光能量分布方面存在明显优势，能够有效降低该类问题求解难度. 参考文献：
　　［1］ ZHAO Jianhua, LI Xiangyang, LIU Hua, et al. Thermal effect in laserHgCdTe interaction［J］. Int Soci Opt Eng, 1998(3436): 961969.
　　［2］ 王玉恒, 刘峰, 陈林柱, 等. 强激光辐照下充压壳体损伤特征温度分析［J］. 固体火箭技术, 2008, 31(5): 512516.
　　WANG Yuheng, LIU Feng, CHEN Linzhu, et al. Analysis of characteristic temperature on damage of pressurized metal shell irradiated by intense laser beam［J］. J Solid Rocket Tech, 2008, 31(5): 512516.
　　［3］ 刘峰, 吴振森, 王玉恒, 等. 重复频率激光辐照柱壳的加热效率［J］. 中国激光, 2006, 33(4): 461466.
　　LIU Feng, WU Zhensen, WANG Yuheng, et al. Heating efficiency of repetitive frequency intensity laser irradiated cylinder［J］. Chin J Lasers, 2006, 33(4): 461466.
　　［4］ 刘全喜, 齐文宗, 郝秋龙, 等. 激光辐照光伏型光电探测器热效应的有限元分析［J］. 应用光学, 2007, 28(3): 275279.
　　LIU Quanxi, QI Wenzon, HAO Qiulong, et al. Finite element analysis of thermal effect of photovoltaic detector irradiated by laser［J］. J Appl Opt, 2007, 28(3): 275279.
　　［5］ 唐芳, 牛燕雄, 张雏, 等. 脉冲激光辐照皮肤组织的热效应数值模拟研究［J］. 激光与红外, 2007, 37(12): 12591261.
　　TANG Fang, NIU Yanxiong, ZHANG Chu, et al. Numerical analysis of the temperature field in skintissue by high power Nd: YAG pulses laser［J］. Laser & Infrared, 2007, 37(12): 12591261.
　　［6］ 孙承伟. 激光辐照效应［M］. 北京: 国防工业出版社, 2002: 2878.
　　［7］ GRUHLKE M. Computer aided thermal analysis of a technology demonstration satellite(NPSAT1)NPSSP03001［R］. Montery: Naval Postgraduate School, 2003.
　　［8］ FIALA P, KADLECOVA E. Progressive methods in the numerical modeling by the finite elements［C］//AEE’05 Proc 4th WSEAS Int Conf Applications Electr Eng. Stevens Point, Wisconsin, USA: World Scientific and Engineering Academy and Society (WSEAS), 2005: 182185.
　　［9］ ZIEMELIS J. Problems of secondary sources for hybrid ray tracing and FEM simulation［J］. Elektronika IR Elektrotechnika, 2003, 47(5): 6568.
　　［10］ ZHANG Shenjin, ZHOU Shouhuan, TANG Xiaojun, et al. Investigation of laser diode facepumped high average power heat capacity laser［J］. Chin Optics Lett, 2006, 4(11): 658660.
　　［11］ 杨自强. 你也需要蒙特卡罗方法——一个得心应手的工具［J］. 数理统计与管理, 2007, 26(1)： 178188.
　　YANG Ziqiang. You need Monte Carlo method: a handy tool［J］. Application Stat & Manage, 2007, 26(1): 178 188.
　　［12］ 王勖成. 有限单元法［M］. 北京：清华大学出版社, 2003: 441467.
　　［13］ LEWIS R W, MORGAN K, THOMAS H R, et al. The finite element method in heat transfer analysis［M］. Hoboken: John Wiley & Sons, 1996: 132.
　　［14］ 张文斌. 2A12铝合金搅拌摩擦焊温度场及塑性材料流动的数值分析［D］. 济南: 山东大学, 2010.
　　［15］ 左永平. 2A12, TC4微观组织及表面状态对红外发射率的影响研究［D］. 南京: 南京航空航天大学, 2008.
　　(编辑 陈锋杰)