Key Notes on Viscous Hyperbolic SPH Equation of Motion of Hydrodynamic Model

Page last edited 4,331 days ago
From LaTeX Notes Series
Jump to navigation Jump to search

Key Notes on Viscous Hyperbolic SPH Equation of Motion of Hydrodynamic Model

考虑体积粘滞系数和剪切粘滞系数在双曲线坐标系中SPH粒子自由度的流体的相对论流体力学方程

本文档除了包括推导,疑惑外,做读书重点的记录

双曲线坐标系及其度规

这里首先列出在一般的曲线坐标下,协变导数的定义和一些基本的关系

<math>\begin{align}

&\Gamma^i_{k\ell}=\frac{1}{2}g^{im} \left(\frac{\partial g_{mk}}{\partial x^\ell} + \frac{\partial g_{m\ell}}{\partial x^k} - \frac{\partial g_{k\ell}}{\partial x^m} \right) = {1 \over 2} g^{im} (g_{mk,\ell} + g_{m\ell,k} - g_{k\ell,m}) \\ &\Gamma^i_{ki}=\frac{1}{2}g^{im}\frac{\partial g_{im}}{\partial x^k}={1 \over 2g}\frac{\partial g}{\partial x^k}=\frac{\partial \log \sqrt{|g|}}{\partial x^k}\\ &D_{\mu}\phi\equiv(\phi)_{;\mu}=\partial_{\mu}\phi\\ &D_{\mu}V_{\nu}\equiv(V_{\nu})_{;\mu}=\partial_{\mu}V_{\nu}-\Gamma^{\alpha}_{\mu\nu}V_{\alpha}\\ &D_{\mu}V^{\nu}\equiv(V^{\nu})_{;\mu}=\partial_{\mu}V^{\nu}+\Gamma^{\nu}_{\mu\alpha}V^{\alpha}\\ &D_{\mu}V_{\alpha\beta}\equiv(V_{\alpha\beta})_{;\mu}=\partial_{\mu}V_{\alpha\beta}-\Gamma^{\lambda}_{\alpha\mu}V_{\lambda\beta}-\Gamma^{\lambda}_{\beta\mu}V_{\lambda\alpha}\\ &D_{\mu}V^{\mu}\equiv(V^{\mu})_{;\mu}=\frac{1}{\sqrt{-g}}\partial_{\mu}(\sqrt{-g}V^{\mu})\\ &D_{\mu}V^{\mu\nu}\equiv(V^{\mu\nu})_{;\mu}=\frac{1}{\sqrt{-g}}\partial_{\mu}(\sqrt{-g}V^{\mu\nu})+\Gamma^{\nu}_{\alpha\mu}V^{\alpha\mu} \end{align}</math>

下面的具体讨论都将在双曲线坐标系中展开

<math>\begin{align}

&\tau=\sqrt{t^2-z^2}\\ &\eta=\frac{1}{2}\ln\left(\frac{t+z}{t-z}\right) \end{align}</math>

该坐标系下度规的具体形式为

<math>\begin{align}

g_{\mu\nu}=\begin{pmatrix} 1 & 0 & 0 & 0\\ 0 & -1 & 0 & 0\\ 0 & 0 & -1 & 0\\ 0 & 0 & 0 & -{\tau^2} \end{pmatrix} \end{align}</math>

<math>\begin{align}

\sqrt{-g}=\tau \end{align}</math>

显然度规的逆为

<math>\begin{align}

g^{\mu\nu}=\begin{pmatrix} 1 & 0 & 0 & 0\\ 0 & -1 & 0 & 0\\ 0 & 0 & -1 & 0\\ 0 & 0 & 0 & -\frac{1}{\tau^2} \end{pmatrix} \end{align}</math>

利用上述度规的具体形式,我们得到,克里斯多夫符号

<math>\begin{align}

\Gamma^i_{k\ell}=\frac{1}{2}g^{im} \left(\frac{\partial g_{mk}}{\partial x^\ell} + \frac{\partial g_{m\ell}}{\partial x^k} - \frac{\partial g_{k\ell}}{\partial x^m} \right) \end{align}</math>

仅仅三项不为零,它们是

<math>\begin{align}

\Gamma^0_{33}=\tau, \Gamma^3_{30}=\Gamma^3_{03}=\frac{1}{\tau} \end{align}</math>

下面我们开始讨论两阶的相对论粘滞流体运动方程在双曲线坐标系中的形式.

SPH自由度

实际应用中,我们需要把方程表达为SPH自由度.按SPH的物理涵义,对任何强度物理量 <math>A</math> 在空间任何点 <math>j</math> 实验室坐标系的值由该物理量在SPH粒子处 <math>i </math>实验室坐标系的取值(依赖于状态方程)以及守恒荷 <math>\nu</math> (对应密度 <math>\sigma</math> )确定下来

<math>\begin{align}

A_j \equiv\sum_i \left(\frac{A}{\sigma^*}\right)_i\nu_iW_{ji} =\sum_i \left(\frac{A}{\sigma\tau\gamma}\right)_i\nu_iW_{ji} \end{align}</math>

注意到一个重要的约定,如果 <math>A</math> 本身也是一个守恒荷 <math>B</math> 在实验室坐标系中的密度 <math>b^*</math> ,那么由于

<math>\begin{align}

&A=b^* \\ &\left(\frac{A}{\sigma^*}\right)_i=\left(\frac{b^*}{\sigma^*}\right)_i=\left(\frac{b}{\sigma}\right)_i \\ &A_j=\sum_i \left(\frac{A}{\sigma^*}\right)_i\nu_iW_{ji}=\sum_i b_iW_{ji} \end{align}</math>

而对于一般的物理量,如压强,则没有上面的简化形式,运动方程中涉及对物理量的空间偏导,从而也必须表达为SPH自由度

<math>\begin{align}

(\partial A)_j\equiv\sum_i \left(\frac{A}{\sigma^*}\right)_i\nu_i\nabla_jW_{ji} \end{align}</math>

对于本模型,实际上只有两项需要写为SPH自由度,其余都仅仅与体系状态方程有关.第一项与压强有关,是 <math>\partial_i {\tilde P}</math> .这里为了得到所需的对称形式,我们必须引入一个貌似多余步骤.这实际上从一个侧面说明了SPH方法不是第一性原理,一些物理量的表达没有唯一性,这时必须引入一些准则,比如表达式必须是对称的(两点之间的相互作用满足牛顿第三定律)

<math>\begin{align}

&(\partial \tilde P)_j= \partial \left(\frac{\tilde P}{\phi} \phi\right)_j = \phi_j\partial \left(\frac{\tilde P}{\phi} \right)_j+\left(\frac{\tilde P}{\phi}\right)_j\partial \left( \phi\right)_j \\ &= \phi_j \sum_i \left(\frac{\tilde P}{\phi\sigma^*}\right)_i\nu_i\nabla_jW_{ji}+\left(\frac{\tilde P}{\phi}\right)_j \sum_i \left(\frac{\phi}{\sigma^*}\right)_i\nu_i\nabla_jW_{ji}\\ &= \sigma^*_j \sum_i \left(\frac{\tilde P}{\sigma^*\sigma^*}\right)_i\nu_i\nabla_jW_{ji}+\left(\frac{\tilde P}{\sigma^*}\right)_j \sum_i \left(\frac{\sigma^*}{\sigma^*}\right)_i\nu_i\nabla_jW_{ji}\\ &=\sum_i \sigma^*_j\frac{\tilde P_i}{{\sigma^*_i}^2}\nu_i\nabla_jW_{ji}+\sum_i \frac{\tilde P_j}{\sigma^*_j}\nu_i\nabla_jW_{ji}\\ &= \sum_i \nu_i\sigma^*_j\left(\frac{\tilde P_i}{{\sigma^*_i}^2}+ \frac{\tilde P_j}{{\sigma^*_j}^2}\right)\nabla_jW_{ji} \end{align}</math>

其中第一行最后引入一个多余的任意的<math> \phi </math>,第二行为代入物理量空间偏导的定义,之后取特例 <math>\phi=\sigma^*</math> ,整理后得

<math>\begin{align}

(\partial {\tilde P})_i = \sum_j \nu_j \sigma_i^* \left(\frac{{\tilde P}_i}{\sigma_i^{*2}}+\frac{{\tilde P}_j}{\sigma_j^{*2}}\right) \bigtriangledown_i W_{ij} \end{align}</math>

另外一个运动方程中必须表达为SPH自由度的项是 <math>\partial_i v^i </math>.它涉及速度从而是对时间的导数(之后再对空间求偏导),与其他物理量(比如重子数密度)不同,我们并不直接计算,而是利用下面的步骤来达到目的.对于守恒荷密度 <math>\sigma</math> 对应的守恒流 <math>\sigma u</math> ,按定义

<math>\begin{align}

\left(\sigma u^{\mu}\right)_{;\mu} =0 \end{align}</math>

等式左边可以表达为

<math>\begin{align}

&\partial_{\mu} (\sigma^* \frac{u^{\mu}}{u^0}) =\partial_{\tau}\sigma^* + \partial_i (\sigma^* v^i) = 0 \\ &\partial_{\tau}\sigma^* + \sigma^* \partial_i v^i = 0 \end{align}</math>

其中第二项是我们需要的结果,而左边可以用SPH自由度表达,这时因为

<math>\begin{align}

\sigma^*_j =\sum_i {\nu}_i W_{ji} \end{align}</math>

对于固定的空间坐标点<math> j</math> 对时间求偏导,我们得到

<math>\begin{align}

\partial_{\tau}\sigma^*_j = \sum_i {\nu}_i \partial_{\tau}W_{ji}(|r_j-r_i|) =\sum_i{\nu}_i \left[\dot{r_j}\cdot \nabla_jW_{ji}(|r_j-r_i|)+\dot{r_i}\cdot\nabla_iW_{ji}(|r_j-r_i|)\right] =\sum_i {\nu}_i (\dot{r_j}-\dot{r_i}) \cdot \nabla_jW_{ji}(|r_j-r_i|)) \end{align}</math>

其中利用

<math>\begin{align}

\nabla_iW_{ji}(|r_j-r_i|)=-\nabla_jW_{ji}(|r_j-r_i|) \end{align}</math>

出于懒惰关系,其中位置矢量并没有加矢量箭头,其定义准确的写开为

<math>\begin{align}

&W_{ji}\equiv W_{ji}(|r_j-r_i|)\equiv W_{ji}(|\vec r_j-\vec r_i|) \\ &\nabla_jW_{ji}\equiv \nabla_jW_{ji}(|r_j-r_i|)\equiv \frac{\partial W_{ji}(|\vec r_j-\vec r_i|)}{\partial \vec r_j} \\ &\nabla_iW_{ji}\equiv \nabla_iW_{ji}(|r_j-r_i|)\equiv \frac{\partial W_{ji}(|\vec r_j-\vec r_i|)}{\partial \vec r_i} \end{align}</math>

最后得到

<math>\begin{align}

\left(\partial_i v^i\right)_j = - \frac{1}{\sigma^*}\sum_k \nu_k (\dot r_j - \dot r_k) \cdot \nabla_j W_{jk} \end{align}</math>

记号和基本关系式

这里我们定义一些记号,和下面需要反复用到的一些常用关系和基本表达式.首先我们引入对前面引入的记号的定义做一个修正.我们把大写的<math>D</math>保留给协变导数,而其他形式仍旧理解为普通的微分.公式中涉及的<math>\tau</math>大多数情况下是时间坐标而非固有时.符号<math>\partial</math>单独出现的时候,一般都是指代协变导数,因为这是之前的在平直空间中公式的自然的推广.注意到下面的定义是上面习惯定义的特例和冲突,我们定义如下

<math>\begin{align}

\frac{D}{D\tau}\equiv u^{\mu}D_{\mu} \end{align}</math>

注意到这个记号是非常容易产生混淆的,它的引入是为了与其他文献中的记号保持一致:但是其中等式左边的<math>\tau</math>被 解释为固有时而非时间坐标.从而,我们有下面的替换关系,

<math>\begin{align}

&D\tau\longleftrightarrow \frac{1}{\gamma}d\tau\\ &\frac{D}{d\tau}\longleftrightarrow \frac{1}{\gamma}u^{\mu}D_\mu \end{align}</math>

在其他情况下<math>\tau</math>都指代时间坐标,比如

<math>\begin{align}

d_\tau f \equiv \frac{df}{d\tau} =\frac{1}{\gamma}u^\mu\partial_\mu f \end{align}</math>

利用上述协变导数的记号和速度投影算子的定义 <math>\Delta^{\mu\nu}\equiv g^{\mu\nu}-u^\mu u^\nu </math> 我们可以把协变导数分开为平行于速度的部分和垂直与速度的部分,即

<math>\begin{align}

D^\mu=g^{\mu\nu}D_{\nu}=\Delta^{\mu\nu}D_\nu+u^\mu\frac{D}{D\tau}\equiv D^\mu_\perp+u^\mu\frac{D}{D\tau} \end{align}</math>

利用<math>\Delta_{\nu\beta}u^\beta=0,u_\beta u^\beta=1</math>以及部分积分,对度规协变导数为零<math>D_\mu g_{\alpha\beta}=0</math>不难直接证明

<math>\begin{align}

\Delta_{\nu\beta}D^\alpha u^\beta=-u^\beta D^\alpha\Delta_{\nu\beta}=D^\alpha u_\nu \end{align}</math>

从而有

<math>\begin{align}

\Delta_{\mu\alpha}\Delta_{\nu\beta}D^{\alpha} u^{\beta}=\Delta_{\mu\alpha}D^{\alpha} u_{\nu}=D_{\mu\perp} u_{\nu}=D_{\mu} u_{\nu}-u_{\mu}\frac{Du_{\nu}}{D\tau} \end{align}</math>

类似的带入投影算子定义以及剪切粘滞系数垂直于速度的性质<math>\pi^{\alpha\beta}u_alpha=\pi^{\alpha\beta}u_beta=0</math>可以直接证明

<math>\begin{align}

&\pi^{\alpha\beta}\frac{D}{D\tau}\Delta_{\alpha\beta}=0 \end{align}</math>

我们在这里给出一个使用频率很高的表达式,注意它仅在我们的双曲线坐标对应的度规下正确

<math>\begin{align}

\Gamma^\lambda_{\mu\nu}u_\lambda=\Gamma^0_{\mu\nu}\gamma+\Gamma^3_{\mu\nu}u_3 =\tau\delta_{\mu 3}\delta_{\nu 3}\gamma+\frac{u_3}{\tau}(\delta_{\mu 0}\delta_{\nu 3}+\delta_{\mu 3}\delta_{\nu 0}) \end{align}</math>

用它可以立刻证明下面的关系

<math>\begin{align}

u^{\mu}\Gamma^{\lambda}_{\mu\nu}u_{\lambda} =-\frac{1}{\tau^3}(u_3)^2\delta_{\nu 0} = -{\tau}(u^3)^2\delta_{\nu 0} \end{align}</math>

接下来我们引入一些记号的定义,和常用的表达式,这些结果都将在运动方程的推导中被反复使用.

<math>\begin{align}

&\partial \cdot u \equiv \left(u^{\mu}\right)_{;\mu} \\ &= \frac{1}{\tau} \partial_{\mu} \left(\tau u^{\mu}\right) \\ &= \frac{1}{\tau} \left(\partial_{\tau} \left(\tau \gamma\right)+\partial_i\left(\tau\gamma v^i \right)\right) \\ &= \frac{1}{\tau} \left(\partial_{\tau} \left(\tau \gamma\right)+v^i \partial_i\left(\tau\gamma \right)\right)+\frac{1}{\tau}(\tau\gamma)\partial_i v^i \\ &= \frac{1}{\tau} d_{\tau} (\tau\gamma) + \gamma \partial_i v^i \\ &= \frac{\gamma}{\tau} + d_{\tau}\gamma+\gamma \partial_i v^{i} \end{align}</math>

比如守恒荷对时间的微分就可以用上述表达式来得到

<math>\begin{align}

&0= D_\mu(\sigma u^{\mu})=\sigma(\partial\cdot u)+\gamma d_{\tau}\sigma \\ &d_\tau\sigma=\frac{d\sigma}{d\tau}=-\frac{\sigma}{\gamma}(\partial\cdot u) \end{align}</math>

我们可以进一步把方程右边逐项表达成上面选定的变量

<math>\begin{align}

&d_{\tau}\gamma = \frac{d}{d\tau}\left(1-g_{ij}u^iu^j\right)^{1/2} \\ &= -\frac{g_{ij}u^i}{\gamma}\frac{du^j}{d\tau}-\frac{u^iu^j}{2\gamma}\frac{dg_{ij}}{d\tau} \\ &= -\frac{u_j}{\gamma}\frac{du^j}{d\tau}+\frac{\left(u^3\right)^2\tau}{\gamma} \end{align}</math>

顺便注意到

<math>\begin{align}

&u_\mu\partial^\mu u_0=u_\mu\partial^\mu\gamma=\gamma d_\tau\gamma\\ &\frac{du_\mu}{d\tau}=\frac{d\gamma}{d\tau}\delta_{\mu 0}+\frac{du_i}{d\tau}\delta_{\mu i} \end{align}</math>

接着我们计算速度的空间导数.由

<math>\begin{align}

&\partial_k(u^\mu u_\mu)=0\\ &0=u^\mu\partial_ku_\mu=\gamma\partial_k\gamma+u^i\partial_ku_i\\ &\partial_k\gamma=-v^i\partial_ku_i\\ &0=\partial_k\gamma+v^i\partial_ku_i=\partial_k\gamma+v^iv_i\partial_k\gamma+u^i\partial_kv_i=\partial_k\gamma+\frac{(1-\gamma^2)}{\gamma^2}\partial_k\gamma+u^i\partial_kv_i\\ &\partial_k\gamma=-{\gamma^2}u^i\partial_kv_i\\ \end{align}</math>

上面第二个关系式虽然正确,但是数值误差大,故我们采用第一个关系式,即计算四速度的空间导数.

速度空间部分的协变固有时微分为

<math>\begin{align}

&u\cdot \partial u_{i} \equiv \gamma \frac{du_i}{d\tau} - u^{\mu} u_{\lambda} \Gamma^{\lambda}_{i\mu} \\ &= \gamma \frac{du_i}{d\tau} \\ &= - \gamma\frac{du^1}{d\tau}\delta_{i1}- \gamma\frac{du^2}{d\tau}\delta_{i2}- \left(\gamma\tau^2\frac{du^3}{d\tau}+2\tau\gamma u^3 \right)\delta_{i3} \end{align}</math>

<math>\begin{align}

\gamma\frac{du_i}{d\tau}=\gamma\frac{du_1}{d\tau}\delta_{i1}+\gamma\frac{du_2}{d\tau}\delta_{i2}+\gamma\frac{du_3}{d\tau}\delta_{i3} \end{align}</math>

其中第二步等式是利用了

<math>\begin{align}

u^{\mu} u_{\lambda} \Gamma^{\lambda}_{i\mu} = 0 \end{align}</math>

注意到在上式的书写中偏导符号并没有歧义,因为第二步等式相当于说

<math>\begin{align}

u\cdot \partial u_i \equiv u_{\mu}{(u_i)_;}^{\mu} = u_{\mu}{(u_i)_,}{^{\mu}} \end{align}</math>

在现有的双曲线坐标下,可以用度规的具体形式证明

<math>\begin{align}

\frac{\partial}{\partial x^{\beta}}\equiv\partial_\beta =g_{\mu\alpha}\frac{\partial}{\partial x_\alpha}=g_{\mu\alpha}\partial^\alpha \end{align}</math>

以及

<math>\begin{align}

\partial^\beta = g^{\mu\alpha}\partial_\alpha \end{align}</math>

值得提起注意,对于一般的度规,普通导数不是协变矢量,故而它们的上下标不再可以方便的调换,因为

<math>\begin{align}

&\partial^\mu \equiv \frac{\partial}{\partial x_\mu}=\frac{\partial}{\partial(g_{\mu\alpha}x^\alpha)}=\frac{\partial x^\beta}{\partial(g_{\mu\alpha}x^\alpha)}\frac{\partial}{\partial x^\beta} \\ &\partial_\beta \equiv \frac{\partial}{\partial x^\beta}=\frac{\partial(g_{\mu\alpha}x^\alpha)}{\partial x^\beta}\frac{\partial}{\partial x_\mu} =(g_{\mu\beta}+g_{\mu\alpha,\beta}x^\alpha)\frac{\partial}{\partial x_\mu} \end{align}</math>

从而在一般情况下

<math>\begin{align}

\partial_\beta \equiv \frac{\partial}{\partial x^\beta}\ne g_{\mu\beta}\frac{\partial}{\partial x_\mu}=g_{\mu\beta}\partial^\alpha \end{align}</math>

我们给出证明如下,对时间分量

<math>\begin{align}

\alpha=0=\tau \end{align}</math>

由于<math>x^0=x_0</math>,我们有

<math>\begin{align}

&\partial_0=\frac{\partial}{\partial x^0}=\frac{\partial}{\partial x_0}=g_{0\mu}\partial^\mu \\ &\partial^0=g^{0\mu}\partial_\mu \\ \end{align}</math>

最后一步利用了度规的对角化的具体形式.

对空间分量,因为度规的具体形式不涉及空间坐标,即<math>g_{\mu\alpha,i}=0</math>,我们有

<math>\begin{align}

&\partial_i \equiv \frac{\partial}{\partial x^i}=\left(\frac{\partial(g_{\mu\alpha} x^\alpha)}{\partial x^i}\right)\frac{\partial}{\partial x_\mu} =\left(g_{\mu i}+g_{\mu\alpha,i} x^\alpha \right)\frac{\partial}{\partial x_\mu}=g_{\mu i}\frac{\partial}{\partial x_\mu}=g_{i\mu}\partial^\mu\\ &\partial^i=g^{i\mu}\partial_\mu \end{align}</math>

综合两者即得

<math>\begin{align}

&\partial_\beta\equiv \frac{\partial}{\partial x^\beta}=g_{\mu\beta}\frac{\partial}{\partial x_\mu}=g_{\mu\beta}\partial^\alpha \\ &\partial^\beta= g^{\mu\alpha}\partial_\alpha \end{align}</math>

用类似的方法,我们可以证明一个下面计算中需要用到的关系

<math>\begin{align}

&T^{\mu\nu}\partial^\alpha(g_{\beta\nu}g_{\alpha\mu})\\ &=T^{\mu\nu}\partial^0(g_{\beta\nu}g_{0\mu})\\ &=T^{\mu\nu}\partial_0(g_{\beta\nu}g_{0\mu})\\ &=T^{\mu\nu}\delta_{\mu 0}\partial_0(g_{\beta\nu})\\ &=T^{\mu\nu}\delta_{\mu 0}\delta_{\beta 3}\delta_{\nu 3}\partial_0(-\tau^2)\\ &=T^{03}\delta_{\beta 3}\partial_0(-\tau^2)\\ &=T^{03}\delta_{\beta 3}(-2\tau)\\ &=\frac{2}{\tau}T^{03}g_{\beta 3} \end{align}</math>

这里我们把上面的常用关系总结如下

<math>\begin{align}

&D^\mu=g^{\mu\nu}D_{\nu}=\Delta^{\mu\nu}D_\nu+u^\mu\frac{D}{D\tau}\equiv D^\mu_\perp+u^\mu\frac{D}{D\tau}\\ &\Delta_{\nu\beta}D^\alpha u^\beta=D^\alpha u_\nu\\ &\Delta_{\nu\beta}D^\alpha u^\beta=-u^\beta D^\alpha\Delta_{\nu\beta}=D^\alpha u_\nu\\ &\Gamma^\lambda_{\mu\nu}u_\lambda =\tau\delta_{\mu 3}\delta_{\nu 3}\gamma+\frac{u_3}{\tau}(\delta_{\mu 0}\delta_{\nu 3}+\delta_{\mu 3}\delta_{\nu 0})\\ &u^{\mu}\Gamma^{\lambda}_{\mu\nu}u_{\lambda} = -{\tau}(u^3)^2\delta_{\nu 0}\\ &\partial_\beta\equiv \frac{\partial}{\partial x^\beta}=g_{\mu\beta}\partial^\alpha \\ &\partial^\beta= g^{\mu\alpha}\partial_\alpha \\ &u\cdot \partial u_i = u_{\mu}{(u_i)_;}^{\mu} = u_{\mu}{(u_i)_,}{^{\mu}} \end{align}</math>

把上面的记号和基本表达式总结如下

<math>\begin{align}

&\theta\equiv Du\equiv\partial \cdot u = \left(u^{\mu}\right)_{;\mu} = \frac{1}{\tau} \partial_{\mu} \left(\tau u^{\mu}\right) =\frac{\gamma}{\tau}+\partial_{\mu}u^{\mu} =\frac{\gamma}{\tau} + d_{\tau}\gamma+\gamma \partial_i v^{i} \\ &d_\tau\sigma=\frac{d\sigma}{d\tau}=-\frac{\sigma}{\gamma}(\partial\cdot u) \frac{du_\mu}{d\tau}=\frac{d\gamma}{d\tau}\delta_{\mu 0}+\frac{du_i}{d\tau}\delta_{\mu i}\\ &\frac{du_\mu}{d\tau}=\frac{d\gamma}{d\tau}\delta_{\mu 0}+\frac{du_i}{d\tau}\delta_{\mu i}\\ &d_{\tau}\gamma = \frac{d}{d\tau}\left(1-g_{ij}u^iu^j\right)^{1/2} = -\frac{g_{ij}u^i}{\gamma}\frac{du^j}{d\tau}-\frac{u^iu^j}{2\gamma}\frac{dg_{ij}}{d\tau} = -\frac{u_j}{\gamma}\frac{du^j}{d\tau}+\frac{\left(u^3\right)^2\tau}{\gamma}\\ &\partial_k\gamma=-v^i\partial_ku_i\\ &u\cdot \partial u_{i} = \gamma \frac{du_i}{d\tau} - u^{\mu} u_{\lambda} \Gamma^{\lambda}_{i\mu} = \gamma \frac{du_i}{d\tau} = - \gamma\frac{du^1}{d\tau}\delta_{i1}- \gamma\frac{du^2}{d\tau}\delta_{i2}- \left(\gamma\tau^2\frac{du^3}{d\tau}+2\tau\gamma u^3 \right)\delta_{i3} \\ &\gamma\frac{du_i}{d\tau}=\gamma\frac{du_1}{d\tau}\delta_{i1}+\gamma\frac{du_2}{d\tau}\delta_{i2}+\gamma\frac{du_3}{d\tau}\delta_{i3}\\ \end{align}</math>

自由度数目和方程数目

体系的运动方程共有1+3+1+5个,一个能量守恒,三个动量守恒,一个体粘滞系数,五个剪切粘滞系数.我们从粘滞系数方程开始讨论,然后讨论能动张量守恒的方程.最后我们总结给出所有的方程和对应的待求解的变量.原则上,所有对空间的导数我们按SPH模型进行计算,所有变量对时间坐标的导数,是我们方程需要求解的对象.如果一个张量分量同时具有时空指标.那些时空指标是哑变量,那么我们必须对它的时间指标部分求和空间指标部分求和分布对待,因为前者原则上将表达为待求解变量对时间坐标的全微分(当然会涉及到对时间坐标的偏微分和对时间坐标的全微分的转化),而后者原则上将转换为已知量.如果那些时空指标并不求和,那么原则上可以通过待求解变量的合理选取,仅仅考察对空间指标的偏导数.

剪切粘滞系数由5个独立的变量构成.在推导它们满足的方程之前,我们先把所有的剪切系数用五个独立的变量表达出来.下面我们通过表达式,把5个独立的变量选取为除去33以外的5个纯空间指标张量.由于剪切粘滞系数是对称的,所以有10个独立的分量.利用和速度正交的性质,我们去掉4个变量,我们有

<math>\begin{align}

&0=u_\mu\pi^{\mu\nu}=\gamma\pi^{0\nu}+u_i\pi^{i\nu}\\ &\pi^{0\nu}=-\frac{u_j}{\gamma}\pi^{j\nu}\\ &\pi^{0i}=-\frac{u_j}{\gamma}\pi^{ji}\\ &\pi^{00}=-\frac{u_j}{\gamma}\pi^{j0}=-\frac{u_j}{\gamma}\pi^{0j}=\frac{u_iu_j}{\gamma^2}\pi^{ij} \end{align}</math>

这样所有的含时间指标的剪切粘滞系数都化为空间指标的剪切系数的函数.最后我们利用无迹条件再消去一个独立的剪切系数,我们得到一个空间指标的剪切系数的方程

<math>\begin{align}

&0=g_{\mu\nu}\pi^{\mu\nu}=\pi^{00}-\pi^{11}-\pi^{22}-\tau^2\pi^{33}\\ &\pi^{00}=\pi^{11}+\pi^{22}+\tau^2\pi^{33}=\frac{u_iu_j}{\gamma^2}\pi^{ij} \end{align}</math>

由此我们可以把其中的一个空间指标剪切系数,表达成其他系数的函数,比如

<math>\begin{align}

&\left(1-\frac{u_1^2}{\gamma^2}\right)\pi^{11}+\left(1-\frac{u_2^2}{\gamma^2}\right)\pi^{22}+\left(1-\frac{u_3^2}{\gamma^2}\right)\tau^2\pi^{33}=\sum_{i\ne j}\frac{u_iu_j}{\gamma^2}\pi^{ij}\\ &\pi^{33}=\left(1-\frac{u_3^2}{\gamma^2}\right)^{-1}\frac{1}{\tau^2}\left[\sum_{i\ne j}\frac{u_iu_j}{\gamma^2}\pi^{ij}-\left(1-\frac{u_1^2}{\gamma^2}\right)\pi^{11}-\left(1-\frac{u_2^2}{\gamma^2}\right)\pi^{22} \right] \end{align}</math>

体粘滞系数方程

运动方程涉及的关于体粘滞系数的变量是<math>\frac{\Pi}{\sigma}</math>. 这样的选取是方便的,因为两个密度量相除,体积项消去,因为分母上是守恒的常数,所以直接得到某个SPH粒子的总体粘滞量的变化.

<math>\begin{align}

\widetilde{\Pi}=-\int_{\tau_0}^{\tau} d\tau'G(\tau,\tau')\frac{\zeta'\theta'}{\sigma'}+e^{-(\tau-\tau_0)/\tau_R}\widetilde{\Pi}_0 \end{align}</math>

利用格林函数的需要满足的条件

<math>\begin{align}

&G(\tau,\tau;\tau_R) \equiv \frac{1}{\tau_R}\\ &\frac{DG(\tau,\tau';\tau_R)}{D\tau}=-\frac{G(\tau,\tau';\tau_R)}{\tau_R} \end{align}</math>

得到具体形式

<math>\begin{align}

G(\tau,\tau';\tau_R) \equiv \frac{1}{\tau_R}e^{-(\tau-\tau')/\tau_R} \end{align}</math>

由此得到体粘滞系数满足的方程

<math>\begin{align}

\dot{\widetilde{\Pi}}=-\frac{\widetilde{\Pi}}{\tau_R}-\frac{\zeta\theta}{\sigma\tau_R} \end{align}</math>

进行如下的推广.按Gabriel的笔记,我们将对故有时的导数记为四速度与协变导数的缩并.(这实际上是一个模型,比如我们可以取为对时间坐标的导数而非对固有时的协变导数,这样会多出一个因子<math>\gamma</math>.)由于体粘滞为标量,故这样定义的量对应对固有时的微分.

<math>\begin{align}

\dot{\widetilde{\Pi}}\equiv\frac{D\widetilde{\Pi}}{D\tau} \end{align}</math>

同时我们将体系的膨胀速率替换为协变散度

<math>\begin{align}

\theta\equiv\partial\cdot u \end{align}</math>

由此我们得到我们需要的关于体粘滞系数的方程.一方面,我们有,

<math>\begin{align}

\frac{D\widetilde{\Pi}}{D\tau}=u^\mu D_\mu\widetilde{\Pi}=u^\mu \partial_\mu\widetilde{\Pi}=\gamma d_\tau\left(\frac{\Pi}{\sigma} \right) \end{align}</math>

另一方面

<math>\begin{align}

d_{\tau}\left(\frac{\Pi}{\sigma}\right) = \frac{d_{\tau}\Pi}{\sigma}-\frac{\Pi}{\sigma^2}d_{\tau}\sigma = \frac{d_{\tau}\Pi}{\sigma} + \frac{\Pi}{\sigma^2}\left(\frac{\sigma}{\gamma}\partial \cdot u \right) = \frac{d_{\tau}\Pi}{\sigma} + \frac{\Pi}{\sigma}\left(\frac{1}{\gamma}\partial \cdot u \right) \end{align}</math>

其中利用关系

<math>\begin{align}

0 = \left(\sigma u^{\mu}\right)_{;\mu} = u \cdot \partial \sigma + \sigma \partial \cdot u \end{align}</math>

从而

<math>\begin{align}

d_\tau\sigma=\frac{d\sigma}{d\tau}=-\frac{\sigma}{\gamma}(\partial\cdot u) \end{align}</math>

联合两者,我们得到需要的方程

<math>\begin{align}

d_{\tau}\left(\frac{\Pi}{\sigma}\right) = \frac{d_{\tau}\Pi}{\sigma} + \frac{\Pi}{\sigma^2}\left(\frac{\sigma}{\gamma}\partial \cdot u\right) =\frac{1}{\gamma}\left(-\frac{\widetilde{\Pi}}{\tau_R}-\frac{\zeta\theta}{\sigma\tau_R} \right) =-\frac{1}{\gamma\tau_R\sigma}\left(\Pi+{\zeta\partial\cdot u}\right) \end{align}</math>

从中,我们可以得到另一个被用到的关系

<math>\begin{align}

{d_{\tau}\Pi}=-\frac{1}{\gamma\tau_R}\left(\Pi+{\left(\zeta+{\tau_R\Pi} \right)\partial\cdot u}\right) \end{align}</math>

注意上面书写中采用的一个默认的约定.当我们把一个表达式写为一些记号和基本表达式时,我们停止继续带入的过程.否则任何公式都将无限的增加其篇幅,同时导致很多不必要的誊写错误.

剪切粘滞系数方程

剪切粘滞系数满足类似与体粘滞系数的方程.剪切粘滞系数的积分形式考虑与守恒荷密度的比对格林函数的积分

<math>\begin{align}

\widetilde{\pi}^{\mu\nu}=\int^\tau_{\tau_0}d\tau'G(\tau,\tau')\frac{\eta\sigma^{\mu\nu}(\tau')}{\sigma}+e^{-(\tau-\tau_0)/\tau_{R_\pi}}\widetilde{\pi}^{\mu\nu}_0 \end{align}</math>

其中利用了得到张量对称无迹部分的投影算符,

<math>\begin{align}

&\sigma^{\mu\nu}\equiv \Delta^{\mu\nu\alpha\beta}D_\alpha u_\beta\\ &\Delta^{\mu\nu\alpha\beta} \equiv \frac{1}{2}\left[\Delta^{\mu\alpha}\Delta^{\nu\beta}+\Delta^{\mu\beta}\Delta^{\nu\alpha}-\frac{2}{3}\Delta^{\mu\nu}\Delta^{\alpha\beta}\right] \end{align}</math>

其实等式的两边其实都受到该投影算符的作用,这样做是基于Curie原理,即能动张量流的无迹部分不受到力(热力学量的梯度)的有迹的影响.显然它满足下面的方程

<math>\begin{align}

&\frac{D\widetilde{\pi}^{\mu\nu}}{D\tau}=-\frac{\widetilde{\pi}^{\mu\nu}}{\tau_R}+\frac{\eta\sigma^{\mu\nu}}{\sigma\tau_R}\\ &\frac{D\widetilde{\pi}^{\mu\nu}}{D\tau}=-\frac{{\pi}^{\mu\nu}}{\sigma\tau_R}+\frac{\eta\sigma^{\mu\nu}}{\sigma\tau_R} \end{align}</math>

方程的左边可以进一步化简为

<math>\begin{align}

&\frac{D\widetilde{\pi}^{\mu\nu}}{D\tau} =u_\alpha D^\alpha\widetilde{\pi}^{\mu\nu} =u_\alpha D^\alpha\left(\frac{\pi^{\mu\nu}}{\sigma}\right)\\ &=\frac{1}{\sigma}u_\alpha D^\alpha{\pi^{\mu\nu}}-\frac{\pi^{\mu\nu}}{\sigma^2}u_\alpha\partial^\alpha{\sigma}\\ &=\frac{1}{\sigma}u_\alpha D^\alpha{\pi^{\mu\nu}}+\frac{\pi^{\mu\nu}}{\sigma^2}\left({\sigma}\partial \cdot u \right)\\ &=\frac{1}{\sigma}u_\alpha D^\alpha{\pi^{\mu\nu}}+\frac{\pi^{\mu\nu}}{\sigma}\left(\partial \cdot u \right) \end{align}</math>

由此可得

<math>\begin{align}

u_\alpha D^\alpha{\pi^{\mu\nu}}=-\frac{{\pi}^{\mu\nu}}{\tau_R}+\frac{\eta\sigma^{\mu\nu}}{\tau_R}-{\pi^{\mu\nu}}\left(\partial \cdot u \right) \end{align}</math>

方程的右边是无迹对称的,所以我们可以用投影算符作用于方程左边,这样做是为了减少数值计算中的误差,我们得到

<math>\begin{align}

\Delta_{\mu\nu\alpha\beta}u_\lambda D^\lambda{\pi^{\alpha\beta}}=-\frac{{\pi}_{\mu\nu}}{\tau_R}+\frac{\eta\sigma_{\mu\nu}}{\tau_R}-{\pi_{\mu\nu}}\left(\partial \cdot u \right) \end{align}</math>

下面我们逐项进行计算.目标是把协变导数和投影算子都替换为普通导数.

<math>\begin{align}

&\eta\sigma_{\mu\nu}=\eta\frac{1}{2}\left[\Delta_{\mu\alpha}\Delta_{\nu\beta}+\Delta_{\mu\beta}\Delta_{\nu\alpha}-\frac{2}{3}\Delta_{\mu\nu}\Delta_{\alpha\beta}\right]D^\alpha u^\beta\\ &=\frac{\eta}{2}\Delta_{\mu\alpha}\Delta_{\nu\beta}\left[D^\alpha u^\beta+D^\beta u^\alpha\right]-\frac{\eta}{3}\Delta_{\mu\nu}\Delta_{\alpha\beta}D^\alpha u^\beta\\ &=\frac{\eta}{2}[\partial_\mu u_\nu+\partial_\nu u_\mu]-\frac{\eta\gamma}{2}\left[u_\mu\frac{du_\nu}{d\tau}+u_\nu\frac{du_\mu}{d\tau}\right]-\frac{\eta}{3}\Delta_{\mu\nu}(\partial\cdot u)-\eta\gamma\tau\delta_{\nu 3}\delta_{\mu 3}-\eta\frac{u_3}{\tau}[\delta_{\mu 0}\delta_{\nu 3}+\delta_{\mu 3}\delta_{\nu 0}]-\frac{\eta(u_3)^2}{2\tau^3}(u_\mu\delta_{\nu 0}+u_\nu\delta_{\mu 0}) \end{align}</math>

其中利用关系

<math>\begin{align}

u_{\mu}u^{\alpha}\Gamma^{\lambda}_{\alpha\nu}u_{\lambda} =-u_\mu\frac{1}{\tau^3}(u_3)^2\delta_{\nu 0} = -u_\mu{\tau}(u^3)^2\delta_{\nu 0} \end{align}</math>

在这里我们不需要过于为上面涉及的四速度对任意坐标的偏导而担心,因为指标并不涉及求和,我们最终可以利用剪切粘滞系数的依赖关系,在实际计算中仅仅计算四速度空间部分的坐标空间导数.

另一方面,剪切粘滞系数等式左边可以写为

<math>\begin{align}

\tau_R\Delta_{\mu\nu\alpha\beta}\frac{D\pi^{\alpha\beta}}{D\tau} =\tau_R\frac{D\pi_{\mu\nu}}{D\tau}-\tau_R\pi^{\alpha\beta}\frac{D\Delta_{\mu\nu\alpha\beta}}{D\tau} \end{align}</math>

其中

<math>\begin{align}

&-\tau_R\pi^{\alpha\beta}\frac{D\Delta_{\mu\nu\alpha\beta}}{D\tau}\\ &=-\frac{1}{2}\tau_R\pi^{\alpha\beta}\frac{D}{D\tau}\left[\Delta^{\mu\alpha}\Delta^{\nu\beta}+\Delta^{\mu\beta}\Delta^{\nu\alpha}-\frac{2}{3}\Delta^{\mu\nu}\Delta^{\alpha\beta}\right]\\ &=\tau_R(u_\mu\pi^\alpha_\nu+u_\nu\pi^\alpha_\mu)\frac{Du_\alpha}{D\tau}\\ &=\tau_R\gamma(u_\mu\pi^\alpha_\nu+u_\nu\pi^\alpha_\mu)\frac{du_\alpha}{d\tau}+\tau_R(u_\mu\pi_{\nu 0}+u_\nu\pi_{\mu 0})\frac{(u_3)^2}{\tau^3}\\ &=\tau_R\gamma(u_\mu\pi^0_\nu+u_\nu\pi^0_\mu)\frac{d\gamma}{d\tau}+\tau_R\gamma(u_\mu\pi^j_\nu+u_\nu\pi^j_\mu)\frac{du_j}{d\tau}+\tau_R(u_\mu\pi_{\nu 0}+u_\nu\pi_{\mu 0})\frac{(u_3)^2}{\tau^3}\\ &=\tau_R(\gamma u_\mu\pi^j_\nu+\gamma u_\nu\pi^j_\mu-u_\mu\pi^0_\nu u^j-u_\nu\pi^0_\mu u^j)\frac{du_j}{d\tau} \end{align}</math>

最后一步利用了对<math>\gamma</math>因子导数的明显形式,原则上不需要,但是这样可以化简结果的形式

<math>\begin{align}

\tau_R\frac{D\pi_{\mu\nu}}{D\tau} =\gamma\tau_R\frac{d\pi_{\mu\nu}}{d\tau}+\frac{\tau_R}{\tau^3}u_3(\delta_{\mu 0}\pi_{\nu 3}+\delta_{\nu 0}\pi_{\mu 3})+\frac{\tau_R}{\tau}u_3(\delta_{\mu 3}\pi_{\nu 0}+\delta_{\nu 3}\pi_{\mu 0}) -\frac{\gamma\tau_R}{\tau}(\delta_{\mu 3}\pi_{\nu 3}+\delta_{\nu 3}\pi_{\mu 3}) \end{align}</math>

两项贡献之和为

<math>\begin{align}

&\tau_R\Delta_{\mu\nu\alpha\beta}\frac{D\pi^{\alpha\beta}}{D\tau}\\ &=\tau_R\gamma(u_\mu\pi^0_\nu+u_\nu\pi^0_\mu)\frac{d\gamma}{d\tau}+\tau_R\gamma(u_\mu\pi^j_\nu+u_\nu\pi^j_\mu)\frac{du_j}{d\tau}+\tau_R(u_\mu\pi_{\nu 0}+u_\nu\pi_{\mu 0})\frac{(u_3)^2}{\tau^3}\\ &+\gamma\tau_R\frac{d\pi_{\mu\nu}}{d\tau}+\frac{\tau_R}{\tau^3}u_3(\delta_{\mu 0}\pi_{\nu 3}+\delta_{\nu 0}\pi_{\mu 3}) +\frac{\tau_R}{\tau}u_3(\delta_{\mu 3}\pi_{\nu 0}+\delta_{\nu 3}\pi_{\mu 0}) -\frac{\gamma\tau_R}{\tau}(\delta_{\mu 3}\pi_{\nu 3}+\delta_{\nu 3}\pi_{\mu 3}) \end{align}</math>

综合以上得到关于剪切粘滞的方程

<math>\begin{align}

&\tau_R(\gamma u_\mu\pi^j_\nu+\gamma u_\nu\pi^j_\mu-u_\mu\pi^0_\nu u^j-u_\nu\pi^0_\mu u^j)\frac{du_j}{d\tau} +\gamma\tau_R\frac{d\pi_{\mu\nu}}{d\tau}+\frac{\tau_R}{\tau^3}u_3(\delta_{\mu 0}\pi_{\nu 3}+\delta_{\nu 0}\pi_{\mu 3})\\ &+\frac{\tau_R}{\tau}u_3(\delta_{\mu 3}\pi_{\nu 0}+\delta_{\nu 3}\pi_{\mu 0}) -\frac{\gamma\tau_R}{\tau}(\delta_{\mu 3}\pi_{\nu 3}+\delta_{\nu 3}\pi_{\mu 3})+\pi_{\mu\nu}\\ &=+\frac{\eta}{2}[\partial_\mu u_\nu+\partial_\nu u_\mu]-\frac{\eta\gamma}{2}\left[u_\mu\frac{du_\nu}{d\tau}+u_\nu\frac{du_\mu}{d\tau}\right]-\frac{\eta}{3}\Delta_{\mu\nu}(\partial\cdot u)\\ &-\eta\gamma\tau\delta_{\nu 3}\delta_{\mu 3}-\eta\frac{u_3}{\tau}[\delta_{\mu 0}\delta_{\nu 3}+\delta_{\mu 3}\delta_{\nu 0}]-\frac{\eta(u_3)^2}{2\tau^3}(u_\mu\delta_{\nu 0}+u_\nu\delta_{\mu 0})-\tau_R\pi_{\mu\nu}(\partial\cdot u) \end{align}</math>

如果我们以 <math>\begin{align} \frac{\pi_{\mu\nu}}{\sigma} \end{align}</math> 作为方程的待解变量,不需要太多的额外劳动,因为考虑了投影算子的方程的左边可以写为

<math>\begin{align}

\tau_R\Delta_{\mu\nu\alpha\beta}\frac{D\widetilde{\pi}^{\alpha\beta}}{D\tau} =\tau_R\frac{D\widetilde{\pi}_{\mu\nu}}{D\tau}-\frac{\tau_R\pi^{\alpha\beta}}{\sigma}\frac{D\Delta_{\mu\nu\alpha\beta}}{D\tau}\\ \end{align}</math>

从而

<math>\begin{align}

&\frac{D\widetilde{\pi}_{\mu\nu}}{D\tau}=\frac{\pi^{\alpha\beta}}{\sigma}\frac{D\Delta_{\mu\nu\alpha\beta}}{D\tau}-\frac{{\pi}^{\mu\nu}}{\sigma\tau_R}+\frac{\eta\sigma^{\mu\nu}}{\sigma\tau_R}\\ &\tau_R\frac{D\widetilde{\pi}_{\mu\nu}}{D\tau}=\frac{\tau_R\pi^{\alpha\beta}}{\sigma}\frac{D\Delta_{\mu\nu\alpha\beta}}{D\tau}-\frac{{\pi}^{\mu\nu}}{\sigma}+\frac{\eta\sigma^{\mu\nu}}{\sigma} \end{align}</math>

和之前的方程相比,方程的右边增加了一项.利用前面的结果,我们需要把方程的左边写成普通导数的形式,借鉴前面的结果,我们直接写下

<math>\begin{align}

\tau_R\frac{D\widetilde{\pi}_{\mu\nu}}{D\tau} =\gamma\tau_R\frac{d\widetilde{\pi}_{\mu\nu}}{d\tau}+\frac{\tau_R}{\sigma\tau^3}u_3(\delta_{\mu 0}\pi_{\nu 3}+\delta_{\nu 0}\pi_{\mu 3})+\frac{\tau_R}{\sigma\tau}u_3(\delta_{\mu 3}\pi_{\nu 0}+\delta_{\nu 3}\pi_{\mu 0}) -\frac{\gamma\tau_R}{\sigma\tau}(\delta_{\mu 3}\pi_{\nu 3}+\delta_{\nu 3}\pi_{\mu 3}) \end{align}</math>

对应的关于单位守恒荷的剪切粘滞系数的运动方程为

<math>\begin{align}

&\frac{\tau_R}{\sigma}(\gamma u_\mu\pi^j_\nu+\gamma u_\nu\pi^j_\mu-u_\mu\pi^0_\nu u^j-u_\nu\pi^0_\mu u^j)\frac{du_j}{d\tau} +\gamma\tau_R\frac{d\widetilde{\pi}_{\mu\nu}}{d\tau}+\frac{\tau_R}{\sigma\tau^3}u_3(\delta_{\mu 0}\pi_{\nu 3}+\delta_{\nu 0}\pi_{\mu 3})\\ &+\frac{\tau_R}{\sigma\tau}u_3(\delta_{\mu 3}\pi_{\nu 0}+\delta_{\nu 3}\pi_{\mu 0}) -\frac{\gamma\tau_R}{\sigma\tau}(\delta_{\mu 3}\pi_{\nu 3}+\delta_{\nu 3}\pi_{\mu 3})+\frac{\pi_{\mu\nu}}{\sigma}\\ &=\frac{\eta}{2\sigma}[\partial_\mu u_\nu+\partial_\nu u_\mu]-\frac{\eta\gamma}{2\sigma}\left[u_\mu\frac{du_\nu}{d\tau}+u_\nu\frac{du_\mu}{d\tau}\right]-\frac{\eta}{3\sigma}\Delta_{\mu\nu}(\partial\cdot u)\\ &-\frac{\eta\gamma\tau}{\sigma}\delta_{\nu 3}\delta_{\mu 3}-\eta\frac{u_3}{\sigma\tau}[\delta_{\mu 0}\delta_{\nu 3}+\delta_{\mu 3}\delta_{\nu 0}]-\frac{\eta(u_3)^2}{2\sigma\tau^3}(u_\mu\delta_{\nu 0}+u_\nu\delta_{\mu 0}) \end{align}</math>

为了之后的需要,我们记下下面的关系

<math>\begin{align}

&\frac{d\widetilde{\pi}_{\mu\nu}}{d\tau}=\frac{d\left(\frac{\pi_{\mu\nu}}{\sigma} \right)}{d\tau}=\frac{1}{\sigma}\frac{d{\pi}_{\mu\nu}}{d\tau}-\frac{{\pi}_{\mu\nu}}{\sigma^2}d_\tau\sigma=\frac{1}{\sigma}\frac{d{\pi}_{\mu\nu}}{d\tau}+\frac{{\pi}_{\mu\nu}}{\sigma\gamma}(\partial\cdot u)\\ &\frac{d{\pi}_{\mu\nu}}{d\tau}=\sigma\frac{d\widetilde{\pi}_{\mu\nu}}{d\tau}-\frac{{\pi}_{\mu\nu}}{\gamma}(\partial\cdot u) \end{align}</math>

能动张量守恒方程,空间部分

这部分方程包含了速度对时间坐标的导数.故在移项以后,方程可以被读作"加速度由力决定"的形式. 可以把能动张量守恒写为下面的形式

<math>\begin{align}

\frac{1}{\tau}\partial^\alpha(\tau T_{\alpha\beta})+g_{\beta 0}\tau T^{33}=0 \end{align}</math>

注意这个形式仅仅出现在Gabriel的笔记中,在推导的其余部分没有被用来得到运动方程,但是其时间分量被用来验证方程的解的自洽性(但是因为时间分量的方程独立与空间部分运动方程,故对理想流体意义不大,在粘滞流体,它包括了能动张量守恒在平行速度分量上的一部分信息).证明过程很大程度上依赖于度规的具体形式.

为了得到能动张量守恒空间部分对应的运动方程,我们利用能动张量的具体形式

<math>\begin{align}

T^{\mu\nu}=(\varepsilon+P+\Pi)u^\mu u^\nu-(P+\Pi)g^{\mu\nu}+\pi^{\mu\nu} \end{align}</math>

从能动张量守恒的协变形式出发,能动张量守恒空间部分的三个分量对应的方程可以写为

<math>\begin{align}

\sigma\gamma\frac{d}{d\tau}\left(\frac{\varepsilon+P+\Pi}{\sigma}u_i\right)+\frac{1}{\tau}\partial^\mu(\tau\pi_{\mu i})=\partial_i(P+\Pi) \end{align}</math>

上面方程的左边第一项可以进一步展开

<math>\begin{align}

&\sigma\gamma\frac{d}{d\tau}\left(\frac{\varepsilon+P+\Pi}{\sigma}u_i\right)\\ &=\gamma(\varepsilon+P+\Pi)\frac{du_i}{d\tau}-u_i(\varepsilon+P)\frac{\gamma}{\sigma}\frac{d\sigma}{d\tau}+u_i{\gamma}\frac{d(\varepsilon+P)}{d\tau}+u_i\gamma\frac{d\widetilde{\Pi}}{d\tau}\\ &=\gamma(\varepsilon+P+\Pi)\frac{du_i}{d\tau}+u_i(\varepsilon+P)(\partial\cdot u)+u_i{\gamma}\frac{d(\varepsilon+P)}{d\tau}+u_i\gamma\frac{d\widetilde{\Pi}}{d\tau}\\ \end{align}</math>

上面等号右边的表达式,除了第三项,根据前面的结果,我们都可以直接得到.我们把第三项的导数部分写为

<math>\begin{align}

\frac{d(\varepsilon+P)}{d\tau} =\frac{d(\varepsilon+P)}{d\varepsilon}\frac{d\varepsilon}{d\tau} =\left(1+\frac{dP}{d\varepsilon}\right)\frac{d\varepsilon}{d\tau} =\left(1+\frac{dP}{d\varepsilon}\right)\frac{1}{\gamma}u^\mu D_\mu\varepsilon =\left(1+\frac{dP}{d\varepsilon}\right)\frac{1}{\gamma}\left({T\sigma\gamma}d_\tau\left(\frac{s}{\sigma}\right)-(\varepsilon+P)(\partial\cdot u)\right) \end{align}</math>

最后一步利用热力学关系

<math>\begin{align}

&u \cdot \partial {\varepsilon} = u \cdot \partial {\varepsilon}+\omega \partial \cdot u - \omega \partial \cdot u \\ &= T u \cdot \partial s + T s \partial \cdot u + \sum_p \mu_p u \cdot \partial n_p + \sum_p \mu_p n_p \partial \cdot u - \omega \partial\cdot u \\ &= T \partial \cdot (su) + \sum_p \mu_p \partial \cdot (n_pu) - \omega \partial \cdot u \\ &= T \partial \cdot (su) - \omega \partial \cdot u \\ &={T \sigma\gamma}d_{\tau}\left(\frac{s}{\sigma}\right) - \omega \partial \cdot u \end{align}</math>

其中用到守恒流满足的关系

<math>\begin{align}

\sum_p \mu_p \partial \cdot (n_p u)=0 \end{align}</math>

还涉及热力学第一定律

<math>\begin{align}

&\partial \varepsilon=\partial(\frac{E}{V})=\frac{V\partial E-E\partial V}{V^2} \\ &\partial E =T\partial S -P \partial V +\sum_p \mu \partial N_p\\ &\partial S = \partial (sV)=V\partial s+s\partial V\\ &\sum_p \mu \partial N_p=\sum_p \mu V\partial n_p+\sum_p \mu n_p \partial V\\ &\partial \varepsilon =T\partial s+\sum_p\mu_p \partial n_p \end{align}</math>

其中关于压强的项自然的被抵消.最后有用到广延量的性质

<math>\begin{align}

\omega=\varepsilon+p=Ts-p+\sum_p \mu \partial n_p+p=Ts+\sum_p \mu \partial n_p \end{align}</math>

方程左边的第二项为,涉及到把对时间的偏导化为对时间坐标的全微分,对空间的偏导可以保留

<math>\begin{align}

&\frac{1}{\tau}\partial^\mu(\tau\pi_{\mu i})=\partial^0\pi_{0i}+\partial^j\pi_{ji}+\frac{1}{\tau}\pi_{\mu i}\partial^\mu\tau=\frac{d\pi_{0i}}{d\tau}-v^j\partial_j\pi_{0i}+\partial^j\pi_{ji}+\frac{\pi_{0i}}{\tau}\\ &=\frac{d\pi_{0i}}{d\tau}-v_j\partial^j\pi_{0i}+\partial^j\pi_{ji}+\frac{\pi_{0i}}{\tau} =\sigma\frac{d\widetilde{\pi}_{0i}}{d\tau}-v_j\partial^j\pi_{0i}+\partial^j\pi_{ji}+\frac{\pi_{0i}}{\tau}-\frac{{\pi}_{0i}}{\gamma}(\partial\cdot u)\\ \end{align}</math>

把上述结果综合起来得到

<math>\begin{align}

&\partial_i(P+\Pi)= \gamma(\varepsilon+P+\Pi)\frac{du_i}{d\tau}+u_i(\varepsilon+P)(\partial\cdot u)+u_i{\gamma}\left(1+\frac{dP}{d\varepsilon}\right)\frac{1}{\gamma}\left({T\sigma\gamma}d_\tau\left(\frac{s}{\sigma}\right)-(\varepsilon+P)(\partial\cdot u)\right)\\ &+u_i\gamma\frac{d\widetilde{\Pi}}{d\tau}+\sigma\frac{d\widetilde{\pi}_{0i}}{d\tau}-v_j\partial^j\pi_{0i}+\partial^j\pi_{ji}+\frac{\pi_{0i}}{\tau}-\frac{{\pi}_{0i}}{\gamma}(\partial\cdot u) \end{align}</math>

能动张量守恒方程,与速度平行部分,熵流方程

我们不采用能动张量守恒方程的时间部分作为我们的最后一个方程,而是采用能动张量守恒方程与速度平行的部分.我们知道,与速度内积后得到的表达式物理上决定体系的熵流.对于非理想流体,熵流并不守恒. 能动方程在速度方向上的投影为

<math>\begin{align}

&0=u_\mu D_\nu T^{\mu\nu}\\ &=u_\mu D_\nu((\varepsilon+P+\Pi)u^\mu u^\nu-(P+\Pi)g^{\mu\nu}+\pi^{\mu\nu})\\ &=\frac{D\varepsilon}{D\tau}+(\varepsilon+P+\Pi)\theta+u_\mu D_\nu\pi^{\mu\nu}\\ &=\frac{D\varepsilon}{D\tau}+(\varepsilon+P+\Pi)\theta-\pi^{\mu\nu}D_\nu u_\mu\\ &=\frac{D\varepsilon}{D\tau}+(\varepsilon+P+\Pi)\theta-\pi^{\mu\nu}\sigma_{\mu\nu}\\ \end{align}</math>

其中利用了(在任意度规下)把一个协变矢量按平行与垂直于速度进行分解的方法

<math>\begin{align}

&D_{\mu}u_{\nu}=u_{\mu}\frac{Du_{\nu}}{D\tau}+D_{\mu\perp}u_{\nu} \\ &=u_{\mu}\frac{Du_{\nu}}{D\tau}+\frac{1}{2}(D_{\mu\perp}u_{\nu}+D_{\nu\perp}u_{\mu})+\frac{1}{2}(D_{\mu\perp}u_{\nu}-D_{\nu\perp}u_{\mu})\\ &=u_{\mu}\frac{Du_{\nu}}{D\tau}+\sigma_{\mu\nu}+\frac{1}{3}\Delta_{\mu\nu}\theta+\Omega_{\mu\nu}\\ &\sigma_{\mu\nu}\equiv\frac{1}{2}(\nabla_{\mu}u_{\nu}+\nabla_{\nu}u_{\mu}-\frac{2}{3}\Delta_{\mu\nu}\theta)\\ &\theta\equiv D_{\mu}u^{\mu}\\ &\Omega_{\mu\nu}\equiv\frac{1}{2}(D_{\mu\perp}u_{\nu}-D_{\nu\perp}u_{\mu})\\ &A^{<\mu\nu>}\equiv\Delta^{\mu\nu\alpha\beta}A_{\alpha\beta}\\ &\nabla^{<\mu}u^{\nu>}\equiv 2\sigma^{\mu\nu}=\Delta^{\mu\nu\alpha\beta}D_\alpha u_\beta=2\nabla^{(\mu}u^{\nu)}-\frac{2}{3}\Delta^{\mu\nu}\theta\\ &\nabla^{(\mu}u^{\nu)}\equiv D_{\mu\perp}u^{\nu}+D_{\nu\perp}u^{\mu} \end{align}</math>

利用热力学关系以及粒子守恒流的条件

<math>\begin{align}

&d\varepsilon=Tds+\mu_P dn_P\\ &\frac{D\varepsilon}{D\tau}=T\frac{Ds}{D\tau}+\mu_P\frac{Dn_P}{D\tau}\\ &\frac{D\varepsilon}{D\tau}=u^\mu D_\mu\epsilon=Tu^\mu D_\mu s+\mu_P u^\mu D_\mu n_P =TD_\mu(su^\mu)+\mu_PD_\mu(n_Pu^\mu)-(Ts+\mu_Pn_P)\theta\\ &\varepsilon=-P+Ts+\mu_Pn_P \end{align}</math>

我们得到

<math>\begin{align}

TD_\mu(su^\mu)=\Pi(-\theta)+\pi^{\mu\nu}\sigma_{\mu\nu}=\Pi(-(\partial\cdot u))+\pi^{\mu\nu}\sigma_{\mu\nu} \end{align}</math>

下面我们把它改写为待解变量对时间的微商的方程的形式.比较复杂的是涉及求和的最后一项,我们不得不把它写开,把时间指标的求和算出来

<math>\begin{align}

&\pi^{\mu\nu}\sigma_{\mu\nu}=\pi^{\mu\nu}D_\mu u_\nu=\pi^{\mu\nu}\partial_\mu u_\nu-\tau\gamma\pi^{33}-\frac{2u_3}{\tau}\pi^{03}\\ &=\left(\pi^{0j}-\pi^{00}\frac{u^j}{\gamma}\right)\frac{du_j}{d\tau}-\frac{(u_3)^2}{\gamma\tau^3}\pi^{00}+(\pi^{ij}+\pi^{00}v^iv^j-\pi^{i0}v^j-\pi^{0j}v^i)\partial_iu_j-\tau\gamma\pi^{33}-\frac{2u_3}{\tau}\pi^{03}\\ &=\left(\pi_0^j-\pi_{00}\frac{u^j}{\gamma}\right)\frac{du_j}{d\tau}+(\pi^{ij}+\pi^{00}v^iv^j-\pi^{i0}v^j-\pi^{0j}v^i)\partial_iu_j-\frac{\gamma}{\tau^3}\pi_{33}+\frac{2u_3}{\tau^3}\pi_{03}-\frac{(u_3)^2}{\gamma\tau^3}\pi_{00} \end{align}</math>

最后由于,

<math>\begin{align}

T D_\mu(su^\mu)=T\sigma\gamma d_\tau(\frac{s}{\sigma})=\Pi(-\theta)+\pi^{\mu\nu}\sigma_{\mu\nu}=\Pi(-(\partial\cdot u))+\pi^{\mu\nu}\sigma_{\mu\nu} =\Pi(-(\partial\cdot u))+\pi^{\mu\nu}\sigma_{\mu\nu} \end{align}</math>

为了今后的需要,我们写下

<math>\begin{align}

\frac{ds}{d\tau}=\sigma{d_\tau(\frac{s}{\sigma})}-\frac{s}{\gamma}(\partial\cdot u) \end{align}</math>

我们写下熵流方程最后表达式

<math>\begin{align}

\\d_\tau(\frac{s}{\sigma}) =\frac{1}{T\sigma\gamma}\Pi(-(\partial\cdot u))+\frac{1}{T\sigma\gamma} \left[\left(\pi_0^j-\pi_{00}\frac{u^j}{\gamma}\right)\frac{du_j}{d\tau}+(\pi^{ij}+\pi^{00}v^iv^j-\pi^{i0}v^j-\pi^{0j}v^i)\partial_iu_j-\frac{\gamma}{\tau^3}\pi_{33}+\frac{2u_3}{\tau^3}\pi_{03}-\frac{(u_3)^2}{\gamma\tau^3}\pi_{00}\right] \end{align}</math>

方程清单

我们在这里把所有的运动方程和待求解变量总结如下.体系的待求解变量为

<math>\begin{align}

u^i,\frac{s}{\sigma},\frac{\Pi}{\sigma},\frac{\pi^{\mu\nu}}{\sigma} \end{align}</math>

其中剪切系数的5个独立变量选取为除去33指标的纯空间指标系数,其余系数都可以用这5个系数表达如下

<math>\begin{align}

&\pi^{33}=\left(1-\frac{u_3^2}{\gamma^2}\right)^{-1}\frac{1}{\tau^2}\left[\sum_{i\ne j}\frac{u_iu_j}{\gamma^2}\pi^{ij}-\left(1-\frac{u_1^2}{\gamma^2}\right)\pi^{11}-\left(1-\frac{u_2^2}{\gamma^2}\right)\pi^{22} \right]\\ &\pi^{0i}=-\frac{u_j}{\gamma}\pi^{ji}\\ &\pi^{00}=\pi^{11}+\pi^{22}+\tau^2\pi^{33}=\frac{u_iu_j}{\gamma^2}\pi^{ij} \end{align}</math>

上述待求解变量满足的运动方程为

<math>\begin{align}

d_{\tau}\left(\frac{\Pi}{\sigma}\right) = \frac{d_{\tau}\Pi}{\sigma} + \frac{\Pi}{\sigma^2}\left(\frac{\sigma}{\gamma}\partial \cdot u\right) =\frac{1}{\gamma}\left(-\frac{\widetilde{\Pi}}{\tau_R}-\frac{\zeta\theta}{\sigma\tau_R} \right) =-\frac{1}{\gamma\tau_R\sigma}\left(\Pi+{\zeta\partial\cdot u}\right) \end{align}</math>

<math>\begin{align}

&\frac{\tau_R}{\sigma}(\gamma u_\mu\pi^j_\nu+\gamma u_\nu\pi^j_\mu-u_\mu\pi^0_\nu u^j-u_\nu\pi^0_\mu u^j)\frac{du_j}{d\tau} +\gamma\tau_R\frac{d\widetilde{\pi}_{\mu\nu}}{d\tau}+\frac{\tau_R}{\sigma\tau^3}u_3(\delta_{\mu 0}\pi_{\nu 3}+\delta_{\nu 0}\pi_{\mu 3})\\ &+\frac{\tau_R}{\sigma\tau}u_3(\delta_{\mu 3}\pi_{\nu 0}+\delta_{\nu 3}\pi_{\mu 0}) -\frac{\gamma\tau_R}{\sigma\tau}(\delta_{\mu 3}\pi_{\nu 3}+\delta_{\nu 3}\pi_{\mu 3})+\frac{\pi_{\mu\nu}}{\sigma}\\ &=\frac{\eta}{2\sigma}[\partial_\mu u_\nu+\partial_\nu u_\mu]-\frac{\eta\gamma}{2\sigma}\left[u_\mu\frac{du_\nu}{d\tau}+u_\nu\frac{du_\mu}{d\tau}\right]-\frac{\eta}{3\sigma}\Delta_{\mu\nu}(\partial\cdot u)\\ &-\frac{\eta\gamma\tau}{\sigma}\delta_{\nu 3}\delta_{\mu 3}-\eta\frac{u_3}{\sigma\tau}[\delta_{\mu 0}\delta_{\nu 3}+\delta_{\mu 3}\delta_{\nu 0}]-\frac{\eta(u_3)^2}{2\sigma\tau^3}(u_\mu\delta_{\nu 0}+u_\nu\delta_{\mu 0}) \end{align}</math>

<math>\begin{align}

&\partial_i(P+\Pi)= \gamma(\varepsilon+P+\Pi)\frac{du_i}{d\tau}+u_i(\varepsilon+P)(\partial\cdot u)+u_i{\gamma}\left(1+\frac{dP}{d\varepsilon}\right)\frac{1}{\gamma}\left({T\sigma\gamma}d_\tau\left(\frac{s}{\sigma}\right)-(\varepsilon+P)(\partial\cdot u)\right)\\ &+u_i\gamma\frac{d\widetilde{\Pi}}{d\tau}+\sigma\frac{d\widetilde{\pi}_{0i}}{d\tau}-v_j\partial^j\pi_{0i}+\partial^j\pi_{ji}+\frac{\pi_{0i}}{\tau}-\frac{{\pi}_{0i}}{\gamma}(\partial\cdot u) \end{align}</math>

<math>\begin{align}

\\d_\tau(\frac{s}{\sigma}) =\frac{1}{T\sigma\gamma}\Pi(-(\partial\cdot u))+\frac{1}{T\sigma\gamma} \left[\left(\pi_0^j-\pi_{00}\frac{u^j}{\gamma}\right)\frac{du_j}{d\tau}+(\pi^{ij}+\pi^{00}v^iv^j-\pi^{i0}v^j-\pi^{0j}v^i)\partial_iu_j-\frac{\gamma}{\tau^3}\pi_{33}+\frac{2u_3}{\tau^3}\pi_{03}-\frac{(u_3)^2}{\gamma\tau^3}\pi_{00}\right] \end{align}</math>

方程中的记号涉及到下面关系

<math>\begin{align}

&\theta\equiv Du\equiv\partial \cdot u = \left(u^{\mu}\right)_{;\mu} = \frac{1}{\tau} \partial_{\mu} \left(\tau u^{\mu}\right) =\frac{\gamma}{\tau}+\partial_{\mu}u^{\mu} =\frac{\gamma}{\tau} + d_{\tau}\gamma+\gamma \partial_i v^{i} \\ &d_\tau\sigma=\frac{d\sigma}{d\tau}=-\frac{\sigma}{\gamma}(\partial\cdot u) \frac{du_\mu}{d\tau}=\frac{d\gamma}{d\tau}\delta_{\mu 0}+\frac{du_i}{d\tau}\delta_{\mu i}\\ &\frac{du_\mu}{d\tau}=\frac{d\gamma}{d\tau}\delta_{\mu 0}+\frac{du_i}{d\tau}\delta_{\mu i}\\ &d_{\tau}\gamma = \frac{d}{d\tau}\left(1-g_{ij}u^iu^j\right)^{1/2} = -\frac{g_{ij}u^i}{\gamma}\frac{du^j}{d\tau}-\frac{u^iu^j}{2\gamma}\frac{dg_{ij}}{d\tau} = -\frac{u_j}{\gamma}\frac{du^j}{d\tau}+\frac{\left(u^3\right)^2\tau}{\gamma}\\ &\partial_k\gamma=-v^i\partial_ku_i\\ &u\cdot \partial u_{i} = \gamma \frac{du_i}{d\tau} - u^{\mu} u_{\lambda} \Gamma^{\lambda}_{i\mu} = \gamma \frac{du_i}{d\tau} = - \gamma\frac{du^1}{d\tau}\delta_{i1}- \gamma\frac{du^2}{d\tau}\delta_{i2}- \left(\gamma\tau^2\frac{du^3}{d\tau}+2\tau\gamma u^3 \right)\delta_{i3} \\ &\gamma\frac{du_i}{d\tau}=\gamma\frac{du_1}{d\tau}\delta_{i1}+\gamma\frac{du_2}{d\tau}\delta_{i2}+\gamma\frac{du_3}{d\tau}\delta_{i3}\\ \end{align}</math>

如需要,我们有下面的相互转化关系

<math>\begin{align}

{d_{\tau}\Pi}=-\frac{1}{\gamma\tau_R}\left(\Pi+{\left(\zeta+{\tau_R\Pi} \right)\partial\cdot u}\right) \end{align}</math>

<math>\begin{align}

\frac{d{\pi}_{\mu\nu}}{d\tau}=\sigma\frac{d\widetilde{\pi}_{\mu\nu}}{d\tau}-\frac{{\pi}_{\mu\nu}}{\gamma}(\partial\cdot u) \end{align}</math>

<math>\begin{align}

\frac{ds}{d\tau}=\sigma{d_\tau(\frac{s}{\sigma})}-\frac{s}{\gamma}(\partial\cdot u) \end{align}</math>