渦電流問題は磁界が広範囲に広がるため、有限要素法では扱いにくいのですが、 円柱導体の場合は境界条件が簡単で、 ベクトルポテンシャルを使って解いてみましょう。 境界条件は、導体外周のベクトルポテンシャルが任意の一定値になることです。

まず、導体内部のベクトルポテンシャルの微分方程式は次のようになります。

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

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

  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

  または、

  Rac/Rdc = Real(E/I)/Rdc = Real(-j*ω*A/(-j*ω*σ∫A*ds))/Rdc

// skin effect on solid circular conductor
// 2006-02-24 Kouichi Hirabayashi
real f1=1e5;		// frequency
real a1=0.5e-3;		// conductor radius
real m1=4e-7*pi;	// permeability
real g1=5.8e7;		// conductivity
complex k1=2*pi*f1*m1*g1*1.0i;
border b0(t=0,2*pi) { x=a1*cos(t); y=a1*sin(t); label=1; }
mesh Th = buildmesh(b0(400));
plot(Th,wait=1);
fespace Vh(Th,P1);

Vh<complex> u,v;
solve ec(u,v) = int2d(Th)(u*v*k1
	+(dx(u)*dx(v)+dy(u)*dy(v)))
	+ on(1,u=1);
Vh ur=real(u), ui=imag(u);
plot(ur,wait=1, fill=true);
plot(ui,wait=1, fill=true);

// heat loss
real w1=2*pi*f1;
Vh w=g1*w1*w1*(real(u)*real(u) + imag(u)*imag(u));
plot(w,wait=1, value=true, fill=true);

// energy method
cout << "energy method:" << endl;
Vh q=real(u)*real(u) + imag(u)*imag(u);
complex i1=int2d(Th)(u);
Vh p=real(dx(u))*real(dx(u)) + imag(dx(u))*imag(dx(u))
	+ real(dy(u))*real(dy(u)) + imag(dy(u))* imag(dy(u));
real x1=pi*a1*a1*int2d(Th)(q)/(real(i1)*real(i1)+imag(i1)*imag(i1));
real y1=pi*a1*a1*int2d(Th)(p)/(real(i1)*real(i1)+imag(i1)*imag(i1))
	/2/pi/f1/m1/g1;
cout << "Zac/Rdc = " << x1 << " + j*" << y1 << endl; 

// voltage/current method
cout << "voltage/current method:" << endl;
complex e=-w1*1.0i;			// V/m
complex i=-w1*g1*1.0i*int2d(Th)(u);	// total current
real r=1.0/(pi*a1*a1*g1);		// DC resistance
complex z=e/i/r;			// Zac/Rdc
cout << "Zac/Rdc = " << real(z) << " + j*" << imag(z) << endl;

周波数が高くなると、表皮効果が大きくなって、 電流が導体周辺に集中しますから、 導体全面を均等に分割すると、大きな誤差がでます。

FreeFem++ の場合は、adaptmesh()という強力なツールがありますので、 もう少し高い周波数で、この関数を試してみましょう。 adaptmesh() を実行する都度変化する分割の仕方に注目してください。

// skin effect on solid circular conductor
// 2006-04-10 Kouichi Hirabayashi
real f1=1e7;		// frequency
real a1=0.5e-3;		// conductor radius
real m1=4e-7*pi;	// permeability
real g1=5.8e7;		// conductivity
complex k1=2*pi*f1*m1*g1*1.0i;
border b0(t=0,2*pi) { x=a1*cos(t); y=a1*sin(t); label=1; }
mesh Th = buildmesh(b0(200));
plot(Th,wait=1);
fespace Vh(Th,P1);
real w1=2*pi*f1;

Vh<complex> u,v;
Vh w;
problem ec(u,v) = int2d(Th)(u*v*k1
	+(dx(u)*dx(v)+dy(u)*dy(v)))
	+ on(1,u=1);
real error=0.1;
for (int i=0;i<3;i++) {
  ec;
  w=g1*w1*w1*(real(u)*real(u) + imag(u)*imag(u));
  Th = adaptmesh(Th,w,err=error);
  plot(Th,wait=1);
  error/=5;
};
Vh ur=real(u), ui=imag(u);
plot(ur,wait=1, fill=true);
plot(ui,wait=1, fill=true);

// heat loss
w=g1*w1*w1*(real(u)*real(u) + imag(u)*imag(u));
plot(w,wait=1, value=true, fill=true);

// energy method
cout << "energy method:" << endl;
Vh q=real(u)*real(u) + imag(u)*imag(u);
complex i1=int2d(Th)(u);
Vh p=real(dx(u))*real(dx(u)) + imag(dx(u))*imag(dx(u))
	+ real(dy(u))*real(dy(u)) + imag(dy(u))* imag(dy(u));
real x1=pi*a1*a1*int2d(Th)(q)/(real(i1)*real(i1)+imag(i1)*imag(i1));
real y1=pi*a1*a1*int2d(Th)(p)/(real(i1)*real(i1)+imag(i1)*imag(i1))
	/2/pi/f1/m1/g1;
cout << "Zac/Rdc = " << x1 << " + j*" << y1 << endl; 

// voltage/current method
cout << "voltage/current method:" << endl;
complex e=-w1*1.0i;			// V/m
complex i=-w1*g1*1.0i*int2d(Th)(u);	// total current
real r=1.0/(pi*a1*a1*g1);		// DC resistance
complex z=e/i/r;			// Zac/Rdc
cout << "Zac/Rdc = " << real(z) << " + j*" << imag(z) << endl;

さらに、最適化を行えば、もっと高い周波数でも計算できますが、 高い周波数については、 表面から表皮深さの３倍とか５倍までの深さだけ残して、 電流がほとんど存在しない内側を削り、 円筒形の導体にしてしまうほうが良いでしょう。

円筒導体の場合は、円筒内表面の内側には磁束がありませんから、ベクトルポテンシャ ルは 0 で、円筒内面の境界条件を下記のように追加します。

real a2=0.4e-3;		// inner surface
..
border b0(t=0,2*pi) { x=a1*cos(t); y=a1*sin(t); label=1; }
border b1(t=0,2*pi) { x=a2*cos(-t); y=a2*sin(-t); label=2; }
mesh Th = buildmesh(b0(200)+b1(200));
..
problem ec(u,v) = int2d(Th)(u*v*k1
	+(dx(u)*dx(v)+dy(u)*dy(v)))
	+ on(1,u=1)
	+ on(2,u=0);
..

円筒の内側磁界の大きさは、円筒内部を埋めた円柱の磁界と、円筒内面と同じ半径を持 つ円柱に同じ大きさの電流を逆方向に流したときの磁界を加算する方法で求めるのが簡 単ですが、円筒外表面の中心と円筒内表面中心との距離に比例した一様(!)磁界になる という面白い結果になります。

なお、同軸ケーブルの外部導体の場合は、ベクトルポテンシャルが外部導体表面で 0 になりますから、境界条件が逆になって、

..
problem ec(u,v) = int2d(Th)(u*v*k1
	+(dx(u)*dx(v)+dy(u)*dy(v)))
	+ on(1,u=1)
	+ on(2,u=0);
..