通量格式与黎曼问题求解器


李东岳


1. 引言

可压缩求解器由于可能存在间断问题一直令人向往。针对可压缩控制方程,尤其是欧拉方程的求解,存在若干不同的方法。1959年,Godunov提出了基于有限体积法的Godunov方法,其需要对网格单元面处的黎曼问题求精确解,然后获得时间推进后做体平均。但由于Godunov方法需要直接求解每个网格面处的非线性方程,导致计算资源要求极高。另一种方法是对网格单元面处做数值处理,将其演变为线性方程并认为在时间步内是不变的,即近似黎曼求解器。还有一种方法是基于流体分子的传输特性构造矢通量分裂格式。后两种方法的共性都在于求解网格单元面心的通量。区别在于基于的物理角度不同。失通量格式基于分子传输的迎风特性,将网格单元体心的分子运动理解为不同的传输方向(如迎风和下风),进而构建通量。近似黎曼求解器通过数值方法,构造包含网格单元面的控制体,对其求解获得黎曼问题的近似解。本文主要讨论这俩种方法。更详细的介绍可以参考Toro编写的教材。由于目前主流CFD软件均使用非结构网格,因此与教材不同的是,本文基于非结构网格进行推导。

在介绍不同的通量格式之前,有必要简单介绍一下双曲系统。典型的双曲系统即无粘无导热欧拉方程。相对于其他系统(如椭圆或抛物),双曲系统对边界条件的要求最为严格。下面的讨论可以更好的理解双曲系统的特点:考虑普通的低速流,仅仅存在一个真实的传输速度(但在算法里一般认为不可压缩流体的音速是无穷大的),大部分变量通过传输速度进行传递。对于双曲系统,不同的特征变量在不同的特征值下进行传输。因此存在多个近似的传输速度。

2. 近似黎曼求解器

下图为近似黎曼求解器考虑的非结构网格单元,中间的灰色的面表示非结构网格的面。$\bfS_f$表示面矢量。左右两边对应的为own与nei网格单元。在$t$时刻(已知时间步),网格体心处$\psi_\own$与$\psi_\nei$已知。虚线内的体积表示波速$a_\own^f,a_\nei^f$(可以理解为双曲系统的最大/小特征值)的影响区域。由于双曲方程有界的特征值,因此此影响区域$\Delta V'$可求。$\psi_\own^f,\psi_\nei^f$表示自网格中心向面的插值,其定义在红色面以及虚线处。同时,对于own网格单元内红线之外的区域,有体平均值$\psi_\own$,红线内的体平均值为$\psi_\own^f$。nei网格单元同理。

如果打不开图像,请右键在新标签页打开图像后刷新几次

首先,对于$\Delta V$区域,参考《无痛苦N-S方程笔记》第2.4.2节的推导,有下述普适性守恒法则:

\begin{equation}\label{1} \begin{split} \int_{\Delta V}\psi^{t+\dt}\rd V &= \int_{\Delta V}\psi\rd V -\int_t^{t+\dt}\left(\psi_\own^f(t)\bfU_\own^f(t)\cdot\bfS_\own^f+\psi_\nei^f(t)\bfU_\nei^f(t)\cdot\bfS_\nei^f\right)\rd t \end{split} \end{equation}

其中$\bfU$表示速度,$\psi^{t+\dt}$为待求变量,$\psi$为当前时间步的已知变量,$\psi_\own^f(t)\bfU_\own^f(t)$与$\psi_\nei^f(t)\bfU_\nei^f(t)$分别对应红色区域$\Delta V$左右两边面上的变量,本质上其为一个随时间变化的量,但由于此区域位于红色区域之外,因此在$\dt$时间步长内保持恒定。因此方程\eqref{1}可以简化为:

\begin{equation}\label{12} \begin{split} \int_{\Delta V}\psi^{t+\dt}\rd V &=\psi\Delta V -\left(\psi_\own^f\bfU_\own^f\cdot\bfS_\own^f+\psi_\nei^f\bfU_\nei^f\cdot\bfS_\nei^f\right)\Delta t \end{split} \end{equation}
同时有
\begin{equation}\label{2} \begin{split} \psi\Delta V &= \psi_\nei^f(\Delta V_\nei+\Delta V_\nei') + \psi_\own^f(\Delta V_\own+\Delta V_\own') \end{split} \end{equation}
\begin{equation}\label{2S} \bfS_\nei^f\approx\bfS_\own^f\approx\bfS_f, \end{equation}
\begin{equation}\label{1211} \begin{split} \int_{\Delta V}\psi^{t+\dt}\rd V &=\psi_\nei^f(\Delta V_\nei+\Delta V_\nei') + \psi_\own^f(\Delta V_\own+\Delta V_\own') \\\\&-\left(\psi_\own^f\bfU_\own^f\cdot\bfS_{f}+\psi_\nei^f\bfU_\nei^f\cdot\bfS_{f}\right)\Delta t \end{split} \end{equation}

现在考虑$\Delta V'_\own$区域,参考方程\eqref{1},对其做积分有:

\begin{equation}\label{11} \begin{split} \int_{\Delta V'_\own}\psi^{t+\dt}\rd V &= \Delta V'_\own\psi_\own^f -\int_t^{t+\dt}\left(\psi_\own^f\bfU_\own^f\cdot\bfS_f+\psi_f(t)\bfU_f(t)\cdot\bfS_f\right)\rd t \end{split} \end{equation}

由于$\psi_\own^f\bfU_\own^f$在$\dt$时间间隔内保持恒定,有:

\begin{equation}\label{111} \begin{split} \int_{\Delta V'_\own}\psi^{t+\dt}\rd V &= \Delta V'_\own\psi_\own^f -\psi_\own^f\bfU_\own^f\cdot\bfS_f\dt -\int_t^{t+\dt}\psi_f(t)\bfU_f(t)\cdot\bfS_f\rd t \end{split} \end{equation}

其中$\int_t^{t+\dt}\psi_f(t)\bfU_f(t)\cdot\bfS_f\rd t$即为未知的通量函数。方程\eqref{12}与\eqref{111}为普适性守恒法则。不同的方法,在于将方程\eqref{111}进行封闭。

2.1. HLL方法

HLL近似黎曼求解器为Harten、Lax、van Leer于1983年提出的用于求解间断问题的数值方法。在HLL近似黎曼求解器中,将方程\eqref{111}中左侧的$\int_{\Delta V'_\own}\psi^{t+\dt}\rd V$转变为已知项,因此封闭进而通量可求。首先,对方程\eqref{1}左侧进行分组:

\begin{equation}\label{3} \begin{split} \int_{\Delta V}\psi^{t+\dt}\rd V &= \int_{\Delta V_\own}\psi^{t+\dt}\rd V + \int_{\Delta V'}\psi^{t+\dt}\rd V + \int_{\Delta V_\nei}\psi^{t+\dt}\rd V \\\\ &= \psi_\own^f\Delta V_\own +\int_{\Delta V'}\psi^{t+\dt}\rd V +\psi_\nei^f\Delta V_\nei \end{split} \end{equation}
将方程\eqref{3}代入到\eqref{1211}有:
\begin{equation}\label{4} \begin{split} \int_{\Delta V'}\psi^{t+\dt}\rd V &= \psi_\nei^f\Delta V_\nei' + \psi_\own^f\Delta V_\own' \\\\ &-\left(\psi_\own^f\bfU_\own^f\cdot\bfS_{f}+\psi_\nei^f\bfU_\nei^f\cdot\bfS_{f}\right)\Delta t \end{split} \end{equation}
同时有波速影响区域的计算公式:
\begin{equation}\label{5} \begin{split} \Delta V_\nei'=\Delta t|\bfS_f| a_\nei^f, \Delta V_\own'=\Delta t|\bfS_f| a_\own^f \end{split} \end{equation}
将方程\eqref{5}代入到\eqref{4}有:
\begin{equation}\label{6} \begin{split} \int_{\Delta V'}\psi^{t+\dt}\rd V = \left(\psi_\nei^f|\bfS_f| a_\nei^f + \psi_\own^f|\bfS_f| a_\own^f-\psi_\own^f\bfU_\own^f\cdot\bfS_f-\psi_\nei^f\bfU_\nei^f\cdot\bfS_f\right)\Delta t \end{split} \end{equation}

对方程\eqref{6}左右除以$\Delta V'=\Delta V_\nei'+\Delta V_\own'$有:

\begin{equation}\label{7} \begin{split} \frac{1}{\Delta V'}\int_{\Delta V'}\psi^{t+\dt}\rd V &= \frac{\psi_\nei^f|\bfS_f| a_\nei^f + \psi_\own^f|\bfS_f| a_\own^f-\psi_\own^f\bfU_\own^f\cdot\bfS_f-\psi_\nei^f\bfU_\nei^f\cdot\bfS_f} {|\bfS_f| a_\own^f+|\bfS_f| a_\nei^f} \\\\ &= \frac{\psi_\nei^f a_\nei^f + \psi_\own^f a_\own^f-\psi_\own^f\bfU_\own^f\cdot\bfn_f-\psi_\nei^f\bfU_\nei^f\cdot\bfn_f} {a_\own^f+ a_\nei^f} \end{split} \end{equation}

方程\eqref{7}右侧在获得波速$a$之后即可求出,可见,方程\eqref{7}左侧为一个常数。同时,方程\eqref{7}左侧表示在$\psi^{t+\dt}$的值,现将这个值表示为$\psi_{HLL}$,有:

\begin{equation}\label{8} \begin{split} \psi_{HLL} &= \frac{\psi_\nei^f a_\nei^f + \psi_\own^f a_\own^f-\psi_\own^f\bfU_\own^f\cdot\bfn_f-\psi_\nei^f\bfU_\nei^f\cdot\bfn_f} {a_\own^f+ a_\nei^f} \end{split} \end{equation}
$\psi_{HLL}$表示整个$\Delta V'$区域内的$\psi$值,因此,网格面上的值必然也为$\psi_{HLL}$。进一步的,Davis通过下述方程定义$a$:
\begin{equation}\label{10} \begin{split} a_\own^f=\min (\bfU_\own^f\cdot\bfS_f-c_\own^f|\bfS_f|, \bfU_\nei^f\cdot\bfS_f-c_\nei^f|\bfS_f|) \\\\ a_\nei^f=\max (\bfU_\own^f\cdot\bfS_f+c_\own^f|\bfS_f|, \bfU_\nei^f\cdot\bfS_f+c_\nei^f|\bfS_f|) \end{split} \end{equation}
其中$c$为音速。同时,在$\Delta V'_\own$内,有$\psi^{t+\dt}=\psi_{HLL}$,将其代入到方程\eqref{111}中有:
\begin{equation}\label{14} \begin{split} \int_t^{t+\dt}\psi(t)_f\bfU_f(t)\cdot\bfS_f\rd t &= \Delta V'_\own\left(\psi_\own^f-\psi_{HLL}\right) -\psi_\own^f\bfU_\own^f\cdot\bfS_f\dt \end{split} \end{equation}
将方程\eqref{8}代入到方程\eqref{14}中有:
\begin{equation}\label{15} \begin{split} &\frac{1}{\dt}\int_t^{t+\dt}\psi(t)_f\bfU_f(t)\cdot\bfS_f\rd t \\\\ &= \frac{\Delta V'_\own}{\dt}\left(\psi_\own^f-\frac{\psi_\nei^f a_\nei^f + \psi_\own^f a_\own^f-\psi_\own^f\bfU_\own^f\cdot\bfn_f-\psi_\nei^f\bfU_\nei^f\cdot\bfn_f} {a_\own^f+ a_\nei^f}\right) -\psi_\own^f\bfU_\own^f\cdot\bfS_f \\\\ &= \frac{ a_\nei^f a_\own^f(\psi_\own^f |\bfS_f|-\psi_\nei^f|\bfS_f|) +\psi_\nei^f\bfU_\nei^f\cdot\bfS_f a_\own^f -\psi_\own^f\bfU_\own^f\cdot\bfS_f a_\nei^f } {a_\own^f+a_\nei^f} \end{split} \end{equation}
方程\eqref{15}即通量的计算公式。需要注意的是,可能存在特征均向nei或者own网格传递的情况,在这种情况下
\begin{equation}\label{16} F= \left\{\begin{matrix} \psi_\own^f\bfU_\own^f\cdot\bfS_f, & a_\own>0 \\ \frac{ a_\nei^f a_\own^f(\psi_\own^f |\bfS_f|-\psi_\nei^f|\bfS_f|) +\psi_\nei^f\bfU_\nei^f\cdot\bfS_f a_\own^f -\psi_\own^f\bfU_\own^f\cdot\bfS_f a_\nei^f } {a_\own^f+a_\nei^f}, & \mathrm{else} \\ \psi_\nei^f\bfU_\nei^f\cdot\bfS_f, & a_\nei<0 \end{matrix}\right. \end{equation}

方程\eqref{16}即最终的HLL通量。

2.2. Roe方法
Roe方法不同于Hll方法。

3. 矢通量分裂格式
3.1. AUSM方法
AUSM格式全称为Advection Upstream Splitting Method,其为Liou and Steffen于1993年提出的一种通量数值计算方法。不同于上述的HLL近似黎曼求解器,AUSM隶属于一种矢通量分裂方法,其物理特性不基于对间断控制体做积分,而是从气体动理学粒子运动的角度对流向作区分进而组建通量。
如果打不开图像,请右键在新标签页打开图像后刷新几次
对于上述的网格系统,欲求解中间灰色面的通量,矢通量分裂考虑own网格与nei网格内不同的流动方向。其中左侧向右侧流动的$F_\own$,与右侧向左流动的$F_\nei$对中间的面通量有影响。因此,矢通量格式的通量必然是$F_\own$与$F_\nei$的函数。其中定义在$f,\own$的变量表示基于own网格向面$f$进行插值,定义在$f,\nei$的变量表示基于nei网格向面$f$进行插值。在这里可以调用不同的插值格式。例如对于迎风格式,$\psi_{f,own}=\psi_\own$。各种矢通量格式通过不同的方法计算$F_\own$与$F_\nei$。考虑守恒变量的欧拉方程:
\begin{equation}\label{166} \begin{split} \frac{\p \rho}{\p t}+ \nabla\cdot(\rho\bfU)=0 \\\\ \frac{\partial \rho_\mathbf{U}}{\partial t}+\nabla \cdot (\rho_\mathbf{U}\mathbf{U})+\nabla p=0 \\\\ \frac{\p E}{\p t}+\nabla\cdot((E+p)\bfU)=0 \end{split} \end{equation}
其中$\rho$表示密度,$\rho_\bfU$表示动量,$E=1/2\rho\bfU^2+\rho e$,其中$e$表示比内能。有焓值:
\begin{equation}\label{17} H=\frac{E+p}{\rho} \end{equation}
有能量方程:
\begin{equation}\label{18} \frac{\p E}{\p t}+\nabla\cdot(\rho_H\bfU)=0 \end{equation}
现只考虑除了时间项之外的项的离散形式为:
\begin{equation}\label{19} \begin{split} \nabla\cdot(\rho\bfU)&\rightarrow \sum \bfU_f\cdot\bfS_f\rho_f \\\\ \nabla \cdot (\rho_\mathbf{U}\mathbf{U})+\nabla p &\rightarrow \sum \bfU_f\cdot\bfS_f{\rho_\bfU}_f+ \sum p_f\bfS_f \\\\ \nabla\cdot(\rho_H\bfU) &\rightarrow \sum \bfU_f\cdot\bfS_f{\rho_H}_f \end{split} \end{equation}
有面心马赫数的定义为
\begin{equation}\label{20} M_f=\frac{\bfU_f\cdot\bfS_f}{c_f|\bfS_f|}=\frac{\bfU_f\cdot\bfn_f}{c_f} \end{equation}
则方程\eqref{19}可以写为:
\begin{equation}\label{21} \begin{split} \sum \bfU_f\cdot\bfS_f\rho_f&=\sum M_fc_f\rho_f|\bfS_f| \\\\ \sum \bfU_f\cdot\bfS_f{\rho_\bfU}_f+ \sum p_f\bfS_f &=\sum M_fc_f{\rho_\bfU}_f|\bfS_f|+ \sum p_f\bfS_f \\\\ \sum \bfU_f\cdot\bfS_f{\rho_H}_f&=\sum M_fc_f{\rho_H}_f|\bfS_f| \end{split} \end{equation}
方程\eqref{21}中的面变量可以定义为:
\begin{equation}\label{22} \begin{split} c_f\rho_f=\left\{\begin{matrix} \rho_\own^f c_\own^f, & M_f \geqslant 0 \\ \rho_\nei^f c_\nei^f, & M_f < 0 \end{matrix}\right. \\\\ c_f{\rho_\bfU}_f=\left\{\begin{matrix} {\rho_\bfU}_\own^f c_\own^f, & M_f \geqslant 0 \\ {\rho_\bfU}_\nei^f c_\nei^f, & M_f < 0 \end{matrix}\right. \\\\ c_f{\rho_H}_f=\left\{\begin{matrix} {\rho_H}_\own^f c_\own^f, & M_f \geqslant 0 \\ {\rho_H}_\nei^f c_\nei^f, & M_f < 0 \end{matrix}\right. \end{split} \end{equation}
方程\eqref{22}与\eqref{23}相结合即常规的通量计算方法。在AUSM中,面心马赫数通过下式进行计算:
\begin{equation}\label{23} M_f=M_\own^f+M_\nei^f \end{equation}
\begin{equation}\label{24} \begin{split} M_\own^f=\left\{\begin{matrix} 0.25(M_\own^f+1)^2, & |M_\own^f|\leqslant 1\\ 0.5(M_\own^f+|M_\own^f|), & |M_\own^f| > 1 \end{matrix}\right. \\\\ M_\nei^f=\left\{\begin{matrix} -0.25(M_\nei^f-1)^2, & |M_\nei^f|\leqslant 1\\ 0.5(M_\nei^f-|M_\nei^f|), & |M_\nei^f| > 1 \end{matrix}\right. \end{split} \end{equation}
\begin{equation}\label{25} M_\own^f=\frac{\bfU_\own^f\cdot\bfn_f}{c_\own^f}, M_\nei^f=\frac{\bfU_\nei^f\cdot\bfn_f}{c_\nei^f} \end{equation}
进一步的, AUSM认为对流与压力是两种物理性质不同的机制。因此压力与其他项需要分开处理。网格单元面上的压力处理方法为:
\begin{equation}\label{26} p_f=p_\own^f+p_\nei^f \end{equation}
\begin{equation}\label{27} \begin{split} p_\own^f=\left\{\begin{matrix} 0.25p_\own^f(M_\own^f+1)^2(2-M_\own^f), & |M_\own^f|\leqslant 1\\ 0.5p_\own^f\frac{M_\own^f+|M_\own^f|}{M_\own^f}, & |M_\own^f| > 1 \end{matrix}\right. \\\\ p_\nei^f=\left\{\begin{matrix} 0.25p_\nei^f(M_\nei^f-1)^2(2+M_\nei^f), & |M_\nei^f|\leqslant 1\\ 0.5p_\nei^f\frac{M_\nei^f-|M_\nei^f|}{M_\nei^f}, & |M_\nei^f| > 1 \end{matrix}\right. \end{split} \end{equation}

将方程\eqref{21}、\eqref{22}、\eqref{23}、\eqref{26}代入到\eqref{19}有AUSM格式定义的通量。

2.2. AUSM+方法

待更新

东岳流体 2014 - 2020
勘误、讨论、补充内容请前往CFD中文网