chemFoam

概述

基本介绍

1. 定义：化学问题求解器，设计用于单细胞案例，提供与其他化学求解器的比较，使用单细胞网格，并从初始条件创建字段。

求解策略

$$
\frac {\partial(\rho Y)}{\partial t}= \omega_i(Y_i,T)
$$

$$
h=u_0+\frac p \rho+\int_0^t\frac q \rho d\tau=u+\frac p \rho
$$

$$
p = \sum p _ { i } = \sum n _ { i } \frac { R T } { V } = \sum \frac { m _ { i } } { M _ { i } } \frac { R T } { V } = \sum \frac { y _ { i } } { M _ { i } } \underbrace { \frac { m _ { t o t } } { V } } R T
$$

1)求解化学:目的是得到所涉及的每种物质的反应速率和化学反应产生的热量

2)求解物种方程:求新时间步长的物种浓度

3)求解能量方程:这里我们得到了新的时间步长的焓(除以热力学和温度)

4)求解压力方程:计算Rspecific，p和rho需要计算焓

invW += Y[i][0]/specieData[i].W();
Rspecific[0] = 1000.0*constant::physicoChemical::R.value()*invW;
p[0] = rho0*Rspecific[0]*thermo.T()[0];
rho[0] = rho0;

知识点

• argList::noParallel()：禁止并行
• 量纲：dimensionedScalar(“Ydefault”, dimless, 1)

//- Construct given a name, a value and its dimensionSet.
//dimensioned(const word&, const dimensionSet&, const Type);
//extern const dimensionSet dimless;

常见的量纲

const dimensionSet dimMass(1, 0, 0, 0, 0, 0, 0);
const dimensionSet dimLength(0, 1, 0, 0, 0, 0, 0);
const dimensionSet dimTime(0, 0, 1, 0, 0, 0, 0);
const dimensionSet dimTemperature(0, 0, 0, 1, 0, 0, 0);
const dimensionSet dimMoles(0, 0, 0, 0, 1, 0, 0);
const dimensionSet dimCurrent(0, 0, 0, 0, 0, 1, 0);
const dimensionSet dimLuminousIntensity(0, 0, 0, 0, 0, 0, 1);

const dimensionSet dimArea(sqr(dimLength));
const dimensionSet dimVolume(pow3(dimLength));
const dimensionSet dimVol(dimVolume);

const dimensionSet dimVelocity(dimLength/dimTime);
const dimensionSet dimAcceleration(dimVelocity/dimTime);

const dimensionSet dimDensity(dimMass/dimVolume);
const dimensionSet dimForce(dimMass*dimAcceleration);
const dimensionSet dimEnergy(dimForce*dimLength);
const dimensionSet dimPower(dimEnergy/dimTime);

const dimensionSet dimPressure(dimForce/dimArea);
const dimensionSet dimGasConstant(dimEnergy/dimMass/dimTemperature);
const dimensionSet dimSpecificHeatCapacity(dimGasConstant);
const dimensionSet dimViscosity(dimArea/dimTime);

• refCast：类型转换模板函数

template<class To, class From>
inline To& refCast(From& r)
{
    try
    {
        return dynamic_cast<To&>(r);
    }
    catch (std::bad_cast)
    {
        FatalErrorIn("refCast<To>(From&)")
            << "Attempt to cast type " << r.type()
            << " to type " << To::typeName
            << abort(FatalError);

        return dynamic_cast<To&>(r);	//此部分在c++基础中介绍
    }
}

• forAll:

#define forAll(list, i) \
      for (Foam::label i=0; i<(list).size(); i++)

• 总焓=对温度压力敏感的焓＋化学反应生成焓

//Sensible enthalpy [J/kg]. 
   /*
   Foam::species::thermo<Thermo, Type>::Hs(const scalar p, const scalar T) const
  {
      return this->hs(p, T)/this->W();
  }
  
  Hs = Ha(p,T) - Hc; //Hs:Sensible enthalpy; Ha:absolute enthalpy; Hc:chemical enthalpy 
  
  */

待解决

• polyPatch(name, size, start, index, bm, patchType)：如何构成边界面的
• T.write()
• Sh：源项的计算

代码详读

Application
    chemFoam

Description
    Solver for chemistry problems
    - designed for use on single cell cases to provide comparison against
      other chemistry solvers
    - single cell mesh created on-the-fly
    - fields created on the fly from the initial conditions

\*---------------------------------------------------------------------------*/

#include "fvCFD.H"						//有限体积法头文件
#include "psiReactionThermo.H"			//基于可压缩性$\phi$的反应混合热物理模型库
#include "turbulenceModel.H"			//湍流模型
#include "psiChemistryModel.H"			//可压缩热力学化学模型
#include "chemistrySolver.H"			//求解化学反应的虚基类
#include "OFstream.H"					//输出文件stream
#include "thermoPhysicsTypes.H"			//热物理模型的类别定义
#include "basicMultiComponentMixture.H"	//多组分混合物。提供质量分数场和辅助函数列表，以查询											混合物成分。
#include "cellModeller.H"				//单元格模型的静态集合，以及查找它们的方法。

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

int main(int argc, char *argv[])
{
    argList::noParallel();		//禁止并行

    #include "setRootCase.H"
    #include "createTime.H"
------------------------------createSingleCellMesh.H-------------------------
    //#include "createSingleCellMesh.H"	
    //建立一个1*1*1的体单元并划分一个hex网格
    Info<< "Constructing single cell mesh" << nl << endl;
    
    labelList owner(6, label(0));	//数组[6]，每个赋值0，不是特别理解
    labelList neighbour(0);
    
    pointField points(8);
    points[0] = vector(0, 0, 0);
    points[1] = vector(1, 0, 0);
    points[2] = vector(1, 1, 0);
    points[3] = vector(0, 1, 0);
    points[4] = vector(0, 0, 1);
    points[5] = vector(1, 0, 1);
    points[6] = vector(1, 1, 1);
    points[7] = vector(0, 1, 1);
    
    const cellModel& hexa = *(cellModeller::lookup("hex"));
    //cellModeller::lookup: Look up a model by name and return a (const     cellModeller) pointer to the model or NULL
    faceList faces = hexa.modelFaces();	//返回模型的面列表
    
    fvMesh mesh
    (
        IOobject
        (
            fvMesh::defaultRegion,
            runTime.timeName(),
            runTime,
            IOobject::READ_IF_PRESENT
        ),
        xferMove<Field<vector> >(points),	//传输列表中的参数
//Construct by transferring the contents of the \a arg
/*
template<class T>
inline Foam::Xfer<T> Foam::xferMove(T& t)
{
    return Foam::Xfer<T>(t, true);
}

template<class T>
inline Foam::Xfer<T>::Xfer(T& t, bool allowTransfer)
:
    ptr_(new T)
{
    if (allowTransfer)
    {
        ptr_->transfer(t);
    }
    else
    {
        ptr_->operator=(t);
    }
}

template<class T>
void Foam::List<T>::transfer(List<T>& a)
{
    if (this->v_) delete[] this->v_;
    this->size_ = a.size_;	//列表中元素个数
    this->v_ = a.v_;	//列表中矢量

    a.size_ = 0;
    a.v_ = 0;
}
*/
    faces.xfer(),
/*
template<class T>
inline Foam::Xfer<Foam::List<T> > Foam::List<T>::xfer()
{
    return xferMove(*this);
}
*/
        owner.xfer(),
        neighbour.xfer()
    );
    
    List<polyPatch*> patches(1);	//创建一个边界
    
    patches[0] = new emptyPolyPatch //用于二维几何的空的前后面
    (
        "boundary",
        6,
        0,
        0,
        mesh.boundaryMesh(),
        emptyPolyPatch::typeName
    );
    
/*
//- Construct from components
emptyPolyPatch
(
    const word& name,
    const label size,
    const label start,
    const label index,
    const polyBoundaryMesh& bm,
    const word& patchType
);
:
	polyPatch(name, size, start, index, bm, patchType)
{}
*/
    mesh.addFvPatches(patches);
------------------------------createSingleCellMesh.H-------------------------
    
------------------------------createFields.H---------------------------------
    //#include "createFields.H"
    if (mesh.nCells() != 1)
    {
        FatalErrorIn(args.executable())
            << "Solver only applicable to single cell cases"
            << exit(FatalError);
    }

    Info<< "Reading initial conditions.\n" << endl;
    IOdictionary initialConditions	
    (
        IOobject
        (
            "initialConditions",
            runTime.constant(),
            runTime,
            IOobject::MUST_READ_IF_MODIFIED,
            IOobject::NO_WRITE
        )
    );

    scalar p0 = readScalar(initialConditions.lookup("p"));	//初始压力
    scalar T0 = readScalar(initialConditions.lookup("T"));	//初始温度

------------------------------createBaseFields.H-----------------------------       
    #include "createBaseFields.H"
// write base thermo fields - not registered since will be re-read by
// thermo package

Info<< "Creating base fields for time " << runTime.timeName() << endl;
{
    volScalarField Ydefault
    (
        IOobject
        (
            "Ydefault",
            runTime.timeName(),
            mesh,
            IOobject::READ_IF_PRESENT,
            IOobject::NO_WRITE,
            false
        ),
        mesh,
        dimensionedScalar("Ydefault", dimless, 1)
        //- Construct given a name, a value and its dimensionSet.
        //dimensioned(const word&, const dimensionSet&, const Type);
        //extern const dimensionSet dimless;        
    );

    Ydefault.write();	//输出Ydefault

    volScalarField p
    (
        IOobject
        (
            "p",
            runTime.timeName(),
            mesh,
            IOobject::READ_IF_PRESENT,
            IOobject::NO_WRITE,
            false
        ),
        mesh,
        dimensionedScalar("p", dimPressure, p0)
    );

    p.write();

    volScalarField T
    (
        IOobject
        (
            "T",
            runTime.timeName(),
            mesh,
            IOobject::READ_IF_PRESENT,
            IOobject::NO_WRITE,
            false
        ),
        mesh,
        dimensionedScalar("T", dimTemperature, T0)
    );

    T.write();
}
------------------------------createBaseFields.H----------------------------- 
    
    Info<< nl << "Reading thermophysicalProperties" << endl;
    autoPtr<psiChemistryModel> pChemistry(psiChemistryModel::New(mesh));
//psiChemistryModel.H:Chemistry model for compressibility-based thermodynamics
/*
Foam::autoPtr<Foam::psiChemistryModel> Foam::psiChemistryModel::New
(
    const fvMesh& mesh
)
{
    return basicChemistryModel::New<psiChemistryModel>(mesh);
}

//- Generic New for each of the related chemistry model
template<class Thermo>
static autoPtr<Thermo> New(const fvMesh&);

Foam::basicChemistryModel::basicChemistryModel(const fvMesh& mesh)
:
    IOdictionary
    (
        IOobject
        (
            "chemistryProperties",
            mesh.time().constant(),
            mesh,
            IOobject::MUST_READ_IF_MODIFIED,
            IOobject::NO_WRITE
        )
    ),
    mesh_(mesh),
    chemistry_(lookup("chemistry")),
    deltaTChemIni_(readScalar(lookup("initialChemicalTimeStep"))),
    deltaTChem_
    (
        IOobject
        (
            "deltaTChem",
            mesh.time().constant(),
            mesh,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        mesh,
        dimensionedScalar("deltaTChem0", dimTime, deltaTChemIni_)
    )
{}
*/    
        
    psiChemistryModel& chemistry = pChemistry();
    scalar dtChem = refCast<const psiChemistryModel>(chemistry).deltaTChem()[0];
    //初始时刻化学反应的时间步
    //refCast：类型转换模板函数

    psiReactionThermo& thermo = chemistry.thermo();
    thermo.validate(args.executable(), "h");

    basicMultiComponentMixture& composition = thermo.composition();
    //- Return the composition of the multi-component mixture
    //- 物种名+质量分数
    PtrList<volScalarField>& Y = composition.Y();

    volScalarField rho
    (
        IOobject
        (
            "rho",
            runTime.timeName(),
            runTime,
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
        thermo.rho()	//return p_*psi_;
    );

    volScalarField& p = thermo.p();	//return p_

    volScalarField Rspecific	//普适气体常数
    (
        IOobject
        (
            "Rspecific",
            runTime.timeName(),
            runTime,
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
        mesh,
        dimensionedScalar
        (
            "zero",
            dimensionSet(dimEnergy/dimMass/dimTemperature),
            0.0
        )
    );

    volVectorField U	//讲道理，感觉也没用到速度场呀，建了没吊用
    (
        IOobject
        (
            "U",
            runTime.timeName(),
            runTime,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        mesh,
        dimensionedVector("zero", dimVelocity, vector::zero),
        p.boundaryField().types()
    );
------------------------------createPhi.H------------------------------------
    //#include "createPhi.H"
    #ifndef createPhi_H
    #define createPhi_H

    Info<< "Reading/calculating face flux field phi\n" << endl;

    surfaceScalarField phi
    (
        IOobject
        (
            "phi",
            runTime.timeName(),
            mesh,
            IOobject::READ_IF_PRESENT,
            IOobject::AUTO_WRITE
        ),
        linearInterpolate(U) & mesh.Sf()	//线性插值
    );

    #endif
------------------------------createPhi.H------------------------------------

    Info << "Creating turbulence model.\n" << endl;
    autoPtr<compressible::turbulenceModel> turbulence
    (
        compressible::turbulenceModel::New
        (
            rho,
            U,
            phi,
            thermo
        )
    );
/*
autoPtr<turbulenceModel> turbulenceModel::New
(
    const volScalarField& rho,
    const volVectorField& U,
    const surfaceScalarField& phi,
    const fluidThermo& thermophysicalModel,
    const word& turbulenceModelName
)
{
    // get model name, but do not register the dictionary
    // otherwise it is registered in the database twice
    const word modelType
    (
        IOdictionary
        (
            IOobject
            (
                "turbulenceProperties",
                U.time().constant(),
                U.db(),
                IOobject::MUST_READ_IF_MODIFIED,
                IOobject::NO_WRITE,
                false
            )
        ).lookup("simulationType")
    );

    Info<< "Selecting turbulence model type " << modelType << endl;

    turbulenceModelConstructorTable::iterator cstrIter =
        turbulenceModelConstructorTablePtr_->find(modelType);

    if (cstrIter == turbulenceModelConstructorTablePtr_->end())
    {
        FatalErrorIn
        (
            "turbulenceModel::New(const volScalarField&, "
            "const volVectorField&, const surfaceScalarField&, "
            "fluidThermo&, const word&)"
        )   << "Unknown turbulenceModel type "
            << modelType << nl << nl
            << "Valid turbulenceModel types:" << endl
            << turbulenceModelConstructorTablePtr_->sortedToc()
            << exit(FatalError);
    }

    return autoPtr<turbulenceModel>
    (
        cstrIter()(rho, U, phi, thermophysicalModel, turbulenceModelName)
    );
}
*/

    OFstream post(args.path()/"chemFoam.out");
    post<< "# Time" << token::TAB << "Temperature [K]" << token::TAB
        << "Pressure [Pa]" << endl;
------------------------------createFields.H---------------------------------
------------------------------readInitialConditions.H------------------------   
    //#include "readInitialConditions.H"
    word constProp(initialConditions.lookup("constantProperty"));
    if ((constProp != "pressure") && (constProp != "volume"))
    {
        FatalError << "in initialConditions, unknown constantProperty type "
            << constProp << nl << " Valid types are: pressure volume."
            << abort(FatalError);
    }

    word fractionBasis(initialConditions.lookup("fractionBasis"));
    if ((fractionBasis != "mass") && (fractionBasis != "mole"))
    {
        FatalError << "in initialConditions, unknown fractionBasis type " << nl
            << "Valid types are: mass or mole."
            << fractionBasis << abort(FatalError);
    }

    label nSpecie = Y.size();	//PtrList<volScalarField>& Y = composition.Y();
    PtrList<gasHThermoPhysics> specieData(Y.size());
/*
typedef
     sutherlandTransport
     <
         species::thermo
         <
             janafThermo
             <
                 perfectGas<specie>
             >,
             sensibleEnthalpy
         >
     > gasHThermoPhysics;
*/        
    forAll(specieData, i)
    {
        specieData.set
        //inline autoPtr<T> set(const label, const autoPtr<T>&);
        //Set element. Return old element (can be NULL).
        //No checks on new element.
        (
            i,
            new gasHThermoPhysics
            (
                dynamic_cast<const reactingMixture<gasHThermoPhysics>&>
                    (thermo).speciesData()[i]
            )
        );
    }

    scalarList Y0(nSpecie, 0.0);
    scalarList X0(nSpecie, 0.0);

    dictionary fractions(initialConditions.subDict("fractions")); //subDict注意下
    if (fractionBasis == "mole")
    {
        forAll(Y, i)
        {
            const word& name = Y[i].name();
            if (fractions.found(name))
            {
                X0[i] = readScalar(fractions.lookup(name));
            }
        }

        scalar mw = 0.0;
        const scalar mTot = sum(X0);	//显示摩尔数值加和求摩尔分数
        forAll(Y, i)
        {
            X0[i] /= mTot;	//计算摩尔分数
            mw += specieData[i].W()*X0[i];	//W():摩尔质量，单位kg/kmol
        }

        forAll(Y, i)
        {
            Y0[i] = X0[i]*specieData[i].W()/mw;
        }
    }
    else  // mass fraction
    {
        forAll(Y, i)
        {
            const word& name = Y[i].name();
            if (fractions.found(name))
            {
                Y0[i] = readScalar(fractions.lookup(name));
            }
        }

        scalar invW = 0.0;
        const scalar mTot = sum(Y0);
        forAll(Y, i)
        {
            Y0[i] /= mTot;
            invW += Y0[i]/specieData[i].W();
        }
        const scalar mw = 1.0/invW;

        forAll(Y, i)
        {
            X0[i] = Y0[i]*mw/specieData[i].W();
        }
    }

    scalar h0 = 0.0;
    forAll(Y, i)
    {
        Y[i] = Y0[i];
        //Y[i]:Return the mass-fraction field for a specie given by index
        h0 += Y0[i]*specieData[i].Hs(p[0], T0);
    }
	//Sensible enthalpy [J/kg]. 
    /*
    Foam::species::thermo<Thermo, Type>::Hs(const scalar p, const scalar T) const
   {
       return this->hs(p, T)/this->W();
   }
   
   Hs = Ha(p,T) - Hc; //Hs:Sensible enthalpy; Ha:absolute enthalpy; Hc:chemical enthalpy 
   
   */        

    thermo.he() = dimensionedScalar("h", dimEnergy/dimMass, h0);
    //Enthalpy/Internal energy [J/kg]
    thermo.correct();
    //更新特性

    rho = thermo.rho();
    scalar rho0 = rho[0];
    scalar u0 = h0 - p0/rho0;
    scalar R0 = p0/(rho0*T0);
    Rspecific[0] = R0;	
    //Rspecific[0] = 1000.0*constant::physicoChemical::R.value()*invW;

    scalar integratedHeat = 0.0;	//the heat supplied to the system. 

    Info << constProp << " will be held constant." << nl
        << " p   = " << p[0] << " [Pa]" << nl
        << " T   = " << thermo.T()[0] << " [K] " << nl
        << " rho = " << rho[0] << " [kg/m3]" << nl
        << endl;
------------------------------readInitialConditions.H------------------------      
   
    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

    Info<< "\nStarting time loop\n" << endl;
    while (runTime.run())
    {
------------------------------readControls.H----------------------------------       
        //#include "readControls.H"
        if (runTime.controlDict().lookupOrDefault("suppressSolverInfo", false))
        {
            lduMatrix::debug = 0;
        }
    
        Switch adjustTimeStep(runTime.controlDict().lookup("adjustTimeStep"));
/*
        //- The various text representations for a switch value.
        //  These also correspond to the entries in names.
        enum switchType
        {
            FALSE   = 0,  TRUE   = 1,
            OFF     = 2,  ON     = 3,
            NO      = 4,  YES    = 5,
            NO_1    = 6,  YES_1  = 7,
            FALSE_1 = 8,  TRUE_1 = 9,
            NONE    = 10, PLACEHOLDER = 11,
            INVALID
        };
*/    
        scalar maxDeltaT(readScalar(runTime.controlDict().lookup("maxDeltaT")));     ------------------------------readControls.H----------------------------------
            
------------------------------setDeltaT.H------------------------------------- 
        //#include "setDeltaT.H"
        if (adjustTimeStep)
        {
            runTime.setDeltaT(min(dtChem, maxDeltaT));
            Info<< "deltaT = " <<  runTime.deltaT().value() << endl;
        }      
------------------------------setDeltaT.H-------------------------------------
        
        runTime++;
        Info<< "Time = " << runTime.timeName() << nl << endl;
------------------------------solveChemistry.H--------------------------------
        //#include "solveChemistry.H"
        dtChem = chemistry.solve(runTime.deltaT().value());
        //solve:
        //- Solve the reaction system for the given time step
        //  and return the characteristic time
/*
 Foam::dimensionedScalar Foam::TimeState::deltaT() const
 {
     return dimensionedScalar("deltaT", dimTime, deltaT_);
 }
*/        
        scalar Sh = chemistry.Sh()()[0]/rho[0];
        //- Return source for enthalpy equation [kg/m/s3]
        integratedHeat += Sh*runTime.deltaT().value();
------------------------------solveChemistry.H--------------------------------  
    
------------------------------YEqn.H------------------------------------------                 //#include "YEqn.H"       
        {
            forAll(Y, specieI)
            {
                volScalarField& Yi = Y[specieI];

                solve
                (
                    fvm::ddt(rho, Yi) - chemistry.RR(specieI),
                    //- Return access to chemical source terms [kg/m3/s]
                    mesh.solver("Yi")
                );
            }
        }
/*
一、solve the matrix using Gaussian elimination with pivoting,
//  returning the solution in the source
template<class Type>
void solve(scalarSquareMatrix& matrix, List<Type>& source);
   
二、RR:basicChemistryModel.H
/*
//- Return const access to chemical source terms [kg/m3/s]
virtual const DimensionedField<scalar, volMesh>& RR
(
      const label i
) const = 0;

三、
const Foam::dictionary& Foam::solution::solver(const word& name) const
{
    if (debug)
    {
        InfoIn("solution::solver(const word&)")
            << "Lookup solver for " << name << endl;
    }

    return solvers_.subDict(name);
}
*/     
------------------------------YEqn.H------------------------------------------
    
------------------------------hEqn.H------------------------------------------
        //#include "hEqn.H"
        {
            volScalarField& h = thermo.he();
        
            if (constProp == "volume")
            {
                h[0] = u0 + p[0]/rho[0] + integratedHeat;
            }
            else
            {
                h[0] = h0 + integratedHeat;
            }
        
            thermo.correct();	//根据焓更新物性
        }
------------------------------hEqn.H------------------------------------------
    
------------------------------pEqn.H------------------------------------------    
        //#include "pEqn.H"
        {
            rho = thermo.rho();
            if (constProp == "volume")
            {
                scalar invW = 0.0;
                forAll(Y, i)
                {
                    invW += Y[i][0]/specieData[i].W();
                }
        
                Rspecific[0] = 1000.0*constant::physicoChemical::R.value()*invW;
        
                p[0] = rho0*Rspecific[0]*thermo.T()[0];
                rho[0] = rho0;
            }
        }
------------------------------pEqn.H------------------------------------------       

------------------------------output.H---------------------------------------- 
        //#include "output.H"
        runTime.write();
    
        Info<< "Sh = " << Sh
            << ", T = " << thermo.T()[0]
            << ", p = " << thermo.p()[0]
            << ", " << Y[0].name() << " = " << Y[0][0]
            << endl;
    
        post<< runTime.value() << token::TAB << thermo.T()[0] << token::TAB
            << thermo.p()[0] << endl;
------------------------------output.H----------------------------------------   

        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
            << nl << endl;
    }

    Info << "Number of steps = " << runTime.timeIndex() << endl;
    Info << "End" << nl << endl;

    return(0);
}


// ************************************************************************* //