很多初学者第一次看到 Dommel 算法时,会被“等效电导”“历史电流源”“伴随模型”这些词吓住。尤其是动态元件突然被改写成一个导纳并联一个电流源,看起来像是电路仿真里凭空长出来的一套新理论。
这篇文章我们先不从这个式子开始。我们先回到最普通、最熟悉的一阶微分方程:
dtdy=f(t)或者写成积分器的样子:
y(t)=∫f(t)dt只要把这个积分器用梯形法则离散,再把式子里的变量换成电路里的电压、电流,就会自然走到 Dommel 形式。所谓历史电流源,并不是新物理概念,而是普通微分方程离散以后留下的上一拍信息。
前两篇文章的铺垫已经很充足了,这里也是快快的回顾一下,熟悉的读者可以跳过。
从最熟悉的积分器开始
设有一个最简单的积分问题:
dtdy=x(t)也就是:
y(t)=∫x(t)dt在连续域里,它对应的传递函数就是:
X(s)Y(s)=s1这就是仿真软件里最常见的积分算子。输入是x,输出是y,连续时间下的含义非常朴素:输出是输入对时间的累积。
现在进入离散仿真。设仿真步长为T,第k步的输入、输出分别为xk和yk。用梯形法则近似一步内的积分面积:
yk−yk−1=2T(xk+xk−1)整理一下:
yk=yk−1+2T(xk+xk−1)这就是梯形积分器最核心的递推式。它只做了一件事:当前输出等于上一拍输出,再加上当前输入和上一拍输入形成的梯形面积。
从 S 到 Z:积分器的离散形式
对上面的递推式做z变换:
Y(z)−z−1Y(z)=2TX(z)+2Tz−1X(z)所以:
Y(z)(1−z−1)=2T(1+z−1)X(z)得到离散积分器:
X(z)Y(z)=2T1−z−11+z−1如果连续域积分器是1/s,那么梯形法则给出的对应关系就是:
s1≈2T1−z−11+z−1也可以反过来写成常见的s到z替换:
s≈T21+z−11−z−1到这里为止,我们完全没有讲电路,也没有讲 Dommel。我们只是把一个普通积分器离散化了。
但这一步已经足够了。因为电感、电容的本质,就是把电压或电流放进这个积分器里。
电感只是一个电流积分器
先把电感看成一个二端元件。规定电感两端电压为u,支路电流为i,电流参考方向和电压参考方向保持一致。
电感可以看成电流积分器
电感的连续时间关系是:
u=Ldtdi整理成积分器形式:
dtdi=L1u也就是说,电感电流i是电感两端电压u/L的积分:
i(t)=i(0)+∫L1u(t)dt这和刚才的y′=x完全一样,只是变量名换了:
| 普通积分器 | 电感 |
|---|
| y | i |
| x | u/L |
| y′=x | di/dt=u/L |
套用梯形积分器:
ik−ik−1=2LT(uk+uk−1)整理为当前步未知量加历史项:
ik=2LTuk+(ik−1+2LTuk−1)令:
GL=2LT上一拍折算出来的历史项记为:
IL,his=ik−1+GLuk−1于是得到:
ik=GLuk+IL,his这个式子已经开始有 Dommel 的味道了。但先别急着给它命名。我们只看它的结构:
- 第一项含有当前未知电压uk。
- 历史项只含有上一拍的ik−1和uk−1,当前步开始前已经知道。
这就是电感的 Dommel 形式。如果不了解数值积分算法,初学时看到这个式子会觉得很突然;但从积分器角度看,它其实只是把前面的梯形积分公式重新整理了一下。
普通积分器的梯形递推式是:
yk=yk−1+2T(xk+xk−1)对电感来说,只不过是把:
yk→ik,xk→Luk代进去:
ik=ik−1+2T(Luk+Luk−1)再把含有当前电压uk的项放到前面,把上一拍已经知道的量放到括号里:
ik=2LTuk+(ik−1+2LTuk−1)所以所谓GL和IL,his,不是额外发明出来的新东西,而只是给这两部分重新起了名字:
GL=2LT,IL,his=ik−1+GLuk−1也就是说,ik就是积分器里的yk,uk/L就是积分器里的xk。Dommel 形式只是把“当前未知量”和“上一拍历史量”分开放,方便后面进入电路求解器。
电容只是一个电压积分器
电容可以看成电压积分器电容的连续时间关系是:
iC=CdtduC整理成积分器形式:
dtduC=C1iC也就是说,电容电压uC是电容电流iC/C的积分:
uC(t)=uC(0)+∫C1iC(t)dt它和刚才的普通积分器也是同一个形式:
| 普通积分器 | 电容 |
|---|
| y | uC |
| x | iC/C |
| y′=x | duC/dt=iC/C |
所以第一步仍然只是套用同一个梯形积分器。普通积分器是:
yk=yk−1+2T(xk+xk−1)对电容来说,只是把:
yk→uk,xk→Cik代进去:
uk−uk−1=2CT(ik+ik−1)注意,这一步得到的是“当前电压怎么由电流积分得到”。它和普通积分器完全同源,没有任何额外技巧。
但是在电路求解里,我们更常希望支路关系写成“当前电流由当前电压和历史量决定”的形式。于是把上式重新整理,把当前电流ik单独解出来:
ik=T2Cuk−T2Cuk−1−ik−1令:
GC=T2C历史项为:
IC,his=−GCuk−1−ik−1于是:
ik=GCuk+IC,his所以电容这里的uk就是积分器里的yk,ik/C就是积分器里的xk。它看起来比电感多绕了一步,只是因为电容的天然积分结果是电压,而我们最后又把式子反过来整理成电流形式。
换句话说,电容的 Dommel 形式也不是新物理模型。它只是把梯形积分递推式里的当前项和历史项分开命名:
- GCuk:当前电压对应的当前电流项。
- IC,his:上一拍电压和上一拍电流折算出来的历史项。
R 和 G:没有积分器,所以没有历史项
电阻和电导没有微分方程,也没有积分器。
电阻:
uk=Rik所以:
ik=R1uk令:
GR=R1得到:
ik=GRuk如果元件本来就是电导G:
ik=Guk所以 R 和 G 没有历史源:
Ihis=0这也很好理解:电阻不储能,不记忆上一拍状态。
Dommel 形式:把物理元件翻译成网络求解器的语言
现在把前面的结果放在一起,就可以给它一个统一名字了:
| 元件 | 连续关系 | 积分器视角 | 梯形离散后整理形式 | 等效导纳 | 历史项 |
|---|
| R | u=Ri | 无积分器 | ik=GRuk | GR=1/R | 0 |
| G | i=Gu | 无积分器 | ik=Guk | G | 0 |
| L | u=Ldi/dt | i积分u/L | ik=GLuk+IL,his | GL=T/(2L) | ik−1+GLuk−1 |
| C | i=Cdu/dt | u积分i/C | ik=GCuk+IC,his | GC=2C/T | −GCuk−1−ik−1 |
所谓 Dommel 形式,就是把离散递推式整理成:
ik=Gequk+Ihis这个式子有两层含义。第一层是物理形式:元件仍然是原来的电阻、电导、电感和电容。电阻和电导是瞬时代数关系;电感储存磁场能量,电流来自电压积分;电容储存电场能量,电压来自电流积分。Dommel 形式没有改变这些物理关系。
第二层是计算形式:在某一个离散时刻k,梯形积分已经把微分关系变成了代数递推关系。把当前未知量和上一拍已知量分开以后,动态元件在求解器眼里就临时变成了一个等效导纳并联一个历史电流源:
- Gequk含有当前未知电压,要进入当前步的节点导纳矩阵。
- Ihis只含上一拍的电压、电流或内部状态,在当前步开始前已经算好,可以进入右端电流注入。
左边就是 Dommel 形式最常用的物理图像:一个等效导纳并联一个历史电流源,也就是诺顿等效形式。右边是同一个端口关系的戴维南等效形式:历史电压源串联等效电阻。两种画法对外端口等效,只是节点导纳法装配时更偏爱左边这种并联形式。
所以Ihis不是神秘电流源,也不是凭空加进电路里的新物理元件。它就是梯形积分器递推式里那些“上一拍”的东西,只是被挪到了电路方程右边。
笔者初学的时候感觉,“历史电流源”这个名字其实有一点误导。因为历史项里不一定只有电流,也可能有上一拍电压,甚至有程序内部保存的其他状态变量。比如上面的电感历史项里有ik−1和uk−1,电容历史项里也同时出现了uk−1和ik−1。从数学角度看,算法并不关心它们原来叫电压还是电流;只要它们来自上一拍,在当前步开始前已经知道,就都属于历史量。
之所以最后把它画成“历史电流源”,是因为我们选择了诺顿等效的写法。观察左边的并联等效图,Gequk是等效导纳支路从上节点流向下节点的电流,Ihis是电流源支路从上节点流向下节点的电流。两条并联支路的电流相加,就得到这个端口从节点k流向节点m的总支路电流:
ik=Gequk+Ihis如果把同一个端口关系改画成戴维南等效,就可以写成历史电压源串联等效电阻的形式:
uk=Reqik+Uhis其中:
Req=Geq1,Uhis=GeqIhis这里把戴维南图中的历史电压源定义为上端为正、下端为负,因此它在端口电压uk中表现为加号。如果后续画图时把电压源极性反过来,只需要把公式里的+Uhis改成−Uhis即可。
所以“历史电流源”更准确地说,是“由历史变量折算出来、在诺顿等效中表现为电流源的那一项”。它的来源是历史变量集合,表现形式才是电流源。
电路仿真喜欢这种写法,是因为节点导纳法天然要解当前节点电压。只要支路电流都能写成Gequk+Ihis,程序就可以把Geq装配进导纳矩阵,把Ihis装配进右端向量。动态元件的物理记忆仍然保留在历史项里,但当前步求解器面对的是一组普通代数方程。这就是“伴随模型”的工程意义:不是改写物理本质,而是把微分方程的递推关系翻译成网络求解器容易处理的语言。
例子一:R+L 串联支路
先看 R 和 L 串联。支路电压为:
u=Ri+Ldtdi这里动态变量仍然是支路电流i。把电感电压写成:
uL=u−Ri再代入电感的梯形积分式:
ik−ik−1=2LT[(uk−Rik)+(uk−1−Rik−1)]整理当前项:
(R+T2L)ik=uk+uk−1+(T2L−R)ik−1于是:
ik=GRLuk+IRL,his其中:
GRL=R+T2L1IRL,his=GRL[uk−1+(T2L−R)ik−1]这个例子说明,串联电阻并没有改变套路。它只是让等效导纳的分母里多了一个R,历史项里多了一些上一拍电压、电流。
例子二:R+C 串联支路
R 和 C 串联时:
u=Ri+uC电容电压满足梯形积分:
uC,k=uC,k−1+2CT(ik+ik−1)代入支路电压:
uk=Rik+uC,k−1+2CT(ik+ik−1)整理:
(R+2CT)ik=uk−uC,k−1−2CTik−1于是:
ik=GRCuk+IRC,his其中:
GRC=R+2CT1IRC,his=−GRC(uC,k−1+2CTik−1)如果程序不想显式保存uC,也可以用上一拍的支路关系uC,k−1=uk−1−Rik−1替换。保存哪个状态,是实现选择;离散化思想没有变。
例子三:R+L+C 串联支路
R、L、C 串联时:
u=Ri+uL+uC电感的梯形形式可以写成:
uL,k=T2L(ik−ik−1)−uL,k−1电容的梯形形式可以写成:
uC,k=uC,k−1+2CT(ik+ik−1)代入总电压:
uk=Rik+T2L(ik−ik−1)−uL,k−1+uC,k−1+2CT(ik+ik−1)把当前电流ik收集起来:
(R+T2L+2CT)ik=uk−[uC,k−1−uL,k−1+(2CT−T2L)ik−1]于是:
ik=GRLCuk+IRLC,his其中:
GRLC=R+T2L+2CT1IRLC,his=−GRLC[uC,k−1−uL,k−1+(2CT−T2L)ik−1]这个式子看起来长了一点,但它仍然只是在做同一件事:先用梯形法则离散积分器,再把当前未知量和历史已知量分开。
不过 RLC 串联支路比前面的 RL、RC 更容易让人误解。它有两个储能元件,因此程序内部必须保存足够的状态量。外部端口电压只告诉我们:
u=uR+uL+uC它并不能单独告诉我们电感上分了多少电压、电容上分了多少电压。上面的历史项里出现了uL,k−1和uC,k−1,正是在提醒这一点:把支路最后压缩成一端口 Dommel 形式,并不等于把内部状态消掉了。程序仍然需要知道上一拍电感和电容的能量分布,或者保存与它等价的一组状态变量。
如果没有额外说明,仿真程序通常会采用零初值:电感初始电流为0,电容初始电压为0。在这个串联支路的写法里,第一步计算时就相当于把历史量取成初始量,例如i0、uC,0,以及这个推导里显式出现的uL,0。若用户给了非零初值,就不能再把这些历史量默认为零,而要把它们写进第一拍的历史项里。更常见的物理初值是“电感电流初值”和“电容电压初值”;至于是否显式保存uL,还是用其他等价状态把它算出来,是程序实现上的选择。
所以这里真正要注意的是:一端口 Dommel 形式只是对外部求解器说“这一拍我可以写成一个等效导纳加一个已知源”。它没有说支路内部只剩一个数。对于含有多个储能元件的支路,历史项背后可能是一组状态,而不是一个简单的上一拍端口电压或端口电流。
实际程序里,RLC 串联支路也可以拆成多个伴随模型并引入内部节点,不一定非要压缩成一个一端口公式。本文这样写,是为了让读者看到:复杂一点的支路也没有逃离普通微分方程数值解。
小结
这篇文章真正想拆掉的是 Dommel 形式的神秘外衣。
- 从y′=x出发,梯形法则给出离散积分器。
- 这个离散积分器对应s≈T21+z−11−z−1。
- 电感是“电压积分成电流”,电容是“电流积分成电压”。
- 梯形离散以后,当前未知量进入等效导纳,上一拍信息进入历史项。
- 最后整理出的Ihis,名字里虽然有“电流”,但数学上本质是“来自上一拍的历史量”或“当前步已经显式已知的量”。在诺顿等效里,它被折算成一个电流源来表示;换成别的等效形式,它也可以表现成历史电压源或其他等价写法。
所以,如果已经理解了微分方程的数值通用解法,再理解 Dommel 方法并不难。它只是把数学上的离散递推,翻译成了电路网络求解器喜欢的语言。