diff --git a/src/finiteVolume/cfdTools/general/porousMedia/porousZone.C b/src/finiteVolume/cfdTools/general/porousMedia/porousZone.C index fb137ca7b..3ce82f8d2 100644 --- a/src/finiteVolume/cfdTools/general/porousMedia/porousZone.C +++ b/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_]; - const vectorField& cellsCellCenter = mesh_.cellCentres(); + // Get current zone + 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_); - label firstCell(firstCell_); - label nVerticalCells(nVerticalCells_); - vector auxUnitVector(auxUnitVector_); + // Get inlet cells + labelList cellsAuxInlet(nCellsAuxInlet_, -1); - labelList cellsAuxInlet(nCellsAuxInlet, -1); - cellsAuxInlet[0] = firstCell; + if (nCellsAuxInlet_ > 0) + { + 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]] - - cellsCellCenter[cellNb[ii]] + cellCentres[cellsAuxInlet[i]] + - cellCentres[cellNb[ii]] )/ mag ( - cellsCellCenter[cellsAuxInlet[i]] - - cellsCellCenter[cellNb[ii]] + cellCentres[cellsAuxInlet[i]] + - cellCentres[cellNb[ii]] ) ) - & auxUnitVector + & auxUnitVector_ ) < 0.5 ) { - cellsAuxInlet[newcounter] = cellNb[ii]; - ++newcounter; + cellsAuxInlet[newCounter] = cellNb[ii]; + ++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]] - - cellsCellCenter[cellNb[iii]] + cellCentres[cellsOrdered[i]] + - cellCentres[cellNb[iii]] )/ mag ( - cellsCellCenter[cellsOrdered[i]] - - cellsCellCenter[cellNb[iii]] + cellCentres[cellsOrdered[i]] + - 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) diff --git a/src/finiteVolume/cfdTools/general/porousMedia/porousZone.H b/src/finiteVolume/cfdTools/general/porousMedia/porousZone.H index 30c0951bf..bfd2208a6 100644 --- a/src/finiteVolume/cfdTools/general/porousMedia/porousZone.H +++ b/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] - Tref = const //[K] + Maux = Mass flow of aux fluid // [kg/s] + 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 - scalar Taux_; + //- Heat to be eliminated + 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 - label nCellsAuxInlet_; + //- Number of inlet cells for aux fluid + label nCellsAuxInlet_; - //- First inlet cell for aux fluid - label firstCell_; + //- First inlet cell for aux fluid + label firstCell_; - //- Inlet direction for auxiliary fluid - vector auxUnitVector_; + //- Inlet direction for auxiliary fluid + vector auxUnitVector_; - //- Number of vertical cells for auxiliary fluid - label nVerticalCells_; + //- Number of vertical cells for auxiliary fluid + label nVerticalCells_; // Private Member Functions diff --git a/src/finiteVolume/cfdTools/general/porousMedia/porousZoneTemplates.C b/src/finiteVolume/cfdTools/general/porousMedia/porousZoneTemplates.C index f8ea6f152..c2ff057aa 100644 --- a/src/finiteVolume/cfdTools/general/porousMedia/porousZoneTemplates.C +++ b/src/finiteVolume/cfdTools/general/porousMedia/porousZoneTemplates.C @@ -33,12 +33,12 @@ void Foam::porousZone::modifyDdt(fvMatrix& 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.source()[cells[i]] *= porosity_; + m.diag()[zoneCells[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); - scalarList dTaux(nMacro, 0.0); - scalar QSum(0.0); + scalarField Tmacro(nMacro, scalar(0)); + scalarField dTaux(nMacro, scalar(0)); + scalar QSum = 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_); + // Mass flow rate for individual channels + const scalar qmAuxi = Maux_/nCellsAuxInlet_; - 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) - Qauxi[cells[i]] = Qcell; + // Heat in each macro (cell) + 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()); - Tmacro[0] = T_aux; + reduce(dTaux, sumOp()); + + Tmacro[0] = Taux_; + forAll (Tmacro, i) { if (i > 0) { - Tmacro[i] = Tmacro[i-1] - dTaux[i-1]; - } - - if ((i > 0) && (i % nVerticalCells_ == 0)) - { - Tmacro[i] = T_aux; + if (i % nVerticalCells_ == 0) + { + Tmacro[i] = Taux_; + } + else + { + Tmacro[i] = Tmacro[i - 1] - dTaux[i - 1]; + } } } reduce(QSum, sumOp()); - 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]]; - Tauxi[cells[i]] = Tauxi_old[cells[i]] - // upwind scheme for Aux fluid (Tmacro = inlet temp) - + Taux_relax*(Tmacro[macro] - Tauxi_old[cells[i]]); + const int macro = Macro[zoneCells[i]]; + + Tauxi[zoneCells[i]] = Tauxi_old[zoneCells[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; } diff --git a/src/finiteVolume/cfdTools/general/porousMedia/porousZones.H b/src/finiteVolume/cfdTools/general/porousMedia/porousZones.H index 29f7bb584..f7bf2c116 100644 --- a/src/finiteVolume/cfdTools/general/porousMedia/porousZones.H +++ b/src/finiteVolume/cfdTools/general/porousMedia/porousZones.H @@ -47,8 +47,10 @@ Description } heatTransfer { - Qtot 50; //[J/s] - Tref 273.15; //[K] + Maux = Mass flow of aux fluid // [kg/s] + 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,