趣味で計算流砂水理

趣味で計算流砂水理 Computational Sediment Hydraulics for Fun Learning

数値計算とか河川工学とかプログラミングのことを書いています

MENU

Section Quasi 2D Uniform Flow model case01

スポンサーリンク

またまた更新間隔が空いてしまいました.結構ひどい風邪を引いてしまい,週の半分も欠勤してしまうなどあり,なかなか手が回りませんでした.

Section Quasi 2D Uniform Flow modelの開発

前回打合せ時の件です.命名はとりあえず上記の通りにしました.先に結論を言うと結構長丁場になりそうです. 関連するデータは全てDBに入れました.

再現対象となる実験結果

冨永先生の以下の論文を対象としました.

長方形断面開水路流の三次元乱流構造に関する実験的研究

今回は矩形水路でアスペクト比1:8の断面における主流の面内分布を対象としました.論文中Fig4の一番上です.

支配方程式

  • 今回は定常,等流としているため,x方向の運動方程式からx微分項およびv,z成分をキャンセルしました.

 {
\frac{\partial u}{\partial t}+\frac{\partial uu}{\partial x}+\frac{\partial uv}{\partial y}+\frac{\partial uw}{\partial z} = X-\frac{1}{\rho }\frac{\partial p}{\partial x}+\nu_t \left( \frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2}+\frac{\partial^2 u}{\partial z^2} \right)
}

 {
0 = X+\nu_t \left(\frac{\partial^2 u}{\partial y^2}+\frac{\partial^2 u}{\partial z^2} \right)
}

離散化

  • 拡散項しか残らないのでシンプルに中央差分です.u(j,k)を左辺に持ってくるとこんな感じになりました.

{
0 = g I_b + \nu_t \left(\frac{u_{j+1,k}-2u_{j,k}+u_{j-1,k}}{\Delta y^2}+\frac{u_{j,k+1}-2u_{j,k}+u_{j,k+1}}{\Delta z^2}\right)
} {
u_{j,k} = \frac{g I_b\frac{\Delta y^2\Delta z^2}{\nu_t} +  \Delta z^2(u_{j+1,k}+u_{j-1,k})+\Delta y^2(u_{j,k+1}+u_{j,k+1}) }{2\Delta y^2+2\Delta z^2}
}

境界条件

  • 多分ここがポイントです.最初は壁面はノンスリップ条件,水面はスリップ条件としました.
  • 今後はここを変更していく形になります.

数値計算方法

  • 差分方法は中央差分,収束計算はシンプルなrelaxation法です.ここはアイデアを下さい.
  • ソースコード 支配方程式の部分だけ公開します.予定通りPython3.4です.境界条件を考慮するとすごく綺麗に書けました.
  • なお,ソースコードの内のerrorは収束判定です.収束判定は10の-12乗にしてますが30秒位で終わります.
  • ぜひそちらの環境で実行してみてください.
def NumericalSimulation(errorcheck=False):
    error = np.zeros((ny,nz))
    for j in range(0,ny):
        for k in range(0,nz):
            if j == 0 : # Bed Boundary: Non-Slip Condtion
                a1 = 3
                diffy = un[j+1][k]
            elif j == ny-1 : # Bed Boundary: Non-Slip Condtion
                a1 = 3
                diffy = un[j-1][k]
            else: # Normal
                a1 = 2
                diffy = un[j+1][k]+un[j-1][k]

            if k == 0 : # Bed Boundary: Non-Slip Condtion
                a2 = 3
                diffz = un[j][k+1]
            elif k == nz-1 : # WaterSurface Boundary: Slip Condtion
                a2 = 1
                diffz = un[j][k-1]
            else: # Normal
                a2 = 2
                diffz = un[j][k+1]+un[j][k-1]

            if errorcheck :
                error[j][k] =np.abs(9.8*ib*dy**2*dz**2/nyut+dz**2*(diffy-a1*un[j][k])+dy**2*(diffz-a2*un[j][k]))
            else :
                un[j][k] =(9.8*ib*dy**2*dz**2/nyut + dz**2*diffy + dy**2*diffz ) / ( a1*dz**2 + a2*dy**2 )
    if errorcheck :
        errormax = np.max(error)
    else :
        errormax = 0.0
    return errormax

結論

  • こんな感じです.論文のものと見比べて下さい.冨永先生のものを今度トレースしておきます.

f:id:SedimentHydraulics:20151025154046p:plain

  • ノンスリップ条件だと底面付近の流速が小さくなりすぎる.
  • 併せて,断面平均流速が小さくなりすぎる.
  • よって,もう少し滑らす必要があるため,次回は対数則でチャレンジします.

先週の半ばに急激にアクセスが伸びてますが何かあったのでしょうか?