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


李东岳


1. 引言

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

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

2. 近似黎曼求解器
2.1. HLL方法
HLL近似黎曼求解器为Harten、Lax、van Leer于1983年提出的用于求解间断问题的数值方法。求解任意一个方程,都可以通过高斯积分将对流项/扩散项转换为变量通量加和的形式进行计算。在HLL近似黎曼求解器中,通过对包含网格面的控制体做积分,可以直接获得网格面的通量值。获得通量之后即可对控制方程进行离散求解。下图为HLL近似黎曼求解器考虑的非结构网格单元,中间的灰色的面表示非结构网格的面。$\bfS_f$表示面矢量。左右两边对应的为own与nei网格单元。在$t$时刻(已知时间步),own网格内的值全部为$\psi_\own$,nei网格内的值全部为$\psi_\nei$。红色虚线表示波速$a_\own,a_\nei$(可以理解为双曲系统的最大/小特征值)的影响位置。红色虚线区域$\Delta V'$表示在$t+\dt$时刻受影响的区域。$\Delta V'$在获得$a_\own,a_\nei$后即可求出。红色虚线外的变量未受影响,因此$\Delta V_\nei$体积内依然为已知时间步的$\psi_\nei$,$\Delta V_\own$体积内依然为已知时间步的$\psi_\own$。$\Delta V$为考虑的积分区域。需要注意的是,$\Delta V$只要大于$\Delta V'$,其选取并不影响结果,最大的$\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^{t}\rd V -\int_t^{t+\dt}\left(\psi_\own(t)\bfU_\own(t)\cdot\bfS_f+\psi_\nei(t)\bfU_\nei(t)\cdot\bfS_f\right)\rd t \end{split} \end{equation}
其中$\bfS_f$表示灰色的面矢量,$\bfU$表示速度,$\psi^{t+\dt}$为待求变量,$\psi^t$为当前时间步的已知变量,$\psi_\own^t\bfU_\own^t$与$\psi_\nei^t\bfU_\nei^t$分别对应绿色的面的变量。此区域位于红色区域之外,因此在$\dt$时间步长内保持恒定。进一步整理方程\eqref{1}有:
\begin{equation}\label{12} \begin{split} \int_{\Delta V}\psi^{t+\dt}\rd V &=\int_{\Delta V}\psi^{t}\rd V -\left(\psi_\own^{t}\bfU_\own^t\cdot\bfS_f+\psi_\nei^{t}\bfU_\nei^t\cdot\bfS_f\right)\Delta t \end{split} \end{equation}
其中
\begin{equation}\label{2} \begin{split} \int_{\Delta V}\psi^{t}\rd V = \psi_\nei^{t}(\Delta V_\nei+\Delta V_\nei') + \psi_\own^{t}(\Delta V_\own+\Delta V_\own') \end{split} \end{equation}
对方程\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^{t}\Delta V_\own +\int_{\Delta V'}\psi^{t+\dt}\rd V +\psi_\nei^{t}\Delta V_\nei \end{split} \end{equation}
将方程\eqref{3},\eqref{2}代入到\eqref{12}有:
\begin{equation}\label{4} \begin{split} \int_{\Delta V'}\psi^{t+\dt}\rd V &= \psi_\nei^{t}\Delta V_\nei' + \psi_\own^{t}\Delta V_\own' \\\\ &-\left(\psi_\own^{t}\bfU_\own^t\cdot\bfS_f+\psi_\nei^{t}\bfU_\nei^t\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, \Delta V_\own'=\Delta t|\bfS_f| a_\own \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^{t}|\bfS_f| a_\nei + \psi_\own^{t}|\bfS_f| a_\own-\psi_\own^{t}\bfU_\own^t\cdot\bfS_f-\psi_\nei^{t}\bfU_\nei^t\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^{t}|\bfS_f| a_\nei + \psi_\own^{t}|\bfS_f| a_\own-\psi_\own^{t}\bfU_\own^t\cdot\bfS_f-\psi_\nei^{t}\bfU_\nei^t\cdot\bfS_f} {|\bfS_f| a_\own+|\bfS_f| a_\nei} \\\\ &= \frac{\psi_\nei^{t} a_\nei + \psi_\own^{t}a_\own-\psi_\own^{t}\bfU_\own^t\cdot\bfn_f-\psi_\nei^{t}\bfU_\nei^t\cdot\bfn_f} {a_\own+a_\nei} \end{split} \end{equation}
方程\eqref{7}右侧在获得波速$a$之后即可求出,可见,方程\eqref{7}左侧为一个常数。同时,方程\eqref{7}左侧表示在$t+\dt$时刻下的$\psi$值,现将这个值表示为$\psi_{HLL}$,有:
\begin{equation}\label{8} \begin{split} \psi_{HLL} &= \frac{\psi_\nei^{t} a_\nei + \psi_\own^{t}a_\own-\psi_\own^{t}\bfU_\own^t\cdot\bfn_f-\psi_\nei^{t}\bfU_\nei^t\cdot\bfn_f} {a_\own+a_\nei} \end{split} \end{equation}
$\psi_{HLL}$表示整个$\Delta V'$区域内的$\psi$值,因此,网格面上的值必然也为$\psi_{HLL}$。进一步的,Davis通过下述方程定义$a$:
\begin{equation}\label{10} \begin{split} a_\own=\min (\bfU_\own^t\cdot\bfS_f-c_\own|\bfS_f|, \bfU_\nei^t\cdot\bfS_f-c_\nei|\bfS_f|), \\\\ a_\nei=\max (\bfU_\own^t\cdot\bfS_f+c_\own|\bfS_f|, \bfU_\nei^t\cdot\bfS_f+c_\nei|\bfS_f|) \end{split} \end{equation}
其中$c$为音速。

至此,面$f$上的$\psi^{t+\dt}$已经可以通过当前时间步的值获得。但$\bfU^{t+\dt}$依然未知,因此通量并不封闭。现在考虑$\Delta V'_\own$区域,参考方程\eqref{1},对其做积分有:
\begin{equation}\label{11} \begin{split} \int_{\Delta V'_\own}\psi^{t+\dt}\rd V &= \int_{\Delta V'_\own}\psi^{t}\rd V -\int_t^{t+\dt}\left(\psi_\own\bfU_\own\cdot\bfS_f+\psi_f(t)\bfU_f(t)\cdot\bfS_f\right)\rd t \\\\ &= \Delta V'_\own\psi^{t}_\own -\int_t^{t+\dt}\left(\psi_\own\bfU_\own\cdot\bfS_f+\psi_f(t)\bfU_f(t)\cdot\bfS_f\right)\rd t \end{split} \end{equation}
其中$\psi_\own\bfU_\own$在$\dt$时间间隔内保持恒定,有:
\begin{equation}\label{123} \begin{split} \int_t^{t+\dt}\psi_f(t)\bfU_f(t)\cdot\bfS_f\rd t &= \Delta V'_\own\psi^{t}_\own -\psi_\own^t\bfU_\own^t\cdot\bfS_f\dt - \int_{\Delta V'_\own}\psi^{t+\dt}\rd V \end{split} \end{equation}
同时,在$\Delta V'_\own$内,有$\psi^{t+\dt}=\psi^{t+\dt}_{HLL}$,代入到方程\eqref{123}中有:
\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^{t}_\own-\psi^{t+\dt}_{HLL}\right) -\psi_\own^t\bfU_\own^t\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^{t}_\own-\frac{\psi_\nei^{t} a_\nei + \psi_\own^{t}a_\own-\psi_\own^{t}\bfU_\own\cdot\bfn_f-\psi_\nei^{t}\bfU_\nei\cdot\bfn_f} {a_\own+a_\nei}\right) -\psi_\own^t\bfU_\own^t\cdot\bfS_f \\\\ &= \frac{ a_\nei a_\own(\psi^{t}_\own |\bfS_f|-\psi_\nei^{t}|\bfS_f|) +\psi_\nei^{t}\bfU_\nei\cdot\bfS_f a_\own -\psi_\own^t\bfU_\own^t\cdot\bfS_f a_\nei } {a_\own+a_\nei} \end{split} \end{equation}
方程\eqref{15}即通量的计算公式。需要注意的是,可能存在特征均向nei或者own网格传递的情况,在这种情况下
\begin{equation}\label{16} F= \left\{\begin{matrix} \psi^t_\own\bfU_\own\cdot\bfS_f, & a_\own>0 \\ \frac{ a_\nei a_\own(\psi^{t}_\own |\bfS_f|-\psi_\nei^{t}|\bfS_f|) +\psi_\nei^{t}\bfU_\nei\cdot\bfS_f a_\own -\psi_\own^t\bfU_\own^t\cdot\bfS_f a_\nei } {a_\own+a_\nei}, & \mathrm{else} \\ \psi^t_\nei\bfU_\nei\cdot\bfS_f, & a_\nei<0 \end{matrix}\right. \end{equation}

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

2.2. Osher方法

待更新

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_{f,\own} c_{f,\own}, & M_f \geqslant 0 \\ \rho_{f,\nei} c_{f,\nei}, & M_f < 0 \end{matrix}\right. \\\\ c_f{\rho_\bfU}_f=\left\{\begin{matrix} {\rho_\bfU}_{f,\own} c_{f,\own}, & M_f \geqslant 0 \\ {\rho_\bfU}_{f,\nei} c_{f,\nei}, & M_f < 0 \end{matrix}\right. \\\\ c_f{\rho_H}_f=\left\{\begin{matrix} {\rho_H}_{f,\own} c_{f,\own}, & M_f \geqslant 0 \\ {\rho_H}_{f,\nei} c_{f,\nei}, & M_f < 0 \end{matrix}\right. \end{split} \end{equation}
方程\eqref{22}与\eqref{23}相结合即常规的通量计算方法。在AUSM中,面心马赫数通过下式进行计算:
\begin{equation}\label{23} M_f=M_{f,\own}^++M_{f,\nei}^- \end{equation}
\begin{equation}\label{24} \begin{split} M^+_{f,\own}=\left\{\begin{matrix} 0.25(M_{f,\own}+1)^2, & |M_{f,\own}|\leqslant 1\\ 0.5(M_{f,\own}+|M_{f,\own}|), & |M_{f,\own}| > 1 \end{matrix}\right. \\\\ M^-_{f,\nei}=\left\{\begin{matrix} -0.25(M_{f,\nei}-1)^2, & |M_{f,\nei}|\leqslant 1\\ 0.5(M_{f,\nei}-|M_{f,\nei}|), & |M_{f,\nei}| > 1 \end{matrix}\right. \end{split} \end{equation}
\begin{equation}\label{25} M_{f,\own}=\frac{\bfU_{f,\own}\cdot\bfn_f}{c_{f,\own}}, M_{f,\nei}=\frac{\bfU_{f,\nei}\cdot\bfn_f}{c_{f,\nei}} \end{equation}
进一步的, AUSM认为对流与压力是两种物理性质不同的机制。因此压力与其他项需要分开处理。网格单元面上的压力处理方法为:
\begin{equation}\label{26} p_f=p_{f,\own}^++p_{f,\nei}^- \end{equation}
\begin{equation}\label{27} \begin{split} p^+_{f,\own}=\left\{\begin{matrix} 0.25p_{f,\own}(M_{f,\own}+1)^2(2-M_{f,\own}), & |M_{f,\own}|\leqslant 1\\ 0.5p_{f,\own}\frac{M_{f,\own}+|M_{f,\own}|}{M_{f,\own}}, & |M_{f,\own}| > 1 \end{matrix}\right. \\\\ p^-_{f,\nei}=\left\{\begin{matrix} 0.25p_{f,\nei}(M_{f,\nei}-1)^2(2+M_{f,\nei}), & |M_{f,\nei}|\leqslant 1\\ 0.5p_{f,\nei}\frac{M_{f,\nei}-|M_{f,\nei}|}{M_{f,\nei}}, & |M_{f,\nei}| > 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中文网