Dual stream heat exchanger for porous media

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

@@ -77,16 +77,18 @@ Foam::porousZone::porousZone
     C1_(0),
     D_("D", dimensionSet(0, -2, 0, 0, 0), tensor::zero),
     F_("F", dimensionSet(0, -1, 0, 0, 0), tensor::zero),
+    rhoPri_(0),
+    Cpri_(0),
     Maux_(0),
+    Qaux_(0),
+    Caux_(0),
     Taux_(0),
-    caux_(0),
     Qepsilon_(0),
-    Tauxrelax_(0),
+    TauxRelax_(1),
     nCellsAuxInlet_(0),
     firstCell_(0),
     auxUnitVector_(vector::zero),
     nVerticalCells_(0)
-
 {
     Info<< "Creating porous zone: " << name_ << endl;
 
@@ -184,7 +186,7 @@ Foam::porousZone::porousZone
     // provide some feedback for the user
     // writeDict(Info, false);
 
-    // it is an error not to define anything
+    // It is an error not to define anything
     if
     (
         C0_ <= VSMALL
@@ -206,32 +208,41 @@ Foam::porousZone::porousZone
     if (const dictionary* dictPtr = dict_.subDictPtr("heatTransfer"))
     {
         Info<< "Reading porous heatTransfer: Maux and Taux" << nl;
+        dictPtr->lookup("rhoPri") >> rhoPri_;
+        dictPtr->lookup("Cpri") >> Cpri_;
+
         dictPtr->lookup("Maux") >> Maux_;
+        dictPtr->lookup("Qaux") >> Qaux_;
+        dictPtr->lookup("Caux") >> Caux_;
         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);
-    }
 
+        dictPtr->lookup("TauxRelax") >> TauxRelax_;
+
+        Info<< "Maux = " << Maux_
+            << ", Taux = " << Taux_
+            << ", Caux = " << Caux_
+            << ", nCellsAuxInlet = " << nCellsAuxInlet_
+            << ", Qepsilon = " << Qepsilon_
+            << endl;
+
+        if (mag(auxUnitVector_) < SMALL)
+        {
+            FatalIOErrorInFunction(dict)
+                << "auxUnitVector for heat transfer zone "
+                << name_ << " has zero length"
+                << exit(FatalIOError);
+        }
+        else
+        {
+            auxUnitVector_ /= mag(auxUnitVector_);
+        }
+    }
 }
 
 
@@ -330,6 +341,7 @@ void Foam::porousZone::addResistance
     }
 
     bool compressible = false;
+
     if (UEqn.dimensions() == dimensionSet(1, 1, -2, 0, 0))
     {
         compressible = true;
@@ -473,35 +485,60 @@ void Foam::porousZone::macroCellOrder
     Info<< "Creating cellsOrdered list for porous heat transfer zone "
         << name_ << nl << endl;
 
-    const labelList& cells = mesh_.cellZones()[cellZoneID_];
+    // Get current zone
-    const vectorField& cellsCellCenter = mesh_.cellCentres();
+    const cellZone& curZone = mesh_.cellZones()[cellZoneID_];
+
+    // Get cell centres
+    const vectorField& cellCentres = mesh_.C().internalField();
+
+    // Get cell-cell addressing
+    const labelListList& cc = mesh_.cellCells();
 
     scalarField& Tauxi = Taux.internalField();
 
-    label nCellsAuxInlet(nCellsAuxInlet_);
+    // Get inlet cells
-    label firstCell(firstCell_);
+    labelList cellsAuxInlet(nCellsAuxInlet_, -1);
-    label nVerticalCells(nVerticalCells_);
-    vector auxUnitVector(auxUnitVector_);
 
-    labelList cellsAuxInlet(nCellsAuxInlet, -1);
+    if (nCellsAuxInlet_ > 0)
-    cellsAuxInlet[0] = firstCell;
+    {
+        cellsAuxInlet[0] = firstCell_;
+    }
+    else
+    {
+        WarningInFunction
+            << "Number of aux inlet cells is set to zero: " << nCellsAuxInlet_
+            << "Reconsider the definition of macro cell order"
+            << endl;
+    }
+
+    // Creating horizontal list of cells, cellsAuxInlet[nCellsAuxInlet_]
+    label newCounter = 1;
+
+    // Note
+    // Terrible search algorithms
+    // Will work only in serial.  Rewrite required
+    // HJ, 17/Jul/2019
 
-    // - Creating horizontal list of cells, cellsAuxInlet[nCellsAuxInlet]
-    label newcounter = 1;
     forAll (cellsAuxInlet, i)
     {
-        const labelList& cellNb = mesh_.cellCells(cellsAuxInlet[i]);
+        const labelList& cellNb = cc[cellsAuxInlet[i]];
 
         forAll (cellNb, ii)
         {
             bool isInsideNb = false;
+
             forAll (cellsAuxInlet, iii)
             {
-                if (cellNb[ii] == cellsAuxInlet[iii]) isInsideNb = true;
+                if (cellNb[ii] == cellsAuxInlet[iii])
+                {
+                    isInsideNb = true;
+                    break;
+                }
             }
+
             if (!isInsideNb)
             {
-                if (mesh_.cellZones()[cellZoneID_].whichCell(cellNb[ii]) != -1)
+                if (curZone.whichCell(cellNb[ii]) != -1)
                 {
                     if
                     (
@@ -509,21 +546,21 @@ void Foam::porousZone::macroCellOrder
                         (
                             (
                                 (
-                                    cellsCellCenter[cellsAuxInlet[i]]
+                                    cellCentres[cellsAuxInlet[i]]
-                                  - cellsCellCenter[cellNb[ii]]
+                                  - cellCentres[cellNb[ii]]
                                 )/
                                 mag
                                 (
-                                    cellsCellCenter[cellsAuxInlet[i]]
+                                    cellCentres[cellsAuxInlet[i]]
-                                  - cellsCellCenter[cellNb[ii]]
+                                  - cellCentres[cellNb[ii]]
                                 )
                             )
-                          & auxUnitVector
+                          & auxUnitVector_
                         ) < 0.5
                     )
                     {
-                        cellsAuxInlet[newcounter] = cellNb[ii];
+                        cellsAuxInlet[newCounter] = cellNb[ii];
-                        ++newcounter;
+                        ++newCounter;
                     }
                 }
             }
@@ -532,42 +569,45 @@ void Foam::porousZone::macroCellOrder
 
     scalarField& Macroi = Macro.internalField();
 
-    labelList cellsOrdered(cells.size(), -1);
+    labelList cellsOrdered(curZone.size(), -1);
 
     // Writing horizontal cells into list cellsOrdered[cells.size()]
     forAll (cellsAuxInlet, i)
     {
-        cellsOrdered[i*nVerticalCells] = 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))
+        if ((i > 1) && (i % nVerticalCells_ == 0))
         {
             ++counter;
         }
 
-        const labelList& cellNb = mesh_.cellCells(cellsOrdered[i]);
+        const labelList& cellNb =cc[cellsOrdered[i]];
 
         forAll (cellNb, iii)
         {
             bool isInsideNb = false;
+
             forAll (cellsOrdered, iiii)
             {
                 if (cellNb[iii] == cellsOrdered[iiii])
                 {
                     isInsideNb = true;
+                    break;
                 }
             }
 
             if (!isInsideNb)
             {
-                if (mesh_.cellZones()[cellZoneID_].whichCell(cellNb[iii]) != -1)
+                if (curZone.whichCell(cellNb[iii]) != -1)
                 {
                     if
                     (
@@ -575,16 +615,16 @@ void Foam::porousZone::macroCellOrder
                         (
                             (
                                 (
-                                    cellsCellCenter[cellsOrdered[i]]
+                                    cellCentres[cellsOrdered[i]]
-                                  - cellsCellCenter[cellNb[iii]]
+                                  - cellCentres[cellNb[iii]]
                                 )/
                                 mag
                                 (
-                                    cellsCellCenter[cellsOrdered[i]]
+                                    cellCentres[cellsOrdered[i]]
-                                  - cellsCellCenter[cellNb[iii]]
+                                  - cellCentres[cellNb[iii]]
                                 )
                             )
-                          & auxUnitVector
+                          & auxUnitVector_
                         ) > 0.85
                     )
                     {
@@ -602,9 +642,11 @@ void Foam::porousZone::macroCellOrder
     forAll (cellsOrdered, i)
     {
         const labelList& cellFaces = mesh_.cells()[cellsOrdered[i]];
-        forAll (cellFaces,ii)
+
+        forAll (cellFaces, ii)
         {
             label faceI = cellFaces[ii];
+
             if (mesh_.isInternalFace(faceI))
             {
                 if (mesh_.faceOwner()[faceI] == cellsOrdered[i])
@@ -626,15 +668,10 @@ void Foam::porousZone::macroCellOrder
             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]] +=
@@ -645,6 +682,7 @@ void Foam::porousZone::macroCellOrder
     }
 }
 
+
 void Foam::porousZone::writeDict(Ostream& os, bool subDict) const
 {
     if (subDict)

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

@@ -39,13 +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
+    In the porous zone a heat source using the dual stream heat exchanger
     approach:
 
     heatTransfer
     @f[
-        Qtot = const            //[J/s]
+        Maux = Mass flow of aux fluid                         // [kg/s]
-        Tref = const            //[K]
+        Qaux = Heat to be eliminated                          // [W]
+        Caux = Specific heat capacity of aux fluid            // [J / kg K]
+        Taux = initial guess for aux fluid inlet temperature  // [K]
     @f]
 
     Darcy-Forchheimer (@e d and @e f parameters)
@@ -138,35 +140,44 @@ class porousZone
         //- Forchheimer coefficient
         dimensionedTensor F_;
 
+        //- Density of primary fluid
+        scalar rhoPri_;
+
+        //- Specific heat capacity of primary fluid
+        scalar Cpri_;
+    
         //- Mass flow of aux fluid
         scalar Maux_;
 
-        //- Inlet temperature of aux fluid
+        //- Heat to be eliminated
-        scalar Taux_;
+        scalar Qaux_;
 
         //- Specific heat capacity of coolant
-        scalar caux_;
+        scalar Caux_;
+
+        //- Inlet temperature of aux fluid
+        mutable scalar Taux_;
 
         //- Epsilon for Single Effectiveness HX
         scalar Qepsilon_;
 
         //- Relaxation factor for Taux
-        scalar Tauxrelax_;
+        scalar TauxRelax_;
 
 
         // Geometrical description of the porous heat exchanger
 
-        //- Number of inlet cells for aux fluid
+            //- Number of inlet cells for aux fluid
-        label nCellsAuxInlet_;
+            label nCellsAuxInlet_;
 
-        //- First inlet cell for aux fluid
+            //- First inlet cell for aux fluid
-        label firstCell_;
+            label firstCell_;
 
-        //- Inlet direction for auxiliary fluid
+            //- Inlet direction for auxiliary fluid
-        vector auxUnitVector_;
+            vector auxUnitVector_;
 
-        //- Number of vertical cells for auxiliary fluid
+            //- Number of vertical cells for auxiliary fluid
-        label nVerticalCells_;
+            label nVerticalCells_;
 
 
     // Private Member Functions

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

@@ -33,12 +33,12 @@ void Foam::porousZone::modifyDdt(fvMatrix<Type>& m) const
 {
     if (porosity_ < 1)
     {
-        const labelList& cells = mesh_.cellZones()[cellZoneID_];
+        const labelList& zoneCells = mesh_.cellZones()[cellZoneID_];
 
-        forAll(cells, i)
+        forAll (zoneCells, i)
         {
-            m.diag()[cells[i]]   *= porosity_;
+            m.diag()[zoneCells[i]] *= porosity_;
-            m.source()[cells[i]] *= porosity_;
+            m.source()[zoneCells[i]] *= porosity_;
         }
     }
 }
@@ -150,68 +150,86 @@ void Foam::porousZone::addHeatSource
     const scalarField dT = Tauxi - T;
 
     const label nMacro(nCellsAuxInlet_*nVerticalCells_);
-    scalarList Tmacro(nMacro, 0.0);
+    scalarField Tmacro(nMacro, scalar(0));
-    scalarList dTaux(nMacro, 0.0);
+    scalarField dTaux(nMacro, scalar(0));
-    scalar QSum(0.0);
+    scalar QSum = 0;
 
-    const scalar Taux_relax(Tauxrelax_);
+    // Mass flow rate for individual channels
-    const scalar c_aux(caux_);
+    const scalar qmAuxi = Maux_/nCellsAuxInlet_;
-    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_];
+    // Get zone cells
+    const labelList& zoneCells = mesh_.cellZones()[cellZoneID_];
 
-    forAll(cells, i)
+    forAll (zoneCells, i)
     {
-        scalar Qcell = Qeps*rho_pri*c_pri*posFlux[cells[i]]*dT[cells[i]];
+        scalar Qcell =
+            Qepsilon_*rhoPri_*Cpri_*posFlux[zoneCells[i]]*dT[zoneCells[i]];
 
-        // heat in each macro(cell)
+        // Heat in each macro (cell)
-        Qauxi[cells[i]] = Qcell;
+        Qauxi[zoneCells[i]] = Qcell;
 
         // make an int out of a macro
-        const int macro = Macro[cells[i]];
+        const int macro = Macro[zoneCells[i]];
 
         // deltaTaux in each macro(cell)
-        dTaux[macro] = Qcell/(c_aux*qm_auxi);
+        dTaux[macro] = Qcell/(Caux_*qmAuxi);
 
         // adding Heat to equation
-        QSource[cells[i]] += Qcell/(rho_pri*c_pri);
+        QSource[zoneCells[i]] += Qcell/(rhoPri_*Cpri_);
 
         // summing for total heat of HX
         QSum += Qcell;
     }
 
-    reduce(dTaux, sumOp<scalarList>());
+    reduce(dTaux, sumOp<scalarField>());
-    Tmacro[0] = T_aux;
+
+    Tmacro[0] = Taux_;
+
     forAll (Tmacro, i)
     {
         if (i > 0)
         {
-            Tmacro[i] = Tmacro[i-1] - dTaux[i-1];
+            if (i % nVerticalCells_ == 0)
-        }
+            {
-
+                Tmacro[i] = Taux_;
-        if ((i > 0) && (i % nVerticalCells_ == 0))
+            }
-        {
+            else
-            Tmacro[i] = T_aux;
+            {
+                Tmacro[i] = Tmacro[i - 1] - dTaux[i - 1];
+            }
         }
     }
 
     reduce(QSum, sumOp<scalar>());
-    Info << "Heat exchanger: " << name_ << endl;
-    Info << "Q = " << QSum << endl;
-    Info << "deltaT = " << QSum/(qm_aux*c_aux) << endl;
 
-    forAll(cells, i)
+    // Adjust Taux_ to match the specified transferred heat
+
+    scalar deltaTAux = (Qaux_ - QSum)/(Maux_*Caux_);
+
+    Info<< "Heat exchanger: " << name_
+        << ": Q = " << QSum
+        << " Taux = " << Taux_
+        << " delta Taux = " << deltaTAux
+        << endl;
+
+    // Relax the heat source
+    forAll (zoneCells, i)
     {
-        const int macro = Macro[cells[i]];
+        const int macro = Macro[zoneCells[i]];
-        Tauxi[cells[i]] = Tauxi_old[cells[i]]
+
-            // upwind scheme for Aux fluid (Tmacro = inlet temp)
+        Tauxi[zoneCells[i]] = Tauxi_old[zoneCells[i]]
-          + Taux_relax*(Tmacro[macro] - Tauxi_old[cells[i]]);
+          + TauxRelax_*(Tmacro[macro] - Tauxi_old[zoneCells[i]]);
     }
+
+    // Note: need better relaxation to speed up convergence close
+    // to the matching value, when deltaTAux -> 0
+    // HJ, 17/Jul/2019
+    
+    // Limit deltaTAux to 10% of Taux_
+    deltaTAux = sign(deltaTAux)*Foam::min(0.1*Taux_, mag(deltaTAux));
+    
+    // Update Taux for next iteration
+    Taux_ += TauxRelax_*deltaTAux;
 }
 
 

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

@@ -47,8 +47,10 @@ Description
         }
         heatTransfer
         {
-            Qtot      50;        //[J/s]
+            Maux = Mass flow of aux fluid                         // [kg/s]
-            Tref      273.15;    //[K]
+            Qaux = Heat to be eliminated                          // [W]
+            Caux = Specific heat capacity of aux fluid            // [J / kg K]
+            Taux = initial guess for aux fluid inlet temperature  // [K]
         }
     }
     )
@@ -168,7 +170,7 @@ public:
             const volScalarField& posFlux
         ) const;
 
-        //- Order cells for Dual Stream model
+        //- Order cells for the Dual Stream model
         void macroCellOrder
         (
             volScalarField& Taux,