Section Quasi 2D Uniform Flow model case01 - 趣味で計算流砂水理 Computational Sediment Hydraulics for Fun Learning

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

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