Copyright (c) 2020, The MathWorks, Inc.
Navier-Stokes 方程式を数値的に解くシリーズ、第２回目です。
- 非圧縮性 Navier-Stokes 方程式の数値解法１：導入編
- 非圧縮性 Navier-Stokes 方程式の数値解法２：拡散項の陰解法
解析対象については前回の記事（非圧縮性 Navier-Stokes 方程式の数値解法１：導入編）を参照してください。
導入編ではまずキャビティ流れを解いてみることに注力して、かなり荒い離散化を行いました。結果として、時間ステップサイズによっても変わりますが、ちょうどよい Reynolds 数（Re = 500程度）以外では発散する結果になっています。
今回はそこを解決するために以下を実施します。
- 拡散項について陰解法（２次精度クランク・ニコルソン法）を実装
- 対流項について２段階の陽解法（２次精度アダムス・バッシュフォース法）を実装
- 高レイノルズ数でのシミュレーション（下図 Re = 5,000）
実際のところルンゲ・クッタ法の方が一般的かと思いますのでそちらも実装しています。文末 Appendix 参照。
**メモ：**拡散項への陰解法の導入だけで事足りるかと期待したんですが、高レイノルズ数では発散しちゃいました。諦めて移送項も改良することにします。
- MATLAB R2019b
- Signal Processing Toolbox (推奨*)
*) ポワソン方程式を解く部分に Signal Processing Toolbox を使っていますが、直接法や反復法にすれば MATLAB 本体だけでOKです。
この辺はノイマンの安定解析やクーラン数（Courant number）の話。CFD の教科書には必ず出てくるネタですね。ここでは要点だけ触れておくと、拡散方程式のオイラー陽解法・中央差分における安定性条件は
であり、グリッドサイズを 1/2 倍にすると時間ステップ（
また同様に移送方程式についてはクーラン数（c は移送速度）
が絡んできます。ただこちらはグリッドサイズを 1/2 倍にすると時間ステップ（
メモ：ルンゲ・クッタ法を使うとクーラン数 1.2 くらいでも発散はしない印象（あくまで観測ベース）
前回は 1 次精度のオイラーの陽解法で時間積分したので、離散化した Navier-Stokes 方程式を行列で表現した結果、
そして行列で表現すると
こんな感じでした。
ここで、拡散項にクランク・ニコルソン法、対流項にアダムスバッシュフォース法を適用すると
こうなります。式が長くなりましたが基本的な構造は同じ。拡散項に次の時間の速度場（
前回と同じく行列で表現すると
ここで
となります。
また LU 分解してみます（Perot 1993 [1]）
これを厳密に解くのは現実的ではないので、通常
メモ：ちなみに
処理を書き下すと陰解法による部分段階法（Fractional Step Method) の出来上がり。
前回と違う部分は第一式だけ。改めて明記すると、
です。
速度勾配（
メモ：ここで要注目ポイントは
実際のコードは GitHub の方を見て頂くとして、ここでは要点だけ紹介します。ポイントは
- 拡散項に関わる行列の確認
- 境界条件の取り扱い
- 連立方程式の解き方
の３つ。前回のオイラー法との違いは仮の速度場を算出するここだけ。
コードとしてはこんな感じ（updateVelocityField_CNAB.m 参照）
今回
% Implicit treatment for xy direction
b = u(2:end-1,2:end-1) + dt*(-3*Nu+Nu_old)/2 + dt/(2*Re)*(Lux+Luy+Lubc);
xu = getIntermediateU_xyCNAB(u, b, dt, Re, Nx, Ny, dx, dy);
b = v(2:end-1,2:end-1) + dt*(-3*Nv+Nv_old)/2 + dt/(2*Re)*(Lvx+Lvy+Lvbc);
xv = getIntermediateV_xyCNAB(v, b, dt, Re, Nx, Ny, dx, dy);
u(2:end-1,2:end-1) = xu;
v(2:end-1,2:end-1) = xv;
ちなみに前回（updateVelocityField_Euler.m 参照**）**
u(2:end-1,2:end-1) = u(2:end-1,2:end-1) - dt*Nu + dt/Re*(Lux+Luy);
v(2:end-1,2:end-1) = v(2:end-1,2:end-1) - dt*Nv + dt/Re*(Lvx+Lvy);
getIntermediateU_xyCNAB.m と getIntermediateV_xyCNAB.m にすべてが詰まっています。CNAB は Crank-Nicolson, Adams-Bashforth の略です。
行列
ということで、２次精度であれば局所的には三重対角行列ですが、x, y の両方向に対して陰解法とするとなると、もう少し広がりのあるスパース行列となります。行列のサイズは
と要素の数的には
・・まぁまぁ大変そうですね。
例えば境界層流れなどで一次元方向のみにグリッドを細かくするような場合は、その方向の拡散成分にだけ陰解法を用いる (Simens et al. 2009 [2]) というケースもあります。その場合は大規模な連立方程式とはいえ、３重帯行列（２次精度の場合）なので直接法でも十分速いので、計算速度が気になる場合にはこの辺を妥協するのもアリです。
境界条件に由来する
ここでオレンジ色の矢印で表される仮想速度は以下のように計算されます。例えば、
境界を挟んだ２つの速度の平均値（2次精度の内挿）が境界値となるという想定です。
この計算に領域内の速度（ここでは
と展開され、y 方向の
また x 方向の拡散項は左側の境界付近において
であり、仮想速度は必要ないので x 方向の
速度
小さいグリッド数で、実際にどんな行列か見てみます。getL4u.m は
clear, close all
addpath('../functions/')
Lx = 1; Ly = 1;
nx = 5; ny = 5;
dx = Lx/nx; dy = Ly/ny;
maskU = false(nx+1,ny+2);
maskU(2:end-1,2:end-1) = true;
L4u = getL4u(nx,ny,dx,dy,maskU);
spy(L4u);
スパースですね。参考まで、冒頭の数値を確認。
full(L4u(1:nx,1:nx))ans = 5x5
-125.0000 25.0000 0 0 25.0000
25.0000 -125.0000 25.0000 0 0
0 25.0000 -125.0000 25.0000 0
0 0 25.0000 -125.0000 0
25.0000 0 0 0 -100.0000
メモ：ここで maskU は速度
maskUmaskU = 6x7 の logical 配列
0 0 0 0 0 0 0
0 1 1 1 1 1 0
0 1 1 1 1 1 0
0 1 1 1 1 1 0
0 1 1 1 1 1 0
0 0 0 0 0 0 0
こんな行列。境界（0）で囲まれているだけのシンプルなものですね。
連立方程式
上で見た通り
dt = 0.01; Re = 100;% お試し設定
A4u = eye(size(L4u),'like',L4u)-dt/(2*Re)*L4u; % A matrix for u
spy(A4u)
MATLAB の関数は入力がスパース行列であればスパース行列用のアルゴリズムを（勝手に）使ってくれます。
r = rand(nx-1,ny); % 試しに右辺を乱数で
xu = A4u\r(:);もし時間ステップ dt とレイノルズ数 Re が変わらなければ、事前に
decomposition オブジェクトは、分解手法を入力行列のプロパティに基づいて自動的に選択してくれるので、実装方法は以下の通り一行追加するだけです。
A4u = decomposition(A4u);
xu1 = A4u\r(:);行列分解を毎回計算すると本末転倒なので、ここは初回だけ計算するように persistent 変数（永続変数）として定義するのもミソ。
最後に反復法。ここにも 11 個のアルゴリズムが用意されています。選択肢がありすぎてむしろ困ります。
ここはひとまず対称行列も正定行列も要求しない cgs（共役傾斜二乗法）で。
A4u = eye(size(L4u),'like',L4u)-dt/(2*Re)*L4u;
xu2 = cgs(A4u,r(:));cgs は、相対残差 5.4e-11 をもつ解に 反復 2 で収束しました。
前処理条件をうまく使えば高速化できそうですが、ここでも「時間ステップ dt とレイノルズ数 Re が変わらなければ」という条件が必要です。
反復法の良いところは行列を明示的に用意する必要はないところ。
xu4 = cgs(@(x) operatorAu_CNAB(x,dt,Re,nx,ny,dx,dy),r(:));cgs は、相対残差 5.4e-11 をもつ解に 反復 2 で収束しました。
みたいに
Ax = x - dt/(2*Re)*(xbig(1:end-2,2:end-1)-2*xbig(2:end-1,2:end-1)+xbig(3:end,2:end-1))/dx^2; % nx-1 * ny Ax = Ax - dt/(2*Re)*(xbig(2:end-1,1:end-2)-2*xbig(2:end-1,2:end-1)+xbig(2:end-1,3:end))/dy^2; % nx-1 * ny Ax = Ax(:);
という風に差分式のまま書き出せるので簡単です。
さて、以上の手法を使って
ここでは少し動きを出すために、 3000 ステップ回す途中で領域上部の速度を反転させています。ここまでの速度場を更新する計算処理は updateVelocityField_CNAB.m として関数化しておきます。
注意：R2019b では recordGIF = true 設定で以下を実行するとエラーが発生します。GIFを作成する場合には script_AnimateVelocityField_part2.m （以下と同じ内容です）を実行してください。
可視化設定**：**矢印は quiver 関数で描けるんですが、矢印が多いと見づらいので visRate 間隔で間引いて表示させます。また毎ステップ可視化すると GIF が重くなるので、recordRate 毎に表示更新して GIF に書き込むようにします。
visRate = 4; % downsample rate of the data for quiver
recordGIF = false; % GIF 作成する場合は true に変更
recordRate = 20;
filename = 'animation_part2.gif'; % Specify the output file name環境設定
Re = 5000; % Reynolds number
nt = 3000; % max time steps
Lx = 1; Ly = 1; % domain size
Nx = 80; Ny = 80; % Number of grids
dt = 0.003; % time step;
% Grid size (Equispaced)
dx = Lx/Nx;
dy = Ly/Ny;
% Coordinate of each grid (cell center)
xce = ((1:Nx)-0.5)*dx;
yce = ((1:Ny)-0.5)*dy;
% Coordinate of each grid (cell corner)
xco = (0:Nx)*dx;
yco = (0:Ny)*dy;速度場の初期化
figure
u = zeros(Nx+1,Ny+2); % velocity in x direction (u)
v = zeros(Nx+2,Ny+1); % velocity in y direction (v)
uce = (u(1:end-1,2:end-1)+u(2:end,2:end-1))/2; % u at cell center
vce = (v(2:end-1,1:end-1)+v(2:end-1,2:end))/2; % v at cell center等高線図
figure
[Xce,Yce] = meshgrid(xce,yce); % cell centerの座標グリッド
[~,h_abs] = contourf(Xce',Yce',sqrt(uce.^2+vce.^2)); % 等高線図警告: 等高線図は ZData が定数の場合はレンダリングされません
h_abs.LevelList = linspace(0,0.6,12);
hold on速度場（矢印）
% 表示用にデータを間引きます（d = downsampled）
xced = xce(1:visRate:end);
yced = yce(1:visRate:end);
[Xced,Yced] = meshgrid(xced, yced);
uced = uce(1:visRate:end,1:visRate:end);
vced = vce(1:visRate:end,1:visRate:end);
h_quiver = quiver(Xced',Yced',uced,vced,3,'Color',[1,1,1]);
hold off
xlim([0 Lx]); ylim([0 Ly]);おまけで領域上部の速度（境界条件）も矢印で表示しておきます。
harrow = annotation('textarrow',[0.3 0.7],[0.96 0.96],"LineWidth",2);余計なものは消しておきましょう。
haxes = gca;
haxes.XTick = [];
haxes.YTick = [];
haxes.CLim = [0,0.6];シミュレーション開始
initialFrame = true;
for ii = 1:nt
bctop = 1; % 境界上部の速度 u
if ii > 1000
bctop = -1;
harrow.X = [0.7, 0.3]; % 矢印の向きも反転
end
% 速度場更新（コサイン変換使用）
[u,v] = updateVelocityField_CNAB_bctop(u,v,Nx,Ny,dx,dy,Re,dt,bctop,'dct');
% [u,v] = updateVelocityField_RK3_bctop(u,v,Nx,Ny,dx,dy,Re,dt,bctop,'dct');
% 描画は recordRate 毎に実施
if mod(ii,recordRate) == 0
% get velocity at the cell center (for visualization)
uce = (u(1:end-1,2:end-1)+u(2:end,2:end-1))/2; % u at cell center
vce = (v(2:end-1,1:end-1)+v(2:end-1,2:end))/2; % v at cell center
% update plot (downsample)
h_quiver.UData = uce(1:visRate:end,1:visRate:end);
h_quiver.VData = vce(1:visRate:end,1:visRate:end);
h_abs.ZData = sqrt(uce.^2+vce.^2);
drawnow
if recordGIF
frame = getframe(gcf); %#ok<UNRCH> % Figure 画面をムービーフレーム（構造体）としてキャプチャ
tmp = frame2im(frame); % 画像に変更
[A,map] = rgb2ind(tmp,256); % RGB -> インデックス画像に
if initialFrame
imwrite(A,map,filename,'gif','LoopCount',Inf,'DelayTime',0.1);
initialFrame = false;
else
imwrite(A,map,filename,'gif','WriteMode','append','DelayTime',0.1);% 画像をアペンド
end
end
end
end
また精度の検証は先延ばしにしますが、高 Reynolds 数がシミュレーションできるようになった模様です。
今回はクランク・ニコルソン法とアダムス・バッシュフォース法を紹介しました。陰解法ではあるものの、高レイノルズ数ではエネルギーの散逸力が足りず時間ステップを小さくしないと発散してしまうケースがみられました。
もし試してみて「発散する！」とか「時間ステップが小さすぎてシミュレーションが進まない・・」ということであればおまけで3段階のルンゲ・クッタも実装してますのでそちらを試してください。
[u,v] = updateVelocityField_CNAB_bctop(u,v,Nx,Ny,dx,dy,Re,dt,bctop,'dct');
の部分を
[u,v] = updateVelocityField_RK3_bctop(u,v,Nx,Ny,dx,dy,Re,dt,bctop,'dct');
に変えてください。
次回「非圧縮性 Navier-Stokes 方程式の数値解法３：境界条件で遊ぶ」では、境界条件に自由度を持たせてキャビティ流れ以外もシミュレーションできるようにします。今度こそ精度検証も行います。
興味のある方いらっしゃいましたら、LGTM で応援してください :)
Spalart et al. (1991) [3] で紹介されているルンゲ・クッタ法（３段階の半陰解法）も実装してみました。updateVelocityField_RK3.m を見てください。
dt 時間ステップ進むのに３倍の計算量を使いますが、安定する分 dt を大きくできるメリットがあります。具体的な式と係数は以下の通り。
ここで
です。
[1] B. Perot, An analysis of the fractional step method. J. Comp. Physics, 108: 51-58, 1993
[2] Simens, M.P., Jim´enez, J., Hoyas, S. and Mizuno, Y. 2009 A high-resolution code for turbulent boundary layers. Journal of Computational Physics 228 (11), 4218-4231.
[3] Spalart, P.R., Moser, R.D. and Rogers, M.M. 1991 Spectral methods for the Navier-Stokes equations with one infinite and two periodic directions. Journal of Computational Physics 96 (2), 297-324.