FreeFEMによるケーブル(電線)の電気設計 (3) - 直流インダクタンスの計算

1. ケーブルのインダクタンス

ケーブルのインダクタンスはキャパシタンスと対になって特性インピーダンスや 位相定数のほとんどを決める重要な特性の一つですが、同軸ケーブルや ツインリードのように、電磁界がともに存在する空間では、キャパシタンスと インダクタンスの間に比誘電率を媒介にした関係

  v = 1/sqrt(L*C) = c/sqrt(εs)
  ここに
	v = 電磁波の伝搬速度 (m/s)
	L = ケーブルのインダクタンス (H/m)
	C = ケーブルのキャパシタンス (F/m)
	c = 真空中の光速 (m/s)
	εs = 比誘電率
が成立ちますから、キャパシタンスが決まればインダクタンスも決まり、 単独のインダクタンス計算は必要ありません。

しかし、シールド付きのツイステッドペアのように、電界が存在しないシールドの外に 磁界が単独で存在する場合は、上記の関係が成立しませんから、 キャパシタンスとインダクタンスを別々に計算しないと、特性インピーダンスや 速度係数が求まりません。

インダクタンスの計算で面倒なのは、いかなる導体材料と絶縁材料を使うにせよ、 原理的に周波数特性があることです。 インダクタンスは回路空間に蓄積されるエネルギに比例しますが、 高周波になると表皮効果等の渦電流により導体内部の磁界が減少しますから、 インダクタンスも周波数の増加とともに減少します。

この渦電流の計算はかなり面倒ですが、回路のインダクタンスを便宜的に

  L = Le + Li

  ここに
	L = ケーブルのインダクタンス (H/m)
	Le = 外部インダクタンス (H/m)
	Li = 内部インダクタンス (H/m)
と分けて考えると、周波数特性の輪郭が掴みやすくなります。

ここで、「外部インダクタンス」(external inductance)というのは導体外部の空間に 蓄積される磁気エネルギに起因したインダクタンス、 「内部インダクタンス」(internal inductance)というのは導体内部の空間に蓄積される 磁気エネルギに起因したインダクタンスですから、 直流のインダクタンスは外部インダクタンスと内部インダクタンスの合計になりますし、 表皮効果により導体内部の電流と磁界のほとんどが存在しなくなる 高周波のインダクタンスはほぼ外部インダクタンスのみになります。 この中間の数 MHz までの交流では、インダクタンスはこの中間の値になります。

一般に 10 MHz 程度以上の周波数では、導体外部の空間の面積に比べて、電流が存在 する導体内部の面積は極めて少なくなりますから、インダクタンスの値はほとんど 外部インダクタンスと等しくなりますが、外部インダクタンスには近接効果による わずかな周波数特性を別にすれば大きな周波数による変化がありませんから、 伝送波形を遅延させる機能だけになり、伝送波形を歪ませることはありません。

ところが、内部インダクタンスには周波数特性があって、一般に周波数の平方根に 反比例して減少しますから、インダクタンス全体としては極わずかな 内部インダクタンスが伝送波形に大きな歪みを与えることになります。

2. 直流に於ける自己インダクタンスの計算

以下、インダクタンス計算の基本になる直流の(自己)インダクタンスを FreeFEMで計算する方法を解説します。

ケーブルの構造が長さ方向に一様だとすると、導体に流れる電流によって生ずる静 磁界のベクトルポテンシャル A はケーブルの軸方向の成分 (Az) だけを持 ち、導体の内部と外部のすべての領域で下記の偏微分方程式が成立ちます。(注1)

  (d/dx)((1/μ)(d/dx)(A) + (d/dy)((1/μ)*(d/dy)(A) = -j                  (2)

  ここに
	A = 磁界のベクトルポテンシャル (weber/m)
	μ = 透磁率 (H/m)
	  = 4e-7*π (真空、空気、銅、アルミなどのとき)
	j = 電流密度 (A/m^2)
FreeFEMの記法なら次のようになります。

  pdu(v) -laplace(v)/m = j                                             (3)

  ここに
	v = ベクトルポテンシャルの z 軸方向成分 (weber/m)
	m = 透磁率 (H/m)

導体を含めた空間に蓄積される磁気エネルギは

  B = rot(A)                                                           (4)

  ここに
	B = 磁束密度 (weber/m^2)
から磁束密度を計算して、
  W = (1/2)*∫(B*H)*dU
    = (1/2)*(1/μ)*∫abs(B)^2*dU                                       (5)

  ここに
	W = 磁界のエネルギ (joul/m^3)
	H = 磁界の強さ (ampare-turn/m)
	abs(x) = ベクトル x の絶対値
から計算するのが普通ですが、FreeFEMの記法では、
  w=0.5/m*(dx(A)*dx(A)+dy(A)*dy(A))                                    (6)
になります。

空間の全エネルギを求める場合は、微分演算を含まない下記の方法で計算すると 精度がよくなります。

  W = (1/2)*∫(A*j)*dU                                                 (7)

  ここに
	W = 磁界のエネルギ (joul/m^3)
以上の方法で空間に蓄積される磁気エネルギがわかれば、後はエネルギ、電流とイ ンダクタンスの関係

W = L*I^2/2 (8) ここに W = 磁界のエネルギ (joul) L = インダクタンス (H/m) I = 電流 (A) を使って、下記のように(自己)インダクタンス (L) を計算することができます。

  L = 2*W/I^2                                                          (9)

2.1 円柱導体の内部インダクタンス

最初に、最も簡単なケースとして、円柱導体の内部インダクタンスを考えます。 導体の全電流を 1 A として、導体表面にぴったりはりついた無限に薄い円筒導体に -1 A の電流を流すことにすれば、 導体表面でベクトルポテンシャルが 0 という境界条件が使えますので、 円柱導体の内部インダクタンスを計算するためのFreeFEMのプログラムは 次のようになります。

導体表面に張り付けた円筒導体は自分自身の内部に磁束を生じませんから、内部の 円柱導体に対しては何の影響もないことに注意してください。 エネルギ計算は (6) 式と (7) 式の両方で行っています。

/* inductance calculation - internal inductance of solid round wire */
a1:=0.4e-3;        /* inner conductor raius */
u1:=4e-7*pi;	/* permeability */

border(1,0,2*pi,100) { x:=a1*cos(t); y:=a1*sin(t) };
buildmesh(1000);

i1:=1;	/* total current */
j1:=i1/pi/a1/a1;	/* current density (total current == 1) */
/* vector potential */
solve(u) {
  onbdy(1) u = 0;
  pde(u) -laplace(u)/u1 = j1;
};

plot (u);

/* magnetic energy */
w=0.5*u*j1;
plot (w);

l1:=2*intt[w]/i1/i1;
printf("Li = %g (H/m, A-i method)\n", l1);

w=0.5/u1*(dx(u)*dx(u)+dy(u)*dy(u));
l1:=2*intt[w]/i1/i1;
printf("Li = %g (H/m, B-H method)\n", l1);
内部インダクタンスとして 4.99e-8 (H/m) が得られますが、これは解析的計算から 得られる正しい値 5.00e-8 (H/m) と比べて 2/1000 の精度です。(注3)

軸対称ですから解析領域を断面の一部に限定すれば、少ない分割で高い精度が得られ るばずで、例えば、断面の 1/8 に限定すれば次のようになります。


/* DC inductance calculation */
/* solid round conductor */
a1:=0.4e-3;	/* conductor raius */
u1:=4e-7*pi;	/* permeability */
t1:=pi/4;	/* computing 0 - pi/4 */

n1:=50;
n2:=floor(n1*a1*t1/a1);

border(1,0,n1+n2+n1,n1+n2+n1+1) {
 if(t<n1) { x:=a1*t/n1; y:= 0; ib:=1 }
 else if(t<=n1+n2) { t2:=(t-n1)*t1/n2; x:=a1*cos(t2); y:=a1*sin(t2); ib:=2 }
 else { t2:=a1-a1*(t-n1-n2)/n1; x:=t2*cos(t1); y:=t2*sin(t1); ib:=1 }
};

buildmesh(1000,1);

i1:=1;	/* total current */
j1:=i1/pi/a1/a1;	/* current density (total current == 1) */
j=j1*(region==0)-j1*(region==2);
/* vector potential */
solve(u) {
  onbdy(2) u = 0;
  pde(u) -laplace(u)/u1 = j;
};

plot (u);

/* magnetic energy */
w=0.5*u*j;
l1:=2*intt[w]/i1/i1*2*pi/t1;
printf("Li = %g (H/m, A-i method)\n", l1);

w=0.5/u1*(dx(u)*dx(u)+dy(u)*dy(u));
l1:=2*intt[w]/i1/i1*2*pi/t1;
printf("Li = %g (H/m, B-H method)\n", l1);
精度は 2/10000 になりました。

2.2 円筒導体の内部インダクタンス

円筒導体の場合も断面の 1/8 を解析することにすれば、次のようになります。単独 円筒とか同軸ケーブルの内部導体として使われる導体の場合は円柱導体と同じよう に、電流の帰路を導体表面に密着した円筒にすれば、導体表面のベクトルポテンシ ャルが 0 という簡単な境界条件になります。


/* DC inductance calculation */
/* tubular conductor */
a1:=2.45e-3;	/* inner radius */
a2:=2.80e-3;	/* outer radius */
u1:=4e-7*pi;	/* permeability */
t1:=pi/4;	/* computing 0 - pi/4 */

n1:=10;
n2:=floor(n1*a1*t1/(a2-a1));
n3:=floor(n1*a2*t1/(a2-a1));

border(1,0,n1+n3+n1+n2,n1+n3+n1+n2+1) {
 if(t<=n2) { t2:=t1-t1*t/n2; x:=a1*cos(t2); y:=a1*sin(t2); ib:=2 }
 else if(t<n2+n1) { x:=a1+(t-n2)*(a2-a1)/n1; y:= 0; ib:=1 }
 else if(t<=n2+n1+n3) {
  t2:=(t-n2-n1)*t1/n3; x:=a2*cos(t2); y:=a2*sin(t2); ib:=3
 }
 else {
  t2:=a2-(a2-a1)*(t-n2-n1-n3)/n1; x:=t2*cos(t1); y:=t2*sin(t1); ib:=1
 }
};

buildmesh(1000,1);

i1:=1;	/* total current */
j1:=i1/pi/(a2*a2-a1*a1);	/* current density (total current == 1) */
/* vector potential */
solve(u) {
  onbdy(3) u = 0;
/*  onbdy(2) u = 0;*/
  pde(u) -laplace(u)/u1 = j1;
};

plot (u);

/* magnetic energy */
w=0.5*u*j1*2*pi/t1;

l1:=2*intt[w]/i1/i1;
printf("Li1 = %g\n", l1);
この計算で得られた内部インダクタンス 8.30e-9 は解析的計算による正しい値の 2/1000 の誤差になります。(注4)

円筒導体を同軸ケーブルの外部導体として使う場合は、電流の帰路を円筒の内部に 密着した円柱にしなければなりませんが、そうすると、円筒の内部表面のベクトル ポテンシャルが一定という境界条件になります。円筒の外表面では、

  (d/dn)(A) = 0
すなわち、FreeFEMの記法では
  onbdy(ib) dnu(u) = 0;
という自然境界条件になります。つまり、ベクトルポテンシャルの法線成分は 0 です。

/* DC inductance calculation */
/* tubular outer conductor */
a1:=2.45e-3;	/* inner radius */
a2:=2.80e-3;	/* outer radius */
u1:=4e-7*pi;	/* permeability */
t1:=pi/4;	/* computing 0 - pi/4 */

n1:=10;
n2:=floor(n1*a1*t1/(a2-a1));
n3:=floor(n1*a2*t1/(a2-a1));

border(1,0,n1+n3+n1+n2,n1+n3+n1+n2+1) {
 if(t<=n2) { t2:=t1-t1*t/n2; x:=a1*cos(t2); y:=a1*sin(t2); ib:=2 }
 else if(t<n2+n1) { x:=a1+(t-n2)*(a2-a1)/n1; y:= 0; ib:=1 }
 else if(t<=n2+n1+n3) {
  t2:=(t-n2-n1)*t1/n3; x:=a2*cos(t2); y:=a2*sin(t2); ib:=3
 }
 else {
  t2:=a2-(a2-a1)*(t-n2-n1-n3)/n1; x:=t2*cos(t1); y:=t2*sin(t1); ib:=1
 }
};

buildmesh(1000,1);

i1:=1;	/* total current */
j1:=i1/pi/(a2*a2-a1*a1);	/* current density (total current == 1) */
/* vector potential */
solve(u) {
  onbdy(2) u = 0;
  pde(u) -laplace(u)/u1 = j1;
};

plot (u);

/* magnetic energy */
w=0.5*u*j1*2*pi/t1;

l1:=2*intt[w]/i1/i1;
printf("Li1 = %g\n", l1);
この計算で得られる内部インダクタンスの値 9.48e-9 は正しい値の 3/1000 の誤差 です。

2.3 同軸ケーブルのインダクタンス

以上を組み合わせて、同軸ケーブルの直流に於けるインダクタンスを計算すると 次のようになります。電流密度の符号が内部導体と外部導体で逆になることに注意 してください。誘電体部分では 0 になります。


/* DC inductance calculation */
/* coaxial cable: 5C-2V */
a1:=0.4e-3;	/* inner conductor raius */
a2:=2.45e-3;	/* insulation radius */
a3:=2.80e-3;	/* outer conductor radius */
u1:=4e-7*pi;	/* permeability */
t1:=pi/8;	/* computing 0 - pi/4 */

n1:=20;
n2:=floor(n1*(a2-a1)/a1);
n3:=floor(n1*(a3-a2)/a1);
n4:=floor(n1*a3*t1/a1);
n5:=floor(n1*a1*t1/a1);
n6:=floor(n1*a2*t1/a1);

border(1,0,n1+n5+n1,n1+n5+n1+1) {
 if(t<n1) { x:=a1*t/n1; y:= 0; ib:=1 }
 else if(t<=n1+n5) { t2:=(t-n1)*t1/n5; x:=a1*cos(t2); y:=a1*sin(t2); ib:=0 }
 else { t2:=a1-a1*(t-n1-n5)/n1; x:=t2*cos(t1); y:=t2*sin(t1); ib:=2 }
};
border(1,0,n2+n6+n2,n2+n6+n2+1) {
 if(t<n2) { x:=a1+(a2-a1)*t/n2; y:= 0; ib:=1 }
 else if(t<=n2+n6) { t2:=(t-n2)*t1/n6; x:=a2*cos(t2); y:=a2*sin(t2); ib:=0 }
 else { t2:=a2-(a2-a1)*(t-n2-n6)/n2; x:=t2*cos(t1); y:=t2*sin(t1); ib:=2 }
};
border(1,0,n3+n4+n3,n3+n4+n3+1) {
 if(t<n3) { x:=a2+(a3-a2)*t/n3; y:= 0; ib:=1 }
 else if(t<=n3+n4) { t2:=(t-n3)*t1/n4; x:=a3*cos(t2); y:=a3*sin(t2); ib:=3 }
 else { t2:=a3-(a3-a2)*(t-n3-n4)/n3; x:=t2*cos(t1); y:=t2*sin(t1); ib:=2 }
};

buildmesh(20000,1);

i1:=1;	/* total current */
j1:=i1/pi/a1/a1;	/* current density (inner conductor) */
j2:=i1/pi/(a3*a3-a2*a2);	/* current density (outer conductor) */
j=j1*(region==0)-j2*(region==2);
/* vector potential */
solve(u) {
  onbdy(3) u = 0;
  pde(u) -laplace(u)/u1 = j;
};

plot (u);

/* magnetic energy */
w=0.5/u1*(dx(u)*dx(u)+dy(u)*dy(u))*2*pi/t1;

printf("B-H method:\n");
l1:=2*intt(0)[w]/i1/i1;
printf("Li1 = %g H/m\n", l1);
l2:=2*intt(1)[w]/i1/i1;
printf("Le = %g H/m\n", l2);
l3:=2*intt(2)[w]/i1/i1;
printf("Li2 = %g H/m\n", l3);
l0:=l1+l2+l3;
printf("L = %g H/m\n", l0);

printf("A-i method:\n");
w=0.5*u*j*2*pi/t1;
/*l0:=2*intt(0,2)[w]/i1/i1;*/
l0:=2*intt[w]/i1/i1;
printf("L = %g H/m\n", l0);
対称性を利用して 1/16 の領域で計算していますが、解析的計算による 4.21982e-07 と比較して、誤差は 1 % 未満です。

2.4 撚線導体同軸ケーブルのインダクタンス

中心導体が7本撚で外部導体が円筒の場合は次のようになります。対称性を利用して、 断面の 1/12 を解析しています。


/* DC inductance calculation */
/* stranded conductor 3D-2V */
a1:=0.16e-3;	/* center conductor radius */
a2:=1.50e-3;	/* outer conductor radius */
u1:=4e-7*pi;	/* permeability */
t1:=pi/6;	/* computing 0 - pi/4 */

n1:=20;
n2:=floor(n1*2*a1/a1);
n3:=floor(n1*(a2-3*a1)/a1);
n4:=floor(n1*a2*t1/a1);
n5:=floor(n1*(a2-2*a1)/a1);
n6:=floor(n1*a1/a1);
n7:=n1;
n8:=floor(n1*a1*t1/a1);
n9:=floor(n1*a1*pi/a1);

border(1,0,n1+n8+n7,n1+n8+n7+1) {
 if(t<n1) { x:=a1*t/n1; y:= 0; ib:=1 }
 else if(t<=n1+n8) {
  t2:=(t-n1)*t1/n8; x:=a1*cos(t2); y:=a1*sin(t2); ib:=2
 }
 else {
  t2:=a1-(t-n1-n8)*a1/n7; x:=t2*cos(t1); y:=t2*sin(t1); ib:=1
 }
};
border(1,0,n2+n3+n4+n5+n6,n2+n3+n4+n5+n6+1) {
 if(t<n2) { x:=a1+2*a1*t/n2; y:=0; ib:=1 }
 else if(t<n2+n3) { x:=3*a1+(a2-3*a1)*(t-n2)/n3; y:=0; ib:=1 }
 else if(t<=n2+n3+n4) {
  t2:=t1*(t-n2-n3)/n4; x:=a2*cos(t2); y:=a2*sin(t2); ib:=3
 }
 else if(t<n2+n3+n4+n5) {
  t2:=a2-(a2-2*a1)*(t-n2-n3-n4)/n5; x:=t2*cos(t1); y:=t2*sin(t1); ib:=1
 }
 else {
  t2:=2*a1-a1*(t-n2-n3-n4-n5)/n6; x:=t2*cos(t1); y:=t2*sin(t1); ib:=1
 }
};
border(0,0,n9,n9+1) { t2:=pi*t/n9; x:=2*a1+a1*cos(t2); y:=a1*sin(t2); ib:=4 };

buildmesh(9000,1);

i1:=1;	/* total current */
j1:=i1/pi/a1/a1/7;	/* current density (inner conductor) */
j=j1*(region==0)+j1*(region==2);
/* vector potential */
solve(u) {
  onbdy(3) u = 0;
  pde(u) -laplace(u)/u1 = j;
};

plot (u);

/* magnetic energy - (1/2/u1)*abs(rot(A))^2 */
w=0.5/u1*(dx(u)*dx(u)+dy(u)*dy(u))*2*pi/t1;

l1:=2*intt(0)[w]/i1/i1;
printf("Li = %g\n", l1);
l2:=2*intt(1)[w]/i1/i1;
printf("Le = %g\n", l2);
l0:=l1+l2;
printf("L = %g\n", l0);
内部インダクタンスと外部インダクタンスの合計を、内部導体の外接円を半径とする 円柱導体を使った下記のケースと比較すると、8 % 程度小さくなることがわかります。 つまり、外接円を円柱に置き換えた近似計算の誤差は 8 % 程度です。キャパシタン スの計算で使われるMeyersの実効径の補正はインダクタンスについては無意味 で、間違ってこれを使うと 13 % の過大評価になることに注意してください。

/* DC inductance calculation */
/* coaxial cable: 5C-2V */
a1:=0.48e-3;	/* inner conductor raius */
a2:=1.50e-3;	/* insulation radius */
u1:=4e-7*pi;	/* permeability */
t1:=pi/4;	/* computing 0 - pi/4 */

n1:=30;
n2:=floor(n1*(a2-a1)/a1);
n3:=floor(n1*a1*t1/a1);
n4:=floor(n1*a2*t1/a1);

border(1,0,n1+n3+n1,n1+n3+n1+1) {
 if(t<n1) { x:=a1*t/n1; y:= 0; ib:=1 }
 else if(t<=n1+n3) { t2:=(t-n1)*t1/n3; x:=a1*cos(t2); y:=a1*sin(t2); ib:=2 }
 else { t2:=a1-a1*(t-n1-n3)/n1; x:=t2*cos(t1); y:=t2*sin(t1); ib:=1 }
};
border(1,0,n2+n4+n2,n2+n4+n2+1) {
 if(t<n2) { x:=a1+(a2-a1)*t/n2; y:= 0; ib:=1 }
 else if(t<=n2+n4) { t2:=(t-n2)*t1/n4; x:=a2*cos(t2); y:=a2*sin(t2); ib:=3 }
 else { t2:=a2-(a2-a1)*(t-n2-n4)/n2; x:=t2*cos(t1); y:=t2*sin(t1); ib:=1 }
};

buildmesh(8000,1);

i1:=1;	/* total current */
j1:=i1/pi/a1/a1;	/* current density (inner conductor) */
j=j1*(region==0);
/* vector potential */
solve(u) {
 onbdy(3) u = 0;
 pde(u) -laplace(u)/u1 = j;
};

plot (u);

/* mugnetic energy - (1/2/u1)*abs(rot(A))^2 */
w=0.5/u1*(dx(u)*dx(u)+dy(u)*dy(u))*2*pi/t1;

l1:=2*intt(0)[w]/i1/i1;
printf("Li1 = %g\n", l1);
l2:=2*intt(1)[w]/i1/i1;
printf("Le = %g\n", l2);
l0:=l1+l2;
printf("L = %g\n", l0);

2.5 横巻円筒同軸ケーブルのインダクタンス

オーディオケーブルの分野で革命を起こした2497の構造で横巻円筒同軸ケーブル のインダクタンスを計算すると、次のようになります。


/* served coaxial cable */
a1:=2.03e-3;
a2:=2.78e-3;
a3:=4.0e-3;
b1:=0.13e-3;
u1:=4e-7*pi;	/* permeability */
m0:=44;
m1:=11;
m2:=15;
t1:=2*pi*m1/m0;
n0:=30;
n1:=floor(a1/2/pi/b1*n0);
n2:=floor((a2-a1)/2/pi/b1*n0);
n3:=floor((a3-a2)/2/pi/b1*n0);
n4:=floor(a3*t1/2/pi/b1*n0);

border(1,0,n1+n2+n3+n4+n3+n2+n1,n1+n2+n3+n4+n3+n2+n1+1) {
 if(t<=n1) { x:=t*a1/n1; y:=0; ib:=1 }
 else if(t<=n1+n2) { x:=a1+(t-n1)*(a2-a1)/n2; y:=0; ib:=1 }
 else if(t<=n1+n2+n3) { x:=a2+(t-n1-n2)*(a3-a2)/n3; y:=0; ib:=1 }
 else if(t<n1+n2+n3+n4)
  { t2:=(t-n1-n2-n3)*t1/n4; x:=a3*cos(t2); y:=a3*sin(t2); ib:=2 }
 else if(t<=n1+n2+n3+n4+n3)
  { t2:=a3-(t-n1-n2-n3-n4)*(a3-a2)/n3; x:=t2*cos(t1); y:=t2*sin(t1); ib:=1 }
 else if(t<=n1+n2+n3+n4+n3+n2)
  { t2:=a2-(t-n1-n2-n3-n4-n3)*(a2-a1)/n2; x:=t2*cos(t1); y:=t2*sin(t1); ib:=1 }
 else
  { t2:=a1-(t-n1-n2-n3-n4-n3-n2)*a1/n1; x:=t2*cos(t1); y:=t2*sin(t1); ib:=1 };
};
for(i1:=0;i1<m1;i1++) {
 t2:=t1/m1*i1+t1/m1/2;
 border(0,0,2*pi,n0) {
  x1:=(a1+b1)*cos(t2); y1:=(a1+b1)*sin(t2);
	x:=x1+b1*cos(t); y:=y1+b1*sin(t); ib:=0
 };
};
for(i1:=0;i1<m2;i1++) {
 t2:=t1/m2*i1+t1/m2/2;
 border(0,0,2*pi,n0) {
  x1:=(a2+b1)*cos(t2); y1:=(a2+b1)*sin(t2);
	x:=x1+b1*cos(t); y:=y1+b1*sin(t); ib:=0
 };
};
buildmesh(50000,1);

i1:=1;	/* total current */
j1:=i1/pi/b1/b1/m0;	/* current density */
j2:=j1*m1/m2;	/* current density */
j=j1*one(region>0&®ion<=m1)-j2*one(region>m1);

/* vector potential */
solve(u) {
 onbdy(2) u = 0;
 pde(u) -laplace(u)/u1 = j;
};
plot(u);

w=0.5/u1*(dx(u)*dx(u)+dy(u)*dy(u));
plot(w);

w=0.5*u*j;
l1:=2*intt[w]/i1/i1*2*pi/t1;
printf("Li = %g (H/m, A-i method)\n", l1);

w=0.5/u1*(dx(u)*dx(u)+dy(u)*dy(u));
l1:=2*intt[w]/i1/i1*2*pi/t1;
printf("Li = %g (H/m, B-H method)\n", l1);

横巻素線の外径と同じ厚さの円筒導体の場合と比べてみると、 1 % 以下の差しかないようです。

2.6 横巻円筒導体の内部インダクタンス

横巻円筒導体の内部インダクタンスを計算してみると、例えば、次のようになって、 ほぼ同じ厚さの円筒導体の内部インダクタンスと等しいことがわかります。


/* DC inductance calculation */
/* served cylindrical conductor */
a1:=2.3e-3;	/* core radius */
a2:=0.09e-3;	/* served conductor radius */
n0:=floor(2*pi*(a1+a2)/(2*a2));	/* number of strucnds */
u1:=4e-7*pi;	/* permeability */
t1:=asin(a2/(a1+a2));
t2:=acos(a2/(a1+a2));
t3:=pi-t2;
a3:=sqrt(a1*a1+2*a1*a2);
a4:=a1+2*a2-a3;
a5:=a1+2*a2;

n2:=20;
n1:=floor(n2*a1/a2);
n3:=n2;
n4:=floor(n2*a5*t1/a2);
n5:=floor(n2*a4/a2);
n6:=floor(n2*a3/a2);
n7:=floor(n2*a2*t3/a2);
n8:=floor(n2*a2*t2/a2);

/*
border(1,0,n7+n8+n2+n3,n7+n8+n2+n3+1) {
 if(t<n7) { t4:=t*t3/n7; x:=a1+a2+a2*cos(t4); y:=a2*sin(t4); ib:=0 }
 else if(t<n7+n8) {
  t4:=t3+(t-n7)*t2/n8; x:=a1+a2+a2*cos(t4); y:=a2*sin(t4); ib:=0 }
 else if(t<n7+n8+n2) { x:=a1+(t-n7-n8)*a2/n2; y:=0; ib:=1 }
 else { x:=a1+a2+(t-n7-n8-n2)*a2/n2; y:=0; ib:=1 };
};
*/
border(1,0,n2+n3+n7+n8,n2+n3+n7+n8+1) {
 if(t<=n2+n3) { x:=a1+t*a2/n2; y:=0; ib:=1 }
 else if(t<n2+n3+n7) {
  t4:=(t-n2-n3)*t3/n7; x:=a1+a2+a2*cos(t4); y:=a2*sin(t4); ib:=0 }
 else {
  t4:=t3+(t-n2-n3-n7)*t2/n8; x:=a1+a2+a2*cos(t4); y:=a2*sin(t4); ib:=0 }
};
border(2,0,n4+n5,n4+n5+1) {
 if(t<n4) { t4:=t*t1/n4; x:=a5*cos(t4); y:=a5*sin(t4); ib:=2 }
 else { t4:=a5-(t-n4)*a4/n5; x:=t4*cos(t1); y:=t4*sin(t1); ib:=3 };
};
border(3,0,n6+n1,n6+n1+1) {
 if(t<n6) { t4:=a3-t*a3/n6; x:=t4*cos(t1); y:=t4*sin(t1); ib:=3 }
 else { x:=(t-n6)*a1/n1; y:=0; ib:=1 };
};
buildmesh(9000,1);

i1:=1;	/* total current */
i2:=i1/pi/a2/a2/n0;	/* current density (inner conductor) */
i=i2*(region==0);
/* vector potential */
solve(u) {
  onbdy(2) u = 0;
  pde(u) -laplace(u)/u1 = i;
};

plot (u);

/* magnetic energy - (1/2/u1)*abs(rot(A))^2 */
w=0.5/u1*(dx(u)*dx(u)+dy(u)*dy(u))*2*pi/t1;
plot (w);

l1:=2*intt(0)[w]/i1/i1;
printf("Li = %g\n", l1);
l2:=2*intt(1,2)[w]/i1/i1;
printf("Le = %g\n", l2);
l0:=l1+l2;
printf("L = %g\n", l0);

2.7 円柱導体ツインリードのインダクタンス

ツインリード(平行線)の場合は磁界が無限遠まで広がっていますし、シールド付きの 場合も銅やアルミの比透磁率は 1 ですから、電界と違って磁界はシールドの外まで 広がるため、インダクタンスについては開領域の計算が必要です。

開領域の計算方法として最も簡単なのは、磁界が十分弱くなるまで離れた 人為的境界を作る方法ですが、例えば、次のようになります。


/* DC inductance calculation */
/* twin lead */
a1:=0.225e-3;	/* conductor radius */
b1:=0.50e-3;	/* conductor separation */
b2:=b1*6;	/* the end of world */
u1:=4e-7*pi;	/* permeability */
t1:=pi/2;

n1:=10;
n2:=floor(n1*pi*a1/a1);
n3:=floor(n1*(b2-b1-a1)/a1);
n4:=floor(n1*b2*t1/a1);
n5:=floor(n1*b2/a1);
n6:=floor(n1*(b1-a1)/a1);

border(1,0,2*n1+n2,2*n1+n2+1) {
 if(t<=2*n1) { x:=b1-a1+t*a1/n1; y:=0; ib:=1 }
 else { t2:=(t-2*n1)*pi/n2; x:=b1+a1*cos(t2); y:=a1*sin(t2); ib:=0 };
};
border(2,0,n3+n4+n5+n6,n3+n4+n5+n6+1) {
 if(t<n3) { x:=b1+a1+t*(b2-b1-a1)/n3; y:=0; ib:=1 }
 else if(t<n3+n4) { t2:=(t-n3)*pi/2/n4; x:=b2*cos(t2); y:=b2*sin(t2); ib:=2 }
 else if(t<=n3+n4+n5) { x:=0; y:=b2-(t-n3-n4)*b2/n5; ib:=3 }
 else { x:=(t-n3-n4-n5)*(b1-a1)/n6; y:=0; ib:=1 };
};
buildmesh(50000,1);

i1:=1;	/* total current */
i2:=i1/pi/a1/a1;	/* current density (inner conductor) */
i=i2*(region==0);
/* vector potential */
solve(u) {
  onbdy(2,3) u = 0;
  pde(u) -laplace(u)/u1 = i;
};

plot (u);

/* magnetic energy - (1/2/u1)*abs(rot(A))^2 */
w=0.5/u1*(dx(u)*dx(u)+dy(u)*dy(u))*2*pi/t1;
plot (w);

l1:=2*intt(0)[w]/i1/i1;
printf("Li = %g\n", l1);
l2:=2*intt(1)[w]/i1/i1;
printf("Le = %g\n", l2);
l0:=l1+l2;
printf("L = %g\n", l0);
対称性を利用して 1/4 の領域で計算していますが、解析的計算で得られる値 7.056e-7 に対して 10 % の誤差です。解析しなければならない面積が広いですから、 計算量のかなり大きくなります。

もう少し効率のよい方法としては、境界条件が計算できるような仮想境界を作って、

  W = (1/2)*∫(B*H)*dU + (1/2)*∫(HxA)*dS                             (10)
  ここに
	W = 磁界のエネルギ (joul/m^3)
	B = 磁束密度 (weber/m^2)
	H = 磁界の強さ (ampare-turn/m)
	A = 磁界のベクトルポテンシャル (weber/m)
で計算する方法があります。最初の積分は境界内部の磁気エネルギ、2番目の積分は 境界外部の磁気エネルギを意味します。(HxA)は磁界と磁気ベクトルポテンシャルの ベクトル積ですが、H は z 軸に垂直、A は z 軸に平行ですから、(HxA) は z 軸と H に垂直な方向になり、その境界とのスカラ積を境界にそって線積分することに なります。空間の全エネルギは (7) 式でも計算できますが、導体内部という狭い範囲 でしか積分しないため、精度は低下します。

この仮想境界が使えるためには、境界上のベクトルポテンシャルが決まらなければ なりませんが、円柱導体外部のベクトルポテンシャルは下記のように簡単に計算でき ます(注6)ので、すべての導体についてこのベクトルポテンシャルを合計すれば、 境界面のベクトルポテンシャルが決まり、境界条件が得られます。

  A = μI/(2*π)*(1/2 + log(r/a))                                     (11)
  ここに
	μ = 透磁率 (H/m)
	     4e-7*π (真空、空気、銅、アルミなどのとき)
	I = 円柱導体の全電流 (A)
	r = 円柱導体中心からの距離 (m)
	a = 円柱導体の半径 (m)
この方法で上記と同じ問題を計算すると、下記のようになります。


/* DC inductance calculation */
/* twin lead */
a1:=0.225e-3;	/* conductor radius */
b1:=0.50e-3;	/* conductor separation */
c1:=b1*2;	/* ficsious boundary */
c2:=a1*2;	/* ficsious boundary */
u1:=4e-7*pi;	/* permeability */
t1:=pi/2;

n1:=30;
n2:=floor(n1*pi*a1/a1);
n3:=floor(n1*(c1-b1-a1)/a1);
n4:=floor(n1*c2/a1);
n5:=floor(n1*c1/a1);
n6:=floor(n1*c2/a1);
n7:=floor(n1*(b1-a1)/a1);

border(1,0,2*n1+n2,2*n1+n2+1) {
 if(t<=2*n1) { x:=b1-a1+t*a1/n1; y:=0; ib:=1 }
 else { t2:=(t-2*n1)*pi/n2; x:=b1+a1*cos(t2); y:=a1*sin(t2); ib:=0 };
};
border(2,0,n3+n4+n5+n6+n7,n3+n4+n5+n6+n7+1) {
 if(t<n3) { x:=b1+a1+t*(c1-b1-a1)/n3; y:=0; ib:=1 }
 else if(t<n3+n4) { x:=c1; y:=(t-n3)*c2/n4; ib:=2 }
 else if(t<n3+n4+n5) { x:=c1-(t-n3-n4)*c1/n5; y:=c2; ib:=3 }
 else if(t<=n3+n4+n5+n6) { x:=0; y:=c2-(t-n3-n4-n5)*c2/n6; ib:=4 }
 else { x:=(t-n3-n4-n5-n6)*(b1-a1)/n7; y:=0; ib:=1 };
};
buildmesh(20000,1);

i1:=1;	/* total current */
j2:=i1/pi/a1/a1;	/* current density (inner conductor) */
j=j2*(region==0);

/* vector potential */
solve(u) {
 onbdy(2,3) u = u1*i1/4/pi*log((x*x+y*y+b1*b1+2*b1*x)/(x*x+y*y+b1*b1-2*b1*x)); 
 onbdy(4) u = 0;
 pde(u) -laplace(u)/u1 = j;
};

plot (u);

/* magnetic energy - (1/2/u1)*abs(rot(A))^2 */
w=0.5/u1*(dx(u)*dx(u)+dy(u)*dy(u));
b=-0.5*dx(u)/u1*u;	/* H×A・dS at x = const */
c=-0.5*dy(u)/u1*u;	/* H×A・dS at y = const */

/*printf("x(%g),y(%g)\n",int(2)[b],int(3)[c]);*/

l1:=2*intt(0)[w]/i1/i1*2*pi/t1;
printf("Li = %g\n", l1);
l2:=2*intt(1)[w]/i1/i1*2*pi/t1;
printf("Le = %g\n", l2);
l3:=2*(int(2)[b]+int(3)[c])/i1/i1*2*pi/t1;
printf("Ls = %g\n", l3);
printf("L = %g (B-H)\n", l1+l2+l3);

w=0.5*u*j;
l0:=2*intt(0)[w]/i1/i1*2*pi/t1;
printf("L = %g (A-i)\n", l0);

printf("Lavarage=%g\n", (l1+l2+l3+l0)/2);

この場合は、解析的計算で得られる値と 0.9 % しが違いません。計算量も激減して います。 (10) 式の B-H を使う方法と (7) 式の A-i を使う方法の誤差は逆向きにでるようで、 この2つの平均値を使うと精度がよくなって、誤差は 0.03 % になります。

3. 直流に於ける相互インダクタンスの計算

3.1 フラットケーブルの相互インダクタンス

相互インダクタンスの計算方法としては、2つの回路を並列接続したときの 合成インダクタンスの計算を元にする方法と、錯交磁束を元にする方法がありますが、 前者は理論的に正確である半面、ごくわずかな相互インダクタンスを 2つのインダクタンスの差から求めるため誤差が大きくなります。 後者は導体内部の磁束の不均一を無視する近似計算になりますが、計算誤差が少なく、 一般には後者が無難です。

まず、最初に前者の方法で、4芯のフラットケーブルを2回路並列に接続したときの 合成インダクタンスを計算してみると、下記のようになります。 対称性を利用して、右半分の上側だけ、つまり 1/4 の領域を解析しています。


/* DC inductance calculation */
/* flat cable - parallel connection */
a1:=0.19e-3;	/* conductor radius */
s1:=1.27e-3;
s2:=s1+2*s1;
c1:=s2*1.5;	/* ficsious boundary */
c2:=a1*3;	/* ficsious boundary */
u1:=4e-7*pi;	/* permeability */
t1:=pi/2;
s3:=(s2+s1)/2;

n1:=15;
n2:=floor(n1*pi*a1/a1);
n3:=floor(n1*(s1-a1)/a1);
n4:=floor(n1*2*(s1-a1)/a1);
n5:=floor(n1*(c1-s2-a1)/a1);
n6:=floor(n1*c2/a1);
n7:=floor(n1*c1/a1);
border(1,0,n3+2*n1+n4+2*n1+n5+n6+n7+n6,n3+2*n1+n4+2*n1+n5+n6+n7+n6+1) {
 if(t<=n3) { x:=t*(s1-a1)/n3; y:=0 }
 else if(t<=n3+2*n1) { x:=s1-a1+(t-n3)*a1/n1; y:=0 }
 else if(t<=n3+2*n1+n4) { x:=s1+a1+(t-n3-2*n1)*(s2-s1-2*a1)/n4; y:=0 }
 else if(t<=n3+2*n1+n4+2*n1) {
  x:=s2-a1+(t-n3-2*n1-n4)*a1/n1; y:=0 }
 else if(t<=n3+2*n1+n4+2*n1+n5) {
  x:=s2+a1+(t-n3-2*n1-n4-2*n1)*(c1-s2-a1)/n5; y:=0 }
 else if(t<n3+2*n1+n4+2*n1+n5+n6) {
  x:=c1; y:=(t-n3-2*n1-n4-2*n1-n5)*c2/n6; ib:=2 }
 else if(t<n3+2*n1+n4+2*n1+n5+n6+n7) {
  x:=c1-(t-n3-2*n1-n4-2*n1-n5-n6)*c1/n7; y:=c2; ib:=3 }
 else {
  x:=0; y:=c2-(t-n3-2*n1-n4-2*n1-n5-n6-n7)*c2/n6; ib:=4 }
};
for(i1:=0;i1<2;i1++) {
 s4:=s1+2*s1*i1;
 border(1,0,n2,n2+1) { t2:=t*pi/n2; x:=s4+a1*cos(t2); y:=a1*sin(t2); ib:=0 };
};
buildmesh(20000,1);

i1:=1;	/* total current */
j1:=i1/pi/a1/a1/2;	/* current density */

j=j1*(region==1)-j1*(region==2);

/* vector potential */
solve(u) {
/*
 onbdy(2) u = u1*i1/4/pi*log((x*x+y*y+s1*s1-2*s1*x)/(x*x+y*y+s1*s1+2*s1*x));
*/
/*
 onbdy(2) u = u1*i1/4/pi
  *log(((x+s3)*(x+s3)+y*y+s1*s1-2*s1*(x+s3))
   /((x+s3)*(x+s3)+y*y+s1*s1+2*s1*(x+s3)));
*/
/*
 onbdy(2) u = u1*i1/4/pi
  *log(((x-s3)*(x-s3)+y*y+s1*s1-2*s1*(x-s3))
   /((x-s3)*(x-s3)+y*y+s1*s1+2*s1*(x-s3)));
 pde(u) -laplace(u)/u1 = j;
*/
/*
 onbdy(2,3) u = u1*i1/4/pi
  *log(((x+s3)*(x+s3)+y*y+s1*s1+2*s1*(x+s3))
   /((x+s3)*(x+s3)+y*y+s1*s1-2*s1*(x+s3))
  *((x-s3)*(x-s3)+y*y+s1*s1-2*s1*(x-s3))
   /((x-s3)*(x-s3)+y*y+s1*s1+2*s1*(x-s3)));
*/
/* onbdy(4) u=u1*i1/4/pi*log((s2*s2+y*y)/(s1*s1+y*y));*/

 onbdy(2,3) u = u1*i1/4/pi/2
  *log((x*x+y*y+s1*s1+2*s1*x)
   /(x*x+y*y+s1*s1-2*s1*x)
  *(x*x+y*y+s2*s2-2*s2*x)
   /(x*x+y*y+s2*s2+2*s2*x));
 onbdy(4) u = 0;
 pde(u) -laplace(u)/u1 = j;
};

plot (u);

/* magnetic energy - (1/2/u1)*abs(rot(A))^2 */
w=0.5/u1*(dx(u)*dx(u)+dy(u)*dy(u));
b=-0.5*dx(u)/u1*u;	/* H×A・dS at x = const */
c=-0.5*dy(u)/u1*u;	/* H×A・dS at y = const */
printf("x(%g),y(%g)\n",int(2,4)[b],int(3)[c]);


l1:=2*intt(1,2,3)[w]/i1/i1*2*pi/t1
 +2*intt(4)[w]/i1/i1*2*pi/t1;
printf("Li = %g\n", l1);
l2:=2*intt(0)[w]/i1/i1*2*pi/t1;
printf("Le = %g\n", l2);
l3:=2*(int(2,4)[b]+int(3)[c])/i1/i1*2*pi/t1;
printf("Ls = %g (B-H)\n", l3);
printf("L = %g (B-H)\n", l1+l2+l3);

w=0.5*u*j;
l0:=2*(intt(1,2)[w]+intt(3,4)[w])/i1/i1*2*pi/t1;
printf("L = %g (A-i)\n", l0);


printf("Lavarage=%g\n", (l1+l2+l3+l0)/2);

解析的計算で得られる値 5.11e-7 と比べると 4 % の誤差ですが、 合成インダクタンスが得られれば、
  M = 2*Lp - L
  ここに
	M = 相互インダクタンス (H/m)
	Lp = 2回路を並列接続したときの合成インダクタンス (H/m)
	L = 個々の回路のインダクタンス (H/m)
から計算できます。

ただし、M そのものの誤差は引き算の結果、55 % に拡大します。

錯交磁束を使う方法は、

  M12 = Φ12/I2
  ここに
	M12 = 回路 2 から回路 1 への相互インダクタンス (H/m)
	Φ12 = 回路 2 による磁束のうち回路 1 に錯交するもの (weber)
	I2 = 回路 2 の電流 (A)
を使うもので、磁束は下記の方法で求めるのが簡単です。
  Φ = ∫B・dS
     = ∫μH・dS
     = ∫A・dl
最初の2つは面積分、最後の1つは1周積分です。たくさんのデータを使う面積分を 使おうとすると微分が必要になりますし、直接ベクトルポテンシャルを使おうとすると 一部の領域だけで線積分しなければならず、一長一短ですが、前期の問題を この方法で計算すると、例えば下記のようになります。

/* DC inductance calculation */
/* flat cable */
a1:=0.19e-3;	/* conductor radius */
s1:=1.27e-3;
s2:=s1+2*s1;
c1:=s2*1.5;	/* the end of world */
c2:=a1*3;	/* the end of world */
u1:=4e-7*pi;	/* permeability */
t1:=pi;
s3:=(s2+s1)/2;

n1:=10;
n2:=floor(n1*pi*a1/a1);
n3:=floor(n1*(s1-a1)/a1);
n4:=floor(n1*2*(s1-a1)/a1);
n5:=floor(n1*(c1-s2-a1)/a1);
n6:=floor(n1*c2/a1);
n7:=floor(n1*2*c1/a1);
border(1,0,n3+2*n1+n4+2*n1+n5+n6+n7+n6+n5+2*n1+n4+2*n1+n3,
  n3+2*n1+n4+2*n1+n5+n6+n7+n6+n5+2*n1+n4+2*n1+n3+1) {
 if(t<=n3) { x:=t*(s1-a1)/n3; y:=0; ib:=1 }
 else if(t<=n3+n1) { x:=s1-a1+(t-n3)*a1/n1; y:=0; ib:=1 }
 else if(t<=n3+2*n1) { x:=s1+(t-n3-n1)*a1/n1; y:=0; ib:=5 }
 else if(t<=n3+2*n1+n4) { x:=s1+a1+(t-n3-2*n1)*(s2-s1-2*a1)/n4; y:=0; ib:=5 }
 else if(t<=n3+2*n1+n4+n1) {
  x:=s2-a1+(t-n3-2*n1-n4)*a1/n1; y:=0; ib:=1 }
 else if(t<=n3+2*n1+n4+2*n1+n5) {
  x:=s2+a1+(t-n3-2*n1-n4-2*n1)*(c1-s2-a1)/n5; y:=0; ib:=1 }
 else if(t<n3+2*n1+n4+2*n1+n5+n6) {
  x:=c1; y:=(t-n3-2*n1-n4-2*n1-n5)*c2/n6; ib:=2 }
 else if(t<n3+2*n1+n4+2*n1+n5+n6+n7) {
  x:=c1-(t-n3-2*n1-n4-2*n1-n5-n6)*2*c1/n7; y:=c2; ib:=3 }
 else if(t<n3+2*n1+n4+2*n1+n5+n6+n7+n6) {
  x:=-c1; y:=c2-(t-n3-2*n1-n4-2*n1-n5-n6-n7)*c2/n6; ib:=4 }
 else if(t<=n3+2*n1+n4+2*n1+n5+n6+n7+n6+n5) {
  x:=-c1+(t-n3-2*n1-n4-2*n1-n5-n6-n7-n6)*(c1-s2-a1)/n5; y:=0; ib:=1 }
 else if(t<=n3+2*n1+n4+2*n1+n5+n6+n7+n6+n5+2*n1) {
  x:=-s2-a1+(t-n3-2*n1-n4-2*n1-n5-n6-n7-n6-n5)*a1/n1; y:=-; ib:=1 }
 else if(t<=n3+2*n1+n4+2*n1+n5+n6+n7+n6+n5+2*n1+n4) {
  x:=-s2+a1+(t-n3-2*n1-n4-2*n1-n5-n6-n7-n6-n5-2*n1)*(s2-s1-2*a1)/n4;
  y:=0; ib:=1 }
 else if(t<=n3+2*n1+n4+2*n1+n5+n6+n7+n6+n5+2*n1+n4+2*n1) {
  x:=-s1-a1+(t-n3-2*n1-n4-2*n1-n5-n6-n7-n6-n5-2*n1-n4)*a1/n1;
  y:=0; ib:=1 }
 else {
  x:=-s1+a1+(t-n3-2*n1-n4-2*n1-n5-n6-n7-n6-n5-2*n1-n4-2*n1)*(s1-a1)/n3;
  y:=0; ib:=1 };
};
for(i1:=0;i1<4;i1++) {
 s4:=-s2+2*s1*i1;
 border(1,0,n2,n2+1) { t2:=t*pi/n2; x:=s4+a1*cos(t2); y:=a1*sin(t2); ib:=0 };
};
buildmesh(20000,1);

i1:=1;	/* total current */
j1:=i1/pi/a1/a1;	/* current density */

j=j1*(region==1)-j1*(region==2);
/*
j=j1*(region==2)-j1*(region==3);
j=j1*(region==1)-j1*(region==2);
j=j1*(region==3)-j1*(region==4);
j=j1*(region==1)-j1*(region==2)+j1*(region==3)-j1*(region==4);
*/

/* vector potential */
solve(u) {
/*
 onbdy(2) u = u1*i1/4/pi*log((x*x+y*y+s1*s1-2*s1*x)/(x*x+y*y+s1*s1+2*s1*x));
*/

 onbdy(2,3,4) u = u1*i1/4/pi
  *log(((x+s3)*(x+s3)+y*y+s1*s1-2*s1*(x+s3))
   /((x+s3)*(x+s3)+y*y+s1*s1+2*s1*(x+s3)));

/*
 onbdy(2) u = u1*i1/4/pi
  *log(((x-s3)*(x-s3)+y*y+s1*s1-2*s1*(x-s3))
   /((x-s3)*(x-s3)+y*y+s1*s1+2*s1*(x-s3)));
 pde(u) -laplace(u)/u1 = j;
*/
/*
 onbdy(2) u = u1*i1/4/pi
  *log(((x+s3)*(x+s3)+y*y+s1*s1-2*s1*(x+s3))
   /((x+s3)*(x+s3)+y*y+s1*s1+2*s1*(x+s3))
  *((x-s3)*(x-s3)+y*y+s1*s1-2*s1*(x-s3))
   /((x-s3)*(x-s3)+y*y+s1*s1+2*s1*(x-s3)));
*/
 pde(u) -laplace(u)/u1 = j;
};

plot (u);

/* magnetic energy - (1/2/u1)*abs(rot(A))^2 */
w=0.5/u1*(dx(u)*dx(u)+dy(u)*dy(u));
b=0.5*dx(u)/u1*u;	/* H×A・dS at x = const */
c=-0.5*dy(u)/u1*u;	/* H×A・dS at y = const */
printf("x(%g),y(%g)\n",int(2,4)[b],int(3)[c]);

l1:=2*intt(1,2)[w]/i1/i1*2*pi/t1;
printf("Li = %g\n", l1);
l2:=2*intt(0,3,4)[w]/i1/i1*2*pi/t1;
printf("Le = %g\n", l2);
l3:=2*(int(2,4)[b]+int(3)[c])/i1/i1*2*pi/t1;
printf("Ls = %g\n", l3);
printf("L = %g (B-H)\n", l1+l2+l3);

w=0.5*u*j;
l0:=2*intt(1,2)[w]/i1/i1*2*pi/t1;
printf("L = %g (A-i)\n", l0);

printf("Lavarage=%g\n", (l1+l2+l3+l0)/2);

f=-dx(u);	/* flux at y=0 */
m1:=int(5)[f];
printf("m1(%g)\n",m1);
精度は 4 % と大幅に向上しました。

なお、直流の相互インダクタンスについては、ほとんどの場合、 2組の往復導体 1-2 と 3-4 の間の相互インダクタンスを求める簡単な公式

  M = μ/(2*π)*log(R14*R23/(R13*R24))
  ここに
	M = 2組の往復導体間の相互インダクタンス (H/m)
	R14 = 導体 1 と導体 4 の距離
	R23 = 導体 2 と導体 3 の距離
	R13 = 導体 1 と導体 3 の距離
	R24 = 導体 2 と導体 4 の距離
で間に合うのが普通で、これは錯交磁束を求める方法で簡単に証明できます。

注1 - 静磁界の偏微分方程式

静磁界問題を解く方法としては

  1. 磁気ベクトルポテンシャルを使う方法
  2. 磁位を使う方法
  3. 磁界を未知数として使う方法
の3つがよく使われますが、FreeFEMを使う場合は、境界条件とエネルギ計 算の簡単さから、ベクトルポテンシャルを使うのが有利で、以下、この方法に限定 して説明します。

Maxwell の方程式

  rot(H) = J + (d/dt)(D)
  rot(E) = -(d/dt)(B)
  div(B) = 0
  div(D) = ρ

  B = μ*H
  D = ε*E
  J = σ*E

  ここに
	H = 磁界の強さ
	J = 電流密度
	E = 電界の強さ
	B = 電束密度
	ρ = 電荷密度
	μ = 透磁率
	ε = 誘電率
	σ = 導電率
に加えて、下記の定義によるベクトルポテンシャル(magnetic vector potantial)
  B = rot(A)

  ここに
	A = ベクトルポテンシャル
を考えますが、A に任意のベクトルを加えても B の値が決まりませんので、制約 条件として、次の「クーロンゲージ」(Coulomb convention)
  div(A) = 0
を使います。2次元場では常にクーロンゲージが満たされます。
  rot(B) = rot(rot(A))
	= grad(div(A) - laplace(A)
ですから
  rot(B) = -laplace(A)
になり、一方
  rot(B) = μ*rot(H) = -μ*(J + (d/dt)(D))
ですから
  laplace(A) = -μ*(J + (d/dt)D))
になります。静磁界の場合は
  (d/dt)(D) = 0
ですから、結局のところ
  laplace(A) = -μ*J
を解けばよいことになります。

注2 - 静磁界のエネルギ計算

エネルギの計算は

  B = rot(A)                                                           (4)

  ここに
	B = 磁束密度 (weber/m^2)
から磁束密度を計算して、
  W = (1/2)*∫(B*H)*dS
    = (1/2)*(1/μ)*∫abs(B)^2*dS                                       (5)

  ここに
	W = 磁界のエネルギ (joul/m^3)
	abs(x) = ベクトル x の絶対値
を使うこともできますが、微分操作が必要になりますので、精度が低下します。実際 に下記のような方法で計算して比較してください。
  /* magnetic energy - (1/2/u1)*abs(rot(A))^2 */
  w=0.5/u1*(dx(u)*dx(u)+dy(u)*dy(u));
  plot (w);

  l1:=2*intt[w]/i1/i1;
  printf("Li = %g (H/m)\n", l1);

注3 - 円柱導体の内部インダクタンス

円柱導体の内部インダクタンスの計算は簡単で、Ampere の貫流式

  Hr = (1/(2*π))*I/r

  ここに
	Hr = 導体中心から r 離れた場所の磁界の強さ (H/m)
	I = 距離 r 内部の全電流 (A)
	r = 導体中心からの距離 (m)
から導体内部の 磁界
  H = I*r/(2*π*a^2)

  ここに
	H = 磁界の強さ (amp-turn/m)
	I = 導体の全電流 (A)
	a = 導体半径 (m)
	r = 導体中心からの距離 (m)
から単位長さあたりに磁気エネルギを求め
       a
  U = ∫(1/2)*μ*H^2*2*π*r*dr
      0
                        a
    = I^2*μ/(4*π*a^4)∫r^3*dr
                       0
    = μ*I^2/(16*π)

電流、インダクタンス、エネルギの関係

  U = (1/2)*Li*I^2
を考慮すれば、
  Li = μ/(8*π)  (H/m)
     = 0.5e-6  (H/m) .. 銅、アルミ等の場合
が得られます。導体のサイズに無関係であることに注意してください。

注4 - 円筒導体の内部インダクタンス

円柱の場合と同じ方針で簡単に計算できますが、同軸ケーブルの中心導体として使 う場合は、円筒内部の磁界が

  H = (1/(2*π*r))*(r^2-a^2)/(b^2-a^2)

  ここに
	H = 円筒内部の磁界の強さ (A-turn/m)
	r = 円筒中心からの距離 (m)
	a = 円筒の内径 (m)
	b = 円筒の外径
になりますから、円筒内部の磁気エネルギは
       b
  U = ∫(1/2)*μ*H^2*2*π*r*dr
      a
    = (μ/(4*π*(b^2-a^2))*I^2*(a^4/(b^2-a^2)*log(b/a)-(3*a^2-b^2)/4)
になりますから、内部インダクタンスは、
    Li = (μ/(2*π*(b^2-a^2))*(a^4/(b^2-a^2)*log(b/a)-(3*a^2-b^2)/4)
になります。積分は簡単です。

円筒を同軸ケーブルの外部導体として使う場合は、

H = (1/(2*π*r))*(b^2-r^2)/(b^2-a^2) になりますから、内部インダクタンスは、

    Li = (μ/(2*π*(b^2-a^2))*(b^4/(b^2-a^2)*log(b/a)-(3*b^2-a^2)/4)
と少し増えます。

注5 - 同軸ケーブルのインダクタンス

同軸ケーブルの直流に於けるインダクタンスは、前記の

  1. 単独円柱導体の内部インダクタンス
  2. 同軸外部導体として使った場合の円筒導体の内部インダクタンス
  3. 外部インダクタンス
の合計になりますが、外部インダクタンスは誘電体部分の磁気エネルギ
       b
  U = ∫(1/2)*μ*H^2*2*π*r*dr
      a
       b
    = ∫(1/2)*μ*(I/(2*π*r))^2*2*π*r*dr
      a
    = (μ/(4*π))*log(b/a)
から
  Le = (μ/(2*π))*log(b/a)
     = 2e-7*log(b/a) .. 銅、アルミの場合
と簡単に計算できます。

注6 - 導体外部の空間に於けるベクトルポテンシャルの計算

導体外部の空間に於けるベクトルポテンシャルは (11) 式で計算できますが、 導体半径の等しい往復導体を対にして、電流双極子として扱うと、下記のように、 導体半径と 1/2 という定数消去できて計算が簡単になります。

  A = μI/(4*π)*log((x^2+y^2+s^2-2*s*x)/(x^2+y^2+s^2+2*s*x))
  ここに
        A = 電流双極子による磁気ベクトルポテンシャル (weber/m)
        μ = 透磁率 (H/m)
        I = 電流 (A)
        x = 2導体の中心を原点としたときの計算点の x 座標 (m)
        y = 2導体の中心を原点としたときの計算点の y 座標 (m)
        s = 2導体の距離の 1/2 (m)
この計算は余弦定理を使うと簡単です。

平林 浩一, 2000-07-17