Commit 0c646223 authored by Hrvoje Jasak's avatar Hrvoje Jasak Committed by Andrew Heather
Browse files

Finite area port, Hrvoje Jasak

- with sphereSurfactantFoam and sphereTransport test case
parent 4c81ee20
liquidFilmFoam.C
EXE = $(FOAM_APPBIN)/liquidFilmFoam
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
{
// 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));
}
{
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;
}
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.schemesDict().setFluxRequired("h");
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);
/*---------------------------------------------------------------------------*\
========= |
\\ / 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 <http://www.gnu.org/licenses/>.
Application
liquidFilmFoam
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"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
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 (runTime.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 = "
<< scalar(runTime.elapsedCpuTime())
<< " s\n" << endl << endl;
}
return 0;
}
// ************************************************************************* //
int nCorr = 0;
if (aMesh.solutionDict().found("nOuterCorrectors"))
{
nCorr =
readInt(aMesh.solutionDict().lookup("nOuterCorrectors"));
}
else
{
FatalErrorIn(args.executable())
<< "Cannot find number of correctors"
<< abort(FatalError);
}
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")
);
/*---------------------------------------------------------------------------*\
========= |
\\ / 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 <http://www.gnu.org/licenses/>.
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;
// ************************************************************************* //
surfactantFoam.C
EXE = $(FOAM_USER_APPBIN)/sphereSurfactantFoam
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
Info<< "Reading field Cs" << endl;
areaScalarField Cs
(
IOobject
(
"Cs",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
aMesh
);
dimensioned<scalar> 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<scalar> 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<vector>("Us", dimVelocity, vector::zero)
);
dimensioned<scalar> 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
),