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/utilities/finiteArea/checkFaMesh/Make/files b/applications/utilities/finiteArea/checkFaMesh/Make/files new file mode 100644 index 0000000000..14a2b16396 --- /dev/null +++ b/applications/utilities/finiteArea/checkFaMesh/Make/files @@ -0,0 +1,3 @@ +checkFaMesh.C + +EXE = $(FOAM_APPBIN)/checkFaMesh diff --git a/applications/utilities/finiteArea/checkFaMesh/Make/options b/applications/utilities/finiteArea/checkFaMesh/Make/options new file mode 100644 index 0000000000..7ab09dc4b6 --- /dev/null +++ b/applications/utilities/finiteArea/checkFaMesh/Make/options @@ -0,0 +1,9 @@ +EXE_INC = \ + -I$(LIB_SRC)/finiteArea/lnInclude \ + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/meshTools/lnInclude \ + +EXE_LIBS = \ + -lfiniteVolume \ + -lmeshTools \ + -lfiniteArea diff --git a/applications/utilities/finiteArea/checkFaMesh/checkFaMesh.C b/applications/utilities/finiteArea/checkFaMesh/checkFaMesh.C new file mode 100644 index 0000000000..27b8990e5e --- /dev/null +++ b/applications/utilities/finiteArea/checkFaMesh/checkFaMesh.C @@ -0,0 +1,83 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 + makeFaMesh + +Description + Check a Finite Area mesh + +Author + Zeljko Tukovic, FAMENA + Hrvoje Jasak, Wikki Ltd. + +\*---------------------------------------------------------------------------*/ + +#include "fvCFD.H" +#include "faCFD.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +using namespace Foam; + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +int main(int argc, char *argv[]) +{ +# include "addRegionOption.H" + +# include "setRootCase.H" +# include "createTime.H" +# include "createNamedMesh.H" +# include "createFaMesh.H" + + Info<< "Time = " << runTime.timeName() << nl << endl; + + // General mesh statistics + Info<< "Number of points: " << aMesh.nPoints() << nl + << "Number of internal edges: " << aMesh.nInternalEdges() << nl + << "Number of edges: " << aMesh.nEdges() << nl + << "Number of faces: " << aMesh.nFaces() << nl + << endl; + + // Check geometry + Info<< "Face area: min = " << min(aMesh.S().field()) + << " max = " << max(aMesh.S().field()) << nl + << "Internal edge length: min = " + << min(aMesh.magLe().internalField()) << nl + << " max = " << max(aMesh.magLe().internalField()) << nl + << "Edge length: min = " + << min(aMesh.magLe()).value() << nl + << " max = " << max(aMesh.magLe()).value() << nl + << "Face area normals: min = " << min(aMesh.faceAreaNormals().field()) + << " max = " << max(aMesh.faceAreaNormals().field()) << nl + << endl; + + + return(0); +} + + +// ************************************************************************* // diff --git a/applications/utilities/finiteArea/makeFaMesh/Make/files b/applications/utilities/finiteArea/makeFaMesh/Make/files new file mode 100644 index 0000000000..fbd0f40cd0 --- /dev/null +++ b/applications/utilities/finiteArea/makeFaMesh/Make/files @@ -0,0 +1,3 @@ +makeFaMesh.C + +EXE = $(FOAM_APPBIN)/makeFaMesh diff --git a/applications/utilities/finiteArea/makeFaMesh/Make/options b/applications/utilities/finiteArea/makeFaMesh/Make/options new file mode 100644 index 0000000000..98ba4f7225 --- /dev/null +++ b/applications/utilities/finiteArea/makeFaMesh/Make/options @@ -0,0 +1,8 @@ +EXE_INC = \ + -I$(LIB_SRC)/finiteArea/lnInclude \ + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/cfdTools/general/lnInclude + +EXE_LIBS = \ + -lfiniteArea \ + -lfiniteVolume diff --git a/applications/utilities/finiteArea/makeFaMesh/makeFaMesh.C b/applications/utilities/finiteArea/makeFaMesh/makeFaMesh.C new file mode 100644 index 0000000000..91b499c09b --- /dev/null +++ b/applications/utilities/finiteArea/makeFaMesh/makeFaMesh.C @@ -0,0 +1,355 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 + makeFaMesh + +Description + A mesh generator for finite area mesh. + +Author + Zeljko Tukovic, FAMENA + Hrvoje Jasak, Wikki Ltd. + +\*---------------------------------------------------------------------------*/ + +#include "objectRegistry.H" +#include "Time.H" +#include "argList.H" +#include "OSspecific.H" +#include "faMesh.H" +#include "fvMesh.H" + +using namespace Foam; + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +class faPatchData +{ +public: + word name_; + word type_; + dictionary dict_; + label ownPolyPatchID_; + label ngbPolyPatchID_; + labelList edgeLabels_; + faPatchData() + : + name_(word::null), + type_(word::null), + ownPolyPatchID_(-1), + ngbPolyPatchID_(-1) + {} +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +int main(int argc, char *argv[]) +{ + #include "addRegionOption.H" + argList::noParallel(); + + #include "setRootCase.H" + #include "createTime.H" + #include "createNamedMesh.H" + + // Reading faMeshDefinition dictionary + IOdictionary faMeshDefinition + ( + IOobject + ( + "faMeshDefinition", + runTime.constant(), + "faMesh", + mesh, + IOobject::MUST_READ, + IOobject::NO_WRITE + ) + ); + + wordList polyMeshPatches + ( + faMeshDefinition.lookup("polyMeshPatches") + ); + + const dictionary& bndDict = faMeshDefinition.subDict("boundary"); + + const wordList faPatchNames(bndDict.toc()); + + List faPatches(faPatchNames.size()+1); + + forAll(faPatchNames, patchI) + { + const dictionary& curPatchDict = bndDict.subDict(faPatchNames[patchI]); + + faPatches[patchI].name_ = faPatchNames[patchI]; + + faPatches[patchI].type_ = word(curPatchDict.lookup("type")); + + const word ownName = curPatchDict.lookup("ownerPolyPatch"); + + faPatches[patchI].ownPolyPatchID_ = + mesh.boundaryMesh().findPatchID(ownName); + + if (faPatches[patchI].ownPolyPatchID_ < 0) + { + FatalErrorIn("makeFaMesh:") + << "neighbourPolyPatch " << ownName << " does not exist" + << exit(FatalError); + } + + const word neiName = curPatchDict.lookup("neighbourPolyPatch"); + + faPatches[patchI].ngbPolyPatchID_ = + mesh.boundaryMesh().findPatchID(neiName); + + if (faPatches[patchI].ngbPolyPatchID_ < 0) + { + FatalErrorIn("makeFaMesh:") + << "neighbourPolyPatch " << neiName << " does not exist" + << exit(FatalError); + } + } + + // Setting faceLabels list size + label size = 0; + + labelList patchIDs(polyMeshPatches.size(), -1); + + forAll(polyMeshPatches, patchI) + { + patchIDs[patchI] = + mesh.boundaryMesh().findPatchID(polyMeshPatches[patchI]); + + if (patchIDs[patchI] < 0) + { + FatalErrorIn("makeFaMesh:") + << "Patch " << polyMeshPatches[patchI] << " does not exist" + << exit(FatalError); + } + + size += mesh.boundaryMesh()[patchIDs[patchI]].size(); + } + + labelList faceLabels(size, -1); + + sort(patchIDs); + + + // Filling of faceLabels list + label faceI = -1; + + forAll(polyMeshPatches, patchI) + { + label start = mesh.boundaryMesh()[patchIDs[patchI]].start(); + + label size = mesh.boundaryMesh()[patchIDs[patchI]].size(); + + for (label i = 0; i < size; ++i) + { + faceLabels[++faceI] = start + i; + } + } + + // Creating faMesh + Info << "Create faMesh ... "; + + faMesh areaMesh + ( + mesh, + faceLabels + ); + Info << "Done" << endl; + + + // Determination of faPatch ID for each boundary edge. + // Result is in the bndEdgeFaPatchIDs list + const indirectPrimitivePatch& patch = areaMesh.patch(); + + labelList faceCells(faceLabels.size(), -1); + + forAll(faceCells, faceI) + { + label faceID = faceLabels[faceI]; + + faceCells[faceI] = mesh.faceOwner()[faceID]; + } + + labelList meshEdges = + patch.meshEdges + ( + mesh.edges(), + mesh.cellEdges(), + faceCells + ); + + const labelListList& edgeFaces = mesh.edgeFaces(); + + const label nTotalEdges = patch.nEdges(); + const label nInternalEdges = patch.nInternalEdges(); + + labelList bndEdgeFaPatchIDs(nTotalEdges - nInternalEdges, -1); + + for (label edgeI = nInternalEdges; edgeI < nTotalEdges; ++edgeI) + { + label curMeshEdge = meshEdges[edgeI]; + + labelList curEdgePatchIDs(2, label(-1)); + + label patchI = -1; + + forAll(edgeFaces[curMeshEdge], faceI) + { + label curFace = edgeFaces[curMeshEdge][faceI]; + + label curPatchID = mesh.boundaryMesh().whichPatch(curFace); + + if (curPatchID != -1) + { + curEdgePatchIDs[++patchI] = curPatchID; + } + } + + for (label pI = 0; pI < faPatches.size() - 1; ++pI) + { + if + ( + ( + curEdgePatchIDs[0] == faPatches[pI].ownPolyPatchID_ + && curEdgePatchIDs[1] == faPatches[pI].ngbPolyPatchID_ + ) + || + ( + curEdgePatchIDs[1] == faPatches[pI].ownPolyPatchID_ + && curEdgePatchIDs[0] == faPatches[pI].ngbPolyPatchID_ + ) + ) + { + bndEdgeFaPatchIDs[edgeI - nInternalEdges] = pI; + break; + } + } + } + + + // Set edgeLabels for each faPatch + for (label pI=0; pI<(faPatches.size()-1); ++pI) + { + SLList