1次元移流分散（有限差分法、Python）

GitHubで水理公式集例題集2024 のコードを眺めていた時に、地下水の不飽和浸透や移流分散計算がUPされているのを見かけました。

EXCEL版 と Fortran版 があり、懐かしいなあと思いつつ「例題 1.24 水平一次元帯水層を対象とした溶質輸送過程の計算（数値解）」の Fortran コードを Python に移植してみました。

で、結果が全く同じだったことには少し驚かされました。

以下、コードに沿ったFTCS（Forward-Time Centered-Space）による離散化とスクリプトです。手元に例題集が無いため、表現は異なるかもしれません。スクリプトではスライスを利用することで計算が1行となり、見やすくなりました。実際の地下では流速や流量が一定ではないので、ここまで簡素にはできませんが、まずはここからなのでしょう。

1. 支配方程式

∂c/∂t + v·∂c/∂x = D·∂²c/∂x² ……(1)

ここで、

c : 濃度 (mg/L)

t : 時間 (day)

x : 距離 (m)

v : 移流速度 (m/day)

D : 分散係数 (m²/day)

2. 差分近似

2.1 時間微分の前進差分

時刻 t における時間微分を前進差分で近似すると、

∂c/∂t ≈ (c_i^(t+1) - c_i^t) / Δt ……(2)

ここで、

c_i^t : 格子点 x_i、時刻 t における濃度

Δt : 時間ステップ

2.2 空間2次微分の中心差分（拡散項）

格子点 x_i における空間2次微分を中心差分で近似すると、

∂²c/∂x² ≈ (c_(i+1)^t - 2·c_i^t + c_(i-1)^t) / Δx² ……(3)

ここで、

Δx : 空間格子間隔

2.3 空間1次微分の後退差分（移流項）

v > 0（下流方向の流れ）を仮定し、上流差分（upwind差分）を適用すると、

∂c/∂x ≈ (c_i^t - c_(i-1)^t) / Δx ……(4)

3. 差分方程式の導出

式(1)に式(2)、(3)、(4)を代入すると、

(c_i^(t+1) - c_i^t) / Δt + v·(c_i^t - c_(i-1)^t) / Δx

= D·(c_(i+1)^t - 2·c_i^t + c_(i-1)^t) / Δx² ……(5)

式(5)の両辺に Δt を掛け、c_i^(t+1) について解くと、

c_i^(t+1) = c_i^t

+ (D·Δt/Δx²)·(c_(i+1)^t - 2·c_i^t + c_(i-1)^t)

- (v·Δt/Δx)·(c_i^t - c_(i-1)^t) ……(6)

ここで、

r_diff = D·Δt / Δx² （拡散数）

r_adv = v·Δt / Δx （クーラン数）

を用いると、式(6)は以下となる。

c_i^(t+1) = c_i^t

+ r_diff·(c_(i+1)^t - 2·c_i^t + c_(i-1)^t)

- r_adv·(c_i^t - c_(i-1)^t) ……(7)

import numpy as np
import csv

# +++++++++++++++++++++++++++++++++++++++++++++++
# 1次元移流分散方程式の数値解析
# 有限差分法：FTCS（Forward-Time Centered-Space）
# ∂c/∂t + v ∂c/∂x = D ∂²c/∂x²
# +++++++++++++++++++++++++++++++++++++++++++++++

def reidai082():
    # ----- パラメータ設定 -----
    xmax  = 30.0        # 解析領域の長さ (m)
    dx    = 1.0         # 差分格子間隔 (m)
    dx2 = dx ** 2     # 差分格子間隔の二乗
    imax = int(xmax / dx) + 1    # 下流境界の格子点番号

    xo  = 10.0   # 濃度観測地点（m）
    ii1 = int(xo / dx)   # Python 0-based インデックス

    vel   = 0.1         # 実流速 (m/day)
    al    = 1.0         # 縦分散長 (m)
    tau   = 1.0         # 屈曲率 (-)
    dm    = 1.0e-5      # 分子拡散係数 (cm^2/s)
    # 分散係数 D (m2/day)
    # dm を m2/day に換算するために 1 cm2/s = 8.64 m2/day
    D     = al * abs(vel) + tau * dm * 8.64

    bc    = 1.0         # 上流境界の濃度 (mg/L)

    tend  = 300.0       # 計算時間 (day)
    dt    = 1.0         # 時間ステップ (day)
    nmax  = int(tend / dt)   # 計算時間のステップ数

    tjikan = 1.0        # 結果をファイルに出力する時間間隔（day）
    iout1 = int(tjikan / dt)   # 結果をファイルに出力する時間間隔のステップ数

    # ----- グリッド設定 -----
    x    = np.linspace(0, xmax, imax)

    # ----- 配列の用意 -----
    c  = np.zeros(imax)  # 新しい時刻の濃度
    co = np.zeros(imax)  # １ステップ前の濃度

    # ----- 係数の用意 -----
    r_diff = D * dt / dx2
    r_adv  = vel * dt / dx

    # 初期条件：上流境界だけ bc
    co[0] = bc

    # ----- 出力ファイル準備 -----
    out_filename = 'output.csv'
    with open(out_filename, mode='w', newline='') as f:
        writer = csv.writer(f)
        writer.writerow(['elapsed-time(days)', 'concentration(mg/L)'])
        # t=0 のデータ
        writer.writerow([f"{0.0:.1f}", f"{co[ii1]:.7f}"])

        # ----- 時間ループ -----
        time = 0.0
        iout  = 0
        for n in range(1, nmax+1):
            time += dt
            iout += 1

            # ---- 境界条件 ----
            co[0]      = bc      # 上流
            co[-1]     = 0.0     # 下流

            # ---- 内部格子の更新 (FTCS) ----
            # c[i] = co[i] + D*dt/dx^2*(co[i+1]-2*co[i]+co[i-1])
            #                  - vel*dt/dx*(co[i]-co[i-1])
            # スライスで一括計算
            # co[1:-1]：インデックス 1～N−2 の内部点
            # co[2:] ：インデックス 2～N−1
            # co[:-2] ：インデックス 0～N−3

            c[1:-1] = (co[1:-1]
                       + r_diff * (co[2:] - 2*co[1:-1] + co[:-2])
                       - r_adv  * (co[1:-1] - co[:-2]))

            # 更新
            co[:] = c[:]

            # ---- 出力 ----
            if iout == iout1:
                print(f"{time:6.1f} days")
                writer.writerow([f"{time:.1f}", f"{co[ii1]:.7f}"])
                iout = 0