Sample Scripts from GB Books
GB012:  2次元の粘性流体

Case3, 4ではチャネル入口部において一定の圧力をかけた場合を想定し解析を行いましたが、ここでは流速値が一定という条件を課したときの流れについて計算を行ってみます。

まずタイトルを設定します。
  TITLE
    'Uniform Velocity of Injection'    { vfluid01f.pde }

次に演算精度に関するセレクタをセットします。Case3, 4に比べて演算精度を落としている点に注意してください。
  SELECT
    Errlim = 3*1e-3

従属変数を規定します。
  VARIABLES
    vx
    vy
    p                             { Pressure minus ambient }

関連するパラメータや数式を定義します。Case3, 4の場合とは全く異なる粘性係数値が使用されている点に注意してください。vx0 の値はReynolds数が小さな値を取るように考慮して設定する必要があります。Cの値は条件に合わせて試行錯誤が必要です。
  DEFINITIONS
    Lx = 1.0  Ly = 1.0
    coef = 0.5                    { Constriction coefficient }
    visc = 1.0                    { Viscosity }
    vx0 = 1e-5                    { Input velocity }
    dens = 1e3                    { Mass density }
    Re = dens*vx0*2*Ly/visc       { Reynolds number }
    v = vector(vx, vy)  vm = magnitude(v)  { Speed }
    unit_x = vector(1, 0)         { Unit vector fields }
    unit_y = vector(0, 1)
    nx = normal(unit_x)           { Direction cosines }
    ny = normal(unit_y)
    natp = visc*(nx*div(grad(vx)) + ny*div(grad(vy)))
                                  { Natural boundary condition for p }
    C = 1e3

方程式を定義します。
  EQUATIONS
    vx:  dx(p) - visc*div(grad(vx)) = 0
    vy:  dy(p) - visc*div(grad(vy)) = 0
    p:   div(grad(p)) - C*visc/Ly^2*div(v) = 0

BOUNDARIESセクションでは境界形状の規定と同時に境界条件を設定します。 圧力 p については基本的にその微分値（自然境界条件）で規定しますが、出口部についてはDirichlet型で p = 0 という指定を行っています。現実の圧力という意味では得られた p の値に大気圧を加える必要があります。
  BOUNDARIES
    Region 1
      Start 'outer' (0, Ly)
        Value(vx) = vx0 Natural(vy) = 0 Natural(p) = natp  { In }
          Line to (0, -Ly)
        Value(vx) = 0 Value(vy) = 0 Natural(p) = natp      { Wall }
          Line to (Lx, -Ly) to (2*Lx, -Ly*coef) to (3*Lx, -Ly*coef)
        Natural(vx) = 0 Natural(vy) = 0 Value(p) = 0       { Out }
          Line to (3*Lx, Ly*coef)
        Value(vx) = 0 Value(vy) = 0 Natural(p) = natp      { Wall }
          Line to (2*Lx, Ly*coef) to (Lx, Ly) to Close

最後に出力すべき情報を規定します。
  PLOTS
    Grid(x, y)
    Elevation(vx) from (0, -Ly) to (0, Ly)
    Elevation(vx, 0.1*dy(vx)) from (3*Lx, -Ly*coef) to (3*Lx, Ly*coef)
    Elevation(p) on 'outer'
    Vector(v) norm Report(Re)
    Contour(vx)
    Contour(vy)
    Contour(vm)
    Contour(p)
    Contour(div(v))
    Contour(curl(v)) painted

  END

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

(2) Elevation(vx) from (0, -Ly) to (0, Ly)
チャネル入口部における流速 vx の値をプロットしたものです。1e-5 という値がチャネル壁に接する部分で 0 に変化する形となっている点に注意してください。

(3) Elevation(vx, 0.1*dy(vx)) from (3*Lx, -Ly*coef) to (3*Lx, Ly*coef)
チャネル出口部分において vx とその導関数の値をプロットしたものです。導関数のグラフが直線であることから、vx の値の変化は2次関数に従っていることがわかります。

(4) Elevation(p) on 'outer'
圧力 p の値が外部境界上でどう変化しているかをプロットしたものです。チャネル出口部分では設定通り p = 0 となっていることが、一方、入口部分では p の値が大きく変化していることがわかります。

(5) Vector(v) norm Report(Re)
流速 v のベクトルプロットです。Reynolds数は 0.02 とレポートされています。

(6) Contour(vx)
流速 vx の値をプロットした等高線図です 。チャネル幅が細まったところでは流速が vx0 の3倍以上になっていることがわかります（チャネル中央部）。

(7) Contour(vy)
流速 vy の値をプロットした等高線図です。

(8) Contour(vm)
流速ベクトル v の絶対値に関する等高線図です。

(9) Contour(p)
圧力 p の値についての等高線図です。チャネル出口部では 0 であるのに対し、入口部分の壁に接する部分で最大の圧力がかかっています。

(10) Contour(div(v))
チャネル入口部の上下端では vx の値が vx0 から 0 に不連続に変化するという境界条件のため計算誤差が大きくなっていますが、その部分を除くと div(v) の値はほぼ 0 となっていることがわかります。

(11) Contour(curl(v)) painted
curl(v) の値の色塗り等高線図です。チャネル幅が狭まった部分で回転が加わっています。