1.   引　言

非金属导电材料具备质量轻和导电性好的特点, 广泛应用于工业和军事等领域作为设备制造材料, 例如飞机蒙皮使用的复合材料就是典型的非金属导电材料. 在飞机等复杂目标中, 大量的传输线铺设在非金属导电材料上, 实现不同设备之间的数据通信. 当这类目标处于复杂电磁环境中时, 空间强电磁干扰源将通过目标内的传输线耦合产生强电流信号, 该信号必然流入传输线端接电路, 对电路上的敏感元件造成干扰或破坏. 因此, 实现对非金属导电材料上传输线的电磁耦合建模, 是分析飞机等电大尺寸目标电磁环境效应的前提.

由传输线理论可知, 在空间电磁场的激励下, 传输线与其邻近的非金属导电材料之间将形成电流回路. 为了便于分析, 将非金属导电材料等效为具有一定介电常数和电导率的有限尺寸大小和厚度的有耗介质层. 全波算法是模拟有耗介质层上传输线电磁耦合的最直接方法, 其中时域有限差分(FDTD)方法^[1-4]是应用非常广泛的时域算法. 但是, 采用FDTD方法进行模拟, 需要对有耗介质层和传输线精细结构直接建模, 剖分所需网格量较大, 势必造成计算效率低. 因此, 国内外学者基于传输线方程理论, 开展了大量的数值算法研究, 实现传输线电磁耦合的高效计算. 这类算法的核心思想是在避免对传输线直接建模的前提下, 获得与全波算法具有相同精度的计算结果. 其中, 主流的算法有3类: Beam-Liu-Tesche (BLT)方程^[5-11]、FDTD-SPICE (simulation program with integrated circuit emphasis)算法^[12-17]和FDTD-TL (transmission line)算法^[18-20]. BLT方程是将传输线看作管道, 激励源和负载看作节点, 构建节点电压电流经管道传播的关系矩阵, 进而求解矩阵以获得负载上的电压和电流响应. 但是, BLT方程是一种频域方法, 当入射波为宽频带信号时, 计算效率不高. FDTD-SPICE算法是一种时域方法, 其首先使用传输线理论建立传输线的SPICE等效电路模型^[21], 然后通过FDTD方法模拟传输线的激励场并引入SPICE软件作为激励源, 最后, 采用SPICE软件仿真得到传输线端接负载上的电压和电流响应. 但是, 该算法在推导SPICE等效电路模型时需要用到大量的理论推导, 且算法中负载瞬态响应与传输线激励场需分开计算, 因此计算效率不高. FDTD-TL算法是本文的前期研究成果, 首先使用传输线方程构建空间电磁场作用传输线的电磁耦合模型, 然后通过FDTD方法模拟传输线周围空间的电磁场分布, 并在FDTD的每个时间步上引入到传输线方程作为等效分布源项, 最后采用FDTD的中心差分格式离散传输线方程^[22], 迭代求解得到传输线和端接负载上的瞬态响应. 相较于其他算法, 该算法实现了空间电磁场辐射与传输线瞬态响应的同步计算. 然而, 这类算法针对的研究对象均是理想地和实际地面上的传输线, 而对于有耗介质层上传输线的研究还未开展, 其原因是现有的传输线方程不适用于有耗介质层上传输线的电磁耦合分析.

因此, 本文首先基于经典传输线方程, 推导了适用于有耗介质层上传输线电磁耦合分析的修正传输线方程. 然后, 结合FDTD-TL算法和插值技术, 提出了一种高效的时域混合算法, 实现有耗介质层上多导体传输线瞬态响应的快速计算. 最后, 通过相应计算实例的数值模拟, 验证了时域混合算法的正确性和高效性.

2.   时域混合算法理论

2.1   有耗介质层上传输线电磁耦合分析的修正传输线方程推导

时域混合算法的核心是建立适用于有耗介质层上传输线电磁耦合分析的传输线方程. 经典的传输线方程包含两个方程, 称之为第一和第二电报方程, 其均由Maxwell方程组推导得到. 下面以经典传输线方程为原型, 详细介绍修正传输线方程的推导过程.

建立第一电报方程时, 需在传输线与接地面零电势位置之间选取合适的闭合回路. 当接地面为金属地时, 零电势位于金属地表面. 接地面为实际地面时, 零电势为无穷远处. 适用于金属地和实际地面上传输线电磁耦合分析的传输线方程及其推导过程, 已在文献[23]中给出. 然而, 有耗介质层是具有一定介电常数和电导率的有限厚度的损耗介质, 已有的传输线方程无法满足电磁波作用有耗介质层上传输线的电磁耦合建模需求.

当空间电磁场进入有耗介质层之后, 随着深度的增加而不断衰减, 因此可以将有耗介质层下表面近似看成零电势位置. 按照图1中的闭合回路, 应用电磁感应定律$\nabla \times {{E}} = - {\rm{j}}\omega {\mu _0}{{H}}$并按照积分形式展开得到:

图 1  闭合回路和闭合曲面的选取

Fig. 1.  Selections of closed loop and surface.

下载: 全尺寸图片 幻灯片

其中, d和h分别表示有耗介质层的厚度与传输线架设的高度; $E_y$, $E_z$和$H_x$均为总场. 考虑到传输线为良导体, 导体表面切向电场为零, 而且电磁场到达有耗介质层下表面时衰减到近似为零, 即$E_y\left( {y, d + h} \right) = E_y\left( {y, 0} \right) = 0$. 此时, 将(1)式对y求导, 可得

y处传输线与有耗介质层下表面之间的电压$V\left( y \right)$可表示为

结合(3)式, 并将(2)式中的磁场$H{}_x$分解为入射磁场$H_x^{{\rm{inc}}}$和散射磁场$H_x^{{\rm{sca}}}$的叠加, 则(2)式可表示为

这里, 传输线的散射磁场$H_x^{{\rm{sca}}}$在闭合回路上产生的磁通$\varphi$与传输线上的电流I通过传输线的单位长度电感L建立联系, 表示为$\varphi = {\mu _0}\displaystyle\int_0^{d + h} {H_x^{{\rm{sca}}}{\rm{d}}z} =$$LI\left( y \right)$. 而$H_x^{{\rm{inc}}}$的积分项通过电磁感应定律表示为

基于此, (4)式可以表示为

其中,

(5)式即为修正以后的第一电报方程. 这里$E_y^{{\rm{inc}}}\left( y \right)$和$E_z^{{\rm{inc}}}\left( y \right)$为入射电场分量.

建立第二电报方程时, 需要在传输线上选取合适的闭合曲面, 如图1所示. 根据安培环路定律$\nabla \!\times\! {{H}} \!=\! {\rm{j}}\omega \varepsilon {{E}} \!+ {{J}}$, 方程两边同时进行闭合面积分, 并应用奥氏公式和旋度的散度恒等于零的结论可得

其中, E为总电场, H为总磁场, J为传输线上的电流密度. 将电场E分解为入射电场E^inc和散射电场E^sca的叠加, 可得

对(9)式右边的两项应用高斯定律可得${\rm{j}}\omega \displaystyle\mathop{{\int\!\!\!\!\!\int}\mkern-21mu \bigcirc}\nolimits_S {\varepsilon {{{E}}^{{\rm{inc}}}} \cdot {\rm{d}}{{S}}} = {\rm{0}}$和${\rm{j}}\omega \displaystyle\mathop{{\int\!\!\!\!\!\int}\mkern-21 mu \bigcirc}\nolimits_S {\varepsilon {{{E}}^{{\rm{sca}}}} \cdot {\rm{d}}{{S}}} = {\rm{j}}\omega q\left( y \right)$, 其中$q\left( y \right)$为传输线沿线的电荷密度.

电荷密度$q\left( y \right)$与散射电压${V^{{\rm{sca}}}}\left( y \right)$和传输线单位长度电容C之间满足关系: $q\left( y \right) = C{V^{{\rm{sca}}}}\left( y \right)$. 而散射电压可由总电压和入射电场表示为

将(10)式代入(9)式, 并将两边对y求导, 可得

(11)式即为修正以后的第二电报方程.

将(5)和(11)式由频域转换为时域, 即可获得适用于有耗介质层上传输线电磁耦合分析的时域传输线方程. 该方程可扩展到多导体传输线的情况, 表示为

式中, ${{V}}\left( {y, t} \right)$和${{I}}\left( {y, t} \right)$分别表示多导线上的电压和电流矢量; L和${{C}}$分别为多导线的单位长度电感和电容分布参数矩阵; ${{{V}}_{\rm{F}}}\left( {y, t} \right)$和${{{I}}_{\rm{F}}}\left( {y, t} \right)$分别为传输线方程的等效分布电压源和电流源项,

其中, ${{{E}}_{\rm{T}}}\left( {y, t} \right)$和${{{E}}_{\rm{L}}}\left( {y, t} \right)$由入射电场计算得到, 计算公式为

修正传输线方程的方程结构形式与经典传输线方程^[18]保持一致, 两者最大的区别在于等效分布源项的计算. 修正传输线方程的${{{E}}_{\rm{T}}}\left( {y, t} \right)$项表示有耗介质层下表面与传输线位置之间的垂直入射电场分量$E_z^{{\rm{inc}}}$的沿线积分, ${{{E}}_{\rm{L}}}\left( {y, t} \right)$为传输线位置的切向入射电场分量$e_y^{{\rm{ex}}}$与有耗介质层下表面的切向电场分量之差. 换言之, 修正传输线方程考虑了有耗介质层内部的电场对传输线电磁耦合的影响.

2.2   有耗介质层上多导线电磁耦合的FDTD求解

修正传输线方程建立了空间电磁场与有耗介质层上多导线之间的电磁耦合关系, 通过求解传输线方程, 即可获得多导线上的瞬态响应. 需要说明的是, 传输线方程中的传输线单位长度分布参数和等效分布源项均为未知量, 其计算精度决定了传输线方程的建模准确度. 因此, 在求解传输线方程之前, 需要准确计算有耗介质层上多导体传输线的单位长度分布参数和等效分布源项.

2.2.1   多导线单位长度分布参数的计算

多导线单位长度电感参数矩阵L可由经验公式${L_{ii}} \!= \!{\mu _0}\ln ( {2{h_i}/{r_i}} )/(2{\text{π}})$和${L_{ij}} \!= \!{\mu _0}\ln ({1 \!+\! 4{h_i}{h_j}/d_{ij}^2})\!/$$(4{\text{π}})$计算得到, 其中i和j分别表示第i根和j根传输线, ${L_{ii}}$和${L_{ij}}$分别表示第i根传输线的自电感以及与第j根传输线的互电感; ${h_i}$, ${h_j}$和${d_{ij}}$分别表示第i根和第j根传输线的高度以及两根传输线之间的距离, 如图2所示. 电容参数矩阵C由公式${{C}} = {\mu _0}{\varepsilon _0}{{{L}}^{ - 1}}$计算得到.

图 2  多导体传输线的横截面几何结构

Fig. 2.  Cross section geometry of multi-conductor transmission lines.

下载: 全尺寸图片 幻灯片

2.2.2   多导体传输线等效分布源项的计算

多导体传输线等效分布源项由传输线周围的空间电磁场获得, 而空间电磁场由FDTD方法计算得到^[24,25]. 由(16)和(17)式可以发现, 等效分布源项只与入射电场有关, 而与多导线的散射电场无关, 这是因为: 一方面, 导线的横截面较小, 散射场较弱; 另一方面, 根据镜像原理, 有耗介质层对入射波的反射会抵消导线的部分散射场. 因此, 采用FDTD方法模拟多导线周围空间电磁场时, 只需对有耗介质层进行网格剖分, 而无需对多导线直接建模. 但是, 多导线的高度和间距可为任意值, 即多导线未必落在FDTD网格的棱边上. 因此, 多导线等效分布源项所需的电场分量需要采用插值技术由相邻棱边上的FDTD电场分量计算得到, 如图3所示.

图 3  多导线沿线和垂直电场分量的插值示意图

Fig. 3.  Interpolation schemes of the electric fields along and perpendicular to the multi-conductor transmission lines.

下载: 全尺寸图片 幻灯片

对于多导线沿线的电场分量${E'_y}$, 插值计算公式为

对于多导线的垂直电场分量${E'_z}$, 插值计算公式为

其中, $\alpha$和$\beta$表示传输线在FDTD网格中的比例因子.

将计算得到的传输线激励场在FDTD的每个时间步进上引入到传输线方程作为等效分布源项. 建立好传输线方程之后, 采用FDTD方法的中心差分格式进行离散, 获得传输线电压和电流的FDTD迭代求解公式, 从而求解得到多导线及其端接负载上的电压和电流响应. 具体的迭代求解公式可由文献[15]获得.

对于多导线端接负载上的电压V₀和V_N不满足中心差分格式, 需分别采用前向差分和后向差分进行离散. 假定多导线按照FDTD网格划分成N段, 如图4所示, 负载上的电压迭代公式表示为

其中, $\Delta y$和$\Delta t$分别表示FDTD差分离散传输线方程所需空间步长和时间步长; ${{{R}}_{\rm{S}}}$和${{{R}}_{\rm{L}}}$分别表示多导线始端和终端负载阻抗矩阵.

图 4  传输线的FDTD网格划分

Fig. 4.  FDTD grid division of transmission lines.

下载: 全尺寸图片 幻灯片

3.   数值仿真与分析

采用时域混合算法对有耗介质层上单导体传输线和多导体传输线的电磁耦合进行数值模拟, 并与商业电磁仿真软件CST的计算结果进行对比, 来验证算法的正确性和高效性.

算例1　有耗介质层上单导线的电磁耦合模型如图5所示, 有耗介质层大小为0.2 m × 0.4 m, 厚度为0.01 m, 相对介电常数为10, 电导率为20 S/m. 单导线长度为20 cm, 高度为1.9 cm, 端接负载分别为50和100 Ω. 入射波为高斯脉冲垂直照射单导线, 幅度为1000 V/m, 脉宽为2 ns. 为了保证计算精度, 时域混合算法选用的网格大小为5 mm. 在计算空间电磁场分布时, 选用各向异性介质完全匹配层(UPML)截断边界, 入射波距离多导线的高度为4个空间网格大小. 图6给出了时域混合算法与CST微波工作室计算得到的负载R₂上的电压响应对比曲线. 可以看出, 两种方法的计算结果振荡周期保持一致, 且幅值吻合度非常高. 表1列出了两种算法计算所需内存和时间的对比, 可以看出, 时域混合算法相较于CST, 节省了47%左右的计算时间, 是因为时域混合算法无需对单导线直接建模. 这里需要说明的是, CST软件虽然提供线缆工作室模拟传输线的电磁耦合, 但是只适用于接地面为金属体的情况.

图 5  有耗介质层上单导线的电磁耦合模型

Fig. 5.  Coupling model of single transmission line on the lossy dielectric layer.

下载: 全尺寸图片 幻灯片

图 6  负载R₂上的电压响应

Fig. 6.  Voltages on the load R₂ computed by the two methods

下载: 全尺寸图片 幻灯片

表 1  两种方法计算算例1时所需内存和时间对比

Table 1.  Memories and computation time needed by the two methods for the first example.

方法	内存/MB	计算时间/min
CST	63	19
时域混合算法	34	10

下载: 导出CSV 

| 显示表格

算例2　有耗介质层上多导体传输线的电磁耦合模型见图7, 有耗介质层的大小为0.4 m × 0.7 m, 厚度为0.01 m, 相对介电常数为10, 电导率为50 S/m. 5根导线平行放置在有耗介质层上, 长度为0.5 m, 高度为1.1 cm, 间距为4 mm, 半径为1 mm. 始端负载R₁—R₅均为50 Ω, 终端负载R₆—R₁₀均为100 Ω. 入射波类型和算法选用的网格大小与算例1的相同.

图 7  有耗介质层上多导体传输线的电磁耦合模型

Fig. 7.  Coupling model of multi-conductor transmission lines on the lossy dielectric layer.

下载: 全尺寸图片 幻灯片

首先, 入射波角度设置为θ = 180°, ϕ = 90°和α = 180°, 即垂直照射多导线. 采用时域混合算法与电磁仿真软件CST计算得到负载R₁和R₇上的电压响应对比曲线, 如图8所示. 可以看出, 两种算法的计算结果基本保持一致.

图 8  入射波垂直照射下的多导线端接负载的电压响应 (a)负载R₁上的电压; (b)负载R₇上的电压

Fig. 8.  Voltages on the terminal loads of multi-conductor transmission lines under the condition of ambient wave perpendicular to the multi-conductor transmission lines: (a) Voltages on R₁; (b) voltages on R₇.

下载: 全尺寸图片 幻灯片

然后, 考虑入射波斜照射的情况, 将入射角度设置为θ = 135°, ϕ = 45°和α = 180°, 同样采用两种方法计算得到负载R₁和R₇上的电压响应对比曲线, 如图9所示. 可以看出, 在入射波斜照射的情况下, 两种算法的计算结果仍能保证很好的吻合度. 另外, 表2列出了两种算法计算所需内存和时间的对比. 相较于算例1, 时域混合算法比CST节省了更多的计算时间, 是因为导线数量增加, 导致CST剖分所需网格量增多.

表 2  两种方法计算算例2时所需内存和时间对比

Table 2.  Memories and computation time needed by the two methods for the second example.

方法	内存/MB	计算时间/min
CST	217	64
时域混合算法	78	15

下载: 导出CSV 

| 显示表格

图 9  入射波斜照射下的多导线端接负载的电压响应　(a)负载R₁上的电压; (b)负载R₇上的电压

Fig. 9.  Voltages on the terminal loads of multi-conductor transmission lines under the condition of ambient wave oblique to the multi-conductor transmission lines: (a) Voltages on R₁; (b) voltages on R₇.

下载: 全尺寸图片 幻灯片

4.   结　论

经典传输线方程不适用于电磁波作用有耗介质层上传输线的电磁耦合问题. 因此, 基于经典传输线方程, 推导了适用于有耗介质层上多导体传输线电磁耦合分析的修正传输线方程. 然后, 结合前期研究的FDTD-TL算法和相应的插值技术, 提出了一种高效的时域混合算法, 实现修正传输线方程的快速求解, 获得多导体传输线及其端接负载上的瞬态响应. 该时域混合算法避免了对多导线精细结构的直接建模, 并实现了空间电磁场辐射与多导线瞬态响应的同步计算. 通过相应计算实例的数值模拟, 验证了时域混合算法能够与CST微波工作室的全波仿真保持相同的计算精度. 当使用相同网格大小时, CST软件所需内存和计算时间较时域混合算法有优势, 但由于在时域混合算法中精细的传输线结构无需剖分网格, 因此, 在获得相同精度的条件下, 时域混合算法相较于CST软件具有显著的优势.