Sample Scripts from GB Books
GB008:  2次元静磁場

マックスウェルの方程式の一つは

であるわけですが、静磁場ということでは単に

を考えれば良いことになります。ここに H は磁場の強度、J は電流密度、D は電束密度を意味します。これを直交座標系での成分で表すと次のようになります。

今、磁場としてx成分とy成分のみを持つケースを考えることにするとJx, Jyが消えて

のみが残る形となります。一方、磁束密度を B としたとき、ベクトルポテンシャル A(x, y, z) を考えることができて

と表現できます。ただし(x, y)平面に平行なベクトル場 H と B を考えようとしているわけですから、Ax, Ay 共に0という条件が付く形となります。B = μH に配慮すると（μ は透磁率）、数式(5)から

が導かれます。このHx, Hy を数式(4)に代入すると Az に関する2階の偏微分方程式

が導出されてきます。FlexPDEに対して指定する方程式はもちろん(7)の形でも良いのですが、DEFINITIONSセクションで(6)の関係を規定してやれば、EQUATIONSセクションで設定する方程式は(4)の形でも構いません。ポテンシャル Az が求まればベクトル場 H と B は(6)によって算出できます。

右の図の灰色の円形は電線の垂直断面を表したものです。電流密度は一様で Jz = 1 としたときの磁場の様子を円形のドメイン上で計算してみます。ただし電線の長さは十分に長いものとします。また透磁率としては電線周囲、電線内（銅製）、共に真空の透磁率 μ₀ を使用します。

この場合、電線の内外で厳密解が知られているため、FlexPDEによる計算結果との対比も行うことにします。

まずタイトルを設定します。
  TITLE
    'Field around a Wire'    { magnetics01a.pde }

次に演算精度に関するセレクタをセットします。
  SELECT
    Errlim = 3e-4

従属変数としてはベクトルポテンシャル A のz成分を使用します。
  VARIABLES
    Az                       { Magnetic vector potential }

偏微分方程式の定義に先立ち、パラメータ類をSI単位系で定義します。なお、B のx成分、y成分は B_x, B_y と表記されている点に注意してください。これは予約語 by とのコンフリクトを避けるためです。
  DEFINITIONS                { SI units }
    r0 = 0.2  r1 = 1.0
    rad = sqrt(x^2 + y^2)    { Radius }
    mu0 = 4*PI*1e-7          { Permeability of vacuum }
    mu = mu0                 { Permeability }
    Jz0 = 1.0                { Jz in wire }
    B_x = dy(Az)  B_y = -dx(Az)
    B = Vector(B_x, B_y)  Bm = magnitude(B)
                             { Magnetic flux density }
    Hx = B_x/mu  Hy = B_y/mu  H = B/mu  Hm = Bm/mu
                             { Magnetic field strength }
    Jz                       { Declared only }
    Bm_ex                    { Exact solution }

次に方程式を上記数式(4)の形で定義します。
  EQUATIONS
    dx(Hy) - dy(Hx) = Jz     { 2nd order PDE in Az }

境界の形状と境界条件を定義します。なお、上記パラメータJz, Bm_ex（厳密解）についてはリージョンごとに値を設定します。
  BOUNDARIES
    Region 1  Jz = 0
              Bm_ex = mu*Jz0*r0^2/(2*rad)
      Start 'outer' (-r1, 0)
        Value(Az) = 0 Arc(Center = 0,0) Angle = 360

    Region 'wire'  Jz = Jz0
                   Bm_ex = mu*Jz0*rad/2
      Start (-r0, 0) Arc(Center = 0,0) Angle = 360

ここで外周上で Value(Az) = 0 という境界条件を与えた意味について補足しておきます。境界上でも

という関係式は成り立っているわけですが、その成分表示式

においてx座標として接ベクトル t を、y座標として法線ベクトル n を取ることにすれば、

という式が導かれます。ここに τ は接線方向の座標変数を意味します。今、外周上における Az の値として定数値を設定したわけですが、このことは最後の式より Bn = 0 という条件を設定したことと等価であることがわかります。

最後に出力すべき情報を規定します。
  PLOTS
    Grid(x, y)
    Elevation(normal(B)/Bm) on 'outer'
    Surface(Az)
    Vector(B) norm zoom(-2*r0, -2*r0, 4*r0, 4*r0)
    Contour(Bm)  Surface(Bm)
    Contour(Bm_ex)
    Contour(Bm - Bm_ex) Report(Globalmax(Bm))
    Elevation(Bm, Bm_ex) from (-r1, 0) to (r1, 0)

  END

(1) Grid(x, y)
FlexPDEによって生成されたメッシュ構成を示しています。メッシュ再構成は1回行われ、電線との境界周辺のメッシュが細分化されています。

(2) Elevation(normal(B)/Bm) on 'outer'
外周上で磁束密度ベクトル B の法線成分が0に近い値になっていることを確認するためのプロットです。

(3) Surface(Az)
ベクトルポテンシャル Az の関数形を曲面図の形で表示したものです。

(4) Vector(B) norm zoom(-2*r0, -2*r0, 4*r0, 4*r0)
電線の近傍部分における磁束密度 B のベクトルプロットを示したものです。フィールドの向きはどこでも動径に直交する形となっています。またその強さは色で表現されていますが、電線の表面上で最大となっていることがわかります。

(5) Contour(Bm)
磁束密度ベクトル B の絶対値に関する等高線図です。電線表面上で最大となっています。

(6) Surface(Bm)
プロット(5)と同じ内容を曲面図の形で表示したものです。電線内部に逆円錐の形状の曲面が隠れているのですが、それについてはプロット(9)を参照ください。

(7) Contour(Bm_ex)
ここで考察している単純なケースについては厳密解が知られています。次の図は厳密解を与える数式をもとに磁束密度ベクトル B の絶対値に関する等高線図を描いたものです。実質(5)と変わりませんが、差分についてはプロット(8)も参照ください。

厳密解は電線内外の円周上で

という周回積分を行うことによって求められます（アンペールの法則）。電線外の場合は

に注意すると

が厳密解を与える数式となります。一方、電線内の点に対しては

という式を使わなくてはならないので、厳密解の数式は次のようになります。

(8) Contour(Bm - Bm_ex) Report(Globalmax(Bm))
FlexPDEによって算出された Bm の値と理論値との差を等高線図の形でプロットしたものです。Bm の最大値はReport文からの出力で 1.30e-7 であることから、その相対誤差は 4.30e-9/1.30e-7, すなわち3.3%以下におさまっていることがわかります。

Note: SELECTセクションで指定した Errlim はポテンシャル Az に対するものである点に注意。磁束密度ベクトル B は Az に対する偏微分操作で得られるため、数値計算に伴う誤差は大きくなっても不思議ではありません。

(9) Elevation(Bm, Bm_ex) from (-r1, 0) to (r1, 0)
右下に示されている横断路に沿って Bm の値とその理論値とをプロットしたものです。2本の曲線が描かれているわけですが、実質重なってしまっているため、1本のグラフのように見えています。積分値にもほとんど差がありません。