一种二维结构散射场的分析方法及系统

1.本发明属于纳米光子学领域，更具体地，涉及一种二维结构散射场的分析方法及系统。


背景技术：

2.有限元方法是计算电磁学中一种常用的数值计算手段，其网格划分的思想使得该方法能够求解任意复杂的结构，因此被广泛应用在光学器件的设计中。从1940年代诞生至今，随着计算机计算能力的大幅度提升，有限元也一步一步发展成熟且形成了比较通用的商业软件，例如comsol，ansys等。有限元计算散射问题需要使用合适的边界条件。其中完美匹配层可以较好的模拟无限大空间，但需要使用额外的空间，增大计算量，降低计算速度。散射边界通常情况下也模拟无限大空间，但在某些情况下并不能得到正确的结果。
3.从海洋学到天文学，从量子力学到高能物理，包括生物医学研究等方向，都广泛的存在光的散射现象；再具体到实验数据的分析、场的探测、遥感观测都需要对光的散射现象进行细致的分析。在对光散射现象研究中，多极子理论是极为重要的一种手段，可以对产生光散射的源、近场、远场进行定量分析。
4.2019年，陈伟锦等人在physical review letters上发表论文“singularities and poincare indexes of electromagnetic multipoles”，研究了电磁多极子的拓扑性质。文中需要得到散射场的多极子系数从而进一步分析拓扑特性，采用多极子系数对产生光散射的源、近场、远场进行定量分析。然而他们求解多极子系数的方法比较繁琐：先使用comsol计算出散射场，再使用多极子展开法求出多极子系数。此外，针对求解复杂结构散射场的多极子系数，市面上并没有直接且简洁的计算方法。因此，目前需要一种能够直接计算出散射场多极子系数，且具有较好性能的方法。


技术实现要素：

5.针对现有技术的缺陷，本发明的目的在于提供一种二维结构散射场的分析方法及系统，旨在解决现有求解二维结构散射场多极子系数的方法求解过程复杂，没有直接且简洁的计算方法，不便于用于对散射场进行定量分析的问题。
6.为实现上述目的，第一方面，本发明提供了一种二维结构散射场的分析方法，包括如下步骤：
7.确定二维结构散射场的标量波动方程；所述标量波动方程反映二维结构所在环境的总电磁场、二维结构的相对磁导率、二维结构的相对介电常数之间的关系；所述二维结构为二维各向同性介质；
8.将所述标量波动方程与第一测试函数相乘，并在预设区域内进行二重积分，并结合格林公式得到所述标量波动方程的弱形式；所述预设区域为圆形；利用磁场边界连续性条件，将预设区域内部磁场和预设区域外部磁场联系起来，以将预设区域之外的电场引入到所述弱形式，得到引入外电场的弱形式；
9.将电场边界连续性条件与第二测试函数相乘，得到电场在预设区域边界处的电场连续性方程；
10.将所述预设区域划分为多个三角网格，并将各个三角网格内的电场用形函数展开，将所述预设区域之外的外电场用电多极子展开；
11.采用伽辽金法，将第一测试函数取遍所有的形函数，第二测试函数取遍所有的电多极子，分别更新所述引入外电场的弱形式和电场连续性方程；
12.将以形函数展开得到的三角网格内电场和以电多极子展开得到的外电场代入更新后的弱形式和电场连续性方程，得到全空间内电场矩阵形式对应的线性方程；
13.对所述线性方程求解得到散射场的电多极子系数，以便利用所述电多极子系数对散射场进行定量分析。
14.在一个可选的示例中，所述标量波动方程为：
[0015][0016]
其中，ez为总电磁场，μr为二维结构的相对磁导率，∈r为二维结构的相对介电常数，k0为真空波数，
[0017]
在一个可选的示例中，所述标量波动方程的弱形式的具体求解过程如下：
[0018]
选取任意连续的第一测试函数w与标量波动方程相乘并在计算区域d内二重积分：
[0019][0020]
其中ds是二重积分的面积元素；
[0021]
使用格林公式对所述二重积分进行进一步处理得到：
[0022][0023]
其中，为区域d的边界，逆时针方向，n为边界单位外法向量，u取为a取为w，b取为ez，dl为曲线积分的积分元素；
[0024]
将总电磁场写成由入射电场和散射电场之和的形式并带入上述处理后的公式，其中为散射电场，为入射电场，可得到标量波动方程的弱形式：
[0025][0026]
在一个可选的示例中，将有限元求解区域内的散射电场用表示，外部区域散射
电场用表示；利用磁场边界连续性条件，标量波动方程的弱形式转化得到引入外电场的弱形式：
[0027][0028]
在一个可选的示例中，对电场边界条件选取任意连续的测试函数wk相乘并在边界上做曲线积分，得到电场连续性方程：
[0029][0030]
其中，为边界单位外法向量，
[0031]
将区域d划分为多个小三角网格，每个小三角网格内形函数为上标(e)表示小网格的形函数，下标z表示标量场，ja取1，2，3；小三角网格内散射电场使用形函数展开开为展开系数；
[0032]
将区域d之外的外部区域散射电场使用电多极子展开，部区域散射电场使用电多极子展开，为电多极子，其中，k为波数，zm(kr)为汉克尔函数，m为阶数，r为柱坐标下半径，φ为柱坐标系下夹角，jb取-k到k整数，k为选取的任意整数，用于决定多极子数目；e
zmnjb
为电多极子对应的系数。
[0033]
在一个可选的示例中，使用伽辽金法，w取遍所有小三角网格的wk取遍所有使用的电多极子此处的下标ia、ib分别与ja、jb取值范围相同，但代表测试函数，则将引入外电场的弱形式和电场连续性方程分别更新为：
[0034][0035][0036]
再将以小网格形函数和多极子展开的形式带入上式，可得每个小三角
内的矩阵形式：
[0037][0038][0039]
其中：每个小网格形函数系数，e
mn
＝[e
zmnjb
]电多极子系数；矩阵az的(ia，ja)项向量bz的ia项矩阵c的(ia，jb)项矩阵d的(ib，ja)项矩阵d的(ib，ja)项矩阵e的(ib，ja)项
[0040]
求解上述矩阵形式对应的线性方程式得到散射场的电多极子系数e
mn
。
[0041]
第二方面，本发明提供了一种二维结构散射场的分析系统，包括：
[0042]
波动方程确定单元，用于确定二维结构散射场的标量波动方程；所述标量波动方程反映二维结构所在环境的总电磁场、二维结构的相对磁导率、二维结构的相对介电常数之间的关系；所述二维结构为二维各向同性介质；
[0043]
弱形式确定单元，用于将所述标量波动方程与第一测试函数相乘，并在预设区域内进行二重积分，并结合格林公式得到所述标量波动方程的弱形式；以及利用磁场边界连续性条件，将预设区域内部磁场和预设区域外部磁场联系起来，以将预设区域之外的电场引入到所述弱形式，得到引入外电场的弱形式；所述预设区域为圆形；
[0044]
电场连续性方程确定单元，用于将电场边界连续性条件与第二测试函数相乘，得到电场在预设区域边界处的电场连续性方程；
[0045]
线性方程确定单元，用于将所述预设区域划分为多个三角网格，并将各个三角网格内的电场用形函数展开，将所述预设区域之外的外电场用电多极子展开；采用伽辽金法，将第一测试函数取遍所有的形函数，第二测试函数取遍所有的电多极子，分别更新所述引入外电场的弱形式和电场连续性方程；以及将以形函数展开得到的三角网格内电场和以电多极子展开得到的外电场代入更新后的弱形式和电场连续性方程，得到全空间内电场矩阵形式对应的线性方程；
[0046]
散射场分析单元，用于对所述线性方程求解得到散射场的电多极子系数，以便利用所述电多极子系数对散射场进行定量分析。
[0047]
在一个可选的示例中，波动方程确定单元确定的所述标量波动方程为：
[0048]
[0049]
其中，ez为总电磁场，μr为二维结构的相对磁导率，∈r为二维结构的相对介电常数，k0为真空波数，
[0050]
在一个可选的示例中，所述弱形式确定单元确定标量波动方程的弱形式的具体求解过程如下：
[0051]
选取任意连续的第一测试函数w与标量波动方程相乘并在计算区域d内二重积分：
[0052][0053]
其中ds是二重积分的面积元素；
[0054]
使用格林公式对所述二重积分进行进一步处理得到：
[0055][0056]
其中，为区域d的边界，逆时针方向，n为边界单位外法向量，u取为a取为w，b取为ez，dl为曲线积分的积分元素；
[0057]
将总电磁场写成由入射电场和散射电场之和的形式并带入上述处理后的公式，其中为散射电场，为入射电场，可得到标量波动方程的弱形式：
[0058][0059]
将有限元求解区域内的散射电场用表示，外部区域散射电场用表示；利用磁场边界连续性条件，标量波动方程的弱形式转化得到引入外电场的弱形式：
[0060][0061]
在一个可选的示例中，所述电场连续性方程确定单元对电场边界条件选取任意连续的测试函数wk相乘并在边界上做曲线积分，得到电场连续性方程：
[0062][0063]
其中，为边界单位外法向量，
[0064]
所述线性方程确定单元，将区域d划分为多个小三角网格，每个小三角网格内形函
数为上标(e)表示小网格的形函数，下标z表示标量场，ja取1，2，3；小三角网格内散射电场使用形函数展开射电场使用形函数展开为展开系数；将区域d之外的外部区域散射电场使用电多极子展开，部区域散射电场使用电多极子展开，为电多极子，为电多极子，其中，k为波数，zm(kr)为汉克尔函数，m为阶数，r为柱坐标下半径，φ为柱坐标系下夹角，jb取-k到k整数，k为选取的任意整数，用于决定多极子数目；e
zmnjb
为电多极子对应的系数；使用伽辽金法，w取遍所有小三角网格的wk取遍所有使用的电多极子此处的下标ia、ib分别与ja、jb取值范围相同，但代表测试函数，则将引入外电场的弱形式和电场连续性方程分别更新为：
[0065][0066][0067]
再将以小网格形函数和多极子展开的形式带入上式，可得每个小三角内的矩阵形式：
[0068][0069][0070]
其中：每个小网格形函数系数，e
mn
＝[e
zmnjb
]电多极子系数；矩阵az的(ia，ja)项向量bz的ia项矩阵c的(ia，jb)项矩阵d的(ib，ja)项矩阵d的(ib，ja)项
矩阵e的(ib，ja)项
[0071]
所述散射场分析单元，求解上述矩阵形式对应的线性方程式得到散射场的电多极子系数e
mn
。
[0072]
总体而言，通过本发明所构思的以上技术方案与现有技术相比，具有以下有益效果：
[0073]
本发明提供一种二维结构散射场的分析方法及系统，利用完备正交的多极子展开外部区域散射场，相较传统计算多极子系数的方法，能够简单便捷地直接获取多极子系数。此外与传统散射边界条件相比，该分析方法计算结果更为准确；与完美匹配层相比，该分析方法计算区域更小，计算效率更高。
附图说明
[0074]
图1是本发明实施例提供的二维结构散射场分析方法流程图。
[0075]
图2是本发明实施例提供的算法的流程图。
[0076]
图3是本发明实施例提供的计算实例的结构示意图。
[0077]
图4(a)是本发明实施例计算出来的散射场。
[0078]
图4(b)是comsol使用散射边界条件计算出来的散射场。
[0079]
图4(c)是comsol使用完美匹配层计算出来的散射场。
[0080]
图5是本发明实施例计算出来的多极子系数与由comsol散射场经多极子展开得到的多极子系数对比图。
[0081]
图6是本发明实施例提供的二维结构散射场分析系统架构图。
具体实施方式
[0082]
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。
[0083]
本发明涉及纳米光子学领域。更具体地，涉及一种有限元法与多极子相结合的混合数值计算方法，来计算各向同性介质的二维结构散射问题中标量散射场的方法，并且能够直接得到散射场的多极子系数。多极子理论可使用求得的多极子系数对产生光散射的源、近场、远场进行定量分析。本发明所要解决的问题是设计出一种计算各向同性介质的二维结构标量散射问题且能得到散射场多极子系数的算法。
[0084]
本发明的方案是用一个圆包裹目标二维结构，圆内部使用有限元法求解，而圆外部使用多极子展开，内部和外部通过连续性边界条件连接。本发明首先通过选取任意连续函数与标量散射问题的波动方程相乘并在求解区域做二重积分，再结合格林公式得到有限元的弱形式。然后利用磁场边界连续性条件将内部电场与外部电场联系起来，同时消除方程中的磁场未知量。将求解区域划分为一个个三角网格，在每个网格内弱形式可以写成刚度矩阵和载荷向量的形式。再通过网格局部编码与全局编码的关系，组装出整个求解区域的刚度矩阵和载荷向量。求解刚度矩阵和载荷向量构成的线性方程组即可得到散射场及其多极子系数。
[0085]
为了验证本发明计算方法的正确性，本发明在数值仿真软件comsol分别使用完美匹配层和散射边界计算同样的散射问题，再在matlab中使用comsol的接口计算出多极子系数。结果显示本发明计算出的散射场和多极子系数与comsol完全一致，从而验证了本发明方法的正确性。
[0086]
目前市面上没有能够直接得到散射场多极子系数的手段，需要通过多步繁琐计算，因此本发明的算法填补了该缺口，为纳米光子学的发展做出一定的贡献。此外本发明的方法比有限元中使用的完美匹配层更快，比使用的散射边界条件更准。
[0087]
图1是本发明实施例提供的二维结构散射场分析方法流程图。如图1所示，包括如下步骤：
[0088]
s101，确定二维结构散射场的标量波动方程；所述标量波动方程反映二维结构所在环境的总电磁场、二维结构的相对磁导率、二维结构的相对介电常数之间的关系；所述二维结构为二维各向同性介质；
[0089]
s102，将所述标量波动方程与第一测试函数相乘，并在预设区域内进行二重积分，并结合格林公式得到所述标量波动方程的弱形式；所述预设区域为圆形；
[0090]
s103，利用磁场边界连续性条件，将预设区域内部磁场和预设区域外部磁场联系起来，以将预设区域之外的电场引入到所述弱形式，得到引入外电场的弱形式；
[0091]
s104，将电场边界连续性条件与第二测试函数相乘，得到电场在预设区域边界处的电场连续性方程；
[0092]
s105，将所述预设区域划分为多个三角网格，并将各个三角网格内的电场用形函数展开，将所述预设区域之外的外电场用电多极子展开；
[0093]
s106，采用伽辽金法，将第一测试函数取遍所有的形函数，第二测试函数取遍所有的电多极子，分别更新所述引入外电场的弱形式和电场连续性方程；
[0094]
s107，将以形函数展开得到的三角网格内电场和以电多极子展开得到的外电场代入更新后的弱形式和电场连续性方程，得到全空间内电场矩阵形式对应的线性方程；
[0095]
s108，对所述线性方程求解得到散射场的电多极子系数，以便利用所述电多极子系数对散射场进行定量分析。
[0096]
具体地，散射问题中，二维各向同性介质结构对应的标量波动方程为：
[0097][0098]
其中ez为总场，μr为二维结构的相对磁导率，∈r为二维结构的相对介电常数，k0为真空波数，
[0099]
选取任意连续的测试函数w与公式(1)相乘并在计算区域d内二重积分：
[0100][0101]
使用第一标量格林公式
为面积分区域d的边界，逆时针方向，n为边界单位外法向量。u取为a取为w，b取为ez，得到：
[0102][0103]
将总场写成由入射场和散射场之和的形式并带入上式，其中为散射场，为入射场，可得到：
[0104][0105]
有限元求解区域内的散射场用表示，外部区域散射场用表示。
[0106]
利用磁场边界连续性条件及磁场与电场关系式因而有其中hs为求解区域内的磁场，为外部区域的磁场，n为边界单位外法线向量，h代表磁场，e代表电场，c为常数。则式(2)可化为：
[0107][0108]
同样，对电场边界条件选取任意连续的测试函数wk相乘并在边界上做曲线积分：
[0109][0110]
其中为边界单位外法向量，
[0111]
待求区域d划分为一个个小三角网格，每个小三角网格内形函数为上标(e)表示小网格的形函数，下标z表示标量场，ja取1，2，3。小三角网格内散射场可以使用形函数展开形函数展开为展开系数。
[0112]
外部区域散射场使用电多极子展开，外部区域散射场使用电多极子展开，为电多极子，其中，k为波数，zm(kr)为汉克尔函数，m为阶数，r为柱坐标下半径，φ为柱坐标系下夹角，jb取-k到k整数，k为选取的任意整
数，用于决定多极子数目。e
zmnjb
为电多极子对应的系数。
[0113]
使用伽辽金法galerkin法，w取遍所有小三角网格的wk取遍所有使用的电多极子此处的ia、ib与上述的ja、jb取值范围相同，但代表测试函数，则式(3)和式(4)变为：
[0114][0115][0116]
再将以小网格形函数和多极子展开的形式带入上式，可得每个小三角内的矩阵形式：
[0117][0118][0119]
其中：每个小网格形函数系数，e
mn
＝[e
zmnjb
]电多极子系数；矩阵az的(ia，ja)项向量bz的ia项矩阵c的(ia，jb)项矩阵d的(ib，ja)项矩阵d的(ib，ja)项矩阵e的(ib，ja)项
[0120]
求解线性方程式(5)即可得到散射场及散射场的多极子系数e
mn
。
[0121]
上述步骤可总结为流程图图2，图中矩量法弱形式代指上述步骤中电场连续性边界条件得到的曲线积分方程。
[0122]
以图3所示二维结构为实例，其中有限元求解区域为空气包裹着圆形散射体，散射体的半径为λ为真空波长，分别使用本发明计算、使用comsol散射边界条件计算、使用
comsol完美匹配层计算，其中comsol完美匹配层计算结果可认为是真实结果。本发明计算得到的散射场如图4(a)所示，comsol散射边界条件计算结果如图4(b)所示，comsol使用完美匹配层计算结果如图4(c)所示。对比可知，本发明计算结果比comsol散射边界条件计算结果更为准确，计算区域相比comsol完美匹配层更小。图4(a)-图4(c)的坐标为空间坐标的横轴和纵轴。
[0123]
随后使用comsol完美匹配层计算结果计算多极子系数与本发明计算的多极子系数对比，如图5所示，在一定误差范围内，本发明也能够正确计算出多极子系数。
[0124]
图6为本发明实施例提供的二维结构散射场的分析系统架构图。如图6所示，包括：
[0125]
波动方程确定单元610，用于确定二维结构散射场的标量波动方程；所述标量波动方程反映二维结构所在环境的总电磁场、二维结构的相对磁导率、二维结构的相对介电常数之间的关系；所述二维结构为二维各向同性介质；
[0126]
弱形式确定单元620，用于将所述标量波动方程与第一测试函数相乘，并在预设区域内进行二重积分，并结合格林公式得到所述标量波动方程的弱形式；以及利用磁场边界连续性条件，将预设区域内部磁场和预设区域外部磁场联系起来，以将预设区域之外的电场引入到所述弱形式，得到引入外电场的弱形式；所述预设区域为圆形；
[0127]
电场连续性方程确定单元630，用于将电场边界连续性条件与第二测试函数相乘，得到电场在预设区域边界处的电场连续性方程；
[0128]
线性方程确定单元640，用于将所述预设区域划分为多个三角网格，并将各个三角网格内的电场用形函数展开，将所述预设区域之外的外电场用电多极子展开；采用伽辽金法，将第一测试函数取遍所有的形函数，第二测试函数取遍所有的电多极子，分别更新所述引入外电场的弱形式和电场连续性方程；以及将以形函数展开得到的三角网格内电场和以电多极子展开得到的外电场代入更新后的弱形式和电场连续性方程，得到全空间内电场矩阵形式对应的线性方程；
[0129]
散射场分析单元650，用于对所述线性方程求解得到散射场的电多极子系数，以便利用所述电多极子系数对散射场进行定量分析。
[0130]
可以理解的是，图6中各个单元的详细功能实现可参见前述方法实施例的介绍，在此不做赘述。
[0131]
本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。