投稿2019/01/11 14:35
編集2019/01/12 08:29永久磁石がある点に作る磁束密度をC言語にて計算したいと考えています。
計算がうまくいかず、どこに間違いがあるのか見つけられていないので、
電磁気ないしC言語に見識のある方のお力添えを頂きたく存じます。
表面磁場 M = 0.24 T
磁石長さ L = 1 mm (y軸方向)
磁石幅 W = 1 mm (z軸方向)
の磁石がある点に作る磁場を求めます。
磁石表面がy-z表面上にあるとして、長さ方向(y軸方向)にn分割してdy、幅方向(z軸方向)にm分割してdzとし、
それぞれの領域が作る磁場をyとzに関する重積分で重畳して求めようとしています。
この時、磁石の表面磁束密度を M で一様だと仮定すると、分割された各要素の
持つ磁極の強さは dydzM [Wb] となります。
したがって、磁場のクーロンの法則より、各要素から r [m]離れた点に磁石の要素が作る磁場dHは、
dH = (1/4πµo) * {(dydzM) / r^2}
磁束密度dB [T]に直して、
dB = (1/4π) * {(dydzM) / r^2)}
ここで、
要素の座標が ( xs0 , ys0 , zs0 ) (添え字はソース点の意味でsです。読みにくくてすみません)
磁束密度を求める点が ( x , y , z )
だとすると、2点間の距離 r = √{(x-xs0)^2 + (y-ys0)^2 + (z-zs0)^2} なので
dB = (1/4π) * [(dydzM) /{(x-xs0)^2 + (y-ys0)^2 + (z-zs0)^2}]
さらに、x軸方向の磁束密度を求めたいため
dB = (1/4π) * [(dydzM) /{(x-xs0)^2 + (y-ys0)^2 + (z-zs0)^2}] * [(x-xs0) / √{(x-xs0)^2 + (y-ys0)^2 + (z-zs0)^2}]
としてそれぞれ求め、これを
y軸方向に 0 → 10^(-3) [m]
z軸方向に 0 → 10^(-3) [m]
の範囲で積分して計算しています。
ここまでで間違いはございますでしょうか。
簡単な静磁場の計算ですが、得手でなくお恥ずかしい限りです。
積分には台形公式を用いています。以下、実際のソースコードになります。
磁束密度を求める点は ( 10^(-3) , 0 , 0 )としています(磁石表面の原点から垂直方向に1 mm の位置)
このプログラムで計算すると、エラーは出ないのですが、解が異常に小さく出たり、
分割数を変えると解の大きさが1桁単位で変化したりと、明らかに間違っており、
おそらくどこか基礎的なことで躓いていることは予想できるのですが、
具体的にどこが間違っているのかが分かりません。
見識ある方のお力添えを頂きたく存じます。
lang
1#define _USE_MATH_DEFINES 2#include<stdio.h> 3#include<math.h> 4#define f(ys0,zs0) ((x-xs0)*M)/(4*M_PI)/((x-xs0)*(x-xs0)+(y-ys0)*(y-ys0)+(z-zs0)*(z-zs0))/sqrt((x-xs0)*(x-xs0)+(y-ys0)*(y-ys0)+(z-zs0)*(z-zs0)) 5 6int main(void) { 7 int n, m, i; 8 double L1, L2, W1, W2, M, xs0, xs1, x, ys0, ys1, y, zs0, zs1, z, dy, dz, dSy1, dSy2, Sy1, Sy2, dV, V; 9 10 //要素分割数 11 n = 1000; 12 m = 1000; 13 //-- パラメータの設定 14 L1 = 0.000;//yの下限値(磁石の長さ) 15 L2 = 0.001;//yの上限値 16 W1 = 0.000;//zの下限値(磁石の幅) 17 W2 = 0.001;//zの上限値 18 M = 0.240;//磁石の磁界強度 [T] = [Wb/m^2] 19 xs0 = 0.000;//ソース点 20 x = 0.001;//目標点 21 y = 0.000;//目標点 22 z = 0.000;//目標点 23 //-- 24 //微小量の定義 25 dy = (L2 - L1) / n; 26 dz = (W2 - W1) / m; 27 //初期化 28 V = 0; 29 30 //-- 積分 31 for (i = 1; i < (n-1); i++) { 32 //y軸ソース位置の指定 33 ys0 = L1 + dy * i; 34 ys1 = ys0 + dy; 35 //ループ毎の初期化 36 Sy1 = 0; 37 Sy2 = 0; 38 for (i = 1; i < (m-1); i++) { 39 //z軸ソース位置の指定 40 zs0 = W1 + dz * i; 41 zs1 = zs0 + dz; 42 //y軸方向の微小面積計算 43 dSy1 = (f(ys0, zs0) + f(ys0, zs1)) * dz / 2; 44 dSy2 = (f(ys1, zs0) + f(ys1, zs1)) * dz / 2; 45 printf("dSy1 = %lf\ndSy2 = %lf\n", dSy1, dSy2); 46 //y軸方向の面積総和 47 Sy1 += dSy1; 48 Sy2 += dSy2; 49 } 50 //z軸方向の微小体積積分 51 dV = (Sy1 + Sy2)*dy / 2; 52 //z軸方向の体積総和 53 V += dV; 54 } 55 //-- 積分終了 56 57 printf("磁束密度 %lf [T]\n", V); 58 getchar(); 59}
分割数の変更、各パラメータの値の変更、台形公式の見直し、簡単な関数での重積分結果と比較など
Microsoft Visual Studio 2017 を使用
回答1件
0
ベストアンサー
forループがどっちもi
投稿2019/01/11 23:50
総合スコア13521
あなたの回答
tips
プレビュー
Q&A
