ENH: laminarFlameSpeed: added RaviPetersen flame speed for Hydrogen

src/thermophysicalModels/laminarFlameSpeed/Make/files

@@ -3,5 +3,6 @@ laminarFlameSpeed/laminarFlameSpeedNew.C
 constant/constant.C
 Gulders/Gulders.C
 GuldersEGR/GuldersEGR.C
+RaviPetersen/RaviPetersen.C
 
 LIB = $(FOAM_LIBBIN)/liblaminarFlameSpeedModels

src/thermophysicalModels/laminarFlameSpeed/RaviPetersen/RaviPetersen.C

@@ -0,0 +1,365 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "RaviPetersen.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace laminarFlameSpeedModels
+{
+    defineTypeNameAndDebug(RaviPetersen, 0);
+
+    addToRunTimeSelectionTable
+    (
+        laminarFlameSpeed,
+        RaviPetersen,
+        dictionary
+    );
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::laminarFlameSpeedModels::RaviPetersen::RaviPetersen
+(
+    const dictionary& dict,
+    const psiuReactionThermo& ct
+)
+:
+    laminarFlameSpeed(dict, ct),
+    coeffsDict_(dict.subDict(typeName + "Coeffs").subDict(fuel_)),
+    pPoints_(coeffsDict_.lookup("pPoints")),
+    EqRPoints_(coeffsDict_.lookup("EqRPoints")),
+    alpha_(coeffsDict_.lookup("alpha")),
+    beta_(coeffsDict_.lookup("beta")),
+    TRef_(readScalar(coeffsDict_.lookup("TRef")))
+{
+    checkPointsMonotonicity("equivalenceRatio", EqRPoints_);
+    checkPointsMonotonicity("pressure", pPoints_);
+    checkCoefficientArrayShape("alpha", alpha_);
+    checkCoefficientArrayShape("beta", beta_);
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::laminarFlameSpeedModels::RaviPetersen::~RaviPetersen()
+{}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::laminarFlameSpeedModels::RaviPetersen::checkPointsMonotonicity
+(
+    const word& name,
+    const List<scalar>& x
+) const
+{
+    for (label i = 1; i < x.size(); i ++)
+    {
+        if (x[i] <= x[i-1])
+        {
+            FatalIOErrorIn
+            (
+                "laminarFlameSpeedModels::RaviPetersen::checkPointsMonotonicity"
+                "("
+                    "const word&, "
+                    "const List<scalar>&"
+                ") const",
+                coeffsDict_
+            )   << "Data points for the " << name
+                << " do not increase monotonically" << endl
+                << exit(FatalIOError);
+        }
+    }
+}
+
+
+void Foam::laminarFlameSpeedModels::RaviPetersen::checkCoefficientArrayShape
+(
+    const word& name,
+    const List<List<List<scalar> > >& x
+) const
+{
+    bool ok = true;
+
+    ok &= x.size() == EqRPoints_.size() - 1;
+
+    forAll(x, i)
+    {
+        ok &= x[i].size() == pPoints_.size();
+
+        forAll(x[i], j)
+        {
+            ok &= x[i][j].size() == x[i][0].size();
+        }
+    }
+
+    if (!ok)
+    {
+        FatalIOErrorIn
+        (
+            "laminarFlameSpeedModels::RaviPetersen::checkCoefficientArrayShape"
+            "("
+                "const word&, "
+                "const List<List<List<scalar> > >&"
+            ") const",
+            coeffsDict_
+        )   << "Inconsistent size of " << name << " coefficients array" << endl
+            << exit(FatalIOError);
+    }
+}
+
+
+inline bool Foam::laminarFlameSpeedModels::RaviPetersen::interval
+(
+    const List<scalar>& xPoints,
+    const scalar x,
+    label& xIndex,
+    scalar& xXi,
+    scalar& xLim
+) const
+{
+    if (x < xPoints.first())
+    {
+        xIndex = 0;
+        xXi = 0.0;
+        xLim = xPoints.first();
+        return false;
+    }
+
+    else if (x > xPoints.last())
+    {
+        xIndex = xPoints.size() - 2;
+        xXi = 1.0;
+        xLim = xPoints.last();
+        return false;
+    }
+
+    for (xIndex = 0; x > xPoints[xIndex+1]; xIndex ++)
+    {
+        // increment xIndex until xPoints[xIndex] < x < xPoints[xIndex+1]
+    }
+
+    xXi = (x - xPoints[xIndex])/(xPoints[xIndex+1] - xPoints[xIndex]);
+    xLim = x;
+
+    return true;
+}
+
+
+inline Foam::scalar Foam::laminarFlameSpeedModels::RaviPetersen::polynomial
+(
+    const List<scalar>& coeffs,
+    const scalar x
+) const
+{
+    scalar xPow = 1.0;
+    scalar y = 0.0;
+    forAll(coeffs, i)
+    {
+        y += coeffs[i]*xPow;
+        xPow *= x;
+    }
+    return y;
+}
+
+
+inline Foam::scalar Foam::laminarFlameSpeedModels::RaviPetersen::dPolynomial
+(
+    const List<scalar>& coeffs,
+    const scalar x
+) const
+{
+    scalar xPow = 1.0;
+    scalar y = 0.0;
+    for (label i = 1; i < coeffs.size(); i++)
+    {
+        y += i*coeffs[i]*xPow;
+        xPow *= x;
+    }
+    return y;
+}
+
+
+inline Foam::scalar Foam::laminarFlameSpeedModels::RaviPetersen::THatPowB
+(
+    const label EqRIndex,
+    const label pIndex,
+    const scalar EqR,
+    const scalar Tu
+) const
+{
+    return pow
+    (
+        Tu/TRef_,                                                     
+        polynomial(beta_[EqRIndex][pIndex],EqR)                       
+    );
+}
+
+
+inline Foam::scalar
+Foam::laminarFlameSpeedModels::RaviPetersen::correlationInRange
+(
+    const label EqRIndex,
+    const label pIndex,
+    const scalar EqR,
+    const scalar Tu
+) const
+{
+    // standard correlation
+    return
+        polynomial(alpha_[EqRIndex][pIndex],EqR)
+       *THatPowB(EqRIndex, pIndex, EqR, Tu);
+}
+ 
+
+inline Foam::scalar
+Foam::laminarFlameSpeedModels::RaviPetersen::correlationOutOfRange
+(
+    const label EqRIndex,
+    const label pIndex,
+    const scalar EqR,
+    const scalar EqRLim,
+    const scalar Tu
+) const
+{
+    scalar A = polynomial(alpha_[EqRIndex][pIndex], EqRLim);
+    scalar dA = dPolynomial(alpha_[EqRIndex][pIndex], EqRLim);
+    scalar dB = dPolynomial(beta_[EqRIndex][pIndex], EqRLim);
+    scalar TB = THatPowB(EqRIndex, pIndex, EqRLim, Tu);
+
+    // linear extrapolation from the bounds of the correlation
+    return max(TB*(A + (dA + A*log(Tu/TRef_)*dB)*(EqR - EqRLim)), 0.0);
+}
+
+
+inline Foam::scalar Foam::laminarFlameSpeedModels::RaviPetersen::speed
+(
+    const scalar EqR,
+    const scalar p,
+    const scalar Tu
+) const
+{
+    scalar Su = 0, s;
+
+    label EqRIndex, pIndex;
+    scalar EqRXi, pXi;
+    scalar EqRLim, pLim;
+    bool EqRInRange;
+
+    EqRInRange = interval(EqRPoints_, EqR, EqRIndex, EqRXi, EqRLim);
+
+    interval(pPoints_, p, pIndex, pXi, pLim);
+
+    for (label pI = 0; pI < 2; pI ++)
+    {
+        if (EqRInRange)
+        {
+            s = correlationInRange(EqRIndex, pIndex + pI, EqR, Tu);
+        }
+        else
+        {
+            s = correlationOutOfRange(EqRIndex, pIndex + pI, EqR, EqRLim, Tu);
+        }
+
+        Su += (1 - pXi)*s;
+        pXi = 1 - pXi;
+    }
+
+    return Su;
+}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::volScalarField>
+Foam::laminarFlameSpeedModels::RaviPetersen::operator()() const
+{
+    const volScalarField& p = psiuReactionThermo_.p();
+    const volScalarField& Tu = psiuReactionThermo_.Tu();
+
+    volScalarField EqR
+    (
+        IOobject
+        (
+            "EqR",
+            p.time().timeName(),
+            p.db(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        p.mesh(),
+        dimensionedScalar("EqR", dimless, 0.0)
+    );
+
+    if (psiuReactionThermo_.composition().contains("ft"))
+    {
+        const volScalarField& ft = psiuReactionThermo_.composition().Y("ft");
+
+        EqR =
+            dimensionedScalar
+            (
+                psiuReactionThermo_.lookup("stoichiometricAirFuelMassRatio")
+            )*ft/max(1 - ft, SMALL);
+    }
+    else
+    {
+        EqR = equivalenceRatio_;
+    }
+
+    tmp<volScalarField> tSu0
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "Su0",
+                p.time().timeName(),
+                p.db(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            p.mesh(),
+            dimensionedScalar("Su0", dimVelocity, 0.0)
+        )
+    );
+
+    volScalarField& Su0 = tSu0();
+
+    forAll(Su0, cellI)
+    {
+        Su0[cellI] = speed(EqR[cellI], p[cellI], Tu[cellI]);
+    }
+
+    return tSu0;
+}
+
+
+// ************************************************************************* //

src/thermophysicalModels/laminarFlameSpeed/RaviPetersen/RaviPetersen.H

@@ -0,0 +1,203 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::laminarFlameSpeedModels::RaviPetersen
+
+Description
+    Laminar flame speed obtained from Ravi and Petersen's correlation.
+
+    The correlation for the laminar flame speed \f$Su\f$ is of the following
+    form:
+    \f[
+        Su = \left( \sum \alpha_i \phi^i \right) 
+        \left( \frac{T}{T_{ref}} \right)^{\left( \sum \beta_j \phi^j \right)} 
+    \f]
+
+    Where \f$\phi\f$ is the equivalence ratio, and \f$\alpha\f$ and \f$\beta\f$
+    are polynomial coefficients given for a number of pressure and equivalence
+    ratio points.
+
+SourceFiles
+    RaviPetersen.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef RaviPetersen_H
+#define RaviPetersen_H
+
+#include "laminarFlameSpeed.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace laminarFlameSpeedModels
+{
+
+/*---------------------------------------------------------------------------*\
+                        Class RaviPetersen Declaration
+\*---------------------------------------------------------------------------*/
+
+class RaviPetersen
+:
+    public laminarFlameSpeed
+{
+    // Private Data
+
+        dictionary coeffsDict_;
+
+        //- Correlation pressure values
+        List<scalar> pPoints_;
+
+        //- Correlation equivalence ratios
+        List<scalar> EqRPoints_;
+
+        //- Correlation alpha coefficients
+        List<List<List<scalar> > > alpha_;
+
+        //- Correlation beta coefficients
+        List<List<List<scalar> > > beta_;
+
+        //- Reference temperature
+        scalar TRef_;
+
+
+    // Private Member Functions
+
+        //- Check that input points are ordered
+        void checkPointsMonotonicity
+        (
+            const word& name,
+            const List<scalar>& x
+        ) const;
+
+        //- Check that the coefficient arrays are of the correct shape
+        void checkCoefficientArrayShape
+        (
+            const word& name,
+            const List<List<List<scalar> > >& x
+        ) const;
+
+        //- Find and interpolate a value in the data point arrays
+        inline bool interval
+        (
+            const List<scalar>& xPoints,
+            const scalar x,
+            label& xIndex,
+            scalar& xXi,
+            scalar& xLim
+        ) const;
+
+        //- Evaluate a polynomial
+        inline scalar polynomial
+        (
+            const List<scalar>& coeffs,
+            const scalar x
+        ) const;
+
+        //- Evaluate a polynomial differential
+        inline scalar dPolynomial
+        (
+            const List<scalar>& coeffs,
+            const scalar x
+        ) const;
+
+        //- Calculate normalised temperature to the power of the B polynomial
+        inline scalar THatPowB
+        (
+            const label EqRIndex,
+            const label pIndex,
+            const scalar EqR,
+            const scalar Tu
+        ) const;
+
+        //- Return the flame speed within the correlation range
+        inline scalar correlationInRange
+        (
+            const label EqRIndex,
+            const label pIndex,
+            const scalar EqR,
+            const scalar Tu
+        ) const;
+
+        //- Extrapolate the flame speed correlation outside its range
+        inline scalar correlationOutOfRange
+        (
+            const label EqRIndex,
+            const label pIndex,
+            const scalar EqR,
+            const scalar EqRLim,
+            const scalar Tu
+        ) const;
+
+        //- Return the laminar flame speed [m/s]
+        inline scalar speed
+        (
+            const scalar EqR,
+            const scalar p,
+            const scalar Tu
+        ) const;
+
+
+        //- Construct as copy (not implemented)
+        RaviPetersen(const RaviPetersen&);
+        void operator=(const RaviPetersen&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("RaviPetersen");
+
+    // Constructors
+
+        //- Construct from dictionary and psiuReactionThermo
+        RaviPetersen
+        (
+            const dictionary&,
+            const psiuReactionThermo&
+        );
+
+
+    //- Destructor
+    virtual ~RaviPetersen();
+
+
+    // Member functions
+
+        //- Return the laminar flame speed [m/s]
+        tmp<volScalarField> operator()() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End laminarFlameSpeedModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

tutorials/combustion/XiFoam/ras/Allclean

@@ -5,7 +5,7 @@ cd ${0%/*} || exit 1    # run from this directory
 . $WM_PROJECT_DIR/bin/tools/CleanFunctions
 
 keepCases="moriyoshiHomogeneous"
-loseCases="moriyoshiHomogeneousPart2"
+loseCases="moriyoshiHomogeneousPart2 moriyoshiHomogeneousHydrogen"
 
 for caseName in $keepCases
 do

tutorials/combustion/XiFoam/ras/Allrun

@@ -6,21 +6,25 @@ cd ${0%/*} || exit 1    # run from this directory
 
 setControlDict()
 {
-    controlDict="system/controlDict"
-
     sed \
-    -e s/"\(deltaT[ \t]*\) 5e-06;"/"\1 1e-05;"/g \
+    -e "s/\(deltaT[ \t]*\) 5e-06;/\1 1e-05;/g" \
-    -e s/"\(endTime[ \t]*\) 0.005;"/"\1 0.015;"/g \
+    -e "s/\(endTime[ \t]*\) 0.005;/\1 0.015;/g" \
-    -e s/"\(writeInterval[ \t]*\) 10;"/"\1 50;"/g \
+    -e "s/\(writeInterval[ \t]*\) 10;/\1 50;/g" \
-    $controlDict > temp.$$
+    -i system/controlDict
-    mv temp.$$ $controlDict
 }
 
+setCombustionProperties()
+{
+    sed \
+    -e "s/\(laminarFlameSpeedCorrelation[ \t]*\) Gulders;/\1 RaviPetersen;/g" \
+    -e "s/\(fuel[ \t]*\) Propane;/\1 HydrogenInAir;/g" \
+    -i constant/combustionProperties
+}
 
 # Do moriyoshiHomogeneous
 ( cd moriyoshiHomogeneous && foamRunTutorials )
 
-# Clone case
+# Clone case for second phase
 cloneCase moriyoshiHomogeneous moriyoshiHomogeneousPart2
 
 # Modify and execute
@@ -32,4 +36,19 @@ cloneCase moriyoshiHomogeneous moriyoshiHomogeneousPart2
     runApplication `getApplication`
 )
 
+# Clone case for hydrogen
+cloneCase moriyoshiHomogeneous moriyoshiHomogeneousHydrogen
+
+# Modify and execute
+(
+    cd moriyoshiHomogeneousHydrogen || exit
+
+    setCombustionProperties
+    mv constant/thermophysicalProperties \
+       constant/thermophysicalProperties.propane
+    mv constant/thermophysicalProperties.hydrogen \
+       constant/thermophysicalProperties
+    runApplication `getApplication`
+)
+
 # ----------------------------------------------------------------- end-of-file

tutorials/combustion/XiFoam/ras/moriyoshiHomogeneous/constant/combustionProperties

@@ -68,6 +68,36 @@ GuldersCoeffs
     }
 }
 
+RaviPetersenCoeffs
+{
+    HydrogenInAir
+    {
+        TRef            320;
+        pPoints         ( 1.0e05 5.0e05 1.0e06 2.0e06 3.0e06 );
+        EqRPoints       ( 0.5 2.0 5.0 );
+        alpha           ( ( (-0.03   -2.347  9.984  -6.734  1.361)
+                            ( 1.61   -9.708 19.026 -11.117  2.098)
+                            ( 2.329 -12.287 21.317 -11.973  2.207)
+                            ( 2.593 -12.813 20.815 -11.471  2.095)
+                            ( 2.728 -13.164 20.794 -11.418  2.086) )
+                          ( ( 3.558   0.162 -0.247   0.0253 0    )
+                            ( 4.818  -0.872 -0.053   0.0138 0    )
+                            ( 3.789  -0.312 -0.208   0.028  0    )
+                            ( 4.925  -1.841  0.211  -0.0059 0    )
+                            ( 4.505  -1.906  0.259  -0.0105 0    ) ) );
+        beta            ( ( ( 5.07  -6.42  3.87  -0.767)
+                            ( 5.52  -6.73  3.88  -0.728)
+                            ( 5.76  -6.92  3.92  -0.715)
+                            ( 6.02  -7.44  4.37  -0.825)
+                            ( 7.84 -11.55  7.14  -1.399) )
+                          ( ( 1.405  0.053 0.022  0    )
+                            ( 1.091  0.317 0      0    )
+                            ( 1.64  -0.03  0.07   0    )
+                            ( 0.84   0.56  0      0    )
+                            ( 0.81   0.64  0      0    ) ) );
+    }
+}
+
 ignite          yes;
 
 ignitionSites

tutorials/combustion/XiFoam/ras/moriyoshiHomogeneous/constant/thermophysicalProperties.hydrogen

@@ -0,0 +1,76 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      thermophysicalProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+thermoType
+{
+    type            heheuPsiThermo;
+    mixture         homogeneousMixture;
+    transport       sutherland;
+    thermo          janaf;
+    equationOfState perfectGas;
+    specie          specie;
+    energy          absoluteEnthalpy;
+}
+
+stoichiometricAirFuelMassRatio stoichiometricAirFuelMassRatio [ 0 0 0 0 0 0 0 ] 34.074;
+
+reactants
+{
+    specie
+    {
+        nMoles          24.8095;
+        molWeight       16.0243;
+    }
+    thermodynamics
+    {
+        Tlow            200;
+        Thigh           5000;
+        Tcommon         1000;
+        highCpCoeffs    ( 3.02082 0.00104314 -2.88613e-07 4.20369e-11 -2.37182e-15 -902.964 2.3064 );
+        lowCpCoeffs     ( 2.99138 0.00343493 -8.43792e-06 9.57755e-09 -3.75097e-12 -987.16 1.95123 );
+    }
+    transport
+    {
+        As              1.67212e-06;
+        Ts              170.672;
+    }
+}
+
+products
+{
+    specie
+    {
+        nMoles          1;
+        molWeight       17.9973;
+    }
+    thermodynamics
+    {
+        Tlow            200;
+        Thigh           5000;
+        Tcommon         1000;
+        highCpCoeffs    ( 2.879 0.00161934 -4.61257e-07 6.41382e-11 -3.3855e-15 -8023.54 4.11691 );
+        lowCpCoeffs     ( 3.3506 0.00176018 -4.28718e-06 5.63372e-09 -2.35948e-12 -8211.42 1.36387 );
+    }
+    transport
+    {
+        As              1.67212e-06;
+        Ts              170.672;
+    }
+}
+
+
+// ************************************************************************* //