一维浅水方程

· · 个人记录

顾名思义，浅水方程只有在水很“浅”的时候成立，这个假设意味着我们假设水在垂直方向没有速度。

那么我们就要求出两个函数：h(x,t) 和 v(x,t)，分别表示 t 时间时 x 处水位和水平速度。

众所周知，在搞这种物理模型的时候我们要找守恒量。

物理中有能量守恒、动量守恒、角动量守恒，看上去能量和角动量没什么用，所以只有动量守恒再加上水的质量守恒。

先假设水的密度为 \rho=1，重力加速度为 g，以右为正方向。

考虑分析 [x_0,x_1] 内水量变化，也就是说 \dfrac{\text d}{\text dt}\int_{x_0}^{x_1}h\text dx。

显然，我们需要考察在 x_0,x_1 处流入或流出的水的总量。

先看 x_0，水以 v(x_0,t) 的速度流入，且高度为 h(x_0,t)，所以 x_0 处水流入的速度为 h(x_0,t)v(x_0,t)。

再看 x_1，水以 v(x_1,t) 的速度流出，且高度为 h(x_1,t)，所以 x_1 处水流出的速度为 h(x_1,t)v(x_1,t)。

所以，我们可以得到：

\frac{\text d}{\text dt}\int_{x_0}^{x_1}h(x,t)\text dx=h(x_0,t)v(x_0,t)-h(x_1,t)v(x_1,t)

这个形式被称作积分形式。

常识告诉我们，一般来说微分形式更简单。

\frac{\text d}{\text dt}\int_{x_0}^{x_1}h(x,t)\text dx=\int_{x_0}^{x_1}\frac\partial{\partial t}h\text dx h(x_0,t)v(x_0,t)-h(x_1,t)v(x_1,t)=\int_{x_0}^{x1}-\frac\partial{\partial x}h(x,t)v(x,t)\text dx

所以：

\int_{x_0}^{x_1}\frac\partial{\partial t}h(x,t)+\frac\partial{\partial x}h(x,t)v(x,t)\text dx=0 \frac\partial{\partial t}h(x,t)+\frac\partial{\partial x}h(x,t)v(x,t)=0

再次分析 [x_0,x_1] 区间内动量的变化，即 \dfrac{\text d}{\text dt}\int_{x_0}^{x_1}h(x,t)v(x,t)\text dt。

众所周知，这次在 x_0 处流入的水有两个动量来源：一是流入的水的动量，二是边界外水的压力。

先分析 x_0：

显然，来源一所带来的动量为 h(x_0,t)v^2(x_0,t)。

考虑求出压力的总和，众所周知压力是压强的积分，且压强的表达式为 \rho gh，在我们的设定中，\rho=1，h=h(x_0,t)-y，所以 x_0 左侧的水对于 x_0 内水体的压力为：

\int_0^{h(x_0,t)}\rho gh\text dy=g\int_0^{h(x_0,t)}[h(x_0,t)-y]\text dy=\frac{gh^2(x_0,t)}2

所以 x_0 处对于区间内水体动量的贡献为：h(x_0,t)v^2(x_0,t)+\dfrac{gh^2(x_0,t)}2。

所以我们可以得到：

\frac{\text d}{\text dt}\int_{x_0}^{x_1}h(x,t)v(x,t)\text dx=h(x_0,t)v^2(x_0,t)+\frac{gh^2(x_0,t)}2-h(x_1,t)v^2(x_1,t)-\frac{gh^2(x_1,t)}2

这就是积分形式，通过一通化简我们可以得到微分形式：

\frac\partial{\partial t}h(x,t)v(x,t)+\frac\partial{\partial t}\left(h(x,t)v^2(x,t)+\frac{gh^2(x,t)}2\right)=0

这样，我们就得到了两个方程：

微分形式：

\begin{cases}\dfrac\partial{\partial t}h(x,t)+\dfrac\partial{\partial x}h(x,t)v(x,t)=0\\\dfrac\partial{\partial t}h(x,t)v(x,t)+\dfrac\partial{\partial t}\left(h(x,t)v^2(x,t)+\dfrac{gh^2(x,t)}2\right)=0\end{cases}

积分形式：

\begin{cases}\dfrac{\text d}{\text dt}\int_{x_0}^{x_1}h(x,t)\text dx=h(x_0,t)v(x_0,t)-h(x_1,t)v(x_1,t)\\\dfrac{\text d}{\text dt}\int_{x_0}^{x_1}h(x,t)v(x,t)\text dx=h(x_0,t)v^2(x_0,t)+\dfrac{gh^2(x_0,t)}2-h(x_1,t)v^2(x_1,t)-\dfrac{gh^2(x_1,t)}2\end{cases}

然后就是解方程组了，这里我选择有限体积法进行数值解。

首先，假设我们要在 [0,1] 区间内求解浅水方程，我们先把这个区间分成 n 段：[x_0,x_1],[x_1,x_2],[x_2,x_3],\dots,[x_{n-1},x_n]。显然 x_0=0,x_n=1，这些区间的长度可以不一样，不过为简便起见，设其均为 \Delta x。即对于所有的 1\le i\le n，都有 x_i-x_{i-1}=\Delta x。

然后，对于每个小区间，我们存储守恒量（在 1D SWE 中就是质量和动量）的平均值，也就是：

M_i(t)=\frac1{\Delta x}\int_{x_i}^{x_{i+1}}h(x,t)\text dx P_i(t)=\frac1{\Delta x}\int_{x_i}^{x_{i+1}}h(x,t)v(x,t)\text dx

代入积分形式的方程组中，可得：

\begin{cases}\dfrac{\text d}{\text dt}\Delta xM_i(t)=h(x_i,t)v(x_i,t)-h(x_{i+1},t)v(x_{i+1},t)\\\dfrac{\text d}{\text dt}\Delta xP_i(t)=h(x_i,t)v^2(x_i,t)+\dfrac{gh^2(x_i,t)}2-h(x_{i+1},t)v^2(x_{i+1},t)-\dfrac{gh^2(x_{i+1},t)}2\end{cases}

但是有一个小问题：我们不知道 h(x,t) 的值。

显然，对于 x_i<x<x_{i+1} 的所有 x，我们都可以近似认为 h(x,t)=M_i(t),v(x,t)=\dfrac{P_i(t)}{M_i(t)}。

但是在边界处怎么办呢？

这里我使用 Lax-Friedrich 通量，即：

h(x_i,t)v(x_i,t)=\frac{P_{i-1}(t)+P_i(t)}2+c\frac{M_{i-1}(t)-M_i(t)}2

其中 c 为波速，这里我们直接取常数即可。

还有就是 M_{-1},P_{-1},P_n,M_n 怎么办。

有两种边界条件：

  1. 假设这是一个水箱，也就是说边界可以“反射”水，那么 M_{-1}=M_{0},P_{-1}=-P_0,M_n=M_{n-1},P_n=-P_{n-1}。

  2. 假设边界是循环的，也就是说 0 和 1 实际上是拼在一起的，那么 M_{-1}=M_{n-1},P_{-1}=P_{n-1},M_n=M_0,M_n=P_0。

现在，我们就有：

\begin{cases}\dfrac{\text d}{\text dt}M_i(t)=\dfrac{h(x_i,t)v(x_i,t)-h(x_{i+1},t)v(x_{i+1},t)}{\Delta x}\\\dfrac{\text d}{\text dt}P_i(t)=\dfrac{h(x_i,t)v^2(x_i,t)-h(x_{i+1},t)v^2(x_{i+1},t)}{\Delta x}+\dfrac{gh^2(x_i,t)-gh^2(x_{i+1},t)}{2\Delta x}\end{cases}

其中：

h(x_i,t)v(x_i,t)=\frac{P_{i-1}(t)+P_i(t)}2+c\frac{M_{i-1}(t)-M_i(t)}2 h(x_i,t)v^2(x_i,t)+\frac{gh^2(x_i,t)}2=\frac{\dfrac{P_{i-1}^2(t)}{M_{i-1}(t)}+\dfrac{P_i^2(t)}{M_i(t)}}2+\frac{gM_{i-1}^2(t)+gM_i^2(t)}4+c\frac{P_{i-1}(t)-P_i(t)}2

然后我们就要对其进行数值积分。

数值积分的要求是：给定一个带有初值的微分方程：

y'(t)=f(t,y(t)),y(t_0)=y_0

求数值解。

一个显然的想法是：

y(t+\Delta t)\approx y(t)+\Delta tf(t,y(t))

这个方法被称作欧拉法。

但是更常用的是四阶 Runge-Kutta 法，又叫 RK4，其公式为：

y(t+\Delta t)=y(t)+\frac{\Delta t}6(k_1+2k_2+2k_3+k_4)

其中：

\begin{cases}k_1=f(y(t),t)\\k_2=f\left(y(t)+\Delta t\dfrac{k_1}2,t+\dfrac{\Delta t}2\right)\\k_3=f\left(y(t)+\Delta t\dfrac{k_2}2,t+\dfrac{\Delta t}2\right)\\k_4=f(y(t)+\Delta tk_3,t+\Delta t)\end{cases}

下面是欧拉法的代码：

import turtle
from math import sin, cos, pi
n = 100
dx = 1 / n
g = 1
c = 0.2
M = [0.0 for i in range(n)]
P = [0.0 for i in range(n)]
for i in range(n):
    if i < n // 2:
        M[i] = 2
    else:
        M[i] = 0.1
'''
for i in range(n):
    M[i] = sin(i / n * 1 * pi) / 5 + 0.3
'''
'''
for i in range(n):
    M[i] = 0.1 + i / n
'''
'''
for i in range(n):
    if i < n // 2:
        M[i] = sin(i / n * 2 * pi) / 5 + 0.3
    else:
        M[i] = 0.3
'''
turtle.tracer(False)
def draw():
    turtle.clear()
    turtle.penup()
    turtle.goto(0, 0)
    turtle.pendown()
    turtle.pencolor((0, 0, 1))
    for i in range(n):
        turtle.goto(i / n * 300, M[i] * 100)
    turtle.update()
dt = 0.0001
def next_time():
    old_M = []
    for i in range(n):
        old_M.append(M[i])
    old_P = []
    for i in range(n):
        old_P.append(P[i])
    for i in range(n):
        if i - 1 == -1:
            hvi = (-old_P[0] + old_P[i]) / 2 + c * (old_M[0] - old_M[i]) / 2
            #hvi = (old_P[n - 1] + old_P[i]) / 2 + c * (old_M[n - 1] - old_M[i]) / 2
        else:
            hvi = (old_P[i - 1] + old_P[i]) / 2 + c * (old_M[i - 1] - old_M[i]) / 2
        if i + 1 == n:
            hvi1 = (old_P[i] - old_P[n - 1]) / 2 + c * (old_M[i] - old_M[n - 1]) / 2
            #hvi1 = (old_P[i] + old_P[0]) / 2 + c * (old_M[i] - old_M[0]) / 2
        else:
            hvi1 = (old_P[i] + old_P[i + 1]) / 2 + c * (old_M[i] - old_M[i + 1]) / 2
        M[i] += (hvi - hvi1) / dx * dt
        if i - 1 == -1:
            hv2i = (old_P[0] * old_P[0] / old_M[0] + old_P[i] * old_P[i] / old_M[i]) / 2 + g * (old_M[0] * old_M[0] + old_M[i] * old_M[i]) / 4 + c * (-old_P[0] - old_P[i]) / 2
            #hv2i = (old_P[n - 1] * old_P[n - 1] / old_M[n - 1] + old_P[i] * old_P[i] / old_M[i]) / 2 + g * (old_M[n - 1] * old_M[n - 1] + old_M[i] * old_M[i]) / 4 + c * (old_P[n - 1] - old_P[i]) / 2
        else:
            hv2i = (old_P[i - 1] * old_P[i - 1] / old_M[i - 1] + old_P[i] * old_P[i] / old_M[i]) / 2 + g * (old_M[i - 1] * old_M[i - 1] + old_M[i] * old_M[i]) / 4 + c * (old_P[i - 1] - old_P[i]) / 2
        if i + 1 == n:
            hv2i1 = (old_P[i] * old_P[i] / old_M[i] + old_P[n - 1] * old_P[n - 1] / old_M[n - 1]) / 2 + g * (old_M[i] * old_M[i] + old_M[n - 1] * old_M[n - 1]) / 4 + c * (old_P[i] + old_P[n - 1]) / 2
            #hv2i1 = (old_P[i] * old_P[i] / old_M[i] + old_P[0] * old_P[0] / old_M[0]) / 2 + g * (old_M[i] * old_M[i] + old_M[0] * old_M[0]) / 4 + c * (old_P[i] + old_P[0]) / 2
        else:
            hv2i1 = (old_P[i] * old_P[i] / old_M[i] + old_P[i + 1] * old_P[i + 1] / old_M[i + 1]) / 2 + g * (old_M[i] * old_M[i] + old_M[i + 1] * old_M[i + 1]) / 4 + c * (old_P[i] - old_P[i + 1]) / 2
        P[i] += (hv2i - hv2i1) / dx * dt
while True:
    for i in range(100):
        next_time()
    draw()

RK4 的代码：

import turtle
from time import sleep
from math import sin, cos, pi, exp
n = 100
dx = 1 / n
g = 1
c = 0.2
y = [[0.0 for i in range(n)], [0.0 for i in range(n)]]
# y[0] = M
# y[1] = P
for i in range(n):
    y[0][i] = exp(-64 * (i / n - 0.25) * (i / n - 0.25)) * 0.1 + 0.1
'''
for i in range(n):
    y[0][i] = 0.1
for i in range(40, 60):
    y[0][i] = 1
'''
'''
for i in range(n):
    if i < n // 2:
        y[0][i] = 2
    else:
        y[0][i] = 0.1
'''
'''
for i in range(n):
    y[0][i] = sin(i / n * 2 * pi) + 1.1
'''
'''
for i in range(n):
    y[0][i] = 0.1 + i / n
'''
'''
for i in range(n):
    if i < n // 2:
        y[0][i] = sin(i / n * 2 * pi) / 5 + 0.3
    else:
        y[0][i] = 0.3
'''
turtle.tracer(False)
def draw():
    turtle.clear()
    turtle.penup()
    turtle.goto(0, 0)
    turtle.pendown()
    turtle.pencolor((0, 0, 1))
    for i in range(n):
        turtle.goto(i / n * 300, y[0][i] * 250)
    turtle.update()
dt = 0.01
def dydt():
    M = y[0]
    P = y[1]
    dydt = [[0.0 for i in range(n)], [0.0 for i in range(n)]]
    for i in range(n):
        if i - 1 == -1:
            hvi = (-P[0] + P[i]) / 2 + c * (M[0] - M[i]) / 2
            #hvi = (P[n - 1] + P[i]) / 2 + c * (M[n - 1] - M[i]) / 2
        else:
            hvi = (P[i - 1] + P[i]) / 2 + c * (M[i - 1] - M[i]) / 2
        if i + 1 == n:
            hvi1 = (P[i] - P[n - 1]) / 2 + c * (M[i] - M[n - 1]) / 2
            #hvi1 = (P[i] + P[0]) / 2 + c * (M[i] - M[0]) / 2
        else:
            hvi1 = (P[i] + P[i + 1]) / 2 + c * (M[i] - M[i + 1]) / 2
        dydt[0][i] = (hvi - hvi1) / dx
        if i - 1 == -1:
            hv2i = (P[0] * P[0] / M[0] + P[i] * P[i] / M[i]) / 2 + g * (M[0] * M[0] + M[i] * M[i]) / 4 + c * (-P[0] - P[i]) / 2
            #hv2i = (P[n - 1] * P[n - 1] / M[n - 1] + P[i] * P[i] / M[i]) / 2 + g * (M[n - 1] * M[n - 1] + M[i] * M[i]) / 4 + c * (P[n - 1] - P[i]) / 2
        else:
            hv2i = (P[i - 1] * P[i - 1] / M[i - 1] + P[i] * P[i] / M[i]) / 2 + g * (M[i - 1] * M[i - 1] + M[i] * M[i]) / 4 + c * (P[i - 1] - P[i]) / 2
        if i + 1 == n:
            hv2i1 = (P[i] * P[i] / M[i] + P[n - 1] * P[n - 1] / M[n - 1]) / 2 + g * (M[i] * M[i] + M[n - 1] * M[n - 1]) / 4 + c * (P[i] + P[n - 1]) / 2
            #hv2i1 = (P[i] * P[i] / M[i] + P[0] * P[0] / M[0]) / 2 + g * (M[i] * M[i] + M[0] * M[0]) / 4 + c * (P[i] + P[0]) / 2
        else:
            hv2i1 = (P[i] * P[i] / M[i] + P[i + 1] * P[i + 1] / M[i + 1]) / 2 + g * (M[i] * M[i] + M[i + 1] * M[i + 1]) / 4 + c * (P[i] - P[i + 1]) / 2
        dydt[1][i] = (hv2i - hv2i1) / dx
    return dydt
def next_time():
    k1 = dydt()
    for i in range(n):
        y[0][i] += k1[0][i] * dt / 2
        y[1][i] += k1[1][i] * dt / 2
    k2 = dydt()
    for i in range(n):
        y[0][i] -= k1[0][i] * dt / 2
        y[1][i] -= k1[1][i] * dt / 2
        y[0][i] += k2[0][i] * dt / 2
        y[1][i] += k2[1][i] * dt / 2
    k3 = dydt()
    for i in range(n):
        y[0][i] -= k2[0][i] * dt / 2
        y[1][i] -= k2[1][i] * dt / 2
        y[0][i] += k3[0][i] * dt
        y[1][i] += k3[1][i] * dt
    k4 = dydt()
    for i in range(n):
        y[0][i] -= k3[0][i] * dt
        y[1][i] -= k3[1][i] * dt
    for i in range(n):
        y[0][i] += dt * (k1[0][i] + 2 * k2[0][i] + 2 * k3[0][i] + k4[0][i]) / 6
        if abs(y[0][i]) < 0.01:
            y[0][i] = -0.01
        y[1][i] += dt * (k1[1][i] + 2 * k2[1][i] + 2 * k3[1][i] + k4[1][i]) / 6
while True:
    for i in range(1):
        next_time()
    draw()
    sleep(0.02)