《基于时域多分辨算法的非球形气溶胶散射特性仿真模拟-胡帅.pdf》由会员分享，可在线阅读，更多相关《基于时域多分辨算法的非球形气溶胶散射特性仿真模拟-胡帅.pdf（18页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、基于时域多分辨算法的非球形气溶胶散射特性仿真模拟胡帅高太长李浩杨波江志东陈鸣李书磊Simulating scattering properties of nonspherical aerosol particles using multiresolution timedomainmethodHu Shuai Gao Tai-Chang Li Hao Yang Bo Jiang Zhi-Dong Chen Ming Li Shu-Lei引用信息Citation: Acta Physica Sinica , 66, 044207 (2017) DOI: 10.7498/aps.66.044207在
2、线阅读View online: http:/dx.doi.org/10.7498/aps.66.044207当期内容View table of contents: http:/ you may be interested in一种利用布里渊增益谱边带解调提高布里渊光时域反射系统测温精度的方法Temperature measurement accuracy enhancement in the Brillouin optical time domain reflectometry sys-tem using the sideband of Brillouin gain spectrum demo
3、dulation物理学报.2016, 65(24): 244203 http:/dx.doi.org/10.7498/aps.65.244203非球形椭球粒子参数变化对光偏振特性的影响Research of the influence of non-spherical ellipsoid particle parameter variation on polarization charac-teristic of light物理学报.2016, 65(6): 064205 http:/dx.doi.org/10.7498/aps.65.064205双Pearcey光束的构建及数学机理研究Con
4、struction of Bi-Pearcey beams and their mathematical mechanism物理学报.2016, 65(21): 214208 http:/dx.doi.org/10.7498/aps.65.214208沙尘大气电磁波多重散射及衰减Multiple scattering and attenuation for electromagnetic wave propagation in sand and dust atmosphere物理学报.2016, 65(9): 094205 http:/dx.doi.org/10.7498/aps.65.094
5、205Pearcey光束簇的实验产生和光学结构研究Generation of a family of Pearcey beams and their optical structure物理学报.2015, 64(23): 234205 http:/dx.doi.org/10.7498/aps.64.234205万方数据物理学报Acta Phys. Sin. Vol. 66, No. 4 (2017) 044207基于时域多分辨算法的非球形气溶胶散射特性仿真模拟胡帅1)高太长1)2)y李浩1)杨波2)江志东3)陈鸣1)李书磊1)1)(解放军理工大学气象海洋学院,南京211101)2)(解放军理工
6、大学,电磁环境效应与电光工程国家级重点实验室,南京210007)3)(海军航空工程学院,青岛266041)(2016年9月19日收到; 2016年11月30日收到修改稿)非球形气溶胶的散射特性是影响辐射传输模拟准确性的重要因素.为实现非球形、非均质气溶胶散射特性的模拟,基于MRTD (multi-resolution time-domain)方法建立了一个新的气溶胶散射模型.采用MRTD技术实现了近场电磁场的计算;考虑气溶胶的特殊性,推导了基于体积积分方法的近远场外推方法,实现了粒子散射振幅矩阵和穆勒矩阵的仿真;构建了粒子吸收和消光截面的计算模型,实现了粒子积分散射特性的高精度模拟.将MRTD
7、散射模型的结果与Mie理论、T矩阵法进行了对比,验证了模型的准确性;讨论了空间网格粗细对模拟精度的影响,并定量分析了模型的运行效率.结果表明, MRTD散射模型的相函数模拟误差小于8%,其中前向散射方向小于4%;当粒径与入射光波长相当时,消光和散射效率因子的相对误差小于0.1%;空间网格粗细对模拟精度影响显著,当粒子尺度参数小于20时,在相同模拟精度要求下,所需网格尺寸随尺度参数呈先增大后减小的特征.关键词:非球形气溶胶,散射特性,时域多分辨算法,穆勒矩阵PACS: 42.68.Mj, 42.68.Ay, 95.30.Jx, 42.25.Bs DOI: 10.7498/aps.66.04420
8、71引言非球形气溶胶散射特性的不确定性一直是制约辐射传输和气候模拟精度的重要因素1;2.对比多年的IPCC报告可发现3,气溶胶直接辐射强迫的评估精度虽然有所提高,但仍存在较大的不确定性(0.90.1 W/m2),其中气溶胶和卷云中冰晶等非球形粒子的散射特性就是制约其评估精确性的一个重要原因4;5.气溶胶光学遥感及大气校正等均需要以辐射传输模式为基础6;7,而无论是标量还是矢量辐射传输模式,均需要气溶胶散射特性作为输入;目前考虑多次散射的辐射传输求解模型如DISORT (discrete-ordinate-method),RT3/PolRadtran及SHDOM (spherical harmo
9、nicsdiscrete ordinate method)等已得到广泛检验8;9,但作为模式输入参数的气溶胶散射特性,如穆勒矩阵和散射相函数等10,仍主要通过Mie散射计算得到,并未充分考虑气溶胶及冰晶的非球形效应,而实际非球形粒子与等效球形粒子的散射特性存在较大差异11;12,这将影响遥感正演模型的准确性,进而影响反演精度,因此准确获取非球形气溶胶的散射特性,提高正演模型模拟精度已成为发展大气遥感技术的关键环节之一13;14.为模拟非球形气溶胶的散射特性,目前已经建立了许多气溶胶散射模型,大致可分为两大类:一类是近似计算模型,主要针对粒子尺度参数较为极端的情形(当入射光波长远大于粒径或远小国
10、家自然科学基金(批准号: 41575025, 41575024)资助的课题.通信作者. E-mail: 2017中国物理学会Chinese Physical Society http:/044207-1万方数据物理学报Acta Phys. Sin. Vol. 66, No. 4 (2017) 044207于粒径的情形),主要的模型包括Rayleigh近似、Rayleigh-Gans-Stevenson近似(RGA)、异常衍射近似和几何光学近似GOA (geometric optics ap-proximation)等15,其中GOA方法广泛用于大尺度冰晶粒子的散射特性计算14;16;第二类是
11、精确求解模型,这类模型通过在一定边界条件下数值求解Maxwell方程组或亥姆霍兹方程,获取模拟粒子电磁散射特性15.按照计算原理,精确求解模型可大致分为三类:基于基函数展开的散射模型、基于体积积分方程的散射模型及基于微元法的散射模型.基于基函数展开的散射模型的特点是基于矢量波函数对电磁入射场、内场及散射场进行展开,然后求解微分方程或积分方程建立入射场与散射场展开系数间的关系,进而实现对散射过程的求解.这一类模型主要包括T矩阵法(T ma-trix method)17, EBCM (extended boundary con-dition method)18, SVM (separation o
12、f variablesmethod)19和PMM (point matching method)20等.这类方法具有计算精度高、速度快、内存消耗小的优点,但当粒子尺度参数较大、形状较为极端或复折射率较大时,易出现数值求解的不稳定;散射场及内部场采用基函数展开时,展开系数存在截断误差;此外,受求解边界条件的限制,这类模型只能模拟一些形状简单的粒子,如椭球粒子、圆柱粒子等.体积积分方程法是通过数值求解亥姆霍兹方程的格林函数形式解来模拟散射场,典型代表是矩量法MoM (method of moments)21, DDA (dis-crete dipole approximation)22;23和FI
13、EM (fred-holm integral equation method)积分方程法. DDA是MoM的改进,其本质差异在于自身感应电场的计算,可适用于任意形状、非均质的气溶胶粒子,但是模拟过程涉及大型矩阵方程的求解,存在不稳定性;同时,计算精度取决于偶极子的数目,计算精度与内存消耗、收敛速度存在突出矛盾.因此,该方法主要适合尺度参数较小的粒子的散射特性模拟,且计算精度不高. FIEM积分方程法的优点是积分过程仅和粒子尺度参数和形状有关,而与入射光的入射角和偏振状态无关,但只适用于尺度和复折射率较小的规则粒子,且模型计算量很大15.基于微元法思想的散射模型主要包括时域有限差分法(nited
14、ierence time domain, FDTD)15;16和有限元法(nite element method, FEM)24.其中FDTD已广泛应用于气溶胶散射特性模拟16,由于该方法避免了积分方程法的奇异核问题,因此更适用于复杂形状、非均质气溶胶的光散射特性模拟.但与DDA相似,该方法存在计算精度与计算量的矛盾.为解决上述问题, Liu等25将时域伪谱法(pseudospectral time domain, PSTD)引入气溶胶散射过程模拟中,由于该方法取相对较少的采样点便可获得较高模拟精度,因此可显著降低内存消耗与运行时间.目前, Liu等26已采用该方法实现了尺度参数达200的粒子
15、散射特性模拟,但受原理限制,该方法只能采用纯散射场技术引入入射波,因此在计算不同形状粒子的散射特性时,相应程序也要进行修改,过程繁琐. MRTD是矩量法和FDTD的结合体27;28,是Maxwell方程的高阶差分展开,具有良好数值色散特性,也广泛用于大目标的散射特性的模拟29.基于此,本文将其引入非球形气溶胶的散射特性的模拟中,结合研究对象和学科需求的特殊性,建立了一套新的非球形气溶胶散射模型.本文思路如下:首先介绍了一些物理概念和MRTD散射模型的基本框架,然后简单介绍了基于MRTD近场计算方案,包括基本离散迭代格式、吸收边界设计及总场散射场技术;针对气溶胶的特殊性,推导了基于体积积分的近远
16、场外推方案和粒子散射特性计算方法;最后对模型的准确性和性能进行了分析和验证.2基本概念2.1 Maxwell旋度方程组及电磁参数转换电磁场传播规律可用Maxwell旋度方程组描述.无源条件下,该方程的矢量表达形式如下:H = Et + E; (1)E = Ht mH; (2)上式中, E和H分别表示电场和磁场的强度矢量;表示介质的介电常数; 表示磁导系数; 和 m分别表示介质的电导率和磁导率.在大气科学领域,粒子的光学属性通常采用复折射率m表示,其实部和虚部分别表征了粒子散射效应和吸收效应的强弱.粒子复折射率m与复介电常数c存在以下转换关系:c = 0m2 = j !; (3)044207-2
17、万方数据物理学报Acta Phys. Sin. Vol. 66, No. 4 (2017) 044207若假设复折射率m = mr jmi,则有关系式 = 0(m2r m2i), = 2!mrmi0成立.采用以上关系便可将粒子的复折射率转化为电磁散射模拟所需的电磁参数.2.2穆勒矩阵和Stokes矢量穆勒矩阵F(rsca;rinc)是辐射传输模式的重要参数,它不仅表征了散射能量的空间分布,同时也定量描述了入射光Stokes矢量Sinc和散射光Stokes矢量Ssc转换关系.穆勒矩阵F的定义如(4)式所示30,其中F11也称为散射相函数.Ssc(rsca)= 1k2r2F(rsca;rinc)S
18、inc= 1k2r2F11 F12 F13 F14F21 F22 F23 F24F31 F32 F33 F34F41 F42 F43 F44IincQincUincVinc; (4)式中rsca;rinc分别表示散射光观测方向以及入射光方向, k表示入射光波数, r为观测点距离散射体中心的长度;偏振光Stokes矢量定义为S = (I;Q;U;V)T,其中I表示光强, Q为水平与垂直偏振分量之差, U为45方向与135方向偏振分量之差, V表示左旋与右旋偏振分量之差.3 MTRD非球形气溶胶散射模型3.1模型的基本框架MRTD非球形气溶胶散射模型包括三个模块:电磁场近场计算模块、近远场外推模块
19、及气溶胶散射特性计算模块.电磁场近场计算模块主要实现包含散射体空间的电磁场模拟,其中关键技术包括Maxwell方程离散迭代格式、吸收边界的设计以及基于总场散射场技术的入射波引入;近远场外推模块的功能是将近场电磁场变换至频域,并实现近远场场变换;气溶胶散射特性计算模块主要实现粒子吸收、消光和散射截面及穆勒矩阵等参数的模拟.由于气溶胶一般为电介质,本文基于介质中的Maxwell方程组推导了基于体积积分远场变换方法;在粒子散射特性计算方面,本文从大气辐射学需求出发,推导了粒子的吸收截面和消光截面计算模型.气溶胶散射特性的计算粒子消光、吸收和散射截面粒子穆勒散射矩阵输入量设计入射光偏振特性设计入射光入
20、射方向设计粒子尺度、形状及复折射率设置MRTD 非球形气溶胶散射模型散射场近远场外推近场的振幅和相位提取（频域）基于体积积分原理的三维时谐场外推电磁场近场的计算迭代网格的设计与稳定性分析Maxwell方程的离散与迭代格式建立总场散射场技术（入射波引入）ADE-PML吸收边界的设计与实现图1 MTRD非球形气溶胶散射模型基本框架Fig. 1. The basic frame of MRTD scattering model for nonspherical aerosol particles.MRTD散射模型的计算域是在直角坐标系下进行离散的,电磁场空间离散格式采用Yee元胞.由于散射体形状不规
21、则,在不同介质的交界处,单个Yee元胞内可能存在多种介质混合的现象,若直接对其进行截断,则可造成较大的“阶梯误差”4,因此需要对元胞内电磁参数进行等效平均.本模型采用Maxwell-Garnett公式进行计算4,具体计算方法如(5)式所示.r(i;j;k)1r(i;j;k) + 2= 1xyzcell(i;j;k)r(x;y;z)1r(x;y;z) + 2dxdydz; (5)式中, r(x;y;z)表示位置(x;y;z)处的相对介电常数,表示复介电常数c与0之比, r是元胞的平均相对介电常数值.044207-3万方数据物理学报Acta Phys. Sin. Vol. 66, No. 4 (2
22、017) 0442073.2基本迭代格式的建立和稳定性条件在直角坐标系下,分别采用Daubechies尺度函数i(x)与矩形脉冲函数hn(t)在时间和空间上对电磁场分量进行展开31:Ex(r;t) =+1i;j;k;n= 1Ex;ni+1/2;j;ki+1/2(x)j(y)k(z)hn(t); (6)Ey(r;t) =+1i;j;k;n= 1Ey;ni;j+1/2;ki(x)j+1/2(y)k(z)hn(t); (7)Ez(r;t) =+1i;j;k;n= 1Ez;ni;j;k+1/2i(x)j(y)k+1/2(z)hn(t); (8)Hx(r;t) =+1i;j;k;n= 1Hx;n+1/2
23、i;j+1/2;k+1/2i(x)j+1/2(y)k+1/2(z)hn+1/2(t); (9)Hy(r;t) =+1i;j;k;n= 1Hy;n+1/2i;j+1/2;k+1/2i+1/2(x)j(y)k+1/2(z)hn+1/2(t); (10)Hz(r;t) =+1i;j;k;n= 1Hz;n+1/2i+1/2;j+1/2;ki+1/2(x)j+1/2(y)k(z)hn+1/2(t); (11)式中r=(x;y;z)为场分量的位置向量; Ex;ni+1/2;j;k,Ey;ni;j+1/2;k, Ez;ni;j;k+1/2为电场分量的展开系数;Hx;n+1/2i;j+1/2;k+1/2, H
24、y;n+1/2i;j+1/2;k+1/2和Hz;n+1/2i+1/2;j+1/2;k表示磁场分量的展开系数; i, j, k分别表示场分量的空间坐标; n表示时间迭代步数. hn(t)满足正交性,其定义如下:hn(t) = h(t/tn); (12)h(t) =1 |t| 1/2:(13)Daubechies尺度函数i(x)定义如下:i(x) = (x/xi+M1); (14)式中M1为Daubechies函数(x)的一阶矩.尺度函数i(x)不仅满足正交性,同时还近似满足位移内插特性,即i(x) i;0, i;i为狄拉克函数.简单而不失一般性,以电场分量Ex为例,对Maxwell离散迭代格式进
25、行简单说明.将(6)式入Maxwell旋度方程,可得+1i;j;k;n= 1Ex;ni+1/2;j;ki+1/2(x)j(y)k(z)hn(t)t+1i;j;k;n= 1Ex;ni+1/2;j;ki+1/2(x)j(y)k(z)hn(t)=+1i;j;k;n= 1Hz;n+1/2i+1/2;j+1/2;ki+1/2(x)j+1/2(y)y k(z)hn+1/2(t)+1i;j;k;n= 1Hy;n+1/2i+1/2;j;k+1/2i+1/2(x)j(y)k+1/2(z)z hn+1/2(t): (15)采用Galerkin法则,在方程两侧同时乘以i+1/2(x)j(y)k(z)hn+1/2(t
26、),并求取对应项的内积,结合函数hn(t)和i(x)的正交性, (15)式可化简为xyz+1i;j;k= 1(Ex;n+1i+1/2;j;k Ex;ni+1/2;j;k) i;i j;j k;k+ xyz2+1i;j;k= 1(Ex;n+1i+1/2;j;k +Ex;ni+1/2;j;k) i;i j;j k;k= txz+1i;j;k= 1Hz;n+1/2i+1/2;j+1/2;k i;i k;kLs 1l= Lsa(l) j+l;j044207-4万方数据物理学报Acta Phys. Sin. Vol. 66, No. 4 (2017) 044207txy+1i;j;k= 1Hy;n+1/
27、2i+1/2;j;k+1/2 i;i j;jLs 1l= Lsa(l) k+l;k; (16)式中, a(l)称为连接系数, l的变化范围为LsLs 1, Ls = 2N 1称为基函数的有效支撑尺寸,N为小波函数消失矩,文中取N = 2. a(l)可由尺度函数在傅里叶频域的内积计算得到:a(l) dj+1/2(x)dx ;j l(x)= 10!|(!)|2 sin(!(l+ 1/2)d!; (17)式中|(!)|为尺度函数的傅里叶变换.整理并化简(16)式可得电场分量的迭代步进方程:Ex;n+1i+1/2;j;k= CA(m)Ex;ni+1/2;j;k +CB(m)Ls 1l= Lsa(l)(
28、Hz;n+1/2i+0:5;j+l+0:5;k/yHy;n+1/2i+0:5;j;k+l+0:5/z); (18)式中, m = (i+ 1/2, j, k)为场分量的网格离散坐标, CA和CB的计算方法如下:CA(m) = 2 t2+ t;CB(m) = 2t2+ t: (19)由于函数i(x)的位移内插特性,空间r0的电磁场值就等于该处的场展开系数,也即有Ex(r0;t0) = Ex;ni+1/2;j;k成立,因此可避免对场分量的重构,这是采用Daubechies尺度函数作为基函数展开的优势所在.同理,电场分量Ey, Ez和磁场分量的MRTD迭代步进方程也有类似形式.由于MRTD是Maxw
29、ell方程的一种显式差分格式,因此为保证迭代过程的收敛性,时间差分间隔t与网格差分步长必须满足Courant稳定性条件,其中对于基于Daubechies尺度函数的MRTD算法,时间迭代步长t必须满足以下关系式31 :t 1cLs 1l=0|a(l)|1(x)2 +1(y)2 +1(z)2;(20)式中, c为光速; x, y和z分别为x, y和z方向的离散步长.3.3 ADE-PML吸收边界采用基于辅助微分方程的PML边界(per-fectly matched layer with auxiliary dierentialequation, ADE-PML)作为吸收边界.本文在刘亚文等32提出
30、的二维ADE-PML边界基础上,将其推广至三维情形.下面仍以Ex分量为例,介绍吸收边界的处理过程.在频域情形下,坐标伸缩的Maxwell旋度方程组如下:j!E + E = s H; (21)j! H mH = s E; (22)其中s = 1sxxx+1syyy+1szzz, si (i = x, y或z)为坐标伸缩因子,具体形式如下:si = i + i i +j!0; (23)求取si的倒数可得1si= 1 i+ ii +j!i,其中i = i, i = 0 2i, i = 2i i + i i, i和 i为非负实数,且 i 132.将其代入(21)式Ex分量的计算公式中,并引进辅助变量
31、exy和 exz可得:1yHzy 1zHyz + exy exz= j!Ex + Ex; (24)exy = yy +j!yyHz;exz = zz +j!zzHy: (25)进一步将上式改写为时域形式:1yHzy 1zHyz + exy exz= tEx + Ex; (26)y t exy + y exy = y yHz;z t exz + z exz = z zHy: (27)参照3.2节的处理方法,采用Galerkin法对三个偏微分方程进行离散,可得场分量Ex;n+1i+1/2;j;k的迭代步进方程:044207-5万方数据物理学报Acta Phys. Sin. Vol. 66, No.
32、 4 (2017) 044207Ex;n+1i+1/2;j;k = CAi+1/2;j;kEx;ni+1/2;j;k +CBi+1/2;j;k 1yyLs 1l= Lsa(l)Hz;n+1/2i+1/2;j+l+1/2;k 1 zzHy;n+1/2i+1/2;j;k+l+1/2+CBi+1/2;j;k(x;n+1exy;i+1/2;j;k x;n+1exz;i+1/2;j;k); (28)x;n+1exy;i+1/2;j;k = by x;nexy;i+1/2;j;k + ayyLs 1l= Lsa(l)Hz;n+1/2i+1/2;j+l+1/2;k; (29)x;n+1exz;i+1/2;j
33、;k = bz x;nexz;i+1/2;j;k + azzLs 1l= Lsa(l)Hy;n+1/2i+1/2;j;k+l+1/2: (30)以上各式中, ai和bi(i = x;y;z)满足以下公式:bi = 2i it2i + it; ai =2 it2i + it: (31)ADE-PML吸收边界通常采用PEC边界进行截断.在边界处,由于MRTD迭代过程涉及边界外的场分量,因此为保证边界附近场分量的正常迭代,需要采用镜像原理对边界外的场分量进行拓展.3.4 MRTD的总场-散射场技术与FDTD模型类似,采用总场散射场技术引入平面入射波.在该方案中,计算域被分为总场区和散射场区,总场区电
34、磁分量为散射场Esca和Hsca与入射场Einc和Hinc之和(如(32)式),散射场区仅包含散射场Esca和Hsca 33.E = Einc +Esca; H = Hinc +Hsca: (32)由MRTD的步进迭代方程可知,在计算连接边界附近节点的场分量时,需要同时使用散射场和总场区的切向分量,因此需要通过引入入射波对相应计算节点进行订正,考虑各场分量迭代公式的相似性,下面仅以电场分量Ex为例进行简单说明.图2所示为Ex分量的总场和散射场节点订正区域,其中图2(a)为散射场需要进行订正的节点区域,共4个;图2(b)为总场区需要订正的节点区域,共8个.考虑各区域订正方法的基本一致,仅以个散射
35、场区域1和总场区域5作为例子进行解释.区域1的范围为: i Imin;Imax 1, j Jmin;Jmax, k Kmax + 1;Kmax + Ls 1.由于该区域为散射场区,在迭代计算过程中,所涉及的总场节点需要扣除入射场的值,因此,该区域内各元胞Ex的计算方法如下:sEx;n+1i+1/2;j;k = sEx;n+1i+1/2;j;k +tyLs 1l= Lsa(l)sHz;n+1/2i+1/2;j+l+1/2;k tz Ls 1l=Kmax ka(l)sHy;n+1/2i+1/2;j;k+l+1/2ImaxIminJmaxJminKminKminKmax KmaxJmaxJminIm
36、axImin123456789101112567891011121234xyz(a) 散射场节点订正区间分布(b) 总场节点订正区间分布图2 Ex分量的总场和散射场节点订正区域Fig. 2. Correction area for the Yee cell distributed in scattering eld and toal eld.044207-6万方数据物理学报Acta Phys. Sin. Vol. 66, No. 4 (2017) 044207+Kmax k 1l= Lsa(l)(tolHy;n+1/2i+1/2;j;k+l+1/2 incHy;n+1/2i+1/2;j;k+l
37、+1/2): (33)区域5的范围为i Imin;Imax 1, j Jmin Ls+ 1;Jmin 1, k Kmax Ls+ 1;Kmax;由于该区域为总场区,计算过程中所涉及的散射场节点需要增加入射场的值.因此,该区域内各节点的计算公式可修正为:tolEx;n+1i+1/2;j;k = tolEx;n+1i+1/2;j;k +ty Ls 1l= Jmin ja(l)tolHz;n+1/2i+1/2;j+l+1/2;k +Ls 1l= Jmin ja(l)(sHz;n+1/2i+1/2;j+l+1/2;k+ incHz;n+1/2i+1/2;j+l+1/2;k) tz Ls 1l=Kmax
38、 ka(l)tolHy;n+1/2i+1/2;j;k+l+1/2+Kmax k 1l= Lsa(l)(sHy;n+1/2i+1/2;j;k+l+1/2 + incHy;n+1/2i+1/2;j;k+l+1/2): (34)以上各式中,下标s, tol和inc分别表示该场分量为散射场、总场和入射场分量.连接边界附近各节点的入射波切向分量采用一维MRTD投影至三维总场边界的方案计算.该方法基本流程与FDTD相似,其中一维MRTD和三维MRTD的时间步长必须相同.3.5基于体积积分的远场计算模型气溶胶和冰晶一般为非导磁电介质,因此在介质内,存在电极化矢量P(r),而磁极化矢量M(r) = 0.取电磁
39、波随时间的变化因子为exp(j!t),则在频域中,介质中的Maxwell方程可改写为34H = j!0(E +P(r);E = j! 0H; (35)上式中P(r) = E(r) = (r)/0 1)E(r), 为电极化率.对(35)式第一个方程两端取旋度,并将第二个方程代入可得E= (j! 0H) = !20 0(E +P(r); (36)采用矢量恒等式E = (E) E进一步化简可得到( +!20 0)E= (!20 0I +)P(r); (37)式中为拉普拉斯算子,由数学知识,该方程的格林函数形式解应满足(38)式和(39)式.E(r) = Einc(r) +VG(r;r)(!20 0I
40、 +)P(r)d3r; (38)G(r;r) +!20 0G(r;r)= (rr): (39)其中格林函数的具体形式如下:G(r;r) = exp(jk|rr|)4 |rr| ; (40)其中k = !0 0.将格林函数的具体形式代入,远场的电场矢量可写为E(r) = Einc(r) +Vexp(jk|rr|)|rr|(k2I +)P(r)d3r: (41)在远场,取近似|r r| = rr er,并将上式写成频域形式,可得散射场的表达形式如下:Esca(r) = E(r)Einc(r)= k2 exp(jkr)4 rV(r)0 1E(r)er er E(r)exp(jker r)d3r: (
41、42)由于MRTD计算区域是离散的, Esca需要采用数值计算方案求解,本文采用以下方法计算:Exsca(r) = k2 exp(jkr)4 rijk (i;j;k)0 1044207-7万方数据物理学报Acta Phys. Sin. Vol. 66, No. 4 (2017) 044207 Ex(i;j;k)sin obcosob Tsexp(jkx ix+jky jy+jkz kz)xyz; (43)Eysca(r) = k2 exp(jkr)4 rijk (i;j;k)0 1 Ey(i;j;k)sin obsinob Tsexp(jkx ix+jky jy+jkz kz)xyz; (44
42、)Ezsca(r) = k2 exp(jkr)4 rijk (i;j;k)0 1 Ez(i;j;k)cos ob Tsexp(jkx ix+jky jy+jkz kz)xyz; (45)式中( ob;ob)为观测方位角; 为Yee元胞的平均介电常数; Ex, Ey和 Ez为电场分量在坐标(i;j;k)的平均值(采用相邻两个Yee元胞平均);k = ker = (kx, ky, kz)为波矢量. Ts为中间变量,如下式所示:Ts = sin obcosob Ex(i;j;k)+sin obsinob Ey(i;j;k)+cosob Ez(i;j;k): (46)进一步将直角坐标系下的场分量变换至
43、球坐标系下,可得Esca(r)E sca(r)Ersca(r) = UTExsca(r)Eysca(r)Ezsca(r)=sini cosi 0cos i cosi cos i sini sin isin i cosi sin i sini cos iExsca(r)Eysca(r)Ezsca(r): (47)值得注意的是, Ersca(r)计算结果应为0,在计算过程,可用该结论验证计算过程的准确性,一般情况下也可省略其计算.3.6粒子散射特性的计算3.6.1散射振幅矩阵与穆勒矩阵的计算散射振幅矩阵S(rsca;rinc)是描述入射电场和散射电场变换关系的重要参数,其定义如下35:E?scaE
44、sca = expj(kr+z)jkr S(rsca;rinc)E?incEinc= expj(kr+z)jkrS1 S4S3 S2E?incEinc; (48)式中, E?inc和Einc表示垂直和平行于散射平面的入射光电场分量, E?sca和Esca是散射光的电场分量,两者满足E?sca = Esca和Esca = E sca 15.根据(48)式,散射振幅矩阵可通过以下方式得到:分别取入射光偏振方向垂直和平行于散射平面,入射光电场强度E?inc = Einc = E0,对应分别可得到两情形下的散射光的电场矢量为(E?sca;?;Esca;?)T和(E?sca;Esca;)T,基于矩阵运算
45、则有:E?sca;? E?sca;Esca;? Esca;= expj(kr+z)jkrS1 S4S3 S2E?inc 00 Einc; (49)S(rsca;rinc)= jkrE0expj(kr+z)E?sca;? E?sca;Esca;? Esca;: (50)通过(50)式便可获得散射振幅矩阵,在此基础上,通过穆勒矩阵F与散射振幅矩阵S的转换关系,便可计算穆勒矩阵,详细请参见文献15, 30.3.6.2吸收截面的推导吸收截面表征的是粒子吸收电磁波的能力,实质是损耗电磁能而将其转化为热能.记电磁波损耗功率功率为Pab,则Pab可采用下式表示:Pab = ReSns S(r)d2r= ReVS(r)d3r= 12ReV(E(r)044207-8万方数据物理学报Acta Phys. Sin. Vol. 66, No. 4 (2017) 044207H (r)d3r; (51)上式中, ns为积分面的外法向矢量, S(r)为电磁场的坡印廷矢量.采用矢量恒等式 (E(r) H (r) = H (r) E(r) E(r) H (r)对体积积分部分进行化简可得(E(r)H (r)= j!(r)E(r)E (r) (r)H(r)H (r) (r)E(r)E (r): (52)代入(51)式,可得电磁波损耗功率Pab的形式为Pab = 12V(r)E(r)E (r)d3r: