Skip to content

Commit

Permalink
Adapted NavierStokes to new narrow band
Browse files Browse the repository at this point in the history
  • Loading branch information
tmaric committed May 20, 2023
1 parent e65085a commit bb5cd0c
Show file tree
Hide file tree
Showing 15 changed files with 396 additions and 74 deletions.
Original file line number Diff line number Diff line change
@@ -0,0 +1,6 @@
fvScalarMatrix YoungLaplaceEqn
(
fvm::laplacian(p_rgh) == fvc::div(fSigma->faceSurfaceTensionForce()*mesh.magSf())
);
YoungLaplaceEqn.solve();

2 changes: 2 additions & 0 deletions applications/solvers/leiaLevelSetTwoPhaseFoam/alphaEqn.H
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,8 @@

phaseInd->calcPhaseIndicator(alpha1, psi);

narrowBand->calc();

// Making U absolute again after advection step
if (mesh.moving())
{
Expand Down
5 changes: 5 additions & 0 deletions applications/solvers/leiaLevelSetTwoPhaseFoam/createFields.H
Original file line number Diff line number Diff line change
Expand Up @@ -147,6 +147,11 @@ volScalarField psi
// Leve Set re-distancing model
autoPtr<redistancer> redist = redistancer::New(mesh);

// Level Set narrowBand model
autoPtr<NarrowBand> narrowBand = NarrowBand::New(mesh, psi);
narrowBand->calc();
narrowBand->write();

// Level Set phase-indicator model
autoPtr<phaseIndicator> phaseInd = phaseIndicator::New(mesh);

Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -62,6 +62,7 @@ Description
#include "phaseIndicator.H"
#include "redistancer.H"
#include "surfaceTensionForce.H"
#include "NarrowBand.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

Expand Down Expand Up @@ -89,6 +90,9 @@ int main(int argc, char *argv[])

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl;

#include "YoungLaplaceEqn.H"
p_rgh.write();

while (runTime.run())
{
Expand Down
2 changes: 2 additions & 0 deletions src/leiaLevelSet/Make/files
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,9 @@ redistancer/redistancer.C
redistancer/lsqRedistancer.C

surfaceTensionForce/surfaceTensionForce.C
surfaceTensionForce/constantCurvatureSurfTension.C
surfaceTensionForce/traceGradGeoNormalSnGradAlpha.C
surfaceTensionForce/traceGradGradPsiSnGradAlpha.C

SDPLSSource/SDPLSSource.C
SDPLSSource/SDPLS_R.C
Expand Down
56 changes: 0 additions & 56 deletions src/leiaLevelSet/phaseIndicator/geometricPhaseIndicator.C
Original file line number Diff line number Diff line change
Expand Up @@ -54,20 +54,6 @@ addToRunTimeSelectionTable(phaseIndicator, geometricPhaseIndicator, Mesh);
geometricPhaseIndicator::geometricPhaseIndicator(const fvMesh& mesh)
:
phaseIndicator(mesh),
// narrowBandTmp_(new volScalarField
// (
// IOobject
// (
// "narrowBand",
// mesh.time().timeName(),
// mesh,
// IOobject::NO_READ,
// IOobject::AUTO_WRITE
// ),
// mesh,
// dimensionedScalar("narrowBand", dimless, 0)
// )
// ),
narrowBand_(mesh.lookupObject<volScalarField>("NarrowBand")),
ncTmp_(new volVectorField
(
Expand Down Expand Up @@ -121,48 +107,6 @@ void geometricPhaseIndicator::calcPhaseIndicator

// Approach: Lookup existing narrowBand
const volScalarField& narrowBand = narrowBand_;

// // Approch: Calc own narrowBand
// volScalarField& narrowBand = narrowBandTmp_.ref();
// const auto& own = mesh.owner();
// const auto& nei = mesh.neighbour();
// narrowBand = dimensionedScalar("narrowBand", narrowBand.dimensions(), 0.);

// // Select the cells in the narrow band using face-connectivity.
// forAll(own, faceI)
// {
// // FIXME: If a cell is clipped by the interface at one of its corners,
// // and the surrounding cell centers do not switch the iso-sign, it
// // may not be picked up by this. May cause numerical inconsistency. TM.
// if (psi[own[faceI]] * psi[nei[faceI]] <= 0)
// {
// narrowBand[own[faceI]] = 1;
// narrowBand[nei[faceI]] = 1;
// }
// }
// // Set narrow band values across coupled process boundaries.
// const auto& psiBdryField = psi.boundaryField(); // needed for Nc LLSQ contribs
// const auto& patches = mesh.boundary(); // needed for Nc LLSQ contribs
// const auto& faceOwner = mesh.faceOwner();
// forAll(psiBdryField, patchI)
// {
// const fvPatch& patch = patches[patchI];
// if (isA<coupledFvPatch>(patch)) // coupled patch
// {
// const auto& psiPatchField = psiBdryField[patchI];
// auto psiPatchNeiFieldTmp =
// psiPatchField.patchNeighbourField();
// const auto& psiPatchNeiField = psiPatchNeiFieldTmp();
// forAll(psiPatchNeiField, faceI)
// {
// label faceJ = faceI + patch.start(); // Global face label.
// if (psi[faceOwner[faceJ]] * psiPatchNeiField[faceI] <= 0)
// {
// narrowBand[faceOwner[faceJ]] = 1;
// }
// }
// }
// }

// Linear Least-Squares Approximation of \Psi
// \psi(x,y,z)= nc_x x + nc_y y + nc_z z + dc
Expand Down
5 changes: 0 additions & 5 deletions src/leiaLevelSet/phaseIndicator/geometricPhaseIndicator.H
Original file line number Diff line number Diff line change
Expand Up @@ -64,7 +64,6 @@ class geometricPhaseIndicator
protected:

// Narrow-band marker
// tmp<volScalarField> narrowBandTmp_;
const volScalarField& narrowBand_;

// Linear Least-Squares Approximation of \Psi
Expand Down Expand Up @@ -93,10 +92,6 @@ public:
const volScalarField& psi
);

// tmp<volScalarField> narrowBand() const
// {
// return narrowBandTmp_;
// }
};

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Expand Down
Original file line number Diff line number Diff line change
@@ -0,0 +1,75 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2023 Tomislav Maric, TU Darmstadt
-------------------------------------------------------------------------------
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/>.
\*---------------------------------------------------------------------------*/

#include "fvPatchFieldsFwd.H"
#include "primitiveFieldsFwd.H"
#include "surfaceTensionForce.H"
#include "constantCurvatureSurfaceTension.H"
#include "addToRunTimeSelectionTable.H"
#include "volFields.H"
#include "surfaceFields.H"
#include "fvcGrad.H"
#include "fvcSnGrad.H"
#include "surfaceInterpolate.H"
#include "volFieldsFwd.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

namespace Foam
{

defineTypeNameAndDebug(constantCurvatureSurfaceTension, false);
addToRunTimeSelectionTable(surfaceTensionForce, constantCurvatureSurfaceTension, Mesh);

// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
constantCurvatureSurfaceTension::constantCurvatureSurfaceTension(const fvMesh& mesh)
:
surfaceTensionForce(mesh),
fvSolutionDict_(mesh_),
levelSetDict_(fvSolutionDict_.subDict("levelSet")),
surfTensionDict_(levelSetDict_.subDict("surfaceTensionForce")),
curvature_
(
"curvature",
pow(dimLength, -1),
surfTensionDict_.get<scalar>("curvature")
),
alpha_(mesh_.lookupObject<volScalarField>(surfTensionDict_.getOrDefault<word>("alpha", "alpha.dispersed")))
{}

// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //

tmp<surfaceScalarField> constantCurvatureSurfaceTension::faceSurfaceTensionForce() const
{
return sigma_ * curvature_ * fvc::snGrad(alpha_);
}

} // End namespace Foam

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

// ************************************************************************* //
Original file line number Diff line number Diff line change
@@ -0,0 +1,91 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2023 Tomislav Maric, TU Darmstadt
-------------------------------------------------------------------------------
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/>.
Class
Foam::constantCurvatureSurfaceTension
Description
A constant-curvature surface-tension force for testing if the pressure-
velocity coupling algorithm is force-balanced.
SourceFiles
constantCurvatureSurfaceTension.C
\*---------------------------------------------------------------------------*/

#ifndef constantCurvatureSurfaceTension_H
#define constantCurvatureSurfaceTension_H

#include "IOdictionary.H"
#include "surfaceTensionForce.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

namespace Foam
{

/*---------------------------------------------------------------------------*\
Class constantCurvatureSurfaceTension Declaration
\*---------------------------------------------------------------------------*/

class constantCurvatureSurfaceTension
:
public surfaceTensionForce
{
protected:

const fvSolution& fvSolutionDict_;
const dictionary& levelSetDict_;
const dictionary& surfTensionDict_;

// Volume fraction
const dimensionedScalar curvature_;
const volScalarField& alpha_;

public:

// Static Data Members
TypeName ("constantCurvatureSurfaceTension");

// Constructors
constantCurvatureSurfaceTension(const fvMesh& mesh);

//- Destructor
virtual ~constantCurvatureSurfaceTension() = default;

// Member Functions
virtual tmp<surfaceScalarField> faceSurfaceTensionForce() const;
};

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

} // End namespace Foam

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

#endif

// ************************************************************************* //
19 changes: 18 additions & 1 deletion src/leiaLevelSet/surfaceTensionForce/surfaceTensionForce.C
Original file line number Diff line number Diff line change
Expand Up @@ -67,7 +67,24 @@ Foam::surfaceTensionForce::New(const fvMesh& mesh)
surfaceTensionForce::surfaceTensionForce(const fvMesh& mesh)
:
mesh_(mesh),
runTime_(mesh.time())
runTime_(mesh.time()),
transportProperties_
(
IOobject
(
"transportProperties",
"constant",
runTime_,
IOobject::MUST_READ,
IOobject::NO_WRITE
)
),
sigma_
(
"sigma",
dimForce / dimLength,
transportProperties_.get<scalar>("sigma")
)
{}

// ************************************************************************* //
Expand Down
4 changes: 4 additions & 0 deletions src/leiaLevelSet/surfaceTensionForce/surfaceTensionForce.H
Original file line number Diff line number Diff line change
Expand Up @@ -62,6 +62,10 @@ protected:
const fvMesh& mesh_;

const Time& runTime_;

// Surface tension coefficient.
IOdictionary transportProperties_;
dimensionedScalar sigma_;

public:

Expand Down
Loading

0 comments on commit bb5cd0c

Please sign in to comment.