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

2次元の流れでも粘性を考慮する場合にはNavier-Stokesの方程式を用いて計算を行う必要があります。関連する方程式の導出過程についてはPDF資料

• GB012a.pdf

を参照してください。

ここでは右図に示すように平行な壁面にはさまれた流体の動きについて考察します。x方向に向かって右側の壁は固定壁ですが、左側の壁は一定速度vx0で移動しているものとします。なお重力場の向きはz軸方向である点に注意してください。

1mのオーダのドメイン、1mm/secといった速度のケースでReynolds数

として小さな値を得るため、粘性係数 η の値が非常に大きな仮想的な流体をここでは想定しています。

まずタイトルを設定します。
  TITLE
    'Flow Due to a Moving Wall'    { vfluid01a.pde }

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

従属変数を定義します。ポテンシャル流れの場合には一つのスカラー変数によって流れが表現できましたが、N-S方程式を使用する場合には速度ベクトルの個々の成分（2次元の場合には vx と vy）と圧力 p を従属変数とする連立偏微分方程式を解く形になります。
  VARIABLES
    vx  vy  p

関連するパラメータや数式を定義します。単位としてはSI単位系を使用します。
GB011では p0 = 1e5 といった圧力（ambient pressure）をかけてきました。これは角の部分で圧力が負値とならないようにするための措置であったわけですが、N-S方程式においては p の微分項しか入ってこないため、p0 の指定は省略できます。ただし解の中で p の値が負となるような事態が発生した場合には適宜ambient pressureの値を指定するようにしてください。
  DEFINITIONS
    Lx = 1.0  Ly = 1.0
    vx0 = 1e-3                    { Speed of moving wall }
    visc = 1e4                    { Viscosity (η) }
    dens = 1e3                    { Mass density }
    Re = dens*vx0*2*Ly/visc       { Reynolds number }
    v = vector(vx, vy)  vm = magnitude(v)  { Speed }

資料 GB012a.pdf 中に記載されている方程式(12), (13)を定義します。その際、従属変数との対応付けを行う必要があるので注意してください。なお、方程式(13)中の発散項についてはここでは省略します（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セクションでは境界形状の規定と同時に境界条件を設定します。vx の値についてはドメイン下端で vx = 0、ドメイン上端で vx = vx0 という設定を行っている点に注意してください。またドメイン左端と右端で Natural(vx) = 0 という指定を行っていますが、これは

という条件を表しています。これに対し圧力については、流れが壁の移動によってもたらされることから、ドメイン左端と右端における圧力 p の値は 0 としてあります。一方、ドメイン上端と下端における圧力については GB012a.pdf 中の数式(14)に基づく自然境界条件を指定しています。ただし法線ベクトルはドメインの中から外へ向う方向に取られるため、上端境界上では n_y = 1, 下端境界上では n_y = -1 となる点に注意が必要です。
  BOUNDARIES
    Region 1
      Start 'outer' (-Lx, Ly)
        Natural(vx) = 0 Value(vy) = 0 Value(p) = 0
          Line to (-Lx, -Ly)                                { Left }
        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)                                  { Right }
        Value(vx) = vx0 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(vy)
    Contour(p)
    Vector(v) norm
    Contour(div(v))
    Contour(curl(v)) painted

  END

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

(2) Elevation(vx, vy) on 'outer' Report(Re)
vx, vy の値を境界上でプロットしたもので、境界条件で指定した内容が反映されていることが確認できます。なおReynolds数を計算していますが、2e-4 という小さな値となっている点に注意してください。

(3) Contour(vx)
x方向の流体速度 vx に関する等高線プロットです。ドメイン境界上端における流速 vx0 がリニアに減少して行き、下端では値が 0 となっていることがわかります。

(4) Contour(vy)
y方向の流体速度 vy に関する等高線プロットです。ドメイン全域で vy = 0 となっており、流れがx方向の層流であることがわかります。

(5) Contour(p)
圧力 p に関する等高線図です。ドメイン全域で p = 0 という結果になっています。

(6) Vector(v) norm
流体速度 v に関するベクトルプロットです。normを指定しているため、速度の大きさはカラーによって表現されています。

(7) Contour(div(v))
div(v) の値をプロットした等高線図です。Scale = E-16 という表示からもわかるように、全域で div(v) = 0 という条件が満足されていることがわかります。

(8) Contour(curl(v)) painted
curl(v) の値をプロットしたもので、ドメイン全域で -5e-4 という一定値となっていることがわかります。回転の定義式

をドメイン境界上に適用して得られる値

と矛盾しない点に注意してください。