（动力工程及工程热物理专业论文）基于有限体积法的超短脉冲激光瞬态辐射研究.pdf

国内图书分类号： tk124 国际图书分类号： 621.1.016.4 工学硕士学位论文 基于有限体积法的超短脉冲激光 瞬态辐射研究 硕 士 研 究 生： 王大林 导师： 阮立明教授 申请学位级别： 工学硕士 学 科 、 专 业： 动力工程及工程热物理 所 在 单 位： 能源科学与工程学院 答 辩 日 期： 2009 年 7 月 授予学位单位： 哈尔滨工业大学 classified index： tk124 u.d.c.： 621.1.016.4 dissertation for the master degree in engineering research of the ultra-short-pulse laser transient radiation based on finite volume method candidate： wang dalin supervisor： prof. ruan liming academic degree applied for： master of engineering specialty： engineering thermophysics unit： school of energy science finite volume method; computational efficiency; ultra-short-pulse light; participate medium - iii - 哈尔滨工业大学工学硕士学位论文 目录目录 摘 要.i abstract. ii 第 1 章 绪论.1 1.1 课题背景及意义.1 1.2 国内外研究现状.3 1.2.1 瞬态辐射.3 1.2.2 有限体积法.4 1.3 主要研究内容.5 第 2 章 求解瞬态辐射传输的数值方法.7 2.1 离散坐标法.7 2.2 有限元法.8 2.2.1 最小二乘有限元法.9 2.2.2 立体角的离散策略.10 2.3 有限体积法.11 2.3.1 用有限体积法对瞬态辐射方程的离散.11 2.3.2 入射短脉冲激光的处理策略.15 2.4 本章小结.17 第 3 章 一维介质瞬态传输的有限体积法.18 3.1 一维平板瞬态辐射传输的离散方程.18 3.2 计算的流程.21 3.3 一维均匀介质的计算.23 3.3.1 程序的验证及算法对比.23 3.3.2 均匀介质中物性变化对透射率和反射率的影响.25 3.3.3 脉冲宽度对透射率和反射率的影响.27 3.3.4 不同入射角度对结果的影响.29 3.4 一维非均匀介质的计算.32 3.4.1 一维两层介质.32 3.4.2 一维三层介质.33 3.4.3 一维连续分布型不均匀介质.34 3.5 本章小结.36 - iv - 哈尔滨工业大学工学硕士学位论文 第 4 章 超短脉冲激光在多维介质内的传输.37 4.1 二维参与性介质内的离散方程.37 4.2 二维均匀介质的计算.37 4.3 二维非均匀介质计算.40 4.3.1 尺寸不同的散射核.40 4.3.2 含有两个散射核的情况.43 4.4 三维介质的计算.46 4.4.1 三维介质的离散方程.46 4.4.2 三维均匀介质的计算.46 4.5 本章小结.48 结 论.50 参考文献.51 攻读学位期间发表的学术论文.56 哈尔滨工业大学硕士学位论文原创性声明.57 哈尔滨工业大学硕士学位论文使用授权书.58 哈尔滨工业大学硕士学位涉密论文管理.59 致 谢.60 - v - 哈尔滨工业大学工学硕士学位论文 第第1章 绪论章 绪论 1.1 课题背景及意义课题背景及意义 热辐射（thermal radiation） ，作为能量和信息传递的一种方式，近年来在 军用和民用等高新技术领域都得到了广泛的发展和应用1, 2。其应用的深度和 广度不断加大，而且与其他学科的交叉也日益增多。辐射是一种电磁波，它的 传播速度是光速，在通常情况下，都是假设辐射传输是瞬间达到稳态的。众所 周知，辐射是一种电磁波，它的传播速度是光速。在通常情况下，都是假设辐 射传输是瞬间达到稳态的。然而，在工程应用中，如果得到热流密度在 s内随时间变化，那么热辐射在介质中的传输就变成了一个瞬态问 题。当光学参与性介质受到短脉冲激光辐射时这种情况就会发生，短脉冲的脉 冲宽度范围在 9 10 15 10 9 10s。短脉冲辐射源可以产生间隔非常短的信号。信号持 续的时间和辐射源的时间是相同数量级的。短暂的时间信号取决于介质的特 性，那么信号就会反映出介质的信息。随着超快激光（ultra-fast laser）的快 速发展及微尺度系统 15 10 3的研究，半透明介质内的瞬态辐射传输受到了广泛的关 注。短脉冲激光对金属的加工4, 5、医学成像技术6、无损伤激光治疗技术7、 粒子的探测研究8、红外遥感技术9等许多方面的进展表明：瞬态辐射的传输 过程需要理论意义上的深刻理解和精确的数值模拟。例如，短脉冲激光医学成 像技术发展的关键之一就是要有一个有效而准确的前向计算。如果没有一个好 的前向计算，那么反向的分析就不可能提供高精度的结果。除了上述新兴技术 的应用背景之外，在传统的辐射传输领域中也需要准确高效的瞬态辐射数值模 型。例如，在天体物理辐射问题10的研究中，由于空间尺度的特性决定了必须 考虑辐射传输随时间的变化。 瞬态辐射传输最具活力的应用就是关于生物光学成像的研究。由于涉及光 学工程、生物医学、计算数学、辐射传输等多个学科领域的专业知识，光学成 像技术的发展一直受到许多学者的关注。其诱人的应用背景，也成为瞬态辐射 研究最直接的动力11。在这项技术中，近红外光（nir）被当作探测光源照射 到生物体的表面，通过生物体的散射、吸收后，光学信号被通过一些设备观测 记录。然后，这些信号被当作输入信息，通过计算机反演得到生物体内部的光 - 1 - 哈尔滨工业大学工学硕士学位论文 学属性分布。通过这样的技术，人们可以监测隐藏在生物体内的一些异常或者 变化，如组织含氧量的变化12、脑出血的检测13、肿瘤的初期发展14等。与 传统的x射线ct、核磁共振（mri） 、正电子成像（pet）等成像技术相比，这 种通过非电离的近红外光成像的技术对人体几乎无害。由于它不需要对放射源 进行维护，因而更加廉价。此外，如果应用近红外的超快激光作其发射源，这 项技术的分辨率也比传统的成像技术提高数倍15。因此，光学成像技术被认为 是对传统的成像技术更安全更廉价的一个补充和替换选择。近期，短脉冲激光 入射在半透明介质中的瞬态辐射传输吸引了许多学者的关注。 另外，瞬态辐射信号作为检测信号能够提供很多在稳态辐射下得不到的信 息，所以很多学者对瞬态辐射传递问题进行了研究，发展了不少的数值计算方 法。 有限体积法自诞生之日起就得到了迅速而持续的发展，这一方法已经在力 学、传热学、电磁学、声学等领域中得到了广泛的应用。作为一个成熟的数值 方法，它有着如下的优点： （1）有限体积法的基本思路易于理解，并能得出直接的物理解释 （2）在有限体积法中，插值函数只用于计算控制体积的积分，得出离散方 程之后，便可忘掉插值函数；如果需要的话，可以对微分方程中不同项采 取不同的插值函数； （3）限体积法得出的离散方程，要求因变量的积分守恒对任意一组控制体 积都得到满足，对整个计算区域，自然也得到满足； （4）易于发展与其他模拟兼容的耦合分析。 尽管如此，这一方法在辐射传输中的应用仍然存在着许多令人关注的问 题，一些关键技术难点仍没有完全解决，如对不规则区域的适应问题，数值误 差的控制，各向异性散射问题的处理等等。因此，求解辐射传输的有限体积法 仍有待于从理论和应用的角度进行深入的研究。 - 2 - 哈尔滨工业大学工学硕士学位论文 1.2 国内外研究现状国内外研究现状 1.2.1 瞬态辐射瞬态辐射 早在1988年，patterson和chance16提出了建立在瞬态辐射传输的p1模型基 础上的漫射近似求解方法。他们将这一方法用于研究短脉冲激光与光学厚介质 的交互作用，并与蒙特卡洛法的计算结果进行了对比。1999年kumar和mitra17 提出很多关于要求解瞬态辐射传递方程（trte）的问题。这包括：混浊介质 中的短脉冲激光18 ，19，材料加工20，21及活体组织22中的脉冲激光等。由于描 述瞬态辐射波传递的微分方程是双曲型的偏微分方程，辐射强度同时是时间， 空间，方向，波长四种参数的函数。因此，与稳态辐射模拟相比，瞬态辐射传 递的模拟要复杂的多，它也成为近十年来光学和计算辐射学的研究热点。近十 年中，已经有许多理论研究和数值方法被提出。关于求解瞬态辐射传递方程的 方法有：积分方程（ie）法23-27，蒙特卡洛（mc）方法28-30，球谐函数法31, 32，离散坐标法（dom）31, 33，和有限体积法（fvm）34-38。 1997 年，mitra和kumat39究一维平板瞬态辐射问题时，比较了四种方法 即：球形函数法，热流近似法，离散坐标法，和直接数值积分法的精确性。比 较了球形函数法（p-n） ，热流近似法和离散坐标法（dom）的结果。离散坐 标法（dom）被广泛的用于参与性介质中稳态辐射传递。最近几年，离散坐 标法也用于求解瞬态辐射问题，在一维问题40，二维问题41，三维问题31中 的应用都有报道。mitra等42用p1 近似法来模拟矩形封闭体内的瞬态辐射传 递。guo和kumar43提出了一种解决平板介质中的瞬态辐射传递的数值方法 辐射元法（radiative element） 。tan和hus44用离散矩形体积法和yix法来模 拟立方封闭体内的瞬态辐射传递，给出了均匀性和非均匀性介质中辐射热流 （radiative flux）和附加辐射（incident radiation）的演变。2000 年，tan等23提 出了瞬态辐射的积分方法。他们的方法是建立在积分型辐射传递方程的基础上 的一种拉格朗日观点的方法，具有精度高的特点，但是这种方法在处理多维问 题的各向异性散射时比较复杂。sakami30和他的同事用monte-carlo方法对瞬 态辐射传递进行了研究。对于多维复杂形体内的瞬态辐射问题，离散坐标法和 有限体积法是比较流行的求解方法。但这两种方法在处理复杂边界条件时，仍 不够灵活。2006 年，run和an等45用最小二乘变分原理的有限元法对多维介质 - 3 - 哈尔滨工业大学工学硕士学位论文 瞬态辐射传输进行模拟，并提出最小二乘有限元法求解瞬态辐射传输的数学模 型。2006 年an提出了用最小二乘有限元法来研究瞬态辐射传递问题46。wu47 等用改进的差分格式和混合的p1-p3模型分析了一维介质中的瞬态辐射问题。另 外 ， 还 有 一 些 其 他 的 方 法 如 ： 改 良 特 征 线 法 （ modified method of characteristics）48、小波分析法(wavelet)49。 1.2.2 有限体积法有限体积法 有限体积法由raithby和chui50首先提出，并将其应用于直角坐标系及柱 坐标系中的辐射问题的求解。chai51将有限体积法应用到多维辐射传输的求解 中。他的基本思想是对辐射传递方程的角度及空间离散均采用有限体积法离 散。有限体积法对角度离散比离散坐标法更自由，而且容易处理空间上不规则 几何形状。由于其本身的守恒特性，有限体积法可以减少射线效应的影响。 2002年chai37利用了有限体积法来计算在二维不规则封闭体内瞬态辐射传 递，采用全隐式格式离散瞬态项。并且利用界面热流及附加的热辐射验证了运 用于吸收散射各项异性散射介质模型中的公式的正确性。目前还没有关于用 do法或fvm法来求不规则几何体的瞬态辐射问题，本文作者是第一个把有限 体积法用于不规则几何体的瞬态辐射问题。chai36提出了用有限容积法来研究 一维板介质瞬态辐射传递的问题。研究了辐射漫射和平行光的影响，吸收和各 向同性散射，连续脉冲辐射都，还提出了脉冲宽度在单脉冲辐射和多脉冲辐射 中的作用。运用了全隐的时间差分格式来离散trte中的瞬态项。用step和 clam格式用来作为空间离散格式。他的研究仅限于介质中的入射光瞬态变 量。用设备探测出的边界处的射和反射信号在他的文章中并没有涉及。也就是 他的研究专注于方形阶梯脉冲，并且没有考虑散射反照率的影响。 chai38构建了瞬态辐射传递的三维模型，利用了他和其同事文献51, 52中的 有限体积法的瞬态形式。对step和clam这两种格式用来检测通过吸收、发射 各向同性散射介质后辐射能量的演变进行了比较。结果表明，clam格式得出 的结果比step格式要精确。并分析了射线效应和假散射。 2006年，mishra53等推导出了单脉冲平行阶梯脉冲公式，比较了离散传 递法（dtm） ，离散坐标法（dom） ，有限体积法（fvm）三种方法。三种方 - 4 - 哈尔滨工业大学工学硕士学位论文 法得到的结果都很接近。fvm方法是很常用的方法并且对任何几何型都有较好 的适应性。不像其他两种方法，它又很小的射线效应。 有限体积法法中的离散和近似可避免的要引起两种误差：假散射（false scattering）和射线效应（ray effect） 。自从rathiby54和cai55提出了fvm和 dom中的假散射和射线效应的问题以来，许多学者对这两种误差进行了研 究。rathiby54认为，这两种误差分别是由于空间离散和立体角离散引起的， 并且它们之间是相互影响的。由于方向离散引起的误差试图将辐射热流集中在 离散方向的中心线上，而由空间离散引起的误差通常导致了辐射强度的假散射 （假漫射）效应，其作用是将传递方向上的强度分布向邻近网格漫射，因此这 两种误差的相互作用有趋于抵消的趋势。此外，rathiby54还指出：将辐射模 型的误差全部的归结于射线效应和假漫射的相互作用是不完整的。实际上，误 差的产生于求解方法、空间离散格式、立体角离散方法的选择有着很大的关 系。而且，未加迎风修正的高阶空间离散格式也会对辐射传输产生较大的误 差。 求解辐射传递的方法很多，但没有哪一种方法在所有方面优于其它的方 法。因此人们从来就没有放弃去探索一个计算精度高、求解速度快，能方便的 与流体流动、燃烧反应等物理现象的模拟方便融合的，高效的辐射换热计算方 法的努力。而且，近期的一些研究在已有的这些方法的基础上，更加关注数值 误差的减小，求解效率的提高，以及处理不规则边界及与其他模拟相互融合等 问题。 1.3 主要研究内容主要研究内容 综合考虑求解精度与效率，本文将开展求解瞬态辐射传输的有限体积方法 研究。利用有体积方法对角度和空间进行离散，将入射激光脉冲在介质内引起 的辐射传输分为平行辐射和漫射辐射两部分，采用适当的时间差分格式，建立 求解瞬态辐射传递的有限体积法模型，其具体内容将包括以下几个方面： （1） 综述有限体积求解法及瞬态辐射传输研究的应用背景、应用领域及应用 前景； （2） 优化并完善现有的有限体积法程序，建立求解瞬态辐射传递方程的有限 - 5 - 哈尔滨工业大学工学硕士学位论文 体积法数值模型，并编制相应的软件； （3） 针对若干一维、多维算例对比分析瞬态辐射的有限体积法求解程序与原 有的有限元法求解程序间的计算效率、计算精度及误差控制等方面的差 异； （4） 研究脉冲激光辐射通过参与性介质的传递以及分析在边界处得到的瞬态 信号，分析其时间及空间的分布规律。 - 6 - 哈尔滨工业大学工学硕士学位论文 第第2章 求解瞬态辐射传输的数值方法章 求解瞬态辐射传输的数值方法 近期，随着超快激光的飞速发展，半透明介质内的瞬态辐射传输受到了广 泛的关注。求解瞬态辐射问题的主要方法和稳态问题类似，主要有离散坐标法 （dom） ，有限元法（fem） ，有限体积法（fvm）等。本章从瞬态辐射传输 方程出发，对上述方法予以简单介绍，重点对有限体积发的原理和瞬态辐射传 递方程的有限体积法离散进行阐述。 直角坐标系下，描述参与性介质内瞬态辐射传输的方程可以表示为如下的 形式： s asab 4 ( ,)( ,) () ( ,)()( , ) ( , ) 4 tt ttt c ts d = +=+ i ri r i ri r,i r , , (2-1) 其中 是介质中光的传播速度，其它符号与稳态辐射传输方程相同。由上 式可见，与稳态辐射传输相比，瞬态辐射传输方程增加了对时间的导数项，辐 射强度也同时成为位置、方向、波长及时间的函数。从数学的角度看，描述瞬 态辐射传输的方程也由稳态的拟线性微分-积分方程，变为具有波动性质的双 曲型微分-积分方程。 c 2.1 离散坐标法离散坐标法 离散坐标法(dom)最早在 1965 年由 love 等将其引入到一维辐射换热的求 解中。与稳态辐射传输方程的求解类似的，我们采用离散坐标法来处理散射积 分项，将方程(2-1)离散为 m 个独立方向上的辐射传输方程。对于特定的方向 ，它沿坐标轴mx，的方向余弦可以表示为yz m ， m ， m ，是入射 方向和散射方向上的散射向函数。这样，在直角坐标系下，瞬态辐射传 输方程可以表示为如下形式： m m mm asab 1 () 4 mmmm m mmm m s mmmm m iiii iiiw c txyz = += + (2-2) 如果定义下面的简化系数，并考虑方形短脉冲入射条件，则方程(2-2)的加权余 量形式可以表示为方程(2-3)的形式： - 7 - 哈尔滨工业大学工学硕士学位论文 s as 4 mm m bw =+ (2-3) s a b 1, 4 m mm m m mmm siiw = s =+ + (2-4) ( * inp0 e()(), 4 ls sih llh lll ) = (2-5) 这里是衰减系数，l是当前位置与入射边界的距离， 0 是脉冲入射方 向， * p l和分别定义为：， * p l * lc= t * pp lct=， p t是脉冲时间宽度。 in i 为入射脉冲 强度，h 表示 heaviside 函数。 离散坐标法的一个严重缺点就是假散射（false scattering）问题，这是空 间坐标离散误差而得到的，类似于在 cfd 计算结果中产生的数值漫射 （numerical diffusion） 。如果一束平行光穿过一个封闭区域，用离散坐标法来 求解，那么光束会随着离起始点的距离的增加而逐渐的变宽。即使不存在真实 的散射，这种辐射强度的非物理拖尾效应，就是所谓的假散射。可以通过更优 异的控制体积的网格来减小假散射的影响。 dom 方法另一个更严重的缺点就是射线效应（ray effect） ，这是由于角 度离散而产生的。考虑一个非常小的封闭区域（体积或表面积） ，具有很大的 发射率。这个区域的辐射强度会被从此区域带到各个离散立体角的方向上。远 离该发射区域的这些射线可能变得很远，而远离某些控制体积，或者表面区域 没有接收到来自高发射率区域的能流，从而产生了非物理性的结果。显然，射 线效应可以通过增大控制体和划分的表面区域来减小。因此，当使用好的空间 网格来减小假散射作用，这又会伴随着增大该方法的阶数。 一般采用阶梯格式（step scheme）不会出现这种非物理性的结果，但是 阶梯格式相对于其它的格式（如菱形格式）误差会大些。 2.2 有限元法有限元法 有限单元法作为数值分析方法的一个重要特点是利用在每一个单元内假设 的近似函数来分片地表示全部求解域上待求的位置场函数。单元内的近似函数 - 8 - 哈尔滨工业大学工学硕士学位论文 通常由未知场函数或其导数在单元各个节点上的数值和其插值函数来表达。因 此，未知场函数或其导数在各个节点的数值就成为新的未知量(自由度)，从而 使一个连续的无限自由度问题变成离散的有限自由度问题，一经求解出这些未 知量，就可以通过插值函数计算各个单元内场函数的近似值，从而得到整个求 解域上的近似解。显然随着单元数目的增加，即单元尺寸的缩小，或者随着单 元自由度的增加及插值函数精度的提高，近似解的精度将不断得到改进。在有 限单元法中，由于对单元作了积分计算，就充分估计了不同节点参数的贡献， 从而克服了有限差分法或有限体积法中不考虑单元本身特性的缺点。有限元法 由于其独有的特性成为目前用于工程计算的热门方法之一。 2.2.1 最小二乘有限元法最小二乘有限元法 如果取权函数 1 ()( n j i jj wl n aln aa = = ) ii a i ，并且使得微分方程 “弱”形式的余量在最小二乘意义下的积分达到最小值时的加权余量法称为 最小二乘有限元法46。其实质是使函数： 2 1 ( )() n ii i i aln a d = = (2-6) 取最小值，即要求0 i i a = ，(1,2, )in=l。 方程(2-2)的加权余量形式可以写为： () mmmm m mmml v iiii bis wdv c txyz + 0= (2-7) 这里v是求解域，是有限元权函数。通过应用最小二乘有限元离散技术， 方向上的辐射强度可以写为下面的矩阵形式： l w m 1 mm nn kiki + =+ r (2-8) 这里 1 m n i + 是时刻的辐射强度，它能够由第时刻的辐射强度值1n+n m n i求出， k是刚度矩阵，是右端项，它们的详细表达式如下： r - 9 - 哈尔滨工业大学工学硕士学位论文 *() *() iii immmi jjje e jmmmj tb xyz kd tb xyz + = + xdydz (2-9) (1)*() *() iii immmi jjje e jmmmj tb xyz kd tb xyz + = + xdydz (2-10) 1 *()*() jjj nnijmmmj e e rtsstbdxdydz xyz + =+ (2-11) 这里是时间格式差分因子，可以取 1，0.5，2 3，或者 0，分别代表全隐 式，半隐式，伽略金和显格式时间差分。时间步长定义为* tc t= 。用同样的 方法也可以写出其他方向上辐射强度的求解方程，然后应用与离散坐标法类似 的技术，就可以求得介质内不同时刻的辐射强度，进而求出边界上的透射和反 射信号。 2.2.2 立体角的离散策略立体角的离散策略 对于立体角离散，与稳态辐射传输相同，可以采用sn型和天顶角-圆周角 型两种离散方式。 an46采用了dom+isw型的立体角离散策略，即将一个独立附加的，具有 无限小积分权重的立体角方向加入到dom的求解中。这个附加的立体角方向 可以是任意的方向，因此可以模拟平行光的入射方向，将其标记为m+1 方向， 它的积分权重设为 10-10。如果入射光平行于求解域的x方向，则相应的方向余 旋 1m + ， 1m + 和 1m + 可以分别设为 1.0，0.0 和 0.0。将任意方向的入射脉冲用 该方法处理，但在辐射传输方程中添加入射脉冲所引起的源项，并与其他方 向散射引起的源项耦合在一起进行求解。同其它的立体角积分格式相比，这种 方法非常简单，也不会违反s s n格式的矩条件。 - 10 - 哈尔滨工业大学工学硕士学位论文 2.3 有限体积法有限体积法 该方法由raithby和chui51首先提出，并将其应用于直角坐标系及柱坐标系 中辐射问题的求解。chai52将有限体积法应用到多维辐射传输的求解中。它的 基本思想是对辐射传递方程的角度及空间离散均采用有限体积法离散。有限体 积法对角度的离散比离散坐标法更自由，而且容易处理空间上不规则几何形 状。由于其本身的守恒特性，有限体积法可以减少射线效应的影响，保证总的 能量平衡。 2.3.1 用有限体积法对瞬态辐射方程的离散用有限体积法对瞬态辐射方程的离散 2.3.1.1 短脉冲激光作用于参与性介质的瞬态辐射方程短脉冲激光作用于参与性介质的瞬态辐射方程 为书写简便起见瞬态辐射传递方程(2-1)写成： () a b 4 , 4 s ii ik iid c ts += + (2-12) 上式中的 i 是位置 ，角度，时间t的函数。对于入射激光脉冲在介质内的传 输问题，介质中的辐射强度可以分两部分 r c i ， d i ，前者是平行辐射强度，它 是与激光脉冲的入射方向同向的，并且它的最大值等于激光入射介质界面时的 强度，后者 d i 是漫射辐射强度。即介质中辐射强度为： cd iii=+ (2-13) 入射界面处，还没发生或者说刚要发生漫射，不考虑介质表面折射影 响，此时 0 x = cin ii=，。平行光辐射强度 d 0i = c i 满足： c c di i ds = (2-14) 把(2-13)，(2-14)式带入(2-12)中，并且有 c 0 i t = ，整理得： dd dcdd ii t issis c ts += += + (2-15) 其中， - 11 - 哈尔滨工业大学工学硕士学位论文 () s cc 4 , 4 si d = (2-16) () s da bd 4 , 4 sk iid =+ (2-17) x y z 0 in i 0 图 2-1 三维几何模型 相函数采用线性散射相函数如下： () , 1coscosa = + (2-18) 式中，、 分别为、方向与 x轴夹角，即为所对应的天顶角，如图 2- 1。系数可以为正值或负值，分别对应前向散射和后向散射，若，则 ，即为各向同性散射介质。 a0a = () , =1 t 2.3.1.2 有限体积法的离散有限体积法的离散 引用长度量纲时间，c为光速。方程 * tc=(2-15)变为： dd d * ii t is st += (2-19) 时间格式按照全隐格式向后差分上式变为： * *ddd dt * ( )()( ) ( )( ) itittdit * its t dst += (2-20) - 12 - 哈尔滨工业大学工学硕士学位论文 上式简写为53： * *d dtd ( ) ( )( )() dit * bitbs tcitt ds +=+ (2-21) 其中， * * 1 t b t = + ， * 1 c t = + ， * t为时间步长。 x y z w e n s t b p s p e t m r m s n e n v x y z t n （a）控制体各表面中心处的积分点 (b) 控制体各表面单位外法矢量 j n 图 2-2 直角坐标系下的控制体 在某一方向 m ，在微元立体角 m 内对(2-21)式积分有： * *d dtd ( ) ( )( )() mmm m mmm it * bdit dbstcitt s d + =+ (2-22) 在三维直角坐标系中上式可写为53： * *ddd d * td ( )( )( ) ( ) ( )() mmm mmmmm xyz mm ititit bdddit xyz bstcitt + m =+ (2-23) 其中 m x d、 m y d 、 m z d为权重 m j d () m m d = n s ? (2-24) 其中，n为控制体界面上的外法线方向，为辐射传输方向。三维空间如 m s图 2- 2中有： - 13 - 哈尔滨工业大学工学硕士学位论文 cossincossinsin m =+sijk (2-25) ,jx jy jz j nnn=+nijk (2-26) j n表示第 j 面的外法线方向，方向权值为： , 22 , 0.50.25 sin 2sin 2 22 sinsincoscos 222 0.5sinsin 22 mmm i mmmm y jz j mm x j d nn n = + + + 2 (2-27) 22 22 sin mm mm mm m mm m dd d + = = 2sinsin 2 m m m = (2-28) 微元体界面上的辐射强度值与微元体中心的辐射强度存在某种关联，即构成某 种空间的差分格式： () d,px d,exd,w 1 mmm if ifi=+ () d,py d,nyd,s 1 mmm if ifi=+ () d,pz d,tzd,b 1 mmm if ifi=+ (2-29) x f ， y f， z f 为差分因子，为控制体中心节点， ， ， pews ntb分别 为控制体的各个界面。当p1f=时即为阶梯格式，0.5f=时为菱形格式。由 (2-29)式带入式(2-23)消去下游界面的辐射强度项，得控制点，方向为m的辐 射强度 p 53： ( ) ( )( )( ) * d,xuyzyd,yuxzzd,zuxyp * , xyzxyz mmmmmmm x m d p m mmm ditf fditf fditf fs it x dddf f f b + = + + (2-30) 式中： () * ptpxy zd,pxy z m mmmm c v svsf f fittf f f b = + - 14 - 哈尔滨工业大学工学硕士学位论文 () () ss t,pcdc0a bd,p , 44 mmmmm sssik ii m =+= + 其中，是激光脉冲入射方向。 0 (2-30)式中表示点 d,xu impx轴方向上游控制体 界面的方向的辐射强度，表示点轴方向上游控制体界面的方向的 辐射强度，表示点轴方向上游控制体界面的方向的辐射强度，它们 与有下述关系： m m d,yu ipym d,zu impzm m j d d,xex d,xu d,xwx i, i i, mm m mm d d = 0 0 = = = (2-34) 把(2-32)式代入(2-16)中得： ()()()( *s cinp ,t, 4 x srti eh tlh tl ) 0 = (2-35) 在实际情况中，入射脉冲经常为时间和空间上高斯型的激光，如对于具有 半峰值时域宽度(fwhm)为 p t的高斯型激光，其辐射强度的解析关系式可以写 为： m+1 cin () l p i ef t =i (2-36) 这里( ) p f t是关于的形状函数。 p t - 16 - 哈尔滨工业大学工学硕士学位论文 2.3.2.2 脉冲链的平行光辐射强度脉冲链的平行光辐射强度 对于方形脉冲链，如图 2-4所示，独立方向的平行光辐射强度 m+1 c i有下列 计算公式： ()()()()( m+1* cinppp 1 *2tt1 *2t l ii eh tlnh tln = ) 0 (2-37) 2 p t in i c i t p t 图 2-4 方形脉冲链的时域特征 式中表示脉冲链的个数，在n图 2-4中3n =。 2.4 本章小结本章小结 本章主要介绍了求解短脉冲激光瞬态辐射问题的经典方法。阐述了这些