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

今度は 平行板コンデンサ内、及び周囲の場について考察します。右の図は2枚の金属板の間に誘電体がはさまった形の平行板コンデンサの断面図です。上側のプレートには1.0e4Vの電位をかけ、下側のプレートは接地するものとします。2mx2m のドメイン周囲は絶縁境界であるとして計算を実行します。

なおプレートの横幅は0.5m、プレート間の幅は0.2m、誘電体の比誘電率は7.0と設定します。

基本形は dielectrics01a.pde と変わりません。
  TITLE
    'Parallel Plate Capacitor'    { dielectrics01b.pde }

  SELECT
    Errlim = 3e-4

  VARIABLES
    U                     { Electric potential }

偏微分方程式の定義に先立ち、パラメータ類をSI単位系で定義します。この段階ではデフォルト値の意味で ε = ε₀ と設定していますが、誘電体のリージョン内においては ε = 7.0ε₀ という形でパラメータ値の変更を行います。なお、結果の検証用に平行板コンデンサの容量を計算する数式 cap_ex を定義しておきます。この数式はElevationプロット中で使用されます。
  DEFINITIONS             { SI units }
    L = 1  xx = 0.25  d0 = 0.1  dd0 = 0.4*d0
    U0 = 1.0e4            { Voltage }
    eps0 = 8.854e-12      { Permittivity of vacuum }
    eps = eps0            { Permittivity, default }
    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 }
    cap_ex = 7.0*eps0*(2*xx)/(2*d0)  { Theoretical capacitance }

次に方程式を定義します。 解析対象のドメインからは上下の金属板はexcludeされている点に注意してください。このドメイン内において真電荷は存在しないので右辺は0となります。
  EQUATIONS
    div(D) = 0            { No total charge }

境界の形状と境界条件を設定します。上下の金属板は解析対象ドメインから除外されます。 一方、誘電体部分は別個のリージョンとして定義し、固有の誘電率を設定します。
  BOUNDARIES
    Region 1
      Start 'outer' (-L, -L)  Natural(U) = 0  { Outer boundary }
        Line to (L, -L) to (L, L) to (-L, L) to Close
      Start 'upper' (-xx, d0)  Value(U) = U0  { Upper plate }
        Line to (xx, d0) to (xx, d0 + dd0) to (-xx, d0 + dd0) to Close
      Start 'lower' (-xx, -d0 - dd0)  Value(U) = 0.0
                                              { Lower plate, ground }
        Line to (xx, -d0 - dd0) to (xx, -d0) to (-xx, -d0) to Close

    Region 'glass' eps = 7.0*eps0    { Dielectric }
      Start (-xx, -d0)  Line to (xx, -d0) to (xx, d0)
                             to (-xx, d0) to Close

最後に出力すべき情報を規定します。
  PLOTS
    Grid(x, y)
    Contour(U)
    Vector(D) norm
    Contour(Dm)  Surface(Dm)
    Contour(Dm) log
    Contour(Dm) painted on 'glass'
    Elevation(D_y) from (-L, 0) to (L, 0)
      Report(cap_ex) as 'Theoretical capacitance'

  END

(1) Grid(x, y)
指定された演算精度を達成するためにメッシュの細分化が2度ほど実施されたことが“Grid#3”という表示からわかります。

(2) Contour(U)
等電位線のプロットです。金属板にはさまれた誘電体内では等電位線がほぼ平行に走っていることが確認できます。なお、周囲は絶縁境界なので等電位線は境界と直交する形となっています。

(3) Vector(D) norm
電束密度 D のベクトルプロットです。上の金属板から出た電気力線はすべて下の金属板に流れ込んでいる様子が示されています。

(4) Contour(Dm)
電束密度ベクトル D の絶対値を等高線図の形でプロットしたものです。ドメイン内の大半の領域で|D|の値が一定であるために、等高線では表現しにくい関数形となっています。

(5) Surface(Dm)
電束密度ベクトル D の絶対値を曲面図の形でプロットしたものです。誘電体内、及びその周囲で平面に近い形状となっているため等高線での表現には不向きと いうことがわかります。

(6) Contour(Dm) log
電束密度ベクトル D の絶対値に関する等高線図ですが、対数スケールでプロットしたものです。これでも誘電体内には等高線がほとんど現れてきません。

(7) Contour(Dm) painted on 'glass'
誘電体内にしぼって|D|の等高線図を描いたものです。外部と接する部分を除くとほぼ平坦な形状となっていることがわかります。

(8) Elevation(D_y) from (-L, 0) to (L, 0) Report(cap_ex)
電束密度ベクトル D のy成分 D_y の値がx軸上でどう変化しているかをプロットしたものです。積分値である -1.68e-6 という値については少々説明を加える必要があるでしょう。

下の金属板を含む周回路 Γ₁, Γ₂, Γ₃, Γ₄ に沿って D の法線成分を積分することを考えます（下図参照）。この場合、Γ₂, Γ₃, Γ₄ 上では法線成分が0となるため、積分に寄与するのは Γ₁ 上での積分のみとなります。

この周回路をz軸方向に単位長のばした直方体の表面を S としたとき、ガウスの定理と数式(1)より

が言えます（Q_F は金属板単位長当りの真電荷量）。一方、

より、FlexPDEが出力してきた -1.68e-6 という値が Q_F を意味することになります。Q = CU という等式にこの Q_F の値と U = 1.0e4 を代入するとコンデンサの容量が 1.68e-10 と算出されます。
一方、Report文で平行板コンデンサ容量の理論値を出力させていますが、その値は 1.55e-10 ということでFlexPDEの計算値とは若干異なります。これは後者が誘電体左右端での場の乱れを想定していないことによるものです。