基于相移降阶模型周期结构的三维电磁场仿真模拟方法

摘要

本发明公开了一种基于相移降阶模型周期结构的三维电磁场仿真模拟方法，通过在单个相移点下将有限元法获得的原本征方程进行Taylor级数展开，缩减为维数更低的本征方程，这一过程也就是建立降阶模型。一旦维数更低的本征方程即降阶模型建立后，就能非常迅速地计算多个相移点下周期结构的本征方程，从而获得宽频带内周期结构的高频特性。利用本发明提出的基于相移降阶模型周期结构的三维电磁场仿真模拟方法可以在一个相移点的情形下对周期结构进行有限元本征分析，能够快速精确地获得周期结构在整个带宽内的高频特性。

说明书

技术领域

本发明属于三维电磁场数值求解技术领域，涉及一种基于相移降阶模型周期结构的三维 电磁场仿真模拟方法。

背景技术

周期结构在微波管中应用非常广泛，包括波导截面的周期性变化，波导周期加载膜片， 周期填充介质等。在微波管中，普遍采用周期结构作为器件的高频电路，形成电子注与高频 场相互作用进行能量交换以实现微波振荡或放大的场所。周期结构的高频特性：包括色散特 性、阻抗特性与衰减特性，直接影响器件的工作频率、频带宽度、换能效率和输出功率，以 及其他一系列整管性能。高精度地获得周期结构的高频特性有着极其重要的意义。

随着计算机技术的发展，利用有限元法对微波管周期结构进行三维本征分析，从而获得 其高频特性，已经成为微波管周期结构研究和设计中最常用和最有效的方法之一。这类方法 的主要过程是用准周期边界条件将微波管周期结构进行截断得到一个周期长度的计算区域， 然后对该计算区域进行三维网格离散，采用有限元法可以将周期结构的电磁问题转换成一个 大型广义本征值方程，对求解该本征值方程所得到的本征值和本征向量进行一系列的后处理 就能最终获得周期结构的高频特性。

目前，周期结构有限元本征分析算法通常在多个不同的相移点下对周期结构进行本征分 析，得到各个相移点对应的高频特性。然后利用插值或拟合等数值技术得到整个带宽内的高 频特性曲线。对具有宽频带特性的周期结构，这些离散相移点的数量必须足够多，才能保证 插值或拟合的精度。这直接导致了周期结构高频特性的整体仿真效率很低，从而严重影响高 性能的周期结构的设计。

发明内容

本发明的目的是克服现有技术中存在的缺陷，提供一种基于相移降阶模型周期结构的三 维电磁场仿真模拟方法，该方法通过在单个相移点下将有限元法获得的原本征方程进行 Taylor级数展开，缩减为维数更低的本征方程，这一过程也就是建立降阶模型。一旦维数更 低的本征方程即降阶模型建立后，就能非常迅速地计算多个相移点下周期结构的本征方程， 从而获得宽频带内周期结构的高频特性。

其技术方案为：

一种基于相移降阶模型周期结构的三维电磁场仿真模拟方法，包括以下步骤：

A.建立周期结构内的电磁场边值问题，通过有限元法的标准变分原理得到电磁场边值 问题的泛函方程；

B.采用四面体网格剖分求解域，并保证周期边界主面和从面上的网格匹配；

C.选择基函数，将电场强度矢量在所有网格内用基函数展开，并运用里兹方法得到周 期结构有限元本征分析中的广义本征值方程；

D.将步骤C所得到的广义本征值方程中的所有变量用Taylor级数展开；

E.获得步骤D中本征向量的Taylor展开系数构成的子空间K；

F.根据步骤E获得的子空间K，建立周期结构有限元本征分析中的相移降阶模型；

G.利用步骤F中获得的周期结构的相移降阶模型，可以快速计算周期结构在给定相移 范围内的频率以及对应的电场，通过进一步后处理可得周期结构高频特性。

进一步优选，所述步骤A中，首先，根据麦克斯韦方程组与周期结构的边界条件与导 体属性得到周期结构内电磁场的边值问题，如下：

$\left(\begin{array}{cc}▿\times {\mu }_r^{-1}▿\times E-k_0^2{ϵ}_{r}E=0& \mathrm{in\Omega }\\ \hat{n}\times (E_s\times \hat{n})=\hat{n}\times (E_m\times \hat{n})e^{-\mathrm{j\beta L}}& \mathrm{on}{\Gamma }_{\mathrm{PBC}}\\ \hat{n}\times E=0& \mathrm{on}{\Gamma }_{\mathrm{PEC}}\end{array}\right)---\left(1\right)$

公式(1)中第一个式子为频域矢量波动方程，是周期结构有限元仿真中的主方程；其 中，Ω为周期结构的仿真区域空间范围，即为公式(1)的求解域，是矢性偏微分算子符 号，μ_r为求解域Ω中介质的相对磁导率，E为求解域Ω的电场强度矢量，k₀为自由空间波 数，ε_r为求解域Ω中介质的相对介电常数；

公式(1)中第二个式子为准周期边界条件，其中，Γ_PBC表示准周期边界，由主面及从 面组成，其中主面定义为沿着周期结构周期性的方向上仿真区域的初始端面，从面定义为沿 着周期结构周期性的方向上仿真区域的最后端面，主面与从面的距离为一个空间周期，为 边界的外法向单位矢量；E_m和E_s分别表示周期边界主面和从面上的电场；j为虚数单位符 号；β为相位常数；L为周期长度；βL即为一个空间周期对应的相移Φ，准周期边界 条件的物理意义是在周期边界从面上的电磁场和主面上的电磁场，相差一个复数相位系数 e^-jβL；

公式(1)中第三个式子为理想导体的电壁边界条件，其中，Γ_PEC表示电壁边界， 理想导体的电壁边界条件的物理意义是理想导体的切向电场为零；

周期边界Γ_PBC和电壁边界Γ_PEC组成了求解域Ω的外边界；

从周期结构内电磁场的边值问题，即从公式(1)出发，通过有限元法的标准变分原理得 到电磁场边值问题的泛函方程F(E)，即公式(2)：

$F\left(E\right)=\frac{1}{2}{\int }_{\Omega }[{(▿\times E)}^{*}\frac{1}{{\mu }_r}(▿\times E)-k_0^2{ϵ}_r{E}^{*}·E]\mathrm{d\Omega }---\left(2\right)$

公式(2)中，上标*表示对物理量取共轭，dΩ表示三维体积分的微元，使泛函方程即公 式(2)取极小值，并且满足公式(1)中的第二个方程的电场函数E即为周期结构内电磁场边 值问题即公式(1)的解。

进一步优选，所述步骤B中采用四面体网格剖分求解域是有限元方法，剖分后的求解域 被人为分割为多个三维四面体网格，从而将连续的几何结构空间转化为离散的网格空间。

进一步优选，所述步骤C中选择合适的基函数N_m，将公式(2)和公式(1)中第二个方 程中的电场E在所有网格内用基函数N_m展开，即

$E=\underset{m}{\Sigma }{x_m}^t{N}_m^t---\left(3\right)$

公式(3)中，下标m取值从0到M_T，M_T为网格中所有基函数的个数，上标t将电场展 开系数x_m和基函数N_m按照所在区域进行区分，t∈{I，M，S}，I表示基函数N_i在除周期 边界外的计算区域内，M表示基函数N_m在周期边界主面上，S表示基函数N_m在周期边界 从面上；

将所有网格内的电场E用基函数N_m展开后，代入泛函方程，即公式(2)，运用里兹方 法，并利用周期边界条件将展开系数用替代，得到周期结构有限元本征分析中的广 义本征值方程；

最后得到周期结构有限元本征分析中的广义本征值方程为：

$\mathrm{Ax}=\mathrm{\lambda Bx}=k_0^2\mathrm{Bx}---\left(4\right)$

公式(4)中，为广义本征值方程的本征值，k₀＝ω/c为自由空间波数，ω为角频率，c 为自由空间光速；x为广义本征值方程的本征向量，也即公式(3)中插值系数组成的向量，A 和B为大型稀疏M×M维有限元矩阵，M等于所采用基函数N_m的总数，有限元矩阵A和 B的每一项A_ij和B_ij可以分别由公式(5)和公式(6)计算可得：

$A_{\mathrm{ij}}={\int }_{\Omega }{W}_i^{*}(▿\times N_i)·\frac{1}{{\mu }_r}(▿\times N_j)W_j\mathrm{d\Omega }---\left(5\right)$

$B_{\mathrm{ij}}={\int }_{\Omega }{W}_i^{*}{N}_i·{ϵ}_r{N}_j{W}_j\mathrm{d\Omega }---\left(6\right)$

公式(5)和公式(6)中，N_i和N_j为插值基函数，dΩ表示三维体积分的微元，Ω为周期结 构的仿真区域空间范围，∫_Ω(·)dΩ表示在仿真区域空间范围内的体积分，W_i与W_j为周期边界 条件因子，如果插值基函数为从面上的基函数，则对应的周期边界条件因子W＝e^-jβL，否则 W＝1；

从公式(5)和公式(6)可以知道，有限元矩阵A和B是相移Φ＝βL的函数。

进一步优选，所述步骤D中从步骤C得到了周期结构有限元本征分析的广义本征值方 程，

Ax＝λBx    (7)

将公式(7)中的所有变量分别采用Taylor级数展开，得到：

$\left(\begin{array}{c}A=A_0+A_1(u-u_e)+{A_2(u-u_e)}^2+A_3{(u-u_e)}^3+...\\ B=B_0+B_1(u-u_e)+B_2{(u-u_e)}^2+B_3{(u-u_e)}^3+...\\ \lambda ={\lambda }_0+{\lambda }_1(u-u_e)+{\lambda }_2{(u-u_e)}^2+{\lambda }_3{(u-u_e)}^3+...\\ x=x_0+x_1(u-u_e)+x_2{(u-u_e)}^2+x_3{(u-u_e)}^3+...\end{array}\right)---\left(8\right)$

公式(8)中的A₀、A₁、A₂…，B₀、B₁、B₂…，λ₀、λ₁、λ₂…，x₀、x₁、x₂…为采用Taylor 级数展开后的多项式的展开系数，当对相移Φ＝βL进行Taylor级数展开时，u＝Φ＝βL，u_e为 展开相移点；

有限元矩阵A和B中的元素分别由公式(5)和公式(6)给出，是相移u＝βL的函数， 因此有限元矩阵A和B关于相移u＝βL的Taylor级数展开系数A₀、A₁、A₂…，B₀、B₁、B₂… 可以直接求得；

将公式(8)代入至公式(7)中展开后，对比系数即可以得到：

A₀x₀＝λ₀B₀x₀    (9)

(A₀-λ₀B₀)x₁＝λ₁B₀x₀-(A₁-λ₀B₁)x₀    (10)

(A₀-λ₀B₀)x₂＝λ₂B₀x₀-(A₁-λ₀B₁)x₁-(A₂-λ₀B₂)x₀+λ₁(B₀x₁+B₁x₀)    (11)。

进一步优选，所述步骤E中在展开点处求解广义本征方程公式(9)，可以获得x₀以及λ₀， 然后在公式(10)中，两边同时乘以x₀的厄密向量后，公式(10)的左端就变为零，于是通 过计算就可以得到λ₁，再将λ₁代入公式(10)求解关于矩阵(A₀-λ₀B₀)的方程可以得到x₁， 接着针对公式(11)按照相同的方式递归求解，最终可以获得由本征向量展开系数构成的子 空间K＝[x₀，x₁，…，x_N]。

进一步优选，所述步骤F中对K矩阵进行奇异值分解获得正交归一化矩阵Q，将公式 (7)两端同时乘以Q的厄密矩阵Q^H，并利用正交归一化矩阵Q满足的Q^HQ＝I，I为单位 矩阵，可以得到相移降阶模型

$\hat{A}\left(u\right)\hat{x}\left(u\right)=\lambda \left(u\right)\hat{B}\left(u\right)---\left(12\right)$

其中

$\hat{A}\left(u\right)=Q^{H}A\left(u\right)Q$

$\hat{B}\left(u\right)=Q^{H}B\left(u\right)Q---\left(13\right)$

$\hat{x}\left(u\right)=Q^{H}x$

由于矩阵Q为M×N维矩阵，Q^H为N×M矩阵，A(u)和B(u)都为M×M维矩阵，由 公式(13)可知，和为N×N矩阵，N是公式(8)中Taylor级数的阶数，而M 为步骤C中所采用基函数N_m的总数，N远远小于M，广义本征值方程的相移降阶模型即 公式(12)的维数，和原广义本征值方程即公式(7)的维数相比较非常小，利用广义本征 值方程的相移降阶模型就可以在改变相移点u下快速获得相应的本征值以及本征向量了，实 现周期结构的快速扫相移本征分析。

进一步优选，所述步骤G中在指定的相移范围内取一个确定的相位值u＝Φ＝βL，代入 相移降阶模型即公式(12)并求解该广义本征方程，可得本征值λ(u)和本征向量，利用 公式(14)(15)得到周期结构在指定相移下的波数k₀和对应的电场基函数插值系数组成的 本征向量x，

$k_0=\sqrt{\lambda \left(u\right)}---\left(14\right)$

$x=Q\hat{x}\left(u\right)---\left(15\right)$

根据指定相移对应的波数k₀即可获得角频率ω，根据对应的电场基函数插值系数组 成的本征向量x，利用公式(3)得到电场，最后通过后处理可以进一步得到周期结构 的高频特性。

与现有技术相比，本发明的有益效果：

利用本发明提出的基于相移降阶模型周期结构的三维电磁场仿真模拟方法可以在 一个相移点的情形下对周期结构进行有限元本征分析，能够快速精确地获得周期结构在 整个带宽内的高频特性。

附图说明

图1是本发明基于相移降阶模型周期结构的三维电磁场仿真模拟方法的流程图。

具体实施方式

下面结合附图和具体实施例来详细描述本发明的技术方案。

参照图1，一种基于相移降阶模型周期结构的三维电磁场仿真模拟方法，包括以下步骤：

A.建立周期结构内的电磁场边值问题，通过有限元法的标准变分原理得到电磁场边值 问题的泛函方程。

首先，根据麦克斯韦方程组与周期结构的边界条件与导体属性得到周期结构内电磁 场的边值问题，如下：

$\left(\begin{array}{cc}▿\times {\mu }_r^{-1}▿\times E-k_0^2{ϵ}_{r}E=0& \mathrm{in\Omega }\\ \hat{n}\times (E_s\times \hat{n})=\hat{n}\times (E_m\times \hat{n})e^{-\mathrm{j\beta L}}& \mathrm{on}{\Gamma }_{\mathrm{PBC}}\\ \hat{n}\times E=0& \mathrm{on}{\Gamma }_{\mathrm{PEC}}\end{array}\right)---\left(1\right)$

公式(1)中第一个式子为频域矢量波动方程，是周期结构有限元仿真中的主方程；其 中，Ω为周期结构的仿真区域空间范围，即为公式(1)的求解域。是矢性偏微分算子符 号，μ_r为求解域Ω中介质的相对磁导率，E为求解域Ω的电场强度矢量，k₀为自由空间波 数，ε_r为求解域Ω中介质的相对介电常数。

公式(1)中第二个式子为准周期边界条件，其中，Γ_PBC表示准周期边界，由主面及从 面组成，其中主面定义为沿着周期结构周期性的方向上仿真区域的初始端面，从面定义为沿 着周期结构周期性的方向上仿真区域的最后端面。主面与从面的距离为一个空间周期。为 边界的外法向单位矢量；E_m和E_s分别表示周期边界主面和从面上的电场；j为虚数单位符 号；β为相位常数；L为周期长度；βL即为一个空间周期对应的相移Φ。准周期边界 条件的物理意义是在周期边界从面上的电磁场和主面上的电磁场，相差一个复数相位系数 e^-jβL；

公式(1)中第三个式子为理想导体的电壁边界条件，其中，Γ_PEC表示电壁边界。 理想导体的电壁边界条件的物理意义是理想导体的切向电场为零。

周期边界Γ_PBC和电壁边界Γ_PEC组成了求解域Ω的外边界；

从周期结构内电磁场的边值问题，即从公式(1)出发，通过有限元法的标准变分原理得 到电磁场边值问题的泛函方程F(E)，即公式(2)：

$F\left(E\right)=\frac{1}{2}{\int }_{\Omega }[{(▿\times E)}^{*}\frac{1}{{\mu }_r}(▿\times E)-k_0^2{ϵ}_r{E}^{*}·E]\mathrm{d\Omega }---\left(2\right)$

公式(2)中，上标*表示对物理量取共轭，dΩ表示三维体积分的微元。使泛函方程即公 式(2)取极小值，并且满足公式(1)中的第二个方程的电场函数E即为周期结构内电磁场边 值问题即公式(1)的解。

B.采用四面体网格剖分求解域，并保证周期边界主面和从面上的网格匹配。

采用四面体网格剖分求解域是有限元方法中的一种公知过程，因此本步骤不再详细描 述。剖分后的求解域被人为分割为多个三维四面体网格，从而将连续的几何结构空间转化为 离散的网格空间。

C.选择基函数，将电场强度矢量在所有网格内用基函数展开，并运用里兹方法得到周 期结构有限元本征分析中的广义本征值方程。

选择合适的基函数N_m，将公式(2)和公式(1)中第二个方程中的电场E在所有网格内 用基函数N_m展开，即

$E=\underset{m}{\Sigma }{x_m}^t{N}_m^t---\left(3\right)$

公式(3)中，下标m取值从0到M_T，M_T为网格中所有基函数的个数。上标t将电场展 开系数x_m和基函数N_m按照所在区域进行区分。t∈{I，M，S}，I表示基函数N_i在除周期 边界外的计算区域内，M表示基函数N_m在周期边界主面上，S表示基函数N_m在周期边界 从面上。

将所有网格内的电场E用基函数N_m展开后，代入泛函方程，即公式(2)，运用里兹方 法，并利用周期边界条件将展开系数用替代，得到周期结构有限元本征分析中的广 义本征值方程。该过程为电磁场有限元分析中众所周知的过程，这里不再赘述。

最后得到周期结构有限元本征分析中的广义本征值方程为：

$\mathrm{Ax}=\mathrm{\lambda Bx}=k_0^2\mathrm{Bx}---\left(4\right)$

公式(4)中，为广义本征值方程的本征值，k₀＝ω/c为自由空间波数，ω为角频率，c 为自由空间光速；x为广义本征值方程的本征向量，也即公式(3)中插值系数组成的向量。A 和B为大型稀疏M×M维有限元矩阵，M等于所采用基函数N_m的总数。有限元矩阵A和 B的每一项A_ij和B_ij可以分别由公式(5)和公式(6)计算可得：

$A_{\mathrm{ij}}={\int }_{\Omega }{W}_i^{*}(▿\times N_i)·\frac{1}{{\mu }_r}(▿\times N_j)W_j\mathrm{d\Omega }---\left(5\right)$

$B_{\mathrm{ij}}={\int }_{\Omega }{W}_i^{*}{N}_i·{ϵ}_r{N}_j{W}_j\mathrm{d\Omega }---\left(6\right)$

公式(5)和公式(6)中，N_i和N_j为插值基函数，dΩ表示三维体积分的微元，Ω为周期结 构的仿真区域空间范围，∫_Ω(·)dΩ表示在仿真区域空间范围内的体积分。W_i与W_j为周期边界 条件因子。如果插值基函数为从面上的基函数，则对应的周期边界条件因子W＝e^-jβL，否则 W＝1。其他符号的含义与公式(1)同。

从公式(5)和公式(6)可以知道，有限元矩阵A和B是相移Φ＝βL的函数。

D.将步骤C所得到的广义本征值方程中的所有变量用Taylor级数展开。

步骤C得到了周期结构有限元本征分析的广义本征值方程，

Ax＝λBx    (7)

将公式(7)中的所有变量分别采用Taylor级数展开，得到：

$\left(\begin{array}{c}A=A_0+A_1(u-u_e)+{A_2(u-u_e)}^2+A_3{(u-u_e)}^3+...\\ B=B_0+B_1(u-u_e)+B_2{(u-u_e)}^2+B_3{(u-u_e)}^3+...\\ \lambda ={\lambda }_0+{\lambda }_1(u-u_e)+{\lambda }_2{(u-u_e)}^2+{\lambda }_3{(u-u_e)}^3+...\\ x=x_0+x_1(u-u_e)+x_2{(u-u_e)}^2+x_3{(u-u_e)}^3+...\end{array}\right)---\left(8\right)$

公式(8)中的A₀、A₁、A₂…，B₀、B₁、B₂…，λ₀、λ₁、λ₂…，x₀、x₁、x₂…为采用Taylor 级数展开后的多项式的展开系数。当对相移Φ＝βL进行Taylor级数展开时，u＝Φ＝βL，u_e为 展开相移点。

有限元矩阵A和B中的元素分别由公式(5)和公式(6)给出，是相移u＝βL的函数， 因此有限元矩阵A和B关于相移u＝βL的Taylor级数展开系数A₀、A₁、A₂…，B₀、B₁、B₂… 可以直接求得。

将公式(8)代入至公式(7)中展开后，对比系数即可以得到：

A₀x₀＝λ₀B₀x₀    (9)

(A₀-λ₀B₀)x₁＝λ₁B₀x₀-(A₁-λ₀B₁)x₀    (10)

(A₀-λ₀B₀)x₂＝λ₂B₀x₀-(A₁-λ₀B₁)x₁-(A₂-λ₀B₂)x₀+λ₁(B₀x₁+B₁x₀)      (11)

.

.

.

E.获得步骤D中本征向量的Taylor展开系数构成的子空间K。

在展开点处求解广义本征方程公式(9)，可以获得x₀以及λ₀。然后在公式(10)中， 两边同时乘以x₀的厄密向量后，公式(10)的左端就变为零，于是通过计算就可以得到λ₁。 再将λ₁代入公式(10)求解关于矩阵(A₀-λ₀B₀)的方程可以得到x₁。接着针对公式(11)按 照相同的方式递归求解，最终可以获得由本征向量展开系数构成的子空间K＝[x₀，x₁，…，x_N]。

F.根据步骤E获得的子空间K，建立周期结构有限元本征分析中的相移降阶模型。

对K矩阵进行奇异值分解获得正交归一化矩阵Q。将公式(7)两端同时乘以Q的厄密 矩阵Q^H，并利用正交归一化矩阵Q满足的Q^HQ＝I，I为单位矩阵，可以得到相移降阶模型

$\hat{A}\left(u\right)\hat{x}\left(u\right)=\lambda \left(u\right)\hat{B}\left(u\right)---\left(12\right)$

其中

$\hat{A}\left(u\right)=Q^{H}A\left(u\right)Q$

$\hat{B}\left(u\right)=Q^{H}B\left(u\right)Q---\left(13\right)$

$\hat{x}\left(u\right)=Q^{H}x$

由于矩阵Q为M×N维矩阵，Q^H为N×M矩阵，A(u)和B(u)都为M×M维矩阵，由 公式(13)可知，和为N×N矩阵。N是公式(8)中Taylor级数的阶数，而M 为步骤C中所采用基函数N_m的总数。N远远小于M。因此广义本征值方程的相移降阶模 型即公式(12)的维数，和原广义本征值方程即公式(7)的维数相比较非常小。利用广义 本征值方程的相移降阶模型就可以在改变相移点u下快速获得相应的本征值以及本征向量 了，实现周期结构的快速扫相移本征分析。

G.利用步骤G中获得的周期结构的相移降阶模型，可以快速计算周期结构在给定相移 范围内的波数k₀以及对应的电场基函数插值系数组成的本征向量x，通过进一步后处理可得 周期结构的高频特性。

在指定的相移范围内取一个确定的相位值u＝Φ＝βL，代入相移降阶模型即公式(12)并 求解该广义本征方程，可得本征值λ(u)和本征向量，利用公式(14)(15)得到周期结 构在指定相移下的波数k₀和对应的电场基函数插值系数组成的本征向量x。

$k_0=\sqrt{\lambda \left(u\right)}---\left(14\right)$

$x=Q\hat{x}\left(u\right)---\left(15\right)$

根据指定相移对应的波数k₀即可获得角频率ω。根据对应的电场基函数插值系数组 成的本征向量x，利用公式(3)得到电场。最后通过后处理可以进一步得到周期结构 的高频特性。该过程为本领域的公知过程，因此不再详细描述。

以上所述，仅为本发明较佳的具体实施方式，本发明的保护范围不限于此，任何熟悉本 技术领域的技术人员在本发明披露的技术范围内，可显而易见地得到的技术方案的简单变化 或等效替换均落入本发明的保护范围内。