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

前ケースと全く同じ問題ですが、今回は div(v) が 0 となるような解を誘導すべく、数式(13)（資料 GB012a.pdf 参照）中におけるパラメータCに0以外の値をセットして計算を行わせます。ただしCに対する適正値を見出すためには少々試行錯誤が必要となります。

前ケースのスクリプト vfluid01c.pde に対する変分のみを記すことにします。
  TITLE
    'Constricted Channel with Divergence Term'    { vfluid01d.pde }

方程式を定義します。変更が加わるのは第3式のみです。ここでは発散項に対して C = 1e3 という重みをかけています。
  EQUATIONS
    vx:  dx(p) - visc*div(grad(vx)) = 0
    vy:  dy(p) - visc*div(grad(vy)) = 0
    p:   div(grad(p)) - 1e3*visc/Ly^2*div(v) = 0

PLOTSセクションに対しては次の変更を加えます。最後の1文は追加です。
  PLOTS
    ...
    Contour(div(v))
    ...
    Elevation(natp) on 'outer'

  END

(1) Grid(x, y)
FlexPDEによって生成されたメッシュ構成を示しています 。メッシュの再構成は2回で済んでいます。

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

(3) Contour(vx) Report(Re)
x方向の流速 vx に関する等高線図です。 チャネル中央部における流速を見た場合、出口側での流速は入口側でのそれの約2倍となっていることがわかります。またReynolds数は 1.57e-4 という小さな値であることがレポートされています。

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

(5) Vector(v) norm
流速 v のベクトルプロットです。チャネル境界の近傍でベクトルが壁面に平行となっている点も改善点の一つです。

(6) Contour(p)
圧力 p の値をプロットした等高線図です 。曲線の形状がチャネルの中心線を軸に上下対称となっている点も改善点の一つです。

(7) Contour(div(v))
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プロットです。下部に示されている積分値 5.25e-4 に注意してください。これはこの断面を単位時間に通過した流量に対応する数値です。

(10) Elevation(vx) from (2.5*Lx, -Ly*coef) to (2.5*Lx, Ly*coef)
チャネル幅が半分になった部分における vx のelevationプロットです。下部に示されている積分値 5.23e-4 はプロット(9)の結果と良く一致しています。

(11) Elevation(natp) on 'outer'
境界上で

の値をプロットしたものです。関数式 natp には2階の微分項が含まれているためグラフは階段状の形状を呈しています。