移流拡散方程式 - 明後日の方向にただ全力

ラプラス方程式がとてもお手軽だったので、調子にのって移流拡散方程式を解いてみます。式は

${\frac{\partial \phi}{\partial t}+u_j\frac{\partial \phi}{\partial x_j}-k\frac{\partial^2 \phi}{\partial x^2} = 0}$

ですね。大気中とか海中の放射能汚染シミュレータとかも原理はこれなんだろうか？大気中に放射能が流れに運ばれながら拡散していくので、たぶんこれかな。

これまでと同様、重み付け残差法で弱形式化して、ガラーキン法で離散化すると、いつも通り連立方程式の形になるので、あとはこれを解く。

…と、解が振動して計算できない。まじか。

原因はおそらく移流項のせい。ガラーキン法は数学的には差分法で言うところの中心差分に等しい。こいつはレイノルズ数ていうか、ペクレ数が高くなるとすぐ振動し始める厄介なやつで、メッシュのサイズをかなり小さくしないとまともに計算できそうにない。差分でコード書いてたときは、一次風上やら三次風上、kkスキームにQUICK法、CIP法まで試して遊んでましたが、有限要素法ではどうしたものか。

…結局差分と同じで、風上化を考えます。ただ人工粘性を付加するだけだと精度が落ちるので、重み関数の形を変えて、移流項に対し、流線方向に人工粘性が加わるように変換してやればいいらしい。SUPG法です。ちょっと名前が恰好よくていい感じです。

平たく言えば、重み関数の形が従来の重み関数＋安定化項という風に変わるだけ。離散化してやれば、最終的にはガラーキン法で導出した行列方程式に安定化項が一個付け足された形になるので、それを解いたら終わり。

ただ安定化パラメータのとり方にいろいろ工夫があるらしく、この部分はちょっと真面目に論文か教科書読んだ方がいいかな…。とりあえず本の通りに計算してみれば完成。
これもラプラス方程式のコードができてれば、付け足せばすぐできるので、大した時間もかからず終了。非定常計算も同様、実装は簡単。

このへんも教科書通り。安定化パラメータの導出とかは、この本にわりと詳しく書いてあった。