Added thermal support in the porous zone. Vanja Skuric

src/finiteVolume/cfdTools/general/porousMedia/porousZone.C

@@ -76,7 +76,17 @@ Foam::porousZone::porousZone
     C0_(0),
     C1_(0),
     D_("D", dimensionSet(0, -2, 0, 0, 0), tensor::zero),
-    F_("F", dimensionSet(0, -1, 0, 0, 0), tensor::zero)
+    F_("F", dimensionSet(0, -1, 0, 0, 0), tensor::zero),
+    Maux_(0),
+    Taux_(0),
+    caux_(0),
+    Qepsilon_(0),
+    Tauxrelax_(0),
+    nCellsAuxInlet_(0),
+    firstCell_(0),
+    auxUnitVector_(vector::zero),
+    nVerticalCells_(0)
+
 {
     Info<< "Creating porous zone: " << name_ << endl;
 
@@ -191,6 +201,37 @@ Foam::porousZone::porousZone
                "nor Darcy-Forchheimer law (d/f) specified"
             << exit(FatalIOError);
     }
+
+    // Heat rate value and temperature of heat exchanger
+    if (const dictionary* dictPtr = dict_.subDictPtr("heatTransfer"))
+    {
+        Info<< "Reading porous heatTransfer: Maux and Taux" << nl;
+        dictPtr->lookup("Maux") >> Maux_;
+        dictPtr->lookup("Taux") >> Taux_;
+        dictPtr->lookup("caux") >> caux_;
+        dictPtr->lookup("Qepsilon") >> Qepsilon_;
+        dictPtr->lookup("firstCell") >> firstCell_;
+        dictPtr->lookup("auxUnitVector") >> auxUnitVector_;
+        dictPtr->lookup("nVerticalCells") >> nVerticalCells_;
+        dictPtr->lookup("nCellsAuxInlet") >> nCellsAuxInlet_;
+        dictPtr->lookup("TauxRelax") >> Tauxrelax_;
+        Info<< "Maux = " << Maux_ <<
+        " ,Taux = " << Taux_ <<
+        " ,caux = " << caux_ <<
+        " ,nCellsAuxInlet = " << nCellsAuxInlet_ <<
+        " ,Qepsilon = " << Qepsilon_ << nl;
+    }
+    else
+    {
+        FatalIOErrorIn
+        (
+            "Foam::porousZone::porousZone"
+            "(const fvMesh&, const word&, const dictionary&)",
+            dict_
+        )   << "\"heatTransfer\" dictionary not specified"
+            << exit(FatalIOError);
+    }
+
 }
 
 
@@ -360,6 +401,245 @@ void Foam::porousZone::addResistance
 }
 
 
+void Foam::porousZone::addHeatResistance
+(
+    fvScalarMatrix& hTEqn,
+    const volScalarField& T,
+    volScalarField& Taux,
+    volScalarField& Qaux,
+    const volVectorField& U,
+    const volScalarField& Macro,
+    const volScalarField& posFlux
+) const
+{
+    if (cellZoneID_ == -1 && mag(Maux_) < SMALL)
+    {
+        return;
+    }
+
+    scalarField& QSource = hTEqn.source();
+    const scalarField& Ti = T.internalField();
+    scalarField& Tauxi = Taux.internalField();
+    scalarField& Qauxi = Qaux.internalField();
+    const vectorField& Ui = U.internalField();
+    const scalarField& Macroi = Macro.internalField();
+    const scalarField& posFluxi = posFlux.internalField();
+
+
+
+    if (hTEqn.dimensions() == dimensionSet(1, -2, -3, 0, 0))
+    {
+        addHeatSource
+        (
+            Macroi,
+            posFluxi,
+            QSource,
+            Qauxi,
+            Ti,
+            Tauxi,
+            Ui,
+            mesh_.lookupObject<volScalarField>("rho")
+        );
+    }
+    else if (hTEqn.dimensions() == dimensionSet(0, 3, -1, 1, 0))
+    {
+        addHeatSource
+        (
+            Macroi,
+            posFluxi,
+            QSource,
+            Qauxi,
+            Ti,
+            Tauxi,
+            Ui,
+            geometricOneField()
+        );
+    }
+    else
+    {
+        Info<< "No hEqn or TEqn, exiting" << nl;
+        return;
+    }
+}
+
+void Foam::porousZone::macroCellOrder
+(
+    volScalarField& Taux,
+    volScalarField& Macro,
+    volScalarField& posFlux,
+    const surfaceScalarField& phi
+) const
+{
+    Info << "Creating cellsOrdered list " << nl << endl;
+
+    const labelList& cells = mesh_.cellZones()[cellZoneID_];
+    const vectorField& cellsCellCenter = mesh_.cellCentres();
+
+    scalarField& Tauxi = Taux.internalField();
+
+    label nCellsAuxInlet(nCellsAuxInlet_);
+    label firstCell(firstCell_);
+    label nVerticalCells(nVerticalCells_);
+    vector auxUnitVector(auxUnitVector_);
+
+    labelList cellsAuxInlet(nCellsAuxInlet, -1);
+    cellsAuxInlet[0] = firstCell;
+
+    // - Creating horizontal list of cells, cellsAuxInlet[nCellsAuxInlet]
+    label newcounter = 1;
+    forAll (cellsAuxInlet, i)
+    {
+        const labelList& cellNb = mesh_.cellCells(cellsAuxInlet[i]);
+
+        forAll (cellNb, ii)
+        {
+            bool isInsideNb = false;
+            forAll (cellsAuxInlet, iii)
+            {
+                if (cellNb[ii] == cellsAuxInlet[iii]) isInsideNb = true;
+            }
+            if (!isInsideNb)
+            {
+                if (mesh_.cellZones()[cellZoneID_].whichCell(cellNb[ii]) != -1)
+                {
+                    if
+                    (
+                        mag
+                        (
+                            (
+                                (
+                                    cellsCellCenter[cellsAuxInlet[i]]
+                                  - cellsCellCenter[cellNb[ii]]
+                                )/
+                                mag
+                                (
+                                    cellsCellCenter[cellsAuxInlet[i]]
+                                  - cellsCellCenter[cellNb[ii]]
+                                )
+                            )
+                          & auxUnitVector
+                        ) < 0.5
+                    )
+                    {
+                        cellsAuxInlet[newcounter] = cellNb[ii];
+                        ++newcounter;
+                    }
+                }
+            }
+        }
+    }
+
+    scalarField& Macroi = Macro.internalField();
+
+    labelList cellsOrdered(cells.size(), -1);
+
+    // Writing horizontal cells into list cellsOrdered[cells.size()]
+    forAll (cellsAuxInlet, i)
+    {
+        cellsOrdered[i*nVerticalCells] = cellsAuxInlet[i];
+    }
+
+    // Creating cellsOrdered list of ordered horizontal cells
+    label counter = 1;
+    forAll (cellsOrdered, i)
+    {
+        Macroi[cellsOrdered[i]] = i;
+        Tauxi[cellsOrdered[i]] = Taux_;
+        
+        if ((i > 1) && (i % nVerticalCells == 0))
+        {
+            ++counter;
+        }
+
+        const labelList& cellNb = mesh_.cellCells(cellsOrdered[i]);
+
+        forAll (cellNb, iii)
+        {
+            bool isInsideNb = false;
+            forAll (cellsOrdered, iiii)
+            {
+                if (cellNb[iii] == cellsOrdered[iiii])
+                {
+                    isInsideNb = true;
+                }
+            }
+
+            if (!isInsideNb)
+            {
+                if (mesh_.cellZones()[cellZoneID_].whichCell(cellNb[iii]) != -1)
+                {
+                    if
+                    (
+                        mag
+                        (
+                            (
+                                (
+                                    cellsCellCenter[cellsOrdered[i]]
+                                  - cellsCellCenter[cellNb[iii]]
+                                )/
+                                mag
+                                (
+                                    cellsCellCenter[cellsOrdered[i]]
+                                  - cellsCellCenter[cellNb[iii]]
+                                )
+                            )
+                          & auxUnitVector
+                        ) > 0.85
+                    )
+                    {
+                        cellsOrdered[counter] = cellNb[iii];
+                        ++counter;
+                    }
+                }
+            }
+        }
+    }
+
+    scalarField& posFluxi = posFlux.internalField();
+    // - Calculating mass flow through each cell
+    forAll (cellsOrdered, i)
+    {
+        const labelList& cellFaces = mesh_.cells()[cellsOrdered[i]];
+        forAll (cellFaces,ii)
+        {
+            label faceI = cellFaces[ii];
+            if (mesh_.isInternalFace(faceI))
+            {
+                if (mesh_.faceOwner()[faceI] == cellsOrdered[i])
+                {
+                    if (phi.internalField()[faceI] > 0.0)
+                    {
+                        posFluxi[cellsOrdered[i]] += phi.internalField()[faceI];
+                    }
+                }
+                else
+                {
+                    if (phi.internalField()[faceI] < 0.0)
+                    {
+                        posFluxi[cellsOrdered[i]] += -phi.internalField()[faceI];
+                    }
+                }
+            }
+            else
+            {
+                const label patchI = mesh_.boundaryMesh().whichPatch(faceI);
+                const label faceIL = mesh_.boundaryMesh()[patchI].whichFace(faceI);
+
+                if (patchI < 0)
+                {
+                    Info << "patchI < 0 " << endl;
+                    return;
+                }
+
+                if (phi.boundaryField()[patchI][faceIL] > 0.0)
+                {
+                    posFluxi[cellsOrdered[i]] += phi.boundaryField()[patchI][faceIL];
+                }
+            }
+        }
+    }
+}
+
 void Foam::porousZone::writeDict(Ostream& os, bool subDict) const
 {
     if (subDict)
@@ -402,6 +682,13 @@ void Foam::porousZone::writeDict(Ostream& os, bool subDict) const
         dictPtr->write(os);
     }
 
+    // Heat transfer inputs
+    if (const dictionary* dictPtr = dict_.subDictPtr("heatTransfer"))
+    {
+        os << indent << "heatTransfer";
+        dictPtr->write(os);
+    }
+
     os << decrIndent << indent << token::END_BLOCK << endl;
 }
 

src/finiteVolume/cfdTools/general/porousMedia/porousZone.H

@@ -39,6 +39,15 @@ Description
         S = - \rho C_0 |U|^{(C_1 - 1)/2} U
     @f]
 
+    In the porous zone a heat source using the single stream heat exchanger
+    approach:
+
+    heatTransfer
+    @f[
+        Qtot = const            //[J/s]
+        Tref = const            //[K]
+    @f]
+
     Darcy-Forchheimer (@e d and @e f parameters)
     @f[
         S = - (\mu \, d + \frac{\rho |U|}{2} \, f) U
@@ -129,6 +138,36 @@ class porousZone
         //- Forchheimer coefficient
         dimensionedTensor F_;
 
+        //- Mass flow of aux fluid
+        scalar Maux_;
+
+        //- Inlet temperature of aux fluid
+        scalar Taux_;
+
+        //- Specific heat capacity of coolant
+        scalar caux_;
+
+        //- Epsilon for Single Effectiveness HX
+        scalar Qepsilon_;
+
+        //- Relaxation factor for Taux
+        scalar Tauxrelax_;
+
+
+        // Geometrical description of the porous heat exchanger
+
+        //- Number of inlet cells for aux fluid
+        label nCellsAuxInlet_;
+
+        //- First inlet cell for aux fluid
+        label firstCell_;
+
+        //- Inlet direction for auxiliary fluid
+        vector auxUnitVector_;
+
+        //- Number of vertical cells for auxiliary fluid
+        label nVerticalCells_;
+
 
     // Private Member Functions
 
@@ -159,7 +198,6 @@ class porousZone
             const vectorField& U
         ) const;
 
-
         //- Power-law resistance
         template<class RhoFieldType>
         void addPowerLawResistance
@@ -181,6 +219,20 @@ class porousZone
             const vectorField& U
         ) const;
 
+        //- Heat source
+        template<class RhoFieldType>
+        void addHeatSource
+        (
+            const scalarField& Macro,
+            const scalarField& posFlux,
+            scalarField& QSource,
+            scalarField& Qaux,
+            const scalarField& T,
+            scalarField& Taux,
+            const vectorField& U,
+            const RhoFieldType& rho
+        ) const;
+
 
         //- Disallow default bitwise copy construct
         porousZone(const porousZone&);
@@ -314,6 +366,27 @@ public:
             bool correctAUprocBC = true
         ) const;
 
+        //- Add the heat source in enthalpy equation
+        void addHeatResistance
+        (
+            fvScalarMatrix& hTEqn,
+            const volScalarField& T,
+            volScalarField& Taux,
+            volScalarField& Qaux,
+            const volVectorField& U,
+            const volScalarField& Macro,
+            const volScalarField& posFlux
+        ) const;
+
+        //- Order cells for Dual Stream model
+        void macroCellOrder
+        (
+            volScalarField& Taux,
+            volScalarField& Macro,
+            volScalarField& posFlux,
+            const surfaceScalarField& phi
+        ) const;
+
         //- Write the porousZone dictionary
         virtual void writeDict(Ostream&, bool subDict = true) const;
 

src/finiteVolume/cfdTools/general/porousMedia/porousZoneTemplates.C

@@ -132,4 +132,72 @@ void Foam::porousZone::addViscousInertialResistance
 }
 
 
+template<class RhoFieldType>
+void Foam::porousZone::addHeatSource
+(
+    const scalarField& Macro,
+    const scalarField& posFlux,
+    scalarField& QSource,
+    scalarField& Qauxi,
+    const scalarField& T,
+    scalarField& Tauxi,
+    const vectorField& U,
+    const RhoFieldType& rho
+) const
+{
+    // - Dual stream model (Single Effectiveness Model)
+    const scalarField Tauxi_old = Tauxi; // create storePrevIter()
+    const scalarField dT = Tauxi - T;
+
+
+    const label nMacro(nCellsAuxInlet_*nVerticalCells_);
+    scalarList Tmacro(nMacro, 0.0);
+    scalarList dTaux(nMacro, 0.0);
+    scalar QSum(0.0);
+
+    const scalar Taux_relax(Tauxrelax_);
+    const scalar c_aux(caux_);
+    const scalar c_pri(1009.0);
+    const scalar qm_aux(Maux_);
+    const scalar qm_auxi = qm_aux/nCellsAuxInlet_;
+    const scalar T_aux(Taux_);
+    const scalar rho_pri(1.1021);
+    const scalar Qeps(Qepsilon_);
+
+    const labelList& cells = mesh_.cellZones()[cellZoneID_];
+
+    forAll(cells, i)
+    {
+        scalar Qcell = Qeps*rho_pri*c_pri*posFlux[cells[i]]*dT[cells[i]];
+
+        Qauxi[cells[i]] = Qcell;                            // heat in each macro(cell)
+        const int macro = Macro[cells[i]];
+        dTaux[macro] = Qcell/(c_aux*qm_auxi);             // deltaTaux in each macro(cell)
+        QSource[cells[i]] += Qcell/(rho_pri*c_pri);       // adding Heat to equation
+        QSum += Qcell;                                    // summing for total heat of HX
+    }
+
+    reduce(dTaux, sumOp<scalarList>());
+    Tmacro[0] = T_aux;
+    forAll (Tmacro, i)
+    {
+        if (i > 0) Tmacro[i] = Tmacro[i-1] - dTaux[i-1];
+        if ((i > 0) && (i % nVerticalCells_ == 0)) Tmacro[i] = T_aux;
+    }
+
+    reduce(QSum, sumOp<scalar>());
+    Info << "Heat exchanger: " << name_ << endl;
+    Info << "Q = " << QSum << endl;
+    Info << "deltaT = " << QSum/(qm_aux*c_aux) << endl;
+
+    forAll(cells, i)
+    {
+        const int macro = Macro[cells[i]];
+        Tauxi[cells[i]] = Tauxi_old[cells[i]]
+            // upwind scheme for Aux fluid (Tmacro = inlet temp)
+          + Taux_relax*(Tmacro[macro] - Tauxi_old[cells[i]]);
+    }
+}
+
+
 // ************************************************************************* //

src/finiteVolume/cfdTools/general/porousMedia/porousZones.C

@@ -88,6 +88,38 @@ void Foam::porousZones::addResistance
 }
 
 
+void Foam::porousZones::addHeatResistance
+(
+    fvScalarMatrix& hTEqn,
+    const volScalarField& T,
+    volScalarField& Taux,
+    volScalarField& Qaux,
+    const volVectorField& U,
+    const volScalarField& Macro,
+    const volScalarField& posFlux
+) const
+{
+    forAll(*this, i)
+    {
+        operator[](i).addHeatResistance(hTEqn, T, Taux, Qaux, U, Macro, posFlux);
+    }
+}
+
+//- Order cells for Dual Stream model
+void Foam::porousZones::macroCellOrder
+(
+    volScalarField& Taux,
+    volScalarField& Macro,
+    volScalarField& posFlux,
+    const surfaceScalarField& phi
+) const
+{
+    forAll(*this, i)
+    {
+        operator[](i).macroCellOrder(Taux, Macro, posFlux, phi);
+    }
+}
+
 bool Foam::porousZones::readData(Istream& is)
 {
     clear();

src/finiteVolume/cfdTools/general/porousMedia/porousZones.H

@@ -45,6 +45,11 @@ Description
             d   d [0 -2 0 0 0]  (-1000 -1000 0.50753e+08);
             f   f [0 -1 0 0 0]  (-1000 -1000 12.83);
         }
+        heatTransfer
+        {
+            Qtot      50;        //[J/s]
+            Tref      273.15;    //[K]
+        }
     }
     )
     @endverbatim
@@ -150,6 +155,28 @@ public:
             volTensorField& AU
         ) const;
 
+        //- Add the heat source contribution to enalphy equation (energy
+        // equation)
+        void addHeatResistance
+        (
+            fvScalarMatrix& hTEqn,
+            const volScalarField& T,
+            volScalarField& Taux,
+            volScalarField& Qaux,
+            const volVectorField& U,
+            const volScalarField& Macro,
+            const volScalarField& posFlux
+        ) const;
+
+        //- Order cells for Dual Stream model
+        void macroCellOrder
+        (
+            volScalarField& Taux,
+            volScalarField& Macro,
+            volScalarField& posFlux,
+            const surfaceScalarField& phi
+        ) const;
+
         //- read modified data
         virtual bool readData(Istream&);