【CFD/Python】CIP法【CFDプログラミング講座　#6　移流項編】

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

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

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

重要な項目は下記の通りです。 二次風上差分とQUICKは風上差分を基とした手法です。参照点を増やすことで、高精度化できることがわかりました。ただし、安定性は悪化する結果となりました。

• 二次風上差分は、風上差分の参照点を増やすことで精度を上げた手法である
• 二次風上差分は、３つの参照点から、２つの1/2点を作成し、最後の計算点を計算する
• QUICK法は、４つの参照点から最後の計算点を得る手法である
• QUICK法は、４つの参照点から、２つの1/2点を作成し、最後の計算点を計算する
• QUICK法は、1/2点の計算に３つの計算点を使用しているので、非線形補完できる

参照点が増えると良い結果が得られる傾向がわかっていただけたと思います。基本的には参照点を増やせば精度が上がるので、マシンパワーが許す限り高精度の離散化スキームを使うと良いです。

ただ、良い結果を得る方法は多項式補完だけではありません。代表的な手法として、今回CIP法について説明します。

CIP法

CIP（Constrained Interpolation Profile scheme）法は、計算点の値とその傾きを使用した手法です。

CIP法を学ぶ前に、まずは多項式の課題を振り返りましょう。

多項式による問題

多項式による高精度化は、どうしても参照点を増やす必要がありました。二次風上差分は３点、QUICK法は４点の参照点を必要としました。

参照点が増えることの問題は、境界点における計算です。今回は、境界に影響のない計算を行っているので特に問題ないですが、実際の計算では境界は重要な問題です。

例えば、流入条件や流出条件、周期境界などをどのように表すかが多項式化すると難しくなります。また、多項式化すると周囲に無駄な点ができて、計算を圧迫します。

多項式が進めば、より遠くの点を参照するようになります。前回の二次風上差分とQUICKで振動が発生したのを見ると、多項式だけでは計算不安定による問題は解決しなさそうです。

CIP法の理論

それに対して、CIP法は計算点の値とその傾き（微分値）から計算を行うので、少ない点から多くの情報を得ることができます。

下記の図のイメージです。$ u_i $ を計算する際に、 $ u_{i-1} $ と $ u_{i+1} $ の値と傾きを使用して計算します。

CIP法は理論的には非常にシンプルですが、計算は難しいです。そのまま式を使用すれば使えるので、細かく理解する必要はありません。全体像として捉えておくと良いでしょう。

CIP法の計算式

CIP法では、下記の基本形関数を使います。

$$ F(\xi) = a \xi^3 + b \xi^2 + c \xi + d $$

a,b,c,dは$ \xi $の係数です。後で、$u_i$や $u_i$ の微分値である$g_i$によって表します。

$ \xi $ は１ステップの移動距離を表します。図で示すと下記のようになります。

$ \xi $ は１ステップの移動距離なので、位置 $ x_{i-1} $ と $ x_{i} $ の差分、もしくは 流速 $ c $ と時間刻み幅 $ \Delta t $ の乗算として表せます。ここでは、どちらも既知の値である $ -c \Delta t $ を使用します。

線形流れにおいて$ \xi $ は定数なので、 $ -c \Delta t $ は最後に代入してやれば良いです。そのため、ここからはわかりやすさも考えて、 $ \xi $を使用します。

$ F(\xi) $ を使用して、物性値uを輸送する式は下記のようになります。

$$ u_i = a\xi^3 + b\xi^2 + g_i \xi + u_i $$

ポイントは、cが$u_i$の微分である$ g_i$になったことと、最後のdが$u_i$になったことです。

$u_i$の計算のためには、$g_i$の計算が必要であることがわかりました。そのため、物性値の微分値gの計算を見てみましょう。

$$ g_i = 3a \xi^2 + 2b\xi + g_i $$

$F(\xi)$ が $ \xi $で微分されていることがわかりますね。こちらでも$u_i$の計算と同様に、cの代わりに$g_i$が出てきました。

あとはaとbの計算ですが、下記のように定義されています。cとdも振り返りを兼ねてまとめておきます。

$$ a=\frac{g_i + g_{i-1}} {\Delta x^2} – \frac{2(u_i – u_{i-1})}{\Delta x^3} $$

$$ b = \frac{3 (u_{i-1} – u_i)}{\Delta x^2} + \frac{2g_i + g_{i-1}}{ \Delta x} $$

$$ c = g_i = \frac{\partial u(\xi)}{\partial \xi} $$

$$ d = u_i $$

全て$u_i$と$g_i$の関数となっており、前ステップの $u_i$と$g_i$ を使用して計算します。

CIP法の計算フロー

必要な計算式は上記で全て出てきました。計算フローについて一旦振り返りましょう。

１．定数 $ \xi $を計算する

２．前ステップの$u_i$ と $g_i$ から a,b,c,d を計算する

３．微分値である$ g_i$ を計算する

４．物性値 $ u_i $ を計算する

線形流れの場合は、流速が定数になるので計算が非常に簡単です。ポイントとしては、a,b,c,dの計算において前ステップの値を使用するので、初期値を入れる必要があることです。$g_i$ の初期値を設定するのを忘れないようにしましょう。

CIP法プログラム

CIP法のコードは下記の通りです。

import numpy 
import matplotlib.pyplot as plt
%matplotlib inline


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

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

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

plt.show()

出力結果は下記のようになります。

CIP法を使えば、今までの多項式と比べて非常に良い結果が得られます。

急激に値の変化する流れにおいて、傾きの情報がいかに重要化がわかっていただけるでしょう。

ライブラリ

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

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
xi = -c * dt

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

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

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

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

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

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

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

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

xiは $ \xi $ のことです。1ステップの移動距離を表すので、$ -c \Delta t$ の値を代入します。

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

gはuの微分値（傾き）です。初期値はゼロで大丈夫です。

移流方程式（CIP法）

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

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

for n in range(nt):  
    u_old = u.copy() 
    g_old = g.copy() 
    
    a = ((g_old[1:-1] + g_old[:-2]) / (dx**2) 
         - 2 * (u_old[1:-1] - u_old[:-2]) / (dx**3))
    b = (3 * (u_old[:-2] - u_old[1:-1]) / (dx**2) 
         + (2 * g_old[1:-1] + g_old[:-2]) / dx)
    u[1:-1] = (a * (xi**3) + b * (xi**2) 
               + g_old[1:-1] * xi + u_old[1:-1])
    g[1:-1] = (3 * a * (xi**2) + 2 * b * xi 
               + g_old[1:-1])
    
    ・・・

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

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

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

次の行からは、移流方程式をCIP法で書いた式です。まずは、a,bの計算のために、下記の式をプログラム化しています。

$$ a=\frac{g_i + g_{i-1}} {\Delta x^2} – \frac{2(u_i – u_{i-1})}{\Delta x^3} $$

$$ b = \frac{3 (u_{i-1} – u_i)}{\Delta x^2} + \frac{2g_i + g_{i-1}}{ \Delta x} $$

cは$g_i$、dは$u_i$を使用するので、これでa,b,c,dが全て揃ったことになります。

次に、$u_i$と$g_i$の計算を行います。下記の式がプログラムにかかれています。

$$ u_i = a\xi^3 + b\xi^2 + g_i \xi + u_i $$

$$ g_i = 3a \xi^2 + 2b\xi + g_i $$

今回は、u_oldとg_oldとして古い値を使用しているので、$u_i$と$g_i$の順番はどちらでも構いません。違和感がある場合は、$g_i$を先に計算しても良いです。

配列の範囲が1:-1となっているのは、領域の両端の２つの節点では移流方程式を計算していないからです。両端の節点は参照点としてだけ使用します。CIP法では、参照点として両端の一つしか使わないため、非常に無駄の少ない計算になります。

グラフ描画

グラフ描画には、初めにインポートしたグラフ描画ライブラリの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でまとめて表示しています。

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

少しだけ振動が見られますが、多項式の離散化スキームに比べると圧倒的に良い結果が得られていることがわかるでしょう。

おわりに

今回は CIP（Constrained Interpolation Profile scheme）法について解説しました。

CIP法は、左右の計算点の値とその傾きを使用した手法でした。重要なポイントは下記のとおりです。

• CIP法は、少ない参照点から多くの情報が得られる
• CIP法は傾きの情報があるため、急激な値変化に対する安定性が高い
• CIP法では１ステップの移動距離$\xi$を使用する
• CIP法では、全てのセルにおける物性値の微分値$g_i$を計算する

理論的な理解としては簡単な内容ですが、実際の計算フローやプログラムとしてはやや複雑です。

CIP法は結果を見れば分かる通り、非常に優れた手法です。しかし、全てが解決するわけではありません。当然課題もあります。

例えば、CIP法は微分値を使うため、二次元にしたときに複雑な式になります。これはフラクショナルステップ法などによって簡易化が行われます。しかし、フラクショナルステップ法は精度低下の原因にもなるので、このあたりは労力と精度の優先度の比較が必要です。

また、CIP法は全てのセルの微分値の情報を記憶するので、メモリを圧迫します。精度を上げるためなので仕方ないですが、大規模計算においては足を引っ張る原因にもなり得ます。

他にもCIP法は体積保存されないという問題もあったりします。

これらはデメリットの一部ですが、デメリットを差し引いてもCIP法は非常に優れた手法です。そのため、今回可視化で体験したオーバーシュートを抑えるRCIP（Rational CIP）法という手法や、拘束条件を追加することで体積を保存するCIP-CSL法などがあります。

このような改良手法を使うことで、CIP法のデメリットは多少抑えられます。離散化スキームの選択肢としては非常に優れたものになるでしょう。

次回

次回はTVDスキームについて解説します。１次元の移流項では最後の解説です。

TVDスキームは、今までの多項式を良い感じに変更することで、精度と安定を両立した手法です。そのため、CIP法のような変化球というよりは、王道の手法をうまくツール的に使用した感じです。

最も使われている手法といっても過言ではないでしょう。実際、OpenFOAMなどの流体ソフトウェアにも入っています。

下記からどうぞ。

youtubeもやってます

youtubeでも解説しています。動画のほうが理解しやすいかたは、こちらもどうぞ。