有限体積法の基礎：リーマンソルバ【MATLABコード付き】

はじめに

　前回は1次元非粘性Bugers方程式について有限体積法で離散化し、以下の物理量\(\bar{q}_{i}\)に関する時間発展式を得ました。

しかし、数値流束\(\tilde{f}^n_{i+\frac{1}{2}},\tilde{f}^n_{i-\frac{1}{2}}\)の計算に必要なセル界面での物理量\(q_{i+\frac{1}{2}},q_{i-\frac{1}{2}}\)の具体的な中身については述べていませんでした。そこで今回は\(\tilde{f}^n_{i+\frac{1}{2}},\tilde{f}^n_{i-\frac{1}{2}}\)と離散変数\(\bar{q}_{i}\)を結びつける関係式「リーマンソルバ」についてまとめたいと思います。

有限体積法の基礎：差分法との比較

はじめに 　流体の方程式の離散化には差分法や有限体積法が用いられます。特に有限体積法は多くのソルバへ実装されており、シミュレーションを行う上でその理論的背景を知っておくことは重要です。そこで今回から有限体積法の基礎についてまとめていきます。...

taku-naka.site

2023.05.27

セル界面での物理量の計算

　式\((2),(3)\)の数値流束の計算のため、まずは界面物理量\(q_{i+\frac{1}{2}},q_{i-\frac{1}{2}}\)の計算方法から考えます。

差分法では格子点の値を用いた簡単な補間からセル界面の物理量\(q_{i+\frac{1}{2}},q_{i-\frac{1}{2}}\)を定義しました（過去の記事参照）。一方、有限体積法ではセル内での物理量の分布を考えることでこれを定義します。

以下の通り、有限体積法では物理量の空間分布\(q(x)\)を考え、そのセル内空間平均\(\bar{\cal{q}}_{i}\)を離散変数としました。

空間分布\(q(x)\)の形が決まっていれば、当然セル界面での物理量\(q_{i+\frac{1}{2}}=q(x_{i+\frac{1}{2}})\)を求めることができます。しかし、\(q(x)\)は連続量であり本来知りたい解そのものなので、\(x\)の関数としてその中身を具体的に決めることはできません。

そこで以下の図のように\(q(x)\)をセルごとに定義された無数の補間関数\(Q_i(x)\)で近似することでその範囲を分割し、なおかつ\(Q_i(x)\)を離散変数\(\bar{\cal{q}}_{i}\)の関数として表現することで、補間関数の形そのものを方程式系に組み込みます。

式\((4)\)はもともと単一セル内部の積分ですので補間関数\(Q_i(x)\)を用いて以下のように書いても問題ないはずです。

このような補間関数\(Q_i(x)\)が決まれば、界面の物理量\(q_{i+\frac{1}{2}}=Q_i(x_{i+\frac{1}{2}})\)が求まるというわけです。

ただし一点ややこしい点があります。補間関数\(Q_i(x)\)はセルごとに定義されるため、界面での物理量は基本的に左右で異なる値を持ってしまう点です。例えばセル界面\(i+\frac{1}{2}\)の物理量は\({}^R\! q_{i+\frac{1}{2}},{}^L\! q_{i+\frac{1}{2}}\)の二種類で定義されます。

数値流束を計算する際にこの２種類のどちらをどのように用いるのか疑問に思えてきますが、この辺りを決めるのが後に出てくるリーマンソルバです。

では\({}^R\! q_{i-\frac{1}{2}}\)を例に実際に界面物理量を求めてみましょう。まずは補間関数\(Q_i(x)\)を離散変数\(\bar{\cal{q}}_{i}\)を用いて表現します。\(Q_i(x)\)を左側界面\(x_{i-\frac{1}{2}}\)においてテイラー展開すると以下の通りです。

両辺をセル\(i\)内で積分し\(\Delta x_{i}\)で割ることで体積平均\(\bar{\cal{q}}_{i}\)を用いて表現します。

ここで積分は\(A=x-x_{i-\frac{1}{2}}\)とおき置換積分で計算します。また\(\Delta x_{i}=x_{i+\frac{1}{2}}-x_{i-\frac{1}{2}}\)と置くことにします。

数値流束の計算で必要なのは界面での物理量\({}^R\! q_{i-\frac{1}{2}}\)なので、これについてまとめます。

以上でセル界面\(x_{i-\frac{1}{2}}\)の右側の物理量\({}^R\! q_{i-\frac{1}{2}}\)を求めることが出来ました。同様の式展開は同じセル\(i\)の\(x_{i+\frac{1}{2}}\)側の界面である\({}^L\! q_{i+\frac{1}{2}}\)でも成立します。

また\({}^R\! q_{i+\frac{1}{2}},{}^L\! q_{i-\frac{1}{2}}\)も同様に考えて以下のように書けます。

ここで補間関数の微分項\(Q’,Q^{\prime\prime}\)が未知のままですが、ひとまず今回は空間１次精度で考えることとし高次の項は無視してしまうことにします。（この部分を考慮すると高次精度の数値流束を計算できますがそれはまた今度…。）

以上よりセル界面\(x_{i-\frac{1}{2}},x_{i+\frac{1}{2}}\)における左右の物理量\({}^L\! q_{i-\frac{1}{2}},{}^R\! q_{i-\frac{1}{2}},{}^L\! q_{i+\frac{1}{2}},{}^R\! q_{i+\frac{1}{2}}\)が求まりました。1次精度で考える限り、界面での物理量は属するセルの離散値そのものの値となることが分かります。

Roeリーマンソルバ

　界面における物理量が得られたので、次はこれらを用いて数値流束\(\tilde{f}_{i+\frac{1}{2}},\tilde{f}_{i-\frac{1}{2}}\)を計算します。1次元Bugers方程式における流束は以下のような式で与えられました。（過去の記事参照）

また、上述の通り有限体積法において界面物理量\(q_{i+\frac{1}{2}}\)はセルごとに異なる値を取り\({}^R\! q_{i+\frac{1}{2}},{}^L\! q_{i+\frac{1}{2}}\)の2種類で表現されました。

ではこれらを用いてどのように数値流束を表現すればよいのでしょうか？

このような数値流束\(\tilde{f}_{i+\frac{1}{2}}\)と\({}^R\! q_{i+\frac{1}{2}},{}^L\! q_{i+\frac{1}{2}}\)の関係を与える式を「リーマンソルバ」と呼び、代表的なものとして以下のRoe（ロー）リーマンソルバが挙げられます。

このままだとこの式の意味するところは分かりずらいので少し変形します。式\((13)\)を用いると式\((15)\)は以下のように書けます。

上式より\(\alpha_{i+\frac{1}{2}}\)はセル界面\({i+\frac{1}{2}}\)における左右の物理量の平均値となっていることが分かります。さらに式\((12)\)を代入すると以下のように書けます。

上式より空間１次精度においては\(\alpha_{i+\frac{1}{2}}\)はその両隣セルの離散変数の平均値となることが分かります。例えば\(q\)が流速だとすると\(\alpha\)は左右セルの平均流速を意味することとなります。

上式を用いると式\((14)\)は次のように書けます。

実際に場合分けしてみて考えます。まず\(\alpha_{i+\frac{1}{2}}=\frac{1}{2}\left( \bar{q}_{i+1}+\bar{q}_{i}\right)>0\)の場合です。

\(\alpha_{i+\frac{1}{2}}=\frac{1}{2}\left( \bar{q}_{i+1}+\bar{q}_{i}\right)<0\)の場合も同様に書けます。

以上をまとめると以下の通りです。

\begin{equation}
\tilde{f}_{i+\frac{1}{2}}=
\begin{cases}
f(\bar{q}_{i})=\frac{1}{2}\bar{q}_{i}^2 & \alpha_{i+\frac{1}{2}}=\frac{1}{2}\left( \bar{q}_{i+1}+\bar{q}_{i}\right)>0の場合\\
f(\bar{q}_{i+1})=\frac{1}{2}\bar{q}_{i+1}^2 & \alpha_{i+\frac{1}{2}}=\frac{1}{2}\left( \bar{q}_{i+1}+\bar{q}_{i}\right)<0の場合 \\
\frac{1}{2}\left(f(\bar{q}_{i})+f(\bar{q}_{i+1})\right) =\frac{1}{4}\left(\bar{q}_{i}^2+\bar{q}_{i+1}^2\right)& \alpha_{i+\frac{1}{2}}=\frac{1}{2}\left( \bar{q}_{i+1}+\bar{q}_{i}\right)=0の場合 \\
\end{cases}
\tag{21}
\end{equation}

上式は「注目するセル界面の両側セルの物理量の平均値から移流方向を定め、その風上側セルの離散変数を数値流束の計算に採用する」ということを意味しています。これは過去記事に記載した風上差分法による離散化と同じです。

言い換えると有限体積法における空間１次精度のRoeリーマンソルバは差分法における１次精度風上差分法と同等ということが出来ます。

\(\tilde{f}_{i+\frac{1}{2}}\)も同様に考えて

\begin{equation}
\tilde{f}_{i-\frac{1}{2}}=
\begin{cases}
f(\bar{q}_{i-1})=\frac{1}{2}\bar{q}_{i-1}^2 & \alpha_{i+\frac{1}{2}}=\frac{1}{2}\left( \bar{q}_{i}+\bar{q}_{i-1}\right)>0の場合\\
f(\bar{q}_{i})=\frac{1}{2}\bar{q}_{i}^2 & \alpha_{i+\frac{1}{2}}=\frac{1}{2}\left( \bar{q}_{i}+\bar{q}_{i-1}\right)<0の場合 \\
\frac{1}{2}\left(f(\bar{q}_{i-1})+f(\bar{q}_{i})\right) =\frac{1}{4}\left(\bar{q}_{i-1}^2+\bar{q}_{i}^2\right)& \alpha_{i+\frac{1}{2}}=\frac{1}{2}\left( \bar{q}_{i}+\bar{q}_{i-1}\right)=0の場合 \\
\end{cases}
\tag{22}
\end{equation}

以上を式\((1)\)に代入することで、\(\bar{q}_{i}\)の時間発展を計算することが出来ます。なおここでは分かりやすさのため場合分けして書きましたが、実際に実装する際は式\((18)\)の形で一本の式としてまとめるのが良いでしょう。

MATLABコード

　以下に1次元非粘性Bugers方程式を空間1次精度のRoeリーマンソルバで求解するMATLABコードを示します。githubからも見ることが出来ます。

% 初期化
clearvars;

% パラメーター
i_max = 100;    % 格子セル数
XL = -1.0;      % 計算領域左端の座標
XR = 1.0;       % 計算領域右端の座標
a = 1.0;        % 線形移流方程式の移流速度
tstop = 0.5;    % 計算停止時刻


% 配列定義
i = 0;                      % セル番号； (i番目セルの左境界の番号はi、右辺境界の番号はi+1とする)
x = zeros(i_max + 1, 1);    % セル境界の座標
u = zeros(i_max, 1);        % セル平均値 （数値解）
ue = zeros(i_max, 1);       % セル平均値 （厳密解）
ul = zeros(i_max + 1, 1);   % セル境界左側の変数値
ur = zeros(i_max + 1, 1);   % セル境界右側の変数値
f = zeros(i_max + 1, 1);    % セル境界の流束
n = 0;                      % 時間ステップ
t = 0;                      % 計算時間


% 線形移流方程式の初期値を選ぶ
sw1 = 1; %"滑らかな分布（正弦波）：１、"不連続な分布（矩形波）：２

dx = (XR - XL) / (i_max - 4.0);% 格子間隔 計算領域外に二つずつ余分なセルを準備。周期的境界条件向けの設定
dt = 0.2 * dx; % 時間刻み

x(1) = XL - 2.0 * dx;% 計算領域外の２セルも考慮した座標を振る。

[x, u] = initc(i_max, x, dx, sw1, u);% 計算格子，時間刻み，初期条件を設定する

ue = exact(sw1, i_max, ue, x, t, dx);%　厳密解を求める

% メインループ
while t <= tstop

    n = n + 1;
    t = t + dt;

    [ul, ur] = reconstruction_pc(i_max, u, ul, ur);% 空間再構築

    f = riemann_roe(i_max, f, ul, ur);  % リーマンソルバー

    u = update(i_max, u, dt, dx, f); % 時間積分

    u = bc(i_max, u); % 境界条件

    ue = exact(sw1, i_max, ue, x, t, dx);%　厳密解を求める

    fprintf("n=%d, t=%f \n", n, t);

    % 動画保存
    if n == 1
        plotconfig(x(1 : end - 1), ue, u,  t)
        v = VideoWriter('u.mp4','MPEG-4');
        v.FrameRate = 40;
        open(v);
    else
        plotconfig(x(1 : end - 1), ue, u,  t)
        frame = getframe(gcf);
        writeVideo(v,frame);
    end

end

% 動画ファイルを閉じる
close(v);

%% 以下ローカル関数

function[f] = f_flux(x)

f = 0.5*x*x;

end


function [x, u] = initc(i_max, x, dx, sw1, u)

for  i = 2 : i_max + 1
    x(i) = x(i - 1) + dx;% 格子点の座標
end

switch  sw1

    case 1  % 初期の変数値（滑らかな分布）

        for  i = 1 : i_max
            u(i) = 0.5*(1.1 + sin(2.* pi *(x(i)-x(3))));% 三番目の要素（XL、XRの座標値）を基準に考える。
        end

    case 2 % 初期の変数値（不連続な分布）

        for  i = 1 : i_max
            u(i) = 0.1;
        end
        for i = i_max / 2 - 10 : i_max / 2 + 10
            u(i) = 1.0;
        end

    otherwise

        disp("初期値を正しく選択されていません。");

end

end

function [ue] = exact(sw1, i_max, ue, x, t, dx)

switch sw1

    case 1 % 初期の変数値（滑らかな分布）

        for i = 1 : i_max

            c = 2 * pi;
            f = ue(i) - 0.5 * (1.1 + sin(c * (x(i)-x(3) - ue(i) * t)));
            df = 1.0 + 0.5 * c * cos(c *(x(i) - x(3) - ue(i) * t)) * t;
            count = 0;
            while abs(f) >= 1.0e-6
                count = count + 1;
                ue(i) = ue(i) - f / df; % ニュートン法の計算式*/
                f = ue(i) - 0.5 * (1.1 + sin(c*((x(i) - x(3) - ue(i) * t))));
                df = 1.0 + 0.5 * c * cos(c * ((x(i) - x(3) - ue(i) * t))) * t;
                if count > 10000
                    disp("厳密解が収束しないので反復を打ち切ります。");
                    break
                end
            end

        end


    case 2 % 初期の変数値（不連続な分布）

        xc = - dx * 10 + t;
        xl = - dx * 10 + 0.1 * t;
        xr = dx * 10 + 0.55 * t;

        % 周期境界条件　
        if  xl > 1.0
            xl = -2.0 + xl;
        end

        % 周期境界条件　
        if xr > 1.0
            xr = -2.0 + xr;
        end

        % 周期境界条件　
        if xc > 1.0
            xc = -2.0 + xc;
        end

        for i = 1 : i_max
            if x(i) <= xl
                ue(i) = 0.1;
            end
            if x(i) >= xl && x(i) <= xc
                ue(i) =(x(i) - xl) / (xc - xl) * 0.9 + 0.1;
            end
            if x(i) >= xc && x(i) <= xr
                ue(i) = 1.0;
            end
            if x(i) >= xr
                ue(i) = 0.1;
            end
        end

end

end

function [ul, ur] = reconstruction_pc(i_max, u, ul, ur)

for i = 2 : i_max - 2
    ul(i + 1) = u(i); % セル境界(i+1/2)左側の値
    ur(i + 1) = u(i + 1); % セル境界(i+1/2)右側の値
end

end

function [f] = riemann_roe(i_max, f, ul, ur) % 流束を計算する

for i = 3 : i_max - 1
    alpha_12 = 0.5 * (ur(i) + ul(i));
    f(i) = 1.0 / 2.0 * (f_flux(ul(i)) + f_flux(ur(i))) - 1.0 / 2.0 * abs(alpha_12) * (ur(i) - ul(i));
end

end

function [u] = update(i_max, u, dt, dx, f)

for  i = 3 : i_max - 2
    u(i) = u(i) - dt / dx * (f(i + 1) - f(i));% 計算変数を更新する
end

end

function [u] = bc(i_max, u) % 周期境界条件

u(1) = u(i_max - 3); % 計算領域左端の境界条件
u(2) = u(i_max - 2); % 計算領域左端の境界条件
u(i_max - 1) = u(3); % 計算領域右端の境界条件
u(i_max) = u(4); % 計算領域右端の境界条件

end

function [] = plotconfig(x, ue, u, t)

plot(x, ue, x, u)
%plot(x, u)

title(['time = ', num2str(t, '%.3f')]);
set(gca, 'FontName', 'Times New Roman', 'FontSize', 12);
axis equal;
axis tight;
axis on;
fig=gcf;
fig.Color='white';
ylim([0 1.1]);
xlabel('x')
ylabel('u')

% 凡例
legend({'exact', '1st order Roe'},'Location','southwest','FontSize', 10)

% 新しいプロットの時、軸設定を保持したまま前のグラフィックスオブジェクトを消去
set(gca,'nextplot','replacechildren');

end

実行結果

　以下に実行結果を示します。Roeスキームで得られた解は1次風上と同等なので拡散成分が強いですが、厳密解を再現していることが分かります。

おわりに

　今回は有限体積法により１次元非粘性Bugers方程式を解きました。リーマンソルバを用いて界面ベースで離散化していくのが有限体積法特有の考え方なのでしょう。次は高次精度解法についてまとめていく予定です。

参考文献

[1] 肖 鋒, “数値流体解析の基礎- Visual C++とgnuplotによる圧縮性・非圧縮性流体解析 -“, コロナ社, 2020