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

誘電体（dielectrics; 電場を加えたときに分極を生じる物質）を対象とする場合には電束密度（electric displacement）ベクトル

を含む関係式を扱う必要があります。ここに ε _r は比誘電率（relative permittivity）を意味します。この場合のマックスウェルの方程式は

で与えられることになります。ここに ρ_v は真電荷（free charge; 自由電荷）密度（分極から生じたものではない電荷密度）を表すパラメータです。ここで扱う問題においては ρ_v が0のケースのみを考えるので、FlexPDEに与える方程式は div(D)=0 となります。
数式(1)はポテンシャル U に関するポアッソン方程式

に展開されることになります。

右の図は同軸ケーブルの断面を示したものです（単位としてはSI単位系を使用）。外側のシールドの電位を0、中心線の電位を1として、その中間に位置するポリマー部分における場を調べることにします。なお、中心線の半径は 2e-3 m、ポリマーの比誘電率は 2.3 とします。
なお、電束密度ベクトルを D としたとき、そのx成分、y成分をDx, Dyと表現したくなりますが、FlexPDEの場合、これらは微分演算子を表す予約語（dx, dy）です。このため以下のスクリプト中では D_x, D_y という表現が使用されている点に注意してください。

まずタイトルを設定します。
  TITLE
    'Coaxial Cable'       { dielectrics01a.pde }

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

従属変数を定義します。
  VARIABLES
    U                     { Electric potential }

偏微分方程式の定義に先立ち、パラメータ類をSI単位系で定義します。なお、cap_ex という数式は円筒型コンデンサの容量を与える理論式であり、Elevationプロット中で使用されます。
  DEFINITIONS             { SI units }
    r1 = 2e-3  r2 = 1e-2  U1 = 1.0
    eps0 = 8.854e-12      { Permittivity of vacuum }
    eps = 2.3*eps0        { Permittivity of polymer }
    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 = 2*PI*eps/ln(r2/r1)  { Theoretical capacitance }

次に方程式を定義します。真電荷が存在しないケースを想定しているので右辺は0となります。
  EQUATIONS
    div(D) = 0            { No volume charge }

境界の形状と境界条件（共にDirichlet型）を設定します。内部境界内は解析対象ドメインから除外されます。
  BOUNDARIES
    Region 1
      Start 'outer' (r2, 0)
        Value(U) = 0   Arc(Center = 0,0) Angle = 360
      Start 'inner' (r1, 0)
        Value(U) = U1  Arc(Center = 0,0) Angle = 360  { Cut-out }

最後に出力すべき情報を規定します。
  PLOTS
    Grid(x, y)
    Contour(U)
    Vector(D) norm
    Contour(D_x)  Contour(D_y)
    Contour(Dm)  Surface(Dm)
    Elevation(Dm) on 'inner'
      Report(cap_ex) as 'Theoretical capacitance'
    Elevation(Dm) on 'outer'

  END

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

(2) Contour(U)
等電位線のプロットです。内周境界で U = 1、外周境界で U = 0 となっています。

(3) Vector(D) norm
電束密度 D のベクトル場をプロットしたものです。中心から放射状に広がるベクトル場が描かれています。normを指定しているため、ベクトル場の大きさは色で判断することになります。

(4) Contour(D_x)
電束密度ベクトル D のx成分 D_x の値を等高線図の形でプロットしたものです。

(5) Contour(D_y)
電束密度ベクトル D のy成分 D_y の値を等高線図の形でプロットしたものです。

(6) Contour(Dm)
電束密度ベクトル D の絶対値を等高線図の形でプロットしたものです。ポテンシャルの勾配が大きい内周部で大きな値となっています。

(7) Surface(Dm)
電束密度ベクトル D の絶対値を曲面図の形でプロットしたものです。

(8) Elevation(Dm) on 'inner' Report(cap_ex)
内周境界に沿って電束密度ベクトル D の絶対値がどう変化しているかをグラフ化したものです。概ね 6.3e-9 前後の値となっており、それを境界上で積分した値が 7.92e-11 と算出されています。

円筒形のコンデンサの容量 C は

で与えられます（スクリプト上のcap_ex）。この値はReport文の形で出力されており、7.95e-11 という値となっています。今、内側と外側の電極間の電位差は1Vであるので、同軸ケーブル単位長当りの真電荷量の理論値は Q_F = CU = C = 7.95e-11 と算出されます。

ここで、内周境界である半径 r1 の円筒状曲面（長さは単位長1） S についてガウスの定理を適用すると、上記数式(1)に注意して

が導かれます。一方、ベクトル D は円周に対して垂直ですから

として計算できるわけですが、この値は 7.92e-11 としてプロット下部に出力されています。FlexPDEによる|D|の計算値にはかなりの変動が見られるわけですが、それでも理論値 7.95e-11 に対しては誤差 0.4% 以内に収まっています。

(9) Elevation(Dm) on 'outer'
外周境界に沿って電束密度ベクトル D の絶対値がどう変化しているかをグラフ化したものです。概ね 1.263e-9 前後の値となっており、それを境界上で積分した値は 7.94e-11 と算出されています。

ここで示したスクリプト上では誘電体の誘電率（eps）として一定値を用いましたが、例えば中心からの半径によって誘電率が変化するケース等についてもわずかな変更で対処 できます。