Sample Scripts from GB Books
GB011:  2次元の完全流体

ここでは側壁が平行なチャネル内に断面が円形の障害物を置いたときの流れの様子を分析します。ただし障害物の中心の位置は側壁から等距離にあるものとします。なお重力場の方向はz軸方向である点に注意してください。境界条件としては入口側で流速 vx0、圧力 p0 を、出口側で流速 vx0 を指定します。

まずタイトルを設定します。
  TITLE
    'Obstacle across a Straight Channel'    { pfluid01b.pde }

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

従属変数を定義します。ポテンシャル流れの場合には一つのスカラー変数（速度ポテンシャル）のみによって流れがすべて規定されます。
  VARIABLES
    phi                           { Velocity potential }

関連するパラメータや数式を定義します。正方形ドメインの各辺は2m、円の半径は0.2mです。
  DEFINITIONS                     { SI units }
    Lx = 1.0  Ly = 1.0  a = 0.2
    vx0 = 5.0                  { x-component of velocity at left end }
    p0 = 1e5                      { Atmospheric pressure at left end }
    dens = 1e3                    { Mass density }
    vx = dx(phi)  vy = dy(phi)    { Velocity components }
    v = vector(vx, vy)  vm = sqrt(vx^2 + vy^2)
    p = p0 + 0.5*dens*(vx0^2 - vm^2)  { Pressure }
    p1 = p0 + 0.5*dens*vx0^2      { Pressure at stagnation point }

ポテンシャル流れを規定するラプラス方程式を記述します。
  EQUATIONS
    div(grad(phi)) = 0

BOUNDARIESセクションでは境界形状の規定と同時に境界条件を設定します。外部境界上での境界条件がすべて微係数指定のNatural型であるため、これでは解が一意に決まらないため、Point value文によって値を固定している点に注意してください。障害物との境界上では Natural(phi) = 0 という指定になります。
  BOUNDARIES
    Region 1
      Start 'outer' (-Lx, Ly) Point value(phi) = 0
        Natural(phi) = -vx0 Line to (-Lx, -Ly)  { In }
        Natural(phi) = 0    Line to (Lx, -Ly)   { Wall }
        Natural(phi) = vx0  Line to (Lx, Ly)    { Out }
        Natural(phi) = 0    Line to Close       { Wall }

      Start 'obstacle' (a, 0)     { Cut-out }
        Natural(phi) = 0 Arc(Center = 0,0) Angle = 360

最後に出力すべき情報を規定します。
  PLOTS
    Grid(x, y)
    Contour(vm) painted
    Contour(p) painted
    Elevation(p) on 'obstacle' Report(p1)
    Vector(v) norm
    Vector(v) norm zoom(-3*a/2, -3*a/2, 3*a, 3*a)

  END

(1) Grid(x, y)
FlexPDEによって自動生成されたメッシュ構造を示しています。円柱周囲でメッシュの細分化が行われています。

(2) Contour(vm) painted
速度ベクトルの絶対値 |v| に関する等高線プロットです。円周の上端、下端部で流速が最大に、一方、円周の左端、右端では流速が 0 になっていることがわかります（淀み点）。

(3) Contour(p) painted
圧力 p に関する等高線図です。左右上下共に対称である点に注意してください。円柱に対しては上流側から働く力と下流側から働く力とが釣り合った状態となっているため、流れの方向に押す効果は存在しないことになります。

(4) Elevation(p) on 'obstacle' Report(p1)
障害物の境界上での圧力の値をプロットしたもので、上下左右の対称性が明解に示されています。左右端において p の値は p0（1e5）よりも大きな 1.25e5 という値を示しています。これはベルヌーイの定理から導かれる

という式において v = 0（淀み点）を代入したときの値（Report文出力参照）に等しいことが確認できます。逆に上端、下端部においては流速が大きいため、圧力は p0 よりはるかに小さな値となっています。

(5) Vector(v) norm
ドメイン全域での流線の様子を示すベクトルプロットです。速度の大きさは色で表現されています。

(6) Vector(v) norm zoom(-3*a/2, -3*a/2, 3*a, 3*a)
障害物の近傍部分における流線の様子を拡大表示したものです。