diff --git a/applications/solvers/combustion/chemFoam/solveChemistry.H b/applications/solvers/combustion/chemFoam/solveChemistry.H index daf56309bd..cb4258c31d 100644 --- a/applications/solvers/combustion/chemFoam/solveChemistry.H +++ b/applications/solvers/combustion/chemFoam/solveChemistry.H @@ -1,7 +1,3 @@ - dtChem = chemistry.solve - ( - runTime.value() - runTime.deltaT().value(), - runTime.deltaT().value() - ); + dtChem = chemistry.solve(runTime.deltaT().value()); scalar Sh = chemistry.Sh()()[0]/rho[0]; integratedHeat += Sh*runTime.deltaT().value(); diff --git a/applications/solvers/combustion/reactingFoam/Allwmake b/applications/solvers/combustion/reactingFoam/Allwmake index 0512e534f3..b60b5b6250 100755 --- a/applications/solvers/combustion/reactingFoam/Allwmake +++ b/applications/solvers/combustion/reactingFoam/Allwmake @@ -5,5 +5,6 @@ set -x wmake wmake rhoReactingFoam wmake rhoReactingBuoyantFoam +wmake LTSReactingFoam # ----------------------------------------------------------------- end-of-file diff --git a/applications/solvers/combustion/reactingFoam/EEqn.H b/applications/solvers/combustion/reactingFoam/EEqn.H index 56ddf867e7..9267c9a9be 100644 --- a/applications/solvers/combustion/reactingFoam/EEqn.H +++ b/applications/solvers/combustion/reactingFoam/EEqn.H @@ -16,7 +16,6 @@ : -dpdt ) - fvm::laplacian(turbulence->alphaEff(), he) -// - fvm::laplacian(turbulence->muEff(), he) // unit lewis no. == reaction->Sh() + fvOptions(rho, he) diff --git a/applications/solvers/combustion/reactingFoam/LTSReactingFoam/LTSReactingFoam.C b/applications/solvers/combustion/reactingFoam/LTSReactingFoam/LTSReactingFoam.C new file mode 100644 index 0000000000..0b34210865 --- /dev/null +++ b/applications/solvers/combustion/reactingFoam/LTSReactingFoam/LTSReactingFoam.C @@ -0,0 +1,107 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2013 OpenFOAM Foundation + \\/ 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 3 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, see . + +Application + LTSReactingFoam + +Description + Local time stepping (LTS) solver for steady, compressible, laminar or + turbulent reacting and non-reacting flow. + +\*---------------------------------------------------------------------------*/ + +#include "fvCFD.H" +#include "turbulenceModel.H" +#include "psiCombustionModel.H" +#include "multivariateScheme.H" +#include "pimpleControl.H" +#include "fvIOoptionList.H" +#include "fvcSmooth.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +int main(int argc, char *argv[]) +{ + #include "setRootCase.H" + #include "createTime.H" + #include "createMesh.H" + + pimpleControl pimple(mesh); + + #include "readGravitationalAcceleration.H" + #include "createFields.H" + #include "createFvOptions.H" + #include "initContinuityErrs.H" + #include "readTimeControls.H" + #include "compressibleCourantNo.H" + #include "setInitialrDeltaT.H" + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + + Info<< "\nStarting time loop\n" << endl; + + while (runTime.run()) + { + #include "readTimeControls.H" + + runTime++; + + Info<< "Time = " << runTime.timeName() << nl << endl; + + #include "setrDeltaT.H" + + #include "rhoEqn.H" + + // --- Pressure-velocity PIMPLE corrector loop + while (pimple.loop()) + { + #include "UEqn.H" + #include "YEqn.H" + #include "EEqn.H" + + // --- Pressure corrector loop + while (pimple.correct()) + { + #include "pEqn.H" + } + + if (pimple.turbCorr()) + { + turbulence->correct(); + } + } + + runTime.write(); + + Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s" + << " ClockTime = " << runTime.elapsedClockTime() << " s" + << nl << endl; + } + + Info<< "End\n" << endl; + + return(0); +} + + +// ************************************************************************* // diff --git a/applications/solvers/combustion/reactingFoam/LTSReactingFoam/Make/files b/applications/solvers/combustion/reactingFoam/LTSReactingFoam/Make/files new file mode 100644 index 0000000000..0b81e7a032 --- /dev/null +++ b/applications/solvers/combustion/reactingFoam/LTSReactingFoam/Make/files @@ -0,0 +1,3 @@ +LTSReactingFoam.C + +EXE = $(FOAM_APPBIN)/LTSReactingFoam diff --git a/applications/solvers/combustion/reactingFoam/LTSReactingFoam/Make/options b/applications/solvers/combustion/reactingFoam/LTSReactingFoam/Make/options new file mode 100644 index 0000000000..abac9d96a4 --- /dev/null +++ b/applications/solvers/combustion/reactingFoam/LTSReactingFoam/Make/options @@ -0,0 +1,28 @@ +EXE_INC = -ggdb3 \ + -I.. \ + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/fvOptions/lnInclude \ + -I$(LIB_SRC)/meshTools/lnInclude \ + -I$(LIB_SRC)/sampling/lnInclude \ + -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \ + -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \ + -I$(LIB_SRC)/ODE/lnInclude \ + -I$(LIB_SRC)/combustionModels/lnInclude + +EXE_LIBS = \ + -lfiniteVolume \ + -lfvOptions \ + -lmeshTools \ + -lsampling \ + -lcompressibleTurbulenceModel \ + -lcompressibleRASModels \ + -lcompressibleLESModels \ + -lreactionThermophysicalModels \ + -lspecie \ + -lfluidThermophysicalModels \ + -lchemistryModel \ + -lODE \ + -lcombustionModels diff --git a/applications/solvers/combustion/reactingFoam/LTSReactingFoam/readTimeControls.H b/applications/solvers/combustion/reactingFoam/LTSReactingFoam/readTimeControls.H new file mode 100644 index 0000000000..133f5b8165 --- /dev/null +++ b/applications/solvers/combustion/reactingFoam/LTSReactingFoam/readTimeControls.H @@ -0,0 +1,48 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2013 OpenFOAM Foundation + \\/ 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 3 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, see . + +\*---------------------------------------------------------------------------*/ + +// Maximum flow Courant number +scalar maxCo(readScalar(pimple.dict().lookup("maxCo"))); + +// Maximum time scale +scalar maxDeltaT(pimple.dict().lookupOrDefault("maxDeltaT", GREAT)); + +// Smoothing parameter (0-1) when smoothing iterations > 0 +scalar rDeltaTSmoothingCoeff +( + pimple.dict().lookupOrDefault("rDeltaTSmoothingCoeff", 1) +); + +// Damping coefficient (1-0) +scalar rDeltaTDampingCoeff +( + pimple.dict().lookupOrDefault("rDeltaTDampingCoeff", 1) +); + +// Maximum change in cell temperature per iteration (relative to previous value) +scalar alphaTemp(pimple.dict().lookupOrDefault("alphaTemp", 0.05)); + + +// ************************************************************************* // diff --git a/applications/solvers/combustion/reactingFoam/LTSReactingFoam/setInitialrDeltaT.H b/applications/solvers/combustion/reactingFoam/LTSReactingFoam/setInitialrDeltaT.H new file mode 100644 index 0000000000..c445eaee71 --- /dev/null +++ b/applications/solvers/combustion/reactingFoam/LTSReactingFoam/setInitialrDeltaT.H @@ -0,0 +1,14 @@ + volScalarField rDeltaT + ( + IOobject + ( + "rDeltaT", + runTime.timeName(), + mesh, + IOobject::READ_IF_PRESENT, + IOobject::AUTO_WRITE + ), + mesh, + dimensionedScalar("one", dimless/dimTime, 1), + zeroGradientFvPatchScalarField::typeName + ); diff --git a/applications/solvers/combustion/reactingFoam/LTSReactingFoam/setrDeltaT.H b/applications/solvers/combustion/reactingFoam/LTSReactingFoam/setrDeltaT.H new file mode 100644 index 0000000000..ce1caf1de5 --- /dev/null +++ b/applications/solvers/combustion/reactingFoam/LTSReactingFoam/setrDeltaT.H @@ -0,0 +1,97 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2013 OpenFOAM Foundation + \\/ 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 3 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, see . + +\*---------------------------------------------------------------------------*/ + +{ + Info<< "Time scales min/max:" << endl; + + // Cache old reciprocal time scale field + volScalarField rDeltaT0("rDeltaT0", rDeltaT); + + // Flow time scale + { + rDeltaT.dimensionedInternalField() = + ( + fvc::surfaceSum(mag(phi))().dimensionedInternalField() + /((2*maxCo)*mesh.V()*rho.dimensionedInternalField()) + ); + + // Limit the largest time scale + rDeltaT.max(1/maxDeltaT); + + Info<< " Flow = " + << gMin(1/rDeltaT.internalField()) << ", " + << gMax(1/rDeltaT.internalField()) << endl; + } + + // Reaction source time scale + if (alphaTemp < 1.0) + { + volScalarField::DimensionedInternalField rDeltaTT + ( + mag(reaction->Sh())/(alphaTemp*rho*thermo.Cp()*T) + ); + + Info<< " Temperature = " + << gMin(1/(rDeltaTT.field() + VSMALL)) << ", " + << gMax(1/(rDeltaTT.field() + VSMALL)) << endl; + + rDeltaT.dimensionedInternalField() = max + ( + rDeltaT.dimensionedInternalField(), + rDeltaTT + ); + } + + // Update tho boundary values of the reciprocal time-step + rDeltaT.correctBoundaryConditions(); + + // Spatially smooth the time scale field + if (rDeltaTSmoothingCoeff < 1.0) + { + fvc::smooth(rDeltaT, rDeltaTSmoothingCoeff); + } + + // Limit rate of change of time scale + // - reduce as much as required + // - only increase at a fraction of old time scale + if + ( + rDeltaTDampingCoeff < 1.0 + && runTime.timeIndex() > runTime.startTimeIndex() + 1 + ) + { + rDeltaT = max + ( + rDeltaT, + (scalar(1.0) - rDeltaTDampingCoeff)*rDeltaT0 + ); + } + + Info<< " Overall = " + << gMin(1/rDeltaT.internalField()) + << ", " << gMax(1/rDeltaT.internalField()) << endl; +} + +// ************************************************************************* // diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H index 14ba97d1c2..fd892bb932 100644 --- a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H +++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H @@ -115,5 +115,5 @@ K = 0.5*magSqr(U); if (thermo.dpdt()) { - dpdt = fvc::ddt(p); + dpdt = fvc::ddt(p) - fvc::div(fvc::meshPhi(rho, U), p); } diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/Allwclean b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/Allwclean new file mode 100755 index 0000000000..8098c05160 --- /dev/null +++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/Allwclean @@ -0,0 +1,8 @@ +#!/bin/sh +cd ${0%/*} || exit 1 # run from this directory +set -x + +wclean libso multiphaseMixtureThermo +wclean + +# ----------------------------------------------------------------- end-of-file diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/Allwmake b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/Allwmake new file mode 100755 index 0000000000..dd002ee069 --- /dev/null +++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/Allwmake @@ -0,0 +1,8 @@ +#!/bin/sh +cd ${0%/*} || exit 1 # run from this directory +set -x + +wmake libso multiphaseMixtureThermo +wmake + +# ----------------------------------------------------------------- end-of-file diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/Make/files b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/Make/files new file mode 100644 index 0000000000..8e6a4cb0de --- /dev/null +++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/Make/files @@ -0,0 +1,3 @@ +compressibleMultiphaseInterFoam.C + +EXE = $(FOAM_APPBIN)/compressibleMultiphaseInterFoam diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/Make/options b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/Make/options new file mode 100644 index 0000000000..0f47e6979e --- /dev/null +++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/Make/options @@ -0,0 +1,17 @@ +EXE_INC = \ + -ImultiphaseMixtureThermo/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \ + -I$(LIB_SRC)/transportModels/twoPhaseMixture/lnInclude \ + -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \ + -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \ + -I$(LIB_SRC)/finiteVolume/lnInclude + +EXE_LIBS = \ + -lmultiphaseMixtureThermo \ + -lfluidThermophysicalModels \ + -lspecie \ + -linterfaceProperties \ + -lcompressibleTurbulenceModel \ + -lcompressibleRASModels \ + -lcompressibleLESModels \ + -lfiniteVolume diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/TEqn.H b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/TEqn.H new file mode 100644 index 0000000000..a7ee3a7bb4 --- /dev/null +++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/TEqn.H @@ -0,0 +1,17 @@ +{ + fvScalarMatrix TEqn + ( + fvm::ddt(rho, T) + + fvm::div(multiphaseProperties.rhoPhi(), T) + - fvm::laplacian(multiphaseProperties.alphaEff(turbulence->mut()), T) + + ( + fvc::div(fvc::absolute(phi, U), p) + + fvc::ddt(rho, K) + fvc::div(multiphaseProperties.rhoPhi(), K) + )*multiphaseProperties.rCv() + ); + + TEqn.relax(); + TEqn.solve(); + + multiphaseProperties.correct(); +} diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/UEqn.H b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/UEqn.H new file mode 100644 index 0000000000..38cfebde6f --- /dev/null +++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/UEqn.H @@ -0,0 +1,27 @@ + fvVectorMatrix UEqn + ( + fvm::ddt(rho, U) + + fvm::div(multiphaseProperties.rhoPhi(), U) + + turbulence->divDevRhoReff(U) + ); + + UEqn.relax(); + + if (pimple.momentumPredictor()) + { + solve + ( + UEqn + == + fvc::reconstruct + ( + ( + multiphaseProperties.surfaceTensionForce() + - ghf*fvc::snGrad(rho) + - fvc::snGrad(p_rgh) + ) * mesh.magSf() + ) + ); + + K = 0.5*magSqr(U); + } diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/alphaCourantNo.H b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/alphaCourantNo.H new file mode 100644 index 0000000000..1fd68ffb79 --- /dev/null +++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/alphaCourantNo.H @@ -0,0 +1,57 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2013 OpenFOAM Foundation + \\/ 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 3 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, see . + +Global + CourantNo + +Description + Calculates and outputs the mean and maximum Courant Numbers. + +\*---------------------------------------------------------------------------*/ + +scalar maxAlphaCo +( + readScalar(runTime.controlDict().lookup("maxAlphaCo")) +); + +scalar alphaCoNum = 0.0; +scalar meanAlphaCoNum = 0.0; + +if (mesh.nInternalFaces()) +{ + scalarField sumPhi + ( + multiphaseProperties.nearInterface()().internalField() + * fvc::surfaceSum(mag(phi))().internalField() + ); + + alphaCoNum = 0.5*gMax(sumPhi/mesh.V().field())*runTime.deltaTValue(); + + meanAlphaCoNum = + 0.5*(gSum(sumPhi)/gSum(mesh.V().field()))*runTime.deltaTValue(); +} + +Info<< "Interface Courant Number mean: " << meanAlphaCoNum + << " max: " << alphaCoNum << endl; + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/compressibleMultiphaseInterFoam.C b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/compressibleMultiphaseInterFoam.C new file mode 100644 index 0000000000..09018b603f --- /dev/null +++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/compressibleMultiphaseInterFoam.C @@ -0,0 +1,111 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2013 OpenFOAM Foundation + \\/ 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 3 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, see . + +Application + compressibleMultiphaseInterFoam + +Description + Solver for n compressible, non-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. + +\*---------------------------------------------------------------------------*/ + +#include "fvCFD.H" +#include "multiphaseMixtureThermo.H" +#include "turbulenceModel.H" +#include "pimpleControl.H" +#include "fixedFluxPressureFvPatchScalarField.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +int main(int argc, char *argv[]) +{ + #include "setRootCase.H" + #include "createTime.H" + #include "createMesh.H" + #include "readGravitationalAcceleration.H" + + pimpleControl pimple(mesh); + + #include "readTimeControls.H" + #include "initContinuityErrs.H" + #include "createFields.H" + #include "CourantNo.H" + #include "setInitialDeltaT.H" + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + + Info<< "\nStarting time loop\n" << endl; + + while (runTime.run()) + { + #include "readTimeControls.H" + #include "CourantNo.H" + #include "alphaCourantNo.H" + #include "setDeltaT.H" + + runTime++; + + Info<< "Time = " << runTime.timeName() << nl << endl; + + // --- Pressure-velocity PIMPLE corrector loop + while (pimple.loop()) + { + multiphaseProperties.solve(); + + solve(fvm::ddt(rho) + fvc::div(multiphaseProperties.rhoPhi())); + + #include "UEqn.H" + #include "TEqn.H" + + // --- Pressure corrector loop + while (pimple.correct()) + { + #include "pEqn.H" + } + + if (pimple.turbCorr()) + { + turbulence->correct(); + } + } + + runTime.write(); + + Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s" + << " ClockTime = " << runTime.elapsedClockTime() << " s" + << nl << endl; + } + + Info<< "End\n" << endl; + + return 0; +} + + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/createFields.H b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/createFields.H new file mode 100644 index 0000000000..5c7f723a87 --- /dev/null +++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/createFields.H @@ -0,0 +1,71 @@ + Info<< "Reading field p_rgh\n" << endl; + volScalarField p_rgh + ( + IOobject + ( + "p_rgh", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh + ); + + Info<< "Reading field U\n" << endl; + volVectorField U + ( + IOobject + ( + "U", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh + ); + + #include "createPhi.H" + + Info<< "Constructing multiphaseMixtureThermo\n" << endl; + multiphaseMixtureThermo multiphaseProperties(U, phi); + + volScalarField& p = multiphaseProperties.p(); + volScalarField& T = multiphaseProperties.T(); + + volScalarField rho + ( + IOobject + ( + "rho", + runTime.timeName(), + mesh, + IOobject::READ_IF_PRESENT + ), + multiphaseProperties.rho() + ); + + Info<< max(rho) << min(rho); + + dimensionedScalar pMin(multiphaseProperties.lookup("pMin")); + + + Info<< "Calculating field g.h\n" << endl; + volScalarField gh("gh", g & mesh.C()); + surfaceScalarField ghf("ghf", g & mesh.Cf()); + + // Construct compressible turbulence model + autoPtr turbulence + ( + compressible::turbulenceModel::New + ( + rho, + U, + multiphaseProperties.rhoPhi(), + multiphaseProperties + ) + ); + + Info<< "Creating field kinetic energy K\n" << endl; + volScalarField K("K", 0.5*magSqr(U)); diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/Make/files b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/Make/files new file mode 100644 index 0000000000..93e6eb9e3d --- /dev/null +++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/Make/files @@ -0,0 +1,5 @@ +phaseModel/phaseModel.C +alphaContactAngle/alphaContactAngleFvPatchScalarField.C +multiphaseMixtureThermo.C + +LIB = $(FOAM_LIBBIN)/libmultiphaseMixtureThermo diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/Make/options b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/Make/options new file mode 100644 index 0000000000..eab8cce15d --- /dev/null +++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/Make/options @@ -0,0 +1,8 @@ +EXE_INC = \ + -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \ + -I$(LIB_SRC)/finiteVolume/lnInclude + +LIB_LIBS = \ + -lfluidThermophysicalModels \ + -lspecie \ + -lfiniteVolume diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/alphaContactAngle/alphaContactAngleFvPatchScalarField.C b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/alphaContactAngle/alphaContactAngleFvPatchScalarField.C new file mode 100644 index 0000000000..6872ae0321 --- /dev/null +++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/alphaContactAngle/alphaContactAngleFvPatchScalarField.C @@ -0,0 +1,146 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2013 OpenFOAM Foundation + \\/ 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 3 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, see . + +\*---------------------------------------------------------------------------*/ + +#include "alphaContactAngleFvPatchScalarField.H" +#include "addToRunTimeSelectionTable.H" +#include "fvPatchFieldMapper.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +alphaContactAngleFvPatchScalarField::interfaceThetaProps::interfaceThetaProps +( + Istream& is +) +: + theta0_(readScalar(is)), + uTheta_(readScalar(is)), + thetaA_(readScalar(is)), + thetaR_(readScalar(is)) +{} + + +Istream& operator>> +( + Istream& is, + alphaContactAngleFvPatchScalarField::interfaceThetaProps& tp +) +{ + is >> tp.theta0_ >> tp.uTheta_ >> tp.thetaA_ >> tp.thetaR_; + return is; +} + + +Ostream& operator<< +( + Ostream& os, + const alphaContactAngleFvPatchScalarField::interfaceThetaProps& tp +) +{ + os << tp.theta0_ << token::SPACE + << tp.uTheta_ << token::SPACE + << tp.thetaA_ << token::SPACE + << tp.thetaR_; + + return os; +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +alphaContactAngleFvPatchScalarField::alphaContactAngleFvPatchScalarField +( + const fvPatch& p, + const DimensionedField& iF +) +: + zeroGradientFvPatchScalarField(p, iF) +{} + + +alphaContactAngleFvPatchScalarField::alphaContactAngleFvPatchScalarField +( + const alphaContactAngleFvPatchScalarField& gcpsf, + const fvPatch& p, + const DimensionedField& iF, + const fvPatchFieldMapper& mapper +) +: + zeroGradientFvPatchScalarField(gcpsf, p, iF, mapper), + thetaProps_(gcpsf.thetaProps_) +{} + + +alphaContactAngleFvPatchScalarField::alphaContactAngleFvPatchScalarField +( + const fvPatch& p, + const DimensionedField& iF, + const dictionary& dict +) +: + zeroGradientFvPatchScalarField(p, iF), + thetaProps_(dict.lookup("thetaProperties")) +{ + evaluate(); +} + + +alphaContactAngleFvPatchScalarField::alphaContactAngleFvPatchScalarField +( + const alphaContactAngleFvPatchScalarField& gcpsf, + const DimensionedField& iF +) +: + zeroGradientFvPatchScalarField(gcpsf, iF), + thetaProps_(gcpsf.thetaProps_) +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +void alphaContactAngleFvPatchScalarField::write(Ostream& os) const +{ + fvPatchScalarField::write(os); + os.writeKeyword("thetaProperties") + << thetaProps_ << token::END_STATEMENT << nl; + writeEntry("value", os); +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +makePatchTypeField +( + fvPatchScalarField, + alphaContactAngleFvPatchScalarField +); + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/alphaContactAngle/alphaContactAngleFvPatchScalarField.H b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/alphaContactAngle/alphaContactAngleFvPatchScalarField.H new file mode 100644 index 0000000000..53ce61d24f --- /dev/null +++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/alphaContactAngle/alphaContactAngleFvPatchScalarField.H @@ -0,0 +1,215 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2013 OpenFOAM Foundation + \\/ 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 3 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, see . + +Class + Foam::alphaContactAngleFvPatchScalarField + +Description + Contact-angle boundary condition for multi-phase interface-capturing + simulations. Used in conjuction with multiphaseMixture. + +SourceFiles + alphaContactAngleFvPatchScalarField.C + +\*---------------------------------------------------------------------------*/ + +#ifndef alphaContactAngleFvPatchScalarField_H +#define alphaContactAngleFvPatchScalarField_H + +#include "zeroGradientFvPatchFields.H" +#include "multiphaseMixtureThermo.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class alphaContactAngleFvPatch Declaration +\*---------------------------------------------------------------------------*/ + +class alphaContactAngleFvPatchScalarField +: + public zeroGradientFvPatchScalarField +{ +public: + + class interfaceThetaProps + { + //- Equilibrium contact angle + scalar theta0_; + + //- Dynamic contact angle velocity scale + scalar uTheta_; + + //- Limiting advancing contact angle + scalar thetaA_; + + //- Limiting receeding contact angle + scalar thetaR_; + + + public: + + // Constructors + interfaceThetaProps() + {} + + interfaceThetaProps(Istream&); + + + // Member functions + + //- Return the equilibrium contact angle theta0 + scalar theta0(bool matched=true) const + { + if (matched) return theta0_; + else return 180.0 - theta0_; + } + + //- Return the dynamic contact angle velocity scale + scalar uTheta() const + { + return uTheta_; + } + + //- Return the limiting advancing contact angle + scalar thetaA(bool matched=true) const + { + if (matched) return thetaA_; + else return 180.0 - thetaA_; + } + + //- Return the limiting receeding contact angle + scalar thetaR(bool matched=true) const + { + if (matched) return thetaR_; + else return 180.0 - thetaR_; + } + + + // IO functions + + friend Istream& operator>>(Istream&, interfaceThetaProps&); + friend Ostream& operator<<(Ostream&, const interfaceThetaProps&); + }; + + typedef HashTable + < + interfaceThetaProps, + multiphaseMixtureThermo::interfacePair, + multiphaseMixtureThermo::interfacePair::hash + > thetaPropsTable; + + +private: + + // Private data + + thetaPropsTable thetaProps_; + + +public: + + //- Runtime type information + TypeName("alphaContactAngle"); + + + // Constructors + + //- Construct from patch and internal field + alphaContactAngleFvPatchScalarField + ( + const fvPatch&, + const DimensionedField& + ); + + //- Construct from patch, internal field and dictionary + alphaContactAngleFvPatchScalarField + ( + const fvPatch&, + const DimensionedField&, + const dictionary& + ); + + //- Construct by mapping given alphaContactAngleFvPatchScalarField + // onto a new patch + alphaContactAngleFvPatchScalarField + ( + const alphaContactAngleFvPatchScalarField&, + const fvPatch&, + const DimensionedField&, + const fvPatchFieldMapper& + ); + + //- Construct and return a clone + virtual tmp clone() const + { + return tmp + ( + new alphaContactAngleFvPatchScalarField(*this) + ); + } + + //- Construct as copy setting internal field reference + alphaContactAngleFvPatchScalarField + ( + const alphaContactAngleFvPatchScalarField&, + const DimensionedField& + ); + + //- Construct and return a clone setting internal field reference + virtual tmp clone + ( + const DimensionedField& iF + ) const + { + return tmp + ( + new alphaContactAngleFvPatchScalarField(*this, iF) + ); + } + + + // Member functions + + //- Return the contact angle properties + const thetaPropsTable& thetaProps() const + { + return thetaProps_; + } + + //- Write + virtual void write(Ostream&) const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.C b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.C new file mode 100644 index 0000000000..1639550564 --- /dev/null +++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.C @@ -0,0 +1,1129 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2013 OpenFOAM Foundation + \\/ 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 3 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, see . + +\*---------------------------------------------------------------------------*/ + +#include "multiphaseMixtureThermo.H" +#include "alphaContactAngleFvPatchScalarField.H" +#include "Time.H" +#include "subCycle.H" +#include "MULES.H" +#include "fvcDiv.H" +#include "fvcGrad.H" +#include "fvcSnGrad.H" +#include "fvcFlux.H" +#include "fvcMeshPhi.H" +#include "surfaceInterpolate.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ + defineTypeNameAndDebug(multiphaseMixtureThermo, 0); +} + + +const Foam::scalar Foam::multiphaseMixtureThermo::convertToRad = + Foam::constant::mathematical::pi/180.0; + + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +void Foam::multiphaseMixtureThermo::calcAlphas() +{ + scalar level = 0.0; + alphas_ == 0.0; + + forAllIter(PtrDictionary, phases_, phase) + { + alphas_ += level*phase(); + level += 1.0; + } + + alphas_.correctBoundaryConditions(); +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::multiphaseMixtureThermo::multiphaseMixtureThermo +( + const volVectorField& U, + const surfaceScalarField& phi +) +: + psiThermo(U.mesh(), word::null), + phases_(lookup("phases"), phaseModel::iNew(p_, T_)), + + mesh_(U.mesh()), + U_(U), + phi_(phi), + + rhoPhi_ + ( + IOobject + ( + "rho*phi", + mesh_.time().timeName(), + mesh_, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh_, + dimensionedScalar("rho*phi", dimMass/dimTime, 0.0) + ), + + alphas_ + ( + IOobject + ( + "alphas", + mesh_.time().timeName(), + mesh_, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + mesh_, + dimensionedScalar("alphas", dimless, 0.0), + zeroGradientFvPatchScalarField::typeName + ), + + sigmas_(lookup("sigmas")), + dimSigma_(1, 0, -2, 0, 0), + deltaN_ + ( + "deltaN", + 1e-8/pow(average(mesh_.V()), 1.0/3.0) + ) +{ + calcAlphas(); + alphas_.write(); + correct(); +} + + +// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * // + +void Foam::multiphaseMixtureThermo::correct() +{ + forAllIter(PtrDictionary, phases_, phasei) + { + phasei().correct(); + } + + PtrDictionary::iterator phasei = phases_.begin(); + + psi_ = phasei()*phasei().thermo().psi(); + mu_ = phasei()*phasei().thermo().mu(); + alpha_ = phasei()*phasei().thermo().alpha(); + + for (++phasei; phasei != phases_.end(); ++phasei) + { + psi_ += phasei()*phasei().thermo().psi(); + mu_ += phasei()*phasei().thermo().mu(); + alpha_ += phasei()*phasei().thermo().alpha(); + } +} + + +void Foam::multiphaseMixtureThermo::correctRho(const volScalarField& dp) +{ + forAllIter(PtrDictionary, phases_, phasei) + { + phasei().thermo().rho() += phasei().thermo().psi()*dp; + } +} + + +bool Foam::multiphaseMixtureThermo::incompressible() const +{ + bool ico = true; + + forAllConstIter(PtrDictionary, phases_, phase) + { + ico &= phase().thermo().incompressible(); + } + + return ico; +} + + +bool Foam::multiphaseMixtureThermo::isochoric() const +{ + bool iso = true; + + forAllConstIter(PtrDictionary, phases_, phase) + { + iso &= phase().thermo().incompressible(); + } + + return iso; +} + + +Foam::tmp Foam::multiphaseMixtureThermo::he +( + const volScalarField& p, + const volScalarField& T +) const +{ + PtrDictionary::const_iterator phasei = phases_.begin(); + + tmp the(phasei()*phasei().thermo().he(p, T)); + + for (++phasei; phasei != phases_.end(); ++phasei) + { + the() += phasei()*phasei().thermo().he(p, T); + } + + return the; +} + + +Foam::tmp Foam::multiphaseMixtureThermo::he +( + const scalarField& p, + const scalarField& T, + const labelList& cells +) const +{ + PtrDictionary::const_iterator phasei = phases_.begin(); + + tmp the + ( + scalarField(phasei(), cells)*phasei().thermo().he(p, T, cells) + ); + + for (++phasei; phasei != phases_.end(); ++phasei) + { + the() += scalarField(phasei(), cells)*phasei().thermo().he(p, T, cells); + } + + return the; +} + + +Foam::tmp Foam::multiphaseMixtureThermo::he +( + const scalarField& p, + const scalarField& T, + const label patchi +) const +{ + PtrDictionary::const_iterator phasei = phases_.begin(); + + tmp the + ( + phasei().boundaryField()[patchi]*phasei().thermo().he(p, T, patchi) + ); + + for (++phasei; phasei != phases_.end(); ++phasei) + { + the() += + phasei().boundaryField()[patchi]*phasei().thermo().he(p, T, patchi); + } + + return the; +} + + +Foam::tmp Foam::multiphaseMixtureThermo::hc() const +{ + PtrDictionary::const_iterator phasei = phases_.begin(); + + tmp thc(phasei()*phasei().thermo().hc()); + + for (++phasei; phasei != phases_.end(); ++phasei) + { + thc() += phasei()*phasei().thermo().hc(); + } + + return thc; +} + + +Foam::tmp Foam::multiphaseMixtureThermo::THE +( + const scalarField& h, + const scalarField& p, + const scalarField& T0, + const labelList& cells +) const +{ + notImplemented("multiphaseMixtureThermo::THE(...)"); + return T0; +} + + +Foam::tmp Foam::multiphaseMixtureThermo::THE +( + const scalarField& h, + const scalarField& p, + const scalarField& T0, + const label patchi +) const +{ + notImplemented("multiphaseMixtureThermo::THE(...)"); + return T0; +} + + +Foam::tmp Foam::multiphaseMixtureThermo::rho() const +{ + PtrDictionary::const_iterator phasei = phases_.begin(); + + tmp trho(phasei()*phasei().thermo().rho()); + + for (++phasei; phasei != phases_.end(); ++phasei) + { + trho() += phasei()*phasei().thermo().rho(); + } + + return trho; +} + + +Foam::tmp Foam::multiphaseMixtureThermo::Cp() const +{ + PtrDictionary::const_iterator phasei = phases_.begin(); + + tmp tCp(phasei()*phasei().thermo().Cp()); + + for (++phasei; phasei != phases_.end(); ++phasei) + { + tCp() += phasei()*phasei().thermo().Cp(); + } + + return tCp; +} + + +Foam::tmp Foam::multiphaseMixtureThermo::Cp +( + const scalarField& p, + const scalarField& T, + const label patchi +) const +{ + PtrDictionary::const_iterator phasei = phases_.begin(); + + tmp tCp + ( + phasei().boundaryField()[patchi]*phasei().thermo().Cp(p, T, patchi) + ); + + for (++phasei; phasei != phases_.end(); ++phasei) + { + tCp() += + phasei().boundaryField()[patchi]*phasei().thermo().Cp(p, T, patchi); + } + + return tCp; +} + + +Foam::tmp Foam::multiphaseMixtureThermo::Cv() const +{ + PtrDictionary::const_iterator phasei = phases_.begin(); + + tmp tCv(phasei()*phasei().thermo().Cv()); + + for (++phasei; phasei != phases_.end(); ++phasei) + { + tCv() += phasei()*phasei().thermo().Cv(); + } + + return tCv; +} + + +Foam::tmp Foam::multiphaseMixtureThermo::Cv +( + const scalarField& p, + const scalarField& T, + const label patchi +) const +{ + PtrDictionary::const_iterator phasei = phases_.begin(); + + tmp tCv + ( + phasei().boundaryField()[patchi]*phasei().thermo().Cv(p, T, patchi) + ); + + for (++phasei; phasei != phases_.end(); ++phasei) + { + tCv() += + phasei().boundaryField()[patchi]*phasei().thermo().Cv(p, T, patchi); + } + + return tCv; +} + + +Foam::tmp Foam::multiphaseMixtureThermo::gamma() const +{ + PtrDictionary::const_iterator phasei = phases_.begin(); + + tmp tgamma(phasei()*phasei().thermo().gamma()); + + for (++phasei; phasei != phases_.end(); ++phasei) + { + tgamma() += phasei()*phasei().thermo().gamma(); + } + + return tgamma; +} + + +Foam::tmp Foam::multiphaseMixtureThermo::gamma +( + const scalarField& p, + const scalarField& T, + const label patchi +) const +{ + PtrDictionary::const_iterator phasei = phases_.begin(); + + tmp tgamma + ( + phasei().boundaryField()[patchi]*phasei().thermo().gamma(p, T, patchi) + ); + + for (++phasei; phasei != phases_.end(); ++phasei) + { + tgamma() += + phasei().boundaryField()[patchi] + *phasei().thermo().gamma(p, T, patchi); + } + + return tgamma; +} + + +Foam::tmp Foam::multiphaseMixtureThermo::Cpv() const +{ + PtrDictionary::const_iterator phasei = phases_.begin(); + + tmp tCpv(phasei()*phasei().thermo().Cpv()); + + for (++phasei; phasei != phases_.end(); ++phasei) + { + tCpv() += phasei()*phasei().thermo().Cpv(); + } + + return tCpv; +} + + +Foam::tmp Foam::multiphaseMixtureThermo::Cpv +( + const scalarField& p, + const scalarField& T, + const label patchi +) const +{ + PtrDictionary::const_iterator phasei = phases_.begin(); + + tmp tCpv + ( + phasei().boundaryField()[patchi]*phasei().thermo().Cpv(p, T, patchi) + ); + + for (++phasei; phasei != phases_.end(); ++phasei) + { + tCpv() += + phasei().boundaryField()[patchi] + *phasei().thermo().Cpv(p, T, patchi); + } + + return tCpv; +} + + +Foam::tmp Foam::multiphaseMixtureThermo::CpByCpv() const +{ + PtrDictionary::const_iterator phasei = phases_.begin(); + + tmp tCpByCpv(phasei()*phasei().thermo().CpByCpv()); + + for (++phasei; phasei != phases_.end(); ++phasei) + { + tCpByCpv() += phasei()*phasei().thermo().CpByCpv(); + } + + return tCpByCpv; +} + + +Foam::tmp Foam::multiphaseMixtureThermo::CpByCpv +( + const scalarField& p, + const scalarField& T, + const label patchi +) const +{ + PtrDictionary::const_iterator phasei = phases_.begin(); + + tmp tCpByCpv + ( + phasei().boundaryField()[patchi]*phasei().thermo().CpByCpv(p, T, patchi) + ); + + for (++phasei; phasei != phases_.end(); ++phasei) + { + tCpByCpv() += + phasei().boundaryField()[patchi] + *phasei().thermo().CpByCpv(p, T, patchi); + } + + return tCpByCpv; +} + + +Foam::tmp Foam::multiphaseMixtureThermo::kappa() const +{ + PtrDictionary::const_iterator phasei = phases_.begin(); + + tmp tkappa(phasei()*phasei().thermo().kappa()); + + for (++phasei; phasei != phases_.end(); ++phasei) + { + tkappa() += phasei()*phasei().thermo().kappa(); + } + + return tkappa; +} + + +Foam::tmp Foam::multiphaseMixtureThermo::kappa +( + const label patchi +) const +{ + PtrDictionary::const_iterator phasei = phases_.begin(); + + tmp tkappa + ( + phasei().boundaryField()[patchi]*phasei().thermo().kappa(patchi) + ); + + for (++phasei; phasei != phases_.end(); ++phasei) + { + tkappa() += + phasei().boundaryField()[patchi]*phasei().thermo().kappa(patchi); + } + + return tkappa; +} + + +Foam::tmp Foam::multiphaseMixtureThermo::kappaEff +( + const volScalarField& alphat +) const +{ + PtrDictionary::const_iterator phasei = phases_.begin(); + + tmp tkappaEff(phasei()*phasei().thermo().kappaEff(alphat)); + + for (++phasei; phasei != phases_.end(); ++phasei) + { + tkappaEff() += phasei()*phasei().thermo().kappaEff(alphat); + } + + return tkappaEff; +} + + +Foam::tmp Foam::multiphaseMixtureThermo::kappaEff +( + const scalarField& alphat, + const label patchi +) const +{ + PtrDictionary::const_iterator phasei = phases_.begin(); + + tmp tkappaEff + ( + phasei().boundaryField()[patchi] + *phasei().thermo().kappaEff(alphat, patchi) + ); + + for (++phasei; phasei != phases_.end(); ++phasei) + { + tkappaEff() += + phasei().boundaryField()[patchi] + *phasei().thermo().kappaEff(alphat, patchi); + } + + return tkappaEff; +} + + +Foam::tmp Foam::multiphaseMixtureThermo::alphaEff +( + const volScalarField& alphat +) const +{ + PtrDictionary::const_iterator phasei = phases_.begin(); + + tmp talphaEff(phasei()*phasei().thermo().alphaEff(alphat)); + + for (++phasei; phasei != phases_.end(); ++phasei) + { + talphaEff() += phasei()*phasei().thermo().alphaEff(alphat); + } + + return talphaEff; +} + + +Foam::tmp Foam::multiphaseMixtureThermo::alphaEff +( + const scalarField& alphat, + const label patchi +) const +{ + PtrDictionary::const_iterator phasei = phases_.begin(); + + tmp talphaEff + ( + phasei().boundaryField()[patchi] + *phasei().thermo().alphaEff(alphat, patchi) + ); + + for (++phasei; phasei != phases_.end(); ++phasei) + { + talphaEff() += + phasei().boundaryField()[patchi] + *phasei().thermo().alphaEff(alphat, patchi); + } + + return talphaEff; +} + + +Foam::tmp Foam::multiphaseMixtureThermo::rCv() const +{ + PtrDictionary::const_iterator phasei = phases_.begin(); + + tmp trCv(phasei()/phasei().thermo().Cv()); + + for (++phasei; phasei != phases_.end(); ++phasei) + { + trCv() += phasei()/phasei().thermo().Cv(); + } + + return trCv; +} + + +Foam::tmp +Foam::multiphaseMixtureThermo::surfaceTensionForce() const +{ + tmp tstf + ( + new surfaceScalarField + ( + IOobject + ( + "surfaceTensionForce", + mesh_.time().timeName(), + mesh_ + ), + mesh_, + dimensionedScalar + ( + "surfaceTensionForce", + dimensionSet(1, -2, -2, 0, 0), + 0.0 + ) + ) + ); + + surfaceScalarField& stf = tstf(); + + forAllConstIter(PtrDictionary, phases_, phase1) + { + const phaseModel& alpha1 = phase1(); + + PtrDictionary::const_iterator phase2 = phase1; + ++phase2; + + for (; phase2 != phases_.end(); ++phase2) + { + const phaseModel& alpha2 = phase2(); + + sigmaTable::const_iterator sigma = + sigmas_.find(interfacePair(alpha1, alpha2)); + + if (sigma == sigmas_.end()) + { + FatalErrorIn + ( + "multiphaseMixtureThermo::surfaceTensionForce() const" + ) << "Cannot find interface " << interfacePair(alpha1, alpha2) + << " in list of sigma values" + << exit(FatalError); + } + + stf += dimensionedScalar("sigma", dimSigma_, sigma()) + *fvc::interpolate(K(alpha1, alpha2))* + ( + fvc::interpolate(alpha2)*fvc::snGrad(alpha1) + - fvc::interpolate(alpha1)*fvc::snGrad(alpha2) + ); + } + } + + return tstf; +} + + +void Foam::multiphaseMixtureThermo::solve() +{ + const Time& runTime = mesh_.time(); + + const dictionary& pimpleDict = mesh_.solutionDict().subDict("PIMPLE"); + + label nAlphaSubCycles(readLabel(pimpleDict.lookup("nAlphaSubCycles"))); + + scalar cAlpha(readScalar(pimpleDict.lookup("cAlpha"))); + + + volScalarField& alpha = phases_.first(); + + if (nAlphaSubCycles > 1) + { + surfaceScalarField rhoPhiSum(0.0*rhoPhi_); + dimensionedScalar totalDeltaT = runTime.deltaT(); + + for + ( + subCycle alphaSubCycle(alpha, nAlphaSubCycles); + !(++alphaSubCycle).end(); + ) + { + solveAlphas(cAlpha); + rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi_; + } + + rhoPhi_ = rhoPhiSum; + } + else + { + solveAlphas(cAlpha); + } +} + + +Foam::tmp Foam::multiphaseMixtureThermo::nHatfv +( + const volScalarField& alpha1, + const volScalarField& alpha2 +) const +{ + /* + // Cell gradient of alpha + volVectorField gradAlpha = + alpha2*fvc::grad(alpha1) - alpha1*fvc::grad(alpha2); + + // Interpolated face-gradient of alpha + surfaceVectorField gradAlphaf = fvc::interpolate(gradAlpha); + */ + + surfaceVectorField gradAlphaf + ( + fvc::interpolate(alpha2)*fvc::interpolate(fvc::grad(alpha1)) + - fvc::interpolate(alpha1)*fvc::interpolate(fvc::grad(alpha2)) + ); + + // Face unit interface normal + return gradAlphaf/(mag(gradAlphaf) + deltaN_); +} + + +Foam::tmp Foam::multiphaseMixtureThermo::nHatf +( + const volScalarField& alpha1, + const volScalarField& alpha2 +) const +{ + // Face unit interface normal flux + return nHatfv(alpha1, alpha2) & mesh_.Sf(); +} + + +// 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::multiphaseMixtureThermo::correctContactAngle +( + const phaseModel& alpha1, + const phaseModel& alpha2, + surfaceVectorField::GeometricBoundaryField& nHatb +) const +{ + const volScalarField::GeometricBoundaryField& gbf + = alpha1.boundaryField(); + + const fvBoundaryMesh& boundary = mesh_.boundary(); + + forAll(boundary, patchi) + { + if (isA(gbf[patchi])) + { + const alphaContactAngleFvPatchScalarField& acap = + refCast(gbf[patchi]); + + vectorField& nHatPatch = nHatb[patchi]; + + vectorField AfHatPatch + ( + mesh_.Sf().boundaryField()[patchi] + /mesh_.magSf().boundaryField()[patchi] + ); + + alphaContactAngleFvPatchScalarField::thetaPropsTable:: + const_iterator tp = + acap.thetaProps().find(interfacePair(alpha1, alpha2)); + + if (tp == acap.thetaProps().end()) + { + FatalErrorIn + ( + "multiphaseMixtureThermo::correctContactAngle" + "(const phaseModel& alpha1, const phaseModel& alpha2, " + "fvPatchVectorFieldField& nHatb) const" + ) << "Cannot find interface " << interfacePair(alpha1, alpha2) + << "\n in table of theta properties for patch " + << acap.patch().name() + << exit(FatalError); + } + + bool matched = (tp.key().first() == alpha1.name()); + + scalar theta0 = convertToRad*tp().theta0(matched); + scalarField theta(boundary[patchi].size(), theta0); + + scalar uTheta = tp().uTheta(); + + // Calculate the dynamic contact angle if required + if (uTheta > SMALL) + { + scalar thetaA = convertToRad*tp().thetaA(matched); + scalar thetaR = convertToRad*tp().thetaR(matched); + + // Calculated the component of the velocity parallel to the wall + vectorField Uwall + ( + U_.boundaryField()[patchi].patchInternalField() + - U_.boundaryField()[patchi] + ); + Uwall -= (AfHatPatch & Uwall)*AfHatPatch; + + // Find the direction of the interface parallel to the wall + vectorField nWall + ( + nHatPatch - (AfHatPatch & nHatPatch)*AfHatPatch + ); + + // Normalise nWall + nWall /= (mag(nWall) + SMALL); + + // Calculate Uwall resolved normal to the interface parallel to + // the interface + scalarField uwall(nWall & Uwall); + + theta += (thetaA - thetaR)*tanh(uwall/uTheta); + } + + + // Reset nHatPatch to correspond to the contact angle + + scalarField a12(nHatPatch & AfHatPatch); + + scalarField b1(cos(theta)); + + scalarField b2(nHatPatch.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); + + nHatPatch = a*AfHatPatch + b*nHatPatch; + + nHatPatch /= (mag(nHatPatch) + deltaN_.value()); + } + } +} + + +Foam::tmp Foam::multiphaseMixtureThermo::K +( + const phaseModel& alpha1, + const phaseModel& alpha2 +) const +{ + tmp tnHatfv = nHatfv(alpha1, alpha2); + + correctContactAngle(alpha1, alpha2, tnHatfv().boundaryField()); + + // Simple expression for curvature + return -fvc::div(tnHatfv & mesh_.Sf()); +} + + +Foam::tmp +Foam::multiphaseMixtureThermo::nearInterface() const +{ + tmp tnearInt + ( + new volScalarField + ( + IOobject + ( + "nearInterface", + mesh_.time().timeName(), + mesh_ + ), + mesh_, + dimensionedScalar("nearInterface", dimless, 0.0) + ) + ); + + forAllConstIter(PtrDictionary, phases_, phase) + { + tnearInt() = max(tnearInt(), pos(phase() - 0.01)*pos(0.99 - phase())); + } + + return tnearInt; +} + + +void Foam::multiphaseMixtureThermo::solveAlphas +( + const scalar cAlpha +) +{ + static label nSolves=-1; + nSolves++; + + word alphaScheme("div(phi,alpha)"); + word alpharScheme("div(phirb,alpha)"); + + surfaceScalarField phic(mag(phi_/mesh_.magSf())); + phic = min(cAlpha*phic, max(phic)); + + PtrList phiAlphaCorrs(phases_.size()); + int phasei = 0; + + forAllIter(PtrDictionary, phases_, phase) + { + phaseModel& alpha = phase(); + + phiAlphaCorrs.set + ( + phasei, + new surfaceScalarField + ( + fvc::flux + ( + phi_, + alpha, + alphaScheme + ) + ) + ); + + surfaceScalarField& phiAlphaCorr = phiAlphaCorrs[phasei]; + + forAllIter(PtrDictionary, phases_, phase2) + { + phaseModel& alpha2 = phase2(); + + if (&alpha2 == &alpha) continue; + + surfaceScalarField phir(phic*nHatf(alpha, alpha2)); + + phiAlphaCorr += fvc::flux + ( + -fvc::flux(-phir, alpha2, alpharScheme), + alpha, + alpharScheme + ); + } + + MULES::limit + ( + geometricOneField(), + alpha, + phi_, + phiAlphaCorr, + zeroField(), + zeroField(), + 1, + 0, + 3, + true + ); + + phasei++; + } + + MULES::limitSum(phiAlphaCorrs); + + rhoPhi_ = dimensionedScalar("0", dimensionSet(1, 0, -1, 0, 0), 0); + + volScalarField sumAlpha + ( + IOobject + ( + "sumAlpha", + mesh_.time().timeName(), + mesh_ + ), + mesh_, + dimensionedScalar("sumAlpha", dimless, 0) + ); + + + volScalarField divU(fvc::div(fvc::absolute(phi_, U_))); + + + phasei = 0; + + forAllIter(PtrDictionary, phases_, phase) + { + phaseModel& alpha = phase(); + + surfaceScalarField& phiAlpha = phiAlphaCorrs[phasei]; + phiAlpha += upwind(mesh_, phi_).flux(alpha); + + volScalarField::DimensionedInternalField Sp + ( + IOobject + ( + "Sp", + mesh_.time().timeName(), + mesh_ + ), + mesh_, + dimensionedScalar("Sp", alpha.dgdt().dimensions(), 0.0) + ); + + volScalarField::DimensionedInternalField Su + ( + IOobject + ( + "Su", + mesh_.time().timeName(), + mesh_ + ), + // Divergence term is handled explicitly to be + // consistent with the explicit transport solution + divU*min(alpha, scalar(1)) + ); + + { + const scalarField& dgdt = alpha.dgdt(); + + forAll(dgdt, celli) + { + if (dgdt[celli] < 0.0 && alpha[celli] > 0.0) + { + Sp[celli] += dgdt[celli]*alpha[celli]; + Su[celli] -= dgdt[celli]*alpha[celli]; + } + else if (dgdt[celli] > 0.0 && alpha[celli] < 1.0) + { + Sp[celli] -= dgdt[celli]*(1.0 - alpha[celli]); + } + } + } + + forAllConstIter(PtrDictionary, phases_, phase2) + { + const phaseModel& alpha2 = phase2(); + + if (&alpha2 == &alpha) continue; + + const scalarField& dgdt2 = alpha2.dgdt(); + + forAll(dgdt2, celli) + { + if (dgdt2[celli] > 0.0 && alpha2[celli] < 1.0) + { + Sp[celli] -= dgdt2[celli]*(1.0 - alpha2[celli]); + Su[celli] += dgdt2[celli]*alpha[celli]; + } + else if (dgdt2[celli] < 0.0 && alpha2[celli] > 0.0) + { + Sp[celli] += dgdt2[celli]*alpha2[celli]; + } + } + } + + MULES::explicitSolve + ( + geometricOneField(), + alpha, + phiAlpha, + Sp, + Su + ); + + rhoPhi_ += fvc::interpolate(alpha.thermo().rho())*phiAlpha; + + Info<< alpha.name() << " volume fraction, min, max = " + << alpha.weightedAverage(mesh_.V()).value() + << ' ' << min(alpha).value() + << ' ' << max(alpha).value() + << endl; + + sumAlpha += alpha; + + phasei++; + } + + Info<< "Phase-sum volume fraction, min, max = " + << sumAlpha.weightedAverage(mesh_.V()).value() + << ' ' << min(sumAlpha).value() + << ' ' << max(sumAlpha).value() + << endl; + + calcAlphas(); +} + + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.H b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.H new file mode 100644 index 0000000000..546048894f --- /dev/null +++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/multiphaseMixtureThermo.H @@ -0,0 +1,434 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2013 OpenFOAM Foundation + \\/ 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 3 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, see . + +Class + Foam::multiphaseMixtureThermo + +Description + +SourceFiles + multiphaseMixtureThermo.C + +\*---------------------------------------------------------------------------*/ + +#ifndef multiphaseMixtureThermo_H +#define multiphaseMixtureThermo_H + +#include "phaseModel.H" +#include "PtrDictionary.H" +#include "volFields.H" +#include "surfaceFields.H" +#include "rhoThermo.H" +#include "psiThermo.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class multiphaseMixtureThermo Declaration +\*---------------------------------------------------------------------------*/ + +class multiphaseMixtureThermo +: + public psiThermo +{ +public: + + class interfacePair + : + public Pair + { + public: + + class hash + : + public Hash + { + public: + + hash() + {} + + label operator()(const interfacePair& key) const + { + return word::hash()(key.first()) + word::hash()(key.second()); + } + }; + + + // Constructors + + interfacePair() + {} + + interfacePair(const word& alpha1Name, const word& alpha2Name) + : + Pair(alpha1Name, alpha2Name) + {} + + interfacePair(const phaseModel& alpha1, const phaseModel& alpha2) + : + Pair(alpha1.name(), alpha2.name()) + {} + + + // Friend Operators + + friend bool operator== + ( + const interfacePair& a, + const interfacePair& b + ) + { + return + ( + ((a.first() == b.first()) && (a.second() == b.second())) + || ((a.first() == b.second()) && (a.second() == b.first())) + ); + } + + friend bool operator!= + ( + const interfacePair& a, + const interfacePair& b + ) + { + return (!(a == b)); + } + }; + + +private: + + // Private data + + //- Dictionary of phases + PtrDictionary phases_; + + const fvMesh& mesh_; + const volVectorField& U_; + const surfaceScalarField& phi_; + + surfaceScalarField rhoPhi_; + + volScalarField alphas_; + + typedef HashTable + sigmaTable; + + sigmaTable sigmas_; + dimensionSet dimSigma_; + + //- Stabilisation for normalisation of the interface normal + const dimensionedScalar deltaN_; + + //- Conversion factor for degrees into radians + static const scalar convertToRad; + + + // Private member functions + + void calcAlphas(); + + void solveAlphas(const scalar cAlpha); + + tmp nHatfv + ( + const volScalarField& alpha1, + const volScalarField& alpha2 + ) const; + + tmp nHatf + ( + const volScalarField& alpha1, + const volScalarField& alpha2 + ) const; + + void correctContactAngle + ( + const phaseModel& alpha1, + const phaseModel& alpha2, + surfaceVectorField::GeometricBoundaryField& nHatb + ) const; + + tmp K + ( + const phaseModel& alpha1, + const phaseModel& alpha2 + ) const; + + +public: + + //- Runtime type information + TypeName("multiphaseMixtureThermo"); + + + // Constructors + + //- Construct from components + multiphaseMixtureThermo + ( + const volVectorField& U, + const surfaceScalarField& phi + ); + + + //- Destructor + virtual ~multiphaseMixtureThermo() + {} + + + // Member Functions + + //- Return the phases + const PtrDictionary& phases() const + { + return phases_; + } + + //- Return non-const access to the phases + PtrDictionary& phases() + { + return phases_; + } + + //- Return the velocity + const volVectorField& U() const + { + return U_; + } + + //- Return the volumetric flux + const surfaceScalarField& phi() const + { + return phi_; + } + + const surfaceScalarField& rhoPhi() const + { + return rhoPhi_; + } + + //- Update properties + virtual void correct(); + + //- Update densities for given pressure change + void correctRho(const volScalarField& dp); + + //- Return true if the equation of state is incompressible + // i.e. rho != f(p) + virtual bool incompressible() const; + + //- Return true if the equation of state is isochoric + // i.e. rho = const + virtual bool isochoric() const; + + + // Access to thermodynamic state variables + + //- Enthalpy/Internal energy [J/kg] + // Non-const access allowed for transport equations + virtual volScalarField& he() + { + notImplemented("multiphaseMixtureThermo::he()"); + return phases_[0]->thermo().he(); + } + + //- Enthalpy/Internal energy [J/kg] + virtual const volScalarField& he() const + { + notImplemented("multiphaseMixtureThermo::he() const"); + return phases_[0]->thermo().he(); + } + + //- Enthalpy/Internal energy + // for given pressure and temperature [J/kg] + virtual tmp he + ( + const volScalarField& p, + const volScalarField& T + ) const; + + //- Enthalpy/Internal energy for cell-set [J/kg] + virtual tmp he + ( + const scalarField& p, + const scalarField& T, + const labelList& cells + ) const; + + //- Enthalpy/Internal energy for patch [J/kg] + virtual tmp he + ( + const scalarField& p, + const scalarField& T, + const label patchi + ) const; + + //- Chemical enthalpy [J/kg] + virtual tmp hc() const; + + //- Temperature from enthalpy/internal energy for cell-set + virtual tmp THE + ( + const scalarField& h, + const scalarField& p, + const scalarField& T0, // starting temperature + const labelList& cells + ) const; + + //- Temperature from enthalpy/internal energy for patch + virtual tmp THE + ( + const scalarField& h, + const scalarField& p, + const scalarField& T0, // starting temperature + const label patchi + ) const; + + + // Fields derived from thermodynamic state variables + + //- Density [kg/m^3] + virtual tmp rho() const; + + //- Heat capacity at constant pressure [J/kg/K] + virtual tmp Cp() const; + + //- Heat capacity at constant pressure for patch [J/kg/K] + virtual tmp Cp + ( + const scalarField& p, + const scalarField& T, + const label patchi + ) const; + + //- Heat capacity at constant volume [J/kg/K] + virtual tmp Cv() const; + + //- Heat capacity at constant volume for patch [J/kg/K] + virtual tmp Cv + ( + const scalarField& p, + const scalarField& T, + const label patchi + ) const; + + //- gamma = Cp/Cv [] + virtual tmp gamma() const; + + //- gamma = Cp/Cv for patch [] + virtual tmp gamma + ( + const scalarField& p, + const scalarField& T, + const label patchi + ) const; + + //- Heat capacity at constant pressure/volume [J/kg/K] + virtual tmp Cpv() const; + + //- Heat capacity at constant pressure/volume for patch [J/kg/K] + virtual tmp Cpv + ( + const scalarField& p, + const scalarField& T, + const label patchi + ) const; + + //- Heat capacity ratio [] + virtual tmp CpByCpv() const; + + //- Heat capacity ratio for patch [] + virtual tmp CpByCpv + ( + const scalarField& p, + const scalarField& T, + const label patchi + ) const; + + + // Fields derived from transport state variables + + //- Thermal diffusivity for temperature of mixture [J/m/s/K] + virtual tmp kappa() const; + + //- Thermal diffusivity of mixture for patch [J/m/s/K] + virtual tmp kappa + ( + const label patchi + ) const; + + //- Effective thermal diffusivity of mixture [J/m/s/K] + virtual tmp kappaEff + ( + const volScalarField& alphat + ) const; + + //- Effective thermal diffusivity of mixture for patch [J/m/s/K] + virtual tmp kappaEff + ( + const scalarField& alphat, + const label patchi + ) const; + + //- Effective thermal diffusivity of mixture [J/m/s/K] + virtual tmp alphaEff + ( + const volScalarField& alphat + ) const; + + //- Effective thermal diffusivity of mixture for patch [J/m/s/K] + virtual tmp alphaEff + ( + const scalarField& alphat, + const label patchi + ) const; + + + //- Return the phase-averaged reciprocal Cv + tmp rCv() const; + + tmp surfaceTensionForce() const; + + //- Indicator of the proximity of the interface + // Field values are 1 near and 0 away for the interface. + tmp nearInterface() const; + + //- Solve for the mixture phase-fractions + void solve(); +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/phaseModel/phaseModel.C b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/phaseModel/phaseModel.C new file mode 100644 index 0000000000..1559f25a48 --- /dev/null +++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/phaseModel/phaseModel.C @@ -0,0 +1,95 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2013 OpenFOAM Foundation + \\/ 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 3 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, see . + +\*---------------------------------------------------------------------------*/ + +#include "phaseModel.H" + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::phaseModel::phaseModel +( + const word& phaseName, + const volScalarField& p, + const volScalarField& T +) +: + volScalarField + ( + IOobject + ( + IOobject::groupName("alpha", phaseName), + p.mesh().time().timeName(), + p.mesh(), + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + p.mesh() + ), + name_(phaseName), + p_(p), + T_(T), + thermo_(NULL), + dgdt_ + ( + IOobject + ( + IOobject::groupName("dgdt", phaseName), + p.mesh().time().timeName(), + p.mesh(), + IOobject::READ_IF_PRESENT, + IOobject::AUTO_WRITE + ), + p.mesh(), + dimensionedScalar("0", dimless/dimTime, 0) + ) +{ + { + volScalarField Tp(IOobject::groupName("T", phaseName), T); + Tp.write(); + } + + thermo_ = rhoThermo::New(p.mesh(), phaseName); + thermo_->validate(phaseName, "e"); + + correct(); +} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +Foam::autoPtr Foam::phaseModel::clone() const +{ + notImplemented("phaseModel::clone() const"); + return autoPtr(NULL); +} + + +void Foam::phaseModel::correct() +{ + thermo_->he() = thermo_->he(p_, T_); + thermo_->correct(); +} + + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/phaseModel/phaseModel.H b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/phaseModel/phaseModel.H new file mode 100644 index 0000000000..66d0ac8d63 --- /dev/null +++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/multiphaseMixtureThermo/phaseModel/phaseModel.H @@ -0,0 +1,156 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2013 OpenFOAM Foundation + \\/ 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 3 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, see . + +Class + Foam::phaseModel + +Description + Single incompressible phase derived from the phase-fraction. + Used as part of the multiPhaseMixture for interface-capturing multi-phase + simulations. + +SourceFiles + phaseModel.C + +\*---------------------------------------------------------------------------*/ + +#ifndef phaseModel_H +#define phaseModel_H + +#include "rhoThermo.H" +#include "volFields.H" +#include "dictionaryEntry.H" + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class phaseModel Declaration +\*---------------------------------------------------------------------------*/ + +class phaseModel +: + public volScalarField +{ + // Private data + + word name_; + const volScalarField& p_; + const volScalarField& T_; + autoPtr thermo_; + volScalarField dgdt_; + + +public: + + // Constructors + + //- Construct from components + phaseModel + ( + const word& phaseName, + const volScalarField& p, + const volScalarField& T + ); + + //- Return clone + autoPtr clone() const; + + //- Return a pointer to a new phaseModel created on freestore + // from Istream + class iNew + { + const volScalarField& p_; + const volScalarField& T_; + + public: + + iNew + ( + const volScalarField& p, + const volScalarField& T + ) + : + p_(p), + T_(T) + {} + + autoPtr operator()(Istream& is) const + { + return autoPtr(new phaseModel(is, p_, T_)); + } + }; + + + // Member Functions + + const word& name() const + { + return name_; + } + + const word& keyword() const + { + return name(); + } + + //- Return const-access to phase rhoThermo + const rhoThermo& thermo() const + { + return thermo_(); + } + + //- Return access to phase rhoThermo + rhoThermo& thermo() + { + return thermo_(); + } + + //- Return const-access to phase divergence + const volScalarField& dgdt() const + { + return dgdt_; + } + + //- Return access to phase divergence + volScalarField& dgdt() + { + return dgdt_; + } + + void correct(); +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/pEqn.H b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/pEqn.H new file mode 100644 index 0000000000..e42a54aabd --- /dev/null +++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/pEqn.H @@ -0,0 +1,142 @@ +{ + volScalarField rAU("rAU", 1.0/UEqn.A()); + surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU)); + + volVectorField HbyA("HbyA", U); + HbyA = rAU*UEqn.H(); + + surfaceScalarField phiHbyA + ( + "phiHbyA", + (fvc::interpolate(HbyA) & mesh.Sf()) + + fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, phi) + ); + + surfaceScalarField phig + ( + ( + multiphaseProperties.surfaceTensionForce() + - ghf*fvc::snGrad(rho) + )*rAUf*mesh.magSf() + ); + + phiHbyA += phig; + + // Update the fixedFluxPressure BCs to ensure flux consistency + setSnGrad + ( + p_rgh.boundaryField(), + ( + phiHbyA.boundaryField() + - (mesh.Sf().boundaryField() & U.boundaryField()) + )/(mesh.magSf().boundaryField()*rAUf.boundaryField()) + ); + + PtrList p_rghEqnComps(multiphaseProperties.phases().size()); + + label phasei = 0; + forAllConstIter + ( + PtrDictionary, + multiphaseProperties.phases(), + phase + ) + { + const rhoThermo& thermo = phase().thermo(); + const volScalarField& rho = thermo.rho()(); + + p_rghEqnComps.set + ( + phasei, + ( + fvc::ddt(rho) + thermo.psi()*correction(fvm::ddt(p_rgh)) + + fvc::div(phi, rho) - fvc::Sp(fvc::div(phi), rho) + ).ptr() + ); + + phasei++; + } + + // Cache p_rgh prior to solve for density update + volScalarField p_rgh_0(p_rgh); + + while (pimple.correctNonOrthogonal()) + { + fvScalarMatrix p_rghEqnIncomp + ( + fvc::div(phiHbyA) + - fvm::laplacian(rAUf, p_rgh) + ); + + tmp p_rghEqnComp; + + phasei = 0; + forAllConstIter + ( + PtrDictionary, + multiphaseProperties.phases(), + phase + ) + { + tmp hmm + ( + (max(phase(), scalar(0))/phase().thermo().rho()) + *p_rghEqnComps[phasei] + ); + + if (phasei == 0) + { + p_rghEqnComp = hmm; + } + else + { + p_rghEqnComp() += hmm; + } + + phasei++; + } + + solve + ( + p_rghEqnComp + + p_rghEqnIncomp, + mesh.solver(p_rgh.select(pimple.finalInnerIter())) + ); + + if (pimple.finalNonOrthogonalIter()) + { + // p = max(p_rgh + multiphaseProperties.rho()*gh, pMin); + // p_rgh = p - multiphaseProperties.rho()*gh; + + phasei = 0; + forAllIter + ( + PtrDictionary, + multiphaseProperties.phases(), + phase + ) + { + phase().dgdt() = + pos(phase()) + *(p_rghEqnComps[phasei] & p_rgh)/phase().thermo().rho(); + } + + phi = phiHbyA + p_rghEqnIncomp.flux(); + + U = HbyA + + rAU*fvc::reconstruct((phig + p_rghEqnIncomp.flux())/rAUf); + U.correctBoundaryConditions(); + } + } + + p = max(p_rgh + multiphaseProperties.rho()*gh, pMin); + + // Update densities from change in p_rgh + multiphaseProperties.correctRho(p_rgh - p_rgh_0); + rho = multiphaseProperties.rho(); + + K = 0.5*magSqr(U); + + Info<< "max(U) " << max(mag(U)).value() << endl; + Info<< "min(p_rgh) " << min(p_rgh).value() << endl; +} diff --git a/applications/solvers/multiphase/interFoam/LTSInterFoam/setrDeltaT.H b/applications/solvers/multiphase/interFoam/LTSInterFoam/setrDeltaT.H index f419153506..3f8ba982cf 100644 --- a/applications/solvers/multiphase/interFoam/LTSInterFoam/setrDeltaT.H +++ b/applications/solvers/multiphase/interFoam/LTSInterFoam/setrDeltaT.H @@ -121,9 +121,11 @@ && runTime.timeIndex() > runTime.startTimeIndex() + 1 ) { - rDeltaT = - rDeltaT0 - *max(rDeltaT/rDeltaT0, scalar(1.0) - rDeltaTDampingCoeff); + rDeltaT = max + ( + rDeltaT, + (scalar(1.0) - rDeltaTDampingCoeff)*rDeltaT0 + ); Info<< "Damped flow time scale min/max = " << gMin(1/rDeltaT.internalField()) diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C index f347de82a0..995d07d7e4 100644 --- a/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C +++ b/applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C @@ -126,11 +126,25 @@ void Foam::multiphaseSystem::solveAlphas() } // Ensure that the flux at inflow BCs is preserved - phiAlphaCorr.boundaryField() = min - ( - phase1.phi().boundaryField()*alpha1.boundaryField(), - phiAlphaCorr.boundaryField() - ); + forAll(phiAlphaCorr.boundaryField(), patchi) + { + fvsPatchScalarField& phiAlphaCorrp = + phiAlphaCorr.boundaryField()[patchi]; + + if (!phiAlphaCorrp.coupled()) + { + const scalarField& phi1p = phase1.phi().boundaryField()[patchi]; + const scalarField& alpha1p = alpha1.boundaryField()[patchi]; + + forAll(phiAlphaCorrp, facei) + { + if (phi1p[facei] < 0) + { + phiAlphaCorrp[facei] = alpha1p[facei]*phi1p[facei]; + } + } + } + } MULES::limit ( diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/alphaEqn.H b/applications/solvers/multiphase/twoPhaseEulerFoam/alphaEqn.H index 8ce175b85c..be420a36e1 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/alphaEqn.H +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/alphaEqn.H @@ -2,6 +2,8 @@ word alphaScheme("div(phi," + alpha1.name() + ')'); word alpharScheme("div(phir," + alpha1.name() + ')'); + alpha1.correctBoundaryConditions(); + surfaceScalarField phic("phic", phi); surfaceScalarField phir("phir", phi1 - phi2); @@ -93,11 +95,25 @@ ); // Ensure that the flux at inflow BCs is preserved - alphaPhic1.boundaryField() = min - ( - phi1.boundaryField()*alpha1.boundaryField(), - alphaPhic1.boundaryField() - ); + forAll(alphaPhic1.boundaryField(), patchi) + { + fvsPatchScalarField& alphaPhic1p = + alphaPhic1.boundaryField()[patchi]; + + if (!alphaPhic1p.coupled()) + { + const scalarField& phi1p = phi1.boundaryField()[patchi]; + const scalarField& alpha1p = alpha1.boundaryField()[patchi]; + + forAll(alphaPhic1p, facei) + { + if (phi1p[facei] < 0) + { + alphaPhic1p[facei] = alpha1p[facei]*phi1p[facei]; + } + } + } + } MULES::explicitSolve ( diff --git a/applications/test/PtrList/Test-PtrList.C b/applications/test/PtrList/Test-PtrList.C index 9241460a59..664dd5b501 100644 --- a/applications/test/PtrList/Test-PtrList.C +++ b/applications/test/PtrList/Test-PtrList.C @@ -32,6 +32,8 @@ Description #include "scalar.H" #include "IOstreams.H" #include "PtrList.H" +#include "plane.H" +#include "DynamicList.H" using namespace Foam; @@ -76,6 +78,7 @@ int main(int argc, char *argv[]) { PtrList list1(10); PtrList list2(15); + PtrList listApp; forAll(list1, i) { @@ -87,9 +90,14 @@ int main(int argc, char *argv[]) list2.set(i, new Scalar(10 + 1.3*i)); } + for (label i = 0; i < 5; ++i) + { + listApp.append(new Scalar(1.3*i)); + } Info<<"list1: " << list1 << endl; Info<<"list2: " << list2 << endl; + Info<<"listApp: " << listApp << endl; Info<<"indirectly delete some items via set(.., 0) :" << endl; for (label i = 0; i < 3; i++) @@ -115,6 +123,13 @@ int main(int argc, char *argv[]) Info<<"list2: " << list2 << endl; Info<<"list3: " << list3 << endl; + PtrList planes; + planes.append(new plane(vector::one, vector::one)); + planes.append(new plane(vector(1,2,3), vector::one)); + + forAll(planes, p) + Info<< "plane " << planes[p] << endl; + Info<< nl << "Done." << endl; return 0; } diff --git a/applications/test/UniformField/Make/files b/applications/test/UniformField/Make/files new file mode 100644 index 0000000000..3e268d0c3e --- /dev/null +++ b/applications/test/UniformField/Make/files @@ -0,0 +1,3 @@ +Test-UniformField.C + +EXE = $(FOAM_USER_APPBIN)/Test-UniformField diff --git a/applications/test/UniformField/Make/options b/applications/test/UniformField/Make/options new file mode 100644 index 0000000000..6a9e9810b3 --- /dev/null +++ b/applications/test/UniformField/Make/options @@ -0,0 +1,2 @@ +/* EXE_INC = -I$(LIB_SRC)/cfdTools/include */ +/* EXE_LIBS = -lfiniteVolume */ diff --git a/applications/test/UniformField/Test-UniformField.C b/applications/test/UniformField/Test-UniformField.C new file mode 100644 index 0000000000..88053afdba --- /dev/null +++ b/applications/test/UniformField/Test-UniformField.C @@ -0,0 +1,15 @@ +#include "UniformField.H" +#include "vector.H" +#include "IOstreams.H" + +using namespace Foam; + +int main() +{ + UniformField uf1(13.1); + UniformField uf2(vector(1, 2, 3)); + + Info<< "uf1 = " << uf1[22] << "; uf2 = " << uf2[1002] << endl; + + return 0; +} diff --git a/applications/test/fieldMapping/Make/files b/applications/test/fieldMapping/Make/files new file mode 100644 index 0000000000..6e0b3b5164 --- /dev/null +++ b/applications/test/fieldMapping/Make/files @@ -0,0 +1,7 @@ +/* +fvMeshGeometry.C +fvMesh.C +*/ +Test-fieldMapping.C + +EXE = $(FOAM_USER_APPBIN)/Test-fieldMapping diff --git a/applications/test/fieldMapping/Make/options b/applications/test/fieldMapping/Make/options new file mode 100644 index 0000000000..8c55d77648 --- /dev/null +++ b/applications/test/fieldMapping/Make/options @@ -0,0 +1,9 @@ +EXE_INC = \ + -DFULLDEBUG -g -O0 \ + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/dynamicMesh/lnInclude \ + -I$(LIB_SRC)/meshTools/lnInclude + +EXE_LIBS = \ + -lfiniteVolume \ + -ldynamicMesh diff --git a/applications/test/fieldMapping/README.txt b/applications/test/fieldMapping/README.txt new file mode 100644 index 0000000000..2c2ad0222e --- /dev/null +++ b/applications/test/fieldMapping/README.txt @@ -0,0 +1,8 @@ +Test application for volField and surfaceField mapping with topology +changes. + +Run + + pipe1D/Allrun + +to compile and map a few fields diff --git a/applications/test/fieldMapping/Test-fieldMapping.C b/applications/test/fieldMapping/Test-fieldMapping.C new file mode 100644 index 0000000000..eece86105e --- /dev/null +++ b/applications/test/fieldMapping/Test-fieldMapping.C @@ -0,0 +1,335 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2013 OpenFOAM Foundation + \\/ 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 3 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, see . + +Application + Test-fieldMapping + +Description + Test app for mapping of fields. + +\*---------------------------------------------------------------------------*/ + +#include "argList.H" +#include "fvMesh.H" +#include "volFields.H" +#include "Time.H" +#include "OFstream.H" +#include "meshTools.H" +#include "removeFaces.H" +#include "mapPolyMesh.H" +#include "polyTopoChange.H" +#include "fvcDiv.H" +#include "Random.H" + +using namespace Foam; + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +bool notEqual(const scalar s1, const scalar s2, const scalar tol) +{ + return mag(s1-s2) > tol; +} + + +// Main program: + +int main(int argc, char *argv[]) +{ +# include "addTimeOptions.H" + argList::validArgs.append("inflate (true|false)"); +# include "setRootCase.H" +# include "createTime.H" +# include "createMesh.H" + + const Switch inflate(args.args()[1]); + + if (inflate) + { + Info<< "Deleting cells using inflation/deflation" << nl << endl; + } + else + { + Info<< "Deleting cells, introducing points at new position" << nl + << endl; + } + + + Random rndGen(0); + + + + // Test mapping + // ------------ + // Mapping is volume averaged + + + // 1. uniform field stays uniform + volScalarField one + ( + IOobject + ( + "one", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + mesh, + dimensionedScalar("one", dimless, 1.0), + zeroGradientFvPatchScalarField::typeName + ); + Info<< "Writing one field " + << one.name() << " in " << runTime.timeName() << endl; + one.write(); + + + // 2. linear profile gets preserved + volScalarField ccX + ( + IOobject + ( + "ccX", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + mesh.C().component(0) + ); + Info<< "Writing x component of cell centres to " + << ccX.name() + << " in " << runTime.timeName() << endl; + ccX.write(); + + + // Uniform surface field + surfaceScalarField surfaceOne + ( + IOobject + ( + "surfaceOne", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + mesh, + dimensionedScalar("one", dimless, 1.0), + calculatedFvsPatchScalarField::typeName + ); + Info<< "Writing surface one field " + << surfaceOne.name() << " in " << runTime.timeName() << endl; + surfaceOne.write(); + + + + + // Force allocation of V. Important for any mesh changes since otherwise + // old time volumes are not stored + const scalar totalVol = gSum(mesh.V()); + + // Face removal engine. No checking for not merging boundary faces. + removeFaces faceRemover(mesh, GREAT); + + + while (runTime.loop()) + { + Info<< "Time = " << runTime.timeName() << nl << endl; + + if (!mesh.nInternalFaces()) + { + break; + } + + // Remove face + label candidateFaceI = rndGen.integer(0, mesh.nInternalFaces()-1); + Info<< "Wanting to delete face " << mesh.faceCentres()[candidateFaceI] + << nl << endl; + + labelList candidates(1, candidateFaceI); + + + // Get compatible set of faces and connected sets of cells. + labelList cellRegion; + labelList cellRegionMaster; + labelList facesToRemove; + + faceRemover.compatibleRemoves + ( + candidates, + cellRegion, + cellRegionMaster, + facesToRemove + ); + + // Topo changes container + polyTopoChange meshMod(mesh); + + // Insert mesh refinement into polyTopoChange. + faceRemover.setRefinement + ( + facesToRemove, + cellRegion, + cellRegionMaster, + meshMod + ); + + // Change mesh and inflate + Info<< "Actually changing mesh" << nl << endl; + autoPtr morphMap = meshMod.changeMesh(mesh, inflate); + + Info<< "Mapping fields" << nl << endl; + mesh.updateMesh(morphMap); + + // Move mesh (since morphing does not do this) + if (morphMap().hasMotionPoints()) + { + Info<< "Moving mesh" << nl << endl; + mesh.movePoints(morphMap().preMotionPoints()); + } + + // Update numbering of cells/vertices. + faceRemover.updateMesh(morphMap); + + + Info<< "Writing fields" << nl << endl; + runTime.write(); + + + // Check mesh volume conservation + if (mesh.moving()) + { + #include "volContinuity.H" + } + else + { + if (mesh.V().size() != mesh.nCells()) + { + FatalErrorIn(args.executable()) + << "Volume not mapped. V:" << mesh.V().size() + << " nCells:" << mesh.nCells() + << exit(FatalError); + } + + const scalar newVol = gSum(mesh.V()); + Info<< "Initial volume = " << totalVol + << " New volume = " << newVol + << endl; + + if (mag(newVol-totalVol)/totalVol > 1e-10) + { + FatalErrorIn(args.executable()) + << "Volume loss: old volume:" << totalVol + << " new volume:" << newVol + << exit(FatalError); + } + else + { + Info<< "Volume check OK" << nl << endl; + } + } + + + // Check constant profile + { + const scalar max = gMax(one); + const scalar min = gMin(one); + + Info<< "Uniform one field min = " << min + << " max = " << max << endl; + + if (notEqual(max, 1.0, 1e-10) || notEqual(min, 1.0, 1e-10)) + { + FatalErrorIn(args.executable()) + << "Uniform volVectorField not preserved." + << " Min and max should both be 1.0. min:" << min + << " max:" << max + << exit(FatalError); + } + else + { + Info<< "Uniform field mapping check OK" << nl << endl; + } + } + + // Check linear profile + { + const scalarField diff = ccX-mesh.C().component(0); + + const scalar max = gMax(diff); + const scalar min = gMin(diff); + + Info<< "Linear profile field min = " << min + << " max = " << max << endl; + + if (notEqual(max, 0.0, 1e-10) || notEqual(min, 0.0, 1e-10)) + { + FatalErrorIn(args.executable()) + << "Linear profile not preserved." + << " Min and max should both be 0.0. min:" << min + << " max:" << max + << exit(FatalError); + } + else + { + Info<< "Linear profile mapping check OK" << nl << endl; + } + } + + // Check face field mapping + if (surfaceOne.size()) + { + const scalar max = gMax(surfaceOne.internalField()); + const scalar min = gMin(surfaceOne.internalField()); + + Info<< "Uniform surface field min = " << min + << " max = " << max << endl; + + if (notEqual(max, 1.0, 1e-10) || notEqual(min, 1.0, 1e-10)) + { + FatalErrorIn(args.executable()) + << "Uniform surfaceScalarField not preserved." + << " Min and max should both be 1.0. min:" << min + << " max:" << max + << exit(FatalError); + } + else + { + Info<< "Uniform surfaceScalarField mapping check OK" << nl + << endl; + } + } + + + Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s" + << " ClockTime = " << runTime.elapsedClockTime() << " s" + << nl << endl; + } + + Info<< "End\n" << endl; + + return 0; +} + + +// ************************************************************************* // diff --git a/applications/test/fieldMapping/pipe1D/Allrun b/applications/test/fieldMapping/pipe1D/Allrun new file mode 100755 index 0000000000..575c24f73d --- /dev/null +++ b/applications/test/fieldMapping/pipe1D/Allrun @@ -0,0 +1,23 @@ +#!/bin/sh +cd ${0%/*} || exit 1 # run from this directory + +# Source tutorial run functions +. $WM_PROJECT_DIR/bin/tools/RunFunctions + +# Get application name +application=Test-fieldMapping + +# Compile +wmake .. + +runApplication blockMesh + +# Run with inflation +runApplication $application true +mv "log.$application" "log.$application-inflate" + +# Run without inflation +runApplication $application false + + +# ----------------------------------------------------------------- end-of-file diff --git a/applications/test/fieldMapping/pipe1D/constant/manualDecomposition b/applications/test/fieldMapping/pipe1D/constant/manualDecomposition new file mode 100644 index 0000000000..8f08c2826f --- /dev/null +++ b/applications/test/fieldMapping/pipe1D/constant/manualDecomposition @@ -0,0 +1,424 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: dev | +| \\ / A nd | Web: www.OpenFOAM.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class labelList; + location "constant"; + object cellDecomposition; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + + +400 +( +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +1 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +) + + +// ************************************************************************* // diff --git a/applications/test/fieldMapping/pipe1D/constant/polyMesh/blockMeshDict b/applications/test/fieldMapping/pipe1D/constant/polyMesh/blockMeshDict new file mode 100644 index 0000000000..14bfd1bc80 --- /dev/null +++ b/applications/test/fieldMapping/pipe1D/constant/polyMesh/blockMeshDict @@ -0,0 +1,61 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: dev | +| \\ / A nd | Web: www.OpenFOAM.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + object blockMeshDict; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +convertToMeters 1; + +vertices +( + (0 0 0) + (1 0 0) + (1 0.1 0) + (0 0.1 0) + (0 0 0.1) + (1 0 0.1) + (1 0.1 0.1) + (0 0.1 0.1) +); + +blocks +( + hex (0 1 2 3 4 5 6 7) (20 1 1) simpleGrading (10 1 1) +); + +edges +( +); + +boundary +( + inlet + { + type patch; + faces + ( + (0 4 7 3) + ); + } + + outlet + { + type patch; + faces + ( + (2 6 5 1) + ); + } +); + +// ************************************************************************* // diff --git a/applications/test/fieldMapping/pipe1D/constant/polyMesh/boundary b/applications/test/fieldMapping/pipe1D/constant/polyMesh/boundary new file mode 100644 index 0000000000..a284128db2 --- /dev/null +++ b/applications/test/fieldMapping/pipe1D/constant/polyMesh/boundary @@ -0,0 +1,41 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: dev | +| \\ / A nd | Web: www.OpenFOAM.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class polyBoundaryMesh; + location "constant/polyMesh"; + object boundary; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +3 +( + inlet + { + type patch; + nFaces 1; + startFace 19; + } + outlet + { + type patch; + nFaces 1; + startFace 20; + } + defaultFaces + { + type empty; + inGroups 1(empty); + nFaces 80; + startFace 21; + } +) + +// ************************************************************************* // diff --git a/applications/test/fieldMapping/pipe1D/constant/transportProperties b/applications/test/fieldMapping/pipe1D/constant/transportProperties new file mode 100644 index 0000000000..fa1c1ca0b1 --- /dev/null +++ b/applications/test/fieldMapping/pipe1D/constant/transportProperties @@ -0,0 +1,21 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: dev | +| \\ / A nd | Web: www.OpenFOAM.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + location "constant"; + object transportProperties; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +nu nu [ 0 2 -1 0 0 0 0 ] 0.01; + + +// ************************************************************************* // diff --git a/applications/test/fieldMapping/pipe1D/system/controlDict b/applications/test/fieldMapping/pipe1D/system/controlDict new file mode 100644 index 0000000000..ff7a82a791 --- /dev/null +++ b/applications/test/fieldMapping/pipe1D/system/controlDict @@ -0,0 +1,56 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: dev | +| \\ / A nd | Web: www.OpenFOAM.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + location "system"; + object controlDict; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +DebugSwitches +{ +// primitiveMesh 1; +// polyMesh 1; +// fvMesh 1; +} + +application Test-fieldMapping; + +startFrom startTime; + +startTime 0; + +stopAt endTime; + +endTime 100; + +deltaT 1; + +writeControl timeStep; + +writeInterval 1; + +purgeWrite 0; + +writeFormat ascii; + +writePrecision 10; + +writeCompression off; + +timeFormat general; + +timePrecision 6; + +runTimeModifiable true; + + +// ************************************************************************* // diff --git a/applications/test/fieldMapping/pipe1D/system/fvSchemes b/applications/test/fieldMapping/pipe1D/system/fvSchemes new file mode 100644 index 0000000000..dede0a6cba --- /dev/null +++ b/applications/test/fieldMapping/pipe1D/system/fvSchemes @@ -0,0 +1,57 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: dev | +| \\ / A nd | Web: www.OpenFOAM.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + location "system"; + object fvSchemes; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +ddtSchemes +{ + default Euler; +} + +gradSchemes +{ + default Gauss linear; + grad(p) Gauss linear; +} + +divSchemes +{ + default none; + div(phi,U) Gauss linear; +} + +laplacianSchemes +{ + default Gauss linear orthogonal; +} + +interpolationSchemes +{ + default linear; +} + +snGradSchemes +{ + default orthogonal; +} + +fluxRequired +{ + default no; + p ; +} + + +// ************************************************************************* // diff --git a/applications/test/fieldMapping/pipe1D/system/fvSolution b/applications/test/fieldMapping/pipe1D/system/fvSolution new file mode 100644 index 0000000000..cc4750f16c --- /dev/null +++ b/applications/test/fieldMapping/pipe1D/system/fvSolution @@ -0,0 +1,46 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: dev | +| \\ / A nd | Web: www.OpenFOAM.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + location "system"; + object fvSolution; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +solvers +{ + p + { + solver PCG; + preconditioner DIC; + tolerance 1e-06; + relTol 0; + } + + U + { + solver PBiCG; + preconditioner DILU; + tolerance 1e-05; + relTol 0; + } +} + +PISO +{ + nCorrectors 2; + nNonOrthogonalCorrectors 0; + pRefCell 0; + pRefValue 0; +} + + +// ************************************************************************* // diff --git a/applications/test/hexRef8/Make/files b/applications/test/hexRef8/Make/files new file mode 100644 index 0000000000..3106abb006 --- /dev/null +++ b/applications/test/hexRef8/Make/files @@ -0,0 +1,3 @@ +Test-hexRef8.C + +EXE = $(FOAM_USER_APPBIN)/Test-hexRef8 diff --git a/applications/test/hexRef8/Make/options b/applications/test/hexRef8/Make/options new file mode 100644 index 0000000000..8c55d77648 --- /dev/null +++ b/applications/test/hexRef8/Make/options @@ -0,0 +1,9 @@ +EXE_INC = \ + -DFULLDEBUG -g -O0 \ + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/dynamicMesh/lnInclude \ + -I$(LIB_SRC)/meshTools/lnInclude + +EXE_LIBS = \ + -lfiniteVolume \ + -ldynamicMesh diff --git a/applications/test/hexRef8/README.txt b/applications/test/hexRef8/README.txt new file mode 100644 index 0000000000..1ffaf7a7d5 --- /dev/null +++ b/applications/test/hexRef8/README.txt @@ -0,0 +1,11 @@ +Test application for volField and surfaceField mapping with +refinement/unrefinement + +Run + + block/Allrun + +to compile and map a few fields. Note that hexRef8 cannot be +run in inflation mode - there is the problem of getting +a set of faces to sweep that is consistent for a cell and +all its neighbours. diff --git a/applications/test/hexRef8/Test-hexRef8.C b/applications/test/hexRef8/Test-hexRef8.C new file mode 100644 index 0000000000..67477aae44 --- /dev/null +++ b/applications/test/hexRef8/Test-hexRef8.C @@ -0,0 +1,371 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2013 OpenFOAM Foundation + \\/ 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 3 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, see . + +Application + Test-hexRef8 + +Description + Test app for refinement and unrefinement. Runs a few iterations refining + and unrefining. + +\*---------------------------------------------------------------------------*/ + +#include "argList.H" +#include "Time.H" +#include "volFields.H" +#include "surfaceFields.H" +#include "hexRef8.H" +#include "mapPolyMesh.H" +#include "polyTopoChange.H" +#include "Random.H" +#include "zeroGradientFvPatchFields.H" +#include "fvcDiv.H" + +using namespace Foam; + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +bool notEqual(const scalar s1, const scalar s2, const scalar tol) +{ + return mag(s1-s2) > tol; +} + + +// Main program: +int main(int argc, char *argv[]) +{ +# include "addTimeOptions.H" + argList::validArgs.append("inflate (true|false)"); +# include "setRootCase.H" +# include "createTime.H" +# include "createMesh.H" + + + const Switch inflate(args.args()[1]); + + if (inflate) + { + Info<< "Splitting/deleting cells using inflation/deflation" << nl + << endl; + } + else + { + Info<< "Splitting/deleting cells, introducing points at new position" + << nl << endl; + } + + + Random rndGen(0); + + + // Force generation of V() + (void)mesh.V(); + + + // Test mapping + // ------------ + + // 1. uniform field stays uniform + volScalarField one + ( + IOobject + ( + "one", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + mesh, + dimensionedScalar("one", dimless, 1.0), + zeroGradientFvPatchScalarField::typeName + ); + Info<< "Writing one field " + << one.name() << " in " << runTime.timeName() << endl; + one.write(); + + + // 2. linear profile gets preserved + volScalarField ccX + ( + IOobject + ( + "ccX", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + mesh.C().component(0) + ); + Info<< "Writing x component of cell centres to " + << ccX.name() + << " in " << runTime.timeName() << endl; + ccX.write(); + + + // Uniform surface field + surfaceScalarField surfaceOne + ( + IOobject + ( + "surfaceOne", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + mesh, + dimensionedScalar("one", dimless, 1.0), + calculatedFvsPatchScalarField::typeName + ); + Info<< "Writing surface one field " + << surfaceOne.name() << " in " << runTime.timeName() << endl; + surfaceOne.write(); + + + // Force allocation of V. Important for any mesh changes since otherwise + // old time volumes are not stored + const scalar totalVol = gSum(mesh.V()); + + + // Construct refiner. Read initial cell and point levels. + hexRef8 meshCutter(mesh); + + + while (runTime.loop()) + { + Info<< "Time = " << runTime.timeName() << nl << endl; + + if (mesh.globalData().nTotalCells() == 0) + { + break; + } + + + + // Mesh changing engine. + polyTopoChange meshMod(mesh); + + if (rndGen.bit()) + { + // Refine + label nRefine = mesh.nCells()/20; + DynamicList