• 返回主页
  • OpenFOAM中的div与snGrad操作符

    作者:邱小平
    编辑:赵一铭

    OpenFOAM的方便之处之一是利用C++的类模板和函数重载等技术定义了很多各种离散操作符,如divlaplaciangrad等等。利用这些操作符,很容易就能对偏微分方程进行离散,并构建起线性方程组。但是,这些操作符真正执行的运算,却需要结合有限体积方法的本质来理解一番才能真正掌握。下面对OpenFOAM中的divsnGrad操作符进行一点解读。

    div操作符的本质

    div操作符表面看,是计算散度的,实际上,在OpenFAOM中,div操作符的作用是加和,比如说$ \nabla \cdot (\bfU \bfU)$,在OpenFOAM中表示为fvm::div(phi,U),这段代码真正执行的是$\sum_f \bfU_f \phi_f$运算,即将每个网格包含的面上的流率与速度乘积,然后再加起来。再比如twoPhaseEulerFoam的UaEqn方程有一项是fvm::sp(fvc::div(phia),Ua),其对应的公式是$\bfU_a(\nabla \cdot \bfU_a)$。为什么是这样呢?以下试图对背后的原理做一点解释。

    单相流的动量守恒方程的微分形式如下:\begin{equation}\frac{\partial \bfU}{\partial t}+\nabla\cdot (\bfU \bfU)+\nabla \cdot \dev(-\nu_{eff} (\nabla \bfU+(\nabla \bfU)^\rT))=-\nabla p +Q\end{equation}为了使用有限体积方法,需要将动量方程写成积分形式,即\begin{equation}\int_V\left [\frac{\partial \bfU}{\partial t}+\nabla\cdot (\bfU\bfU)+\nabla \cdot \dev(-\nu_{eff} (\nabla \bfU+(\nabla \bfU)^T))\right ]\rd V=\int_V\left [ -\nabla p +Q \right ]\rd V\end{equation}这里只分析$\int_V \nabla\cdot (\bfU\bfU) \rd V$这一项,利用高斯定理,可以将这一个体积分,转化成对包围该体积微元的表面的面积分:$\oint_{\partial V} (\bfU \bfU)\cdot \rd \bfS$,其中 $\rd \bfS$ 是面积微元矢量。又因为实际操作中,一个体积微元总是由有限的几个面组成的,所以有:\begin{equation}\oint_{\partial V} (\bfU\bfU)\cdot \rd\bfS=\sum_f\int_{\bfS_f}(\bfU\bfU)\cdot \rd\bfS_f=\sum_f(\bfU\bfU)_f\cdot \bfS_f\end{equation}其中\begin{equation}(\bfU\bfU)_f=\frac{1}{\mag(\bfS_f)}\int_{\bfS_f}(\bfU\bfU)\cdot \rd\bfS_f\end{equation}这里做一个近似:\begin{equation}(\bfU\bfU)_f\approx(\bfU_f\bfU_f)\end{equation}于是得到\begin{equation}\sum_f(\bfU\bfU)_f\cdot \bfS_f\approx \sum_f(\bfU_f \bfU_f)\cdot \bfS_f\end{equation}运用张量计算规则$[\mathbf{uv}\cdot\mathbf{w}]=\mathbf{u}(\mathbf{v} \cdot \mathbf{w})$,得到\begin{equation}\sum_f(\bfU_f\bfU_f)\cdot \bfS_f=\sum_f\bfU_f(\bfU_f\cdot \bfS_f)=\sum_f\bfU_f\phi_f\end{equation}综上,OpenFOAM中的代码fvm::div(UU)对应的公式其实是$\int_V \nabla\cdot (\bfU\bfU) \rd V$,而根据推导,在有限体积方法中 $\int_V \nabla\cdot (\bfU\bfU) \rd V=\sum_f(\bfU_f\phi_f)$,所以,fvm::div(UU)实质上进行的运算是,把网格当作体积微元,将网格中心的速度$\bfU$插值到包围该网格的所有面上的面心上得到 $\bfU_f$,并计算每个面上的速度通量$\phi_f$,然后返回每一个面上的速度与通量乘积的加和。$\bfU_f$的计算需要用到本网格与邻近网格的速度值,离散格式的作用就体现在如果利用本网格与邻近网格的速度得到面上的速度。

    再来看上面提到的另一项$\bfU_a(\nabla \cdot \bfU_a)$,根据上面类似的推导\begin{equation}\int_V\left [ \bfU_a(\nabla \cdot \bfU_a)\right ]\rd V=\bfU_a\int_V(\nabla \cdot \bfU_a)\rd V\end{equation}注意,这里之所以能这样变换,是因为在一个体积微元$\rd V$内,$\bfU_a$是常数。根据高斯定理\begin{equation}\int_V(\nabla \cdot \bfU_a)\rd V=\oint_{\partial V}\bfU_a\cdot \rd \bfS=\sum_f(\bfU_a)_f\cdot \bfS_f=\sum_f(\phi_a)_f\end{equation}于是得到\begin{equation}\int_V\left [ \bfU_a(\nabla \cdot \bfU_a)\right ]\rd V=\bfU_a\sum_f(\phi_a)_f\end{equation}这一项在OpenFOAM中的表达是fvm::sp(fvc::div(phia),Ua),含义就很明显了,这一项相当于是一个系数与需要求解的量$\bfU_a$的乘积,所以被当作隐式的源项来处里。注意我先前以为把$(\nabla \cdot \bfU_a)$当作显式处理是人为简化的结果,其实不然,这是自然而然的结果。而在这里也可以看出,$\sum_f(\phi_a)_f$对应的代码是fvc::div(phia),也印证了上面的观点,即div操作符本质上是在作加和运算。

    gradsnGrad

    在twoPhaseEulerFoam中,$\frac{\nabla \alpha}{\alpha}$在两个不同的地方用了两种不同的表示。\begin{equation}\label{a}\nabla \cdot \left \{ [\nu_{eff,a} \frac{\nabla \alpha}{\alpha}][\bfU_a] \right \}\end{equation}对应的代码是:fvm::div(phiRa,Ua),其中

    1
    2
    phiRa=-fvc::interpolate(nuEffa)*mesh.magSf()*fvc::snGrad(alpha)
         /fvc::interpolate(alpha + scalar(0.001));
    

    而另一项\begin{equation}\label{b}[\frac{\nabla \alpha }{\alpha}]\cdot{Rca}\end{equation}对应的代码是

    1
    fvc::grad(alpha)/fvc::average(alpha + scalar(0.001)) & Rca
    

    下面是我对这个的理解:

    对于方程$\eqref{a}$其处理方法跟$\nabla \cdot(\bfU_a\bfU_a)$是一样的,因为$ \frac{\nabla \alpha}{\alpha}$也是一个矢量。phiRa相当于是$\left [\nu_{eff,a} \frac{\nabla \alpha}{\alpha} \right ]$这个矢量的面通量,类比于phiaphia的定义是linearInterpolate(Ua) & mesh.Sf(),而phiRa理论上应该也可以定义成类似于linearInterpolate(gradalpha) & mesh.Sf(),前提是要先定义一个 volVectorField gradalpha=-nuEffa*fvc::grad(alpha)/alpha。但是考虑到界面通量本质上是先将一个volVectorField插值到面上,并乘以面积矢量,对于 $\frac{\nabla \alpha}{\alpha}$,完全可以直接求出每个面上的$\frac{\nabla \alpha}{\alpha}$,然后乘以面积矢量,而不需要先建立体中心的$\frac{\nabla \alpha}{\alpha}$再插值到面上。snGrad就是这样一个用来求面上的梯度量的函数,它求解面上的梯度时,采用如下公式:\begin{equation}(\nabla \phi)_f=\frac{\phi_N-\phi_P}{|\mathbf{d}|}\end{equation}即用相邻网格的值减去面所属网格的值,再除以两个网格中心矢量的模。注意这个处理对于正交网格是精确的,对于非正交网格,只是一个近似处理。而且要注意,这样求得的面上的梯度值已经是一个标量了。再回到phiRa的定义,既然fvc::snGrad(alpha)已经是标量了,那只要再乘以面积矢量的模mesh.magSf(),便是界面上的通量了。fvc::interpolate(nuEffa)fvc::interpolate(alpha + scalar(0.001))则分别是将volField插值到面。

    再来看方程$\eqref{b}$,这一项是被当作显式的源项来处理,所以,我们需要得到的是一个volVcetorField。所以这里将$\nabla \alpha$处理成volVectorField,再与作为volTensorField的Rca进行点乘,得到的就是 volVcetorField。因此,这里的$\frac{\nabla \alpha }{\alpha}$对应的代码是fvc::grad(alpha)/fvc::average(alpha + scalar(0.001))。注意这里的fvc::average()函数的定义

     1
     2
     3
     4
     5
     6
     7
     8
     9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    template<class Type>
    tmp<GeometricField<Type, fvPatchField, volMesh> >
    average
    (
        const GeometricField<Type, fvsPatchField, surfaceMesh>& ssf
    )
    {
        const fvMesh& mesh = ssf.mesh();
    
        tmp<GeometricField<Type, fvPatchField, volMesh> > taverage
        (
            new GeometricField<Type, fvPatchField, volMesh>
            (
                IOobject
                (
                    "average("+ssf.name()+')',
                    ssf.instance(),
                    mesh,
                    IOobject::NO_READ,
                    IOobject::NO_WRITE
                ),
                mesh,
                ssf.dimensions()
            )
        );
    
        GeometricField<Type, fvPatchField, volMesh>& av = taverage();
    
        av.internalField() =
        (
            surfaceSum(mesh.magSf()*ssf)/surfaceSum(mesh.magSf())
        )().internalField();
    

    可以看出,fvc::average()函数的返回类型是tmp<GeometricField<Type, fvPatchField, volMesh> >,对应alpha,返回类型就是tmp<GeometricField<scalar, fvPatchField, volMesh> >,即tmp<volScalarField>average函数的核心定义见代码29-32行,可见average采用的是面积加权平均,即\begin{equation}\overline{\phi}=\frac{\sum_f\phi_f*|\bfS_f|}{\sum_f |\bfS_f|}\end{equation}如果给average的参数类型是surfaceField,那么直接计算平均值后返回volField,如果给的参数的volField,那就先将volField插值到面,计算平均值后再返回volField


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