SPICEによるケーブルの過渡応答の計算

業界標準の汎用回路解析プログラムSpice

  ftp://ic.eecs.berkeley.edu/pub/Spice3/sp3f4.kit.tar.Z
  ftp://ic.eecs.berkeley.edu/pub/Spice3/sp3f4.patch.tar.Z
を使ってケーブルの過渡応答を計算するための手法を解説します。

1. 集中定数回路で扱える場合

ケーブルの伝送特性計算で直面する最初の問題が

  1. ケーブルを集中定数回路で扱うことができるか
  2. 分布定数で扱わなければならないのか
の判定ですが、これは、ケーブルの長さが短いとき、つまり、
  電磁波がケーブルを往復する時間 < パルスの立上り時間       (1)
が成り立つときは、ケーブルを集中定数回路で扱うことができます。つまり、ケー ブル全体を1図のような L, C, R, G の1つのセグメントで表すことができます。


1図 ケーブルの集中定数モデル

L, C, R, G はそれぞれケーブル全体のインダクタンス(H)、キャパシタンス(F)、 抵抗(Ohm)、コンダクタンス(G)ですが、下記のような周波数特性があります。ωは 角周波数、δは絶縁材料の誘電損失角です。

  1. R は導体の渦電流により、高周波で sqrt(ω) に比例する
  2. G = ω*C*tan(δ)
  3. L は導体の内部インダクタンスが周波数の増加とともに 0 に近付く
  4. C は絶縁体によって周波数特性や温度特性の大きなものがある
このうち特に注意が必要なのは R と G で、どちらも大きな周波数特性があります が、G については誘電体損失の少ない絶縁材料なら他のパラメータに比べて無視で きるほど小さくなるのが普通で、通常は省略可能。例えば、EIA-232(V.28) 19200bps の伝送回路で 22 AWG の PVC 絶縁多芯ケーブルを長さ 10 m で使った回 路だと次のようになります。ドライバの出力インピーダンスは半導体メーカのデー タシートから計算するか、実測します。また、多芯ケーブルの構造、回路に対する 導体割り当て、送受端に接続されるドライバとレシーバによって C と L が変わり ますが、典型的な値を使いました。

V.28 CIRCUIT EXAMPLE
VIN 1 0 PULSE(-9 +9 1US 0.6US 0.6US 52.3MS 104.6MS)
RG 1 2 75
L1 2 3 3.5UH
R1 3 4 0.75
C1 4 0 1.6NF
R2 4 5 0.75
L2 5 6 3.5UH
RL 6 0 3K
.TRAN 100P 10U
.PLOT V(1) V(6)

ケーブル全体の L, R を半分に分けて、C を真中に入れたのは、少しでも分布定数 線路に近づけようという意図です。

この例では R1, R2 がドライバの出力インピーダンスやレシーバの入力抵抗に比べ て小さいため無視できて、例えば R1 と R2 を削除しても結果はほとんど変わりま せんから、R の周波数特性を考えずに済みますが、R を無視できないときは、次の 二つのアプローチがあります。

  1. 最高ないし主な周波数成分に対する値を使う
  2. 導体の周波数特性を近似する回路を付加る
最初の方法は、例えば
  f = 0.35 / Tr                   (2)

  ここに、f = 周波数 (Hz)
          Tr = パルスの立上り時間 (s)
の周波数に対する交流抵抗を R として使うもので、これ以下の周波数については 過大評価になりますが、「これ以上悪くならない」という評価はできます。もう 少し低い周波数に対する値を使えば、多少過小評価にはなりますが、かなり実際 に近い結果が得られます。

二つ目の方法については、後で説明します。

2. 分布定数回路で扱わなければならない場合

条件 (1) が成り立たないときは分布定数回路として扱うことになって、上記の回 路 (L, C, R, G のセグメント)を従属接続した2図のようなモデルが基本になりま す。


2図 ケーブルの分布定数モデル

この場合はケーブル全体をいくつのセグメントに分割するかが問題になりますが、 損失の少ない線路で分割(セグメントの数)が少なすぎると本来無限小セグメントな ら存在しないはずのリンギングが出ますし、かといって分割を増やすと計算時間や メモリ消費が急速に増加して計算不能になります。

無難なセグメント分割の目安は、次のとおりです。

  1. 損失の少ない線路 (R << ω*L)
    1. 1つのセグメントの遅延時間(sqrt(L*C)*l) <= パルスの立上り時間の 1/10
    2. 1つのセグメントの直列抵抗(R*l) <= 特性インピーダンスの 1/10
  2. 損失の多い線路 (R >> ω*L)
    1. 3から5セクションで充分
「損失の多い線路」というのは LSI の内部配線等の場合で、PCB(プリント配線板) やケーブルは「損失の少ない線路」になります。次の典型的な線路定数で R/ω*L の値を比べてみるとよくわかるでしょう。
Case    L0       C0      R0       Z0    T0     k   R0/(ω*L)(b)
      (nH/cm)  (pF/cm) (Ohm/cm) (Ohm) (ps/cm) 
--------------------------------------------------------------
PCB     5       1       0.050(a) 70.7 70.7    4.5   0.0023
MCM(d)  5       1       5        70.0 70.0    4.5   0.23
Chip    2.2     2     500        32.9 65.8    3.9  52.5
Coax(c) 2.6     0.97    0.0001   52.0 50.6    2.3   0.000005
--------------------------------------------------------------
(a) ω=2*π*f=2*π*0.35/Tr=4.4e9, Tr=0.5ns
(b) 4-mil width, 1-oz Cu
(c) RG-8/U
(d) MultiChip Module

例えば、ポリエチレン絶縁の同軸ケーブルの遅延時間は 3e8/sqrt(2.3) = 5.0nS/m ですから、立ち上がり 2 nS のパルスで使うとすると、4 cm 程度のセグメントに しなければなりません。IC の内部配線やプリント基板のパターン程度ですと楽で すが、100 m とか数 100 m のケーブルになると非常に大きな計算になって、 SPICEを使うのは実用的ではありません。

計算量を減らすには素子数を減らすのが一番ですが、ケーブルの場合はG が無視できるのが普通で、多くの場合、下記のようなT-sectionが使われます。


3図 T-section

例えば、RG-8/U同軸ケーブル 10 m で立ち上がり 25 ns、周期 200 ns の 方形波を伝送する場合だと次のようになります。

COAXIAL CABLE EXAMPLE
* RG-8/U, 1m
.SUBCKT XRG8U1M 10 30
L1 10 15 0.132UH
R1 15 20 0.086
C  20 0  97PF
R2 20 25 0.086
L2 25 30 0.132UH
.ENDS
* RG-8/U, 10m
.SUBCKT XRG8U10M 10 20
X10M0 10 11 XRG8U1M
X10M1 11 12 XRG8U1M
X10M2 12 13 XRG8U1M
X10M3 13 14 XRG8U1M
X10M4 14 15 XRG8U1M
X10M5 15 16 XRG8U1M
X10M6 16 17 XRG8U1M
X10M7 17 18 XRG8U1M
X10M8 18 19 XRG8U1M
X10M9 19 20 XRG8U1M
.ENDS

VIN   1  0  PULSE(0 1 0 25NS 25NS 100NS 200NS)
RS    1  10 52
XLINE 10 90 XRG8U10M
RT    90 0  52

.TRAN 1N 250N
.PLOT V(10) v(90)
.END
少しリップルが出ますが、特性インピーダンスで終端された線路ですから、本来こ のリップルはないはずで、この原因は有限分割とSPICEの計算方法の組合せ に起因するものです。分割を細かくしても完全になくすことはできません。

なお、上記の例では、T セクションを素直にSUBCIRCUIT化していますが、 SPICET(Lossless Transmission Line)モデルを使って、下記のよ うに書くこともできます。

.SUBCKT XRG8U1M 10 30
TL1 10 0 20 0 Z0=52 TD=5.06N
R1  20 30 0.172
.ENDS
また、SPICE3以降で導入されたLTRA(Lossy Transmission Line) モデルを使えば、次のようにもっと簡単になります。
COAXIAL CABLE EXAMPLE
* RG-8/U, 100 m
VIN   1  0  PULSE(0 1 0 25NS 25NS 100NS 200NS)
RS    1  10 52
O1    10 0 90 0 RG8U
RT    90 0  52
.MODEL RG8U LTRA R=0.0018 L=2.64NH C=0.97PF G=0 LEN=10
.TRAN 1N 250N
.PLOT V(10) v(90)
.END

先に述べたように、このモデルには、導体の渦電流による R と L の周波数特性が 組み込まれていませんから、シミュレートできるのは、

  1. 時間遅れ
  2. 反射
  3. 直流的な減衰
だけです。つまり減衰と位相歪みが問題にならず、反射だけが問題になるような回 路ということになりますが、これは単距離の高速伝送が典型的なケースで、 「SCSI ケーブル」や「GP-IB ケーブル」などがこの条件になります。

例えば、実際に市販されていた特性インピーダンスが 50 Ohm 程度の安い Single-EndedSCSIケーブルで何が起きるか調べて見ましょう。直流的な減衰 は無視できますから、SPICET(Transmission)モデルという無損失 線路を使います。

全長 6 m のSCSIケーブルの一端でバスを駆動し、その隣 0.5 m と他端に 接続されたデバイスの入力波形を調べると、次のようになります。

SCSI CABLE EXAMPLE
* driver
VIN 1 0 PULSE(0 +5 10NS 15NS 15NS 300NS 600NS)
RS 1 2 10K
* 3=collector, 2=base, 0=emitter
Q1 3 2 0 BJT
* terminator
R1 3 10 220
R2 3 0  330
* low cost SCSI cable
T1 3 0 5 0 Z0=53 TD=3.2NS
T2 5 0 9 0 Z0=53 TD=34.3NS
* terminator
R3 9 10 220
R4 9 0  330
* SCSI power
V1 10 0  4.8V
.MODEL BJT NPN
.IC V(3)=2.9 V(5)=2.9 V(9)=2.9
.TRAN 5N 600N
.PLOT V(3) V(5) V(9)
SCSIバスでは、 が論理の判定基準ですから、かなり遅いデバイス以外は使いものにならないことが わかります。

同期転送でどうなるかは、信号源を下記のように指定するとわかりますが、とても 使用に耐える状態ではありません。

VIN 1 0 PULSE(0 +5 10NS 15NS 15NS 90NS 180NS)

反射の原因はケーブルの特性インピーダンスと終端抵抗のミスマッチですから、終端 抵抗のインピーダンスを下げれば改善されるのは明らかで、アクティブ・ターミネー タの 110 Ohm と交換すると、次のように多少改善されます。

SCSI CABLE EXAMPLE
* driver
VIN 1 0 PULSE(0 +5 10NS 15NS 15NS 90NS 180NS)
*VIN 1 0 PULSE(0 +5 10NS 15NS 15NS 300NS 600NS)
RS 1 2 10K
* 3=collector, 2=base, 0=emitter
Q1 3 2 0 BJT
* Active terminator
R1 3 10 110
* line
T1 3 0 5 0 Z0=53 TD=3.2NS
T2 5 0 9 0 Z0=53 TD=34.3NS
* Active Terminator
R2 9 10 110
V1 10 0  2.85V
.MODEL BJT NPN
.IC V(3)=2.85 V(5)=2.85 V(9)=2.85
.TRAN 5N 600N
.PLOT V(3) V(5) V(9)

SCSI-1では「望ましい」としていた 90 Ohm 以上の特性インピーダンスが SCSI-2では必須になりましたが、90 Ohm になると下記のようにテブナン 終端でもほとんど問題がなくなります。

SCSI CABLE EXAMPLE
* driver
VIN 1 0 PULSE(0 +5 10NS 15NS 15NS 300NS 600NS)
RS 1 2 10K
* 3=collector, 2=base, 0=emitter
Q1 3 2 0 BJT
* terminator
R1 3 10 220
R2 3 0  330
* low cost SCSI cable
T1 3 0 5 0 Z0=90 TD=3.2NS
T2 5 0 9 0 Z0=90 TD=34.3NS
* terminator
R3 9 10 220
R4 9 0  330
V1 10 0  4.8V
.MODEL BJT NPN
.IC V(3)=2.85 V(5)=2.85 V(9)=2.85
.TRAN 5N 600N
.PLOT V(3) V(5) V(9)

もちろん、コストは上昇しますが、

T1 3 0 5 0 Z0=132 TD=3.2NS
T2 5 0 9 0 Z0=132 TD=34.3NS
と 132 Ohm の特性インピーダンスを持ったSCSIケーブルを使えば無反射に なります。

3. 導体抵抗とインダクタンスの周波数特性

前記のモデルはすべて R と L を定数として扱っていますが、実際には大きな周波数 特性があって、パルス波形に次のような大きな影響を与えます。

  1. R は sqrt(ω) に比例するため波形が鈍る
  2. L の内部インダクタンスが 1/sqrt(ω) に比例するため位相歪みを生ずる
L の変化は極めてわずかで、ネットワークアナライザでは測定できませんが、波形に 影響するのはその微分的変化ですから、これを無視すると正しい結果が得られません。

R の周波数特性をモデル化しなければならないのは高周波トランスを扱う場合など も同じですが、これらの周波数特性をSPICEで扱う場合は、R と L のネット ワーク回路で実現するのが普通で、いくつかのモデルが作れますが、 例えば下記の回路が使えます。


4図 導体の R と L の周波数特性を近似するための回路モデル

同軸ケーブルの中心導体の場合は下記のようになることが解析的に証明できます。 (注4)

  Rn = Rdc  (Ohm/m)
  Ln = μ/(jn^2*π)  (H/m)

  ここに、Rdc = 導体の直流抵抗 (Ohm/m)
          μ = 透磁率 (H/m)
          jn = 1次の Bessel 関数の 0 点

jn の最初の 20 個は次のようになります。

 n    jn
----------
 1 3.83171
 2 7.01559
 3 10.17347
 4 13.32369
 5 16.47063
 6 19.61586
 7 22.76008
 8 25.90367
 9 29.04683
10 32.48968
11 35.33231
12 38.47477
13 41.61709
14 44.75932
15 47.90146
16 51.04354
17 54.18555
18 57.32753
19 60.46946
20 63.61136

すべての定数が直流抵抗だけで決まり、導体の内部インダクタンスはサイズに関係 に一定という面白いモデルですが、これは内部インダクタンスがサイズに無関係な ことが原因です。


5図 R, L の周波数特性を考慮したケーブルモデル

4図のモデルは導体による損失と導体の内部インダクタンスを含んでいますから、 この後に外部インダクタンス、キャパシタンス、コンダクタンスを追加すれば、 R と L の周波数特性を考慮した、ケーブルモデルができます。

導体の R, L をモデル化したとき、同軸ケーブルの内部導体、つまり宇宙に単独で 置かれた無限長の円形断面導体というかなり特殊な状況を考えていましたが、実際 には、同軸ケーブルの外部導体の特性もほとんど同じですし、ツイステッド・ペア の場合も近接効果や周囲導体による渦電流の増加があまり大きくないため、通常の ほとんどのケーブルで、かなりの精度が得られます。

このモデルでケーブルが長くなると素子数が増えて計算量が爆発的に増えますが、 ケーブルが長いときは反射波も減衰しますので、反射が問題にならないことが多く、 ケーブル全体を4図の単一セグメントにしてしまうことができます。

また、5図の G が無視できる場合は右端の L, C は

  Td = 1/sqrt(L*C)

  ここに、Td = 遅延時間 (s)
          L = ケーブルの外部インダクタンス (H)
          C = ケーブルのキャパシタンス (F)
という無歪みの純粋な遅延を表しますので、この部分をSPICET モデルに置き換えた下図のモデルで計算量を減らすことができます。


6図 導体の R, L による歪みに無損失線路を組み合わせたケーブルモデル

この方法で長さ 1000 m の RG-8/U の過渡特性を計算すると次のようにな ります。

RG-8/U EXAMPLE
* RG-8/U, 1000 m
.SUBCKT XRG8U1000M 10 31
R0  10 11 9.84OHM
R1  11 12 9.84OHM
L1  11 12 27.2UH
R2  12 13 9.84OHM
L2  12 13 8.13UH
R3  13 14 9.84OHM
L3  13 14 3.86UH
R4  14 15 9.84OHM
L4  14 15 2.25UH
R5  15 16 9.84OHM
L5  15 16 1.47UH
R6  16 17 9.84OHM
L6  16 17 1.03UH
R7  17 18 9.84OHM
L7  17 18 0.07U2NH
R8  18 19 9.84OHM
L8  18 19 0.0596UH
R9  19 20 9.84OHM
L9  19 20 0.0474UH
R10 20 21 9.84OHM
L10 20 21 0.0379UH
R11 21 22 9.84OHM
L11 21 22 0.0320UH
R12 22 23 9.84OHM
L12 22 23 0.0270UH
R13 23 24 9.84OHM
L13 23 24 0.0231UH
R14 24 25 9.84OHM
L14 24 25 0.0200UH
R15 25 26 9.84OHM
L15 25 26 0.0174UH
R16 26 27 9.84OHM
L16 26 27 0.0154UH
R17 27 28 9.84OHM
L17 27 28 0.0136UH
R18 28 29 9.84OHM
L18 28 29 0.0122UH
R19 29 30 9.84OHM
L19 29 30 0.0109UH
R20 30 31 9.84OHM
L20 30 31 0.00988UH
.ENDS

VIN   1  0  PULSE(0 1 0 25NS 25NS 100NS 200NS)
RS    1  10 52
XLINE 10 80 XRG8U1000M
T1    80 0 90 0 Z0=52 TD=5.06U
RT    90 0  52

.TRAN 1N 6U 5U
.PLOT V(10) v(90)
.END

導体の内部インピーダンスを展開した L1, L2, L3, .. をどこまで計算に使うかの 判断は、 (2) 式の周波数の目安から求めたωを使って

  ω*Ln << Rdc を満足する項を捨てる
という基準によることになりますが、この導体モデルは素直でわかりやすい半面、 収束が悪いのが欠点で、立上りの早いパルスに適用する場合は非常にたくさんの 素子を必要とします。

市販の文献には出ていませんが、より収束の速いモデルを作ることも可能で、 いずれ改めてこの問題に触れようと思います。 (注5)

4. 表皮効果を R-C 線路に変換する方法

汎用性はありませんが、導体の表皮効果を R-C 線路に変換できる場合があります。

伝送線路の伝送線路の基本的な特性は、下記のようになりますが、

  Z0 = sqrt((R + j*ω*L)/(G + j*ω*C))
  γ = sqrt((R + j*ω*L)*(G + j*ω*C))
  ここに
	Z0 = 特性インピーダンス (Ohm)
	γ = 伝搬定数
	   = α*l + j*β*l
	α = 減衰定数 (neper/m)
	β = 位相定数 (rad/m)
	l = 線路の長さ (m)
	j = sqrt(-1)
周波数が十分高く、G << ω*C, R << ω*L が成り立つ場合は、

  Z0 〜 sqrt(L/C)
  γ = sqrt((R + j*ω*L)*(G + j*ω*C))
     〜 sqrt((j*ω*L + R)*j*ω*C)
     = j*ω*sqrt(L*C)*sqrt(1 + R/(j*ω*L))
     〜 j*ω*sqrt(L*C)*(1 + R/(j*ω*L)/2)
     = R/2*sqrt(C/L) + j*ω*sqrt(L*C)
     = R/2/Z0 + j*ω*sqrt(L*C)
になります。

周波数が高く、表皮の厚さ (skin depth) に比べて、導体の曲率が無視できる場合は、 R 〜 k*sqrt(j*ω) が成り立ちますから、

  γ = k*sqrt(ω)/2*sqrt(2)/Z0 + j*k*sqrt(ω)/2/sqrt(2)/Z0 + j*ω*sqrt(LC)

一方、G = L = 0 の場合は、

  Z0 = sqrt(R/(j*ω*C)
  γ = sqrt(j*ω*R*C)
     = ω*R*C/sqrt(2) + j*ω*R*C/sqrt(2)
これを高周波で表皮効果のある場合と比べると、 高周波で導体の表皮効果を考慮したケーブルの伝搬定数は、 皮効果を考慮したケーブルの伝搬定数 は、 R と C から構成される分布定数線路と同型であることがわかります。 つまり、
  R = 2*sqrt(2)*Z0*α/sqrt(ω)
とすれば、 下記のような等価回路が作れます。


7図 表皮効果を R-C 線路に変換したケーブルモデル

理想的な同軸ケーブルでは、

  α = Rs*(1/(2*π*a + 2*π*b)/(2*Z0)
  ここに
	Rs = sqrt(ω*μ/σ)
	σ = 導電率 (S/m)
	σ = 導電率 (S/m)
	a = 内部導体外径 (m)
	b = 外部導体内径 (m)
ですから、

  R = sqrt(μ/σ)*(1/(2*π*a + 2*π*b)
で等価回路の定数が決まります。

例えば、RG-174 同軸ケーブル 100 m のステップ応答なら、 次のように計算できます。

Step response of RG-174 (100m)
VIN 1 0 PULSE(-9 +9 1US 0.6US 0.6US 52.3MS 104.6MS)
RG 1 2 50
T1 2 0 3 0 Z0=50 TD=220N
U1 3 4 0 URCMOD L=100
RL 4 0 50
.MODEL URCMOD URC RPERL=2.37e-4 CPERL=100e-12
.TRAN 10P 10U
.PLOT V(1) V(4)
この方法は簡単ですが、 等価回路の特性インピーダンスが現実のケーブルとは全く違っていて、 減衰特性だけが互換ですから、 無反射の場合しか使えないことに注意してください。 例えば、この例で V(2) は正しい結果を与えません。 つまり、 表皮効果による損失を R-C 線路で表現しようとすると、 反射波の存在する環境では、 似ても似付かない波形になってしまいます。

5. 注と参考事項

注1 f = 0.35/Tr の根拠

多くの文献に出ている「常識」ですが、疑問を感じる方は 「パズル」のページ を御覧ください。

注2 SPICE以外の解法

SPICEでは実用的に解けない、5図のようなセググメントをたくさん従属 接続したモデルでも、自分でプログラムを書けば、例えば Backward-Euler 法で 能率良く解くことができます。

注3 RG-8/U 同軸ケーブルの概要

  特性インピーダンス = 52 Ohm
  キャパシタンス = 97 PF/m
  直流抵抗 = 0.10 Ohm/m
  減衰 = 19 dB/km
  遅延時間 = 5.06 ns/m
  中心導体 = 7/.724 mm
  絶縁体径 = 7.2 (ポリエチレン)
  外部導体 = 編組
  ジャケット径 = 10.3 mm (PVC)

注4 円筒導体の内部インピーダンスの等価回路を求める

4図の等価回路が正しいことは、例えば次の方法で証明できます。

まず、円筒座標 (r-z) で表した導体内部の Maxwell 方程式が次のようになること はよく知られています。(d/dr) は半径方向の偏微分演算子です。

  (1/r)(d/dr)(r*H(r)) = σ*E           (a)

  ここに、H(r) = 半径 r における磁界 (H/m)
          σ = 導電率 (S/m)
          σ = 導電率 (S/m)
          E = 長さ(z 軸)方向k電界 (V/m)\
一方、好運に恵まれた上で少し考えると、半径 r、時刻 t に於ける磁界の時間的 変化が
           ∞
  H(r,t) = Σ(An*exp(-Bn*t)*J1(jn*r/a)))    (b)
          n=1

  ここに、J1() = 1次の第1種 Bessel 関数
          J2() = 2次の第1種 Bessel 関数
          An = I/(J2(jn)*π*a*jn)
          Bn = (jn)^2/(μ*σ*a*a)
          jn = 1次の Bessel 関数の n 番目の 0 点
          μ = 透磁率 (H/m)
          σ = 導電率 (S/m)
          σ = 導電率 (S/m)
          a = 導体半径 (m)
          I = 導体の全電流 (A)
と Bessel 関数に展開すれば (a) が簡単に積分できることに気づきますので、 E のステップ応答が次のように求まります。
                      ∞
  E(a,t) = Rdc*I*(1 + Σ(exp(-Rdc/Ln*t))
                      n=1

  ここに、Rdc = 1/(σ*π*a*a)
          Ln = μ/(jn^2*π)
これが4図の等価回路になりますが、内部インダクタンスの合計を計算してみると、
        ∞
  Lin = Σ(Ln) = μ/(8*π)
        n=1

  ここに、Lin = 円筒導体の内部インダクタンス (H/m)
になって、電磁気学の結果と一致することを確認できます。

最後に、 SPICEの概要と参考書を付けておきます。

注5 - 4図より効率的な回路モデル

私が在職中(非公開で)使っていた方法は 伝送線路の SPICE モデルで解説しました。 (2015-08)

平林浩一, 1999-02