浅水方程式

By M.Sato

August 5, 2026

1. 概要

水平方向に広く、鉛直方向に狭い流れは浅水流近似を使ってモデル化することができる。浅水流の解析には浅水方程式が用いられる。

浅水方程式では、鉛直方向の速度成分が水平成分に比べて十分小さいと仮定し、水深平均した水平方向の流速および水面高さを計算する。

ここでは、浅水方程式の導出過程を示し、 OpenFOAMの標準ソルバーshallowWaterFoamを用いて解析を行ってみた。

2. 浅水方程式の導出

2.1. 非圧縮性流体の基礎式

非圧縮性流体の連続方程式および運動量方程式は次式で表される。流れ場の水平スケールが鉛直スケールよりも十分大きい場合、鉛直方向の速度成分は水平成分に比べて小さいと考えられる。浅水方程式はこれらの基礎方程式を水深方向に積分することで導かれる。

$$ \begin{align} \frac{\partial u_i}{\partial x_i} &= 0 \notag \\ \frac{\partial u}{\partial t} + \frac{\partial u u_j}{\partial x_j} &= -\frac{\partial }{\partial x} \left(\frac{p}{\rho}+gz\right)+ \frac{\partial }{\partial x_j} \left(\frac{\tau_{xj}}{\rho}\right) \notag \\ \frac{\partial v}{\partial t} + \frac{\partial v u_j}{\partial x_j} &= -\frac{\partial }{\partial y} \left(\frac{p}{\rho}+gz\right)+ \frac{\partial }{\partial x_j} \left(\frac{\tau_{yj}}{\rho}\right) \notag \\ \frac{\partial w}{\partial t} + \frac{\partial w u_j}{\partial x_j} &= -\frac{\partial }{\partial z} \left(\frac{p}{\rho}+gz\right)+ \frac{\partial }{\partial x_j} \left(\frac{\tau_{zj}}{\rho}\right) \notag \end{align} $$

2.2. ライプニッツの積分則

浅水方程式の導出にはライプニッツの積分則を用いる。ライプニッツの積分則は次式で表される。

$$ \begin{align} \frac{d}{dx} \int_{a(x)}^{b(x)} f(x,y) dy = \int_{a(x)}^{b(x)} \frac{\partial f(x,y)}{\partial x} dy + f(x,b(x)) \frac{db}{dx} - f(x,a(x)) \frac{da}{dx} \notag \end{align} $$

2.3. 連続方程式の水深積分

図1 — 変数の定義

変数の定義は上図のようにする。ここで、$z_s$は自由表面位置、$z_b$は底面位置である。水深は$h=z_s-z_b$と定義される。$i,k$はそれぞれx,z方向の単位ベクトルである。

非圧縮性流体の連続方程式を水深方向に積分する。

$$ \begin{align} \int_{z_b}^{z_s} \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} + \frac{\partial w}{\partial z} dz &= 0 \notag \end{align} $$

各項は次式のように変形できる。ただし、ここではライプニッツの積分則を利用した。

$$ \begin{align} \int_{z_b}^{z_s} \frac{\partial u}{\partial x} dz &= \frac{\partial }{\partial x} \int_{z_b}^{z_s} u dz - u(z_s) \frac{\partial z_s}{\partial x} + u(z_b) \frac{\partial z_b}{\partial x} \notag \\ \int_{z_b}^{z_s} \frac{\partial v}{\partial y} dz &= \frac{\partial }{\partial y} \int_{z_b}^{z_s} v dz - v(z_s) \frac{\partial z_s}{\partial y} + v(z_b) \frac{\partial z_b}{\partial y} \notag \\ \int_{z_b}^{z_s} \frac{\partial w}{\partial z} dz &= w(z_s) - w(z_b) \notag \end{align} $$

結局これらをまとめると次式のようになる。

$$ \begin{align} \int_{z_b}^{z_s} \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} + \frac{\partial w}{\partial z} dz &= \frac{\partial }{\partial x} \int_{z_b}^{z_s} u dz +\frac{\partial }{\partial y} \int_{z_b}^{z_s} v dz \notag \\ &+w(z_s) -u(z_s) \frac{\partial z_s}{\partial x} - v(z_s) \frac{\partial z_s}{\partial y} \notag \\ &-w(z_b) + u(z_b) \frac{\partial z_b}{\partial x} + v(z_b) \frac{\partial z_b}{\partial y} \notag \end{align} $$

次に境界において成立する条件について考える。自由表面上の流体粒子は自由表面に沿って移動することになる。従って次式が成り立つ。

$$ \begin{align} \left( \mathbf{u}(z_s)-\frac{\partial z_s}{\partial t} \mathbf{k} \right) \cdot \nabla(z-z_s(x,y)) = 0 \notag \end{align} $$

成分ごとに展開すると次式のようになる。

$$ \begin{align} w(z_s) = \frac{\partial z_s}{\partial t} + u(z_s) \frac{\partial z_s}{\partial x} + v(z_s) \frac{\partial z_s}{\partial y} \notag \end{align} $$

同様に水底面上の流体粒子は水底面に沿って移動する。従って次式が成り立つ。

$$ \begin{align} \left( \mathbf{u}(z_b)-\frac{\partial z_b}{\partial t} \mathbf{k} \right) \cdot \nabla(z-z_b(x,y)) = 0 \notag \end{align} $$

$$ \begin{align} w(z_b) = \frac{\partial z_b}{\partial t} + u(z_b) \frac{\partial z_b}{\partial x} + v(z_b) \frac{\partial z_b}{\partial y} \notag \end{align} $$

ここで断面平均流速$U,V$を次式で定義する。

$$ \begin{align} U = \frac{1}{h} \int_{z_b}^{z_s} u dz, \quad V = \frac{1}{h} \int_{z_b}^{z_s} v dz \notag \end{align} $$

これらの関係式を用いると、最終的に連続方程式の水深積分は次式のように変形できる。

$$ \begin{align} \frac{\partial h}{\partial t} + \frac{\partial (Uh)}{\partial x} + \frac{\partial (Vh)}{\partial y} = 0 \notag \end{align} $$

2.4. 運動量方程式の水深積分

z方向の運動量方程式において、鉛直方向の速度、加速度成分が水平方向成分に比べて十分小さいと仮定すると、次式が成り立つ。

$$ \begin{align} 0 &= -\frac{\partial }{\partial z} \left(\frac{p}{\rho}+gz\right) \notag \end{align} $$

$z=z_s$で$p=0$とすると、圧力を静水圧近似で表すことができる。

$$ \begin{align} p &= \rho g (z_s - z) \notag \end{align} $$

x方向の運動量方程式の左辺項を水深方向に積分する。

$$ \begin{align} \int_{z_b}^{z_s} \frac{\partial u}{\partial t} dz &= \frac{\partial }{\partial t} \int_{z_b}^{z_s} u dz -u(z_s) \frac{\partial z_s}{\partial t} + u(z_b) \frac{\partial z_b}{\partial t} \notag \\ \int_{z_b}^{z_s} \frac{\partial uu}{\partial x} dz &= \frac{\partial }{\partial x} \int_{z_b}^{z_s} uu dz -u(z_s)u(z_s) \frac{\partial z_s}{\partial x} + u(z_b)u(z_b) \frac{\partial z_b}{\partial x} \notag \\ \int_{z_b}^{z_s} \frac{\partial uv}{\partial y} dz &= \frac{\partial }{\partial y} \int_{z_b}^{z_s} uv dz -u(z_s)v(z_s) \frac{\partial z_s}{\partial y} + u(z_b)v(z_b) \frac{\partial z_b}{\partial y} \notag \\ \int_{z_b}^{z_s} \frac{\partial uw}{\partial z} dz &= u(z_s)w(z_s) - u(z_b)w(z_b) \notag \end{align} $$

これらの項をまとめると次式のようになる。

$$ \begin{align} \int_{z_b}^{z_s} \left( \frac{\partial u}{\partial t} + \frac{\partial u u_j}{\partial x_j} \right) dz &= \frac{\partial }{\partial t} \int_{z_b}^{z_s} u dz +\frac{\partial }{\partial x} \int_{z_b}^{z_s} uu dz +\frac{\partial }{\partial y} \int_{z_b}^{z_s} uv dz \notag \\ &-u(z_s) \left( \frac{\partial z_s}{\partial t} +u(z_s) \frac{\partial z_s}{\partial x} +v(z_s) \frac{\partial z_s}{\partial y} -w(z_s) \right) \notag \\ &+u(z_b) \left( \frac{\partial z_b}{\partial t} +u(z_b) \frac{\partial z_b}{\partial x} +v(z_b) \frac{\partial z_b}{\partial y} -w(z_b) \right) \notag \end{align} $$

ここで流速の二次項を断面平均流速$U,V$を用いて次式のように定義する。

$$ \begin{align} \beta U U = \frac{1}{h} \int_{z_b}^{z_s} uu dz, \quad \beta U V = \frac{1}{h} \int_{z_b}^{z_s} uv dz \notag \end{align} $$

$\beta$は運動量補正係数である。自由表面および底面の境界条件の関係式を用いると、最終的にx方向の運動量方程式の左辺項の水深積分は次式のように表される。（ここでは$\beta=1$と近似）

$$ \begin{align} \int_{z_b}^{z_s} \left( \frac{\partial u}{\partial t} + \frac{\partial u u_j}{\partial x_j} \right) dz &= \frac{\partial Uh}{\partial t} + \frac{\partial U U h}{\partial x} + \frac{\partial U V h}{\partial y} \notag \end{align} $$

x方向の運動量方程式の右辺の圧力項を水深方向に積分すると次式のようになる。

$$ \begin{align} -g \int_{z_b}^{z_s} \frac{\partial z_s}{\partial x} dz &= -g \left( \frac{\partial }{\partial x} \int_{z_b}^{z_s} z_s dz -z_s \frac{\partial z_s}{\partial x} + z_s \frac{\partial z_b}{\partial x} \right) \notag \\ &= -g \left( \frac{\partial }{\partial x} z_s(z_s-z_b) -z_s \frac{\partial z_s}{\partial x} + z_s \frac{\partial z_b}{\partial x} \right) \notag \\ &= -g (z_s-z_b) \frac{\partial z_s }{\partial x} = -gh \frac{\partial z_s}{\partial x} \notag \end{align} $$

x方向の運動量方程式の右辺の粘性応力項を水深方向に積分する。

$$ \begin{align} \int_{z_b}^{z_s} \frac{\partial }{\partial x}\left( \frac{\tau_{xx}}{\rho}\right) dz &= \frac{\partial }{\partial x} \int_{z_b}^{z_s} \frac{\tau_{xx}}{\rho} dz - \frac{\tau_{xx}(z_s)}{\rho} \frac{\partial z_s}{\partial x} + \frac{\tau_{xx}(z_b)}{\rho} \frac{\partial z_b}{\partial x} \notag \\ \int_{z_b}^{z_s} \frac{\partial }{\partial y}\left( \frac{\tau_{xy}}{\rho}\right) dz &= \frac{\partial }{\partial y} \int_{z_b}^{z_s} \frac{\tau_{xy}}{\rho} dz - \frac{\tau_{xy}(z_s)}{\rho} \frac{\partial z_s}{\partial y} + \frac{\tau_{xy}(z_b)}{\rho} \frac{\partial z_b}{\partial y} \notag \\ \int_{z_b}^{z_s} \frac{\partial }{\partial z}\left( \frac{\tau_{xz}}{\rho}\right) dz &= \frac{\tau_{xz}(z_s)}{\rho} - \frac{\tau_{xz}(z_b)}{\rho} \notag \end{align} $$

これらの項をまとめると次式のようになる。

$$ \begin{align} \int_{z_b}^{z_s} \frac{\partial }{\partial x_j} \left( \frac{\tau_{xj}}{\rho} \right) dz &= \frac{\partial }{\partial x} \int_{z_b}^{z_s} \frac{\tau_{xx}}{\rho} dz +\frac{\partial }{\partial y} \int_{z_b}^{z_s} \frac{\tau_{xy}}{\rho} dz \notag \\ &-\frac{1}{\rho} \left( \tau_{xx}(z_s) \frac{\partial z_s}{\partial x} + \tau_{xy}(z_s) \frac{\partial z_s}{\partial y} - \tau_{xz}(z_s) \right) \notag \\ &+\frac{1}{\rho} \left( \tau_{xx}(z_b) \frac{\partial z_b}{\partial x} + \tau_{xy}(z_b) \frac{\partial z_b}{\partial y} - \tau_{xz}(z_b) \right) \notag \end{align} $$

ここで、自由表面形状は $f(x,y)=z-z_s(x,y)=0$と表されることから、自由表面上の法線ベクトルは次式で表される。

$$ \begin{align} \mathbf{n} = \frac{\nabla f}{|\nabla f|} &=(-\frac{\partial z_s}{\partial x}, -\frac{\partial z_s}{\partial y}, 1)/\sqrt{1+\left(\frac{\partial z_s}{\partial x}\right)^2+\left(\frac{\partial z_s}{\partial y}\right)^2} \notag \\ &\simeq(-\frac{\partial z_s}{\partial x}, -\frac{\partial z_s}{\partial y}, 1) \notag \end{align} $$

応力テンソルを$T$とすると法線ベクトル$\mathbf{n}$の作用面に働く応力$\mathbf{t}$は次式で表される。

$$ \begin{align} \mathbf{t} = \begin{pmatrix} t_x \\ t_y \\ t_z \end{pmatrix} =\mathbf{T} \cdot \mathbf{n} = \begin{pmatrix} \tau_{xx} & \tau_{xy} & \tau_{xz} \\ \tau_{yx} & \tau_{yy} & \tau_{yz} \\ \tau_{zx} & \tau_{zy} & \tau_{zz} \end{pmatrix} \begin{pmatrix} -\frac{\partial z_s}{\partial x} \\ -\frac{\partial z_s}{\partial y} \\ {1} \end{pmatrix} = \begin{pmatrix} -\tau_{xx} \frac{\partial z_s}{\partial x} - \tau_{xy} \frac{\partial z_s}{\partial y} + \tau_{xz} \\ -\tau_{yx} \frac{\partial z_s}{\partial x} - \tau_{yy} \frac{\partial z_s}{\partial y} + \tau_{yz} \\ -\tau_{zx} \frac{\partial z_s}{\partial x} - \tau_{zy} \frac{\partial z_s}{\partial y} + \tau_{zz} \end{pmatrix} \notag \end{align} $$

応力の断面平均値を次式で定義する。

$$ \begin{align} \bar{\tau}_{xx} = \frac{1}{h} \int_{z_b}^{z_s} \tau_{xx} dz, \quad \bar{\tau}_{xy} = \frac{1}{h} \int_{z_b}^{z_s} \tau_{xy} dz \notag \end{align} $$

これらの関係式を用いると、最終的にx方向の運動量方程式の水深積分は次式のように表される。

$$ \begin{align} &\frac{\partial Uh}{\partial t} + \frac{\partial U U h}{\partial x} + \frac{\partial U V h}{\partial y} = \notag \\ &-gh \frac{\partial z_s}{\partial x} + \frac{\partial }{\partial x} \left( \frac{h \bar{\tau}_{xx}}{\rho} \right) + \frac{\partial }{\partial y} \left( \frac{h \bar{\tau}_{xy}}{\rho} \right) + \frac{t_x(z_s)}{\rho} - \frac{t_x(z_b)}{\rho} \notag \end{align} $$

y方向についても同様に運動量方程式を水深積分した方程式を導くことができる。

$$ \begin{align} &\frac{\partial Vh}{\partial t} + \frac{\partial U V h}{\partial x} + \frac{\partial V V h}{\partial y} = \notag \\ &-gh \frac{\partial z_s}{\partial y} + \frac{\partial }{\partial x} \left( \frac{h \bar{\tau}_{yx}}{\rho} \right) + \frac{\partial }{\partial y} \left( \frac{h \bar{\tau}_{yy}}{\rho} \right) + \frac{t_y(z_s)}{\rho} - \frac{t_y(z_b)}{\rho} \notag \end{align} $$

以上で導いた連続方程式の水深積分式、およびx方向・y方向の運動量方程式の水深積分式が、浅水流を計算する方程式である。

3. OpenFOAMにおける浅水方程式

OpenFOAMでは浅水方程式を解くために標準ソルバー shallowWaterFoamが利用可能である。上記の導出では示していないが、流体にコリオリ力が働く場合も考慮可能である。

ただし、粘性応力に起因する項$\bar{\tau}_{xx},\bar{\tau}_{xy}$、および境界作用力$t(z_s),t(z_b)$などは組み込まれていない。そのため底面摩擦力などを考慮したい場合にはソルバーをカスタマイズする必要がある。

4. 解析の安定性

浅水方程式は双曲型の微分方程式であるため、計算を精度よく安定して行うためには時間刻み$\Delta t$を許容値以下に設定する必要がある。

簡単のため、$z_b=0$とし、流れがx方向のみの一次元場を考える。この場合、解くべき連立微分方程式は以下のようになる。

$$ \begin{align} \frac{\partial}{\partial t} \begin{bmatrix} h \\ U \\ V \end{bmatrix} + \begin{bmatrix} U & h & 0 \\ g & U & 0 \\ 0 & 0 & U \end{bmatrix} \frac{\partial}{\partial x} \begin{bmatrix} h \\ U \\ V \end{bmatrix} =0 \notag \end{align} $$

係数行列の固有値を計算すると、$U,U\pm\sqrt{gh}$が得られる。 $\Delta t$をメッシュの刻み幅とこれらの特性流速値の比から求められる時間刻み幅よりも小さくしておく必要がある。