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

2. 圧力による流れ

今度は圧力によって流れが生ずるケースについて考察することにします。チャネルの形状はCase1の場合と同一ですが、入口側の圧力を p = δp, 出口側の圧力を p = 0 と設定することにします。

この単純な問題に関しては次のような解析解が知られています。
　　
ここに l はチャネル長を、2w はチャネル幅を意味します。FlexPDEによる計算結果をこの厳密解と比較するコードも以下のスクリプト中には含めてあります。

2.1 Problem descriptor [ vfluid01b.pde ]

まずタイトルを設定します。
TITLE
'Pressure-Driven Flow through a Channel' { vfluid01b.pde }

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

従属変数の設定は基本的にはCase1の場合と同一ですが、vy に対してはThresholdというクローズを設けています。FlexPDEの演算ロジックはErrlimで指定された相対誤差精度を達成すべく機能しますが、たまたまこの問題では vy = 0 が解となるため相対誤差の判定が有効に機能せず、メッシュ生成が際限なく繰返されるという問題が発生します。Thresholdクローズによる絶対誤差の明示はこの問題を回避するための措置と言えます。
VARIABLES
    vx
    vy(Threshold = 1e-5)
    p

関連するパラメータや数式を定義します。単位としてはSI単位系を使用します。
Case1ではReynolds数の計算（資料 GB012a.pdf 中の数式(5)）に vx0 という定数値を用いることができましたが、本ケースではそれが利用できないため、演算結果の中から vx の最大値を抽出（globalmax）、その値を用いてReynolds数を算出するというアプローチを取ります。
DEFINITIONS
    Lx = 1.0 Ly = 1.0
    visc = 1e4                                { Viscosity }
    delp = 100                                { Driving pressure }
    vx_ex = delp/(2*Lx)/(2*visc)*(Ly^2 - y^2) { Exact solution }
    dens = 1e3                                { Mass density }
    Re = dens*globalmax(vx)*2*Ly/visc         { Reynolds number }
    v = vector(vx, vy) vm = magnitude(v)     { Speed }

資料 GB012a.pdf 中に記載されている方程式(12), (13)を定義します。なお、方程式(13)中の発散項についてはCase1の場合と同様省略します（C = 0）。
EQUATIONS { Tagged by dominant variable } { For vanishing Re }
    vx: dx(p) - visc*div(grad(vx)) = 0
    vy: dy(p) - visc*div(grad(vy)) = 0
    p:   div(grad(p)) = 0    { Divergence term neglected }

BOUNDARIESセクションでは境界形状の規定と同時に境界条件を設定します。入口側では p = δp, 出口側では p = 0 という指定を行っています。入口、出口部におけるNatural(vx) = 0 という指定は
　　
を意味しています。またドメイン上端、下端におけるNatural(p)の指定は資料 GB012a.pdf 中の数式(14)に基づくものです。
BOUNDARIES
    Region 1
      Start 'outer' (-Lx, Ly)
        Natural(vx) = 0 Value(vy) = 0 Value(p) = delp
          Line to (-Lx, -Ly)                   { In }
        Value(vx) = 0 Value(vy) = 0 Natural(p) = -visc*div(grad(vy))
          Line to (Lx, -Ly)                    { Lower }
        Natural(vx) = 0 Value(vy) = 0 Value(p) = 0
          Line to (Lx, Ly)                     { Out }
        Value(vx) = 0 Value(vy) = 0 Natural(p) = visc*div(grad(vy))
          Line to Close                        { Upper }

最後に出力すべき情報を規定します。
PLOTS
    Grid(x, y)
    Elevation(vx, vy) on 'outer' Report(Re)
    Contour(vx)
    Contour(vx - vx_ex) Report(globalmax(vx))
    Contour(vy)
    Contour(p)
    Vector(v) norm
    Contour(div(v))
    Contour(curl(v)) painted

END

2.2 実行結果

(1) Grid(x, y)
FlexPDEによって生成されたメッシュ構成を示しています。問題が単純であるためメッシュの再構成は行われていません。

(2) Elevation(vx, vy) on 'outer' Report(Re)
vx, vy の値を境界上でプロットしたものです。vx の値はチャネル壁に接している部分で 0、チャネル中央部で最大となっています。変化の形状は放物線状ですが、これは厳密解の数式に符合しています。一方、vy の値は境界上で一律 0 です。なおReynolds数は 5e-4 という小さな値であることが確認できます。

(3) Contour(vx)
x方向の流速 vx に関する等高線図です。チャネル壁に接している部分で 0、チャネル中央部で最大となっていることがわかります。

(4) Contour(vx - vx_ex) Report(globalmax(vx))
FlexPDEが算出した vx の値を理論値と比較したものです。Scale = E-11という表示に見られるように厳密解との良い一致が見られます。なお vx の最大値は 2.5e-3 である点に注意してください。