有限要素法によるケーブル(電線)の電気設計 (1)

1960 年代に航空宇宙工学の構造解析で開発された有限要素法は、その後、偏微分方程式の数値解析の標準的な方法として確立され、今ではあらゆる技術分野の常識として、多量の文献がありますが、あまり勉強しなかった人には理解が困難ですから、社内で絶えず必要になる問題に絞って、有限要素法による初歩的なケーブル設計の実務を解説し、ある程度

  分からなくても使える

ことを目的に、そのまま使える豊富な例題を用意することにしました。

有限要素法のプログラムは実にたくさんありますが、ここではFreeFEMというフランスの大学で開発され、私がケーブル設計で使いやすいように少し拡張した版を使うことにします。オリジナルのFreeFEMではprintf()以外にも、for()文やインクリメント等の演算子、fmod()等の関数がありませんので、複雑なケーブルではかなり長いプログラムになってしまいます。

FreeFEMの全体像と構文の詳細についてはFreeFEMのドキュメントを参照してください。日本語の解説も用意してあります。

FreeFEMはborder()文で指定された形状の頂点座標を計算し、 Delauney-Boronoi三角分割を行った結果を表示し、solve()文で指定された偏微分方程式と境界条件を使って、問題を解くようになっています。このコースで用意された例題は、三角分割の結果と、偏微分方程式を解いて得られたポテンシャル分布の等高線を表示した後、ケーブルの特性を計算して表示するようになっていますので、グラフィック表示の後は、マウスの左釦をして先に進めてください。もし、マウスの操作を省いて、結果だけを表示したければ、プログラムの始めの方に、

  nowait;

を追加します。ただ、三角分割と等ポテンシャル線図をチェックしないと、とんでもない間違いを見逃す危険があります。

1. キャパシタンス、特性インピーダンス、速度係数の計算法

ケーブルの電気特性で最も基本になるのがキャパシタンス(capacitance)で、これが計算できると、次のようにして、高周波の特性インピーダンス(characteristic impedance) と速度係数(velocity ratio)が求まります。

  速度係数 = c*sqrt(C0/C)
  特性インピーダンス = 1/(c*sqrt(C0*C))
  ここに、
    c = 光速 (2.99792458e8)
    C = キャパシタンス
    C0 = 誘電体部分を真空にしたときのキャパシタンス

つまり、キャパシタンスの計算ができれば、特性インピーダンスと速度係数は簡単に分かるのですが、キャパシタンスが解析的に計算できるのは、同軸ケーブル等、かなり簡単な形状に限られ、一般には解析的な計算はできません。有限要素法とか差分法といった数値計算を使うことになります。

数値計算の場合は、スカラ・ポテンシャル(電位)を使った、２次元のポアソン方程式 (Poisson's eqation)

  (d/dx+d/dy)(e*(d/dx+d/dy)(V)) = q

  ここに、
    V(x,y) = 電位 (V)
    e(x,y) = 誘電率 (F/m)
    q(x,y) = 電荷 (C)

を実際の境界条件のもとで解いて、２つの導体間の電位差と電界のエネルギから、キャパシタンスを計算します。

もちろん、導体表面の電界強度と誘電率から、導体の表面電荷を求め、それを積分した導体の全電荷と導体間電位差から、注(0)の方法で、キャパシタンスを求めることもできますが、有限要素法の場合はエネルギの誤差はポテンシャルの誤差より少ないという際だった性格があって、エネルギを使うかほうが高い精度が得られます。 (1)

  C = 2*U/(V*V)  (F/m)

  ここに、
    C = キャパシタンス (F/m)
    U = 電界のエネルギ (joul/m^3)
    V = 電位差 (V)

エネルギは電位傾度(電界)の自乗を全領域に渡って積分すれば算出できますから、次のように計算します。

  U = ∫(e(x,y)*(1/2)*(((d/dx)V(x,y))^2 + ((d/dy)V(x,y))^2)))*dS

積分は導体以外の全領域が対象です。

開領域の場合は、解析対象外の空間のエネルギも考慮しなければなりませんが、これは、外部境界の周積分を使って、下記のように計算します。

  U = ∫(e(x,y)*(1/2)*(((d/dx)V(x,y))^2 + ((d/dy)V(x,y))^2)))*dS
	- ∫e(x,y)*(1/2)*V*dV/dn*dl

最初の項は閉領域と同じ断面の面積分ですが、第２項は外部境界上の線積分です。

なお、キャパシタンスの計算で電位分布を求めなければならないのは誘電体内部だけで、電荷が存在する導体の電位をディリクレ(Direcret)型境界条件として指定しますから、電荷を考慮する必要はなく、

  (d/dx+d/dy)(e*(d/dx+d/dy)*V) = 0

というラプラス方程式(Laplace's equation)を解くだけの簡単な仕事になり、 FreeFEMでは、laplace()演算子だけで間に合います。

以下、実例をもとに説明しますが、キャパシタンス、速度係数、特性インピーダンスは幾何学的に相似なら同じ値になりますから、寸法は mm 単位にします。m 単位でFreeFEMのスクリプトを書くと、値が小さくなるため、同じ場所ではないのに同じ場所と見られることがあります。

2. 同軸ケーブル

2.1. 円筒導体

最も簡単な同軸ケーブルから始めます。下記のFreeFEMスクリプトを引数にFreeFEMを起動すると、5C-2V(JIS C 3501) のキャパシタンスを計算します。1e7/(4*pi*2.998e8^2) は真空の誘電率です。この例では、誘電体が一つしかありませんので、誘電体を真空にしても電位分布は変わらず、微分方程式の求解も１度で済みます。


/* coaxial cable: 5C-2V */
a1:=0.4;	/* inner conductor raius */
a2:=2.45;	/* outer conductor radius */
eps:=2.3;	/* dielectric constant */

border(1,0,2*pi,120) { x:=a2*cos(t); y:=a2*sin(t) };
border(2,0,2*pi,20) { x:=a1*cos(-t); y:=a1*sin(-t) };
buildmesh(5000);

v1:=0.0; v2:=1.0;
solve(v) {
 onbdy(1) v = v1;
 onbdy(2) v = v2;
 pde(v) -laplace(v) = 0;
};
plot(v);

u=0.5*(dx(v)*dx(v)+dy(v)*dy(v));	/* electro static energy for air */
u0:=intt[u];
u1=eps*u0;	/* electro static energy for eps */
v:=sqrt(u0/u1);	/* velocity ratio */
c:=u1*17.707;	/* pF/m */
z:=1/(2.998e-4*v*c);	/* characteristic impedance */
printf("C = %.1f (pF/m)\n", c);
printf("Vr = %.3f\n", v);
printf("Z0 = %.1f\n", z);

キャパシタンスを計算するときの 17.707 という数字は、面積分 u0 に真空の誘電率

  eps0 = 1e7/(4*pi*(c^2))

を掛けてエネルギに換算し、pF 単位にしたものです。自分でも計算して確認してください。

同軸ケーブルのキャパシタンスについては次のような解析的な厳密解が得られます。

  C = 2*pi*eps/log(a2/a1) = 55.63e-12/log(a2/a1)  (F/m)

  ここに
    pi = 3.149265..
    eps = 誘電率 (真空の誘電率は 1e7/(4*pi*c*c))
    a2 = 外部導体内径 (m)
    a1 = 内部導体外径 (m)

計算結果は 69.8 pF/m、厳密解は 70.597 pF/m ですから、この程度の分解でも 1.1 % の精度が得られます。border文の４番目の引数でメッシュの細かさが調整できますので、いろいろ試してください。

JIS C 3501ではキャパシタンスが 67+-3 pF/m になっていますが、この中心値は少し低すぎます。外部導体が編組ですから、空隙によって誘電率は下るはずで、厳密解より小さくなることはあり得ません。

2.2. 偏芯

30 % 偏芯したらどうなるかを調べてみましょう。e1は偏芯率(0<=e1<1) です。


/* coaxial cable: 5C-2V eccentric */
a1:=0.4;	/* inner conductor raius */
a2:=2.45;	/* outer conductor radius */
eps:=2.3;	/* dielectric constant */
e:=0.30;	/* eccentric factor */

border(1,0,2*pi,120) { x:=a2*cos(t); y:=a2*sin(t) };
if (e>=1.0) then e:=0.99;
d1:=(a2-a1)*e;
border(2,0,2*pi,20) { x:=a1*cos(-t)+d1; y:=a1*sin(-t) };
buildmesh(5000);

v1:=0.0; v2:=1.0;
solve(v) {
 onbdy(1) v = v1;
 onbdy(2) v = v2;
 pde(v) -laplace(v) = 0;
};
plot(v);

u=0.5*(dx(v)*dx(v)+dy(v)*dy(v));	/* electro static energy for air */
u0:=intt[u];
u1=eps*u0;	/* electro static energy for eps */
v:=sqrt(u0/u1);	/* velocity ratio */
c:=u1 * 17.707;	/* pF/m */
z:=1/(2.998e-4*v*c);	/* characteristic impedance */
printf("C = %.1f (pF/m)\n", c);
printf("Vr = %.3f\n", v);
printf("Z0 = %.1f\n", z);

キャパシタンスで 3.7 % の増加、特性インピーダンスで 1.1 % の減少ですから、この２つの特性が偏芯に対して鈍感なことがわかります。

2.3. 撚線導体

8D-2Vの場合は内部導体が 7/0.18 という７本撚構造になっていて、その外接円の半径 1.2e-3 を使って上記の計算をすると、107 pF/m という約 5 % 大きな値が得られます。これは撚線外周のへこみを考慮していないためで、 MayerのEffective diameter factorを使って、内部導体半径を 1.2e-3*0.939 とすると、101.6 pF/m という非常に良い値になります。

撚線のような形状について厳密な解析解を求めることはできず、近似解でさえかなり困難な仕事になりますが、有限要素法だと簡単に扱うことができて、例えば、次のようになります。撚線の外形のプログラミングの方法はいくつか考えられますが、この方法が一番すっきりします。n1の値を変えることで任意の素線数に対応できるようにプログラムしていることに注意してください。


/* coaxial cable: 8D-2V */
a1:=0.4;	/* inner conductor raius */
a2:=3.9;	/* outer conductor radius */
eps:=2.3;	/* dielectric constant */
n1:=6;		/* number of strands on outer layer */

/* outer conductor */
border(1,0,2*pi,60) { x:=a2*cos(t); y:=a2*sin(t) };
/* inner conductor */
r1:=a1/sin(pi/n1);
border(2,0,2*pi,80) {
 t1:=t-fmod(t,2*pi/n1);
 t3:=t1+pi/n1;
 t4:=t3+(t-t3)*(n1/2+1);
 x1:=r1*cos(-t3); y1:=r1*sin(-t3);
 x:=x1+a1*cos(-t4); y:=y1+a1*sin(-t4);
};
buildmesh(5000);

v1:=0.0; v2:=1.0;
solve(v) {
 onbdy(1) v = v1;
 onbdy(2) v = v2;
 pde(v) -laplace(v) = 0;
};
plot(v);

u=0.5*(dx(v)*dx(v)+dy(v)*dy(v));	/* electro static energy for air */
u0:=intt[u];
u1=eps*u0;	/* electro static energy for eps */
vr:=sqrt(u0/u1);	/* velocity ratio */
c:=u1*17.707;	/* pF/m */
z:=1/(2.998e-4*vr*c);	/* characteristic impedance */
printf("C = %.1f (pF/m)\n", c);
printf("Vr = %.3f\n", vr);
printf("Z0 = %.1f Ohm\n", z);

結果は 102.3 pF/m で、Mayer補正した円筒の計算結果と 0.7 % の差になります。JISでは 100+-4 になっていますが、この中心値も少し小さすぎだと思います。

2.4. 横巻外部導体

外部導体が横巻の場合も解析的には解けませんが、有限要素法だとまともに計算することができます。例えば、2422 の構造なら、次のようになります。


/* 2422 (1/0.20*1.3*?/0.1*2.0) */
a1:=0.1;	/* inner conductor raius */
a2:=0.05;	/* outer conductor radius */
eps:=2.3;	/* dielectric constant */

/* outer conductor */
r1:=0.65+a2; n1:=floor(pi/asin(a2/r1)+0.5); r1:=a2/sin(pi/n1);
printf("r1 = %g, r1-a2 =%g\n", r1, r1-a2);
border(1,0,2*pi,300) {
 t1:=t-fmod(t,2*pi/n1);
 t3:=t1+pi/n1;
 t4:=pi+t3-(t-t3)*(n1/2-1);
 x1:=r1*cos(t3); y1:=r1*sin(t3);
 x:=x1+a2*cos(t4); y:=y1+a2*sin(t4);
};
/* inner conductor */
border(2,0,2*pi,30) { x:=a1*cos(-t); y:=a1*sin(-t) };
buildmesh(5000);

v1:=0.0; v2:=1.0;
solve(v) {
 onbdy(1) v = v1;
 onbdy(2) v = v2;
 pde(v) -laplace(v) = 0;
};
plot(v);

u=0.5*(dx(v)*dx(v)+dy(v)*dy(v));	/* electro static energy for air */
u0:=intt[u];
u1=eps*u0;	/* electro static energy for eps */
vr:=sqrt(u0/u1);	/* velocity ratio */
c:=u1*17.707;	/* pF/m */
z:=1/(2.998e-4*vr*c);	/* characteristic impedance */
printf("C = %.1f (pF/m)\n", c);
printf("Vr = %.3f\n", vr);
printf("Z0 = %.1f Ohm\n", z);

結果は 67.6 pF/m になりますが、これは絶縁体径に等しい円筒外部導体として計算したときと同じです。つまり、横巻外部導体の場合は、導体素線の内周部に接する円筒導体に置き換えてもほとんど同じ値が得られ、Meyersの補正係数 (2)のような細工が要らないことを意味します。

2.5. ２重絶縁

発泡ポリエチレンの上に機械的保護を目的とした「スキン」(skin)をかぶせた構造を考えます。例えば、次のようになります。FreeFEMは材料定数を指定するための領域(region)番号を自動的に割り振りますが、予期しない割り当てになることがありますので、三角分解の結果を見てからintt()[]の積分領域を書くのが無難です。


/* coaxial cable - Celler PE with PE skin */
a1:=1.0; a2:=3.0; a3:=2.99;
e1:=2.3; e2:= 1.5;

border(1,0,2*pi,200) { x:=a2*cos(t); y:=a2*sin(t) };
border(0,0,2*pi,200) { x:=a3*cos(t); y:=a3*sin(t) };
border(2,0,2*pi,100) { x:=a1*cos(-t); y:=a1*sin(-t) };
buildmesh(5000);

e=one(region==0)*e1 + one(region>0)*e2;
v1:=0.0; v2:=1.0;
solve(v) {
 onbdy(1) v = v1;
 onbdy(2) v = v2;
 pde(v) -laplace(v)*e = 0;
};
plot(v);

/* electro static energy */
u=0.5*(dx(v)*dx(v)+dy(v)*dy(v));
u1:=e1*intt(0)[u]+e2*(intt(1)[u]+intt(2)[u]);
c:=u1*17.707;	/* pF/m */

/* electro static energy for vacuum */
solve(v) {
 onbdy(1) v = v1;
 onbdy(2) v = v2;
 pde(v) -laplace(v) = 0;
};
nowait;
plot(v);

u=0.5*(dx(v)*dx(v)+dy(v)*dy(v));
u0:=intt[u];
vr:=sqrt(u0/u1);	/* velocity ratio */
c:=u1*17.707;	/* pF/m */
z:=1/(2.99e-4*vr*c);	/* characteristic impedance */

printf("C = %.1f pF\n", c);
printf("Vr = %.3f\n", vr);
printf("Z0 = %.1f Ohm\n", z);

誘電率が場所によって違うときは、誘電体の堺に於ける「境界条件」(bounday condition)

  1) E (電界)の接線成分は常に連続
  2) D (電束)の法線成分は界面にある電荷に対応する量だけ不連続である

を考慮しなければなりませんが、この場合は電荷が存在しないため、D の法線成分も連続になります。解析的な計算ではこの２種類の境界条件を意識的に組み込む必要がありますが、有限要素法では、これが自動的に満足され、取扱が極めて簡単です。

もとろん、エネルギ計算の積分では場所による誘電率の違いを考慮しなければなりませんし、誘電率の分布が変わると、電位分布も変わりますから、速度係数と特性インピーダンスを求めるには、実際の電位分布と誘電体を真空で置き換えたときの電位分布を別々に解く必要があります。

この結果をスキンのない場合と比べると、この程度のスキンはごくわずかな影響しか与えないことがわかります。また、絶縁体外周のスキンの厚さを多少厚くしても、特性はほとんど変わりません。特性に大きな影響を与えるのは、導体近くの絶縁体です。

なお、この例のように、異種の同心円筒誘電体を被覆した同軸ケーブルのキャパシタンスは２つの同軸ケーブルの直列接続として、解析的に解け、この場合は 1.9 % の誤差になります。誤差を減らしたい場合はスキン部分のメッシュを増やすために、

  border(0,...) {...};

の構文をスキン層の内部に追加して、スキン層の要素数を増やします。

3. ２芯シールド線

3.1. 円形シールド

現実的てはありませんが、最初に完全な円で構成した２芯シールド線を考えます。例えば、次のようになります。絶縁体がなければ解析的な近似式が得られますが、この程度の構造でも解析的な計算はかなり面倒になります。


/* shielded twin cable (rigid) */
a1:=0.27*0.939;	/* conductor radius */
a2:=0.68;	/* insulation radius */
e1:=2.3;	/* dielectric constant of insulation */
e2:=1.0;	/* dielectric constant of air gap */

x1:=a2;
a3:=a2*2;
/* shield */
border(1,0,2*pi,100) { x:=a3*cos(t); y:=a3*sin(t) };
/* insulation */
border(0,0,2*pi,60) { x:=a2*cos(t)+x1; y:=a2*sin(t) };
border(0,0,2*pi,60) { x:=a2*cos(t)-x1; y:=a2*sin(t) };
/* conductor */
border(2,0,2*pi,50) { x:=a1*cos(-t)+x1; y:=a1*sin(-t) };
border(3,0,2*pi,50) { x:=a1*cos(-t)-x1; y:=a1*sin(-t) };
buildmesh(5000,1);

v1:=0.0; v2:=0.5;
e=e1*one(region>0) + e2*one(region==0);
solve(v) {
 onbdy(1) v = v1;
 onbdy(2) v = v2;
 onbdy(3) v = -v2;
 pde(v) -laplace(v)*e = 0;
};
plot(v);

/* electro static energy */
u=0.5*(dx(v)*dx(v)+dy(v)*dy(v));
u1:=e2*intt(0)[u]+e1*(intt(1)[u]+intt(2)[u]+intt(3)[u]+intt(4)[u]);
c:=u1*17.707; /* pF/m */

solve(v) {
 onbdy(1) v = v1;
 onbdy(2) v = v2;
 onbdy(3) v = -v2;
 pde(v) -laplace(v) = 0;
};
plot(v);

u=0.5*(dx(v)*dx(v)+dy(v)*dy(v));
u0:=intt[u];
vr:=sqrt(u0/u1); /* velocity ratio */
c:=u1 * 17.707; /* pF/m */
z:=1/(2.99e-4*vr*c);     /* characteristic impedance */

printf("C = %.1f pF/m\n", c);
printf("Vr = %.3f\n", vr);
printf("Z0 = %.1f Ohm\n", z);

幾何学的には 64.3 % 発泡分の空隙がありますが、得られた速度係数 0.729 に対応する等価誘電率は 1.88 (1/(0.729^2)) ですから、28 % 発泡に相当し、この空隙は誘電率を下げるのに有効とは言えません。理由は空隙部分が電界密度の低い導体から離れた場所にあるためです。

実際の製品では、シールドの隙間や導電率の差から、この計算とは、かなり違った値が得られますが、その理由は、この後の「ドレインワイヤ」の項で説明します。

3.2. 最小周長構造のシールド

パイプで作ったもの以外はシールドが円形断面になりませんから、編組、横巻、テープによるシールドといった、もっと実用的な構造を考えると、次のようになります。これらの構造では、周長が最小になることに注意してください。形状決定には「最小作用原理」が働きます。


/* shielded twin cable (flexible) */
a1:=0.27*0.939;	/* conductor radius */
a2:=0.68;	/* insulation radius */
e1:=2.3;	/* dielectric constant of insulation */
e2:=1.0;	/* dielectric constant of air gap */

x1:=a2;
a3:=a2*1.005;
/* shield */
border(1,0,2*pi,100) {
 if (t<pi/2) then { t1:=2*(t-pi/4); x:=a3*cos(t1)+a3; y:=a3*sin(t1) }
 else if (t<=pi) then { t1:=a3-(t-pi/2)*4*a3/pi; x:=t1; y:=a3 }
 else if (t<3*pi/2) then { t1:=2*t-3*pi/2; x:=a3*cos(t1)-a3; y:=a3*sin(t1) }
 else { t1:=(t-3*pi/2)*4*a3/pi-a3; x:=t1; y:=-a3 }
};
/* insulation */
border(0,0,2*pi,60) { x:=a2*cos(t)+x1; y:=a2*sin(t) };
border(0,0,2*pi,60) { x:=a2*cos(t)-x1; y:=a2*sin(t) };
/* conductor */
border(2,0,2*pi,50) { x:=a1*cos(-t)+x1; y:=a1*sin(-t) };
border(3,0,2*pi,50) { x:=a1*cos(-t)-x1; y:=a1*sin(-t) };
buildmesh(5000,1);

v1:=0.0; v2:=0.5; v3:=-0.5;
e=e1*one(region>0) + e2*one(region==0);
solve(v) {
 onbdy(1) v = v1;
 onbdy(2) v = v2;
 onbdy(3) v = v3;
 pde(v) -laplace(v)*e = 0;
};
plot(v);

/* electro static energy for eps */
u=0.5*(dx(v)*dx(v)+dy(v)*dy(v));
u1:=e2*intt(0)[u]+e1*(intt(1)[u]+intt(2)[u]+intt(3)[u]+intt(4)[u]);
c:=u1 * 17.707; /* pF/m */

solve(v) {
 onbdy(1) v = v1;
 onbdy(2) v = v2;
 onbdy(3) v = v3;
 pde(v) -laplace(v) = 0;
};
plot(v);

u=0.5*(dx(v)*dx(v)+dy(v)*dy(v));
u0:=intt[u];
vr:=sqrt(u0/u1); /* velocity ratio */
z:=1/(2.99e-4*vr*c);     /* characteristic impedance */

printf("C = %.1f pF\n", c);
printf("Vr = %.3f\n", vr);
printf("Z0 = %.1f\n", z);

この例でも、空隙は 50 % ありますが、速度係数から求めた等価誘電率は 2.17 ですから、8 % 発泡相当になって、実効誘電率にはあまり効かないことがわかります。

前記の結果は、シールドの電位を 0 にして、２つの導体に絶対値が等しく符号が反対の電位を与えた計算ですから、「平衡伝送」(ballanced circuit transmition) で使うという条件での「実効キャパシタンス」(effective capacitance) を求めたことになります。

「不平衡伝送」(unballanced circuit transmition) の場合は境界条件を、

 v1:=0.0; v2:=0.0; v3:=1.0;

と書き換えて、導体の１つをシールドと同じ 0 V、他の導体を 1 V にします。等電位線の変化をよく観察してください。

平衡伝送と不平衡伝送では特性インピーダンスも違ってきますが、ケーブルに接続される回路の不平衡度を問題にする場合は、次のパラメータもよく使われます。

  Normal Mode 特性インピーダンス - 平衡伝送の特性インピーダンス
  Common Mode 特性インピーダンス - 不平衡伝送の特性インピーダンス

コモンモード特性インピーダンスは、２つの導体を短絡し、それらとシールドの間で不平衡伝送を行ったときの特性インピーダンスで、意図的にはこういう使い方はしませんが、伝送回路のアンバランスによって生まれる、コモンモード・ノイズがこのモードで伝送されますので、この特性の評価が必要になります。

3.3. 部分容量

利用者の使いかた分からない場合、シールドを含めた導体が３つ以上になると、すべての導体間の「部分容量」(partial capacitance)を求めなければなりません。導体数と部分容量数には次の関係があり、かなりの数になりますが、簡単な構造では幾何学的対称性により少なくなることがあります。

  部分容量の数 = n!/(2*(n-2)!

  ここに n = 導体数

!は階乗です。

不規則な構造で、導体数が増えると部分容量の数も加速度的に増えて、収拾がつかなくなり、規格値の表示としては、特定の導体とシールドを含めた他の導体全てとのキャパシタンスを求めて参考値とするといったやりかたになったりします。

導体 1 と導体 2 の部分容量を求める場合は、次のようにします。

  C12 = (Ca + Cb - Cc)/2
  C10 = Ca - C12
  C20 = Cb - C12

  ここに
    C12 = 導体 1 と導体 2 の部分容量
    C10 = 導体 1 とシールドの部分容量
    C20 = 導体 2 とシールドの部分容量
    Ca = 導体 1 とシールドを含む他の全ての導体とのキャパシタンス
    Cb = 導体 2 とシールドを含む他の全ての導体とのキャパシタンス
    Cc = 短絡した導体 1 と導体 2 と、他の全ての導体とのキャパシタンス

つまり、Ca, Cb, Cc を測定しては C12, C10, C20 を計算することを繰り返せば、全ての部分容量が求まります。

２芯シールド線の場合は、対称性を利用して、次のように計算するほうが無駄がありません。

  C10 = C20 = Cc/2
  C12 = Cn - C10/2

  ここに
    C12 = 導体 1 と導体 2 の部分容量
    C10 = 導体 1 とシールドの部分容量
    C20 = 導体 2 とシールドの部分容量
    Cc = コモンモードのキャパシタンス
	(短絡した導体 1 と導体 2 とシールドのキャパシタンス)
    Cn = ノーマルモードのキャパシタンス
	(導体 1 と導体 2 のキャパシタンス)

この方式による計算例は次のとおりです。


/* shielded twin cable */
a1:=0.27*0.939;	/* conductor radius */
a2:=0.68;	/* insulation radius */
e1:=2.3;	/* dielectric constant of insulation */
e2:=1.0;	/* dielectric constant of air gap */

a3:=a2*1.005;
/* shield */
border(1,0,2*pi,100) {
 if (t<pi/2) then { t1:=2*(t-pi/4); x:=a3*cos(t1)+a3; y:=a3*sin(t1) }
 else if (t<=pi) then { t1:=a3-(t-pi/2)*4*a3/pi; x:=t1; y:=a3 }
 else if (t<3*pi/2) then { t1:=2*t-3*pi/2; x:=a3*cos(t1)-a3; y:=a3*sin(t1) }
 else { t1:=(t-3*pi/2)*4*a3/pi-a3; x:=t1; y:=-a3 }
};
/* insulation */
x1:=a2;
border(0,0,2*pi,60) { x:=a2*cos(t)+x1; y:=a2*sin(t) };
border(0,0,2*pi,60) { x:=a2*cos(t)-x1; y:=a2*sin(t) };
/* conductor */
border(2,0,2*pi,50) { x:=a1*cos(-t)+x1; y:=a1*sin(-t) };
border(3,0,2*pi,50) { x:=a1*cos(-t)-x1; y:=a1*sin(-t) };
buildmesh(5000,1);

/* normal mode */
v1:=0.0; v2:=0.5; v3:=-0.5;
e=e1*one(region>0) + e2*one(region==0);
solve(v) {
 onbdy(1) v = v1;
 onbdy(2) v = v2;
 onbdy(3) v = v3;
 pde(v) -laplace(v)*e = 0;
};
plot(v);

/* electro static energy */
u=0.5*(dx(v)*dx(v)+dy(v)*dy(v));
u1:=e2*intt(0)[u]+e1*(intt(1)[u]+intt(2)[u]+intt(3)[u]+intt(4)[u]);
c1:=u1 * 17.707; /* pF/m */

/* electro static energy with vacuum */
solve(v) {
 onbdy(1) v = v1;
 onbdy(2) v = v2;
 onbdy(3) v = v3;
 pde(v) -laplace(v) = 0;
};
plot(v);

u=0.5*(dx(v)*dx(v)+dy(v)*dy(v));
u01:=intt[u];
vr1:=sqrt(u01/u1); /* velocity ratio */

/* common mode */
v1:=0.0; v2:=1.0; v3:=1.0;
solve(v) {
 onbdy(1) v = v1;
 onbdy(2) v = v2;
 onbdy(3) v = v3;
 pde(v) -laplace(v)*e = 0;
};
plot(v);

/* electro static energy */
u=0.5*(dx(v)*dx(v)+dy(v)*dy(v));
u2:=e2*intt(0)[u]+e1*(intt(1)[u]+intt(2)[u]+intt(3)[u]+intt(4)[u]);
c2:=u2 * 17.707; /* pF/m */

/* electro static energy with vacuum */
solve(v) {
 onbdy(1) v = v1;
 onbdy(2) v = v2;
 onbdy(3) v = v3;
 pde(v) -laplace(v) = 0;
};
plot(v);

u=0.5*(dx(v)*dx(v)+dy(v)*dy(v));
u02:=intt[u];
vr2:=sqrt(u02/u2); /* velocity ratio */

c10:=c2/2;
c12:=c1 - c10/2;

printf("Partial Capacitance:\n");
printf(" C10 = C20 = %.1f pF/m\n", c10);
printf(" C12 = %.1f pF/m\n", c12);

z1:=1/(2.99e-4*vr1*c1);     /* normal mode characteristic impedance */
z2:=1/(2.99e-4*vr2*c2);     /* common mode characteristic impedance */
printf("Normal Mode:\n");
printf(" Characteristic Impedance = %.1f\n", z1);
printf(" Velocity ratio = %.3f\n", vr1);
printf("Common Mode:\n");
printf(" Characteristic Impedance = %.1f\n", z2);
printf(" Velocity ratio = %.3f\n", vr2);

3.4. ドレイン・ワイヤ

シールドの端末処理を容易にするためのドレイン・ワイヤ(drain wire)を入れることがありますが、その影響を解析する例です。要素が増えて、region 番号の割当が変わりますので、誘電率の指定と積分に注意が必要です。


/* shielded twin cable with drain wire */
a1:=0.27*0.939;	/* conductor radius */
a2:=0.68;	/* insulation radius */
a3:=0.27;	/* drain radius */
e1:=2.3;	/* dielectric constant of insulation */
e2:=1.0;	/* dielectric constant of air gap */

k1:=1.01; xk:=a2*k1; yk:=0; rk:=-a2*k1; xl:=0;
yl:=sqrt(2*a2*a3+a3*a3)*k1; rl:=-a3*k1;
rlk:=rl-rk; ylk:=yl-yk; xlk:=xl-xk;
w1:=xlk*xlk+ylk*ylk; w2:=sqrt(w1-rlk*rlk); w1:=1/w1;
a9:=(-rlk*xlk-ylk*w2)*w1;
b9:=(-rlk*ylk+xlk*w2)*w1;
c9:=-(rk+a9*xk+b9*yk);
w1:=1/(1+a9*xk/c9+b9*yk/c9)/c9; w2:=rk*rk*w1; x1:=xk-a9*w2; y1:=yk-b9*w2;
w1:=1/(1+a9*xl/c9+b9*yl/c9)/c9; w2:=rl*rl*w1; x2:=xl-a9*w2; y2:=yl-b9*w2;
ta:=atan2(y1,x1); tb:=atan2(y2,x2);
/* shield */
border(1,0,2*pi,100) {
 if (t<=pi/4+ta) then {
  t1:=t-pi/4; w1:=xk*cos(t1); p1:=w1+sqrt(w1*w1-xk*xk+rk*rk);
 }
 else if (t<pi/4+tb) then {
  t1:=t-pi/4; p1:=-c9/(a9*cos(t1)+b9*sin(t1));
 }
 else if (t<=5*pi/4-tb) then {
  t1:=t-pi/4; w1:=yl*sin(t1); p1:=w1+sqrt(w1*w1-yl*yl+rl*rl);
 }
 else if (t<5*pi/4-ta) then {
  t1:=t-pi/4; p1:=-c9/(-a9*cos(t1)+b9*sin(t1));
 }
 else if (t<=3*pi/2) then {
  t1:=t-pi/4; w1:=xk*cos(t1); p1:=-w1+sqrt(w1*w1-xk*xk+rk*rk);
 }
 else {
  t1:=t-pi/4; p1:=rk/sin(t1);
 };
 x:=p1*cos(t1); y:=p1*sin(t1)
};
/* drain */
y1:=sqrt(2*a2*a3+a3*a3);
border(4,0,2*pi,50) { x:=a3*cos(-t); y:=a3*sin(-t)+y1 };
/* insulation */
x1:=a2;
border(0,0,2*pi,60) { x:=a2*cos(t)+x1; y:=a2*sin(t) };
border(0,0,2*pi,60) { x:=a2*cos(t)-x1; y:=a2*sin(t) };
/* conductor */
border(2,0,2*pi,50) { x:=a1*cos(-t)+x1; y:=a1*sin(-t) };
border(3,0,2*pi,50) { x:=a1*cos(-t)-x1; y:=a1*sin(-t) };
buildmesh(5000,1);

/* normal mode */
v1:=0.0; v2:=0.5; v3:=-0.5; v4:=0.0;
e=e1*one(region>1) + e2*one(region<=1);
solve(v) {
 onbdy(1) v = v1;
 onbdy(2) v = v2;
 onbdy(3) v = v3;
 onbdy(4) v = v4;
 pde(v) -laplace(v)*e = 0;
};
plot(v);

/* electro static energy */
u=0.5*(dx(v)*dx(v)+dy(v)*dy(v));
u1:=e2*(intt(0)[u]+intt(1)[u])+e1*(intt(2)[u]+intt(3)[u]+intt(4)[u]+intt(5)[u]);
c1:=u1 * 17.707; /* pF/m */

/* electro static energy with vacuum */
solve(v) {
 onbdy(1) v = v1;
 onbdy(2) v = v2;
 onbdy(3) v = v3;
 onbdy(4) v = v4;
 pde(v) -laplace(v) = 0;
};
plot(v);

u=0.5*(dx(v)*dx(v)+dy(v)*dy(v));
u01:=intt[u];
vr1:=sqrt(u01/u1); /* velocity ratio */

/* common mode */
v1:=0.0; v2:=1.0; v3:=1.0; v4:=0.0;
solve(v) {
 onbdy(1) v = v1;
 onbdy(2) v = v2;
 onbdy(3) v = v3;
 onbdy(4) v = v4;
 pde(v) -laplace(v)*e = 0;
};
plot(v);

/* electro static energy */
u=0.5*(dx(v)*dx(v)+dy(v)*dy(v));
u2:=e2*(intt(0)[u]+intt(1)[u])+e1*(intt(2)[u]+intt(3)[u]+intt(4)[u]+intt(5)[u]);
c2:=u2 * 17.707; /* pF/m */

/* electro static energy with vacuum */
solve(v) {
 onbdy(1) v = v1;
 onbdy(2) v = v2;
 onbdy(3) v = v3;
 onbdy(4) v = v4;
 pde(v) -laplace(v) = 0;
};
plot(v);

u=0.5*(dx(v)*dx(v)+dy(v)*dy(v));
u02:=intt[u];
vr2:=sqrt(u02/u2); /* velocity ratio */

c10:=c2/2;
c12:=c1 - c10/2;

printf("Partial Capacitance:\n");
printf(" C10 = C20 = %.1f pF/m\n", c10);
printf(" C12 = %.1f pF/m\n", c12);

z1:=1/(2.99e-4*vr1*c1);     /* normal mode characteristic impedance */
z2:=1/(2.99e-4*vr2*c2);     /* common mode characteristic impedance */
printf("Normal Mode:\n");
printf(" Characteristic Impedance = %.1f\n", z1);
printf(" Velocity ratio = %.3f\n", vr1);
printf("Common Mode:\n");
printf(" Characteristic Impedance = %.1f\n", z2);
printf(" Velocity ratio = %.3f\n", vr2);

この程度のドレインワイヤだと、キャパシタンスで 4 % の増加になります。特性インピーダンスはキャパシタンスの平方根に反比例しますから、ほぼこの 1/2 の増加になります。速度係数は空隙が多少増えるため、1 % 程度の減少になります。

実際の製品では、

編組や横巻の隙間が原因で、キャパシタンスが増える (15 % 程度が多い)
シールドが完全導体でないため、磁界がシールドの外まで広がり、インダクタンスが増える (35 % 程度のケースが多い)
その結果、速度係数、特性インピーダンスも大きくなる。

ため、例えば、キャパシタンスで 15 % 程度減少、特性インピーダンスで 40 % 程度増加することがあります。

3.5. モノフィラメントのフィラー

モノフィラメント(monofilament)のフィラー(filler - 介在)を使った２芯シールド線を解析します。プログラムは、後で３芯以上のケーブルでも n1 を書き換えるだけで、同じ形状記述が使えるように工夫しています。


/* shielded twin cable with monofilament filler */
a1:=0.31;	/* conductor radius */
a2:=0.8;	/* insulation radius */
e1:=2.3;	/* dielectric constant of insulation */
e2:=2.3;	/* dielectric constant of filler */
e3:=1.0;	/* dielectric constant of air gap */

n1:=2;		/* number of conductors */
r1:=a2/sin(pi/n1);	/* pitch diameter */
/* outer filler */
w1:=1+1/sin(pi/n1)-1/tan(pi/n1); a3:=w1*w1/(1+w1)/2*a2; r3:=r1+a2-a3;
b1:=1;	/* boundary number */
/* overall shield */
k1:=1.01; xk:=r1*k1; yk:=0; rk:=-a2*k1;
t1:=pi/n1; xl:=r3*cos(t1)*k1; yl:=r3*sin(t1)*k1; rl:=-a3*k1;
rlk:=rl-rk; ylk:=yl-yk; xlk:=xl-xk;
w1:=xlk*xlk+ylk*ylk; w2:=sqrt(w1-rlk*rlk); w1:=1/w1;
a9:=(-rlk*xlk-ylk*w2)*w1;
b9:=(-rlk*ylk+xlk*w2)*w1;
c9:=-(rk+a9*xk+b9*yk);
/* a9*x + b9*y + c9 = 0 <- common tangent to insulation and filler */
w1:=1/(1+a9*xk/c9+b9*yk/c9)/c9; w2:=rk*rk*w1; x1:=xk-a9*w2; y1:=yk-b9*w2;
w1:=1/(1+a9*xl/c9+b9*yl/c9)/c9; w2:=rl*rl*w1; x2:=xl-a9*w2; y2:=yl-b9*w2;
/* (x1,y1), (x2,y2) <- inter section between common tangent and two circle */
ta:=atan2(y1,x1); tb:=atan2(y2,x2); tc:=2*pi/n1-tb; td:=2*pi/n1-ta;
t0:=2*pi/n1; t1:=acos(-a9); t2:=2*pi/n1-t1;
border(b1++,0,2*pi,60) {
 t3:=floor(t/t0)*t0; t4:=fmod(t,t0);
 x1:=r1*cos(t3); y1:=r1*sin(t3);
 t5:=t3+pi/n1;
 x2:=r3*cos(t5); y2:=r3*sin(t5);
 p2:=t-ta;
 if(t4<2*ta) then {
  w1:=x1*cos(p2)+y1*sin(p2); p1:=w1+sqrt(w1*w1-x1*x1-y1*y1+a2*a2);
 }
 else if (t4<=tb+ta) then {
  w1:=t1+t3; p1:=c9/cos(p2-w1);
 }
 else if (t4<tc+ta) then {
  w1:=x2*cos(p2)+y2*sin(p2);
  p1:=w1+sqrt(w1*w1-x2*x2-y2*y2+a3*a3);
 }
 else {
  w1:=t2+t3; p1:=c9/cos(p2-w1);
 };
 p1:=p1*1.005;
 x:=p1*cos(p2); y:=p1*sin(p2);
};
for (i1:=0;i1<n1;i1++) {
 t3:=2*pi/n1*i1;
 /* filler */
 t4:=t3+pi/n1;
 x1:=r3*cos(t4); y1:=r3*sin(t4); 
 border(0,0,2*pi,20) { x:=a3*cos(t)+x1; y:=a3*sin(t)+y1 };
 /* insulation */
 x1:=r1*cos(t3); y1:=r1*sin(t3);
 border(0,0,2*pi,20) { x:=a2*cos(t)+x1; y:=a2*sin(t)+y1 };
 /* conductor */
 border(b1++,0,2*pi,20) { x:=a1*cos(-t)+x1; y:=a1*sin(-t)+y1 };
};
buildmesh(5000,1);

/* normal mode */
v1:=0.0; v2:=0.5; v3:=-0.5; v4:=0.0;
e=e1*one(region>2) + e2*one((0<region)and(region<3)) + e3*one(region==0);
solve(v) {
 onbdy(1) v = v1;
 onbdy(2) v = v2;
 onbdy(3) v = v3;
 onbdy(4) v = v4;
 pde(v) -laplace(v)*e = 0;
};
plot(v);

/* electro static energy */
u=0.5*(dx(v)*dx(v)+dy(v)*dy(v));
u1:=e3*intt(0)[u]+e2*(intt(1)[u]+intt(4)[u])
 +e1*(intt(2)[u]+intt(3)[u]+intt(5)[u]+intt(6)[u]);
c1:=u1 * 17.707; /* pF/m */

/* electro static energy with vacuum */
solve(v) {
 onbdy(1) v = v1;
 onbdy(2) v = v2;
 onbdy(3) v = v3;
 pde(v) -laplace(v) = 0;
};
plot(v);

u=0.5*(dx(v)*dx(v)+dy(v)*dy(v));
u01:=intt[u];
vr1:=sqrt(u01/u1); /* velocity ratio */

/* common mode */
v1:=0.0; v2:=1.0; v3:=1.0; v4:=0.0;
solve(v) {
 onbdy(1) v = v1;
 onbdy(2) v = v2;
 onbdy(3) v = v3;
 onbdy(4) v = v4;
 pde(v) -laplace(v)*e = 0;
};
plot(v);

/* electro static energy */
u=0.5*(dx(v)*dx(v)+dy(v)*dy(v));
u2:=e3*intt(0)[u]+e2*(intt(1)[u]+intt(4)[u])
 +e1*(intt(2)[u]+intt(3)[u]+intt(5)[u]+intt(6)[u]);
c2:=u2 * 17.707; /* pF/m */

/* electro static energy with vacuum */
solve(v) {
 onbdy(1) v = v1;
 onbdy(2) v = v2;
 onbdy(3) v = v3;
 onbdy(4) v = v4;
 pde(v) -laplace(v) = 0;
};
plot(v);

u=0.5*(dx(v)*dx(v)+dy(v)*dy(v));
u02:=intt[u];
vr2:=sqrt(u02/u2); /* velocity ratio */

c10:=c2/2;
c12:=c1 - c10/2;

printf("Partial Capacitance:\n");
printf(" C10 = C20 = %.1f pF/m\n", c10);
printf(" C12 = %.1f pF/m\n", c12);

z1:=1/(2.99e-4*vr1*c1);     /* normal mode characteristic impedance */
z2:=1/(2.99e-4*vr2*c2);     /* common mode characteristic impedance */
printf("Normal Mode:\n");
printf(" Characteristic Impedance = %.1f\n", z1);
printf(" Velocity ratio = %.3f\n", vr1);
printf("Common Mode:\n");
printf(" Characteristic Impedance = %.1f\n", z2);
printf(" Velocity ratio = %.3f\n", vr2);

モノフィラメントの材質を発泡ポリエチレンに変えても、特性インピーダンスはほとんど変わりません。かけるコストの割に効果が少ないのは、モノフィラメントの位置が電位傾度の大きい導体から離れているためです。逆に、絶縁体を発泡ポリエチレンに変えると、特性インピーダンスは 22 % も増えます。

4. ３芯シールド・ケーブル

３芯シールド・ケーブルの場合は対称度が高いため、部分容量は導体間と導体シールド間の２つしかなく、解析も簡単です。例えば、

  c1 = 全ての導体を短絡したものとシールド間のキャパシタンス
  c2 = 一つの導体と残りの２導体とシールドの間のキャパシタンス

を求めて、

  k10 = k20 = k30 = c1 / 3
  k12 = k23 = k31 = (c2 - c1)/2

  ここに、
    kab = a-b 間の部分容量 (1, 2, 3 は導体、0 はシールド)

と計算することができます。


/* shielded 3 core cable with monofilament filler */
a1:=0.31;	/* conductor radius */
a2:=0.8;	/* insulation radius */
e1:=2.3;	/* dielectric constant of insulation */
e2:=1.0;	/* dielectric constant of air gap */

n1:=3;		/* number of conductors */
r1:=a2/sin(pi/n1);	/* pitch diameter */
b1:=1;	/* boundary number */
/* overall shield */
k1:=1.01; xk:=r1*k1; yk:=0; rk:=-a2*k1;
t1:=2*pi/n1; xl:=r1*cos(t1)*k1; yl:=r1*sin(t1)*k1; rl:=-a2*k1;
rlk:=rl-rk; ylk:=yl-yk; xlk:=xl-xk;
w1:=xlk*xlk+ylk*ylk; w2:=sqrt(w1-rlk*rlk); w1:=1/w1;
a9:=(-rlk*xlk-ylk*w2)*w1;
b9:=(-rlk*ylk+xlk*w2)*w1;
c9:=-(rk+a9*xk+b9*yk);
/* a9*x + b9*y + c9 = 0 <- common tangent to insulation and filler */
w1:=1/(1+a9*xk/c9+b9*yk/c9)/c9; w2:=rk*rk*w1; x1:=xk-a9*w2; y1:=yk-b9*w2;
w1:=1/(1+a9*xl/c9+b9*yl/c9)/c9; w2:=rl*rl*w1; x2:=xl-a9*w2; y2:=yl-b9*w2;
/* (x1,y1), (x2,y2) <- inter section between common tangent and two circle */
ta:=atan2(y1,x1); tb:=atan2(y2,x2); tc:=2*pi/n1-tb; td:=2*pi/n1-ta;
t0:=2*pi/n1; t1:=acos(-a9); t2:=2*pi/n1-t1;
border(b1++,0,2*pi,60) {
 t3:=floor(t/t0)*t0; t4:=fmod(t,t0);
 x1:=r1*cos(t3); y1:=r1*sin(t3);
 p2:=t-ta;
 if(t4<=2*ta) then {
  w1:=x1*cos(p2)+y1*sin(p2); p1:=w1+sqrt(w1*w1-x1*x1-y1*y1+a2*a2);
 }
 else {
  w1:=t2+t3; p1:=c9/cos(p2-w1);
 };
 p1:=p1*1.005;
 x:=p1*cos(p2); y:=p1*sin(p2);
};
for (i1:=0;i1<n1;i1++) {
 t3:=2*pi/n1*i1;
 /* insulation */
 x1:=r1*cos(t3); y1:=r1*sin(t3);
 border(0,0,2*pi,20) { x:=a2*cos(t)+x1; y:=a2*sin(t)+y1 };
 /* conductor */
 border(b1++,0,2*pi,20) { x:=a1*cos(-t)+x1; y:=a1*sin(-t)+y1 };
};
buildmesh(5000,1);

/* k10+k20+k30 */
v1:=0.0; v2:=1.0;
e=e1*one(region>0) + e2*one(region==0);
solve(v) {
 onbdy(1) v = v1;
 onbdy(2) v = v2;
 onbdy(3) v = v2;
 onbdy(4) v = v2;
 pde(v) -laplace(v)*e = 0;
};
plot(v);
/* electro static energy */
u=0.5*(dx(v)*dx(v)+dy(v)*dy(v));
u1:= e1*(intt(1)[u]+intt(2)[u]+intt(3)[u]+intt(4)[u]+intt(5)[u]+intt(6)[u])
	+e2*intt(0)[u];
c1:=u1 * 17.707; /* pF/m */

/* k10+k12+k13 */
v1:=0.0; v2:=1.0;
solve(v) {
 onbdy(1) v = v1;
 onbdy(2) v = v2;
 onbdy(3) v = v1;
 onbdy(4) v = v1;
 pde(v) -laplace(v)*e = 0;
};
plot(v);
/* electro static energy */
u=0.5*(dx(v)*dx(v)+dy(v)*dy(v));
u2:= e1*(intt(1)[u]+intt(2)[u]+intt(3)[u]+intt(4)[u]+intt(5)[u]+intt(6)[u])
	+e2*intt(0)[u];
c2:=u2 * 17.707; /* pF/m */

k10:=c1/3;
k12:=(c2-k10)/2;

printf("Partial Capacitance:\n");
printf(" k10 = k20 = k30 = %.1f pF/m\n", k10);
printf(" k12 = k23 = k31 = %.1f pF/m\n", k12);

5. ４芯ケーブル

４芯ケーブルの場合は、３種類の部分容量ができますが、

  c1 = 全ての導体とシールド間のキャパシタンス
  c2 = 対角線になる導体２本と残りのすべての間のキャパシタンス
  c3 = １つの導体と残りの全ての間のキャパシタンス

から

  k10 = c1/4
  k12 = (c2 - 2*k10)/4
  k13 = c3 - k10 - 2*k12

で求めるのが簡単です。


/* shielded quad cable with monofilament filler */
a1:=0.31;	/* conductor radius */
a2:=0.8;	/* insulation radius */
e1:=2.3;	/* dielectric constant of insulation */
e2:=1.0;	/* dielectric constant of air gap */

n1:=4;		/* number of conductors */
r1:=a2/sin(pi/n1);	/* pitch diameter */
b1:=1;	/* boundary number */
/* overall shield */
k1:=1.01; xk:=r1*k1; yk:=0; rk:=-a2*k1;
t1:=2*pi/n1; xl:=r1*cos(t1)*k1; yl:=r1*sin(t1)*k1; rl:=-a2*k1;
rlk:=rl-rk; ylk:=yl-yk; xlk:=xl-xk;
w1:=xlk*xlk+ylk*ylk; w2:=sqrt(w1-rlk*rlk); w1:=1/w1;
a9:=(-rlk*xlk-ylk*w2)*w1;
b9:=(-rlk*ylk+xlk*w2)*w1;
c9:=-(rk+a9*xk+b9*yk);
/* a9*x + b9*y + c9 = 0 <- common tangent to insulation and filler */
w1:=1/(1+a9*xk/c9+b9*yk/c9)/c9; w2:=rk*rk*w1; x1:=xk-a9*w2; y1:=yk-b9*w2;
w1:=1/(1+a9*xl/c9+b9*yl/c9)/c9; w2:=rl*rl*w1; x2:=xl-a9*w2; y2:=yl-b9*w2;
/* (x1,y1), (x2,y2) <- inter section between common tangent and two circle */
ta:=atan2(y1,x1); tb:=atan2(y2,x2); tc:=2*pi/n1-tb; td:=2*pi/n1-ta;
t0:=2*pi/n1; t1:=acos(-a9); t2:=2*pi/n1-t1;
border(b1++,0,2*pi,60) {
 t3:=floor(t/t0)*t0; t4:=fmod(t,t0);
 x1:=r1*cos(t3); y1:=r1*sin(t3);
 p2:=t-ta;
 if(t4<=2*ta) then {
  w1:=x1*cos(p2)+y1*sin(p2); p1:=w1+sqrt(w1*w1-x1*x1-y1*y1+a2*a2);
 }
 else {
  w1:=t2+t3; p1:=c9/cos(p2-w1);
 };
 p1:=p1*1.005;
 x:=p1*cos(p2); y:=p1*sin(p2);
};
for (i1:=0;i1<n1;i1++) {
 t3:=2*pi/n1*i1;
 /* insulation */
 x1:=r1*cos(t3); y1:=r1*sin(t3);
 border(0,0,2*pi,20) { x:=a2*cos(t)+x1; y:=a2*sin(t)+y1 };
 /* conductor */
 border(b1++,0,2*pi,20) { x:=a1*cos(-t)+x1; y:=a1*sin(-t)+y1 };
};
buildmesh(5000,1);

/* k10+k20+k30+k40 */
v1:=0.0; v2:=1.0;
e=e1*one(region>0) + e2*one(region==0);
solve(v) {
 onbdy(1) v = v1;
 onbdy(2) v = v2;
 onbdy(3) v = v2;
 onbdy(4) v = v2;
 onbdy(5) v = v2;
 pde(v) -laplace(v)*e = 0;
};
plot(v);
/* electro static energy */
u=0.5*(dx(v)*dx(v)+dy(v)*dy(v));
u1:= e1*(intt(1)[u]+intt(2)[u]+intt(3)[u]+intt(5)[u]+intt(6)[u]
	+intt(7)[u]+intt(8)[u]) + e2*intt(0)[u];
c1:=u1 * 17.707; /* pF/m */

/* 2*k10+4*k12 */
v1:=0.0; v2:=1.0;
solve(v) {
 onbdy(1) v = v1;
 onbdy(2) v = v2;
 onbdy(3) v = v1;
 onbdy(4) v = v2;
 onbdy(5) v = v1;
 pde(v) -laplace(v)*e = 0;
};
plot(v);
/* electro static energy */
u=0.5*(dx(v)*dx(v)+dy(v)*dy(v));
u2:= e1*(intt(1)[u]+intt(2)[u]+intt(3)[u]+intt(5)[u]+intt(6)[u]
	+intt(7)[u]+intt(8)[u]) + e2*intt(0)[u];
c2:=u2 * 17.707; /* pF/m */

/* k10+2*k12+k13 */
v1:=0.0; v2:=1.0;
solve(v) {
 onbdy(1) v = v1;
 onbdy(2) v = v2;
 onbdy(3) v = v1;
 onbdy(4) v = v1;
 onbdy(5) v = v1;
 pde(v) -laplace(v)*e = 0;
};
plot(v);
/* electro static energy */
u=0.5*(dx(v)*dx(v)+dy(v)*dy(v));
u3:= e1*(intt(1)[u]+intt(2)[u]+intt(3)[u]+intt(5)[u]+intt(6)[u]
	+intt(7)[u]+intt(8)[u]) + e2*intt(0)[u];
c3:=u3 * 17.707; /* pF/m */

k10:=c1/4;
k12:=(c2-2*k10)/4;
k13:=c3-k10-2*k12;

printf("Partial Capacitance:\n");
printf(" k10 = k20 = k30 = k40= %.1f pF/m\n", k10);
printf(" k12 = k23 = k31 = k41= %.1f pF/m\n", k12);
printf(" k13 = k24 = %.1f pF/m\n", k13);

この他、２芯ケーブルと同じようにして、カッド(quad)の対角線に位置する導体で平衡伝送するときの特性等も計算できます。

6. 多芯ケーブル

シールドされた７対ケーブルの中心に置かれた対の特性を求めてみます。


/* shielded 7 pair cable */
a1:=0.191;	/* conductor raius */
a2:=0.475;	/* insulation radius */
a3:=0.95;
a4:=2*a2*1.01;
a5:=6*a2*1.02;
eps:=3.1;	/* dielectric constant */

b1:=1;		/* boundary number */

/* overall shield */
border(b1++,0,2*pi,80) { x:=a5*cos(t); y:=a5*sin(t) };

/* pairs */
t2:=pi/4;
x2:=a2*cos(t2); y2:=a2*sin(t2);
/* border(0,0,2*pi,20) { x:=a4*cos(t); y:=a4*sin(t) }; */
border(0,0,2*pi,20) { x:=a2*cos(t)+x2; y:=a2*sin(t)+y2 };
border(b1++,0,2*pi,20) { x:=a1*cos(-t)+x2; y:=a1*sin(-t)+y2 };
border(0,0,2*pi,20) { x:=a2*cos(t)-x2; y:=a2*sin(t)-y2 };
border(b1++,0,2*pi,20) { x:=a1*cos(-t)-x2; y:=a1*sin(-t)-y2 };

r1:=1.01*a4/sin(pi/6);
for(i1:=0;i1<6;i1++) {
 t1:=2*pi/6*i1;
 x1:=r1*cos(t1); y1:=r1*sin(t1);
 x2:=a2*cos(t2); y2:=a2*sin(t2);
/* border(0,0,2*pi,20) { x:=a4*cos(t)+x1; y:=a4*sin(t)+y1 }; */
 border(0,0,2*pi,20) { x:=a2*cos(t)+x1+x2; y:=a2*sin(t)+y1+y2 };
 border(b1++,0,2*pi,20) { x:=a1*cos(-t)+x1+x2; y:=a1*sin(-t)+y1+y2 };
 border(0,0,2*pi,20) { x:=a2*cos(t)+x1-x2; y:=a2*sin(t)+y1-y2 };
 border(b1++,0,2*pi,20) { x:=a1*cos(-t)+x1-x2; y:=a1*sin(-t)+y1-y2 };
};

buildmesh(5000);

v0:=0; v1:=0.5; v2:=-0.5;
e=1*one(region<1) + eps*one(region>0);
solve(v) {
 onbdy(1) v = v0;
 onbdy(2) v = v1; onbdy(3) v = v2;
 onbdy(4) v = v0; onbdy(5) v = v0;
 onbdy(6) v = v0; onbdy(7) v = v0;
 onbdy(8) v = v0; onbdy(9) v = v0;
 onbdy(10) v = v0; onbdy(11) v = v0;
 onbdy(12) v = v0; onbdy(13) v = v0;
 onbdy(14) v = v0; onbdy(15) v = v0;
 pde(v) -laplace(v)*e = 0;
};
plot(v);
u=0.5*(dx(v)*dx(v)+dy(v)*dy(v));
u0:=intt(0)[u];		/* enegy at air gap */
u1=1*u0+eps*(intt[u]-u0);	/* energy at insulation */

solve(v) {
 onbdy(1) v = v0;
 onbdy(2) v = v1; onbdy(3) v = v2;
 onbdy(4) v = v0; onbdy(5) v = v0;
 onbdy(6) v = v0; onbdy(7) v = v0;
 onbdy(8) v = v0; onbdy(9) v = v0;
 onbdy(10) v = v0; onbdy(11) v = v0;
 onbdy(12) v = v0; onbdy(13) v = v0;
 onbdy(14) v = v0; onbdy(15) v = v0;
 pde(v) -laplace(v) = 0;
};
plot(v);
u=0.5*(dx(v)*dx(v)+dy(v)*dy(v));
u0:=intt[u];		/* enegy at air gap */

v:=sqrt(u0/u1);	/* velocity ratio */
c:=u1 * 17.707;	/* pF/m */
z:=1/(2.998e-4*v*c);	/* characteristic impedance */
printf("v = %.2f\n", v);
printf("c = %.1f pF\n", c);
printf("Z0 = %.1f Ohm\n", z);

対の位置を入れ換えて、外周の対の特性を求めると、次のようになって、特性インピーダンスが 5 % 程度下るのがわかります。近くにシールドという大きな導体があるためです。


/* shielded 7 pair cable */
a1:=0.191;	/* conductor raius */
a2:=0.475;	/* insulation radius */
a3:=0.95;
a4:=2*a2*1.01;
a5:=6*a2*1.02;
eps:=3.1;	/* dielectric constant */

b1:=1;		/* boundary number */

/* overall shield */
border(b1++,0,2*pi,80) { x:=a5*cos(t); y:=a5*sin(t) };

/* pairs */
t2:=pi/4;

r1:=1.01*a4/sin(pi/6);
for(i1:=0;i1<6;i1++) {
 t1:=2*pi/6*i1;
 x1:=r1*cos(t1); y1:=r1*sin(t1);
 x2:=a2*cos(t2); y2:=a2*sin(t2);
/* border(0,0,2*pi,20) { x:=a4*cos(t)+x1; y:=a4*sin(t)+y1 }; */
 border(0,0,2*pi,20) { x:=a2*cos(t)+x1+x2; y:=a2*sin(t)+y1+y2 };
 border(b1++,0,2*pi,20) { x:=a1*cos(-t)+x1+x2; y:=a1*sin(-t)+y1+y2 };
 border(0,0,2*pi,20) { x:=a2*cos(t)+x1-x2; y:=a2*sin(t)+y1-y2 };
 border(b1++,0,2*pi,20) { x:=a1*cos(-t)+x1-x2; y:=a1*sin(-t)+y1-y2 };
};

x2:=a2*cos(t2); y2:=a2*sin(t2);
/* border(0,0,2*pi,20) { x:=a4*cos(t); y:=a4*sin(t) }; */
border(0,0,2*pi,20) { x:=a2*cos(t)+x2; y:=a2*sin(t)+y2 };
border(b1++,0,2*pi,20) { x:=a1*cos(-t)+x2; y:=a1*sin(-t)+y2 };
border(0,0,2*pi,20) { x:=a2*cos(t)-x2; y:=a2*sin(t)-y2 };
border(b1++,0,2*pi,20) { x:=a1*cos(-t)-x2; y:=a1*sin(-t)-y2 };

buildmesh(5000);

v0:=0; v1:=0.5; v2:=-0.5;
e=1*one(region<1) + eps*one(region>0);
solve(v) {
 onbdy(1) v = v0;
 onbdy(2) v = v1; onbdy(3) v = v2;
 onbdy(4) v = v0; onbdy(5) v = v0;
 onbdy(6) v = v0; onbdy(7) v = v0;
 onbdy(8) v = v0; onbdy(9) v = v0;
 onbdy(10) v = v0; onbdy(11) v = v0;
 onbdy(12) v = v0; onbdy(13) v = v0;
 onbdy(14) v = v0; onbdy(15) v = v0;
 pde(v) -laplace(v)*e = 0;
};
plot(v);
u=0.5*(dx(v)*dx(v)+dy(v)*dy(v));
u0:=intt(0)[u];		/* enegy at air gap */
u1=1*u0+eps*(intt[u]-u0);	/* energy at insulation */

solve(v) {
 onbdy(1) v = v0;
 onbdy(2) v = v1; onbdy(3) v = v2;
 onbdy(4) v = v0; onbdy(5) v = v0;
 onbdy(6) v = v0; onbdy(7) v = v0;
 onbdy(8) v = v0; onbdy(9) v = v0;
 onbdy(10) v = v0; onbdy(11) v = v0;
 onbdy(12) v = v0; onbdy(13) v = v0;
 onbdy(14) v = v0; onbdy(15) v = v0;
 pde(v) -laplace(v) = 0;
};
plot(v);
u=0.5*(dx(v)*dx(v)+dy(v)*dy(v));
u0:=intt[u];		/* enegy at air gap */

v:=sqrt(u0/u1);	/* velocity ratio */
c:=u1 * 17.707;	/* pF/m */
z:=1/(2.998e-4*v*c);	/* characteristic impedance */
printf("v = %.2f\n", v);
printf("c = %.1f pF\n", c);
printf("v = %.1f Ohm\n", z);

7. 対称性の利用

有限要素法を含む一般の数値解析では、メッシュを細かくすると精度が上がりますが、計算量とメモリ消費量は加速度的 (3)に増加しますから、解析すべき領域が対称な場合は、

  1) 自然境界条件 (対称性からポテンシャルの法線成分が 0)
  2) ディリクレ境界条件 (対称性からポテンシャルが決まる)

場所を見付けて、繰り返しの最小単位だけを解析するようにすると精度と時間が稼げます。

このコースでは分かりやすいように、ほとんどの例で、断面全体を解析対称にしていますが、精度と計算速度をあげるには、可能な限り対称性を利用して、解析範囲をしぼらなければなりません。

例えば、前記の２芯のシールドケーブルなら、下記のようにすれば 1/4 の領域で済みます。エネルギの計算では 4 倍しなければならないことを忘れないでください。


/* shielded twin cable (1/4) */
a1:=0.27*0.939;	/* conductor radius */
a2:=0.68;	/* insulation radius */
e1:=2.3;	/* dielectric constant of insulation */
e2:=1.0;	/* dielectric constant of air gap */

x1:=a2;
border(2,0,pi,40) { t1:=pi-t; x:=a1*cos(t1)+a2; y:=a1*sin(t1) };
border(4,0,pi/2,50) { x:=a2*cos(t)+a2; y:=a2*sin(t) };
border(1,0,a2-a1,20) { x:=t; y:=0 };
border(3,a2+a1,2*a2,20) { x:=t; y:=0 };
border(5,0,a2,50) { x:=a2-t; y:=a2 };
border(6,0,a2,50) { x:=0; y:=a2-t };
border(0,pi/2,pi,50) { x:=a2*cos(t)+a2; y:=a2*sin(t) };
buildmesh(5000,1);

v1:=0.0; v2:=0.5;
e=e1*one(region>0) + e2*one(region==0);
solve(v) {
 onbdy(1) dnu(v) = 0;
 onbdy(2) v = v2;
 onbdy(3) dnu(v) = 0;
 onbdy(4) v = v1;
 onbdy(5) v = v1;
 onbdy(6) v = v1;
 pde(v) -laplace(v)*e = 0;
};
plot(v);

/* electro static energy for eps */
u=0.5*(dx(v)*dx(v)+dy(v)*dy(v));
u1:=e2*(intt(4)[u]+intt(5)[u])
 +e1*(intt(1)[u]+intt(2)[u]+intt(3)[u]+intt(6)[u]);
c:=4*u1*17.707; /* pF/m */

solve(v) {
 onbdy(1) dnu(v) = 0;
 onbdy(2) v = v2;
 onbdy(3) dnu(v) = 0;
 onbdy(4) v = v1;
 onbdy(5) v = v1;
 onbdy(6) v = v1;
 pde(v) -laplace(v) = 0;
};
plot(v);

u=0.5*(dx(v)*dx(v)+dy(v)*dy(v));
u0:=intt[u];
vr:=sqrt(u0/u1); /* velocity ratio */
z:=1/(2.99e-4*vr*c);     /* characteristic impedance */

printf("C = %.1f pF\n", c);
printf("Vr = %.3f\n", vr);
printf("Z0 = %.1f\n", z);

７芯撚内部導体、円筒外部導体の同軸ケーブルなら、わずか 1/12 の領域を解析すれば十分で、例えば、次のようになります。

/* coaxial cable: 8D-2V (1/12) */
a1:=0.4;	/* inner conductor raius */
a2:=3.9;	/* outer conductor radius */
eps:=2.3;	/* dielectric constant */
n1:=6;		/* number of strands on outer layer */

t1:=pi/n1; r1:=a1/sin(t1);
/* outer conductor */
border(1,0,t1,60) { x:=a2*cos(t); y:=a2*sin(t) };
/* inner conductor */
t2:=pi/2+t1;
border(2,0,t2,20) { t3:=t2-t; x:=a1*cos(t3)+r1; y:=a1*sin(t3); };
border(3,r1+a1,a2,20) { x:=t; y:=0 };
r2:=r1*cos(pi/6);
border(4,r2,a2,20) { p1:=a2-t+r2; x:=p1*cos(t1); y:=p1*sin(t1) };
buildmesh(5000);

v1:=0.0; v2:=1.0;
solve(v) {
 onbdy(1) v = v1;
 onbdy(2) v = v2;
 onbdy(3) dnu(v) = 0;
 onbdy(4) dnu(v) = 0;
 pde(v) -laplace(v) = 0;
};
plot(v);

u=0.5*(dx(v)*dx(v)+dy(v)*dy(v));	/* electro static energy for air */
u0:=intt[u];
u1=eps*u0;	/* electro static energy for eps */
vr:=sqrt(u0/u1);	/* velocity ratio */
c:=2*n1*u1*17.707;	/* pF/m */
z:=1/(2.998e-4*vr*c);	/* characteristic impedance */
printf("C = %.1f (pF/m)\n", c);
printf("Vr = %.3f\n", vr);
printf("Z0 = %.1f Ohm\n", z);

8. 開領域問題

シールドなしの対(twin-lead, twisted pair)とか、テレビの平衡フィーダのように宇宙全体に電界が広がってしまうような構造を数値解析で解くには、通常、次の２つの方法が使われます。

ポテンシャルが 0 とみなせるような十分離れた場所にシールドを置く
電位が解析的に求められる境界を使う

要素数が増えるのが欠点ですが、プログラミングは前者が簡単です。問題は人為的なシールドをどの程度離れたところに置くかですが、平衡ケーブルの場合は空間のエネルギ密度が導体群からの距離の自乗に反比例して小さくなりますので、数倍から１０倍離れた場所にすれば、実用的な精度が得られます。また、シールの距離をいくつか変えて計算した結果から、距離が無限大になったときの値を外挿するといった方法もあります。 (4)

後者の方法は、電位の計算を、

  V = -Σ(q/(2*π*ε)*log(r/R))
  ここに、
	V = 空間の電位 (V)
	q = 導体の線電荷密度 (Q/m)
	ε = 空間の誘電率 (F/m)
	r = 導体からの距離 (m)
	R = 電位 0 の基準点からの距離 (m)

で求めますが、ラプラス方程式を解かないと導体上の電荷密度もわかりませんから、境界上の電位分布の近似値を初期値にして、ラプラス方程式を解き、そこから得られる電荷分布を元に、境界上の電位分布を修正するといった方法で、正しい値に近付きます。

こうして、導体を含む部分空間の電位分布がわかれば、下記のいずれかの方法で、キャパシタンスの計算ができます。

エネルギを使う方法
下記のように、全空間のエネルギと導体の電位差からキャパシタンスを求めます。
```
  U = ∫(ε*2)*((dV/dx)^2+(dV/dy)^2)*ds - ∫(ε*2)*V*(dv/dn)*dl
  C = U/(2*V^2)
  ここに、
	U = 全空間のエネルギ (J)
	ε = 空間の誘電率 (F/m)
	∫ds = 全領域を対象にした面積分
	dv/dn = 境界の法線方向微分 (電界の法線方向成分) (V/m)
	∫dl = 境界上の１周積分
	C = キャパシタンス (F/m)
	V = 導体間の電位差 (V)
```
この第２項の積分 ∫ε*(dv/dn)*dl は境界内部のすべてを導体で置き換えたときの、導体の表面電荷密度になっていることに注意してください。この電荷密度を q とすれば、よく知られたエネルギ計算式 U = Q*V/2 と一致します。つまり、境界領域外部空間のエネルギを領域上の電荷に換算しているわけです。

電荷を使う方法
下記のように、導体の電荷と電位差からキャパシタンスを求めます。

  Q = ∫ε*(dv/dn)*dl
  C = Q / V
  ここに、
	Q = 導体上の電荷 (Q)
	ε = 導体表面に接する空間の誘電率 (F/m)
	dv/dn = 境界の法線方向微分 (電界の法線方向成分) (V/m)
	∫dl = 導体表面上の１周積分
	C = キャパシタンス (F/m)
	V = 導体間の電位差 (V)

8.1. VHF フィーダ

JIS C 3330 の1VHFを解析すると次のようになります。


/* VHF feeder code (JIS C 3330 1VHF) */
pi2:=2*pi; a1:=0.345; a2:=0.75; s1:=(9 - 1.5)/2; h1:=1/2; a3:=20;
eps:=2.3;
y1:=h1; x1:=sqrt(a2*a2 - h1*h1);
t1:=atan(y1/x1);
x1:=s1-x1;
border(1,0,pi2,60) { x:=a3*cos(t); y:=a3*sin(t) };
border(0,0,4*pi-4*t1+6,100) {
 if (t<=2*pi-2*t1) then {
	p1:=t+t1; x:=a2*cos(p1)-s1; y:=a2*sin(p1)
 }
 else if (t<=2*pi-2*t1+3) then {
	p1:=t-(2*pi-2*t1);
	x:=-x1+2*x1/3*p1; y:=-h1
 }
 else if (t<=4*pi-4*t1+3) then {
	p1:=t-(pi2-2*t1+3)+t1-pi;
	x:=a2*cos(p1)+s1; y:=a2*sin(p1)
 }
 else {
	p1:=t-(4*pi-4*t1+3);
	x:=x1-2*x1/3*p1; y:=h1
 }
};
border(2,0,pi2,20) { x:=a1*cos(-t)+s1; y:=a1*sin(-t) };
border(3,0,pi2,20) { x:=a1*cos(-t)-s1; y:=a1*sin(-t) };

buildmesh(10000);

v1:=0; v2:=0.5;
e=one(region>0)*eps+one(region==0)*1;
solve(v) {
 onbdy(1) v = v1;
 onbdy(2) v = v2;
 onbdy(3) v = -v2;
 pde(v) -laplace(v)*e = 0;
};
plot(v);
u=0.5*(dx(v)*dx(v)+dy(v)*dy(v));	/* electro static energy for air */
u1:=intt(0)[u]+eps*(intt(1)[u]+intt(2)[u]+intt(3)[u]+intt(4)[u]);

solve(v) {
 onbdy(1) v = v1;
 onbdy(2) v = v2;
 onbdy(3) v = -v2;
 pde(v) -laplace(v) = 0;
};
plot(v);
u=0.5*(dx(v)*dx(v)+dy(v)*dy(v));	/* electro static energy for air */
u0:=intt[u];

c:=u1*17.707;	/* pF/m */
vr:=sqrt(u0/u1);	/* velocity ratio */
z:=1/(2.998e-4*vr*c);   /* characteristic impedance */
printf("C = %.1f (pF/m)\n", c);
printf("Vr = %.3f\n", vr);
printf("Z0 = %.1f Ohm\n", z);

8.2. UHF フィーダ (1)

JIS C 3330のUHF-Mフィーダの解析例です。


/* UHF feeder code (JIS C 3330 UHF-M) */
pi2:=2*pi; a1:=0.435; a2:=1.8; a3:=2.75; a4:=20;
eps1:=1.5; eps2:=2.3;
h1:=sqrt(a3*a3-a2*a2);
t1:=atan(h1/a2);
border(1,0,pi2,60) { x:=a4*cos(t); y:=a4*sin(t) };
border(0,0,2,40) {
 if (t<=1) then {
	p1:=t*2*(pi-t1)+t1; x:=a3*cos(p1)-a2; y:=a3*sin(p1)
 }
 else {
	p1:=(t-1)*2*(pi-t1)-(pi-t1); x:=a3*cos(p1)+a2; y:=a3*sin(p1)
 }
};
border(0,0,pi2,60) { x:=a2*cos(t)+a2; y:=a2*sin(t) };
border(0,0,pi2,60) { x:=a2*cos(t)-a2; y:=a2*sin(t) };
border(2,0,pi2,20) { x:=a1*cos(-t)+a2; y:=a1*sin(-t) };
border(3,0,pi2,20) { x:=a1*cos(-t)-a2; y:=a1*sin(-t) };

buildmesh(10000,1);

v1:=0; v2:=0.5;
e=one(region>0)*eps+one(region==0)*1;
solve(v) {
 onbdy(1) v = v1;
 onbdy(2) v = v2;
 onbdy(3) v = -v2;
 pde(v) -laplace(v)*e = 0;
};
plot(v);

u=0.5*(dx(v)*dx(v)+dy(v)*dy(v));	/* electro static energy for air */
u1:=intt(0)[u]+eps2*intt(1)[u]
 +eps1*(intt(2)[u]+intt(3)[u]+intt(4)[u]+intt(5)[u]);

solve(v) {
 onbdy(1) v = v1;
 onbdy(2) v = v2;
 onbdy(3) v = -v2;
 pde(v) -laplace(v) = 0;
};
plot(v);
u=0.5*(dx(v)*dx(v)+dy(v)*dy(v));	/* electro static energy for air */
u0:=intt[u];

c:=u1*17.707;	/* pF/m */
vr:=sqrt(u0/u1);	/* velocity ratio */
z:=1/(2.998e-4*vr*c);	/* characteristic impedance */
printf("C = %.1f (pF/m)\n", c);
printf("Vr = %.3f\n", vr);
printf("Z0 = %.1f Ohm\n", z);

8.3. UHF フィーダ (2)

JIS C 3330のもうひとつの UHF フィーダの解析例です。

/* UHF feeder code (JIS C 3330 UHF) */
pi2:=2*pi; a1:=0.435; a2:=1.8; a3:=2.75; a4:=20;
eps1:=1.5; eps2:=2.3;
t1:=atan(a3/a2);
border(1,0,pi2,80) { x:=a4*cos(t); y:=a4*sin(t) };
border(0,-t1,pi2-t1,80) {
 if (t<=t1) then {
	w1:=a2*cos(t);p1:=w1+sqrt(w1*w1-a2*a2+a3*a3);
	x:=p1*cos(t); y:=p1*sin(t)
 }
 else if (t<=pi-t1) then {
	p1:=a3/sin(t);
	x:=p1*cos(t); y:=p1*sin(t)
 }
 else if (t<=pi+t1) then {
	w1:=a2*cos(t); p1:=-w1+sqrt(w1*w1-a2*a2+a3*a3);
	x:=p1*cos(t); y:=p1*sin(t)
 }
 else {
	p1:=-a3/sin(t);
	x:=p1*cos(t); y:=p1*sin(t)
 }
};
border(0,0,pi2,60) { x:=a2*cos(t)+a2; y:=a2*sin(t) };
border(0,0,pi2,60) { x:=a2*cos(t)-a2; y:=a2*sin(t) };
border(2,0,pi2,20) { x:=a1*cos(-t)+a2; y:=a1*sin(-t) };
border(3,0,pi2,20) { x:=a1*cos(-t)-a2; y:=a1*sin(-t) };

buildmesh(10000,1);

v1:=0; v2:=0.5;
e=one(region>0)*3+one(region==0)*1;
solve(v) {
 onbdy(1) v = v1;
 onbdy(2) v = v2;
 onbdy(3) v = -v2;
 pde(v) -laplace(v)*e = 0;
};
plot(v);

u=0.5*(dx(v)*dx(v)+dy(v)*dy(v));	/* electro static energy for air */
u1:=intt(0)[u]+eps2*intt(1)[u]
 +eps1*(intt(2)[u]+intt(3)[u]+intt(4)[u]+intt(5)[u]);

solve(v) {
 onbdy(1) v = v1;
 onbdy(2) v = v2;
 onbdy(3) v = -v2;
 pde(v) -laplace(v) = 0;
};
plot(v);

u=0.5*(dx(v)*dx(v)+dy(v)*dy(v));	/* electro static energy for air */
u0:=intt[u];

c:=u1*17.707;	/* pF/m */
vr:=sqrt(u0/u1);	/* velocity ratio */
z:=1/(2.998e-4*vr*c);	/* characteristic impedance */
printf("C = %.1f (pF/m)\n", c);
printf("Vr = %.3f\n", vr);
printf("Z0 = %.1f Ohm\n", z);

8.4. 開領域を直接取り扱う方法

要素数を減らしたい場合や解析領域全体のポテンシャルの精度が必要な場合は、境界領域の電位を導体上の電荷から解析的に計算します。これは、全ての導体を含む境界上のポテンシャル(電位)が導体上の電荷密度と導体までの距離で決まるという条件を利用するのですが、 FreeFEMで解く場合は、解析すべき領域を含む境界のディリクレ条件として、

  V = 1/(2*pi*eps)*(q1*log(r1) + q2*log(r2) + .. )
  ここに、
    q1, q2, .. = 線電荷密度 (q/m)
    eps = 誘電率 (F/m)
    r1, r2, .. = 境界上の点 (x,y) から線電荷までの距離 (m)

を使って、逐次近似法で解くことになります。

例えば、往復２導体(対)の中心座標をそれぞれ (-a/2,0), (a/2,0)、線電荷密度をそれぞれ -q, q とすると、

  V = q/(4*pi*eps)*log(((x+a)^2+y^2)/((x-a)^2+y^2));

になります。

境界が対から十分はなれていて、a<<sqrt(x^2+y^2) が成り立つ場合は、

  V = q/(4*pi*eps)*a*x*/(x^2+y^2)

と近似できます。

現実の導体は線とちがって広がりを持ちますから、導体中心に線電荷を仮定した計算では誤差が出ますが、これは、電気映像法で、導体半径 a、中心間距離 d の場合、中心間距離 2*sqrt(d*d-a*a) に置かれた線電荷がつくる電界と同じになることが分かっていて、補正も簡単です。ただ、一般には、補正なしで計算しても結果はほとんど変わりません。

静電エネルギは無限遠の空間に広がっていますが、その一部しか考えていませんので、導体の表面電荷を使う方法が考えやすいかもしれませんが、境界外部のエネルギは、境界上の情報だけから、下記のように計算できますから、これに領域内部の静電エネルギを加えれば、全空間の静電エネルギが求まります。

  -∫(ε/2)*En*ds
  ここに、
	ε = 空間の誘電率 (F/m)
	En = 境界の電界の法線方向成分 (V/m)
	∫*ds = 境界上の周積分

この方針で裸導線のペアを解析すると、例えば、次のとおりになります。電荷として適当な初期値を仮定してLaplace方程式を解き、その結果から電荷を求めて初期値にフィードバックするというサイクルを繰り返し、誤差が小さくなったら終了します。収束過程を見るため毎回表示を止めていますが、止めたくない場合は、nowait命令を入れます。


/* open boundary - bare wire pair */
a1:=1;		/* conductor radius */
d1:=2.5;	/* (distance between conductor)/2 */
x0:=6; y0:=3;	/* boundary size */
/* bundary */
border(1,-x0,x0,60) { x:=t; y:=-y0 };
border(1,-y0,y0,60) { x:=x0; y:=t };
border(1,-x0,x0,40) { x:=-t; y:=y0 };
border(1,-y0,y0,40) { x:=-x0; y:=-t };
/* conductor */
border(2,0,2*pi,60) { x:=a1*cos(-t)+d1; y:=a1*sin(-t) };
border(3,0,2*pi,60) { x:=a1*cos(-t)-d1; y:=a1*sin(-t) };
buildmesh(10000,1);

q1:=1e-12;	/* electric charge, initial guess */
e0:=1e7/(4*pi*2.998e8^2);	/* dielectric constant for vaccum */
s1:=sqrt(d1*d1-a1*a1);	/* conductor spacing */
v1:=0; v2:=0.5;
for (q0:=999; abs(q0-q1)/q1>1e-2; ) {
 solve(v) {
  onbdy(1) v = q1/4/pi/e0*log(((x+s1)^2+y^2)/((x-s1)^2+y^2));
  onbdy(2) v = v2;
  onbdy(3) v = -v2;
  pde(v) -laplace(v) = 0;
 };
 plot(v);

 /* charge on conductor */
 u=sqrt(dx(v)*dx(v)+dy(v)*dy(v));
 q0:=q1;
 q1:=e0*int(2)[u];	/* Q/m */
/* printf("C = %.1f (pF/m)\n", q*1e12);*/
/* printf("err: %g\n", abs(q0-q)/q);*/
};
printf("C = %.1f (pF/m by charge\n", q1*1e12);
/* energy */
u=dx(v)*dx(v)+dy(v)*dy(v);
u1:=0.5*e0*intt[u];
q=0.5*e0*v*sqrt((dx(v)*nx)^2+(dy(v)*ny)^2);
u2:=int(1)[q];
c1:=2*(u1-u2)/((v2-(-v2))^2)*1e12;
printf("C = %.1f (pF/m) by energy\n", c1);
/* analytical */
c2:=pi*e0/log((2*d1+sqrt((2*d1)^2-4*a1))/(2*a1));
printf("C = %.1f (pF/m) by analisys\n", c2*1e12);

この構造は解析的に解けて、

  C = 2*π*ε/acosh((d^2-2*a^2)/(2*a^2))  (F/m)
    = π*ε/log((d-sqrt(d^2-4*a))/(2*a))  (F/m)

  ここに、
    a = 導体半径 (m)
    d = 導体の中心間距離 (m)

になりますが、電荷を使う方法でも、エネルギを使う方法でも、誤差 5.8 % になります。

絶縁体導体の対で、複数の誘電体が存在する場合は、数電体境界に分極電荷が現れますから、その電荷も考慮して境界の電位を求めなければなりませんが、 FreeFEMの場合は、プログラムを改造しないと無理ですから、分極電荷を無視できるていどに外部境界を広げるのが良いようです。

9. 静電遮蔽効果

実用的な静電シールドは相互キャパシタンスがあまりに小さくて、数値解析は無理ですが、極めて目の荒いものは解析できて、例えば、次のようになります。


/* elctrostatic shield */
a1:=0.245;	/* inner conductor raius */
a2:=0.75;	/* insulation radius */
a3:=0.06;	/* served conductor radius */
a4:=2;		/* outer conductor (electrode) radius */
eps1:=2.3;	/* dielectric constant inisulation */
eps2:=1;	/* dielectric constant in air */
n1:=5;		/* number of served conductor */

border(1,0,2*pi,40) { x:=a1*cos(-t); y:=a1*sin(-t) };
border(0,0,2*pi,50) { x:=a2*cos(t); y:=a2*sin(t) };
r1:=a2+a3;	/* pitch diameter */
for(i1:=0;i1<n1;i1++) {
 t1:=2*pi/n1*i1; x1:=r1*cos(t1); y1:=r1*sin(t1);
 border(2,0,2*pi,20) { x:=a3*cos(-t)+x1; y:=a3*sin(-t)+y1 };
};
/* outer conductor */
border(3,0,2*pi,100) { x:=a4*cos(t); y:=a4*sin(t) };

buildmesh(5000,1);

e=eps1*one(region<2)+eps2*one(region>1);
v1:=0; v2:=1;
solve(v,1) {
 onbdy(1) v = v1;
 onbdy(2) v = v1;
 onbdy(3) v = v2;
 pde(v) -laplace(v)*e = 0;
};
plot(v);
u=0.5*(dx(v)*dx(v)+dy(v)*dy(v));
u1:=eps1*(intt(0)[u]+intt(1)[u])+eps2*(intt[u]-intt(0)[u]-intt(1)[u]);
c1:=u1*17.707;	/* pF/m */
printf("C1 = %.g (pF/m)\n", c1);

solve(v,-1) {
 onbdy(1) v = v1;
 onbdy(2) v = v2;
 onbdy(3) v = v1;
 pde(v) -laplace(v)*e = 0;
};
plot(v);
u=0.5*(dx(v)*dx(v)+dy(v)*dy(v));
u2:=eps1*(intt(0)[u]+intt(1)[u])+eps2*(intt[u]-intt(0)[u]-intt(1)[u]);
c2:=u2*17.707;	/* pF/m */
printf("C2 = %.g (pF/m)\n", c2);

solve(v,-1) {
 onbdy(1) v = v2;
 onbdy(2) v = v1;
 onbdy(3) v = v1;
 pde(v) -laplace(v)*e = 0;
};
plot(v);
u=0.5*(dx(v)*dx(v)+dy(v)*dy(v));
u3:=eps1*(intt(0)[u]+intt(1)[u])+eps2*(intt[u]-intt(0)[u]-intt(1)[u]);
c3:=u3*17.707;	/* pF/m */
printf("C3 = %.g (pF/m)\n", c3);

k31:=(c1-c2+c3)/2;
k12:=c3-k31;
k23:=c2-k12;
printf("K12 = %.1f (pF/m)\n", k12*e-12);
printf("K23 = %.1f (pF/m)\n", k23*e-12);
printf("K31 = %.1f (pF/m)\n", k31*e-12);

同じ構造で電位だけ変えて解きなおしますので、最初の「LU 分解」の結果を保存しています。

10. あとがき

このコースで解説した手法は、現代の有限要素法の技術としては極めて初歩的ですが、コンピュータによる数値解析の威力に目覚め、たとえこの程度の技術でも知らない人と比べたらまるでレベルが違ってしまうことに気づいてほしいという希望から、かなり丁寧な例題を用意して、技術書には書かれていない、私のノーハウを伝えようと努力しました。これを機会にさらに研讃を積むための意欲と学問と技術に対する感動を生むことを願っています。日常に仕事の中にも、新鮮な工夫がほしいものです。

この教材は 1998-08-26 に準備を始めて、1998-09-13 に仕上げました。

11. 補足

(0) 電荷を利用したキャパシタンス計算

空間エネルギを使わずに、導体の表面電荷からキャパシタンスを計算する場合は、次のようになります。

  C = Q / V
  Q = ∫ρ*ds
   = ∫ε*E*ds
     =∫ε*dV/dn*ds					(2)
  ここに、
	C = キャパシタンス (F/m)
	Q = 導体の全表面電荷 (Q/m)
	V = 導体間電位差 (V)
	ρ = 導体の表面電荷密度 (C/m^2)
	ε = 導体表面の誘電率 (F/m)
	E = 導体の表面の電界の強さ (V/m)
	dV/dn = 導体表面に於ける V の法線方向微分 (V/m)
	∫*ds = 導体表面の面積分(２次元なら周積分)

積分は導体の表面の周積分ですから、 FreeFEMの場合は int()[]を使うことになります。

導体表面の電界は法線方向しかありませんから、表面の法線ベクトルを考えずに、単に電界を積分をするだけで済むことに注意してください。

(1) エネルギ誤差

電気工学に置ける有限要素法の偉人Silvesterの名著

  P.P.Silvester & R.L.Ferrari,- Einit Element for Electrical Engineers
	(Cambridge University Press 1983)

にあるGreenの定理を使った証明(2.05)が明解で、真のポテンシャル解 u(x,y) と、指定された境界条件で正確に 0 になるような微分可能な誤差関数 v(x,y) の組み合わせ u(x,y)+h*v(x,y) (h は誤差を表す小さなスカラ値)で有限要素法の誤差モデルを作ると、エネルギ W(u+h*v) は、

  W(u+h*v) = W(u) + h*h*W(v)

となることを証明している。つまり、エネルギ誤差は h でなく h の自乗に比例するわけで、エネルギ誤差はポテンシャル誤差に比べて十分少ない！

(2) Meyers の実効径

同軸ケーブルの内部導体として撚線を使ったときと、同じキャパシタンスになるような単線の直径がMeyersの「等価外径」で、外接円の直径より小さくなる。

総素線数 1 3 7 12 19 27 37

外層素線数 1 3 6 9 12 15 18

外径=d0=q*d q 1 2.16 3.00 4.16 5.00 6.16 7.00

等価外径=k1*d0 k1 1.000 0.871 0.939 0.957 0.970 0.976 0.980

総素線数	1	3	7	12	19	27	37
外層素線数	1	3	6	9	12	15	18
外径=d0=q*d q	1	2.16	3.00	4.16	5.00	6.16	7.00
等価外径=k1*d0 k1	1.000	0.871	0.939	0.957	0.970	0.976	0.980

(3) 要素(頂点)数と計算量

多くの数値計算の文献で解説されているように、n 元の連立方程式を「LU 分解」 (Gaus decomposition)で解くには、

  UL 分解に約 n^3/3 回
  右辺の計算に約 n^2 回

の乗算が必要になるため、計算すべき頂点数が増えるに従って、計算量が極端に増える。定数項だけが違う場合は右辺だけの計算で済ませるとかなり速くなるため、 FreeFEMではsolve(u, -1)のように既に LU 分解したマトリクスを保存して何度も使える仕組みを用意してある。連立方程式の計算に逆行列を使わず、「LU 分解」(三角分解)を使うのは、これができるため。

(4) 平衡対による空間のエネルギ密度

線電荷密度 q, -q を持つ２本の導体をそれぞれ (a/2,0), (-a/2,0) に配置したときの (x,y) 点の電位 V は、

  V = q/(2*π*ε)
    *log((1+(a/r)^2+0.5*(a/r)*cos(θ))/((1+(a/r)^2-0.5*(a/r)*cos(θ)))
  ここに、
    r = sqrt(x*x+y*y)
    θ = atan2(y,x);

対から離れた空間では a<<r だから log を展開して高次の項を捨てると

  V = q*a/(2*π*ε)*cos(θ)/r = q*a/(2*π*ε)*x/(x*x+y*y)

電界のエネルギ密度 u は、u = (ε/2)*((dV/dx)^2+(dV/dy)^2) を計算して、

  u = (ε/2)/(x^2+y^2)^2 = (ε/2)/r^4

半径 r1 の円から外に存在するエネルギ U は

                             ∞ 2π
  U = (ε/2)*(q/(2*π*ε))^2*∫∫(1/r^4)*dθ*dr
                             r  0
 
    = (q/4)/r1^2

となって、半径の自乗に比例して低下する。

これは私のアイデアによる証明。ケーブル設計に従事していないと、こんなことは気にしないらしいが、こういった計算はReduce等の数式処理プログラムを使うと簡単！

平林浩一, 1998-09-23