1.   引　言

计算流体力学(computational fluid dynamics, CFD)是流体力学、计算数学与计算机技术相结合而发展起来的新兴学科, 它通过数值仿真和可视化处理, 对流体流动和热传导等物理现象进行计算机数值分析和研究. 其中有限差分法^[1-3]是应用最广泛、最成的方法之一.

随着CFD的不断发展, 人们所求解的问题也日益复杂化, 这些问题通常归结为高速和低速流动问题, 高速流动需采用可压缩流动计算方法^[4,5], 而低速流动采用不可压缩流动计算方法^[6-8]. 然而, 这两种方法的特性存在很大差别, 例如, 高超声速往往伴随着如激波一类的强间断存在^[9], 亚声速流动的流场中任何微小的扰动都将影响整个流场^[10], 因此具体数值方法要在特定环境下使用. 针对高速的可压缩流动数值求解问题, 学者们已研究出较成熟的理论, 而面对低速流动时, 无论采取显格式还是隐格式计算, 速度越小对计算误差和收敛速度影响越大, 这也被称为刚性问题^[11,12], 这种刚性给数值方法带来了很大的困难.

同时, 不可压缩流动数值方法在实际工程应用又有很大的需求, 如翼型绕流, 边界层内流动^[13-15]. 在某些问题中, 还会出现高低速流动并存的现象^[16-18]. 为了求解此类问题, 需要发展能够有效求解低速流动的数值方法. 为了求解有很强不可压缩性的低速流动问题, Fasel^[19]提出涡函数-流函数法, Aziz^[20]提出势函数-流函数法, 其本质是引入流函数和涡函数并对控制方程进行变换, 但是壁面处涡量难以处理, 而且推广到三维问题后过于复杂. Halrow和Wetch^[21]提出格子中标记点算法(marker in cell, MAC), Hirt等^[22]则在MAC算法上加以改进提出SOLA算法, 此类方法可以有效地处理含有复杂的自由面或界面的流动, 但是对时间步限制非常严格, 所以求解效率较低. Spalding和Patankar^[23]提出压力联系方程的半隐式法(semi implicit method for pressure linked equations, SIMPLE), 该方法核心在于不断修正压力项, 通过修正后的压力求解速度并满足连续方程, 直到收敛. 为了提高收敛速度, Patankar^[24]又提出SIMPLER (SIMPLE revised)算法, van Doormaal和Raithby^[25]提出SIMPLEC (SIMPLE consistent)算法, Minkowycz等^[26]提出SIMPLEX (SIMPLE extraplation)算法, Issa等^[27]提出PISO (pressure implicit with splitting of operators)算法, 他们都通过加强速度与压力的耦合来达到加速收敛. 但是这类方法占用内存很大, 插值过程繁琐, 大量时间都消耗在修正压力项上, 且编程复杂.

实际上, 不可压缩流体是为了简化某些问题提出的理想模型, 任何流体都是可压缩的. 而且很多时候热量的影响是不可忽略的^[28], 质量、动量和能量方程需要耦合在一起求解. 鉴于可压缩流动的数值方法相对成熟, 因此有研究者提出使用计算可压缩流动的数值方法来计算不可压缩流动问题. Chorin^[29]于1967年提出虚拟压缩法就是基于这种思想, 这也成为了预处理方法的雏形, 该方法另一优势是可将高速与低速流动统一在同一算法框架下求解, 避免了切换数值方法的繁琐和交换信息而产生的精度损失. 最早的预处理矩阵由Turkle^[30]于1986年提出, 其核心作用在于改变时间导数项, 减小方程组特征速度的量级差异, 1993年Choi和Merkle^[31]提出与Chorin-Turkel矩阵相似的适用于Euler方程的预处理矩阵. 最早适用于全速域的预处理矩阵是由Viviand^[32]设计的四参数矩阵, 但在驻点附近刚性依然很强. 针对Navier-Stokes方程, Choi与Merkle^[33]以及Buelow 等^[34]提出了考虑Reynolds数的影响的预处理矩阵, Godfrey等^[35]提出的矩阵将Euler处理矩阵结合了离散方程的黏性/无黏性的块Jacobi矩阵, 但仍然难以降低驻点附近的刚性. Weiss等^[36], Pletcher和Chen^[37]提出的矩阵形式简单, 且降低了驻点附近刚性, 被广泛的应用于FLUENT, OVERFLOW等商业软件. Dailey和Pletcher^[38]提另一种矩阵并将其与多重网格法结合. 且预处理方法不仅可以与有限差分法搭配, 也可以与有限体积法^[39]、间断有限元法^[40]搭配, 此外也可以与双时间步技术结合用于计算非定常问题^[41].

本文提出一种新的预处理矩阵, 该矩阵引入了两个关于Mach数的函数b和M_L, 使得方程组的刚性进一步降低, 在确保(甚至能提高)精度和稳定性前提下加速了收敛速度, 并且在高速与低速流动之间过渡更加光滑(其中M_L在有黏和无黏情况下拥有不同的形式). 通过数值实验, 初步证明了本文方法的可靠性, 并与其他方法相对比, 体现本文方法的优势.

2.   基本概念

2.1   物理量的无量纲化

为了方便实现数值实验的相似模拟, 本文对变量采取了如下无量纲化处理:

$$\begin{split} &x=\frac{{x}^{*}}{L},~\;y=\frac{{y}^{*}}{L},~\;z=\frac{{z}^{*}}{L},~\overline{t}=\frac{{t}^{*}}{L/{c}_{\infty }}, \\ &u=\frac{{u}^{*}}{{c}_{\infty }},~\;v=\frac{{v}^{*}}{{c}_{\infty }},\;w=\frac{{w}^{*}}{{c}_{\infty }},~\;g=\frac{{g}^{*}}{{g}_{\infty }} , \\ &p=\frac{{p}^{*}}{{\rho }_{\infty }{c}_{\infty }{}^{2}},\;~\rho =\frac{{\rho }^{*}}{{\rho }_{\infty }},~\;E=\frac{{e}^{*}}{{\rho }_{\infty }{c}_{\infty }{}^{2}}, \\ &c=\frac{{c}^{*}}{{c}_{\infty }},\;~T=\frac{{T}^{*}}{{T}_{\infty }},~\;\mu =\frac{{\mu }^{*}}{{\mu }_{\infty }}, \end{split}$$

其中, x, y, z代表空间直角坐标系的三个坐标分量; L代表参考长度; t代表时间; u, v, w分别代表三个坐标轴方向上的速度; g代表重力加速度; ρ代表密度; E为单位体积流体的总能量; c代表声速; T代表温度; μ为黏性系数; 上标“ * ”代表有量纲量, 下标“ ∞ ”代表参考值; 下文2.2节中的Re为来流Reynolds数, $Re = \rho_\infty c_\infty L/\mu_\infty$.

2.2   流体力学控制方程组

Navier-Stokes方程组是流体动力学中最重要的基本方程, 它是黏性流体微团应用牛顿第二定律得到的运动微分方程. 以三维可压缩流动的情形为例, 忽略外加热源和彻体力的无量纲化方程:

$$\frac{{\partial {{\boldsymbol{Q}}}}}{{\partial t}} + \frac{{\partial {{\boldsymbol{E}}}}}{{\partial x}} + \frac{{\partial {{\boldsymbol{F}}}}}{{\partial y}} + \frac{{\partial {{\boldsymbol{G}}}}}{{\partial z}} = \frac{1}{{Re}}\left( {\frac{{\partial {{{\boldsymbol{E}}}_\nu }}}{{\partial x}} + \frac{{\partial {{{\boldsymbol{F}}}_\nu }}}{{\partial y}} + \frac{{\partial {{{\boldsymbol{G}}}_\nu }}}{{\partial z}}} \right) \text{, }$$

(1)

其中

$$\begin{split} &{{E}}=\frac{p}{\gamma -1}+\frac{\rho {\left|U\right|}^{2}}{2}\text{ }, \\ &{{\boldsymbol{Q}}}={\left(\rho ,\rho u,\rho v,\rho w,E\right)}^{\text{T}}\text{ }, \\ &{{\boldsymbol{E}}}={\left(\rho u,\rho {u}^{2}+p,\rho uv,\rho uw,\left(E+p\right)u\right)}^{\text{T}}\text{ }, \\ &{{\boldsymbol{F}}}={\left(\rho v,\rho uv,\rho {v}^{2}+p,\rho vw,\left(E+p\right)v\right)}^{\text{T}}\text{ }, \\ &{{{\boldsymbol{{G}}}}}={\left(\rho w,\rho uw,\rho vw,\rho {w}^{2}+p,\left(E+p\right)w\right)}^{\text{T}}\text{ }, \\ &{{{\boldsymbol{E}}}}_{\nu }=\left(\begin{array}{c}0\\ {\tau }_{xx}\\ {\tau }_{xy}\\ {\tau }_{xz}\\ u{\tau }_{xx}+v{\tau }_{xy}+w{\tau }_{xz}-{q}_{x}\end{array}\right)\text{ }, \\ &{{{\boldsymbol{F}}}}_{\nu }=\left(\begin{array}{c}0\\ {\tau }_{yx}\\ {\tau }_{yy}\\ {\tau }_{yz}\\ u{\tau }_{yx}+v{\tau }_{yy}+w{\tau }_{yz}-{q}_{y}\end{array}\right)\text{ }, \\ &{{{\boldsymbol{G}}}}_{\nu }=\left(\begin{array}{c}0\\ {\tau }_{zx}\\ {\tau }_{zy}\\ {\tau }_{zz}\\ u{\tau }_{zx}+v{\tau }_{zy}+w{\tau }_{zz}-{q}_{z}\end{array}\right)\text{ }. \end{split}$$

黏性应力项与热传导流通量分别为

$$\begin{split} &{\tau }_{ik}=\mu \left[\left(\frac{\partial {u}_{i}}{\partial {x}_{k}}+\frac{\partial {u}_{k}}{\partial {x}_{i}}\right)-\frac{2}{3}{\delta }_{ik}\frac{\partial {u}_{j}}{\partial {x}_{j}}\right]\text{ }, \\ &{q}_{j}=\frac{-\mu }{\left(\gamma -1\right)M{a}_{\infty }^{2}Pr}\frac{\partial T}{\partial {x}_{j}}\text{ }. \end{split}$$

其中当i = j时, δ_ij = 1, 如果i ≠ j, δ_ij = 0; Pr是Prandtl数; U代表速度矢量, 对于无黏流动, (1)式的右端为零. 在实际的应用中, 经常会用到曲线坐标系形式的方程组, 具体如下:

$$\frac{{\partial {{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{{\boldsymbol{Q}}} }}}}{{\partial \tau }} + \frac{{\partial {{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{{\boldsymbol{E}}} }}}}{{\partial \xi }} + \frac{{\partial {{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{{\boldsymbol{F}}} }}}}{{\partial \eta }} + \frac{{\partial {{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{{\boldsymbol{G}}} }}}}{{\partial \zeta }} = \frac{1}{{Re}}\left( {\frac{{\partial {{{{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{{\boldsymbol{E}}} }}}_\nu }}}{{\partial \xi }} + \frac{{\partial {{{{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{{\boldsymbol{F}}} }}}_\nu }}}{{\partial \eta }} + \frac{{\partial {{{{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{{\boldsymbol{G}}} }}}_\nu }}}{{\partial \zeta }}} \right) .$$

(2)

这里$\left( {\xi , \eta , \zeta } \right)$为曲线坐标系.

$$\begin{split} &J=\frac{\partial \left(\xi ,\eta ,\zeta \right)}{\partial \left(x,y,z\right)},\;{{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{{\boldsymbol{Q}}} }}=\frac{1}{J}{\boldsymbol{Q}}\text{ }, \\ &{{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{{\boldsymbol{E}}} }}=\frac{1}{J}\left({\xi }_{t}{\boldsymbol{Q}}+{\xi }_{x}{\boldsymbol{E}}+{\xi }_{y}{\boldsymbol{F}}+{\xi }_{z}{\boldsymbol{G}}\right)\text{ }, \\ &{{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{{\boldsymbol{F}}} }}=\frac{1}{J}\left({\eta }_{t}{\boldsymbol{Q}}+{\eta }_{x}{\boldsymbol{E}}+{\eta }_{y}{\boldsymbol{F}}+{\eta }_{z}{\boldsymbol{G}}\right)\text{ }, \\ &{{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{{\boldsymbol{G}}} }}=\frac{1}{J}\left({\zeta }_{t}{\boldsymbol{Q}}+{\zeta }_{x}{\boldsymbol{E}}+{\zeta }_{y}{\boldsymbol{F}}+{\zeta }_{z}{\boldsymbol{G}}\right)\text{ }, \\ &{{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{{\boldsymbol{E}}} }}_{\nu }=\frac{1}{J}\left({\xi }_{x}{{\boldsymbol{E}}}_{\nu }+{\xi }_{y}{{\boldsymbol{F}}}_{\nu }+{\xi }_{z}{{\boldsymbol{G}}}_{\nu }\right)\text{ }, \\ &{{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{{\boldsymbol{F}}} }}_{\nu }=\frac{1}{J}\left({\eta }_{x}{{\boldsymbol{E}}}_{\nu }+{\eta }_{y}{{\boldsymbol{F}}}_{\nu }+{\eta }_{z}{{\boldsymbol{G}}}_{\nu }\right)\text{ }, \\ &{{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{{\boldsymbol{G}}} }}_{\nu }=\frac{1}{J}\left({\zeta }_{x}{{\boldsymbol{E}}}_{\nu }+{\zeta }_{y}{{\boldsymbol{F}}}_{\nu }+{\zeta }_{z}{{\boldsymbol{G}}}_{\nu }\right)\text{ }. \end{split}$$

本文采用的是方程组(2), 并取τ = t, 后文将省略该方程组中变量的上标“$\frown$ ”.

3.   方程组的刚性问题

3.1   预处理方法

使用与时间相关的可压缩流动数值方法来求解不可压缩流动时, 方程组的特征根量级差异过大(也称刚性过大), 这会严重影响计算的收敛性. 基于虚拟压缩法原理的预处理法可以有效地解决这一问题, 即在时间导数项前乘上个预处理矩阵, 改变方程组的特征系统. 虽然改变了方程组形式, 但对于定常问题而言, 经过足够长的时间后取得的收敛解不会受此影响. 为了更好地计算黏性问题, 通常在带预处理矩阵的控制方程中使用原始变量${{\boldsymbol{Q}}_T} = {\left( {p, u, v, T} \right)^{\text{T}}}$, 下面以二维带预处理矩阵的有黏性控制方程为例:

$${\boldsymbol{\varGamma}}\frac{{\partial {{\boldsymbol{Q}}_T}}}{{\partial t}} + \frac{{\partial {\boldsymbol{E}}}}{{\partial \xi }} + \frac{{\partial {\boldsymbol{F}}}}{{\partial \eta }} = \frac{1}{{Re}}\left( {\frac{{\partial {{\boldsymbol{E}}_\nu }}}{{\partial \xi }} + \frac{{\partial {{\boldsymbol{F}}_\nu }}}{{\partial \eta }}} \right) \text{, }$$

(3)

其中各变量与通量意义同前文所述. Γ是预处理矩阵, 其形式有多种, 常见的形式有Weiss-Smith矩阵^[36]:

$${{\boldsymbol{\varGamma}}_{{\text{WS}}}} = \left( {\begin{array}{*{20}{c}} \varTheta &0&0&{{\rho _T}} \\ {\varTheta u}&\rho &0&{{\rho _T}u} \\ {\varTheta v}&0&\rho &{{\rho _T}v} \\ {\varTheta H - 1}&{\rho u}&{\rho v}&{{\rho _T}q} \end{array}} \right) \text{, }$$

(4)

其中

$$\begin{split} & q = \frac{{{u^2} + {v^2}}}{2}, ~~ H = \frac{{{c^2}}}{{\gamma - 1}} + q, \hfill \\ \end{split}$$

$$\begin{split} & \varTheta = \frac{1}{{U}_{\text{r}}^{{2}}} + \frac{\gamma -1}{{c}^{2}},\\ & {U}_{\text{r}} = \min\left\{\max\left\{\left|{\boldsymbol{U}}\right|,{k}_{0}\left|{{\boldsymbol{U}}}_{\infty }\right|\right\},c\right\} . \end{split}$$

(5)

k₀通常取1, 如果是有黏方程, 还要限制:

$${U_{\text{r}}} = \max \left\{ {{U_{\text{r}}},\frac{\mu }{{\rho \Delta L}}} \right\} .$$

(6)

Weiss与Smith^[36]构造参考速度U_r, 是为了既能在Ma较大时, 让预处理矩阵还原成非预处理状态, 又能Ma较小时, 减小方程组刚性. 此外, 对于黏性流动而言, 还需要限制U_r不小于局部扩散速度μ/(ρΔL), 这种限制在扩散作用为主导且网格间距较小的情况下是有必要的. 对预处理后的特征根的讨论详见3.2节.

Pletcher-Chen矩阵^[37]:

$${{\boldsymbol{\varGamma}} _{{\text{PC}}}} = \left( {\begin{array}{*{20}{c}} {{\varOmega _{{\text{PC}}}}}&0&0&{{\rho _T}} \\ {{\varOmega _{{\text{PC}}}}u}&\rho &0&{{\rho _T}u} \\ {{\varOmega _{{\text{PC}}}}v}&0&\rho &{{\rho _T}v} \\ {{\varOmega _{{\text{PC}}}}H - 1}&{\rho u}&{\rho v}&{{\rho _T}q} \end{array}} \right){\text{ ,}}$$

(7)

其中

$${\varOmega _{{\text{PC}}}} = \frac{{\gamma \left( {\gamma - 1} \right)M{a^2} + 1}}{{\gamma {{\left| {\boldsymbol{U}} \right|}^2}}} = \frac{1}{{\gamma {{\left| {\boldsymbol{U}} \right|}^2}}} + \frac{{\gamma - 1}}{{{c^2}}} .$$

(8)

Choi-Merkle矩阵^[33]:

$${{\boldsymbol{\varGamma}} _{{\text{CM}}}} = \left( {\begin{array}{*{20}{c}} {\dfrac{1}{{\beta Ma_{\text{r}}^{\text{2}}}}}&0&0&0 \\ {\dfrac{u}{{\beta Ma_{\text{r}}^{\text{2}}}}}&\rho &0&0 \\ {\dfrac{v}{{\beta Ma_{\text{r}}^{\text{2}}}}}&0&\rho &0 \\ {\dfrac{H}{{\beta Ma_{\text{r}}^{\text{2}}}} - 1}&{\rho u}&{\rho v}&{\dfrac{\rho }{{\gamma - 1}}} \end{array}} \right)$$

(9)

其中

$$\begin{split} \beta Ma_{\text{r}}^{\text{2}} =\;& \min \Bigg\{ {c^2},\max \Bigg\{ {c^2}M{a^2}\Bigg( 1 \\ &+ \frac{{1 - M{a_0}^2}}{{M{a_0}^4}}M{a^2} \Bigg),{k_2}^2{{\left| {{U_\infty }} \right|}^2} \Bigg\} \Bigg\}{\text{ }}.\end{split}$$

(10)

Ma₀是参考Mach数, k₂是常数, 通过大量实验发现, Ma₀取0.5, k₂在Euler方程时取0.5, 在NS方程时取1, 这样有利于长时间运算的稳定性. 该方案是Choi与Turkle^[31]提出的, 其作用是为避免驻点附近预处理矩阵的奇异性, 对βMa_r²实施一定的人为控制.

3.2   特征系统分析

如果不预处理, 则方程(3)应该为

$$\frac{{\partial {\boldsymbol{Q}}}}{{\partial {{\boldsymbol{Q}}_T}}}\frac{{\partial {{\boldsymbol{Q}}_T}}}{{\partial t}} + \frac{{\partial {\boldsymbol{E}}}}{{\partial \xi }} + \frac{{\partial {\boldsymbol{F}}}}{{\partial \eta }} = \frac{1}{{Re}}\left( {\frac{{\partial {{\boldsymbol{E}}_\nu }}}{{\partial \xi }} + \frac{{\partial {{\boldsymbol{F}}_\nu }}}{{\partial \eta }}} \right) .$$

(11)

将(11)式等号左端全部替换成以Q_T为变量的形式:

$$\begin{split} &\frac{{\partial {\boldsymbol{Q}}}}{{\partial {{\boldsymbol{Q}}_T}}}\frac{{\partial {{\boldsymbol{Q}}_T}}}{{\partial t}} + \frac{{\partial {\boldsymbol{E}}}}{{\partial {{\boldsymbol{Q}}_T}}}\frac{{\partial {{\boldsymbol{Q}}_T}}}{{\partial \xi }} + \frac{{\partial {\boldsymbol{F}}}}{{\partial {{\boldsymbol{Q}}_T}}}\frac{{\partial {{\boldsymbol{Q}}_T}}}{{\partial \eta }} \\ =\;& \frac{1}{{Re}}\left( {\frac{{\partial {{\boldsymbol{E}}_\nu }}}{{\partial \xi }} + \frac{{\partial {{\boldsymbol{F}}_\nu }}}{{\partial \eta }}} \right) .\end{split}$$

(12)

即

$$\begin{split} &\frac{{\partial {{\boldsymbol{Q}}_T}}}{{\partial t}} + {P_0}^{ - 1}A\frac{{\partial {{\boldsymbol{Q}}_T}}}{{\partial \xi }} + {P_0}^{ - 1}B\frac{{\partial {{\boldsymbol{Q}}_T}}}{{\partial \eta }} \\ =\;& \frac{{{P_0}^{ - 1}}}{{Re}}\left( {\frac{{\partial {{\boldsymbol{E}}_\nu }}}{{\partial \xi }} + \frac{{\partial {{\boldsymbol{F}}_\nu }}}{{\partial \eta }}} \right) , \end{split}\tag{12a}$$

其中守恒变量Q到原始变量Q_T的Jacobi变换矩阵为

$${{\boldsymbol{P}}_0} = \left( {\begin{array}{*{20}{c}} {{\rho _p}}&0&0&{{\rho _T}} \\ {{\rho _p}u}&\rho &0&{{\rho _T}u} \\ {{\rho _p}v}&0&\rho &{{\rho _T}v} \\ {{\rho _p}H - 1}&{\rho u}&{\rho v}&{{\rho _T}q} \end{array}} \right) .$$

(13)

下面以直角坐标系的 ξ 方向为例分析预处理前后方程组的特征系统:

$${\boldsymbol{A}} = \left( {\begin{array}{*{20}{c}} {\dfrac{{\gamma U}}{T}}&{\rho {\xi _x}}&{\rho {\xi _y}}&{ - \dfrac{{\rho U}}{T}} \\ {\dfrac{{\gamma uU}}{T} + {\xi _x}}&{\rho \left( {U + {\xi _x}u} \right)}&{\rho u{\xi _y}}&{ - \dfrac{{\rho uU}}{T}} \\ {\dfrac{{\gamma vU}}{T} + {\xi _y}}&{\rho v{\xi _x}}&{\rho \left( {U + {\xi _y}v} \right)}&{ - \dfrac{{\rho vU}}{T}} \\ {\dfrac{{\gamma UH}}{T}}&{\rho \left( {{\xi _x}H + uU} \right)}&{\rho \left( {{\xi _y}H + vU} \right)}&{ - \dfrac{{\rho Uq}}{T}} \end{array}} \right) \text{, }$$

(14)

其中

$$\begin{split} & U={\xi }_{x}u+{\xi }_{y}v,~~ q= ({{u}^{2}+{v}^{2}})/{2},\\ & H=\frac{\gamma p}{\left(\gamma -1\right)\rho }+q , \end{split}\tag{15a}$$

$$h={\lambda }_{\mathrm{max}}\Delta t,\;{\lambda }_{\mathrm{min}}\Delta t=\frac{h}{{C}_{N}} .\tag{15b}$$

可以求得$P_0^{ - 1}A$的特征根为

$${\lambda }_{1,2}=U,\;{\lambda }_{3,4}=U\pm {C}_{\xi } .$$

(16)

特别地, 当$\xi =x, \eta =y$时有

$${\lambda }_{1,2}=u,\;{\lambda }_{3,4}=u\pm c . \tag{16a}$$

条件数通常用来衡量方程组的刚性大小, 当u → 0时:

$${C_N} = \frac{{\left| {{\lambda _{\max }}} \right|}}{{\left| {{\lambda _{\min }}} \right|}} = \frac{{\left| u \right| + c}}{{\left| u \right|}} = 1 + \frac{1}{{M{a_x}}} .$$

(17)

假如快波传输了1个网格的长度, 则有

$$h={\lambda }_{\mathrm{max}}\Delta t,\;{\lambda }_{\mathrm{min}}\Delta t=\frac{h}{{C}_{N}} .$$

(18)

可见, 当条件数非常大时, 相同时间内快慢波传播的距离会相差很多数量级, 这给方程的求解造成很大困难. 若使用预处理矩阵替换P₀, 则由Choi^[31], Pletcher^[37]和Weiss^[36]等提出的预处理法得到的特征根都有相同的形式:

$$\frac{{1 + f}}{2}u \pm \sqrt {{{\left( {\frac{{1 + f}}{2}u} \right)}^2} + \left( {{c^2} - {u^2}} \right)f} \text{, }\tag{19a}$$

其中f是关于Ma的函数, 虽然每种预处理法都不一样, 但该函数都是改变特征根的关键所在. 以Weiss-Smith形式为例:

$$\begin{split} &{\lambda }_{1,2}=u,\\ &{\lambda }_{3,4}=\frac{\left(1+M{a}_{\text{WS}}{}^{2}\right)u}{2}\\ &\quad\pm \frac{1}{2}\sqrt{{\left(1-M{a}_{\text{WS}}{}^{2}\right)}^{2}{u}^{2} + 4M{a}_{\text{WS}}{}^{2}{c}^{2}} , \end{split}\tag{19b}$$

其中U_r同(5)式:

$$M{a_{{\text{WS}}}} = \frac{{{U_{\text{r}}}}}{c} .$$

(20)

当Mach数很小时, 特征根依然能保持在相近的量级范围内, 刚性得到降低, 条件数的取值范围被限制在1 < C_N ≤ 2内. 为避免驻点附近预处理矩阵的奇异性, Turkle^[30]对预处理参数引入了(10)式的处理方法. 但是, 将经(10)式处理后的特征根看作Ma的函数时, 该函数在Ma₀处连续但不可导, 这会导致当速度接近Ma₀c时, 控制方程从预处理到关闭预处理的过度不够光滑, 可能导致误差放大.

对此, 本文设计了一种新的预处理矩阵, 使得方程组在预处理与非预处理之间的过度更加光滑, 并进一步降低了特征系统的条件数, 从而增加收敛效率, 并且拥有良好的精度和稳定性.

3.3   新预处理矩阵的提出

Turkle^[30]曾在Wiess-Smith矩阵^[36]基础上加以改进, 增加一个参数从而提高预处理矩阵的灵活性, 其矩阵具体如下:

$${{\boldsymbol{\varGamma}} _{{\text{Tur}}}} = \left( {\begin{array}{*{20}{c}} {{\varOmega _{{\text{Tur}}}}}&0&0&{\dfrac{{{C_p}}}{T}\left( {\dfrac{\delta }{{{c^2}}} - \dfrac{\rho }{{{C_p}}}} \right)} \\ {{\varOmega _{{\text{Tur}}}}u}&\rho &0&{u\dfrac{{{C_p}}}{T}\left( {\dfrac{\delta }{{{c^2}}} - \dfrac{\rho }{{{C_p}}}} \right)} \\ {{\varOmega _{{\text{Tur}}}}v}&0&\rho &{v\dfrac{{{C_p}}}{T}\left( {\dfrac{\delta }{{{c^2}}} - \dfrac{\rho }{{{C_p}}}} \right)} \\ {{\varOmega _{{\text{Tur}}}}H - 1}&m&n&{\dfrac{\delta }{{{c^2}}} + \dfrac{{M{a^2}}}{2} \left( {\dfrac{\delta }{{{c^2}}} - \dfrac{\rho }{{{C_p}}}} \right)} \end{array}} \right) \text{, }$$

(21)

式中

$${\varOmega _{{\text{Tur}}}} = \frac{1}{{\beta M{a_{\text{r}}}^2}} - \frac{1}{{\gamma p}}\left( {\frac{\delta }{{{c^2}}} - \frac{\rho }{{{C_p}}}} \right) \text{, }$$

(22)

其中 δ 为根据求解需要选定的0—1之间的常数, βMa_r²同(10)式. 本文受Turkel引入调节参数 δ 的启发, 引入调节函数b₁, b₂和 Ω_L, 这些函数会随着流场中某些物理量的变化而变化, 从而可实现对预处理矩阵动态调控, 具体形式如下:

$${{\boldsymbol{\varGamma}} _{{\text{Liu}}}} = \left( {\begin{array}{*{20}{c}} {{\varOmega _L}}&0&0&{{\rho _T}} \\ {\left( {{\varOmega _L} - \dfrac{{{b_1}}}{{{c^2}{M_L}^2}}} \right)u}&\rho &0&{{\rho _T}u} \\ {\left( {{\varOmega _L} - \dfrac{{{b_2}}}{{{c^2}{M_L}^2}}} \right)v}&0&\rho &{{\rho _T}v} \\ {{\varOmega _L}H - 1}&{\rho u}&{\rho v}&{{\rho _T}\dfrac{{{u^2} + {v^2}}}{2}} \end{array}} \right)$$

(23)

其中

$${\varOmega _L} = \frac{1}{{{c^2}{M_L}^2}} + \frac{{\gamma - 1}}{{{c^2}}} \text{, } \tag{24a}$$

$$M_L =\min\{\max\{Ma,\varepsilon_m\},1\}.\tag{24b}$$

ε_m是防止M_L = 0的一个小量, 本文取ε_m = 10^–6. 若Reynolds数较小, 则需限制其不小于局部扩散速度, 同时既要避免驻点附近的奇异性, 又要保证M_L函数的光滑性, 还要保证后文中伪声速$C'_x$的光滑性, 转化为数学分析语言:

1) M_L函数在(0, 1)上为单调递增的至少C₁连续函数, 且$M_L'$(1) = 0.

2) 当x → 0时, M_L(x) → 0; M_L(1) = 1.

3) 若M_L ≤ μ/(ρc△L), 则M_L = μ/(ρc△L); 若M_L ≥ 1时, 则M_L = 1.

根据以上原则构造出如(24c)式的M_L函数, 其中ε_m同(24b)式:

$${M_L} = \max \left\{ {\min \left\{ {\max \left\{ {Ma\sqrt {1 + M{a^2} - M{a^4}} ,{\varepsilon _m}} \right\},1} \right\},\frac{\mu }{{\rho c\Delta L}}} \right\} . \tag{24c}$$

b₁, b₂是与Ma相关的函数, 可以相等也可以不相等, 在本文中取b₁ = b₂ = b, 具体的有

$$b = 2B - M_L^2 - 1 . \tag{25a}$$

为了在低速时进一步降低方程组刚性, 在高速时还原成原方程, 需使得

$$\mathop {\lim }\limits_{Ma \to 0} \frac{{\left| {{\lambda _{\max }}} \right|}}{{\left| {{\lambda _{\min }}} \right|}} = 1 \text{, }$$

$$Ma \geqslant m{a_0} \gt 1 时,\;\frac{\left|{\lambda }_{\mathrm{max}}\right|}{\left|{\boldsymbol{U}}\right|}=1+\frac{1}{Ma},\; \frac{\left|{\lambda }_{\mathrm{min}}\right|}{\left|{\boldsymbol{U}}\right|}=1-\frac{1}{Ma}$$

故B函数应满足:

1) 当Ma = 0时, B = 0; 当Ma ≥ ma₀时, B = 1.

2) 当0< Ma ≤ma₀时, B是关于Ma 的严格单调递增的至少C₂连续函数.

3) B'(0) = B''(0) = B'(ma₀) = B''(ma₀) = 0.

4) 在全Mach数范围内, 需要满足λ₃ > λ_1,2 > |λ₄|.

5) 在Ma < ma₀时, 要求λ₃/λ_1,2, |λ₄|/λ_1,2分别是关于Ma的单调递增, 递减函数.

6) 在Ma > ma₀时, 要求λ₃/λ_1,2, |λ₄|/λ_1,2分别逼近未经预处理时的λ₃/λ_1,2, |λ₄|/λ_1,2.

7) 使B函数在左右端点Ma = 0和Ma = ma₀附近应保持变化缓慢, 从而抑制速度u在低Mach数和临界点处由于预处理作用而出现的剧变.

8) 与(24c)式的M_L匹配, 使得后文中伪声速$C_x'$具备良好的光滑性.

按照以上8条要求, 本文构造了如(25b)式的B函数:

$$B{\text{ = }}\max \left\{ {\frac{{M{a^{{n_1} + 1}}/m{a_0}}}{{M{a^{{n_1}}} - Ma \cdot m{a_0}^{{n_1} - 2}\left( {m{a_0} - Ma} \right) + {C_L}{{\left( {m{a_0} - Ma} \right)}^{{n_2}}}}},1} \right\} \text{, } \tag{25b}$$

其中ma₀是事先给定的常数; n₁, n₂是正整数; C_L是常数. 为了使B满足上述5)—7)性质, 且尽可能使特征根差异最小, 并使经与处理后的特征根满足λ₃ > λ_1,2 > |λ₄|, 计算无黏和黏性很小的流动时, 选择:

$$n_1=0,\;n_2=4,\;C_L=0.5,\;ma_0=1.73 \tag{26a}$$

如果黏性不可忽略, 可以令

$$n_1=1,\;n_2=3,\;C_L=1,\;ma_0=1.73 \tag{26b}$$

函数B的图像可见图1(d). 通过简单运算可知该矩阵的行列式为

$$\left| {{{\boldsymbol{\varGamma}} _{{\text{Liu}}}}} \right| = \frac{{{\rho ^3}\left[ {1 + b\left( {\gamma - 1} \right)M{a^2}} \right]}}{{\left( {\gamma - 1} \right){c^2}{M_L}^{\text{2}}}} \gt 0 .$$

(27)

可见, 该矩阵恒为可逆矩阵. 此时的预处理矩阵在计算低速问题时效果较好, 如果用于跨声速和全速域问题的计算, 可以选取ma₀ = Ma_max (此最大Mach数是实验前的估计值).

3.4   经新预处理法后的特征系统

3.4.1   特征值和方程组刚性

下面来定量分析该处理法得到的特征系统. 在二维曲线坐标系$\left( {\xi , \eta } \right)$中, 以$\xi$方向为例, 求出经Γ_Liu预处理后Euler方程的特征根:

$${\lambda }_{1,2}=U,\;{\lambda }_{3,4}={U}^{\prime }\pm {C}_{\xi' }{} ,$$

(28)

其中伪速度和伪声速分别为

$$U' = \frac{{{M_L}^{\text{2}} + b + 1}}{2}U = BU ,\tag{29a}$$

$${C_\xi' } = \sqrt {{{\left( {BU} \right)}^2} + {M_L}^2\left( {{C_\xi }^2 - {U^2}} \right)} .\tag{29b}$$

$U, {C_\xi }$意义同(16)式, 预处理参数选用(26a)式, 且U在(0, ma₀U)上时, ${\lambda _3}$为关于Ma的单调递增函数, 与未经预处理时特征根随速度U增加而变化的趋势一致. 经过Γ_Liu预处理后的特征根受两个函数B和M_L调控, 形式上比Weiss, Pletcher和Choi等得到的特征根(19a)式拥有更多自由度. 当Ma → 0时有

$$\underset{Ma\to 0}{\mathrm{lim}}\frac{{\lambda }_{3}}{{\lambda }_{1,2}}=1,\;\underset{Ma\to 0}{\mathrm{lim}}\frac{{\lambda }_{4}}{{\lambda }_{1,2}}=-1 .$$

(30)

此时所有特征根在不考虑正负号的情况下为等价无穷小, 即在Ma很小的时候, 条件数趋于1. 先证明对于任意Ma > ε_m, 有λ₃ > λ_1,2 > |λ₄|.

当Ma≤1时:

$$\begin{split} \;& \frac{{{\lambda _3}}}{{{\lambda _{1,2}}}} = B + \sqrt {{B^2} + {M_L}^2\left( {\frac{1}{{M{a^2}}} - 1} \right)}\\ =\;& B + \sqrt {{B^2} + [ {1 + M{a^2}( {1 - M{a^2}} )} ]( {1 - M{a^2}} )} \gt 1 . \end{split}$$

当1 < Ma ≤ ma₀时:

$$\begin{split} \frac{{{\lambda _3}}}{{{\lambda _{1,2}}}} =\;& B + \sqrt {{B^2} + {M_L}^2\left( {\frac{1}{{M{a^2}}} - 1} \right)} \\ =\;& B + \sqrt {{B^2} + \left( {\frac{1}{{M{a^2}}} - 1} \right)} \gt 1. \end{split}$$

当Ma > ma₀时:

$$\begin{split} \frac{{{\lambda _3}}}{{{\lambda _{1,2}}}} =\;& B + \sqrt {{B^2} + {M_L}^2\left( {\frac{1}{{M{a^2}}} - 1} \right)} \\ =\;& 1 + \sqrt {1 + \left( {\frac{1}{{M{a^2}}} - 1} \right)} \\ = \;& 1 + \frac{1}{{Ma}} \gt 1. \end{split}$$

由此可知对于任意Ma > ε_m有λ₃ > λ_1,2, 同理可证λ_1,2 > |λ₄|, 这也与未经预处理的特征根的大小关系一致, 图1(e)更直观地展示了预处理后的特征根大小关系. 如果$m{a_0} \gt \sqrt 3$则会出现Ma₀>ε_m , 使得λ₃ < λ_1,2, 因此需要限制$m{a_0} \leqslant \sqrt 3$, 故上文选取ma₀ = 1.73. λ₃, λ_1,2恒为正值. 当Ma < 1时, λ₄ < 0; 当Ma > 1时, λ₄ > 0, 与未经预处理时特征根的符号相同. 特殊地, 选取直角坐标系, 可以求出其特征系统的条件数满足:

$$\begin{split}{C_N} =\;& \frac{{\left| {{\lambda _{\max }}} \right|}}{{\left| {{\lambda _{\min }}} \right|}} = \frac{{{\lambda _3}}}{{\left| {{\lambda _4}} \right|}}\\ =\;& \dfrac{{{M_L}^2\left( {\dfrac{1}{{M{a^2}}} - 1} \right)}}{{{{\left[ {\sqrt {{B^2} + {M_L}^2\left( {\dfrac{1}{{M{a^2}}} - 1} \right)} - B} \right]}^2}}} . \end{split}$$

当Ma > ma₀时, 预处理矩阵退化成P₀, 方程恢复成原控制方程. 图1为控制方程经Γ_WS, Γ_PC, Γ_CM, Γ_D以及Γ_Liu预处理后系统的条件数随Ma的发展趋势. 从图1可见, Γ_D在$M{a_\infty }$附近效果良好, 在驻点处刚性很大, 经Γ_Liu处理后的系统刚性小于其他方法, 且特征根随着Mach数从附近增大到Ma > ma₀时的变化过程是光滑的, 当Ma超过ma₀后, 伪速度和伪声速都还原成真实速度和声速. 由于

$$\Delta t = \frac{{{\text{CFL}}\min \left\{ {\Delta x,\Delta y} \right\}}}{{\left| {{\lambda _{\max }}} \right|}} \text{, }$$

(31)

其中CFL数是差分法中判断收敛的条件数, 当CFL数固定时, $\left| {{\lambda _{\max }}} \right|$减小可以增大时间步长, 从而达到加速效果.

3.4.2   特征根的连续性

下面证明特征根${\lambda _i}$关于$\left( {\xi , \eta } \right)$的连续性, 以直角坐标系的x方向为例, 注意到(28)式中$U', {\text{ }}{C_x'}$是关于Ma以1和ma₀为分界点的分段函数, 只需要证明分界点处连续性. 经无量纲化后, 不考虑y方向的速度影响, 则有Ma = U = u, 从而$U', {C_x'}$ 均可视为Ma (或u)的函数, 那么有

$$\left\{ \begin{aligned} & \frac{{\partial U'}}{{\partial x}} = \frac{{{\text{d}}U'}}{{{\text{d}}u}}\frac{{\partial u}}{{\partial x}}{ {,}} \hfill \\ & \frac{{{\partial ^2}U'}}{{\partial {x^2}}} = \frac{{{\text{d}}U'}}{{{\text{d}}u}}\frac{{{\partial ^2}u}}{{\partial {x^2}}} + \frac{{{{\text{d}}^2}U'}}{{{\text{d}}{u^2}}}{\left( {\frac{{\partial u}}{{\partial x}}} \right)^2}{\text{ }}{\text{. }} \end{aligned} \right.$$

(32)

当u存在关于x的二阶连续偏导数时, 只需验证$U'$关于u的可导性. 对(29a)式和(29b)式分别在(0, 1), (1, ma₀), (ma₀, +∞)上关于u求导, 可得

$$\frac{\text{d}{U}^{\prime }}{\text{d}u}=\frac{\text{d}B}{\text{d}u}u+B,\frac{{\text{d}}^{2}{U}^{\prime }}{\text{d}{u}^{2}}=\frac{{\text{d}}^{2}B}{\text{d}{u}^{2}}u+2\frac{\text{d}B}{\text{d}u} .$$

(33)

而B在(0, +∞)上存在关于u的二阶连续导函数, 那么易证:

$$\begin{split} \frac{\text{d}{U}^{\prime }}{\text{d}u}\Big|{}_{u\to {1}^-}=\;&\frac{\text{d}{U}^{\prime }}{\text{d}u} \Big|{}_{u\to {1}^+},\frac{{\text{d}}^{2}{U}^{\prime }}{\text{d}{u}^{2}}\Big|{}_{u\to m{a}_{0}{}^-}\\ =\;&\frac{{\text{d}}^{2}{U}^{\prime }}{\text{d}{u}^{2}}\Big|{}_{u\to m{a}_{0}{}^+} . \end{split}$$

(34)

即$U'$存在关于x的二阶连续偏导数, 同理可证${C_x'}$也存在关于x的二阶连续偏导数. 即经Γ_Liu预处理后得到的, 只要u存在关于x的二阶连续导数, ${\lambda _i}$亦存在关于x的二阶连续导数. 从而克服了Γ_WS, Γ_PC和Γ_CM预处理后的特征根, 在预处理和非预处理间过渡不够光滑的缺陷, 图1(a)—(c)也可反映出这点.

4.   数值求解

4.1   空间离散

采用Roe的通量差分格式作为Riemann求解器^[42], 下面以$\xi$方向介绍该方法, 其他方向与之同理. 带预处理的Roe格式可以表示为

$${{{\boldsymbol{E}}}_{i + \tfrac{1}{2}}} = \frac{1}{2}{\left[ {{{\boldsymbol{E}}}\left( {{{{\boldsymbol{Q}}}^{\text{L}}}} \right) + {{\boldsymbol{E}}}\left( {{{{\boldsymbol{Q}}}^{\text{R}}}} \right) - \tilde A\left( {{{{\boldsymbol{Q}}}^{\text{R}}} - {{{\boldsymbol{Q}}}^{\text{L}}}} \right)} \right]_{i + \tfrac{1}{2}}} ,$$

(35)

其中

$$\tilde A ={\boldsymbol{ \varGamma}} {M_{\boldsymbol{\varGamma}}}\left| {{{{\boldsymbol{\varLambda }}}_{\boldsymbol{\varGamma}}}} \right|M_{\boldsymbol{\varGamma}}^{ - 1}.$$

(36)

${{\boldsymbol{\varLambda}} _{\boldsymbol{\varGamma }}}$是Γ_LiuA的特征根, 即(28)式组成的对角矩阵, A的意义同(14)式, M_Γ是将其化为对角矩阵的变换矩阵, 即

$${{{\boldsymbol{\varLambda}} }_{\boldsymbol{\varGamma}}} = M_{\boldsymbol{\varGamma}}^{ - 1}{{\boldsymbol{\varGamma}}_{{\text{Liu}}}}A{M_{\boldsymbol{\varGamma}}} .$$

(37)

此处$\left| {{{\boldsymbol{\varLambda}} _{\boldsymbol{\varGamma}} }} \right|$表示将该矩阵所有元素取绝对值后组成的矩阵. 高精度的通量E(Q)^L和E(Q)^R采用5阶WENO-JS格式来求得, 此过程可详见文献[43]. 对于黏性项的离散, 本文采取一种4阶紧致的中心差分格式, 具体可见文献[44].

4.2   时间推进

为了与空间离散的精度相匹配, 采取三阶显式Runge-Kutta法做时间推进^[45]:

$$\left\{ \begin{aligned} &{\boldsymbol{Q}}_j^{n\left( 1 \right)} = {\boldsymbol{Q}}_j^n + \Delta t{L_h}\left( {{\boldsymbol{Q}}_j^n} \right){\text{ ,}} \hfill \\ & {\boldsymbol{Q}}_j^{n\left( 2 \right)} = \frac{3}{4}{\boldsymbol{Q}}_j^n + \frac{1}{4}{\boldsymbol{Q}}_j^{n\left( 1 \right)} + \frac{{\Delta t}}{4}{L_h}\left( {{\boldsymbol{Q}}_j^{n\left( 1 \right)}} \right){\text{ ,}} \hfill \\ &{\boldsymbol{Q}}_j^{n + 1} = \frac{1}{3}{\boldsymbol{Q}}_j^n + \frac{2}{3}{\boldsymbol{Q}}_j^{n\left( 2 \right)} + \frac{{2\Delta t}}{3}{L_h}\left( {{\boldsymbol{Q}}_j^{n\left( 2 \right)}} \right){\text{ ,}} \end{aligned} \right.$$

(38)

其中${L_h}\left( {{\boldsymbol{Q}}_j^n} \right)$为${\boldsymbol{Q}}_j^n$关于空间的离散格式.

4.3   非定常问题的时间推进

本文仅是对新预处理法的合理性和可靠性进行检验, 故选取较为简单的定常流动来验证. 如用于非定常流动计算, 预处理会破坏时间精度, 因此需要做特殊处理. 这里介绍一种应用较广泛的双时间步法, 其本质是在控制方程中填入一项关于伪时间τ的导数项, 并对其做预处理, 即

$${\boldsymbol{\varGamma}}\frac{{\partial {{\boldsymbol{Q}}_T}}}{{\partial \tau }} + \frac{{\partial {\boldsymbol{Q}}}}{{\partial t}} + \frac{{\partial {\boldsymbol{E}}}}{{\partial x}} + \frac{{\partial {\boldsymbol{F}}}}{{\partial y}} = \frac{1}{{Re}}\left( {\frac{{\partial {{\boldsymbol{E}}_\nu }}}{{\partial x}} + \frac{{\partial {{\boldsymbol{F}}_\nu }}}{{\partial y}}} \right){\text{ }}{\text{. }}$$

其中 τ 和t分别是伪时间和真实物理时间, 其余物理量同(3)式所述. 当 τ → ∞时, 上式因左边的第一项消失而复原为方程. 分别对 τ 和t离散, 该方法包含两个循环, 即伪时间 τ 的内循环和物理时间t的外循环, 在每个物理时间步当作定常问题用伪时间推进求解, 无需时间精确的伪时间内迭代可以采用预处理加速收敛. 当解的残差变化稳定后, 即达到当前物理时间步的真实解, 随即物理时间t外循环推进至下一步, 然后进行伪时间的内迭代并重复上述步骤, 最终得到欲求时刻的非定常解.

5.   数值模拟验证

下面选取一维有精确解的方程组, 平板层流边界层, 方腔顶盖驱动, 凸鼓包管道流动等定常问题来验证本文方法的可靠性、收敛性和稳定性. 由于Γ_D在远离Ma_∞的效果不如其他方法, 文献[46]的实验说明了使用Weiss & Smith和Choi & Merkle矩阵的效果接近, 故本文选取Weiss & Smith, Pletcher & Chen和本文提出的预处理矩阵进行实验, 并对比得到的结果. 算例5.1采用matlab2019编程, 算例5.2—5.4采用fortran90编程.

5.1   精度测试

选取有定常精确解的一维Euler方程组, 以${{\boldsymbol{Q}}_T} = {\left( {p, u, T} \right)^{\text{T}}}$为变量:

$$\frac{{\partial {{\boldsymbol{Q}}}}}{{\partial t}} + \frac{{\partial {{\boldsymbol{Q}}}}}{{\partial {{\boldsymbol{W}}}}}\frac{{\partial {{\boldsymbol{E}}}}}{{\partial x}} = {{\boldsymbol{S}}} \text{, }$$

(39)

$$\begin{split} &{{\boldsymbol{Q}}}={\left(p,u,T\right)}^{\text{T}},\;{{\boldsymbol{W}}}={\left(\rho ,\rho u,E\right)}^{\text{T}},\\ &{{\boldsymbol{E}}}={\left(\rho u,\rho {u}^{2}+p,\left(\rho +E\right)u\right)}^{\text{T}},\;{{\boldsymbol{S}}}={\left({s}_{1},{s}_{2},{s}_{3}\right)}^{\text{T}},\end{split}$$

$$\left\{ \begin{aligned} & p = \left( {1 + 0.2{{\text{e}}^{ - t}}} \right)/\gamma,\hfill \\ &u = \left( {1 + k - {{\text{e}}^{ - t}}} \right)\left( {0.5 + kx} \right), \hfill \\ & T = 1, \end{aligned} \right.$$

(40)

$$\left\{\begin{aligned} &{s}_{1}=k\left(\frac{{\text{e}}^{-t}}{5}+1\right)\left(k-{\text{e}}^{-t}+1\right)-\frac{{\text{e}}^{-t}}{5}\text{ }\text{, }\\ &{s}_{2}=\left(kx+\frac{1}{2}\right)\left\{\left(k-{\text{e}}^{-t}+1\right)+\left(\frac{{\text{e}}^{-t}}{5}+1\right)\left[{\text{e}}^{-t}+2k\left(k-{\text{e}}^{-t}+1\right)\left(k-{\text{e}}^{-t}\right)\right]\right\},\\ &{s}_{3}=k\left(k-{\text{e}}^{-t}+1\right)\left\{{\left(kx+\frac{1}{2}\right)}^{2}\left[\left(k+\frac{1}{2}\right){\left(k-{\text{e}}^{-t}+1\right)}^{2}+\frac{{\text{e}}^{-t}}{k}\right]+\frac{1}{\gamma \left(\gamma -1\right)}+\frac{\left( {{\text{e}}^{-t}}/{5}+1\right)}{{\gamma}}\right\}. \end{aligned}\right.$$

(41)

其计算域为0 < x < 1, 按照精确解给出源项及边界条件, 设网格数N = 100, 选取CFL = 1. 注意到t足够大后, p, u和T在固定的x处均趋于常数, 可认为是定常问题. 选取k = 0.5, 0.04, 0.01, 0.005, 0.001, 分别使用本文预处理法和无预处理来计算该问题, 并计算到t = 100, 150, 200时刻, 得到的结果与精确解对比如表1所列, 预处理的效果与未进行预处理的对比如图2所示.

由表1可见, 当k = 0.5时, 速度u会经历从小于1 (无量纲化)到大于ma₀的跨越, 使用本文的预处理法可以在保证精度的同时实现从预处理还原到不处理的连续过渡; 在相同的计算时间t = 100下, 使用预处理方法比未处理得到的数值解精度高出约1—2个量级. 增加一倍的计算时间, 到t = 200时刻未用预处理的解的精度才勉强达到t = 100时刻使用预处理的精度. 当k进一步减小时, 预处理的效果更加明显, 从图2(c)—(f)可以清晰地看出, 预处理比不处理收敛速度更快. 这表明本文的预处理方法能够在保证精度的同时, 提高计算的效率.

5.2   平板层流边界层

一长为L的平板静置在水平面上, 从左侧的与平板平行的均匀来流通过其上表面, 以板长为参考长度, 以来流处的声速为参考速度进行无量纲化, 给定来流速度u_∞ = 0.1, 基于板长的Reynolds数Re_∞ = 10⁶, 设定来流处压强和温度以及出口的压强, 计算网格为180×50, 其中x方向为均匀网格, y方向随着y减小而加密, 贴近平板上表面的第1层网格厚度为0.002, CFL = 5. 计算到收敛后, 选取x = 0.8处的速度u随y的分布, 并与Bulsius准解析解对比, 对比结果如图3所示, 其中

$$R{e}_{x}=\frac{R{e}_{\infty }x}{L},\;\eta =\frac{y\sqrt{R{e}_{x}}}{x} .$$

(42)

本文预处理法的收敛速度同Weiss & Smith和Pletcher & Chen方法的对比如图4所示. 可见, 本文的预处理方法得到的数值解与Blasius解基本一致. 由图4可见, 本文方法的收敛效率略好于Weiss & Smith和Pletcher & Chen方法.

5.3   空腔顶盖驱动

此算例描述的是一边长为1 (此算例中变量均无量纲化)的正方形空腔, 初始Mach数Ma₀ = 1, 顶盖以u₀ = 1速度向右运动, Reynolds数分别为100, 1000, 10000. 计算网格为60×60成中心对称的非均匀网格, 并在壁面附近加密, CFL = 100, 分别使用本文预处理法和Weiss & Smith和Pletcher & Chen方法计算此问题. 图5—图7是计算到收敛后得到的空腔内部流场线图, 可以观察到, 在空腔的几何中心附近产生一个大涡, 在底部两角附近产生两个较小的涡结构. 随着Reynolds数的增大, 大涡中心逐渐向空腔几何中心移动, 底部小涡逐渐增大, 当Reynolds数为10000时, 左上角附近也产生一个涡结构. 图8是本文水平与垂直中心线上速度与使用多重网格法的计算结果^[47]的对比, 图9是三种预处理法的收敛速度对比.

通过图8和图9可知, 本文提出的方法精确度略高于Weiss & Smith和Pletcher & Chen方法(这两种方法的精度对比可参考文献[46]). 本文方法的收敛速度也更快一些, 且随着黏性效应增加(Reynolds数减小)效果越明显, 这表明本文格式可以实现对低速有黏问题的高精度高效率的模拟.

5.4   带凸鼓包的二维管道流动

此算例^[48]目的在于验证本文提出的预处理方法在曲线坐标系下的计算性能, 故选择下表面带鼓包的二维管道流动问题, 该算例模型和流场相对简单, 便于对结果进行分析. 具体模型如下: 管道的长度为3 (单位均无量纲化), 高度为1, 凸鼓包位于下表面的中心位置, 长度为1, 其尺寸为y = 0.1sin²πx (1 < x < 2). 计算网格选取160 × 60, 并在下表面处加密, 具体见图10. 来流Mach数选取为Ma₀ = 0.01, 0.05, 0.1和0.2, 设p₀ = 1, T₀ = 1, 取CFL = 5, 设置出口反压与来流值相等, 整个流场无黏, 故在壁面处的法向速度为0.

图11为三种方法的残差收敛曲线对比. Ma₀ = 0.05时三种方法得到的流场Mach数等值线分布见图12(a)—(c). 图12(d)为在加密网格(600×200)下使用Weiss & Smith法计算得到的流场, 并将其作为准精确解. 在低速流动下, 该流场的Mach等值线接近轴对称. 可见本文方法的计算结果与准精确解基本一致, 相比另外两种方法精度更高, 尤其第3—5条等值线的位置. 另外本文方法的收敛效率也略高于其他方法, 虽然Ma₀ = 0.2时由于刚性并不是很大, 加速收敛的效果比Ma₀更低时明显, 但最终的残差均小于其他预处理法, 这说明本文格式在保证精度的同时, 拥有更良好的稳定性.

6.   结　论

针对低速流动问题, 本文基于虚拟压缩法的思想, 提出一种以当地Mach数、来流Mach数、密度、温度、每个维度的分速度和Reynolds数为变量的预处理矩阵, 相比Weiss, Choi, Dailey, Pletcher和Turkle等提出的矩阵, 把控制方程的刚性在Ma<0.5时进一步降低到1左右, 使得特征根之间的量级差异进一步减小, 从而提升了求解效率, 同时改进了从预处理到关闭预处理的过渡不够光滑的不足.

展开数值实验验证了新方法的可靠性. 通过一维算例验证了新方法的准确性, 将预处理和不处理的结果加以对比, 证明了预处理方法可加速计算低速流动问题的收敛. 通过平板层流边界层和黏性方腔流动问题进一步验证了本文方法在二维低速流动问题中可以获得高精度数值解, 并且拥有更高效的收敛性, 随着黏性作用越强, 这种高效性越明显. 通过二维鼓包管道流动问题, 验证了本文方法在曲线坐标系下仍能保持良好的精度和收敛性, 而且具备更好的稳定性. 目前本文重点在于改进算法, 故选取算例的结构比较简单, 未来将会尝试模拟复杂网格, 非定常流动或引入湍流模型.

感谢中国科学院力学研究所高温气体动力学国家重点实验室的李诗尧博士, 天津大学水利工程仿真与安全国家重点实验室的陈嘉禹硕士为本文提供的计算资源; 感谢天津大学数学学院2018级程启豪同学, 南开大学数学学院学习部2018级左晨琛、2019级的戚红利、2020级赵振岐等同学为本文的矩阵特征系统在曲线坐标系下数学性质以及特定几何结构从笛卡尔坐标系下到曲线坐标下变换涉及的复分析和拓扑问题的讨论; 感谢国防科技大学空天科学学院田正雨教授, 南京航空航天大学航空学院陈红全教授, 南京航空航天大学数学系张燕博士为本文提出建议.