zhyang-dev/CFDEMcoupling-PFM
1
0
Fork 0
mirror of https://github.com/ParticulateFlow/CFDEMcoupling-PFM.git synced 2025-12-08 06:37:44 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

Updating solvers to OF 4.x.

Browse Source
This commit is contained in:
Thomas Lichtenegger
2016-07-21 12:16:55 +02:00
parent 0a2f1c793a
commit d30f828fe5
7 changed files with 189 additions and 332 deletions
Download Patch File Download Diff File Expand all files Collapse all files

71
applications/solvers/cfdemSolverPiso/cfdemSolverPiso.C
View File

@ -38,6 +38,7 @@ Description
#include "singlePhaseTransportModel.H"
#include "turbulentTransportModel.H"
#include "pisoControl.H"
#include "fvOptions.H"
#include "cfdemCloud.H"
#include "implicitCouple.H"
@ -54,6 +55,7 @@ int main(int argc, char *argv[])
#include "createMesh.H"
#include "createControl.H"
#include "createFields.H"
#include "createFvOptions.H"
#include "initContinuityErrs.H"
// create cfdemCloud
@ -101,78 +103,19 @@ int main(int argc, char *argv[])
// Pressure-velocity PISO corrector
{
// Momentum predictor
fvVectorMatrix UEqn
(
fvm::ddt(voidfraction,U) - fvm::Sp(fvc::ddt(voidfraction),U)
+ fvm::div(phi,U) - fvm::Sp(fvc::div(phi),U)
// + turbulence->divDevReff(U)
+ particleCloud.divVoidfractionTau(U, voidfraction)
==
- fvm::Sp(Ksl/rho,U)
);
UEqn.relax();
if (piso.momentumPredictor() && (modelType=="B" || modelType=="Bfull"))
solve(UEqn == - fvc::grad(p) + Ksl/rho*Us);
else if (piso.momentumPredictor())
solve(UEqn == - voidfraction*fvc::grad(p) + Ksl/rho*Us);
#include "UEqn.H"
// --- PISO loop
while (piso.correct())
{
volScalarField rUA = 1.0/UEqn.A();
surfaceScalarField rUAf("(1|A(U))", fvc::interpolate(rUA));
volScalarField rUAvoidfraction("(voidfraction2|A(U))",rUA*voidfraction);
surfaceScalarField rUAfvoidfraction("(voidfraction2|A(U)F)", fvc::interpolate(rUAvoidfraction));
U = rUA*UEqn.H();
phi = ( fvc::interpolate(U*voidfraction) & mesh.Sf() )
+ rUAfvoidfraction*fvc::ddtCorr(U, phi);
surfaceScalarField phiS(fvc::interpolate(Us*voidfraction) & mesh.Sf());
surfaceScalarField phiGes = phi + rUAf*(fvc::interpolate(Ksl/rho) * phiS);
if (modelType=="A")
rUAvoidfraction = volScalarField("(voidfraction2|A(U))",rUA*voidfraction*voidfraction);
// Non-orthogonal pressure corrector loop
while (piso.correctNonOrthogonal())
{
// Pressure corrector
fvScalarMatrix pEqn
(
fvm::laplacian(rUAvoidfraction, p) == fvc::div(phiGes) + particleCloud.ddtVoidfraction()
);
pEqn.setReference(pRefCell, pRefValue);
pEqn.solve(mesh.solver(p.select(piso.finalInnerIter())));
if (piso.finalNonOrthogonalIter())
{
phiGes -= pEqn.flux();
phi = phiGes;
}
} // end non-orthogonal corrector loop
#include "continuityErrorPhiPU.H"
if (modelType=="B" || modelType=="Bfull")
U -= rUA*fvc::grad(p) - Ksl/rho*Us*rUA;
else
U -= voidfraction*rUA*fvc::grad(p) - Ksl/rho*Us*rUA;
U.correctBoundaryConditions();
} // end piso loop
#include "pEqn.H"
}
}
laminarTransport.correct();
turbulence->correct();
}// end solveFlow
}
else
{
Info << "skipping flow solution." << endl;

2
applications/solvers/cfdemSolverPiso/createFields.H
View File

@ -122,3 +122,5 @@ surfaceScalarField phi
(
incompressible::turbulenceModel::New(U, phi, laminarTransport)
);
#include "createMRF.H"

114
applications/solvers/cfdemSolverPisoMS/cfdemSolverPisoMS.C
View File

@ -25,12 +25,12 @@ License
along with CFDEMcoupling. If not, see <http://www.gnu.org/licenses/>.
Application
cfdemSolverPisoMS
cfdemSolverPiso
Description
Transient solver for incompressible flow.
Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected.
The code is an evolution of the solver pisoFoam in OpenFOAM(R) 1.6,
The code is an evolution of the solver pisoFoam in OpenFOAM(R) 1.6,
where additional functionality for CFD-DEM coupling is added.
\*---------------------------------------------------------------------------*/
@ -38,11 +38,13 @@ Description
#include "singlePhaseTransportModel.H"
#include "turbulentTransportModel.H"
#include "pisoControl.H"
#include "fvOptions.H"
#include "cfdemCloudMS.H"
#include "implicitCouple.H"
#include "clockModel.H"
#include "smoothingModel.H"
#include "forceModel.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -53,6 +55,7 @@ int main(int argc, char *argv[])
#include "createMesh.H"
#include "createControl.H"
#include "createFields.H"
#include "createFvOptions.H"
#include "initContinuityErrs.H"
// create cfdemCloud
@ -61,11 +64,10 @@ int main(int argc, char *argv[])
#include "checkModelType.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl;
while (runTime.loop())
{
Info<< "\nStarting time loop\n" << endl;
particleCloud.clockM().start(1,"Global");
particleCloud.clockM().start(1,"Global");
Info<< "Time = " << runTime.timeName() << nl << endl;
@ -73,91 +75,51 @@ int main(int argc, char *argv[])
// do particle stuff
particleCloud.clockM().start(2,"Coupling");
particleCloud.evolve(voidfraction,Us,U);
bool hasEvolved = particleCloud.evolve(voidfraction,Us,U);
if(hasEvolved)
{
particleCloud.smoothingM().smoothen(particleCloud.forceM(0).impParticleForces());
}
Info << "update Ksl.internalField()" << endl;
Ksl = particleCloud.momCoupleM(0).impMomSource();
particleCloud.smoothingM().smoothen(Ksl);
Ksl.correctBoundaryConditions();
//Force Checks
vector fTotal(0,0,0);
vector fImpTotal = sum(mesh.V()*Ksl.internalField()*(Us.internalField()-U.internalField())).value();
reduce(fImpTotal, sumOp<vector>());
Info << "TotalForceExp: " << fTotal << endl;
Info << "TotalForceImp: " << fImpTotal << endl;
#include "solverDebugInfo.H"
particleCloud.clockM().stop("Coupling");
particleCloud.clockM().start(26,"Flow");
// Pressure-velocity PISO corrector
if(particleCloud.solveFlow())
{
// Momentum predictor
fvVectorMatrix UEqn
(
fvm::ddt(voidfraction,U) //particleCloud.ddtVoidfractionU(U,voidfraction) //
+ fvm::div(phi, U)
// + turbulence->divDevReff(U)
+ particleCloud.divVoidfractionTau(U, voidfraction)
==
- fvm::Sp(Ksl/rho,U)
);
if (modelType=="B")
UEqn == - fvc::grad(p) + Ksl/rho*Us;
else
UEqn == - voidfraction*fvc::grad(p) + Ksl/rho*Us;
UEqn.relax();
if (piso.momentumPredictor())
solve(UEqn);
// --- PISO loop
while (piso.correct())
// Pressure-velocity PISO corrector
{
volScalarField rUA = 1.0/UEqn.A();
// Momentum predictor
#include "UEqn.H"
surfaceScalarField rUAf("(1|A(U))", fvc::interpolate(rUA));
volScalarField rUAvoidfraction("(voidfraction2|A(U))",rUA*voidfraction);
// --- PISO loop
U = rUA*UEqn.H();
phi = (fvc::interpolate(U*voidfraction) & mesh.Sf() );
//+ fvc::ddtPhiCorr(rUAvoidfraction, U, phi);
surfaceScalarField phiS(fvc::interpolate(Us*voidfraction) & mesh.Sf());
surfaceScalarField phiGes = phi + rUAf*(fvc::interpolate(Ksl/rho) * phiS);
if (modelType=="A")
rUAvoidfraction = volScalarField("(voidfraction2|A(U))",rUA*voidfraction*voidfraction);
// Non-orthogonal pressure corrector loop
while (piso.correctNonOrthogonal())
while (piso.correct())
{
// Pressure corrector
fvScalarMatrix pEqn
(
fvm::laplacian(rUAvoidfraction, p) == fvc::div(phiGes) + particleCloud.ddtVoidfraction()
);
pEqn.setReference(pRefCell, pRefValue);
#include "pEqn.H"
}
}
pEqn.solve(mesh.solver(p.select(piso.finalInnerIter())));
if (piso.finalNonOrthogonalIter())
{
phiGes -= pEqn.flux();
}
} // end non-orthogonal corrector loop
#include "continuityErrorPhiPU.H"
if (modelType=="B")
U -= rUA*fvc::grad(p) - Ksl/rho*Us*rUA;
else
U -= voidfraction*rUA*fvc::grad(p) - Ksl/rho*Us*rUA;
U.correctBoundaryConditions();
} // end piso loop
laminarTransport.correct();
turbulence->correct();
}
else
{
Info << "skipping flow solution." << endl;
}
turbulence->correct();
runTime.write();
@ -170,7 +132,7 @@ int main(int argc, char *argv[])
}
Info<< "End\n" << endl;
return 0;
}

1
applications/solvers/cfdemSolverPisoScalar/Make/options
View File

@ -7,6 +7,7 @@ EXE_INC = \
-I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
-I$(LIB_SRC)/transportModels \
-I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
-I ../cfdemSolverPiso \
-I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/lnInclude \
-I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/cfdTools \

103
applications/solvers/cfdemSolverPisoScalar/cfdemSolverPisoScalar.C
View File

@ -25,7 +25,7 @@ License
along with CFDEMcoupling. If not, see <http://www.gnu.org/licenses/>.
Application
cfdemSolverPisoScalar
cfdemSolverPiso
Description
Transient solver for incompressible flow.
@ -38,9 +38,11 @@ Description
#include "singlePhaseTransportModel.H"
#include "turbulentTransportModel.H"
#include "pisoControl.H"
#include "fvOptions.H"
#include "cfdemCloud.H"
#include "implicitCouple.H"
#include "clockModel.H"
#include "smoothingModel.H"
#include "forceModel.H"
@ -53,6 +55,7 @@ int main(int argc, char *argv[])
#include "createMesh.H"
#include "createControl.H"
#include "createFields.H"
#include "createFvOptions.H"
#include "initContinuityErrs.H"
// create cfdemCloud
@ -64,11 +67,14 @@ int main(int argc, char *argv[])
Info<< "\nStarting time loop\n" << endl;
while (runTime.loop())
{
particleCloud.clockM().start(1,"Global");
Info<< "Time = " << runTime.timeName() << nl << endl;
#include "CourantNo.H"
// do particle stuff
particleCloud.clockM().start(2,"Coupling");
bool hasEvolved = particleCloud.evolve(voidfraction,Us,U);
if(hasEvolved)
@ -80,100 +86,38 @@ int main(int argc, char *argv[])
Ksl = particleCloud.momCoupleM(0).impMomSource();
Ksl.correctBoundaryConditions();
//Force Checks
vector fTotal(0,0,0);
vector fImpTotal = sum(mesh.V()*Ksl.internalField()*(Us.internalField()-U.internalField())).value();
reduce(fImpTotal, sumOp<vector>());
Info << "TotalForceExp: " << fTotal << endl;
Info << "TotalForceImp: " << fImpTotal << endl;
#include "solverDebugInfo.H"
particleCloud.clockM().stop("Coupling");
// get scalar source from DEM
particleCloud.forceM(1).manipulateScalarField(Tsource);
Tsource.correctBoundaryConditions();
// solve scalar transport equation
fvScalarMatrix TEqn
(
fvm::ddt(voidfraction,T) - fvm::Sp(fvc::ddt(voidfraction),T)
+ fvm::div(phi, T) - fvm::Sp(fvc::div(phi),T)
- fvm::laplacian(DT*voidfraction, T)
==
Tsource
);
TEqn.relax();
TEqn.solve();
particleCloud.clockM().start(26,"Flow");
#include "TEqn.H"
if(particleCloud.solveFlow())
{
// Pressure-velocity PISO corrector
{
// Momentum predictor
fvVectorMatrix UEqn
(
fvm::ddt(voidfraction,U) - fvm::Sp(fvc::ddt(voidfraction),U)
+ fvm::div(phi,U) - fvm::Sp(fvc::div(phi),U)
// + turbulence->divDevReff(U)
+ particleCloud.divVoidfractionTau(U, voidfraction)
==
- fvm::Sp(Ksl/rho,U)
);
UEqn.relax();
if (piso.momentumPredictor() && (modelType=="B" || modelType=="Bfull"))
solve(UEqn == - fvc::grad(p) + Ksl/rho*Us);
else if (piso.momentumPredictor())
solve(UEqn == - voidfraction*fvc::grad(p) + Ksl/rho*Us);
#include "UEqn.H"
// --- PISO loop
while (piso.correct())
{
volScalarField rUA = 1.0/UEqn.A();
surfaceScalarField rUAf("(1|A(U))", fvc::interpolate(rUA));
volScalarField rUAvoidfraction("(voidfraction2|A(U))",rUA*voidfraction);
surfaceScalarField rUAfvoidfraction("(voidfraction2|A(U)F)", fvc::interpolate(rUAvoidfraction));
U = rUA*UEqn.H();
phi = ( fvc::interpolate(U*voidfraction) & mesh.Sf() )
+ rUAfvoidfraction*fvc::ddtCorr(U, phi);
surfaceScalarField phiS(fvc::interpolate(Us*voidfraction) & mesh.Sf());
surfaceScalarField phiGes = phi + rUAf*(fvc::interpolate(Ksl/rho) * phiS);
if (modelType=="A")
rUAvoidfraction = volScalarField("(voidfraction2|A(U))",rUA*voidfraction*voidfraction);
while (piso.correctNonOrthogonal())
{
// Pressure corrector
fvScalarMatrix pEqn
(
fvm::laplacian(rUAvoidfraction, p) == fvc::div(phiGes) + particleCloud.ddtVoidfraction()
);
pEqn.setReference(pRefCell, pRefValue);
pEqn.solve(mesh.solver(p.select(piso.finalInnerIter())));
if (piso.finalNonOrthogonalIter())
{
phiGes -= pEqn.flux();
phi = phiGes;
}
} // end non-orthogonal corrector loop
#include "continuityErrorPhiPU.H"
if (modelType=="B" || modelType=="Bfull")
U -= rUA*fvc::grad(p) - Ksl/rho*Us*rUA;
else
U -= voidfraction*rUA*fvc::grad(p) - Ksl/rho*Us*rUA;
U.correctBoundaryConditions();
} // end piso loop
#include "pEqn.H"
}
}
laminarTransport.correct();
turbulence->correct();
}// end solveFlow
}
else
{
Info << "skipping flow solution." << endl;
@ -184,6 +128,9 @@ int main(int argc, char *argv[])
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
particleCloud.clockM().stop("Flow");
particleCloud.clockM().stop("Global");
}
Info<< "End\n" << endl;

228
applications/solvers/cfdemSolverPisoScalar/createFields.H
View File

@ -1,36 +1,36 @@
Info<< "Reading field p\n" << endl;
volScalarField p
(
IOobject
(
"p",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
Info<< "Reading physical velocity field U" << endl;
Info<< "Note: only if voidfraction at boundary is 1, U is superficial velocity!!!\n" << endl;
volVectorField U
(
IOobject
(
"U",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
//========================
// drag law modelling
//========================
Info<< "Reading field p\n" << endl;
volScalarField p
(
IOobject
(
"p",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
Info<< "Reading physical velocity field U" << endl;
Info<< "Note: only if voidfraction at boundary is 1, U is superficial velocity!!!\n" << endl;
volVectorField U
(
IOobject
(
"U",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
//========================
// drag law modelling
//========================
Info<< "\nReading momentum exchange field Ksl\n" << endl;
volScalarField Ksl
(
@ -44,8 +44,8 @@
),
mesh
//dimensionedScalar("0", dimensionSet(0, 0, -1, 0, 0), 1.0)
);
);
Info<< "\nReading voidfraction field voidfraction = (Vgas/Vparticle)\n" << endl;
volScalarField voidfraction
(
@ -58,8 +58,8 @@
IOobject::AUTO_WRITE
),
mesh
);
);
Info<< "\nCreating density field rho\n" << endl;
volScalarField rho
(
@ -71,27 +71,27 @@
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
mesh,
mesh,
dimensionedScalar("0", dimensionSet(1, -3, 0, 0, 0), 1.0)
);
Info<< "Reading particle velocity field Us\n" << endl;
volVectorField Us
(
IOobject
(
"Us",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
//========================
// scalar field modelling
//========================
);
Info<< "Reading particle velocity field Us\n" << endl;
volVectorField Us
(
IOobject
(
"Us",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
//========================
// scalar field modelling
//========================
Info<< "\nCreating dummy density field rho = 1\n" << endl;
volScalarField T
(
@ -103,10 +103,10 @@
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh//,
mesh//,
//dimensionedScalar("0", dimensionSet(0, 0, -1, 1, 0), 273.15)
);
);
Info<< "\nCreating fluid-particle heat flux field\n" << endl;
volScalarField Tsource
(
@ -118,57 +118,59 @@
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh//,
mesh//,
//dimensionedScalar("0", dimensionSet(0, 0, -1, 1, 0), 0.0)
);
IOdictionary transportProperties
(
IOobject
(
"transportProperties",
runTime.constant(),
mesh,
IOobject::MUST_READ,
IOobject::NO_WRITE
)
);
dimensionedScalar DT
(
transportProperties.lookup("DT")
);
//========================
//# include "createPhi.H"
#ifndef createPhi_H
#define createPhi_H
Info<< "Reading/calculating face flux field phi\n" << endl;
surfaceScalarField phi
(
IOobject
(
"phi",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
linearInterpolate(U*voidfraction) & mesh.Sf()
);
#endif
label pRefCell = 0;
scalar pRefValue = 0.0;
setRefCell(p, mesh.solutionDict().subDict("PISO"), pRefCell, pRefValue);
singlePhaseTransportModel laminarTransport(U, phi);
autoPtr<incompressible::turbulenceModel> turbulence
(
incompressible::turbulenceModel::New(U, phi, laminarTransport)
);
);
IOdictionary transportProperties
(
IOobject
(
"transportProperties",
runTime.constant(),
mesh,
IOobject::MUST_READ,
IOobject::NO_WRITE
)
);
dimensionedScalar DT
(
transportProperties.lookup("DT")
);
//========================
//# include "createPhi.H"
#ifndef createPhi_H
#define createPhi_H
Info<< "Reading/calculating face flux field phi\n" << endl;
surfaceScalarField phi
(
IOobject
(
"phi",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
linearInterpolate(U*voidfraction) & mesh.Sf()
);
#endif
label pRefCell = 0;
scalar pRefValue = 0.0;
setRefCell(p, mesh.solutionDict().subDict("PISO"), pRefCell, pRefValue);
singlePhaseTransportModel laminarTransport(U, phi);
autoPtr<incompressible::turbulenceModel> turbulence
(
incompressible::turbulenceModel::New(U, phi, laminarTransport)
);
#include "createMRF.H"

2
src/lagrangian/cfdemParticle/cfdTools/continuityErrorPhiPU.H
View File

@ -34,7 +34,7 @@ Description
\*---------------------------------------------------------------------------*/
{
volScalarField contErr( fvc::div(phiGes) + fvc::ddt(voidfraction) );
volScalarField contErr( fvc::div(phi) + fvc::ddt(voidfraction) );
scalar sumLocalContErr = runTime.deltaTValue()*
mag(contErr)().weightedAverage(mesh.V()).value();