马驰， 李江涛， 张伯昱， 严泽鑫， 孙子涵

(西安交通大学 电力设备电气绝缘国家重点实验室，陕西 西安 710049)

0 引 言

在对电力系统电磁暂态过程(electromagnetic transient process,EMTP)进行分析时，往往会涉及埋地电缆高频下的串联阻抗矩阵、并联导纳矩阵计算。现有大部分的电磁暂态分析程序，如电磁暂态程序(electromagnetic transient programs - the alternative transients program, EMTP-ATP)、电力系统计算机辅助设计软件(power systems computer aided design, PSCAD)等，其对电缆阻抗、导纳参数的计算通常使用准静态电磁场条件[1]。在低频条件下，土壤中的位移电流远小于传导电流，准静态电磁场条件不会引起较大误差；但随着频率升高到0.1MHz左右，土壤中位移电流大小的数量级与传导电流可比，此时传统的解析算法精度开始下降。此外，传统的解析算法在计算电缆阻抗、导纳参数时往往只考虑集肤效应，却无法将邻近效应带来的影响考虑在内[2]。

针对埋地电缆的参数计算，许多学者对此进行了研究。Wedephol等提出了同轴电缆阻抗、导纳矩阵参数的解析算法[3]，Ametani在此基础上将电缆参数的计算扩展到含缆芯、屏蔽、铠装三导体模型，并将算法应用于EMTP-ATP的CABLE CONSTANTS子程序中[4]。但是，上述方法都是建立在准静态电磁场中，忽略了位移电流及邻近效应，使得高频下参数的准确性受到影响；Ametani等采用对导体分区的方法计算含邻近效应的电缆参数矩阵[5-6]，然而导体分区的过程十分繁琐，且会带来大量的时间成本，在实际工程中难以应用。文献[7]利用表面导纳算子对埋地电缆的串联阻抗矩阵进行计算，但是其计算过程工作量大、模型适应性差，难以适用于实际工程计算中。

在对埋地电缆的串联阻抗及并联导纳矩阵进行计算时，尤其是在高频条件下，大地返回阻抗的计算至关重要。Pollaczek[8]第一个推导出埋地导线大地返回自、互阻抗的精确公式。Pollaczek公式中的无穷积分项不易计算，Wedephol[3]用截断无穷级数的方法对Pollaczek公式进行计算。Srivallipuranandan[9]以Romberg外推法为基础对Pollaczek公式中的无穷积分项直接进行数值积分，但因被积函数的剧烈震荡特性使得直接数值积分并不容易。Ametani[4]使用Carson公式近似计算大地返回阻抗。上述计算大地返回阻抗的解析公式的问题在于，随着频率的增加，积分项需要更长的时间才能收敛，导致截断误差；当频率接近1 MHz时，积分中的某些项会出现剧烈震荡的现象。Petrache[10]提出了一种简单的接地阻抗对数近似法，该法未考虑电缆的埋深，这可能会导致结果的误差。

本文采用有限元方法求解波动方程，对埋在有损土壤中任意截面的多导体电缆进行电磁模态分析。目前，已经有许多学者利用有限元法对电缆的参数进行计算[11-14]，但是这些方法都是建立在准静态电场条件下，其对位移电流的忽略可能造成高频下计算结果的误差。本文将电缆二维截面上电磁场的分布以及各种模式下传播常数的求解问题转化为对波动方程特征值的求解。在求得二维截面上电磁场的分布后，对电场及磁场进行积分，即可得到电缆中导体的电压、电流大小，从而通过进一步计算得到埋地电缆的串联阻抗及并联导纳矩阵。

1 电磁模态分析理论

对埋地电缆进行电磁模态分析，由于相域中电缆各导体存在电磁耦合，无法直接计算单一的传播常数，因此需要利用相模变换将电缆导体的电压、电流关系从相域转化为模域[15]，从而实现解耦。解耦后的电缆会存在不同的工作模式，对不同工作模式下电缆参数的计算需要至于时变电磁场中，以考虑高频下土壤中位移电流带来的影响。

1.1 相模变换理论

对含有n个导体的电缆系统，其存在n个独立的模式[16]。

首先研究单根电缆，以图1所示的双导体同轴电缆(导体为缆芯、护层)为例，假设电缆内的正弦电磁波沿z轴方向传播，其时域电波方程如下：

(1)

(2)

图1 双导体同轴电缆截面图Fig.1 Cross section of a two-conductor coaxial cable

方程中U=[U1U2]T表示电缆中导体的电压向量，I=[I1I2]T表示电缆中导体的电流向量，Z表示两导体的串联阻抗矩阵，Y表示两导体的并联导纳矩阵。由式(1)、式(2)可以推导出电缆内电磁波沿z轴方向传播的二维波动方程：

(3)

因为矩阵Z和矩阵Y都是对称阵，故式(3)中的矩阵YZ为矩阵ZY的转置，二者具有相同的特征值。现假设存在可逆矩阵S、Q，将相域中的电压电流向量U、I转换为模域中的电压电流向量Um、Im，即

(4)

若以下关系成立，则矩阵S、Q分别实现了对矩阵ZY、YZ的相似对角化：

(5)

式中：λ为矩阵ZY、YZ的特征值构成的相似对角阵，S、Q为特征值对应特征向量构成的矩阵。

由以上分析可知，通过对矩阵ZY、YZ的相似对角化处理，即可相域中的电压、电流向量U、I转换为模域中的电压电流向量Um、Im，两个导体组成的电缆系统即对应于模域中的两个模式。因为矩阵ZY变换为了对角阵，模域中的各模态电压、电流向量间不存在相互耦合关系，第k(k=1,2)个模式的传播常数γk即对应于矩阵λ中第k个特征值λk的平方根值

(6)

式中：αk对应于第k个模式的衰减常数；而βk对应于第k个模式的相位常数。

随着频率逐渐升高，在高频下，变换矩阵S中元素的虚部变得很小，通常可以忽略不计。通过分析相域电流I和模域电流Im之间的关系，可以得到第k个模式下对应的相域电流的分布。以双导体电缆为例，其变换矩阵Q为2阶方阵，将矩阵Q进行列分块，则式(4)中的第二个公式可以变形为

(7)

令Im1=1，Im2=0，可以得到第1个模式下电缆中导体的相域电流分布；同理令Im1=0，Im2=1，可以得到第2个模式下电缆中导体的相域电流分布。通过上述分析可知，变换矩阵Q的每一列都对应于相应模式下的相域电流分布。

对于双导体同轴电缆而言，其在高频下一般存在两种工作模式：接地模式和同轴模式。接地模式表示电流不经过缆芯，直接由护层和大地构成回路;同轴模式表示电流从缆芯流入，从护层返回，从而构成通路。两种模式对应的相域电流分布如图2所示。

图2 双导体同轴电缆不同模式下相域电流分布Fig.2 Phase domain current distribution of a two-conductor coaxial cable under different modes

相模变换理论可以推广至多根埋地电缆的模型中。假设现有三根双导体同轴电缆呈水平排布，则整个电缆系统中的导体总数为6，在高频下存在6个不同的模式：3种同轴模式、两种内护层模式和一种接地模式。记三根导体分别为A、B、C，则这6种模式对应的相域电流分布如图3所示。

图3 3根双导体同轴电缆不同模式下相域电流分布Fig.3 Phase domain current distribution of three two-conductor coaxial cables under different modes

1.2 不同模式电磁场的计算

在时变电磁场中，通常使用麦克斯韦方程组对电磁波的传播特性进行研究。假设所要分析的空间中无自由电荷，则麦克斯韦方程组可以表示如下：

(8)

麦克斯韦方程组中的第一个方程为全定律电流，全电流可以变形为

(9)

故只要导体满足似稳条件ωε/σ<<1，即可忽略位移电流。一般金属导体的电导率σ≥107S·m，故在电磁暂态计算关注的频率范围内，金属内的位移电流均可忽略。对于土壤，其电导率很低，在高频下土壤中位移电流的大小往往不可忽略。对于某些导电性很差的土壤，从50 kHz开始位移电流的大小就不能忽略，属于电磁暂态模拟关注的范围。

因此，在利用电磁模态分析对电缆参数进行计算时，可以对电缆中缆芯、护套等导体层中的位移电流进行忽略，仅仅计及土壤中的位移电流分布即可。

在三维坐标空间中，假设电缆中的电磁波沿着z轴方向传播，由全电流定律和电磁感应定律，可以得到如下方程：

(10)

(11)

γ=αz+jβ。

(12)

式中：kz表示面外波数；k0表示真空中的波数；α表示电磁波沿z轴方向传播时的衰减常数；β表示电磁波沿z轴方向传播时的相位常数。利用有限元法，即可将二维截面上电场、磁场的求解转化为对方程特征值问题的求解。图4给出了双导体同轴电缆在接地模式和同轴模式下的电场强度分布云图。

图4 双导体同轴电缆不同模式电场强度云图Fig.4 Cloud diagram of electric field intensity of two-conductor coaxial cable under different modes

在利用有限元法对电缆中的电磁场进行模拟时，需要设置一个范围很大的截断空间，在截断空间的边界上设置表面电流密度为0。截断空间半径的确定方法如下：

对于渗透进土壤的交流电流，其在土壤中的散流程度可以用大地渗透深度进行表示，即为

(13)

式中：μ0表示真空中的磁导率；σ表示土壤的电导率；ρearth表示土壤的电阻率。

在设置截断空间来模拟空气及大地时，至少应将截断半径设置为3～5倍的大地渗透深度。对电阻率为50 Ω·m的土壤来说，其在工频下的透入深度约为0.5 km，意味着在使用有限元法建立工频下的土壤模型时，其截断半径至少为1.5 km，形成的巨大网格将极大的增加计算时间。因此，虽然使用电磁模态分析可以计算低频下的电缆参数，但却会带来巨大的时间成本，且低频下已有解析公式的计算精度已经足够。

2 高频下电缆电气参数的计算

在高频下，利用有限元法对地下电缆进行电磁模态分析，可以得到各种模式下电缆的传播常数、电磁场分布。根据求得的传播常数，可以计算各模式下电磁波的衰减常数和波速；根据各模式下电磁场的分布，经过一系列变换，可以计算埋地电缆的串联阻抗矩阵及并联导纳矩阵。

2.1 各模式衰减常数、波速及相对波速的计算

建立埋地电缆模型后，可以根据有效折射率对想要求得的模式进行搜索，并在此过程中求得各模式下电磁波的传播常数。由式(12)可知，求解器求得的复传播常数的实部αz，即是电磁波在该模式下的衰减常数，表征电磁波在传播方向上每单位长度的振幅和能量衰减程度。根据电磁波理论，波速可以由下述公式求得：

(14)

式中：ω表示角频率；β表示相位常数(对应于传播常数的虚部)。

在电磁暂态分析和有关电缆的故障诊断、定位中，经常实际波速相对于真空波速的比值，即相对波速，可以表示为

(15)

式中：c表示真空中的光速；β表示真空中波的相位常数(也可理解为真空波数k0)，可以用关系式β=k0=ω/c进行求解。

2.2 电缆串联阻抗、并联导纳矩阵的计算

假设电缆中电磁波沿z轴方向传播，电缆中导体的电压向量U和电流向量I可以作如下改写：

(16)

(17)

将式(16)、式(17)代入时域电报方程(1)、方程(2)，可以得到以下关系：

(18)

(19)

对于不同的模式，电磁波具有不同的传播常数γ，可以列写不同模式下的式(18)、式(19)。对于双导体同轴电缆，两个模式下的传播常数已经求出，只要知道对应模式下导体的电压、电流值，即可以通过下述关系求出电缆的串联阻抗和并联导纳矩阵：

(20)

(21)

上文通过有限元法求解了不同模式下电缆内的电场分布和磁场分布。电缆的埋设图如图5所示。图中O点为电压参考点，不同参考点的选取可能对结果产生轻微影响。在高频条件下，导体的表面可以看做等位面，因此沿线l1和l2作电场强度的线积分，即可得到两个导体的电压值；沿面S1、S2作电流密度的面积分，即可得到两个导体的电流值：

(22)

(23)

图5 地下双导体同轴电缆截面图Fig.5 Cross section of an underground two-conductor coaxial cable

3 案例分析

以单根双导体同轴电缆为例进行分析：首先，在电磁暂态仿真软件EMTP-ATP中计算其利用传统解析法得出的串联阻抗和并联导纳矩阵；其次，利用电磁模态法计算频率f≥0.5 MHz时的串联阻抗和并联导纳矩阵，并将上述两者的计算结果进行对比；最后，在EMTP-ATP中建立JMarti模型进行时域仿真，并利用电磁模态法的计算结果替换f≥0.5 MHz时的输入参数文件，将传统解析法和电磁模态法得到的时域波形进行对比。

3.1 单根双导体同轴电缆

假设有一双导体同轴电缆，其结构如图5所示。电缆的几何半径设为：a=2 cm，b=3.75 cm，c=4 cm，d=5 cm；缆芯的电导率σc=9×107S/m，护层的电导率σS=4×107S/m；模型中导体、绝缘、空气以及土壤的相对磁导率全部设置为1；内绝缘的相对介电常数εr1=3，外绝缘的相对介电常数εr2=3，导体、空气及土壤的相对介电常数全部设置为1。如图5所示，空气的电导率σ0=0，土壤的电导率σg=0.02 S/m。电缆的埋深h分别设置为0.1、0.25和1 m，以便对比不同埋深下电缆的高频参数特性。

图6、图7分别给出了不同埋深下，衰减常数以及相对波速随频率变化的关系图。图中，模式1表示单根双导体同轴电缆的接地模式，模式2表示单根双导体同轴电缆的同轴模式。

图6 单根埋地电缆衰减常数随频率变化关系图Fig.6 Variation of the attenuation constant with frequency of a two-conductor coaxial cable

由图6、图7可以发现，对于单根双导体同轴电缆而言，相较于同轴模式，接地模式拥有更高的衰减程度，其相对波速却更低。Ametani[15]分析了不同模式下电磁波传播路径的阻抗和电容的影响，对其进行了解释：对接地模式而言，电流从护层进入，流经大地后返回护层，回路中包含数值较大的大地返回阻抗，且外绝缘与大地间的电容较大，使得电磁波传播过程中有较多电流泄入大地，导致电磁波以很快的速度衰减且波速较慢；对同轴模式而言，其与大地间的电磁联系不强，故其具有更低的衰减常数和更高的波速。

图7 单根埋地电缆相对波速随频率变化关系图Fig.7 Variation of the relative velocity with frequency of a two-conductor coaxial cable

图8给出了不同埋深下，电缆的串联阻抗中的元素随频率变化的关系图。由于缆芯与护层的互阻抗、护层的自阻抗与缆芯的自阻抗值十分相似，故图8中只给出了元素R11、L11随频率变化的关系图。

图8 单根埋地电缆串联阻抗矩阵随频率变化关系图Fig.8 Variation of the series impedance matrix with frequency of a two-conductor coaxial cable

由图8可以看出，当频率升至0.5 MHz时，传统解析法与电磁模态法的计算结果开始呈现明显差异。相较于传统的解析算法，利用电磁模态法计算得到的串联阻抗矩阵元素的电阻更大，电感更小；电缆的埋深越大，自电阻及互电阻越大，自电感和互电感越小。

图9给出了不同埋深下，电缆的并联导纳中的元素随频率变化的关系图。本文建立的电缆中不考虑半导体层的影响，故并联导纳矩阵中的并联电导部分可以忽略；对于考虑半导体层影响的电缆，则需考虑并联电导的影响。

图9 单根埋地电缆并联导纳矩阵随频率变化关系图Fig.9 Variation of the shunt admittance with frequency of a two-conductor coaxial cable

由图9可以看出，对于缆芯和护层间的电容，两种方法的计算结果完全相同；对于护层的自电容，由于其包含护层对地的电容，由于电磁模态法考虑了高频下的位移电流的影响，故当频率升至0.5 MHz时两种计算方法的结果开始呈现明显差异。

3.2 时域仿真

为了进一步说明电磁模态法与传统解析法对电磁暂态计算造成的影响，故需在时域中探究两种方法计算得到暂态波形的差异。相较于同轴模式，接地模式的电磁暂态过程与大地间的联系更为紧密，故选取接地模式作为时域暂态仿真的电路，如图10所示。

图10 接地模式的电路模型Fig.10 Circuit model of grounding mode

图中，电缆的截面形状如图5所示，电缆的电气参数设置与3.1中的双导体同轴电缆的设置相同。电缆长度为1 km，埋深为0.25 m。

首先在EMTP-ATP中采用频率相关的JMarti模型，对线路分别施加1.2/50微秒的雷电波、1.2/5微秒的雷电波以0.1/1微秒的VFTO，得到输出波形；接着用电磁模态法计算得到的串联阻抗和并联导纳矩阵替换原EMTP-ATP输入文件中频率大于等于500 kHz的部分，得到输出波形，并将结果与EMTP-ATP中直接使用解析法得到的结果进行对比。图11给出了不同激励下使用传统解析法和电磁模态法输出的暂态过电压波形。

由图11可知，在接地模式下，当主频率较高、波前较陡的过电压波侵入电缆时，相较于传统的解析算法，将电磁模态法得到的串联阻抗和并联导纳矩阵应用于暂态计算时，其波形的传播速度更快，二者在极短的暂态时间内存在一些差异。当土壤电阻率越高、大地的相对介电常数越大时，大地中位移电流的影响越大，两种方法得到的暂态波形的差异越明显。

4 结 论

在高频条件下(f≥0.5 MHz)，大地中的位移电流会对埋地电缆串联阻抗及并联导纳矩阵的计算产生较大影响，原有准静态电磁场下的解析算法对位移电流的忽略可能会使结果产生较大偏差，具体的偏差与电缆的结构、埋深、土壤电阻率等因素有关；在时域下的仿真说明了高频下串联阻抗及并联导纳矩阵的差异可能会导致电磁暂态计算结果的差异。位移电流的影响直接决定了时域中输出波形的差异大小，若大地电阻率越高、介电常数越大，则位移电流的影响越大，传统解析法和电磁模态法的时域仿真结果的差异越大。