diff --git a/applications/solvers/cfdemSolverPimple/UEqn.H b/applications/solvers/cfdemSolverPimple/UEqn.H index d72c8c23..f449301e 100644 --- a/applications/solvers/cfdemSolverPimple/UEqn.H +++ b/applications/solvers/cfdemSolverPimple/UEqn.H @@ -4,7 +4,7 @@ tmp tUEqn ( fvm::ddt(voidfraction,U) - fvm::Sp(fvc::ddt(voidfraction),U) + fvm::div(phi,U) - fvm::Sp(fvc::div(phi),U) - + particleCloud.divVoidfractionTau(U, voidfraction) + + particleCloud.divVoidfractionTau(U, voidfraction) // in case of doing "periodic box" simulations viscous term can be commented and the solver can be compiled // again. However, the effect of this term on the results is negligible (about 1-2%). - fOther/rho == fvOptions(U) diff --git a/applications/solvers/cfdemSolverPimple/cfdemSolverPimple.C b/applications/solvers/cfdemSolverPimple/cfdemSolverPimple.C index 1b8199f3..c77f896d 100644 --- a/applications/solvers/cfdemSolverPimple/cfdemSolverPimple.C +++ b/applications/solvers/cfdemSolverPimple/cfdemSolverPimple.C @@ -58,6 +58,12 @@ int main(int argc, char *argv[]) cfdemCloud particleCloud(mesh); #include "checkModelType.H" + // switch for periodic box simulations + Switch periodicBoxSwitch + ( + pimple.dict().lookupOrDefault("periodicBoxSwitch", false) + ); + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // Info<< "\nStarting time loop\n" << endl; while (runTime.loop()) @@ -117,6 +123,11 @@ int main(int argc, char *argv[]) Info << "skipping flow solution." << endl; } + if (periodicBoxSwitch) + { + #include "periodicBoxProperties.H" + } + runTime.write(); Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s" diff --git a/applications/solvers/cfdemSolverPimple/createFields.H b/applications/solvers/cfdemSolverPimple/createFields.H index 4d3fd024..d68c901a 100644 --- a/applications/solvers/cfdemSolverPimple/createFields.H +++ b/applications/solvers/cfdemSolverPimple/createFields.H @@ -139,6 +139,30 @@ linearInterpolate(U) & mesh.Sf() ); +//Periodic box + volScalarField unity + ( + IOobject + ( + "unity", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh, + dimensionedScalar("unity", dimless, 1.0) + ); + +// gravity vector for periodic boundary conditions + dimensionedVector gN + ( + "gN", + dimensionSet(0,1,-2,0,0,0,0), + vector(0,0,-9.81) + ); + + label pRefCell = 0; scalar pRefValue = 0.0; diff --git a/applications/solvers/cfdemSolverPimple/periodicBoxProperties.H b/applications/solvers/cfdemSolverPimple/periodicBoxProperties.H new file mode 100644 index 00000000..8ef5fecb --- /dev/null +++ b/applications/solvers/cfdemSolverPimple/periodicBoxProperties.H @@ -0,0 +1,65 @@ +// gravity direction should always be -z! +dimensionedScalar volume = fvc::domainIntegrate(unity); + +Info<< "particle_ENSTROPHY: " + << ( + fvc::domainIntegrate( 0.5*magSqr(fvc::curl(Us))) + /volume + ).value() + << endl; + +Info<< "air_ENSTROPHY: " + << ( + fvc::domainIntegrate( 0.5*magSqr(fvc::curl(U))) + /volume + ).value() + << endl; + +Info<< "slip_velocity: " + << - (( + fvc::domainIntegrate(voidfraction*(U&gN)).value() + /fvc::domainIntegrate(voidfraction*mag(gN)).value() + ) + - ( + fvc::domainIntegrate((1.0-voidfraction)*(Us&gN)).value() + /fvc::domainIntegrate((1.0-voidfraction)*mag(gN)).value() + )) + << endl; + +dimensionedVector alpha1Us = fvc::domainIntegrate((1.0-voidfraction)*(Us))/fvc::domainIntegrate((1.0-voidfraction)); +dimensionedVector alpha2U = fvc::domainIntegrate(voidfraction*(U))/fvc::domainIntegrate(voidfraction); +dimensionedScalar alpha1M = fvc::domainIntegrate((1.0-voidfraction))/volume; +dimensionedScalar alpha2M = scalar(1.0) - alpha1M; + +Info<< "TKE gas: " + << 0.5 + *( + fvc::domainIntegrate(voidfraction*(U&U)).value() + /fvc::domainIntegrate(voidfraction).value() + ) + - 0.5 + *( + alpha2U.value() + &alpha2U.value() + ) + << endl; + +Info<< "TKE solid: " + << 0.5 + *( + fvc::domainIntegrate((1.0-voidfraction)*(Us&Us)).value() + /fvc::domainIntegrate(1.0-voidfraction).value() + ) + - 0.5 + *( + alpha1Us.value() + &alpha1Us.value() + ) + << endl; + +Info<< "PhiP2: " + << fvc::domainIntegrate((1.0-voidfraction)*(1.0-voidfraction)).value() + /fvc::domainIntegrate(unity).value() + - alpha1M.value()*alpha1M.value() + << endl; +