• 返回主页
  • boundaryFoam解析

    Reviewed by:
    求解器类型:稳态不可压缩1D湍流求解器。

    boundaryFoam为OpenFOAM中一个比较特殊的求解器。其为一个多网格单元的1D求解器。在应用中,我们可以用此求解器来求解充分发展的管道流(类似泊肃叶流)来获得最快速的收敛解。例如:模拟一个充分长的2D管道,要获得充分发展的湍流,则需要此管道足够的长,这样会大大增加网格数量并大大增加计算时间。但通过boundaryFoam,仅仅需要几十个网格就可以完成求解。

    此求解器的特点如下:

    1.使用OpenFOAM处理多网格单元的1D问题;
    2.1D问题下,速度矢量仅具有一个分量,$\mathbf{U}_x$,其他方向人工置零。且$\mathbf{U}_x$在y方向非均匀分布;
    3.1D问题下,y方向的压力梯度为0,即:$\mathrm{d}p/\mathrm{d}y=0$;
    4.无对流项,因此植入的速度方程为:$\frac{\mathrm{d}p}{\mathrm{d}x}=\frac{\mathrm{d}}{\mathrm{d}y}(ν \frac{\mathrm{d}u}{\mathrm{d}y})+$湍流项;
    5.压力x分量和速度的耦合采用迭代方法求解,速度型线满足连续性方程;

    首先进入createFields.H文件(略去部分内容):

        dimensionedVector Ubar(transportProperties.lookup("Ubar"));//在transportProperties中定义一个矢量Ubar,这个Ubar是管道中的平均速度(实际上为抛物线型
    )。注意,此处我们的求解器定义为1D求解器,因此Ubar的形式应该为:(5,0,0)这样的形式,即仅有一个分量。
    vector flowDirection = (Ubar/mag(Ubar)).value(); tensor flowMask = sqr(flowDirection); dimensionedVector gradP ( "gradP", dimensionSet(0, 1, -2, 0, 0), vector::zero );

    其中的flowDirection为定义一个流动的矢量,不管Ubar为任何值,flowDirection只有3种情况: \begin{equation} (1,0,0),(0,1,0),(0,0,1) \end{equation} 然后定义flowMask张量,具有9个元素。同样,在我们的1D情况下,也只有3种情况:

    $ \left( \begin{matrix} 1 & 0 & 0\\ 0 & 0 & 0\\ 0 & 0 & 0 \end{matrix} \right) , \left( \begin{matrix} 0 & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 & 0 \end{matrix} \right) , \left( \begin{matrix} 0 & 0 & 0\\ 0 & 0 & 0\\ 0 & 0 & 1 \end{matrix} \right) $


    下面进入主程序boundaryFoam.C(略去部分内容):

    while (runTime.loop())
        {
            Info<< "Time = " << runTime.timeName() << nl << endl;
    
            fvVectorMatrix divR(turbulence->divDevReff(U));//定义湍流项
            divR.source() = flowMask & divR.source();//在有了离散divR的矩阵系数后,因为OpenFOAM在求解1D问题的时候,
    需要人工把第二个维度的分量去掉(这并不同于OpenFOAM处理2D和3D问题的方式)
    fvVectorMatrix UEqn ( divR == gradP& divR.source();//组建方程 );

    需要提及的是,在boundaryFoam中省略了对流项:fvm::div(phi, U)。因为在1D的情况下,$\mathbf{U}_y=0$且$\frac{\partial \mathbf{U}_x}{\partial x}=0$。

            UEqn.relax();
    
            UEqn.solve();
    
            // Correct driving force for a constant volume flow rate
            dimensionedVector UbarStar = flowMask & U.weightedAverage(mesh.V());//步骤1
    
            U += (Ubar - UbarStar);//步骤2
            gradP += (Ubar - UbarStar)/(1.0/UEqn.A())().weightedAverage(mesh.V());//步骤3
    
            turbulence->correct();//步骤4
    
            Info<< "Uncorrected Ubar = " << (flowDirection & UbarStar.value())
                << ", pressure gradient = " << (flowDirection & gradP.value())
                << endl;
    
    在求解速度的时候 ,需要迭代求解,其可以表示为:

    1.对速度方程使用已知的gradP进行求解,我们有$\mathbf{U}'$;
    2.$\mathbf{U}'$和平均速度不同,计算$\mathbf{U}'$和Ubar的差值,对$\mathbf{U}'$进行修正:
    3.修正后我们有$\mathbf{U}''$,使用$\mathbf{U}''$对压力进行修正:
    4.湍流动能,湍流耗散方程计算;
    5.我们有符合$\mathbf{U}''$的压力梯度,使用此压力梯度,重复第一步循环计算;

    最终收敛的时候,我们的平均修正速度会和定义的Ubar相同,压力梯度和速度也符合NS方程。
    如何引用:
    李东岳. boundaryFoam解析[OL]. http://dyfluid.com/boundaryFoam.html, 2016.04.08

    东岳流体
    info@dyfluid.com