From b6f44b333852dfd77c2ff62f25b445cfdcdc8bd2 Mon Sep 17 00:00:00 2001 From: asanaz Date: Mon, 17 Aug 2020 14:12:26 +0200 Subject: [PATCH] recurrence solver for modal decomposition of flow fleids --- .../solvers/rctfSpeciesTransport/CEq.H | 19 ++ .../solvers/rctfSpeciesTransport/Make/files | 3 + .../solvers/rctfSpeciesTransport/Make/options | 27 ++ .../rctfSpeciesTransport/continuityErrCalc.H | 17 ++ .../rctfSpeciesTransport/createFields.H | 281 ++++++++++++++++++ .../rctfSpeciesTransport/createRecBase.H | 66 ++++ .../rctfSpeciesTransport.C | 131 ++++++++ .../solvers/rctfSpeciesTransport/readFields.H | 40 +++ etc/solver-list.txt | 1 + 9 files changed, 585 insertions(+) create mode 100755 applications/solvers/rctfSpeciesTransport/CEq.H create mode 100755 applications/solvers/rctfSpeciesTransport/Make/files create mode 100755 applications/solvers/rctfSpeciesTransport/Make/options create mode 100644 applications/solvers/rctfSpeciesTransport/continuityErrCalc.H create mode 100755 applications/solvers/rctfSpeciesTransport/createFields.H create mode 100644 applications/solvers/rctfSpeciesTransport/createRecBase.H create mode 100755 applications/solvers/rctfSpeciesTransport/rctfSpeciesTransport.C create mode 100755 applications/solvers/rctfSpeciesTransport/readFields.H diff --git a/applications/solvers/rctfSpeciesTransport/CEq.H b/applications/solvers/rctfSpeciesTransport/CEq.H new file mode 100755 index 00000000..53265ae9 --- /dev/null +++ b/applications/solvers/rctfSpeciesTransport/CEq.H @@ -0,0 +1,19 @@ + +volScalarField alphaEff("alphaEff", turbulence->nu()/Sc + alphat); + +CEqn = +( + fvm::ddt(C) + + fvm::div(phiRec, C) + - fvm::laplacian(alphaEff, C) + == + fvOptions(C) +); + +CEqn.relax(relaxCoeff); + +fvOptions.constrain(CEqn); + +CEqn.solve(); + +fvOptions.correct(C); diff --git a/applications/solvers/rctfSpeciesTransport/Make/files b/applications/solvers/rctfSpeciesTransport/Make/files new file mode 100755 index 00000000..e803d38c --- /dev/null +++ b/applications/solvers/rctfSpeciesTransport/Make/files @@ -0,0 +1,3 @@ +rctfSpeciesTransport.C + +EXE=$(CFDEM_APP_DIR)/rctfSpeciesTransport diff --git a/applications/solvers/rctfSpeciesTransport/Make/options b/applications/solvers/rctfSpeciesTransport/Make/options new file mode 100755 index 00000000..686f676c --- /dev/null +++ b/applications/solvers/rctfSpeciesTransport/Make/options @@ -0,0 +1,27 @@ +include $(CFDEM_ADD_LIBS_DIR)/additionalLibs + +EXE_INC = \ + -I$(CFDEM_OFVERSION_DIR) \ + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/meshTools/lnInclude \ + -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \ + -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \ + -I$(LIB_SRC)/transportModels \ + -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \ + -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/lnInclude \ + -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/cfdTools \ + -I$(CFDEM_SRC_DIR)/recurrence/lnInclude \ + -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/derived/cfdemCloudRec \ + +EXE_LIBS = \ + -L$(CFDEM_LIB_DIR)\ + -lrecurrence \ + -lturbulenceModels \ + -lincompressibleTurbulenceModels \ + -lincompressibleTransportModels \ + -lfiniteVolume \ + -lmeshTools \ + -lfvOptions \ + -l$(CFDEM_LIB_NAME) \ + $(CFDEM_ADD_LIB_PATHS) \ + $(CFDEM_ADD_LIBS) diff --git a/applications/solvers/rctfSpeciesTransport/continuityErrCalc.H b/applications/solvers/rctfSpeciesTransport/continuityErrCalc.H new file mode 100644 index 00000000..a0b5b991 --- /dev/null +++ b/applications/solvers/rctfSpeciesTransport/continuityErrCalc.H @@ -0,0 +1,17 @@ +// calculate the continuity error according to phiRec + +{ + volScalarField contErr(fvc::div(phiRec)); + + scalar sumLocalContErr = runTime.deltaTValue()* + mag(contErr)().weightedAverage(mesh.V()).value(); + + scalar globalContErr = runTime.deltaTValue()* + contErr.weightedAverage(mesh.V()).value(); + cumulativeContErr += globalContErr; + + Info<< "time step continuity errors : sum local = " << sumLocalContErr + << ", global = " << globalContErr + << ", cumulative = " << cumulativeContErr + << endl; +} diff --git a/applications/solvers/rctfSpeciesTransport/createFields.H b/applications/solvers/rctfSpeciesTransport/createFields.H new file mode 100755 index 00000000..3b6ef0e5 --- /dev/null +++ b/applications/solvers/rctfSpeciesTransport/createFields.H @@ -0,0 +1,281 @@ +//creating the fields according to the recurrence dictionary + +IOdictionary recProperties_ +( + IOobject + ( + "recProperties0", + runTime.constant(), + mesh, + IOobject::MUST_READ, + IOobject::NO_WRITE + ) +); + +List fieldsDict_(recProperties_.lookup("fieldsPairs")); + +wordList fieldNames(fieldsDict_.size()); + +for(int i = 0; i < fieldsDict_.size(); i++) +{ + + fieldNames[i]= fieldsDict_[i][0]; +} + +Info<< "\n list of the fields: \n" << fieldNames << endl; + +//reading coherent velocity field name +label k = findIndex(fieldNames,"coh_velocity"); + +if (k < 0) +{ + FatalError <<"\n No field is defiened for the coherent velocity\n" << abort(FatalError); +} +const word Ucoh_pair = fieldsDict_[k][1]; + + volVectorField UcohRec + ( + IOobject + ( + Ucoh_pair, + runTime.timeName(), + mesh, + IOobject::READ_IF_PRESENT, + IOobject::AUTO_WRITE + ), + mesh, + dimensionedVector("zero",dimensionSet(0, 1, -1, 0, 0),vector::zero) + ); + +//reading incoherent velocity field name +k = findIndex(fieldNames,"inc_velocity"); + +if (k < 0) +{ + FatalError <<"\n No field is defiened for the incoherent velocity\n" << abort(FatalError); +} +const word Uinc_pair = fieldsDict_[k][1]; + + volVectorField UincRec + ( + IOobject + ( + Uinc_pair, + runTime.timeName(), + mesh, + IOobject::READ_IF_PRESENT, + IOobject::AUTO_WRITE + ), + mesh, + dimensionedVector("zero",dimensionSet(0, 1, -1, 0, 0),vector::zero) + ); + +//reading coherent turb kinetic energy field name +k = findIndex(fieldNames,"kSGS_coh"); +if (k < 0) +{ + FatalError <<"\n No field is defiened for the coherent subgrid-scale turbulent kinetic energy\n" << abort(FatalError); +} +const word kSGScoh_pair = fieldsDict_[k][1]; + + volScalarField kcohRec + ( + IOobject + ( + kSGScoh_pair, + runTime.timeName(), + mesh, + IOobject::READ_IF_PRESENT, + IOobject::AUTO_WRITE + ), + mesh, + dimensionedScalar("zero",dimensionSet(0, 2, -2, 0, 0),0.0) + ); + +//reading incoherent turb kinetic energy field name +k = findIndex(fieldNames,"kSGS_inc"); +if (k < 0) +{ + FatalError <<"\n No field is defiened for the coherent subgrid-scale turbulent kinetic energy\n" << abort(FatalError); +} +const word kSGSinc_pair = fieldsDict_[k][1]; + + volScalarField kincRec + ( + IOobject + ( + kSGSinc_pair, + runTime.timeName(), + mesh, + IOobject::READ_IF_PRESENT, + IOobject::AUTO_WRITE + ), + mesh, + dimensionedScalar("zero",dimensionSet(0, 2, -2, 0, 0),0.0) + ); + +// calculated fields +Info<< "\nCreating cell volume field\n" << endl; + +volScalarField delta +( + IOobject + ( + "delta", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh, + dimensionedScalar("delta", dimLength, 0.0) +); + +delta.primitiveFieldRef()=pow(mesh.V(),1.0/3.0); +delta.write(); + + +volVectorField URec +( + IOobject + ( + "URec", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh +); +Info<< "\nCreating turb kinetic energy field\n" << endl; + +volScalarField kRec +( + IOobject + ( + "kRec", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + kcohRec+kincRec +); + +// check if there is any negative values +forAll(kRec, cellI) +{ + if (kRec[cellI] < SMALL) + { + kRec[cellI] = 0.0; + } + +} + +const fvPatchList& patches = mesh.boundary(); +forAll(patches, patchI) +{ + kRec.boundaryFieldRef()[patchI] = 0.0; +} + + +kRec.write(); + +Info<< "\nCreating turb viscosity field\n" << endl; +volScalarField nutRec +( + IOobject + ( + "nutRec", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + sqrt(kRec)*delta*0.094 +); + +nutRec.write(); + + +Info<< "Calculating face flux field phiRec\n" << endl; +surfaceScalarField phiRec + ( + IOobject + ( + "phiRec", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + linearInterpolate(URec) & mesh.Sf() + ); + + phiRec.write(); + + singlePhaseTransportModel laminarTransport(URec, phiRec); + + autoPtr turbulence + ( + incompressible::turbulenceModel::New(URec, phiRec, laminarTransport) + ); + + dimensionedScalar Sc("Sc", dimless, laminarTransport); + dimensionedScalar Sct("Sct", dimless, laminarTransport); + +volScalarField alphat +( + IOobject + ( + "alphat", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + nutRec/Sct +); + + // create the scalar field + Info<< "Creating scalar transport field\n" << endl; + + volScalarField C + ( + IOobject + ( + "C", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh + ); + + fvScalarMatrix CEqn(C, dimless*dimVolume/(dimTime)); + + scalar relaxCoeff(0.0); + + Info<< "reading clockProperties\n" << endl; + + IOdictionary clockProperties + ( + IOobject + ( + "clockProperties", + mesh.time().constant(), + mesh, + IOobject::MUST_READ, + IOobject::NO_WRITE + ) + ); + + autoPtr myClock + ( + clockModel::New + ( + clockProperties, + mesh.time() + ) + ); diff --git a/applications/solvers/rctfSpeciesTransport/createRecBase.H b/applications/solvers/rctfSpeciesTransport/createRecBase.H new file mode 100644 index 00000000..79eb7be7 --- /dev/null +++ b/applications/solvers/rctfSpeciesTransport/createRecBase.H @@ -0,0 +1,66 @@ +// check which recProperties dicts are present, read them in and construct a PtrList of recBases +// names for dicts can be "recProperties" or "recPropertiesN" where N in {0, 1, ...} + +#include "error.H" + +word dictName = "recProperties"; +wordList recPropertiesList(0); +PtrList recBases(0); +label maxDictNumber = 100; + +{ + IOdictionary recPropDict + ( + IOobject + ( + dictName, + mesh.time().constant(), + mesh, + IOobject::READ_IF_PRESENT, + IOobject::NO_WRITE + ) + ); + + if (recPropDict.headerOk()) + { + recPropertiesList.append(dictName); + } +} + + +for (label counter = 0; counter < maxDictNumber; counter++) +{ + word dictNameIter = dictName + Foam::name(counter); + IOdictionary recPropDict + ( + IOobject + ( + dictNameIter, + mesh.time().constant(), + mesh, + IOobject::READ_IF_PRESENT, + IOobject::NO_WRITE + ) + ); + + if (recPropDict.headerOk()) + { + recPropertiesList.append(dictNameIter); + } +} + +if (recPropertiesList.size() == 0) +{ + FatalError << "no recProperties dicts found" << endl; +} + +else +{ + Info << "found " << recPropertiesList.size() << " dicts with names " << recPropertiesList << endl; +} + + +for (label counter = 0; counter < recPropertiesList.size(); counter++) +{ + recBases.append( new recBase(mesh, recPropertiesList[counter])); +} diff --git a/applications/solvers/rctfSpeciesTransport/rctfSpeciesTransport.C b/applications/solvers/rctfSpeciesTransport/rctfSpeciesTransport.C new file mode 100755 index 00000000..061a0fb2 --- /dev/null +++ b/applications/solvers/rctfSpeciesTransport/rctfSpeciesTransport.C @@ -0,0 +1,131 @@ +/*---------------------------------------------------------------------------*\ + CFDEMcoupling academic - Open Source CFD-DEM coupling + + Contributing authors: + Thomas Lichtenegger, Gerhard Holzinger, Sanaz Abbasi + Copyright (C) 2015- Johannes Kepler University, Linz +------------------------------------------------------------------------------- +License + This file is part of CFDEMcoupling academic. + + CFDEMcoupling academic 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. + + CFDEMcoupling academic 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 CFDEMcoupling academic. If not, see . + +Application + Turbulent Transport Recurrence Solver for modal decomposition + +Description + Solves a transport equation for a passive scalar on a single-phase solution + for a solver based on recurrence statistics + +Rules + Solution data to compute the recurrence statistics from, needs to + reside in $CASE_ROOT/dataBase(0...N) + Time step data in the first dataBase needs to be evenly spaced in time + A list of indices for the corresponding incoherent fields to coherent ones + should be provided. + +\*---------------------------------------------------------------------------*/ + +#include "fvCFD.H" +#include "singlePhaseTransportModel.H" +#include "turbulentTransportModel.H" +#include "fvOptions.H" + +#include "recBase.H" +#include "recModel.H" + +#include "clockModel.H" + +#include "objectRegistry.H" +#include "VectorSpace.H" +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +int main(int argc, char *argv[]) +{ + #include "postProcess.H" + #include "setRootCase.H" + #include "createTime.H" + #include "createMesh.H" + #include "createControl.H" + #include "createFields.H" + #include "createFvOptions.H" + scalar cumulativeContErr = 0; + + //create recBases according to a list of recProperties + #include "createRecBase.H" + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + + + label recTimeIndex(0); + label currTimeIndex(0); + + scalar recTimeStep_=recBases[0].recM().recTimeStep(); + labelPairList incPairTimeIndex_(0); + + IFstream pairFile("incIndexPairList"); + pairFile >> incPairTimeIndex_; + + while (runTime.run()) + { + + myClock().start(1,"Global"); + runTime++; + + myClock().start(11,"Total"); + + Info<< "Time = " << runTime.timeName() << nl << endl; + + myClock().start(2,"fieldUpdate"); + + if ( runTime.timeOutputValue() - (recTimeIndex+1)*recTimeStep_ + 1.0e-5 > 0.0 ) + { + Info << "Updating fields at run time " << runTime.timeOutputValue() + << " corresponding to recurrence time " << (recTimeIndex+1)*recTimeStep_ << ".\n" << endl; + recBases[0].updateRecFields(); + #include "readFields.H" + + recTimeIndex++; + } + + myClock().stop("fieldUpdate"); + + #include "continuityErrCalc.H" + + myClock().start(3,"speciesEqn"); + #include "CEq.H" + myClock().stop("speciesEqn"); + + myClock().stop("Total"); + + runTime.write(); + + Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s" + << " ClockTime = " << runTime.elapsedClockTime() << " s" + << nl << endl; + + myClock().stop("Global"); + + } + + + myClock().evalPar(); + myClock().normHist(); + + Info<< "End\n" << endl; + + return 0; +} + + +// ************************************************************************* // diff --git a/applications/solvers/rctfSpeciesTransport/readFields.H b/applications/solvers/rctfSpeciesTransport/readFields.H new file mode 100755 index 00000000..a5ee930e --- /dev/null +++ b/applications/solvers/rctfSpeciesTransport/readFields.H @@ -0,0 +1,40 @@ + +currTimeIndex = recBases[0].recM().currentTimeIndex(); + +Info << "current Time Index = " << currTimeIndex << endl; + +recBases[0].recM().exportVolVectorField(Ucoh_pair,UcohRec); +recBases[0].recM().exportVolScalarField(kSGScoh_pair,kcohRec); + + +label incTimeIndex = incPairTimeIndex_[currTimeIndex][1]; + +Info << " incoherent pair Time Index = " << incTimeIndex << endl; + +UincRec = recBases[1].recM().exportVolVectorField(Uinc_pair,incTimeIndex); +kincRec = recBases[1].recM().exportVolScalarField(kSGSinc_pair,incTimeIndex); + +kRec = kcohRec+kincRec; + +forAll(kRec, cellI) +{ + if (kRec[cellI] < SMALL) + { + kRec[cellI] = 0.0; + } +} + +const fvPatchList& patches = mesh.boundary(); +forAll(patches, patchI) +{ + kRec.boundaryFieldRef()[patchI] = 0.0; +} + +URec = UcohRec + UincRec; +phiRec = linearInterpolate(URec) & mesh.Sf(); + +nutRec = sqrt(kRec)*delta*0.094; + +alphat = nutRec/Sct; +alphat.correctBoundaryConditions(); + diff --git a/etc/solver-list.txt b/etc/solver-list.txt index 8391f45e..5eeca1d2 100644 --- a/etc/solver-list.txt +++ b/etc/solver-list.txt @@ -12,3 +12,4 @@ cfdemSolverIB/dir cfdemSolverPisoScalar/dir cfdemSolverRhoPimpleChem/dir cfdemSolverMultiphase/dir +rctfSpeciesTransport/dir