SPH法の理論と実装 ~理論編1~ - Qiita ですっ飛ばしたところを少し詳しく書きます。自分用のメモも兼ねて。
より正確に知りたい人は、Rosswog (2009) や Springel (2010) を眺めたほうがいいです。
記号
| 記号 | 説明 |
|---|---|
| $\boldsymbol{r}$ | 位置 |
| $\boldsymbol{v}$ | 速度 |
| $m$ | 質量 |
| $p$ | 圧力 |
| $\rho$ | 密度 |
| $u$ | 単位質量当たり内部エネルギー |
| $L$ | ラグランジアン |
| $W(r,h)$ | カーネル関数 |
| $h$ | smoothing length |
| $d$ | 次元 |
導出
ラグランジアン
$$
L = \sum_j \left( m_j \frac{\boldsymbol{v}^2_j}{2} - m_j u_j \right)
$$
に対して、Euler-Lagrange 方程式
\frac{d}{dt} \left( \frac{\partial L}{\partial \dot{\boldsymbol{r}}_i} \right) - \frac{\partial L}{\partial \boldsymbol{r}_i} = 0
を考えます。第1項は簡単で、
\frac{d}{dt} \left( \frac{\partial L}{\partial \dot{\boldsymbol{r}}_i} \right) = m_i \frac{d \boldsymbol{v}_i}{dt}
です。
第2項は
\frac{\partial L}{\partial \boldsymbol{r}_i} = - \sum_j m_j \frac{\partial u_j}{\partial \boldsymbol{r}_i} = - \sum_j m_j \frac{p_j}{\rho_j^2} \frac{\partial \rho_j}{\partial \boldsymbol{r}_i}
となります。等エントロピーを仮定し、
$$
du = -p dV = -p d\left( \frac{1}{\rho} \right) = \frac{p}{\rho^2} d\rho
$$
の関係式 (熱力学第一法則)を使用しています。SPH法の密度の計算式
\rho_j = \sum_k m_k W_{jk}\\
W_{jk} \equiv W(|\boldsymbol{r}_j - \boldsymbol{r}_k|, h)
から
\frac{\partial \rho_j}{\partial \boldsymbol{r}_i} =
\begin{cases}
\sum_k m_k \nabla_i W_{ik} & (i = j)\\
m_i \nabla_i W_{ij} & (i \neq j)
\end{cases}
となり、
\frac{\partial L}{\partial \boldsymbol{r}_i} = - m_i \frac{p_i}{\rho_i^2} \sum_k m_k \nabla_i W_{ik} - \sum_{j \neq i} m_j \frac{p_j}{\rho_j^2} m_i \nabla_i W_{ij}
と変形できます。添字を付け替えると
\frac{\partial L}{\partial \boldsymbol{r}_i} = - m_i \sum_j m_j \left( \frac{p_i}{\rho_i^2} + \frac{p_j}{\rho_j^2} \right) \nabla_i W_{ij}
のように整理できます。ここで、 $ \nabla_i W_{ii} = 0 $であることを使っています。
第1項と第2項を合わせて、運動方程式
m_i \frac{d \boldsymbol{v}_i}{dt} = - m_i \sum_j m_j \left( \frac{p_i}{\rho_i^2} + \frac{p_j}{\rho_j^2} \right) \nabla_i W_{ij}
が導出できます。
可変smoothing lengthの場合の導出
smoothing length が粒子ごとに与えられる場合を考えます。密度の計算式は
\rho_j = \sum_k m_k W_{jk}(h_j)\\
W_{jk}(h_j) \equiv W(|\boldsymbol{r}_j - \boldsymbol{r}_k|, h_j)
を使います。このとき、密度の空間微分は、smoothing lengthが変数であることを考慮して
\frac{\partial \rho_j}{\partial \boldsymbol{r}_i} = \sum_k m_k \left[ \nabla_i W_{jk}(h_j) + \frac{\partial W_{jk}(h_j)}{\partial h_j} \frac{\partial h_j}{\partial \rho_j} \frac{\partial \rho_j}{\partial \boldsymbol{r}_i} \right]
となります。smoothing length の条件として
$$
\rho_j h_j^d = {\rm const.}
$$
を与えると、以下のように整理できます。
\frac{\partial \rho_j}{\partial \boldsymbol{r}_i} = f_j \times
\begin{cases}
\sum_k m_k \nabla_i W_{ik}(h_i) & (i = j)\\
m_i \nabla_i W_{ij}(h_j) & (i \neq j)
\end{cases}\\
f_j = \left(1 + \frac{h_j}{\rho_j d} \sum_k m_k \frac{\partial W_{jk}(h_j)}{\partial h_j} \right)^{-1}
これを使うと、Euler-Lagrange方程式の第2項は
\frac{\partial L}{\partial \boldsymbol{r}_i} = - m_i \sum_j m_j \left( f_i \frac{p_i}{\rho_i^2} \nabla_i W_{ij}(h_i) + f_j \frac{p_j}{\rho_j^2} \nabla_i W_{ij}(h_j) \right)
となるので、運動方程式が
m_i \frac{d \boldsymbol{v}_i}{dt} = - m_i \sum_j m_j \left( f_i \frac{p_i}{\rho_i^2} \nabla_i W_{ij}(h_i) + f_j \frac{p_j}{\rho_j^2} \nabla_i W_{ij}(h_j) \right)
に変わります。
参考文献
- Springel, V. (2010). Smoothed Particle Hydrodynamics in Astrophysics. Annual Review of Astronomy and Astrophysics, 48(1), 391–430. https://doi.org/10.1146/annurev-astro-081309-130914
- Rosswog, S. (2009). Astrophysical smooth particle hydrodynamics. New Astronomy Reviews, 53(4–6), 78–104. https://doi.org/10.1016/j.newar.2009.08.007