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

ここで紹介するスクリプトはGunnar Backstrom氏の承諾のもと、書籍 “Simple Fields of Physics by Finite Element Analysis” に記されている多数のFlexPDE適用事例 の中からその一部を紹介するものです。

PDF版 (849KB)

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

  

(1)

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

  

(2)

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

  

(3)

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

  

(4)

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

  

(5)

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

  

(6)

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

  

(7)

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

1. 電線周辺の磁場

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

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

1.1 Problem descriptor [ magnetics01a.pde ]

まずタイトルを設定します。
  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.2 実行結果

(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)

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

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

  

(9)

(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本のグラフのように見えています。積分値にもほとんど差がありません。

次へ

page_top_icon