基于时域谱元法的微波无源电路电磁热一体化分析方法

摘要

本发明公开了一种基于谱元法的微波无源电路电磁热一体化分析方法该方法首先采用时域谱元法求解麦克斯韦方程组，求出电路在高功率脉冲作用下瞬时的电场和磁场分布，得出当前时刻的电磁损耗。假设模型内部电磁损耗全部转化为热能，将所得热能带入热传导方程，得到当前时刻各点的温度分布情况。利用介电参数随温度变化的关系式得到下一时刻材料的电特性参数并再次计算电磁场方程，得出电磁损耗。如此反复循环，直到完成预定加热时间。电磁热一体化分析可以清楚的得到滤波器在不同脉冲的作用下，模型内部温度随时间变化的分布情况，建模灵活，剖分方便，形成的矩阵具有良好的稀疏性，求解效率较高。

说明书

技术领域

本发明属于微波无源器件的电磁热一体化分析，特别是针对SIW滤波器设计的数值分析方法。 

背景技术

上世纪70年代以来，由于计算机技术的飞速发展，计算电磁学获得了迅猛发展。相对于解析方法，数值方法能够解决很大一类计算巨大，结构复杂而解析方法难以或无法得到精确结果的问题。针对于微波无源电路，如滤波器、微带线的热效应作用过程，可以单独考虑为电磁求解过程和热模型求解过程。即分别对应于麦克斯韦电场波动方程和热传导方程的求解。研究微波热模型就是选择合适的数值分析方法求解分析这两大过程。对于电磁模型数值求解，目前主流方法按时频域分为时域方法(如时域有限差分法)和频域方法(如有限元，矩量法等)等。在应用上述方法分析微波热模型中的电磁问题时，有限元，矩量法，时域有限差分法(FDTD)等当前热门数值算法都可以使用，但考虑到微波热模型中电参数为时变函数和热模型时间延续性的特点，采用时域方法更为合适。一般FDTD，FEM等方法更为普遍。然而由于FDTD的Yee网格特性决定了其无法模拟结构复杂的模型。FEM应用到时域时每个时间步都涉及到对线性方程组的求解，计算量非常庞大，很浪费时间。 

发明内容

本发明的目的在于提供一种微波无源电路的电磁热一体化分析，同时考虑介质损耗和导体损耗（忽略辐射损耗）从而实现快速得到器件内部温度分布的方法。 

实现本发明目的的技术解决方案为：一种微波无源器件电热一体化分析方法，步骤表述如下： 

第一步，建立求解模型及网格剖分。并采用曲六面体对模型进行整体剖分，得到模型的结构信息，包括每个曲六面体单元的结点编号和坐标等。剖分网格的尺寸大于满足精度所需的剖分尺寸。 

第二步，从麦克斯韦方程方程组出发，对等式采用伽辽金法测试，强加边界条件， 得到电场波动方程，进而求解得到各节点t_n时刻的电场及电流分布。 

第三步，由以上步骤得到的电场及电流分布得出各点的导体损耗及介质损耗（忽略辐射损耗）。在求解介质损耗时用等效电导率σ'代替σ，即σ'＝σ+2πfε″。其中，ε″为复介电常数的虚部，f为工作频率。然后通过公式计算得到单位体积内介质损耗所转化的热能。求解导体损耗的时候，首先通过麦克斯韦方程由磁场分布计算得到电流密度分布。由于良导体中电流分布主要集中在导体表面。假设微波电路所敷金属贴片厚度即为电磁波进入金属的趋肤深度，计算得到表面电阻。由得到导体损耗。 

第四步，建立微波无源电路的热传导方程，将电磁损耗作为热源项代入该方程中，求解得到各节点温度分布； 

第五步，判断t_n时刻是否为预定加热时刻，若是预定加热时间则停止计算，若不是是预定加热时间则更新电场波动方程中的介电常数和电导率，计算t_n+1的电场分布及电流密度分布。重复步骤二、三、四，五步，如此反复循环，直到完成预定加热时间； 

本发明与现有技术相比，其显著优点：（1）SETD采用曲六面体剖分，建模灵活，剖分方便，使用用特定的正交多项式作为基函数，随着多项式阶数的提高，计算误差将呈指数下降。（2）对于电磁场场计算空间，得到的质量矩阵为块对角阵；可用块对角矩阵的求逆方法事先求出质量矩阵的逆，使整个方程的求解变为显式，降低计算量.（3）微波器件的电特性与热特性是一体化分析的，没有将其割裂开来。而且下一步可以将电热之间的相互影响联系起来。（4）热传导方程形成的矩阵方程性态较好，直接求逆相当方便，使用直接解法可快速求得器件内部的温度分布。 

附图说明

图1是SIW带通滤波器的模型图。 

图2是本发明计算得到带通滤波器在通过高斯调制脉冲时温度分布图。 

具体实施方式

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

一、传导方程简单推导 

在热媒质内(x,y,z)点，t时刻热量传输满足导热偏微分方程，即热传导方程(HTE)， 方程的形式为： 

${\rho }_m{c}_m\frac{\partial T(x,y,z)}{\partial t}=k_t{▿}^{2}T(x,y,z)-V_S(T(x,y,z)-T_a)+P_d(x,y,z,t)---\left(1.1\right)$

式中，T[℃]为物体的瞬态温度；t[s]为过程进行的时间；k_t[W/(m·℃)]为材料的导热系数；ρ_m[kg/m³]为材料的密度；c_m[J/(kg·℃)]为材料的定压比热；V_S[W/(m³·℃)]为冷却流的热容流积；T_a[℃]为冷却流的温度；P_d[W/m³]为内部热源的功率密度。 

为简化问题，假定ρ_m、c_m和k_t均不是时间、空间和温度的函数，且不考虑冷却流，在直角坐标系下，式(1)改写为 

$\frac{\partial T}{\partial t}=D_t((\frac{{\partial }^{2}T}{{\partial x}^2}+\frac{{\partial }^{2}T}{{\partial y}^2}+\frac{{\partial }^{2}T}{{\partial z}^2})+\frac{P_d}{k_t})---\left(1.2\right)$

上式中 

$D_t=\frac{k_t}{{\rho }_m{c}_m}---\left(1.3\right)$

D_t定义为热扩散率。需要说明的是，导热偏微分方程是描述导热过程共性的通用表达式，适用于任何导热过程。其中热源项P_d我们仅考虑介质损耗和导体损耗，忽略辐射损耗，将在第三节详细说明。 

传导方程的差分形式存在数值稳定性条件。为了防止随着迭代次数的增加而导致场值结果的溢出，在求解热模型时，时间步长和物体单元的网格大小必须要满足一定的准则。具体表达式为： 

$\mathrm{\Delta t}\le \frac{1}{{2D}_t}{(\frac{1}{{\mathrm{\Delta x}}^2}+\frac{1}{{\mathrm{\Delta y}}^2}+\frac{1}{{\mathrm{\Delta z}}^2})}^{-1}---\left(1.4\right)$

需要指出的是，当同时求解麦克斯韦方程和热传导方程时，必须选择最小时间步的值。 

在计算热方程时，需要边界上的温度场值，同样可以采用类似电磁场中Mur吸收边界的方法，利用边界附近网格的温度场去构建边界面上的温度场值。为了得到固体导热偏微分方程的唯一解，必须附加边界条件和初始条件，统称为定解条件，这里介绍初始条件和三类边界条件。 

初始条件已知是过程开始时物体整个区域中所具有的温度为已知值，用公式表示为 

T|_t＝0＝T₀                (1.5) 

式中，T₀[℃]为一已知常数，表示物体初温是均匀的。 

第一类边界条件是指物体边界上的温度函数为已知，用公式表示为 

T|_Γ＝T_w                 (1.6) 

式中，Γ为物体边界，其方向是逆时针方向；T_w[℃]为已知壁面的温度（常数）。 

第二类边界条件是指物体边界上的热流密度为已知，用公式表示为 

$-k\frac{\partial T}{\partial n}{|}_{\Gamma }=q_0---\left(1.7\right)$

式中，q₀[W/m²]为已知热流密度（常数）；当q₀＝0时，称为绝热边界条件。绝热边界条件本身不能唯一确定温度场，需要别的边界条件或初始条件配合下才能确定温度场。 

第三类边界条件是指与物体相接触的环境温度T_f和热对流系数h为已知，用公式表示为 

$k\frac{\partial T}{\partial n}{|}_{\Gamma }=-h(T-T_f)---\left(1.8\right)$

式中，h与T_f可以是常数，也可以是某种随时间和位置而变化的函数，如果h和T_f不是常数，则在数值计算中经常分段取其平均值作为常数。 

二、时域谱元法的简单介绍 

谱元法结合了有限元法和谱方法，求解过程与有限元法类似，其主要步骤是：(1)首先把计算的区域分成许多子域（单元），每个子域由若干节点组成；(2)在每个子域上把近似解表示成截断的切比雪夫正交多项式或勒让德多项式展开；(3)用伽辽金法对波动方程进行测试从而形成矩阵方程；(4)求解方程。它的基本思想是在每一个单元上Gauss-Lobatto-Legendre(GLL)基函数离散。在一维标准参考单元ξ∈[-1,1]中，N阶GLL基函数为 

${{\phi }_j}^{\left(N\right)}\left(\xi \right)=\frac{-1}{N(N+1)L_N\left({\xi }_j\right)}\frac{(1-{\xi }^2)L_N^{\prime }\left(\xi \right)}{(\xi -{\xi }_j)}---\left(2.1\right)$

其中，j＝0,1,…N，L_N(ξ)是N阶勒让德多项式，L'_N(ξ)是其导数。GLL积分点的选取为方程式(1-ξ_j²)L_N’(ξ_j)＝0的根，它们均定义在式ξ∈[-1,1]内。因而温度T可以表示成： 

利用变分原理和Galerkin法，将热传导方程在离散空间中转换成其残量形式： 

R_i＝∫_ΩN_irdΩ＝0           (2.3) 

其中N_i为加权函数，残量r为 

$r=-\frac{k_t}{{\rho }_m{c}_m}{▿}^{2}T+\frac{\partial T}{\partial t}-\frac{P_d}{{\rho }_m{c}_m}---\left(2.4\right)$

将(2.2)式和(2.4)式代入(2.3)式并利用矢量恒等式及 后得到 

在不同媒质交界面上，温度场满足一类特定的关系。以第三类边界条件为例，温度场满足以下关系： 

$\frac{\partial T}{\partial n}=\frac{h}{k_t}(T_{\mathrm{air}}-T_{\mathrm{surf}})---\left(2.6\right)$

其中，h为对流传热系数，一般取h＝10[W/(m²·℃)]；n为边界面的单位法向矢量，T_air为媒质外部温度即环境温度，T_surf为媒质的表面温度。将(2.6)式代入(2.5)式得到， 

其紧凑格式为： 

$\left[S\right]T_j+\left[T\right]\frac{{\partial T}_j}{\partial t}+\left[R\right]T_j=\left[\mathrm{Rq}\right]+\left[F\right]---\left(2.8\right)$

其中： 

${\left[S\right]}_{\mathrm{ij}}=\frac{k_t}{{\rho }_m{c}_m}\int \int \int ▿N_i·▿N_j\mathrm{dv}---\left(2.9\right)$

[T]_ij＝∫∫∫N_i·N_jdv            (2.10) 

${\left[R\right]}_{\mathrm{ij}}=\frac{h}{{\rho }_m{c}_m}{\int \int N}_i·N_j\mathrm{ds}---\left(2.11\right)$

${\left[\mathrm{Rq}\right]}_i=\frac{h{T}_{\mathrm{air}}}{{\rho }_m{c}_m}\int \int N_i\mathrm{ds}---\left(2.12\right)$

${\left[F\right]}_i=\frac{P_d}{{\rho }_m{c}_m}\int \int \int N_i\mathrm{dv}---\left(2.13\right)$

采用前向差分格式得： 

[T]Tⁿ⁺¹＝([T]-Δt([S]+[R]))Tⁿ+Δt[Rq]+Δt[F]       (2.14) 

其中，T为质量矩阵，S为刚度矩阵，R边界积分矩阵，F为源矢量矩阵，R_q由边界条件所引起的边界源矢量矩阵。由于单元基函数的正交特性，我们所得到的质量矩阵为对角阵或块对角阵，这样我们可利用块对角矩阵求逆的方法事先求出质量矩阵的逆，使得被求解方程为显式方程，降低了计算量，提高了计算效率。这是谱元法的一大重要特性。 

三、热源项分析 

微带线的辐射损耗相对于介质损耗和导体损耗比较小决定了辐射损耗可以忽略不计。在SETD公式推导中，当只记入介质损耗时，此时用等效电导率σ′代替σ： 

σ'＝σ+ωε″                   (3.1) 

其中，σ为电导率，ε″为复介电常数的虚部；ε、σ'分别对应待求电场值处的复介电常数实部和等效电导率。在求解电磁问题时，假定介质中的场分布经过激励源加载后的短暂瞬时过程后，达到一个相对稳定的分布。一旦得到电场值分布，据此就可以计算出媒质损耗电磁能以及转化到热模型中的热能。每个单元处的单位体积损耗能转化为热能功率的计算式如下： 

$P={\sigma }^{\prime }\overline{E^2}={\mathrm{\omega ϵ}}_0{{ϵ}_r}^{\prime \prime }\overline{E^2}---\left(3.2\right)$

上式中，是电场E的均方根值，ε_r″是复相对介电系数的虚部。 

在SETD公式推导中，只记入导体损耗时，由麦克斯韦第二方程：求解得到介质基片表面磁场分布，再由麦克斯韦第一方程：求解得到电流然而电磁波在良导体中的衰减很快，在传播很短的一段距离之后就几乎衰减完了。因此，良导体中的电磁波局限于导体表面附近的区域，即导体的趋肤效应。导体的趋肤深度非常小，以至于可以认为电流仅分布于导体表面很薄的一层内，也即介质基片表面一层。这与恒定电流或低频电流均匀分布于导体横截面上的情况不同，在高频时，导体的实际载流截面减小了，因而导体的高频电阻大于直流或低频电阻，它与电导率和趋肤深度有关，可表示成即厚度为δ的导体每平方米的电阻，称为导体的表面电阻率，简称表面电阻。最后由公式得出导体损耗。 

四、电磁—热耦合的基本理论 

数值方法模拟微波加热问题的总体思想是：首先，在t_n时刻根据初始的电特性条件，输入一定的微波功率，计算出各点的微波耗散功率，在模拟仿真时，假设电磁损耗能全部转化为热能。然后，将得到的热能代入到热传导方程中，得到各点的温度。根据物体各点处不同的温度分布，利用介电参数随温度变化的关系式得到下一时刻t_n+1的介质材料的电特性参数并再次重复计算电磁场方程，计算耗散功率，如此反复循环，直到完成预定的加热时间。在热场和电磁场的求解过程中采用相同的剖分网格单元。 

以图1真分析的带通滤波器为例：带通滤波器的通带范围为12GHz-12.5GHz。介 质基板采用Rogers RT5880，板厚0.508mm，具体尺寸由于涉密不方便给出。 

在进行热模拟仿真时需要具体知道微带电路的各部分性能参数：电导率，相对介电常数等。一般情况下，这些性能参数都会受到随周围的环境温度，压强的影响，这些参数都比较难以精确的获得，这也是数值模拟过程和实际升温过程存在差异的重要原因。 

表1，表2列出了一般情况下Rogers RT5880介质基板的性能参数，及覆铜片的的性能参数： 

表1介质基板的性能参数 

表2铜的性能参数 

将周围环境温度设为20℃，由于微波电路本身接触空气，故将周围表面设为第三类边界条件即为热对流形式。在微波电路的输入端注入幅度为10MV/m的调制高斯脉冲，中心频率设为12.5GHz，脉冲相关量设为9。图2为脉冲通过微波电路0.03us时的瞬态热分布图，可以清楚的看到，能量聚集在输入端口，这是由于调制高斯脉冲的频谱范围覆盖比较广，只有处在通带范围内的一部分频谱分量通过微波电路，到达输出端口。