diff --git a/.gitmodules b/.gitmodules new file mode 100644 index 0000000000..ee2b053a21 --- /dev/null +++ b/.gitmodules @@ -0,0 +1,3 @@ +[submodule "cfmesh"] + path = modules/cfmesh + url = https://develop.openfoam.com/Community/integration-cfmesh.git diff --git a/Allwmake b/Allwmake index b9cfa95a5e..ed3d03e789 100755 --- a/Allwmake +++ b/Allwmake @@ -43,6 +43,15 @@ echo "Compile OpenFOAM applications" echo applications/Allwmake $targetType $* +# Additional components/modules +if [ -d "$WM_PROJECT_DIR/modules" ] +then + echo "========================================" + echo "Compile OpenFOAM modules" + echo + (cd $WM_PROJECT_DIR/modules 2>/dev/null && wmake -all) +fi + # Some summary information echo date "+%Y-%m-%d %H:%M:%S %z" 2>/dev/null || echo "date is unknown" diff --git a/applications/solvers/basic/potentialFoam/overPotentialFoam/Make/files b/applications/solvers/basic/potentialFoam/overPotentialFoam/Make/files new file mode 100644 index 0000000000..511066afa5 --- /dev/null +++ b/applications/solvers/basic/potentialFoam/overPotentialFoam/Make/files @@ -0,0 +1,3 @@ +potentialFoam.C + +EXE = $(FOAM_APPBIN)/overPotentialFoam diff --git a/applications/solvers/basic/potentialFoam/overPotentialFoam/Make/options b/applications/solvers/basic/potentialFoam/overPotentialFoam/Make/options new file mode 100644 index 0000000000..155894029a --- /dev/null +++ b/applications/solvers/basic/potentialFoam/overPotentialFoam/Make/options @@ -0,0 +1,12 @@ +EXE_INC = \ + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/dynamicFvMesh/lnInclude \ + -I$(LIB_SRC)/overset/lnInclude \ + -I$(LIB_SRC)/meshTools/lnInclude \ + -I$(LIB_SRC)/sampling/lnInclude + +EXE_LIBS = \ + -lfiniteVolume \ + -lmeshTools \ + -lsampling \ + -loverset diff --git a/applications/solvers/basic/potentialFoam/overPotentialFoam/createControls.H b/applications/solvers/basic/potentialFoam/overPotentialFoam/createControls.H new file mode 100644 index 0000000000..7015cae02e --- /dev/null +++ b/applications/solvers/basic/potentialFoam/overPotentialFoam/createControls.H @@ -0,0 +1,9 @@ +const dictionary& potentialFlow +( + mesh.solutionDict().subDict("potentialFlow") +); + +const int nNonOrthCorr +( + potentialFlow.lookupOrDefault("nNonOrthogonalCorrectors", 0) +); diff --git a/applications/solvers/basic/potentialFoam/overPotentialFoam/createFields.H b/applications/solvers/basic/potentialFoam/overPotentialFoam/createFields.H new file mode 100644 index 0000000000..b7daa4fe47 --- /dev/null +++ b/applications/solvers/basic/potentialFoam/overPotentialFoam/createFields.H @@ -0,0 +1,146 @@ +Info<< "Reading velocity field U\n" << endl; +volVectorField U +( + IOobject + ( + "U", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh +); + +// Initialise the velocity internal field to zero +U = dimensionedVector("0", U.dimensions(), Zero); + +surfaceScalarField phi +( + IOobject + ( + "phi", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + fvc::flux(U) +); + +if (args.optionFound("initialiseUBCs")) +{ + U.correctBoundaryConditions(); + phi = fvc::flux(U); +} + + +// Construct a pressure field +// If it is available read it otherwise construct from the velocity BCs +// converting fixed-value BCs to zero-gradient and vice versa. +word pName("p"); + +// Update name of the pressure field from the command-line option +args.optionReadIfPresent("pName", pName); + +// Infer the pressure BCs from the velocity +wordList pBCTypes +( + U.boundaryField().size(), + fixedValueFvPatchScalarField::typeName +); + +forAll(U.boundaryField(), patchi) +{ + if (U.boundaryField()[patchi].fixesValue()) + { + pBCTypes[patchi] = zeroGradientFvPatchScalarField::typeName; + } +} + +// Note that registerObject is false for the pressure field. The pressure +// field in this solver doesn't have a physical value during the solution. +// It shouldn't be looked up and used by sub models or boundary conditions. +Info<< "Constructing pressure field " << pName << nl << endl; +volScalarField p +( + IOobject + ( + pName, + runTime.timeName(), + mesh, + IOobject::READ_IF_PRESENT, + IOobject::NO_WRITE, + false + ), + mesh, + dimensionedScalar(pName, sqr(dimVelocity), 0), + pBCTypes +); + +// Infer the velocity potential BCs from the pressure +wordList PhiBCTypes +( + p.boundaryField().size(), + zeroGradientFvPatchScalarField::typeName +); + +forAll(p.boundaryField(), patchi) +{ + if (p.boundaryField()[patchi].fixesValue()) + { + PhiBCTypes[patchi] = fixedValueFvPatchScalarField::typeName; + } +} + +Info<< "Constructing velocity potential field Phi\n" << endl; +volScalarField Phi +( + IOobject + ( + "Phi", + runTime.timeName(), + mesh, + IOobject::READ_IF_PRESENT, + IOobject::NO_WRITE + ), + mesh, + dimensionedScalar("Phi", dimLength*dimVelocity, 0), + PhiBCTypes +); + +label PhiRefCell = 0; +scalar PhiRefValue = 0; +setRefCell +( + Phi, + potentialFlow.dict(), + PhiRefCell, + PhiRefValue +); +mesh.setFluxRequired(Phi.name()); + +#include "createMRF.H" + +// Add overset specific interpolations +{ + dictionary oversetDict; + oversetDict.add("Phi", true); + oversetDict.add("U", true); + + const_cast + ( + mesh.schemesDict() + ).add + ( + "oversetInterpolationRequired", + oversetDict, + true + ); +} + +// Mask field for zeroing out contributions on hole cells +#include "createCellMask.H" + +// Create bool field with interpolated cells +#include "createInterpolatedCells.H" diff --git a/applications/solvers/basic/potentialFoam/overPotentialFoam/potentialFoam.C b/applications/solvers/basic/potentialFoam/overPotentialFoam/potentialFoam.C new file mode 100644 index 0000000000..7d340fbc36 --- /dev/null +++ b/applications/solvers/basic/potentialFoam/overPotentialFoam/potentialFoam.C @@ -0,0 +1,260 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2017 OpenCFD Ltd + \\/ 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 + potentialFoam + +Group + grpBasicSolvers + +Description + Potential flow solver which solves for the velocity potential, to + calculate the flux-field, from which the velocity field is obtained by + reconstructing the flux. + + \heading Solver details + The potential flow solution is typically employed to generate initial fields + for full Navier-Stokes codes. The flow is evolved using the equation: + + \f[ + \laplacian \Phi = \div(\vec{U}) + \f] + + Where: + \vartable + \Phi | Velocity potential [m2/s] + \vec{U} | Velocity [m/s] + \endvartable + + The corresponding pressure field could be calculated from the divergence + of the Euler equation: + + \f[ + \laplacian p + \div(\div(\vec{U}\otimes\vec{U})) = 0 + \f] + + but this generates excessive pressure variation in regions of large + velocity gradient normal to the flow direction. A better option is to + calculate the pressure field corresponding to velocity variation along the + stream-lines: + + \f[ + \laplacian p + \div(\vec{F}\cdot\div(\vec{U}\otimes\vec{U})) = 0 + \f] + where the flow direction tensor \f$\vec{F}\f$ is obtained from + \f[ + \vec{F} = \hat{\vec{U}}\otimes\hat{\vec{U}} + \f] + + \heading Required fields + \plaintable + U | Velocity [m/s] + \endplaintable + + \heading Optional fields + \plaintable + p | Kinematic pressure [m2/s2] + Phi | Velocity potential [m2/s] + | Generated from p (if present) or U if not present + \endplaintable + + \heading Options + \plaintable + -writep | write the Euler pressure + -writePhi | Write the final velocity potential + -initialiseUBCs | Update the velocity boundaries before solving for Phi + \endplaintable + +\*---------------------------------------------------------------------------*/ + +#include "fvCFD.H" +#include "pisoControl.H" +#include "dynamicFvMesh.H" +#include "cellCellStencilObject.H" +#include "localMin.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +int main(int argc, char *argv[]) +{ + argList::addOption + ( + "pName", + "pName", + "Name of the pressure field" + ); + + argList::addBoolOption + ( + "initialiseUBCs", + "Initialise U boundary conditions" + ); + + argList::addBoolOption + ( + "writePhi", + "Write the final velocity potential field" + ); + + argList::addBoolOption + ( + "writep", + "Calculate and write the Euler pressure field" + ); + + argList::addBoolOption + ( + "withFunctionObjects", + "execute functionObjects" + ); + + #include "setRootCase.H" + #include "createTime.H" + #include "createNamedDynamicFvMesh.H" + + pisoControl potentialFlow(mesh, "potentialFlow"); + + #include "createFields.H" + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + + Info<< nl << "Calculating potential flow" << endl; + + mesh.update(); + + surfaceScalarField faceMask(localMin(mesh).interpolate(cellMask)); + + // Since solver contains no time loop it would never execute + // function objects so do it ourselves + runTime.functionObjects().start(); + + MRF.makeRelative(phi); + adjustPhi(phi, U, p); + + // Non-orthogonal velocity potential corrector loop + while (potentialFlow.correct()) + { + phi = fvc::flux(U); + + while (potentialFlow.correctNonOrthogonal()) + { + fvScalarMatrix PhiEqn + ( + fvm::laplacian(faceMask, Phi) + == + fvc::div(phi) + ); + + PhiEqn.setReference(PhiRefCell, PhiRefValue); + PhiEqn.solve(); + + if (potentialFlow.finalNonOrthogonalIter()) + { + phi -= PhiEqn.flux(); + } + } + + MRF.makeAbsolute(phi); + + Info<< "Continuity error = " + << mag(fvc::div(phi))().weightedAverage(mesh.V()).value() + << endl; + + U = fvc::reconstruct(phi); + U.correctBoundaryConditions(); + + Info<< "Interpolated velocity error = " + << (sqrt(sum(sqr(fvc::flux(U) - phi)))/sum(mesh.magSf())).value() + << endl; + } + + // Write U and phi + U.write(); + phi.write(); + + // Optionally write Phi + if (args.optionFound("writePhi")) + { + Phi.write(); + } + + // Calculate the pressure field from the Euler equation + if (args.optionFound("writep")) + { + Info<< nl << "Calculating approximate pressure field" << endl; + + label pRefCell = 0; + scalar pRefValue = 0.0; + setRefCell + ( + p, + potentialFlow.dict(), + pRefCell, + pRefValue + ); + + // Calculate the flow-direction filter tensor + volScalarField magSqrU(magSqr(U)); + volSymmTensorField F(sqr(U)/(magSqrU + SMALL*average(magSqrU))); + + // Calculate the divergence of the flow-direction filtered div(U*U) + // Filtering with the flow-direction generates a more reasonable + // pressure distribution in regions of high velocity gradient in the + // direction of the flow + volScalarField divDivUU + ( + fvc::div + ( + F & fvc::div(phi, U), + "div(div(phi,U))" + ) + ); + + // Solve a Poisson equation for the approximate pressure + while (potentialFlow.correctNonOrthogonal()) + { + fvScalarMatrix pEqn + ( + fvm::laplacian(p) + divDivUU + ); + + pEqn.setReference(pRefCell, pRefValue); + pEqn.solve(); + } + + p.write(); + } + + runTime.functionObjects().end(); + + Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s" + << " ClockTime = " << runTime.elapsedClockTime() << " s" + << nl << endl; + + Info<< "End\n" << endl; + + return 0; +} + + +// ************************************************************************* // diff --git a/applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/pEqn.H b/applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/pEqn.H index 5f1fc7451c..ba9cbf752c 100644 --- a/applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/pEqn.H +++ b/applications/solvers/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam/pEqn.H @@ -7,8 +7,8 @@ surfaceScalarField faceMask(localMin(mesh).interpolate(cellMask)); volScalarField rAU(1.0/UEqn.A()); surfaceScalarField rhorAUf("rhorAUf", faceMask*fvc::interpolate(rho*rAU)); -volVectorField HbyA("HbyA", constrainHbyA(rAU*UEqn.H(), U, p)); -//mesh.interpolate(HbyA); +volVectorField HbyA("HbyA", U); +HbyA = constrainHbyA(rAU*UEqn.H(), U, p); if (pimple.nCorrPISO() <= 1) { diff --git a/applications/solvers/compressible/rhoSimpleFoam/overRhoSimpleFoam/Make/files b/applications/solvers/compressible/rhoSimpleFoam/overRhoSimpleFoam/Make/files new file mode 100644 index 0000000000..0519e8b961 --- /dev/null +++ b/applications/solvers/compressible/rhoSimpleFoam/overRhoSimpleFoam/Make/files @@ -0,0 +1,3 @@ +rhoSimpleFoam.C + +EXE = $(FOAM_APPBIN)/overRhoSimpleFoam diff --git a/applications/solvers/compressible/rhoSimpleFoam/overRhoSimpleFoam/Make/options b/applications/solvers/compressible/rhoSimpleFoam/overRhoSimpleFoam/Make/options new file mode 100644 index 0000000000..e027221b0c --- /dev/null +++ b/applications/solvers/compressible/rhoSimpleFoam/overRhoSimpleFoam/Make/options @@ -0,0 +1,27 @@ +EXE_INC = \ + -I.. \ + -I$(LIB_SRC)/transportModels/compressible/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \ + -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \ + -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \ + -I$(LIB_SRC)/finiteVolume/cfdTools \ + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/sampling/lnInclude \ + -I$(LIB_SRC)/dynamicFvMesh/lnInclude \ + -I$(LIB_SRC)/dynamicMesh/lnInclude \ + -I$(LIB_SRC)/meshTools/lnInclude \ + -I$(LIB_SRC)/overset/lnInclude + +EXE_LIBS = \ + -lcompressibleTransportModels \ + -lfluidThermophysicalModels \ + -lspecie \ + -lturbulenceModels \ + -lcompressibleTurbulenceModels \ + -lfiniteVolume \ + -lsampling \ + -lmeshTools \ + -lfvOptions \ + -loverset \ + -ldynamicFvMesh \ + -ltopoChangerFvMesh diff --git a/applications/solvers/compressible/rhoSimpleFoam/overRhoSimpleFoam/UEqn.H b/applications/solvers/compressible/rhoSimpleFoam/overRhoSimpleFoam/UEqn.H new file mode 100644 index 0000000000..9d4afec58e --- /dev/null +++ b/applications/solvers/compressible/rhoSimpleFoam/overRhoSimpleFoam/UEqn.H @@ -0,0 +1,23 @@ + // Solve the Momentum equation + MRF.correctBoundaryVelocity(U); + + tmp tUEqn + ( + fvm::div(phi, U) + + MRF.DDt(rho, U) + + turbulence->divDevRhoReff(U) + == + fvOptions(rho, U) + ); + fvVectorMatrix& UEqn = tUEqn.ref(); + + UEqn.relax(); + + fvOptions.constrain(UEqn); + + if (simple.momentumPredictor()) + { + solve(UEqn == -cellMask*fvc::grad(p)); + } + + fvOptions.correct(U); diff --git a/applications/solvers/compressible/rhoSimpleFoam/overRhoSimpleFoam/createFieldRefs.H b/applications/solvers/compressible/rhoSimpleFoam/overRhoSimpleFoam/createFieldRefs.H new file mode 100644 index 0000000000..502b3b4230 --- /dev/null +++ b/applications/solvers/compressible/rhoSimpleFoam/overRhoSimpleFoam/createFieldRefs.H @@ -0,0 +1 @@ +const volScalarField& psi = thermo.psi(); diff --git a/applications/solvers/compressible/rhoSimpleFoam/overRhoSimpleFoam/createFields.H b/applications/solvers/compressible/rhoSimpleFoam/overRhoSimpleFoam/createFields.H new file mode 100644 index 0000000000..2544f3e559 --- /dev/null +++ b/applications/solvers/compressible/rhoSimpleFoam/overRhoSimpleFoam/createFields.H @@ -0,0 +1,91 @@ +Info<< "Reading thermophysical properties\n" << endl; + +autoPtr pThermo +( + fluidThermo::New(mesh) +); +fluidThermo& thermo = pThermo(); +thermo.validate(args.executable(), "h", "e"); + +volScalarField& p = thermo.p(); + +volScalarField rho +( + IOobject + ( + "rho", + runTime.timeName(), + mesh, + IOobject::READ_IF_PRESENT, + IOobject::AUTO_WRITE + ), + thermo.rho() +); + +Info<< "Reading field U\n" << endl; +volVectorField U +( + IOobject + ( + "U", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh +); + +#include "compressibleCreatePhi.H" + +pressureControl pressureControl(p, rho, simple.dict()); + +mesh.setFluxRequired(p.name()); + +Info<< "Creating turbulence model\n" << endl; +autoPtr turbulence +( + compressible::turbulenceModel::New + ( + rho, + U, + phi, + thermo + ) +); + +dimensionedScalar initialMass = fvc::domainIntegrate(rho); + +#include "createMRF.H" + +//- Overset specific + +// Add solver-specific interpolations +{ + dictionary oversetDict; + oversetDict.add("U", true); + oversetDict.add("p", true); + oversetDict.add("HbyA", true); + oversetDict.add("grad(p)", true); + oversetDict.add("rho", true); + + const_cast + ( + mesh.schemesDict() + ).add + ( + "oversetInterpolationRequired", + oversetDict, + true + ); +} + +// Mask field for zeroing out contributions on hole cells +#include "createCellMask.H" + +#include "createInterpolatedCells.H" + +bool adjustFringe +( + simple.dict().lookupOrDefault("oversetAdjustPhi", false) +); diff --git a/applications/solvers/compressible/rhoSimpleFoam/overRhoSimpleFoam/createUpdatedDynamicFvMesh.H b/applications/solvers/compressible/rhoSimpleFoam/overRhoSimpleFoam/createUpdatedDynamicFvMesh.H new file mode 100644 index 0000000000..db58f90828 --- /dev/null +++ b/applications/solvers/compressible/rhoSimpleFoam/overRhoSimpleFoam/createUpdatedDynamicFvMesh.H @@ -0,0 +1,22 @@ + Info<< "Create dynamic mesh for time = " + << runTime.timeName() << nl << endl; + + autoPtr meshPtr + ( + dynamicFvMesh::New + ( + IOobject + ( + dynamicFvMesh::defaultRegion, + runTime.timeName(), + runTime, + IOobject::MUST_READ + ) + ) + ); + + dynamicFvMesh& mesh = meshPtr(); + + // Calculate initial mesh-to-mesh mapping. Note that this should be + // done under the hood, e.g. as a MeshObject + mesh.update(); diff --git a/applications/solvers/compressible/rhoSimpleFoam/overRhoSimpleFoam/pEqn.H b/applications/solvers/compressible/rhoSimpleFoam/overRhoSimpleFoam/pEqn.H new file mode 100644 index 0000000000..b4e670bee7 --- /dev/null +++ b/applications/solvers/compressible/rhoSimpleFoam/overRhoSimpleFoam/pEqn.H @@ -0,0 +1,123 @@ +{ + surfaceScalarField faceMask(localMin(mesh).interpolate(cellMask)); + + volScalarField rAU(1.0/UEqn.A()); + surfaceScalarField rhorAUf("rhorAUf", faceMask*fvc::interpolate(rho*rAU)); + volVectorField HbyA("HbyA", U); + HbyA = constrainHbyA(cellMask*rAU*UEqn.H(), U, p); + tUEqn.clear(); + + bool closedVolume = false; + + surfaceScalarField phiHbyA("phiHbyA", fvc::interpolate(rho)*fvc::flux(HbyA)); + MRF.makeRelative(fvc::interpolate(rho), phiHbyA); + + // Update the pressure BCs to ensure flux consistency + constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF); + + if (simple.transonic()) + { + surfaceScalarField phid + ( + "phid", + (fvc::interpolate(psi)/fvc::interpolate(rho))*phiHbyA + ); + + phiHbyA -= fvc::interpolate(psi*p)*phiHbyA/fvc::interpolate(rho); + + while (simple.correctNonOrthogonal()) + { + fvScalarMatrix pEqn + ( + fvc::div(phiHbyA) + + fvm::div(phid, p) + - fvm::laplacian(rhorAUf, p) + == + fvOptions(psi, p, rho.name()) + ); + + // Relax the pressure equation to ensure diagonal-dominance + pEqn.relax(); + + pEqn.setReference + ( + pressureControl.refCell(), + pressureControl.refValue() + ); + + pEqn.solve(); + + if (simple.finalNonOrthogonalIter()) + { + phi = phiHbyA + pEqn.flux(); + } + } + } + else + { + closedVolume = adjustPhi(phiHbyA, U, p); + + if (adjustFringe) + { + oversetAdjustPhi(phiHbyA, U); + } + + while (simple.correctNonOrthogonal()) + { + fvScalarMatrix pEqn + ( + fvc::div(phiHbyA) + - fvm::laplacian(rhorAUf, p) + == + fvOptions(psi, p, rho.name()) + ); + + pEqn.setReference + ( + pressureControl.refCell(), + pressureControl.refValue() + ); + + pEqn.solve(); + + if (simple.finalNonOrthogonalIter()) + { + phi = phiHbyA + pEqn.flux(); + } + } + } + + #include "incompressible/continuityErrs.H" + + // Explicitly relax pressure for momentum corrector + p.relax(); + + volVectorField gradP(fvc::grad(p)); + + U = HbyA - rAU*cellMask*gradP; + U.correctBoundaryConditions(); + fvOptions.correct(U); + + + bool pLimited = pressureControl.limit(p); + + // For closed-volume cases adjust the pressure and density levels + // to obey overall mass continuity + if (closedVolume) + { + p += (initialMass - fvc::domainIntegrate(psi*p)) + /fvc::domainIntegrate(psi); + } + + if (pLimited || closedVolume) + { + p.correctBoundaryConditions(); + } + + rho = thermo.rho(); + + if (!simple.transonic()) + { + rho.relax(); + } +} diff --git a/applications/solvers/compressible/rhoSimpleFoam/overRhoSimpleFoam/rhoSimpleFoam.C b/applications/solvers/compressible/rhoSimpleFoam/overRhoSimpleFoam/rhoSimpleFoam.C new file mode 100644 index 0000000000..d2aced6375 --- /dev/null +++ b/applications/solvers/compressible/rhoSimpleFoam/overRhoSimpleFoam/rhoSimpleFoam.C @@ -0,0 +1,92 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2017 OpenCFD Ltd + \\/ 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 + overRhoSimpleFoam + +Group + grpCompressibleSolvers + +Description + Overset steady-state solver for turbulent flow of compressible fluids. + +\*---------------------------------------------------------------------------*/ + +#include "fvCFD.H" +#include "dynamicFvMesh.H" +#include "fluidThermo.H" +#include "turbulentFluidThermoModel.H" +#include "simpleControl.H" +#include "pressureControl.H" +#include "fvOptions.H" +#include "cellCellStencilObject.H" +#include "localMin.H" +#include "oversetAdjustPhi.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +int main(int argc, char *argv[]) +{ + #define CREATE_MESH createUpdatedDynamicFvMesh.H + #include "postProcess.H" + + #include "setRootCase.H" + #include "createTime.H" + #include "createUpdatedDynamicFvMesh.H" + #include "createControl.H" + #include "createFields.H" + #include "createFieldRefs.H" + #include "createFvOptions.H" + #include "initContinuityErrs.H" + + turbulence->validate(); + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + + Info<< "\nStarting time loop\n" << endl; + + while (simple.loop()) + { + Info<< "Time = " << runTime.timeName() << nl << endl; + + // Pressure-velocity SIMPLE corrector + #include "UEqn.H" + #include "EEqn.H" + #include "pEqn.H" + + 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/doc/solver.dox b/applications/solvers/doc/solver.dox index 4c829de162..56bc8f1fd7 100644 --- a/applications/solvers/doc/solver.dox +++ b/applications/solvers/doc/solver.dox @@ -3,7 +3,7 @@ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | \\ / A nd | Copyright (C) 2015 OpenFOAM Foundation - \\/ M anipulation | + \\/ M anipulation | Copyright (C) 2017 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -39,6 +39,7 @@ The available solvers are grouped into the following categories: - \ref grpLagrangianSolvers - \ref grpMultiphaseSolvers - \ref grpStressAnalysisSolvers + - \ref grpFiniteAreaSolvers \*---------------------------------------------------------------------------*/ diff --git a/applications/solvers/doc/solversDoc.H b/applications/solvers/doc/solversDoc.H index da58bc0ce3..c0a01e5d76 100644 --- a/applications/solvers/doc/solversDoc.H +++ b/applications/solvers/doc/solversDoc.H @@ -3,7 +3,7 @@ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | \\ / A nd | Copyright (C) 2015 OpenFOAM Foundation - \\/ M anipulation | + \\/ M anipulation | Copyright (C) 2017 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -34,4 +34,10 @@ License This group contains moving mesh solvers solvers @} +\defgroup grpFiniteAreaSolvers Finite area solvers +@{ + \ingroup grpSolvers + This group contains finite area solvers +@} + \*---------------------------------------------------------------------------*/ diff --git a/applications/solvers/finiteArea/liquidFilmFoam/Make/files b/applications/solvers/finiteArea/liquidFilmFoam/Make/files new file mode 100755 index 0000000000..39ee7f75e9 --- /dev/null +++ b/applications/solvers/finiteArea/liquidFilmFoam/Make/files @@ -0,0 +1,3 @@ +liquidFilmFoam.C + +EXE = $(FOAM_APPBIN)/liquidFilmFoam diff --git a/applications/solvers/finiteArea/liquidFilmFoam/Make/options b/applications/solvers/finiteArea/liquidFilmFoam/Make/options new file mode 100755 index 0000000000..02549914ff --- /dev/null +++ b/applications/solvers/finiteArea/liquidFilmFoam/Make/options @@ -0,0 +1,10 @@ +EXE_INC = \ + -I$(LIB_SRC)/finiteArea/lnInclude \ + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/meshTools/lnInclude \ + -I$(LIB_SRC)/cfdTools/general/lnInclude + +EXE_LIBS = \ + -lfiniteArea \ + -lfiniteVolume \ + -lmeshTools diff --git a/applications/solvers/finiteArea/liquidFilmFoam/calcFrictionFactor.H b/applications/solvers/finiteArea/liquidFilmFoam/calcFrictionFactor.H new file mode 100644 index 0000000000..701b4635a1 --- /dev/null +++ b/applications/solvers/finiteArea/liquidFilmFoam/calcFrictionFactor.H @@ -0,0 +1,7 @@ +{ + // Stabilisation of friction factor calculation + // Friction factor is defined with standard gravity + frictionFactor.primitiveFieldRef() = + mag(2*9.81*sqr(manningField.primitiveField())/ + pow(mag(h.primitiveField()) + 1e-7, 1.0/3.0)); +} diff --git a/applications/solvers/finiteArea/liquidFilmFoam/capillaryCourantNo.H b/applications/solvers/finiteArea/liquidFilmFoam/capillaryCourantNo.H new file mode 100644 index 0000000000..8e66fd171f --- /dev/null +++ b/applications/solvers/finiteArea/liquidFilmFoam/capillaryCourantNo.H @@ -0,0 +1,13 @@ +{ + scalar CoNumSigma = max + ( + sqrt + ( + 2*M_PI*sigma*sqr(aMesh.edgeInterpolation::deltaCoeffs()) + *aMesh.edgeInterpolation::deltaCoeffs() + /rhol + ) + ).value()*runTime.deltaT().value(); + + Info<< "Max Capillary Courant Number = " << CoNumSigma << '\n' << endl; +} diff --git a/applications/solvers/finiteArea/liquidFilmFoam/createFaFields.H b/applications/solvers/finiteArea/liquidFilmFoam/createFaFields.H new file mode 100644 index 0000000000..5c760788c9 --- /dev/null +++ b/applications/solvers/finiteArea/liquidFilmFoam/createFaFields.H @@ -0,0 +1,158 @@ + Info<< "Reading field h" << endl; + areaScalarField h + ( + IOobject + ( + "h", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + aMesh + ); + + + Info<< "Reading field Us" << endl; + areaVectorField Us + ( + IOobject + ( + "Us", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + aMesh + ); + + + edgeScalarField phis + ( + IOobject + ( + "phis", + runTime.timeName(), + mesh, + IOobject::READ_IF_PRESENT, + IOobject::AUTO_WRITE + ), + fac::interpolate(Us) & aMesh.Le() + ); + + + edgeScalarField phi2s + ( + IOobject + ( + "phi2s", + runTime.timeName(), + mesh, + IOobject::READ_IF_PRESENT, + IOobject::AUTO_WRITE + ), + fac::interpolate(h*Us) & aMesh.Le() + ); + + + const areaVectorField& Ns = aMesh.faceAreaNormals(); + areaVectorField Gs(g - Ns*(Ns & g)); + areaScalarField Gn(mag(g - Gs)); + + // Mass source + areaScalarField Sm + ( + IOobject + ( + "Sm", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + aMesh, + dimensionedScalar("Sm", dimLength/dimTime, 0) + ); + + // Mass sink + areaScalarField Sd + ( + IOobject + ( + "Sd", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + aMesh, + dimensionedScalar("Sd", dimLength/dimTime, 0) + ); + + areaVectorField Ug + ( + IOobject + ( + "Sg", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + aMesh, + dimensionedVector("Ug", dimVelocity, vector::zero) + ); + + + // Surface pressure + areaScalarField ps + ( + IOobject + ( + "ps", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + rhol*Gn*h - sigma*fac::laplacian(h) + ); + + // Friction factor + areaScalarField dotProduct + ( + aMesh.faceAreaNormals() & (g/mag(g)) + ); + + Info<< "View factor: min = " << min(dotProduct.internalField()) + << " max = " << max(dotProduct.internalField()) << endl; + + areaScalarField manningField + ( + IOobject + ( + "manningField", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + aMesh + ); + + areaScalarField frictionFactor + ( + IOobject + ( + "frictionFactor", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + aMesh, + dimensionedScalar("one", dimless, 0.01) + ); + + aMesh.setFluxRequired("h"); diff --git a/applications/solvers/finiteArea/liquidFilmFoam/createFvFields.H b/applications/solvers/finiteArea/liquidFilmFoam/createFvFields.H new file mode 100644 index 0000000000..788b155588 --- /dev/null +++ b/applications/solvers/finiteArea/liquidFilmFoam/createFvFields.H @@ -0,0 +1,31 @@ +volVectorField U +( + IOobject + ( + "U", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + mesh, + dimensionedVector("0", dimVelocity, vector::zero) +); + + +volScalarField H +( + IOobject + ( + "H", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + mesh, + dimensionedScalar("0", dimLength, 0) +); + +// Create volume-to surface mapping object +volSurfaceMapping vsm(aMesh); diff --git a/applications/solvers/finiteArea/liquidFilmFoam/liquidFilmFoam.C b/applications/solvers/finiteArea/liquidFilmFoam/liquidFilmFoam.C new file mode 100644 index 0000000000..1f1f8e389d --- /dev/null +++ b/applications/solvers/finiteArea/liquidFilmFoam/liquidFilmFoam.C @@ -0,0 +1,160 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | + \\/ M anipulation | +------------------------------------------------------------------------------- + | Copyright (C) 2016-2017 Wikki Ltd +------------------------------------------------------------------------------- +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 + liquidFilmFoam + +Group + grpFiniteAreaSolvers + +Description + Transient solver for incompressible, laminar flow of Newtonian fluids in + liquid film formulation. + +Author + Zeljko Tukovic, FMENA + Hrvoje Jasak, Wikki Ltd. + +\*---------------------------------------------------------------------------*/ + +#include "fvCFD.H" +#include "faCFD.H" +#include "loopControl.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +int main(int argc, char *argv[]) +{ + #include "setRootCase.H" + #include "createTime.H" + #include "createMesh.H" + #include "createFaMesh.H" + #include "readGravitationalAcceleration.H" + #include "readTransportProperties.H" + #include "createFaFields.H" + #include "createFvFields.H" + #include "createTimeControls.H" + + Info<< "\nStarting time loop\n" << endl; + + while (runTime.run()) + { + #include "readSolutionControls.H" + #include "readTimeControls.H" + #include "surfaceCourantNo.H" + #include "capillaryCourantNo.H" + #include "setDeltaT.H" + + runTime++; + + Info<< "Time = " << runTime.timeName() << nl << endl; + + while (iters.loop()) + { + phi2s = fac::interpolate(h)*phis; + + #include "calcFrictionFactor.H" + + faVectorMatrix UsEqn + ( + fam::ddt(h, Us) + + fam::div(phi2s, Us) + + fam::Sp(0.0125*frictionFactor*mag(Us), Us) + == + Gs*h + - fam::Sp(Sd, Us) + ); + + UsEqn.relax(); + solve(UsEqn == - fac::grad(ps*h)/rhol + ps*fac::grad(h)/rhol); + + areaScalarField UsA(UsEqn.A()); + + Us = UsEqn.H()/UsA; + Us.correctBoundaryConditions(); + + phis = + (fac::interpolate(Us) & aMesh.Le()) + - fac::interpolate(1.0/(rhol*UsA))*fac::lnGrad(ps*h)*aMesh.magLe() + + fac::interpolate(ps/(rhol*UsA))*fac::lnGrad(h)*aMesh.magLe(); + + faScalarMatrix hEqn + ( + fam::ddt(h) + + fam::div(phis, h) + == + Sm + - fam::Sp + ( + Sd/(h + dimensionedScalar("small", dimLength, SMALL)), + h + ) + ); + + hEqn.relax(); + hEqn.solve(); + + phi2s = hEqn.flux(); + + // Bound h + h.primitiveFieldRef() = max + ( + max + ( + h.primitiveField(), + fac::average(max(h, h0))().primitiveField() + *pos(h0.value() - h.primitiveField()) + ), + h0.value() + ); + + ps = rhol*Gn*h - sigma*fac::laplacian(h); + ps.correctBoundaryConditions(); + + Us -= (1.0/(rhol*UsA))*fac::grad(ps*h) + - (ps/(rhol*UsA))*fac::grad(h); + Us.correctBoundaryConditions(); + } + + if (runTime.outputTime()) + { + vsm.mapToVolume(h, H.boundaryFieldRef()); + vsm.mapToVolume(Us, U.boundaryFieldRef()); + + runTime.write(); + } + + Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s" + << " ClockTime = " << runTime.elapsedClockTime() << " s" + << nl << endl; + } + + Info<< "End\n" << endl; + + return 0; +} + + +// ************************************************************************* // diff --git a/applications/solvers/finiteArea/liquidFilmFoam/readSolutionControls.H b/applications/solvers/finiteArea/liquidFilmFoam/readSolutionControls.H new file mode 100644 index 0000000000..3d0c97261c --- /dev/null +++ b/applications/solvers/finiteArea/liquidFilmFoam/readSolutionControls.H @@ -0,0 +1 @@ +loopControl iters(runTime, aMesh.solutionDict(), "solution"); diff --git a/applications/solvers/finiteArea/liquidFilmFoam/readTransportProperties.H b/applications/solvers/finiteArea/liquidFilmFoam/readTransportProperties.H new file mode 100644 index 0000000000..2f2c99ffda --- /dev/null +++ b/applications/solvers/finiteArea/liquidFilmFoam/readTransportProperties.H @@ -0,0 +1,41 @@ +IOdictionary transportProperties +( + IOobject + ( + "transportProperties", + runTime.constant(), + mesh, + IOobject::MUST_READ, + IOobject::NO_WRITE + ) +); + +dimensionedScalar mug +( + transportProperties.lookup("mug") +); + +dimensionedScalar mul +( + transportProperties.lookup("mul") +); + +dimensionedScalar sigma +( + transportProperties.lookup("sigma") +); + +dimensionedScalar rhol +( + transportProperties.lookup("rhol") +); + +dimensionedScalar rhog +( + transportProperties.lookup("rhog") +); + +dimensionedScalar h0 +( + transportProperties.lookup("h0") +); diff --git a/applications/solvers/finiteArea/liquidFilmFoam/surfaceCourantNo.H b/applications/solvers/finiteArea/liquidFilmFoam/surfaceCourantNo.H new file mode 100644 index 0000000000..9e58e23282 --- /dev/null +++ b/applications/solvers/finiteArea/liquidFilmFoam/surfaceCourantNo.H @@ -0,0 +1,63 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | + \\/ M anipulation | +------------------------------------------------------------------------------- + | Copyright (C) 2016-2017 Wikki Ltd +------------------------------------------------------------------------------- +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 + surfaceCourantNo + +Author + Hrvoje Jasak, Wikki Ltd. + +Description + Calculates and outputs the mean and maximum Courant Numbers for the + Finite Area method. + +\*---------------------------------------------------------------------------*/ + +scalar CoNum = 0.0; +scalar meanCoNum = 0.0; +scalar velMag = 0.0; + +if (aMesh.nInternalEdges()) +{ + edgeScalarField SfUfbyDelta + ( + aMesh.edgeInterpolation::deltaCoeffs()*mag(phis) + ); + + CoNum = max(SfUfbyDelta/aMesh.magLe()) + .value()*runTime.deltaT().value(); + + meanCoNum = (sum(SfUfbyDelta)/sum(aMesh.magLe())) + .value()*runTime.deltaT().value(); + + velMag = max(mag(phis)/aMesh.magLe()).value(); +} + +Info<< "Courant Number mean: " << meanCoNum + << " max: " << CoNum + << " velocity magnitude: " << velMag << endl; + + +// ************************************************************************* // diff --git a/applications/solvers/finiteArea/sphereSurfactantFoam/Make/files b/applications/solvers/finiteArea/sphereSurfactantFoam/Make/files new file mode 100644 index 0000000000..574a497031 --- /dev/null +++ b/applications/solvers/finiteArea/sphereSurfactantFoam/Make/files @@ -0,0 +1,3 @@ +surfactantFoam.C + +EXE = $(FOAM_USER_APPBIN)/sphereSurfactantFoam diff --git a/applications/solvers/finiteArea/sphereSurfactantFoam/Make/options b/applications/solvers/finiteArea/sphereSurfactantFoam/Make/options new file mode 100644 index 0000000000..02549914ff --- /dev/null +++ b/applications/solvers/finiteArea/sphereSurfactantFoam/Make/options @@ -0,0 +1,10 @@ +EXE_INC = \ + -I$(LIB_SRC)/finiteArea/lnInclude \ + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/meshTools/lnInclude \ + -I$(LIB_SRC)/cfdTools/general/lnInclude + +EXE_LIBS = \ + -lfiniteArea \ + -lfiniteVolume \ + -lmeshTools diff --git a/applications/solvers/finiteArea/sphereSurfactantFoam/createFaFields.H b/applications/solvers/finiteArea/sphereSurfactantFoam/createFaFields.H new file mode 100644 index 0000000000..4ec416f95b --- /dev/null +++ b/applications/solvers/finiteArea/sphereSurfactantFoam/createFaFields.H @@ -0,0 +1,78 @@ +Info<< "Reading field Cs" << endl; +areaScalarField Cs +( + IOobject + ( + "Cs", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + aMesh +); + +dimensioned Cs0 +( + "Cs0", + dimensionSet(1, -2, 0, 0, 0, 0, 0), + 1.0 +); + +const areaVectorField& R = aMesh.areaCentres(); + +Cs = Cs0*(1.0 + R.component(vector::X)/mag(R)); + + +dimensioned Ds +( + "Ds", + dimensionSet(0, 2, -1, 0, 0, 0, 0), + 1.0 +); + + +areaVectorField Us +( + IOobject + ( + "Us", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + aMesh, + dimensioned("Us", dimVelocity, vector::zero) +); + +dimensioned Uinf("Uinf", dimVelocity, 1.0); + +forAll (Us, faceI) +{ + Us[faceI].x() = + Uinf.value()*(0.25*(3.0 + sqr(R[faceI].x()/mag(R[faceI]))) - 1.0); + + Us[faceI].y() = + Uinf.value()*0.25*R[faceI].x()*R[faceI].y()/sqr(mag(R[faceI])); + + Us[faceI].z() = + Uinf.value()*0.25*R[faceI].x()*R[faceI].z()/sqr(mag(R[faceI])); +} + +Us -= aMesh.faceAreaNormals()*(aMesh.faceAreaNormals() & Us); + + +edgeScalarField phis +( + IOobject + ( + "phis", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + linearEdgeInterpolate(Us) & aMesh.Le() +); + diff --git a/applications/solvers/finiteArea/sphereSurfactantFoam/createFaMesh.H b/applications/solvers/finiteArea/sphereSurfactantFoam/createFaMesh.H new file mode 100644 index 0000000000..905c8e7405 --- /dev/null +++ b/applications/solvers/finiteArea/sphereSurfactantFoam/createFaMesh.H @@ -0,0 +1,2 @@ + // Create Finite Area mesh + faMesh aMesh(mesh); diff --git a/applications/solvers/finiteArea/sphereSurfactantFoam/createVolFields.H b/applications/solvers/finiteArea/sphereSurfactantFoam/createVolFields.H new file mode 100644 index 0000000000..fed56e0409 --- /dev/null +++ b/applications/solvers/finiteArea/sphereSurfactantFoam/createVolFields.H @@ -0,0 +1,36 @@ + // Create volume-to surface mapping object + volSurfaceMapping vsm(aMesh); + + volScalarField Cvf + ( + IOobject + ( + "Cvf", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + mesh, + dimensionedScalar("0", dimless/dimLength, 0) + ); + + vsm.mapToVolume(Cs, Cvf.boundaryFieldRef()); + Cvf.write(); + + volVectorField U + ( + IOobject + ( + "U", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + mesh, + dimensionedVector("zero", dimVelocity, vector::zero) + ); + + vsm.mapToVolume(Us, U.boundaryFieldRef()); + U.write(); diff --git a/applications/solvers/finiteArea/sphereSurfactantFoam/surfactantFoam.C b/applications/solvers/finiteArea/sphereSurfactantFoam/surfactantFoam.C new file mode 100644 index 0000000000..893d46284d --- /dev/null +++ b/applications/solvers/finiteArea/sphereSurfactantFoam/surfactantFoam.C @@ -0,0 +1,91 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | + \\/ M anipulation | +------------------------------------------------------------------------------- + | Copyright (C) 2016-2017 Wikki Ltd +------------------------------------------------------------------------------- +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 + surfactantFoam for sphere transport + +Group + grpFiniteAreaSolvers + +Description + Passive scalar transport equation solver on a sphere + +\*---------------------------------------------------------------------------*/ + +#include "fvCFD.H" +#include "faCFD.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +int main(int argc, char *argv[]) +{ + #include "setRootCase.H" + #include "createTime.H" + #include "createMesh.H" + + #include "createFaMesh.H" + #include "createFaFields.H" + #include "createVolFields.H" + + Info<< "Total mass of surfactant: " + << sum(Cs.internalField()*aMesh.S()) << endl; + + Info<< "\nStarting time loop\n" << endl; + + while (runTime.loop()) + { + Info<< "Time = " << runTime.value() << endl; + + faScalarMatrix CsEqn + ( + fam::ddt(Cs) + + fam::div(phis, Cs) + - fam::laplacian(Ds, Cs) + ); + + CsEqn.solve(); + + if (runTime.writeTime()) + { + vsm.mapToVolume(Cs, Cvf.boundaryFieldRef()); + + runTime.write(); + } + + Info<< "Total mass of surfactant: " + << sum(Cs.internalField()*aMesh.S()) << endl; + + Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s" + << " ClockTime = " << runTime.elapsedClockTime() << " s" + << nl << endl; + } + + Info<< "End\n" << endl; + + return 0; +} + + +// ************************************************************************* // diff --git a/applications/solvers/finiteArea/surfactantFoam/Make/files b/applications/solvers/finiteArea/surfactantFoam/Make/files new file mode 100644 index 0000000000..b4086fd825 --- /dev/null +++ b/applications/solvers/finiteArea/surfactantFoam/Make/files @@ -0,0 +1,3 @@ +surfactantFoam.C + +EXE = $(FOAM_APPBIN)/surfactantFoam diff --git a/applications/solvers/finiteArea/surfactantFoam/Make/options b/applications/solvers/finiteArea/surfactantFoam/Make/options new file mode 100644 index 0000000000..02549914ff --- /dev/null +++ b/applications/solvers/finiteArea/surfactantFoam/Make/options @@ -0,0 +1,10 @@ +EXE_INC = \ + -I$(LIB_SRC)/finiteArea/lnInclude \ + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/meshTools/lnInclude \ + -I$(LIB_SRC)/cfdTools/general/lnInclude + +EXE_LIBS = \ + -lfiniteArea \ + -lfiniteVolume \ + -lmeshTools diff --git a/applications/solvers/finiteArea/surfactantFoam/createFaFields.H b/applications/solvers/finiteArea/surfactantFoam/createFaFields.H new file mode 100644 index 0000000000..a285b21cee --- /dev/null +++ b/applications/solvers/finiteArea/surfactantFoam/createFaFields.H @@ -0,0 +1,63 @@ +Info<< "Reading field Cs" << endl; +areaScalarField Cs +( + IOobject + ( + "Cs", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + aMesh +); + +Info<< "Reading transportProperties\n" << endl; + +IOdictionary transportProperties +( + IOobject + ( + "transportProperties", + runTime.constant(), + mesh, + IOobject::MUST_READ, + IOobject::NO_WRITE + ) +); + + +Info<< "Reading diffusivity D\n" << endl; + +dimensionedScalar Ds +( + transportProperties.lookup("Ds") +); + +areaVectorField Us +( + IOobject + ( + "Us", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::NO_WRITE + ), + aMesh +); + + +edgeScalarField phis +( + IOobject + ( + "phis", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + linearEdgeInterpolate(Us) & aMesh.Le() +); + diff --git a/applications/solvers/finiteArea/surfactantFoam/createFaMesh.H b/applications/solvers/finiteArea/surfactantFoam/createFaMesh.H new file mode 100644 index 0000000000..905c8e7405 --- /dev/null +++ b/applications/solvers/finiteArea/surfactantFoam/createFaMesh.H @@ -0,0 +1,2 @@ + // Create Finite Area mesh + faMesh aMesh(mesh); diff --git a/applications/solvers/finiteArea/surfactantFoam/createVolFields.H b/applications/solvers/finiteArea/surfactantFoam/createVolFields.H new file mode 100644 index 0000000000..fed56e0409 --- /dev/null +++ b/applications/solvers/finiteArea/surfactantFoam/createVolFields.H @@ -0,0 +1,36 @@ + // Create volume-to surface mapping object + volSurfaceMapping vsm(aMesh); + + volScalarField Cvf + ( + IOobject + ( + "Cvf", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + mesh, + dimensionedScalar("0", dimless/dimLength, 0) + ); + + vsm.mapToVolume(Cs, Cvf.boundaryFieldRef()); + Cvf.write(); + + volVectorField U + ( + IOobject + ( + "U", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + mesh, + dimensionedVector("zero", dimVelocity, vector::zero) + ); + + vsm.mapToVolume(Us, U.boundaryFieldRef()); + U.write(); diff --git a/applications/solvers/finiteArea/surfactantFoam/surfactantFoam.C b/applications/solvers/finiteArea/surfactantFoam/surfactantFoam.C new file mode 100644 index 0000000000..630ca8e668 --- /dev/null +++ b/applications/solvers/finiteArea/surfactantFoam/surfactantFoam.C @@ -0,0 +1,114 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | + \\/ M anipulation | +------------------------------------------------------------------------------- + | Copyright (C) 2016-2017 Wikki Ltd +------------------------------------------------------------------------------- +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 + surfactantFoam + +Group + grpFiniteAreaSolvers + +Description + Passive scalar transport equation solver. + + \heading Solver details + The equation is given by: + + \f[ + \ddt{Cs} + \div \left(\vec{U} Cs\right) - \div \left(D_T \grad Cs \right) + = 0 + \f] + + Where: + \vartable + Cs | Passive scalar + Ds | Diffusion coefficient + \endvartable + + \heading Required fields + \plaintable + Cs | Passive scalar + Us | Velocity [m/s] + \endplaintable + +Author + Zeljko Tukovic, FMENA + Hrvoje Jasak, Wikki Ltd. + +\*---------------------------------------------------------------------------*/ + +#include "fvCFD.H" +#include "faCFD.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +int main(int argc, char *argv[]) +{ + #include "setRootCase.H" + #include "createTime.H" + #include "createMesh.H" + #include "createFaMesh.H" + #include "createFaFields.H" + #include "createVolFields.H" + + Info<< "Total mass of surfactant: " + << sum(Cs.internalField()*aMesh.S()) << endl; + + Info<< "\nStarting time loop\n" << endl; + + while (runTime.loop()) + { + Info<< "Time = " << runTime.value() << endl; + + faScalarMatrix CsEqn + ( + fam::ddt(Cs) + + fam::div(phis, Cs) + - fam::laplacian(Ds, Cs) + ); + + CsEqn.solve(); + + if (runTime.writeTime()) + { + vsm.mapToVolume(Cs, Cvf.boundaryFieldRef()); + + runTime.write(); + } + + Info<< "Total mass of surfactant: " + << sum(Cs.internalField()*aMesh.S()) << endl; + + Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s" + << " ClockTime = " << runTime.elapsedClockTime() << " s" + << nl << endl; + } + + Info<< "End\n" << endl; + + return 0; +} + + +// ************************************************************************* // diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C index 155b08a0f8..80e7e78fbb 100644 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C @@ -3,7 +3,7 @@ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation - \\/ M anipulation | + \\/ M anipulation | Copyright (C) 2017 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -47,6 +47,7 @@ Description #include "radiationModel.H" #include "fvOptions.H" #include "coordinateSystem.H" +#include "loopControl.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -89,11 +90,10 @@ int main(int argc, char *argv[]) } } - // --- PIMPLE loop - for (int oCorr=0; oCorr 1) + { + loopControl looping(runTime, pimple, "energyCoupling"); + + while (looping.loop()) + { + Info<< nl << looping << nl; + + forAll(fluidRegions, i) + { + Info<< "\nSolving for fluid region " + << fluidRegions[i].name() << endl; + #include "setRegionFluidFields.H" + #include "readFluidMultiRegionPIMPLEControls.H" + frozenFlow = true; + #include "solveFluid.H" + } + + forAll(solidRegions, i) + { + Info<< "\nSolving for solid region " + << solidRegions[i].name() << endl; + #include "setRegionSolidFields.H" + #include "readSolidMultiRegionPIMPLEControls.H" + #include "solveSolid.H" + } + } + } } runTime.write(); diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/chtMultiRegionSimpleFoam.C b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/chtMultiRegionSimpleFoam.C index ea89befe12..20aa2d59cd 100644 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/chtMultiRegionSimpleFoam.C +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/chtMultiRegionSimpleFoam.C @@ -3,7 +3,7 @@ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation - \\/ M anipulation | + \\/ M anipulation | Copyright (C) 2017 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -42,6 +42,7 @@ Description #include "radiationModel.H" #include "fvOptions.H" #include "coordinateSystem.H" +#include "loopControl.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -57,7 +58,6 @@ int main(int argc, char *argv[]) #include "createFields.H" #include "initContinuityErrs.H" - while (runTime.loop()) { Info<< "Time = " << runTime.timeName() << nl << endl; @@ -80,6 +80,35 @@ int main(int argc, char *argv[]) #include "solveSolid.H" } + // Additional loops for energy solution only + { + loopControl looping(runTime, "SIMPLE", "energyCoupling"); + + while (looping.loop()) + { + Info<< nl << looping << nl; + + forAll(fluidRegions, i) + { + Info<< "\nSolving for fluid region " + << fluidRegions[i].name() << endl; + #include "setRegionFluidFields.H" + #include "readFluidMultiRegionSIMPLEControls.H" + frozenFlow = true; + #include "solveFluid.H" + } + + forAll(solidRegions, i) + { + Info<< "\nSolving for solid region " + << solidRegions[i].name() << endl; + #include "setRegionSolidFields.H" + #include "readSolidMultiRegionSIMPLEControls.H" + #include "solveSolid.H" + } + } + } + runTime.write(); Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s" diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/setRegionFluidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/setRegionFluidFields.H index 6cde0ec06b..119c2ad7ea 100644 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/setRegionFluidFields.H +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/setRegionFluidFields.H @@ -12,6 +12,14 @@ volScalarField& p = thermo.p(); const volScalarField& psi = thermo.psi(); + volScalarField& p_rgh = p_rghFluid[i]; + + const dimensionedVector& g = gFluid[i]; + const volScalarField& gh = ghFluid[i]; + const surfaceScalarField& ghf = ghfFluid[i]; + + radiation::radiationModel& rad = radiation[i]; + IOMRFZoneList& MRF = MRFfluid[i]; fv::options& fvOptions = fluidFvOptions[i]; @@ -22,14 +30,7 @@ initialMassFluid[i] ); - radiation::radiationModel& rad = radiation[i]; + bool frozenFlow = frozenFlowFluid[i]; const label pRefCell = pRefCellFluid[i]; const scalar pRefValue = pRefValueFluid[i]; - const bool frozenFlow = frozenFlowFluid[i]; - - volScalarField& p_rgh = p_rghFluid[i]; - - const dimensionedVector& g = gFluid[i]; - const volScalarField& gh = ghFluid[i]; - const surfaceScalarField& ghf = ghfFluid[i]; diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/solid/solveSolid.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/solid/solveSolid.H index 0600a1c650..2b6306e5c2 100644 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/solid/solveSolid.H +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/solid/solveSolid.H @@ -1,5 +1,5 @@ { - for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) + for (int nonOrth=0; nonOrth<=nNonOrthCorr; ++nonOrth) { fvScalarMatrix hEqn ( @@ -20,9 +20,9 @@ fvOptions.correct(h); } + + thermo.correct(); + + Info<< "Min/max T:" << min(thermo.T()).value() << ' ' + << max(thermo.T()).value() << endl; } - -thermo.correct(); - -Info<< "Min/max T:" << min(thermo.T()).value() << ' ' - << max(thermo.T()).value() << endl; diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H index 87e372fcae..15aba4cb27 100644 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H @@ -19,8 +19,8 @@ List frozenFlowFluid(fluidRegions.size(), false); PtrList MRFfluid(fluidRegions.size()); PtrList fluidFvOptions(fluidRegions.size()); -List