有限要素法によるケーブル(電線)の電気設計 (4) - 渦電流問題

高周波ケーブルの減衰と単純な遅延以外の位相特性のほとんどを決めるのが、 導体に生ずる渦電流によるエネルギ損失と導体内部インダクタンスの変化ですが、 次の３つの機構があります。

• 表皮効果 - 導体中の電流分布が導体表面に偏る
• 近接効果 - 他の導体に流れる電流との相互インダクタンスによる電流の偏り
• シールド等、他の導体に生ずる渦電流の反作用

そして、この渦電流の計算は、なかなか難しいのですが、 通常TEM 波(Transverse Electro-magnetic Wave - 電磁界の進行方向成分がない) の伝送線路として使われる、 導電性のよい導体と損失の少ない絶縁体の組み合わせについては、 完全導体と損失のない絶縁体の組み合わせと同じように、 静電界と同じ断面の電磁界と 軸方向に進む進行波に分離することができて (注１)、 導体に流れる交流(正弦波)電流の分布を決める、 ベクトルポテンシャルの微分方程式はμの空間的変化がないとき次のようになります。

  -(d/dx+d/dy)(d/dx+d/dy)(A)) +j*ω*μ*σ*A + ω^2*μ*ε*A = 0     (1)

  ここに
        A = ベクトルポテンシャル (weber/m)
        j = sqrt(-1)
        ω = 角速度 (rad/s)
           = 2*π*f
        f = 周波数 (Hz)
        μ = 透磁率 (H/m)
        σ = 導電率 (S/m)
        σ = 導電率 (S/m)
           = 4e7*π (真空、空気、銅、アルミ等の場合)
        ε = 誘電率 (F/m)
           = 1e7/4/π/c^2*εs
        c = 真空中の光速 (2.998e8)
        εs = 比誘電率

この方程式を有限要素法で解こうとすると、 境界条件が必要になりますが、 ベクトルポテンシャルを決めるのは電流で、

  J = -σ*(d/dt)(A)
    = - j*σ*ω*A                                                         (2)
  A = (μ/(2*π))*∫J*log(r)*ds                                          (3)
  ここに、
	J = 電流密度 (A/m^2)
	r = 電流密度 J の場所から A を求める場所までの距離 (m)
	(d/dt)(A) = A の時間微分
	∫f*ds = f の面積分

しかし、まさに電流分布を知るために A を求めるわけですから、 鶏と卵の関係になって、うまくゆきません。 渦電流問題の難しさは、この境界条件が簡単に得られないことにあって、 偏微分方程式を扱う有限要素法は、 どちらかと言うと不向きです。 この点、境界条件の要らない積分方程式で扱うと、 この問題を避けられますから、 渦電流問題には積分方程式のほうが有利です。

偏微分方程式で必要な境界条件が容易に求まるのは、 次の場合です。

1. 単独円筒や同軸ケーブルのように軸対象の導体系
2. 完全導体や高周波で、電流が導体表面にしか存在しない場合

幸い、渦電流が大きな問題になるのは、高周波で、 表皮効果が大きくなった場合ですから、 この場合は、後者の表面電流だけを考えれば済む問題になって、 下記の手順で、渦電流による導体のインピーダンスを計算することができます。

まず、 スカラ・ポテンシャル(電位)を使った、 ラプラス方程式 (Laplace equation)

  -(d/dx+d/dy)(ε*(d/dx+d/dy)(V)) = 0                                    (2)

  ここに、
    V(x,y) = 電位 (V)
    ε(x,y) = 誘電率 (F/m)

次ぎに、以下の手順で導体表面の A の値を求めます。

  E = sqrt((dV/dx)^2 + (dV/dy)^2)                                        (3)
  q = ε*E                                                               (4)
  i = q*v                                                                (5)
  v = 1/sqrt(μ*ε)                                                      (6) 
  A = j*i/(ω*σ)                                                        (7)

  ここに、
	E = 導体表面の電界強度 (V/m)
	q = 導体表面の電荷密度 (C/m^2)
	ε = 導体表面の誘電体の誘電率 (F/m)
	i = 導体表面の電流密度 (A/m^2)
	v = 導体表面の進行波の位相速度 (m/s)
	μ = 導体表面の誘電体の透磁率 (H/m)
	A = 導体表面のベクトルポテンシャル (weber/m)
	σ = 導体表面の導電率 (S/m)
	σ = 導体表面の導電率 (S/m)
	j = sqrt(-1)

こうして、ベクトルポテンシャルが得られると、 導体抵抗と導体内部インダクタンスは下記の方法で計算できます。 (注２)

  Rac/Rdc = ∫abs(A)^2*ds/(∫A*ds)^2*S                                   (8)
  ω*Lac/Rdc = (1/(ω*μ*σ)*∫((dA/dx)^2 + (dA/dy)^2)*ds*S              (9)

  ここに、
	Rac = 周波数 ω/(2*π) に於ける交流抵抗 (Ohm)
	Rdc = 直流抵抗 (Ohm)
	    = 1/(S*σ)
	Lac = 周波数 ω/(2*π) に於ける導体内部インダクタンス (H)
	abs(z) = z の絶対値
	       = sqrt(x^2+y^2)  (z == x*i*y)
	∫*ds = 導体断面の面積分
	S = 導体断面積 (m^2)

ここで、Rac、Lac のいずれも、 ベクトルポテンシャルの偏りで決まることに注意してください。

以下、最初に軸対象問題を扱い、 次に、任意の形状の高周波に於ける渦電流問題を扱います。

1. 同軸ケーブル

同軸ケーブルの場合は、軸対象ですから、 導体の内外周上のベクトルポテンシャルが一定になって、 境界条件は極めて簡単で(注７)、 内部導体外周と外部導体内周のベクトルポテンシャルの比が、 計算できますので、それを使うと、次のようになります。

なお、同軸ケーブルの内部導体は孤立単線導体と同じ、 同軸ケーブルの外部導体は孤立した円筒導体と同じであることに注意してください。

5D-2V 同軸ケーブルの寸法で、上記の計算を実行すると、 次ぎのようになります。

/* skin effect of coaxial cable (5D-2V) */
complex;
f1:=1e5;	/* frequency */
a1:=0.7e-3;	/* radius of innner conductor */
a2:=2.4e-3;	/* inner rudius of outer conductor */
a3:=2.75e-3;	/* outer rudius of outer conductor */
u1:=4e-7*pi;	/* permeability */
g1:=5.8e7;	/* conductivity */
k1:=I*2*pi*f1*u1*g1;
e0:=1e7/(4*pi*2.998e8^2);
e1:=2.3*e0;
w1:=2*pi*f1;
k2:=w1*w1*u1*e0;
border(1,0,2*pi,100) { x:=a1*cos(t); y:=a1*sin(t) };
border(2,0,2*pi,300) { x:=a2*cos(t); y:=a2*sin(t) };
border(3,0,2*pi,300) { x:=a3*cos(t); y:=a3*sin(t) };
buildmesh(20000);

b1:=1;
b2:=b1*(1+2*log(a2/a1));
b3:=0;

k=k1*(one(region==0)+one(region==2))+k2*one(region==1);

/* vector potential */
solve(u) {
  onbdy(1) u = b1;	/* vector potential at surface */
  onbdy(2) u = b2;	/* vector potential at surface */
  onbdy(3) u = b3;	/* vector potential at surface */
  pde(u) -laplace(u) + id(u)*k = 0;
};

w=sqrt(Re(u)*Re(u)+Im(u)*Im(u))*1e12;	/* 1e12 = scaling for plot */
plot(w);

/* A method */
printf("Inner conductor:\n");
w=Re(u)*Re(u)+Im(u)*Im(u);
x1:=intt(0)[u];
x2:=Re(x1)*Re(x1)+Im(x1)*Im(x1);
s1:=pi*a1*a1;
x3:=intt(0)[w]/x2*s1;
wx=dx(u);
wy=dy(u);
w=Re(wx)*Re(wx)+Im(wx)*Im(wx)+Re(wy)*Re(wy)+Im(wy)*Im(wy);
x4:=intt(0)[w];
x5:=x4/x2*s1/(2*pi*f1*u1*g1);
printf("Zac/Rdc = %g + j*%g\n", x3, x5);
printf("Rac = %g, w*Lac = %g\n", x3/(s1*g1), x5/(s1*g1));

printf("Outer conductor:\n");
w=Re(u)*Re(u)+Im(u)*Im(u);
x1:=intt(2)[u];
x2:=Re(x1)*Re(x1)+Im(x1)*Im(x1);
s2:=pi*(a3*a3-a2*a2);
x3:=intt(2)[w]/x2*s2;
wx=dx(u);
wy=dy(u);
w=Re(wx)*Re(wx)+Im(wx)*Im(wx)+Re(wy)*Re(wy)+Im(wy)*Im(wy);
x4:=intt(2)[w];
x5:=x4/x2*s2/(2*pi*f1*u1*g1);
printf("Zac/Rdc = %g + j*%g\n", x3, x5);
printf("Rac = %g, w*Lac = %g\n", x3/(s2*g1), x5/(s2*g1));

かなりの計算量にはなりますが、 比較的低い周波数で、かなり細かく分割すれば、実用的な精度が得られます。 ただ、厚さの薄い外部導体については、肝心な半径方向の分割が荒くなりますし、 メッシュのほとんどが、本来不要な絶縁体部分に食われるため、 微分値を必要とする、インダクタンス計算では、 特に大きな誤差が出ます。

対称性を利用して、領域を絞っても、あまり改善されません。

/* skin effect of coaxial cable (5D-2V) */
complex;
f1:=1e5;	/* frequency */
a1:=0.7e-3;	/* radius of innner conductor */
a2:=2.4e-3;	/* inner rudius of outer conductor */
a3:=2.75e-3;	/* outer rudius of outer conductor */
u1:=4e-7*pi;	/* permeability */
g1:=5.8e7;	/* conductivity */
k1:=I*2*pi*f1*u1*g1;
e0:=1e7/(4*pi*2.998e8^2);
e1:=2.3*e0;
w1:=2*pi*f1;
k2:=w1*w1*u1*e0;
t1:=pi/4;
n3:=15;
n2:=floor(n3*(a2-a1)/(a3-a2)+0.5);
n1:=floor(n3*a1/(a3-a2)+0.5);
n4:=floor(n1*t1+0.5);
n5:=floor((n1+n2)*t1+0.5);
n6:=floor((n1+n2+n3)*t1+0.5);
border(1,0,n1+n4+n1,n1+n4+n1+1) {
 if (t<=n1) { x:=a1*t/n1; y:=0; ib:=4 }
 else if (t<n1+n4) { t0:=t1*(t-n1)/n4; x:=a1*cos(t0); y:=a1*sin(t0); ib:=1  }
 else { a0:=a1-(t-n1-n4)/n1*a1; x:=a0*cos(t1); y:=a0*sin(t1); ib:=5 };
};
border(2,0,n2+n5+n2,n2+n5+n2+1) {
 if (t<=n2) { x:=a1+(a2-a1)*t/n2; y:=0; ib:=4 }
 else if (t<n2+n5) { t0:=t1*(t-n2)/n5; x:=a2*cos(t0); y:=a2*sin(t0); ib:=2 }
 else { a0:=a2-(a2-a1)*(t-n2-n5)/n2; x:=a0*cos(t1); y:=a0*sin(t1); ib:=5 };
};
border(3,0,n3+n6+n3,n3+n6+n3+1) {
 if (t<=n3) { x:=a2+(a3-a2)*t/n3; y:=0; ib:=4 }
 else if (t<n3+n6) { t0:=t1*(t-n3)/n6; x:=a3*cos(t0); y:=a3*sin(t0); ib:=3 }
 else { a0:=a3-(a3-a2)*(t-n3-n6)/n3; x:=a0*cos(t1); y:=a0*sin(t1); ib:=5 };
};
buildmesh(20000);

b1:=1;
b2:=b1*(1+2*log(a2/a1));
b3:=0;

k=k1*(one(region==0)+one(region==2))+k2*one(region==1);

/* vector potential */
solve(u) {
  onbdy(1) u = b1;	/* vector potential at surface */
  onbdy(2) u = b2;	/* vector potential at surface */
  onbdy(3) u = b3;	/* vector potential at surface */
  pde(u) -laplace(u) + id(u)*k = 0;
};

w=sqrt(Re(u)*Re(u)+Im(u)*Im(u))*1e12;	/* 1e12 = scaling for plot */
plot(w);

/* A method */
printf("Inner conductor:\n");
w=Re(u)*Re(u)+Im(u)*Im(u);
x1:=intt(0)[u];
x2:=Re(x1)*Re(x1)+Im(x1)*Im(x1);
s1:=pi*a1*a1*t1/(2*pi);
x3:=intt(0)[w]/x2*s1;
wx=dx(u);
wy=dy(u);
w=Re(wx)*Re(wx)+Im(wx)*Im(wx)+Re(wy)*Re(wy)+Im(wy)*Im(wy);
x4:=intt(0)[w];
x5:=x4/x2*s1/(2*pi*f1*u1*g1);
printf("Zac/Rdc = %g + j*%g\n", x3, x5);
printf("Rac = %g, w*Lac = %g\n", x3/(s1*g1), x5/(s1*g1));

printf("Outer conductor:\n");
w=Re(u)*Re(u)+Im(u)*Im(u);
x1:=intt(2)[u];
x2:=Re(x1)*Re(x1)+Im(x1)*Im(x1);
s2:=pi*(a3*a3-a2*a2)*t1/(2*pi);
x3:=intt(2)[w]/x2*s2;
wx=dx(u);
wy=dy(u);
w=Re(wx)*Re(wx)+Im(wx)*Im(wx)+Re(wy)*Re(wy)+Im(wy)*Im(wy);
x4:=intt(2)[w];
x5:=x4/x2*s2/(2*pi*f1*u1*g1);
printf("Zac/Rdc = %g + j*%g\n", x3, x5);
printf("Rac = %g, w*Lac = %g\n", x3/(s2*g1), x5/(s2*g1));

さらに、周波数が高くなって、表皮効果が大きくなると、 メッシュのほとんどに電流が存在しなくなりますから、 まるで計算できなくなります。

1.1 同軸ケーブルの内部導体

そこで、内部導体と外部導体を別々に計算することにして、 本来計算しなくて良い絶縁体部分を対象領域から外します。

まず、内部導体については、 導体外周のベクトルポテンシャルが導体の電流によって決まる一定値になりますが、 その値そのもは、(8), (9) の計算では不要ですから、 何でもよくて、ここでは 1 にしておきます。

/* skin effect on solid circular conductor */
complex;
f1:=1e5;	/* frequency */
a1:=0.7e-3;	/* conductor radius */
u1:=4e-7*pi;	/* permeability */
g1:=5.8e7;	/* conductivity */
k1:=I*2*pi*f1*u1*g1;
border(1,0,2*pi,80) { x:=a1*cos(t); y:=a1*sin(t) };
buildmesh(8000);

/* vector potential */
solve(u) {
  onbdy(1) u = 1;	/* vector potential at surface */
  pde(u) -laplace(u) + id(u)*k1 = 0;
};

/* A method */
printf("Inner conductor:\n");
w=Re(u)*Re(u)+Im(u)*Im(u);
x1:=intt[u];
x2:=Re(x1)*Re(x1)+Im(x1)*Im(x1);
s1:=pi*a1*a1;
x3:=intt[w]/x2*s1;
wx=dx(u);
wy=dy(u);
w=Re(wx)*Re(wx)+Im(wx)*Im(wx)+Re(wy)*Re(wy)+Im(wy)*Im(wy);
x4:=intt[w];
x5:=x4/x2*s1/(2*pi*f1*u1*g1);
printf("Zac/Rdc = %g + j*%g\n", x3, x5);
printf("Rac = %g, w*Lac = %g\n", x3/(s1*g1), x5/(s1*g1));

w1:=2*pi*f1;
w=g1*w1*w1*(Re(u)*Re(u)+Im(u)*Im(u));	/* heat loss */
plot(w);	/* energy loss */

printf("energy method:\n");
q=Re(u)*Re(u)+Im(u)*Im(u);
i1:=intt[u];
printf("i1 = %g, %g\n", Re(i1), Im(i1));
p=Re(dx(u))*Re(dx(u))+Im(dx(u))*Im(dx(u))
	+ Re(dy(u))*Re(dy(u))+Im(dy(u))*Im(dy(u));
printf("p(%g), q(%g)\n", intt[p], intt[q]);
printf("Zac/Rdc = %g + j*%g\n",
	pi*a1*a1*intt[q]/(Re(i1)*Re(i1)+Im(i1)*Im(i1)),
	pi*a1*a1*intt[p]/(Re(i1)*Re(i1)+Im(i1)*Im(i1))/2/pi/f1/u1/g1);

printf("Rac = %g (Ohm/m)\n",
	intt[q]/(Re(i1)*Re(i1)+Im(i1)*Im(i1))/g1);
printf("Lac = %g (uH/m)\n",
	intt[p]/(Re(i1)*Re(i1)+Im(i1)*Im(i1))/(2*pi*f1)^2/u1/g1^2*1e6);

printf("\nvoltage/current method:\n");
e:=-I*w1*1;	/* V/m */
i:=-I*w1*g1*intt[u];	/* total current */
r:=1/(pi*a1*a1*g1);	/* DC registance */
z:=e/i/r;	/* Zac/Rdc */
printf("Zac/Rdc = %g + j*%g\n", Re(z), Im(z));

計算領域が狭いので、比較的低い周波数なら、 この程度の分割でも、0.5 % 程度の誤差になります。

さらに精度をよくするには、 対称性を利用して、領域を絞り、半径方向の分割を増やすと共に、 高い周波数に対応するため、 導体半径が表皮深さの３倍を越えるときは、 電流がほとんど存在しない導体内部を切捨てるようにすれば、 すべての周波数で誤差 1 % 以下になります。

/* skin effect of solid conductor */
complex;
f1:=1e6;	/* frequency */
a1:=0.7e-3;	/* conductor outer radius */
a2:=0;
m1:=4e-7*pi;	/* permeability */
g1:=5.8e7;	/* conductivity */
k1:=I*2*pi*f1*m1*g1;
s1:=sqrt(2/2/pi/f1/m1/g1);	/* skin depth */
t1:=pi/16;
printf("a/s = %g\n", a1/s1);
n2:=50;

if (a1/s1>3.0) {	/* high frequency */
 a2:=a1-3*s1;
 n1:=floor(n2*t1*a1/(a1-a2));
 n3:=floor(n2*t1*a2/(a1-a2));
 border(1,0,n1+n2+n3+n2,n1+n2+n3+n2+1) {
  if (t<=n1) { t0:=t1*t/n1; x:=a1*cos(t0); y:=a1*sin(t0); ib:=1 }
  else if (t<n1+n2) { t0:=a1-(a1-a2)*(t-n1)/n2;
   x:=t0*cos(t1); y:=t0*sin(t1); ib:=3 }
  else if (t<n1+n2+n3) { t0:=t1-t1*(t-n1-n2)/n3;
   x:=a2*cos(t0); y:=a2*sin(t0); ib:=2 }
  else { t0:=a2+(a1-a2)*(t-n1-n2-n3)/n2; x:=t0; y:=0; ib:=3 };
 };
}
else {
 n3:=floor(n2*t1);
 border(1,0,n3+n2+n2,n3+n2+n2+1) {
  if (t<=n3) { t0:=t1*t/n3; x:=a1*cos(t0); y:=a1*sin(t0); ib:=1 }
  else if (t<n3+n2) { a0:=a1-a1*(t-n3)/n2; x:=a0*cos(t1); y:=a0*sin(t1); ib:=2 }
  else { x:=a1*(t-n3-n2)/n2; y:=0; ib:=2 };
 };
};

buildmesh(10000);

/* vector potential */
solve(u) {
 onbdy(1) u = 1;       /* vector potential at surface */
 pde(u) -laplace(u) + id(u)*k1 = 0;
};

w1:=2*pi*f1;
w=g1*w1*w1*(Re(u)*Re(u)+Im(u)*Im(u));	/* heat loss */
plot(w);	/* energy loss */
/* energy method */
q=(Re(u)*Re(u)+Im(u)*Im(u))*2*pi/t1;
i1:=intt[u]*2*pi/t1;
p=(Re(dx(u))*Re(dx(u))+Im(dx(u))*Im(dx(u))
	+ Re(dy(u))*Re(dy(u))+Im(dy(u))*Im(dy(u)))*2*pi/t1;
printf("Zac/Rdc = %g + j*%g\n",
	pi*a1*a1*intt[q]/(Re(i1)*Re(i1)+Im(i1)*Im(i1)),
	pi*a1*a1*intt[p]/(Re(i1)*Re(i1)+Im(i1)*Im(i1))/2/pi/f1/m1/g1);

printf("Rac = %g (Ohm/m)\n",
	intt[q]/(Re(i1)*Re(i1)+Im(i1)*Im(i1))/g1);
printf("Lac = %g (uH/m)\n",
	intt[p]/(Re(i1)*Re(i1)+Im(i1)*Im(i1))/(2*pi*f1)^2/m1/g1^2*1e6);
/* V-I method */
e:=-I*w1*1;	/* V/m */
i:=-I*w1*g1*intt[u]*2*pi/t1;	/* total current */
r:=1/(pi*(a1*a1-a2*a2)*g1);	/* DC registance */
z:=e/i/r;	/* Zac/Rdc */
printf("Zac/Rdc = %.2f + j*%.2f\n", Re(z), Im(z));

1.2 同軸ケーブルの外部導体

次ぎに、外部導体の場合は、外部導体の外では内部の電流の合計が 0 になりますから、 ベクトルポテンシャルも 0 になりますので、 外周でベクトルポテンシャルが 0、 内周で一定値になります。

/* skin effect on solid circular conductor */
complex;
f1:=1e5;	/* frequency */
a1:=2.4e-3;	/* outer conductor inner radius */
a2:=2.75e-3;	/* outer conductor outer radius */
u1:=4e-7*pi;	/* permeability */
g1:=5.8e7;	/* conductivity */
k1:=I*2*pi*f1*u1*g1;
border(1,0,2*pi,800) { x:=a2*cos(t); y:=a2*sin(t) };
border(2,0,2*pi,900) { x:=a1*cos(-t); y:=a1*sin(-t) };
buildmesh(20000);

/* vector potential */
solve(u) {
  onbdy(1) u = 0;	/* vector potential at outer surface */
  onbdy(2) u = 1;	/* vector potential at inner surface */
  pde(u) -laplace(u) + id(u)*k1 = 0;
};

/* A method */
w=Re(u)*Re(u)+Im(u)*Im(u);
plot(w);
x1:=intt[u];
x2:=Re(x1)*Re(x1)+Im(x1)*Im(x1);
s1:=pi*(a2*a2-a1*a1);
x3:=intt[w]/x2*s1;
wx=dx(u);
wy=dy(u);
w=Re(wx)*Re(wx)+Im(wx)*Im(wx)+Re(wy)*Re(wy)+Im(wy)*Im(wy);
x4:=intt[w];
x5:=x4/x2*s1/(2*pi*f1*u1*g1);
printf("Zac/Rdc = %g + j*%g\n", x3, x5);
printf("Rac = %g, w*Lac = %g\n", x3/(s1*g1), x5/(s1*g1));

printf("\nvoltage/current method:\n");
w1:=2*pi*f1;
e:=-I*w1*1;	/* V/m */
i:=-I*w1*g1*intt[u];	/* total current */
r:=1/(s1*g1);	/* DC registance */
z:=e/i/r;	/* Zac/Rdc */
printf("Zac/Rdc = %g + j*%g\n", Re(z), Im(z));

このように全形状を使うと誤差が大きいのですが、 これは、厚さの薄い円筒構造では、 本来細かな分割を必要としない外周に比べて、 電流密度が変化する半径方向の分割が荒くなってしまうためです。

そこで、先程の円柱導体と同じように、 無駄な円周方向の分割を減らして、 計算資源を半径方向に回し、 さらに、表皮深さに比べて厚さが大きいときは、 電流が存在しない領域を削るようにすれば、 すべての周波数で誤差 1 % 未満になります。

/* skin effect on tubular conductor */
complex;
f1:=1e6;	/* frequency */
a1:=2.75e-3;	/* conductor outer radius */
a2:=2.4e-3;	/* conductor inner radius */
a20:=a2;

m1:=4e-7*pi;	/* permeability */
g1:=5.8e7;	/* conductivity */
k1:=I*2*pi*f1*m1*g1;
s1:=sqrt(2/2/pi/f1/m1/g1);	/* skin depth */
t1:=pi/16;
printf("a/s = %g\n", a1/s1);
n2:=50;

if ((a1-a2)/s1>3.0) {	/* high frequency */
 a2:=a1-5*s1;
};
n1:=floor(n2*t1*a1/(a1-a2));
n3:=floor(n2*t1*a2/(a1-a2));
border(1,0,n1+n2+n3+n2,n1+n2+n3+n2+1) {
 if (t<=n1) { t0:=t1*t/n1; x:=a1*cos(t0); y:=a1*sin(t0); ib:=1 }
 else if (t<n1+n2) { t0:=a1-(a1-a2)*(t-n1)/n2;
  x:=t0*cos(t1); y:=t0*sin(t1); ib:=3 }
 else if (t<n1+n2+n3) { t0:=t1-t1*(t-n1-n2)/n3;
  x:=a2*cos(t0); y:=a2*sin(t0); ib:=2 }
 else { t0:=a2+(a1-a2)*(t-n1-n2-n3)/n2; x:=t0; y:=0; ib:=3 };
};

buildmesh(8000);

/* vector potential */
solve(u) {
  onbdy(1) u = 0;	/* vector potential at surface */
  onbdy(2) u = 1;	/* vector potential at surface */
  pde(u) -laplace(u) + id(u)*k1 = 0;
};

w1:=2*pi*f1;
w=g1*w1*w1*(Re(u)*Re(u)+Im(u)*Im(u));	/* heat loss */
plot(w);	/* energy loss */

q=(Re(u)*Re(u)+Im(u)*Im(u))*2*pi/t1;
i1:=intt[u]*2*pi/t1;
p=(Re(dx(u))*Re(dx(u))+Im(dx(u))*Im(dx(u))
	+ Re(dy(u))*Re(dy(u))+Im(dy(u))*Im(dy(u)))*2*pi/t1;
printf("Zac/Rdc = %g + j*%g\n",
	pi*(a1*a1-a20*a20)*intt[q]/(Re(i1)*Re(i1)+Im(i1)*Im(i1)),
	pi*(a1*a1-a20*a20)*intt[p]/(Re(i1)*Re(i1)+Im(i1)*Im(i1))/2/pi/f1/m1/g1);

printf("Rac = %g (Ohm/m)\n",
	intt[q]/(Re(i1)*Re(i1)+Im(i1)*Im(i1))/g1);
printf("Lac = %g (uH/m)\n",
	intt[p]/(Re(i1)*Re(i1)+Im(i1)*Im(i1))/(2*pi*f1)^2/m1/g1^2*1e6);
/* V-I method */
e:=-I*w1*1;	/* V/m */
i:=-I*w1*g1*intt[u]*2*pi/t1;	/* total current */
r:=1/(pi*(a1*a1-a20*a20)*g1);	/* DC registance */
z:=e/i/r;	/* Zac/Rdc */
printf("Zac/Rdc = %.2f + j*%.2f\n", Re(z), Im(z));

同軸ケーブルの内部導体が円筒の場合は、 導体内周のベクトルポテンシャルが 0、 導体外周のベクトルポテンシャルが一定値になりますので、 境界条件を逆にするだけで、 同じプログラムが使えます。

2. その他のケーブルの概算

同軸ケーブル以外は、 FreeFEM内で一部の導体だけを選択して解析する方法がありませんので、 実用的な計算方法がありませんが、 一般に、近接効果や他の導体の反作用は表皮効果に比べて、 そう大きなものではありませんし、 周波数が低いときは、近接効果も小さいので、 個々の導体の表皮効果だけを計算して、 それを経験値を元にして、経験値で適当に割増しする程度でも、 数 % から 15 % 程度の誤差で真の値を予測することができます。

高周波については、 事項で述べるようなTEM波の特性を活かした解法があって、 渦電流の完全な解析ができます。

3. 高周波ケーブル

周波数が高くなって、 導体サイズが表皮深さ(skin depth)に比べて、充分大きくなると、 導体の交流抵抗と内部インダクタンスは、 導体の表面電流密度と表皮抵抗(surface registance)だけで計算することができます。 (注４)

そして、伝送線路のTEM波の場合、 導体断面の電界は同じ構造の静電界と同じLaplace方程式を見たしますから、 静電界問題として、導体表面の電界強度が計算できて、 そこから導体表面の電荷密度が得られます。 一方、導体表面の(見掛けの)電荷の移動も速度は、 導体表面の電磁波の位相速度を同じですから、 この２つをもとに、導体の表面電流が計算できることになって、 高周波ケーブルの渦電流問題は、 取り扱いの容易な、静電界問題に帰着できることになります。 (注５)

この手法による計算手順は、次のようになります。(注６)

1. 導体系の静電界に於ける電位 V(x,y) を求める。
2. 表面電荷密度を次式で求める。

     ρ = ε*dV/dn					(3)
     ここに、
   	ρ = 導体の表面電荷密度 (C/m^2)
   	ε = 導体表面の誘電率 (F/m)
   	dV/dn = 導体表面に於ける V の法線方向微分 (V/m)
   

3. 表面電流密度(surface currecnt density)を次式で計算する。

     i = ρ / sqrt(ε*μ)					(4)
     ここに、
   	i = 導体の表面電流密度 (A/m)
   	ε = 導体表面の誘電率 (F/m)
   	μ = 導体表面の透磁率 (H/m)
   

4. 個々の導体について、電力損失を次式で計算する。

     P = ∫Rs*i^2*dl (W)				(5)
     ここに、
   	Rs = sqrt(ω*μ/(2*ρ))
   	   = 表皮抵抗 (Ω/m^2)
   	∫dl = 境界(表面)上の周積分
   

5. 個々の導体について、全電流を次式で計算する。

     I = ∫i*dl						(6)
     ここに、
   	I = 全電流 (A)
   

6. 個々の導体について、導体抵抗を次式で計算する。

     R = P / I^2						(7)
     ここに、
   	R = 導体抵抗 (Ω/m)
   

7. 個々の導体について、導体内部インダクタンスを次式で計算する。

     Li = R / ω						(8)
     ここに、
   	L1 = 導体内部インダクタンス (H/m)
   	ω = 2*π*f
   	   = 角周波数 (rad/s)
   	f = 周波数 (Hz)
   

  ∫i^2*dl >= (∫i*dl)^2

導体の高周波抵抗が得られれば、 特性インピーダンスの値から、 下記のように減衰定数が計算できます。

  α = Rac / 2 / Z0
  ここに、
	α = 減衰定数 (neper/m)
	Rac = 導体の高周波抵抗 (Ω/m - 往復)
	Z0 = 特性インピーダンス (Ω - 正確には特性インピーダンスの実部)

以下、解析解がわかっている、同軸ケーブルとツインリードの構造について、 上記のアルゴリズムを検証してみましょう。

3.1 同軸ケーブル

最初に、表皮効果だけを考えれば済む同軸構造を扱うことにして、 内部導体径 0.8 mm、外部導体内径 4.9 mm、 ポリエチレン絶縁の同軸ケーブル(5C-2V の寸法)で、 100 MHz の交流抵抗と減衰定数を求めました。

ここでは、断面全体を解析していますが、 対称性を利用して、断面の一部を解析すると、 対称軸に直交する境界(自然境界)の分割をうまくやらないと、 要素数が劇的に増えて要素分割を打ち切られ、 かえって精度が低下することがあります。

/* coaxial cable: 5C-2V */
a1:=0.4e-3;     /* inner conductor raius */
a2:=2.45e-3;   /* outer conductor radius */
eps:=2.3;       /* dielectric constant */
e0:=8.854e-12;	/* permittivity */
e1:=eps*e0;
m1:=4e-7*pi;	/* permeability */
g1:=5.8e7;	/* conductivity */
f1:=100e6;	/* frequency */

n1:=50;
border(1,0,2*pi,6*n1) { x:=a2*cos(t); y:=a2*sin(t) };
border(2,0,2*pi,n1) { x:=a1*cos(-t); y:=a1*sin(-t) };
buildmesh(8000);

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

/* Characteristic impedance, Velocity ratio */
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\n", z);

/* Rac at high frequency */

e=sqrt(dx(v)^2+dy(v)^2);
plot(e);
w1:=2*pi*f1;
rs:=sqrt(2*pi*f1*m1/2/g1);	/* surface resistance */
j=-sqrt(e1/m1)*sqrt(dx(v)^2+dy(v)^2);	/* surface current density */
j2=(e1/m1)*(dx(v)^2+dy(v)^2);
i1:=int(1)[j];		/* total current in outer conductor */
i2:=int(2)[j];		/* total current in inner conductor */
p1:=int(1)[j2]*rs;	/* power loss in outer conductor */
p2:=int(2)[j2]*rs;	/* power loss in inner conductor */
r1:=p1/i1/i1;		/* outer conductor ac-resistance */
r2:=p2/i2/i2;		/* inner conductor ac-resistance */
l1:=r1/w1;		/* inner conductor internal inductance */
l2:=r2/w1;		/* outer conductor internal inductance */
printf("Rac = %g + %g \n", r2, r1);
printf("Attenuation = %g neper/m (%g Hz)\n", (r1+r2)/z/2, f1);	/* attenuation */

中心導体の交流抵抗を解析的計算から得られた値と比較した結果は、 次のようなります。 導体厚 t は内部導体半径の 1/2 として t/δ を表示してあります。

周波数(Hz)	理論値	計算値	t/δ	誤差(%)
1k	0.0343016	0.00328525	0.191	940
10k	0.0343963	0.0103889	0.605	230
100k	0.0421711	0.0328525	1.914	28
1M	0.112903	0.103889	6.052	8.6
10M	0.33701	0.328525	19.1	2.6
100M	1.0467	1.03889	60.5	0.75
1G	3.29126	3.28525	191	0.18

この計算方法の前提は t/δ >> 1 ですが、 10 MHz 程度でもかなりの精度が得られ、 100 MHz 以上なら充分な精度が得られることがわかります。 また、解析的計算による高周波円筒導体の近似式、

  Rac = Rs / (2 * π * a)
  ここに、
	Rac = 円形導体の交流抵抗 (Ω/m)
	π = 3.14159265..
	a = 導体半径(m - 同軸ケーブルの外部導体の場合は内半径)

特性インピーダンス Z0、速度係数 vr は静電ポテンシャルの計算から得られますから、 高周波特性のすべてが、有限要素法で簡単に解けることになります。

3.2 ツインリード

表皮効果以外に近接効果が無視できない、 導体径 0.5 mm、導体間隔 1.0 mm のツインリードで、 100 MHz の場合を計算します。

/* open boundary - bare wire pair */

a1:=0.25e-3;		/* conductor radius */
d1:=1.0e-3;	/* (distance between conductor)/2 */
x0:=10*a1; y0:=5*a1;	/* boundary size */
e0:=1e7/(4*pi*2.998e8^2);	/* dielectric constant for vaccum */
m0:=4e-7*pi;	/* permeability */
g1:=5.8e7;	/* frequency */
f1:=100e6;	/* frequency */

/* bundary */
border(1,-x0,x0,40) { x:=t; y:=-y0 };
border(1,-y0,y0,40) { 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);

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

 /* charge on conductor */
 e=sqrt(dx(v)*dx(v)+dy(v)*dy(v));
 q0:=q;
 q:=e0*int(2)[e];	/* Q/m */
/* printf("C = %.1f (pF/m)\n", q*1e12);*/
/* printf("err: %g\n", abs(q0-q)/q);*/
};
printf("C = %.1f (pF/m)\n", q*1e12);

w1:=2*pi*f1;
rs:=sqrt(2*pi*f1*m0/2/g1);	/* surface resistance */
j=-sqrt(e0/m0)*sqrt(dx(v)^2+dy(v)^2);	/* surface current density */
j2=(e0/m0)*(dx(v)^2+dy(v)^2);
i1:=int(2)[j];		/* total current in outer conductor */
i2:=int(3)[j];		/* total current in inner conductor */
p1:=int(2)[j2]*rs;	/* power loss in outer conductor */
p2:=int(3)[j2]*rs;	/* power loss in inner conductor */
r1:=p1/i1/i1;		/* outer conductor ac-resistance */
r2:=p2/i2/i2;		/* inner conductor ac-resistance */
l1:=r1/w1;		/* inner conductor internal inductance */
l2:=r2/w1;		/* outer conductor internal inductance */
printf("Rac = %g + %g\n", r2, r1);

  外径 0.5 mm の表皮効果を考慮した交流抵抗 = 1.68308 Ω/m
  近接効果(係数) = 1.841

つまり、前記のアルゴリズムから、 表皮効果と近接効果の両方を考慮した交流抵抗が求まることがわかります。

4. 注

4.1 注１- TEM 波の特徴

  John Carson,- The Rigorous and Approximate Theories of Electrical
	Transmission Along Wires, B.S.T.J., Oct. 1926

  太田浩一,- 電磁気学 II
	(丸善) ISBN 4-621-04805-8

  E = -j*ω*A
  I = σ*E
  abs(Iac)^2*Rac = ∫E*I*ds
  Iac = ∫I*ds
  Rdc = 1/(σ*S)

  ここに、
	E = 軸方向電界強度 (V/m)
	i = sqrt(-1)
	ω = 2*π*f
	f = 周波数 (Hz)
	A = 軸方向ベクトルポテンシャル (weber/m)
	I = 軸方向電流密度 (A/m^2)
	σ = 導電率 (S/m)
	σ = 導電率 (S/m)
	∫*ds = 導体断面の面積分
	Iac = 導体全体の交流電流 (A)
	Rdc = 導体の直流抵抗 (Ohm/m)
	S = 導体断面積 (m^2)

  Rac/Rdc = ∫abs(A)*ds/abs(∫A*ds)^2*S

  Wm = (1/2)*Lac*abs(Iac)^2
     = ∫(B*H/2)*ds
     = (1/2*μ)*∫[(dA/dx)^2 + (dA/dy)^2]*ds
  ここに、
	Wm = 導体内部の磁気エネルギ (joul/m^3)
	Lac = 交流の導体内部インダクタンス (H/m)
	B = 磁束密度 (weber/m)
	H = 磁界強度 (ampere-turn/m)

  ω*Lac/Rdc = (1/(ω*μ*σ)*∫((dA/dx)^2 + (dA/dy)^2)*ds*S

 P.Silvester,- Network Analog Solution of Skin and Proximity Effect Problems
	IEEE Transuctions on Power Apparatus and Systems, Vol.PAS-86, No.2
	pp 241-247, Febtuary, 1967

  A = (μ*i/(4*π))*r^2/a^2              導体内部
    = (μ*i/(2*π))*(1/2 + log(r/a))     導体外部

  ここに、
	A = ベクトルポテンシャル (weber/m) (軸方向)
	μ = 透磁率 (H/m)
	i = 導体の全電流 (A)
	r = 導体中心軸からの距離 (m)
	a = 導体半径 (m)

  Aa = μ*i/(4*π)
  Ab = μ*i/(4*π)*(1/2 + log(b/a))
     = Aa*(1 + 2*log(b/a))
  ここに、
	Aa = 内部導体外周のベクトルポテンシャル (weber/m)
	Ab = 外部導体内周のベクトルポテンシャル (weber/m)
	b = 外部導体内径 (m)

外部導体外周では、それよち内側の全電流が 0 ですから、 ベクトルポテンシャルは 0 になります。

4.4 注４ - 表皮の深さと表面インピーダンス

  δ = sqrt(2/(ω*μ*σ))
  ここに、
	δ = 表皮深さ(skin depth) (m)
	ω = 2*π*f (rad/s)
	f = 周波数 (Hz)
	μ = 透磁率 (H/m)
	σ = 導電率 (S/m)
	σ = 導電率 (S/m)

  V = Zs * I
  Zs = Rs + j*Xs
     = Rs * (i + j)
  Rs = 1 / (σ * δ)
  ここに、
	V = 導体表面の電位差 (V/m)
	I = 導体内部の全電流 (A)
	Zs = 表面インピーダンス (Ω/m)
	Rs = 表面抵抗 (Ω/m)
	Xs = 表面リアクタンス (Ω/m)
	σ = 導体の導電率 (S/m)
	σ = 導体の導電率 (S/m)
	δ = 表皮深さ (m)

  Li = Xs/ω
     = Rs/ω (H/m)

4.5 注５ - 導体を流れる電流

  I = j*ω*ε*∫En*ds
  ここに、
	I = 導体を流れる電流 (A)
	ω = 2*π*f (rad/s)
	f = 周波数 (Hz)
	ε = 誘電率 (F/m)
	En = 導体の表面(法線方向)の電界の強さ
	∫ds = 導体の表面上の１周積分

  V.Belevitch,- Theory of the Proximity Effect in Multiwire Cables
        (Philips Res. Repts, 32, 16-43, 1977)

4.6 注６ - 計算例

  E = V/(r*log(b/a))
  ここに、
        E = 中心軸から r 離れた場所の電界の強さ (V/m)
        a = 内部導体半径 (m)
        b = 外部導体半径 (m)
        V = 内部導体の電位 (V .. 外部導体は 0 V)

  ρ = ε*E
     = ε*V/(a*log(b/a))

  i = ρ / sqrt(ε*μ)
    = sqrt(ε/μ)*V/(a*log(b/a))

  P = ∫Rs*i^2*dl
    = Rs*2*π*a*(ε/μ)*V^2/(a*log(b/a))^2

  I = ∫j*ds
    = 2*π*a*sqrt(ε/μ)*V/(a*log(b/a))

  R = P / I^2
    = Rs/(2*π*a)

  V/I = (1/2*π)*sqrt(μ/ε)*log(b/a)

4.7 注７ - 同軸ケーブルのベクトルポテンシャル

  Aa = μ*i/(4*π)
  Ab = μ*i/(4*π)*(1/2 + log(b/a))
     = Aa*(1 + 2*log(b/a))
  ここに、
	Aa = 内部導体外周のベクトルポテンシャル (weber/m)
	Ab = 外部導体内周のベクトルポテンシャル (weber/m)
	b = 外部導体内径 (m)