coalChemistryFoam part 1: main function

源代码：coalChemistryFoam.C
为了提高看代码的效率，先把英文的粘过来，之后抽空补充中文翻译

主程序

先来看coalChemistryFoam.C
首先给出了程序的简要说明
程序名：coalChemistryFoam
类别：grpLagrangianSolvers
简介：Transient solver for compressible, turbulent flow, with coal and limestone particle clouds, an energy source, and combustion.

头文件

#include "fvCFD.H"
#include "turbulenceModel.H"
#include "basicThermoCloud.H"
#include "coalCloud.H"
#include "psiCombustionModel.H"
#include "fvIOoptionList.H"
#include "radiationModel.H"
#include "SLGThermo.H"
#include "pimpleControl.H"

以上是v2.3.1，在v6和V1806中增加了如下头文件：

#include "pressureControl.H"
#include "localEulerDdtScheme.H"
#include "fvcSmooth.H"

main function

  • addCheckCaseOptions.H includes predefined arguments to check case set-up.
  • setRootCase.H will create an object from argList class to register and output informations according the case settings. argList is a very basic class defined in argList.H, and its function is to extract command arguments and options from the inputs.
  • createTime.H will create the object runtime from the Time class. This class is used to control time loops during the calculation.
  • createMesh.H will create the object mesh from the fvMesh class. fvMesh is a multiple inheritence class, and it inherits from polyMesh class and other classes related to finite volume method.
  • createControl.H In createControl.H three different solver controllers are defined (SIMPLE, PISO and PIMPLE). In coalChemistryFoam, PIMPLE algorithm is used. These algorithm control class will supply convergence information checks for the calculation loop. The implementation of pressure correction algorithm is usually in the pEqn.H file.
  • creatTimeControls.H Three variables adjustTimeStep, maxCo and maxDeltaT will be created and initialized from the controlDict
  • creatFields.H In this file, the variables are created as fields data object. Two pointers, i.e. combustion and turbulence, are defined by autoPrt template type. The pointer created by New function is using for the ”Run Time Selection(RTS)” method, it allows to select the submodel class used in template after compiling of the source code. We will look into more detail in this file in the following sections.

After the object and variable declaration is almost done, the turbulent model is checked by the following code:

turbulence->validate();

Before entering the time loop, set initial time step and calculate Courant number. In general, as the courant number changes from small to large, the convergence speed gradually increases, but the stability gradually decreases.

if (!LTS)
    {
        #include "compressibleCourantNo.H"
        #include "setInitialDeltaT.H"
    }

 Info<< "\nStarting time loop\n" << endl;

    while (runTime.run())
    {
        #include "readTimeControls.H"

        if (LTS)
        {
            #include "setRDeltaT.H"
        }
        else
        {
            #include "compressibleCourantNo.H"
            #include "setDeltaT.H"
        }

        runTime++;

        Info<< "Time = " << runTime.timeName() << nl << endl;

LTS is short for local time stepping. When the reactions and radiation model is activated, LTS is usually used. It needs to be specified according to the solver itself, so we can find the file ”setDeltaT.H” in the solver’s directory.

        rhoEffLagrangian = coalParcels.rhoEff() + limestoneParcels.rhoEff();
        pDyn = 0.5*rho*magSqr(U);

        coalParcels.evolve();

        limestoneParcels.evolve();

Then, the solver calculates and updates particles’ parameter by calling the evolve() function in the parcel class, the source terms from the particle to the gas phase equations will also be calculated. rhoEffLagrangian and pDyn is not used in the calculation. They are only calculated for post-process.
coalParcels is a coalCould class, as defined in createClouds.H and envolve() is a function of Class ThermoCloud.H

Info<< "\nConstructing coal cloud" << endl;
coalCloud coalParcels
(
    "coalCloud1",
    rho,
    U,
    g,
    slgThermo
);

Info<< "\nConstructing limestone cloud" << endl;
basicThermoCloud limestoneParcels
(
    "limestoneCloud1",
    rho,
    U,
    g,
    slgThermo
);

Next, solve the governing eqautions for the gas phase, and correct the pressure using PIMPLE algorithm.

        #include "rhoEqn.H"

        // --- Pressure-velocity PIMPLE corrector loop
        while (pimple.loop())
        {
            #include "UEqn.H"
            #include "YEqn.H"
            #include "EEqn.H"

            // --- Pressure corrector loop
            while (pimple.correct())
            {
                #include "pEqn.H"
            }

Then, solve the turbulence based on pimple algorithm.

            if (pimple.turbCorr())
            {
                turbulence->correct();
            }
        }

At last, the density of gas phase is updated according to the "equation of state", and write data before entering the next time step.

        rho = thermo.rho();

        runTime.write();

        runTime.printExecutionTime(Info);
    }

    Info<< "End\n" << endl;

    return 0;
}

方程求解

下面以大涡模拟中气相控制方程为例

质量守恒方程

质量守恒方程

rhoEqn.H

{
    fvScalarMatrix rhoEqn
    (
        fvm::ddt(rho)
      + fvc::div(phi)
      ==
        coalParcels.Srho(rho)
      + fvOptions(rho)
    );

    rhoEqn.solve();

    fvOptions.correct(rho);
}

【OF中的方程离散】微分方程是通过离散从而转化为代数方程组进行求解。对于一个代数方程组主要有两部分：系数矩阵A，右边值b。构造成类似Ax＝b的代数方程形式。方程的离散有显式和隐式，隐式离散得到方程稀疏矩阵A的，显式离散得到右边的值b。一旦通过某种离散形式得到A和b就可以求解线性方程组了。
OpenFOAM的所有的显式离散都在fvc空间中，所以以fvc开头的都为显式，比如fvc::div(phi)。隐式离散在fvm空间中，所有的隐式离散都是以fvm开头，比如fvm::ddt, fvm::div(phi,U)等。所以，以fvm开头的项声称A，以fvc开头的项生成b。
对于质量守恒方程，fvm::ddt(rho)表示ddt(rho)是待求解的未知量，fvc::div(phi)表示div(phi)是已知量，即上一时刻的值。
这里可压缩流中，phi定义为rho*U，见compressibleCreatePhi_H：

surfaceScalarField phi
(
    IOobject
    (
        "phi",
        runTime.timeName(),
        mesh,
        IOobject::READ_IF_PRESENT,
        IOobject::AUTO_WRITE
    ),
    linearInterpolate(rho*U) & mesh.Sf()
);
  • coalParcels.Srho(rho)表示质量源项，Srho是ReactingCloud.H类的函数
  • fvOptions(rho)表示附加源项，OpenFOAM 2.2.0: fvOptions，在system/fvOptions中设置，可用作点火源

动量方程

动量守恒方程

UEqn.H

    MRF.correctBoundaryVelocity(U);

    fvVectorMatrix UEqn
    (
        fvm::ddt(rho, U) + fvm::div(phi, U)
      + MRF.DDt(rho, U)
      + turbulence->divDevRhoReff(U)
     ==
        rho()*g
      + coalParcels.SU(U)
      + limestoneParcels.SU(U)
      + fvOptions(rho, U)
    );

    UEqn.relax();

    fvOptions.constrain(UEqn);

    if (pimple.momentumPredictor())
    {
        solve(UEqn == -fvc::grad(p));

        fvOptions.correct(U);
        K = 0.5*magSqr(U);
    }
  • MRF is for rotating reference frame, MRF.DDt will return the contribution from the rotating reference frame effects.
  • turbulence是creatFields.H中定义的湍流模型：
    Info<< "Creating turbulence model\n" << endl;
    autoPtr<compressible::turbulenceModel> turbulence
    (
        compressible::turbulenceModel::New
        (
            rho,
            U,
            phi,
            thermo
        )
    );
  • divDevRhoReff(U) function is determined by the turbulence model. It will return the source term for the momentum equation. 在RANS下如果couplingFactor为0（事实上，在ReynoldsStress.C中已经将其初始化为0， 它将对Reynold应力项的显式计算进行加权，详细信息可以在定义中找到）ReynoldsStress.C，在它的定义为：
            fvc::laplacian
            (
                this->alpha_*rho*this->nut(),
                U,
                "laplacian(nuEff,U)"
            )
          + fvc::div(this->alpha_*rho*R_)
          - fvc::div(this->alpha_*rho*this->nu()*dev2(T(fvc::grad(U))))
          - fvm::laplacian(this->alpha_*rho*this->nuEff(), U)

LES下（Smagorinsky模型），它的定义为linearViscousStress.C

Foam::linearViscousStress<BasicTurbulenceModel>::divDevRhoReff
(
    volVectorField& U
) const
{
    return
    (
      - fvc::div((this->alpha_*this->rho_*this->nuEff())*dev2(T(fvc::grad(U))))
      - fvm::laplacian(this->alpha_*this->rho_*this->nuEff(), U)
    );
}

最后，考虑压力梯度项，利用Pimple循环：

    if (pimple.momentumPredictor())
    {
        solve(UEqn == -fvc::grad(p));

        fvOptions.correct(U);
        K = 0.5*magSqr(U);
    }

完成动量方程求解

组分方程

image.png

YEqn.H

In the YEqn.H file, a convectionScheme, which is named mvConvection, is defined at the beginning. It is created on free store in random-access memory (RAM). When the number of gas species is large, it enhances the computational efficiency.

tmp<fv::convectionScheme<scalar>> mvConvection
(
    fv::convectionScheme<scalar>::New
    (
        mesh,
        fields,
        phi,
        mesh.divScheme("div(phi,Yi_h)")
    )
);

Before calculate the species transport equation, the gas phase combustion is calculated and the heat effects of the reactions is stored in the parameter Qdot().

    combustion->correct();
    Qdot = combustion->Qdot();

A variable Yt is initialized

volScalarField Yt(0.0*Y[0]);

and used in the calculation loop to make sure all the species mass fraction together is 1. The error is made up by the inert gas. The transport equation is defined as:

              fvScalarMatrix YiEqn
            (
                fvm::ddt(rho, Yi)
              + mvConvection->fvmDiv(phi, Yi)
              - fvm::laplacian(turbulence->muEff(), Yi)
              ==
                coalParcels.SYi(i, Yi)
              + combustion->R(Yi)
              + fvOptions(rho, Yi)
            );

能量方程

能量守恒方程

EEqn.H

    fvScalarMatrix EEqn
    (
        fvm::ddt(rho, he) + mvConvection->fvmDiv(phi, he)
      + fvc::ddt(rho, K) + fvc::div(phi, K)
      + (
            he.name() == "e"
          ? fvc::div
            (
                fvc::absolute(phi/fvc::interpolate(rho), U),
                p,
                "div(phiv,p)"
            )
          : -dpdt
        )
      - fvm::laplacian(turbulence->alphaEff(), he)
     ==
        rho*(U&g)
      + Qdot
      + coalParcels.Sh(he)
      + limestoneParcels.Sh(he)
      + radiation->Sh(thermo, he)
      + fvOptions(rho, he)
    );

压力修正

pimple的基本思想是：将每个时间步长内用simple稳态算法求解（也就是将每个时间步内看成稳态流动），时间步长的步进用piso算法来完成。
在有限容积离散中，时间项的离散仍采用的差分格式，这样做可以得到某个时间点的流场信息，而非某个时间步长的内的平均值。采用传统的piso算法求解变化较快的流动的时候，需要的时间步长较小（因为相邻两个时间点的流场不能差别太大，否则会发散），这样会造成求解的某种流动需要的耗时过长。 pimple算法将每个时间步长内看成一种稳态的流动（采用亚松驰来解决相邻两个时间段变化大的情况），当按照稳态的求解器求解到一定的时候，则采用标准的piso做最后一步求解。

子模型

OpenFOAM通过RTS定义子模型，OpenFOAM中的RunTime Selection（RTS）机制——TypeName

createFields.H

createFields主要是用来加载各种子模型的，它包含了若干.H文件：

#include "createRDeltaT.H" //为RTS模型创建时间步长
#include "readGravitationalAcceleration.H"//重力加速度
#include "compressibleCreatePhi.H"//可压缩流
#include "createDpdt.H"
#include "createK.H"
#include "createMRF.H"
#include "createClouds.H"//多相流模型
#include "createRadiationModel.H"//辐射模型
#include "createFvOptions.H"

参考文献

  1. Zhang, J.: Modifying coalChemistryFoam for dense gas-solid simulation. In Proceedings of CFD with OpenSource Software, 2018, Edited by Nilsson. H., http://dx.doi.org/10.17196/OS_CFD#YEAR_2018
  2. 苏军伟，OpenFOAM中的神奇方程定义方式的背后. http://blog.sina.com.cn/s/blog_5fdfa7e60100d83c.html
  3. 苏军伟，pimple算法简述. http://blog.sina.com.cn/s/blog_5fdfa7e60100fbb1.html
  4. 王伟，OpenFOAM中的RunTime Selection（RTS）机制——TypeName，http://youmengtian.coding.me/2017/03/25/runtimeselection/
最后编辑于 ：
©著作权归作者所有,转载或内容合作请联系作者
  • 人面猴
    序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
    沈念sama阅读 216,496评论 6赞 501
  • 死咒
    序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
    沈念sama阅读 92,407评论 3赞 392
  • 救了他两次的神仙让他今天三更去死
    文/潘晓璐 我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
    开封第一讲书人阅读 162,632评论 0赞 353
  • 道士缉凶录：失踪的卖姜人
    文/不坏的土叔 我叫张陵，是天一观的道长。 经常有香客问我，道长，这世上最难降的妖魔是什么？ 我笑而不...
    开封第一讲书人阅读 58,180评论 1赞 292
  • 港岛之恋（遗憾婚礼）
    正文 为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
    茶点故事阅读 67,198评论 6赞 388
  • 恶毒庶女顶嫁案：这布局不是一般人想出来的
    文/花漫 我一把揭开白布。 她就那样静静地躺着，像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上，一...
    开封第一讲书人阅读 51,165评论 1赞 299
  • 城市分裂传说
    那天，我揣着相机与录音，去河边找鬼。 笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。 我是一名探鬼主播，决...
    沈念sama阅读 40,052评论 3赞 418
  • 双鸳鸯连环套：你想象不到人心有多黑
    文/苍兰香墨 我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
    开封第一讲书人阅读 38,910评论 0赞 274
  • 万荣杀人案实录
    序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
    沈念sama阅读 45,324评论 1赞 310
  • 护林员之死
    正文 独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 37,542评论 2赞 332
  • 白月光启示录
    正文 我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 39,711评论 1赞 348
  • 活死人
    序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
    沈念sama阅读 35,424评论 5赞 343
  • 日本核电站爆炸内幕
    正文 年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
    茶点故事阅读 41,017评论 3赞 326
  • 男人毒药：我在死后第九天来索命
    文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 31,668评论 0赞 22
  • 一桩弑父案，背后竟有这般阴谋
    文/苍兰香墨 我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,823评论 1赞 269
  • 情欲美人皮
    我被黑心中介骗来泰国打工， 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。 一个月前我还...
    沈念sama阅读 47,722评论 2赞 368
  • 代替公主和亲
    正文 我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
    茶点故事阅读 44,611评论 2赞 353

推荐阅读更多精彩内容

  • rljs
    rljs by sennchi Timeline of History Part One The Cognitiv...
    sennchi阅读 7,322评论 0赞 10
  • 2003-2014考研真题
    **2014真题Directions:Read the following text. Choose the be...
    又是夜半惊坐起阅读 9,476评论 0赞 23
  • 2018-06-11
    mean to add the formatted="false" attribute?.[ 46% 47325/...
    ProZoom阅读 2,696评论 0赞 3
  • 起信论7
    2016-11-21 05:20:11 原文地址：起信论7_yc791022_新浪博客 复次，此真如者，依言说分别...
    露电梦幻泡影阅读 722评论 0赞 0
  • 在雨中奔跑的感觉
    暴雨后接着就是细雨绵绵，梅雨季节里云儿爬的很低，一挥衣袖伸手就能撸到云朵，不知道是晚炊烟云还是从天空中飞翔过来的云...
    象浦阅读 279评论 0赞 1