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

2次元における流体の流れを分析するに際し、ここでは

• 非圧縮性の完全流体（粘性がない流体）

• 渦を持たない流れ

という理想的なケースを想定することにします。この場合には流れはポテンシャル流れとなり、速度ベクトルはスカラーのポテンシャル関数 Φ を用いて

のように表現することができます。一方、連続の式

（ρ は密度、v は速度ベクトル）は非圧縮性流体の場合、ρ =一定となるため、

に帰着されます。従って (1) を (3) に代入することによって次のラプラス方程式を得ることができます。

この型の流れの場合、渦ベクトル ω のz成分を計算してみると

となることから渦は存在しないことがわかります。一方、独立変数 p （圧力）についてはエネルギー保存則に相当するベルヌーイの定理

を使って算出することができます。ただし数式(5)では重力の働く方向を z 軸方向と仮定しています。

ここでは右図のよう に先で太さが狭まったチャネルにおける水平な流れを分析します。重力の働く方向はz軸方向であり、ここで解析対象とする2次元の流れには影響が及ばないため、ベルヌーイの定理における gz の項は無視します。

方程式(4)と(5)は連立させる必要はないため、EQUATIONSセクションで指定するのは(4)のみとし、圧力 p を計算する数式はDEFINITIONSセクション中に記述するというアプローチを取ります。

境界条件としてはチャネル左端 x = 0 におけるx方向の流速 vx0 = 3.0 と圧力 p0 = 1e5 を規定します（単位はSI系）。従属変数が Φ であることから境界条件の記述の中での vx0 の指定にはNatural文が使用されることになります。一方 x = 3 におけるチャネル右端においてはValue文により Φ = 0 という定数値を指定しています。定数値としては別に 0 でなくても良いわけですが、これによってこの境界上で

が保証されることになる点に注意してください。その他の境界上では Natural(Φ)=0 という設定になります。

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

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

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

関連するパラメータや数式を定義します。単位としてはSI単位系を使用します。p に関する数式は数式(5)に基づくものですが、2次元の流れなので gz の項は無視できます。div(v), curl(v) という数式は検証用に設定しています。
  DEFINITIONS                     { SI units }
    Lx = 1  Ly = 1
    coef = 0.5                    { Constriction coefficient }
    vx0 = 3.0                     { Velocity at input end }
    p0 = 1e5                      { Atmospheric pressure }
    dens = 1e3                    { Mass density }
    vx = dx(phi)  vy = dy(phi)    { Velocity components }
    v = vector(vx, vy)  vm = sqrt(vx^2 + vy^2)  { Speed }
    p = p0 + 1/2*dens*(vx0^2 - vm^2)            { Pressure }
    div_v = dx(vx) + dy(vy)       { Divergence, or div(v) }
    curl_z = dx(vy) - dy(vx)      { Vorticity, or curl(v) }

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

BOUNDARIESセクションでは境界形状の規定と同時に境界条件を設定します。チャネル左端部では流速規定型のため、自然境界条件 Natural(phi) = -vx0 という指定になります。ドメインから外に向かう方向が + である点に注意してください。一方、チャネル右端部ではValue文によるDirichlet型の境界条件設定となっています。その他の境界では Natural(phi) = 0 という指定を行っていますが、これらの境界上では法線方向の流れの成分が 0 であることが約束されます。
  BOUNDARIES
    Region 1
      Start 'outer' (0, Ly)
        Natural(phi) = -vx0 Line to (0, -Ly)         { In }
        Natural(phi) = 0    Line to (Lx, -Ly)
                            to (2*Lx, -Ly*coef) to (3*Lx, -Ly*coef)
        Value(phi) = 0      Line to (3*Lx, Ly*coef)  { Out }
        Natural(phi) = 0    Line to (2*Lx, Ly*coef) to (Lx, Ly)
                            to Close

最後に出力すべき情報を規定します。
  PLOTS
    Grid(x, y)
    Contour(phi)
    Vector(v) norm
    Contour(vm) painted
    Contour(p) painted
    Contour(p) zoom(1.5*Lx, 0, Lx, Ly)
    Surface(p) zoom(1.5*Lx, 0, Lx, Ly)
    Elevation(vm) on 'outer'      { Verify boundary conditions }
    Contour(div_v)
    Contour(curl_z)

  END

(1) Grid(x, y)
FlexPDEによって生成されたメッシュ構成を示しています 。メッシュ再構成は3回行われ、角の部分におけるメッシュ密度が高くなっています。

(2) Contour(phi)
速度ポテンシャル Φ に関する等高線図を示したものです。チャネルの側壁に対しては直交するパターンとなっており、また右端における Φ の値は指定値通り 0 となっています。

(3) Vector(v) norm
流体の速度場 v を表すベクトルプロットです。normを指定しているため速度の大きさはカラーによって表現されています。プロット(2)の曲線群とは直交する形の流線が描かれています。チャネル出口における流速は入口におけるそれの概ね倍となっていることがカラーコードから読み取れます。

(4) Contour(vm) painted
速度ベクトル v の絶対値についての色塗り等高線図です。完全流体ゆえ側壁の近傍であっても速度が落ちることはありません。|v| の最大、最小値が境界上のどこで起きているかについてはプロット(8)を参照ください。

(5) Contour(p) painted
圧力 p に関する色塗り等高線図です。入口部と出口部での圧力の違いはそれほど大きくはありません。

(6) Contour(p) zoom(1.5*Lx, 0, Lx, Ly)
圧力 p に関する等高線図を (x, y) = (2, 0.5) の近辺で拡大したものです。流速 |v| が最大となる点において圧力は最小となります。

(7) Surface(p) zoom(1.5*Lx, 0, Lx, Ly)
プロット(6)と同じ内容を3次元曲面の形で表示したものです。

(8) Elevation(vm) on 'outer'
速度ベクトルの大きさ |v| が外周上でどう変化しているかをプロットしたものです。最大値、最小値がどの箇所で起るかが明確に示されています。また入口部、出口部における流速がそれぞれ 3.0, 6.0 であることがプロットから読み取れます。

(9) Contour(div_v)
ドメイン全域にわたり div(v) の値を計算したものです。ほぼ全域で 0 となっており、連続の式(3)を満足していることがわかります。

(10) Contour(curl_z)
ドメイン全域にわたり curl(v) （rot(v)）の値を計算したものです。全域で 0 となっており、渦のない流れであることが確認できます。