1960 年代に航空宇宙工学の構造解析で開発された有限要素法は、その後、偏微分 方程式の数値解析の標準的な方法として確立され、今ではあらゆる技術分野の常識 として、多量の文献がありますが、あまり勉強しなかった人には理解が困難ですか ら、社内で絶えず必要になる問題に絞って、有限要素法による初歩的なケーブル設 計の実務を解説し、ある程度
分からなくても使えることを目的に、そのまま使える豊富な例題を用意することにしました。
有限要素法のプログラムは実にたくさんありますが、ここではFreeFEMとい うフランスの大学で開発され、私がケーブル設計で使いやすいように少し拡張した 版を使うことにします。オリジナルのFreeFEMではprintf()以外に も、for()文やインクリメント等の演算子、fmod()等の関数がありま せんので、複雑なケーブルではかなり長いプログラムになってしまいます。
FreeFEMの全体像と構文の詳細についてはFreeFEMのドキュメントを 参照してください。日本語の解説も用意してあります。
FreeFEMはborder()文で指定された形状の頂点座標を計算し、 Delauney-Boronoi三角分割を行った結果を表示し、solve()文で指定 された偏微分方程式と境界条件を使って、問題を解くようになっています。このコ ースで用意された例題は、三角分割の結果と、偏微分方程式を解いて得られたポテ ンシャル分布の等高線を表示した後、ケーブルの特性を計算して表示するようにな っていますので、グラフィック表示の後は、マウスの左釦をして先に進めてくださ い。もし、マウスの操作を省いて、結果だけを表示したければ、プログラムの始め の方に、
nowait;を追加します。ただ、三角分割と等ポテンシャル線図をチェックしないと、とんで もない間違いを見逃す危険があります。
ケーブルの電気特性で最も基本になるのがキャパシタンス(capacitance)で、これ が計算できると、次のようにして、高周波の特性インピーダンス(characteristic impedance) と速度係数(velocity ratio)が求まります。
速度係数 = c*sqrt(C0/C) 特性インピーダンス = 1/(c*sqrt(C0*C)) ここに、 c = 光速 (2.99792458e8) C = キャパシタンス C0 = 誘電体部分を真空にしたときのキャパシタンスつまり、キャパシタンスの計算ができれば、特性インピーダンスと速度係数は簡単 に分かるのですが、キャパシタンスが解析的に計算できるのは、同軸ケーブル等、 かなり簡単な形状に限られ、一般には解析的な計算はできません。有限要素法とか 差分法といった数値計算を使うことになります。
数値計算の場合は、スカラ・ポテンシャル(電位)を使った、2次元のポアソン方 程式 (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)を実際の境界条件のもとで解いて、2つの導体間の電位差と電界のエネルギから、 キャパシタンスを計算します。
もちろん、導体表面の電界強度と誘電率から、導体の表面電荷を求め、 それを積分した導体の全電荷と導体間電位差から、 注(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最初の項は閉領域と同じ断面の面積分ですが、 第2項は外部境界上の線積分です。
なお、キャパシタンスの計算で電位分布を求めなければならないのは誘電体内部だ けで、電荷が存在する導体の電位をディリクレ(Direcret)型境界条件として指定し ますから、電荷を考慮する必要はなく、
(d/dx+d/dy)(e*(d/dx+d/dy)*V) = 0というラプラス方程式(Laplace's equation)を解くだけの簡単な仕事になり、 FreeFEMでは、laplace()演算子だけで間に合います。
以下、実例をもとに説明しますが、キャパシタンス、速度係数、特性インピーダン スは幾何学的に相似なら同じ値になりますから、寸法は mm 単位にします。m 単位 でFreeFEMのスクリプトを書くと、値が小さくなるため、同じ場所ではない のに同じ場所と見られることがあります。
最も簡単な同軸ケーブルから始めます。下記のFreeFEMスクリプトを引数 にFreeFEMを起動すると、5C-2V(JIS C 3501) のキャパシタンスを 計算します。1e7/(4*pi*2.998e8^2) は真空の誘電率です。この例では、誘電体が 一つしかありませんので、誘電体を真空にしても電位分布は変わらず、微分方程 式の求解も1度で済みます。
/* 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文の4番目の引数でメッシュの細か さが調整できますので、いろいろ試してください。
JIS C 3501ではキャパシタンスが 67+-3 pF/m になっていますが、この 中心値は少し低すぎます。外部導体が編組ですから、空隙によって誘電率は下る はずで、厳密解より小さくなることはあり得ません。
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つの特性が偏芯に対して鈍感なことがわかります。
8D-2Vの場合は内部導体が 7/0.18 という7本撚構造になっていて、その 外接円の半径 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 になっていますが、この中心値も少し小さすぎだと 思います。
外部導体が横巻の場合も解析的には解けませんが、有限要素法だとまともに計算す ることができます。例えば、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)のような細工が要らないことを意味します。
発泡ポリエチレンの上に機械的保護を目的とした「スキン」(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 の法線成分 も連続になります。解析的な計算ではこの2種類の境界条件を意識的に組み込む 必要がありますが、有限要素法では、これが 自動的に満足され、取扱が極めて簡単です。
もとろん、エネルギ計算の積分では場所による誘電率の違いを考慮しなければなり ませんし、誘電率の分布が変わると、電位分布も変わりますから、速度係数と特性 インピーダンスを求めるには、実際の電位分布と誘電体を真空で置き換えたときの 電位分布を別々に解く必要があります。
この結果をスキンのない場合と比べると、この程度のスキンはごくわずかな影響し か与えないことがわかります。 また、絶縁体外周のスキンの厚さを多少厚くしても、 特性はほとんど変わりません。 特性に大きな影響を与えるのは、 導体近くの絶縁体です。
なお、この例のように、異種の同心円筒誘電体を被覆した同軸ケーブルのキャパシ タンスは2つの同軸ケーブルの直列接続として、解析的に解け、この場合は 1.9 % の誤差になります。誤差を減らしたい場合はスキン部分のメッシュを増やすために、
border(0,...) {...};の構文をスキン層の内部に追加して、スキン層の要素数を増やします。
現実的てはありませんが、最初に完全な円で構成した2芯シールド線を考えます。 例えば、次のようになります。絶縁体がなければ解析的な近似式が得られますが、 この程度の構造でも解析的な計算はかなり面倒になります。
/* 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 % 発泡に相当し、こ の空隙は誘電率を下げるのに有効とは言えません。理由は空隙部分が電界密度 の低い導体から離れた場所にあるためです。
実際の製品では、シールドの隙間や導電率の差から、 この計算とは、かなり違った値が得られますが、 その理由は、この後の「ドレインワイヤ」の項で説明します。
パイプで作ったもの以外はシールドが円形断面になりませんから、編組、横巻、テ ープによるシールドといった、もっと実用的な構造を考えると、次のようになりま す。これらの構造では、周長が最小になることに注意してください。形状決定には 「最小作用原理」が働きます。
/* 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 にして、2つの導体に絶対値が等しく符号が 反対の電位を与えた計算ですから、「平衡伝送」(ballanced circuit transmition) で使うという条件での「実効キャパシタンス」(effective capacitance) を求めた ことになります。
「不平衡伝送」(unballanced circuit transmition) の場合は境界条件を、
v1:=0.0; v2:=0.0; v3:=1.0;と書き換えて、導体の1つをシールドと同じ 0 V、他の導体を 1 V にします。等 電位線の変化をよく観察してください。
平衡伝送と不平衡伝送では特性インピーダンスも違ってきますが、ケーブルに接続 される回路の不平衡度を問題にする場合は、次のパラメータもよく使われます。
Normal Mode 特性インピーダンス - 平衡伝送の特性インピーダンス Common Mode 特性インピーダンス - 不平衡伝送の特性インピーダンスコモンモード特性インピーダンスは、2つの導体を短絡し、それらとシールドの間 で不平衡伝送を行ったときの特性インピーダンスで、意図的にはこういう使い方は しませんが、伝送回路のアンバランスによって生まれる、コモンモード・ノイズが このモードで伝送されますので、この特性の評価が必要になります。
利用者の使いかた分からない場合、シールドを含めた導体が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 を計算することを繰り返せば、 全ての部分容量が求まります。
2芯シールド線の場合は、対称性を利用して、次のように計算するほうが無駄があ りません。
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);
シールドの端末処理を容易にするためのドレイン・ワイヤ(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 % 程度の減少にな ります。
実際の製品では、
モノフィラメント(monofilament)のフィラー(filler - 介在)を使った2芯シール ド線を解析します。プログラムは、後で3芯以上のケーブルでも 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 % も増えます。
3芯シールド・ケーブルの場合は対称度が高いため、部分容量は導体間と導体シー ルド間の2つしかなく、解析も簡単です。例えば、
c1 = 全ての導体を短絡したものとシールド間のキャパシタンス c2 = 一つの導体と残りの2導体とシールドの間のキャパシタンスを求めて、
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);
4芯ケーブルの場合は、3種類の部分容量ができますが、
c1 = 全ての導体とシールド間のキャパシタンス c2 = 対角線になる導体2本と残りのすべての間のキャパシタンス c3 = 1つの導体と残りの全ての間のキャパシタンスから
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);この他、2芯ケーブルと同じようにして、カッド(quad)の対角線に位置する導体で 平衡伝送するときの特性等も計算できます。
シールドされた7対ケーブルの中心に置かれた対の特性を求めてみます。
/* 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);
有限要素法を含む一般の数値解析では、メッシュを細かくすると精度が上がります が、計算量とメモリ消費量は加速度的 (3)に増加しますから、解析すべき領域が 対称な場合は、
1) 自然境界条件 (対称性からポテンシャルの法線成分が 0) 2) ディリクレ境界条件 (対称性からポテンシャルが決まる)場所を見付けて、繰り返しの最小単位だけを解析するようにすると精度と時間が稼 げます。
このコースでは分かりやすいように、ほとんどの例で、断面全体を解析対称にして いますが、精度と計算速度をあげるには、可能な限り対称性を利用して、解析範囲 をしぼらなければなりません。
例えば、前記の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);
7芯撚内部導体、円筒外部導体の同軸ケーブルなら、わずか 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);
シールドなしの対(twin-lead, twisted pair)とか、テレビの平衡フィーダのよう に宇宙全体に電界が広がってしまうような構造を数値解析で解くには、通常、次 の2つの方法が使われます。
後者の方法は、電位の計算を、
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 = 境界上の1周積分 C = キャパシタンス (F/m) V = 導体間の電位差 (V)この第2項の積分 ∫ε*(dv/dn)*dl は境界内部のすべてを導体で置き換えたときの、 導体の表面電荷密度になっていることに注意してください。 この電荷密度を q とすれば、 よく知られたエネルギ計算式 U = Q*V/2 と一致します。 つまり、境界領域外部空間のエネルギを領域上の電荷に換算しているわけです。
Q = ∫ε*(dv/dn)*dl C = Q / V ここに、 Q = 導体上の電荷 (Q) ε = 導体表面に接する空間の誘電率 (F/m) dv/dn = 境界の法線方向微分 (電界の法線方向成分) (V/m) ∫dl = 導体表面上の1周積分 C = キャパシタンス (F/m) V = 導体間の電位差 (V)
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);
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);
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);
要素数を減らしたい場合や解析領域全体のポテンシャルの精度が必要な場合は、 境界領域の電位を導体上の電荷から解析的に計算します。 これは、全ての導体を含む境界上のポテンシャル(電位)が導体 上の電荷密度と導体までの距離で決まるという条件を利用するのですが、 FreeFEMで解く場合は、解析すべき領域を含む境界のディリクレ条件として、
V = 1/(2*pi*eps)*(q1*log(r1) + q2*log(r2) + .. ) ここに、 q1, q2, .. = 線電荷密度 (q/m) eps = 誘電率 (F/m) r1, r2, .. = 境界上の点 (x,y) から線電荷までの距離 (m)を使って、逐次近似法で解くことになります。
例えば、往復2導体(対)の中心座標をそれぞれ (-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の場合は、プログラムを改造しないと無理ですから、分極電荷を無視 できるていどに外部境界を広げるのが良いようです。
実用的な静電シールドは相互キャパシタンスがあまりに小さくて、数値解析は無理 ですが、極めて目の荒いものは解析できて、例えば、次のようになります。
/* 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 分解」の結果を保存 しています。
このコースで解説した手法は、現代の有限要素法の技術としては極めて初歩的です が、コンピュータによる数値解析の威力に目覚め、たとえこの程度の技術でも知ら ない人と比べたらまるでレベルが違ってしまうことに気づいてほしいという希望か ら、かなり丁寧な例題を用意して、技術書には書かれていない、私のノーハウを伝 えようと努力しました。これを機会にさらに研讃を積むための意欲と学問と技術に 対する感動を生むことを願っています。日常に仕事の中にも、新鮮な工夫がほしい ものです。
この教材は 1998-08-26 に準備を始めて、1998-09-13 に仕上げました。
空間エネルギを使わずに、
導体の表面電荷からキャパシタンスを計算する場合は、
次のようになります。
導体表面の電界は法線方向しかありませんから、
表面の法線ベクトルを考えずに、
単に電界を積分をするだけで済むことに注意してください。
電気工学に置ける有限要素法の偉人Silvesterの名著
同軸ケーブルの内部導体として撚線を使ったときと、同じキャパシタンスになる
ような単線の直径がMeyersの「等価外径」で、外接円の直径より小さく
なる。
(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 = 導体表面の面積分(2次元なら周積分)
積分は導体の表面の周積分ですから、
FreeFEMの場合は int()[]を使うことになります。
(1) エネルギ誤差
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 の実効径
総素線数 | 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 |