はじめに

冬休みの自由研究の2つ目として、1次元オイラー方程式を解いてみました。前回と同じく、C++でシミュレータを作成し、Jupyter Notebookを使って可視化しました。

前回の記事：

アルゴリズム

1次元オイラー方程式の保存型表示は次式で与えられます。

\frac{\partial \mathbf{U}}{\partial t} + \frac{\partial \mathbf{F}}{\partial x} = 0,  \tag{1}

ここで保存変数ベクトル$\mathbf{U}$と流束ベクトル$\mathbf{F}$は次式で表されます。

\mathbf{U} = 
\begin{bmatrix}
\rho \\
\rho u \\
\rho E
\end{bmatrix},
\quad
\mathbf{F} =
\begin{bmatrix}
\rho u \\
\rho u^2 + p \\
u \left( \rho E + p \right)
\end{bmatrix}
\tag{2}

ここで$\rho$は密度、$u$は速度、$p$は圧力、$E$は総エネルギーです。この式を離散化しオイラー前進法を適用すると

$$
\mathbf{U}_j^{n+1} = \mathbf{U}_j^n - \frac{\Delta t}{\Delta x} ( \hat{\mathbf{F}}^n _ {j + 1/2} - \hat{\mathbf{F}} _ {j - 1/2}^n ) \tag{3}
$$

となります。ここで$\hat{\mathbf{F}}_{j+1/2}^n$は時間ステップ$n$におけるセル境界$j+1/2$での数値流束ベクトルです。

数値流束ベクトルを計算する手法として、[1]ではリーマンソルバを使う方法と、2段階Lax-Wendroff法を使う方法が解説されています。リーマンソルバを使う場合は、まずセル格子点の値$\mathbf{U}_j^n$からセル境界の左側と右側の補間値$\mathbf{U} _ {j+1/2}^L$・$\mathbf{U} _ {j+1/2}^R$を計算し、リーマンソルバを用いてそれら補間値から数値流束$\hat{\mathbf{F}} _ {j+1/2}^n$を計算します。2段階Lax-Wendroff法を使う場合は、セル境界における保存変数の中間値$\mathbf{U} _ {j+1/2}^*$を求め、その中間値から直接数値流束を計算します。

以上をまとめると、アルゴリズムは以下の通りになります。

初期化
$n = 1 \dots N$
1. 数値流束の計算
  1. リーマンソルバの場合：$\mathbf{U} _ j^n$から$\mathbf{U} _ {j+1/2}^L$・$\mathbf{U} _ {j+1/2}^R$を計算し、リーマンソルバを使って$\hat{\mathbf{F}} _ {j+1/2}$を計算する。
  2. 2段階Lax-Wendroff法の場合：$\mathbf{U} _ {j+1/2}^*$を計算してから$\hat{\mathbf{F}} _ {j+1/2}$を計算する。
2. 時間積分：$\mathbf{U} _ j^n$・$\hat{\mathbf{F}} _ {j+1/2}$から$\mathbf{U} _ j^{n+1}$を求める。
3. 境界条件を適用する。
計算終了

空間再構築およびリーマンソルバについては、肖・長﨑（2020）¹で解説されている手法のうち、以下の手法を実装しました。

空間再構築
- 1次精度
- 2次精度：Lax-Wendroff法
- TVD法
リーマンソルバ
- Steger-Warmingリーマンソルバ
- Roeリーマンソルバ

実装

ここからはシミュレータの実装について解説します。可能な限り各離散化式をそのまま実装できるように、シミュレータでは式(2)の保存変数ベクトルと流束ベクトルを以下のように定義しています。

\mathbf{U} =
\begin{bmatrix}
U_{11} & U_{12} & U_{13} \\
U_{21} & U_{22} & U_{23} \\
 & \vdots & \\
U_{N1} & U_{N2} & U_{N3}
\end{bmatrix}
, \quad
\mathbf{F} =
\begin{bmatrix}
F_{11} & F_{12} & F_{13} \\
F_{21} & F_{22} & F_{23} \\
 & \vdots & \\
F_{N1} & F_{N2} & F_{N3}
\end{bmatrix}

ここで$U_{ij}$・$F_{ij}$は、それぞれセル$i$における保存変数ベクトル・流束ベクトルの$j$成分です。シミュレータ内ではこれらベクトルをEigen::MatrixXdを使って受け渡ししています。このように定義すると、1次元を2・3次元に拡張しても、セル番号をベクトル行番号に変換する方法を追加すれば、アルゴリズムの変更を最小限に抑えることができます。

初期条件は、基本変数ベクトル$\mathbf{V} = \begin{bmatrix} \rho & u & p \end{bmatrix}^T$を離散化したベクトル

\mathbf{V} =
\begin{bmatrix}
V_{11} & V_{12} & V_{13} \\
V_{21} & V_{22} & V_{23} \\
& \vdots  \\
V_{N1} & V_{N2} & V_{N3}
\end{bmatrix}

でシミュレータに与え、シミュレータ内で保存変数ベクトルに変換します。

数値流束を計算するクラスは、リーマンソルバを使うものと2段階Lax-Wendroff法を使うものに分けて以下の様に定義しました。

template <typename SpacialReconstructor, typename RiemannSolver>
class RiemannFluxCalculator {
 public:
  /**
   * @brief Compute flux vector
   * 
   * @param[in] U Conservation variables vector
   * @return Flux vector
   */
  template <typename Derived>
  Eigen::MatrixXd compute(const Eigen::MatrixBase<Derived>& U) const noexcept {
    using Eigen::MatrixXd;
    const MatrixXd Ul = reconstructor_.calc_left(U);
    const MatrixXd Ur = reconstructor_.calc_right(U);
    return solver_.calc_flux(Ul, Ur);
  }

 private:
  SpacialReconstructor reconstructor_;
  RiemannSolver solver_;
};

class LaxWendroffFluxCalculator {
 public:
  /**
   * @brief Compute flux vector
   * 
   * @param[in] U Conservation variables vector
   * @param[in] dt Time step length
   * @return Flux vector
   */
  template <typename Derived>
  Eigen::MatrixXd compute(const Eigen::MatrixBase<Derived>& U,
                          double dt) const noexcept {
    /* ... */
  }
};

時間積分はオイラー陽解法と2次ルンゲ・クッタ法を定義しました。

class ExplicitEulerTimeIntegrator {
 public:
  /**
   * @brief Update variables using the explicit Euler method.
   *
   * @param U[in,out] Conservation variables vector
   * @param dt[in] Time step length
   * @param flux[in] Numerical flux calculator
   * @param boundary[in] Boundary condition
   */
  template <typename Derived, typename FluxCalculator, typename Boundary>
  void update(Eigen::MatrixBase<Derived>& U, double dt,
              const FluxCalculator& flux,
              const Boundary& boundary) const noexcept {
    using Eigen::all, Eigen::seqN, Eigen::MatrixXd;
    const auto i = n_boundary_cells_;
    const auto n = n_domain_cells_;
    // Two-step Lax-Wendroff method
    if constexpr (std::is_same_v<FluxCalculator, LaxWendroffFluxCalculator>) {
      const MatrixXd F = flux.compute(U, dt);
      U(seqN(i, n), all) -= (dt / dx_) * (F.bottomRows(n) - F.topRows(n));
    // Riemann-solver-based flux calculators
    } else {
      const MatrixXd F = flux.compute(U);
      U(seqN(i, n), all) -= (dt / dx_) * (F.bottomRows(n) - F.topRows(n));
    }
    boundary.apply(U);
  }

 private:
  double dx_;             ///> Grid length
  int n_boundary_cells_;  ///> Number of boundary cells
  int n_domain_cells_;    ///> Number of domain cells
};

class RungeKutta2ndOrderTimeIntegrator {
 public:
  /**
   * @brief Update conservation variables.
   *
   * @param U[in,out] Conservation variables vector
   * @param dt[in] Time step length
   * @param flux[in] Numerical flux calculator
   * @param boundary[in] Boundary condition
   */
  template <typename Derived, typename FluxCalculator, typename Boundary>
  void update(Eigen::MatrixBase<Derived>& U, double dt,
              const FluxCalculator& flux,
              const Boundary& boundary) const noexcept {
    static_assert(!std::is_same_v<FluxCalculator, LaxWendroffFluxCalculator>,
                  "LaxWendroffFluxCalculator cannot be used.");

    using Eigen::all, Eigen::seqN, Eigen::MatrixXd;

    const auto i = n_boundary_cells_;
    const auto n = n_domain_cells_;
    const auto ntotal = 2 * n_boundary_cells_ + n_domain_cells_;
    const auto rng = seqN(i, n);

    // First step
    const MatrixXd F1 = flux.compute(U);
    const MatrixXd dU1 = (dt / dx_) * (F1.topRows(n) - F1.bottomRows(n));
    MatrixXd U1 = U(rng, all) + dU1;
    boundary.apply(U1);

    // Second step
    const MatrixXd F2 = flux.compute(U1);
    const MatrixXd dU2 = (dt / dx_) * (F2.topRows(n) - F2.bottomRows(n));
    U(rng, all) += 0.5 * (dU1 + dU2);
    boundary.apply(U);
  }

 private:
  double dx_;             ///> Grid length
  int n_boundary_cells_;  ///> Number of boundary cells
  int n_domain_cells_;    ///> Number of domain cells
};

これらを用いるとシミュレータのメイン関数は以下の様に実装できます。

template <typename FluxCalculator, typename TimeIntegrator>
class EulerEquationSimulator1d {
 public:
  /**
   * @brief Run a simulation case.
   *
   * @tparam Derived
   * @param V Primitive variables vector at the initial condition.
   * @return Primitive Variables vector at the end of time steps.
   */
  template <typename Derived>
  Eigen::MatrixXd run(const Eigen::MatrixBase<Derived>& V) const noexcept {
    using Eigen::MatrixXd;
    MatrixXd U = this->to_conservation_vars(V);
    boundary_.apply(U);
    double t = 0.0;
    int tsteps = 0;

    while (t < tend_) {
      auto dt = timestep_.compute(U);
      if (t + dt > tend_) {
        dt = tend_ - t;
      }
      tsteps += 1;
      t += dt;
      integrator_.update(U, dt, flux_, boundary_);
    }

    return this->to_primitive_vars(U);
  }

 private:
  /* ... */
  NoFlowBoundary boundary_;
  FluxCalculator flux_;
  TimeIntegrator integrator_;
  TimestepLengthCalculator timestep_;
};

結果

肖・長﨑（2020）の5.1.3節にある1次元衝撃波管問題を解き、ほぼ同じ結果が得られました。

環境について

このシミュレータを作成した環境は以下の通りです。

OS: Windows 11 Pro
コンパイラ: MSVC 2022
ビルドシステム: CMake (Visual Studio付属)
エディタ: Visual Studio Code
Python環境: miniconda

Windows上で作成していますが、CMakeとPythonを使っているので、Linux/Mac上でもコンパイルして動かせるはずです。問題が起きたら連絡してください。

肖鋒・長﨑孝夫　2020　数値流体解析の基礎－Visual C++とgnuplotによる圧縮性・非圧縮性流体解析－　コロナ社 ↩