レッスン1.3 — 数値安定性・シミュレーション精度(レッスン1.2からの続き)

前回までのレッスンでは、シミュレーションの始め方(レッスン1.1)やいくつかの画面操作(レッスン1.2)について学びました。さらにレッスン1.2では、デフォルトのシミュレーション設定では非物理的な数値的振動が現れることを確認しました。今回のレッスンでは、数値安定性と精度を高めるために、どのようにgrid.txtを設定するかについて学びます。

一般的にCFD(数値流体力学)においては精度は計算資源(計算時間、メモリー、CPU数)と正の相関があります。多くの場合、計算精度と使える計算資源を上手くバランスさせるためにどちらか、あるいは両方を妥協する必要があります。このことを頭の隅に置き、先に進みましょう。

数値振動の抑制

まず最初にレッスン1.2で明らかになった数値振動について取り組みます。まずメインフォルダにあるデフォルトの2Dチャネル流れ用のgrid.txtを開いて見てください(ダウンロードはこちらから)。レッスン1.1で触れたようにgrid.txtではシミュレーション上必要なパラメーターを指定します。セパレーター(区切り行)以外の全ての行は3つのブロックから構成されています。例えば二つ目のセパレーターの直後には以下のような行があります。

1
01:cmode        0  // Simulation mode, ....

最初のブロック (01:cmode) は各パラメーターの名前になっています(ほとんどの場合、ユーザーガイドの数式中の変数名と類似)。2番目のブロック (“0” 上の例)はそれぞれのパラメーターの値です。そして3番目のブロックはそれぞれのパラメーターに関するコメントで、シミュレーションに直接影響しません。これらのブロックはスペースで仕切られている必要があります。

次にgrid.txt内(上方)で以下の行を見つけてください。

1
2
3
4
02:nx         384  // No. grid points in x
03:ny         128  // No. grid points in y
04:lx       0.015  // Domain x-size
05:ly       0.005  // Domain y-size

3番目のブロック内のコメントが説明しているように、nx及びnyはそれぞれx方向、y方向の格子点数の数であり、lx及びlyはそれぞれx方向、y方向のシミュレーション領域の長さです。単位はSIユニット(メーター)を使います。Flowsquareでは差分法を使って計算するので格子点数は一般的に多ければ多いほど良いです。しかし、格子点数の増加は2乗で計算時間の増加に寄与します。従って、格子点数と計算時間との間で上手いバランスをとる必要があります。また、特別な場合を除いて、一般的にはx方向及びy方向の格子点密度(それぞれlx/nx、ly/ny)は等しくあるべきです。

さて、さらにgrid.txtを見ると、以下のような行があります。

1
2
10:nfil         0  // Interval time steps for filtering
11:wfil         0  // Relaxation parameter for filtering

これらの行は直接今回の問題(数値振動)へ関係します。Filtering(フィルター)とは、流れの解にいわゆる人口粘性を与え、非物理的な小さな振動を取り除く、という操作です。上の2行を下のように書き換え、今までと同様にシミュレーションを行いましょう。このフィルターを使うシミュレーションのケース名は Ch0_filter にします。

1
2
10:nfil         1  // Interval time steps for filtering
11:wfil         1  // Relaxation parameter for filtering

だいたい、600タイムステップ程計算したところで(1分程度)、シミュレーションを一旦停止(ESCキー)し、v速度成分(縦方向速度成分)を表示してみましょう(操作方法を忘れたらレッスン1.2を参照)。前回フィルターなしで行った計算結果(レッスン1.2の図9)と比べると、見事に数値振動が無くなっていることが確認できます(図1)。

図1:フィルター使用時のv速度成分分布。

図1:フィルター使用時のv速度成分分布。

一般的にほとんど全てのシミュレーションにおいてフィルターは必須となるでしょう。しかし、wfil = 1はもしかしたら人口粘性を与えすぎかもしれません。計算が発散せずに行える範囲で出来るだけ小さなwfil(0.01–0.1)を使うのが望ましいです。

計算精度

grid.txt内でシミュレーション精度を指定する以下の行を見つけてください。

1
09:iorder       0  // 0: low order, 1: high order, ....

iorderには0から3までの数字を指定することができ、それぞれの数字は以下の数値計算法を表しています。

  • iorder=0: 低精度計算法;2次精度中心差分+オイラー法(1次精度)
  • iorder=1: 高精度計算法;4次精度中心差分+3次精度ルンゲクッタ法
  • iorder=2: 2次精度中心差分+Lax-Wendroff (中間点)法
  • iorder=3: 4次精度中心差分+Lax-Wendroff (中間点)法

ここでは、以下のようにiorderとして1を指定しましょう。これにより、シミュレーションは高精度計算法を用いて行われます。

1
09:iorder       1  // 0: low order, 1: high order, ....

また、高精度計算を行うために、人口粘性を取り除く必要があります。先ほど変更したフィルター関係の行を以下のように戻してください。

1
2
10:nfil         0  // Interval time steps for filtering
11:wfil         0  // Relaxation parameter for filtering

これで、高精度シミュレーションを実行する準備が整いました。プロジェクト名を“Ch1”とし、シミュレーションを実行しましょう。高精度シミュレーションは少し時間がかかります(4000ステップまでおよそ3-5分)。シミュレーション途中で一旦停止して、様々な物理量のグラフを描くのもいいかもしれません。シミュレーションは4000ステップで停止するはずです。次回のレッスンで、いよいよ高精度シミュレーション結果を解析します。最後まで読んでくださりありがとうございます。