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

今度は流れの中に平板を斜めに配置したときに、平板にかかる力について分析してみます。流れについてはこれまでと同一の条件、すなわち入口での流速 vx0、圧力 p0、出口での流速 vx0 を仮定します。板の流れに対する角度を α 、板の幅と厚みをそれぞれ a = 0.5, d = 0.1 としたとき、板の4隅の座標値は次のように与えられる点に注意してください。なお、ここでは α  = PI/6 のケースについて計算を行います。

まずタイトルを設定します。
  TITLE
    'Drag and Lift on a Plate'    { pfluid01d.pde }

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

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

関連するパラメータや数式を定義します。最後のbrute_forceという式にかかっている 2*y2 という式は y2 - y4 を意味し、板の流れの方向（x方向）への射影面積（z方向には1を取る）に対応します。
  DEFINITIONS                   { SI units }
    Lx = 1.0  Ly = 1.0  a = 0.5*Ly  d = 0.2*a
    vx0 = 5.0                   { x-component of velocity at left end }
    alpha = PI/6                { Angle of attack, radians }
    si = sin(alpha)  co = cos(alpha)
    x1 = -d/2*si - a/2*co       { Corner coordinates }
    y1 = -d/2*co + a/2*si
    x2 = d/2*si - a/2*co
    y2 = d/2*co + a/2*si
    x3 = -x1  y3 = -y1  x4 = -x2  y4 = -y2
    p0 = 1e5                    { Atmospheric pressure at left end }
    dens = 1e3                  { Mass density }
    vx = dx(phi)  vy = dy(phi)  { Velocity components }
    v = vector(vx, vy)  vm = magnitude(v)
    p = p0 + 0.5*dens*(vx0^2 - vm^2)  { Pressure }
    brute_force = p0*2*y2

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

BOUNDARIESセクションにおける設定のしかたはこれまでのケースと 同様です。
  BOUNDARIES
    Region 1
      Start 'outer' (-Lx, Ly) Point value(phi) = 0
        Natural(phi) = -vx0 Line to (-Lx, -Ly)
        Natural(phi) = 0    Line to (Lx, -Ly)
        Natural(phi) = vx0  Line to (Lx, Ly)
        Natural(phi) = 0    Line to Close

      Start 'obstacle' (x1, y1)    { Cut-out }
        Natural(phi) = 0    Line to (x2, y2) to (x3, y3)
                                 to (x4, y4) to Close

最後に出力すべき情報を規定します。
  PLOTS
    Grid(x, y)
    Contour(vm) painted
    Vector(v) norm
    Contour(p) painted
    Elevation(p) from (-Lx, -Ly) to (-Lx, Ly)  { Left }
      Report(brute_force)
    Elevation(p) from (Lx, -Ly) to (Lx, Ly)    { Right }
    Elevation(p) from (-Lx, -Ly) to (Lx, -Ly)  { Bottom }
    Elevation(p) from (-Lx, Ly) to (Lx, Ly)    { Top }

  END

(1) Grid(x, y)
FlexPDEによって自動生成されたメッシュ構造を示しています 。メッシュ再構成が2回行われ、板の両端部でメッシュ密度が高くなっています。

(2) Contour(vm) painted
速度ベクトルの絶対値 |v| に関する等高線プロットです。流速は(x₂, y₂), (x₄, y₄)の部分で最大に、(x₁, y₁), (x₃, y₃)の部分で最小になっていることがわかります。

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

(4) Contour(p) painted
圧力 p に関する等高線図です。流速が早い部分で圧力が低くなっています。平板に対しては時計方向の回転モーメントがかかっていることがわかります。

(5) Elevation(p) from (-Lx, -Ly) to (-Lx, Ly) Report(brute_force)
ドメイン左端における p の値をプロットしたもので、概ね1e5という値を示しています。出力されているbrute_forceという値は計算式からわかるように、流れが平板に及ぼす力の概数を示しています。p の線積分値 199982 という値は次のプロットで参照します。

(6) Elevation(p) from (Lx, -Ly) to (Lx, Ly)
ドメイン右端における p の値をプロットしたもので、概ね1e5という値を示しています。p の線積分値を左端と右端とで比較した場合、その差はbrute_forceのオーダに比べてnegligible smallと言えます。

(7) Elevation(p) from (-Lx, -Ly) to (Lx, -Ly)
ドメイン下端（チャネル右岸）における p の値をプロットしたものです。p の線積分値 198699 という値に注意してください。

(8) Elevation(p) from (-Lx, Ly) to (Lx, Ly)
ドメイン上端（チャネル左岸）における p の値をプロットしたものです。p の線積分値を下端と上端とで比較した場合、その差はbrute_forceのオーダに比べてnegligible smallと言えます 。

プロット(5)-(8)の線積分値の結果が示唆するところは、平板に対してdragやliftというnet forceが全く働いていないということです。実際ポテンシャルフローの場合にはこのような結果になることが解析的に示されています。これは常識と矛盾する結果であるわけですが、その原因は粘性を全く持たない完全流体を仮定したことに起因しています。より現実に即した結果を導くためには粘性に配慮したモデルを考える必要があります。