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

Bob在学习OpenFOAM

目的：求解充分发展槽道流动，进出口设置为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);
      }
   }
}

东岳

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

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

Bob在学习OpenFOAM

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

CjjJoy

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

东岳

@CjjJoy

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

D.Benjamin

@东岳老师您好，我想请教几个问题：

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

东岳

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

不知道，在声明pGrad之后，需要赋值，比如：

pGrad = fvc::grad(p);

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

U ((10 0 0) 0);

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

D.Benjamin

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

东岳

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

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

气固两相流

@d-benjamin 请问按照方案二加入压力梯度后，在算例中应该如何给压力梯度赋值呢？

东岳

他已经赋值了

pGrad = fvc::grad(p);

气固两相流

@东岳

东岳老师，我的错误显示是我没有在0文件夹下建立 pGrad 文件，建立pGrad的文件的话，这个地方应该怎么赋值，而且如果已经在blockMesh中设置周期性边界条件cyclic的话，那么pGrad不就和p一样了依然没有压差么？

东岳

volVectorField pGrad
(
    IOobject
    (
        "pGrad",
        runTime.timeName(),
        mesh,
        //IOobject::MUST_READ,
        IOobject::NO_READ,
        IOobject::AUTO_WRITE
    ),
    mesh,
    Zero
    //vector(1, 0, 0)
);

因为你代码pGrad为必须读取，你做上面的改动。可以把pGrad设置为默认的(0 0 0)。也可以写为vector(1, 0, 0)

D.Benjamin

@气固两相流你好，正如东岳老师所说，只要在0文件夹下有p文件，就会按照给的式子进行赋值了

气固两相流

@东岳非常感谢东岳老师的回答，但是直接对他赋具体值zero或者vector（1 0 0）没有成功，但是fvc::grad( p )是没有问题的

气固两相流

@d-benjamin 嗯嗯，谢谢您的回答，直接赋值fvc::grad ( p)是没有问题的

我是河滩

你好，你在尝试1中加压力梯度源项，为什么对UEqn.H和pEqn.H都做了修改？
simple算法求NS方程是分两步，UEqn.H是预估步，pEqn.H包括校正步。既然在预估步就加入了源项，为什么在校正步pEqn.H还要添加源项呢？
另外，对于尝试3 fvoptions，这种方法的源项是加在simple求解的哪一步中了？是和尝试2相同吗？

D.Benjamin

@我是河滩你好，他只修改了U的方程，并没有修改p方程，只是定义了一个新的压力梯度场。他上面三个尝试应该都是一样的。有一点我不是很明白，就是尝试3中的volumeMode，两个关键字absolute和specific有什么具体的区别。如果东岳老师可以帮忙解答一下就好了。@东岳

我是河滩

@东岳在动量方程如何加入加压力梯度源项？中说：

多相流的的表面张力项

如何把源项加入到压力方程中，和动量方程加源项有何区别？

张RP

@D-Benjamin 请参照 http://www.cfd-china.com/topic/208/fvoptions-rho-yi-与scalarsemiimplicitsource/9

东岳

@我是河滩
加入到压力方程更加无震荡
加入到动量方程更简单

我是河滩

谢谢。

东岳

@我是河滩你做什么方向的啊，

我是河滩

运动物体计算，浸入边界。