动量方程如何加入加压力梯度源项?



  • 目的:求解充分发展槽道流动,进出口设置为cyclic周期性边界。通过动量方程中加入流动的驱动力:压力梯度。
    求解器:OpenFOAM v3.0+
    其他参考:fvOptions 的demo:simpleFoam\channel395,设置的是平均速度Ubar,还不是驱动压力;
    目前遇到的问题:通过下面的这3种方法,加入压力梯度源项后,测试了 kEpsilon 和 kOmega 湍流模型,kEpsilon表现正常,kOmega在槽道中心区域算出诡异的速度分布,分析问题可能出现在两方面:1)加源项的方法,2)kOmega湍流模型本身。具体算例,我会在另一个帖子里说明。

    这篇帖子向大家请教添加源项的方法,下面是我尝试的几种方法,不知道是不是正确,请大家帮忙批评指正 :)

    • 尝试1:
      createField 中定义一个压力梯度向量 pGrad,变量类型为volVectorField
    volVectorField pGrad
    (
        IOobject
        (
            "pGrad",
            runTime.timeName(),
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh
    );
    

    simpleFoam 作如下改动:

    solve(UEqn() == -fvc::grad(p)-pGrad); // UEqn.H 
    
    U = HbyA - rAtU()*fvc::grad(p) - rAtU()*pGrad; // pEqn.H
    
    • 尝试2:
      同样需要在 createField.H 中定义 pGrad,不同的是换一种方法改动
      simpleFoam,在 UEqn 中加入 -pGrad:
        tmp<fvVectorMatrix> UEqn
       (
           fvm::div(phi, U)
         + MRF.DDt(U)
         + turbulence->divDevReff(U)
        ==
           fvOptions(U) - pGrad
       );
    
    • 尝试3:
      通过fvOptions加入源项,该方法参考 cfdonline
      这样给动量方程加入了 pGrad=10 的压力梯度
    momentumSource
    {
       type vectorSemiImplicitSource;
       active on;
       selectionMode all;
    
       vectorSemiImplicitSourceCoeffs
       {
          //volumeMode      absolute; // specific
          volumeMode        specific;
          injectionRateSuSp
          {
             U           ( (10 0 0) 0);
          }
       }
    }
    

  • 网格教授 OpenFOAM教授 管理员

    Hi,我觉得你这3种方法都差不多:expressionless: ,

    另外一种不同的方法是把源项添加到压力方程中,OpenFOAM经常把fvOptions不可包括的源项添加到压力方程中。:expressionless: 例如多相流的的表面张力项



  • @李东岳 谢谢你关注我的问题,我又尝试了槽道层流,上面的源项方法没有问题,个人比较偏向第2种,把压力梯度加入到 UEqn 矩阵中,感觉有点像隐式的求解方法。



  • @李东岳 你好,如何在压力方程中添加源项 包含多相流的的表面张力项?有没有推荐的文档可以学习下?谢谢


  • 网格教授 OpenFOAM教授 管理员

    @CjjJoy

    这个interFoam里面加入到了UEqn,源项添加代码上非常简单啊没有任何难点?你写好了源项放到方程右边不就可以了?



  • @东岳 老师您好,我想请教几个问题:

    1. 楼主加入了一个源项体积矢量pGrad,这个矢量,OpenFOAM是如何识别为压力梯度呢?也就是说OpenFOAM怎么知道pGrad就是▽p呢?楼主上面的场的设置,OpenFOAM应该只认识pGrad是一矢量,而不知道它就是代表着▽p吧?
    2. 楼主尝试的第三种方法,U ((10 0 0) 0);表示什么意思?(10 0 0)不就表示速度矢量值了吗?
      若能够获得老师的回复,我将不胜感激!!

  • 网格教授 OpenFOAM教授 管理员

    也就是说OpenFOAM怎么知道pGrad就是▽p呢?

    不知道,在声明pGrad之后,需要赋值,比如:

    pGrad = fvc::grad(p);
    

    这样才可以。楼主那个应该没写这部分

    U ((10 0 0) 0);

    这不是简单的速度矢量了,这是给injectionRateSuSp关键词赋值



  • @东岳 老师您好,谢谢你的回复,我还想追问一下,赋值的这部分代码应该写到哪个地方?下面我写了两个方案,哪一个方案是正确的?还有,我发现在laplacianFoam的解析中,有两种定义场的方式,如下面的图片所示,定义场"T"的时候是mesh,而定义场"gradTx"的时候,mesh被替换成了gradT.component(vector::x),正是因为对此疑惑,我才有上面两种方案是否合理的疑惑,望老师给与答复,感激不尽!
    0_1537532008564_1.PNG

    0_1537532016263_2.PNG

    0_1537532538546_3.PNG

    0_1537532648432_4.PNG


  • 网格教授 OpenFOAM教授 管理员

    方案二是正确的,方案一mesh下面不能赋值,可以直接写fvc::grad(p)

    同理,下面那个T,如果直接给gradT....,就不需要写mesh了,直接赋值。