加载中…
个人资料
WayneWangY
WayneWangY
  • 博客等级:
  • 博客积分:0
  • 博客访问:91,777
  • 关注人气:155
  • 获赠金笔:0支
  • 赠出金笔:0支
  • 荣誉徽章:
相关博文
推荐博文
谁看过这篇博文
加载中…
正文 字体大小:

OpenFOAM求解器(二):通量与求解器算法

(2015-10-30 17:14:58)
标签:

openfoam

openfoam求解器

openfoam求解器开发

openfoam通量

peqn.flux

分类: OpenFOAM-实用技术
更新2016.05.09
1.
参照http://www.tfd.chalmers.se/~hani/kurser/OS_CFD_2007/rhiechow.pdf
公式(9)是半离散化的动量方程。
公式(13)是压力泊松方程:fvm::laplacian(rAU, p) == fvc::div(phiHbyA)
使用的是phiHbyA来求解压力场,故依据(11)和(9),自然有
 phi = phiHbyA - pEqn.flux();
 U = HbyA - rAU*fvc::grad(p);

fvc::div(phiHbyA)的问题是由于还要做一个积分,可参考
https://drive.google.com/file/d/0B2lpdhG-Zh05Y0ZzOHVLb3lGekk/edit中的(2.64)和

template
tmp >
div
(
    const GeometricField& ssf
)
{
    return tmp >
    (
        new GeometricField
        (
            "div("+ssf.name()+')',
            fvc::surfaceIntegrate(ssf)
        )
    );
}


2.
    For cases which do no have a pressure boundary adjust the balance of
    fluxes to obey continuity.  Return true if the domain is closed.

bool adjustPhi
(
    surfaceScalarField& phi,
    const volVectorField& U,
    volScalarField& p
);
如果算例存在给定压力的边界条件,直接返回false,flag是 p.needReference();
对于给定固定值且不是inletOutletFvPatchVectorField类型的速度边界,计算总的入流量massIn,总的固定出流量fixedMassOut;
对于其他速度边界,继续计算massIn和总的可调节出流量adjustableMassOut;

        if
        (
            magAdjustableMassOut > VSMALL
         && magAdjustableMassOut/totalFlux > SMALL
        )
        {
            massCorr = (massIn - fixedMassOut)/adjustableMassOut;
        }
        else if (mag(fixedMassOut - massIn)/totalFlux > 1e-8)
        {
            FatalErrorIn (这是这个函数唯一有报错的地方)
  
   通俗解释上面的代码:啥时候报错呢?就是可调节流量很小,而流量的不平衡较大时。
如果没有报错的话,就按照比例系数massCorr对所有的可调节流量进行修正。
参考 http://www.cfd-china.com/topic/305/adjustphi的作用是检查边界条件

http://www.cfd-online.com/Forums/openfoam-programming-development/83876-pref-adjustphi.html

3.correctBoundaryConditions();
template<<span style=" color:#808000;">class Type, template<<span style=" color:#808000;">class> class PatchField, class GeoMesh>
void Foam::GeometricField<<span style=" color:#800080;">Type, PatchField, GeoMesh>::
correctBoundaryConditions()
{
    this->setUpToDate();
    storeOldTimes();
    boundaryField_.evaluate();
}
  这里涉及evaluate()函数在每个边界条件中的重载和边界条件的调用,参考http://www.cfd-online.com/Forums/openfoam-programming-development/129271-how-boundary-conditions-called-openfoam-solvers.html
最终的实现在每个边界条件中找。

4.ddtCorr(U, phi) 
位于fvcDdt.H
template<<span style=" color:#808000;">class Type>
tmptypename flux<<span style=" color:#800080;">Type>::type, fvsPatchField, surfaceMesh> >
ddtCorr
(
    const GeometricField<<span style=" color:#800080;">Type, fvPatchField, volMesh>& U,
    const GeometricField
    <</pre>
        typename flux<<span style=" color:#800080;">Type>::type,
        fvsPatchField,
        surfaceMesh
    >& phi
)
{
    return fv::ddtScheme<<span style=" color:#800080;">Type>::New
    (
        U.mesh(),
        U.mesh().ddtScheme("ddt(" + U.name() + ')')
    )().fvcDdtPhiCorr(U, phi); 这个函数位于ddtScheme.H
}

  virtual tmp<<span style="color: rgb(128, 0, 128);">fluxFieldType> fvcDdtPhiCorr
        (
            const GeometricField<<span style=" color:#800080;">Type, fvPatchField, volMesh>& U,
            const fluxFieldType& phi
        ) = 0;

纯虚函数。其中一个实现如下:

template<<span style=" color:#808000;">class Type>
tmp<<span style=" color:#808000;">typename EulerDdtScheme<<span style=" color:#800080;">Type>::fluxFieldType>
EulerDdtScheme<<span style=" color:#800080;">Type>::fvcDdtPhiCorr
(
    const GeometricField<<span style=" color:#800080;">Type, fvPatchField, volMesh>& U,
    const fluxFieldType& phi
)
{
    dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT();

    fluxFieldType phiCorr
    (
        phi.oldTime() - (mesh().Sf() & fvc::interpolate(U.oldTime()))
    );

    return tmp<<span style=" color:#800080;">fluxFieldType>
    (
        new fluxFieldType
        (
            IOobject
            (
                "ddtCorr(" + U.name() + ',' + phi.name() + ')',
                mesh().time().timeName(),
                mesh()
            ),
            this->fvcDdtPhiCoeff(U.oldTime(), phi.oldTime(), phiCorr)
           *rDeltaT*phiCorr
        )
    );
}

5.pEqn.flux()的代码实现

一个代码解读见
http://www.cfd-online.com/Forums/openfoam/114933-description-flux-method.html



以下内容为早期的笔记,不保证正确,仅供参考。

===========================================================


通量phi在of中十分关键,理解并正确使用它可以使得我们在开发自己的求解器时避免很多问题。

通量的一般意义见[2]。本文以分析OF中的通量为突破口来解释OpenFOAM中求解器算法的代码实现。

在porousMultiphaseFoam中[1],求解完压力方程后,通过
phiP = pEqn.flux();
获得压力梯度造成的流量phiP。这里的pEqn.flux()是啥意思呢?

在solvers目录下通过 find ./ -name "*" | xargs grep "pEqn.flux"命令查看一下 pEqn.flux()的使用情况,有如下几种:

phi = phiHbyA - pEqn.flux(); //如icoFoam, simpleFoam, pimpleFoam

phi == pEqn.flux(); //如reactingFoam, rhoReactingFoam
phi = phiHbyA + pEqn.flux();

phic = phiHbyA - pEqn.flux()/alphacf; //如DPMFoam
+rAUc*fvc::reconstruct((phicForces - pEqn.flux()/alphacf)/rAUcf);

phi += (phiGradp + pEqn.flux())/rhof; //如cavitatingFoam

phi = pEqn.flux(); //如sonicFoam

phiP = pEqn.flux(); //如porousMultiphaseFoam

参考Jasak的博士论文或者[3-4],对于不可压系统:
OpenFOAM求解器(二):通量与求解器算法(*)
需要解决对流项的非线性和速度、压力的耦合问题。
通用形式的输运方程:
OpenFOAM求解器(二):通量与求解器算法(3.2)
对流项的离散:
OpenFOAM求解器(二):通量与求解器算法(3.17)
一样地,(*)中对流项的离散是:
OpenFOAM求解器(二):通量与求解器算法
线性化对流项,通量的速度用当前时间步来计算,保留另一个速度为未知量[7]。注意这里的通量F 要满足连续性方程,然后用它来计算系数a_p和a_N。关于线性化的讨论可参见东岳流体技术文档:PISO,SIMPLE 详解

结合东岳流体技术文档:icoFoam 详解,我们来看动量方程的半离散化:
OpenFOAM求解器(二):通量与求解器算法(3.136)
速度可以表示为:
OpenFOAM求解器(二):通量与求解器算法(3.139)
面上的速度由插值得到:
OpenFOAM求解器(二):通量与求解器算法(3.140)
面上的通量:
OpenFOAM求解器(二):通量与求解器算法(3.144)
其中,H(U)包括对流扩散项(输运),瞬态项和除压力梯度之外的所有源项:
OpenFOAM求解器(二):通量与求解器算法(3.137)
离散化的连续性方程:
OpenFOAM求解器(二):通量与求解器算法(3.138)
将(3.140)代入(3.138),得压力方程:
OpenFOAM求解器(二):通量与求解器算法(3.141)
最终得到离散化的不可压缩NS方程:
OpenFOAM求解器(二):通量与求解器算法(3.142;3.143)
其中3.142即[7]中的(2.4)。

对于瞬态PISO 算法,首先是动量预测,求解[7]中(2.4)得到预测速度,即icoFoam中的
  fvVectorMatrix UEqn
        (
            fvm::ddt(U)
          + fvm::div(phi, U)
          - fvm::laplacian(nu, U)  //层流
        );
        solve(UEqn == -fvc::grad(p));
simpleFoam中的动量预测UEqn.H与上面的相同,使用turbulence->divDevReff(U)考虑了湍流,使用fvOptions考虑了源项,并引入了方程的松弛。(这块比较熟了,有任何问题可留言讨论)

// --- PISO loop
       然后是压力求解,
            volScalarField rAU(1.0/UEqn.A()); //A() 返回fvVectorMatrix的中心系数,这里即a_P。This is the famous "operator splitting".

            volVectorField HbyA("HbyA", U);
            HbyA = rAU*UEqn.H(); //H() 返回fvVectorMatrix的H源项,这里即H(U)。HbyA见[7]中(3.1)
            surfaceScalarField phiHbyA
            (
                "phiHbyA",
                (fvc::interpolate(HbyA) & mesh.Sf())
              + fvc::interpolate(rAU)*fvc::ddtCorr(U, phi) //simpleFoam中这里是 surfaceScalarField phiHbyA("phiHbyA", fvc::interpolate(HbyA) & mesh.Sf()); ddtCorr最后分析。
            );
            


            adjustPhi(phiHbyA, U, p); //修改压力为zeroGradient边界条件的通量以确保总体的通量守恒,见[8]

            for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
            {
                fvScalarMatrix pEqn
                (
                    fvm::laplacian(rAU, p) == fvc::div(phiHbyA) //用不满足无散度的通量求解压力方程3.143,[7]中的(3.9)
                );

                pEqn.setReference(pRefCell, pRefValue); //压力参考值
                pEqn.solve();

                if (nonOrth == nNonOrthCorr)
                {
                    phi = phiHbyA - pEqn.flux();  //得到无散度的通量[6],最后分析
                }
            }

最后是显式速度修正(实际上HbyA相对于U滞后,故需要进行迭代,参见Jasak的博士论文第148页的精彩论述):
            U = HbyA - rAU*fvc::grad(p); //3.139,此速度满足连续性方程
            U.correctBoundaryConditions();

pisoFoam.c与上述过程基本一致。

至此,本文还剩下两处没有解释:
  "phiHbyA",(fvc::interpolate(HbyA) & mesh.Sf()) + fvc::interpolate(rAU)*fvc::ddtCorr(U, phi)phi = phiHbyA - pEqn.flux(); 
这两个问题并不简单,有错误请指出,一起讨论。
 fvc::interpolate(rAU)*fvc::ddtCorr(U, phi)的作用与PISO的非正交修正有关,大意可能是利用前一次的速度场插值通量场和通量场的差别来对本次的通量场进行修正,对于结构化,完全正交的网格可能不起作用,可参考[9-10],代码在finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.C,其中函数及之间的调用关系见[12]
pEqn.flux()是由fvScalarMatrix返回通量phi = phiHbyA - pEqn.flux(); 是为了得到无散度的通量[5-6]。具体来说,flux()函数作用于fvScalarMatrix,返回的应该是(fvc::interpolate(rUA*fvc::grad(p)) & mesh.Sf(见3.144),这与[1]中的用法是相符合的,因为在porousMultiphaseFoam中速度本身没有输运(达西定律),也就没有H(U)相关的项(见3.137。在实际实现中,flux()返回的值可能还包含非正交修正相关的项[5]。flux()具体的返回值与方程中的fvm项有关,更具体的讨论见[13][14],可以多看几个求解器来build confidence。
对于本文开头,不同求解器中对pEqn.flux()貌似不同的用法,可以参照[11]来理解,依据具体的求解器进行分析。

reference:
  1. porousMultiphaseFoam详解(一)& OpenFOAM求解器开发入门(一):求解策略
  2. 东岳流体技术文档:何为通量
  3. OpenFOAM guide/The SIMPLE algorithm in OpenFOAM
  4. OpenFOAM guide/The PISO algorithm in OpenFOAM
  5. http://www.cfd-online.com/Forums/openfoam/114933-description-flux-method.html
  6. http://www.cfd-online.com/Forums/openfoam-programming-development/83638-phi-peqn-flux-vs-linearinterpolate-u-mesh-sf.html(关于PISO的讨论)
  7. 东岳流体技术文档:icoFoam 详解
  8. 东岳流体技术文档:potentialFoam 详解
  9. http://www.cfd-online.com/Forums/openfoam-solving/60096-ddtphicorr.html
  10. http://www.cfd-online.com/Forums/openfoam/84210-fvc-ddtphicorr-rua-u-phi.html
  11. http://www.cfd-online.com/Forums/openfoam-programming-development/89144-wrong-operator-transonic-pressure-equation.html
  12. http://www.cfd-online.com/Forums/openfoam-programming-development/70444-programming-questions.html
  13. http://www.cfd-online.com/Forums/openfoam-solving/109440-peqn-flux.html
  14. http://www.cfd-online.com/Forums/openfoam-programming-development/68027-different-implementation-solving-peqn.html
  15. http://www.cfd-online.com/Forums/openfoam-solving/90435-interpolation-pressure.html
  16. http://www.mathematik.uni-dortmund.de/~kuzmin/cfdintro/lecture11.pdf
注:本文公式直接取自Jasak的博士论文,并保留了标号以供查询。


敬请期待,OpenFOAM求解器开发入门(三):组分输运



有CFD或者of相关的问题,欢迎联系一起讨论,邮箱:wangyan14@mails.tsinghua.edu.cn。
请使用单位邮箱。
所需信息:
姓名
具体单位
职务职称
问题的具体描述,可以的话请带上相关的文献和代码等附件

0

阅读 评论 收藏 转载 喜欢 打印举报/Report
  • 评论加载中,请稍候...
发评论

    发评论

    以上网友发言只代表其个人观点,不代表新浪网的观点或立场。

      

    新浪BLOG意见反馈留言板 电话:4000520066 提示音后按1键(按当地市话标准计费) 欢迎批评指正

    新浪简介 | About Sina | 广告服务 | 联系我们 | 招聘信息 | 网站律师 | SINA English | 会员注册 | 产品答疑

    新浪公司 版权所有