CN114781202A - 一种包含非线性色散媒质电磁问题的仿真方法 - Google Patents

Abstract

本发明公开了一种包含非线性色散媒质电磁问题的仿真方法，包括：根据器件媒质的几何尺寸，建立媒质的几何模型，采用六面体单元对该模型进行剖分，获取器件的几何信息并且设置激励源的位置信息；根据并行参数的核数对已经剖分完成的几何模型进行区域划分；采用时域不连续伽辽金有限元方法建立微分方程组，对媒质的几何模型中任意色散模型和非线性色散媒质模型求解。本发明以时域不连续伽辽金有限元为基础，结合局部时间步进和并行技术，可以实现复杂媒质的建模，以及对非线性色散媒质的电磁特性分析；由于时域不连续伽辽金有限元法有高精度及建模灵活的特点，可以直接得到逆矩阵，由于局部时间步进和并行技术的加入，能极大降低计算时间。

Description

一种包含非线性色散媒质电磁问题的仿真方法

技术领域

本发明属于非线性问题的电磁特性仿真技术，特别是一种包含非线性色散媒质电磁问题的仿真方法。

背景技术

由于激光器的出现，非线性媒质效应逐步映入人们的视线当中，其极化强度不再与电场成线性关系，而是成高次幂的形式，随着入射波幅值的增大，其非线性效应也随之增大。非线性媒质在光学中的应用前景广泛，比如光学开关、四波混频(FWM)、放大器、调制器，在光纤通信中可以形成形状和速度均不变的光孤子波，有助于实现无损通信。对于非线性媒质问题，一般可以计算Maxwell方程组，得到非线性媒质的电磁特性。可以使用辅助微分方法和拟线性求解非线性媒质问题，从而避免了牛顿迭代，其中拟线性求解要求前一时刻与后一时刻的时间间隔不能变化太快，否则容易降低求解精度，所以时间迭代步长会选取得很小，随着非线性的增强，步长选取会越来越小，但是每个时刻求解一步即可，且矩阵为快对角，为后面并行提供可行性，而牛顿迭代的时间步长选取较为宽松，但是求解容易不收敛。

非线性媒质是一个重要的研究领域，在2015年，S.Yan和金建铭利用牛顿迭代方法结合FDTD分析了非线性铁磁问题，并且引入了非均匀时间步长技术。同年，Zhu B等人将非线性电磁参数视为分段线性参数，在求解速度上有了提高，但是分段线性参数构造不具有通用性。2019年，伊朗科技大学的Moradi M等人采用基于交替隐式法的FDTD分析非线性色散媒质的三次谐波现象，其算法时间步进选取突破了传统FDTD稳定性条件限制，但当时间步稍大时，会导致计算结果不精确。

发明内容

本发明的目的在于提供一种包含非线性色散媒质电磁问题的仿真方法，提高了计算准确度，且能够大大缩短计算时间，降低计算成本。

实现本发明目的的技术方案为：一种包含非线性色散媒质电磁问题的仿真方法，包括步骤：

第一步，根据器件媒质的几何尺寸建立媒质的几何模型，采用六面体单元对几何模型进行剖分，得到器件媒质的几何信息并且设置激励源的位置信息以及施加激励源；

第二步，调用并行参数，根据调用的核数对已经剖分完成的几何模型进行区域划分；

第三步，采用时域不连续伽辽金有限元方法建立微分方程组，即将选定的基函数对未知的电场、磁场变量近似展开，然后分别代入偏微分方程组中，通过伽辽金测试，得到矩阵方程组，根据色散媒质特点，选取色散媒质参数，根据任意色散模型进行描述，用辅助微分方程法对其求解，根据非线性媒质特点，选取求解方案；

第四步，引入局部时间步进技术，将整个求解区域根据线性与非线性区域划分，将上一时刻大时间步下的解插值到小时间步中，求解出对应小时间步下的非线性方程解，再将小时间步求得的解带入大时间步下的方程，求解出全部电磁场的解。

本发明与现有技术相比，其显著优点为：

(1)本发明提及的方法是一种包含非线性媒质的高效不连续伽辽金有限元方法，在以往对非线性媒质问题进行模拟仿真时，其思路是用牛顿迭代直接对非线性项求解，影响全局时间步进，本发明方法用线性拟合方法，且引入局部时间步进技术和并行技术，对其进行加速；

(2)采用时域不连续伽辽金有限元(DGTD)作为仿真所用的方法，它使用曲六面体剖分，建模灵活，剖分方便，随着多项式阶数的提高，计算误差将呈指数下降；同时DGTD生成的质量矩阵为块对角矩阵，可以直接求逆计算，能够大大缩短计算时间，降低计算成本。

附图说明

图1为引入局部时间步进技术求解过程示意图。

图2为非线性色散结构示意图。

图3为观察点处归一化功率谱密度和参考文献对比图。

具体实施方式

本发明提出的方法位基于时域不连续伽辽金有限元方法的对包含非线性色散媒质的高效数值仿真方法，该方法将局部时间步进技术和并行技术，引入到时域不连续伽辽金有限元中，非线性媒质包含克尔非线性模型和拉曼非线性模型，色散模型包含德拜、德鲁德、洛伦兹，对非线性色散模型仿真提供支持。

本发明提出了一种包含非线性色散媒质的高效不连续伽辽金有限元方法。包括步骤如下：

第一步，根据媒质的几何尺寸，建立媒质的几何模型，采用六面体单元对该模型进行剖分，得到目标的几何信息并且设置激励源的位置信息；

第二步，调用并行参数，根据调用的核数对已经剖分完成的几何模型进行区域划分；

第三步，采用时域不连续伽辽金有限元方法建立微分方程组。首先，通过伽辽金对麦克斯韦方程组两端进行测试，然后用选定的基函数对方程中未知的电场、磁场变量近似展开，得到矩阵方程组。最后，媒质是由色散媒质和非线性媒质共同组成(媒质包含任意色散媒质和非线性媒质)。对于含色散媒质，运用任意色散模型(GDM)进行表征，用辅助微分方程法对其求解，对于含非线性色散媒质，用辅助微分方程法对其求解；

第四步，引入局部时间步进技术，将整个求解区域根据线性与非线性区域划分，将上一时刻大时间下的解插值到小时间中，求解出对应小时间步下的非线性方程解，再将小时间步求得的解带入大时间步方程，求解出全部电磁场的解；

所述第三步中，给出媒质的麦克斯韦方程组在频域的公式为：

其中E表示电场，H表示磁场，P_Kerr和P_Raman为非线性极化矢量。

所述第三步中，含色散媒质时，给出任意色散模型的表达式：

式中，ε₀和ε_∞表示真空和频率无穷大时的相对介电常数，任意色散模型的参数包括a_0,p，a_1,p，b_0,p，b_1,p和极点数M，表1给出了典型的任意色散模型参数。

表1典型色散模型的GDM参数

通过辅助微分方程化简式(1.2)中线性色散部分，令

令

将式(1.3)进行傅里叶逆变换：

对方程(1.4)进行伽辽金测试并对未知量进行展开，消去方程两边共同项得：

其中J′_p为J_p的展开系数，e为E的展开系数，引入变量γⁿ＝∫_tJ_p，并对式(1.4)进行中心差分得：

其中

所述第三步中，含非线性媒质的求解具体为：

其中已典型三阶项构成的非线性模型为例：克尔非线性(Kerr)模型和拉曼非线性(Raman)模型。

考虑到三阶非线性的极化公式如下：

式中

是三阶磁化率张量，对公式(1.7)使用Born-Oppenheimer近似，公式简化为：

公式中g(t)是非线性因果响应函数。

对于克尔非线性模型，公式(1.8)中非线性因果响应函数表达为：

g(t)＝αδ(t) (1.9)

其中δ(t)是狄拉克三角函数(Dirac delta function)，用来模拟瞬时克尔非共振跃迁，α是一个范围在0到1的常数，用来表示克尔非线性和拉曼非线性相对强度。

对于拉曼非线性模型，公式(1.8)中非线性因果项响应函数表达式为：

g(t)＝(1-α)g_Raman(t) (1.11)

其中g_Raman为拉曼非线性响应，拉曼非线性极化表达为：

引入辅助变量

其中

将公式(1.13)进行傅里叶变换，得到表达式为：

式中F{·}代表傅里叶变换，频谱响应函数表达式如下：

其中

将式(1.16)代入式(1.15)可以得到

将式(1.18)左右两端乘上

因为不连续伽辽金时域方法需要时域信息，所以对其进行傅里叶逆变换，得到ADE方程的表示形式为：

对式(1.19)以时间步n进行中心差分，得到S的离散表达式：

将式(1.13)代入式(1.12)式可以得到：

P_Raman(t)＝ε₀S(t)E(t) (1.21)

所述第四步中，局部时间步进技术的计算过程，以求解(t+Δt_l)时刻为例：

将求解区域用四面体和六面体进行混合离散，利用稳定性条件确定电磁场区域的时间步长Δt_l，其为大时间步单元。

由于非线性媒质中非线性项采用上一时刻近似，需要取较小的时间步进来保证解的精确性，且往往会远小于线性区域的时间步长，本专利中取Δt_d＝Δt_l/M(M为正整数)，其为小时间步单元。

通过不同的时间步进进行未知场值的更新，求解基于任意高阶导数的DGTD构建非线性色散媒质结构的系统方程：

其中矩阵

其中v为待计算的电场和磁场向量，Q是离散形式中左端方程的系数，K为右端电场和磁场向量前面的系数矩阵，N为基函数，Z为本体的本征阻抗，Y为本征导纳，S和e_n分别是计算范围的外表面和外法向，ε为介电常数，μ为磁导率，不属于本单元的体通过“+”来标记。

小时间单元区域的时间迭代公式可以表达为：

其中v为待计算的电场和磁场向量，Q是离散形式中左端方程的系数，K为右端电场和磁场向量前面的系数矩阵，m为ADER的阶数。式(1.22)中第二项需要用到大单元的相关插值信息，由于相邻大单元中不存在中间时刻t+(n-1)Δt_l/M，因此无法直接获得未知场值v(t+(n-1)Δt_l/M)。那么，式(1.23)中的表达将不能继续求解，需要将该项表示成泰勒级数展开的原始形式

通过相邻的大时间步单元计算出来的变量信息插值得：

时间步为Δt_l的电磁场区域场值的更新格式同理可得：

给出局部时间步技术中，两个区域大小时间步的更新示意图如图1所示。

实施例

以图2所示结构器件为例，对本发明的具体步骤作进一步详细描述。

媒质中包含线性洛伦兹极化和非线性极化，其中非线性中包含克尔非线性和拉曼非线性。具体参数如下，线性Lorentz极化参数：ε_s＝5.25，ε_∞＝2.25，ω_p＝4×10¹⁴sec^-1，δ_p＝2×10⁹sec^-1。克尔和拉曼极化参数：

α＝0.7，τ₁＝12.2fs，τ₂＝32fs。其中线性洛伦兹色散部分使用GDM进行描述。非线性色散模型如图2所示，模型x方向和y方向使用PBC边界条件，z方向两端使用ABC边界条件，其具体操作步骤如下：

第一步，根据图2所示包含结构的几何尺寸，利用Ansys软件对其进行建模，采用六面体单元对结构进行剖分，即得到结构的未知量个数以及节点的坐标信息，对于结构的边界信息，可依据模型的尺寸，将x和y方向设为周期边界条件，z方向设为吸收边界条件来处理，并且设置激励源的位置信息，模型采用两个载波频率的调制高斯平面波，其表达式为：

其中f₁＝300THz，f₂＝400THz，t₀＝86.66fs，τ＝28.89fs

第二步，调用并行参数，根据调用的核数对已经剖分完成的几何模型进行区域划分；

第三步，根据建立的非线性色散方程组(媒质的麦克斯韦方程组)，采用DGTD建立微分方程组，首先对方程两端进行伽辽金测试，接着将选定的六面体基函数对未知的电场、磁场变量近似展开，得到DGTD半离散方程。

第四步，首先计算小时间步进区域中非线性的场值信息，并将大时间步进单元信息插值到小时间步进单元，接着更新大时间步进单元场值信息。

第五步，对计算出的电场值做归一化傅里叶变换。

根据本发明所述方法对如图1所示构器件进行仿真，观察点处归一化功率谱密度与参考文献对比见图3，表2计算效率比，图3和表2的对比结果表明，本发明仿真效果好，能够大大缩短计算时间，降低计算成本。

表2计算效率比较

Claims (10)

1.一种包含非线性色散媒质电磁问题的仿真方法，其特征在于，包括步骤：

根据器件媒质的几何尺寸，建立媒质的几何模型，采用六面体单元对该模型进行剖分，获取器件的几何信息并且设置激励源的位置信息；

根据并行参数的核数对已经剖分完成的几何模型进行区域划分；

采用时域不连续伽辽金有限元方法建立微分方程组，对媒质的几何模型中任意色散模型和非线性色散媒质模型求解；

引入局部时间步进方法，求解出全部电磁场的解。

2.根据权利要求1所述的一种包含非线性色散媒质电磁问题的仿真方法，其特征在于，所述采用时域不连续伽辽金有限元方法建立微分方程组，对媒质的几何模型中任意色散模型和非线性色散媒质模型求解具体包括：

通过伽辽金对媒质的麦克斯韦方程组进行测试，用选定的基函数对未知的电场、磁场变量近似展开，获取矩阵方程组；并通过辅助微分方程法对任意色散模型求解，对于非线性色散媒质模型，根据非线性媒质特点对其求解。

3.根据权利要求2所述的一种包含非线性色散媒质电磁问题的仿真方法，其特征在于，所述任意色散模型为：

其中，ε₀和ε_∞表示真空和频率无穷大时的相对介电常数，a_0,p，a_1,p，b_0,p，b_1,p为任意色散模型的参数，M为极点数。

4.根据权利要求3所述的一种包含非线性色散媒质电磁问题的仿真方法，其特征在于，所述媒质的麦克斯韦方程组为:

其中E表示电场，H表示磁场，P_Kerr和P_Raman为非线性极化矢量。

5.根据权利要求1所述的一种包含非线性色散媒质电磁问题的仿真方法，其特征在于，所述通过辅助微分方程法对任意色散模型求解具体包括：

通过辅助微分方程化简任意色散模型中线性色散部分，令

令

将式(1.3)进行傅里叶逆变换得到：

对方程(1.4)进行伽辽金测试并对未知量进行展开，消去方程两边共同项得：

其中J′_p为J_p的展开系数，e为E的展开系数，引入变量γⁿ＝∫_tJ_p，并对式(1.4)进行中心差分得：

其中

6.根据权利要求2所述的一种包含非线性色散媒质电磁问题的仿真方法，其特征在于，所述非线性媒质包括克尔非线性模型和拉曼非线性模型。

7.根据权利要求6所述的一种包含非线性色散媒质电磁问题的仿真方法，其特征在于，所述对于非线性色散媒质模型，根据非线性媒质特点对其求解具体包括：

所述非线性色散媒质的三阶非线性的极化公式为：

其中

是三阶磁化率张量，对公式(1.7)使用Born-Oppenheimer近似简化为：

其中g(t)是非线性因果响应函数；

对于克尔非线性模型，公式(1.8)中非线性因果响应函数为：

g(t)＝αδ(t) (1.9)

其中δ(t)是狄拉克三角函数，α是一个范围在0到1的常数，克尔非线性极化为：

对于拉曼非线性模型，公式(1.8)中非线性因果项响应函数表达式为：

g(t)＝(1-α)g_Raman(t) (1.11)

其中g_Raman为拉曼非线性响应，拉曼非线性极化为：

引入辅助变量

其中

将公式(1.13)进行傅里叶变换，得到表达式为：

式中F{·}代表傅里叶变换，频谱响应函数表达式如下：

其中

将式(1.16)代入式(1.15)可以得到

将式(1.18)左右两端乘上

对其进行傅里叶逆变换，得到ADE方程的表示形式为：

对式(1.19)以时间步n进行中心差分，得到S的离散表达式：

将式(1.13)代入式(1.12)式可以得到：

P_Raman(t)＝ε₀S(t)E(t) (1.21)

8.根据权利要求7所述的一种包含非线性色散媒质电磁问题的仿真方法，其特征在于，所述入局部时间步进方法，求解出全部电磁场的解具体为：将整个求解区域根据线性与非线性区域划分，将上一时刻大时间步单元下的解插值到小时间步单元中，求解出对应小时间步单元下的非线性方程解，再将小时间步单元下求得的解带入大时间步单元下的方程，求解出全部电磁场的解。

9.根据权利要求8所述的一种包含非线性色散媒质电磁问题的仿真方法，其特征在于，将整个求解区域根据线性与非线性区域划分，将上一时刻大时间步单元下的解插值到小时间步单元中，求解出对应小时间步单元下的非线性方程解，再将小时间步单元下求得的解带入大时间步单元下的方程，求解出全部电磁场的解具体为：

将求解区域用六面体进行离散，利用稳定性条件确定电磁场区域的时间步长Δt_l，其为大时间步单元，另小时间步单元为Δt_d＝Δt_l/M，M为正整数；

小时间单元区域的时间迭代公式为：

其中：

其中v为待计算的电场和磁场向量，m为ADER的阶数，；N为基函数，Z为本体的本征阻抗，Y为本征导纳，S和e_n分别是计算范围的外表面和外法向，ε为介电常数，μ为磁导率，不属于本单元的体通过“+”来标记；式(1.22)表示成泰勒级数展开的原始形式：

通过相邻的大时间步单元计算出来的变量信息插值得：

同理时间步为Δt_l的大时间步区域场值的更新格式为：

10.根据权利要求1～9任一所述的一种包含非线性色散媒质电磁问题的仿真方法，其特征在于，所述媒质通过Ansys软件建模，并行参数的核数为8核。

Images