Sample Scripts from GB Books
GB007:  2次元誘電体

今度は 右図に示されるように、平行な電場の中に断面が円形のガラス棒を置いたときの場の変化について考察します。この場合にポイントとなるのはドメインの左側境界と右側境界において

すなわちポテンシャルの法線微分値をNatural文によって指定することになるわけですが、方程式が
div(D) = 0
という形式で与えられているため、U には ε がかかっていると考える必要があります。このため境界条件で法線微分値を指定する際にも ε がかかった形で指定することになるのでご注意ください。これは部分積分からもたらされる要請と言えます。

一つ前の dielectrics01b.pde をベースにスクリプトを用意します。
  TITLE
    'Glass Rod across a Parallel Field'    { dielectrics01c.pde }

  SELECT
    Errlim = 3e-5

  VARIABLES
    U                     { Electric potential }

偏微分方程式の定義に先立ち、パラメータ類をSI単位系で定義します。この段階では ε  の値は設定されていませんが、リージョン定義の中で具体的な値を指定しています。
  DEFINITIONS             { SI units }
    L = 1  r0 = 0.2
    eps0 = 8.854e-12      { Permittivity of vacuum }
    epsr2 = 7.0           { Relative permittivity of glass }
    eps                   { Declared only }
    Ex0 = 1e4             { Far field }
    Ex = -dx(U)  Ey = -dy(U)  E = -grad(U)  Em = magnitude(E)
    D_x = eps*Ex  D_y = eps*Ey  D = eps*E  Dm = magnitude(D)
                          { Dielectric displacement }

次に方程式を定義します。 解析対象のドメイン内には真電荷は存在しないので右辺は0となります。
  EQUATIONS
    div(D) = 0            { div(eps*(-grad(U)) = 0 }

境界の形状と境界条件を設定します。法線微分値のみの指定では従属変数値が一意に決まらないため（GB001参照）、下辺上でPoint valueの指定を入れている点に注意してください。また左辺の法線微分値はドメインから見て内向きとなるため負号を付ける必要があります。ここでは'domain'と'glass'という2つのリージョンを設定、その各々で別個の誘電率を規定していますが、これらリージョン名は境界上での積分計算を行う際に重要な意味を持ってきます。
  BOUNDARIES
    Region 'domain' eps = eps0    { Vacuum outside }
      Start 'outer' (-L, -L)
        Natural(U) = 0         Line to (0, -L) Point value(U) = 0
                               Line to (L, -L)
        Natural(U) = eps*Ex0   Line to (L, L)
        Natural(U) = 0         Line to (-L, L)
        Natural(U) = -eps*Ex0  Line to Close

    Region 'glass' eps = epsr2*eps0    { Glass }
      Start 'circle' (r0, 0) Arc(Center = 0,0) Angle = 360

最後に出力すべき情報を規定します。
  PLOTS
    Grid(x, y)
    Contour(U)
    Elevation(U) on 'outer'
    Elevation(normal(D)) on 'outer' Report(eps0*Ex0)
    Vector(D) norm    Vector(D) on 'glass'
    Elevation(tangential(E)) on 'circle' on Region 'domain'
    Elevation(tangential(E)) on 'circle' on Region 'glass'
    Elevation(normal(D)) on 'circle' on Region 'domain'
    Elevation(normal(D)) on 'circle' on Region 'glass'
    Elevation(normal(E)) on 'circle' on Region 'domain'
    Elevation(normal(E)) on 'circle' on Region 'glass'

  END

(1) Grid(x, y)
FlexPDEによって自動生成された有限要素メッシュの構成を示しています。

(2) Contour(U)
等電位線を表す U に関する等高線図です。

(3) Elevation(U) on 'outer'
外周境界上で U の値がどう変化しているかを示しています。左右の境界上における電位はガラス棒の存在によりわずかながら変化している点にご注意ください。

(4) Elevation(normal(D)) on 'outer' Report(eps0*Ex0)
電束密度ベクトル D の外周上での法線成分の値をプロットしたものです。左右の境界上で取るべき値はReport文によりプロット下部に出力されています。Natural文で指定した境界条件が正しく反映されていることが確認できます。

(5) Vector(D) norm
電束密度 D のベクトルプロットです。左から右への流れが表現されています。

(6) Vector(D) on 'glass'
誘電体内での D に関するベクトルプロットです。全体として平行なベクトル場となっており、またベクトルの大きさもほぼ一様であることが右の表からわかります。

(7) Elevation(tangential(E)) on 'circle' on Region 'domain'
内部境界 'circle' 上で電場ベクトル E の接線方向成分の値をプロットしたものです。ただし境界線はリージョン 'domain' とリージョン 'glass' 双方に属するため、そのどちらに視点を置くかによって符号が変わってきます。(8)は 'glass' 側でプロットしているためグラフは反転しています。

(8) Elevation(tangential(E)) on 'circle' on Region 'glass'
内部境界 'circle' 上で電場ベクトル E の接線方向成分の値をリージョン 'glass' 側からプロットしたものです。(7)の符号を反転させれば(8)と重なったものとなります。これは
電場の接線成分は境界面において連続である
という理論上の要請に基づくものである点に注意してください。

(9) Elevation(normal(D)) on 'circle' on Region 'domain'
内部境界 'circle' 上で電束密度ベクトル D の法線方向成分の値をリージョン 'domain' 側からプロットしたものです。

(10) Elevation(normal(D)) on 'circle' on Region 'glass'
内部境界 'circle' 上で電束密度ベクトル D の法線方向成分の値をリージョン 'glass' 側からプロットしたものです。(9)の符号を反転させれば(10)と重なったものとなります。これは
電束密度の法線成分は境界面において連続である
という理論上の要請に基づくものである点に注意してください。

(11) Elevation(normal(E)) on 'circle' on Region 'domain'
同じ法線成分であっても電場ベクトル E の場合には境界の両側で連続とはなりません。(12)とはスケールが異なる点に注意してください。

(12) Elevation(normal(E)) on 'circle' on Region 'glass'
電場ベクトル E の法線成分の値をリージョン 'glass' 側からプロットしたものです。符号を反転させてもスケールが異なるため、(11)のグラフとは重なりません。