数学とか語学とか楽しいよね

フランス語、ドイツ語、ロシア語、アラビア語、オランダ語、英語、スペイン語、ラテン語とか数学とか数値計算(有限要素法、有限体積法、差分法、格子ボルツマン法、数理最適化、C++コード付き)とか勉強したことをまとめます。右のカテゴリーから興味のある記事を探してください。最近はクラシックの名演も紹介しています。noteにも書いています。https://note.mu/kumag_min

【浅水流方程式】1次元浅水流方程式の導出(勾配と摩擦なし)

【浅水流方程式】サイトマップ(ここから関連記事が探せます)


1次元浅水流方程式(Shallow Water Equation, SWE)の導出(勾配と摩擦なし)をします。取りあえず導出すべき1次元浅水流方程式の基本形を書いておきます。要は慣れです。慣れると怖くなくなります。

水平な長方形の一様断面を考えることにします。今回は勾配と摩擦は無しにしましょう。浅水流方程式においては非圧縮性が仮定されています(密度が時間・場所変化しない)。まず、水深と流速で記述した場合は

 \displaystyle
\begin{eqnarray} 
&\frac{\partial h}{\partial t} + \frac{\partial (hu)}{\partial x} = 0 \\
&\frac{\partial (hu)}{\partial t} + \frac{\partial}{\partial x} \left( hu^2 + \frac{1}{2}gh^2 \right) = 0
\end{eqnarray}

となります。ここで  t は時間、 x は1次元座標、 h は水深、 u は流速、 g は重力加速度です。上の式が連続式で、下の式が運動量方程式です。次に水深と流量を用いて記述する場合です。水深と流量の関係

 \displaystyle
\begin{eqnarray} 
&q = uh
\end{eqnarray}

を用いて上式を変形すると

 \displaystyle
\begin{eqnarray} 
&\frac{\partial h}{\partial t} + \frac{\partial q}{\partial x} = 0 \\
&\frac{\partial q}{\partial t} + \frac{\partial}{\partial x} \left( \frac{q^2}{h} + \frac{1}{2}gh^2 \right) = 0
\end{eqnarray}

となります。こちらもどうぞ。


さて変数は以下の図のように配置されているとします。

f:id:mutsumunemitsutan:20181218222147j:plain:w400

手前から奥へ向けて( x 方向へ)水が流れています。各地点の流速  u(x) と水深  h(x) が与えられています。また、水路の幅は  B とします。一様断面なので  B が一定であることに注意してください。区間  [x, x+\Delta x] を考えます。

まず連続式からいきましょう。要するにこの式は質量保存則をあらわしています。区間  [x, x+\Delta x] における質量の保存を考えます。保存則は

考えている範囲の時間変化量 =
単位時間あたりに入ってくる量 - 単位時間あたりに出ていく量

という関係式から導かれます。物理においてはだいたいこのパターンなのでよく覚えておいてください。まず、区間  [x, x+\Delta x] における質量の時間変化は、水の密度を  \rho として(非圧縮性を仮定しているので定数扱い)

 \displaystyle
\begin{eqnarray} 
&\frac{\partial}{\partial t} \left( h(x) B \Delta x \rho \right) \\
\end{eqnarray}

と書けます。次に  x から単位時間あたりに入ってくる質量は

 \displaystyle
\begin{eqnarray} 
&h(x) B u(x) \rho\\
\end{eqnarray}

であり、 x+\Delta x から単位時間あたりに出ていく質量は

 \displaystyle
\begin{eqnarray} 
&h(x+\Delta x) B u(x+\Delta x) \rho \\
\end{eqnarray}

と書けます。これは  x において以下の図のような水の塊が流れ込んでくるからです。

f:id:mutsumunemitsutan:20181218224124g:plain:w200

高さが  h(x) で、 x 方向の奥行が  u(x) の水の塊です。 x 方向の奥行が  u(x) になるのは、単位時間で考えているからです。流速 (m/s) × 1 (s) = 距離 (m) のようになっています。 x+\Delta x からは高さが  h(x+\Delta x) で、 x 方向の奥行が  u(x+\Delta x) の水の塊が出ていきます。さて、これで保存則の要素が揃いました。関係式に代入すると

 \displaystyle
\begin{eqnarray} 
&\frac{\partial}{\partial t} \left( h(x) B \Delta x \rho \right) = h(x) B u(x) \rho - h(x+\Delta x) B u(x+\Delta x) \rho \\
\end{eqnarray}

となります。定数である  B, \rho で両辺を割ってから右辺を左辺に移項し、 \Delta x で割ると

 \displaystyle
\begin{eqnarray} 
&\frac{\partial h(x)}{\partial t} + \frac{h(x+\Delta x) u(x+\Delta x)-h(x) u(x)}{\Delta x} = 0 \\
\end{eqnarray}

を得ます。 \Delta x \to 0 の極限をとると、微分の定義より

 \displaystyle
\begin{eqnarray} 
&\frac{\partial h}{\partial t} + \frac{\partial (hu)}{\partial x} = 0 \\
\end{eqnarray}

が得られます。これが連続式です!

次に運動量方程式を導出します。運動量方程式は運動量保存則のことです。



区間  [x, x+\Delta x] における質量の保存を考えます。

参考
E. F. Toro, Shock-Capturing Methods for Free-Surface Shallow Flows, 第2章
R. J. Leveque, Finite Volume Methods for Hyperbolic Problems, 第13章
R. Szymkiewicz, Numerical Modeling in Open Channel Hydraulics, 第1章

Shock-Capturing Methods for Free-Surface Shallow Flows

Shock-Capturing Methods for Free-Surface Shallow Flows

Finite Volume Methods for Hyperbolic Problems (Cambridge Texts in Applied Mathematics)

Finite Volume Methods for Hyperbolic Problems (Cambridge Texts in Applied Mathematics)