系数矩阵和通量分裂

摘要

Steger Warming 分裂的精粹在于对无粘Jacobian矩阵分裂。而LUSGS推进也需要利用$A^{\pm}$矩阵。
这部分内容从数学上来讲很容易，就是求无粘通量$F$对守恒变量$Q$的Jacobian矩阵$A=\frac{\partial F}{\partial Q}$，然后对此矩阵相似对角化$A=R \Lambda L$，得到左右特征矩阵和对角特征值。对对角特征值进行迎风分裂$\Lambda=\Lambda^+ +\Lambda^-$，最后乘以左右特征矩阵得到迎风分裂的系数矩阵$A^{\pm}=R \Lambda^{\pm} L$。迎风分裂的系数矩阵乘以守恒变量就是SW迎风通量$F^{\pm}=A^{\pm}Q$。

但是推导起来还是复杂的很。思路是这样子的，$A$很容易直接逐项求无粘通量$F$对守恒变量$Q$的偏导数得到，但是该矩阵比较复杂，不容易对角化。曲线救国的方案是求无粘通量$F$对原始变量$W$的偏导数矩阵$M=\frac{\partial F}{\partial W}$,$M$比较容易对角化，而且$M$和$A$是相似的。通过守恒变量$Q$对原始变量$W$的偏导数矩阵联系。求得$M$的对角化后很容易进而求解$A$的对角化（实在是太烦了）。

这里我省去了繁杂的推导过程，直接给出系数矩阵$A$、左右特征矩阵$R$和$L$、迎风分裂的系数矩阵$A^{\pm}$以及SW迎风通量$F^{\pm}$。

公式

三维Euler方程

设无粘三维Euler方程为

$$\frac{\partial Q}{\partial t} + \frac{\partial \tilde{F}_k}{\partial x_k} =0，\quad x_k=x,y,z$$
曲线坐标系中为
$$\frac{\partial JQ}{\partial t} + \frac{\partial J F_k}{\partial x^k} =0,\quad x^k=\xi,\eta,\zeta$$
其中

$$Q = \left[ \begin{array}{c}
{\mathbf{\rho}_{\bf{s}}}\\
\rho {\bf{u}}\\
\rho {e^v}\\
\rho {e^e}
\end{array} \right],
F_k = \left[ \begin{array}{c}
\mathbf{\rho}_\bf{s}\theta \\
\rho \theta \bf{u} + \mathcal{p} \mathbf{k}\\
\rho \theta {e^v}\\
\rho \theta H
\end{array} \right]$$

其中$\mathbf{k}$已经归一化，$\theta=\mathbf{k} \cdot \mathbf{u}$

系数矩阵

系数矩阵$A_k=\frac{F_k}{\partial Q}$具体形式为：

$$A_k = \left( \begin{array}{cccc}
\theta (\delta_{sj} - Y_s) & {Y_s}{\bf{k}}^T & 0 & 0\\
k_x \beta_s - \theta {\bf{u}} & \theta {\bf{I}} - \varphi {\bf{k}}{\bf{u}}^T + {\bf{u}} {\bf{k}}^T & - \varphi {\bf{k}} & \varphi {\bf{k}}\\
{- {e^v}\theta} & e^v {\bf{k}}^T & \theta & 0\\
(\beta_s - H)\theta &H {\bf{k}}^T - \varphi \theta {\bf{u}}^T & -\varphi \theta & (\varphi+1)\theta
\end{array} \right)$$

其中
$\varphi = \frac{\partial p}{\partial (\rho e^e)} =\frac{R_u}{c_V^{tr}\bar M}$
${\beta_s} = \frac{\partial p}{\partial {\rho_s}} = \frac{R_u}{M_s}T - \varphi c_{V,s}^{tr}T + \varphi \left[ \frac{1}{2}\left( u^2 + v^2+ w^2 \right) - e_s^0 \right]$

系数矩阵的对角化

$$A=R \Lambda L$$

$$\Lambda = diag \left[{\bf{\theta}_s} ,\theta ,\theta ,\theta + c,\theta ,\theta - c\right]$$

找到与${\bf{k}} = [k_x,k_y,k_z]^T$垂直的的两个向量${\bf{b}}_i = [b_1,b_2,b_2]^T$与其协变矢量${\bf{b}}^j = [b^1,b^2,b^2]^T$，使 ${\bf{k}}\cdot {\bf{b}}_i = 0,{\bf{b}}_i\cdot {\bf{b}}^j=\delta_i^j$,那么

$$R = \left( \begin{array}
\delta_{sp} & 0 & 0 & Y_s & 0 & Y_s\\
{\bf{u}} & c{\bf{b}}_{\bf{i}}^{(1)} & c{\bf{b}}_{\bf{i}}^{(2)} & {\bf{u}} + c{\bf{k}} & 0 & {\bf{u}} - c{\bf{k}}\\
0 & 0 & 0 & e^v & 1 & e^v\\
{\bf{u}} \cdot {\bf{u}} - \beta_p/ \varphi & c(\bf{u} \cdot \bf{b}_{i}^{(1)}) & c(\bf{u} \cdot \bf{b}_i^{(2)}) & H + \theta c & 1 & H - \theta c
\end{array} \right)$$

$$L_k = \frac{1}{c^2}\left( {\begin{array}{cccc}
(c^2 \delta_{qj} - Y_q \beta_j) & \varphi {\bf{u}}^T Y_q & \varphi Y_q & -\varphi Y_q\\
-c({\bf{b}}^{j(1)} \cdot {\bf{u}}) & c{\bf{b}}^{j(1)T} & 0 & 0\\
-c({\bf{b}}^{j(2)} \cdot {\bf{u}}) & c{\bf{b}}^{j(2)T} & 0 & 0\\
\frac{\beta_j - \theta c}{2} & \frac{c {\bf{k}^T} - \varphi {\bf{u}}^T}{2} & -\frac{\varphi }{2} & \frac{\varphi}{2}\\
-e^v \beta_j & \varphi {\bf{u}} e^v & c^2 + \varphi e^v & -\varphi e^v\\
\frac{\beta_j + \theta c}{2} &-\frac{c {\bf{k}^T} + \varphi {\bf{u}}^T}{2} & -\frac{\varphi }{2} & \frac{\varphi}{2}
\end{array}} \right)$$

系数矩阵的分裂

按照特征值$\lambda$的正负，可以将系数矩阵分裂为
$$A^{\pm}= A^+ +A^-$$
用$\tilde{A}$代表$A^+$或者$A^-$，那么
$$\tilde{A}= R\tilde{\Lambda}L$$
将$\tilde{A}$写成如下形式：
$$\tilde{A}=[\tilde{A}_{1j},\tilde{\bf{A}}_2,\tilde{A}_3,\tilde{A}_4]$$

$$\tilde A_{1,j} = \left[ \begin{array}{c}
{\tilde \lambda_1} {\delta_{sj}} + {Y_s}({J_1}{\beta_j} - {J_2}\theta )\\
(\beta_j {\bf{u}} -\theta c^2+{\bf{k}}) {J_1}+ (\beta_j \bf{k}- \theta\bf{u}) J_2 \\
{e^v}({J_1}{\beta_j} - {J_2}\theta )\\
\theta ({J_2}{\beta_s} - {J_1}\theta {c^2}) + H({J_1}{\beta_j} - {J_2}\theta )
\end{array} \right]$$

$$\tilde{\bf{A}}_2 = \left[ \begin{array}{c}
Y_s ({\bf{k}}^T J_2 - \varphi {\bf{u}}^T J_1)\\
{\bf{u}} ({\bf{k}}^T J_2 - \varphi {\bf{u}}^T J_1) + {\tilde \lambda }_1 {\bf{I}} + {\bf{k}} ({\bf{k}}^T J_1 c^2 - \varphi {\bf{u}}^T J_2)\\
e^v ({\bf{k}}^T J_2 - \varphi {\bf{u}}^T J_1)\\
\theta ({\bf{k}}^T J_1 c^2 - \varphi {\bf{u}}^T J_2) + H( {\bf{k}}^T J_2 - \varphi {\bf{u}}^T J_1)
\end{array} \right]$$

$$\tilde{A}_3 = - \left[ \begin{array}{c}
Y_s\varphi J_1\\
{\bf{u}} \varphi J_1 + {\bf{k}} \varphi J_2\\
e^v \varphi J_1 - {\tilde \lambda }_1\\
\theta \varphi J_2 + H \varphi {J_1}
\end{array} \right]$$

$$\tilde{A}_4 = \left[ \begin{array}{c}
{Y_s} \varphi {J_1}\\
{\bf{u}} \varphi {J_1} + {\bf{k}} \varphi {J_2}\\
{e^v} \varphi {J_1}\\
{\tilde \lambda }_1 + \theta \varphi {J_2} + H\varphi {J_1}
\end{array} \right]$$

其中$J_1=\frac{\tilde\lambda_2+\tilde\lambda_3-2\tilde\lambda_1}{2c^2},J_2=\frac{\tilde\lambda_2-\tilde\lambda_3}{2c}$

正负通量分裂

$$F^{\pm}=A^{\pm}Q=R \Lambda^{\pm} LQ $$

$$F^{\pm} = \left[ \begin{array}{c}
{\bf{\rho}}_{\bf{s}} dd\\
\rho {\bf{u}} dd + \frac{\lambda_2 - \lambda_3}{2\gamma}\rho c{\bf{k}}\\
\rho e^v dd\\
\rho e^e dd + \frac{p}{2\gamma }(\lambda_2 + \lambda_3 - 2\lambda_1) + \frac{\lambda_2 - \lambda_3}{2\gamma} \rho c \theta
\end{array} \right]$$

其中$dd = \frac{2\lambda_1 (\gamma - 1) + \lambda_2 + \lambda_3}{2\gamma}$
作为自洽性的证明，当不对特征值分裂时，该结果退化为通量$F$的表达式。

最后

1. 我现在觉得把守恒变量按照这个顺序排列比较好$Q = \left[ \begin{array}{c}
   {\mathbf{\rho}_{\bf{s}}}\\
   \rho {e^v}\\
   \rho {\bf{u}}\\
   \rho {e^e}
   \end{array} \right]$.这样方便源项偏导数的计算。

2. 本文中的公式实在是太繁琐了。但是由于热化学非平衡流动控制方程中所述的原因，以及周禹的博士论文中公式存在印刷错误，我还是花了大量时间进行了这种似乎没有意义的繁杂矩阵计算。