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

これはGB011で扱った問題とほぼ同一の問題設定ではありますが、粘性の影響を考慮した方程式を使用する点が本質的に異なります。

ただし最初に発散項を含まない形（C = 0、資料 GB012a.pdf 中の数式(13)参照）での求解を試みます。結果的に誤った解が誘導されてしまう点に注意してください。

まずタイトルを設定します。
  TITLE
    'Viscous Flow through a Constricted Channel'    { vfluid01c.pde }

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

従属変数を規定します。
  VARIABLES
    vx  vy  p

関連するパラメータや数式を定義します。単位としてはSI単位系を使用します。
ここで任意境界面上における外向き単位法線ベクトルを n としたとき、そのx成分、y成分である n_x, n_y の算出方法について触れておきます。n の偏角を θ とすると

一方、x方向、y方向の単位ベクトルを e_x, e_y としたとき、これらに normal 演算子を施すことによって各々の n 方向の成分を求められるわけですが、

という関係式が成り立ちます。従って n_x, n_y は

のように算出できることがわかります。
  DEFINITIONS
    Lx = 1.0  Ly = 1.0
    coef = 0.5                    { Constriction coefficient }
    visc = 1e4                    { Viscosity }
    delp = 100                    { Driving pressure }
    dens = 1e3                    { Mass density }
    Re = dens*globalmax(vx)*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 }

資料 GB012a.pdf 中に記載されている方程式(12), (13)を定義します。なお、方程式(13)中の発散項については取り敢えず無視して（C = 0）計算を行ってみます。
  EQUATIONS
    vx:  dx(p) - visc*div(grad(vx)) = 0
    vy:  dy(p) - visc*div(grad(vy)) = 0
    p:   div(grad(p)) = 0

BOUNDARIESセクションでは境界形状の規定と同時に境界条件を設定します。入口側では p = δp, 出口側では p = 0 という指定を行っています。入口、出口部におけるNatural(vx) = 0 という指定は

を意味しています。またドメイン上端、下端におけるNatural(p)の指定は資料 GB012a.pdf 中の数式(14)に基づくものです。
  BOUNDARIES
    Region 1
      Start 'outer' (0, Ly)
        Natural(vx) = 0 Value(vy) = 0 Value(p) = delp    { In }
          Line to (0, -Ly)
        Value(vx) = 0 Value(vy) = 0 Natural(p) = natp    { Lower }
          Line to (Lx, -Ly) to (2*Lx, -Ly*coef) to (3*Lx, -Ly*coef)
        Natural(vx) = 0 Value(vy) = 0 Value(p) = 0       { Out }
          Line to (3*Lx, Ly*coef)
        Value(vx) = 0 Value(vy) = 0 Natural(p) = natp    { Upper }
          Line to (2*Lx, Ly*coef) to (Lx, Ly) to Close

最後に出力すべき情報を規定します。
  PLOTS
    Grid(x, y)
    Elevation(nx, ny) on 'outer' as 'direction cosines'
    Contour(vx) Report(Re)
    Contour(vm)
    Vector(v) norm
    Contour(p)
    Contour(div(v)) painted
    Contour(curl(v)) painted
    Elevation(vx) from (0.5*Lx, -Ly) to (0.5*Lx, Ly)
    Elevation(vx) from (2.5*Lx, -Ly*coef) to (2.5*Lx, Ly*coef)

  END

(1) Grid(x, y)
FlexPDEによって生成されたメッシュ構成を示しています 。メッシュの再構成が3回行われているため計算には10秒前後を要しますが、一応所定の演算精度を持つ解に到達することができます。

(2) Elevation(nx, ny) on 'outer' as 'direction cosines'
境界に沿って n_x, n_y の値がどう変化するかをプロットしたものです。

(3) Contour(vx) Report(Re)
x方向の流速 vx に関する等高線図です。出口側での流速が入口側よりも小さいという奇異な結果になっている点に注意してください。

(4) Contour(vm)
流速ベクトル v の絶対値に関する等高線図です。実質プロット(3)と変わりません。

(5) Vector(v) norm
流速 v のベクトルプロットです。

(6) Contour(p)
圧力 p の値をプロットした等高線図です。

(7) Contour(div(v)) painted
div(v)の値をプロットしたものです。非圧縮性の定常流の場合、これは全域で 0 でなくてはならないはずです。ここで得られた解が連続の式に違反していることがわかります。

(8) Contour(curl(v)) painted
curl(v)の値に関する等高線図です。

(9) Elevation(vx) from (0.5*Lx, -Ly) to (0.5*Lx, Ly)
チャネル幅が広い部分における vx のelevationプロットです。下部に示されている積分値 2.06e-3 に注意してください。これはこの断面を単位時間に通過した流量に対応する数値です。

(10) Elevation(vx) from (2.5*Lx, -Ly*coef) to (2.5*Lx, Ly*coef)
チャネル幅が半分になった部分における vx のelevationプロットです。下部に示されている積分値 2.54e-4 はプロット(9)の結果と明らかに異なっており、質量保存則（連続の式）が満たされていない解であることがこれからもわかります。