CN104636553A - 微波铁氧体元器件的时域谱元仿真方法 - Google Patents

Abstract

本发明公开了一种微波铁氧体元器件的时域谱元仿真方法，该方法使用曲六面体单元对包含铁氧体材料的整个电磁空间进行离散，将具有张量形式的、包含阻尼因子的铁氧体磁导率引入电场时域亥姆赫兹方程，并且以单轴电各向异性完美匹配层作为吸收边界条件，经过伽辽金变换后，在时间的离散上采用中心差分格式，得到时域电场迭代式；利用时域电场迭代式对目标铁氧体元器件进行两次时域仿真，根据微波网络散射参量的定义，确定目标铁氧体元件在工作频段内各个端口的插入损耗、回波损耗和隔离度。本发明具有计算精度高、计算速度快、适用范围广的优点。

Description

微波铁氧体元器件的时域谱元仿真方法

技术领域

本发明属于电磁仿真技术领域，特别是一种微波铁氧体元器件的时域电磁仿真方法。

背景技术

一般情况下，由各向同性材料所组成的无源微波网络，其各个端口间总是满足互易特性的。但是，如果其中使用了各向异性的材料，就能够得到具有非互易特性的无源微波网络。铁氧体就是一种具有磁各向异性的亚铁磁材料。自1949年发现电磁波在铁氧体中发生的铁磁共振现象以来，人们已经对电磁波在铁氧体材料中传播的非互易传输机理展开了深入的研究。研究表明，铁氧体材料的磁各向异性实际上是由外加的恒定磁场所诱导的，该恒定磁场作用于铁氧体中的磁偶极子，产生合成的磁偶极矩，从而使得磁偶极子在该磁场力的控制下作进动。与进动方向一致的极化波与磁偶极矩发生强烈的互作用，而相反方向的极化波与磁偶极矩发生较弱的互作用。对于给定极化方向的电磁波，其通过加载外部恒定偏置磁场的铁氧体材料时，在不同的方向上会产生不同的传播特性。利用铁氧体的这种特性可以设计出各种非互易微波元器件，例如环形器、隔离器、移相器、开关调制器等。

随着理论研究的不断深入，以及现代微波器件内部结构的复杂程度越来越高，使用数值仿真方法来研究和设计微波铁氧体元器件变得越来越现实。作为一种无源的非线性微波器件，对铁氧体元件的仿真既可以使用频域仿真算法，也可以使用时域仿真算法。有限元法是一种典型的能够仿真复杂媒质的频域算法，大型商用软件HFSS中使用频域有限元法实现了对铁氧体媒质的有效仿真。但使用频域有限元法仿真铁氧体元件时存在两个问题。第一个问题是铁氧体材料的介电常数通常很高，需要使用极细的网格单元进行离散，进而形成一个巨大的线性方程组，导致计算规模非常庞大；第二个问题是为了获得全波信息，在频域有限元中需要按频率逐点仿真，即在每个频率点都要计算一组大型的线性方程组，最终导致仿真效率极低。

时域仿真是研究铁氧体元器件的另一种有效手段。通过时域仿真，不仅可以观察微波信号在铁氧体器件中的瞬态变化过程，而且通过快速傅立叶变换还可以获得近乎连续的频率响应曲线，实现真正的全波仿真。针对铁氧体材料，长久以来，研究人员一直关注于时域有限差分法(FDTD)和时域有限元法(FETD)。

FDTD是一种成熟的数值算法，该算法中，电、磁场在每一时间步均为显示迭代，因此FDTD算法的最大优点在于计算速度快。同时，FDTD算法的内存消耗很低，只与未知量的个数相关。但FDTD算法受限于只能采用矩形网格离散，对于结构较为复杂的器件模型，其拟合精度较差，导致最终的计算精度降低。如果使用极细的网格进行结构拟合，不仅会使未知量成倍增长，增加迭代求解的时间，而且网格离散尺寸过小，也会导致色散误差的增加。

与FDTD相反，时域有限元法(FETD)凭借其在建模时对复杂电磁目标的精确拟合特性，正受到越来越多多关注。但是，使用有限元法都要涉及到对一个大型线性方程组的求解，对于电大尺寸电磁问题，其计算量将非常庞大。当有限元法应用到时域时，存在一个无法回避的问题，即在算法的每个时间步都要计算一次这样的方程组，其总的耗时非常庞大，这也正是时域有限元法难以得到推广的主要原因。

发明内容

本发明的目的在于提供一种微波铁氧体元器件的时域谱元仿真方法。

实现本发明目的的技术方案为：一种微波铁氧体元器件的时域谱元仿真方法，包括以下步骤：

步骤1、根据目标铁氧体元器件的结构尺寸，建立实际坐标系下的电磁空间模型；采用曲六面体单元对该电磁空间模型进行离散，依据每个离散单元的几何参数，建立实际坐标系中各个曲六面体单元与参量坐标系下标准立方体单元的几何映射关系；所述参量坐标系是以标准立方体单元的重心为原点，以标准立方体单元的棱边为坐标轴方向的直角坐标系；

步骤2、选定时域谱元法中基函数的阶数，确定参量坐标系下标准立方体单元中高斯积分点的几何参数，根据上一步的几何映射关系，得到各个曲六面体单元中相应的实际高斯积分点的几何参数，建立参量坐标系与实际坐标系中基函数的映射关系；

步骤3、将具有张量形式的、包含阻尼因子的铁氧体磁导率引入电场时域亥姆赫兹方程，并且以单轴电各向异性完美匹配层作为吸收边界条件，经过伽辽金变换后，在时间的离散上采用中心差分格式，得到时域电场迭代式；

步骤4、利用时域电场迭代式对目标铁氧体元器件进行两次时域仿真；第一次仿真时，将铁氧体材质退化为普通材质，使元件内部在材质上保持连续，记录观察点上的时域仿真数据，当时域仿真数据稳定后，通过FFT变换，获得入射波信号的频率响应；第二次仿真时，恢复铁氧体材质，记录观察点上的时域仿真数据，当时域仿真数据稳定后，通过FFT变换，获得前向和后向传输信号的频率响应；

步骤5、利用上一步中两次仿真获得的频率响应，根据微波网络散射参量的定义，确定目标铁氧体元件在工作频段内各个端口的插入损耗、回波损耗和隔离度。

本发明与现有技术相比，其显著优点为：(1)离散模型的拟合度高：本发明采用曲六面体单元对电磁计算空间进行离散，能够很好地拟合各种复杂的空间形状，保证了离散后模型的精确性；(2)内存需求低：本发明使用高阶基函数，因此空间离散尺寸由0.07倍的最小波长，增加到0.3倍的最小波长，导致空间离散后产生的未知量减少，从而使内存需求大大降低；(3)实现简单：本发明中用于描绘铁氧体元器件内部电场分布的通式不仅适用于电、磁各向同性的普通媒质，也适用于具有磁各向异性的铁氧体材料，还适用于具有单轴电各向异性的PML吸收边界，从而避免了在不同区域使用不同的方程的不便，使得程序的实现过程大大简化；(4)求解速度快：本发明所形成的方程组的系数矩阵为块对角阵，其逆矩阵很容易求得，因此时域的电场迭代格式表现为显式，无需求解方程，因此计算速度极快；另外，现有技术中采用电场旋度方程和磁场旋度方程作为控制方程，在每个时间步需要经过两次计算，本发明采用电场亥姆赫兹方程作为控制方程，每个时间步只需经过一次计算，效率提高了一倍；(5)计算结果可信度高：本发明采用了经过严格论证的包含阻尼因子的铁氧体张量磁导率，并通过合理的代数推导，使其与时域谱元法相融合，形成了严密的时域电场迭代格式，因此，该计算结果可信度高，可以为微波铁氧体元器件的电磁特性分析与设计提供重要的数值参考。

下面结合附图对本发明作进一步详细描述。

附图说明

图1为本发明的微波铁氧体元器件的时域谱元(SETD)仿真方法的流程图。

图2为本发明的实施例中Y形铁氧体波导环行器的结构示意图。

图3为图2所述Y形铁氧体波导环行器的电磁计算空间网格离散图。

图4为本发明的时域谱元法中几何映射关系示意图。

图5为三个方向的基函数均为2阶时，GLL积分点在单位正立方体中的分布图。

图6为包含铁氧体材料的第二次仿真中三个端口观察点处信号的时域波形图。

图7为该Y形铁氧体波导环形器的回波损耗S11的曲线图。

图8为该Y形铁氧体波导环形器的隔离度S21的曲线图。

图9为该Y形铁氧体波导环形器的插入损耗S31的曲线图。

具体实施方式

结合图1，一种微波铁氧体元器件的时域谱元仿真方法，包括以下步骤：

步骤1、根据目标铁氧体元器件的结构尺寸，建立实际坐标系下的电磁空间模型；采用曲六面体单元对该电磁空间模型进行离散，依据每个离散单元的几何参数，建立实际坐标系中各个曲六面体单元与参量坐标系下标准立方体单元的几何映射关系；所述参量坐标系是以标准立方体单元的重心为原点，以标准立方体单元的棱边为坐标轴方向的直角坐标系；具体为：

步骤1-1、根据微波铁氧体元器件的结构尺寸，用计算机辅助设计工具建立实际坐标系下的电磁空间模型；所述实际坐标系为三维直角坐标系，反映了目标铁氧体元器件在真实空间中的位置，实际坐标系的原点和三坐标轴可任意设置，目标铁氧体元器件的外加恒定磁场必须与三坐标轴中的某一坐标轴相平行；

步骤1-2、采用20个节点的曲六面体单元对目标模型进行全局离散，得到所有单元的几何参数；所述20个节点为所在曲六面体单元的8个顶点和12条棱边的中点，所述几何参数包括单元的全局编号、每个单元的20个节点的全局编号和所有节点的坐标信息；

步骤1-3、读入目标模型的几何参数，并将每个曲六面体单元中20个节点分别对应到参量坐标系中的标准立方体单元中，建立实际坐标系与参量坐标系的几何映射关系；其中，几何映射关系的表达式为：

$x=\underset{i=1}{\overset{20}{\mathrm{\Σ}}}{P}_i(\mathrm{\ξ},\mathrm{\η},\mathrm{\ζ})x_i$

$y=\underset{i=1}{\overset{20}{\mathrm{\Σ}}}{P}_i(\mathrm{\ξ},\mathrm{\η},\mathrm{\ζ})y_i---\left(1\right)$

$z=\underset{i=1}{\overset{20}{\mathrm{\Σ}}}{P}_i(\mathrm{\ξ},\mathrm{\η},\mathrm{\ζ})z_i$

$\begin{array}{c}{P}_i(\mathrm{\ξ},\mathrm{\η},\mathrm{\ζ})=\frac{{\mathrm{\ξ}}_i^2{\mathrm{\η}}_i^2{\mathrm{\ζ}}_i^2}{8}(1+{\mathrm{\ξ}}_i\mathrm{\ξ})(1+{\mathrm{\η}}_i\mathrm{\η})(1+{\mathrm{\ζ}}_i\mathrm{\ζ})({\mathrm{\ξ}}_i\mathrm{\ξ}+{\mathrm{\η}}_i\mathrm{\η}+{\mathrm{\ζ}}_i\mathrm{\ζ}-2)\\ +\frac{{\mathrm{\η}}_i^2{\mathrm{\ζ}}_i^2}{4}(1-{\mathrm{\ξ}}^2)(1+{\mathrm{\η}}_i\mathrm{\η})(1+{\mathrm{\ζ}}_i\mathrm{\ζ})(1-{\mathrm{\ξ}}_i^2)\\ +\frac{{\mathrm{\ξ}}_i^2{\mathrm{\ζ}}_i^2}{4}(1-{\mathrm{\η}}^2)(1+{\mathrm{\ζ}}_i\mathrm{\ζ})(1+{\mathrm{\ξ}}_i\mathrm{\ξ})(1-{\mathrm{\η}}_i^2)\\ +\frac{{\mathrm{\ξ}}_i^2{\mathrm{\η}}_i^2}{4}(1-{\mathrm{\ζ}}^2)(1+{\mathrm{\ξ}}_i\mathrm{\ξ})(1+{\mathrm{\η}}_i\mathrm{\η})(1-{\mathrm{\ζ}}^2)\end{array},$

其中，(x_i,y_i,z_i)为实际坐标系中第i个节点的坐标，(ξ_i,η_i,ζ_i)为该节点在参量坐标系中的坐标，i表示曲六面体单元中节点的序号，1≤i≤20；(ξ,η,ζ)为参量坐标系中任意一点坐标，(x,y,z)为该点映射到实际直角坐标系中的坐标。

步骤2、选定时域谱元法中基函数的阶数，确定参量坐标系下标准立方体单元中高斯积分点的几何参数，根据上一步的几何映射关系，得到各个曲六面体单元中相应的实际高斯积分点的几何参数，建立参量坐标系与实际坐标系中基函数的映射关系；具体为：时域谱元法中，电场展开基函数采用Gauss-Lobatto-Legendre多项式，Gauss-Lobatto-Legendre多项式为高阶正交多项式，由该多项式组成的三维矢量基函数满足单元分界面上的切向连续性条件，其表达式为：

${\mathrm{\Φ}}_{\mathrm{rst}}^{\mathrm{\ξ}}=\widehat{\mathrm{\ξ}}{\mathrm{\φ}}_r^{\left(N_{\mathrm{\ξ}}\right)}\left(\mathrm{\ξ}\right){\mathrm{\φ}}_s^{\left(N_{\mathrm{\η}}\right)}\left(\mathrm{\η}\right){\mathrm{\φ}}_t^{\left(N_{\mathrm{\ζ}}\right)}\left(\mathrm{\ζ}\right)$

${\mathrm{\Φ}}_{\mathrm{rst}}^{\mathrm{\η}}=\widehat{\mathrm{\η}}{\mathrm{\φ}}_r^{\left(N_{\mathrm{\ξ}}\right)}\left(\mathrm{\ξ}\right){\mathrm{\φ}}_s^{\left(N_{\mathrm{\η}}\right)}\left(\mathrm{\η}\right){\mathrm{\φ}}_t^{\left(N_{\mathrm{\ζ}}\right)}\left(\mathrm{\ζ}\right)---\left(2\right)$

${\mathrm{\Φ}}_{\mathrm{rst}}^{\mathrm{\ζ}}=\widehat{\mathrm{\ζ}}{\mathrm{\φ}}_r^{\left(N_{\mathrm{\ξ}}\right)}\left(\mathrm{\ξ}\right){\mathrm{\φ}}_s^{\left(N_{\mathrm{\η}}\right)}\left(\mathrm{\η}\right){\mathrm{\φ}}_t^{\left(N_{\mathrm{\ζ}}\right)}\left(\mathrm{\ζ}\right)$

这组基函数的定义域为步骤1中的标准立方体单元；其中，分别表示参量坐标系中沿ξ、η、ζ方向的电场展开矢量基函数，分别表示参量坐标系中沿ξ、η、ζ方向的单位向量；另外，分别为ξ轴、η轴和ζ轴上的一维标量基函数，其定义域为(-1，1)，表达式为：

${\mathrm{\φ}}_r^{\left(N_{\mathrm{\ξ}}\right)}\left(\mathrm{\ξ}\right)=\frac{-1}{N_{\mathrm{\ξ}}(N_{\mathrm{\ξ}}+1)L_{N_{\mathrm{\ξ}}}\left({\mathrm{\ξ}}_r\right)}\frac{(1-{\mathrm{\ξ}}^2)L_{N_{\mathrm{\ξ}}}^{\′}\left(\mathrm{\ξ}\right)}{\mathrm{\ξ}{-\mathrm{\ξ}}_r}$

${\mathrm{\φ}}_s^{\left(N_{\mathrm{\η}}\right)}\left(\mathrm{\η}\right)=\frac{-1}{N_{\mathrm{\η}}(N_{\mathrm{\η}}+1)L_{N_{\mathrm{\η\ξ}}}\left({\mathrm{\η}}_s\right)}\frac{(1-{\mathrm{\η}}^2)L_{N_{\mathrm{\η}}}^{\′}\left(\mathrm{\η}\right)}{{\mathrm{\η}-\mathrm{\η\ξ}}_s}---\left(3\right)$

${\mathrm{\φ}}_t^{\left(N_{\mathrm{\ζ}}\right)}\left(\mathrm{\ζ}\right)=\frac{-1}{N_{\mathrm{\ζ\ξ}}(N_{\mathrm{\ζ}}+1)L_{N_{\mathrm{\ζ\ξ}}}\left({\mathrm{\ζ\ξ}}_t\right)}\frac{(1-{\mathrm{\ζ\ξ}}^2)L_{N_{\mathrm{\ζ}}}^{\′}\left(\mathrm{\ζ}\right)}{\mathrm{\ζ\ξ}{-\mathrm{\ζ}}_t}$

其中，N_ξ、N_η、N_ζ分别表示在ξ、η、ζ轴上，一维标量基函数的阶数；ξ_r、η_s、ζ_t中的下标r、s、t分别表示在ξ、η、ζ方向上，一维标量基函数插值积分点的编号，其中，r＝1～N_ξ，s＝1～N_η，t＝1～N_ζ，r、s、t的组合可以唯一地表示单元体内三维矢量基函数的积分点所在的位置；表示在ξ轴上的N_ξ阶Legendre多项式，为其一阶导数；表示在η轴上的N_η阶Legendre多项式，为其一阶导数；表示在ζ轴上的N_ζ阶Legendre多项式，为其一阶导数；

实际坐标系与参量坐标系中基函数的映射关系如式(4)所示：

N_j＝J^-1Φ_j

$\▿\×N_j=\frac{1}{\left|J\right|}{J}^T\▿\×{\mathrm{\Φ}}_j---\left(4\right)$

其中，Φ_j表示参量坐标系中第j个基函数，j＝1～3(N_ξ+1)(N_η+1)(N_ζ+1)；Ν_j为实际直角坐标系中的对应基函数，J为雅克比矩阵，由表示几何映射关系的式(1)可得到J的表达式：

$J=\left[\begin{array}{ccc}\frac{\∂x}{\∂\mathrm{\ξ}}& \frac{\∂y}{\∂\mathrm{\ξ}}& \frac{\∂z}{\∂\mathrm{\ξ}}\\ \frac{\∂x}{\∂\mathrm{\η}}& \frac{\∂y}{\∂\mathrm{\η}}& \frac{\∂z}{\∂\mathrm{\η}}\\ \frac{\∂x}{\∂\mathrm{\ζ}}& \frac{\∂y}{\∂\mathrm{\ζ}}& \frac{\∂z}{\∂\mathrm{\ζ}}\end{array}\right]---\left(5\right)$

通过式(4)、(5)，得到实际坐标系与参量坐标系中基函数的映射关系。

步骤3、将具有张量形式的、包含阻尼因子的铁氧体磁导率引入电场时域亥姆赫兹方程，并且以单轴电各向异性完美匹配层作为吸收边界条件，经过伽辽金变换后，在时间的离散上采用中心差分格式，得到时域电场迭代式；

将具有张量形式的、包含阻尼因子的铁氧体磁导率引入电场时域亥姆赫兹方程，该电场时域亥姆赫兹方程的表达式为：

$\▿\×\frac{1}{{\mathrm{\μ}}_0}{{\stackrel{=}{\mathrm{\μ}}}_r}^{-1}{\stackrel{=}{\mathrm{\Λ}}}^{-1}\▿\×E+{\mathrm{\ϵ}}_0{\mathrm{\ϵ}}_r\stackrel{=}{\mathrm{\Λ}}\frac{{\∂}^{2}E}{\∂t^2}=0---\left(6\right)$

其中μ₀为自由空间磁导率；ε₀为自由空间介电常数；ε_r为材料的相对介电常数；E为矢量电场；为具有张量形式的、单轴电各向异性完美匹配层的本构参数；为具有张量形式的、包含阻尼因子的铁氧体磁导率；

当外加恒定磁场为x轴方向时，的逆矩阵表达式为：

${\stackrel{=}{\mathrm{\μ}}}_r^{-1}=\left(\begin{array}{ccc}1& 0& 0\\ 0& \mathrm{\δ}& -\mathrm{\β}\\ 0& \mathrm{\β}& \mathrm{\δ}\end{array}\right)$

当外加恒定磁场为y轴方向时，的逆矩阵表达式为：

${{\stackrel{\‾}{\stackrel{\‾}{\mathrm{\μ}}}}_r}^{-1}=\left(\begin{array}{ccc}\mathrm{\δ}& 0& \mathrm{\β}\\ 0& 1& 0\\ -\mathrm{\β}& 0& \mathrm{\δ}\end{array}\right)$

当外加恒定磁场为z轴方向时，的逆矩阵表达式为：

${{\stackrel{\‾}{\stackrel{\‾}{\mathrm{\μ}}}}_r}^{-1}=\left(\begin{array}{ccc}\mathrm{\δ}& -\mathrm{\β}& 0\\ \mathrm{\β}& \mathrm{\δ}& 0\\ 0& 0& 1\end{array}\right)$

其中，

$\mathrm{\δ}=\frac{{({\mathrm{\ω}}_0+\mathrm{j\ω\α})}^2-{\mathrm{\ω}}^2+({\mathrm{\ω}}_0+\mathrm{j\ω\α}){\mathrm{\ω}}_m}{{(({\mathrm{\ω}}_0+\mathrm{j\ω\α})+{\mathrm{\ω}}_m)}^2-{\mathrm{\ω}}^2}$

$\mathrm{\β}=\frac{\mathrm{j\ω}{\mathrm{\ω}}_m}{{(({\mathrm{\ω}}_0+\mathrm{j\ω\α})+{\mathrm{\ω}}_m)}^2-{\mathrm{\ω}}^2}$

α为阻尼因子；进动角频率ω₀＝γH₀，γ为旋磁比，H₀为外加恒定磁场的强度；磁化角频率ω_m＝γ4πM_s，M_s为饱和磁化率；

当外加恒定磁场与坐标轴的方向相反时，β表达式取相反数。

描绘铁氧体元器件内部电场分布的时域亥姆赫兹方程的通用表达式为：

$\begin{array}{c}\frac{{\mathrm{\α}}^2+1}{{\mathrm{\μ}}_0}\▿\×\▿\×E+\frac{1}{{\mathrm{\μ}}_0}\▿\×B\▿\×E^{\′}+\frac{2}{{\mathrm{\ϵ}}_0{\mathrm{\μ}}_0}\▿\×R\▿\×E^{\′}+\frac{1}{{{\mathrm{\ϵ}}_0}^2{\mathrm{\μ}}_0}\▿\×R^2\▿\×E^{\′\′}\\ +\frac{1}{{\mathrm{\μ}}_0}\▿\×A\▿\×E^{\′\′}+\mathrm{\ϵ}{({\mathrm{\ω}}_0+{\mathrm{\ω}}_m)}^{2}E+\frac{(1+{\mathrm{\α}}^2){\mathrm{\ϵ}}_r}{{\mathrm{\ϵ}}_0}{K}^{2}E\\ +2\mathrm{\α}({\mathrm{\ω}}_0+{\mathrm{\ω}}_m)\mathrm{\ϵ}\frac{\∂E}{\∂t}+2(1+{\mathrm{\α}}^2){\mathrm{\ϵ}}_{r}K\frac{\∂E}{\∂t}+\mathrm{\ϵ}(1+{\mathrm{\α}}^2)\frac{{\∂}^{2}E}{\∂t^2}=0\end{array}---\left(7\right)$

其中_，E′＝∫_tE,E″＝∫∫_tE

$A=\left(\begin{array}{ccc}{{\mathrm{\ω}}_0}^2+{\mathrm{\ω}}_0{\mathrm{\ω}}_m& & \\ & {({\mathrm{\ω}}_0+{\mathrm{\ω}}_m)}^2& \\ & & {{\mathrm{\ω}}_0}^2+{\mathrm{\ω}}_0{\mathrm{\ω}}_m\end{array}\right)$

$B=\left(\begin{array}{ccc}2\mathrm{\α}{\mathrm{\ω}}_0+\mathrm{\α}{\mathrm{\ω}}_m& 0& -{\mathrm{\ω}}_m\\ 0& 2\mathrm{\α}{\mathrm{\ω}}_0+2\mathrm{\α}{\mathrm{\ω}}_m& 0\\ {\mathrm{\ω}}_m& 0& 2\mathrm{\α}{\mathrm{\ω}}_0+\mathrm{\α}{\mathrm{\ω}}_m\end{array}\right)$

$K=\left(\begin{array}{ccc}{\mathrm{\σ}}_z& & \\ & {\mathrm{\σ}}_z& \\ & & 0\end{array}\right),R=\left(\begin{array}{ccc}0& & \\ & 0& \\ & & {\mathrm{\σ}}_z\end{array}\right)$

其中，σ_z为单轴电各向异性完美匹配层设置在z轴方向时的本构参数，通式(7)既适用于电、磁各向同性的普通媒质，也适用于具有磁各向异性的铁氧体材料，还适用于具有单轴电各向异性完美匹配层的吸收边界；

通式(7)经伽辽金变换，并化简得到

$(S+T_p)e+S_1{e}^{\′}+S_{\mathrm{th}}{e}^{\′\′}+T_q\frac{\∂e}{\∂t}+T\frac{{\∂}^{2}e}{\∂t^2}=0$

其中，

$S=\frac{{\mathrm{\α}}^2+1}{{\mathrm{\μ}}_0}\∫\∫\∫\▿\×N_k\·\▿\×N_j\mathrm{dV}$

$S=\frac{{\mathrm{\α}}^2+1}{{\mathrm{\μ}}_0}\∫\∫\∫\▿\×N_k\·(\frac{1}{{\mathrm{\μ}}_0}B+\frac{2}{{\mathrm{\ϵ}}_0{\mathrm{\μ}}_0}R)\·\▿\×N_j\mathrm{dV}$

$S_{\mathrm{th}}=\∫\∫\∫\▿\×N_k\·(\frac{1}{{\mathrm{\μ}}_0}A+\frac{2}{{{\mathrm{\ϵ}}_0}^2{\mathrm{\μ}}_0}{R}^2)\·\▿\×N_j\mathrm{dV}$

T＝ε(1+α²)∫∫∫N_k·N_jdV

$T_p=\∫\∫\∫N_k\·[\mathrm{\ϵ}{({\mathrm{\ω}}_0+{\mathrm{\ω}}_m)}^2+\frac{(1+{\mathrm{\α}}^2){\mathrm{\ϵ}}_r}{{\mathrm{\ϵ}}_0}{K}^2]\·N_j\mathrm{dV}$

T_q＝∫∫∫N_k·[2α(ω₀+ω_m)ε+2(1+α²)ε_rK]·N_jdV

其中，N_k为实际坐标系中的第k个测试基函数，N_j为实际坐标系中第j个基函数；在谱元法应用中，利用式(4)将其转换到参量坐标系中即可；

将式(8)采用中心差分格式展开，得到

(0.5ΔtT_q+T)eⁿ⁺¹＝(2T-Δt²S-Δt²T_p)eⁿ+(0.5ΔtT_q-T)e^n-1

                                                   (9)

-Δt²S₁e′ⁿ-Δt²S_the″ⁿ

式(9)为微波铁氧体元件内部电场的时间迭代式，其中，Δt为时间步进，n为时间步，其最大取值必须保证经过n个时间步的运算后，信号能够从源的位置传播到观察点，并保持稳定；e^n-1、eⁿ、eⁿ⁺¹分别表示第n-1时间步、第n时间步和第n+1时间步的电场，e′ⁿ、e″ⁿ分别表示第n时间步电场对时间的一次积分和二次积分，其递推关系为：

e′ⁿ＝e′^n-1+eⁿΔt

e″ⁿ＝e″^n-1+e′ⁿΔt。

步骤4、利用时域电场迭代式对目标铁氧体元器件进行两次时域仿真；第一次仿真时，将铁氧体材质退化为普通材质，使元件内部在材质上保持连续，记录观察点上的时域仿真数据，当时域仿真数据稳定后，通过FFT变换，获得入射波信号的频率响应；第二次仿真时，恢复铁氧体材质，记录观察点上的时域仿真数据，当时域仿真数据稳定后，通过FFT变换，获得前向和后向传输信号的频率响应；

步骤5、利用上一步中两次仿真获得的频率响应，根据微波网络散射参量的定义，确定目标铁氧体元件在工作频段内各个端口的插入损耗、回波损耗和隔离度。

下面以一个Y形的铁氧体波导环行器为例，详细阐述本发明的具体实施方式。

实施例1

结合图1，一种微波铁氧体元器件的时域谱元仿真方法，其中微波铁氧体元器件为Y形波导环形器，如图2所示，三段波导的口径均为22.86mm×10.16mm，之间互为120°设置，三段波导的交汇处设置铁氧体柱。铁氧体柱的半径为3.5mm，高为10.16mm；铁氧体材料选用TT1-109，参数为ε_rf＝11.7，4πM_s＝1317G，H₀＝200Oe。该时域谱元仿真方法具体包括以下步骤：

步骤1、建立目标的几何模型：根据微波铁氧体元器件的结构尺寸，用计算机辅助设计工具(如ANSYS)进行建模，并采用20个点的曲六面体单元对目标模型进行离散，得到所有单元的几何参数信息，如图3所示；在三段波导的端口处各设置单轴各向异性完美匹配层(PML)；其余部分为空气，模型最外层的边界为理想电壁。

读入目标模型的几何参数，每个曲六面体有20个节点坐标，前8个节点坐标分别对应曲六面体的8个顶点，后12个节点坐标分别对应曲六面体的12条棱边的中点；按每次读入的节点的顺序，依次将20个节点分别对应到参量坐标系中一个每边都为-1～1的标准立方体中，建立实际坐标系与参量坐标系的几何映射关系，如图4所示；

步骤2、根据基函数的阶数，确定参量坐标系下标准立方体高斯点的几何参数；当三个方向的基函数均为2阶时，GLL积分点与GLL基函数在标准参量立方体单元中的分布如图5所示；

步骤3、将具有张量形式的、包含阻尼因子的铁氧体磁导率引入电场时域亥姆赫兹方程，并且以单轴电各向异性完美匹配层作为吸收边界条件，经过伽辽金变换后，在时间的离散上采用中心差分格式，得到时域电场迭代式；

步骤4、利用上一步所推导出的时域电场迭代通式，对Y形铁氧体环行器做两次时域仿真，并且在3个端口各设置一个观察点，调制高斯脉冲信号由端口1输入；第一次仿真时，将铁氧体材质退化为普通材质，使元件内部在材质上保持连续，记录观察点1上的时域仿真数据；通过这一次仿真，可以得到没有任何干扰情况下的纯的输入信号V_a1(t)，完成预定的m步计算后，通过FFT变换，获得入射波信号的频率响应V_a1(f)；第二次仿真时，恢复铁氧体材质，记录3个观察点上的时域仿真数据V_b1(t)，V_b2(t)，V_b3(t)；完成预定的n步计算后，通过FFT变换，获得前向和后向传输波信号的频率响应V_b1(f)，V_b2(f)，V_b2(f)；图6为第二次仿真时，信号从端口1输入，另两个端口输出，在3个端口观察点处的信号时域波形，显然，端口3的输出信号较之端口2的输出信号大很多，因此端口3为该环形器的信号传输端，端口2为信号隔离端。

步骤5、根据散射参量的定义，确定Y形铁氧体环行器的插入损耗S31、回波损耗S11和隔离度S21：

$S11=-20\log \frac{V_{b1}\left(f\right)-V_{a1}\left(f\right)}{V_{a1}\left(f\right)}\mathrm{dB}$

$S21=-20\log \frac{V_{b2}\left(f\right)}{V_{a1}\left(f\right)}\mathrm{dB}$

$S31=-20\log \frac{V_{b3}\left(f\right)}{V_{a1}\left(f\right)}\mathrm{dB}$

经计算，该Y形铁氧体波导环形器的S11、S21、S31如图7～9所示。显然，该Y形铁氧体波导环形器的工作频率范围是9.8GHz～10.2GHz，此时，信号从端口1输入，从端口3输出，端口2为隔离端口。仿真结果与使用FDTD方法计算的结果基本吻合。

以上所述仅为本发明的优选实施例，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (5)

1.一种微波铁氧体元器件的时域谱元仿真方法，其特征在于，包括以下步骤：

步骤1、根据目标铁氧体元器件的结构尺寸，建立实际坐标系下的电磁空间模型；采用曲六面体单元对该电磁空间模型进行离散，依据每个离散单元的几何参数，建立实际坐标系中各个曲六面体单元与参量坐标系下标准立方体单元的几何映射关系；所述参量坐标系是以标准立方体单元的重心为原点，以标准立方体单元的棱边为坐标轴方向的直角坐标系；

步骤2、选定时域谱元法中基函数的阶数，确定参量坐标系下标准立方体单元中高斯积分点的几何参数，根据上一步的几何映射关系，得到各个曲六面体单元中相应的实际高斯积分点的几何参数，建立参量坐标系与实际坐标系中基函数的映射关系；

步骤3、将具有张量形式的、包含阻尼因子的铁氧体磁导率引入电场时域亥姆赫兹方程，并且以单轴电各向异性完美匹配层作为吸收边界条件，经过伽辽金变换后，在时间的离散上采用中心差分格式，得到时域电场迭代式；

步骤4、利用时域电场迭代式对目标铁氧体元器件进行两次时域仿真；第一次仿真时，将铁氧体材质退化为普通材质，使元件内部在材质上保持连续，记录观察点上的时域仿真数据，当时域仿真数据稳定后，通过FFT变换，获得入射波信号的频率响应；第二次仿真时，恢复铁氧体材质，记录观察点上的时域仿真数据，当时域仿真数据稳定后，通过FFT变换，获得前向和后向传输信号的频率响应；

步骤5、利用上一步中两次仿真获得的频率响应，根据微波网络散射参量的定义，确定目标铁氧体元件在工作频段内各个端口的插入损耗、回波损耗和隔离度。

2.根据权利要求1所述的微波铁氧体元器件的时域谱元仿真方法，其特征在于，步骤1具体为：

步骤1-1、根据微波铁氧体元器件的结构尺寸，用计算机辅助设计工具建立实际坐标系下的电磁空间模型；所述实际坐标系为三维直角坐标系，反映了目标铁氧体元器件在真实空间中的位置，实际坐标系的原点和三坐标轴可任意设置，目标铁氧体元器件的外加恒定磁场必须与三坐标轴中的某一坐标轴相平行；

步骤1-2、采用20个节点的曲六面体单元对目标模型进行全局离散，得到所有单元的几何参数；所述20个节点为所在曲六面体单元的8个顶点和12条棱边的中点，所述几何参数包括单元的全局编号、每个单元的20个节点的全局编号和所有节点的坐标信息；

步骤1-3、读入目标模型的几何参数，并将每个曲六面体单元中20个节点分别对应到参量坐标系中的标准立方体单元中，建立实际坐标系与参量坐标系的几何映射关系；其中，几何映射关系的表达式为：

$x=\underset{i=1}{\overset{20}{\mathrm{\Σ}}}{P}_i(\mathrm{\ξ},\mathrm{\η},\mathrm{\ζ})x_i$

$y=\underset{i=1}{\overset{20}{\mathrm{\Σ}}}{P}_i(\mathrm{\ξ},\mathrm{\η},\mathrm{\ζ})y_i---\left(1\right)$

$z=\underset{i=1}{\overset{20}{\mathrm{\Σ}}}{P}_i(\mathrm{\ξ},\mathrm{\η},\mathrm{\ζ})z_i$

$\begin{array}{c}{P}_i(\mathrm{\ξ},\mathrm{\η},\mathrm{\ζ})=\frac{{\mathrm{\ξ}}_i^2{\mathrm{\η}}_i^2{\mathrm{\ζ}}_i^2}{8}(1+{\mathrm{\ξ}}_i\mathrm{\ξ})(1+{\mathrm{\η}}_i\mathrm{\η})(1+{\mathrm{\ζ}}_i\mathrm{\ζ})({\mathrm{\ξ}}_i\mathrm{\ξ}+{\mathrm{\η}}_i\mathrm{\η}+{\mathrm{\ζ}}_i\mathrm{\ζ}-2)\\ +\frac{{\mathrm{\η}}_i^2{\mathrm{\ζ}}_i^2}{4}(1-{\mathrm{\ξ}}^2)(1+{\mathrm{\η}}_i\mathrm{\η})(1+{\mathrm{\ζ}}_i\mathrm{\ζ})(1-{\mathrm{\ξ}}_i^2)\\ +\frac{{\mathrm{\ξ}}_i^2{\mathrm{\ζ}}_i^2}{4}(1-{\mathrm{\η}}^2)(1+{\mathrm{\ζ}}_i\mathrm{\ζ})(1+{\mathrm{\ξ}}_i\mathrm{\ξ})(1-{\mathrm{\η}}_i^2)\\ +\frac{{\mathrm{\ξ}}_i^2{\mathrm{\η}}_i^2}{4}(1-{\mathrm{\ζ}}^2)(1+{\mathrm{\ξ}}_i\mathrm{\ξ})(1+{\mathrm{\η}}_i\mathrm{\η})(1-{\mathrm{\ζ}}_i^2)\end{array},$

其中，(x_i,y_i,z_i)为实际坐标系中第i个节点的坐标，(ξ_i,η_i,ζ_i)为该节点在参量坐标系中的坐标，i表示曲六面体单元中节点的序号，1≤i≤20；(ξ,η,ζ)为参量坐标系中任意一点坐标，(x,y,z)为该点映射到实际直角坐标系中的坐标。

3.根据权利要求2所述的微波铁氧体元器件的时域谱元仿真方法，其特征在于，步骤2具体包括以下步骤：

时域谱元法中，电场展开基函数采用Gauss-Lobatto-Legendre多项式，Gauss-Lobatto-Legendre多项式为高阶正交多项式，由该多项式组成的三维矢量基函数满足单元分界面上的切向连续性条件，其表达式为：

${\mathrm{\Φ}}_{\mathrm{rst}}^{\mathrm{\ξ}}={\widehat{\mathrm{\ξ}}\mathrm{\φ}}_r^{\left(N_{\mathrm{\ξ}}\right)}\left(\mathrm{\ξ}\right){\mathrm{\φ}}_s^{\left(N_{\mathrm{\η}}\right)}\left(\mathrm{\η}\right){\mathrm{\φ}}_t^{\left(N_{\mathrm{\ζ}}\right)}\left(\mathrm{\ζ}\right)$

${\mathrm{\Φ}}_{\mathrm{rst}}^{\mathrm{\η}}={\widehat{\mathrm{\η}}\mathrm{\φ}}_r^{\left(N_{\mathrm{\ξ}}\right)}\left(\mathrm{\ξ}\right){\mathrm{\φ}}_s^{\left(N_{\mathrm{\η}}\right)}\left(\mathrm{\η}\right){\mathrm{\φ}}_t^{\left(N_{\mathrm{\ζ}}\right)}\left(\mathrm{\ζ}\right)---\left(2\right)$

${\mathrm{\Φ}}_{\mathrm{rst}}^{\mathrm{\ζ}}={\widehat{\mathrm{\ζ}}\mathrm{\φ}}_r^{\left(N_{\mathrm{\ξ}}\right)}\left(\mathrm{\ξ}\right){\mathrm{\φ}}_s^{\left(N_{\mathrm{\η}}\right)}\left(\mathrm{\η}\right){\mathrm{\φ}}_t^{\left(N_{\mathrm{\ζ}}\right)}\left(\mathrm{\ζ}\right)$

这组基函数的定义域为步骤1中的标准立方体单元；其中，分别表示参量坐标系中沿ξ、η、ζ方向的电场展开矢量基函数，分别表示参量坐标系中沿ξ、η、ζ方向的单位向量；另外，分别为ξ轴、η轴和ζ轴上的一维标量基函数，其定义域为(-1，1)，表达式为：

${\mathrm{\φ}}_r^{\left(N_{\mathrm{\ξ}}\right)}\left(\mathrm{\ξ}\right)=\frac{-1}{N_{\mathrm{\ξ}}(N_{\mathrm{\ξ}}+1)L_{N_{\mathrm{\ξ}}}\left({\mathrm{\ξ}}_r\right)}\frac{(1-{\mathrm{\ξ}}^2)L_{N_{\mathrm{\ξ}}}^{\′}\left(\mathrm{\ξ}\right)}{\mathrm{\ξ}-{\mathrm{\ξ}}_r}$

${\mathrm{\φ}}_s^{\left(N_{\mathrm{\η}}\right)}\left(\mathrm{\η}\right)=\frac{-1}{N_{\mathrm{\η\ξ}}(N_{\mathrm{\η}}+1)L_{N_{\mathrm{\η}}}\left({\mathrm{\η}}_s\right)}\frac{(1-{\mathrm{\η\ξ}}^2)L_{N_{\mathrm{\η}}}^{\′}\left(\mathrm{\η}\right)}{\mathrm{\η}-{\mathrm{\η}}_s}---\left(3\right)$

${\mathrm{\φ}}_t^{\left(N_{\mathrm{\ζ}}\right)}\left(\mathrm{\ζ}\right)=\frac{-1}{N_{\mathrm{\ζ}}(N_{\mathrm{\ζ}}+1)L_{N_{\mathrm{\ζ}}}\left({\mathrm{\ζ}}_r\right)}\frac{(1-{\mathrm{\ζ}}^2)L_{N_{\mathrm{\ζ}}}^{\′}\left(\mathrm{\ζ}\right)}{\mathrm{\ζ}-{\mathrm{\ζ}}_t}$

其中，N_ξ、N_η、N_ζ分别表示在ξ、η、ζ方向上，一维标量基函数的阶数；ξ_r、η_s、ζ_t中的下标r、s、t分别表示在ξ、η、ζ方向上，一维标量基函数插值积分点的编号，其中，r＝1～N_ξ，s＝1～N_η，t＝1～N_ζ；表示在ξ轴上的N_ξ阶Legendre多项式，为其一阶导数；表示在η轴上的N_η阶Legendre多项式，为其一阶导数；表示在ζ轴上的N_ζ阶Legendre多项式，为其一阶导数；

实际坐标系与参量坐标系中基函数的映射关系如式(4)所示：

N_j＝J^-1Φ_j

$\▿\×N_j=\frac{1}{\left|j\right|}{J}^T\▿\×{\mathrm{\Φ}}_j---\left(4\right)$

其中，Φ_j表示参量坐标系中第j个基函数，j＝1～3(N_ξ+1)(N_η+1)(N_ζ+1)；Ν_j为实际直角坐标系中的对应基函数，J为雅克比矩阵，由表示几何映射关系的式(1)可得到J的表达式：

$J=\left[\begin{array}{ccc}\frac{\∂x}{\∂\mathrm{\ξ}}& \frac{\∂y}{\∂\mathrm{\ξ}}& \frac{\∂z}{\∂\mathrm{\ξ}}\\ \frac{\∂x}{\∂\mathrm{\η}}& \frac{\∂y}{\∂\mathrm{\η}}& \frac{\∂z}{\∂\mathrm{\η}}\\ \frac{\∂x}{\∂\mathrm{\ζ}}& \frac{\∂y}{\∂\mathrm{\ζ}}& \frac{\∂z}{\∂\mathrm{\ζ}}\end{array}\right]---\left(5\right)$

通过式(4)、(5)，得到实际坐标系与参量坐标系中基函数的映射关系。

4.根据权利要求3所述的微波铁氧体元器件的时域谱元仿真方法，其特征在于，步骤3中将具有张量形式的、包含阻尼因子的铁氧体磁导率引入电场时域亥姆赫兹方程，该电场时域亥姆赫兹方程的表达式为：

$\▿\×\frac{1}{{\mathrm{\μ}}_0}{{\stackrel{=}{\mathrm{\μ}}}_r}^{-1}{\stackrel{=}{\mathrm{\Λ}}}^{-1}\▿\×E+{\mathrm{\ϵ}}_0{\mathrm{\ϵ}}_r\stackrel{=}{\mathrm{\Λ}}\frac{{\∂}^{2}E}{{\∂t}^2}=0---\left(6\right)$

其中μ₀为自由空间磁导率；ε₀为自由空间介电常数；ε_r为材料的相对介电常数；E为矢量电场；为具有张量形式的、单轴电各向异性完美匹配层的本构参数；为具有张量形式的、包含阻尼因子的铁氧体磁导率；

当外加恒定磁场为x轴方向时，的逆矩阵表达式为：

${{\stackrel{=}{\mathrm{\μ}}}_r}^{-1}=\left(\begin{array}{ccc}1& 0& 0\\ 0& \mathrm{\δ}& -\mathrm{\β}\\ 0& \mathrm{\β}& \mathrm{\δ}\end{array}\right)$

当外加恒定磁场为y轴方向时，的逆矩阵表达式为：

${{\stackrel{=}{\mathrm{\μ}}}_r}^{-1}=\left(\begin{array}{ccc}\mathrm{\δ}& 0& \mathrm{\β}\\ 0& 1& 0\\ -\mathrm{\β}& 0& \mathrm{\δ}\end{array}\right)$

当外加恒定磁场为z轴方向时，的逆矩阵表达式为：

${{\stackrel{=}{\mathrm{\μ}}}_r}^{-1}=\left(\begin{array}{ccc}\mathrm{\δ}& -\mathrm{\β}& 0\\ \mathrm{\β}& \mathrm{\δ}& 0\\ 0& 0& 1\end{array}\right)$

其中，

$\mathrm{\δ}=\frac{{({\mathrm{\ω}}_0+\mathrm{j\ω\α})}^2-{\mathrm{\ω}}^2+({\mathrm{\ω}}_0+\mathrm{j\ω\α}){\mathrm{\ω}}_m}{{(({\mathrm{\ω}}_0+\mathrm{j\ω\α})+{\mathrm{\ω}}_m)}^2-{\mathrm{\ω}}^2}$

$\mathrm{\β}=\frac{\mathrm{j\ω}{\mathrm{\ω}}_m}{{(({\mathrm{\ω}}_0+\mathrm{j\ω\α})+{\mathrm{\ω}}_m)}^2-{\mathrm{\ω}}^2}$

α为阻尼因子；进动角频率ω₀＝γH₀，γ为旋磁比，H₀为外加恒定磁场的强度；磁化角频率ω_m＝γ4πM_s，M_s为饱和磁化率；

当外加恒定磁场与坐标轴的方向相反时，β表达式取相反数。

5.根据权利要求4所述的微波铁氧体元器件的时域谱元仿真方法，其特征在于，描绘铁氧体元器件内部电场分布的时域亥姆赫兹方程的通用表达式为：

$\begin{array}{c}\frac{{\mathrm{\α}}^2+1}{{\mathrm{\μ}}_0}\▿\×\▿\×E+\frac{1}{{\mathrm{\μ}}_0}\▿\×B\▿\×E^{\′}+\frac{2}{{\mathrm{\ϵ}}_0{\mathrm{\μ}}_0}\▿\×R\▿\×E^{\′}+\frac{1}{{{\mathrm{\ϵ}}_0}^2{\mathrm{\μ}}_0}\▿\×R^2\▿\×E^{\′\′}\\ +\frac{1}{{\mathrm{\μ}}_0}\▿\×A\▿\×E^{\′\′}+\mathrm{\ϵ}{({\mathrm{\ω}}_0+{\mathrm{\ω}}_m)}^{2}E+\frac{(1+{\mathrm{\α}}^2){\mathrm{\ϵ}}_r}{{\mathrm{\ϵ}}_0}{K}^{2}E\\ +2\mathrm{\α}({\mathrm{\ω}}_0+{\mathrm{\ω}}_m)\mathrm{\ϵ}\frac{\∂E}{\∂t}+2(1+{\mathrm{\α}}^2){\mathrm{\ϵ}}_{r}K\frac{\∂E}{\∂t}+\mathrm{\ϵ}(1+{\mathrm{\α}}^2)\frac{{\∂}^{2}E}{{\∂t}^2}=0\end{array}---\left(7\right)$

其中，E′＝∫_tE,E″＝∫∫_tE

$A=\left(\begin{array}{ccc}{{\mathrm{\ω}}_0}^2+{\mathrm{\ω}}_0{\mathrm{\ω}}_m& & \\ & {({\mathrm{\ω}}_0+{\mathrm{\ω}}_m)}^2& \\ & & {{\mathrm{\ω}}_0}^2+{\mathrm{\ω}}_0{\mathrm{\ω}}_m\end{array}\right)$

$B=\left(\begin{array}{ccc}2{\mathrm{\α\ω}}_0+{\mathrm{\α\ω}}_m& 0& -{\mathrm{\ω}}_m\\ 0& 2{\mathrm{\α\ω}}_0+2{\mathrm{\α\ω}}_m& 0\\ {\mathrm{\ω}}_m& 0& 2{\mathrm{\α\ω}}_0+{\mathrm{\α\ω}}_m\end{array}\right)$

$K=\left(\begin{array}{ccc}{\mathrm{\σ}}_z& & \\ & {\mathrm{\σ}}_z& \\ & & 0\end{array}\right),R=\left(\begin{array}{ccc}0& & \\ & 0& \\ & & {\mathrm{\σ}}_z\end{array}\right)$

其中，σ_z为单轴电各向异性完美匹配层设置在z轴方向时的本构参数，通式(7)既适用于电、磁各向同性的普通媒质，也适用于具有磁各向异性的铁氧体材料，还适用于具有单轴电各向异性完美匹配层的吸收边界；

通式(7)经伽辽金变换，并化简得到

$(S+T_p)e+S_1{e}^{\′}+S_{\mathrm{th}}{e}^{\′\′}+T_q\frac{\∂e}{\∂t}+T\frac{{\∂}^{2}e}{{\∂t}^2}=0---\left(8\right)$

其中，

$S=\frac{{\mathrm{\α}}^2+1}{{\mathrm{\μ}}_0}\∫\∫\∫\▿\×N_k\·\▿\×N_j\mathrm{dV}$

$S_1=\∫\∫\∫\▿\×N_k\·(\frac{1}{{\mathrm{\μ}}_0}B+\frac{2}{{\mathrm{\ϵ}}_0{\mathrm{\μ}}_0}R)\·\▿\×N_j\mathrm{dV}$

$S_{\mathrm{th}}=\∫\∫\∫\▿\×N_k\·(\frac{1}{{\mathrm{\μ}}_0}A+\frac{2}{{{\mathrm{\ϵ}}_0}^2{\mathrm{\μ}}_0}{R}^2)\·\▿\×N_j\mathrm{dV}$

T＝ε(1+α²)∫∫∫N_k·N_jdV

$T_p=\∫\∫\∫N_k\·[\mathrm{\ϵ}{({\mathrm{\ω}}_0+{\mathrm{\ω}}_m)}^2+\frac{(1+{\mathrm{\α}}^2){\mathrm{\ϵ}}_r}{{\mathrm{\ϵ}}_0}{K}^2]\·N_j\mathrm{dV}$

T_q＝∫∫∫N_k·[2α(ω₀+ω_m)ε+2(1+α²)ε_rK]·N_jdV

其中，N_k为实际坐标系中的第k个测试基函数，N_j为实际坐标系中第j个基函数；将式(8)采用中心差分格式展开，得到

(0.5ΔtT_q+T)eⁿ⁺¹＝(2T-Δt²S-Δt²T_p)eⁿ+(0.5ΔtT_q-T)e^n-1

                                 (9)

-Δt²S₁e′ⁿ-Δt²S_the″ⁿ

式(9)为微波铁氧体元件内部电场的时间迭代式，其中，Δt为时间步进，n为时间步，其最大取值必须保证经过n个时间步的运算后，信号能够从源的位置传播到观察点，并保持稳定；e^n-1、eⁿ、eⁿ⁺¹分别表示第n-1时间步、第n时间步和第n+1时间步的电场，e′ⁿ、e″ⁿ分别表示第n时间步电场对时间的一次积分和二次积分，其递推关系为：

e′ⁿ＝e′^n-1+eⁿΔt

e″ⁿ＝e″^n-1+e′ⁿΔt。