Sample Scripts from GB Books
GB009:  3次元の電場

今度は 右図のような円錐台の形状をしたシリコン製の物体が与えられたとします。上面に 1V の電位をかけたときの電導の様子について、FlexPDEを用いて解析を行います。 その際、下面の電位は 0 に、周囲の斜面は絶縁状態に保たれているものとします。ただし、関連する寸法は次の通りです。
r1 = 1cm    r0 = 0.2cm    h = 1cm

この場合も extrusion の考え方で3次元オブジェクトの定義を行うわけですが、その場合、下面lower（下図で青の部分）の設定が少々わかりにくいかも知れません。円柱座標系で動径の長さを r としたとき、斜面を規定する方程式は

従って、lowerという面全体で言えば

という数式で定義できる点にご注意ください。

まずタイトルを設定します。
  TITLE
    'Electrical Conduction in a Cone'    { 3dfields1d.pde }

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

座標系が3次元直交座標系であることを明示します（デフォルトはCartesian2）。
  COORDINATES
    Cartesian3

従属変数を定義します。
  VARIABLES
    U

偏微分方程式の定義に先立ち、パラメータ類 や関係式を定義します。本ケースでは電流を伴うため、電流密度ベクトルJの規定が必要となります。 背景情報については GB003 を参照ください。
  DEFINITIONS
    r0 = 2e-3  r1 = 10e-3  h = 10e-3
    cond = 1.0e-3                   { Conductivity of silicon }
    Ex = -dx(U)  Ey = -dy(U)  Ez = -dz(U)
    E = -grad(U)  Em = magnitude(E)
    Jx = cond*Ex  Jy = cond*Ey  Jz = cond*Ez
    J = cond*E  Jm = magnitude(J)
    rad = sqrt(x^2 + y^2)           { Radius }
    power = Vol_Integral(Jm*Em)     { Dissipation }
    current = Surf_Integral(normal(J), 'top')

3次元の電導を表す偏微分方程式を定義します。背景情報については GB003 を参照ください。
  EQUATIONS
    div(J) = 0

z軸方向へのextrusionを定義します。実際には下面と上面の定義があるのみなので、その間の空間の特性は一様となります。なお、面'lower'の定義文中で使われている rad というのはDEFINITIONSセクションで定義されている通り、円柱座標系での動径を意味します。
  EXTRUSION
    Surface 'lower' z = max(0, h*(rad-r0)/(r1-r0))
                                    { Flat base plus cone }
    Surface 'top' z = h             { Flat upper surface }

BOUNDARIESセクションにおいてbase plane上でのリージョン構成定義と境界条件の設定を行います。Base planeがこの場合平面ではないため、やや変則的です。また半径が r0 の円形リージョンを定義していますが、これは下面'lower'内でこの部分のみ境界条件が異なるためです。
  BOUNDARIES
    Surface 'top' Value(U) = 1.0    { 1 volt applied on top }

    Region 'domain'                 { Total domain }
      Surface 'lower' Natural(U) = 0  { dU/dn=0 on cone }
      Start(r1, 0) Arc(Center=0,0) Angle = 360 to Close
                                    { Extension in (x, y) }

    Region 'cylinder'               { Overwrites total }
      Surface 'lower' Value(U) = 0  { Lower flat face }
      Start (r0, 0) Arc(Center=0,0) Angle = 360 to Close

最後に出力すべき情報を規定します。
  PLOTS
    Grid(x, y, z) Report(power) Report(current)
    Contour(U) on x = 0
    Vector(J) norm on x = 0
    Elevation(tangential(E))
      Line_Integrate from (0, r1/2, h) to (0, 0, 0)  { Voltage }

  END

(1) Grid(x, y, z) Report(power) Report(current)
FlexPDEによって生成された3次元のメッシュ構成を示しています。Report文によって


という2つの積分値が出力されている点に注意ください（S は電位のかかっている上面（'top'）を意味しています）。I₁はボリューム全体で散逸される電力を、I₂は上面から流入する電流の量を表しています。今、電位は 1 に設定してあるので、I₁ + I₂ = 0 というのが理論上の要請で ある点に注意してください。

(2) Contour(U) on x = 0
平面 x = 0 上における等電位線のプロットです 。これらの曲線は円錐斜面とは直交関係にあります。

(3) Vector(J) norm on x = 0
平面 x = 0 上における電流密度ベクトル J のベクトルプロットです。上面から下面に至る流れが表現されています。

(4) Elevation(tangential(E))
      Line_Integrate from (0, r1/2, h) to (0, 0, 0)
上面'top'上の1点から下面に至る直線に沿って電場ベクトル E の接線方向成分をelevationプロットの形で図にしたものです。同時に線積分を計算させていますが、FlexPDEによる計算値は 0.999919 と理論値である 1 （電位差）に十分近い値となっています。