diff --git a/Allwmake b/Allwmake index 749c60fc8..3dcdff3e7 100755 --- a/Allwmake +++ b/Allwmake @@ -1,19 +1,27 @@ #!/bin/sh -set -x +cd ${0%/*} || exit 1 # run from this directory -# run from this directory only -cd ${0%/*} || exit 1 +if [ "$PWD" != "$WM_PROJECT_DIR" ] +then + echo "Error: Current directory is not \$WM_PROJECT_DIR" + echo " The environment variable are not consistent with the installation." + echo " Check the OpenFOAM entries in your dot-files and source them." + exit 1 +fi # wmake is required for subsequent targets -(cd wmake/src && make) +( cd wmake/src && make ) -(cd $WM_THIRD_PARTY_DIR && ./Allwmake) +# build ThirdParty sources +( cd $WM_THIRD_PARTY_DIR && ./Allwmake ) -(cd src && ./Allwmake) - -(cd applications && ./Allwmake) +# build OpenFOAM libraries and applications +src/Allwmake +applications/Allwmake if [ "$1" = doc ] then - (cd doc && ./Allwmake) + doc/Allwmake fi + +# ----------------------------------------------------------------- end-of-file diff --git a/applications/Allwmake b/applications/Allwmake index 7437e4f9b..82a2ec0df 100755 --- a/applications/Allwmake +++ b/applications/Allwmake @@ -1,5 +1,17 @@ #!/bin/sh +cd ${0%/*} || exit 1 # run from this directory + +if [ "$PWD" != "$WM_PROJECT_DIR/applications" ] +then + echo "Error: Current directory in not \$WM_PROJECT_DIR/applications" + echo " The environment variable are not consistent with the installation." + echo " Check the OpenFOAM entries in your dot-files and source them." + exit 1 +fi + set -x -( cd solvers && wmake all ) -( cd utilities && wmake all ) +wmake all solvers +wmake all utilities + +# ----------------------------------------------------------------- end-of-file diff --git a/applications/solvers/combustion/dieselEngineFoam/dieselEngineFoam.C b/applications/solvers/combustion/dieselEngineFoam/dieselEngineFoam.C index bb76f0ce9..77f897135 100644 --- a/applications/solvers/combustion/dieselEngineFoam/dieselEngineFoam.C +++ b/applications/solvers/combustion/dieselEngineFoam/dieselEngineFoam.C @@ -77,12 +77,6 @@ int main(int argc, char *argv[]) Info<< "Crank angle = " << runTime.theta() << " CA-deg" << endl; mesh.move(); - Check that this is unnecessary. HJ - // 1.6.x merge. HJ, 26/Aug/2010 - const_cast - ( - volPointInterpolation::New(mesh) - ).updateMesh(); dieselSpray.evolve(); diff --git a/applications/solvers/engine/sonicTurbDyMEngineFoam/Make/options b/applications/solvers/engine/sonicTurbDyMEngineFoam/Make/options index 1383b9ac5..d463180c5 100644 --- a/applications/solvers/engine/sonicTurbDyMEngineFoam/Make/options +++ b/applications/solvers/engine/sonicTurbDyMEngineFoam/Make/options @@ -4,7 +4,7 @@ EXE_INC = \ -I$(LIB_SRC)/dynamicMesh/lnInclude \ -I$(LIB_SRC)/dynamicFvMesh/lnInclude \ -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \ - -I$(LIB_SRC)/turbulenceModels/RAS + -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel EXE_LIBS = \ -lengine \ @@ -12,6 +12,7 @@ EXE_LIBS = \ -ldynamicMesh \ -ltopoChangerFvMesh \ -lcompressibleRASModels \ + -lcompressibleLESModels \ -lbasicThermophysicalModels \ -lspecie \ -lmeshTools \ diff --git a/applications/solvers/engine/sonicTurbDyMEngineFoam/createFields.H b/applications/solvers/engine/sonicTurbDyMEngineFoam/createFields.H index b87373e87..a55b6a621 100644 --- a/applications/solvers/engine/sonicTurbDyMEngineFoam/createFields.H +++ b/applications/solvers/engine/sonicTurbDyMEngineFoam/createFields.H @@ -2,10 +2,11 @@ Info<< "Reading thermophysical properties\n" << endl; - autoPtr thermo + autoPtr pThermo ( - basicThermo::New(mesh) + basicPsiThermo::New(mesh) ); + basicPsiThermo& thermo = pThermo(); volScalarField rho ( @@ -17,16 +18,16 @@ IOobject::NO_READ, IOobject::AUTO_WRITE ), - thermo->rho() + thermo.rho() ); rho.oldTime(); - volScalarField& p = thermo->p(); + volScalarField& p = thermo.p(); p.oldTime(); - - const volScalarField& psi = thermo->psi(); - const volScalarField& T = thermo->T(); - volScalarField& h = thermo->h(); + + const volScalarField& psi = thermo.psi(); + const volScalarField& T = thermo.T(); + volScalarField& h = thermo.h(); Info<< "\nReading field U\n" << endl; @@ -47,14 +48,14 @@ Info<< "Creating turbulence model\n" << endl; - autoPtr turbulence + autoPtr turbulence ( - compressible::RASModel::New + compressible::turbulenceModel::New ( rho, U, phi, - thermo() + thermo ) ); diff --git a/applications/solvers/engine/sonicTurbDyMEngineFoam/hEqn.H b/applications/solvers/engine/sonicTurbDyMEngineFoam/hEqn.H index 4d65bbc46..530b4ef30 100644 --- a/applications/solvers/engine/sonicTurbDyMEngineFoam/hEqn.H +++ b/applications/solvers/engine/sonicTurbDyMEngineFoam/hEqn.H @@ -11,5 +11,5 @@ hEqn.solve(); - thermo->correct(); + thermo.correct(); } diff --git a/applications/solvers/engine/sonicTurbDyMEngineFoam/pEqn.H b/applications/solvers/engine/sonicTurbDyMEngineFoam/pEqn.H index 623caef9c..7ea343d36 100644 --- a/applications/solvers/engine/sonicTurbDyMEngineFoam/pEqn.H +++ b/applications/solvers/engine/sonicTurbDyMEngineFoam/pEqn.H @@ -1,12 +1,12 @@ { - rho = thermo->rho(); + rho = thermo.rho(); rUA = 1.0/UEqn.A(); H = UEqn.H(); U = rUA*UEqn.H(); phi = fvc::interpolate(rho) - *((fvc::interpolate(U) & mesh.Sf()) - fvc::meshPhi(rho, U)); + *((fvc::interpolate(U) & mesh.Sf()) - fvc::meshPhi(rho, U)); // Store pressure for under-relaxation p.storePrevIter(); @@ -24,7 +24,7 @@ if (nonOrth == nNonOrthCorr) { - phi += pEqn.flux(); + phi += pEqn.flux(); } } @@ -35,8 +35,8 @@ // rho does not carry working boundary conditions and needs to be updated // strictly according to the thermodynamics package // HJ, 22/Aug/2007 - thermo->correct(); - rho = thermo->rho(); + thermo.correct(); + rho = thermo.rho(); DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p); diff --git a/applications/solvers/engine/sonicTurbDyMEngineFoam/sonicTurbDyMEngineFoam.C b/applications/solvers/engine/sonicTurbDyMEngineFoam/sonicTurbDyMEngineFoam.C index 9b44ff56f..9b5f522d4 100644 --- a/applications/solvers/engine/sonicTurbDyMEngineFoam/sonicTurbDyMEngineFoam.C +++ b/applications/solvers/engine/sonicTurbDyMEngineFoam/sonicTurbDyMEngineFoam.C @@ -26,7 +26,7 @@ Application sonicTurbDyMEngineFoam Description - Solver for compressible cold flow in internal combustion engines + Solver for compressible cold flow in internal combustion engines with mesh motion and topological changes. @@ -36,8 +36,8 @@ Description #include "engineTime.H" #include "dynamicFvMesh.H" #include "engineTopoChangerMesh.H" -#include "basicThermo.H" -#include "compressible/RASModel/RASModel.H" +#include "basicPsiThermo.H" +#include "turbulenceModel.H" #include "Switch.H" #include "OFstream.H" @@ -61,8 +61,8 @@ int main(int argc, char *argv[]) // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - thermo->correct(); - + thermo.correct(); + Info << "\nStarting time loop\n" << endl; while (runTime.run()) @@ -72,24 +72,24 @@ int main(int argc, char *argv[]) # include "setDeltaT.H" runTime++; - + Info<< "Crank angle = " << runTime.theta() << " CA-deg" << endl; // make phi relative - + phi += meshFlux; - + bool meshChanged = mesh.update(); if(meshChanged) { - thermo->correct(); + thermo.correct(); # include "checkTotalVolume.H" # include "compressibleCorrectPhi.H" # include "CourantNo.H" } - + meshFlux = fvc::interpolate(rho)*fvc::meshPhi(rho, U); // Make phi absolute @@ -109,11 +109,11 @@ int main(int argc, char *argv[]) } turbulence->correct(); - + # include "logSummary.H" - rho = thermo->rho(); - + rho = thermo.rho(); + runTime.write(); # include "infoDataOutput.H" diff --git a/applications/solvers/engine/turbDyMEngineFoam/Make/options b/applications/solvers/engine/turbDyMEngineFoam/Make/options index 5582f6f44..e58f37d58 100644 --- a/applications/solvers/engine/turbDyMEngineFoam/Make/options +++ b/applications/solvers/engine/turbDyMEngineFoam/Make/options @@ -4,8 +4,9 @@ EXE_INC = \ -I$(LIB_SRC)/dynamicMesh/lnInclude \ -I$(LIB_SRC)/topoChangerFvMesh/lnInclude \ -I$(LIB_SRC)/meshTools/lnInclude \ - -I$(LIB_SRC)/turbulenceModels/RAS \ -I$(LIB_SRC)/transportModels \ + -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \ + -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \ -I$(LIB_SRC)/finiteVolume/lnInclude EXE_LIBS = \ @@ -15,6 +16,7 @@ EXE_LIBS = \ -ltopoChangerFvMesh \ -lmeshTools \ -lincompressibleRASModels \ + -lincompressibleLESModels \ -lincompressibleTransportModels \ -lfiniteVolume \ $(WM_DECOMP_LIBS) diff --git a/applications/solvers/engine/turbDyMEngineFoam/createFields.H b/applications/solvers/engine/turbDyMEngineFoam/createFields.H index c887b7ab2..2db24767d 100644 --- a/applications/solvers/engine/turbDyMEngineFoam/createFields.H +++ b/applications/solvers/engine/turbDyMEngineFoam/createFields.H @@ -39,9 +39,9 @@ singlePhaseTransportModel laminarTransport(U, phi); - autoPtr turbulence + autoPtr turbulence ( - incompressible::RASModel::New(U, phi, laminarTransport) + incompressible::turbulenceModel::New(U, phi, laminarTransport) ); Info<< "Reading field rUA if present\n" << endl; diff --git a/applications/solvers/engine/turbDyMEngineFoam/turbDyMEngineFoam.C b/applications/solvers/engine/turbDyMEngineFoam/turbDyMEngineFoam.C index 70abf3ba3..94ea3df8c 100644 --- a/applications/solvers/engine/turbDyMEngineFoam/turbDyMEngineFoam.C +++ b/applications/solvers/engine/turbDyMEngineFoam/turbDyMEngineFoam.C @@ -32,8 +32,8 @@ Description \*---------------------------------------------------------------------------*/ #include "fvCFD.H" -#include "incompressible/singlePhaseTransportModel/singlePhaseTransportModel.H" -#include "incompressible/RASModel/RASModel.H" +#include "singlePhaseTransportModel.H" +#include "turbulenceModel.H" #include "dynamicFvMesh.H" #include "engineTime.H" @@ -115,7 +115,7 @@ int main(int argc, char *argv[]) { pEqn.solve(mesh.solver(p.name())); } - + if (nonOrth == nNonOrthCorr) { phi -= pEqn.flux(); diff --git a/applications/solvers/heatTransfer/boussinesqBuoyantFoam/boussinesqBuoyantFoam.C b/applications/solvers/heatTransfer/boussinesqBuoyantFoam/boussinesqBuoyantFoam.C index 600236af6..5051eb54d 100644 --- a/applications/solvers/heatTransfer/boussinesqBuoyantFoam/boussinesqBuoyantFoam.C +++ b/applications/solvers/heatTransfer/boussinesqBuoyantFoam/boussinesqBuoyantFoam.C @@ -43,7 +43,7 @@ int main(int argc, char *argv[]) # include "createTime.H" # include "createMesh.H" # include "readTransportProperties.H" -# include "readEnvironmentalProperties.H" +# include "readGravitationalAcceleration.H" # include "createFields.H" # include "initContinuityErrs.H" diff --git a/applications/solvers/heatTransfer/buoyantPisoFoam/buoyantPisoFoam.C b/applications/solvers/heatTransfer/buoyantPisoFoam/buoyantPisoFoam.C index 85f88ce13..7e08febe3 100644 --- a/applications/solvers/heatTransfer/buoyantPisoFoam/buoyantPisoFoam.C +++ b/applications/solvers/heatTransfer/buoyantPisoFoam/buoyantPisoFoam.C @@ -43,15 +43,15 @@ Description int main(int argc, char *argv[]) { - #include "setRootCase.H" - #include "createTime.H" - #include "createMesh.H" - #include "readGravitationalAcceleration.H" - #include "createFields.H" - #include "initContinuityErrs.H" - #include "readTimeControls.H" - #include "compressibleCourantNo.H" - #include "setInitialDeltaT.H" +# include "setRootCase.H" +# include "createTime.H" +# include "createMesh.H" +# include "readGravitationalAcceleration.H" +# include "createFields.H" +# include "initContinuityErrs.H" +# include "readTimeControls.H" +# include "compressibleCourantNo.H" +# include "setInitialDeltaT.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -59,26 +59,26 @@ int main(int argc, char *argv[]) while (runTime.run()) { - #include "readTimeControls.H" - #include "readPISOControls.H" - #include "compressibleCourantNo.H" - #include "setDeltaT.H" +# include "readTimeControls.H" +# include "readPISOControls.H" +# include "compressibleCourantNo.H" +# include "setDeltaT.H" runTime++; Info<< "Time = " << runTime.timeName() << nl << endl; - #include "rhoEqn.H" +# include "rhoEqn.H" - #include "UEqn.H" +# include "UEqn.H" - #include "hEqn.H" +# include "hEqn.H" // --- PISO loop for (int corr=0; corrcorrect(); diff --git a/applications/solvers/multiphase/bubbleFoam/bubbleFoam.C b/applications/solvers/multiphase/bubbleFoam/bubbleFoam.C index 214602005..aa905ae78 100644 --- a/applications/solvers/multiphase/bubbleFoam/bubbleFoam.C +++ b/applications/solvers/multiphase/bubbleFoam/bubbleFoam.C @@ -26,7 +26,7 @@ Application bubbleFoam Description - Solver for a system of 2 incompressible fluid phases with one phase + Solver for a system of 2 incompressible fluid phases with one phase dispersed, e.g. gas bubbles in a liquid. \*---------------------------------------------------------------------------*/ @@ -40,12 +40,11 @@ Description int main(int argc, char *argv[]) { - # include "setRootCase.H" # include "createTime.H" # include "createMesh.H" -# include "readEnvironmentalProperties.H" +# include "readGravitationalAcceleration.H" # include "createFields.H" # include "initContinuityErrs.H" @@ -53,7 +52,7 @@ int main(int argc, char *argv[]) Info<< "\nStarting time loop\n" << endl; - for (runTime++; !runTime.end(); runTime++) + while (runTime.loop()) { Info<< "Time = " << runTime.timeName() << nl << endl; diff --git a/applications/solvers/multiphase/lesCavitatingFoam/CourantNo.H b/applications/solvers/multiphase/cavitatingFoam/CourantNo.H similarity index 100% rename from applications/solvers/multiphase/lesCavitatingFoam/CourantNo.H rename to applications/solvers/multiphase/cavitatingFoam/CourantNo.H diff --git a/applications/solvers/multiphase/cavitatingFoam/Make/files b/applications/solvers/multiphase/cavitatingFoam/Make/files new file mode 100644 index 000000000..832391f03 --- /dev/null +++ b/applications/solvers/multiphase/cavitatingFoam/Make/files @@ -0,0 +1,3 @@ +cavitatingFoam.C + +EXE = $(FOAM_APPBIN)/cavitatingFoam diff --git a/applications/solvers/multiphase/rasCavitatingFoam/Make/options b/applications/solvers/multiphase/cavitatingFoam/Make/options similarity index 81% rename from applications/solvers/multiphase/rasCavitatingFoam/Make/options rename to applications/solvers/multiphase/cavitatingFoam/Make/options index a694ac217..9cb749d63 100644 --- a/applications/solvers/multiphase/rasCavitatingFoam/Make/options +++ b/applications/solvers/multiphase/cavitatingFoam/Make/options @@ -3,12 +3,12 @@ EXE_INC = \ -I$(LIB_SRC)/transportModels \ -I$(LIB_SRC)/transportModels/incompressible/lnInclude \ -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \ - -I$(LIB_SRC)/turbulenceModels/RAS \ + -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \ -I$(LIB_SRC)/thermophysicalModels/barotropicCompressibilityModel/lnInclude EXE_LIBS = \ -lincompressibleTransportModels \ -lincompressibleRASModels \ + -lincompressibleLESModels \ -lfiniteVolume \ -lbarotropicCompressibilityModel - diff --git a/applications/solvers/multiphase/rasCavitatingFoam/UEqn.H b/applications/solvers/multiphase/cavitatingFoam/UEqn.H similarity index 86% rename from applications/solvers/multiphase/rasCavitatingFoam/UEqn.H rename to applications/solvers/multiphase/cavitatingFoam/UEqn.H index 374e410c2..01911faaa 100644 --- a/applications/solvers/multiphase/rasCavitatingFoam/UEqn.H +++ b/applications/solvers/multiphase/cavitatingFoam/UEqn.H @@ -14,7 +14,11 @@ - fvc::div(muEff*(fvc::interpolate(dev(fvc::grad(U))) & mesh.Sf())) ); + UEqn.relax(); + if (momentumPredictor) { solve(UEqn == -fvc::grad(p)); } + + Info<< "max(U) " << max(mag(U)).value() << endl; diff --git a/applications/solvers/multiphase/rasCavitatingFoam/rasCavitatingFoam.C b/applications/solvers/multiphase/cavitatingFoam/cavitatingFoam.C similarity index 70% rename from applications/solvers/multiphase/rasCavitatingFoam/rasCavitatingFoam.C rename to applications/solvers/multiphase/cavitatingFoam/cavitatingFoam.C index 8fa6a7a69..e149fdc59 100644 --- a/applications/solvers/multiphase/rasCavitatingFoam/rasCavitatingFoam.C +++ b/applications/solvers/multiphase/cavitatingFoam/cavitatingFoam.C @@ -23,56 +23,58 @@ License Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA Application - rasCavitatingFoam + cavitatingFoam Description - Transient cavitation code with RAS turbulence. + Transient cavitation code based on the homogeneous equilibrium model + from which the compressibility of the liquid/vapour "mixture" is obtained. + + Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected. \*---------------------------------------------------------------------------*/ #include "fvCFD.H" #include "barotropicCompressibilityModel.H" #include "twoPhaseMixture.H" -#include "incompressible/RASModel/RASModel.H" +#include "turbulenceModel.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // int main(int argc, char *argv[]) { + #include "setRootCase.H" -# include "setRootCase.H" + #include "createTime.H" + #include "createMesh.H" + #include "readThermodynamicProperties.H" + #include "readControls.H" + #include "createFields.H" + #include "initContinuityErrs.H" + #include "compressibleCourantNo.H" + #include "setInitialDeltaT.H" -# include "createTime.H" -# include "createMesh.H" -# include "readThermodynamicProperties.H" -# include "readControls.H" -# include "createFields.H" -# include "initContinuityErrs.H" -# include "compressibleCourantNo.H" -# include "setInitialDeltaT.H" - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // Info<< "\nStarting time loop\n" << endl; while (runTime.run()) { -# include "readControls.H" -# include "CourantNo.H" -# include "setDeltaT.H" + #include "readControls.H" + #include "CourantNo.H" + #include "setDeltaT.H" runTime++; Info<< "Time = " << runTime.timeName() << nl << endl; for (int outerCorr=0; outerCorr turbulence + // Create incompressible turbulence model + autoPtr turbulence ( - incompressible::RASModel::New(U, phiv, twoPhaseProperties) + incompressible::turbulenceModel::New(U, phiv, twoPhaseProperties) ); diff --git a/applications/solvers/multiphase/lesCavitatingFoam/gammaPsi.H b/applications/solvers/multiphase/cavitatingFoam/gammaPsi.H similarity index 100% rename from applications/solvers/multiphase/lesCavitatingFoam/gammaPsi.H rename to applications/solvers/multiphase/cavitatingFoam/gammaPsi.H diff --git a/applications/solvers/multiphase/rasCavitatingFoam/pEqn.H b/applications/solvers/multiphase/cavitatingFoam/pEqn.H similarity index 72% rename from applications/solvers/multiphase/rasCavitatingFoam/pEqn.H rename to applications/solvers/multiphase/cavitatingFoam/pEqn.H index c9382dfc0..bcaedb7a1 100644 --- a/applications/solvers/multiphase/rasCavitatingFoam/pEqn.H +++ b/applications/solvers/multiphase/cavitatingFoam/pEqn.H @@ -1,7 +1,7 @@ { if (nOuterCorr == 1) { - p = + p = ( rho - (1.0 - gamma)*rhol0 @@ -24,7 +24,7 @@ phiv -= phiGradp/rhof; -# include "resetPhivPatches.H" + #include "resetPhivPatches.H" for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) { @@ -37,7 +37,14 @@ - fvm::laplacian(rUAf, p) ); - pEqn.solve(); + if (corr == nCorr-1 && nonOrth == nNonOrthCorr) + { + pEqn.solve(mesh.solver(p.name() + "Final")); + } + else + { + pEqn.solve(mesh.solver(p.name())); + } if (nonOrth == nNonOrthCorr) { @@ -45,9 +52,32 @@ } } - Info<< "max-min p: " << max(p).value() + Info<< "Predicted p max-min : " << max(p).value() << " " << min(p).value() << endl; + rho == max + ( + psi*p + + (1.0 - gamma)*rhol0 + + ((gamma*psiv + (1.0 - gamma)*psil) - psi)*pSat, + rhoMin + ); + + #include "gammaPsi.H" + + p = + ( + rho + - (1.0 - gamma)*rhol0 + - ((gamma*psiv + (1.0 - gamma)*psil) - psi)*pSat + )/psi; + + p.correctBoundaryConditions(); + + Info<< "Phase-change corrected p max-min : " << max(p).value() + << " " << min(p).value() << endl; + + // Correct velocity U = HbyA - rUA*fvc::grad(p); @@ -63,18 +93,4 @@ U.correctBoundaryConditions(); Info<< "max(U) " << max(mag(U)).value() << endl; - - rho == max - ( - psi*p - + (1.0 - gamma)*rhol0 - + ((gamma*psiv + (1.0 - gamma)*psil) - psi)*pSat, - rhoMin - ); - - Info<< "max-min rho: " << max(rho).value() - << " " << min(rho).value() << endl; - -# include "gammaPsi.H" - } diff --git a/applications/solvers/multiphase/lesCavitatingFoam/readControls.H b/applications/solvers/multiphase/cavitatingFoam/readControls.H similarity index 100% rename from applications/solvers/multiphase/lesCavitatingFoam/readControls.H rename to applications/solvers/multiphase/cavitatingFoam/readControls.H diff --git a/applications/solvers/multiphase/lesCavitatingFoam/readThermodynamicProperties.H b/applications/solvers/multiphase/cavitatingFoam/readThermodynamicProperties.H similarity index 100% rename from applications/solvers/multiphase/lesCavitatingFoam/readThermodynamicProperties.H rename to applications/solvers/multiphase/cavitatingFoam/readThermodynamicProperties.H diff --git a/applications/solvers/multiphase/lesCavitatingFoam/resetPhiPatches.H b/applications/solvers/multiphase/cavitatingFoam/resetPhiPatches.H similarity index 100% rename from applications/solvers/multiphase/lesCavitatingFoam/resetPhiPatches.H rename to applications/solvers/multiphase/cavitatingFoam/resetPhiPatches.H diff --git a/applications/solvers/multiphase/lesCavitatingFoam/resetPhivPatches.H b/applications/solvers/multiphase/cavitatingFoam/resetPhivPatches.H similarity index 100% rename from applications/solvers/multiphase/lesCavitatingFoam/resetPhivPatches.H rename to applications/solvers/multiphase/cavitatingFoam/resetPhivPatches.H diff --git a/applications/solvers/multiphase/lesCavitatingFoam/rhoEqn.H b/applications/solvers/multiphase/cavitatingFoam/rhoEqn.H similarity index 100% rename from applications/solvers/multiphase/lesCavitatingFoam/rhoEqn.H rename to applications/solvers/multiphase/cavitatingFoam/rhoEqn.H diff --git a/applications/solvers/multiphase/lesCavitatingFoam/setDeltaT.H b/applications/solvers/multiphase/cavitatingFoam/setDeltaT.H similarity index 100% rename from applications/solvers/multiphase/lesCavitatingFoam/setDeltaT.H rename to applications/solvers/multiphase/cavitatingFoam/setDeltaT.H diff --git a/applications/solvers/multiphase/lesCavitatingFoam/setInitialDeltaT.H b/applications/solvers/multiphase/cavitatingFoam/setInitialDeltaT.H similarity index 100% rename from applications/solvers/multiphase/lesCavitatingFoam/setInitialDeltaT.H rename to applications/solvers/multiphase/cavitatingFoam/setInitialDeltaT.H diff --git a/applications/solvers/multiphase/compressibleInterDyMFoam/Make/files b/applications/solvers/multiphase/compressibleInterDyMFoam/Make/files new file mode 100644 index 000000000..121264b1a --- /dev/null +++ b/applications/solvers/multiphase/compressibleInterDyMFoam/Make/files @@ -0,0 +1,3 @@ +compressibleInterDyMFoam.C + +EXE = $(FOAM_APPBIN)/compressibleInterDyMFoam diff --git a/applications/solvers/multiphase/lesCavitatingFoam/Make/options b/applications/solvers/multiphase/compressibleInterDyMFoam/Make/options similarity index 51% rename from applications/solvers/multiphase/lesCavitatingFoam/Make/options rename to applications/solvers/multiphase/compressibleInterDyMFoam/Make/options index 91833c136..af146aaa3 100644 --- a/applications/solvers/multiphase/lesCavitatingFoam/Make/options +++ b/applications/solvers/multiphase/compressibleInterDyMFoam/Make/options @@ -1,16 +1,20 @@ EXE_INC = \ - -I$(LIB_SRC)/finiteVolume/lnInclude \ -I$(LIB_SRC)/transportModels \ -I$(LIB_SRC)/transportModels/incompressible/lnInclude \ -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \ - -I$(LIB_SRC)/turbulenceModels/LES \ - -I$(LIB_SRC)/turbulenceModels/LES/incompressible/lnInclude \ - -I$(LIB_SRC)/turbulenceModels/LES/LESdeltas/lnInclude \ - -I$(LIB_SRC)/thermophysicalModels/barotropicCompressibilityModel/lnInclude + -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \ + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/dynamicMesh/lnInclude \ + -I$(LIB_SRC)/meshTools/lnInclude \ + -I$(LIB_SRC)/dynamicFvMesh/lnInclude EXE_LIBS = \ + -linterfaceProperties \ -lincompressibleTransportModels \ + -lincompressibleRASModels \ -lincompressibleLESModels \ -lfiniteVolume \ - -lbarotropicCompressibilityModel + -ldynamicMesh \ + -lmeshTools \ + -ldynamicFvMesh diff --git a/applications/solvers/multiphase/lesInterFoam/UEqn.H b/applications/solvers/multiphase/compressibleInterDyMFoam/UEqn.H similarity index 66% rename from applications/solvers/multiphase/lesInterFoam/UEqn.H rename to applications/solvers/multiphase/compressibleInterDyMFoam/UEqn.H index 1bf56f02d..138e58fc7 100644 --- a/applications/solvers/multiphase/lesInterFoam/UEqn.H +++ b/applications/solvers/multiphase/compressibleInterDyMFoam/UEqn.H @@ -2,7 +2,7 @@ ( "muEff", twoPhaseProperties.muf() - + fvc::interpolate(rho*turbulence->nuSgs()) + + fvc::interpolate(rho*turbulence->nut()) ); fvVectorMatrix UEqn @@ -11,9 +11,11 @@ + fvm::div(rhoPhi, U) - fvm::laplacian(muEff, U) - (fvc::grad(U) & fvc::grad(muEff)) - //- fvc::div(muEff*(fvc::interpolate(dev(fvc::grad(U))) & mesh.Sf())) + //- fvc::div(muf*(mesh.Sf() & fvc::interpolate(fvc::grad(U)().T()))) ); + UEqn.relax(); + if (momentumPredictor) { solve @@ -22,10 +24,10 @@ == fvc::reconstruct ( - ( - fvc::interpolate(interface.sigmaK())*fvc::snGrad(gamma) - - ghf*fvc::snGrad(rho) - - fvc::snGrad(pd) + fvc::interpolate(rho)*(g & mesh.Sf()) + + ( + fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1) + - fvc::snGrad(p) ) * mesh.magSf() ) ); diff --git a/applications/solvers/multiphase/compressibleLesInterFoam/alphaEqns.H b/applications/solvers/multiphase/compressibleInterDyMFoam/alphaEqns.H similarity index 95% rename from applications/solvers/multiphase/compressibleLesInterFoam/alphaEqns.H rename to applications/solvers/multiphase/compressibleInterDyMFoam/alphaEqns.H index 819cd0f53..dd704c069 100644 --- a/applications/solvers/multiphase/compressibleLesInterFoam/alphaEqns.H +++ b/applications/solvers/multiphase/compressibleInterDyMFoam/alphaEqns.H @@ -59,7 +59,7 @@ alpharScheme ); - MULES::explicitSolve(oneField(), alpha1, phi, phiAlpha1, Sp, Su, 1, 0); + MULES::explicitSolve(geometricOneField(), alpha1, phi, phiAlpha1, Sp, Su, 1, 0); surfaceScalarField rho1f = fvc::interpolate(rho1); surfaceScalarField rho2f = fvc::interpolate(rho2); diff --git a/applications/solvers/multiphase/compressibleLesInterFoam/alphaEqnsSubCycle.H b/applications/solvers/multiphase/compressibleInterDyMFoam/alphaEqnsSubCycle.H similarity index 87% rename from applications/solvers/multiphase/compressibleLesInterFoam/alphaEqnsSubCycle.H rename to applications/solvers/multiphase/compressibleInterDyMFoam/alphaEqnsSubCycle.H index c52dce969..32a716313 100644 --- a/applications/solvers/multiphase/compressibleLesInterFoam/alphaEqnsSubCycle.H +++ b/applications/solvers/multiphase/compressibleInterDyMFoam/alphaEqnsSubCycle.H @@ -10,9 +10,11 @@ ); surfaceScalarField phic = mag(phi/mesh.magSf()); - phic = min(interface.cGamma()*phic, max(phic)); + phic = min(interface.cAlpha()*phic, max(phic)); + fvc::makeAbsolute(phi, U); volScalarField divU = fvc::div(phi); + fvc::makeRelative(phi, U); if (nAlphaSubCycles > 1) { diff --git a/applications/solvers/multiphase/compressibleInterDyMFoam/compressibleInterDyMFoam.C b/applications/solvers/multiphase/compressibleInterDyMFoam/compressibleInterDyMFoam.C new file mode 100644 index 000000000..c4cfbc48e --- /dev/null +++ b/applications/solvers/multiphase/compressibleInterDyMFoam/compressibleInterDyMFoam.C @@ -0,0 +1,142 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright held by original author + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM 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 2 of the License, or (at your + option) any later version. + + OpenFOAM 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 OpenFOAM; if not, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +Application + compressibleLesInterFoam + +Description + Solver for 2 compressible, isothermal immiscible fluids using a VOF + (volume of fluid) phase-fraction based interface capturing approach, + with optional mesh motion and mesh topology changes including adaptive + re-meshing. + + The momentum and other fluid properties are of the "mixture" and a + single momentum equation is solved. Turbulence modelling is generic, + i.e. laminar, RAS or LES may be selected. + +\*---------------------------------------------------------------------------*/ + +#include "fvCFD.H" +#include "dynamicFvMesh.H" +#include "MULES.H" +#include "subCycle.H" +#include "interfaceProperties.H" +#include "twoPhaseMixture.H" +#include "turbulenceModel.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +int main(int argc, char *argv[]) +{ + #include "setRootCase.H" + #include "createTime.H" + #include "createDynamicFvMesh.H" + #include "readGravitationalAcceleration.H" + #include "readControls.H" + #include "initContinuityErrs.H" + #include "createFields.H" + #include "CourantNo.H" + #include "setInitialDeltaT.H" + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + Info<< "\nStarting time loop\n" << endl; + + while (runTime.run()) + { + #include "readControls.H" + #include "CourantNo.H" + + // Make the fluxes absolute + fvc::makeAbsolute(phi, U); + + #include "setDeltaT.H" + + runTime++; + + Info<< "Time = " << runTime.timeName() << nl << endl; + + { + // Store divU from the previous mesh for the correctPhi + volScalarField divU = fvc::div(phi); + + scalar timeBeforeMeshUpdate = runTime.elapsedCpuTime(); + + // Do any mesh changes + mesh.update(); + + if (mesh.changing()) + { + Info<< "Execution time for mesh.update() = " + << runTime.elapsedCpuTime() - timeBeforeMeshUpdate + << " s" << endl; + } + + if (mesh.changing() && correctPhi) + { + #include "correctPhi.H" + } + } + + // Make the fluxes relative to the mesh motion + fvc::makeRelative(phi, U); + + if (mesh.changing() && checkMeshCourantNo) + { + #include "meshCourantNo.H" + } + + turbulence->correct(); + + // --- Outer-corrector loop + for (int oCorr=0; oCorr turbulence + // Construct incompressible turbulence model + autoPtr turbulence ( - incompressible::LESModel::New(U, phi, twoPhaseProperties) + incompressible::turbulenceModel::New(U, phi, twoPhaseProperties) ); + + + wordList pcorrTypes + ( + p.boundaryField().size(), + zeroGradientFvPatchScalarField::typeName + ); + + for (label i=0; i pEqnComp; + + if (transonic) + { + pEqnComp = + (fvm::ddt(p) + fvm::div(phi, p) - fvm::Sp(fvc::div(phi), p)); + } + else + { + pEqnComp = + (fvm::ddt(p) + fvc::div(phi, p) - fvc::Sp(fvc::div(phi), p)); + } + + + U = rUA*UEqn.H(); + + surfaceScalarField phiU + ( + "phiU", + (fvc::interpolate(U) & mesh.Sf()) + ); + + phi = phiU + + ( + fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)*mesh.magSf() + + fvc::interpolate(rho)*(g & mesh.Sf()) + )*rUAf; + + for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) + { + fvScalarMatrix pEqnIncomp + ( + fvc::div(phi) + - fvm::laplacian(rUAf, p) + ); + + if + ( + oCorr == nOuterCorr-1 + && corr == nCorr-1 + && nonOrth == nNonOrthCorr + ) + { + solve + ( + ( + max(alpha1, scalar(0))*(psi1/rho1) + + max(alpha2, scalar(0))*(psi2/rho2) + ) + *pEqnComp() + + pEqnIncomp, + mesh.solver(p.name() + "Final") + ); + } + else + { + solve + ( + ( + max(alpha1, scalar(0))*(psi1/rho1) + + max(alpha2, scalar(0))*(psi2/rho2) + ) + *pEqnComp() + + pEqnIncomp + ); + } + + if (nonOrth == nNonOrthCorr) + { + dgdt = + (pos(alpha2)*(psi2/rho2) - pos(alpha1)*(psi1/rho1)) + *(pEqnComp & p); + phi += pEqnIncomp.flux(); + } + } + + U += rUA*fvc::reconstruct((phi - phiU)/rUAf); + U.correctBoundaryConditions(); + + p.max(pMin); + + rho1 = rho10 + psi1*p; + rho2 = rho20 + psi2*p; + + Info<< "max(U) " << max(mag(U)).value() << endl; + Info<< "min(p) " << min(p).value() << endl; + + // Make the fluxes relative to the mesh motion + fvc::makeRelative(phi, U); +} diff --git a/applications/solvers/multiphase/compressibleInterDyMFoam/readControls.H b/applications/solvers/multiphase/compressibleInterDyMFoam/readControls.H new file mode 100644 index 000000000..a2e4ef374 --- /dev/null +++ b/applications/solvers/multiphase/compressibleInterDyMFoam/readControls.H @@ -0,0 +1,32 @@ + #include "readPISOControls.H" + #include "readTimeControls.H" + + label nAlphaCorr + ( + readLabel(piso.lookup("nAlphaCorr")) + ); + + label nAlphaSubCycles + ( + readLabel(piso.lookup("nAlphaSubCycles")) + ); + + if (nAlphaSubCycles > 1 && nOuterCorr != 1) + { + FatalErrorIn(args.executable()) + << "Sub-cycling alpha is only allowed for PISO, " + "i.e. when the number of outer-correctors = 1" + << exit(FatalError); + } + + bool correctPhi = true; + if (piso.found("correctPhi")) + { + correctPhi = Switch(piso.lookup("correctPhi")); + } + + bool checkMeshCourantNo = false; + if (piso.found("checkMeshCourantNo")) + { + checkMeshCourantNo = Switch(piso.lookup("checkMeshCourantNo")); + } diff --git a/applications/solvers/multiphase/compressibleInterFoam/Make/files b/applications/solvers/multiphase/compressibleInterFoam/Make/files new file mode 100644 index 000000000..de5437219 --- /dev/null +++ b/applications/solvers/multiphase/compressibleInterFoam/Make/files @@ -0,0 +1,3 @@ +compressibleInterFoam.C + +EXE = $(FOAM_APPBIN)/compressibleInterFoam diff --git a/applications/solvers/multiphase/rasInterFoam/Make/options b/applications/solvers/multiphase/compressibleInterFoam/Make/options similarity index 73% rename from applications/solvers/multiphase/rasInterFoam/Make/options rename to applications/solvers/multiphase/compressibleInterFoam/Make/options index bd701d51d..0e21fc2d6 100644 --- a/applications/solvers/multiphase/rasInterFoam/Make/options +++ b/applications/solvers/multiphase/compressibleInterFoam/Make/options @@ -1,14 +1,13 @@ EXE_INC = \ - -I../interFoam \ -I$(LIB_SRC)/transportModels \ -I$(LIB_SRC)/transportModels/incompressible/lnInclude \ -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \ - -I$(LIB_SRC)/turbulenceModels/RAS \ + -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \ -I$(LIB_SRC)/finiteVolume/lnInclude EXE_LIBS = \ -linterfaceProperties \ -lincompressibleTransportModels \ -lincompressibleRASModels \ - -lfiniteVolume \ - -llduSolvers + -lincompressibleLESModels \ + -lfiniteVolume diff --git a/applications/solvers/multiphase/rasInterFoam/UEqn.H b/applications/solvers/multiphase/compressibleInterFoam/UEqn.H similarity index 84% rename from applications/solvers/multiphase/rasInterFoam/UEqn.H rename to applications/solvers/multiphase/compressibleInterFoam/UEqn.H index e5b553bfe..0b1a9ac02 100644 --- a/applications/solvers/multiphase/rasInterFoam/UEqn.H +++ b/applications/solvers/multiphase/compressibleInterFoam/UEqn.H @@ -24,10 +24,10 @@ == fvc::reconstruct ( - ( - fvc::interpolate(interface.sigmaK())*fvc::snGrad(gamma) - - ghf*fvc::snGrad(rho) - - fvc::snGrad(pd) + fvc::interpolate(rho)*(g & mesh.Sf()) + + ( + fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1) + - fvc::snGrad(p) ) * mesh.magSf() ) ); diff --git a/applications/solvers/multiphase/compressibleInterFoam/alphaEqns.H b/applications/solvers/multiphase/compressibleInterFoam/alphaEqns.H new file mode 100644 index 000000000..dd704c069 --- /dev/null +++ b/applications/solvers/multiphase/compressibleInterFoam/alphaEqns.H @@ -0,0 +1,76 @@ +{ + word alphaScheme("div(phi,alpha)"); + word alpharScheme("div(phirb,alpha)"); + + surfaceScalarField phir = phic*interface.nHatf(); + + for (int gCorr=0; gCorr 0.0 && alpha1[celli] > 0.0) + { + Sp[celli] -= dgdt[celli]*alpha1[celli]; + Su[celli] += dgdt[celli]*alpha1[celli]; + } + else if (dgdt[celli] < 0.0 && alpha1[celli] < 1.0) + { + Sp[celli] += dgdt[celli]*(1.0 - alpha1[celli]); + } + } + + + surfaceScalarField phiAlpha1 = + fvc::flux + ( + phi, + alpha1, + alphaScheme + ) + + fvc::flux + ( + -fvc::flux(-phir, alpha2, alpharScheme), + alpha1, + alpharScheme + ); + + MULES::explicitSolve(geometricOneField(), alpha1, phi, phiAlpha1, Sp, Su, 1, 0); + + surfaceScalarField rho1f = fvc::interpolate(rho1); + surfaceScalarField rho2f = fvc::interpolate(rho2); + rhoPhi = phiAlpha1*(rho1f - rho2f) + phi*rho2f; + + alpha2 = scalar(1) - alpha1; + } + + Info<< "Liquid phase volume fraction = " + << alpha1.weightedAverage(mesh.V()).value() + << " Min(alpha1) = " << min(alpha1).value() + << " Min(alpha2) = " << min(alpha2).value() + << endl; +} diff --git a/applications/solvers/multiphase/compressibleInterFoam/alphaEqnsSubCycle.H b/applications/solvers/multiphase/compressibleInterFoam/alphaEqnsSubCycle.H new file mode 100644 index 000000000..89ba7a4e7 --- /dev/null +++ b/applications/solvers/multiphase/compressibleInterFoam/alphaEqnsSubCycle.H @@ -0,0 +1,43 @@ +{ + label nAlphaCorr + ( + readLabel(piso.lookup("nAlphaCorr")) + ); + + label nAlphaSubCycles + ( + readLabel(piso.lookup("nAlphaSubCycles")) + ); + + surfaceScalarField phic = mag(phi/mesh.magSf()); + phic = min(interface.cAlpha()*phic, max(phic)); + + volScalarField divU = fvc::div(phi); + + if (nAlphaSubCycles > 1) + { + dimensionedScalar totalDeltaT = runTime.deltaT(); + surfaceScalarField rhoPhiSum = 0.0*rhoPhi; + + for + ( + subCycle alphaSubCycle(alpha1, nAlphaSubCycles); + !(++alphaSubCycle).end(); + ) + { + #include "alphaEqns.H" + rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi; + } + + rhoPhi = rhoPhiSum; + } + else + { + #include "alphaEqns.H" + } + + if (oCorr == 0) + { + interface.correct(); + } +} diff --git a/applications/solvers/multiphase/compressibleLesInterFoam/compressibleLesInterFoam.C b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C similarity index 92% rename from applications/solvers/multiphase/compressibleLesInterFoam/compressibleLesInterFoam.C rename to applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C index becf85fa9..0689d1ec0 100644 --- a/applications/solvers/multiphase/compressibleLesInterFoam/compressibleLesInterFoam.C +++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C @@ -23,14 +23,16 @@ License Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA Application - compressibleLesInterFoam + compressibleInterFoam Description Solver for 2 compressible, isothermal immiscible fluids using a VOF (volume of fluid) phase-fraction based interface capturing approach. + The momentum and other fluid properties are of the "mixture" and a single - momentum equation is solved. Turbulence is modelled using a run-time - selectable incompressible LES model. + momentum equation is solved. + + Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected. \*---------------------------------------------------------------------------*/ @@ -39,7 +41,7 @@ Description #include "subCycle.H" #include "interfaceProperties.H" #include "twoPhaseMixture.H" -#include "incompressible/LESModel/LESModel.H" +#include "turbulenceModel.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -48,7 +50,7 @@ int main(int argc, char *argv[]) #include "setRootCase.H" #include "createTime.H" #include "createMesh.H" - #include "readEnvironmentalProperties.H" + #include "readGravitationalAcceleration.H" #include "readControls.H" #include "initContinuityErrs.H" #include "createFields.H" @@ -69,8 +71,6 @@ int main(int argc, char *argv[]) Info<< "Time = " << runTime.timeName() << nl << endl; - turbulence->correct(); - // --- Outer-corrector loop for (int oCorr=0; oCorrcorrect(); + runTime.write(); Info<< "ExecutionTime = " @@ -98,7 +100,7 @@ int main(int argc, char *argv[]) Info<< "End\n" << endl; - return(0); + return 0; } diff --git a/applications/solvers/multiphase/compressibleInterFoam/createFields.H b/applications/solvers/multiphase/compressibleInterFoam/createFields.H new file mode 100644 index 000000000..3e6904d38 --- /dev/null +++ b/applications/solvers/multiphase/compressibleInterFoam/createFields.H @@ -0,0 +1,129 @@ + Info<< "Reading field p\n" << endl; + volScalarField p + ( + IOobject + ( + "p", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh + ); + + Info<< "Reading field alpha1\n" << endl; + volScalarField alpha1 + ( + IOobject + ( + "alpha1", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh + ); + + Info<< "Calculating field alpha1\n" << endl; + volScalarField alpha2("alpha2", scalar(1) - alpha1); + + Info<< "Reading field U\n" << endl; + volVectorField U + ( + IOobject + ( + "U", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh + ); + + #include "createPhi.H" + + + Info<< "Reading transportProperties\n" << endl; + twoPhaseMixture twoPhaseProperties(U, phi); + + dimensionedScalar rho10 + ( + twoPhaseProperties.subDict + ( + twoPhaseProperties.phase1Name() + ).lookup("rho0") + ); + + dimensionedScalar rho20 + ( + twoPhaseProperties.subDict + ( + twoPhaseProperties.phase2Name() + ).lookup("rho0") + ); + + dimensionedScalar psi1 + ( + twoPhaseProperties.subDict + ( + twoPhaseProperties.phase1Name() + ).lookup("psi") + ); + + dimensionedScalar psi2 + ( + twoPhaseProperties.subDict + ( + twoPhaseProperties.phase2Name() + ).lookup("psi") + ); + + dimensionedScalar pMin(twoPhaseProperties.lookup("pMin")); + + volScalarField rho1 = rho10 + psi1*p; + volScalarField rho2 = rho20 + psi2*p; + + volScalarField rho + ( + IOobject + ( + "rho", + runTime.timeName(), + mesh, + IOobject::READ_IF_PRESENT, + IOobject::AUTO_WRITE + ), + alpha1*rho1 + alpha2*rho2 + ); + + + // Mass flux + // Initialisation does not matter because rhoPhi is reset after the + // alpha1 solution before it is used in the U equation. + surfaceScalarField rhoPhi + ( + IOobject + ( + "rho*phi", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + fvc::interpolate(rho)*phi + ); + + volScalarField dgdt = + pos(alpha2)*fvc::div(phi)/max(alpha2, scalar(0.0001)); + + // Construct interface from alpha1 distribution + interfaceProperties interface(alpha1, U, twoPhaseProperties); + + // Construct incompressible turbulence model + autoPtr turbulence + ( + incompressible::turbulenceModel::New(U, phi, twoPhaseProperties) + ); diff --git a/applications/solvers/multiphase/compressibleLesInterFoam/pEqn.H b/applications/solvers/multiphase/compressibleInterFoam/pEqn.H similarity index 60% rename from applications/solvers/multiphase/compressibleLesInterFoam/pEqn.H rename to applications/solvers/multiphase/compressibleInterFoam/pEqn.H index ebf24498a..9d2dc2391 100644 --- a/applications/solvers/multiphase/compressibleLesInterFoam/pEqn.H +++ b/applications/solvers/multiphase/compressibleInterFoam/pEqn.H @@ -2,17 +2,17 @@ volScalarField rUA = 1.0/UEqn.A(); surfaceScalarField rUAf = fvc::interpolate(rUA); - tmp pdEqnComp; + tmp pEqnComp; if (transonic) { - pdEqnComp = - (fvm::ddt(pd) + fvm::div(phi, pd) - fvm::Sp(fvc::div(phi), pd)); + pEqnComp = + (fvm::ddt(p) + fvm::div(phi, p) - fvm::Sp(fvc::div(phi), p)); } else { - pdEqnComp = - (fvm::ddt(pd) + fvc::div(phi, pd) - fvc::Sp(fvc::div(phi), pd)); + pEqnComp = + (fvm::ddt(p) + fvc::div(phi, p) - fvc::Sp(fvc::div(phi), p)); } @@ -26,16 +26,16 @@ phi = phiU + ( - fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1) - - ghf*fvc::snGrad(rho) - )*rUAf*mesh.magSf(); + fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)*mesh.magSf() + + fvc::interpolate(rho)*(g & mesh.Sf()) + )*rUAf; for(int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) { - fvScalarMatrix pdEqnIncomp + fvScalarMatrix pEqnIncomp ( fvc::div(phi) - - fvm::laplacian(rUAf, pd) + - fvm::laplacian(rUAf, p) ); solve @@ -44,31 +44,27 @@ max(alpha1, scalar(0))*(psi1/rho1) + max(alpha2, scalar(0))*(psi2/rho2) ) - *pdEqnComp() - + pdEqnIncomp + *pEqnComp() + + pEqnIncomp ); if (nonOrth == nNonOrthCorr) { dgdt = (pos(alpha2)*(psi2/rho2) - pos(alpha1)*(psi1/rho1)) - *(pdEqnComp & pd); - phi += pdEqnIncomp.flux(); + *(pEqnComp & p); + phi += pEqnIncomp.flux(); } } U += rUA*fvc::reconstruct((phi - phiU)/rUAf); U.correctBoundaryConditions(); - p = max - ( - (pd + gh*(alpha1*rho10 + alpha2*rho20))/(1.0 - gh*(alpha1*psi1 + alpha2*psi2)), - pMin - ); + p.max(pMin); rho1 = rho10 + psi1*p; rho2 = rho20 + psi2*p; Info<< "max(U) " << max(mag(U)).value() << endl; - Info<< "min(pd) " << min(pd).value() << endl; + Info<< "min(p) " << min(p).value() << endl; } diff --git a/applications/solvers/multiphase/compressibleLesInterFoam/readControls.H b/applications/solvers/multiphase/compressibleInterFoam/readControls.H similarity index 100% rename from applications/solvers/multiphase/compressibleLesInterFoam/readControls.H rename to applications/solvers/multiphase/compressibleInterFoam/readControls.H diff --git a/applications/solvers/multiphase/compressibleLesInterFoam/Make/files b/applications/solvers/multiphase/compressibleLesInterFoam/Make/files deleted file mode 100644 index 05dafae8b..000000000 --- a/applications/solvers/multiphase/compressibleLesInterFoam/Make/files +++ /dev/null @@ -1,3 +0,0 @@ -compressibleLesInterFoam.C - -EXE = $(FOAM_APPBIN)/compressibleLesInterFoam diff --git a/applications/solvers/multiphase/compressibleLesInterFoam/UEqn.H b/applications/solvers/multiphase/compressibleLesInterFoam/UEqn.H deleted file mode 100644 index 089591736..000000000 --- a/applications/solvers/multiphase/compressibleLesInterFoam/UEqn.H +++ /dev/null @@ -1,29 +0,0 @@ - surfaceScalarField muf = - twoPhaseProperties.muf() - + fvc::interpolate(rho*turbulence->nuSgs()); - - fvVectorMatrix UEqn - ( - fvm::ddt(rho, U) - + fvm::div(rhoPhi, U) - - fvm::laplacian(muf, U) - - (fvc::grad(U) & fvc::grad(muf)) - //- fvc::div(muf*(mesh.Sf() & fvc::interpolate(fvc::grad(U)().T()))) - ); - - if (momentumPredictor) - { - solve - ( - UEqn - == - fvc::reconstruct - ( - ( - fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1) - - ghf*fvc::snGrad(rho) - - fvc::snGrad(pd) - ) * mesh.magSf() - ) - ); - } diff --git a/applications/solvers/multiphase/interDyMFoam/Make/options b/applications/solvers/multiphase/interDyMFoam/Make/options index 48a32a49d..39b00bbd9 100644 --- a/applications/solvers/multiphase/interDyMFoam/Make/options +++ b/applications/solvers/multiphase/interDyMFoam/Make/options @@ -1,9 +1,9 @@ EXE_INC = \ + -I../interFoam \ -I$(LIB_SRC)/transportModels \ -I$(LIB_SRC)/transportModels/incompressible/lnInclude \ -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \ - -I$(LIB_SRC)/turbulenceModels/RAS \ - -I$(LIB_SRC)/turbulenceModels/RAS/incompressible/lnInclude \ + -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \ -I$(LIB_SRC)/finiteVolume/lnInclude \ -I$(LIB_SRC)/dynamicMesh/lnInclude \ -I$(LIB_SRC)/meshTools/lnInclude \ @@ -13,8 +13,10 @@ EXE_LIBS = \ -linterfaceProperties \ -lincompressibleTransportModels \ -lincompressibleRASModels \ + -lincompressibleLESModels \ -lfiniteVolume \ -ldynamicMesh \ -lmeshTools \ - -ldynamicFvMesh - + -ldynamicFvMesh \ + -ltopoChangerFvMesh \ + -llduSolvers diff --git a/applications/solvers/multiphase/interDyMFoam/UEqn.H b/applications/solvers/multiphase/interDyMFoam/UEqn.H deleted file mode 100644 index 7e4a29468..000000000 --- a/applications/solvers/multiphase/interDyMFoam/UEqn.H +++ /dev/null @@ -1,29 +0,0 @@ - surfaceScalarField muf = twoPhaseProperties.muf(); - - fvVectorMatrix UEqn - ( - fvm::ddt(rho, U) - + fvm::div(rhoPhi, U) - - fvm::laplacian(muf, U) - - (fvc::grad(U) & fvc::grad(muf)) - //- fvc::div(muf*(fvc::interpolate(dev(fvc::grad(U))) & mesh.Sf())) - ); - - UEqn.relax(); - - if (momentumPredictor) - { - solve - ( - UEqn - == - fvc::reconstruct - ( - ( - fvc::interpolate(interface.sigmaK())*fvc::snGrad(gamma) - - ghf*fvc::snGrad(rho) - - fvc::snGrad(pd) - ) * mesh.magSf() - ) - ); - } diff --git a/applications/solvers/multiphase/interDyMFoam/correctPhi.H b/applications/solvers/multiphase/interDyMFoam/correctPhi.H index c975c9b37..88472df9b 100644 --- a/applications/solvers/multiphase/interDyMFoam/correctPhi.H +++ b/applications/solvers/multiphase/interDyMFoam/correctPhi.H @@ -1,4 +1,26 @@ { + if (mesh.changing()) + { + forAll(U.boundaryField(), patchi) + { + if (U.boundaryField()[patchi].fixesValue()) + { + U.boundaryField()[patchi].initEvaluate(); + } + } + + forAll(U.boundaryField(), patchi) + { + if (U.boundaryField()[patchi].fixesValue()) + { + U.boundaryField()[patchi].evaluate(); + + phi.boundaryField()[patchi] = + U.boundaryField()[patchi] & mesh.Sf().boundaryField()[patchi]; + } + } + } + # include "continuityErrs.H" volScalarField pcorr diff --git a/applications/solvers/multiphase/interDyMFoam/createFields.H b/applications/solvers/multiphase/interDyMFoam/createFields.H index c2b7ee49e..afadd3666 100644 --- a/applications/solvers/multiphase/interDyMFoam/createFields.H +++ b/applications/solvers/multiphase/interDyMFoam/createFields.H @@ -12,12 +12,12 @@ mesh ); - Info<< "Reading field gamma\n" << endl; - volScalarField gamma + Info<< "Reading field alpha1\n" << endl; + volScalarField alpha1 ( IOobject ( - "gamma", + "alpha1", runTime.timeName(), mesh, IOobject::MUST_READ, @@ -44,7 +44,7 @@ Info<< "Reading transportProperties\n" << endl; - twoPhaseMixture twoPhaseProperties(U, phi, "gamma"); + twoPhaseMixture twoPhaseProperties(U, phi, "alpha1"); const dimensionedScalar& rho1 = twoPhaseProperties.rho1(); const dimensionedScalar& rho2 = twoPhaseProperties.rho2(); @@ -60,15 +60,15 @@ mesh, IOobject::READ_IF_PRESENT ), - gamma*rho1 + (scalar(1) - gamma)*rho2, - gamma.boundaryField().types() + alpha1*rho1 + (scalar(1) - alpha1)*rho2, + alpha1.boundaryField().types() ); rho.oldTime(); // Mass flux // Initialisation does not matter because rhoPhi is reset after the - // gamma solution before it is used in the U equation. + // alpha1 solution before it is used in the U equation. surfaceScalarField rhoPhi ( IOobject @@ -83,18 +83,22 @@ ); - // Construct interface from gamma distribution - interfaceProperties interface(gamma, U, twoPhaseProperties); + // Construct interface from alpha1 distribution + interfaceProperties interface(alpha1, U, twoPhaseProperties); - // Construct incompressible RAS model - autoPtr turbulence + // Construct incompressible turbulence model + autoPtr turbulence ( - incompressible::RASModel::New(U, phi, twoPhaseProperties) + incompressible::turbulenceModel::New(U, phi, twoPhaseProperties) ); - wordList pcorrTypes(pd.boundaryField().types()); + wordList pcorrTypes + ( + pd.boundaryField().size(), + zeroGradientFvPatchScalarField::typeName + ); - for (label i=0; i 1) -{ - dimensionedScalar totalDeltaT = runTime.deltaT(); - surfaceScalarField rhoPhiSum = 0.0*rhoPhi; - - for - ( - subCycle gammaSubCycle(gamma, nGammaSubCycles); - !(++gammaSubCycle).end(); - ) - { -# include "gammaEqn.H" - rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi; - } - - rhoPhi = rhoPhiSum; -} -else -{ -# include "gammaEqn.H" -} - -interface.correct(); - -rho == gamma*rho1 + (scalar(1) - gamma)*rho2; diff --git a/applications/solvers/multiphase/interDyMFoam/interDyMFoam.C b/applications/solvers/multiphase/interDyMFoam/interDyMFoam.C index 48f7be461..285059a8c 100644 --- a/applications/solvers/multiphase/interDyMFoam/interDyMFoam.C +++ b/applications/solvers/multiphase/interDyMFoam/interDyMFoam.C @@ -39,8 +39,7 @@ Description #include "subCycle.H" #include "interfaceProperties.H" #include "twoPhaseMixture.H" -#include "incompressible/RASModel/RASModel.H" -#include "EulerDdtScheme.H" +#include "turbulenceModel.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -49,7 +48,7 @@ int main(int argc, char *argv[]) #include "setRootCase.H" #include "createTime.H" #include "createDynamicFvMesh.H" - #include "readEnvironmentalProperties.H" + #include "readGravitationalAcceleration.H" #include "readPISOControls.H" #include "initContinuityErrs.H" #include "createFields.H" @@ -106,7 +105,7 @@ int main(int argc, char *argv[]) twoPhaseProperties.correct(); - #include "gammaEqnSubCycle.H" + #include "alphaEqnSubCycle.H" #include "UEqn.H" @@ -139,7 +138,7 @@ int main(int argc, char *argv[]) Info<< "End\n" << endl; - return(0); + return 0; } diff --git a/applications/solvers/multiphase/interDyMFoam/pEqn.H b/applications/solvers/multiphase/interDyMFoam/pEqn.H index 73c53b570..bda17ff52 100644 --- a/applications/solvers/multiphase/interDyMFoam/pEqn.H +++ b/applications/solvers/multiphase/interDyMFoam/pEqn.H @@ -5,17 +5,18 @@ U = rAU*UEqn.H(); surfaceScalarField phiU("phiU", (fvc::interpolate(U) & mesh.Sf())); - phi = phiU + - ( - fvc::interpolate(interface.sigmaK())*fvc::snGrad(gamma) - - ghf*fvc::snGrad(rho) - )*rAUf*mesh.magSf(); - if (pd.needReference()) { adjustPhi(phi, U, pd); } + phi = phiU + + ( + fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1) + - ghf*fvc::snGrad(rho) + )*rAUf*mesh.magSf(); + + for(int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) { fvScalarMatrix pdEqn diff --git a/applications/solvers/multiphase/interFoam/Make/options b/applications/solvers/multiphase/interFoam/Make/options index 376c09beb..cda832118 100644 --- a/applications/solvers/multiphase/interFoam/Make/options +++ b/applications/solvers/multiphase/interFoam/Make/options @@ -2,10 +2,13 @@ EXE_INC = \ -I$(LIB_SRC)/transportModels \ -I$(LIB_SRC)/transportModels/incompressible/lnInclude \ -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \ + -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \ -I$(LIB_SRC)/finiteVolume/lnInclude EXE_LIBS = \ -linterfaceProperties \ -lincompressibleTransportModels \ + -lincompressibleRASModels \ + -lincompressibleLESModels \ -lfiniteVolume \ -llduSolvers diff --git a/applications/solvers/multiphase/interFoam/UBlendingFactor.H b/applications/solvers/multiphase/interFoam/UBlendingFactor.H deleted file mode 100644 index 217a9c397..000000000 --- a/applications/solvers/multiphase/interFoam/UBlendingFactor.H +++ /dev/null @@ -1,6 +0,0 @@ - surfaceScalarField gammaf = fvc::interpolate(gamma); - surfaceScalarField UBlendingFactor - ( - "UBlendingFactor", - sqrt(max(min(4*gammaf*(1.0 - gammaf), 1.0), 0.0)) - ); diff --git a/applications/solvers/multiphase/interFoam/UEqn.H b/applications/solvers/multiphase/interFoam/UEqn.H index 7e4a29468..528e0aaaf 100644 --- a/applications/solvers/multiphase/interFoam/UEqn.H +++ b/applications/solvers/multiphase/interFoam/UEqn.H @@ -1,12 +1,17 @@ - surfaceScalarField muf = twoPhaseProperties.muf(); + surfaceScalarField muEff + ( + "muEff", + twoPhaseProperties.muf() + + fvc::interpolate(rho*turbulence->nut()) + ); fvVectorMatrix UEqn ( fvm::ddt(rho, U) + fvm::div(rhoPhi, U) - - fvm::laplacian(muf, U) - - (fvc::grad(U) & fvc::grad(muf)) - //- fvc::div(muf*(fvc::interpolate(dev(fvc::grad(U))) & mesh.Sf())) + - fvm::laplacian(muEff, U) + - (fvc::grad(U) & fvc::grad(muEff)) + //- fvc::div(muEff*(fvc::interpolate(dev(fvc::grad(U))) & mesh.Sf())) ); UEqn.relax(); @@ -20,7 +25,7 @@ fvc::reconstruct ( ( - fvc::interpolate(interface.sigmaK())*fvc::snGrad(gamma) + fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1) - ghf*fvc::snGrad(rho) - fvc::snGrad(pd) ) * mesh.magSf() diff --git a/applications/solvers/multiphase/interFoam/alphaEqn.H b/applications/solvers/multiphase/interFoam/alphaEqn.H new file mode 100644 index 000000000..0b2fb4ebf --- /dev/null +++ b/applications/solvers/multiphase/interFoam/alphaEqn.H @@ -0,0 +1,35 @@ +{ + word alphaScheme("div(phi,alpha)"); + word alpharScheme("div(phirb,alpha)"); + + surfaceScalarField phic = mag(phi/mesh.magSf()); + phic = min(interface.cAlpha()*phic, max(phic)); + surfaceScalarField phir = phic*interface.nHatf(); + + for (int aCorr=0; aCorr 1) +{ + dimensionedScalar totalDeltaT = runTime.deltaT(); + surfaceScalarField rhoPhiSum = 0.0*rhoPhi; + + for + ( + subCycle alphaSubCycle(alpha1, nAlphaSubCycles); + !(++alphaSubCycle).end(); + ) + { +# include "alphaEqn.H" + rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi; + } + + rhoPhi = rhoPhiSum; +} +else +{ +# include "alphaEqn.H" +} + +interface.correct(); + +rho == alpha1*rho1 + (scalar(1) - alpha1)*rho2; diff --git a/applications/solvers/multiphase/interFoam/correctPhi.H b/applications/solvers/multiphase/interFoam/correctPhi.H index 171e1670f..ed5ccc9a7 100644 --- a/applications/solvers/multiphase/interFoam/correctPhi.H +++ b/applications/solvers/multiphase/interFoam/correctPhi.H @@ -1,7 +1,11 @@ { # include "continuityErrs.H" - wordList pcorrTypes(pd.boundaryField().types()); + wordList pcorrTypes + ( + pd.boundaryField().size(), + zeroGradientFvPatchScalarField::typeName + ); for (label i=0; i turbulence + ( + incompressible::turbulenceModel::New(U, phi, twoPhaseProperties) + ); diff --git a/applications/solvers/multiphase/interFoam/gammaEqn.H b/applications/solvers/multiphase/interFoam/gammaEqn.H deleted file mode 100644 index 8978d1d29..000000000 --- a/applications/solvers/multiphase/interFoam/gammaEqn.H +++ /dev/null @@ -1,35 +0,0 @@ -{ - word gammaScheme("div(phi,gamma)"); - word gammarScheme("div(phirb,gamma)"); - - surfaceScalarField phic = mag(phi/mesh.magSf()); - phic = min(interface.cGamma()*phic, max(phic)); - surfaceScalarField phir = phic*interface.nHatf(); - - for (int gCorr=0; gCorr 1) -{ - dimensionedScalar totalDeltaT = runTime.deltaT(); - surfaceScalarField rhoPhiSum = 0.0*rhoPhi; - - for - ( - subCycle gammaSubCycle(gamma, nGammaSubCycles); - !(++gammaSubCycle).end(); - ) - { -# include "gammaEqn.H" - rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi; - } - - rhoPhi = rhoPhiSum; -} -else -{ -# include "gammaEqn.H" -} - -interface.correct(); - -rho == gamma*rho1 + (scalar(1) - gamma)*rho2; diff --git a/applications/solvers/multiphase/interFoam/interFoam.C b/applications/solvers/multiphase/interFoam/interFoam.C index 8f30c061c..16e4b0558 100644 --- a/applications/solvers/multiphase/interFoam/interFoam.C +++ b/applications/solvers/multiphase/interFoam/interFoam.C @@ -28,9 +28,12 @@ Application Description Solver for 2 incompressible, isothermal immiscible fluids using a VOF (volume of fluid) phase-fraction based interface capturing approach. + The momentum and other fluid properties are of the "mixture" and a single momentum equation is solved. + Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected. + For a two-fluid approach see twoPhaseEulerFoam. \*---------------------------------------------------------------------------*/ @@ -40,6 +43,7 @@ Description #include "subCycle.H" #include "interfaceProperties.H" #include "twoPhaseMixture.H" +#include "turbulenceModel.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -48,7 +52,7 @@ int main(int argc, char *argv[]) #include "setRootCase.H" #include "createTime.H" #include "createMesh.H" - #include "readEnvironmentalProperties.H" + #include "readGravitationalAcceleration.H" #include "readPISOControls.H" #include "initContinuityErrs.H" #include "createFields.H" @@ -74,7 +78,7 @@ int main(int argc, char *argv[]) twoPhaseProperties.correct(); - #include "gammaEqnSubCycle.H" + #include "alphaEqnSubCycle.H" #include "UEqn.H" @@ -107,7 +111,7 @@ int main(int argc, char *argv[]) Info<< "End\n" << endl; - return(0); + return 0; } diff --git a/applications/solvers/multiphase/interFoam/pEqn.H b/applications/solvers/multiphase/interFoam/pEqn.H index 77256a8e2..e93021a46 100644 --- a/applications/solvers/multiphase/interFoam/pEqn.H +++ b/applications/solvers/multiphase/interFoam/pEqn.H @@ -7,16 +7,18 @@ surfaceScalarField phiU ( "phiU", - (fvc::interpolate(U) & mesh.Sf()) + fvc::ddtPhiCorr(rUA, rho, U, phi) + (fvc::interpolate(U) & mesh.Sf()) + + fvc::ddtPhiCorr(rUA, rho, U, phi) ); + adjustPhi(phiU, U, p); + phi = phiU + ( - fvc::interpolate(interface.sigmaK())*fvc::snGrad(gamma) + fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1) - ghf*fvc::snGrad(rho) )*rUAf*mesh.magSf(); - adjustPhi(phi, U, pd); for(int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) { diff --git a/applications/solvers/multiphase/interMixingFoam/Make/files b/applications/solvers/multiphase/interMixingFoam/Make/files new file mode 100644 index 000000000..488cd77e5 --- /dev/null +++ b/applications/solvers/multiphase/interMixingFoam/Make/files @@ -0,0 +1,6 @@ +incompressibleThreePhaseMixture/threePhaseMixture.C +threePhaseInterfaceProperties/threePhaseInterfaceProperties.C +interMixingFoam.C + +EXE = $(FOAM_APPBIN)/interMixingFoam + diff --git a/applications/solvers/multiphase/interMixingFoam/Make/options b/applications/solvers/multiphase/interMixingFoam/Make/options new file mode 100644 index 000000000..d8e4da231 --- /dev/null +++ b/applications/solvers/multiphase/interMixingFoam/Make/options @@ -0,0 +1,17 @@ +INTERFOAM = $(FOAM_SOLVERS)/multiphase/interFoam + +EXE_INC = \ + -I$(INTERFOAM) \ + -IincompressibleThreePhaseMixture \ + -IthreePhaseInterfaceProperties \ + -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \ + -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \ + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/transportModels + +EXE_LIBS = \ + -linterfaceProperties \ + -lincompressibleTransportModels \ + -lincompressibleRASModels \ + -lincompressibleLESModels \ + -lfiniteVolume diff --git a/applications/solvers/multiphase/interMixingFoam/alphaEqns.H b/applications/solvers/multiphase/interMixingFoam/alphaEqns.H new file mode 100644 index 000000000..1fc545d45 --- /dev/null +++ b/applications/solvers/multiphase/interMixingFoam/alphaEqns.H @@ -0,0 +1,164 @@ +{ + word alphaScheme("div(phi,alpha)"); + word alpharScheme("div(phirb,alpha)"); + + surfaceScalarField phir + ( + IOobject + ( + "phir", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + interface.cAlpha()*mag(phi/mesh.magSf())*interface.nHatf() + ); + + for (int gCorr=0; gCorr(mesh, phi).flux(alpha1); + + // Calculate the flux correction for alpha1 + phiAlpha1 -= phiAlpha1BD; + + // Calculate the limiter for alpha1 + MULES::limiter + ( + allLambda, + geometricOneField(), + alpha1, + phiAlpha1BD, + phiAlpha1, + zeroField(), + zeroField(), + 1, + 0, + 3 + ); + + // Create the complete flux for alpha2 + surfaceScalarField phiAlpha2 = + fvc::flux + ( + phi, + alpha2, + alphaScheme + ) + + fvc::flux + ( + -fvc::flux(phir, alpha1, alpharScheme), + alpha2, + alpharScheme + ); + + // Create the bounded (upwind) flux for alpha2 + surfaceScalarField phiAlpha2BD = + upwind(mesh, phi).flux(alpha2); + + // Calculate the flux correction for alpha2 + phiAlpha2 -= phiAlpha2BD; + + // Further limit the limiter for alpha2 + MULES::limiter + ( + allLambda, + geometricOneField(), + alpha2, + phiAlpha2BD, + phiAlpha2, + zeroField(), + zeroField(), + 1, + 0, + 3 + ); + + // Construct the limited fluxes + phiAlpha1 = phiAlpha1BD + lambda*phiAlpha1; + phiAlpha2 = phiAlpha2BD + lambda*phiAlpha2; + + // Solve for alpha1 + solve(fvm::ddt(alpha1) + fvc::div(phiAlpha1)); + + // Create the diffusion coefficients for alpha2<->alpha3 + volScalarField Dc23 = D23*max(alpha3, scalar(0))*pos(alpha2); + volScalarField Dc32 = D23*max(alpha2, scalar(0))*pos(alpha3); + + // Add the diffusive flux for alpha3->alpha2 + phiAlpha2 -= fvc::interpolate(Dc32)*mesh.magSf()*fvc::snGrad(alpha1); + + // Solve for alpha2 + fvScalarMatrix alpha2Eqn + ( + fvm::ddt(alpha2) + + fvc::div(phiAlpha2) + - fvm::laplacian(Dc23 + Dc32, alpha2) + ); + alpha2Eqn.solve(); + + // Construct the complete mass flux + rhoPhi = + phiAlpha1*(rho1 - rho3) + + (phiAlpha2 + alpha2Eqn.flux())*(rho2 - rho3) + + phi*rho3; + + alpha3 = 1.0 - alpha1 - alpha2; + } + + Info<< "Air phase volume fraction = " + << alpha1.weightedAverage(mesh.V()).value() + << " Min(alpha1) = " << min(alpha1).value() + << " Max(alpha1) = " << max(alpha1).value() + << endl; + + Info<< "Liquid phase volume fraction = " + << alpha2.weightedAverage(mesh.V()).value() + << " Min(alpha2) = " << min(alpha2).value() + << " Max(alpha2) = " << max(alpha2).value() + << endl; +} diff --git a/applications/solvers/multiphase/interMixingFoam/alphaEqnsSubCycle.H b/applications/solvers/multiphase/interMixingFoam/alphaEqnsSubCycle.H new file mode 100644 index 000000000..765087a18 --- /dev/null +++ b/applications/solvers/multiphase/interMixingFoam/alphaEqnsSubCycle.H @@ -0,0 +1,43 @@ +label nAlphaCorr +( + readLabel(piso.lookup("nAlphaCorr")) +); + +label nAlphaSubCycles +( + readLabel(piso.lookup("nAlphaSubCycles")) +); + +if (nAlphaSubCycles > 1) +{ + surfaceScalarField rhoPhiSum = 0.0*rhoPhi; + dimensionedScalar totalDeltaT = runTime.deltaT(); + + for + ( + subCycle alphaSubCycle(alpha1, nAlphaSubCycles); + !(++alphaSubCycle).end(); + ) + { +# include "alphaEqns.H" + rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi; + } + + rhoPhi = rhoPhiSum; +} +else +{ +# include "alphaEqns.H" +} + +interface.correct(); + +{ + volScalarField rhoNew = alpha1*rho1 + alpha2*rho2 + alpha3*rho3; + + //solve(fvm::ddt(rho) + fvc::div(rhoPhi)); + //Info<< "density error = " + // << max((mag(rho - rhoNew)/mag(rhoNew))().internalField()) << endl; + + rho == rhoNew; +} diff --git a/applications/solvers/multiphase/lesInterFoam/createFields.H b/applications/solvers/multiphase/interMixingFoam/createFields.H similarity index 58% rename from applications/solvers/multiphase/lesInterFoam/createFields.H rename to applications/solvers/multiphase/interMixingFoam/createFields.H index 2b858bd30..6fdd34175 100644 --- a/applications/solvers/multiphase/lesInterFoam/createFields.H +++ b/applications/solvers/multiphase/interMixingFoam/createFields.H @@ -1,4 +1,4 @@ - Info<< "Reading field pd\n" << endl; + Info<< "Reading field p\n" << endl; volScalarField pd ( IOobject @@ -12,12 +12,12 @@ mesh ); - Info<< "Reading field gamma\n" << endl; - volScalarField gamma + Info<< "Reading field alpha1\n" << endl; + volScalarField alpha1 ( IOobject ( - "gamma", + "alpha1", runTime.timeName(), mesh, IOobject::MUST_READ, @@ -26,6 +26,39 @@ mesh ); + + Info<< "Reading field alpha2\n" << endl; + volScalarField alpha2 + ( + IOobject + ( + "alpha2", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh + ); + + + Info<< "Reading field alpha3\n" << endl; + volScalarField alpha3 + ( + IOobject + ( + "alpha3", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh + ); + + alpha3 == 1.0 - alpha1 - alpha2; + + Info<< "Reading field U\n" << endl; volVectorField U ( @@ -42,12 +75,13 @@ # include "createPhi.H" - Info<< "Reading transportProperties\n" << endl; - twoPhaseMixture twoPhaseProperties(U, phi, "gamma"); + threePhaseMixture threePhaseProperties(U, phi); - const dimensionedScalar& rho1 = twoPhaseProperties.rho1(); - const dimensionedScalar& rho2 = twoPhaseProperties.rho2(); + const dimensionedScalar& rho1 = threePhaseProperties.rho1(); + const dimensionedScalar& rho2 = threePhaseProperties.rho2(); + const dimensionedScalar& rho3 = threePhaseProperties.rho3(); + dimensionedScalar D23(threePhaseProperties.lookup("D23")); // Need to store rho for ddt(rho, U) volScalarField rho @@ -59,14 +93,15 @@ mesh, IOobject::READ_IF_PRESENT ), - gamma*rho1 + (scalar(1) - gamma)*rho2, - gamma.boundaryField().types() + alpha1*rho1 + alpha2*rho2 + alpha3*rho3, + alpha1.boundaryField().types() ); rho.oldTime(); + // Mass flux // Initialisation does not matter because rhoPhi is reset after the - // gamma solution before it is used in the U equation. + // alpha solution before it is used in the U equation. surfaceScalarField rhoPhi ( IOobject @@ -80,7 +115,6 @@ rho1*phi ); - Info<< "Calculating field g.h\n" << endl; volScalarField gh("gh", g & mesh.C()); surfaceScalarField ghf("gh", g & mesh.Cf()); @@ -99,7 +133,6 @@ pd + rho*gh ); - label pdRefCell = 0; scalar pdRefValue = 0.0; setRefCell(pd, mesh.solutionDict().subDict("PISO"), pdRefCell, pdRefValue); @@ -121,11 +154,12 @@ ); } - // Construct interface from gamma distribution - interfaceProperties interface(gamma, U, twoPhaseProperties); + // Construct interface from alpha distribution + threePhaseInterfaceProperties interface(threePhaseProperties); - // Construct LES model - autoPtr turbulence + + // Construct incompressible turbulence model + autoPtr turbulence ( - incompressible::LESModel::New(U, phi, twoPhaseProperties) + incompressible::turbulenceModel::New(U, phi, threePhaseProperties) ); diff --git a/applications/solvers/multiphase/interMixingFoam/incompressibleThreePhaseMixture/threePhaseMixture.C b/applications/solvers/multiphase/interMixingFoam/incompressibleThreePhaseMixture/threePhaseMixture.C new file mode 100644 index 000000000..42c7769bc --- /dev/null +++ b/applications/solvers/multiphase/interMixingFoam/incompressibleThreePhaseMixture/threePhaseMixture.C @@ -0,0 +1,204 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright held by original author + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM 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 2 of the License, or (at your + option) any later version. + + OpenFOAM 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 OpenFOAM; if not, write to the Free Software Foundation, + Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + +Class + threePhaseMixture + +\*---------------------------------------------------------------------------*/ + +#include "threePhaseMixture.H" +#include "addToRunTimeSelectionTable.H" +#include "surfaceFields.H" +#include "fvc.H" + +// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * // + +//- Calculate and return the laminar viscosity +void Foam::threePhaseMixture::calcNu() +{ + // Average kinematic viscosity calculated from dynamic viscosity + nu_ = mu()/(alpha1_*rho1_ + alpha2_*rho2_ + alpha3_*rho3_); +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::threePhaseMixture::threePhaseMixture +( + const volVectorField& U, + const surfaceScalarField& phi +) +: + transportModel(U, phi), + + phase1Name_("phase1"), + phase2Name_("phase2"), + phase3Name_("phase3"), + + nuModel1_ + ( + viscosityModel::New + ( + "nu1", + subDict(phase1Name_), + U, + phi + ) + ), + nuModel2_ + ( + viscosityModel::New + ( + "nu2", + subDict(phase2Name_), + U, + phi + ) + ), + nuModel3_ + ( + viscosityModel::New + ( + "nu3", + subDict(phase2Name_), + U, + phi + ) + ), + + rho1_(nuModel1_->viscosityProperties().lookup("rho")), + rho2_(nuModel2_->viscosityProperties().lookup("rho")), + rho3_(nuModel3_->viscosityProperties().lookup("rho")), + + U_(U), + phi_(phi), + + alpha1_(U_.db().lookupObject ("alpha1")), + alpha2_(U_.db().lookupObject ("alpha2")), + alpha3_(U_.db().lookupObject ("alpha3")), + + nu_ + ( + IOobject + ( + "nu", + U_.time().timeName(), + U_.db() + ), + U_.mesh(), + dimensionedScalar("nu", dimensionSet(0, 2, -1, 0, 0), 0), + calculatedFvPatchScalarField::typeName + ) +{ + calcNu(); +} + + +// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * // + +Foam::tmp Foam::threePhaseMixture::mu() const +{ + return tmp + ( + new volScalarField + ( + "mu", + alpha1_*rho1_*nuModel1_->nu() + + alpha2_*rho2_*nuModel2_->nu() + + alpha3_*rho3_*nuModel3_->nu() + ) + ); +} + + +Foam::tmp Foam::threePhaseMixture::muf() const +{ + surfaceScalarField alpha1f = fvc::interpolate(alpha1_); + surfaceScalarField alpha2f = fvc::interpolate(alpha2_); + surfaceScalarField alpha3f = fvc::interpolate(alpha3_); + + return tmp + ( + new surfaceScalarField + ( + "mu", + alpha1f*rho1_*fvc::interpolate(nuModel1_->nu()) + + alpha2f*rho2_*fvc::interpolate(nuModel2_->nu()) + + alpha3f*rho3_*fvc::interpolate(nuModel3_->nu()) + ) + ); +} + + +Foam::tmp Foam::threePhaseMixture::nuf() const +{ + surfaceScalarField alpha1f = fvc::interpolate(alpha1_); + surfaceScalarField alpha2f = fvc::interpolate(alpha2_); + surfaceScalarField alpha3f = fvc::interpolate(alpha3_); + + return tmp + ( + new surfaceScalarField + ( + "nu", + ( + alpha1f*rho1_*fvc::interpolate(nuModel1_->nu()) + + alpha2f*rho2_*fvc::interpolate(nuModel2_->nu()) + + alpha3f*rho3_*fvc::interpolate(nuModel3_->nu()) + )/(alpha1f*rho1_ + alpha2f*rho2_ + alpha3f*rho3_) + ) + ); +} + + +bool Foam::threePhaseMixture::read() +{ + if (transportModel::read()) + { + if + ( + nuModel1_().read(*this) + && nuModel2_().read(*this) + && nuModel3_().read(*this) + ) + { + nuModel1_->viscosityProperties().lookup("rho") >> rho1_; + nuModel2_->viscosityProperties().lookup("rho") >> rho2_; + nuModel3_->viscosityProperties().lookup("rho") >> rho3_; + + return true; + } + else + { + return false; + } + } + else + { + return false; + } +} + + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/interMixingFoam/incompressibleThreePhaseMixture/threePhaseMixture.H b/applications/solvers/multiphase/interMixingFoam/incompressibleThreePhaseMixture/threePhaseMixture.H new file mode 100644 index 000000000..b55b75a2c --- /dev/null +++ b/applications/solvers/multiphase/interMixingFoam/incompressibleThreePhaseMixture/threePhaseMixture.H @@ -0,0 +1,197 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright held by original author + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM 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 2 of the License, or (at your + option) any later version. + + OpenFOAM 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 OpenFOAM; if not, write to the Free Software Foundation, + Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + +Class + threePhaseMixture + +Description + +SourceFiles + threePhaseMixture.C + +\*---------------------------------------------------------------------------*/ + +#ifndef threePhaseMixture_H +#define threePhaseMixture_H + +#include "incompressible/transportModel/transportModel.H" +#include "incompressible/viscosityModels/viscosityModel/viscosityModel.H" +#include "dimensionedScalar.H" +#include "volFields.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class threePhaseMixture Declaration +\*---------------------------------------------------------------------------*/ + +class threePhaseMixture +: + public transportModel +{ + // Private data + + word phase1Name_; + word phase2Name_; + word phase3Name_; + + autoPtr nuModel1_; + autoPtr nuModel2_; + autoPtr nuModel3_; + + dimensionedScalar rho1_; + dimensionedScalar rho2_; + dimensionedScalar rho3_; + + const volVectorField& U_; + const surfaceScalarField& phi_; + + const volScalarField& alpha1_; + const volScalarField& alpha2_; + const volScalarField& alpha3_; + + volScalarField nu_; + + + // Private Member Functions + + //- Calculate and return the laminar viscosity + void calcNu(); + + +public: + + // Constructors + + //- Construct from components + threePhaseMixture + ( + const volVectorField& U, + const surfaceScalarField& phi + ); + + + // Destructor + + ~threePhaseMixture() + {} + + + // Member Functions + + //- Return const-access to phase1 viscosityModel + const viscosityModel& nuModel1() const + { + return nuModel1_(); + } + + //- Return const-access to phase2 viscosityModel + const viscosityModel& nuModel2() const + { + return nuModel2_(); + } + + //- Return const-access to phase3 viscosityModel + const viscosityModel& nuModel3() const + { + return nuModel3_(); + } + + //- Return const-access to phase1 density + const dimensionedScalar& rho1() const + { + return rho1_; + } + + //- Return const-access to phase2 density + const dimensionedScalar& rho2() const + { + return rho2_; + }; + + //- Return const-access to phase3 density + const dimensionedScalar& rho3() const + { + return rho3_; + }; + + const volScalarField& alpha1() const + { + return alpha1_; + } + + const volScalarField& alpha2() const + { + return alpha2_; + } + + const volScalarField& alpha3() const + { + return alpha3_; + } + + //- Return the velocity + const volVectorField& U() const + { + return U_; + } + + //- Return the dynamic laminar viscosity + tmp mu() const; + + //- Return the face-interpolated dynamic laminar viscosity + tmp muf() const; + + //- Return the kinematic laminar viscosity + const volScalarField& nu() const + { + return nu_; + } + + //- Return the face-interpolated dynamic laminar viscosity + tmp nuf() const; + + //- Correct the laminar viscosity + void correct() + { + calcNu(); + } + + //- Read base transportProperties dictionary + bool read(); +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/lesInterFoam/lesInterFoam.C b/applications/solvers/multiphase/interMixingFoam/interMixingFoam.C similarity index 70% rename from applications/solvers/multiphase/lesInterFoam/lesInterFoam.C rename to applications/solvers/multiphase/interMixingFoam/interMixingFoam.C index 0c0b0e698..3eaf09aa7 100644 --- a/applications/solvers/multiphase/lesInterFoam/lesInterFoam.C +++ b/applications/solvers/multiphase/interMixingFoam/interMixingFoam.C @@ -20,26 +20,23 @@ License You should have received a copy of the GNU General Public License along with OpenFOAM; if not, write to the Free Software Foundation, - Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA Application - lesInterFoam + interMixingFoam Description - Solver for 2 incompressible, isothermal immiscible fluids using a VOF - (volume of fluid) phase-fraction based interface capturing approach. - The momentum and other fluid properties are of the "mixture" and a single - momentum equation is solved. Turbulence is modelled using a run-time - selectable incompressible LES model. + Solver for 3 incompressible fluids, two of which are miscible, + using a VOF method to capture the interface. \*---------------------------------------------------------------------------*/ #include "fvCFD.H" #include "MULES.H" #include "subCycle.H" -#include "interfaceProperties.H" -#include "twoPhaseMixture.H" -#include "incompressible/LESModel/LESModel.H" +#include "threePhaseMixture.H" +#include "threePhaseInterfaceProperties.H" +#include "turbulenceModel.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -48,16 +45,16 @@ int main(int argc, char *argv[]) #include "setRootCase.H" #include "createTime.H" #include "createMesh.H" - #include "readEnvironmentalProperties.H" + #include "readGravitationalAcceleration.H" #include "readPISOControls.H" #include "initContinuityErrs.H" #include "createFields.H" #include "readTimeControls.H" - #include "correctPhi.H" #include "CourantNo.H" #include "setInitialDeltaT.H" + #include "correctPhi.H" -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // Info<< "\nStarting time loop\n" << endl; @@ -69,33 +66,25 @@ int main(int argc, char *argv[]) #include "setDeltaT.H" runTime++; + Info<< "Time = " << runTime.timeName() << nl << endl; - #include "gammaEqnSubCycle.H" + threePhaseProperties.correct(); - turbulence->correct(); + #include "alphaEqnsSubCycle.H" + #define twoPhaseProperties threePhaseProperties #include "UEqn.H" // --- PISO loop - for (int corr=0; corr < nCorr; corr++) + for (int corr=0; corrcorrect(); runTime.write(); @@ -104,7 +93,7 @@ int main(int argc, char *argv[]) << nl << endl; } - Info<< "End\n" << endl; + Info<< "\n end \n"; return(0); } diff --git a/applications/solvers/multiphase/interMixingFoam/threePhaseInterfaceProperties/threePhaseInterfaceProperties.C b/applications/solvers/multiphase/interMixingFoam/threePhaseInterfaceProperties/threePhaseInterfaceProperties.C new file mode 100644 index 000000000..0ea347f0d --- /dev/null +++ b/applications/solvers/multiphase/interMixingFoam/threePhaseInterfaceProperties/threePhaseInterfaceProperties.C @@ -0,0 +1,209 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright held by original author + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM 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 2 of the License, or (at your + option) any later version. + + OpenFOAM 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 OpenFOAM; if not, write to the Free Software Foundation, + Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + +Application + threePhaseInterfaceProperties + +Description + Properties to aid interFoam : + 1. Correct the alpha boundary condition for dynamic contact angle. + 2. Calculate interface curvature. + +\*---------------------------------------------------------------------------*/ + +#include "threePhaseInterfaceProperties.H" +#include "alphaContactAngleFvPatchScalarField.H" +#include "mathematicalConstants.H" +#include "surfaceInterpolate.H" +#include "fvcDiv.H" +#include "fvcGrad.H" + +// * * * * * * * * * * * * * * * Static Member Data * * * * * * * * * * * * // + +const Foam::scalar Foam::threePhaseInterfaceProperties::convertToRad = + Foam::mathematicalConstant::pi/180.0; + + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +// Correction for the boundary condition on the unit normal nHat on +// walls to produce the correct contact angle. + +// The dynamic contact angle is calculated from the component of the +// velocity on the direction of the interface, parallel to the wall. + +void Foam::threePhaseInterfaceProperties::correctContactAngle +( + surfaceVectorField::GeometricBoundaryField& nHatb +) const +{ + const volScalarField::GeometricBoundaryField& alpha1 = + mixture_.alpha1().boundaryField(); + const volScalarField::GeometricBoundaryField& alpha2 = + mixture_.alpha2().boundaryField(); + const volScalarField::GeometricBoundaryField& alpha3 = + mixture_.alpha3().boundaryField(); + const volVectorField::GeometricBoundaryField& U = + mixture_.U().boundaryField(); + + const fvMesh& mesh = mixture_.U().mesh(); + const fvBoundaryMesh& boundary = mesh.boundary(); + + forAll(boundary, patchi) + { + if (isA(alpha1[patchi])) + { + const alphaContactAngleFvPatchScalarField& a2cap = + refCast + (alpha2[patchi]); + + const alphaContactAngleFvPatchScalarField& a3cap = + refCast + (alpha3[patchi]); + + scalarField twoPhaseAlpha2 = max(a2cap, scalar(0)); + scalarField twoPhaseAlpha3 = max(a3cap, scalar(0)); + + scalarField sumTwoPhaseAlpha = + twoPhaseAlpha2 + twoPhaseAlpha3 + SMALL; + + twoPhaseAlpha2 /= sumTwoPhaseAlpha; + twoPhaseAlpha3 /= sumTwoPhaseAlpha; + + fvsPatchVectorField& nHatp = nHatb[patchi]; + + scalarField theta = + convertToRad + *( + twoPhaseAlpha2*(180 - a2cap.theta(U[patchi], nHatp)) + + twoPhaseAlpha3*(180 - a3cap.theta(U[patchi], nHatp)) + ); + + vectorField nf = boundary[patchi].nf(); + + // Reset nHatPatch to correspond to the contact angle + + scalarField a12 = nHatp & nf; + + scalarField b1 = cos(theta); + + scalarField b2(nHatp.size()); + + forAll(b2, facei) + { + b2[facei] = cos(acos(a12[facei]) - theta[facei]); + } + + scalarField det = 1.0 - a12*a12; + + scalarField a = (b1 - a12*b2)/det; + scalarField b = (b2 - a12*b1)/det; + + nHatp = a*nf + b*nHatp; + + nHatp /= (mag(nHatp) + deltaN_.value()); + } + } +} + + +void Foam::threePhaseInterfaceProperties::calculateK() +{ + const volScalarField& alpha1 = mixture_.alpha1(); + + const fvMesh& mesh = alpha1.mesh(); + const surfaceVectorField& Sf = mesh.Sf(); + + // Cell gradient of alpha + volVectorField gradAlpha = fvc::grad(alpha1); + + // Interpolated face-gradient of alpha + surfaceVectorField gradAlphaf = fvc::interpolate(gradAlpha); + + // Face unit interface normal + surfaceVectorField nHatfv = gradAlphaf/(mag(gradAlphaf) + deltaN_); + correctContactAngle(nHatfv.boundaryField()); + + // Face unit interface normal flux + nHatf_ = nHatfv & Sf; + + // Simple expression for curvature + K_ = -fvc::div(nHatf_); + + // Complex expression for curvature. + // Correction is formally zero but numerically non-zero. + //volVectorField nHat = gradAlpha/(mag(gradAlpha) + deltaN_); + //nHat.boundaryField() = nHatfv.boundaryField(); + //K_ = -fvc::div(nHatf_) + (nHat & fvc::grad(nHatfv) & nHat); +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::threePhaseInterfaceProperties::threePhaseInterfaceProperties +( + const threePhaseMixture& mixture +) +: + mixture_(mixture), + cAlpha_ + ( + readScalar + ( + mixture.U().mesh().solutionDict().subDict("PISO").lookup("cAlpha") + ) + ), + sigma12_(mixture.lookup("sigma12")), + sigma13_(mixture.lookup("sigma13")), + + deltaN_ + ( + "deltaN", + 1e-8/pow(average(mixture.U().mesh().V()), 1.0/3.0) + ), + + nHatf_ + ( + ( + fvc::interpolate(fvc::grad(mixture.alpha1())) + /(mag(fvc::interpolate(fvc::grad(mixture.alpha1()))) + deltaN_) + ) & mixture.alpha1().mesh().Sf() + ), + + K_ + ( + IOobject + ( + "K", + mixture.alpha1().time().timeName(), + mixture.alpha1().mesh() + ), + -fvc::div(nHatf_) + ) +{ + calculateK(); +} + + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/interMixingFoam/threePhaseInterfaceProperties/threePhaseInterfaceProperties.H b/applications/solvers/multiphase/interMixingFoam/threePhaseInterfaceProperties/threePhaseInterfaceProperties.H new file mode 100644 index 000000000..b5594b6e9 --- /dev/null +++ b/applications/solvers/multiphase/interMixingFoam/threePhaseInterfaceProperties/threePhaseInterfaceProperties.H @@ -0,0 +1,157 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright held by original author + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM 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 2 of the License, or (at your + option) any later version. + + OpenFOAM 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 OpenFOAM; if not, write to the Free Software Foundation, + Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + +Class + threePhaseInterfaceProperties + +Description + Properties to aid interFoam : + 1. Correct the alpha boundary condition for dynamic contact angle. + 2. Calculate interface curvature. + +SourceFiles + threePhaseInterfaceProperties.C + +\*---------------------------------------------------------------------------*/ + +#ifndef threePhaseInterfaceProperties_H +#define threePhaseInterfaceProperties_H + +#include "threePhaseMixture.H" +#include "surfaceFields.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class threePhaseInterfaceProperties Declaration +\*---------------------------------------------------------------------------*/ + +class threePhaseInterfaceProperties +{ + // Private data + + const threePhaseMixture& mixture_; + + //- Compression coefficient + scalar cAlpha_; + + //- Surface tension 1-2 + dimensionedScalar sigma12_; + + //- Surface tension 1-3 + dimensionedScalar sigma13_; + + //- Stabilisation for normalisation of the interface normal + const dimensionedScalar deltaN_; + + surfaceScalarField nHatf_; + volScalarField K_; + + + // Private Member Functions + + //- Disallow default bitwise copy construct and assignment + threePhaseInterfaceProperties(const threePhaseInterfaceProperties&); + void operator=(const threePhaseInterfaceProperties&); + + //- Correction for the boundary condition on the unit normal nHat on + // walls to produce the correct contact dynamic angle + // calculated from the component of U parallel to the wall + void correctContactAngle + ( + surfaceVectorField::GeometricBoundaryField& nHat + ) const; + + //- Re-calculate the interface curvature + void calculateK(); + + +public: + + //- Conversion factor for degrees into radians + static const scalar convertToRad; + + + // Constructors + + //- Construct from volume fraction field alpha and IOdictionary + threePhaseInterfaceProperties(const threePhaseMixture& mixture); + + + // Member Functions + + scalar cAlpha() const + { + return cAlpha_; + } + + const dimensionedScalar& deltaN() const + { + return deltaN_; + } + + const surfaceScalarField& nHatf() const + { + return nHatf_; + } + + const volScalarField& K() const + { + return K_; + } + + tmp sigma() const + { + volScalarField limitedAlpha2 = max(mixture_.alpha2(), scalar(0)); + volScalarField limitedAlpha3 = max(mixture_.alpha3(), scalar(0)); + + return + (limitedAlpha2*sigma12_ + limitedAlpha3*sigma13_) + /(limitedAlpha2 + limitedAlpha3 + SMALL); + } + + tmp sigmaK() const + { + return sigma()*K_; + } + + void correct() + { + calculateK(); + } +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/Make/options b/applications/solvers/multiphase/interPhaseChangeFoam/Make/options index fa6943e58..c427781c7 100644 --- a/applications/solvers/multiphase/interPhaseChangeFoam/Make/options +++ b/applications/solvers/multiphase/interPhaseChangeFoam/Make/options @@ -2,13 +2,14 @@ EXE_INC = \ -I$(LIB_SRC)/transportModels \ -I$(LIB_SRC)/transportModels/incompressible/lnInclude \ -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \ - -I$(LIB_SRC)/turbulenceModels/LES \ - -I$(LIB_SRC)/turbulenceModels/LES/LESdeltas/lnInclude \ + -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \ -IphaseChangeTwoPhaseMixtures/phaseChangeTwoPhaseMixture \ -I$(LIB_SRC)/finiteVolume/lnInclude EXE_LIBS = \ -linterfaceProperties \ -lincompressibleTransportModels \ + -lincompressibleRASModels \ -lincompressibleLESModels \ - -lfiniteVolume + -lfiniteVolume \ + -llduSolvers diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/UEqn.H b/applications/solvers/multiphase/interPhaseChangeFoam/UEqn.H index fefdeb41c..c59137c7b 100644 --- a/applications/solvers/multiphase/interPhaseChangeFoam/UEqn.H +++ b/applications/solvers/multiphase/interPhaseChangeFoam/UEqn.H @@ -1,15 +1,18 @@ - surfaceScalarField muf = + surfaceScalarField muEff + ( + "muEff", twoPhaseProperties->muf() - + fvc::interpolate(rho*turbulence->nuSgs()); + + fvc::interpolate(rho*turbulence->nut()) + ); fvVectorMatrix UEqn ( fvm::ddt(rho, U) + fvm::div(rhoPhi, U) - fvm::Sp(fvc::ddt(rho) + fvc::div(rhoPhi), U) - - fvm::laplacian(muf, U) - - (fvc::grad(U) & fvc::grad(muf)) - //- fvc::div(muf*(fvc::interpolate(dev2(fvc::grad(U))) & mesh.Sf())) + - fvm::laplacian(muEff, U) + - (fvc::grad(U) & fvc::grad(muEff)) + //- fvc::div(muEff*(fvc::interpolate(dev2(fvc::grad(U))) & mesh.Sf())) ); UEqn.relax(); @@ -23,7 +26,7 @@ fvc::reconstruct ( ( - fvc::interpolate(interface.sigmaK())*fvc::snGrad(gamma) + fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1) - ghf*fvc::snGrad(rho) - fvc::snGrad(pd) ) * mesh.magSf() diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/gammaEqn.H b/applications/solvers/multiphase/interPhaseChangeFoam/alphaEqn.H similarity index 57% rename from applications/solvers/multiphase/interPhaseChangeFoam/gammaEqn.H rename to applications/solvers/multiphase/interPhaseChangeFoam/alphaEqn.H index 018694892..a232056b8 100644 --- a/applications/solvers/multiphase/interPhaseChangeFoam/gammaEqn.H +++ b/applications/solvers/multiphase/interPhaseChangeFoam/alphaEqn.H @@ -1,23 +1,23 @@ { - word gammaScheme("div(phi,gamma)"); - word gammarScheme("div(phirb,gamma)"); + word alphaScheme("div(phi,alpha)"); + word alpharScheme("div(phirb,alpha)"); surfaceScalarField phir("phir", phic*interface.nHatf()); - for (int gCorr=0; gCorr > vDotAlphal = @@ -46,22 +46,22 @@ ), // Divergence term is handled explicitly to be // consistent with the explicit transport solution - divU*gamma + divU*alpha1 + vDotcAlphal ); - //MULES::explicitSolve(gamma, phi, phiGamma, 1, 0); - //MULES::explicitSolve(oneField(), gamma, phi, phiGamma, Sp, Su, 1, 0); - MULES::implicitSolve(oneField(), gamma, phi, phiGamma, Sp, Su, 1, 0); + //MULES::explicitSolve(alpha1, phi, phiAlpha, 1, 0); + //MULES::explicitSolve(geometricOneField(), alpha1, phi, phiAlpha, Sp, Su, 1, 0); + MULES::implicitSolve(geometricOneField(), alpha1, phi, phiAlpha, Sp, Su, 1, 0); rhoPhi += (runTime.deltaT()/totalDeltaT) - *(phiGamma*(rho1 - rho2) + phi*rho2); + *(phiAlpha*(rho1 - rho2) + phi*rho2); } Info<< "Liquid phase volume fraction = " - << gamma.weightedAverage(mesh.V()).value() - << " Min(gamma) = " << min(gamma).value() - << " Max(gamma) = " << max(gamma).value() + << alpha1.weightedAverage(mesh.V()).value() + << " Min(alpha1) = " << min(alpha1).value() + << " Max(alpha1) = " << max(alpha1).value() << endl; } diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/gammaEqnSubCycle.H b/applications/solvers/multiphase/interPhaseChangeFoam/alphaEqnSubCycle.H similarity index 53% rename from applications/solvers/multiphase/interPhaseChangeFoam/gammaEqnSubCycle.H rename to applications/solvers/multiphase/interPhaseChangeFoam/alphaEqnSubCycle.H index c6355d6fb..dd1d82803 100644 --- a/applications/solvers/multiphase/interPhaseChangeFoam/gammaEqnSubCycle.H +++ b/applications/solvers/multiphase/interPhaseChangeFoam/alphaEqnSubCycle.H @@ -11,37 +11,37 @@ surfaceScalarField rhoPhi ); { - label nGammaCorr + label nAlphaCorr ( - readLabel(piso.lookup("nGammaCorr")) + readLabel(piso.lookup("nAlphaCorr")) ); - label nGammaSubCycles + label nAlphaSubCycles ( - readLabel(piso.lookup("nGammaSubCycles")) + readLabel(piso.lookup("nAlphaSubCycles")) ); surfaceScalarField phic = mag(phi/mesh.magSf()); - phic = min(interface.cGamma()*phic, max(phic)); + phic = min(interface.cAlpha()*phic, max(phic)); volScalarField divU = fvc::div(phi); dimensionedScalar totalDeltaT = runTime.deltaT(); - if (nGammaSubCycles > 1) + if (nAlphaSubCycles > 1) { for ( - subCycle gammaSubCycle(gamma, nGammaSubCycles); - !(++gammaSubCycle).end(); + subCycle alphaSubCycle(alpha1, nAlphaSubCycles); + !(++alphaSubCycle).end(); ) { -# include "gammaEqn.H" +# include "alphaEqn.H" } } else { -# include "gammaEqn.H" +# include "alphaEqn.H" } if (nOuterCorr == 1) @@ -49,5 +49,5 @@ surfaceScalarField rhoPhi interface.correct(); } - rho == gamma*rho1 + (scalar(1) - gamma)*rho2; + rho == alpha1*rho1 + (scalar(1) - alpha1)*rho2; } diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/correctPhi.H b/applications/solvers/multiphase/interPhaseChangeFoam/correctPhi.H index 171e1670f..ed5ccc9a7 100644 --- a/applications/solvers/multiphase/interPhaseChangeFoam/correctPhi.H +++ b/applications/solvers/multiphase/interPhaseChangeFoam/correctPhi.H @@ -1,7 +1,11 @@ { # include "continuityErrs.H" - wordList pcorrTypes(pd.boundaryField().types()); + wordList pcorrTypes + ( + pd.boundaryField().size(), + zeroGradientFvPatchScalarField::typeName + ); for (label i=0; i twoPhaseProperties = - phaseChangeTwoPhaseMixture::New(U, phi, "gamma"); + phaseChangeTwoPhaseMixture::New(U, phi, "alpha1"); const dimensionedScalar& rho1 = twoPhaseProperties->rho1(); const dimensionedScalar& rho2 = twoPhaseProperties->rho2(); @@ -60,19 +60,32 @@ mesh, IOobject::READ_IF_PRESENT ), - gamma*rho1 + (scalar(1) - gamma)*rho2, - gamma.boundaryField().types() + alpha1*rho1 + (scalar(1) - alpha1)*rho2, + alpha1.boundaryField().types() ); rho.oldTime(); - label pdRefCell = 0; - scalar pdRefValue = 0.0; - setRefCell(pd, mesh.solutionDict().subDict("PISO"), pdRefCell, pdRefValue); + // Mass flux + // Initialisation does not matter because rhoPhi is reset after the + // alpha1 solution before it is used in the U equation. + surfaceScalarField rhoPhi + ( + IOobject + ( + "rho*phi", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + rho1*phi + ); - Info<< "Calculating field g.h" << endl; + + Info<< "Calculating field g.h\n" << endl; volScalarField gh("gh", g & mesh.C()); - surfaceScalarField ghf("ghf", g & mesh.Cf()); + surfaceScalarField ghf("gh", g & mesh.Cf()); volScalarField p ( @@ -88,11 +101,32 @@ ); - // Construct interface from gamma distribution - interfaceProperties interface(gamma, U, twoPhaseProperties()); + label pdRefCell = 0; + scalar pdRefValue = 0.0; + setRefCell(pd, mesh.solutionDict().subDict("PISO"), pdRefCell, pdRefValue); - // Construct LES model - autoPtr turbulence + scalar pRefValue = 0.0; + + if (pd.needReference()) + { + pRefValue = readScalar + ( + mesh.solutionDict().subDict("PISO").lookup("pRefValue") + ); + + p += dimensionedScalar + ( + "p", + p.dimensions(), + pRefValue - getRefCellValue(p, pdRefCell) + ); + } + + // Construct interface from alpha1 distribution + interfaceProperties interface(alpha1, U, twoPhaseProperties()); + + // Construct incompressible turbulence model + autoPtr turbulence ( - incompressible::LESModel::New(U, phi, twoPhaseProperties()) + incompressible::turbulenceModel::New(U, phi, twoPhaseProperties()) ); diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeFoam.C b/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeFoam.C index caae02a3f..2f44534ce 100644 --- a/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeFoam.C +++ b/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeFoam.C @@ -28,13 +28,17 @@ Application Description Solver for 2 incompressible, isothermal immiscible fluids with phase-change (e.g. cavitation). Uses a VOF (volume of fluid) phase-fraction based - interface capturing approach. The momentum and other fluid properties are - of the "mixture" and a single momentum equation is solved. + interface capturing approach. + + The momentum and other fluid properties are of the "mixture" and a + single momentum equation is solved. The set of phase-change models provided are designed to simulate cavitation but other mechanisms of phase-change are supported within this solver framework. + Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected. + \*---------------------------------------------------------------------------*/ #include "fvCFD.H" @@ -42,23 +46,23 @@ Description #include "subCycle.H" #include "interfaceProperties.H" #include "phaseChangeTwoPhaseMixture.H" -#include "incompressible/LESModel/LESModel.H" +#include "turbulenceModel.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // int main(int argc, char *argv[]) { - #include "setRootCase.H" - #include "createTime.H" - #include "createMesh.H" - #include "readEnvironmentalProperties.H" - #include "readPISOControls.H" - #include "initContinuityErrs.H" - #include "createFields.H" - #include "readTimeControls.H" - #include "correctPhi.H" - #include "CourantNo.H" - #include "setInitialDeltaT.H" +# include "setRootCase.H" +# include "createTime.H" +# include "createMesh.H" +# include "readGravitationalAcceleration.H" +# include "readPISOControls.H" +# include "initContinuityErrs.H" +# include "createFields.H" +# include "readTimeControls.H" +# include "correctPhi.H" +# include "CourantNo.H" +# include "setInitialDeltaT.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -66,18 +70,16 @@ int main(int argc, char *argv[]) while (runTime.run()) { - #include "readPISOControls.H" - #include "readTimeControls.H" - #include "CourantNo.H" - #include "setDeltaT.H" +# include "readPISOControls.H" +# include "readTimeControls.H" +# include "CourantNo.H" +# include "setDeltaT.H" runTime++; Info<< "Time = " << runTime.timeName() << nl << endl; - twoPhaseProperties->correct(); - - #include "gammaEqnSubCycle.H" +# include "alphaEqnSubCycle.H" turbulence->correct(); @@ -89,12 +91,14 @@ int main(int argc, char *argv[]) // --- PISO loop for (int corr=0; corrcorrect(); + runTime.write(); Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s" @@ -104,7 +108,7 @@ int main(int argc, char *argv[]) Info<< "End\n" << endl; - return(0); + return 0; } diff --git a/applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H b/applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H index 3a7f7e4eb..a89d2d08a 100644 --- a/applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H +++ b/applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H @@ -11,13 +11,14 @@ + fvc::ddtPhiCorr(rUA, rho, U, phi) ); + adjustPhi(phiU, U, p); + phi = phiU + ( - fvc::interpolate(interface.sigmaK())*fvc::snGrad(gamma) + fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1) - ghf*fvc::snGrad(rho) )*rUAf*mesh.magSf(); - adjustPhi(phi, U, pd); Pair > vDotP = twoPhaseProperties->vDotP(); const volScalarField& vDotcP = vDotP[0](); diff --git a/applications/solvers/multiphase/lesCavitatingFoam/Make/files b/applications/solvers/multiphase/lesCavitatingFoam/Make/files deleted file mode 100644 index b577562e2..000000000 --- a/applications/solvers/multiphase/lesCavitatingFoam/Make/files +++ /dev/null @@ -1,5 +0,0 @@ -lesCavitatingFoam.C - -devOneEqEddy/devOneEqEddy.C - -EXE = $(FOAM_APPBIN)/lesCavitatingFoam diff --git a/applications/solvers/multiphase/lesCavitatingFoam/UEqn.H b/applications/solvers/multiphase/lesCavitatingFoam/UEqn.H deleted file mode 100644 index 73033ddfa..000000000 --- a/applications/solvers/multiphase/lesCavitatingFoam/UEqn.H +++ /dev/null @@ -1,20 +0,0 @@ - surfaceScalarField muEff - ( - "muEff", - twoPhaseProperties.muf() - + fvc::interpolate(rho*turbulence->nuSgs()) - ); - - fvVectorMatrix UEqn - ( - fvm::ddt(rho, U) - + fvm::div(phi, U) - - fvm::laplacian(muEff, U) - //- (fvc::grad(U) & fvc::grad(muf)) - - fvc::div(muEff*(fvc::interpolate(dev(fvc::grad(U))) & mesh.Sf())) - ); - - if (momentumPredictor) - { - solve(UEqn == -fvc::grad(p)); - } diff --git a/applications/solvers/multiphase/lesCavitatingFoam/createFields.H b/applications/solvers/multiphase/lesCavitatingFoam/createFields.H deleted file mode 100644 index d2e19bcae..000000000 --- a/applications/solvers/multiphase/lesCavitatingFoam/createFields.H +++ /dev/null @@ -1,85 +0,0 @@ - Info<< "Reading field p\n" << endl; - volScalarField p - ( - IOobject - ( - "p", - runTime.timeName(), - mesh, - IOobject::MUST_READ, - IOobject::AUTO_WRITE - ), - mesh - ); - - volScalarField rho - ( - IOobject - ( - "rho", - runTime.timeName(), - mesh, - IOobject::MUST_READ, - IOobject::AUTO_WRITE - ), - mesh - ); - - volScalarField gamma - ( - IOobject - ( - "gamma", - runTime.timeName(), - mesh, - IOobject::NO_READ, - IOobject::AUTO_WRITE - ), - max(min((rho - rholSat)/(rhovSat - rholSat), scalar(1)), scalar(0)) - ); - gamma.oldTime(); - - Info<< "Creating compressibilityModel\n" << endl; - autoPtr psiModel = - barotropicCompressibilityModel::New - ( - thermodynamicProperties, - gamma - ); - - const volScalarField& psi = psiModel->psi(); - - rho == max - ( - psi*p - + (1.0 - gamma)*rhol0 - + ((gamma*psiv + (1.0 - gamma)*psil) - psi)*pSat, - rhoMin - ); - - Info<< "Reading field U\n" << endl; - volVectorField U - ( - IOobject - ( - "U", - runTime.timeName(), - mesh, - IOobject::MUST_READ, - IOobject::AUTO_WRITE - ), - mesh - ); - -# include "createPhiv.H" -# include "compressibleCreatePhi.H" - - Info<< "Reading transportProperties\n" << endl; - - twoPhaseMixture twoPhaseProperties(U, phiv, "gamma"); - - // Create LES model - autoPtr turbulence - ( - incompressible::LESModel::New(U, phiv, twoPhaseProperties) - ); diff --git a/applications/solvers/multiphase/lesCavitatingFoam/devOneEqEddy/devOneEqEddy.C b/applications/solvers/multiphase/lesCavitatingFoam/devOneEqEddy/devOneEqEddy.C deleted file mode 100644 index 682a49a34..000000000 --- a/applications/solvers/multiphase/lesCavitatingFoam/devOneEqEddy/devOneEqEddy.C +++ /dev/null @@ -1,130 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | OpenFOAM: The Open Source CFD Toolbox - \\ / O peration | - \\ / A nd | Copyright held by original author - \\/ M anipulation | -------------------------------------------------------------------------------- -License - This file is part of OpenFOAM. - - OpenFOAM 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 2 of the License, or (at your - option) any later version. - - OpenFOAM 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 OpenFOAM; if not, write to the Free Software Foundation, - Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA - -\*---------------------------------------------------------------------------*/ - -#include "devOneEqEddy.H" -#include "addToRunTimeSelectionTable.H" - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -namespace Foam -{ -namespace incompressible -{ -namespace LESModels -{ - -// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // - -defineTypeNameAndDebug(devOneEqEddy, 0); -addToRunTimeSelectionTable(LESModel, devOneEqEddy, dictionary); - -// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // - -devOneEqEddy::devOneEqEddy -( - const volVectorField& U, - const surfaceScalarField& phi, - transportModel& transport -) -: - LESModel(typeName, U, phi, transport), - GenEddyVisc(U, phi, transport), - - k_ - ( - IOobject - ( - "k", - runTime_.timeName(), - mesh_, - IOobject::MUST_READ, - IOobject::AUTO_WRITE - ), - mesh_ - ), - ck_ - ( - dimensioned::lookupOrAddToDict - ( - "ck", - coeffDict(), - 0.07 - ) - ) -{ - printCoeffs(); -} - - -// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // - -void devOneEqEddy::correct(const tmp& gradU) -{ - GenEddyVisc::correct(gradU); - - //volScalarField G = 2*nuSgs_*magSqr(symm(gradU)); - volScalarField G = 2*nuSgs_*(gradU() && dev(symm(gradU()))); - - solve - ( - fvm::ddt(k_) - + fvm::div(phi(), k_) - - fvm::Sp(fvc::div(phi()), k_) - - fvm::laplacian(DkEff(), k_) - == - G - - fvm::Sp(ce_*sqrt(k_)/delta(), k_) - ); - - bound(k_, k0()); - - nuSgs_ = ck_*sqrt(k_)*delta(); - nuSgs_.correctBoundaryConditions(); -} - - -bool devOneEqEddy::read() -{ - if (GenEddyVisc::read()) - { - ck_.readIfPresent(coeffDict()); - - return true; - } - else - { - return false; - } -} - - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -} // End namespace LESModels -} // End namespace incompressible -} // End namespace Foam - -// ************************************************************************* // diff --git a/applications/solvers/multiphase/lesCavitatingFoam/lesCavitatingFoam.C b/applications/solvers/multiphase/lesCavitatingFoam/lesCavitatingFoam.C deleted file mode 100644 index f41a99923..000000000 --- a/applications/solvers/multiphase/lesCavitatingFoam/lesCavitatingFoam.C +++ /dev/null @@ -1,94 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | OpenFOAM: The Open Source CFD Toolbox - \\ / O peration | - \\ / A nd | Copyright held by original author - \\/ M anipulation | -------------------------------------------------------------------------------- -License - This file is part of OpenFOAM. - - OpenFOAM 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 2 of the License, or (at your - option) any later version. - - OpenFOAM 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 OpenFOAM; if not, write to the Free Software Foundation, - Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA - -Application - lesCavitatingFoam - -Description - Transient cavitation code with LES turbulence. - -\*---------------------------------------------------------------------------*/ - -#include "fvCFD.H" -#include "barotropicCompressibilityModel.H" -#include "twoPhaseMixture.H" -#include "incompressible/LESModel/LESModel.H" - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -int main(int argc, char *argv[]) -{ - -# include "setRootCase.H" - -# include "createTime.H" -# include "createMesh.H" -# include "readThermodynamicProperties.H" -# include "readControls.H" -# include "createFields.H" -# include "initContinuityErrs.H" -# include "compressibleCourantNo.H" -# include "setInitialDeltaT.H" - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - - Info<< "\nStarting time loop\n" << endl; - - while (runTime.run()) - { -# include "readControls.H" -# include "CourantNo.H" -# include "setDeltaT.H" - - runTime++; - Info<< "Time = " << runTime.timeName() << nl << endl; - - turbulence->correct(); - - for (int outerCorr=0; outerCorr turbulence + ( + incompressible::turbulenceModel::New(U, phi, mixture) + ); diff --git a/applications/solvers/multiphase/multiphaseInterFoam/multiphaseInterFoam.C b/applications/solvers/multiphase/multiphaseInterFoam/multiphaseInterFoam.C index 4d7f78da9..54cab3275 100644 --- a/applications/solvers/multiphase/multiphaseInterFoam/multiphaseInterFoam.C +++ b/applications/solvers/multiphase/multiphaseInterFoam/multiphaseInterFoam.C @@ -27,22 +27,24 @@ Application Description Solver for n incompressible fluids which captures the interfaces and - includes surface-tension and contact-angle effects for each. + includes surface-tension and contact-angle effects for each phase. + + Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected. \*---------------------------------------------------------------------------*/ #include "fvCFD.H" #include "multiphaseMixture.H" +#include "turbulenceModel.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // int main(int argc, char *argv[]) { - # include "setRootCase.H" # include "createTime.H" # include "createMesh.H" -# include "readEnvironmentalProperties.H" + #include "readGravitationalAcceleration.H" # include "readPISOControls.H" # include "initContinuityErrs.H" # include "createFields.H" @@ -79,6 +81,8 @@ int main(int argc, char *argv[]) # include "continuityErrs.H" + turbulence->correct(); + runTime.write(); Info<< "ExecutionTime = " @@ -88,7 +92,7 @@ int main(int argc, char *argv[]) Info<< "\n end \n"; - return(0); + return 0; } diff --git a/applications/solvers/multiphase/porousInterFoam/Make/files b/applications/solvers/multiphase/porousInterFoam/Make/files new file mode 100644 index 000000000..dda65f11f --- /dev/null +++ b/applications/solvers/multiphase/porousInterFoam/Make/files @@ -0,0 +1,3 @@ +porousInterFoam.C + +EXE = $(FOAM_APPBIN)/porousInterFoam diff --git a/applications/solvers/multiphase/compressibleLesInterFoam/Make/options b/applications/solvers/multiphase/porousInterFoam/Make/options similarity index 58% rename from applications/solvers/multiphase/compressibleLesInterFoam/Make/options rename to applications/solvers/multiphase/porousInterFoam/Make/options index 45e4e10f6..82185b813 100644 --- a/applications/solvers/multiphase/compressibleLesInterFoam/Make/options +++ b/applications/solvers/multiphase/porousInterFoam/Make/options @@ -1,15 +1,19 @@ INTERFOAM = $(FOAM_SOLVERS)/multiphase/interFoam EXE_INC = \ + -I$(INTERFOAM) \ -I$(LIB_SRC)/transportModels \ -I$(LIB_SRC)/transportModels/incompressible/lnInclude \ -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \ - -I$(LIB_SRC)/turbulenceModels/LES \ - -I$(LIB_SRC)/turbulenceModels/LES/LESdeltas/lnInclude \ - -I$(LIB_SRC)/finiteVolume/lnInclude + -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \ + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/meshTools/lnInclude EXE_LIBS = \ -linterfaceProperties \ -lincompressibleTransportModels \ + -lincompressibleRASModels \ -lincompressibleLESModels \ - -lfiniteVolume + -lfiniteVolume \ + -lmeshTools + diff --git a/applications/solvers/multiphase/porousInterFoam/UEqn.H b/applications/solvers/multiphase/porousInterFoam/UEqn.H new file mode 100644 index 000000000..efb1df5ea --- /dev/null +++ b/applications/solvers/multiphase/porousInterFoam/UEqn.H @@ -0,0 +1,39 @@ + surfaceScalarField muEff + ( + "muEff", + twoPhaseProperties.muf() + + fvc::interpolate(rho*turbulence->nut()) + ); + + // Calculate and cache mu for the porous media + volScalarField mu(twoPhaseProperties.mu()); + + fvVectorMatrix UEqn + ( + pZones.ddt(rho, U) + + fvm::div(rhoPhi, U) + - fvm::laplacian(muEff, U) + - (fvc::grad(U) & fvc::grad(muEff)) + //- fvc::div(muEff*(fvc::interpolate(dev(fvc::grad(U))) & mesh.Sf())) + ); + + UEqn.relax(); + + pZones.addResistance(UEqn); + + if (momentumPredictor) + { + solve + ( + UEqn + == + fvc::reconstruct + ( + fvc::interpolate(rho)*(g & mesh.Sf()) + + ( + fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1) + - fvc::snGrad(p) + ) * mesh.magSf() + ) + ); + } diff --git a/applications/solvers/multiphase/rasInterFoam/createFields.H b/applications/solvers/multiphase/porousInterFoam/createFields.H similarity index 78% rename from applications/solvers/multiphase/rasInterFoam/createFields.H rename to applications/solvers/multiphase/porousInterFoam/createFields.H index 2d932dadc..f441e8b33 100644 --- a/applications/solvers/multiphase/rasInterFoam/createFields.H +++ b/applications/solvers/multiphase/porousInterFoam/createFields.H @@ -12,12 +12,12 @@ mesh ); - Info<< "Reading field gamma\n" << endl; - volScalarField gamma + Info<< "Reading field alpha1\n" << endl; + volScalarField alpha1 ( IOobject ( - "gamma", + "alpha1", runTime.timeName(), mesh, IOobject::MUST_READ, @@ -42,8 +42,9 @@ # include "createPhi.H" + Info<< "Reading transportProperties\n" << endl; - twoPhaseMixture twoPhaseProperties(U, phi, "gamma"); + twoPhaseMixture twoPhaseProperties(U, phi, "alpha1"); const dimensionedScalar& rho1 = twoPhaseProperties.rho1(); const dimensionedScalar& rho2 = twoPhaseProperties.rho2(); @@ -59,15 +60,15 @@ mesh, IOobject::READ_IF_PRESENT ), - gamma*rho1 + (scalar(1) - gamma)*rho2, - gamma.boundaryField().types() + alpha1*rho1 + (scalar(1) - alpha1)*rho2, + alpha1.boundaryField().types() ); rho.oldTime(); // Mass flux // Initialisation does not matter because rhoPhi is reset after the - // gamma solution before it is used in the U equation. + // alpha1 solution before it is used in the U equation. surfaceScalarField rhoPhi ( IOobject @@ -86,7 +87,6 @@ volScalarField gh("gh", g & mesh.C()); surfaceScalarField ghf("gh", g & mesh.Cf()); - volScalarField p ( IOobject @@ -122,11 +122,13 @@ ); } - // Construct interface from gamma distribution - interfaceProperties interface(gamma, U, twoPhaseProperties); + // Construct interface from alpha1 distribution + interfaceProperties interface(alpha1, U, twoPhaseProperties); - // Construct RAS model - autoPtr turbulence + // Construct incompressible turbulence model + autoPtr turbulence ( - incompressible::RASModel::New(U, phi, twoPhaseProperties) + incompressible::turbulenceModel::New(U, phi, twoPhaseProperties) ); + + porousZones pZones(mesh); diff --git a/applications/solvers/multiphase/rasInterFoam/rasInterFoam.C b/applications/solvers/multiphase/porousInterFoam/porousInterFoam.C similarity index 82% rename from applications/solvers/multiphase/rasInterFoam/rasInterFoam.C rename to applications/solvers/multiphase/porousInterFoam/porousInterFoam.C index 4fee888ef..3c7643ecd 100644 --- a/applications/solvers/multiphase/rasInterFoam/rasInterFoam.C +++ b/applications/solvers/multiphase/porousInterFoam/porousInterFoam.C @@ -23,14 +23,19 @@ License Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA Application - rasInterFoam + porousInterFoam Description Solver for 2 incompressible, isothermal immiscible fluids using a VOF (volume of fluid) phase-fraction based interface capturing approach. The momentum and other fluid properties are of the "mixture" and a single - momentum equation is solved. Turbulence is modelled using a run-time - selectable incompressible RAS model. + momentum equation is solved. + + Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected. + + For a two-fluid approach see twoPhaseEulerFoam. + + Explicit handling of porous zones is included. \*---------------------------------------------------------------------------*/ @@ -39,7 +44,8 @@ Description #include "subCycle.H" #include "interfaceProperties.H" #include "twoPhaseMixture.H" -#include "incompressible/RASModel/RASModel.H" +#include "turbulenceModel.H" +#include "porousZones.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -48,7 +54,7 @@ int main(int argc, char *argv[]) # include "setRootCase.H" # include "createTime.H" # include "createMesh.H" -# include "readEnvironmentalProperties.H" +# include "readGravitationalAcceleration.H" # include "readPISOControls.H" # include "initContinuityErrs.H" # include "createFields.H" @@ -57,7 +63,7 @@ int main(int argc, char *argv[]) # include "CourantNo.H" # include "setInitialDeltaT.H" -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // Info<< "\nStarting time loop\n" << endl; @@ -72,30 +78,20 @@ int main(int argc, char *argv[]) Info<< "Time = " << runTime.timeName() << nl << endl; -# include "gammaEqnSubCycle.H" + twoPhaseProperties.correct(); + +# include "alphaEqnSubCycle.H" # include "UEqn.H" // --- PISO loop - for (int corr=0; corr < nCorr; corr++) + for (int corr=0; corrcorrect(); runTime.write(); @@ -107,7 +103,7 @@ int main(int argc, char *argv[]) Info<< "End\n" << endl; - return(0); + return 0; } diff --git a/applications/solvers/multiphase/rasCavitatingFoam/CourantNo.H b/applications/solvers/multiphase/rasCavitatingFoam/CourantNo.H deleted file mode 100644 index b7450c92d..000000000 --- a/applications/solvers/multiphase/rasCavitatingFoam/CourantNo.H +++ /dev/null @@ -1,59 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | OpenFOAM: The Open Source CFD Toolbox - \\ / O peration | - \\ / A nd | Copyright held by original author - \\/ M anipulation | -------------------------------------------------------------------------------- -License - This file is part of OpenFOAM. - - OpenFOAM 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 2 of the License, or (at your - option) any later version. - - OpenFOAM 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 OpenFOAM; if not, write to the Free Software Foundation, - Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA - -Global - CourantNo - -Description - Calculates and outputs the mean and maximum Courant Numbers. - -\*---------------------------------------------------------------------------*/ - -scalar CoNum = 0.0; -scalar meanCoNum = 0.0; -scalar acousticCoNum = 0.0; - -if (mesh.nInternalFaces()) -{ - surfaceScalarField SfUfbyDelta = - mesh.surfaceInterpolation::deltaCoeffs()*mag(phiv); - - CoNum = max(SfUfbyDelta/mesh.magSf()) - .value()*runTime.deltaT().value(); - - meanCoNum = (sum(SfUfbyDelta)/sum(mesh.magSf())) - .value()*runTime.deltaT().value(); - - acousticCoNum = max - ( - mesh.surfaceInterpolation::deltaCoeffs()/sqrt(fvc::interpolate(psi)) - ).value()*runTime.deltaT().value(); -} - -Info<< "phiv Courant Number mean: " << meanCoNum - << " max: " << CoNum - << " acoustic max: " << acousticCoNum - << endl; - -// ************************************************************************* // diff --git a/applications/solvers/multiphase/rasCavitatingFoam/Make/files b/applications/solvers/multiphase/rasCavitatingFoam/Make/files deleted file mode 100644 index 450946379..000000000 --- a/applications/solvers/multiphase/rasCavitatingFoam/Make/files +++ /dev/null @@ -1,3 +0,0 @@ -rasCavitatingFoam.C - -EXE = $(FOAM_APPBIN)/rasCavitatingFoam diff --git a/applications/solvers/multiphase/rasCavitatingFoam/continuityErrs.H b/applications/solvers/multiphase/rasCavitatingFoam/continuityErrs.H deleted file mode 100644 index 6f1622510..000000000 --- a/applications/solvers/multiphase/rasCavitatingFoam/continuityErrs.H +++ /dev/null @@ -1,22 +0,0 @@ -{ - volScalarField thermoRho = psi*p + (1.0 - gamma)*rhol0; - - dimensionedScalar totalMass = fvc::domainIntegrate(rho); - - scalar sumLocalContErr = - ( - fvc::domainIntegrate(mag(rho - thermoRho))/totalMass - ).value(); - - scalar globalContErr = - ( - fvc::domainIntegrate(rho - thermoRho)/totalMass - ).value(); - - cumulativeContErr += globalContErr; - - Info<< "time step continuity errors : sum local = " << sumLocalContErr - << ", global = " << globalContErr - << ", cumulative = " << cumulativeContErr - << endl; -} diff --git a/applications/solvers/multiphase/rasCavitatingFoam/gammaPsi.H b/applications/solvers/multiphase/rasCavitatingFoam/gammaPsi.H deleted file mode 100644 index b259ddd32..000000000 --- a/applications/solvers/multiphase/rasCavitatingFoam/gammaPsi.H +++ /dev/null @@ -1,10 +0,0 @@ -{ - gamma = max(min((rho - rholSat)/(rhovSat - rholSat), scalar(1)), scalar(0)); - - Info<< "max-min gamma: " << max(gamma).value() - << " " << min(gamma).value() << endl; - - psiModel->correct(); - - //Info<< "min a: " << 1.0/sqrt(max(psi)).value() << endl; -} diff --git a/applications/solvers/multiphase/rasCavitatingFoam/readControls.H b/applications/solvers/multiphase/rasCavitatingFoam/readControls.H deleted file mode 100644 index f53e7b9eb..000000000 --- a/applications/solvers/multiphase/rasCavitatingFoam/readControls.H +++ /dev/null @@ -1,9 +0,0 @@ -#include "readTimeControls.H" - -scalar maxAcousticCo -( - readScalar(runTime.controlDict().lookup("maxAcousticCo")) -); - - -#include "readPISOControls.H" diff --git a/applications/solvers/multiphase/rasCavitatingFoam/readThermodynamicProperties.H b/applications/solvers/multiphase/rasCavitatingFoam/readThermodynamicProperties.H deleted file mode 100644 index d3fbb9307..000000000 --- a/applications/solvers/multiphase/rasCavitatingFoam/readThermodynamicProperties.H +++ /dev/null @@ -1,27 +0,0 @@ - Info<< "Reading thermodynamicProperties\n" << endl; - - IOdictionary thermodynamicProperties - ( - IOobject - ( - "thermodynamicProperties", - runTime.constant(), - mesh, - IOobject::MUST_READ, - IOobject::NO_WRITE - ) - ); - - dimensionedScalar psil(thermodynamicProperties.lookup("psil")); - - dimensionedScalar rholSat(thermodynamicProperties.lookup("rholSat")); - - dimensionedScalar psiv(thermodynamicProperties.lookup("psiv")); - - dimensionedScalar pSat(thermodynamicProperties.lookup("pSat")); - - dimensionedScalar rhovSat("rhovSat", psiv*pSat); - - dimensionedScalar rhol0("rhol0", rholSat - pSat*psil); - - dimensionedScalar rhoMin(thermodynamicProperties.lookup("rhoMin")); diff --git a/applications/solvers/multiphase/rasCavitatingFoam/resetPhiPatches.H b/applications/solvers/multiphase/rasCavitatingFoam/resetPhiPatches.H deleted file mode 100644 index e7d0c2f93..000000000 --- a/applications/solvers/multiphase/rasCavitatingFoam/resetPhiPatches.H +++ /dev/null @@ -1,15 +0,0 @@ -fvsPatchScalarFieldField& phiPatches = phi.boundaryField(); -const fvPatchScalarFieldField& rhoPatches = rho.boundaryField(); -const fvPatchVectorFieldField& Upatches = U.boundaryField(); -const fvsPatchVectorFieldField& SfPatches = mesh.Sf().boundaryField(); - -forAll(phiPatches, patchI) -{ - if (phi.boundaryField().types()[patchI] == "calculated") - { - calculatedFvsPatchScalarField& phiPatch = - refCast(phiPatches[patchI]); - - phiPatch == ((rhoPatches[patchI]*Upatches[patchI]) & SfPatches[patchI]); - } -} diff --git a/applications/solvers/multiphase/rasCavitatingFoam/resetPhivPatches.H b/applications/solvers/multiphase/rasCavitatingFoam/resetPhivPatches.H deleted file mode 100644 index 7e8b040bb..000000000 --- a/applications/solvers/multiphase/rasCavitatingFoam/resetPhivPatches.H +++ /dev/null @@ -1,14 +0,0 @@ -surfaceScalarField::GeometricBoundaryField& phivPatches = phiv.boundaryField(); -const volVectorField::GeometricBoundaryField& Upatches = U.boundaryField(); -const surfaceVectorField::GeometricBoundaryField& SfPatches = mesh.Sf().boundaryField(); - -forAll(phivPatches, patchI) -{ - if (phiv.boundaryField().types()[patchI] == "calculated") - { - calculatedFvsPatchScalarField& phivPatch = - refCast(phivPatches[patchI]); - - phivPatch == (Upatches[patchI] & SfPatches[patchI]); - } -} diff --git a/applications/solvers/multiphase/rasCavitatingFoam/rhoEqn.H b/applications/solvers/multiphase/rasCavitatingFoam/rhoEqn.H deleted file mode 100644 index d0bd6e1da..000000000 --- a/applications/solvers/multiphase/rasCavitatingFoam/rhoEqn.H +++ /dev/null @@ -1,16 +0,0 @@ -{ - fvScalarMatrix rhoEqn - ( - fvm::ddt(rho) - + fvm::div(phiv, rho) - ); - - rhoEqn.solve(); - - phi = rhoEqn.flux(); - - Info<< "max-min rho: " << max(rho).value() - << " " << min(rho).value() << endl; - - rho == max(rho, rhoMin); -} diff --git a/applications/solvers/multiphase/rasCavitatingFoam/setDeltaT.H b/applications/solvers/multiphase/rasCavitatingFoam/setDeltaT.H deleted file mode 100644 index 9233804c4..000000000 --- a/applications/solvers/multiphase/rasCavitatingFoam/setDeltaT.H +++ /dev/null @@ -1,54 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | OpenFOAM: The Open Source CFD Toolbox - \\ / O peration | - \\ / A nd | Copyright held by original author - \\/ M anipulation | -------------------------------------------------------------------------------- -License - This file is part of OpenFOAM. - - OpenFOAM 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 2 of the License, or (at your - option) any later version. - - OpenFOAM 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 OpenFOAM; if not, write to the Free Software Foundation, - Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA - -Global - setDeltaT - -Description - Reset the timestep to maintain a constant maximum courant Number. - Reduction of time-step is imediate but increase is damped to avoid - unstable oscillations. - -\*---------------------------------------------------------------------------*/ - -if (adjustTimeStep) -{ - scalar maxDeltaTFact = - min(maxCo/(CoNum + SMALL), maxAcousticCo/(acousticCoNum + SMALL)); - - scalar deltaTFact = min(min(maxDeltaTFact, 1.0 + 0.1*maxDeltaTFact), 1.2); - - runTime.setDeltaT - ( - min - ( - deltaTFact*runTime.deltaT().value(), - maxDeltaT - ) - ); - - Info<< "deltaT = " << runTime.deltaT().value() << endl; -} - -// ************************************************************************* // diff --git a/applications/solvers/multiphase/rasCavitatingFoam/setInitialDeltaT.H b/applications/solvers/multiphase/rasCavitatingFoam/setInitialDeltaT.H deleted file mode 100644 index 61823a001..000000000 --- a/applications/solvers/multiphase/rasCavitatingFoam/setInitialDeltaT.H +++ /dev/null @@ -1,54 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | OpenFOAM: The Open Source CFD Toolbox - \\ / O peration | - \\ / A nd | Copyright held by original author - \\/ M anipulation | -------------------------------------------------------------------------------- -License - This file is part of OpenFOAM. - - OpenFOAM 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 2 of the License, or (at your - option) any later version. - - OpenFOAM 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 OpenFOAM; if not, write to the Free Software Foundation, - Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA - -Global - setInitialDeltaT - -Description - Set the initial timestep corresponding to the timestep adjustment - algorithm in setDeltaT - -\*---------------------------------------------------------------------------*/ - -if (adjustTimeStep) -{ -# include "CourantNo.H" - - if (CoNum > SMALL) - { - scalar maxDeltaTFact = - min(maxCo/(CoNum + SMALL), maxAcousticCo/(acousticCoNum + SMALL)); - - runTime.setDeltaT - ( - min - ( - maxDeltaTFact*runTime.deltaT().value(), - maxDeltaT - ) - ); - } -} - -// ************************************************************************* // diff --git a/applications/solvers/multiphase/rasInterFoam/Make/files b/applications/solvers/multiphase/rasInterFoam/Make/files deleted file mode 100644 index 39b992ad4..000000000 --- a/applications/solvers/multiphase/rasInterFoam/Make/files +++ /dev/null @@ -1,3 +0,0 @@ -rasInterFoam.C - -EXE = $(FOAM_APPBIN)/rasInterFoam diff --git a/applications/solvers/multiphase/settlingFoam/settlingFoam.C b/applications/solvers/multiphase/settlingFoam/settlingFoam.C index b6a464c1d..072e92794 100644 --- a/applications/solvers/multiphase/settlingFoam/settlingFoam.C +++ b/applications/solvers/multiphase/settlingFoam/settlingFoam.C @@ -43,12 +43,11 @@ Description int main(int argc, char *argv[]) { - # include "setRootCase.H" # include "createTime.H" # include "createMesh.H" -# include "readEnvironmentalProperties.H" +# include "readGravitationalAcceleration.H" # include "createFields.H" # include "initContinuityErrs.H" @@ -57,7 +56,7 @@ int main(int argc, char *argv[]) Info<< "\nStarting time loop\n" << endl; - for (runTime++; !runTime.end(); runTime++) + while (runTime.loop()) { Info<< "Time = " << runTime.timeName() << nl << endl; @@ -66,22 +65,24 @@ int main(int argc, char *argv[]) # include "rhoEqn.H" -# include "calcVdj.H" - -# include "UEqn.H" - -# include "alphaEqn.H" - -# include "correctViscosity.H" - - - // --- PISO loop - for (int corr=0; corrnut()) + ); fvVectorMatrix UEqn ( fvm::ddt(rho, U) + fvm::div(rhoPhi, U) - - fvm::laplacian(muf, U) - //- (fvc::grad(U) & fvc::grad(muf)) - - fvc::div(muf*(mesh.Sf() & fvc::interpolate(fvc::grad(U)().T()))) + - fvm::laplacian(muEff, U) + - (fvc::grad(U) & fvc::grad(muEff)) + //- fvc::div(muEff*(fvc::interpolate(dev(fvc::grad(U))) & mesh.Sf())) ); + if (oCorr == nOuterCorr-1) + { + UEqn.relax(1); + } + else + { + UEqn.relax(); + } + if (momentumPredictor) { solve ( UEqn == - -fvc::reconstruct + fvc::reconstruct ( - mesh.magSf()*(fvc::snGrad(pd) + ghf*fvc::snGrad(rho)) - ) + (- ghf*fvc::snGrad(rho) - fvc::snGrad(p_rgh))*mesh.magSf() + ), + mesh.solver(oCorr == nOuterCorr-1 ? "UFinal" : "U") ); } diff --git a/applications/solvers/multiphase/twoLiquidMixingFoam/alphaEqn.H b/applications/solvers/multiphase/twoLiquidMixingFoam/alphaEqn.H new file mode 100644 index 000000000..8194753c8 --- /dev/null +++ b/applications/solvers/multiphase/twoLiquidMixingFoam/alphaEqn.H @@ -0,0 +1,24 @@ +{ + fvScalarMatrix alpha1Eqn + ( + fvm::ddt(alpha1) + + fvm::div(phi, alpha1) + //- fvm::Sp(fvc::div(phi), alpha1) + - fvm::laplacian + ( + Dab + alphatab*turbulence->nut(), alpha1, + "laplacian(Dab,alpha1)" + ) + ); + + alpha1Eqn.solve(); + + rhoPhi = alpha1Eqn.flux()*(rho1 - rho2) + phi*rho2; + rho = alpha1*rho1 + (scalar(1) - alpha1)*rho2; + + Info<< "Phase 1 volume fraction = " + << alpha1.weightedAverage(mesh.V()).value() + << " Min(alpha1) = " << min(alpha1).value() + << " Max(alpha1) = " << max(alpha1).value() + << endl; +} diff --git a/applications/solvers/multiphase/twoLiquidMixingFoam/createFields.H b/applications/solvers/multiphase/twoLiquidMixingFoam/createFields.H index 652d1ad64..7cd78a521 100644 --- a/applications/solvers/multiphase/twoLiquidMixingFoam/createFields.H +++ b/applications/solvers/multiphase/twoLiquidMixingFoam/createFields.H @@ -1,9 +1,9 @@ - Info<< "Reading field pd\n" << endl; - volScalarField pd + Info<< "Reading field p_rgh\n" << endl; + volScalarField p_rgh ( IOobject ( - "pd", + "p_rgh", runTime.timeName(), mesh, IOobject::MUST_READ, @@ -12,12 +12,12 @@ mesh ); - Info<< "Reading field gamma\n" << endl; - volScalarField gamma + Info<< "Reading field alpha1\n" << endl; + volScalarField alpha1 ( IOobject ( - "gamma", + "alpha1", runTime.timeName(), mesh, IOobject::MUST_READ, @@ -40,24 +40,27 @@ mesh ); -# include "createPhi.H" + #include "createPhi.H" Info<< "Reading transportProperties\n" << endl; - twoPhaseMixture twoPhaseProperties(U, phi, "gamma"); - + twoPhaseMixture twoPhaseProperties(U, phi); + const dimensionedScalar& rho1 = twoPhaseProperties.rho1(); const dimensionedScalar& rho2 = twoPhaseProperties.rho2(); dimensionedScalar Dab(twoPhaseProperties.lookup("Dab")); + // Read the reciprocal of the turbulent Schmidt number + dimensionedScalar alphatab(twoPhaseProperties.lookup("alphatab")); + // Need to store rho for ddt(rho, U) - volScalarField rho("rho", gamma*rho1 + (scalar(1) - gamma)*rho2); + volScalarField rho("rho", alpha1*rho1 + (scalar(1) - alpha1)*rho2); rho.oldTime(); // Mass flux // Initialisation does not matter because rhoPhi is reset after the - // gamma solution before it is used in the U equation. + // alpha1 solution before it is used in the U equation. surfaceScalarField rhoPhi ( IOobject @@ -71,11 +74,53 @@ rho1*phi ); + Info<< "Calculating field (g.h)f\n" << endl; + volScalarField gh("gh", g & mesh.C()); + surfaceScalarField ghf = surfaceScalarField("gh", g & mesh.Cf()); - Info<< "Calculating field g.h\n" << endl; - surfaceScalarField ghf("gh", g & mesh.Cf()); + volScalarField p + ( + IOobject + ( + "p", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + p_rgh + rho*gh + ); + + label p_rghRefCell = 0; + scalar p_rghRefValue = 0.0; + setRefCell + ( + p_rgh, + mesh.solutionDict().subDict("PIMPLE"), + p_rghRefCell, + p_rghRefValue + ); + + scalar pRefValue = 0.0; + + if (p_rgh.needReference()) + { + pRefValue = readScalar + ( + mesh.solutionDict().subDict("PIMPLE").lookup("pRefValue") + ); + + p += dimensionedScalar + ( + "p", + p.dimensions(), + pRefValue - getRefCellValue(p, p_rghRefCell) + ); + } - label pdRefCell = 0; - scalar pdRefValue = 0.0; - setRefCell(pd, mesh.solutionDict().subDict("PISO"), pdRefCell, pdRefValue); + // Construct incompressible turbulence model + autoPtr turbulence + ( + incompressible::turbulenceModel::New(U, phi, twoPhaseProperties) + ); diff --git a/applications/solvers/multiphase/twoLiquidMixingFoam/gammaEqn.H b/applications/solvers/multiphase/twoLiquidMixingFoam/gammaEqn.H deleted file mode 100644 index 144339f39..000000000 --- a/applications/solvers/multiphase/twoLiquidMixingFoam/gammaEqn.H +++ /dev/null @@ -1,19 +0,0 @@ -{ - fvScalarMatrix gammaEqn - ( - fvm::ddt(gamma) - + fvm::div(phi, gamma) - - fvm::laplacian(Dab, gamma) - ); - - gammaEqn.solve(); - - rhoPhi = gammaEqn.flux()*(rho1 - rho2) + phi*rho2; - rho = gamma*rho1 + (scalar(1) - gamma)*rho2; - - Info<< "Phase 1 volume fraction = " - << gamma.weightedAverage(mesh.V()).value() - << " Min(gamma) = " << min(gamma).value() - << " Max(gamma) = " << max(gamma).value() - << endl; -} diff --git a/applications/solvers/multiphase/twoLiquidMixingFoam/pEqn.H b/applications/solvers/multiphase/twoLiquidMixingFoam/pEqn.H deleted file mode 100644 index eaa0cef62..000000000 --- a/applications/solvers/multiphase/twoLiquidMixingFoam/pEqn.H +++ /dev/null @@ -1,35 +0,0 @@ -{ - volScalarField rUA = 1.0/UEqn.A(); - surfaceScalarField rUAf = fvc::interpolate(rUA); - - U = rUA*UEqn.H(); - - surfaceScalarField phiU - ( - "phiU", - (fvc::interpolate(U) & mesh.Sf()) + fvc::ddtPhiCorr(rUA, rho, U, phi) - ); - - phi = phiU - ghf*fvc::snGrad(rho)*rUAf*mesh.magSf(); - - adjustPhi(phi, U, pd); - - for(int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) - { - fvScalarMatrix pdEqn - ( - fvm::laplacian(rUAf, pd) == fvc::div(phi) - ); - - pdEqn.setReference(pdRefCell, pdRefValue); - pdEqn.solve(); - - if (nonOrth == nNonOrthCorr) - { - phi -= pdEqn.flux(); - } - } - - U += rUA*fvc::reconstruct((phi - phiU)/rUAf); - U.correctBoundaryConditions(); -} diff --git a/applications/solvers/multiphase/twoLiquidMixingFoam/p_rghEqn.H b/applications/solvers/multiphase/twoLiquidMixingFoam/p_rghEqn.H new file mode 100644 index 000000000..6e12b0478 --- /dev/null +++ b/applications/solvers/multiphase/twoLiquidMixingFoam/p_rghEqn.H @@ -0,0 +1,56 @@ +{ + volScalarField rUA = 1.0/UEqn.A(); + surfaceScalarField rUAf = fvc::interpolate(rUA); + + U = rUA*UEqn.H(); + + surfaceScalarField phiU + ( + "phiU", + (fvc::interpolate(U) & mesh.Sf()) + + fvc::ddtPhiCorr(rUA, rho, U, phi) + ); + + phi = phiU - ghf*fvc::snGrad(rho)*rUAf*mesh.magSf(); + + for(int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) + { + fvScalarMatrix p_rghEqn + ( + fvm::laplacian(rUAf, p_rgh) == fvc::div(phi) + ); + + p_rghEqn.setReference(p_rghRefCell, p_rghRefValue); + + p_rghEqn.solve + ( + mesh.solver + ( + p_rgh.name() + + ((corr == nCorr-1 && nonOrth == nNonOrthCorr) ? "Final" : "") + ) + ); + + if (nonOrth == nNonOrthCorr) + { + phi -= p_rghEqn.flux(); + } + } + + p == p_rgh + rho*gh; + + if (p_rgh.needReference()) + { + p += dimensionedScalar + ( + "p", + p.dimensions(), + pRefValue - getRefCellValue(p, p_rghRefCell) + ); + } + + #include "continuityErrs.H" + + U += rUA*fvc::reconstruct((phi - phiU)/rUAf); + U.correctBoundaryConditions(); +} diff --git a/applications/solvers/multiphase/twoLiquidMixingFoam/twoLiquidMixingFoam.C b/applications/solvers/multiphase/twoLiquidMixingFoam/twoLiquidMixingFoam.C index 6f3174dc3..69caf3deb 100644 --- a/applications/solvers/multiphase/twoLiquidMixingFoam/twoLiquidMixingFoam.C +++ b/applications/solvers/multiphase/twoLiquidMixingFoam/twoLiquidMixingFoam.C @@ -28,47 +28,61 @@ Application Description Solver for mixing 2 incompressible fluids. + Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected. + \*---------------------------------------------------------------------------*/ #include "fvCFD.H" #include "twoPhaseMixture.H" +#include "turbulenceModel.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // int main(int argc, char *argv[]) { - # include "setRootCase.H" # include "createTime.H" # include "createMesh.H" -# include "readEnvironmentalProperties.H" +# include "readGravitationalAcceleration.H" +# include "readPIMPLEControls.H" # include "initContinuityErrs.H" # include "createFields.H" +# include "readTimeControls.H" +# include "CourantNo.H" +# include "setInitialDeltaT.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // Info<< "\nStarting time loop\n" << endl; - for (runTime++; !runTime.end(); runTime++) + while (runTime.run()) { +# include "readPIMPLEControls.H" +# include "readTimeControls.H" +# include "CourantNo.H" +# include "setDeltaT.H" + + runTime++; + Info<< "Time = " << runTime.timeName() << nl << endl; -# include "readPISOControls.H" -# include "CourantNo.H" - - twoPhaseProperties.correct(); - -# include "gammaEqn.H" - -# include "UEqn.H" - - // --- PISO loop - for (int corr=0; corrcorrect(); + } runTime.write(); @@ -79,7 +93,7 @@ int main(int argc, char *argv[]) Info<< "End\n" << endl; - return(0); + return 0; } diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C index b809e491a..292ff8046 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C @@ -47,12 +47,11 @@ Description int main(int argc, char *argv[]) { - # include "setRootCase.H" # include "createTime.H" # include "createMesh.H" -# include "readEnvironmentalProperties.H" +# include "readGravitationalAcceleration.H" # include "createFields.H" # include "readPPProperties.H" # include "initContinuityErrs.H" @@ -108,7 +107,7 @@ int main(int argc, char *argv[]) Info<< "End\n" << endl; - return(0); + return 0; } diff --git a/applications/solvers/newStressAnalysis/Allwmake b/applications/solvers/newStressAnalysis/Allwmake index 61fce76ac..965c9fc2d 100755 --- a/applications/solvers/newStressAnalysis/Allwmake +++ b/applications/solvers/newStressAnalysis/Allwmake @@ -1,4 +1,6 @@ #!/bin/sh +cd ${0%/*} || exit 1 # run from this directory + set -x wmake libso materialModels diff --git a/applications/solvers/newStressAnalysis/materialModels/cohesiveLaws/cohesiveLaw/cohesiveLaw.C b/applications/solvers/newStressAnalysis/materialModels/cohesiveLaws/cohesiveLaw/cohesiveLaw.C index 5a4b452c0..2ea3e3a7d 100644 --- a/applications/solvers/newStressAnalysis/materialModels/cohesiveLaws/cohesiveLaw/cohesiveLaw.C +++ b/applications/solvers/newStressAnalysis/materialModels/cohesiveLaws/cohesiveLaw/cohesiveLaw.C @@ -104,7 +104,7 @@ Foam::cohesiveLaw::~cohesiveLaw() void Foam::cohesiveLaw::writeDict(Ostream& os) const { - os.writeKeyword(type() + "Coeffs") + os.writeKeyword(word(type() + "Coeffs")) << cohesiveLawCoeffs(); } diff --git a/applications/solvers/newStressAnalysis/newContactStressFoam/newContactStressFoam.C b/applications/solvers/newStressAnalysis/newContactStressFoam/newContactStressFoam.C index f51d9ba12..5086cf4b7 100644 --- a/applications/solvers/newStressAnalysis/newContactStressFoam/newContactStressFoam.C +++ b/applications/solvers/newStressAnalysis/newContactStressFoam/newContactStressFoam.C @@ -58,7 +58,7 @@ int main(int argc, char *argv[]) // Force n-sqaured projection // polyPatch::setNSquaredProjection(true); - for (runTime++; !runTime.end(); runTime++) + while (runTime.loop()) { Info<< "Iteration: " << runTime.timeName() << nl << endl; diff --git a/applications/solvers/newStressAnalysis/newStressedFoam/newStressedFoam.C b/applications/solvers/newStressAnalysis/newStressedFoam/newStressedFoam.C index 7c440eeba..68420167a 100644 --- a/applications/solvers/newStressAnalysis/newStressedFoam/newStressedFoam.C +++ b/applications/solvers/newStressAnalysis/newStressedFoam/newStressedFoam.C @@ -52,11 +52,11 @@ int main(int argc, char *argv[]) // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // Info<< "\nCalculating displacement field\n" << endl; - + volTensorField gradU = fvc::grad(U); volScalarField rho = rheology.rho(); - for (runTime++; !runTime.end(); runTime++) + while (runTime.loop()) { Info<< "Iteration: " << runTime.timeName() << nl << endl; diff --git a/applications/solvers/stressAnalysis/icoFsiFoam/createFields.H b/applications/solvers/stressAnalysis/icoFsiFoam/createFields.H index b772ed7d5..0e1e4a701 100644 --- a/applications/solvers/stressAnalysis/icoFsiFoam/createFields.H +++ b/applications/solvers/stressAnalysis/icoFsiFoam/createFields.H @@ -12,13 +12,11 @@ ) ); - dimensionedScalar nu ( transportProperties.lookup("nu") ); - dimensionedScalar rhoFluid ( transportProperties.lookup("rho") @@ -56,11 +54,6 @@ # include "createPhi.H" - -// surfaceScalarField phiNet = phi; -// phiNet.rename("phiNet"); - - label pRefCell = 0; scalar pRefValue = 0.0; setRefCell(p, mesh.solutionDict().subDict("PISO"), pRefCell, pRefValue); diff --git a/applications/solvers/stressAnalysis/icoFsiFoam/createStressMesh.H b/applications/solvers/stressAnalysis/icoFsiFoam/createStressMesh.H index bfcb4eb27..e9a081363 100644 --- a/applications/solvers/stressAnalysis/icoFsiFoam/createStressMesh.H +++ b/applications/solvers/stressAnalysis/icoFsiFoam/createStressMesh.H @@ -9,11 +9,5 @@ ) ); - pointMesh pStressMesh(stressMesh); - - volPointInterpolation cpi - ( - stressMesh, - pStressMesh - ); + const volPointInterpolation& cpi = volPointInterpolation::New(stressMesh); diff --git a/applications/solvers/stressAnalysis/icoFsiFoam/icoFsiFoam.C b/applications/solvers/stressAnalysis/icoFsiFoam/icoFsiFoam.C index d0285f34c..a1a3ee468 100644 --- a/applications/solvers/stressAnalysis/icoFsiFoam/icoFsiFoam.C +++ b/applications/solvers/stressAnalysis/icoFsiFoam/icoFsiFoam.C @@ -26,9 +26,9 @@ Application icoFoam Description - Transient solver for incompressible, laminar flow of Newtonian fluids with + Transient solver for incompressible, laminar flow of Newtonian fluids with mesh motion. Set up as a fake fluid structure interaction solver - + \*---------------------------------------------------------------------------*/ #include "fvCFD.H" diff --git a/applications/utilities/mesh/advanced/combinePatchFaces/combinePatchFaces.C b/applications/utilities/mesh/advanced/combinePatchFaces/combinePatchFaces.C index 99af2c373..3ca753d9a 100644 --- a/applications/utilities/mesh/advanced/combinePatchFaces/combinePatchFaces.C +++ b/applications/utilities/mesh/advanced/combinePatchFaces/combinePatchFaces.C @@ -86,21 +86,25 @@ void checkSnapMesh // Max nonorthogonality allowed scalar maxNonOrtho(readScalar(snapDict.lookup("maxNonOrtho"))); + primitiveMesh::nonOrthThreshold_ = maxNonOrtho; + // Max concaveness allowed. scalar maxConcave(readScalar(snapDict.lookup("maxConcave"))); + primitiveMesh::faceAngleThreshold_ = maxConcave; + // Min volume allowed (factor of minimum cellVolume) scalar relMinVol(readScalar(snapDict.lookup("minVol"))); const scalar minCellVol = min(mesh.cellVolumes()); const scalar minPyrVol = relMinVol*minCellVol; + // Min area scalar minArea(readScalar(snapDict.lookup("minArea"))); - if (maxNonOrtho < 180.0-SMALL) + if (maxNonOrtho < 180.0 - SMALL) { Pout<< "Checking non orthogonality" << endl; label nOldSize = wrongFaces.size(); - mesh.setNonOrthThreshold(maxNonOrtho); mesh.checkFaceOrthogonality(false, &wrongFaces); Pout<< "Detected " << wrongFaces.size() - nOldSize @@ -123,7 +127,8 @@ void checkSnapMesh Pout<< "Checking face angles" << endl; label nOldSize = wrongFaces.size(); - mesh.checkFaceAngles(false, maxConcave, &wrongFaces); + mesh.checkFaceAngles(false, &wrongFaces); + Pout<< "Detected additional " << wrongFaces.size() - nOldSize << " faces with concavity > " << maxConcave << " degrees" << endl; @@ -144,6 +149,7 @@ void checkSnapMesh wrongFaces.insert(faceI); } } + Pout<< "Detected additional " << wrongFaces.size() - nOldSize << " faces with area < " << minArea << " m^2" << endl; } diff --git a/applications/utilities/mesh/advanced/refineHexMesh/refineHexMesh.C b/applications/utilities/mesh/advanced/refineHexMesh/refineHexMesh.C index d9cc48307..39c28bf28 100644 --- a/applications/utilities/mesh/advanced/refineHexMesh/refineHexMesh.C +++ b/applications/utilities/mesh/advanced/refineHexMesh/refineHexMesh.C @@ -164,7 +164,6 @@ int main(int argc, char *argv[]) // Update fields mesh.updateMesh(map); - pMesh.updateMesh(map); // Update numbering of cells/vertices. meshCutter.updateMesh(map); @@ -173,11 +172,11 @@ int main(int argc, char *argv[]) if (map().hasMotionPoints()) { mesh.movePoints(map().preMotionPoints()); - pMesh.movePoints(map().preMotionPoints()); } Pout<< "Refined from " << returnReduce(map().nOldCells(), sumOp