コイルの軸から離れた点の磁場を計算したい その2(ビオ・サバール汎用コード製作)
コイル軸上に無い(オフアクシスの)任意の点の磁場は?ヘルムホルツコイルを作って磁場を発生させたとき、径方向にどれぐらいの成分が出るのか、つまりどれぐらいの領域なら平行磁場と言っていいのか、を簡単に知りたい。解析解はない、ということで、ビオサバールの式を数値的に積分して磁場計算するコードを作りました。昨日書いた通り、外積形式のビオ・サバールの式に座標値を入れれば計算できます。微小電流要素Cの座標を(Cx, Cy, Cz)磁場を計算したい点Pの座標を(Px, Py, Pz)とするとCからPに向くベクトルrの成分は (Px-Cx, Py-Cy, Pz-Cz)微小電流要素を流れる電流ベクトルdIの成分を(Ix, Iy, Iz)とすると、誘起される磁場のベクトルdH=(Hx, Hy, Hz)の各成分は、ビオサバールの式に各成分を入れてやって以下のようになります。ただしまずは汎用ビオサバールの計算コードを作りました。Visual Bacic(.NET)でつくりました。磁場を求めたい点Pの座標(Px,Py,Pz)、微小電流要素Cの座標(Cx, Cy,Cz)、電流ベクトルIの成分(Ix, Iy, Iz)を引数として関数BSelementCalcを呼ぶと、点Pでの磁場成分(Bx, By,Hz)を一次元配列として返します。引数の単位は[m]と[A]で、戻り値は磁場[テスラ]です。Private Function BSelementCalc(Px As Double, Py As Double, Pz As Double, Cx As Double, Cy As Double, Cz As Double, Ix As Double, Iy As Double, Iz As Double) '戻り値の配列 Dim ABB(2) As Double 'ビオサバールの定数部分 Dim AAA As Double = 1 / ((Px - Cx) ^ 2 + (Py - Cy) ^ 2 + (Pz - Cz) ^ 2) ^ 1.5 * 0.0000001 '本当は4πで割るが、どうせ磁場を求めるときに透磁率μ(=4πx1e7)を掛けるので、もう最初からかけたことにして4πは省略 'ビオサバールの外積部分 ABB(0) = AAA * (Iy * (Pz - Cz) - Iz * (Py - Cy)) ABB(1) = AAA * (Iz * (Px - Cx) - Ix * (Pz - Cz)) ABB(2) = AAA * (Ix * (Py - Cy) - Iy * (Px - Cx)) '配列で戻す Return ABBEnd Functionこれを使って、Z軸上に置かれたコイルが任意点に作る磁場を計算したのが下のコード。微小電流要素の座標Cと、電流ベクトルIをコイルの円環に合わせて変化させて、足し合わせます。例として、原点から150[mm]に置かれた半径40[mm]のコイルが、電流500[A]で座標(10,20,200)[mm]に作る磁場(Bx, By, Bz)[gauss]を計算してみます。コイルの円環は180分割(2度刻み)にしてます。 '汎用3次元ビオサバールでコイルの作る磁場を計算 'わかりやすく作ってあるので高速化は無視Private sub BStestCalculation() Dim CoilR As Double = 40/1000 ’コイルの半径[m] Dim CoilZ As Double = 150/1000 'コイルのz座標[m] Dim CoilI As Double =500 ’コイルに流れる電流[A] '磁場を求めたい任意点の座標[m] Dim Px As Double = 10/10000 Dim Py As Double = 20/10000 Dim Pz As Double = 200/10000 '微小電流要素の座標と、電流ベクトル Dim Cx, Cy, Cz As Double Dim Ix, Iy, Iz As Double 'コイルの円周上の回転角 Dim theta As Double Dim ThetaDivide As Integer = 180 'コイル円周の分割数 Dim ttt As Integer '計算された磁場 Dim Bx As Double = 0 Dim By As Double = 0 Dim Bz As Double = 0 '関数戻りを収める配列 磁場[テスラ] Dim resultABB() As Double 'まず、コイルはz軸上に置かれているのでz成分はゼロに固定 Cz = 0 Iz = 0 'コイルの円周を分割した微小要素からの磁場の寄与の和を求める For ttt = 0 To ThetaDivide - 1 ’コイル上の座標 theta = (ttt / ThetaDivide) * 2 * PI Cx = CoilR * Cos(theta) Cy = CoilR * Sin(theta) ’電流ベクトル Ix = -1 * Sin(theta) * CoilI Iy = Cos(theta) * CoilI '求めたい位置のxyz座標、要素のxyz座標 電流ベクトルの成分を引数にBS計算 resultABB = BSelementCalc(Px, Py, Pz, Cx, Cy, Cz, Ix, Iy, Iz) '微小要素の寄与の和を計算 Bx = Bx + resultABB(0) By = By + resultABB(1) Bz = Bz + resultABB(2) Next ttt 'コイルを分割した数だけ倍化されているので割る テスラ → ガウスの変換で1万倍する '足し合わせた微小要素の長さは2π×coilRの円周になるので、コイル全体の寄与としてはそれだけ大きくなるので円周分かけないといけない 微小要素の寄与トータルが必要 Bx =Bx / ThetaDivide * (2 * PI * CoilR) * 10000 By =By / ThetaDivide * (2 * PI * CoilR) * 10000 Bz =Bz / ThetaDivide * (2 * PI * CoilR) * 10000End subこれを実行すると、コイル中心軸とオフアクシスな点Pの磁場はBx=3.14By=6.28Bz=15.65ガウスであると求まります。どの教科書にも載っている式で検証してみます。半径a[m]の円環電流I[A]が中心に作る磁場B[T]は、これにa=1[m]、I=1[A]を計算するとB=6.28e-7[T]になります。上記プログラムも同じ値になります。コイルの座標と電流ベクトルさえ正しく振ってやれば、コイルがどのような位置・どのような配置・どのような傾きに置かれていても、任意点Pの磁場が計算できます。直線であろうが四角であろうが楕円であろうが何でも求まります。これでも十分使えますが、ヘルムホルツコイルを作って磁場を発生させたとき、径方向にどれぐらいの成分が出るのか、つまりどれぐらいの領域なら平行磁場と言っていいのか、を簡単に知りたいんでしょ?というのがもともとの動機なので、コイルのジオメトリを限定して、もうちょっとラクに使えるコードにしました。とりあえず休憩しましょうOVALE オヴァール とろとろプリンタルト 9個入 お取り寄せ/洋菓子/スイーツ/タルト/プリン/ギフト/贈答/詰め合わせ/セット/個包装明日に続く。