【CFD/Python】Lax-Wendroff法【CFDプログラミング講座　#4　移流項編】

この講座ではPythonを使用したCFDプログラミングを解説していきます。最後まで受講して頂くと、差分法による流体解析ができるようになるようにしています。ぜひ最後までお付き合いください。

第２回からは移流項に関する解説をしていきます。移流項は最も難しい部分ですが、特に重要なので先に解説します。

前回の第３回では移流方程式の代表的な離散化スキームである中心差分と風上差分について解説しました。まだ見てない方は下記からどうぞ。

重要な項目は下記の通りです。風上差分と中心差分は簡単ですが安定性と精度に極端なスキームです。これらをどう工夫して精度と安定性を両立した手法を考えるかがここからの課題です。

• 風上差分は数値的に安定するが、誤差が大きい
• 風上差分は誤差により波がなめらかになる
• 中心差分は誤差は少ないが、数値振動が起きる
• 中心差分は計算が不安定なので、実用性はない
• CFDプログラミングは、「変数定義→初期化→実行→可視化」の流れである
• CFDプログラミングには、numpyとmatplotlibを使う
• numpyは配列ライブラリでmatplotlibはグラフ描画ライブラリである

Lax-Wendroff法

Lax-Wendroff法は２次多項式を使うことで数値振動を抑制する手法です。

離散化スキームは高次方程式になればなるほど精度は上がります。なぜでしょうか？

それは、離散化が連続的な微分を線形な形に置き換えているからです。二次・三次と補間関数が複雑になれば、もとの連続的な式も復元しやすくなります。そのため、補完関数は高次であればあるほど精度が上がります。

ただし一方で、高次方程式になるだけでは安定性は必ずしも上がりません。精度の高い中心差分よりも風上差分のほうが安定性が高いのは、誤差として人工粘性が含まれるためです。人工粘性が振動を抑えてくれるのです。

つまり、高次方程式に対して粘性成分を与えることで、精度と安定性を両立した手法を作ることができます。それが Lax-Wendroff法 です。

移流方程式を下記に示します。

$$ \frac{\partial u}{\partial t} + c \frac{\partial u}{\partial x} = 0$$

Lax-Wendroff法 では左辺の$ \frac{\partial u}{\partial x} $ を下記の２次多項式で置き換えます。

$$ u_{x} = a (x – x_i)^2 + b (x – x_i) + c $$

$$ a = \frac{u_{i+1} – 2 u_{i} + u_{i-1}}{ 2 {\Delta x}^2 } $$

$$ b = \frac{ u_{i+1} – u_{i-1}}{ 2 \Delta x} $$

$$ c = u_{i} $$

上記の式で理解しても良いですが、Lax-Wendroff法はもっとわかりやすい書き方があります。

上記の式を展開すると、中心差分とそれ以外の成分に分けることができます。

$$ u^{n+1} = u_{i} – c \Delta t \frac{u_{i+1}-u_{i-1}}{2 \Delta x} + c^2 {\Delta t}^2 \frac{u_{i+1} -2 u_{i} + u_{i-1}}{2 {\Delta x}^2} $$

右辺第３項以外を見てください。中心差分の形になっています。そして、右辺第３項は、二次の拡散を与える項となっています。

つまり、 Lax-Wendroff法 は中心差分に対して２次の粘性を与えた式ということです。中心差分に二次の項を入れているので２次精度を担保できますし、二次の粘性項があるので安定した計算ができるようになっています。

仕組みがわかると多項式も簡単でしょう。このように項が増えることで、精度を上げつつ、粘性を与えて計算を安定させる事ができます。

Lax-Wendroff法プログラム

Lax-Wendroff法は下記のコードで表せます。風上差分や中心差分と異なるのはメインループ内の移流方程式だけですが、一応振り返りも兼ねて全部解説します。

import numpy 
import matplotlib.pyplot as plt
%matplotlib inline


nx = 501
dx = 2 / (nx - 1)
nt = 1000 
dt = .001
c = 1

u = numpy.zeros(nx) 
u[int(0.2 / dx) : int(0.5 / dx + 1)] = 1  

for n in range(nt):  
    
    u_old = u.copy() 
    u[1:-1] = (u_old[1:-1] - c * dt / (2 * dx) * (u_old[2:] - u_old[:-2])
    + c**2 / 2  * (dt / dx)**2 * (u_old[2:] - 2 * u_old[1:-1] + u_old[:-2]))            
    
    if n % 200 == 0:
        plt.plot(numpy.linspace(0, 2, nx), u)

plt.show()

結果は下記のようになります。二次精度を維持しているのに、中心差分よりも振動が抑えられる結果が得られています。これは二次の拡散項のおかげです。ただし、まだ振動が大きく実用とは言えないですね。

ライブラリ

上から順にプログラムを解説していきます。まずは使用するライブラリのインポートからです。

import numpy 
import matplotlib.pyplot as plt
%matplotlib inline

使うライブラリはnumpyとmatplotlibです。import 〇〇と書くことでライブラリをインポートでき、as〇〇と書くことで別名を定義できます。

numpyは配列演算に関するライブラリです。まとめてデータを扱えて計算が高速であることから、数値計算では非常によく使われます。npという省略表記されることも多いです。

pythonにはリストという機能がありますが、今回はnumpyを使用します。numpyでは格納される型が揃っており、他のライブラリとの相性の面でも使い勝手がよいです。

matplotlibはグラフ描画用ライブラリです。計算結果を可視化するために使用します。

今回はmatplotlib内のpylotというメソッドを多用するので、pltという省略名で定義しておきます。これによってplt.〇〇という表記で使用することができ、プログラムがスッキリします。

初期化

ここでは解析の初期条件を与えます。

初期条件として、格子数や時間刻み幅などのパラメータから、どこにどのくらいの速度を与えるかまで決定してやる必要があります。

nx = 501
dx = 2 / (nx - 1)
nt = 1000 
dt = .001 
c = 1

u = numpy.zeros(nx) 
u[int(0.2 / dx) : int(0.5 / dx + 1)] = 1  

nxは格子の節点数です。nxは節点なのでセル数はnx-1個になることは注意してください。今回から説明する差分法では計算点が節点になります。そのため、501点の計算点があると考えてください。

dxは接点間距離です。セルの大きさと言い換えてもらっても良いです。

今回は計算領域の長さを２ｍとするため、このような式となっています。複雑そうに見える場合は下記の節点が少ない場合の例を表した図を見ると、わかりやすいと思います。

節点が両端に位置することになるので、接点間距離はnx-1個の間隔で割ることがわかります。

ntはステップ数です。今回はforループで時間を進めるので、1000回ループ計算を行うことになります。言い換えると、1000回も移流方程式を解くことになります。

dtは時間刻み幅です。１回のステップで進む時間を表します。CFDでは時間刻み幅は非常に重要な要素です。計算の安定化のために調節してやる必要があります。

cは流速（速度）です。線形移流方程式なので、流速は定数になります。

uは輸送される任意の変数です。ナビエ・ストークス方程式なら流速、熱やVOF値の場合もあります。初期値として、0.2~0.5の範囲で1を与えてやります。他は0に設定しているので、 0.2~0.5の範囲だけ1となる矩形波が現れます。

移流方程式（Lax-Wendroff法）

ここからがメインの計算ループになります。

今回は矩形波を時間経過とともに右側に流すプログラムになるので、時間をforループで進めます。

for n in range(nt):  
    
    u_old = u.copy() 
    u[1:-1] = (u_old[1:-1] - c * dt / (2 * dx) * (u_old[2:] - u_old[:-2])
    + c**2 / 2  * (dt / dx)**2 * (u_old[2:] - 2 * u_old[1:-1] + u_old[:-2]))            
    
    ・・・

nステップを0~ ntまで変化させて時間変化を計算します。forループ内に移流方程式が入ります。

u_oldは古い時間ステップにおける輸送変数uの値です。uの値は移流方程式で更新されるため、u_oldを右辺において新しい値はuに代入することとします。

pythonでは配列をまとめて演算するのでこのような一時変数は不要ですが、他の言語にも対応できるように一般的な書き方をしています。

次の行は、移流方程式をLax-Wendroff法で書いた式です。下記の式をプログラム化しています。

$$ u^{n+1} = u_{i} – c \Delta t \frac{u_{i+1}-u_{i-1}}{2 \Delta x} + c^2 {\Delta t}^2 \frac{u_{i+1} -2 u_{i} + u_{i-1}}{2 {\Delta x}^2} $$

配列の範囲が1:-1となっているのは、領域の両端の節点では移流方程式を計算していないからです。両端の節点は参照点としてだけ使用します。

これは、移流方程式をLax-Wendroff法で離散化した際に、計算点の隣の節点を参照するためです。そのため、端の計算点は計算出来ないことになります。

離散化スキームの変更には、この１行を変更すれば良いことがわかります。

グラフ描画

グラフ描画には、初めにインポートしたグラフ描画ライブラリのmatplotlibを使用します。

for n in range(nt):  
    ・・・            
    
    if n % 200 == 0:
        plt.plot(numpy.linspace(0, 2, nx), u)

plt.show()

for文の中にif文があります。ここで、時間ステップが200進むごとに画面出力をしています。

plt.plotはグラフの描画を行います。横軸と縦軸を指定することでグラフを描画できます。

横軸はnumpy.linspaceを使用します。linspaceは指定した範囲を等間隔で分割した配列を作成してくれます。ここでは横軸には0~2までをnx分割し、縦軸はuを使用することにより、縦横で同じ要素数の配列を作成できます。

最後のplt.showはグラフの画面表示です。plotは描画に対してshowは表示です。plotでグラフを重ねて描画し、showでまとめて表示しています。

結果として下記のようなグラフが出力されます。

二次精度を維持しているのに、中心差分よりも振動が抑えられる結果が得られています。これは二次の拡散項のおかげです。ただし、まだ振動が大きいので何とかしたいところです。

おわりに

今回はLax-Wendroff法について解説しました。

Lax-Wendroff法は２次多項式を使った手法であり、多項式の中でももっともシンプルなものです。多項式による高精度化の基本的な考え方が詰まっているので、これからの勉強で非常に役に立つでしょう。

重要なポイントは下記の通りです。

• Lax-Wendroff法は二次多項式を使う手法である
• 二次多項式を使うことで、精度を維持したまま人工粘性を与えられる
• Lax-Wendroff法は中心差分に二次の粘性を加えることで、精度と安定性を両立する
• 高次多項式はLax-Wendroff法の考え方の延長である

特に、精度を維持したまま高次の人工粘性を加えるところがポイントです。高次方程式には共通する考え方となります。

移流項はまだまだ深いです。次回は、二次風上差分とQUICK法について解説します。これらはその名の通り、風上差分を高次にしたものです。Lax-Wendroff法は安定性に欠けていましたが、次回紹介する2つは安定性を重視しています。

結局の所、安定性がないと精度がいくら良くても使えません。ここまでは理論の理解を主にやってきましたが、次回からはいよいよ実用面を意識した手法に触れていきましょう。

Youtubeもやってます

youtubeでも解説しています。ぜひ御覧ください。