返回主页

非结构网格以及通量分裂

由于非结构网格的向下兼容特性,目前主流的CFD求解器均支持非结构网格。但是经典的CFD教材(数值传热学、Computatioanl Methods for Fluid Dynamics等)由于出版年限的原因主要介绍的还是结构网格离散方法。从CFD求解器的角度俩说,非结构网格和结构网格最大的区别在于无规则的网格序号标志。本文力图让用户形成在非结构网格进行离散的思维方式,并简单介绍一下非结构网格下的通量分裂思想。

考虑下面的二维结构网格。结构网格上的网格单元可以用序号进行标识。例如我们可以在$x$方向使用$i$作为标识,$y$方向使用$j$作为标识。对于中心处$(i,j)$,其上下左右的网格单元面可以通过$(i,j+\frac{1}{2})$,$(i,j-\frac{1}{2})$,$(i-\frac{1}{2},j)$,$(i+\frac{1}{2},j)$来标识。在对$i$网格单元进行离散的时候,需要调用面$(i,j+\frac{1}{2})$,面$(i,j-\frac{1}{2})$,面$(i-\frac{1}{2},j)$,面$(i+\frac{1}{2},j)$四点的插值信息。比如如果要对第$256, 29$个网格单元进行离散,那么其周围的面单元标识也可以相应的计算出来,即$256+\frac{1}{2}, 29$等。

考虑下图中的非结构网格。对于非结构网格,这种相对应的标识并不存在。这会带来面变量插值的计算问题。因此在非结构网格中,需要定义面的左右单元对应关系,在下文中省略$y$方向的网格单元序号。如果面$230$左右的两个网格单元标识号为$3$和$99$。这样每个面均被两个网格单元共有。以一阶空间格式举例,在对面$230$上的变量进行插值的时候,只需要寻找网格单元标识号为$3$和$99$上定义的网格体心值即可。

在下面的例子中,考虑一个纯对流方程,对其进行非结构网格上的离散。

对流方程

标量$m$的纯对流方程:

\begin{equation}\label{m} \frac{\p m}{\p t}+\nabla\cdot\left(m\bfU\right)=0 \end{equation}
其中的$\mathbf{U}$为速度矢量。对方程($\ref{m}$)在非结构网格上的$\own$网格点进行欧拉显式时间和空间离散:
\begin{equation}\label{m1} \int \int {\frac{{\partial m}}{{\partial t}}\mathrm{d}V\mathrm{d}t = \left( {m_\own^{n+1} - m_\own^n} \right)\;} {\Delta V_\own} \end{equation}
\begin{equation}\label{m2} \int \int \nabla \cdot \left( m\bfU\right)\mathrm{d}V\mathrm{d}t = \int \int {m\mathbf{U}\mathrm{d}\bfS\mathrm{d}t} = \sum {{m_f^n\bfU_f^n}} \bfS_f\Delta t = \sum m_f^nF^n\Delta t \end{equation}
其中上标$^n$和$^{n+1}$表示为当前和下一个时间步,下标$_f$表示网格单元面上的插值,$\Delta V_\own$表示$\own$网格单元体积,$\bfS_f$表示$\own$网格单元各个面的面矢量,$\Delta t$表示时间步长,$F$为$\own$网格单元各个面的通量。

将方程($\ref{m1}$),($\ref{m2}$)带入到方程($\ref{m}$)有离散的对流方程:

\begin{equation}\label{m4} \left(m_\own^{n+1}-m_\own^n \right)\Delta V_\own =-\sum m_f^nF^n\Delta t \end{equation}
\begin{equation}\label{m41} m_\own^{n+1} =m_\own^n -\frac{\Delta t}{\Delta V_\own}\sum F^nm_f^n \end{equation}
其中的$m_f^n$需要插值而来进行封闭。

普通插值方法

在经过速度压力耦合求解后,方程($\ref{m4}$)中的$F_f^n$为已经求得的量。为了封闭方程($\ref{m4}$),其中的$m_f^n$可以通过插值方法来获取,其可以表示为

\begin{equation}\label{m5} m_f^n=f\left(m_\own,m_\nei\right) \end{equation}
注意,用户在非结构网格下考虑面上的插值一定要转换为$\own$和$\nei$的思想。实际上,在结构网格上方程\eqref{m5}也可以表示为$m_{i+1/2}^n=f\left(m_i,m_{i+1}\right)$,这是因为结构网格上的面标识和网格单元标识是有规律可循的。

回到非结构网格上,对于迎风格式:

\begin{equation}\label{upwind} f\left(m_\own,m_\nei\right)= \left\{\begin{matrix} m_\own &F^n > 0 \\ m_\nei &F^n < 0 \end{matrix}\right. \end{equation}
也即
\begin{equation}\label{upwind2} f\left(m_\own,m_\nei\right)=\max\left(F^n,0\right)m_\own+\min\left(F^n,0\right)m_\nei \end{equation}
将方程($\ref{upwind2}$)带入到方程($\ref{m41}$)有
\begin{equation}\label{upwind3} m_\own^{n+1} =m_\own^n -\frac{\Delta t}{\Delta V_\own}\sum \left(\max\left(F^n,0\right)m_\own+\min\left(F^n,0\right)m_\nei\right) \end{equation}

通量分裂

通量分裂的思想在一些高精度空间格式上存在。例如对于多相颗粒流,可以使用通量分裂(Flux splitting)的思想将面上向左流动的颗粒和向右流动的颗粒进行区分。在空气动力学中,其(Flux-vector splitting)思想也来源于颗粒流,即将流体当做颗粒,可以向左也可以向右流动。目前,开源CFD软件OpenFOAM中的rhoCentralFoam中也引入了通量分裂。

在非结构网格上,通量分裂需要将流入$\own$和流出$\own$的通量进行区分并分别定义。考虑下面的非结构网格单元,在左图中,其中红色的箭头表示面上流入的通量,黑色的箭头表示面上流出的通量。

右图显示的是单独的一个面,其中红色箭头表示变量从面的$\rleft$方向流向$\rright$,黑色箭头表示变量从面的$\rright$方向流向$\rleft$。在通量分裂中,同一个面上的$m_\rleft$和$m_\rright$可以分别定义。在本文中,$\rleft$定义为从$\nei$网格指向$\own$网格,$\rright$定义为从$\own$网格指向$\nei$网格。

我们把面上的通量分裂为$\inlet$通量(入流通量)以及$\outlet$通量(出流通量),其可以定义为:

\begin{equation}\label{f} F^nm_f^n=F_\inlet^n+F_\outlet^n \end{equation}
其中$F_\inlet^n$表示进入的通量,$F_\outlet^n$表示流出的通量。如果连续性方程可以表示为
\begin{equation}\label{f2} \sum F^n=0 \end{equation}
也即
\begin{equation}\label{f3} \sum F_\inlet^n+\sum F_\outlet^n=0 \end{equation}
其表示$\own$网格单元上所有面进入的通量$=$流出的通量。

对于某个面,存在下面的关系:

  • 如果$F>0$,表示流体从 $\own$($\rright$) 流出到 $\nei$($\rleft$);
  • 如果$F<0$,表示流体从 $\nei$($\rleft$) 流入到 $\own$($\rright$);
  • 如果$F>0$,有:

    \begin{equation}\label{f4} F_\inlet=0,F_\outlet=Fm_\rright \end{equation}
    如果$F<0$,有:
    \begin{equation}\label{f5} F_\inlet=Fm_\rleft,F_\outlet=0 \end{equation}
    结合方程($\ref{f4}$)和($\ref{f5}$)有:
    \begin{equation}\label{f6} F_\inlet=\min\left(F^n,0\right)m_\rright \end{equation}
    \begin{equation}\label{f7} F_\outlet=\max\left(F^n,0\right)m_\rleft \end{equation}
    把方程($\ref{f6}$)和($\ref{f7}$)带入到方程($\ref{f}$)有
    \begin{equation}\label{f8} F^n=\max\left(F^n,0\right)m_\rright+\min\left(F^n,0\right)m_\rleft \end{equation}
    考虑一阶迎风格式有
    \begin{equation}\label{f9} m_\rright=m_\own,m_\rleft=m_\nei \end{equation}
    将方程($\ref{f9}$)代入到方程($\ref{f8}$),其与方程($\ref{upwind3}$)相同。

    非结构网格上的面积分

    在计算网格单元通量的时候存在俩种方法:

  • 遍历$N_1$个网格单元标志,并在每个网格单元上对其拥有的$N_2$个面进行加和;
  • 遍历$N_2$个网格单元的面并进行数值操作;

  • 很明显第一种操作需要进行$N_1\times N_2$次,因此在开源CFD软件OpenFOAM中采用的是第二种方式。

    考虑下图中的网格系统,其中黑色数字表示网格标识,其由网格生成软件产生。考虑23号网格单元,其周围的16, 78, 30, 14, 25均表示其相邻的网格单元标号(nei)。红色数字表示网格单元23的面标识。其由网格生成软件产生。

    针对每个面标识,OpenFOAM进一步的定义了每个面的own以及nei。own为面相邻的具有较低标识的网格单元。nei为面相邻的具有较低标识的网格单元(如下图所示)。

    上图中蓝色的箭头表示已知的速度矢量,其具有方向。OpenFOAM进一步定义每个面上的通量如果大于零,则从own指向nei,反之亦然(如下图所示)。

    假定各个面上的通量如下图所示。从图中可以看出流出标量$T$的通量为$9+1=10$,流入的通量为$2+7+1=10$。因此净通量为$0$。下面我们用OpenFOAM遍历面的方式来计算23网格单元的变化。

    首先初始网格23单元的$T=0$,然后遍历面:

  • 面0:$T_{23}=T_{\nei}=0-F_{0}=0-(-1)=1$;
  • 面1:$T_{23}=T_{\own}=1+F_{1}=1+(9)=10$;
  • 面2:$T_{23}=T_{\own}=10+F_{2}=10+(-2)=8$;
  • 面3:$T_{23}=T_{\nei}=8-F_{3}=8-(7)=1$;
  • 面4:$T_{23}=T_{\own}=1+F_{4}=1+(-1)=0$;

  • 最终,23网格单元的净通量为$0$。

    更新历史

    2017.12.05:创立页面

    东岳流体dyfluid.com
    勘误、讨论、补充内容请前往CFD中国