9.5.3 Form of the Matrices 矩阵的形式[cfd-9-5-3]

在原始变量\(\vec{W}_p\)的各种选择中,形式en

Among the various choices for the primitive variables \(\vec{W}_p\), the form

\[\vec{W}_p = [p,\ u,\ v,\ w,\ T]^T \tag{9.56}\]

出现得最多。因此,后面的讨论将局限于这一特定形式的\(\vec{W}_p\)。下面,我们将针对一般流体和完全气体,给出变换矩阵和预条件矩阵,以及特征值和左、右特征向量。en

appears most often. Therefore, we shall restrict the further discussion to this particular form of \(\vec{W}_p\). In the following, we will present the transformation and the preconditioning matrices together with the eigenvalues and the left and right eigenvectors for a general fluid, as well as a perfect gas.

Transformation matrices 变换矩阵

对于一般流体,从守恒变量到原始变量的变换矩阵\(\bar{P}^{-1} = \partial\vec{W}_p/\partial\vec{W}\)由文献[91]给出en

For a general fluid, the transformation matrix from the conservative into the primitive variables \(\bar{P}^{-1} = \partial\vec{W}_p/\partial\vec{W}\) is given by [91]

\[\bar{P}^{-1} = \begin{bmatrix} \dfrac{\rho h_T + \rho_T\left(H - q^2\right)}{a_1} & \dfrac{\rho_T u}{a_1} & \dfrac{\rho_T v}{a_1} & \dfrac{\rho_T w}{a_1} & -\dfrac{\rho_T}{a_1} \\[3ex] -\dfrac{u}{\rho} & \dfrac{1}{\rho} & 0 & 0 & 0 \\[3ex] -\dfrac{v}{\rho} & 0 & \dfrac{1}{\rho} & 0 & 0 \\[3ex] -\dfrac{w}{\rho} & 0 & 0 & \dfrac{1}{\rho} & 0 \\[3ex] \dfrac{1 - \rho_p\left(H - q^2\right) - \rho h_p}{a_1} & -\dfrac{\rho_p u}{a_1} & -\dfrac{\rho_p v}{a_1} & -\dfrac{\rho_p w}{a_1} & \dfrac{\rho_p}{a_1} \end{bmatrix} \tag{9.57}\]

其中en

with

\[\begin{aligned} q^2 &= \|\vec{v}\|_2^2 = u^2 + v^2 + w^2\\ a_1 &= \rho\rho_p h_T + \rho_T\left(1 - \rho h_p\right). \end{aligned} \tag{9.58}\]

从原始变量到守恒变量的变换矩阵\(\bar{P} = \partial\vec{W}/\partial\vec{W}_p\)为[91]en

The transformation matrix from the primitive into the conservative variables \(\bar{P} = \partial\vec{W}/\partial\vec{W}_p\) reads [91]

\[\bar{P} = \begin{bmatrix} \rho_p & 0 & 0 & 0 & \rho_T \\ \rho_p u & \rho & 0 & 0 & \rho_T u \\ \rho_p v & 0 & \rho & 0 & \rho_T v \\ \rho_p w & 0 & 0 & \rho & \rho_T w \\ \rho_p H - 1 - \rho h_p & \rho u & \rho v & \rho w & \rho_T H + \rho h_T \end{bmatrix}. \tag{9.59}\]

式(9.57)-(9.59)中密度和焓对压力及温度的导数可以写成en

Derivatives of the density and of the enthalpy with respect to the pressure and the temperature in Eqs. (9.57)-(9.59) can be written in the form

\[\begin{aligned} \rho_p &= \rho\,\alpha_p\\ \rho_T &= -\rho\,\alpha_T\\ h_p &= \frac{1 - \alpha_T T}{\rho}\\ h_T &= c_p\,, \end{aligned} \tag{9.60}\]

其中\(\alpha_p\)和\(\alpha_T\)分别为定压和定温压缩性系数。声速可以由下式计算en

where \(\alpha_p\) and \(\alpha_T\) are the compressibility coefficients at constant pressure and temperature, respectively. The speed of sound can be computed from

\[c^2 = \frac{\rho h_T}{a_1}\,. \tag{9.61}\]

对于完全气体(见2.4.1小节),式(9.60)中的压缩性系数成为\(\alpha_p = 1/p\)和\(\alpha_T = 1/T\)。此时,式(9.57)中从守恒变量到原始变量的变换矩阵\(\bar{P}^{-1} = \partial\vec{W}_p/\partial\vec{W}\)可以写成en

In the case of a perfect gas (see Subsection 2.4.1), the compressibility coefficients in Eq. (9.60) become \(\alpha_p = 1/p\) and \(\alpha_T = 1/T\). In this case, the transformation matrix from the conservative into the primitive variables \(\bar{P}^{-1} = \partial\vec{W}_p/\partial\vec{W}\) from Eq. (9.57) can be cast into

\[\bar{P}^{-1} = \begin{bmatrix} (\gamma-1)\dfrac{q^2}{2} & (1-\gamma)u & (1-\gamma)v & (1-\gamma)w & \gamma-1 \\[3ex] -\dfrac{u}{\rho} & \dfrac{1}{\rho} & 0 & 0 & 0 \\[3ex] -\dfrac{v}{\rho} & 0 & \dfrac{1}{\rho} & 0 & 0 \\[3ex] -\dfrac{w}{\rho} & 0 & 0 & \dfrac{1}{\rho} & 0 \\[3ex] \dfrac{1}{\rho}\left[\dfrac{\gamma q^2}{2c_p} - T\right] & -\dfrac{\gamma u}{c_p\rho} & -\dfrac{\gamma v}{c_p\rho} & -\dfrac{\gamma w}{c_p\rho} & \dfrac{\gamma}{c_p\rho} \end{bmatrix}. \tag{9.62}\]

对完全气体,从原始变量到守恒变量的变换矩阵\(\bar{P} = \partial\vec{W}/\partial\vec{W}_p\)为en

The transformation matrix from the primitive into the conservative variables \(\bar{P} = \partial\vec{W}/\partial\vec{W}_p\) reads for a perfect gas

\[\bar{P} = \begin{bmatrix} \dfrac{\rho}{p} & 0 & 0 & 0 & -\dfrac{\rho}{T} \\[3ex] \dfrac{\rho u}{p} & \rho & 0 & 0 & -\dfrac{\rho u}{T} \\[3ex] \dfrac{\rho v}{p} & 0 & \rho & 0 & -\dfrac{\rho v}{T} \\[3ex] \dfrac{\rho w}{p} & 0 & 0 & \rho & -\dfrac{\rho w}{T} \\[3ex] \dfrac{\rho E}{p} & \rho u & \rho v & \rho w & -\dfrac{\rho q^2}{2T} \end{bmatrix} \tag{9.63}\]

其中\(q^2\)按式(9.58)定义。en

with \(q^2\) as defined in Eq. (9.58).

Preconditioning matrices 预条件矩阵

在欧拉方程的情形下,预条件矩阵的构造相对容易。van Leer等人[92]提出的构造达到了可达的最低条件数。该方法在文献[79]中有详细讨论。然而,Navier-Stokes方程的预条件则更为复杂,原因在于黏性项会导致复数波速,使预条件系统难以分析。黏性流动最著名的预条件器分别由Choi和Merkle[93]、[94],Turkel[90]、[95]-[97],Lee和van Leer[98]、[99]以及Lee[80],Jorgenson和Pletcher[100],以及Weiss和Smith[101]、[45]提出。应用实例见文献[91]、[102]-[109]。en

The construction of a preconditioning matrix is relatively easy in the case of the Euler equations. The formulation proposed by van Leer at al. [92] achieves the lowest attainable condition number. The methodology was discussed in detail in [79]. However, the preconditioning of the Navier-Stokes equations is more involved. The reason is that the viscous terms lead to complex wave speeds, which makes the preconditioned system difficult to analyse. The most recognised preconditioners for viscous flows were proposed by Choi and Merkle [93], [94], Turkel [90], [95]-[97], Lee and van Leer [98], [99] and Lee [80], Jorgenson and Pletcher [100], and by Weiss and Smith [101], [45], respectively. Examples of applications can be found in Refs. [91], [102]-[109].

Weiss and Smith Preconditioner Weiss和Smith预条件器

Weiss和Smith[101]、[45]提出的预条件矩阵\(\bar{\Gamma}\),在一般流体情形下,其形式与式(9.59)中的\(\bar{P}\)相同,只是把\(\rho_p\)替换为适当的预条件参数(preconditioning parameter)\(\theta\)。预条件矩阵的逆也是如此,即\(\bar{\Gamma}^{-1}\)与式(9.57)中的\(\bar{P}^{-1}\)类似。因此,式(9.58)中的参数\(a_1\)变为en

The preconditioning matrix \(\bar{\Gamma}\) due to Weiss and Smith [101], [45] has, in the case of a general fluid, a form identical to \(\bar{P}\) in Eq. (9.59) with \(\rho_p\) replaced by a suitable preconditioning parameter \(\theta\). The same holds also for the inverse of the preconditioning matrix, i.e., \(\bar{\Gamma}^{-1}\) which resembles \(\bar{P}^{-1}\) from Eq. (9.57). Consequently, the parameter \(a_1\) from Eq. (9.58) is changed into

\[a_1^{\Gamma} = \rho\,\theta h_T + \rho_T\left(1 - \rho h_p\right). \tag{9.64}\]

对于完全气体,预条件矩阵\(\bar{\Gamma}\)可以写成en

In the case of a perfect gas, the preconditioning matrix \(\bar{\Gamma}\) can be written as

\[\bar{\Gamma} = \begin{bmatrix} \theta & 0 & 0 & 0 & -\dfrac{\rho}{T} \\[3ex] \theta u & \rho & 0 & 0 & -\dfrac{\rho u}{T} \\[3ex] \theta v & 0 & \rho & 0 & -\dfrac{\rho v}{T} \\[3ex] \theta w & 0 & 0 & \rho & -\dfrac{\rho w}{T} \\[3ex] \theta H - 1 & \rho u & \rho v & \rho w & -\dfrac{\rho q^2}{2T} \end{bmatrix}, \tag{9.65}\]

其中\(\theta\)仍是预条件参数,稍后定义。预条件矩阵的逆由下式给出en

where \(\theta\) is again the preconditioning parameter, which will be defined later. The inverse of the preconditioning matrix is given by

\[\bar{\Gamma}^{-1} = \begin{bmatrix} a_2\left[\dfrac{c^2}{\gamma-1} - \left(H - q^2\right)\right] & -a_2 u & -a_2 v & -a_2 w & a_2 \\[3ex] -\dfrac{u}{\rho} & \dfrac{1}{\rho} & 0 & 0 & 0 \\[3ex] -\dfrac{v}{\rho} & 0 & \dfrac{1}{\rho} & 0 & 0 \\[3ex] -\dfrac{w}{\rho} & 0 & 0 & \dfrac{1}{\rho} & 0 \\[3ex] a_3\left[1 - \theta\left(H - q^2\right)\right] & -\theta a_3 u & -\theta a_3 v & -\theta a_3 w & \theta a_3 \end{bmatrix} \tag{9.66}\]

其中采用缩写en

with the abbreviations

\[\begin{aligned} a_2 &= (\gamma-1)\,\phi\\ a_3 &= \frac{(\gamma-1)\,\phi T}{\rho}\\ \phi &= \frac{1}{\theta c^2 - (\gamma-1)}\,. \end{aligned} \tag{9.67}\]

式(9.64)-(9.67)中预条件参数\(\theta\)有若干种取法。文献[45]给出了一种基于参考速度\(u_r\)的定义。对一般流体,它为en

Several choices exist for the preconditioning parameter \(\theta\) in Equations (9.64)-(9.67). One definition, which is based on a reference velocity \(u_r\) was provided in Ref. [45]. It reads for a general fluid

\[\theta = \frac{1}{u_r^2} - \frac{\rho_T}{\rho h_T}\,, \tag{9.68}\]

相应地,对完全气体,en

or, correspondingly for a perfect gas,

\[\theta = \frac{1}{u_r^2} + (\gamma-1)\frac{1}{c^2}\,. \tag{9.69}\]

式(9.68)或式(9.69)中的参考速度\(u_r\)由如下关系式求得en

The reference velocity \(u_r\) in Eq. (9.68) or Eq. (9.69) is obtained from the relation

\[u_r = \min\left[\max\left(\|\vec{v}\|_2,\ \frac{\nu}{\Delta h},\ \frac{\kappa}{\Delta h},\ \epsilon\sqrt{\frac{|\Delta p|}{\rho}}\right),\ c\right], \tag{9.70}\]

其中\(\Delta h\)是控制体尺寸的度量,\(\Delta p\)表示相邻控制体之间的压差,\(\epsilon\)是一个小数(\(\approx 10^{-3}\))。从式(9.70)的定义可以看到,参考速度受局部输运速度的限制。\(\nu/\Delta h\)和\(\kappa/\Delta h\)两项在以扩散或热传导为主的边界层内变得重要。压力项的目的是防止\(u_r\)在驻点处趋于零。当\(u_r = c\)时,\(\bar{\Gamma}\)变得与\(\bar{P}\)完全相同,\(\bar{\Gamma}^{-1}\)则变为\(\bar{P}^{-1}\)。因此,正如所期望的那样,预条件在超声速流动时自动关闭。en

where \(\Delta h\) is a measure of the control volume size, the quantity \(\Delta p\) stands for the pressure difference between the adjacent control volumes and \(\epsilon\) is a small number (\(\approx 10^{-3}\)). As we can see from the definition in Eq. (9.70), the reference velocity is bounded by the local transport velocity. The terms \(\nu/\Delta h\) and \(\kappa/\Delta h\) become important in boundary layers with dominant diffusion or heat conduction. The pressure term is intended to prevent \(u_r\) from vanishing at stagnation points. In the case that \(u_r = c\), \(\bar{\Gamma}\) becomes identical to \(\bar{P}\) and \(\bar{\Gamma}^{-1}\) is converted into \(\bar{P}^{-1}\). Hence, the preconditioning is turned off for a supersonic flow as intended.

另一种针对完全气体的做法是令[104]en

Another possibility, which was devised for a perfect gas, is to set [104]

\[\theta = \frac{1}{\beta\gamma R T} = \frac{1}{\beta c^2}\,. \tag{9.71}\]

式(9.71)中的参数\(\beta\)定义为en

The parameter \(\beta\) in Eq. (9.71) is defined as

\[\beta = \frac{M_r^2}{1 + (\gamma-1)M_r^2} \tag{9.72}\]

其中参考马赫数由下式给出en

with the reference Mach number given by

\[M_r^2 = \max\left[\min(M^2,\, 1),\ M^2_{min}\right]\,, \tag{9.73}\]

而\(M\)为当地马赫数(\(M^2 = \|\vec{v}\|_2^2/c^2\))。按照文献[104],参数\(M^2_{min} = K M_\infty^2\)且\(K \approx 3\)(也有人取\(K = 1\)甚至\(K = 0.15\);其确切取值似乎取决于驻点区或边界层内控制体的数目)。容易验证,当\(M \ge 1\)时,参数\(\beta\)等于\(1/\gamma\),矩阵\(\bar{\Gamma}\)与\(\bar{P}\)完全相同。应当指出,式(9.69)与式(9.71)中\(\theta\)的两种定义是等价的。因此,借助式(9.70),参考马赫数也可以取为\(M_r = u_r/c\)。en

and \(M\) being the local Mach number (\(M^2 = \|\vec{v}\|_2^2/c^2\)). According to Ref. [104], parameter \(M^2_{min} = K M_\infty^2\) and \(K \approx 3\) (others choose \(K = 1\) or even \(K = 0.15\); the exact value seems to depend on the number of control volumes in the stagnation region or inside the boundary layer). As it can be easily verified, when \(M \ge 1\) the parameter \(\beta\) equals to \(1/\gamma\) and the matrix \(\bar{\Gamma}\) becomes identical to \(\bar{P}\). It should be mentioned that both definitions of \(\theta\) in Eq. (9.69) and in Eq. (9.71) are equivalent. Hence, the reference Mach number could be determined as \(M_r = u_r/c\) with the help of Eq. (9.70).

Eigenvalues of the Preconditioned System 预条件系统的特征值

预条件系统即式(9.49)的特征值矩阵,即\(\bar{P}\bar{\Gamma}^{-1}\bar{A}_c = \bar{P}\bar{\Gamma}^{-1}\bar{A}_{c,p}\bar{P}^{-1}\),由下式给出en

The matrix of the eigenvalues of the preconditioned system Eq. (9.49), i.e, \(\bar{P}\bar{\Gamma}^{-1}\bar{A}_c = \bar{P}\bar{\Gamma}^{-1}\bar{A}_{c,p}\bar{P}^{-1}\) is given by

\[\bar{\Lambda}_{c,p} = \begin{bmatrix} V & 0 & 0 & 0 & 0 \\ 0 & V & 0 & 0 & 0 \\ 0 & 0 & V & 0 & 0 \\ 0 & 0 & 0 & \dfrac{a_4+1}{2}V + c' & 0 \\[2ex] 0 & 0 & 0 & 0 & \dfrac{a_4+1}{2}V - c' \end{bmatrix}, \tag{9.74}\]

其中\(V = \vec{v}\cdot\vec{n}\)表示逆变速度,而en

where \(V = \vec{v}\cdot\vec{n}\) represents the contravariant velocity, and

\[c' = \frac{1}{2}\sqrt{V^2\left(a_4-1\right)^2 + 4a_5} \tag{9.75}\]

表示修正后的声速。式(9.74)和式(9.75)中的参数\(a_4\)和\(a_5\)在一般情形下为en

denotes the modified speed of sound. The parameters \(a_4\) and \(a_5\) in Eq. (9.74) and Eq. (9.75) read in the general case

\[\begin{aligned} a_4 &= \frac{a_1}{a_1^{\Gamma}}\\ a_5 &= \frac{\rho h_T}{a_1^{\Gamma}}\,, \end{aligned} \tag{9.76}\]

其中\(a_1\)按式(9.58)定义,\(a_1^{\Gamma}\)按式(9.64)定义。在完全气体以及按式(9.69)定义预条件参数\(\theta\)的情形下,参数成为\(a_4 = \phi\)和\(a_5 = \phi c^2\),\(\phi\)由式(9.67)给出。当采用式(9.71)中\(\theta\)的第二种定义时,式(9.76)中的参数分别简化为\(a_4 = M_r^2\)和\(a_5 = M_r^2 c^2\)。可以看到,当\(|M| \rightarrow 0\)时,\(c' \approx (V/2)\sqrt{5}\),因此各特征值如预期那样被均衡。这样,由条件数即式(9.40)所代表的刚性得以降低(条件数为\(C_N \approx 2.6\)),时间推进或迭代求解过程的收敛性得到极大增强。另一方面,当\(|M| \ge 1\)时,\(a_4 = 1\)、\(c' = c\),从而恢复\(\bar{A}_c\)的特征值。en

where \(a_1\) is defined in Eq. (9.58) and \(a_1^{\Gamma}\) in Eq. (9.64). In the case of a perfect gas and the definition of the preconditioning parameter \(\theta\) according to Eq. (9.69), the parameters become \(a_4 = \phi\) and \(a_5 = \phi c^2\), with \(\phi\) given by Eq. (9.67). When using the second definition of \(\theta\) from Eq. (9.71), the parameters in Eq. (9.76) simplify to \(a_4 = M_r^2\) and \(a_5 = M_r^2 c^2\), respectively. As we can see, \(c' \approx (V/2)\sqrt{5}\) for \(|M| \rightarrow 0\) and hence the eigenvalues become equalised as intended. In this way, the stiffness represented by the condition number Eq. (9.40) is reduced (the condition number is \(C_N \approx 2.6\)) and the convergence of the time-stepping or iterative solution process is dramatically enhanced. On the other hand, \(a_4 = 1\) and \(c' = c\) for \(|M| \ge 1\), and thus the eigenvalues of \(\bar{A}_c\) are recovered.

基于式(9.74),预条件系统的谱半径成为en

Based on Eq. (9.74), the spectral radius of the preconditioned system becomes

\[\hat{\Lambda}_{c,p} = \left[\frac{a_4+1}{2}|V| + c'\right]\Delta S \tag{9.77}\]

其中\(\Delta S\)表示面面积。该表达式用于计算式(6.14)、(6.18)、(6.20)或(6.22)中的时间步长。在中心人工耗散的情形下,它也取代式(4.53)中的谱半径(另见式(5.33)、(9.52)、(9.53))。en

with \(\Delta S\) denoting the face area. This expression is employed to compute the time step in Eq. (6.14), (6.18), (6.20) or (6.22). It also replaces the spectral radius in Eq. (4.53) in the case of the central artificial dissipation (see also Eqs. (5.33), (9.52), (9.53)).

Eigenvectors of the Preconditioned System 预条件系统的特征向量

对流通量雅可比矩阵的左、右特征向量(参见附录A.11节)会因预条件而改变。这一点对于基于特征变量的空间离散很重要,例如Roe上风格式(4.3.3小节)以及4.3.4小节介绍的上风TVD格式。en

The left and right eigenvectors of the convective flux Jacobian (cf. Section A.11) will be changed by the preconditioning. This is of importance for spatial discretisations based on characteristic variables like Roe's upwind scheme (Subsection 4.3.3), or the upwind TVD scheme presented in the Subsection 4.3.4.

式(9.47)中预条件新通量雅可比矩阵\(\Gamma^{-1}\bar{A}_{c,p}\)的左特征向量矩阵可写为[91]en

The matrix of the left eigenvectors of the preconditioned new flux Jacobian \(\Gamma^{-1}\bar{A}_{c,p}\) in Eq. (9.47) can be written as [91]

\[\bar{T}_p^{-1} = \begin{bmatrix} -\dfrac{\rho_T n_x}{a_1^{\Gamma}} & 0 & \dfrac{a_6 n_z}{c'} & -\dfrac{a_6 n_y}{c'} & -\dfrac{a_6 n_x}{T} \\[3ex] -\dfrac{\rho_T n_y}{a_1^{\Gamma}} & -\dfrac{a_6 n_z}{c'} & 0 & \dfrac{a_6 n_x}{c'} & -\dfrac{a_6 n_y}{T} \\[3ex] -\dfrac{\rho_T n_z}{a_1^{\Gamma}} & \dfrac{a_6 n_y}{c'} & -\dfrac{a_6 n_x}{c'} & 0 & -\dfrac{a_6 n_z}{T} \\[3ex] \dfrac{1}{2} + a_7 & \dfrac{a_6 n_x}{2c'} & \dfrac{a_6 n_y}{2c'} & \dfrac{a_6 n_z}{2c'} & 0 \\[3ex] \dfrac{1}{2} - a_7 & -\dfrac{a_6 n_x}{2c'} & -\dfrac{a_6 n_y}{2c'} & -\dfrac{a_6 n_z}{2c'} & 0 \end{bmatrix} \tag{9.78}\]

如式(9.55)所示,\(\Gamma^{-1}\bar{A}_{c,p}\)的右特征向量矩阵乘以\(\bar{\Gamma}\)后,由如下公式给出[91]en

The matrix of the right eigenvectors of \(\Gamma^{-1}\bar{A}_{c,p}\) multiplied by \(\bar{\Gamma}\), as indicated in Eq. (9.55), is given by the formula [91]

\[\bar{\Gamma}\bar{T}_p = \frac{a_1^{\Gamma}}{\rho h_T} \begin{bmatrix} -a_8 n_x & -a_8 n_y & -a_8 n_z & 1 & 1 \\[2ex] -a_8 u n_x & -a_8 u n_y - c' n_z & -a_8 u n_z + c' n_y & u + (a_9 + c')n_x & u + (a_9 - c')n_x \\[2ex] -a_8 v n_x + c' n_z & -a_8 v n_y & -a_8 v n_z - c' n_x & v + (a_9 + c')n_y & v + (a_9 - c')n_y \\[2ex] -a_8 w n_x - c' n_y & -a_8 w n_y + c' n_x & -a_8 w n_z & w + (a_9 + c')n_z & w + (a_9 - c')n_z \\[2ex] a_{11} - a_{10} n_x & a_{12} - a_{10} n_y & a_{13} - a_{10} n_z & H + (a_9 + c')V & H + (a_9 - c')V \end{bmatrix}, \tag{9.79}\]

其中\(a_1^{\Gamma}\)按式(9.64)定义。上两式(9.78)和(9.79)中的其余参数为en

where \(a_1^{\Gamma}\) is defined in Eq. (9.64). Further parameters in the above Eqs. (9.78) and (9.79) read

\[\begin{aligned} a_6 &= \rho a_5\\[1ex] a_7 &= \frac{V\left(a_4-1\right)}{4c'}\\[1ex] a_8 &= \frac{\rho_T T}{\rho}\\[1ex] a_9 &= -\frac{1}{2}V\left(a_4-1\right)\\[1ex] a_{10} &= a_8 H + T h_T\\ a_{11} &= c'\left(v n_z - w n_y\right)\\ a_{12} &= c'\left(w n_x - u n_z\right)\\ a_{13} &= c'\left(u n_y - v n_x\right) \end{aligned} \tag{9.80}\]

其中\(a_4\)、\(a_5\)按式(9.76)定义,\(c'\)按式(9.75)定义。en

with \(a_4\), \(a_5\) defined in Eq. (9.76) and \(c'\) in Eq. (9.75), respectively.

在守恒变量下表述的预条件Roe格式(见式(9.54))需要矩阵\(\bar{P}\bar{\Gamma}^{-1}\bar{A}_c = \bar{P}\bar{\Gamma}^{-1}\bar{A}_{c,p}\bar{P}^{-1}\)的特征向量\(\bar{T}_{c,p}\)和\(\bar{T}_{c,p}^{-1}\),它们可以通过修改上述特征向量式(9.78)和(9.79)得到。可以证明,\(\bar{P}\bar{\Gamma}^{-1}\bar{A}_c\)的右特征向量矩阵的列由\(\bar{P}\vec{x}\)构成,其中\(\vec{x}\)是\(\Gamma^{-1}\bar{A}_{c,p}\)的右特征向量。因此,\(\bar{P}\bar{\Gamma}^{-1}\bar{A}_c\)的右、左特征向量矩阵可以表示为en

The eigenvectors \(\bar{T}_{c,p}\) and \(\bar{T}_{c,p}^{-1}\) of the matrix \(\bar{P}\bar{\Gamma}^{-1}\bar{A}_c = \bar{P}\bar{\Gamma}^{-1}\bar{A}_{c,p}\bar{P}^{-1}\), which are required for the preconditioned Roe scheme formulated in the conservative variables (see Eq. (9.54)), can be obtained by a modification of the above eigenvectors Eqs. (9.78) and (9.79). It can be shown that the columns of the right-eigenvector matrix of \(\bar{P}\bar{\Gamma}^{-1}\bar{A}_c\) are composed of \(\bar{P}\vec{x}\), where \(\vec{x}\) are the right eigenvectors of \(\Gamma^{-1}\bar{A}_{c,p}\). Hence, the matrices of the right and left eigenvectors of \(\bar{P}\bar{\Gamma}^{-1}\bar{A}_c\) can be expressed as

\[\begin{aligned} \bar{T}_{c,p} &= \bar{P}\,\bar{T}_p\\ \bar{T}_{c,p}^{-1} &= \bar{T}_p^{-1}\bar{P}^{-1}\,. \end{aligned} \tag{9.81}\]

把上述关系式(9.81)代入守恒变量下预条件Roe格式的公式即式(9.54),我们得到en

Inserting the above relations (9.81) into the formula for the preconditioned Roe scheme in the conservative variables Eq. (9.54), we obtain

\[\begin{aligned} \left(\vec{F}_c\right)_{I+1/2} = \frac{1}{2}\Big[&\vec{F}_c(\vec{W}_R) + \vec{F}_c(\vec{W}_L)\\ &- \left(\bar{\Gamma}\bar{T}_p\left|\bar{\Lambda}_{c,p}\right|\bar{T}_p^{-1}\bar{P}^{-1}\right)_{I+1/2}\left(\vec{W}_R - \vec{W}_L\right)\Big] \end{aligned} \tag{9.82}\]

其中\(\bar{\Gamma}\bar{T}_p\)由式(9.79)给出,\(\bar{T}_p^{-1}\)由式(9.78)给出。en

with \(\bar{\Gamma}\bar{T}_p\) given by Eq. (9.79) and \(\bar{T}_p^{-1}\) by Eq. (9.78), respectively.

将式(9.51)的预条件有限体积格式(其中\(\bar{\Gamma}\)按式(9.65)取)应用于翼型流动的例子见图9.10-9.12。可以观察到,对于来流马赫数0.01,中心格式和Roe上风格式都无法给出正确的解。无预条件的格式不能预测压力分布,因而也不能正确预测升力系数(\(C_L = 0.323\)和\(0.324\),而正确值为\(0.352\))。如图9.12所示,预条件有助于获得正确的解(\(C_L = 0.353\)),并且还显著加速了收敛。en

An example of the application of the preconditioned finite-volume scheme from Eq. (9.51) with \(\bar{\Gamma}\) according to Eq. (9.65) to airfoil flow is presented in Figs. 9.10-9.12. As we can observe, both the central scheme as well as Roe's upwind scheme fail to deliver the correct solution for an inflow Mach number of 0.01. The schemes without preconditioning cannot predict the pressure distribution and hence the lift coefficient (\(C_L = 0.323\) and \(0.324\) versus the correct \(0.352\)). As demonstrated in Fig. 9.12, preconditioning helps to obtain the correct solution (\(C_L = 0.353\)), and it also significantly accelerates the convergence.

图9.10:绕对称Joukowsky翼型(厚度10%)的二维无黏流动

图9.10:绕对称Joukowsky翼型(厚度10%)的二维无黏流动。结构网格,\(M_\infty = 10^{-2}\),\(\alpha = 3^{\circ}\),中心空间离散,显式多级时间推进格式,无预条件。上图为收敛历史,下图为压力系数与精确势流解的比较。图例:上图——convergence(收敛残差,左纵轴为\(\log(\text{res})\))与lift(升力系数,右纵轴),横轴为迭代次数(iteration);下图——Euler solver(欧拉解算器,曲线)与exact solution(精确解,圆点),纵轴为\(-C_p\),横轴为\(x/L\)。

图9.11:绕对称Joukowsky翼型(厚度10%)的二维无黏流动

图9.11:绕对称Joukowsky翼型(厚度10%)的二维无黏流动。结构网格,\(M_\infty = 10^{-2}\),\(\alpha = 3^{\circ}\),二阶Roe上风离散,显式多级时间推进格式,无预条件。上图为收敛历史,下图为压力系数与精确势流解的比较。图例:同图9.10。

图9.12:绕对称Joukowsky翼型(厚度10%)的二维无黏流动

图9.12:绕对称Joukowsky翼型(厚度10%)的二维无黏流动。结构网格,\(M_\infty = 10^{-2}\),\(\alpha = 3^{\circ}\),中心空间离散,显式多级时间推进格式,Weiss-Smith预条件。上图为收敛历史,下图为压力系数与精确势流解的比较。图例:同图9.10。