从支路电压到节点电压法

上一篇文章里,我们一直把一个支路两端的电压写成uku_k

这里的下标kk表示离散时刻,不是节点编号。换句话说,uku_k的意思是“第kk个时间步的支路电压”。这一点很容易和节点电压法里的节点编号混在一起,所以从这篇文章开始,我们把两个概念分开写。

把支路两端电压先记成uu,是一种很适合入门的写法。读者可以先把注意力放在“一个元件或一个支路怎么离散”,不用一上来就面对完整电路网络里的节点编号、矩阵装配和右端向量符号约定。

但是实际电路仿真程序通常不会真的把每条支路的端口电压都当成独立未知量。更常见的做法是节点电压法:先选定参考地,再把每个非参考节点的电压作为未知量。支路两端电压不是新的未知量,而是两个节点电压之差。

也就是说,前面写的:

uu

到了网络求解器里,通常要换成:

uab=uaubu_{ab}=u_a-u_b

其中aabb是节点编号,uau_aubu_b是当前时刻这两个节点对参考地的节点电压。若要显式写离散时刻,可以写成:

uab(k)=ua(k)ub(k)u_{ab}^{(k)}=u_a^{(k)}-u_b^{(k)}

这里上标(k)(k)表示第kk个离散时刻。为了让公式更轻,后文默认都在同一个当前时刻讨论,于是省略时间上标,只写:

uab=uaubu_{ab}=u_a-u_b

从一个 R 支路开始

先看最简单的电阻支路。设电阻连接在节点aa和节点bb之间,支路电流参考方向从aa流向bb

R 支路连接两个节点
R 支路连接两个节点

电阻支路关系是:

i=G(uaub)i=G(u_a-u_b)

其中G=1/RG=1/R

这句话已经很接近节点电压法了,因为右边只剩节点电压。对节点aa来说,这条支路从节点aa流出去的电流是:

ia=G(uaub)i_a=G(u_a-u_b)

对节点bb来说,同一条支路从节点bb流出去的电流是反方向:

ib=G(ubua)i_b=G(u_b-u_a)

把这两行放在一起:

[iaib]=[GGGG][uaub]\begin{bmatrix}i_a\\i_b\end{bmatrix}=\begin{bmatrix}G&-G\\-G&G\end{bmatrix}\begin{bmatrix}u_a\\u_b\end{bmatrix}

这就是最普通的电阻支路对节点导纳矩阵的贡献。

它也解释了为什么电路程序喜欢节点电压法。只要每条支路电流都能写成节点电压的线性组合,程序就可以把每条支路贡献的2×22\times2小矩阵装配到全局矩阵里。支路连接哪两个节点,就把这四个数加到哪两个节点对应的行列上。

L 支路

电阻支路只有当前节点电压,没有历史项。接下来换成电感支路,就会看到 Dommel 形式真正进入节点电压法的样子。

上一篇文章已经推过电感的梯形离散形式:

i=GLu+IL,hisi=G_Lu+I_{L,\mathrm{his}}

其中:

GL=T2LG_L=\frac{T}{2L}IL,his=ik1+GLuk1I_{L,\mathrm{his}}=i_{k-1}+G_Lu_{k-1}

这里的k1k-1仍然表示上一拍离散时刻。现在把这条电感支路接在节点aa和节点bb之间,并规定支路电流参考方向从aa流向bb,支路电压参考方向也取:

u=uaubu=u_a-u_b

于是支路电流就是:

i=GL(uaub)+IL,hisi=G_L(u_a-u_b)+I_{L,\mathrm{his}}

现在要规定清楚正负号。本文采用下面这套约定:

  • 支路电流参考方向从节点aa流向节点bb
  • 对每个节点写方程时,从该节点流出的支路电流记为正
  • 所以对节点aa,这条支路电流是流出节点,贡献为+i+i
  • 对节点bb,同一条支路电流是流入节点,按“流出为正”的约定,贡献为i-i

因此:

ia=GL(uaub)+IL,hisi_a=G_L(u_a-u_b)+I_{L,\mathrm{his}}ib=i=GL(ubua)IL,hisi_b=-i=G_L(u_b-u_a)-I_{L,\mathrm{his}}

把两行放到一起,就得到节点电压法里最常见的形式:

[iaib]=[GLGLGLGL][uaub]+[IL,hisIL,his]\begin{bmatrix}i_a\\i_b\end{bmatrix}=\begin{bmatrix}G_L&-G_L\\-G_L&G_L\end{bmatrix}\begin{bmatrix}u_a\\u_b\end{bmatrix}+\begin{bmatrix}I_{L,\mathrm{his}}\\-I_{L,\mathrm{his}}\end{bmatrix}

这就是截图里那种结构。把电感换成任意已经整理好的 Dommel 支路,只要它能写成:

i=Gequ+Ihisi=G_{\mathrm{eq}}u+I_{\mathrm{his}}

并且支路方向仍然从aa指向bb,就有:

[iaib]=[GeqGeqGeqGeq][uaub]+[IhisIhis]\begin{bmatrix}i_a\\i_b\end{bmatrix}=\begin{bmatrix}G_{\mathrm{eq}}&-G_{\mathrm{eq}}\\-G_{\mathrm{eq}}&G_{\mathrm{eq}}\end{bmatrix}\begin{bmatrix}u_a\\u_b\end{bmatrix}+\begin{bmatrix}I_{\mathrm{his}}\\-I_{\mathrm{his}}\end{bmatrix}

左边的两个量ia,ibi_a,i_b,表示这条支路分别对节点aa、节点bb的“流出电流贡献”。如果程序采用“流入节点为正”或者把所有已知项移到方程右端,历史项向量的符号会跟着变,但物理方向和支路方程本身没有变。

R+L+C 串联支路

现在回到上一篇文章最后的 RLC 串联支路。

支路层面已经得到:

i=GRLCu+IRLC,hisi=G_{RLC}u+I_{RLC,\mathrm{his}}

其中:

GRLC=1R+2LT+T2CG_{RLC}=\frac{1}{R+\frac{2L}{T}+\frac{T}{2C}}IRLC,his=GRLC[uC,k1uL,k1+(T2C2LT)ik1]I_{RLC,\mathrm{his}}=-G_{RLC}\left[u_{C,k-1}-u_{L,k-1}+\left(\frac{T}{2C}-\frac{2L}{T}\right)i_{k-1}\right]

注意这里的k1k-1是离散时刻,不是节点编号。本文统一用a,ba,b表示节点编号,用ua,ubu_a,u_b表示节点电压,避免把“时间步下标”和“节点编号下标”混在一起。

如果 RLC 串联支路连接节点aa和节点bb,并规定支路电流从aa流向bb,那么:

u=uaubu=u_a-u_b

代入支路 Dommel 形式:

i=GRLC(uaub)+IRLC,hisi=G_{RLC}(u_a-u_b)+I_{RLC,\mathrm{his}}

于是它对两个节点的贡献就是:

[iaib]=[GRLCGRLCGRLCGRLC][uaub]+[IRLC,hisIRLC,his]\begin{bmatrix}i_a\\i_b\end{bmatrix}=\begin{bmatrix}G_{RLC}&-G_{RLC}\\-G_{RLC}&G_{RLC}\end{bmatrix}\begin{bmatrix}u_a\\u_b\end{bmatrix}+\begin{bmatrix}I_{RLC,\mathrm{his}}\\-I_{RLC,\mathrm{his}}\end{bmatrix}

这和电阻支路的矩阵结构完全一样。不同的只是:

  • 电阻支路的导纳GG是常数。
  • RLC 串联支路的等效导纳GRLCG_{RLC}来自梯形离散后的伴随模型。
  • RLC 串联支路多了历史项IRLC,hisI_{RLC,\mathrm{his}},它要进入当前步的右端向量。

所以从网络求解器的角度看,Dommel 形式的好处就在这里:动态支路虽然内部有电感、电容和历史状态,但在当前时刻装配矩阵时,它仍然像一个“等效导纳加已知源”的二端支路。

展开 RLC 的内部节点

上面的 RLC 写法,是把整个 R+L+C 串联支路压缩成一个二端支路,只看外部两个节点aabb。这一步对网络求解器很友好,因为它只需要给节点aa和节点bb装配一个2×22\times2小矩阵。

现在把难度加一点:不再把 RLC 压缩成一个整体,而是把 R、L、C 三个元件各自看成一条单独支路。这样串联支路内部会多出两个节点。设节点顺序从左到右为:

acdba\rightarrow c\rightarrow d\rightarrow b

其中:

  • R 支路连接节点aa和节点cc
  • L 支路连接节点cc和节点dd
  • C 支路连接节点dd和节点bb
RLC 串联支路的内部节点电压
RLC 串联支路展开后的内部节点电压

三个支路的电压分别是:

uR=uaucu_R=u_a-u_cuL=ucudu_L=u_c-u_duC=udubu_C=u_d-u_b

支路方向都取从左到右,于是三条支路可以写成:

iR=GR(uauc)i_R=G_R(u_a-u_c)iL=GL(ucud)+IL,hisi_L=G_L(u_c-u_d)+I_{L,\mathrm{his}}iC=GC(udub)+IC,hisi_C=G_C(u_d-u_b)+I_{C,\mathrm{his}}

这里GR=1/RG_R=1/R,而GLG_LGCG_C来自电感和电容各自的 Dommel 伴随模型。历史项的正负号仍然沿用前面的约定:支路方向从左到右,流出节点为正。

把这三条支路分别装配到节点a,c,d,ba,c,d,b上,就得到一个4×44\times4矩阵:

[iaicidib]=[GRGR00GRGR+GLGL00GLGL+GCGC00GCGC][uaucudub]+[0IL,hisIL,his+IC,hisIC,his]\begin{bmatrix}i_a\\i_c\\i_d\\i_b\end{bmatrix}=\begin{bmatrix}G_R&-G_R&0&0\\-G_R&G_R+G_L&-G_L&0\\0&-G_L&G_L+G_C&-G_C\\0&0&-G_C&G_C\end{bmatrix}\begin{bmatrix}u_a\\u_c\\u_d\\u_b\end{bmatrix}+\begin{bmatrix}0\\I_{L,\mathrm{his}}\\-I_{L,\mathrm{his}}+I_{C,\mathrm{his}}\\-I_{C,\mathrm{his}}\end{bmatrix}

这个大矩阵其实没有新东西。它只是把三条二端支路的小矩阵加在一起:

  • R 支路把GRG_R加到a,ca,c对应的四个位置。
  • L 支路把GLG_L加到c,dc,d对应的四个位置,同时把+IL,his+I_{L,\mathrm{his}}加到节点cc,把IL,his-I_{L,\mathrm{his}}加到节点dd
  • C 支路把GCG_C加到d,bd,b对应的四个位置,同时把+IC,his+I_{C,\mathrm{his}}加到节点dd,把IC,his-I_{C,\mathrm{his}}加到节点bb

所以,中间节点dd的历史项看起来是:

IL,his+IC,his-I_{L,\mathrm{his}}+I_{C,\mathrm{his}}

这不是新的物理源,而是因为节点dd同时连着 L 和 C:L 支路的电流方向是从ccdd,所以对dd来说历史电流是流入;C 支路方向是从ddbb,所以对dd来说历史电流是流出。按“流出为正”的约定,两项自然一负一正。

这也解释了为什么程序里常说“装配”。所谓装配,不是重新推一遍全电路方程,而是每条支路只关心自己两端节点,把自己的小矩阵和历史项加到全局矩阵、全局向量对应的位置上。

实战:正弦电压源驱动 RLC 支路

最后做一个很小的实战。为了不把问题一下子推进到改进节点电压法,这里先采用最简单的电路:理想正弦电压源一端接地,另一端接 RLC 串联支路。

接地电压源驱动 RLC 串联支路

这时右端节点是参考地:

ub=0u_b=0

左端节点由电压源强制给定:

ua=vs(t)u_a=v_s(t)

所以 RLC 支路两端电压不是未知量,而是已知量:

u=uaub=vs(t)u=u_a-u_b=v_s(t)

这一步很重要。这个例子不是在求一个含浮接理想电压源的完整节点矩阵,而是在已知端口电压的条件下,计算 RLC 串联支路的动态响应。也就是说,节点电压法里的“参考地”和“已知节点电压”已经先处理好了,下面只需要沿用前面推出来的 RLC Dommel 递推式。

取:

vs(t)=100sin(2π50t)v_s(t)=100\sin(2\pi 50t)

示例参数为:

R=10 Ω,L=20 mH,C=100 μF,T=20 μsR=10\ \Omega,\qquad L=20\ \mathrm{mH},\qquad C=100\ \mu\mathrm{F},\qquad T=20\ \mu\mathrm{s}

计算流程可以按下面理解:

  • 初始时刻默认uC(0)=0u_C(0)=0i(0)=0i(0)=0,因此第一拍的历史量来自零初值。
  • 每一个当前步先由电压源得到uk=vs(tk)u_k=v_s(t_k)
  • 再把上一拍保存的uL,k1u_{L,k-1}uC,k1u_{C,k-1}ik1i_{k-1}合成历史项,求出当前步iki_k
  • 最后用当前步iki_k更新uL,ku_{L,k}uC,ku_{C,k}。下一拍再计算历史项时,用的就是这组刚刚更新过的状态。

所以程序里看起来是在更新uLu_LuCu_C,本质上是在保存电感和电容的内部状态。IhisI_{\mathrm{his}}不是一个必须单独长期保存的神秘变量;它可以在每一拍开始时,由上一拍状态重新算出来。

Python 代码如下,完整脚本保存在同目录的 simulate-rlc-voltage-source.py

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
import numpy as np
import matplotlib.pyplot as plt

R = 10.0
L = 20e-3
C = 100e-6
f = 50.0
U_peak = 100.0
T = 20e-6
t_end = 0.12

t = np.arange(0.0, t_end + T, T)
u = U_peak * np.sin(2.0 * np.pi * f * t)
i = np.zeros_like(t)
u_L = np.zeros_like(t)
u_C = np.zeros_like(t)

G_rlc = 1.0 / (R + 2.0 * L / T + T / (2.0 * C))

for k in range(1, len(t)):
    history = u_C[k - 1] - u_L[k - 1] + (T / (2.0 * C) - 2.0 * L / T) * i[k - 1]
    i[k] = G_rlc * (u[k] - history)
    u_L[k] = 2.0 * L / T * (i[k] - i[k - 1]) - u_L[k - 1]
    u_C[k] = u_C[k - 1] + T / (2.0 * C) * (i[k] + i[k - 1])

fig, ax1 = plt.subplots(figsize=(8.6, 3.8))
ax1.plot(t * 1000.0, u, label="source voltage")
ax1.set_xlabel("time / ms")
ax1.set_ylabel("voltage / V")
ax1.grid(True)

ax2 = ax1.twinx()
ax2.plot(t * 1000.0, i, color="tab:red", label="branch current")
ax2.set_ylabel("current / A")

fig.tight_layout()
fig.savefig("rlc-voltage-source-result.svg", format="svg")

仿真结果如下:

正弦电压源驱动 RLC 串联支路的响应
正弦电压源驱动 RLC 串联支路的响应

左图中蓝色曲线是外加电压源,红色曲线是 RLC 串联支路电流;右图画出同一时刻 R、L、C 三个元件各自承担的电压。由于这里的电容、电感初值都取零,刚开始会有一段暂态;经过几个周波以后,响应逐渐进入 50 Hz 正弦稳态。这个例子虽然很小,但它已经把前面的推导串起来了:先由电压源给出当前步支路电压,再用历史项保存上一拍的电感、电容状态,最后递推出当前步支路电流。

读者的思考

很多资料讲 Dommel 形式,讲到:

i=Gv+Ihisi=Gv+I_{\mathrm{his}}

就停下来了。这个式子当然重要,但真正写程序时,很快会遇到一堆它没有直接回答的问题。

比如,刚刚的例子,求解的电路为什么一定要有参考地?如果一个网络完全浮起来,不接地,节点电压矩阵会发生什么?再比如,理想电压源如果不是一端接地,而是浮接在两个非参考节点之间,还能不能直接塞进普通节点导纳矩阵?

这篇文章先不深究这些问题。本文只做一件事:把单条支路的 Dommel 形式,接到节点电压差和节点矩阵装配上。至于参考地、浮接电压源、改进节点电压法、理想源约束方程、矩阵奇异这些问题,后面单独开文章再慢慢拆。


相关内容

Buy me a coffee~
RLC侠 支付宝支付宝
RLC侠 微信微信