choudendoutaienchu.nb

磁界中の超伝導円柱

超伝導円柱を平行磁界中に置くとマイスナー効果で反磁性を示し、内部に磁界を侵入させない。このときの磁位は次のように与えられる。
　　　　　　φ= cosθ　　　　　　　　　　　　　(1)
第一項は外部からの平行磁界であり、第二項は超伝導体による反磁性の項である。これよりこの磁位を可視化する。さらに磁界を求めて、これを可視化する。

定義

平行磁界=1が半径Rの超伝導円柱にかかっているときの磁位phiは以下のように与えられる

Clear[phi] ; R = 1 ; phi[r_, th_] := r Cos[th] + R^2/r * Cos[th] /; r>R ; phi[r_, th_] := 0 /; r <= R ;

円柱座標系から直交座標系に変換するための定義

th[x_, y_] := Which[x>0, ArcTan[y/x],
&nbs ... p; x == 0, 0] ;

磁位の表示

磁位はポテンシャルなので、高さを使って表現できる。まず三次元表示。

[Graphics:HTMLFiles/index_8.gif]

上の図を見ると原点付近の超伝導内部では磁位は0で一定であり、外部で均一な勾配がついていることが分かる。これを円柱座標系で表したのが式(1)である。式(1)からこの図を想像することができるだろうか。
次のようにParametricPlot3Dを使えば直接表示することもできる

fx = r Cos[th] ; fy = r Sin[th] ; ParametricPlot3D[{fx, fy, phi[r, th]}, {r, 0, 5}, {th, 0, 2 Pi}, PlotRange -> {{-2, 2}, {-2, 2}, {-2, 2}}]

[Graphics:HTMLFiles/index_11.gif]

PlotRangeがないと円盤が出てしまう。

[Graphics:HTMLFiles/index_14.gif]

次に等高線表示。三次元情報は等高線を使って表示することができる。

[Graphics:HTMLFiles/index_17.gif]

次に濃淡表示

[Graphics:HTMLFiles/index_20.gif]

磁界の表示

必要なパッケージを読み込む

まず、円柱座標系で勾配(grad)をとってみる

R = 1 ; Clear[phi2] ; phi2[r_, th_, z_] := +r Cos[th] + R^2/r * Cos[th] /; r>R ; phi2[r_, th_, z_] := 0 /; r <= R ;

${Er[r_, th_, z_], Eth[r_, th_, z_], Ez[r_, th_, z_]} = Grad[-phi2[r, th, z], Cylindrical[r, th, z]]$

${-phi2^(1, 0, 0)[r, th, z], -phi2^(0, 1, 0)[r, th, z]/r, -phi2^(0, 0, 1)[r, th, z]}$

次にそれを直交座標系に変換する。式は次のとおり
Bx = Br cosθ - Bθ sinθ
By= Br sinθ + Bθ cosθ

Ex[x_, y_] := Er[r[x, y], th[x, y], 0] Cos[th[x, y]] - Eth[r[x, y], th[x, y], 0] Sin[th[x, y]] ... Ey[x_, y_] := Er[r[x, y], th[x, y], 0] Sin[th[x, y]] + Eth[r[x, y], th[x, y], 0] Cos[th[x, y]] ;

ベクター表示 PlotVectorFieldでは直交座標系のみ扱うので変換して用いる。
実行には時間がかかるので注意。

[Graphics:HTMLFiles/index_31.gif]

中心がうまくいっていないが、外部で磁界がはじかれていることが分かる。磁位の傾き(gradient)がベクトルの向きに対応することがわかる。
ScaleFunctionをうまく使うとよい。この場合には長さを強制的に0.5にしている。ScaleFunctionを使わないと矢印の長さがそこでのベクトルの大きさを示すことになり、磁位の傾きが大きいところで矢印は長くなる。

次のように直接Br, Bthを計算した結果を入れて、ベクトルをプロットすることもできる。計算速度は速くなる。

R = 1 ; Clear[Br, Bth] Br[r_, th_] := -(1 - R^2/r^2) Cos[th] /; r>R ; Br[r_, th_] := 0 /;r <= R ; Bth[r_, th_] := (1 + R^2/r^2) Sin[th] /; r>R ; Bth[r_, th_] := 0 /; r <= R ;

Bx[x_, y_] := Br[r[x, y], th[x, y]] Cos[th[x, y]] - Bth[r[x, y], th[x, y]] Sin[th[x, y]] ; By[x_, y_] := Br[r[x, y], th[x, y]] Sin[th[x, y]] + Bth[r[x, y], th[x, y]] Cos[th[x, y]] ;

[Graphics:HTMLFiles/index_36.gif]

Created by Mathematica (April 23, 2004)