ChemFoam

Материал из ru.wiki.laduga.ru
Перейти к: навигация, поиск

chemFoam 

   Решатель для химических задач, предназначенный для использования в случаях с одной ячейкой для
   обеспечения сравнения с другими химическими решателями, использующими одну ячейку
   сетка и поля, созданные из начальных условий.

1.Стратегия решения

Этот решатель обеспечивает отличную отправную точку для тех, кто хочет получить первое впечатление о влиянии химических реакций в уравнениях, описывающих эволюцию концентрации видов, а также эволюцию температуры (то есть уравнение энергии). Поскольку вычислительная область состоит только из одной ячейки, единственным механизмом, влияющим на эволюцию концентрации веществ и температуры, являются химические реакции. Решаемая система уравнений имеет следующий вид:


$ \frac{ \partial \rho Y_i}{\partial t} = \dot \omega_i (Y_i,T) $
(1)
$ h = u_0 + \frac{p}{\rho} + \int_0^t \frac{\dot q}{\rho} d\tau $
(2)
$ p = \frac{ \rho R T}{M_{avg}} $
(3)

$ Y_i, \omega_i, T, p, \rho, \dot q, R, M_{avg}, u_0, $ это массовая доля фракции, скорость реакции, температура, давление, плотность, теплота, выделяемая при сгорании, газовая постоянная, средняя молярная масса и начальная внутренняя энергия соответственно. Заметим, что для заданной энтальпии, давления и плотности можно рассчитать температуру.

Как и большинство решателей, присутствующих в OpenFOAM, этот решатель также следует стратегии сегрегированного решения. Это означает, что для каждой интересующей величины решается одно линейное уравнение и связь между уравнениями достигается явными исходными членами. Кратко резюмируя, решение достигается следующим образом:

1) Решение химии: цель состоит в том, чтобы получить скорость реакции для каждого участвующего вида вещества и тепло, реализуемое химической реакцией;

2) Решение видового уравнения: цель - получить видовую концентрацию на новом временном шаге;

3) Решение уравнения энергии: здесь мы получаем энтальпию (по термодинамике также температуру) на новом временном шаге;

4) Решение уравнения давления: требуется вычислить энтальпию.


Исходный код можно найти в chemFoam.C 


/*---------------------------------------------------------------------------*\
 =========                 |
 \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
  \\    /   O peration     |
   \\  /    A nd           | Copyright (C) 2010-2011 OpenCFD Ltd.
    \\/     M anipulation  |
-------------------------------------------------------------------------------
                           | Copyright (C) 2011-2017 OpenFOAM Foundation
-------------------------------------------------------------------------------
Лицензия
   Этот файл является частью OpenFOAM.

   OpenFOAM - бесплатное программное обеспечение: вы можете распространять его и / или изменять в
   соответствии с условиями GNU General Public License, опубликованной Фондом
   свободного программного обеспечения, либо версии 3 Лицензии, либо
   (по вашему выбору) любой более поздней версии.

   OpenFOAM распространяется в надежде, что он будет полезен, но без
   КАКИХ-ЛИБО ГАРАНТИЙ; даже без подразумеваемой гарантии товарной пригодности или
   ПРИГОДНОСТЬ ДЛЯ ОПРЕДЕЛЕННОЙ ЦЕЛИ. Смотрите стандартную общественную лицензию GNU
   для получения более подробной информации.

   Для получения более подробной информации. Вы должны были получить копию лицензии GNU General Public License
   вместе с OpenFOAM. Если нет, то смотрите <http://www.gnu.org/licenses/>.

Приложение
   chemFoam

Группа
   grpCombustionSolvers

Описание
   Решатель для химических задач, предназначенный для использования в случаях с одной ячейкой для
   обеспечения сравнения с другими химическими решателями, использующими одну ячейку
   сетки и поля, созданные из начальных условий.

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

#include "fvCFD.H"
#include "psiReactionThermo.H"
#include "BasicChemistryModel.H"
#include "reactingMixture.H"
#include "chemistrySolver.H"
#include "OFstream.H"
#include "thermoPhysicsTypes.H"
#include "basicSpecieMixture.H"
#include "hexCellFvMesh.H"

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

int main(int argc, char *argv[])
{
    argList::addNote
    (
        "Решатель химических задач для использования в случаях с одной ячейкой"
        " для обеспечения сравнения с другими химическими решателями"
    );

    argList::noParallel();

    #define CREATE_MESH createSingleCellMesh.H
    #define NO_CONTROL
    #include "postProcess.H"

    #include "setRootCaseLists.H"
    #include "createTime.H"
    #include "createSingleCellMesh.H"
    #include "createFields.H"
    #include "createFieldRefs.H"
    #include "readInitialConditions.H"
    #include "createControls.H"

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

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

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

        #include "setDeltaT.H"

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

        #include "solveChemistry.H"
        #include "YEqn.H"
        #include "hEqn.H"
        #include "pEqn.H"

        #include "output.H"

        runTime.printExecutionTime(Info);
    }

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

    return 0;
}



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

2. Химические реакции

2.1 Химические уравнения

Простейшая химическая реакция с участием $ K $ видов веществ в $ I $ реакциях может быть выражена следующим образом (смотрите также https://www3.nd.edu/~powers/ame.60636/chemkin2000.pdf or the book of [1]):

$ \sum_{k=1}^K \nu_{ki}' \chi_k \rightleftharpoons \sum_{k=1}^K \nu_{ki}'' \chi_k \ \ \ i=1...I $
(4)

$ \nu_{ki}', \nu_{ki}'' $ and $ \chi_k $ являются прямым, обратным стехиометрическими коэффициентами и химическим символом вида $ k $, соответственно.

Чтобы сделать значение более понятным для начинающих пользователей CFD для моделирования процесса горения, механизм Зельдовича для окисления N0 (см. https://en.wikipedia.org/wiki/Zeldovich_mechanism или книга [2]) используется в качестве примера:

$ N_2 + O \rightleftharpoons NO + N $
$ N + O_2 \rightleftharpoons NO + O $

На основе приведённого выше примера, прямой $ \nu_{ki}' $ , и обратный стехиометрические коэффициенты $ \nu_{ki}'' $ и химический символ вещества k $ \chi_k $ могут быть записаны следующим образом:

$ \nu_{ki}' = \begin{pmatrix}1 & 1 & 0 & 0 & 0 \\0 & 0 & 0 & 1 & 1 \\\end{pmatrix}, \nu_{ki}'' = \begin{pmatrix}0 & 0 & 1 & 1 & 0 \\0 & 1 & 1 & 0 & 0 \\\end{pmatrix} , \chi_k = \begin{pmatrix}N_2 \\ O \\ NO \\ N \\ O_2 \\\end{pmatrix} $

Реакции задаются в файле constant/thermophysicalProperties , где находится химический считыватель и файл, содержащий реакции и термодинамические свойства: 



/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v1912                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "constant";
    object      thermophysicalProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

thermoType
{
    type            hePsiThermo;
    mixture         reactingMixture;
    transport       sutherland;
    thermo          janaf;
    energy          sensibleEnthalpy;
    equationOfState perfectGas;
    specie          specie;
}



chemistryReader foamChemistryReader;
foamChemistryFile "<constant>/reactions";
foamChemistryThermoFile "<constant>/thermo";

//CHEMKINFile         "<case>/chemkin/chem.inp";
//CHEMKINThermoFile   "<case>/chemkin/therm.dat";
//CHEMKINTransportFile "<case>/chemkin/transportProperties";

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

В приведённом выше примере реакции выглядят следующим образом: 


elements        2(O N);

species         5(O O2 N2 N NO);

reactions
{
    un-named-reaction-0
    {
        type            reversibleArrheniusReaction;
        reaction        "NO^1 + N^1 = N2^1 + O^1";
        A               2.1e+10;
        beta            0;
        Ta              0;
    }
    un-named-reaction-1
    {
        type            reversibleArrheniusReaction;
        reaction        "N^1 + O2^1 = NO^1 + O^1";
        A               5.83e+06;
        beta            1.01;
        Ta              3120;
    }
}


В нашем примере обе реакции являются обратимыми реакциями Аррениуса. A - это per-exponent фактор, Ta - температура активации и beta - показатель температуры при расчёте скорости реакции (см. дополнительные пояснения ниже). Значения взяты из книги [3]. Заметим, что значения констант в OpenFOAM вычисляются в собственной системе единиц измерения: kmol, m3, s, K.

2.2 Скорость химической реакции

Установив химические реакции, которые желательно должны быть разрешёнными, следующим шагом является описание скорости химической реакции, т.е. изменение концентрации отдельных веществ со временем. Скорость реакции $ \omega_i $ можно записать так (смотрите также книгу [4]):

$ \omega_i = k_{fi} \prod_k [X]_k^{\nu_{ki}'} - k_{ri} \prod_k [X]_k^{\nu_{ki}''} $

$ k_{fi} $ и $ k_{ri} $ - прямая и обратная константы скорости реакции i. Обратите внимание, что в OpenFoam по-умолчанию используются значения показателей - стехиометрические коэффициенты: $ \nu_{ki}', \nu_{ki}'' $. Однако, пользователь может установить их явно, добавив значок ^ после химического символа (смотрите также https://openfoamwiki.net/images/8/8a/Christ_OFW6_slides.pdf). $ [X]_k $ - это молярная концентрация вещества k. Чтобы продолжить наш пример, скорость реакции механизма Зельдовича можно записать так:

$ \omega_1 = k_{f1}[NO][N] - k_{r1}[N2][O] $
$ \omega_2 = k_{f2}[N][O2] - k_{r2}[NO][O] $

Расчёт скорости реакции каждой реакции производится в функции omega файла $FOAM_SRC/thermophysicalModels/chemistryModel/chemistryModel/StandardChemistryModel/StandardChemistryModel.C: 



template<class ReactionThermo, class ThermoType>
Foam::scalar Foam::StandardChemistryModel<ReactionThermo, ThermoType>::omega
(
    const Reaction<ThermoType>& R,
    const scalarField& c,
    const scalar T,
    const scalar p,
    scalar& pf,
    scalar& cf,
    label& lRef,
    scalar& pr,
    scalar& cr,
    label& rRef
) const
{
    const scalar kf = R.kf(p, T, c);
    const scalar kr = R.kr(kf, p, T, c);

    pf = 1.0;
    pr = 1.0;

    const label Nl = R.lhs().size();
    const label Nr = R.rhs().size();

    label slRef = 0;
    lRef = R.lhs()[slRef].index;

    pf = kf;
    for (label s = 1; s < Nl; s++)
    {
        const label si = R.lhs()[s].index;

        if (c[si] < c[lRef])
        {
            const scalar exp = R.lhs()[slRef].exponent;
            pf *= pow(max(c[lRef], 0.0), exp);
            lRef = si;
            slRef = s;
        }
        else
        {
            const scalar exp = R.lhs()[s].exponent;
            pf *= pow(max(c[si], 0.0), exp);
        }
    }
    cf = max(c[lRef], 0.0);

    {
        const scalar exp = R.lhs()[slRef].exponent;
        if (exp < 1.0)
        {
            if (cf > SMALL)
            {
                pf *= pow(cf, exp - 1.0);
            }
            else
            {
                pf = 0.0;
            }
        }
        else
        {
            pf *= pow(cf, exp - 1.0);
        }
    }

    label srRef = 0;
    rRef = R.rhs()[srRef].index;

    // Find the matrix element and element position for the rhs
    pr = kr;
    for (label s = 1; s < Nr; s++)
    {
        const label si = R.rhs()[s].index;
        if (c[si] < c[rRef])
        {
            const scalar exp = R.rhs()[srRef].exponent;
            pr *= pow(max(c[rRef], 0.0), exp);
            rRef = si;
            srRef = s;
        }
        else
        {
            const scalar exp = R.rhs()[s].exponent;
            pr *= pow(max(c[si], 0.0), exp);
        }
    }
    cr = max(c[rRef], 0.0);

    {
        const scalar exp = R.rhs()[srRef].exponent;
        if (exp < 1.0)
        {
            if (cr>SMALL)
            {
                pr *= pow(cr, exp - 1.0);
            }
            else
            {
                pr = 0.0;
            }
        }
        else
        {
            pr *= pow(cr, exp - 1.0);
        }
    }

    return pf*cf - pr*cr;
}


Чтобы получить скорость реакции веществ k, скорости реакции каждой реакции должны быть умножены на стехиометрический коэффициент и суммированы с положительным знаком, если вещества - продукты, и с отрицательным знаком, если вещества - реагенты. (см. также e.g. https://en.wikipedia.org/wiki/Reaction_rate). Это может быть записано так (смотрите также https://www3.nd.edu/~powers/ame.60636/chemkin2000.pdf и другие книги здесь):

$ \omega_k = \sum_i \omega_i ({\nu_{ki}'} - {\nu_{ki}''}) $

Для нашего примера реакции мы получаем:

$ \omega_k = \frac{d [X_k]}{dt} = \begin{pmatrix}\frac{d[N_2]}{dt} \\ \frac{d[O]}{dt} \\ \frac{d[NO]}{dt} \\ \frac{d[N]}{dt} \\ \frac{d[O_2]}{dt} \\\end{pmatrix} = \begin{pmatrix} -\omega_1 \\ -\omega_1 + \omega_2 \\ \omega_1 + \omega_2 \\ \omega_1 - \omega_2 \\ \omega_2 \\\end{pmatrix} = \begin{pmatrix} -(k_{f1}[NO][N] - k_{r1}[N2][O]) \\ -(k_{f1}[NO][N] - k_{r1}[N2][O]) + k_{f2}[N][O2] - k_{r2}[NO][O] \\ k_{f1}[NO][N] - k_{r1}[N2][O] + k_{f2}[N][O2] - k_{r2}[NO][O] \\ k_{f1}[NO][N] - k_{r1}[N2][O] - (k_{f2}[N][O2] - k_{r2}[NO][O]) \\ k_{f2}[N][O2] - k_{r2}[NO][O] \\\end{pmatrix} $

Рассматривая приведённое выше уравнение, мы видим, что уравнение, описывающее изменение концентрации веществ, является связанным нелинейным дифференциальным уравнением. Расчёт скорости реакции каждого вещества производится в функции файла $FOAM_SRC/thermophysicalModels/chemistryModel/chemistryModel/StandardChemistryModel/StandardChemistryModel.C: 



template<class ReactionThermo, class ThermoType>
void Foam::StandardChemistryModel<ReactionThermo, ThermoType>::omega
(
    const scalarField& c,
    const scalar T,
    const scalar p,
    scalarField& dcdt
) const
{
    scalar pf, cf, pr, cr;
    label lRef, rRef;

    dcdt = Zero;

    forAll(reactions_, i)
    {
        const Reaction<ThermoType>& R = reactions_[i];

        scalar omegai = omega
        (
            R, c, T, p, pf, cf, lRef, pr, cr, rRef
        );

        forAll(R.lhs(), s)
        {
            const label si = R.lhs()[s].index;
            const scalar sl = R.lhs()[s].stoichCoeff;
            dcdt[si] -= sl*omegai;
        }

        forAll(R.rhs(), s)
        {
            const label si = R.rhs()[s].index;
            const scalar sr = R.rhs()[s].stoichCoeff;
            dcdt[si] += sr*omegai;
        }
    }
}

2.3 Константы скорости химической реакции

Остальные величины, которые надо определить - это константы скорости прямой и обратной реакции $ k_f $ и $ k_r $.

2.3.1 Константы скорости прямой реакции

Константа скорости прямой реакции вычисляется по уравнению Аррениуса следующим образом (смотрите также [5]):

$ k_f = A T^\beta exp(-T_a/T) $

A - это пер-экспонентный фактор, Ta - температура активации и beta - показатель температуры при расчёте скорости реакции. Код soure может быть найден в $FOAM_SRC/thermophysicalModels/specie/reaction/reactionRate/ArrheniusReactionRate/ArrheniusReactionRateI.H: 



inline Foam::scalar Foam::ArrheniusReactionRate::operator()
(
    const scalar p,
    const scalar T,
    const scalarField&
) const
{
    scalar ak = A_;

    if (mag(beta_) > VSMALL)
    {
        ak *= pow(T, beta_);
    }

    if (mag(Ta_) > VSMALL)
    {
        ak *= exp(-Ta_/T);
    }

    return ak;
}


2.3.2 Константы обратной скорости

Обратная скорость 1-й реакции вычисляется следующим образом (см. также https://www3.nd.edu/~powers/ame.60636/chemkin2000.pdf) :

$ k_{ri} = \frac{k_{fi}}{K_{ci}} $

Константа равновесия вычисляется следующим образом:

$ {K_{ci}} = K_{pi} (\frac{p_{atm}}{R T})^{\sum_k \nu_{ki}} $

$ \sum \nu_{ki} = \sum \nu_{ki}' - \sum \nu_{ki}'' $ является суммой стехиометрического коэффициента 1-й реакции. Кроме того


$ K_{pi} = exp \left( \sum_k \nu_{ki}( \frac{S^0_k}{R} - \frac{H^0_k}{RT}) \right) $

$ H^0_k $ и $ S^0_k $ являются молярной энтальпией и энтропией вида k. .

Константы обратной скорости подсчитываются в файле $FOAM_SRC/thermophysicalModels/specie/reaction/Reactions/ReversibleReaction/ReversibleReaction.C: 


template
<
    template<class> class ReactionType,
    class ReactionThermo,
    class ReactionRate
>
Foam::scalar Foam::ReversibleReaction
<
    ReactionType,
    ReactionThermo,
    ReactionRate
>::kr
(
    const scalar kfwd,
    const scalar p,
    const scalar T,
    const scalarField& c
) const
{
    return kfwd/max(this->Kc(p, T), 1e-6);
}


Расчёт Kc можно найти в файле $FOAM_SRC/thermophysicalModels/specie/thermo/thermo/thermoI.H. 


template<class Thermo, template<class> class Type>
inline Foam::scalar
Foam::species::thermo<Thermo, Type>::K(const scalar p, const scalar T) const
{
    scalar arg = -this->Y()*this->G(Pstd, T)/(RR*T);

    if (arg < 600)
    {
        return exp(arg);
    }
    else
    {
        return VGREAT;
    }
}



template<class Thermo, template<class> class Type>
inline Foam::scalar
Foam::species::thermo<Thermo, Type>::Kp(const scalar p, const scalar T) const
{
    return K(p, T);
}



template<class Thermo, template<class> class Type>
inline Foam::scalar
Foam::species::thermo<Thermo, Type>::Kc(const scalar p, const scalar T) const
{
    const scalar nm = this->Y()/this->W();

    if (equal(nm, SMALL))
    {
        return Kp(p, T);
    }
    else
    {
        return Kp(p, T)*pow(Pstd/(RR*T), nm);
    }
}


Не так легко увидеть, как сумма стехиометрических коэффициентов входит в расчёт константы равновесия. Поэтому это объясняется немного более подробно (см. https://www.cfd-online.com/Forums/openfoam/224970-understanding-reverse-rate-calculation-chemical-reactions.html). Для расчёта константы равновесия для каждой реакции задаётся псевдосмесь, в которой термодинамические свойства уравновешиваются стехиометрическим коэффициентом. Это делается в конструкторе объекта реакции путём вызова функции setThermo: 


template<class ReactionThermo>
Foam::Reaction<ReactionThermo>::Reaction
(
    const speciesTable& species,
    const List<specieCoeffs>& lhs,
    const List<specieCoeffs>& rhs,
    const HashPtrTable<ReactionThermo>& thermoDatabase,
    bool initReactionThermo
)
:
    ReactionThermo::thermoType(*thermoDatabase[species[0]]),
    name_("un-named-reaction-" + Foam::name(getNewReactionID())),
    species_(species),
    lhs_(lhs),
    rhs_(rhs)
{
    if (initReactionThermo)
    {
        setThermo(thermoDatabase);
    }
}

template<class ReactionThermo>
void Foam::Reaction<ReactionThermo>::setThermo
(
    const HashPtrTable<ReactionThermo>& thermoDatabase
)
{

    typename ReactionThermo::thermoType rhsThermo
    (
        rhs_[0].stoichCoeff
        *(*thermoDatabase[species_[rhs_[0].index]]).W()
        *(*thermoDatabase[species_[rhs_[0].index]])
    );

    for (label i=1; i<rhs_.size(); ++i)
    {
        rhsThermo +=
            rhs_[i].stoichCoeff
        *(*thermoDatabase[species_[rhs_[i].index]]).W()
        *(*thermoDatabase[species_[rhs_[i].index]]);
    }

    typename ReactionThermo::thermoType lhsThermo
    (
        lhs_[0].stoichCoeff
       *(*thermoDatabase[species_[lhs_[0].index]]).W()
       *(*thermoDatabase[species_[lhs_[0].index]])
    );

    for (label i=1; i<lhs_.size(); ++i)
    {
        lhsThermo +=
            lhs_[i].stoichCoeff
           *(*thermoDatabase[species_[lhs_[i].index]]).W()
           *(*thermoDatabase[species_[lhs_[i].index]]);
    }

    ReactionThermo::thermoType::operator=(lhsThermo == rhsThermo);
}


Чтобы определить, как вычисляется константа равновесия Kc, мы должны взглянуть на операторы +=, * и == определённые в термодинамических пакетах. В приведённом выше перечне строятся виды, - термодинамические свойства сначала умножаются на стехиометрический коэффициент и молярный вес, а затем в правой и левой части химической реакции выполняется суммирование видов, для последнего шага применяется оператор == . В качестве наглядного примера для объяснения этой операции возьмем постоянные свойства термодинамического пакета. Выводы, сделанные здесь,могут быть применены и ко всем другим пакетам. 


template<class EquationOfState>
inline Foam::hConstThermo<EquationOfState> Foam::operator*
(
    const scalar s,
    const hConstThermo<EquationOfState>& ct
)
{
    return hConstThermo<EquationOfState>
    (
        s*static_cast<const EquationOfState&>(ct),
        ct.Cp_,
        ct.Hf_
    );
}

template<class EquationOfState>
inline void Foam::hConstThermo<EquationOfState>::operator+=
(
    const hConstThermo<EquationOfState>& ct
)
{
    scalar Y1 = this->Y();

    EquationOfState::operator+=(ct);

    if (mag(this->Y()) > SMALL)
    {
        Y1 /= this->Y();
        scalar Y2 = ct.Y()/this->Y();

        Cp_ = Y1*Cp_ + Y2*ct.Cp_;
        Hf_ = Y1*Hf_ + Y2*ct.Hf_;
    }
}

template<class EquationOfState>
inline Foam::hConstThermo<EquationOfState> Foam::operator==
(
    const hConstThermo<EquationOfState>& ct1,
    const hConstThermo<EquationOfState>& ct2
)
{
    EquationOfState eofs
    (
        static_cast<const EquationOfState&>(ct1)
     == static_cast<const EquationOfState&>(ct2)
    );

    return hConstThermo<EquationOfState>
    (
        eofs,
        ct2.Y()/eofs.Y()*ct2.Cp_
      - ct1.Y()/eofs.Y()*ct1.Cp_,
        ct2.Y()/eofs.Y()*ct2.Hf_
      - ct1.Y()/eofs.Y()*ct1.Hf_
    );
}


Первый оператор анализирует * оператор. Когда объект hConstThermo мультиплицируется скаляром s, оператор возвращает новый объект, и только атрибут equationOfState мультиплицируется этим скаляром. По этой причине мы должны посмотреть на оператор * в объекте equationOfState. В качестве примера мы рассмотрим уравнение состояния идеального газа perfectGas: 


template<class Specie>
inline Foam::perfectGas<Specie> Foam::operator*
(
    const scalar s,
    const perfectGas<Specie>& pg
)
{
    return perfectGas<Specie>(s*static_cast<const Specie&>(pg));
}


Объект perfectGas возвращается, если атрибут Specie мультиплицируется с этим скаляром. Для атрибута Specie оператор * operator возвращает вид вместе с Y_ attribute мультиплецированным по скаляру s: 


inline specie operator*(const scalar s, const specie& st)
{
    return specie
    (
        s*st.Y_,
        st.molWeight_
    );
}


Таким образом, применение оператора приводит только к изменению атрибута Y_ attribute как в левой, так и в правой части химической реакции:

$ Y' = \nu' W' $
$ Y'' = \nu'' W'' $

Следующий оператор, который мы рассмотрим - это оператор += . Мы видим, что сначала атрибут Y_ attribute сохраняется в переменной Y1. Затем выполняется оператор += уравнения состояния. После этого берется взвешенная сумма термодинамических свойств двух вовлеченных видов.. 


template<class Specie>
inline void Foam::perfectGas<Specie>::operator+=(const perfectGas<Specie>& pg)
{
    Specie::operator+=(pg);
}





inline void specie::operator+=(const specie& st)
{
    const scalar sumY = Y_ + st.Y_;
    if (mag(sumY) > SMALL)
    {
        molWeight_ = sumY/(Y_/molWeight_ + st.Y_/st.molWeight_);
    }

    Y_ = sumY;
}



Легко заметить, что для атрибутов Y и molWeight для левой и правой частей химического уравнения мы получаем:

$ Y' = \sum_k \nu_k' W_k' $
$ Y'' = \sum_k \nu_k'' W_k'' $
$ molWeight' = \frac{\sum_k \nu_k' W_k'}{\sum_k \nu_k'} $
$ molWeight'' = \frac{\sum_k \nu_k'' W_k''}{\sum_k \nu_k''} $

Для термодинамических свойств Hf (химической энтальпии) и удельной теплоёмкости получаем:

$ Hf' = \frac{\sum_k \nu_k' W_k'Hf'_k}{\sum_k \nu_k'W_k'} $
$ Hf'' = \frac{\sum_k \nu_k'' W_k''Hf''_k}{\sum_k \nu_k''W_k''} $
$ cp' = \frac{\sum_k \nu_k' W_k'cp'_k}{\sum_k \nu_k'W_k'} $
$ cp'' = \frac{\sum_k \nu_k'' W_k''cp''_k}{\sum_k \nu_k''W_k''} $

Последний оператор, который мы рассмотрим - это == operator. Также здесь сначала выполняется уравнение состояния == operator. Также, как и для других операторов, уравнение состояния == operator выполняет == видов: 


template<class Specie>
inline Foam::perfectGas<Specie> Foam::operator==
(
    const perfectGas<Specie>& pg1,
    const perfectGas<Specie>& pg2
)
{
    return perfectGas<Specie>
    (
        static_cast<const Specie&>(pg1) == static_cast<const Specie&>(pg2)
    );
}




inline specie operator==(const specie& st1, const specie& st2)
{
    scalar diffY = st2.Y_ - st1.Y_;
    if (mag(diffY) < SMALL)
    {
        diffY = SMALL;
    }

    const scalar diffRW = st2.Y_/st2.molWeight_ - st1.Y_/st1.molWeight_;

    #ifdef __clang__
    // Using intermediate volatile bool to prevent compiler optimising out the
    // if block (above) - CLANG 3.7.1
    volatile const bool valid = (mag(diffRW) > SMALL);
    const scalar molWeight = valid ? diffY/diffRW : GREAT;
    #else
    scalar molWeight = GREAT;
    if (mag(diffRW) > SMALL)
    {
         molWeight = diffY/diffRW;
    }
    #endif

    return specie(diffY, molWeight);
}



Таким образом, после применения == мы получаем для Y и molWeight атрибут вида:

$ Y = \sum_k \nu_k' W_k' - \sum_k \nu_k'' W_k'' $
$ molWeight' = \frac{\sum_k \nu_k' W_k' - \sum_k \nu_k'' W_k''}{\sum_k \nu_k' - \sum_k \nu_k''} $

Для термодинамических свойств мы имеем:

$ cp = \frac{\sum_k \nu_k' W_k'cp_k' - \sum_k \nu_k'' W_k''cp_k''}{\sum_k \nu_k' W_k' - \sum_k \nu_k'' W_k'' } $
$ Hf = \frac{\sum_k \nu_k' W_k'Hf_k' - \sum_k \nu_k'' W_k''Hf_k''}{\sum_k \nu_k' W_k' - \sum_k \nu_k'' W_k'' } $

В конце концов, если мы посмотрим, как вычисляются энтальпия и энтропия для комбинации hconst и perfectGas, то очень скоро получим то же определение константы равновесия $ K_{ci} $ , что и для chemkin: 


template<class Thermo, template<class> class Type>
inline Foam::scalar
Foam::species::thermo<Thermo, Type>::G(const scalar p, const scalar T) const
{
    return this->Ha(p, T) - T*this->S(p, T);
}




template<class EquationOfState>
inline Foam::scalar Foam::hConstThermo<EquationOfState>::Ha
(
    const scalar p, const scalar T
) const
{
    return Hs(p, T) + Hc();
}



template<class EquationOfState>
inline Foam::scalar Foam::hConstThermo<EquationOfState>::Hs
(
    const scalar p, const scalar T
) const
{
    return Cp_*T + EquationOfState::H(p, T);
}



template<class EquationOfState>
inline Foam::scalar Foam::hConstThermo<EquationOfState>::Hc() const
{
    return Hf_;
}




template<class Specie>
inline Foam::scalar Foam::perfectGas<Specie>::S(scalar p, scalar T) const
{
    return -this->R()*log(p/Pstd);
}




template<class Specie>
inline Foam::scalar Foam::perfectGas<Specie>::H(scalar p, scalar T) const
{
    return 0;
}


2.4 Тепло, выделяющееся в результате химических реакций

Тепло, выделяющееся в результате химических реакций, возвращается функцией Qdot() в файле $FOAM_SRC/thermophysicalModels/chemistryModel/chemistryModel/StandardChemistryModel/StandardChemistryModel.C (если выбрана стандартная химическая модель):


template<class ReactionThermo, class ThermoType> Foam::tmp<Foam::volScalarField> Foam::StandardChemistryModel<ReactionThermo, ThermoType>::Qdot() const { 

   tmp<volScalarField> tQdot
   (
       new volScalarField
       (
           IOobject
           (
               "Qdot",
               this->mesh_.time().timeName(),
               this->mesh_,
               IOobject::NO_READ,
               IOobject::NO_WRITE,
               false
           ),
           this->mesh_,
           dimensionedScalar(dimEnergy/dimVolume/dimTime, Zero)
       )
   );

   if (this->chemistry_)
   {
       scalarField& Qdot = tQdot.ref();

       forAll(Y_, i)
       {
           forAll(Qdot, celli)
           {
               const scalar hi = specieThermo_[i].Hc();
               Qdot[celli] -= hi*RR_[i][celli];
           }
       }
   }

   return tQdot;

} </cpp>

Hc является химической энтальпией и зависит от выбранной термодинамической модели.

3. Эволюция видов


{
    forAll(Y, specieI)
    {
        volScalarField& Yi = Y[specieI];

        solve
        (
            fvm::ddt(rho, Yi) - chemistry.RR(specieI),
            mesh.solver("Yi")
        );
    }
}


4. Эволюция энергии

Уравнение энергии выведено из первого закона термодинамики (смотрите также https://en.wikipedia.org/wiki/Enthalpy): Для замкнутых систем изменение внутренней энергии равно количеству тепла, подаваемого в систему. Для открытых систем изменение энтальпии равно теплу, подаваемому в систему. Для первого случая мы имеем:

$ _ u = u_0 + q = u_0 + \int \dot q dt $
(x)

Если мы добавим с обеих сторон давление, делённое на плотность, мы получим:


$ _ u + \frac{p}{\rho} = u_0 + \frac{p}{\rho} + \int \dot q dt \ \ , \ \ h = u + \frac{p}{\rho} $
(x)

$ u, u_0 $ это энергия и внутренняя энергия в момент времени 0. 


{
    volScalarField& h = thermo.he();

    if (constProp == "volume")
    {
        h[0] = u0 + p[0]/rho[0] + integratedHeat;
    }
    else
    {
        h[0] = h0 + integratedHeat;
    }

    thermo.correct();
}


С помощью приведённого выше уравнения можно подсчитать эволюцию энтальпии во времени. Зная функциональную зависимость теплоёмкости при постоянном давлении $ c_p $ , можно установить связь между энтальпией и температурой. Пример того, как это происходит в OpenFOAM с полиномиальной зависимостью теплоемкости от температуры, можно найти здесь https://caefn.com/openfoam/temperature-calculation.

5. Эволюция давления

Для расчёта давления используется модель Дальтона: В ней говорится, что термодинамическое давление может быть вычислено как сумма парциальных давлений отдельных видов (см. https://en.wikipedia.org/wiki/Dalton%27s_law):


$ _ 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_{tot}}{V}}_{\rho_0} R T $
(x)

$ p_i, n_i, m_i, M_i, y_i, m_{tot}, \rho_0 $ это парциальное давление, число молекул, масса, молярная масса, массовая доля вида i, общая масса и начальная плотность соответственно. 


{
    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;
    }
}


Обратите внимание, что константа 1000 в приведенном выше коде исходит из того факта, что в OpenFOAM молярный вес указан в [kg/kmol]. Смотрите файл $FOAM_SRC/thermophysicalModels/specie/specie/specie.H.

Ссылки

  1. ↑ Powers, Joseph. Combustion thermodynamics and dynamics. Cambridge University Press, 2016.
  2. ↑ Powers, Joseph. Combustion thermodynamics and dynamics. Cambridge University Press, 2016.
  3. ↑ Powers, Joseph. Combustion thermodynamics and dynamics. Cambridge University Press, 2016.
  4. ↑ Poinsot, Thierry, and Denis Veynante. Theoretical and numerical combustion. RT Edwards, Inc., 2005.
  5. ↑ Poinsot, Thierry, and Denis Veynante. Theoretical and numerical combustion. RT Edwards, Inc., 2005.
Ссылка на оригинал - https://openfoamwiki.net/index.php/ChemFoam
Источник — «http://ru.wiki.laduga.ru/index.php?title=ChemFoam&oldid=417»