移流分散方程式を解いてみる

背景

ある日、先生からメールが届いた。
要約すると、宿題出すから解いてこい、とのことだった。マジで急にきてビビり散らかしました。
宿題が課されているのは文面から理解できたけど、具体的にどのような回答が期待されているか推し量りかねた（プログラミングが要りそうという見当はついた）。

コーディングするには、まず要件定義から。先生にところに足を運び、具体的に解くべき問題を明らかにした。この時点で回答期限が迫ってるんだよなぁ…

宿題内容

こんな問題
観測孔のある区間に，鉛直深さLでトレーサをT0時間注入し続け，その後はゼロ濃度を保つ．観測孔からX離れた下流側の観測孔にて得られるBTCsを深度ごとに考察せよ．

BTCsは、ブレークスルーカーブのこと。流体力学のアレです。

回答

まずは一次元から

1次元の移流分散方程式は、

$$
D\frac{∂^2C}{∂x^2} - v\frac{∂C}{∂x} = R\frac{∂C}{∂t} \tag{1}
$$

にて与えられる$^{1)}$。ここに、$D~$は拡散係数、$v~$は実流速、$C(x,t)~$は位置$~x~$と時間$~t~$を変数に持つ濃度、$~R~$は遅延係数である。今回は、遅延係数$~R~$は経験的に用いられている$~1~$を採用した。

移流分散方程式は、陽解法を用いて離散化して解いた$^{2)}$。離散化の手順については割愛。

まずは、深度方向を考えずに1次元において得られるブレークスルーカーブを取得する。
以下の条件を仮定する。

• 初期条件：$C(x=0,t=0)=C_0$
• 境界条件：$C(x=x,t=0) , x \neq 0 $

次に、グラフ（BTCs）を描く関数を作成。

graphを描く関数（1,2次元両方OK）

import numpy as np
import matplotlib.pyplot as plt

def graph(C,pos_x=None,pos_z=None,title=None):
  t = np.linspace(0,A-1,A)
  plt.plot(t,C)
  plt.title(title)
  plt.xlabel("Time, $t$")
  if pos_x is not None and pos_z is not None:
    plt.ylabel(f"Concentrate, $C$ | $x={pos_x}$,$z={pos_z}$")
  elif pos_x is not None:
    plt.ylabel(f"Concentrate, $C$ | $x={pos_x}$")
  elif pos_z is not None:
    plt.ylabel(f"Concentrate, $C$ | $z={pos_z}$")
  else:
    plt.ylabel(f"Concentrate, $C$")
  plt.show()

引数は濃度C、BTCsを得たい地点の座標、グラフタイトル。
なお、今後の2次元を想定して座標の引数には pos_z も投入可能にした。

続いて、一次元移流分散方程式を解く関数を定義。前述の通り、離散化して解いてる。

一次元移流分散方程式を計算

def BTC1(C, D=0.3,v=0.1):
  delta_x, delta_t = 1, 1
  C[0,:T0] = C0

  for t in range(0,A-1):
    for x in range(1,A-1):
      # t=tにおける，Cの位置xが可変的な配列．つまり，時間tにおける，濃度分布
      C[x,t+1] = delta_t*(D * (C[x+1,t]- 2*C[x,t] + C[x-1,t]) / delta_x**2 - v * (C[x+1,t]-C[x,t]) / delta_x) + C[x,t]

  return C 

後は、テキトーに（但し各変数のオーダーはある程度自然現象のそれと合わせた）決める。

一次元移流分散方程式を計算（続き）

# ----------------------  variable data ------------------------ #

A = 200 # size of t, x
C = np.zeros((A,A)) #C(x,t)
C[:,0] = 0 # 境界条件
T0 = 3 # 1なら、瞬間注入。Aならずっと定常注入。0ならトレーサーは流さない。
C0 = 10 # 初期濃度
X = 3 # 濃度変化をグラフに出力したい地点
# --------------------------------------------------------------- #

C = BTC1(C, D=0.2,v=0.01) # C:濃度行列, D:分散係数, v;実流速
graph(C[X,:],pos_x=X,title=f"BTCs of X = {X}, Tracer kept flowing to {T0}") # 可視化
plt.show()

てな感じで、描いたBTCsがコレ。

うんうん、急激な立ち上がりがあってロングテールがあるな。良いだろう（感覚派）。

お題である、二次元へ

続いて、二次元を考える。二次元移流分散方程式は、

$$
D_x\frac{∂^2C}{∂x^2} + D_z\frac{∂^2C}{∂z^2} - u\frac{∂C}{∂x} - v\frac{∂C}{∂z}= R\frac{∂C}{∂t} \tag{2}
$$

にて与えられる$^{3)}$。ここに、$D_x$は縦分散係数、$D_z$は横分散係数、$v$は実流速、$C(x,z,t)$は位置$~x,z~$と時間$~t~$を変数に持つ濃度、$R$は遅延係数であり、$1~$とした。

二次元移流分散方程式を解くプログラムでは、次の条件にて二次元におけるブレークスルーカーブを得た。

• トレーサー注入時間$~T_0 = 3$
• 初期濃度$~C_0 = 10$
• 濃度変化をグラフに出力したい地点$~X , Z = ( 3, 4 )$
• 縦分散係数$~D_x = 0.3$
• 横分散係数$~D_z = 0.03$
• 縦方向の速度$~v = 0.0001$
• 横方向の速度$~u = 0$

なお、グラフ化したときにBTCsが見やすいように値を設定しており、値そのものに意味は無い。ただし、縦分散係数と横分散係数と流速の比率はある程度室内実験をしてる友人から聞いた値をもとに決定した。

本質は一次元と一緒。graph を描く関数は同じものを使用してる。

二次元移流分散方程式を解く

def BTC2(C, Dx=0.3,Dz=0,v=0.1,u=0):
  C[0,0,:T0] = C0
  delta_x,delta_z, delta_t = 1, 1, 1

  for t in range(0,A-1):
    for z in range(A-1):
      for x in range(A-1):
        # t=tにおける，Cの位置xが可変的な配列．つまり，時間tにおける，濃度分布
        C[x,z,t+1] = delta_t*(
            Dx * (C[x+1,z,t]- 2*C[x,z,t] + C[x-1,z,t]) / delta_x**2 + Dz * (C[x,z+1,t]- 2*C[x,z,t] + C[x,z-1,t]) / delta_z**2  
            - v * (C[x+1,z,t]-C[x,z,t]) / delta_x - u * (C[x,z+1,t]-C[x,z,t]) / delta_z
            ) + C[x,z,t]
  return C 

# ----------------------  variable data ------------------------ #

A = 200 # size of t, x
C = np.zeros((A,A,A)) #C(x,z,t)
C[:,:,0] = 0 # 境界条件
T0 = 3 # 1なら、瞬間注入。Aならずっと定常注入。0ならトレーサーは流さない。
C0 = 10 # 初期濃度
X , Z = ( 3, 4 ) # 濃度変化をグラフに出力したい地点
Dx = 0.3 # 縦分散係数
Dz = 0.03 # 横分散係数
v  = 0.0001 # x方向の速度
u  = 0 # z方向の速度

# --------------------------------------------------------------- #

C = BTC2(C, Dx=Dx,Dz=Dz,v=v,u=u) # C:濃度行列, D:分散係数, v,u;実流速


# 可視化
graph(C[X,Z,:],pos_x=X,pos_z=Z,title=f"BTCs of (X,Z) = ({X},{Z}), Tracer kept flowing to {T0}")

得られたのがコレ。

考察っぽいことしようぜ

ここで、改めてお題を。

観測孔からX離れた下流側の観測孔にて得られるBTCsを深度ごとに考察せよ．

てことで、深度ごとに（zの値を変えて）得られるBTCsについて考える。
観測孔から距離 X=3 における、深度 Z=3,5,7 を表示した。

見やすいように、複数のグラフを表示する関数を以下に別途定義し、BTCsを描く。

BTCsを並べて描く

def multi_graph(pos_x=3,pos_z=[3,5,7]):
  t = np.linspace(0,A-1,A)
  fig, axes = plt.subplots(1, len(pos_z),sharey=True,figsize=((len(pos_z)+1)*4,len(pos_z)+1))
  axes[0].set_ylabel('Concentrate, $C$',labelpad=15)

  for i,z in enumerate(pos_z):
    axes[i].plot(t,C[X,z,:],label= f'BTCs at (X,Z)=({pos_x},{pos_z[i]})')
    axes[i].set_xlabel('Time, $t$',labelpad=15)
    axes[i].axvline(x=np.argmax(C[X,z,:]),color='cyan',linestyle='--',label="$t_{max}= $"+f"{np.argmax(C[X,z,:])}")
    axes[i].legend(loc="upper center",bbox_to_anchor=(0.5, 1.23))

  plt.show()


multi_graph(pos_x=X, pos_z=[3,5,7]) # 可視化

コレが描かれる。

グラフの青線は各地点におけるBTCsを、水色の点線はピーク濃度における時間$~t~$（グラフの極大値を通ってる線）を表している。

以下、回答にて示した考察の抜粋。来年同じ課題が課された後輩諸君！この記事を見つけてコピペすると既視感あるなってバレるから注意してね。…と言っても、考察というより事実の列挙だけど。

----- ここから -----
得られたBTCsより、深度が大きければ大きいほど、到達する濃度の量は少なくなることがわかる。上のグラフは縦分散係数$~D_x=0.3~$、横分散係数$~D_z~$を縦分散係数の約十分の一である$~0.03~$としている。横方向に比べて縦方向の分散の影響が大きく、また、横方向の実流速も$~0~$としているため、地点$~(0,0)~$にて流したトレーサーの到達量が小さくなるためと考えられる。

また、BTCsに着目すると、深度が大きくなるに連れてピーク濃度値を与える時間$~t~$は大きくなることが分かる。このことからも、深度が深くなるに連れて到着が遅れていることが推察される。
----- ここまで -----

終わりに

ここまでで既に回答期限を超えちゃってたけど←
…まあ、その分キチンと回答できたかな。

唐突な課題が今後も来るかもしれない。つかの間の休息。ま、個人的には Qiita に投稿するネタができてよかった。

参考

1. 藤縄克之: 環境地下水学, 共立出版, pp.179 - 207, 2010.

2. Yuki Kagawa's Page: 偏微分方程式の離散化, 2022/09/05確認.

3. 松田 安弘, 邵 長城: 二次元移流拡散方程式の有限要素法および差分法による新しい陽的定式化 : 誤差解析手法によるアプローチ, 59(559), pp.833 - 839, 1993.

↓ 全体的にめちゃくちゃ参考にさせて頂きました！！m(_ _)m
【Python】流体シミュレーション：線形から非線形へ

↓ 離散化する手順はこちらをさんこうにしました！
【第7回Python流体の数値計算】１次元拡散方程式を差分法で実装する。