定式化 - YuriOku/1D_SPH GitHub Wiki

SPH法で時間発展方程式を定式化するにはいくつかの方法があります。 1D_SPHには以下の2通りの方程式系が実装されています。これらの方法はRosswog (2009)で詳しく導出、解説されています。興味のある方はそちらも参照ください。

Vanilla ice

SPH粒子の時間発展方程式を求める標準的な方法は、Lagrange形式で表現された流体方程式系をカーネル近似によって離散化する方法です。 SPH法ではSPH粒子の運動が流体の流れを表現するため、質量保存則が自動的に満たされ、連続の式は不要です。運動方程式とエネルギー式は

$$ \begin{align*} \frac{dv}{dt} &= -\frac{\nabla P}{\rho}\ \frac{du}{dt} &= -\frac{P}{\rho}\nabla \cdot v \end{align*} $$

です。ここで、 $v$ は速度、 $P$ は圧力、 $\rho$ は密度、 $u$ は単位質量当たり内部エネルギーです。空間微分をカーネル近似を用いて計算します。カーネル近似の下では、物理量 $f$ は

$$ f(x) = \int f(x') W(|x - x'|, h) dx' \cdot\cdot\cdot\cdot\cdot(1) $$

と表されます。ここで$W$ は粒子の広がりを表すカーネル関数です。式(1)を各SPH粒子の体積要素 $dV$ を用いて離散化すると

$$ f_i = \sum_j f_j dV_j W(| x_i - x_j|, h_i) \cdot\cdot\cdot\cdot\cdot(2) $$

となります。体積要素はSPH粒子が持つ物理量 $Z$ とそれに対応する連続量 $y$ から $dV = \frac{Z}{y}$ と求められ、標準的には $Z = m,,y = \rho$ が使われます。詳しくは体積要素のページをご覧ください。 $f$の空間微分は

$$ \nabla_i f_i = \sum_j f_j dV_j \nabla_i W(| x_i - x_j|, h_i) \cdot\cdot\cdot\cdot\cdot(3) $$

と表されます。 $\nabla_i W$ の計算については、カーネル補間法をご覧ください。運動方程式の右辺に式(3)を適用する前に、式が i, j の交換に対して反対称となるように、

$$ \begin{align*} \frac{\nabla P}{\rho} = \frac{P}{\rho^\lambda}\nabla\left( \frac{1}{\rho^{1-\lambda}} \right) + \frac{1}{\rho^{2 - \lambda}}\nabla \left(\frac{P}{\rho^{\lambda - 1}} \right) \end{align*} $$

の関係 (Monaghan 1992)を使います。ここで $\lambda$ は定数です。伝統的には $P/\rho$ の微分から自然に得られる$\lambda = 2$の場合の形が使われていますが、 1D_SPH では $\lambda = 1$ を採用しています。 $\lambda = 1$ とすると、異なる質量のSPH粒子があるときに発生するノイズを抑えることができます (Ritchie & Thomas 2001)。 $\lambda = 1$ とする形式はGASOLINE2 (Wadsley et al. 2017)、MAGMA2 (Rosswog 2020) で採用されています。このときSPH方程式系は

$$ \begin{align*} \frac{d\vec{v_i}}{dt} &= -\frac{1}{m_i}\sum_j Z_i Z_j \left( \frac{P_i + P_j}{y_i y_j}\right) \nabla_i \tilde{W}{ij}\ \frac{du_i}{dt} &= \frac{P_i}{m_i} \sum_j \frac{Z_i}{y_i}\frac{Z_j}{y_j}\vec{v{ij}}\cdot\nabla_i \tilde{W}_{ij} \end{align*} $$

となります。ここで、 $\vec{v_{ij}}=\vec{v_i} - \vec{v_j}$ , $\tilde{W_{ij}} = [W(r_{ij}, h_i)+W(r_{ij}, h_j)]/2$ です。

Lagrangian

SPH方程式を導出するもう一つの方法はLagrangianから出発して変分原理を用いる方法です。 SPH粒子系のLagrangianは

$$ L = \sum_i m_i(\frac{v_i^2}{2} - u_i) \cdot\cdot\cdot\cdot\cdot(4) $$

です。これをEuler-Lagrange方程式

$$ \frac{d}{dt}\left(\frac{\partial L}{\partial \vec{v}_i} \right) - \frac{\partial L}{\partial \vec{r}_i} = 0 \cdot\cdot\cdot\cdot\cdot(5) $$

に代入します。式(5)の第一項は運動方程式の左辺になります。第二項にLagrangianを代入すると

$$ \frac{\partial L}{\partial \vec{r}_i} = -\sum_j m_j\frac{\partial u_j}{\partial \rho_j} \frac{\partial \rho_j}{\partial \vec{r}_i} \cdot\cdot\cdot\cdot\cdot(6) $$

となります。ここで等エントロピー流での熱力学第一法則

$$ du = - PdV = \frac{P}{\rho^2}d\rho \cdot\cdot\cdot\cdot\cdot(7) $$

を使うとEuler-Lagrange方程式は

$$ m_i\frac{d\vec{v}_i}{dt} = -\sum_j m_j\frac{P_j}{\rho_j^2} \frac{\partial \rho_j}{\partial \vec{r}_i} \cdot\cdot\cdot\cdot\cdot(8) $$

となります。 SPH方程式系を密度から定式化する場合 ($Z = m,,y = \rho$の場合)、式(8)に

$$ \rho_j = \sum_k m_k W(r_{jk}, h_j) \cdot\cdot\cdot\cdot\cdot(9) $$

を代入します。すると、式(8)の密度の微分は

$$ \frac{\partial \rho_j}{\partial \vec{r}i} = \sum_k m_k \left(\nabla_i W(r{jk}, h_j) +\frac{\partial W(r_{jk}, h_j)}{\partial h_j} \frac{\partial h_j}{\partial \rho_k}\frac{\partial \rho_k}{\partial \vec{r}i} \right) = f_j \sum_k m_k \nabla_i W(r{jk}, h_j) \cdot\cdot\cdot\cdot\cdot(10) $$

となります。ここで、

$$ f_i = \left(1 - \frac{\partial h_i}{\partial \rho_i} \sum_j m_j \frac{\partial}{\partial h_i} W(r_{ij}, h_i) \right)^{-1} $$

はスムージング長$h$の変化を考慮に入れるための項で "grad-h-term"と呼ばれます。添え字に注意して変形すると、運動方程式

$$ \frac{d \vec{v}i}{dt} = - \sum_j m_j \left( f_i \frac{P_i}{\rho_i^2} \nabla_i W(r{ij}, h_i) + f_j \frac{P_j}{\rho_j^2} \nabla_i W(r_{ij}, h_j) \right) $$

が得られます。エネルギーの時間微分は、熱力学第一法則から

$$ \frac{du_i}{dt} = \frac{P_i}{\rho_i^2}\frac{d \rho_i}{dt} = f_i \frac{P_i}{\rho_i^2} \sum_j m_j \vec{v_{ij}} \cdot \nabla_i W(r_{ij}, h_i) $$

となります。

一般の体積要素の場合では、Euler-Lagrange 方程式(8)に

$$ y_j = \sum_k Z_k W(r_{jk}, h_j) \cdot\cdot\cdot\cdot\cdot(11) $$

を代入して、

$$ \begin{align*} \frac{d \vec{v}i}{dt} &= - \frac{1}{m_i}\sum_j Z_i Z_j \left( f_i \frac{P_i}{y_i^2} \nabla_i W(r{ij}, h_i) + f_j \frac{P_j}{y_j^2} \nabla_i W(r_{ij}, h_j) \right)\ \frac{du_i}{dt} &= f_i \frac{P_i}{m_i y_i^2} \sum_j Z_i Z_j \vec{v_{ij}} \cdot \nabla_i W(r_{ij}, h_i) \end{align*} $$

が得られます。

参考文献

Monaghan, J. J., “Smoothed particle hydrodynamics.”, Annual Review of Astronomy and Astrophysics, vol. 30, pp. 543–574, 1992. doi:10.1146/annurev.aa.30.090192.002551.
Ritchie, B. W. and Thomas, P. A., “Multiphase smoothed-particle hydrodynamics”, Monthly Notices of the Royal Astronomical Society, vol. 323, no. 3, pp. 743–756, 2001. doi:10.1046/j.1365-8711.2001.04268.x.
Rosswog, S., ``Astrophysical smooth particle hydrodynamics'', New Astronomy Reviews, vol. 53, no. 4–6, pp. 78–104, 2009. doi:10.1016/j.newar.2009.08.007.
Rosswog, S., “The Lagrangian hydrodynamics code MAGMA2”, Monthly Notices of the Royal Astronomical Society, vol. 498, no. 3, pp. 4230–4255, 2020. doi:10.1093/mnras/staa2591.
Wadsley, J. W., Keller, B. W., and Quinn, T. R., “Gasoline2: a modern smoothed particle hydrodynamics code”, Monthly Notices of the Royal Astronomical Society, vol. 471, no. 2, pp. 2357–2369, 2017. doi:10.1093/mnras/stx1643.