ENH: New outletMachNumberPressure BC. It sets pressure at outlet

keeping chocked conditions of Mach number.
This BC can work in two modes, chocked or non-chocked. In the
chocked mode the Ma is an input. In the non-chocked mode
the Ma is calculated from model inputs.

src/TurbulenceModels/compressible/Make/files

@@ -13,6 +13,7 @@ $(BCs)/convectiveHeatTransfer/convectiveHeatTransferFvPatchScalarField.C
 $(BCs)/fixedIncidentRadiation/fixedIncidentRadiationFvPatchScalarField.C
 $(BCs)/lumpedMassWallTemperature/lumpedMassWallTemperatureFvPatchScalarField.C
 $(BCs)/outletMappedUniformInletHeatAddition/outletMappedUniformInletHeatAdditionFvPatchField.C
+$(BCs)/outletMachNumberPressure/outletMachNumberPressureFvPatchScalarField.C
 
 turbulentFluidThermoModels/derivedFvPatchFields/wallFunctions/alphatWallFunctions/alphatWallFunction/alphatWallFunctionFvPatchScalarField.C
 turbulentFluidThermoModels/derivedFvPatchFields/wallFunctions/alphatWallFunctions/alphatJayatillekeWallFunction/alphatJayatillekeWallFunctionFvPatchScalarField.C

src/TurbulenceModels/compressible/turbulentFluidThermoModels/derivedFvPatchFields/outletMachNumberPressure/outletMachNumberPressureFvPatchScalarField.C

+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018 OpenCFD Ltd.
+     \\/     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 "outletMachNumberPressureFvPatchScalarField.H"
+#include "addToRunTimeSelectionTable.H"
+#include "fvPatchFieldMapper.H"
+#include "volFields.H"
+#include "surfaceFields.H"
+#include "fluidThermo.H"
+#include "constants.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::outletMachNumberPressureFvPatchScalarField::
+outletMachNumberPressureFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    fixedValueFvPatchScalarField(p, iF),
+    M_(1),
+    pBack_(0.0),
+    c1_(0.0),
+    A1_(0.0),
+    phiName_("phi"),
+    rhoName_("rho"),
+    UName_("U"),
+    choked_(false),
+    relax_(0.0)
+{}
+
+
+Foam::outletMachNumberPressureFvPatchScalarField::
+outletMachNumberPressureFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    fixedValueFvPatchScalarField(p, iF, dict),
+    M_(dict.lookupOrDefault<scalar>("M", 0.0)),
+    pBack_(readScalar(dict.lookup("pBack"))),
+    c1_(dict.lookupOrDefault<scalar>("c1", 0.0)),
+    A1_(dict.lookupOrDefault<scalar>("A1", 0.0)),
+    phiName_(dict.lookupOrDefault<word>("phi", "phi")),
+    rhoName_(dict.lookupOrDefault<word>("rho", "rho")),
+    UName_(dict.lookupOrDefault<word>("U", "U")),
+    choked_(dict.lookup("choked")),
+    relax_(dict.lookupOrDefault<scalar>("relax", 0))
+{
+    Info<< "The outletMachNumber pressure BC is using the flag choked: "
+        <<  choked_  << endl;
+}
+
+
+Foam::outletMachNumberPressureFvPatchScalarField::
+outletMachNumberPressureFvPatchScalarField
+(
+    const outletMachNumberPressureFvPatchScalarField& ptf,
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    fixedValueFvPatchScalarField(ptf, p, iF, mapper),
+    M_(ptf.M_),
+    pBack_(ptf.pBack_),
+    c1_(ptf.c1_),
+    A1_(ptf.A1_),
+    phiName_(ptf.phiName_),
+    rhoName_(ptf.rhoName_),
+    UName_(ptf.UName_),
+    choked_(ptf.choked_),
+    relax_(ptf.relax_)
+{}
+
+
+Foam::outletMachNumberPressureFvPatchScalarField::
+outletMachNumberPressureFvPatchScalarField
+(
+    const outletMachNumberPressureFvPatchScalarField& tppsf
+)
+:
+    fixedValueFvPatchScalarField(tppsf),
+    M_(tppsf.M_),
+    pBack_(tppsf.pBack_),
+    c1_(tppsf.c1_),
+    A1_(tppsf.A1_),
+    phiName_(tppsf.phiName_),
+    rhoName_(tppsf.rhoName_),
+    UName_(tppsf.UName_),
+    choked_(tppsf.choked_),
+    relax_(tppsf.relax_)
+{}
+
+
+Foam::outletMachNumberPressureFvPatchScalarField::
+outletMachNumberPressureFvPatchScalarField
+(
+    const outletMachNumberPressureFvPatchScalarField& tppsf,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    fixedValueFvPatchScalarField(tppsf, iF),
+    M_(tppsf.M_),
+    pBack_(tppsf.pBack_),
+    c1_(tppsf.c1_),
+    A1_(tppsf.A1_),
+    phiName_(tppsf.phiName_),
+    rhoName_(tppsf.rhoName_),
+    UName_(tppsf.UName_),
+    choked_(tppsf.choked_),
+    relax_(tppsf.relax_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::outletMachNumberPressureFvPatchScalarField::updateCoeffs()
+{
+    if (updated())
+    {
+        return;
+    }
+
+     const volScalarField& p =
+        db().lookupObject<volScalarField>
+        (
+            this->internalField().name()
+        );
+
+    const label patchi = patch().index();
+
+    const scalarField pb(p.oldTime().boundaryField()[patchi]);
+
+    const fvsPatchField<scalar>& phi =
+        patch().lookupPatchField<surfaceScalarField, scalar>(phiName_);
+
+    // Calculate the current mass flow rate
+    if (phi.internalField().dimensions() != dimDensity*dimVelocity*dimArea)
+    {
+        FatalErrorInFunction
+            <<"phi is not a mass flux." << exit(FatalError);
+    }
+
+    const fluidThermo* thermoPtr =
+        db().lookupObjectPtr<fluidThermo>(basicThermo::dictName);
+
+    const volVectorField& U = db().lookupObject<volVectorField>(UName_);
+
+    vectorField Ub(U.boundaryField()[patchi]);
+    const vectorField UbOld(U.oldTime().boundaryField()[patchi]);
+
+    // relax U
+    Ub = relax_*UbOld + (1 -relax_)*Ub;
+
+    const scalarField gamma(thermoPtr->gamma()().boundaryField()[patchi]);
+
+    const fvPatchField<scalar>& rho =
+        patch().lookupPatchField<volScalarField, scalar>(rhoName_);
+
+    const scalarField Mb(mag(Ub)/sqrt(gamma*pb/rho));
+
+    const scalarField ptot
+    (
+        pb*(pow(1 + (gamma-1)/2*sqr(gAverage(Mb)), gamma/(gamma-1)))
+    );
+
+    scalarField M(patch().size(), 1.0);
+
+    if (choked_)
+    {
+        if (M_ > 0.0)
+        {
+            M = M_;
+        }
+        else
+        {
+            FatalErrorInFunction <<" Mach number is lower than zero" << endl
+                << "Pelase specify M in the dictionary"
+                << exit(FatalError);
+        }
+    }
+    else
+    {
+
+        if (A1_ == 0.0 || c1_ == 0.0)
+        {
+             FatalErrorInFunction <<" Please enter values for A1 and c1" << endl
+                << exit(FatalError);
+        }
+
+        const scalarField r = pBack_/ptot;
+        const scalar area = gSum(mag(patch().Sf()));
+        M =
+            A1_/(c1_*area)
+           *sqrt(2/(gamma-1)*(pow(r, 2/gamma) - pow(r, (gamma+1)/gamma)));
+
+        forAll(M, i)
+        {
+            if (M[i] < 0 || r[i] >=1)
+            {
+                WarningInFunction <<" Mach number is lower than zero" << endl
+                    << "or pBack/ptot ratio is larger then one"
+                    << endl;
+            }
+        }
+    }
+
+    scalarField pbNew
+    (
+        ptot/(pow(1 + (gamma-1)/2*sqr(gAverage(M)), gamma/(gamma-1)))
+    );
+
+    // relax pressure
+    pbNew = relax_*pb + (1 -relax_)*pbNew;
+
+    operator==(pbNew);
+
+    fixedValueFvPatchScalarField::updateCoeffs();
+}
+
+
+void Foam::outletMachNumberPressureFvPatchScalarField::write(Ostream& os) const
+{
+    fvPatchScalarField::write(os);
+    os.writeEntry("pBack", pBack_);
+    os.writeEntryIfDifferent("c1", 0.0, c1_);
+    os.writeEntryIfDifferent("A1", 0.0, A1_);
+    os.writeEntry("choked", choked_);
+    os.writeEntryIfDifferent("relax", 0.0, relax_);
+
+    os.writeEntryIfDifferent<word>("phi", "phi", phiName_);
+    os.writeEntryIfDifferent<word>("rho", "rho", rhoName_);
+    os.writeEntryIfDifferent<word>("U", "U", UName_);
+    os.writeEntryIfDifferent<scalar>("M", 0.0, M_);
+
+    writeEntry("value", os);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    makePatchTypeField
+    (
+        fvPatchScalarField,
+        outletMachNumberPressureFvPatchScalarField
+    );
+}
+
+// ************************************************************************* //
+

src/TurbulenceModels/compressible/turbulentFluidThermoModels/derivedFvPatchFields/outletMachNumberPressure/outletMachNumberPressureFvPatchScalarField.H

+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018 OpenCFD Ltd.
+     \\/     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::outletMachNumberPressureFvPatchScalarField
+
+Group
+    grpOutletBoundaryConditions
+
+Description
+    This BC set the outlet pressure for compressible flows is such as a
+    choked condition is achieved given an input Mach number.
+    Optionally, if the BC is operated in a non-choked condition (choked false),
+    the constants C1 and A1 are used to calculate the corresponding Mach number at
+    the outlet.
+
+    The static pressure is calculates as :
+
+        p_s = p_tot / [1+(k-1)/2*M^2]^(k/(k-1))
+
+    where
+
+    p_tot: [Pa] mass averaged total pressure on outlet patch
+    k: [-] mass averaged heat capacity ratio on outlet patch
+    M: [-] target Mach number on the outlet patch defining as either choked or non-choked
+
+    For choked conditions:
+
+        M = constant defined in dict
+
+    For non-choked conditions:
+
+        M = A1/(C1*A_outlet)*sqrt(2/(k-1)*[r^(2/k)-r^((k+1)/k)], r=p_back/p_tot
+
+    where
+
+        A_outlet: [m2] area of outlet patch
+        A1: [m2] constant defined in dict
+        C1: [-] constant defined in dict
+        p_back: [Pa] constant defined in dict
+
+
+Usage
+    \table
+        Property        | Description                 | Required | Default value
+        pBack           | Back pressure               | yes
+        M               | Average desired mach number | no
+        C1              | Model input                 | no       | 0.0
+        A1              | Model input                 | no       | 0.0
+        choked          | The outlet is considered chocked | yes
+        relax           | relaxation factor (1 fully relax) | no | 0.0
+    \endtable
+
+    Example of the boundary condition specification:
+    \verbatim
+    <patchName>
+    {
+        type            outletMachNumberPressure;
+        pBack           101325;
+        c1              1;
+        A1              0.008;
+        choked          false;
+        value           uniform 101325;
+    }
+    \endverbatim
+
+SourceFiles
+    outletMachNumberPressureFvPatchScalarField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef outletMachNumberPressureFvPatchScalarField_H
+#define outletMachNumberPressureFvPatchScalarField_H
+
+#include "fixedValueFvPatchFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                      Class plenumPressureFvPatch Declaration
+\*---------------------------------------------------------------------------*/
+
+class outletMachNumberPressureFvPatchScalarField
+:
+    public fixedValueFvPatchScalarField
+{
+    // Private data
+
+        //- Mach number
+        scalar M_;
+
+        //- Back pressure
+        scalar pBack_;
+
+        //- Model constant
+        scalar c1_;
+
+        //- Model constant area input
+        scalar A1_;
+
+        //- The name of the flux field
+        word phiName_;
+
+        //- The name of the rho field
+        word rhoName_;
+
+        //- The name of the U field
+        word UName_;
+
+        //- Choked
+        Switch choked_;
+
+        //- Relaxation factor
+        scalar relax_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("outletMachNumberPressure");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        outletMachNumberPressureFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        outletMachNumberPressureFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given outletMachNumberPressureFvPatchScalarField
+        //  onto a new patch
+        outletMachNumberPressureFvPatchScalarField
+        (
+            const outletMachNumberPressureFvPatchScalarField&,
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct as copy
+        outletMachNumberPressureFvPatchScalarField
+        (
+            const outletMachNumberPressureFvPatchScalarField&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchScalarField> clone() const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new outletMachNumberPressureFvPatchScalarField(*this)
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        outletMachNumberPressureFvPatchScalarField
+        (
+            const outletMachNumberPressureFvPatchScalarField&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct and return a clone setting internal field reference
+        virtual tmp<fvPatchScalarField> clone
+        (
+            const DimensionedField<scalar, volMesh>& iF
+        ) const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new outletMachNumberPressureFvPatchScalarField(*this, iF)
+            );
+        }
+
+
+    // Member functions
+
+        // Evaluation functions
+
+            //- Update the coefficients associated with the patch field
+            virtual void updateCoeffs();
+
+
+        //- Write
+        virtual void write(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

tutorials/combustion/reactingFoam/RAS/chockedNozzle/0/CH4

+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version         2;
+    format          ascii;
+    class           volScalarField;
+    location        "0";
+    object          CH4;
+}
+
+dimensions      [ 0 0 0 0 0 0 0 ];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    inletair
+    {
+        type            fixedValue;
+        value