表皮効果(プログラミング) - ケーブルの渦電流問題 (2)

数値計算は理工学で解析的に得られた数式を実務に活かすだけでなく、 数式の理解にも重要で、 計算するためのアルゴリズムを見出し、 プログラムを書き、 プログラムの実行結果が観測と一致したとき 始めてわかったという実感が得られます。

定常状態の渦電流問題では複素数のかなり複雑な計算が不可欠ですが、 Unixで生まれた K&R から ANSI の c89 までの c では

  typedef struct { double re; double im; } complex;

この解説では、c90 の実装が進んできた最近の c を使うことにして、 動作確認には FreeBSD-8.4 (gcc 4.2.1) と FreeBSD-10.3 (clang 3.4.1) の c を使うことにしました。 (注１)

ただ、前者の場合は csqrt() 以外の複素関数がありませんので、 必要な関数は自分で用意することにします。 (注２)

1. プログラミング例 - 円柱導体の電流分布

表皮効果を実感するのに最も適した計算は円柱導体の電流分布 (導体の渦電流問題 (1) - 表皮効果 の (12) 式)

  i(r)/i0 = I0(k*r)/I0(k*a)
  ここに、
	i(r) = 円柱導体の中心からの距離 r に於ける電流密度 (A/m^2)
	i0 = 円柱導体の表面電流密度 (A/m^2)
	r = 円柱導体の中心からの距離 (0 <= r <= a)
	a = 円柱導体の半径 (m)
	k = sqrt(ω*μ*ρ)
	ω = 2*π*f  .. 角周波数 (rad/s)
	f = 周波数 (Hz)
	μ = 透磁率 (H/m)
	ρ = 導電率 (S/m)
	ρ = 導電率 (S/m)
	I0(x) = 0 階の第１種変形 Bessel 関数

  周波数 絶対値 位相差

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <complex.h>
#include "Complex.h"

// current distribution of solid conductor

int main(argc, argv) char **argv;
{
  double a, f, u, g, r, PI;
  COMPLEX k, z;

  PI = atan2(0.0, -1.0);
  a = 5e-4;		// radius of solid copper conductor
  u = 4e-7 * PI;	// permeability
  g = 5.80e7;		// conductivity
  f = 2.1e5;		// frequency
  for (r = 0.0; r <= a * 1.001; r += a * 0.01) {
	k = csqrt(I * 2.0 * PI * f * u * g);
	z = I0(r * k) / I0(a * k);
	printf("%g %g %g\n", r, cabs(z), carg(z));
//	printf("%g %g %g %g\n", r, cabs(z), creal(z), cimag(z));
  }
  return 0;
}

上記のプログラムファイルの名前を x.c とすれば、 コンパイルの仕方は次のとおりです。

  cc -o x x.c Complex.c -lm

このプログラムの出力をグラフ化した例は次のとおりです。 X 軸が導体中心からの距離 (mm)、 黒線が電流分布の絶対値、 赤線が表面電流との位相差(ラジアン単位)ですが、 中心近くでは電流が表面とはほぼ逆方向に流れていることに注意してください。

電流の実部と虚部を見たいときは、コメント(//)行を使います。

１図 外径 1 mm 銅単線導体内部の電流分布 (周波数 210 kHz)

2. プログラミング例 - 表皮効果

以下、導体の渦電流問題 (1) - 表皮効果 で解説した各種の表皮効果 Zac/Rdc と同軸ケーブル外部導体の伝達インピーダンス を計算するためのプログラム例です。

  sold() - 単独単線導体の Zac/Rdc 周波数特性 - 1.2. (12)
  tube() - 単独円筒導体の Zac/Rdc 周波数特性 - 1.1. (11)
  bimetal() - 単独バイメタル導体の Zac/Rdc 周波数特性 - 1.4. (20)
  coax_outer() - 同軸ケーブル外部導体の Zac/Rdc 周波数特性 - 1.3 (20)
  Zt() - 同軸ケーブル外部導体の伝達インピーダンス(transfer impedance)周波数特性 - 1.3. (22)

プログラムのダウンロード

これらのプログラムの実行結果をグラフ化すると、下記のようになります。

• 単独単線導体の Zac/Rdc 周波数特性
• 単独円筒導体の Zac/Rdc 周波数特性
• 同軸ケーブル外部導体の Zac/Rdc 周波数特性
• 単独バイメタル導体の Zac/Rdc 周波数特性
• 同軸ケーブル外部導体の伝達インピーダンス(transfer impedance)周波数特性

同軸ケーブル外部導体の伝達インピーダンスの実部と虚部をみたいとか、 位相差も見たい場合は、 必要な printf() 文を活かしてください。

コンパイルには c90, c99, c11 の複素数が扱える c compiler が必要です。 FreeBSD-4.8 (gcc 4.2.1) と FreeBSD-10.3 (clang 3.4.1) で動作確認しましたが、 Mac OSX の Xcode や Linux でも使えると思います。

a.c 以外で link している -lamos (amos.a) は

  http://www.intex.tokyo/unix/cbes/cbes.html

このプログラムでは高精度の複素 Bessel 関数は必要ありませんので、 多公式近似と漸近展開の組み合わせで計算してもかまいませんが、 後に述べる近接効果や７芯導体の計算では精度の高い計算が必要になります。

なお、動作確認に使った FreeBSD-4.8 (gcc 4.2.1) では csqrt(), cabs(), carg() 以外の複素関数が用意されていませんが、 FreeBSD-10.3 (clang 3.4.1), MacOSX (clang-800.0.38) では、 cacos(), cacosh(), casin(), casinh(), catan(), catanh(), ccos(), ccosh(), cexp(), conj(), csin(), csinh(), ctan(), ctanh() など、一般的な関数は用意されています。

N. 注

N.1. 注１ - FreeBSD 以外の状況

N.2. 注２ - 複素数の Bessel 関数

N.3. 注３ - 計算結果の可視化

数値計算プログラムの開発過程では、 プログラムを書き換えてはグラフ化し、 それを見てプログラムを修正するという作業を何度も繰り返しますので、 Unixの場合は、 グラフ作成プログラムを数値計算プログラムにパイプライン(pipe line) で接続できるように設計すると、 仕事が非常に速くなります。

  $ 数値計算プログラム | グラフ化プログラム