こんにちはJS2IIUです。Pythonで攻略する流体力学シミュレーション連載の第11回です。
前回、第10回では偏微分方程式（PDE）の数学的・物理的な分類について学びました。定常状態を表す「楕円型」、なめらかな広がりを見せる「放物線型」、そして情報を波として運ぶ「双曲型」。これらの特性を理解することは、シミュレーションの設計図を描く上で非常に重要です。
今回扱うテーマは、流体計算において最も基本的でありながら、同時に最も手強い相手の一つである「移流拡散方程式（Advection-Diffusion Equation）」です。これは、第10回で学んだ「移流（双曲型）」と「拡散（放物線型）」が同時に起こる現象を記述した方程式です。
これまで学んできた中心差分をそのまま適用すると、計算が「振動」してしまうという奇妙な現象に直面します。この問題をいかに解決するか。流体計算の実戦的なテクニックである「上流差分法」について、Pythonでの実装を通じて深く掘り下げていきましょう。
1. 移流拡散方程式とは何か
移流拡散方程式は、流体力学や物質輸送現象の基礎となる偏微分方程式です。物理量（例えば温度、濃度、速度など）が「流れ（移流）」によって運ばれつつ、「分子運動や乱流（拡散）」によって周囲に広がる様子を同時に記述します。
現象のイメージ
川の流れにインクを垂らすと、インクは川下へ運ばれ（移流）、同時に水中にじんで広がります（拡散）。このような「運ばれながら広がる」現象は、空気中の煙や香水、海洋の塩分・温度分布など、自然界の多くの場面で見られます。
数式による定式化
1次元の移流拡散方程式は次のように表されます。
$$
\frac{\partial \phi}{\partial t} + u \frac{\partial \phi}{\partial x} = \nu \frac{\partial^2 \phi}{\partial x^2}
$$
ここで、
- \(\phi\) (ファイ): 輸送される物理量（例：温度、濃度、速度場など）
- \(u\): 流速（移流の速さ、単位は m/s）
- \(\nu\) (ニュー): 拡散係数（物質拡散なら分子拡散係数、流体なら動粘性係数、単位は \(m^2/s\)）
各項の物理的意味
- \(\frac{\partial \phi}{\partial t}\): 時間変化（物理量が時間とともにどう変化するか）
- \(u \frac{\partial \phi}{\partial x}\): 移流項（流れによる物理量の輸送。流速\(u\)の方向に情報が伝播する）
- \(\nu \frac{\partial^2 \phi}{\partial x^2}\): 拡散項（物理量が周囲に広がる効果。濃度勾配があると均一化しようとする）
方程式の成り立ち
この式は、ナビエ・ストークス方程式の「移流項」と「粘性項」を抽出したもので、流体の運動や物質の輸送現象の本質を簡潔に表現しています。移流と拡散が同時に働くことで、物理量の分布は「流れに沿って運ばれつつ、なだらかに広がる」特徴的な振る舞いを示します。
応用例
- 大気や海洋の温度・塩分分布の予測
- 汚染物質の拡散シミュレーション
- 化学反応器内の濃度分布解析
この方程式を数値的に解くことは、現代の流体シミュレーションや環境予測、工学設計の基盤となっています。以降では、この方程式を安定かつ正確に解くための数値手法について詳しく解説します。
例えば、川の流れ（移流）にインクを垂らしたときの様子を想像してください。インクの塊は川下に流されていきますが、同時に水の中ににじんで広がっていきます。あるいは、部屋の中を漂う香水の匂いや、煙突から出た煙の広がりも同様です。
1次元の移流拡散方程式は、以下のように書くことができます。
$$
\frac{\partial \phi}{\partial t} + u \frac{\partial \phi}{\partial x} = \nu \frac{\partial^2 \phi}{\partial x^2}
$$
- \(\phi\) (ファイ): 輸送される物理量（温度、濃度、速度など）。
- \(u\): 流速（移流の速さ）。
- \(\nu\) (ニュー): 拡散係数（または動粘性係数）。
この式は、ナビエ・ストークス方程式の「移流項」と「粘性項」をシンプルに抽出した形になっており、流体計算のエッセンスが詰まった方程式と言えます。
2. 中心差分の罠：数値的振動
第6回で学んだ通り、差分法において「中心差分」は2次精度を持つ非常に優れた近似手法でした。一見すると、移流項 \(u \frac{\partial \phi}{\partial x}\) に対しても中心差分を適用するのが最良の選択に思えます。
しかし、移流が非常に強い（(u) が大きく、(\nu) が小さい）状況、すなわち拡散よりも移流が支配的な場合に中心差分を用いると、計算結果に不自然なギザギザの波（高周波成分）が現れることがあります。これを数値的振動（Numerical Oscillation）と呼びます。
なぜ振動するのか？
この現象の本質は、数値安定性と情報伝播の方向性にあります。中心差分は「右側の情報」と「左側の情報」を対称的に扱うため、物理的な流れの方向（上流→下流）を無視してしまいます。移流項は本来、情報が上流から下流へ一方向に伝わる性質（双曲型PDEの特徴）を持っていますが、中心差分は下流側の値も参照してしまうため、物理的な因果関係と計算上の参照関係が矛盾します。
この安定性を測る指標として、セルレイノルズ数（またはセルペクレ数） \(P_e\) が重要になります。
$$
P_e = \frac{|u| \Delta x}{\nu}
$$
理論的には、この値が 2 を超えると、中心差分による計算は振動しやすくなることが知られています。格子間隔 \(\Delta x\) を非常に細かくすれば回避できますが、それでは計算コストが膨大になってしまいます。
この矛盾は、離散化した差分方程式の安定性解析（例えばフーリエモードによるvon Neumann安定性解析）を行うと明確になります。中心差分を用いた移流項の離散化では、格子上の高周波成分（ジグザグなモード）が減衰せず、むしろ増幅される場合があることが分かります。特に、拡散項が弱い場合（\(\nu\) が小さい）、この不安定性が顕著に現れます。
この安定性の指標として、セルレイノルズ数（Cell Reynolds Number, またはセルペクレ数） \(P_e\) が重要です。
3. 上流差分法の魔法
この数値的振動の問題を根本的に解決するために考案されたのが、上流差分法（Upwind Scheme）です。
上流差分法の本質は、「物理的な情報の伝播方向（流れの向き）に合わせて、数値計算でも上流側の値のみを参照する」ことにあります。これにより、計算上の因果関係と物理現象の因果関係が一致し、数値的な不安定性を抑えることができます。
数式による離散化の違い
例えば、一次元の移流項 \(u \frac{\partial \phi}{\partial x}\) を離散化する場合、
- (u > 0)（右向きの流れ）のときは、後退差分（backward difference）を用います：
$$u \frac{\partial \phi}{\partial x} \approx u \frac{\phi_i – \phi_{i-1}}{\Delta x}$$ - (u < 0)（左向きの流れ）のときは、前進差分（forward difference）を用います：
$$u \frac{\partial \phi}{\partial x} \approx u \frac{\phi_{i+1} – \phi_i}{\Delta x}$$
このように、流速の符号によって参照する格子点の方向を切り替えることで、物理的な情報の伝播方向（上流→下流）と数値計算の参照方向が一致します。
安定性の観点
上流差分法は1次精度ですが、中心差分法と異なり、移流が支配的な場合でも高周波成分の数値的不安定性（振動）が発生しません。これは、上流差分の離散化によって、数値的に「人工的な拡散項（数値粘性）」が付加されるため、ジグザグなモードが抑制されるからです。
また、von Neumann安定性解析を行うと、上流差分法では全ての波数成分に対して減衰効果が働くことが示されます。これにより、たとえセルレイノルズ数 (P_e) が大きくても、計算が爆発したり振動したりすることなく、安定して解を求めることが可能になります。
物理的意味と数値的意味
上流差分法は、物理現象の「因果律」を数値計算に忠実に反映する手法です。流れのある系では、情報は必ず上流から下流へ伝わるため、下流側の値を参照しないことで、物理的な矛盾や数値的な不安定性を回避できます。
このように、上流差分法は「安定性」と「物理的整合性」を両立するための基本的かつ強力な手法として、数値流体力学の基礎となっています。
4. Pythonによる実装：中心差分 vs 上流差分
それでは、PyTorchを用いて、中心差分と上流差分の挙動を比較するシミュレーションを行ってみましょう。
厳しい条件（高いセルレイノルズ数）において、両者がどのような結果を出すかを確認します。
import torch
import matplotlib.pyplot as plt
def simulate_transport(nx=100, method='upwind'):
"""
移流拡散方程式を計算する
method: 'central' (中心差分) or 'upwind' (上流差分)
"""
# パラメータ設定
L = 2.0
dx = L / (nx - 1)
u = 1.0 # 移流速度 (右向き)
nu = 0.005 # 拡散係数 (あえて小さくして移流を強くする)
dt = 0.001
nt = 300 # ステップ数
# セルレイノルズ数の表示
pe = (u * dx) / nu
print(f"Method: {method}, Cell Reynolds Number: {pe:.2f}")
x = torch.linspace(0, L, nx)
# 初期条件: 矩形波
phi = torch.zeros(nx)
phi[(x >= 0.4) & (x <= 0.8)] = 1.0
phi_init = phi.clone()
for n in range(nt):
phi_new = phi.clone()
# 1. 拡散項 (中心差分: 常に安定)
diffusion = nu * (phi[2:] - 2*phi[1:-1] + phi[:-2]) / dx**2
# 2. 移流項 (手法の切り替え)
if method == 'central':
# 中心差分: (phi[i+1] - phi[i-1]) / 2dx
advection = u * (phi[2:] - phi[:-2]) / (2 * dx)
else:
# 上流差分 (u > 0 なので後退差分を使用): (phi[i] - phi[i-1]) / dx
advection = u * (phi[1:-1] - phi[:-2]) / dx
# 更新
phi_new[1:-1] = phi[1:-1] + dt * (diffusion - advection)
phi = phi_new
return x, phi_init, phi
# --- 実行と可視化 ---
nx_test = 100
x, init, res_central = simulate_transport(nx=nx_test, method='central')
_, _, res_upwind = simulate_transport(nx=nx_test, method='upwind')
plt.figure(figsize=(12, 6))
plt.plot(x.numpy(), init.numpy(), 'k--', label='Initial (Rectangular Wave)')
plt.plot(x.numpy(), res_central.numpy(), 'r-', label='Central Difference (Unstable)')
plt.plot(x.numpy(), res_upwind.numpy(), 'b-', label='Upwind Scheme (Stable)')
plt.title("Advection-Diffusion: Central vs Upwind")
plt.xlabel("x")
plt.ylabel("phi")
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()
コードの解説と動作結果
- セルレイノルズ数の設定:
このコードでは、\(u=1.0\), \(dx=0.02\), \(\nu=0.005\) としており、\(P_e = 4.0\) となります。これは、中心差分が不安定になる境界（\(P_e=2\)）を超えた設定です。 - 移流項の計算:
method='central'の場合は、現在の点phi[1:-1]の両隣の差をとっています。一方、method='upwind'では、流速 \(u\) が正であることを利用して、自分自身の点と一つ左（上流）の点の差だけを使っています。 - 結果の観察:
赤い線（中心差分）を見てください。波の前後、特に後ろ側に大きな振動（オーバーシュートとアンダーシュート）が発生していることがわかります。これは物理的にはあり得ない負の値まで取ってしまっています。
一方、青い線（上流差分）は、全く振動することなく波を安定して運んでいます。これが上流差分法の最大のメリットである「安定性」です。
5. 数値拡散：上流差分がもたらす副作用
上流差分法は万能ではありません。安定性を手に入れる代わりに、私たちは大きな代償を払っています。それが数値拡散（Numerical Diffusion）です。
先ほどのグラフの青い線（上流差分）をよく見ると、元の矩形波に比べて、波の角が丸まり、全体的に「ボケて」しまっていることに気づきませんか？
実は、上流差分法で一階微分を近似することは、数学的に見ると、「真の微分」に「余計な二階微分（拡散効果）」をこっそり付け加えた計算を行っているのと等価なのです。この計算によって勝手に付け加えられた偽の拡散成分を「数値拡散」と呼びます。
数値拡散のジレンマ
- 中心差分: 精度は高い（2次精度）が、振動に弱く不安定。
- 上流差分: 非常に安定（1次精度）だが、数値拡散によって解がなまってしまう。
「計算を安定させたいが、解の鮮明さも保ちたい」というジレンマは、数値流体力学の歴史において長年の課題でした。第12回で扱うような、より高度な「高精度スキーム」は、この数値拡散をいかに抑えつつ安定性を確保するかという工夫の結晶なのです。
6. まとめ
第11回では、移流と拡散のせめぎ合いをどう計算するかを学びました。
- 移流拡散方程式は、流体輸送の最も基本的なモデル。
- 中心差分は高精度だが、移流が強いと数値的振動を起こす。
- 上流差分法は、情報の流れる向きを意識することで、圧倒的な安定性を実現する。
- ただし、上流差分には解をなまらせる数値拡散という副作用がある。
物理的な情報の流れ（因果関係）を計算手法に組み込むという「上流差分」の考え方は、複雑な3次元流体シミュレーションにおいても根底に流れる重要な哲学です。
次回、第12回では、この移流の難しさがさらに牙を向く「非線形移流」の世界へ足を踏み入れます。衝撃波という、流体力学における最もダイナミックな現象をどのようにコンピュータで扱うのか。より高度な数値スキームと共に解説していきます。
最後まで読んでいただきありがとうございます。
コメント