非结构网格与通量分裂


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

1. 非结构网格

由于非结构网格向下兼容的特性,目前主流的CFD求解器均支持非结构网格。但是经典的CFD教材由于出版年限的原因主要介绍的还是结构网格离散方法。从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$上定义的网格体心值即可。

2. 对流方程

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

\begin{equation}\label{m} \frac{\p m}{\p t}+\nabla\cdot\left(m\bfU\right)=0 \end{equation}

其中的$\mathbf{U}$为速度矢量。对方程\eqref{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_f^n\Delta t \end{equation}

其中上标$^n$和$^{n+1}$表示为当前和下一个时间步,下标$_f$表示网格单元面上的插值,$\Delta V_\own$表示$\own$网格单元体积,$\bfS_f$表示$\own$网格单元各个面的面矢量,$\Delta t$表示时间步长,$F$为$\own$网格单元各个面的通量。将方程\eqref{m1},\eqref{m2}带入到方程\eqref{m}有离散的对流方程:

\begin{equation}\label{m4} \left(m_\own^{n+1}-m_\own^n \right)\Delta V_\own =-\sum m_f^nF_f^n\Delta t \end{equation}

\begin{equation}\label{m41} m_\own^{n+1} =m_\own^n -\frac{\Delta t}{\Delta V_\own}\sum F_f^nm_f^n \end{equation}

其中$m_\own^n$,$\Delta V_\own$和$\Delta t$均为已知量。$m_\own^{n+1}$为待求量。$m_f^n$需要插值而来进行封闭。

3. 普通插值

在经过速度压力耦合求解后,方程\eqref{m4}中的$F_f^n$(通量)为已经求得的量。为了封闭方程\eqref{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_f^n > 0 \\ m_\nei, &F_f^n < 0 \end{matrix}\right. \end{equation}

也即

\begin{equation}\label{upwind2} f\left(m_\own,m_\nei\right)=\max\left(F_f^n,0\right)m_\own+\min\left(F_f^n,0\right)m_\nei \end{equation}

将方程\eqref{upwind2}带入到\eqref{m41}有

\begin{equation}\label{upwind3} m_\own^{n+1} =m_\own^n -\frac{\Delta t}{\Delta V_\own}\sum \left(\max\left(F_f^n,0\right)m_\own+\min\left(F_f^n,0\right)m_\nei\right) \end{equation}
4. 通量分裂

通量分裂的思想在一些高精度空间格式上存在。例如对于多相颗粒流,可以使用通量分裂的思想将面上向左流动的颗粒和向右流动的颗粒进行区分。空气动力学中的矢通量分裂思想也来源于颗粒流,即将流体当做颗粒,可以向左也可以向右流动。目前,开源CFD软件OpenFOAM中的rhoCentralFoam中也引入了通量分裂的格式。在非结构网格上,通量分裂需要将流入$\own$和流出$\own$的通量进行区分并分别定义。考虑下面的非结构网格单元,在左图中,其中红色的箭头表示面上流入的通量,黑色的箭头表示面上流出的通量。

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

\begin{equation}\label{f}F_f^nm_f^n=F_\inlet^nm_\inlet^n+F_\outlet^nm_\outlet^n \end{equation}

其中$F_\inlet^nm_\inlet^n$表示进入的通量,$F_\outlet^nm_\outlet^n$表示流出的通量。如果连续性方程可以表示为

\begin{equation}\label{f2} \sum F_f^nm_f^n=0 \end{equation}

则有

\begin{equation}\label{f3} \sum F_\inlet^nm_\inlet^n+\sum F_\outlet^nm_\outlet^n=0 \end{equation}

其表示$\own$网格单元上所有面进入的通量和流出的通量相等。对于$\own$网格单元上的某个面,存在下面的关系:

  • $F_\inlet^nm_\inlet^n$表示流入网格单元的量,其最大值为$0$(否则表示流出)。如果$F_\inlet^nm_\inlet^n<0$,表示有流体流入,则

    \begin{equation}\label{f4}F_\inlet^nm_\inlet^n=\min(F_\rright^n,0)m_\rright^n \end{equation}

  • $F_\outlet^nm_\outlet^n$表示流出网格单元的量,其最小值为$0$(否则表示流入)。如果$F_\outlet^nm_\outlet^n>0$,表示有流体流入,则
    \begin{equation}\label{f5}F_\outlet^nm_\outlet^n=\max(F_\rleft^n,0)m_\rleft^n \end{equation}

    把方程\eqref{f4}和\eqref{f5}带入到方程\eqref{f}有

    \begin{equation}\label{f8} \sum\left(\min(F_\rright^n,0)m_\rright^n+\max(F_\rleft^n,0)m_\rleft^n\right)=0 \end{equation}

    考虑一阶迎风格式有

    \begin{equation}\label{f9} m_\rright^n=m_\nei^n,m_\rleft^n=m_\own^n,F_\rleft^n=F_\own^n,F_\rright^n=F_\nei^n \end{equation}

    将方程\eqref{f9}代入到方程\eqref{f8},其与方程\eqref{upwind3}相同。

    更新历史
    2019.02.02小修

    东岳流体®版权所有
    勘误、讨论、补充内容请前往CFD中文网