（工程热物理专业论文）求解矢量辐射传输方程的有限体积法研究.pdf

国内图书分类号：tk124 国际图书分类号：621.1.016 工学硕士学位论文工学硕士学位论文 求解矢量辐射传输方程的有限体积法研究 硕 士 研 究 生 ： 李 俊 导 师： 阮立明教授 申 请 学 位： 工学硕士 学 科 、 专 业 ： 工程热物理 所 在 单 位 ： 能源科学与工程学院 答 辩 日 期 ： 2011 年 6 月 授予学位单位： 哈尔滨工业大学 classified index：tk124 u.d.c.：621.1.016 a dissertation for the degree of m. eng. the finite volume method for vector radiative transfer equation candidate： li jun supervisor： prof. ruan liming academic degree applied for： master of engineering speciality： engineering thermophysics affiliation: school of energy science , ) =-( , , ) ( , , ; , )( , , )( , , ) d ( , , ) ( , , ;,) ( , , ; , ;,) 4 eab s i x y z k x y z i x y zkx y z ix y z s k x y z i x y zx y zdd (2-18) 式中，( , ) 为传播方向，其中cos，为天顶角的余弦，为圆周角， 1,1 ，0,2 。( , , ; , )i x y z 为空间位置( , , )x y z处，方向为( , ) 的 辐射强度；( , , ) b ix y z为黑体辐射强度； a k、 s k、 e k分别为介质的吸收、散 射、衰减系数。( , , ; , ;,)x y z 为散射相函数。 基于 stokes 不连续的叠加准则( the incoherent addition principle )，对于 准单色光的偏振辐射，在转轴取向随机的粒子群介质中的三维稳态矢量辐射 传输方程可以描述为： 21 01 d ( , , ; , ) =-( , , ; , ) ( , , ; , )( , , ; , )( , , ) d ( , , ; , ;,) ( , , ;,) eab i x y z k x y zi x y zkx y zix y z s x y zi x y zdd (2-19) 关于该方程中各项的具体形式可以见参考文献67的附件 a。 关于式(2-19)的完整推导是非常繁琐复杂的过程，很多学者对此做出了 研究，这里就不再复述。1960 年 chandrasekhar68给出了 rayleigh 散射大气 中下的表达式；1982 年 ishimaru 和 cheung69推导了非发射性粒子散射的矢 量辐射传输方程，tsang70推广到了含粒子发射的情况。 与标量形式的辐射传输方程相比，矢量辐射传输方程的表达式上有如下 改变： 1 矢量辐射传输方程中将i替换成斯托克斯矢量( , )ti i q u t； 2 将仅与位置相关的标量衰减系数 e k替换为与位置-方向相关的消 光矩阵 e k； 哈尔滨工业大学工学硕士学位论文 - 14 - 3 将仅与位置相关的标量吸收系数 a k替换为与位置-方向相关的吸 收（发射）矢量 a k； 4 将标量的散射相函数替换为通过旋转所得的散射相矩阵，且 原来标量辐射传输方程中的散射系数 s k也包含在散射相矩阵中。 从物理意义上来看，标量辐射传输方程没有考虑光的偏振，认为各个方 向分量的强度是一样的；矢量辐射传输方程考虑了光的偏振，加入了偏振信 息，由斯托克斯矢量可以看出，由原来的一个标量方程转化成一个含有四个 方程的方程组。 2.5 散射相矩阵散射相矩阵 2.5.1 定义定义 对于一个粒子群，散射相矩阵的定义是通过对这些穆勒散射矩阵取粒子 数和转轴取向取平均， 2 1 ( , ;,)( , )( , ;,; , )pn rmr k (2-20) 式中( , )n r 表示粒径分布，它不仅是粒子半径 r 的函数，还是粒子转轴取向 的函数。符号 表示该参量对粒径和转轴取向取平均。例如，对于球形 粒子，参量 表示 0 ( )n r dr ，式(2-20)的操作是为了确保散射相函数的 单位与矢量传输方程保持一致。 2.5.2 相矩阵的化简相矩阵的化简 一般情况下，散射相矩阵是一个4 4的矩阵，但是对于一些特殊的情 况，这个相矩阵中的某些元素为 0，可以进行化简。 (1) 六元素 六元素矩阵简化成（详细叙述参见 van de hulst16和 hansen，travis72 的文章） ： 哈尔滨工业大学工学硕士学位论文 - 15 - 12 25 34 46 00 00 00 00 pp pp p pp pp (2-21) 该矩阵的使用范围：(1)随机导向轴对称粒子；(2)随机导向非对称粒子，但 一半数的粒子是另外一半的镜像；(3)rayleigh 散射和 mie 散射。 (2) 球形粒子 球形粒子群的散射一般成为 mie 散射。在这种情况下，元素 15 pp， 63 pp，式(2-21)化简成四个元素： 12 21 34 43 00 00 00 00 pp pp p pp pp (2-22) 矩阵中各个元素的含义为73： 1 p 描述入射辐射与出射辐射能量分布的相函数，通常是非零 的。 21 /pp 表示非偏振辐射能入射后散射出辐射能的偏振度。对于垂直 偏振和 rayleigh 散射 2 0p ；对于水平偏振 2 0p 。 31 /pp 表示非水平、非垂直的“斜”偏振，即与散射平面成一定角 度的线偏振辐射。 41 /pp 表示圆偏振度。即与斜偏振辐射能转换成圆偏振的程度。如 果电矢量顺时针旋转，表示 4 0p ，否则 4 0p 。 (3) rayleigh 散射 发生散射的球形粒子远小于波长的散射，称为 rayleigh 散射。散射相矩 阵为： 22 22 1 cos1 cos00 1 cos1 cos003 4002cos0 0002cos p (2-23) 哈尔滨工业大学工学硕士学位论文 - 16 - 其中为散射角。但是由于粒子的各向异性，实验证明，实际上辐射能通常 并不是与(2-23)所给的相矩阵一致。对于各向异性粒子的 rayleigh 散射相矩 阵为： 22 22 1 cos1 cos00 1 cos1 cos003 4002cos0 0002cos 1000 0000 (1) 0000 0000 p (2-24) 其中， 11 2 , 1/21 (2-25) 为去偏振因子，是与波长有关的。一般来说，对于气体分子，取 0.03。关 于各向异性的 rayleigh 散射的研究有很多，这里不做过多介绍。 (4) 旋转对称的导向散射 对于矢量辐射能入射到旋转对称的散射体上（例如圆柱） ，其长轴平行于 水平面，且导向是随机的，tang 和 aydin73给出了 mueller 矩阵表达式， 000 1 000 1 2 12 000 2 0002 p pr m pr r (2-23) 其中 22 12 pss (2-24) * 21 rs s (2-25) j s为散射振幅函数。 哈尔滨工业大学工学硕士学位论文 - 17 - 2.5.3 参考平面的旋转参考平面的旋转 散射相矩阵的定义是在其散射的平面内完成的，而矢量辐射传输方程则 是在其自己的笛卡尔坐标系下定义的。如图 2-3 所示，在原点o发生的散 射，散射平面为 12 opp，入射射线( ,) 的子午线面是 1 poy，出射射线( , ) 的子午线面为 2 poy。 图 2-3 偏振参考面的坐标示意图 如上图所示，子午线平面 1 poy是产生散射源项的入射波 s i的平面。散射相 矩阵是在散射平面 12 opp上定义的。出射波强度 i i是在出射平面 2 poy上定义 的。因此，对于矢量形式的辐射传输方程，偏振量必须旋转角度 1 i至到散射 平面 12 opp上，应用散射相矩阵进行计算，再把计算得到的偏振量旋转出此 平面，旋转过角度 2 i得到旋转过的散射相矩阵，其数学表达形式为： 21 ( , ;,)() ( , ;,) ()li pli (2-26) 对于斯托克斯量( ,) t ii q u v，基本的偏振旋转矩阵为（详情可以 参见 liou71，1980） ， 哈尔滨工业大学工学硕士学位论文 - 18 - 1000 0cos2sin20 ( ) 0sin2cos20 0001 l (2-27) 利用球面三角学（详情参见文献71的附件 a） ， 1 i和 2 i可以由下面的式 子得到： 1 21/22 1/2 cos cos (1 cos)(1) i (2-28) 2 21/22 1/2 cos cos (1 cos)(1) i (2-29) 其中是入射波与反射波所夹的散射角度， 2 1/22 1/2 cos(1)(1)cos() (2-30) 式(2-28)、(2-29)中，当2时，取号；当0时， 取号。此外，evans 和 stephens40(1991)给出等效的表达形式为， 2 1/2 1 21/2 (1)sin() sin (1 cos) i (2-31) 2 1/2 2 21/2 (1)sin() sin (1 cos) i (2-32) 2 1/22 1/2 1 21/2 (1)(1)cos() cos (1 cos) i (2-33) 2 1/22 1/2 2 21/2 (1)(1)cos() cos (1 cos) i (2-34) 这样，以随机导向轴对称粒子为例，利用式(2-26)(2-34)，相矩阵(2-21)通 过旋转后得到的矩阵为： 哈尔滨工业大学工学硕士学位论文 - 19 - 21 12121 512512 2242 312312 512512 2242 312312 4 ()() cos2sin20 cos2 cos2sin2 cos2 cos2sin2 sin2 sin2cos2 sin2 cos2 sin2sin2 sin2 sin2cos2 sin2 cos2cos2 cos2 0sin li pli ppipi piipii pipi piipii piipii pipi piipii p 2416 2sin2ipip (2-35) 相矩阵中的元素 i p可以用 legendre 多项式展开得到 ( ) 0 (cos)(cos) l n l iell l pkp (2-36) 2.5.4 归一化归一化 相矩阵的归一化，主要是要求它的第一个元素满足归一化条件74， 即 11 4 4d (2-37) 上式中的积分项 11 即标量辐射传输方程中的相函数，与方程(2-18)中的 一致。定义散射截面为 11 2 4 1 s cm d k (2-38) 归一化的穆勒矩阵可以定义为 2 4 s m c k (2-39) 将式(2-39)代入式(2-14)中，得到入射和散射斯托克斯量之间的关系，用归 一化的穆勒矩阵表示为 2 22222 11 444 ss siiii c kc imiiidi k rk rr (2-40) 其中 2 s c d r 为与散射相关的离散立体角。 哈尔滨工业大学工学硕士学位论文 - 20 - 2.6 消光矩阵和发射系数向量消光矩阵和发射系数向量 从电磁理论的角度出发，当粒子受到电磁波照射的时候，除了吸收一部 分能量以外，还会由于振荡的电磁场作用而发生极化，从而感应出振荡的电 磁多级子。多级子又会产生电磁振荡，向各个方向发射电磁波。这就形成了 粒子的散射过程。因此我们能够用电磁理论来分析与计算辐射特性。 在分析辐射传输方程中的各项系数时，粒子的复折射率 ri mmim是一 个重要的参数。m的实部 r m与光在介质中的传播速度有关，虚部 i m与介质 吸收能量的多少有关。特别需要注意的是，对于某一特定材料的复折射率是 温度和波长的函数。只要知道了某种物质的复折射率，就可以运用电磁理论 计算它的辐射特性。在描述粒子的尺寸对辐射特性的影响时，常常用尺度参 数来表示，对于球形粒子，是通过球心的截圆周长与波长之比， 2 r kr ，其中2 /k 。 2.6.1 消光矩阵和发射向量消光矩阵和发射向量 一般除了多色散的情况，我们研究的粒子通常是非球形，非轴对称的， 有特定的转轴取向。对于这种一般的情况，矢量辐射传输方程(2-19)中的消 光矩阵可以写成如下的形式： vvhhvvhhvhhvvhhv vvhhvvhhvhhvvhhv hvvhhvvhvvhhhhvv hvvhvhhvvvhhvvhh re mmre mmre mmim mm re mmre mmre mmim mm re mmre mmre mmim mm im mm im mmim mmre mm e ag k k (2-42) 为特征矩阵。上式表明，衰减主要源于吸收，散射主要由粒子引起，吸收 主要由气体引起。上式子中的各项m定义为 2 2 m( , )n rs k (2-43) 上式中s表示散射振幅函数，上一章已经介绍。下标满足这样的法 哈尔滨工业大学工学硕士学位论文 - 21 - 则： vv2vh3hv4hh1 ,ssssssss。符号 表示表示该参量对粒径和转 轴取向取平均。 与相矩阵类似，在某些特定条件下 e k也可以作一定的简化 (1) 对角矩阵 000 000 000 000 e e e e e k k k k k (2-44) 其中 e k是标量辐射传输方程中使用的消光系数。该简化适用于：(1)由随机 导向轴对称粒子组成的各向同性介质（如球体和椭球体） ；(2)随机导向非对 称粒子，但一半数的粒子是另外一半的镜像； (2) 三元素矩阵 ( )( )00 ( )( )00 00( )( ) 00( )( ) ep pe e ec ce kk kk k kk kk (2-45) 其中 e k、 p k、 c k分别是消光系数，线偏振系数，圆偏振系数。为出射辐 射的天顶角余弦。它适用于：(1)在水平面上随机导向的非球形粒子；(2)非 球形轴导向粒子。 利用消光矩阵和相矩阵，发射(吸收)系数向量为 1111 4 1212 4 1313 4 1414 4 e e a e e kp d kp d k kp d kp d (2-46) 边边界条件和表面特性界条件和表面特性 包含表面发射和反射的矢量辐射传输方程的边界条件可以写成： 哈尔滨工业大学工学硕士学位论文 - 22 - 2 0 () ()(1) oo d o bbacdmc f iit iifririvri d (2-47) 其中上标, 和o分别表示从边界向内部，从内部到边界，从边界外部到内 部。 ba i和 c i代表背景辐射和准直辐射。 m i是i 的反射斯托克斯量。 d f代 表漫反射所占的比例， 为 dirac delta 函数。符号 0 ，t，r分别表示 发射向量，传输矩阵以及双向反射分布矩阵（对于包含大气顶部的准直和背 景辐射的气象领域的应用，传输矩阵一般视为单位矩阵） 。符号v为传播方 向与参考边界法线的夹角的余弦。式(2-47)表明，从给定的边界所发出的能 量的斯托克斯参数，是边界发射，背景和准直辐射，以及镜面反射与漫反射 的综合作用的结果。 兰贝特表面兰贝特表面 兰贝特（lambertian）表面，又称作漫射（diffuse）表面，是一种理想 的粗糙表面，其表面的发射与反射与方向无关。其反射辐射为非偏振的，双 向反射分布矩阵为 0 000 0000 0000 0000 r (2-48) 其中 0 为边界的半球反射率。对于兰贝特表面，其发射也是非偏振的，其 发射矢量为， 0 0 0 0 0 (2-49) 其中 0 为边界的半球发射率。其传输矩阵为， 0 000 0000 0000 0000 t t (2-50) 边界的半球反射率满足 000 1t 。 哈尔滨工业大学工学硕士学位论文 - 23 - 菲涅尔表面菲涅尔表面 菲涅尔（fresnel）表面，又称作镜面（specular）表面，是一种理想的 非常光滑的表面，其表面的反射辐射是偏振的且与角度有关，其双向反射 分布矩阵为， 2222 2222 * * 11 00 22 11 00 22 00re()im() 00im()re() vhvh vhvh vhvh vhvh rrrr rrrr r r rr r r rr r (2-51) 由菲涅尔定律可以得到垂直和平行的反射系数分别为， 222 222 coscos1 ( ) coscos1 v mm r mm (2-52) 22 22 coscos1 ( ) coscos1 h m r m (2-53) m 为表面的复折射率。菲涅尔表面的能量发射也是偏振，与角度有关的， 发射矢量为， 22 22 0 1 1 2 1 2 0 0 vh vh rr rr (2-54) 传输矩阵为， 2 2 * 12 21 000 000 00 00 v h r r t gg gg (2-55) 其中，下标*表示该矩阵适用于修正的斯托克斯量 * 1 cos re(1)(1) cos t vh i grr (2-56) 哈尔滨工业大学工学硕士学位论文 - 24 - * 2 cos im(1)(1) cos t vh i grr (2-57) 其中 i 为入射角， t 为折射角。这两个角的关系可以由斯涅尔折射定律 （snells law）得到， , , sinsin r t it r i m m (2-58) , r t m， , r i m分别为入射和出射的介质的复折射率的实部， ,r tr i mm。当大于 临界角 c 时， 12 0gg。其中 1 , sin / cr tr i mm (2-59) 2.7 辐射亮温辐射亮温 根据普朗克定律 1 5 2 exp/1 b c e ct (2-60) 其 中 b e为 光 谱 半 球 黑 体 辐 射 力 ， 842 1 0.59544 10 w m /mc ， 4 2 1.4388 10 m kc 。 辐射强度为： b b e i (2-61) 对于一个给定波长和温度的实际物体，其发射的能量要小于黑体。一 个实际物体的亮温定义是：与该实际物体发射出相同的能量的黑体的温度。 其定义式为： 2 5 1 ( , ) ln1/( , ) b c t ci (2-62) 其中(,)i 为该实际物体的辐 射强度。对于长波 或者高温的情况 ， 2/ 1ct，满足瑞利-琼斯假设，有 1 4 2 b ct e c (2-63) 得到该假设下的的亮温表达式, 哈尔滨工业大学工学硕士学位论文 - 25 - 4 2 , 1 ( , )( , ) b rj c ti c (2-64) 2.8 辐射热流辐射热流 根据 siegel 和 howell75的研究，定义离开一个单位体积的净辐射能量 为辐射热流矢量的散度， y xz q qq q xyz (2-65) 其中净辐射热流矢量为 4 0 n qiv d (2-66) 在式(2-66)中，下标 n 代表 x，y，z 坐标中的任意一个，v 代表辐射强度 i 的方向与坐标轴的夹角的余弦值。i 表示全波长的总的辐射强度，定义式 为， 0 ii d (2-67) 利用(2-66)，(2-67)，式(2-18)代入(2-65 中)， 44 000 44 00 () () (,) 4 eab s qkidki d k iddd (2-68) 对于矢量的辐射传输，将式(2-19)代入式(2-65)可以得到辐射热流的散度形 式，对于第一个斯托克斯参量 i，辐射热流的散度可以写成， 44 111213141 000 44 11121314 00 () 1 () 4 ieeeeab qkikqk uk v dk i d z iz qz uz v ddd (2-69) 另外三个量的辐射热流散度的表达式可以同理推出。 2.9 本本章小结章小结 本章从电磁场理论入门，介绍了偏振特性以及矢量辐射传输相关的基本 知识，比较了标量辐射传输方程与矢量形式的区别与联系。与标量的辐射传 哈尔滨工业大学工学硕士学位论文 - 26 - 输方程不同，矢量辐射传输方程中，用四个 stokes 量来描述辐射，研究它 们在传播过程中的吸收、发射和散射。与标量的辐射传输方程相比，矢量辐 射传输方程能够完全精确地描述含偏振的辐射在介质中的传播，计算结果更 加准确。由于四个参量的引入，它们的求解过程是相互耦合的，因此矢量辐 射传输的问题变得复杂，求解难度也增大了很多。 哈尔滨工业大学工学硕士学位论文 - 27 - 第第3章章 矢量辐射传输方程的有限体积法矢量辐射传输方程的有限体积法 目前国内外学者发展的，对矢量辐射传输求解的数值方法有很多，包括 倍加累加法，球谐函数法，谱元法，逐次散射法，离散坐标法等等。这些算 法各有千秋，但是目前还没有人提出使用有限体积法来求解矢量辐射传输问 题。 有限体积法（简称 fvm）与离散坐标法（简称 dom）比较类似，都是 通过选用某种差分格式将控制体中心的辐射强度与表面的联系起来，通过数 值积分近似散射相函数的积分，能够采用与流场计算协调的网格。此外， fvm 与 dom 最根本的区别在于，fvm 可以保证在每一个单元，每一个立 体角内的能量守恒。目前在辐射传输领域，对于使用 fvm 来求解标量形式 的辐射传输方程的研究已经很成熟9，这对于本文将其拓展到矢量辐射传输 方程的求解也很有帮助。 求解矢量辐射传输方程的矢量有限体积法（为了与求解标量方程的有限 体积法区分开来，称之为 vector finite volume method, vfvm）是基于传统 的 fvm，利用 stokes 矢量以及 mueller 矩阵求解的。通过对空间和角度的 离散，将矢量辐射传输方程转化成含有 4 个 stokes 量的方程组，然后对每 一个控制体单元和控制立体角单元积分得到。 3.1 笛卡尔坐标系下矢量辐射传输方程的有限体积法笛卡尔坐标系下矢量辐射传输方程的有限体积法 基于斯托克斯量的不连续叠加准则( the incoherent addition principle )， 对于准单色光的偏振辐射，在随机的粒子群介质中的三维稳态矢量辐射传输 方程可以描述为： 21 01 d ( , , ; , ) =-( , , ; , ) ( , , ; , )( , , ; , )( , , ) d ( , , ; , ;,) ( , , ;,) eab i x y z k x y zi x y zkx y zix y z s x y zi x y zdd (3-1) 由于( , , ; , )i x y z 为一个包含 4 个元素的列向量( , , )tii q u v，将式(3-1)展 开，可以得到由 4 个方程组成的方程组， 哈尔滨工业大学工学硕士学位论文 - 28 - 111213141 21 11121314 01 = 1 4 eeeeab di kikqk uk vk i ds iquv dd (3-2a) 212223242 21 21222324 01 = 1 4 eeeeab dq kikqkukvk i ds iquv dd (3-2b) 313233343 21 31323334 01 = 1 4 eeeeab du kikqkukvk i ds iquv dd (3-2c) 414243444 21 41424344 01 = 1 4 eeeeab dv kikqkukvk i ds iquv dd (3-2d) 有限体积法的基本思想是保证在每个立体角内其辐射能量守恒。在控制 体积 p v和控制立体角 m 内对辐射传输方程积分，可以得到 p v， m 内辐射 能量守恒方程的有限体积表达式。即在立体角 m 内从控制体 p v各表面进出 的辐射能量之和等于该控制体在立体角 m 内吸收、发射和散射的净辐射 能。 以第一个方程(3-2a)为例，积分得 111213141 21 11121314 01 ()= 1 4 mm cc mmmmmmm cjceeeeab aa m c ida dkikqk uk vk i iquv ddda d sn (3-3) 式子中 m c i 控制体环表面第 m 个立体角的辐射强度， 2 w/(m sr)； c a 控制体环表面积； m s 方向矢量，表示第 m 个立体角的中心； j n 控制体外表面单位法矢量。 假设在控制体积 p v和控制立体角 m 内辐射强度相等，边界节点的值可 哈尔滨工业大学工学硕士学位论文 - 29 - 近似由内节点 p 上的值表示，式(3-3)简化得， ,11,12,13,14,1, , , , , , 11121314 1 4 c m mmmmmm c jc jjeppeppeppeppa pb p j e w s n t b lllllm ppppp l a idkikqkukvki iquvv (3-4) 可以用某种空间差分格式将控制体表面上的辐射强度与控制体中心的辐 射强度关联起来。若使用阶梯格式，即 , d0 mmm c jpj ii (3-5) (3-4)式可写为 , , , =+b mmmmm ppjjp ie w s n t b a ia i (3-6) 下标, , ,je w s n t b表示与控制体 p 相邻的各个控制体的中心节点。 其中， 11, , , , , , max,0 mmm pjjepp j e w s n t b aa dkv (3-7) max,0 mm jjj aa d (3-8) 4 12,13,14, 11121314 b = 1 4 mmmm p peppeppepp lllllm ppppp l t kqkukv iquvv (3-9) 离散的边界条件为， 1111213141112 131411121314 1112131411121314 0 () () (1) i oooooo obbabababacc oo ccdmmmm dd cccciiiii v iit it qt ut vt it q t ut vfr ir qr ur v ff r ir qr ur vr ir qr ur vv (3-10) 另外三个量的推导同理可以求得。 3.2 程序的计算流程程序的计算流程 求解矢量辐射传输方程的有限体积法的基本思想和计算程序流程如下： 哈尔滨工业大学工学硕士学位论文 - 30 - (1) 设定各项基本参数，比如空间网格划分的密度，立体角划分的大小 等。 (2) 根据设定的数据，对天顶角和圆周角进行划分和编号，然后确定空间 立体角的编号。对空间网格进行离散。 (3) 根据模型的条件设置温度场，然后根据模型给定的条件设定 legendre 多项式展开的系数，利用 legendre 展开求得计算所需的散射相矩阵 中的各个元素。 (4) 对所求的参数进行初始化，计算所有点，各个面的方向权值 m j d，准 备开始迭代。 (5) vfvm 计算主程序对于 nl 个空间立体角，ix-2 个节点进行迭代。首 先根据计算的方向确定迭代的顺序，根据所给的模型，利用式(3-10) 设定边界条件的值，然后利用初始化后的数值代入式(3-4)中，求得第 一次迭代，第一个节点，第一个立体角的 i 值，同理求得第一次迭 代，第一个节点，第一个立体角的 q，u，v 值。然后对 nl 个空间 立体角，ix-2 个节点进行循环迭代，求出第一次迭代后，全场的每 一个节点，每一个离散的空间立体角的 i，q，u，v 的数值。 (6) 检验计算出的数值是否满足收敛条件，若不满足，则回到第(5)步， 进行下一次迭代。用该次计算出的 i，q，u，v 数值作为下一次检验 的标准，直到最终迭代收敛为止。 (7) 输出所要求的参数：强度，亮温，热流密度等等。 利用 vfvm 求解矢量辐射传输方程的计算程序流程图如图 3-1 所示： 哈尔滨工业大学工学硕士学位论文 - 31 - 开 始 call grid 角度离散和空间离散 call phasecompute 散射相矩阵中各个元素的计算 其他各项物性参数的计算 call calcr 开始计算 计算i,q,u,v四个量 四个参数是否均是否 满足收敛条件？ 否 call output 输出计算结果（强度，亮温，热流等） 计算结束 是 图 3-1 vfvm 程序计算流程图 哈尔滨工业大学工学硕士学位论文 - 32 - 3.3 计算和分析计算和分析 为了验证所编制程序的准确性，将针对多种条件下，多个算例与文献中 的结果作对比，证明 vfvm 的准确性。 3.3.1 菲涅尔表面的发射模型菲涅尔表面的发射模型 如图 3-2 所示，算例 1 为一个 300k 的无限大的平静水面（视为镜面， 即菲涅尔表面，负折射率为 3.724-2.212i） ，水面上方无气体，求频率为 85.5ghz的单色上行（upwelling）辐射强度。 水面 300 k m = 3.724-2.212i 图 3-2 算例 1 模型示意图 由于没有参与性气体，使得该问题大大简化，方程(2-18)等式右边的数 值全为 0。因此任意点和任意方向的辐射强度取决于边界的发射，在任意的 方向上是一个常量，并不衰减，该问题转换成一个边界问题。根据式(2-47) 得到该问题的边界条件： o b ii (3-11) 由于模型是镜面的，根据式(2-54)知道 o 的后两项为 0，则对于该问题， u，v 的数值恒为 0。利用式(2-64)，将辐射强度写成辐射亮温的形式， bo tt (3-12) 由于模型是一维的，计算结果只是天顶角余弦的函数。天顶角网格划分为 50 份，将计算的结果与 evans 和 stephens40的计算结果作对比。 哈尔滨工业大学工学硕士学位论文 - 33 - 0.00.20.40.60.81.0 60 80 100 120 140 160 180 tb,i(k) vfvm ref.40 (a) ,b i t 0.00.20.40.60.81.0 0 20 40 60 80 100 120 tb,q(k) vfvm ref.40 (a) ,b q t 图 3-3 算例 1 的计算结果 如图 3-3 所示，vfvm 的计算结果与文献40几乎完全吻合，证明了模 型的程序的正确性。该算例也很好地说明了镜面的起偏作用。 哈尔滨工业大学工学硕士学位论文 - 34 - 3.3.2 等温均质吸收发射性气体（无散射）模型等温均质吸收发射性气体（无散射）模型 300 k m = 3.724-2.212i 8 km ka1 = ke = 0.01 km-1 85.5ghz 图 3-4 算例 2 模型示意图 如图 3-4 所示，算例 2 的下表面使用与算例 1 一样的镜反射面，在镜面 上方有一层8 km厚的均质、均温、吸收发射性气体，无散射作用，上表面 为透明边界。气体温度为 300k。气体分子为随机导向的，形状任意，各向 同性的，其消光矩阵为对角阵，且只含有一个参数 e k。由于无散射，发射 向量只有第一项 -1 1 0.01km ae kk，其他项都为 0。evans 和 stephens40指 出，对于这样的问题，u 和 v 都为 0，这样所有的矩阵都简化成2 2，只需 要求解 i 和 q 的值。使用 vfvm 计算的结果与文献67中第五章 test 2 的 结果作对比，计算0 kmy 处的下行（downwelling）辐射亮温 ,b i t，以及 0 kmy 处的频率为 85.5ghz 的单色上行辐射亮温（upwelling） ,b i t、 ,b q t。在高度方向上的计算网格采取均匀划分，数量为 80 份，空间立体角 的划分为50 4nn。如图 3-5 所示，计算结果与参考文献基本一致。 哈尔滨工业大学工学硕士学位论文 - 35 - 0.00.20.40.60.81.0 180 200 220 240 260 280 300 tb,i(k) vfvm ref.67 upwelling at y=8km 0.00.20.40.60.81.0 0 10 20 30 40 50 60 70 tb,q(k) vfvm ref.67 upwelling at y=8km 0.00.20.40.60.81.0 0 50 100 150 200 250 300 tb,i(k) vfvm ref.67 downwelling at y=0km 图 3-5 算例 2 的计算结果 哈尔滨工业大学工学硕士学位论文 - 36 - 3.3.3 双层参与性大气（含散射）模型双层参与性大气（含散射）模型 t =300k 4 km t = 273k t = 245k 4 km 1 0.15224 km er k 1 0.13536 km ei k 0.38175 r 0.98190 i 雨 冰 r p i p 2.7 ba tk 85.5ghz 图 3-6 算例 3 模型示意图 如图 3-6 所示，算例 3 的下表面使用与算例 1 一样的反射镜面，在镜面 上方有一个 8km 厚，双层的参与性大气（含散射） ，上表面为透明边界。其 中下层的大气模型为雨，上层为冰。在双层介质中，温度呈线性分布，上边 界，冰水分界面，以及下边界的温度如图所示。由于大气是平板平行的， 且只有辐射热源，u 和 v 项恒为 0，所以散射相矩阵简化成一个2 2的矩 阵。两层雨滴均为球形粒子、且遵循 marshall-palmer 分布，雨层的降雨速 率为 0.5mm/h，冰层的降雨速率为 2.0mm/h。两层中粒子的最大的粒径为 1cm。上表面为非偏振的漫黑体的背景辐射，其等效的辐射亮温为 2.7k。两 层大气的衰减系数和单散射反照率已经由 evans 和 stephens40给出，如图所 示。由于粒子为球形的，由(2-44)可知，两层的消光矩阵均为对角阵，分别 只含有一个参数 e k，散射相矩阵的计算可以由 legendre 多项式展开得到， 展开系数已经由 evans 和 stephens40计算得到，如表格 3-1 所示。 哈尔滨工业大学工学硕士学位论文 - 37 - 表 3-1 散射相矩阵中的三个元素的 legendre 展开系数40 l 雨 层 冰 层 1 l 2 l 3 l 1 l 2 l 3 l 0 1.000000 -0.378560 0.122149 1.000000 -0.203665 0.706712 1 0.365211 -0.082115 1.482953 1.305650 -0.111353 1.669173 2 0.518055 0.353963 0.356369 0.915656 0.176185 0.856457 3 0.115569 0.077666 0.067922 0.348084 0.097906 0.357357 4 0.032699 0.023756 0.008445 0.131959 0.025961 0.114952 5 0.006058 0.004312 0.001049 0.032286 0.012476 0.031836 6 0.001139 0.000819 0.000071 0.011950 0.001494 0.009323 7 0.000188 0.000133 0.000000 0.001923 0.000916 0.002046 8 0.000030 0.000021 -0.000002 0.000899 0.000030 0.000561 9 0.000005 0.000006 0.000000 0.000085 0.000052 0.000112 10 0.000061 -0.000004 0.000025 11 0.000003 0.000004 0.0