一种基于有限元边界积分法的纳米结构局域表面等离激元共振分析仿真方法

1.本发明涉及纳米结构局域表面等离激元共振分析仿真方法，特别是一种基于有限元边界积分法纳米结构局域表面等离激元共振分析仿真方法。


背景技术：

2.纳米粒子与入射电磁波相互作用产生的独特性质一直以来都是科研的热点问题，随着国内外科研人员不断地耕耘研究，等离激元光子学成为一门逐渐兴起的学科，尤其是当入射光照射在金属纳米粒子上时产生的局域表面等离激元共振(lspr)现象受到越来越多科研工作者的关注。当入射光作用在金属纳米粒子表面时，使得金属纳米粒子表面的自由电子产生集体震荡，并产生沿纳米粒子表面传播的电磁波，这使得纳米粒子对特定波长的入射光吸收与散射显著加强，自由电子发生共振现象，这种共振现象被称为局域表面等离激元共振。局域表面等离激元共振现象已经应用在了生活中的方方面面，例如在生物医疗中利用特定波长下金属纳米粒子发生局域表面等离激元共振产生释放大量的热能来杀死病变的细胞组织，在太阳能电池板中加入金属纳米粒子，利用特定波长下发生局域表面等离激元共振对入射光吸收增强来提高光能转换为电能的效率。
3.对于纳米结构结构局域表面等离激元共振的分析计算方法中，近年来数值分析方法由于适用范围广精确度高等优点变得越来越受欢迎。目前使用较为广泛的数值分析方法主要有时域有限差分法、有限元法和矩量法。其中时域有限差分法和有限元法可以分析内部比较复杂的纳米结构，但是由于不能天然满足索末菲辐射边界条件，因此在分析散射和辐射等开域问题时需要人为地添加截断边界条件，增加了离散区域。矩量法由于格林函数的存在天然满足索末菲辐射边界条件，适合分析开域问题，而且由于没有数值色散误差精确度比时域有限差分和有限元法也要高，但是矩量法由于剖分的是纳米结构的表面不适合分析内部比较复杂的纳米结构。有限元边界积分法保留了有限元和矩量法两种数值方法的优点摒弃了两者的缺点，引入一个虚构边界，边界内部区域的场利用有限元方法描述，边界外部的场利用边界积分描述，然后利用边界两侧场的连续性建立有限元边界积分方程。使用有限元边界积分方法既能分析内部复杂的纳米结构，在边界外部也不需要人为地添加截断边界条件，减少了离散区域。
4.因此发展一种基于有限元边界积分法纳米结构局域表面等离激元共振分析仿真方法具有重要的应用价值。


技术实现要素：

5.本发明的目的在于提供一种基于有限元边界积分法纳米结构局域表面等离激元共振分析仿真方法。该方法精度高，适用范围广，通过精确的电磁建模，可对内部复杂的纳米结构实现局域表面等离激元共振的精确分析。
6.实现本发明的技术解决方案是：一种基于有限元边界积分法纳米结构局域表面等
离激元共振分析仿真方法，包括如下步骤：
7.第1步，基于有限元边界积分法基本理论引入一个包围纳米结构的虚构边界，边界内部的场用有限元方法表达，边界外部的场用边界积分表达。
8.第2步，边界内部区域利用四面体单元进行剖分，并使用whitney基函数对有限元方程进行离散，边界面利用三角形单元进行剖分，并使用rwg基函数对边界积分方程进行离散，通过边界两侧场的连续性获得有限元边界积分方程组。
9.第3步，求解有限元边界积分方程组获得区域内场的系数和边界面上等效电流和磁流的系数，进而求得区域内部的电场和磁场以及边界面上的等效电流和磁流。
10.利用内观法，外观法或者混合算法求解有限元边界积分方程组获得区域内场的系数和边界面上等效电磁流的系数后，通过高斯积分定理对边界内部区域的每个小四面体单元进行循环并结合whitney基函数和内部区域场的系数可以求得内部区域的场，对边界面的小三角形单元进行循环结合rwg基函数和等效电磁流的系数即可求出边界面上的等效电流和等效磁流。
11.第4步，通过边界面上的等效电流和磁流计算边界外部的散射电场和散射磁场并计算不同入射角度下的散射电场强度。
12.根据散射电磁场的计算公式对格林函数取梯度并结合等效电流和等效磁流，利用高斯积分定理对边界面上的每一个小三角形进行循环并累加计算出散射电场和散射磁场。设置不同的入射角度，通过传播矢量的计算公式求得方向坐标，进而代入入射场改变右端向量根据散射电场的计算公式计算出当前入射角度下的散射电场强度。
13.第5步，计算入射场，散射场和总场的坡应廷矢量求得纳米结构的消光截面和局域表面等离激元共振波长。
14.将入射电场和入射磁场的共轭进行叉乘得到入射场的坡印廷矢量，将散射电场和散射磁场的共轭进行叉乘得到散射场的坡印廷矢量，将总电场和总磁场的共轭进行叉乘得到总场的坡印廷矢量，然后根据消光截面的计算公式利用高斯积分定理结合入射场，散射场和总场的坡印廷矢量求得消光截面，分析消光截面共振峰的位置进而获得局域表面等离激元共振波长。
15.第6步，计算纳米结构垂直于入射波方向的截面积进而求得消光系数。
16.利用公式法，积分法或网格法计算不同形态纳米结构垂直于入射波方向的截面积，然后利用消光截面与计算出来的截面积的比值求得消光系数。
17.本发明与现有技术相比，其显著的效果是：
18.本发明所述方法在纳米结构局域表面等离激元共振分析中，保留了有限元和矩量法两种数值方法的优点摒弃了两者的缺点，和其它现有技术相比具有以下两个明显的优势。
19.(1)精度高，未知量少。纳米结构外部的场使用边界积分方法进行分析，和时域有限差分法、有限元等数值方法，物理光学法、几何光学法等高频近似方法相比精确度高，而且由于矩量法天然满足索末菲辐射边界条件不需要人为地添加索末菲辐射边界条件，因此和时域有限差分，有限元等基于微分方程的方法相比离散区域小，未知量少。
20.(2)适用范围广。纳米结构内部的场使用有限元方法进行分析，和矩量法，解析法等方法相比适用范围更广，不仅可以分析内部均匀的纳米结构也可以分析内部复杂的纳米
结构。
附图说明
21.图1是具体实施方式中基于有限元边界积分法纳米结构局域表面等离激元共振分析仿真方法流程图；
22.图2是具体实施方式中有限元边界积分法纳米结构示意图；
23.图3是具体实施方式的测试例1中银纳米球壳模型xoz截面图；
24.图4是具体实施方式的测试例1中本发明所述方法离散银纳米球壳模型示意图；
25.图5是具体实施方式的测试例1中有限元方法离散银纳米球壳模型示意图；
26.图6是具体实施方式的测试例1中银纳米球壳归一化消光截面与comsol仿真结果对比图；
27.图7是具体实施方式的测试例2中银纳米球壳阵列模型xoz截面图；
28.图8是具体实施方式的测试例2中不同间距银纳米球壳阵列归一化消光截面结果图。
具体实施方式
29.为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。此外，下面所描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。
30.图1是本发明基于有限元边界积分法纳米结构局域表面等离激元共振分析仿真方法的整体流程图。如图1所示，该流程主要包括以下的处理步骤：
31.第1步，基于有限元边界积分法基本理论引入一个包围纳米结构的虚构边界，边界内部的场用有限元方法表达，边界外部的场用边界积分表达。
32.图2给出了使用有限元边界积分法时纳米结构模型示意图，其中，s为包围纳米结构的虚构边界，为了使积分区域最小设置s为纳米结构外边界，ε
r
为区域v内纳米结构的介电常数，μ
r
为区域v内纳米结构的磁导率。
33.区域v内的场利用有限元方法表示，在无源的区域v内，电场满足矢量波动方程
[0034][0035]
其中，e(r)表示区域v内的待求电场，ε和μ分别表示自由空间的介电常数和磁导率。令并利用用矢量恒等式和散度定理，得到等效变分问题为
[0036][0037]
通过上式可以看出，区域v内可以包含多种纳米结构，并且介电常数可以为实数也可以为复数。
[0038]
边界s外的场利用边界积分表达，可以建立电场积分方程(efie)。
[0039][0040]
也可以建立磁场积分方程(mfie)。
[0041][0042]
上式中，为边界面上的等效电流，为边界面上的等效磁流，e
i
和分别为入射电场和入射磁场。
[0043]
l和k算子分别为：
[0044][0045]
其中，为格林函数，y(r)＝x(r)
×
n，ω是相对于观察点的立体角，r
′
代表源点位置，r代表场点位置。
[0046]
第2步，边界内部区域利用四面体单元进行剖分，并使用whitney基函数对有限元方程进行离散，边界面利用三角形单元进行剖分，并使用rwg基函数对边界积分方程进行离散，通过边界两侧场的连续性获得有限元边界积分方程组。
[0047]
边界面s的内部区域v使用四面体单元进行剖分，每个小四面体单元上的场都可以利用 whitney基函数表示为：
[0048][0049]
其中，为待求系数，为whitney基函数，e
e
为每个小四面体单元上的场。
[0050]
将(6)代入有限元泛函变分方程(2)式中进行离散得到有限元矩阵方程为：
[0051][0052]
边界面s上使用三角形单元进行剖分，每个小三角形单元上的等效电流和等效磁流可以表示为：
[0053][0054]
其中，就是rwg基函数，和e
i
为待求系数，选取不同的试函数分别对电场积分方程和磁场积分方程进行离散可以得到tenenh方程，如下式所示。关于tenenh方程的内容不是本发明的重点，在这里不再赘述。
[0055]
[0056]
根据边界面两侧场的连续性，可以得到有限元边界积分方程的矩阵形式为：
[0057][0058]
其中，p和q是矩量法生成的阻抗矩阵，它们是稠密矩阵，k和b是有限元生成的阻抗矩阵，它们是稀疏矩阵，并且k是对称矩阵，b是反对称矩阵，b是矩量法生成的右端向量， e
i
,e
s
,是待求系数。
[0059]
第3步，求解有限元边界积分方程组获得区域内场的系数和边界面上等效电流和磁流的系数，进而求得区域内部的电场和磁场以及边界面上的等效电流和磁流。
[0060]
通过内观法，外观法或混合算法求解有限元边界积分方程组获得纳米结构内部区域内场的系数和纳米结构表面等效电流和等效磁流的系数，利用高斯积分定理对边界内部区域的每个小四面体单元进行循环并结合whitney基函数和内部区域场的系数可以求得内部区域的场，对纳米结构表面的小三角形单元进行循环，结合rwg基函数和每个小三角形单元上每条边上等效电流和磁流的系数能够得到当前小三角形单元上的等效电流和磁流，然后进行累加即可得到纳米结构表面的等效电流和等效磁流。
[0061]
第4步，通过边界面上的等效电流和磁流计算边界外部的散射电场和散射磁场并计算不同入射角度下的散射电场强度。
[0062]
散射电场和散射磁场的计算公式：
[0063][0064][0065]
其中，g是格林函数。
[0066]
在求得纳米结构表面的等效电流和等效磁流后，根据散射电磁场的计算公式(11)(12) 式对格林函数取梯度并结合等效电流和等效磁流，利用高斯积分定理对纳米结构表面上的每一个小三角形进行循环并累加可以计算出散射电场和散射磁场。
[0067]
入射波是一个均匀平面波，它的传播矢量为：
[0068][0069]
其中，θ
i
是与z轴正向的夹角，φ
i
是与x轴正向的夹角，通过设置不同的θ
i
和φ
i
可以得到不同的入射方向和方向坐标，进而改变有限元边界积分矩阵方程中的右端向量，计算出不同入射角度下的散射电场强度。
[0070]
第5步，计算入射场，散射场和总场的坡应廷矢量求得纳米结构的消光截面和局域表面等离激元共振波长。
[0071]
入射波是一个均匀平面波，均匀平面波的坡印廷矢量是一个常数，入射光强度为：
[0072][0073]
其中，e0为入射电场强度，η为自由空间波阻抗，散射场的坡应廷矢量为散射电场
和散射磁场共轭的叉乘，散射光强度的计算公式为：
[0074][0075]
其中e
s
表示纳米结构外部散射电场，h
s
表示纳米结构外部散射磁场，利用高斯积分定理对散射场的坡印廷矢量在纳米结构的外表面上做面积分即可得到散射光强度。
[0076]
纳米结构外部的总场为入射场和散射场之和，总场的坡印廷矢量为总电场和总磁场共轭的叉乘，吸收光强度的计算公式为：
[0077][0078]
其中，e表示总电场，h表示总磁场，利用高斯积分定理对总场的坡印廷矢量在纳米结构外表面上做面积分即可得到吸收光强度。
[0079]
消光截面也称为强度衰减截面，表示当电磁波入射到物体上被物体吸收和散射的总的辐射通量，计算公式为：
[0080][0081]
求得消光截面后，通过分析共振峰的强弱和位置可以得到局域表面等离激元共振波长。
[0082]
第6步，计算纳米结构垂直于入射波方向的截面积进而求得消光系数。
[0083]
根据不同形态的纳米结构选择不同的计算垂直于入射波方向的截面积的方法，对于球这种规则形状，可以采用公式法进行计算，对于螺旋体这种形状比较复杂的纳米结构，需要使用积分法或者网格法进行计算。
[0084]
求得纳米结构垂直于入射波方向的截面积后，结合消光截面即可求出消光系数，消光系数的计算公式为：
[0085][0086]
其中，σ
e
表示消光截面，s表示纳米结构垂直于入射波方向的截面积，c
e
为消光系数。
[0087]
测试例1
[0088]
采用本发明所述方法对银纳米球壳进行局域表面等离激元共振分析，银纳米球壳模型的 xoz截面图如图3所示，r0和r1分别为银纳米球壳的内核半径和外壳半径 (r0＝20nm,r1＝25nm)，外层是一个厚度为s的银壳(s＝5nm)，银的介电常数等参数取自 ag(silver)(johnson and christy 0.188
‑
1.937μm)模型，入射光的波长范围400
‑
700nm，内核介质的介电常数为1.44。
[0089]
图4给出了本发明所述方法的模型离散示意图，银纳米球壳表面使用三角形单元进行剖分，内部使用四面体单元进行剖分，剖分密度均为0.025λ。图5给出了使用传统有限元方法计算时的模型离散示意图，银纳米球壳模型、空气和pml均使用四面体单元进行剖分，银纳米球壳的剖分密度为0.025λ，空气和pml的剖分密度为0.05λ，λ表示入射光的波长。表1 给出了两种不同的方法模型离散未知量的对比，可以看到采用本发明所述方法的离散未知量明显少于传统有限元方法。
[0090]
表1
[0091]
数值方法四面体单元数棱边数节点数传统有限元方法345855236710104本发明方法14322386323
[0092]
为了验证本发明所述方法的正确性，图6给出了本发明所述方法计算的归一化消光截面与comsol软件仿真结果的对比图，可以发现本发明所述方法的计算结果与comsol仿真结果吻合较好，验证了本发明所述方法的正确性。并且从图中可以发现当入射光波长等于500nm 时，银纳米球壳发生了局域表面等离激元共振，消光截面出现了大幅度的增大，表明此时银纳米壳对入射光的散射和吸收大幅度增大。银纳米球壳垂直于入射波方向的截面积s＝1962.5nm2，共振波长处消光截面的峰值为σ
emax
＝423.65nm2,在共振波长处消光系数为 s/σ
emax
＝0.216。
[0093]
测试例2
[0094]
采用本发明所述方法对不同间距银纳米球壳阵列进行局部表面等离激元共振分析，图7 给出了银纳米球壳阵列模型的xoz截面图，如图所示r0和r1分别为银纳米球壳的内核半径和外壳半径(r0＝20nm,r1＝25nm)，s为的银壳的厚度(s＝5nm)，d为两个银纳米球壳之间的间距。银的介电常数等参数取自ag(silver)(johnson and christy 0.188
‑
1.937μm)模型，入射光的波长范围400
‑
800nm，内核介质的介电常数为1.44。图8给出了在不同间距d下银纳米球壳阵列归一化消光截面计算结果，从图中可以看出随着银纳米球壳间距的减小，局域表面等离激元共振波长从530nm红移到590nm，并且消光截面共振峰的峰值也在不断增大，说明本发明所述方法能够分析计算内部比较复杂的纳米结构阵列的局域表面等离激元共振现象。