diff --git a/applications/solvers/incompressible/pimpleFoam/Make/options b/applications/solvers/incompressible/pimpleFoam/Make/options index 584c27112d..572e65531f 100644 --- a/applications/solvers/incompressible/pimpleFoam/Make/options +++ b/applications/solvers/incompressible/pimpleFoam/Make/options @@ -7,7 +7,8 @@ EXE_INC = \ -I$(LIB_SRC)/transportModels \ -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \ -I$(LIB_SRC)/dynamicMesh/lnInclude \ - -I$(LIB_SRC)/dynamicFvMesh/lnInclude + -I$(LIB_SRC)/dynamicFvMesh/lnInclude \ + -I$(LIB_SRC)/dynamicFvMesh/interfaceTrackingFvMesh/lnInclude EXE_LIBS = \ -lfiniteVolume \ @@ -20,4 +21,5 @@ EXE_LIBS = \ -ldynamicMesh \ -ldynamicFvMesh \ -ltopoChangerFvMesh \ - -latmosphericModels + -latmosphericModels \ + -linterfaceTrackingFvMesh diff --git a/src/Allwmake b/src/Allwmake index 648deae6d7..147c3c080a 100755 --- a/src/Allwmake +++ b/src/Allwmake @@ -98,6 +98,8 @@ wmake $targetType waveModels wmake $targetType engine conversion/Allwmake $targetType $* + +phaseSystemModels/Allwmake $targetType $* functionObjects/Allwmake $targetType $* wmake $targetType lumpedPointMotion @@ -109,7 +111,10 @@ wmake $targetType semiPermeableBaffle wmake $targetType atmosphericModels wmake $targetType optimisation/adjointOptimisation/adjoint -phaseSystemModels/Allwmake $targetType $* + +# interfaceTracking libs +wmake dynamicFvMesh/interfaceTrackingFvMesh + # Needs access to Turbulence diff --git a/src/dynamicFvMesh/interfaceTrackingFvMesh/Make/files b/src/dynamicFvMesh/interfaceTrackingFvMesh/Make/files new file mode 100644 index 0000000000..4c6c1f79d2 --- /dev/null +++ b/src/dynamicFvMesh/interfaceTrackingFvMesh/Make/files @@ -0,0 +1,9 @@ +interfaceTrackingFvMesh.C +freeSurfacePointDisplacement.C +fvPatchFields/freeSurfacePressure/freeSurfacePressureFvPatchScalarField.C +fvPatchFields/freeSurfaceVelocity/freeSurfaceVelocityFvPatchVectorField.C + +functionObjects/pointHistory/pointHistory.C +functionObjects/writeFreeSurface/writeFreeSurface.C + +LIB = $(FOAM_LIBBIN)/libinterfaceTrackingFvMesh diff --git a/src/dynamicFvMesh/interfaceTrackingFvMesh/Make/options b/src/dynamicFvMesh/interfaceTrackingFvMesh/Make/options new file mode 100644 index 0000000000..e3cea0e094 --- /dev/null +++ b/src/dynamicFvMesh/interfaceTrackingFvMesh/Make/options @@ -0,0 +1,23 @@ +EXE_INC = \ + -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \ + -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \ + -I$(LIB_SRC)/transportModels \ + -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \ + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/finiteArea/lnInclude \ + -I$(LIB_SRC)/meshTools/lnInclude \ + -I$(LIB_SRC)/dynamicFvMesh/lnInclude \ + -I$(LIB_SRC)/dynamicMesh/lnInclude \ + -I$(LIB_SRC)/fvMotionSolver/lnInclude + + +LIB_LIBS = \ + -lturbulenceModels \ + -lincompressibleTurbulenceModels \ + -lincompressibleTransportModels \ + -lfiniteVolume \ + -lfiniteArea \ + -lmeshTools \ + -ldynamicFvMesh \ + -lfvMotionSolvers \ + -ldynamicMesh diff --git a/src/dynamicFvMesh/interfaceTrackingFvMesh/boundaryProcessorFaPatchPoints.H b/src/dynamicFvMesh/interfaceTrackingFvMesh/boundaryProcessorFaPatchPoints.H new file mode 100644 index 0000000000..d83f117875 --- /dev/null +++ b/src/dynamicFvMesh/interfaceTrackingFvMesh/boundaryProcessorFaPatchPoints.H @@ -0,0 +1,44 @@ +// Boundary processor patch points +{ + const labelList& curPointEdges = pointEdges[curPoint]; + + label patchID = -1; + label edgeID = -1; + + for (const label curEdgei : curPointEdges) + { + if (edgeFaces[curEdgei].size() == 1) + { + forAll(aMesh().boundary(), patchI) + { + const labelList& curEdges = aMesh().boundary()[patchI]; + + label index = curEdges.find(curEdgei); + + if (index != -1) + { + if + ( + aMesh().boundary()[patchI].type() + != processorFaPatch::typeName + ) + { + patchID = patchI; + edgeID = index; + break; + } + } + } + } + } + + if (patchID != -1) + { + vector mirrorCtrlPoint = + patchMirrorPoints[patchID][edgeID]; + + lsPoints[curPatchPoint].setSize(lsPoints[curPatchPoint].size() + 1); + lsPoints[curPatchPoint][lsPoints[curPatchPoint].size() - 1] = + mirrorCtrlPoint; + } +} diff --git a/src/dynamicFvMesh/interfaceTrackingFvMesh/freeSurfacePointDisplacement.C b/src/dynamicFvMesh/interfaceTrackingFvMesh/freeSurfacePointDisplacement.C new file mode 100644 index 0000000000..c11c4f91bb --- /dev/null +++ b/src/dynamicFvMesh/interfaceTrackingFvMesh/freeSurfacePointDisplacement.C @@ -0,0 +1,623 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2019 Zeljko Tukovic, FSB Zagreb. +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +\*---------------------------------------------------------------------------*/ + +#include "interfaceTrackingFvMesh.H" +#include "primitivePatchInterpolation.H" +#include "emptyFaPatch.H" +#include "wedgeFaPatch.H" +#include "wallFvPatch.H" +#include "PstreamCombineReduceOps.H" +#include "coordinateSystem.H" +#include "unitConversion.H" +#include "scalarMatrices.H" +#include "tensor2D.H" + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +Foam::tmp +Foam::interfaceTrackingFvMesh::pointDisplacement() +{ + const pointField& points = aMesh().patch().localPoints(); + const labelListList& pointFaces = aMesh().patch().pointFaces(); + + auto tdisplacement = tmp::New(points.size(), Zero); + auto& displacement = tdisplacement.ref(); + + // Calculate displacement of internal points + const vectorField& pointNormals = aMesh().pointAreaNormals(); + const edgeList& edges = aMesh().patch().edges(); + labelList internalPoints = aMesh().internalPoints(); + + for (const label curPoint : internalPoints) + { + const labelList& curPointFaces = pointFaces[curPoint]; + + vectorField lsPoints(curPointFaces.size(), Zero); + + for (label i=0; i patchMirrorPoints(aMesh().boundary().size()); + + // Old faMesh points + vectorField oldPoints(aMesh().nPoints(), Zero); + const labelList& meshPoints = aMesh().patch().meshPoints(); + forAll(oldPoints, pI) + { + oldPoints[pI] = + mesh().oldPoints()[meshPoints[pI]]; + } + + forAll(patchMirrorPoints, patchI) + { + patchMirrorPoints.set + ( + patchI, + new vectorField + ( + aMesh().boundary()[patchI].faPatch::size(), + Zero + ) + ); + + vectorField N + ( + aMesh().boundary()[patchI].ngbPolyPatchFaceNormals() + ); + + const labelList& eFaces = + aMesh().boundary()[patchI].edgeFaces(); + + // Correct N according to specified contact angle + if (contactAnglePtr_) + { + label ngbPolyPatchID = + aMesh().boundary()[patchI].ngbPolyPatchIndex(); + + if (ngbPolyPatchID != -1) + { + if + ( + mesh().boundary()[ngbPolyPatchID].type() + == wallFvPatch::typeName + ) + { + // Info<< aMesh().boundary()[patchI].name() << endl; + + scalar rotAngle = degToRad + ( + gAverage + ( + 90 + - contactAnglePtr_->boundaryField()[patchI] + ) + ); + + const vectorField& pEdgN = + aMesh().edgeAreaNormals().boundaryField()[patchI]; + + vectorField rotationAxis( N^pEdgN ); + + const edgeList::subList patchEdges = + aMesh().boundary()[patchI].patchSlice(aMesh().edges()); + + forAll(rotationAxis, edgeI) + { + vector e = patchEdges[edgeI].vec(oldPoints); + // vector e = patchEdges[edgeI].vec(aMesh().points()); + + // Adjust direction + rotationAxis[edgeI] = + e*(e&rotationAxis[edgeI]) + /mag((e&rotationAxis[edgeI])); + } + rotationAxis /= mag(rotationAxis) + SMALL; + + vectorField rotationAxis2 = rotationAxis; + forAll(rotationAxis2, edgeI) + { + rotationAxis2[edgeI] = + (N[edgeI]^facesDisplacementDir()[eFaces[edgeI]]); + + // Adjust direction + rotationAxis2[edgeI] = + rotationAxis2[edgeI] + *(rotationAxis2[edgeI]&rotationAxis[edgeI]) + /( + mag((rotationAxis2[edgeI]&rotationAxis[edgeI])) + + SMALL + ); + } + rotationAxis2 /= mag(rotationAxis2) + SMALL; + + // Rodrigues' rotation formula + N = N*cos(rotAngle) + + rotationAxis*(rotationAxis & N)*(1 - cos(rotAngle)) + + (rotationAxis^N)*sin(rotAngle); + + N /= mag(N) + SMALL; + + N = (rotationAxis^N); + + N = (N^rotationAxis2); + + N /= mag(N) + SMALL; + + // Info<< N << endl; + } + } + } + + const labelList peFaces = + labelList::subList + ( + aMesh().edgeOwner(), + aMesh().boundary()[patchI].faPatch::size(), + aMesh().boundary()[patchI].start() + ); + + const labelList& pEdges = aMesh().boundary()[patchI]; + + vectorField peCentres(pEdges.size(), Zero); + forAll(peCentres, edgeI) + { + peCentres[edgeI] = + edges[pEdges[edgeI]].centre(points); + } + + vectorField delta + ( + vectorField(controlPoints(), peFaces) + - peCentres + ); + + // Info<< aMesh().boundary()[patchI].name() << endl; + // Info<< vectorField(controlPoints(), peFaces) << endl; + + patchMirrorPoints[patchI] = + peCentres + ((I - 2*N*N)&delta); + + // Info<< patchMirrorPoints[patchI] << endl; + } + + // Calculate displacement of boundary points + labelList boundaryPoints = aMesh().boundaryPoints(); + + const labelListList& edgeFaces = aMesh().patch().edgeFaces(); + const labelListList& pointEdges = aMesh().patch().pointEdges(); + + for (const label curPoint : boundaryPoints) + { + if (motionPointsMask()[curPoint] == 1) + { + // Calculating mirror points + const labelList& curPointEdges = pointEdges[curPoint]; + + vectorField mirrorPoints(2, Zero); + + label counter = -1; + + forAll(curPointEdges, edgeI) + { + label curEdge = curPointEdges[edgeI]; + + if(edgeFaces[curEdge].size() == 1) + { + label patchID = -1; + label edgeID = -1; + forAll(aMesh().boundary(), patchI) + { + const labelList& pEdges = + aMesh().boundary()[patchI]; + label index = pEdges.find(curEdge); + if (index != -1) + { + patchID = patchI; + edgeID = index; + break; + } + } + + mirrorPoints[++counter] = + patchMirrorPoints[patchID][edgeID]; + } + } + + // Calculating LS plane fit + const labelList& curPointFaces = pointFaces[curPoint]; + + vectorField lsPoints + ( + curPointFaces.size() + mirrorPoints.size(), + Zero + ); + + counter = -1; + + for (label i=0; i(aMesh().boundary()[patchI]); + + const labelList& patchPointLabels = + procPatch.pointLabels(); + + FieldField lsPoints(patchPointLabels.size()); + forAll(lsPoints, pointI) + { + lsPoints.set(pointI, new vectorField(0, Zero)); + } + + const labelList& nonGlobalPatchPoints = + procPatch.nonGlobalPatchPoints(); + + forAll(nonGlobalPatchPoints, pointI) + { + label curPatchPoint = + nonGlobalPatchPoints[pointI]; + + label curPoint = + patchPointLabels[curPatchPoint]; + + const labelList& curPointFaces = pointFaces[curPoint]; + + lsPoints[curPatchPoint].setSize(curPointFaces.size()); + + forAll(curPointFaces, faceI) + { + label curFace = curPointFaces[faceI]; + + lsPoints[curPatchPoint][faceI] = controlPoints()[curFace]; + } + + #include "boundaryProcessorFaPatchPoints.H" + } + + // Parallel data exchange + { + OPstream toNeighbProc + ( + Pstream::commsTypes::blocking, + procPatch.neighbProcNo() + ); + + toNeighbProc << lsPoints; + } + + FieldField ngbLsPoints(patchPointLabels.size()); + { + IPstream fromNeighbProc + ( + Pstream::commsTypes::blocking, + procPatch.neighbProcNo() + ); + + fromNeighbProc >> ngbLsPoints; + } + + forAll(nonGlobalPatchPoints, pointI) + { + label curPatchPoint = + nonGlobalPatchPoints[pointI]; + + label curPoint = + patchPointLabels[curPatchPoint]; + + label curNgbPoint = procPatch.neighbPoints()[curPatchPoint]; + + vectorField allLsPoints + ( + lsPoints[curPatchPoint].size() + + ngbLsPoints[curNgbPoint].size(), + Zero + ); + + label counter = -1; + forAll(lsPoints[curPatchPoint], pointI) + { + allLsPoints[++counter] = lsPoints[curPatchPoint][pointI]; + } + forAll(ngbLsPoints[curNgbPoint], pointI) + { + allLsPoints[++counter] = ngbLsPoints[curNgbPoint][pointI]; + } + + vectorField pointAndNormal + ( + lsPlanePointAndNormal + ( + allLsPoints, + points[curPoint], + pointNormals[curPoint] + ) + ); + + vector& P = pointAndNormal[0]; + vector& N = pointAndNormal[1]; + + if (motionPointsMask()[curPoint] != 0) + { + displacement[curPoint] = + pointsDisplacementDir()[curPoint] + *((P - points[curPoint])&N) + /(pointsDisplacementDir()[curPoint]&N); + } + } + } + } + + + // Calculate displacement of global processor patch points + if (aMesh().globalData().nGlobalPoints() > 0) + { + const labelList& spLabels = + aMesh().globalData().sharedPointLabels(); + + const labelList& addr = aMesh().globalData().sharedPointAddr(); + + for (label k=0; k > procLsPoints(Pstream::nProcs()); + + label curSharedPointIndex = addr.find(k); + + if (curSharedPointIndex != -1) + { + label curPoint = spLabels[curSharedPointIndex]; + + const labelList& curPointFaces = pointFaces[curPoint]; + + procLsPoints[Pstream::myProcNo()] = + List(curPointFaces.size()); + + forAll(curPointFaces, faceI) + { + label curFace = curPointFaces[faceI]; + + procLsPoints[Pstream::myProcNo()][faceI] = + controlPoints()[curFace]; + } + } + + Pstream::gatherList(procLsPoints); + Pstream::scatterList(procLsPoints); + + if (curSharedPointIndex != -1) + { + label curPoint = spLabels[curSharedPointIndex]; + + label nAllPoints = 0; + forAll(procLsPoints, procI) + { + nAllPoints += procLsPoints[procI].size(); + } + + vectorField allPoints(nAllPoints, Zero); + + label counter = 0; + forAll(procLsPoints, procI) + { + forAll(procLsPoints[procI], pointI) + { + allPoints[counter++] = + procLsPoints[procI][pointI]; + } + } + + vectorField pointAndNormal + ( + lsPlanePointAndNormal + ( + allPoints, + points[curPoint], + pointNormals[curPoint] + ) + ); + + const vector& P = pointAndNormal[0]; + const vector& N = pointAndNormal[1]; + + displacement[curPoint] = + pointsDisplacementDir()[curPoint] + *((P - points[curPoint])&N) + /(pointsDisplacementDir()[curPoint]&N); + } + } + } + + return tdisplacement; +} + + +Foam::tmp +Foam::interfaceTrackingFvMesh::lsPlanePointAndNormal +( + const vectorField& points, + const vector& origin, + const vector& axis +) const +{ + // LS in local CS + vector dir = (points[0] - origin); + dir -= axis*(axis&dir); + dir /= mag(dir); + coordinateSystem cs("cs", origin, axis, dir); + + vectorField localPoints(cs.localPosition(points)); + vector avgLocalPoint = average(localPoints); + + // scalarField W = 1.0/(mag(points - origin) + SMALL); + scalarField W(points.size(), scalar(1)); + + const label nCoeffs = 2; + scalarRectangularMatrix M(points.size(), nCoeffs, Zero); + + scalar L = 2*max(mag(localPoints-avgLocalPoint)); + for (label i=0; i::New(2); + auto& pointAndNormal = tpointAndNormal.ref(); + + pointAndNormal[0] = p0; + pointAndNormal[1] = n0; + + return tpointAndNormal; +} + + +// ************************************************************************* // diff --git a/src/dynamicFvMesh/interfaceTrackingFvMesh/functionObjects/pointHistory/pointHistory.C b/src/dynamicFvMesh/interfaceTrackingFvMesh/functionObjects/pointHistory/pointHistory.C new file mode 100644 index 0000000000..96470ec5fe --- /dev/null +++ b/src/dynamicFvMesh/interfaceTrackingFvMesh/functionObjects/pointHistory/pointHistory.C @@ -0,0 +1,215 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2019 Zeljko Tukovic, FSB Zagreb. +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +\*----------------------------------------------------------------------------*/ + +#include "pointHistory.H" +#include "addToRunTimeSelectionTable.H" +#include "IOmanip.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ + defineTypeNameAndDebug(pointHistory, 0); + + addToRunTimeSelectionTable + ( + functionObject, + pointHistory, + dictionary + ); +} + + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +bool Foam::pointHistory::writeData() +{ + const fvMesh& mesh = + time_.lookupObject(polyMesh::defaultRegion); + + vector position(Zero); + + if (processor_ == Pstream::myProcNo()) + { + position = mesh.points()[historyPointID_]; + } + + reduce(position, sumOp()); + + if (Pstream::master()) + { + historyFilePtr_() << setprecision(12); + + historyFilePtr_() + << time_.time().value() << tab + << position.x() << tab + << position.y() << tab + << position.z() << endl; + } + + return true; +} + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::pointHistory::pointHistory +( + const word& name, + const Time& runTime, + const dictionary& dict +) +: + functionObject(name), + name_(name), + time_(runTime), + regionName_(polyMesh::defaultRegion), + historyPointID_(-1), + refHistoryPoint_(dict.lookup("refHistoryPoint")), + processor_(-1), + fileName_(dict.get("fileName")), + historyFilePtr_() +{ + Info<< "Creating " << this->name() << " function object." << endl; + + dict.readIfPresent("region", regionName_); + dict.readIfPresent("historyPointID", historyPointID_); + + const fvMesh& mesh = + time_.lookupObject(regionName_); + + const vectorField& points = mesh.points(); + + List minDist(Pstream::nProcs(), GREAT); + + if (historyPointID_ == -1) + { + forAll(points, pointI) + { + scalar dist = mag(refHistoryPoint_ - points[pointI]); + + if (dist < minDist[Pstream::myProcNo()]) + { + minDist[Pstream::myProcNo()] = dist; + historyPointID_ = pointI; + } + } + } + + Pstream::gatherList(minDist); + Pstream::scatterList(minDist); + + processor_ = -1; + scalar min = GREAT; + + forAll(minDist, procI) + { + if (minDist[procI] < min) + { + min = minDist[procI]; + processor_ = procI; + } + } + + if (processor_ == Pstream::myProcNo()) + { + Pout<< "History point ID: " << historyPointID_ << nl + << "History point coordinates: " + << points[historyPointID_] << nl + << "Reference point coordinates: " << refHistoryPoint_ + << endl; + } + + // Create history file if not already created + if (historyFilePtr_.empty()) + { + // File update + if (Pstream::master()) + { + fileName historyDir; + + word startTimeName = + time_.timeName(mesh.time().startTime().value()); + + if (Pstream::parRun()) + { + // Put in undecomposed case (Note: gives problems for + // distributed data running) + historyDir = time_.path()/".."/"history"/startTimeName; + } + else + { + historyDir = time_.path()/"history"/startTimeName; + } + + // Create directory if does not exist. + mkDir(historyDir); + + + // Open new file at start up + + // OStringStream FileName; + // FileName() << "point_" << historyPointID_ << ".dat"; + + historyFilePtr_.reset + ( + new OFstream(historyDir/fileName_) + ); + + // Add headers to output data + if (historyFilePtr_.valid()) + { + historyFilePtr_() + << "# Time" << tab << "X" << tab + << "Y" << tab << "Z"; + + historyFilePtr_() << endl; + } + } + } + + // Write start time data + writeData(); +} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +bool Foam::pointHistory::execute() +{ + return writeData(); +} + + +bool Foam::pointHistory::read(const dictionary& dict) +{ + dict.readIfPresent("region", regionName_); + + return true; +} + + +// ************************************************************************* // diff --git a/src/dynamicFvMesh/interfaceTrackingFvMesh/functionObjects/pointHistory/pointHistory.H b/src/dynamicFvMesh/interfaceTrackingFvMesh/functionObjects/pointHistory/pointHistory.H new file mode 100644 index 0000000000..9eef7d7489 --- /dev/null +++ b/src/dynamicFvMesh/interfaceTrackingFvMesh/functionObjects/pointHistory/pointHistory.H @@ -0,0 +1,135 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2019 Zeljko Tukovic, FSB Zagreb. +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +Class + Foam::pointHistory + +Description + +SourceFiles + pointHistory.C + +\*---------------------------------------------------------------------------*/ + +#ifndef pointHistory_H +#define pointHistory_H + +#include "functionObject.H" +#include "dictionary.H" +#include "fvMesh.H" +#include "OFstream.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class pointHistory Declaration +\*---------------------------------------------------------------------------*/ + +class pointHistory +: + public functionObject +{ + // Private Data + + //- Name + const word name_; + + //- Reference to main object registry + const Time& time_; + + //- Region name + word regionName_; + + //- History point ID + label historyPointID_; + + //- History point + vector refHistoryPoint_; + + //- Processor of history point + label processor_; + + //- Output file name + word fileName_; + + //- Output file stream + autoPtr historyFilePtr_; + + // Private Member Functions + + //- Write data + bool writeData(); + + //- No copy construct + pointHistory(const pointHistory&) = delete; + + //- No copy assignment + void operator=(const pointHistory&) = delete; + +public: + + //- Runtime type information + TypeName("pointHistory"); + + + // Constructors + + //- Construct from components + pointHistory + ( + const word& name, + const Time& runTime, + const dictionary& + ); + + + // Member Functions + + //- execute is called at each ++ or += of the time-loop + virtual bool execute(); + + //- Read and set the function object if its data has changed + virtual bool read(const dictionary& dict); + + //- No-op + virtual bool write() + { + return true; + } +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/dynamicFvMesh/interfaceTrackingFvMesh/functionObjects/writeFreeSurface/writeFreeSurface.C b/src/dynamicFvMesh/interfaceTrackingFvMesh/functionObjects/writeFreeSurface/writeFreeSurface.C new file mode 100644 index 0000000000..2ca76be4fa --- /dev/null +++ b/src/dynamicFvMesh/interfaceTrackingFvMesh/functionObjects/writeFreeSurface/writeFreeSurface.C @@ -0,0 +1,114 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2019 Zeljko Tukovic, FSB Zagreb. +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +\*----------------------------------------------------------------------------*/ + +#include "writeFreeSurface.H" +#include "addToRunTimeSelectionTable.H" +#include "IOmanip.H" +#include "interfaceTrackingFvMesh.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ + defineTypeNameAndDebug(writeFreeSurface, 0); + + addToRunTimeSelectionTable + ( + functionObject, + writeFreeSurface, + dictionary + ); +} + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +bool Foam::writeFreeSurface::writeData() +{ + if (time_.outputTime()) + { + const fvMesh& mesh = + time_.lookupObject(polyMesh::defaultRegion); + + interfaceTrackingFvMesh& itm = + refCast + ( + const_cast + ( + mesh.lookupObject("fvSolution") + ) + ); + + itm.writeVTKControlPoints(); + } + + return true; +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::writeFreeSurface::writeFreeSurface +( + const word& name, + const Time& runTime, + const dictionary& dict +) +: + functionObject(name), + name_(name), + time_(runTime), + regionName_(polyMesh::defaultRegion) +{ + Info<< "Creating " << this->name() << " function object." << endl; + + dict.readIfPresent("region", regionName_); +} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +bool Foam::writeFreeSurface::start() +{ + return writeData(); +} + + +bool Foam::writeFreeSurface::execute() +{ + return writeData(); +} + + +bool Foam::writeFreeSurface::read(const dictionary& dict) +{ + dict.readIfPresent("region", regionName_); + + return true; +} + + +// ************************************************************************* // diff --git a/src/dynamicFvMesh/interfaceTrackingFvMesh/functionObjects/writeFreeSurface/writeFreeSurface.H b/src/dynamicFvMesh/interfaceTrackingFvMesh/functionObjects/writeFreeSurface/writeFreeSurface.H new file mode 100644 index 0000000000..48d6d6a0fb --- /dev/null +++ b/src/dynamicFvMesh/interfaceTrackingFvMesh/functionObjects/writeFreeSurface/writeFreeSurface.H @@ -0,0 +1,125 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2019 Zeljko Tukovic, FSB Zagreb. +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +Class + Foam::writeFreeSurface + +Description + A function object to write the control points on the free surface + +SourceFiles + writeFreeSurface.C + +\*---------------------------------------------------------------------------*/ + +#ifndef writeFreeSurface_H +#define writeFreeSurface_H + +#include "functionObject.H" +#include "dictionary.H" +#include "fvMesh.H" +#include "OFstream.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class writeFreeSurface Declaration +\*---------------------------------------------------------------------------*/ + +class writeFreeSurface +: + public functionObject +{ + // Private Data + + //- Name + const word name_; + + //- Reference to main object registry + const Time& time_; + + //- Region name + word regionName_; + + + // Private Member Functions + + //- Write data + bool writeData(); + + //- No copy construct + writeFreeSurface(const writeFreeSurface&) = delete; + + //- No copy assignment + void operator=(const writeFreeSurface&) = delete; + +public: + + //- Runtime type information + TypeName("writeFreeSurface"); + + + // Constructors + + //- Construct from components + writeFreeSurface + ( + const word& name, + const Time& runTime, + const dictionary& dict + ); + + + // Member Functions + + //- start is called at the start of the time-loop + virtual bool start(); + + //- execute is called at each ++ or += of the time-loop + virtual bool execute(); + + //- Read and set the function object if its data has changed + virtual bool read(const dictionary& dict); + + //- No-op + virtual bool write() + { + return true; + } +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/dynamicFvMesh/interfaceTrackingFvMesh/fvPatchFields/freeSurfacePressure/freeSurfacePressureFvPatchScalarField.C b/src/dynamicFvMesh/interfaceTrackingFvMesh/fvPatchFields/freeSurfacePressure/freeSurfacePressureFvPatchScalarField.C new file mode 100644 index 0000000000..ede6ca8fd2 --- /dev/null +++ b/src/dynamicFvMesh/interfaceTrackingFvMesh/fvPatchFields/freeSurfacePressure/freeSurfacePressureFvPatchScalarField.C @@ -0,0 +1,188 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2013-2016 OpenFOAM Foundation + Copyright (C) 2018 OpenCFD Ltd. + Copyright (C) 2019 Zeljko Tukovic, FSB Zagreb. +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +\*---------------------------------------------------------------------------*/ + +#include "freeSurfacePressureFvPatchScalarField.H" +#include "addToRunTimeSelectionTable.H" +#include "fvPatchFieldMapper.H" +#include "volFields.H" +#include "gravityMeshObject.H" +#include "turbulentTransportModel.H" +#include "interfaceTrackingFvMesh.H" +#include "singlePhaseTransportModel.H" + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::freeSurfacePressureFvPatchScalarField:: +freeSurfacePressureFvPatchScalarField +( + const fvPatch& p, + const DimensionedField& iF +) +: + fixedValueFvPatchScalarField(p, iF), + pa_(p.size(), Zero) +{} + + +Foam::freeSurfacePressureFvPatchScalarField:: +freeSurfacePressureFvPatchScalarField +( + const fvPatch& p, + const DimensionedField& iF, + const dictionary& dict +) +: + fixedValueFvPatchScalarField(p, iF, dict, false), + pa_("pa", dict, p.size()) +{ + if (dict.found("value")) + { + fvPatchScalarField::operator= + ( + scalarField("value", dict, p.size()) + ); + } + else + { + fvPatchField::operator=(pa_); + } +} + + +Foam::freeSurfacePressureFvPatchScalarField:: +freeSurfacePressureFvPatchScalarField +( + const freeSurfacePressureFvPatchScalarField& ptf, + const fvPatch& p, + const DimensionedField& iF, + const fvPatchFieldMapper& mapper +) +: + fixedValueFvPatchScalarField(ptf, p, iF, mapper), + pa_(ptf.pa_, mapper) +{} + + +Foam::freeSurfacePressureFvPatchScalarField:: +freeSurfacePressureFvPatchScalarField +( + const freeSurfacePressureFvPatchScalarField& ptf +) +: + fixedValueFvPatchScalarField(ptf), + pa_(ptf.pa_) +{} + + +Foam::freeSurfacePressureFvPatchScalarField:: +freeSurfacePressureFvPatchScalarField +( + const freeSurfacePressureFvPatchScalarField& ptf, + const DimensionedField& iF +) +: + fixedValueFvPatchScalarField(ptf, iF), + pa_(ptf.pa_) +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +void Foam::freeSurfacePressureFvPatchScalarField::autoMap +( + const fvPatchFieldMapper& m +) +{ + fixedValueFvPatchScalarField::autoMap(m); + pa_.autoMap(m); +} + + +void Foam::freeSurfacePressureFvPatchScalarField::rmap +( + const fvPatchScalarField& ptf, + const labelList& addr +) +{ + fixedValueFvPatchScalarField::rmap(ptf, addr); + + const freeSurfacePressureFvPatchScalarField& tiptf = + refCast(ptf); + + pa_.rmap(tiptf.pa_, addr); +} + + +void Foam::freeSurfacePressureFvPatchScalarField::updateCoeffs() +{ + if (updated()) + { + return; + } + + const fvMesh& mesh = patch().boundaryMesh().mesh(); + + interfaceTrackingFvMesh& itm = + refCast + ( + const_cast + ( + mesh.lookupObject("fvSolution") + ) + ); + + operator== + ( + pa_ + itm.freeSurfacePressureJump() + ); + + fixedValueFvPatchScalarField::updateCoeffs(); +} + + +void Foam::freeSurfacePressureFvPatchScalarField::write(Ostream& os) const +{ + fvPatchScalarField::write(os); + pa_.writeEntry("pa", os); + writeEntry("value", os); +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + makePatchTypeField + ( + fvPatchScalarField, + freeSurfacePressureFvPatchScalarField + ); +} + +// ************************************************************************* // diff --git a/src/dynamicFvMesh/interfaceTrackingFvMesh/fvPatchFields/freeSurfacePressure/freeSurfacePressureFvPatchScalarField.H b/src/dynamicFvMesh/interfaceTrackingFvMesh/fvPatchFields/freeSurfacePressure/freeSurfacePressureFvPatchScalarField.H new file mode 100644 index 0000000000..022d834133 --- /dev/null +++ b/src/dynamicFvMesh/interfaceTrackingFvMesh/fvPatchFields/freeSurfacePressure/freeSurfacePressureFvPatchScalarField.H @@ -0,0 +1,217 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2013-2016 OpenFOAM Foundation + Copyright (C) 2019 Zeljko Tukovic, FSB Zagreb. +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +Class + Foam::freeSurfacePressureFvPatchScalarField + +Group + grpGenericBoundaryConditions + +Description + This boundary condition provides static pressure condition for p_rgh, + calculated as: + + \f[ + p = pa - \vec{g} & \vec{r} + \f] + + where + \vartable + p | Free surface modified pressure + pa | Free surface ambient pressure + g | acceleration due to gravity [m/s^2] + \endtable + +Usage + \table + Property | Description | Required | Default value + pa | static ambient pressure | yes | 0 + \endtable + + Example of the boundary condition specification: + \verbatim + + { + type freeSurfacePressure; + pa uniform 0; + value uniform 0; // optional initial value + } + \endverbatim + +See also + Foam::fixedValueFvPatchScalarField + +SourceFiles + freeSurfacePressureFvPatchScalarField.C + +\*---------------------------------------------------------------------------*/ + +#ifndef freeSurfacePressureFvPatchScalarField_H +#define freeSurfacePressureFvPatchScalarField_H + +#include "fixedValueFvPatchFields.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class freeSurfacePressureFvPatchScalarField Declaration +\*---------------------------------------------------------------------------*/ + +class freeSurfacePressureFvPatchScalarField +: + public fixedValueFvPatchScalarField +{ +protected: + + // Protected data + + //- Ambient pressure + scalarField pa_; + +public: + + //- Runtime type information + TypeName("freeSurfacePressure"); + + + // Constructors + + //- Construct from patch and internal field + freeSurfacePressureFvPatchScalarField + ( + const fvPatch&, + const DimensionedField& + ); + + //- Construct from patch, internal field and dictionary + freeSurfacePressureFvPatchScalarField + ( + const fvPatch&, + const DimensionedField&, + const dictionary& + ); + + //- Construct by mapping given + // freeSurfacePressureFvPatchScalarField onto a new patch + freeSurfacePressureFvPatchScalarField + ( + const freeSurfacePressureFvPatchScalarField&, + const fvPatch&, + const DimensionedField&, + const fvPatchFieldMapper& + ); + + //- Construct as copy + freeSurfacePressureFvPatchScalarField + ( + const freeSurfacePressureFvPatchScalarField& + ); + + //- Construct and return a clone + virtual tmp clone() const + { + return tmp + ( + new freeSurfacePressureFvPatchScalarField(*this) + ); + } + + //- Construct as copy setting internal field reference + freeSurfacePressureFvPatchScalarField + ( + const freeSurfacePressureFvPatchScalarField&, + const DimensionedField& + ); + + //- Construct and return a clone setting internal field reference + virtual tmp clone + ( + const DimensionedField& iF + ) const + { + return tmp + ( + new freeSurfacePressureFvPatchScalarField(*this, iF) + ); + } + + + // Member functions + + // Access + + //- Return the ambient pressure + const scalarField& pa() const + { + return pa_; + } + + //- Return reference to the ambient pressure to allow adjustment + scalarField& pa() + { + return pa_; + } + + + // Mapping functions + + //- Map (and resize as needed) from self given a mapping object + virtual void autoMap + ( + const fvPatchFieldMapper& + ); + + //- Reverse map the given fvPatchField onto this fvPatchField + virtual void rmap + ( + const fvPatchScalarField&, + const labelList& + ); + + + // Evaluation functions + + //- Update the coefficients associated with the patch field + virtual void updateCoeffs(); + + + //- Write + virtual void write(Ostream&) const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/dynamicFvMesh/interfaceTrackingFvMesh/fvPatchFields/freeSurfaceVelocity/freeSurfaceVelocityFvPatchVectorField.C b/src/dynamicFvMesh/interfaceTrackingFvMesh/fvPatchFields/freeSurfaceVelocity/freeSurfaceVelocityFvPatchVectorField.C new file mode 100644 index 0000000000..06c742d7b4 --- /dev/null +++ b/src/dynamicFvMesh/interfaceTrackingFvMesh/fvPatchFields/freeSurfaceVelocity/freeSurfaceVelocityFvPatchVectorField.C @@ -0,0 +1,131 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2019 Zeljko Tukovic, FSB Zagreb. +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +\*---------------------------------------------------------------------------*/ + +#include "freeSurfaceVelocityFvPatchVectorField.H" +#include "addToRunTimeSelectionTable.H" +#include "fvPatchFieldMapper.H" +#include "volFields.H" +#include "surfaceFields.H" +#include "interfaceTrackingFvMesh.H" + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::freeSurfaceVelocityFvPatchVectorField:: +freeSurfaceVelocityFvPatchVectorField +( + const fvPatch& p, + const DimensionedField& iF +) +: + fixedGradientFvPatchVectorField(p, iF) +{} + + +Foam::freeSurfaceVelocityFvPatchVectorField:: +freeSurfaceVelocityFvPatchVectorField +( + const freeSurfaceVelocityFvPatchVectorField& ptf, + const fvPatch& p, + const DimensionedField& iF, + const fvPatchFieldMapper& mapper +) +: + fixedGradientFvPatchVectorField(ptf, p, iF, mapper) +{} + + +Foam::freeSurfaceVelocityFvPatchVectorField:: +freeSurfaceVelocityFvPatchVectorField +( + const fvPatch& p, + const DimensionedField& iF, + const dictionary& dict +) +: + fixedGradientFvPatchVectorField(p, iF) +{ + patchType() = dict.lookupOrDefault("patchType", word::null); + fvPatchVectorField::operator=(patchInternalField()); +} + + +Foam::freeSurfaceVelocityFvPatchVectorField:: +freeSurfaceVelocityFvPatchVectorField +( + const freeSurfaceVelocityFvPatchVectorField& fcvpvf, + const DimensionedField& iF +) +: + fixedGradientFvPatchVectorField(fcvpvf, iF) +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +void Foam::freeSurfaceVelocityFvPatchVectorField::updateCoeffs() +{ + if (updated()) + { + return; + } + + const fvMesh& mesh = patch().boundaryMesh().mesh(); + + interfaceTrackingFvMesh& itm = + refCast + ( + const_cast + ( + mesh.lookupObject("fvSolution") + ) + ); + + gradient() = itm.freeSurfaceSnGradU(); + + fixedGradientFvPatchVectorField::updateCoeffs(); +} + + +void Foam::freeSurfaceVelocityFvPatchVectorField::write(Ostream& os) const +{ + fixedGradientFvPatchVectorField::write(os); + writeEntry("value", os); +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + makePatchTypeField + ( + fvPatchVectorField, + freeSurfaceVelocityFvPatchVectorField + ); +} + +// ************************************************************************* // diff --git a/src/dynamicFvMesh/interfaceTrackingFvMesh/fvPatchFields/freeSurfaceVelocity/freeSurfaceVelocityFvPatchVectorField.H b/src/dynamicFvMesh/interfaceTrackingFvMesh/fvPatchFields/freeSurfaceVelocity/freeSurfaceVelocityFvPatchVectorField.H new file mode 100644 index 0000000000..77c81522de --- /dev/null +++ b/src/dynamicFvMesh/interfaceTrackingFvMesh/fvPatchFields/freeSurfaceVelocity/freeSurfaceVelocityFvPatchVectorField.H @@ -0,0 +1,140 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2019 Zeljko Tukovic, FSB Zagreb. +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +Class + Foam::freeSurfaceVelocityFvPatchVectorField + +Group + grpOutletBoundaryConditions + +Description + This boundary condition provides a velocity outlet boundary condition for + free surface patches. + +SourceFiles + freeSurfaceVelocityFvPatchVectorField.C + +\*---------------------------------------------------------------------------*/ + +#ifndef freeSurfaceVelocityFvPatchVectorField_H +#define freeSurfaceVelocityFvPatchVectorField_H + +#include "fvPatchFields.H" +#include "fixedGradientFvPatchFields.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class freeSurfaceVelocityFvPatchVectorField Declaration +\*---------------------------------------------------------------------------*/ + +class freeSurfaceVelocityFvPatchVectorField +: + public fixedGradientFvPatchVectorField +{ +public: + + //- Runtime type information + TypeName("freeSurfaceVelocity"); + + + // Constructors + + //- Construct from patch and internal field + freeSurfaceVelocityFvPatchVectorField + ( + const fvPatch&, + const DimensionedField& + ); + + //- Construct from patch, internal field and dictionary + freeSurfaceVelocityFvPatchVectorField + ( + const fvPatch&, + const DimensionedField&, + const dictionary& + ); + + //- Construct by mapping given freeSurfaceVelocityFvPatchVectorField + // onto a new patch + freeSurfaceVelocityFvPatchVectorField + ( + const freeSurfaceVelocityFvPatchVectorField&, + const fvPatch&, + const DimensionedField&, + const fvPatchFieldMapper& + ); + + //- Construct and return a clone + virtual tmp clone() const + { + return tmp + ( + new freeSurfaceVelocityFvPatchVectorField(*this) + ); + } + + //- Construct as copy setting internal field reference + freeSurfaceVelocityFvPatchVectorField + ( + const freeSurfaceVelocityFvPatchVectorField&, + const DimensionedField& + ); + + //- Construct and return a clone setting internal field reference + virtual tmp clone + ( + const DimensionedField& iF + ) const + { + return tmp + ( + new freeSurfaceVelocityFvPatchVectorField(*this, iF) + ); + } + + + // Member functions + + //- Update the coefficients associated with the patch field + virtual void updateCoeffs(); + + //- Write + virtual void write(Ostream&) const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/dynamicFvMesh/interfaceTrackingFvMesh/interfaceTrackingFvMesh.C b/src/dynamicFvMesh/interfaceTrackingFvMesh/interfaceTrackingFvMesh.C new file mode 100644 index 0000000000..d588bc3060 --- /dev/null +++ b/src/dynamicFvMesh/interfaceTrackingFvMesh/interfaceTrackingFvMesh.C @@ -0,0 +1,2372 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2019 Zeljko Tukovic, FSB Zagreb. +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +\*---------------------------------------------------------------------------*/ + +#include "interfaceTrackingFvMesh.H" +#include "addToRunTimeSelectionTable.H" +#include "motionSolver.H" +#include "volFields.H" +#include "wedgeFaPatch.H" +#include "wedgeFaPatchFields.H" +#include "slipFaPatchFields.H" +#include "fixedValueFaPatchFields.H" +#include "slipFvPatchFields.H" +#include "symmetryFvPatchFields.H" +#include "wallFvPatch.H" +#include "polyPatchID.H" +#include "fvcMeshPhi.H" +#include "velocityLaplacianFvMotionSolver.H" +#include "EulerDdtScheme.H" +#include "CrankNicolsonDdtScheme.H" +#include "backwardDdtScheme.H" +#include "twoDPointCorrector.H" +#include "gravityMeshObject.H" +#include "turbulentTransportModel.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ + defineTypeNameAndDebug(interfaceTrackingFvMesh, 0); + addToRunTimeSelectionTable + ( + dynamicFvMesh, + interfaceTrackingFvMesh, + IOobject + ); +} + + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +void Foam::interfaceTrackingFvMesh::initializeData() +{ + // Set free surface patch index + { + const word fsPatchName(motion().get("fsPatchName")); + + polyPatchID patch(fsPatchName, this->boundaryMesh()); + + if (!patch.active()) + { + FatalErrorInFunction + << "Patch name " << fsPatchName << " not found." + << abort(FatalError); + } + + fsPatchIndex_ = patch.index(); + } + + // Set point normal correction for finite area mesh + { + boolList& correction = aMesh().correctPatchPointNormals(); + + for (const word& patchName : pointNormalsCorrectionPatches_) + { + label patchID = aMesh().boundary().findPatchID(patchName); + + if (patchID == -1) + { + FatalErrorInFunction + << "Patch name '" << patchName + << "' for point normals correction does not exist" + << abort(FatalError); + } + + correction[patchID] = true; + } + } + + // Read motion direction + if (!normalMotionDir_) + { + motionDir_ = normalised(motion().get("motionDir")); + } + + // Check if contact angle is defined + makeContactAngle(); + + motion().readIfPresent + ( + "nonReflectingFreeSurfacePatches", + nonReflectingFreeSurfacePatches_ + ); +} + + +void Foam::interfaceTrackingFvMesh::makeUs() const +{ + DebugInFunction + << "making free-surface velocity field" << nl; + + if (UsPtr_) + { + FatalErrorInFunction + << "free-surface velocity field already exists" + << abort(FatalError); + } + + wordList patchFieldTypes + ( + aMesh().boundary().size(), + zeroGradientFaPatchVectorField::typeName + ); + + forAll(aMesh().boundary(), patchI) + { + if + ( + aMesh().boundary()[patchI].type() + == wedgeFaPatch::typeName + ) + { + patchFieldTypes[patchI] = + wedgeFaPatchVectorField::typeName; + } + else + { + label ngbPolyPatchID = + aMesh().boundary()[patchI].ngbPolyPatchIndex(); + + if (ngbPolyPatchID != -1) + { + if + ( + mesh().boundary()[ngbPolyPatchID].type() + == wallFvPatch::typeName + ) + { + patchFieldTypes[patchI] = + slipFaPatchVectorField::typeName; + } + } + } + } + + for (const word& patchName : fixedFreeSurfacePatches_) + { + const label fixedPatchID = + aMesh().boundary().findPatchID(patchName); + + if (fixedPatchID == -1) + { + FatalErrorInFunction + << "Wrong faPatch name '" << patchName + << "' in the fixedFreeSurfacePatches list" + << " defined in the dynamicMeshDict dictionary" + << abort(FatalError); + } + + label ngbPolyPatchID = + aMesh().boundary()[fixedPatchID].ngbPolyPatchIndex(); + + if (ngbPolyPatchID != -1) + { + if + ( + mesh().boundary()[ngbPolyPatchID].type() + == wallFvPatch::typeName + ) + { + patchFieldTypes[fixedPatchID] = + fixedValueFaPatchVectorField::typeName; + } + } + } + + UsPtr_ = new areaVectorField + ( + IOobject + ( + "Us", + mesh().time().timeName(), + mesh(), + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + aMesh(), + dimensioned(dimVelocity, Zero), + patchFieldTypes + ); + + for (const word& patchName : fixedFreeSurfacePatches_) + { + const label fixedPatchID = aMesh().boundary().findPatchID(patchName); + + if (fixedPatchID == -1) + { + FatalErrorInFunction + << "Wrong faPatch name '" << patchName + << "' in the fixedFreeSurfacePatches list" + << " defined in the dynamicMeshDict dictionary" + << abort(FatalError); + } + + label ngbPolyPatchID = + aMesh().boundary()[fixedPatchID].ngbPolyPatchIndex(); + + if (ngbPolyPatchID != -1) + { + if + ( + mesh().boundary()[ngbPolyPatchID].type() + == wallFvPatch::typeName + ) + { + UsPtr_->boundaryFieldRef()[fixedPatchID] == Zero; + } + } + } +} + + +void Foam::interfaceTrackingFvMesh::makeFsNetPhi() const +{ + DebugInFunction + << "making free-surface net flux" << nl; + + if (fsNetPhiPtr_) + { + FatalErrorInFunction + << "free-surface net flux already exists" + << abort(FatalError); + } + + fsNetPhiPtr_ = new areaScalarField + ( + IOobject + ( + "fsNetPhi", + mesh().time().timeName(), + mesh(), + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + aMesh(), + dimensionedScalar(dimVelocity*dimArea, Zero) + ); +} + + +void Foam::interfaceTrackingFvMesh::makeControlPoints() +{ + DebugInFunction + << "making control points" << nl; + + if (controlPointsPtr_) + { + FatalErrorInFunction + << "control points already exists" + << abort(FatalError); + } + + IOobject controlPointsHeader + ( + "controlPoints", + mesh().time().timeName(), + mesh(), + IOobject::MUST_READ + ); + + if (controlPointsHeader.typeHeaderOk()) + { + Info<< "Reading control points" << endl; + controlPointsPtr_ = + new vectorIOField + ( + IOobject + ( + "controlPoints", + mesh().time().timeName(), + mesh(), + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ) + ); + } + else + { + Info<< "Creating new control points" << endl; + controlPointsPtr_ = + new vectorIOField + ( + IOobject + ( + "controlPoints", + mesh().time().timeName(), + mesh(), + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + aMesh().areaCentres().internalField() + ); + + initializeControlPointsPosition(); + } +} + + +void Foam::interfaceTrackingFvMesh::makeMotionPointsMask() +{ + DebugInFunction + << "making motion points mask" << nl; + + if (motionPointsMaskPtr_) + { + FatalErrorInFunction + << "motion points mask already exists" + << abort(FatalError); + } + + motionPointsMaskPtr_ = new labelList + ( + mesh().boundaryMesh()[fsPatchIndex()].nPoints(), + 1 + ); + + // Mark free surface boundary points + // that belong to processor patches + forAll(aMesh().boundary(), patchI) + { + if + ( + aMesh().boundary()[patchI].type() + == processorFaPatch::typeName + ) + { + const labelList& patchPoints = + aMesh().boundary()[patchI].pointLabels(); + + forAll(patchPoints, pointI) + { + motionPointsMask()[patchPoints[pointI]] = -1; + } + } + } + + // Mark fixed free surface boundary points + for (const word& patchName : fixedFreeSurfacePatches_) + { + const label fixedPatchID = aMesh().boundary().findPatchID(patchName); + + if (fixedPatchID == -1) + { + FatalErrorInFunction + << "Wrong faPatch name in the fixedFreeSurfacePatches list" + << " defined in the dynamicMeshDict dictionary" + << abort(FatalError); + } + + const labelList& patchPoints = + aMesh().boundary()[fixedPatchID].pointLabels(); + + forAll(patchPoints, pointI) + { + motionPointsMask()[patchPoints[pointI]] = 0; + } + } +} + + +void Foam::interfaceTrackingFvMesh::makeDirections() +{ + DebugInFunction + << "make displacement directions for points and control points" << nl; + + if (pointsDisplacementDirPtr_ || facesDisplacementDirPtr_) + { + FatalErrorInFunction + << "points, control points displacement directions already exist" + << abort(FatalError); + } + + pointsDisplacementDirPtr_ = + new vectorField + ( + mesh().boundaryMesh()[fsPatchIndex()].nPoints(), + Zero + ); + + facesDisplacementDirPtr_ = + new vectorField + ( + mesh().boundaryMesh()[fsPatchIndex()].size(), + Zero + ); + + if (!normalMotionDir()) + { + if (mag(motionDir_) < SMALL) + { + FatalErrorInFunction + << "Zero motion direction" + << abort(FatalError); + } + + facesDisplacementDir() = motionDir_; + pointsDisplacementDir() = motionDir_; + } + + updateDisplacementDirections(); +} + + +void Foam::interfaceTrackingFvMesh::makePhis() +{ + DebugInFunction + << "making free-surface flux" << nl; + + if (phisPtr_) + { + FatalErrorInFunction + << "free-surface flux already exists" + << abort(FatalError); + } + + phisPtr_ = new edgeScalarField + ( + IOobject + ( + "phis", + mesh().time().timeName(), + mesh(), + IOobject::NO_READ, + IOobject::NO_WRITE + ), + linearEdgeInterpolate(Us()) & aMesh().Le() + ); +} + + +void Foam::interfaceTrackingFvMesh::makeSurfactConc() const +{ + DebugInFunction + << "making free-surface surfactant concentration field" << nl; + + if (surfactConcPtr_) + { + FatalErrorInFunction + << "free-surface surfactant concentration field already exists" + << abort(FatalError); + } + + surfactConcPtr_ = new areaScalarField + ( + IOobject + ( + "Cs", + mesh().time().timeName + ( + mesh().time().startTime().value() + ), + // mesh().time().timeName(), + mesh(), + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + aMesh() + ); +} + + +void Foam::interfaceTrackingFvMesh::makeBulkSurfactConc() const +{ + DebugInFunction + << "making volume surfactant concentration field" << nl; + + if (bulkSurfactConcPtr_) + { + FatalErrorInFunction + << "volume surfactant concentration field already exists" + << abort(FatalError); + } + + bulkSurfactConcPtr_ = new volScalarField + ( + IOobject + ( + "C", + mesh().time().timeName + ( + mesh().time().startTime().value() + ), + // mesh().time().timeName(), + mesh(), + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh() + ); + volScalarField& bulkSurfactConc = *bulkSurfactConcPtr_; + + if (mesh().time().timeIndex()-1 == 0) + { + // Initialize uniform volume surfactant concentration + bulkSurfactConc = surfactant().bulkConc(); + bulkSurfactConc.correctBoundaryConditions(); + } +} + + +void Foam::interfaceTrackingFvMesh::makeSurfaceTension() const +{ + DebugInFunction + << "making surface tension field" << nl; + + if (surfaceTensionPtr_) + { + FatalErrorInFunction + << "surface tension field already exists" + << abort(FatalError); + } + + surfaceTensionPtr_ = new areaScalarField + ( + IOobject + ( + "surfaceTension", + mesh().time().timeName(), + mesh(), + IOobject::NO_READ, + IOobject::NO_WRITE + ), + sigma() + surfactant().dSigma(surfactantConcentration())/rho_ + ); +} + + +void Foam::interfaceTrackingFvMesh::makeSurfactant() const +{ + DebugInFunction + << "making surfactant properties" << nl; + + if (surfactantPtr_) + { + FatalErrorInFunction + << "surfactant properties already exists" + << abort(FatalError); + } + + const dictionary& surfactProp = + motion().subDict("surfactantProperties"); + + surfactantPtr_ = new surfactantProperties(surfactProp); +} + + +void Foam::interfaceTrackingFvMesh::makeContactAngle() +{ + DebugInFunction + << "making contact angle field" << nl; + + if (contactAnglePtr_) + { + FatalErrorInFunction + << "contact angle already exists" + << abort(FatalError); + } + + // Check if contactAngle is defined + IOobject contactAngleHeader + ( + "contactAngle", + mesh().time().timeName(), + mesh(), + IOobject::MUST_READ + ); + + if (contactAngleHeader.typeHeaderOk()) + { + Info<< "Reading contact angle field" << endl; + + contactAnglePtr_ = + new areaScalarField + ( + IOobject + ( + "contactAngle", + mesh().time().timeName(), + mesh(), + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + aMesh() + ); + } +} + + +void Foam::interfaceTrackingFvMesh::updateDisplacementDirections() +{ + if (normalMotionDir()) + { + // Update point displacement direction + pointsDisplacementDir() = aMesh().pointAreaNormals(); + + // Correct point displacement direction at contact line + forAll(aMesh().boundary(), patchI) + { + if (contactAnglePtr_) + { + label ngbPolyPatchID = + aMesh().boundary()[patchI].ngbPolyPatchIndex(); + + if (ngbPolyPatchID != -1) + { + if + ( + mesh().boundary()[ngbPolyPatchID].type() + == wallFvPatch::typeName + ) + { + labelList patchPoints = + aMesh().boundary()[patchI].pointLabels(); + + vectorField N + ( + aMesh().boundary()[patchI] + .ngbPolyPatchPointNormals() + ); + + forAll(patchPoints, pointI) + { + pointsDisplacementDir()[patchPoints[pointI]] -= + N[pointI] + *( + N[pointI] + & pointsDisplacementDir()[patchPoints[pointI]] + ); + + pointsDisplacementDir()[patchPoints[pointI]] /= + mag + ( + pointsDisplacementDir() + [ + patchPoints[pointI] + ] + ) + SMALL; + } + } + } + } + } + + // Update face displacement direction + facesDisplacementDir() = + aMesh().faceAreaNormals().internalField(); + + // Correction of control points postion + const vectorField& Cf = aMesh().areaCentres().internalField(); + + controlPoints() = + facesDisplacementDir() + *(facesDisplacementDir()&(controlPoints() - Cf)) + + Cf; + } +} + + +void Foam::interfaceTrackingFvMesh::initializeControlPointsPosition() +{ + { + const faceList& faces = aMesh().faces(); + const pointField& points = aMesh().points(); + + pointField displacement(pointDisplacement()); + scalarField sweptVolCorr(faces.size(), Zero); + correctPointDisplacement(sweptVolCorr, displacement); + + pointField newPoints(points + displacement); + + scalarField sweptVol(faces.size(), Zero); + + forAll(faces, faceI) + { + sweptVol[faceI] = -faces[faceI].sweptVol(points, newPoints); + } + + vectorField faceArea(faces.size(), Zero); + + forAll(faceArea, faceI) + { + faceArea[faceI] = faces[faceI].unitNormal(newPoints); + } + + scalarField deltaH = scalarField(aMesh().nFaces(), Zero); + + forAll(deltaH, faceI) + { + deltaH[faceI] = sweptVol[faceI]/ + ((faceArea[faceI] & facesDisplacementDir()[faceI]) + SMALL); + + if (mag(faceArea[faceI] & facesDisplacementDir()[faceI]) < SMALL) + { + // Info<< (faceArea[faceI] & facesDisplacementDir()[faceI]) + // << ", " << faceArea[faceI] + // << ", " << facesDisplacementDir()[faceI] << endl; + + FatalError + << "Something is wrong with specified motion direction" + << abort(FatalError); + } + } + + for (const word& patchName : fixedFreeSurfacePatches_) + { + label fixedPatchID = aMesh().boundary().findPatchID(patchName); + + if (fixedPatchID == -1) + { + FatalError + << "Wrong faPatch name in the fixedFreeSurfacePatches list" + << " defined in the freeSurfaceProperties dictionary" + << abort(FatalError); + } + + const labelList& eFaces = + aMesh().boundary()[fixedPatchID].edgeFaces(); + + forAll(eFaces, edgeI) + { + deltaH[eFaces[edgeI]] *= 2.0; + } + } + + controlPoints() += facesDisplacementDir()*deltaH; + } +} + + +void Foam::interfaceTrackingFvMesh::smoothFreeSurfaceMesh() +{ + Info<< "Smoothing free surface mesh" << endl; + + controlPoints() = aMesh().areaCentres().internalField(); + + pointField displacement(pointDisplacement()); + + const faceList& faces = aMesh().faces(); + const pointField& points = aMesh().points(); + + pointField newPoints(points + displacement); + + scalarField sweptVol(faces.size(), Zero); + forAll(faces, faceI) + { + sweptVol[faceI] = -faces[faceI].sweptVol(points, newPoints); + } + + vectorField faceArea(faces.size(), Zero); + forAll(faceArea, faceI) + { + faceArea[faceI] = faces[faceI].unitNormal(newPoints); + } + + scalarField deltaHf + ( + sweptVol/(faceArea & facesDisplacementDir()) + ); + + for (const word& patchName : fixedFreeSurfacePatches_) + { + label fixedPatchID = aMesh().boundary().findPatchID(patchName); + + if (fixedPatchID == -1) + { + FatalError + << "Wrong faPatch name fixedFreeSurfacePatches list" + << " defined in the dynamicMeshDict dictionary" + << abort(FatalError); + } + + const labelList& eFaces = + aMesh().boundary()[fixedPatchID].edgeFaces(); + + forAll(eFaces, edgeI) + { + deltaHf[eFaces[edgeI]] *= 2.0; + } + } + + controlPoints() += facesDisplacementDir()*deltaHf; + + displacement = pointDisplacement(); + + velocityMotionSolver& vMotion = + refCast + ( + const_cast(motion()) + ); + + pointVectorField& pointMotionU = vMotion.pointMotionU(); + pointMotionU.primitiveFieldRef() = Zero; + + fixedValuePointPatchVectorField& fsPatchPointMeshU = + refCast + ( + const_cast + ( + pointMotionU.boundaryField()[fsPatchIndex()] + ) + ); + + fsPatchPointMeshU == + displacement/mesh().time().deltaT().value(); + + dynamicMotionSolverFvMesh::update(); +} + + +void Foam::interfaceTrackingFvMesh::updateSurfaceFlux() +{ + Phis() = fac::interpolate(Us()) & aMesh().Le(); +} + + +void Foam::interfaceTrackingFvMesh::updateSurfactantConcentration() +{ + if (!pureFreeSurface()) + { + Info<< "Correct surfactant concentration" << endl << flush; + + updateSurfaceFlux(); + + // Crate and solve the surfactanta transport equation + faScalarMatrix CsEqn + ( + fam::ddt(surfactantConcentration()) + + fam::div(Phis(), surfactantConcentration()) + - fam::laplacian + ( + surfactant().diffusion(), + surfactantConcentration() + ) + ); + + if (surfactant().soluble()) + { + #include "solveBulkSurfactant.H" + + areaScalarField Cb + ( + IOobject + ( + "Cb", + mesh().time().timeName(), + mesh(), + IOobject::NO_READ, + IOobject::NO_WRITE + ), + aMesh(), + dimensionedScalar(dimMoles/dimVolume, Zero), + zeroGradientFaPatchScalarField::typeName + ); + + Cb.ref().field() = + bulkSurfactantConcentration().boundaryField()[fsPatchIndex()]; + Cb.correctBoundaryConditions(); + + CsEqn += + fam::Sp + ( + surfactant().adsorptionCoeff()*Cb + + surfactant().adsorptionCoeff() + *surfactant().desorptionCoeff(), + surfactantConcentration() + ) + - surfactant().adsorptionCoeff() + *Cb*surfactant().saturatedConc(); + } + + CsEqn.solve(); + + // Info<< "Correct surface tension" << endl; + + surfaceTension() = + sigma() + surfactant().dSigma(surfactantConcentration())/rho_; + + if (neg(min(surfaceTension().internalField().field()))) + { + FatalErrorInFunction + << "Surface tension is negative" + << abort(FatalError); + } + } +} + + +Foam::vector Foam::interfaceTrackingFvMesh::totalPressureForce() const +{ + const scalarField& S = aMesh().S(); + + const vectorField& n = aMesh().faceAreaNormals().internalField(); + + const scalarField& P = p().boundaryField()[fsPatchIndex()]; + + vectorField pressureForces(S*P*n); + + return gSum(pressureForces); +} + + +Foam::vector Foam::interfaceTrackingFvMesh::totalViscousForce() const +{ + const auto& turbulence = + mesh().lookupObject("turbulenceProperties"); + + scalarField nu(turbulence.nuEff(fsPatchIndex())); + + // const singlePhaseTransportModel& properties = + // mesh().lookupObject + // ( + // "transportProperties" + // ); + + // dimensionedScalar nu(properties.lookup("nu")); + + const scalarField& S = aMesh().S(); + const vectorField& n = aMesh().faceAreaNormals().internalField(); + + vectorField nGradU + ( + U().boundaryField()[fsPatchIndex()].snGrad() + ); + + vectorField viscousForces + ( + - nu*S + *( + nGradU + + (fac::grad(Us())().internalField()&n) + - (n*fac::div(Us())().internalField()) + ) + ); + + return gSum(viscousForces); +} + + +Foam::vector Foam::interfaceTrackingFvMesh::totalSurfaceTensionForce() const +{ + const scalarField& S = aMesh().S(); + + const vectorField& n = aMesh().faceAreaNormals().internalField(); + + const scalarField& K = aMesh().faceCurvatures().internalField(); + + vectorField surfTensionForces(n.size(), Zero); + + if (pureFreeSurface()) + { + surfTensionForces = + S*sigma().value() + *fac::edgeIntegrate + ( + aMesh().Le()*aMesh().edgeLengthCorrection() + )().internalField(); + } + else + { + surfTensionForces = surfaceTension().internalField().field()*K*S*n; + } + + return gSum(surfTensionForces); +} + + +Foam::scalar Foam::interfaceTrackingFvMesh::maxCourantNumber() +{ + scalar CoNum = 0; + + if (pureFreeSurface()) + { + const scalarField& dE = aMesh().lPN(); + + CoNum = gMax + ( + mesh().time().deltaT().value()/ + sqrt + ( + Foam::pow(dE, 3.0)/2.0/M_PI/(sigma().value() + SMALL) + ) + ); + } + else + { + scalarField sigmaE + ( + linearEdgeInterpolate(surfaceTension())().internalField().field() + + SMALL + ); + + const scalarField& dE = aMesh().lPN(); + + CoNum = gMax + ( + mesh().time().deltaT().value()/ + sqrt + ( + Foam::pow(dE, 3.0)/2.0/M_PI/sigmaE + ) + ); + } + + return CoNum; +} + + +void Foam::interfaceTrackingFvMesh::updateProperties() +{ + const singlePhaseTransportModel& properties = + mesh().lookupObject + ( + "transportProperties" + ); + + rho_ = dimensionedScalar("rho", properties); + + sigma0_ = dimensionedScalar("sigma", properties)/rho_; +} + + +void Foam::interfaceTrackingFvMesh::correctPointDisplacement +( + const scalarField& sweptVolCorr, + vectorField& displacement +) +{ + const labelListList& pFaces = + aMesh().patch().pointFaces(); + + const faceList& faces = + aMesh().patch().localFaces(); + + const pointField& points = + aMesh().patch().localPoints(); + + for (const word& patchName : fixedFreeSurfacePatches_) + { + label fixedPatchID = aMesh().boundary().findPatchID(patchName); + + const labelList& pLabels = + aMesh().boundary()[fixedPatchID].pointLabels(); + + const labelList& eFaces = + aMesh().boundary()[fixedPatchID].edgeFaces(); + + labelHashSet pointSet; + + forAll(eFaces, edgeI) + { + label curFace = eFaces[edgeI]; + + const labelList& curPoints = faces[curFace]; + + forAll(curPoints, pointI) + { + label curPoint = curPoints[pointI]; + label index = pLabels.find(curPoint); + + if (index == -1) + { + if (!pointSet.found(curPoint)) + { + pointSet.insert(curPoint); + } + } + } + } + + labelList corrPoints = pointSet.toc(); + + labelListList corrPointFaces(corrPoints.size()); + + forAll(corrPoints, pointI) + { + label curPoint = corrPoints[pointI]; + + labelHashSet faceSet; + + forAll(pFaces[curPoint], faceI) + { + label curFace = pFaces[curPoint][faceI]; + + label index = eFaces.find(curFace); + + if (index != -1) + { + if (!faceSet.found(curFace)) + { + faceSet.insert(curFace); + } + } + } + + corrPointFaces[pointI] = faceSet.toc(); + } + + forAll(corrPoints, pointI) + { + label curPoint = corrPoints[pointI]; + + scalar curDisp = 0; + + const labelList& curPointFaces = corrPointFaces[pointI]; + + forAll(curPointFaces, faceI) + { + const face& curFace = faces[curPointFaces[faceI]]; + + label ptInFace = curFace.which(curPoint); + label next = curFace.nextLabel(ptInFace); + label prev = curFace.prevLabel(ptInFace); + + vector a = points[next] - points[curPoint]; + vector b = points[prev] - points[curPoint]; + const vector& c = pointsDisplacementDir()[curPoint]; + + curDisp += 2*sweptVolCorr[curPointFaces[faceI]]/((a^b)&c); + } + + curDisp /= curPointFaces.size(); + + displacement[curPoint] = + curDisp*pointsDisplacementDir()[curPoint]; + } + } + + + for (const word& patchName : nonReflectingFreeSurfacePatches_) + { + label nonReflectingPatchID = + aMesh().boundary().findPatchID(patchName); + + const labelList& pLabels = + aMesh().boundary()[nonReflectingPatchID].pointLabels(); + + const labelList& eFaces = + aMesh().boundary()[nonReflectingPatchID].edgeFaces(); + + labelList corrPoints = pLabels; + + labelListList corrPointFaces(corrPoints.size()); + + forAll(corrPoints, pointI) + { + label curPoint = corrPoints[pointI]; + + labelHashSet faceSet; + + forAll(pFaces[curPoint], faceI) + { + label curFace = pFaces[curPoint][faceI]; + + label index = eFaces.find(curFace); + + if (index != -1) + { + if (!faceSet.found(curFace)) + { + faceSet.insert(curFace); + } + } + } + + corrPointFaces[pointI] = faceSet.toc(); + } + + + forAll(corrPoints, pointI) + { + label curPoint = corrPoints[pointI]; + + scalar curDisp = 0; + + const labelList& curPointFaces = corrPointFaces[pointI]; + + forAll(curPointFaces, faceI) + { + const face& curFace = faces[curPointFaces[faceI]]; + + label ptInFace = curFace.which(curPoint); + label next = curFace.nextLabel(ptInFace); + label prev = curFace.prevLabel(ptInFace); + + label p0 = -1; + label p1 = -1; + label p2 = -1; + + if (corrPoints.find(next) == -1) + { + p0 = curPoint; + p1 = next; + p2 = curFace.nextLabel(curFace.which(next)); + } + else + { + p0 = curFace.prevLabel(curFace.which(prev)); + p1 = prev; + p2 = curPoint; + } + + vector a0 = points[p1] - points[p0]; + vector b0 = points[p2] - points[p1]; + vector c0 = displacement[p1]; + + scalar V0 = mag((a0^b0)&c0)/2; + + scalar DV = sweptVolCorr[curPointFaces[faceI]] - V0; + + if (corrPoints.find(prev) != -1) + { + vector a = points[curPoint] - points[prev]; + vector b = + (points[next] + displacement[next]) + - points[curPoint]; + const vector& c = pointsDisplacementDir()[curPoint]; + + curDisp += 2*DV/((a^b)&c); + } + else + { + vector a = points[curPoint] + - (points[prev] + displacement[prev]); + vector b = points[next] - points[curPoint]; + const vector& c = pointsDisplacementDir()[curPoint]; + + curDisp += 2*DV/((a^b)&c); + } + } + + curDisp /= curPointFaces.size(); + + displacement[curPoint] = + curDisp*pointsDisplacementDir()[curPoint]; + } + } +} + + +void Foam::interfaceTrackingFvMesh::correctContactLinePointNormals() +{ + // Correct normals for contact line points + // according to specified contact angle + + vectorField& N = + const_cast + ( + aMesh().pointAreaNormals() + ); + + if (contactAnglePtr_ && correctContactLineNormals()) + { + Info<< "Correcting contact line normals" << endl; + + vectorField oldPoints(aMesh().nPoints(), Zero); + + const labelList& meshPoints = aMesh().patch().meshPoints(); + + forAll(oldPoints, ptI) + { + oldPoints[ptI] = + mesh().oldPoints()[meshPoints[ptI]]; + } + +// # include "createTangentField.H" + areaVectorField tangent + ( + IOobject + ( + "tangent", + mesh().time().timeName(), + mesh(), + IOobject::NO_READ, + IOobject::NO_WRITE + ), + aMesh(), + dimensioned(dimless, Zero) + ); + + if (Pstream::parRun()) + { + const labelListList& edgeFaces = aMesh().patch().edgeFaces(); + const labelListList& pointEdges = aMesh().patch().pointEdges(); + const labelListList& pointFaces = aMesh().patch().pointFaces(); + const edgeList& edges = aMesh().edges(); + + forAll(aMesh().boundary(), patchI) + { + if + ( + aMesh().boundary()[patchI].type() + == processorFaPatch::typeName + ) + { + const processorFaPatch& procPatch = + refCast + ( + aMesh().boundary()[patchI] + ); + + const labelList& patchPointLabels = + procPatch.pointLabels(); + + forAll(patchPointLabels, pointI) + { + label curPoint = patchPointLabels[pointI]; + + // Check if processor point is boundary point + + label patchID = -1; + label edgeID = -1; + + const labelList& curPointEdges = pointEdges[curPoint]; + + forAll(curPointEdges, edgeI) + { + label curEdge = curPointEdges[edgeI]; + + if (edgeFaces[curEdge].size() == 1) + { + forAll(aMesh().boundary(), pI) + { + const labelList& curEdges = + aMesh().boundary()[pI]; + + label index = curEdges.find(curEdge); + + if (index != -1) + { + if + ( + aMesh().boundary()[pI].type() + != processorFaPatch::typeName + ) + { + patchID = pI; + edgeID = index; + break; + } + } + } + } + } + + if (patchID != -1) + { + label curEdge = + aMesh().boundary()[patchID].start() + edgeID; + + vector t = edges[curEdge].vec(oldPoints); + t /= mag(t) + SMALL; + + const labelList& curPointFaces = + pointFaces[curPoint]; + + forAll(curPointFaces, fI) + { + tangent.ref().field()[curPointFaces[fI]] = t; + } + } + } + } + } + + tangent.correctBoundaryConditions(); + } + + forAll(aMesh().boundary(), patchI) + { + label ngbPolyPatchID = + aMesh().boundary()[patchI].ngbPolyPatchIndex(); + + if (ngbPolyPatchID != -1) + { + if + ( + mesh().boundary()[ngbPolyPatchID].type() + == wallFvPatch::typeName + ) + { + const scalar rotAngle = degToRad + ( + gAverage + ( + 90 + - contactAnglePtr_->boundaryField()[patchI] + ) + ); + + vectorField ngbN + ( + aMesh().boundary()[patchI].ngbPolyPatchPointNormals() + ); + + const labelList& patchPoints = + aMesh().boundary()[patchI].pointLabels(); + + vectorField pN(N, patchPoints); + + vectorField rotationAxis(ngbN^pN); + rotationAxis /= mag(rotationAxis) + SMALL; + + + // Calc rotation axis using edge vectors + + const edgeList& edges = aMesh().edges(); + + const labelListList& pointEdges = + aMesh().boundary()[patchI].pointEdges(); + + forAll(pointEdges, pointI) + { + vector rotAx = Zero; + + forAll(pointEdges[pointI], eI) + { + label curEdge = + aMesh().boundary()[patchI].start() + + pointEdges[pointI][eI]; + + vector e = edges[curEdge].vec(oldPoints); + + e *= (e&rotationAxis[pointI]) + /mag(e&rotationAxis[pointI]); + + e /= mag(e) + SMALL; + + rotAx += e; + } + + if (pointEdges[pointI].size() == 1) + { + label curPoint = patchPoints[pointI]; + + const labelListList& ptEdges = + aMesh().patch().pointEdges(); + const labelList& curPointEdges = + ptEdges[curPoint]; + + label procPatchID = -1; + label edgeID = -1; + + const labelListList& edgeFaces = + aMesh().patch().edgeFaces(); + + forAll(curPointEdges, edgeI) + { + label curEdge = curPointEdges[edgeI]; + + if (edgeFaces[curEdge].size() == 1) + { + forAll(aMesh().boundary(), pI) + { + const labelList& curEdges = + aMesh().boundary()[pI]; + + label index = + curEdges.find(curEdge); + + if (index != -1) + { + if + ( + aMesh().boundary()[pI].type() + == processorFaPatch::typeName + ) + { + procPatchID = pI; + edgeID = index; + break; + } + } + } + } + } + + if (procPatchID != -1) + { + vector t = + tangent.boundaryField()[procPatchID] + .patchNeighbourField()()[edgeID]; + + t *= (t&rotationAxis[pointI]) + /(mag(t&rotationAxis[pointI]) + SMALL); + + t /= mag(t) + SMALL; + + rotAx += t; + } + } + + rotationAxis[pointI] = rotAx/(mag(rotAx) + SMALL); + } + + // Rodrigues' rotation formula + ngbN = ngbN*cos(rotAngle) + + rotationAxis*(rotationAxis & ngbN)*(1 - cos(rotAngle)) + + (rotationAxis^ngbN)*sin(rotAngle); + + // Info<< aMesh().boundary()[patchI].name() << endl; + forAll(patchPoints, pointI) + { + N[patchPoints[pointI]] -= + ngbN[pointI]*(ngbN[pointI]&N[patchPoints[pointI]]); + + N[patchPoints[pointI]] /= + mag(N[patchPoints[pointI]]) + SMALL; + + // Info<< N[patchPoints[pointI]] << endl; + } + } + } + } + } +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::interfaceTrackingFvMesh::interfaceTrackingFvMesh(const IOobject& io) +: + dynamicMotionSolverFvMesh(io), + aMeshPtr_(new faMesh(*this)), + fsPatchIndex_(-1), + fixedFreeSurfacePatches_ + ( + motion().lookup("fixedFreeSurfacePatches") + ), + nonReflectingFreeSurfacePatches_(), + pointNormalsCorrectionPatches_ + ( + motion().lookup("pointNormalsCorrectionPatches") + ), + normalMotionDir_ + ( + motion().lookup("normalMotionDir") + ), + motionDir_(Zero), + smoothing_(motion().lookupOrDefault("smoothing", false)), + pureFreeSurface_(motion().lookupOrDefault("pureFreeSurface", true)), + rigidFreeSurface_(motion().lookupOrDefault("rigidFreeSurface", false)), + correctContactLineNormals_ + ( + motion().lookupOrDefault("correctContactLineNormals", false) + ), + sigma0_("zero", dimForce/dimLength/dimDensity, Zero), + rho_("one", dimDensity, 1.0), + timeIndex_(-1), + UsPtr_(nullptr), + controlPointsPtr_(nullptr), + motionPointsMaskPtr_(nullptr), + pointsDisplacementDirPtr_(nullptr), + facesDisplacementDirPtr_(nullptr), + fsNetPhiPtr_(nullptr), + phisPtr_(nullptr), + surfactConcPtr_(nullptr), + bulkSurfactConcPtr_(nullptr), + surfaceTensionPtr_(nullptr), + surfactantPtr_(nullptr), + contactAnglePtr_(nullptr) +{ + initializeData(); +} + + +Foam::interfaceTrackingFvMesh::interfaceTrackingFvMesh +( + const IOobject& io, + pointField&& points, + faceList&& faces, + labelList&& allOwner, + labelList&& allNeighbour, + const bool syncPar +) +: + dynamicMotionSolverFvMesh + ( + io, + std::move(points), + std::move(faces), + std::move(allOwner), + std::move(allNeighbour), + syncPar + ), + aMeshPtr_(new faMesh(*this)), + fsPatchIndex_(-1), + fixedFreeSurfacePatches_ + ( + motion().lookup("fixedFreeSurfacePatches") + ), + nonReflectingFreeSurfacePatches_(), + pointNormalsCorrectionPatches_ + ( + motion().lookup("pointNormalsCorrectionPatches") + ), + normalMotionDir_ + ( + motion().lookup("normalMotionDir") + ), + motionDir_(Zero), + smoothing_(motion().lookupOrDefault("smoothing", false)), + pureFreeSurface_(motion().lookupOrDefault("pureFreeSurface", true)), + sigma0_("zero", dimForce/dimLength/dimDensity, Zero), + rho_("one", dimDensity, 1.0), + timeIndex_(-1), + UsPtr_(nullptr), + controlPointsPtr_(nullptr), + motionPointsMaskPtr_(nullptr), + pointsDisplacementDirPtr_(nullptr), + facesDisplacementDirPtr_(nullptr), + fsNetPhiPtr_(nullptr), + phisPtr_(nullptr), + surfactConcPtr_(nullptr), + bulkSurfactConcPtr_(nullptr), + surfaceTensionPtr_(nullptr), + surfactantPtr_(nullptr), + contactAnglePtr_(nullptr) +{ + initializeData(); +} + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +Foam::interfaceTrackingFvMesh::~interfaceTrackingFvMesh() +{ + deleteDemandDrivenData(UsPtr_); + deleteDemandDrivenData(controlPointsPtr_); + deleteDemandDrivenData(motionPointsMaskPtr_); + deleteDemandDrivenData(pointsDisplacementDirPtr_); + deleteDemandDrivenData(facesDisplacementDirPtr_); + deleteDemandDrivenData(fsNetPhiPtr_); + deleteDemandDrivenData(phisPtr_); + deleteDemandDrivenData(surfactConcPtr_); + deleteDemandDrivenData(bulkSurfactConcPtr_); + deleteDemandDrivenData(surfaceTensionPtr_); + deleteDemandDrivenData(surfactantPtr_); + deleteDemandDrivenData(contactAnglePtr_); +} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +Foam::areaVectorField& Foam::interfaceTrackingFvMesh::Us() +{ + if (!UsPtr_) + { + makeUs(); + } + + return *UsPtr_; +} + + +const Foam::areaVectorField& Foam::interfaceTrackingFvMesh::Us() const +{ + if (!UsPtr_) + { + makeUs(); + } + + return *UsPtr_; +} + + +Foam::areaScalarField& Foam::interfaceTrackingFvMesh::fsNetPhi() +{ + if (!fsNetPhiPtr_) + { + makeFsNetPhi(); + } + + return *fsNetPhiPtr_; +} + + +const Foam::areaScalarField& Foam::interfaceTrackingFvMesh::fsNetPhi() const +{ + if (!fsNetPhiPtr_) + { + makeFsNetPhi(); + } + + return *fsNetPhiPtr_; +} + + +void Foam::interfaceTrackingFvMesh::correctUsBoundaryConditions() +{ + forAll(Us().boundaryField(), patchI) + { + if + ( + Us().boundaryField()[patchI].type() + == calculatedFaPatchVectorField::typeName + ) + { + vectorField& pUs = Us().boundaryFieldRef()[patchI]; + + pUs = Us().boundaryField()[patchI].patchInternalField(); + + label ngbPolyPatchID = + aMesh().boundary()[patchI].ngbPolyPatchIndex(); + + if (ngbPolyPatchID != -1) + { + if + ( + ( + U().boundaryField()[ngbPolyPatchID].type() + == slipFvPatchVectorField::typeName + ) + || + ( + U().boundaryField()[ngbPolyPatchID].type() + == symmetryFvPatchVectorField::typeName + ) + ) + { + vectorField N + ( + aMesh().boundary()[patchI].ngbPolyPatchFaceNormals() + ); + + pUs -= N*(N&pUs); + } + } + } + } + + Us().correctBoundaryConditions(); +} + + +void Foam::interfaceTrackingFvMesh::updateUs() +{ + // Info<< "Update Us" << endl; + + Us().ref().field() = U().boundaryField()[fsPatchIndex()]; + + // // Correct normal component of free-surface velocity + // const vectorField& nA = aMesh().faceAreaNormals().internalField(); + // vectorField UnFs = nA*phi().boundaryField()[fsPatchIndex()] + // /mesh().boundary()[fsPatchIndex()].magSf(); + // Us().ref().field() += UnFs - nA*(nA&Us().internalField()); + + correctUsBoundaryConditions(); +} + + +const Foam::volVectorField& Foam::interfaceTrackingFvMesh::U() const +{ + return *getObjectPtr("U"); +} + + +const Foam::volScalarField& Foam::interfaceTrackingFvMesh::p() const +{ + return *getObjectPtr("p"); +} + + +const Foam::surfaceScalarField& Foam::interfaceTrackingFvMesh::phi() const +{ + return *getObjectPtr("phi"); +} + + +Foam::tmp +Foam::interfaceTrackingFvMesh::freeSurfaceSnGradU() +{ + auto tSnGradU = tmp::New(aMesh().nFaces(), Zero); + auto& SnGradU = tSnGradU.ref(); + + const vectorField& nA = aMesh().faceAreaNormals().internalField(); + + areaScalarField divUs + ( + fac::div(Us()) + - aMesh().faceCurvatures()*(aMesh().faceAreaNormals()&Us()) + ); + + areaTensorField gradUs(fac::grad(Us())); + + // Remove component of gradient normal to surface (area) + const areaVectorField& n = aMesh().faceAreaNormals(); + gradUs -= n*(n & gradUs); + gradUs.correctBoundaryConditions(); + + const turbulenceModel& turbulence = + mesh().lookupObject("turbulenceProperties"); + + scalarField nu(turbulence.nuEff(fsPatchIndex())); + + vectorField tangentialSurfaceTensionForce(nA.size(), Zero); + + if (!pureFreeSurface() && max(nu) > SMALL) + { + tangentialSurfaceTensionForce = + surfaceTensionGrad()().internalField(); + } + + SnGradU = + tangentialSurfaceTensionForce/(nu + SMALL) + - nA*divUs.internalField() + - (gradUs.internalField()&nA); + + return tSnGradU; +} + + +Foam::tmp +Foam::interfaceTrackingFvMesh::freeSurfaceSnGradUn() +{ + auto tSnGradUn = tmp::New(aMesh().nFaces(), Zero); + auto& SnGradUn = tSnGradUn.ref(); + + areaScalarField divUs + ( + fac::div(Us()) + - aMesh().faceCurvatures()*(aMesh().faceAreaNormals()&Us()) + ); + + SnGradUn = -divUs.internalField(); + + return tSnGradUn; +} + + +Foam::tmp +Foam::interfaceTrackingFvMesh::freeSurfacePressureJump() +{ + auto tPressureJump = tmp::New(aMesh().nFaces(), Zero); + auto& pressureJump = tPressureJump.ref(); + + const scalarField& K = aMesh().faceCurvatures().internalField(); + + const uniformDimensionedVectorField& g = + meshObjects::gravity::New(mesh().time()); + + const turbulenceModel& turbulence = + mesh().lookupObject("turbulenceProperties"); + + scalarField nu(turbulence.nuEff(fsPatchIndex())); + + pressureJump = + - (g.value() & mesh().Cf().boundaryField()[fsPatchIndex()]) + + 2.0*nu*freeSurfaceSnGradUn(); + + if (pureFreeSurface()) + { + pressureJump -= sigma().value()*K; + } + else + { + pressureJump -= surfaceTension().internalField()*K; + } + + return tPressureJump; +} + + +Foam::vectorField& Foam::interfaceTrackingFvMesh::controlPoints() +{ + if (!controlPointsPtr_) + { + makeControlPoints(); + } + + return *controlPointsPtr_; +} + + +Foam::labelList& Foam::interfaceTrackingFvMesh::motionPointsMask() +{ + if (!motionPointsMaskPtr_) + { + makeMotionPointsMask(); + } + + return *motionPointsMaskPtr_; +} + + +Foam::vectorField& Foam::interfaceTrackingFvMesh::pointsDisplacementDir() +{ + if (!pointsDisplacementDirPtr_) + { + makeDirections(); + } + + return *pointsDisplacementDirPtr_; +} + + +Foam::vectorField& Foam::interfaceTrackingFvMesh::facesDisplacementDir() +{ + if (!facesDisplacementDirPtr_) + { + makeDirections(); + } + + return *facesDisplacementDirPtr_; +} + + +Foam::edgeScalarField& Foam::interfaceTrackingFvMesh::Phis() +{ + if (!phisPtr_) + { + makePhis(); + } + + return *phisPtr_; +} + +Foam::areaScalarField& +Foam::interfaceTrackingFvMesh::surfactantConcentration() +{ + if (!surfactConcPtr_) + { + makeSurfactConc(); + } + + return *surfactConcPtr_; +} + + +const Foam::areaScalarField& +Foam::interfaceTrackingFvMesh::surfactantConcentration() const +{ + if (!surfactConcPtr_) + { + makeSurfactConc(); + } + + return *surfactConcPtr_; +} + + +Foam::volScalarField& +Foam::interfaceTrackingFvMesh::bulkSurfactantConcentration() +{ + if (!bulkSurfactConcPtr_) + { + makeBulkSurfactConc(); + } + + return *bulkSurfactConcPtr_; +} + + +const Foam::volScalarField& +Foam::interfaceTrackingFvMesh::bulkSurfactantConcentration() const +{ + if (!bulkSurfactConcPtr_) + { + makeBulkSurfactConc(); + } + + return *bulkSurfactConcPtr_; +} + + +Foam::areaScalarField& +Foam::interfaceTrackingFvMesh::surfaceTension() +{ + if (!surfaceTensionPtr_) + { + makeSurfaceTension(); + } + + return *surfaceTensionPtr_; +} + + +const Foam::areaScalarField& +Foam::interfaceTrackingFvMesh::surfaceTension() const +{ + if (!surfaceTensionPtr_) + { + makeSurfaceTension(); + } + + return *surfaceTensionPtr_; +} + + +const Foam::surfactantProperties& +Foam::interfaceTrackingFvMesh::surfactant() const +{ + if (!surfactantPtr_) + { + makeSurfactant(); + } + + return *surfactantPtr_; +} + + +Foam::tmp +Foam::interfaceTrackingFvMesh::surfaceTensionGrad() +{ + auto tgrad = tmp::New + ( + IOobject + ( + "surfaceTensionGrad", + mesh().time().timeName(), + mesh(), + IOobject::NO_READ, + IOobject::NO_WRITE + ), + fac::grad(surfaceTension()) + // (-fac::grad(surfactantConcentration()/rho_)* + // surfactant().surfactR()*surfactant().surfactT()/ + // (1.0 - surfactantConcentration()/ + // surfactant().surfactSaturatedConc()))() + ); + auto& grad = tgrad.ref(); + + // Remove component of gradient normal to surface (area) + const areaVectorField& n = aMesh().faceAreaNormals(); + grad -= n*(n & grad); + grad.correctBoundaryConditions(); + + return tgrad; +} + + +bool Foam::interfaceTrackingFvMesh::update() +{ + if (timeIndex_ != mesh().time().timeIndex()) + { + if (smoothing_ && !rigidFreeSurface_) + { + smoothFreeSurfaceMesh(); + clearControlPoints(); + } + + updateDisplacementDirections(); + + updateProperties(); + + Info<< "Maximal capillary Courant number: " + << maxCourantNumber() << endl; + + const scalarField& K = aMesh().faceCurvatures().internalField(); + + Info<< "Free surface curvature: min = " << gMin(K) + << ", max = " << gMax(K) << ", average = " << gAverage(K) << nl; + + timeIndex_ = mesh().time().timeIndex(); + } + + if (!rigidFreeSurface_) + { + // This is currently relaltive flux + scalarField sweptVolCorr = + phi().boundaryField()[fsPatchIndex()]; + + // Info<< "Free surface flux: sum local = " + // << gSum(mag(sweptVolCorr)) + // << ", global = " << gSum(sweptVolCorr) << endl; + + // if (mesh().moving()) + // { + // sweptVolCorr -= + // fvc::meshPhi(U())().boundaryField()[fsPatchIndex()]; + // } + + Info<< "Free surface continuity error : sum local = " + << gSum(mag(sweptVolCorr)) << ", global = " << gSum(sweptVolCorr) + << endl; + + // For postprocessing + fsNetPhi().ref().field() = sweptVolCorr; + + word ddtScheme + ( + mesh().ddtScheme + ( + "ddt(" + U().name() + ')' + ) + ); + + if + ( + ddtScheme + == fv::CrankNicolsonDdtScheme::typeName + ) + { + sweptVolCorr *= (1.0/2.0)*mesh().time().deltaT().value(); + } + else if + ( + ddtScheme + == fv::EulerDdtScheme::typeName + ) + { + sweptVolCorr *= mesh().time().deltaT().value(); + } + else if + ( + ddtScheme + == fv::backwardDdtScheme::typeName + ) + { + if (mesh().time().timeIndex() == 1) + { + sweptVolCorr *= mesh().time().deltaT().value(); + } + else + { + sweptVolCorr *= (2.0/3.0)*mesh().time().deltaT().value(); + } + } + else + { + FatalErrorInFunction + << "Unsupported temporal differencing scheme : " + << ddtScheme << nl + << abort(FatalError); + } + + const scalarField& Sf = aMesh().S(); + const vectorField& Nf = aMesh().faceAreaNormals().internalField(); + + scalarField deltaHf + ( + sweptVolCorr/(Sf*(Nf & facesDisplacementDir())) + ); + + for (const word& patchName : fixedFreeSurfacePatches_) + { + label fixedPatchID = + aMesh().boundary().findPatchID(patchName); + + if (fixedPatchID == -1) + { + FatalErrorInFunction + << "Wrong faPatch name '" << patchName + << "'in the fixedFreeSurfacePatches list" + << " defined in dynamicMeshDict dictionary" + << abort(FatalError); + } + + const labelList& eFaces = + aMesh().boundary()[fixedPatchID].edgeFaces(); + + forAll(eFaces, edgeI) + { + deltaHf[eFaces[edgeI]] *= 2.0; + } + } + + controlPoints() += facesDisplacementDir()*deltaHf; + + pointField displacement(pointDisplacement()); + correctPointDisplacement(sweptVolCorr, displacement); + + velocityMotionSolver& vMotion = + refCast + ( + const_cast(motion()) + ); + + pointVectorField& pointMotionU = vMotion.pointMotionU(); + pointMotionU.primitiveFieldRef() = Zero; + + fixedValuePointPatchVectorField& fsPatchPointMeshU = + refCast + ( + const_cast + ( + pointMotionU.boundaryField()[fsPatchIndex()] + ) + ); + + fsPatchPointMeshU == + displacement/mesh().time().deltaT().value(); + + dynamicMotionSolverFvMesh::update(); + + correctContactLinePointNormals(); + } + else + { + vectorField displacement + ( + mesh().boundaryMesh()[fsPatchIndex()].nPoints(), + Zero + ); + + velocityMotionSolver& vMotion = + refCast + ( + const_cast(motion()) + ); + + pointVectorField& pointMotionU = vMotion.pointMotionU(); + pointMotionU.primitiveFieldRef() = Zero; + + fixedValuePointPatchVectorField& fsPatchPointMeshU = + refCast + ( + const_cast + ( + pointMotionU.boundaryField()[fsPatchIndex()] + ) + ); + + fsPatchPointMeshU == + displacement/mesh().time().deltaT().value(); + + dynamicMotionSolverFvMesh::update(); + } + + updateUs(); + + updateSurfactantConcentration(); + + return true; +} + + +void Foam::interfaceTrackingFvMesh::writeVTK() const +{ + // Write patch and points into VTK + OFstream mps(mesh().time().timePath()/"freeSurface.vtk"); + + const vectorField& points = aMesh().patch().points(); + const IndirectList& faces = aMesh().patch(); + + mps << "# vtk DataFile Version 2.0" << nl + << mesh().time().timePath()/"freeSurface.vtk" << nl + << "ASCII" << nl + << "DATASET POLYDATA" << nl + << "POINTS " << points.size() << " float" << nl; + + // Write points + List mlpBuffer(3*points.size()); + + label counter = 0; + forAll(points, i) + { + mlpBuffer[counter++] = float(points[i].x()); + mlpBuffer[counter++] = float(points[i].y()); + mlpBuffer[counter++] = float(points[i].z()); + } + + forAll(mlpBuffer, i) + { + mps << mlpBuffer[i] << ' '; + + if (i > 0 && (i % 10) == 0) + { + mps << nl; + } + } + + // Write faces + label nFaceVerts = 0; + + forAll(faces, faceI) + { + nFaceVerts += faces[faceI].size() + 1; + } + labelList mlfBuffer(nFaceVerts); + + counter = 0; + forAll(faces, faceI) + { + const face& f = faces[faceI]; + + mlfBuffer[counter++] = f.size(); + + forAll(f, fpI) + { + mlfBuffer[counter++] = f[fpI]; + } + } + mps << nl; + + mps << "POLYGONS " << faces.size() << ' ' << nFaceVerts << endl; + + forAll(mlfBuffer, i) + { + mps << mlfBuffer[i] << ' '; + + if (i > 0 && (i % 10) == 0) + { + mps << nl; + } + } + mps << nl; + + // aMesh().patch().writeVTK + // ( + // mesh().time().timePath()/"freeSurface", + // aMesh().patch(), + // aMesh().patch().points() + // ); +} + + +void Foam::interfaceTrackingFvMesh::writeVTKControlPoints() +{ + // Write control points into VTK + fileName name(mesh().time().timePath()/"freeSurfaceControlPoints.vtk"); + OFstream mps(name); + + Info<< "Writing free surface control point to " << name << endl; + + mps << "# vtk DataFile Version 2.0" << nl + << name << nl + << "ASCII" << nl + << "DATASET POLYDATA" << nl + << "POINTS " << controlPoints().size() << " float" << nl; + + forAll(controlPoints(), pointI) + { + mps << controlPoints()[pointI].x() << ' ' + << controlPoints()[pointI].y() << ' ' + << controlPoints()[pointI].z() << nl; + } + + // Write vertices + mps << "VERTICES " << controlPoints().size() << ' ' + << controlPoints().size()*2 << nl; + + forAll(controlPoints(), pointI) + { + mps << 1 << ' ' << pointI << nl; + } +} + + +// ************************************************************************* // diff --git a/src/dynamicFvMesh/interfaceTrackingFvMesh/interfaceTrackingFvMesh.H b/src/dynamicFvMesh/interfaceTrackingFvMesh/interfaceTrackingFvMesh.H new file mode 100644 index 0000000000..638ad8cce4 --- /dev/null +++ b/src/dynamicFvMesh/interfaceTrackingFvMesh/interfaceTrackingFvMesh.H @@ -0,0 +1,442 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2019 Zeljko Tukovic, FSB Zagreb. +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +Class + Foam::interfaceTrackingFvMesh + +Description + The interfaceTrackingFvMesh + +SourceFiles + interfaceTrackingFvMesh.C + +\*---------------------------------------------------------------------------*/ + +#ifndef interfaceTrackingFvMesh_H +#define interfaceTrackingFvMesh_H + +#include "dynamicMotionSolverFvMesh.H" +#include "regIOobject.H" +#include "faCFD.H" +#include "volFields.H" +#include "surfaceFields.H" +#include "surfactantProperties.H" +#include "singlePhaseTransportModel.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class interfaceTrackingFvMesh Declaration +\*---------------------------------------------------------------------------*/ + +class interfaceTrackingFvMesh +: + public dynamicMotionSolverFvMesh +{ + // Private Data + + //- Finite area mesh + autoPtr aMeshPtr_; + + //- Free surface patch index + label fsPatchIndex_; + + //- Free surface faPatch-es which do not move + wordList fixedFreeSurfacePatches_; + + //- Free surface faPatch-es where wave shuld not reflect + wordList nonReflectingFreeSurfacePatches_; + + //- Free surface patches for which point normals must be corrected + wordList pointNormalsCorrectionPatches_; + + //- Is it free-surface points displacement direction + // parallel with free-surface point normals + Switch normalMotionDir_; + + //- Free-surface points displacement direction + // if not normal motion direction + vector motionDir_; + + //- Interface smoothing at the begining of time step + Switch smoothing_; + + //- Pure free-surface + Switch pureFreeSurface_; + + //- Rigid free-surface + Switch rigidFreeSurface_; + + //- Correct contact line point normals + Switch correctContactLineNormals_; + + //- Surface tension coefficient of pure free-surface + dimensionedScalar sigma0_; + + //- Fluid density + dimensionedScalar rho_; + + //- Current interface tracking time index + label timeIndex_; + + //- Free-surface velocity field + mutable areaVectorField* UsPtr_; + + //- Points which are attached to the free-surface A side faces + // and which defines the free-surface shape + mutable vectorIOField* controlPointsPtr_; + + //- Field which additionally determines + // the motion of free-surface points + mutable labelList* motionPointsMaskPtr_; + + //- Displacement direction of free-surface points + mutable vectorField* pointsDisplacementDirPtr_; + + //- Displacement direction of free-surface control points + mutable vectorField* facesDisplacementDirPtr_; + + //- Free-surface net flux + mutable areaScalarField* fsNetPhiPtr_; + + //- Free-surface flux + mutable edgeScalarField* phisPtr_; + + //- Free-surface surfactant concetration + mutable areaScalarField* surfactConcPtr_; + + //- Volume surfactant concetration + mutable volScalarField* bulkSurfactConcPtr_; + + //- Surface tension field + mutable areaScalarField* surfaceTensionPtr_; + + //- Surfactant properties + mutable surfactantProperties* surfactantPtr_; + + //- Contact angle + mutable areaScalarField* contactAnglePtr_; + + // Private Member Functions + + // Make demand-driven data + + //- Create free surface velocity field + void makeUs() const; + + //- Create control points + void makeControlPoints(); + + //- Create motion points mask + void makeMotionPointsMask(); + + //- Create points and control point motion direction + void makeDirections(); + + //- Create free surface net flux + void makeFsNetPhi() const; + + //- Create free-surface flux + void makePhis(); + + //- Create surfactant volume concentration field + void makeBulkSurfactConc() const; + + //- Create surfactant concentration field + void makeSurfactConc() const; + + //- Create surface tension field + void makeSurfaceTension() const; + + //- Create surfactant properties + void makeSurfactant() const; + + //- Create contact angle + void makeContactAngle(); + + //- No copy construct + interfaceTrackingFvMesh(const interfaceTrackingFvMesh&) = delete; + + //- No copy assignment + void operator=(const interfaceTrackingFvMesh&) = delete; + + //- Initialize data + void initializeData(); + + //- Update control points end displacement directions + void updateDisplacementDirections(); + + //- Initialize control points position + void initializeControlPointsPosition(); + + //- Calculate free surface points displacement + // for given new control points position + tmp pointDisplacement(); + + //- Calc least sqare plane point and normal + tmp lsPlanePointAndNormal + ( + const vectorField& points, + const vector& origin, + const vector& axis + ) const; + + //- Smooth free-surface mesh + void smoothFreeSurfaceMesh(); + + //- Update free-surface flux + void updateSurfaceFlux(); + + //- Update free-surface surfactant concentration + void updateSurfactantConcentration(); + + //- Calculate total pressure force + vector totalPressureForce() const; + + //- Calculate total viscous force + vector totalViscousForce() const; + + //- Calculate total surface tension force + vector totalSurfaceTensionForce() const; + + //- Maximal surface tension based Courant number + scalar maxCourantNumber(); + + //- Update properties + void updateProperties(); + + //- Correct free surface point normals at contact line + void correctContactLinePointNormals(); + + //- Correct free surface point displacement next to fixed contact line + void correctPointDisplacement + ( + const scalarField& sweptVolCorr, + vectorField& displacement + ); + +public: + + //- Runtime type information + TypeName("interfaceTrackingFvMesh"); + + + // Constructors + + //- Construct from IOobject + interfaceTrackingFvMesh(const IOobject& io); + + //- Construct from components without boundary. + // Boundary is added using addFvPatches() member function + interfaceTrackingFvMesh + ( + const IOobject& io, + pointField&& points, + faceList&& faces, + labelList&& allOwner, + labelList&& allNeighbour, + const bool syncPar = true + ); + + + //- Destructor + ~interfaceTrackingFvMesh(); + + + // Member Functions + + fvMesh& mesh() + { + return *this; + } + + const fvMesh& mesh() const + { + return *this; + } + + //- Return reference to finite area mesh + faMesh& aMesh() + { + return aMeshPtr_(); + } + + //- Return reference to finite area mesh + const faMesh& aMesh() const + { + return aMeshPtr_(); + } + + const label& fsPatchIndex() const + { + return fsPatchIndex_; + } + + //- Pure free-surface + const Switch& pureFreeSurface() const + { + return pureFreeSurface_; + } + + //- Rigid free-surface + const Switch& rigidFreeSurface() const + { + return rigidFreeSurface_; + } + + //- Rigid free-surface + Switch& rigidFreeSurface() + { + return rigidFreeSurface_; + } + + //- Correct contact line point normals + const Switch& correctContactLineNormals() const + { + return correctContactLineNormals_; + } + + //- Correct contact line point normals + Switch& correctContactLineNormals() + { + return correctContactLineNormals_; + } + + //- Surface tension coefficient of pure free-surface + const dimensionedScalar& sigma() const + { + return sigma0_; + } + + //- Return free-surface velocity field + areaVectorField& Us(); + + //- Return free-surface velocity field + const areaVectorField& Us() const; + + //- Return free-surface net flux + const areaScalarField& fsNetPhi() const; + + //- Return free-surface net flux + areaScalarField& fsNetPhi(); + + //- Correct surface velocity boundary conditions + void correctUsBoundaryConditions(); + + //- Update free-surface velocity field + void updateUs(); + + //- Return constant reference to velocity field + const volVectorField& U() const; + + //- Return constant reference to pressure field + const volScalarField& p() const; + + //- Return constant reference to velocity field + const surfaceScalarField& phi() const; + + //- Return free surface normal derivative of velocity + tmp freeSurfaceSnGradU(); + + //- Return free surface normal derivative of normal velocity comp + tmp freeSurfaceSnGradUn(); + + //- Return free surface pressure jump + tmp freeSurfacePressureJump(); + + //- Return control points + vectorField& controlPoints(); + + //- Return reference to motion points mask field + labelList& motionPointsMask(); + + //- Return reference to point displacement direction field + vectorField& pointsDisplacementDir(); + + //- Return reference to control points displacement direction field + vectorField& facesDisplacementDir(); + + //- Motion direction swithc + bool normalMotionDir() const + { + return normalMotionDir_; + } + + //- Return free-surface fluid flux field + edgeScalarField& Phis(); + + //- Return free-surface surfactant concentration field + areaScalarField& surfactantConcentration(); + + //- Return free-surface surfactant concentration field + const areaScalarField& surfactantConcentration() const; + + //- Return volume surfactant concentration field + volScalarField& bulkSurfactantConcentration(); + + //- Return volume surfactant concentration field + const volScalarField& bulkSurfactantConcentration() const; + + //- Return surface tension field + areaScalarField& surfaceTension(); + + //- Return surface tension field + const areaScalarField& surfaceTension() const; + + //- Return surface tension gradient + tmp surfaceTensionGrad(); + + //- Return surfactant properties + const surfactantProperties& surfactant() const; + + //- Update the mesh for both mesh motion and topology change + virtual bool update(); + + //- Clear control points + void clearControlPoints() + { + deleteDemandDrivenData(controlPointsPtr_); + } + + //- Write VTK freeSurface mesh + void writeVTK() const; + + //- Write VTK freeSurface control points + void writeVTKControlPoints(); +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/dynamicFvMesh/interfaceTrackingFvMesh/solveBulkSurfactant.H b/src/dynamicFvMesh/interfaceTrackingFvMesh/solveBulkSurfactant.H new file mode 100644 index 0000000000..0fb27546a3 --- /dev/null +++ b/src/dynamicFvMesh/interfaceTrackingFvMesh/solveBulkSurfactant.H @@ -0,0 +1,46 @@ +{ + volScalarField& C = bulkSurfactantConcentration(); + + const dimensionedScalar& D = surfactant().bulkDiffusion(); + + scalar ka = surfactant().adsorptionCoeff().value(); + scalar kb = surfactant().desorptionCoeff().value(); + scalar CsInf = surfactant().saturatedConc().value(); + + const scalarField& Cs = + surfactantConcentration().internalField(); + + // const scalarField& Cfs = C.boundaryField()[fsPatchIndex()]; + + if + ( + C.boundaryField()[fsPatchIndex()].type() + == fixedGradientFvPatchScalarField::typeName + ) + { + fixedGradientFvPatchScalarField& fsC = + refCast + ( + C.boundaryFieldRef()[fsPatchIndex()] + ); + + fsC.gradient() = (ka*kb*Cs - ka*(CsInf-Cs)*fsC)/D.value(); + } + else + { + FatalErrorInFunction + << "Bulk concentration boundary condition " + << "at the free-surface patch is not " + << fixedGradientFvPatchScalarField::typeName + << exit(FatalError); + } + + fvScalarMatrix CEqn + ( + fvm::ddt(C) + + fvm::div(phi(), C, "div(phi,C)") + - fvm::laplacian(D, C, "laplacian(D,C)") + ); + + CEqn.solve(); +} diff --git a/src/dynamicFvMesh/interfaceTrackingFvMesh/surfactantProperties.H b/src/dynamicFvMesh/interfaceTrackingFvMesh/surfactantProperties.H new file mode 100644 index 0000000000..5d10771bd4 --- /dev/null +++ b/src/dynamicFvMesh/interfaceTrackingFvMesh/surfactantProperties.H @@ -0,0 +1,208 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2019 Zeljko Tukovic, FSB Zagreb. +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +Class + surfactantProperties + +Description + +SourceFiles + surfactantProperties.H + +\*---------------------------------------------------------------------------*/ + +#ifndef SurfactantProperties_H +#define SurfactantProperties_H + +#include "fvCFD.H" +#include "faCFD.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class surfactantProperties Declaration +\*---------------------------------------------------------------------------*/ + +class surfactantProperties +{ + // Private data + + //- Surfactant concentration in the bulk of fluid + dimensionedScalar bulkConc_; + + //- Saturated surfactant concentration on the free-surface + dimensionedScalar saturatedConc_; + + //- Adsorption coefficient of surfactant + dimensionedScalar adsorptionCoeff_; + + //- Desorption coefficient of surfactant + dimensionedScalar desorptionCoeff_; + + //- Diffusion coefficient of surfactant in the bulk of fluid + dimensionedScalar bulkDiffusion_; + + //- Diffusion coefficient of surfactant at the free-surface + dimensionedScalar diffusion_; + + //- Temperature of surfactant at the free-surface + dimensionedScalar T_; + + //- Universal gas constant + dimensionedScalar R_; + + //- Equilibrium surfactant concentration at the free-surface + dimensionedScalar equilibriumConc_; + + //- Is the surfactant soluble? + Switch soluble_; + + +public: + + // Constructors + + explicit surfactantProperties(const dictionary& dict) + : + bulkConc_("bulkConc", dict), + saturatedConc_("saturatedConc", dict), + adsorptionCoeff_("adsorptionCoeff", dict), + desorptionCoeff_("desorptionCoeff", dict), + bulkDiffusion_("bulkDiffusion", dict), + diffusion_("diffusion", dict), + T_("temperature", dict), + R_("R", dimGasConstant*dimMass/dimMoles, 8.3144), + equilibriumConc_ + ( + saturatedConc_ + /( + 1.0 + + desorptionCoeff_ + /bulkConc_ + ) + ), + soluble_(dict.get("soluble")) + {} + + + // Member Functions + + //- Return surfactant concentration in the bulk of fluid + const dimensionedScalar& bulkConc() const + { + return bulkConc_; + } + + //- Return saturated surfactant concentration at the free-surface + const dimensionedScalar& saturatedConc() const + { + return saturatedConc_; + } + + //- Return surfactant adsorption coefficient + const dimensionedScalar& adsorptionCoeff() const + { + return adsorptionCoeff_; + } + + //- Return surfactant desorption coefficient + const dimensionedScalar& desorptionCoeff() const + { + return desorptionCoeff_; + } + + //- Return diffusion coefficient of the surfactant in the bulk of fluid + const dimensionedScalar& bulkDiffusion() const + { + return bulkDiffusion_; + } + + //- Return diffusion coefficient of the surfactant at the free-surface + const dimensionedScalar& diffusion() const + { + return diffusion_; + } + + //- Return surfactant temeprature + const dimensionedScalar& T() const + { + return T_; + } + + //- Return universal gas constant + const dimensionedScalar& R() const + { + return R_; + } + + //- Return equilibrium surfactant concentration at the free-surface + const dimensionedScalar& equilibriumConc() const + { + return equilibriumConc_; + } + + //- Is the surfactant soluble + Switch soluble() const + { + return soluble_; + } + + //- Surface tension change due to presense of surfactants + tmp dSigma + ( + const areaScalarField& surfactConc + ) const + { + return tmp::New + ( + IOobject + ( + "dSigma", + surfactConc.mesh().mesh().time().timeName(), + surfactConc.mesh().mesh(), + IOobject::NO_READ, + IOobject::NO_WRITE + ), + ( + R()*T()*saturatedConc() + * log(1.0 - surfactConc/saturatedConc()) + ) + ); + } +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/dynamicMesh/motionSolvers/motionSolver/motionSolver.H b/src/dynamicMesh/motionSolvers/motionSolver/motionSolver.H index d594af1694..7c3ecb880c 100644 --- a/src/dynamicMesh/motionSolvers/motionSolver/motionSolver.H +++ b/src/dynamicMesh/motionSolvers/motionSolver/motionSolver.H @@ -6,7 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2011-2017 OpenFOAM Foundation - Copyright (C) 2004-2019 OpenCFD Ltd. + Copyright (C) 2016-2019 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. diff --git a/src/finiteVolume/Make/files b/src/finiteVolume/Make/files index eb8b9577e0..f313468299 100644 --- a/src/finiteVolume/Make/files +++ b/src/finiteVolume/Make/files @@ -434,6 +434,7 @@ $(snGradSchemes)/uncorrectedSnGrad/uncorrectedSnGrads.C $(snGradSchemes)/orthogonalSnGrad/orthogonalSnGrads.C $(snGradSchemes)/quadraticFitSnGrad/quadraticFitSnGrads.C $(snGradSchemes)/linearFitSnGrad/linearFitSnGrads.C +$(snGradSchemes)/skewCorrectedSnGrad/skewCorrectedSnGrads.C convectionSchemes = finiteVolume/convectionSchemes $(convectionSchemes)/convectionScheme/convectionSchemes.C diff --git a/src/finiteVolume/finiteVolume/snGradSchemes/skewCorrectedSnGrad/skewCorrectedSnGrad.C b/src/finiteVolume/finiteVolume/snGradSchemes/skewCorrectedSnGrad/skewCorrectedSnGrad.C new file mode 100644 index 0000000000..7138197f69 --- /dev/null +++ b/src/finiteVolume/finiteVolume/snGradSchemes/skewCorrectedSnGrad/skewCorrectedSnGrad.C @@ -0,0 +1,261 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2019 Zeljko Tukovic, FSB Zagreb. +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +\*---------------------------------------------------------------------------*/ + +#include "skewCorrectedSnGrad.H" +#include "volFields.H" +#include "surfaceFields.H" +#include "linear.H" +#include "fvcGrad.H" +#include "gaussGrad.H" + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +template +Foam::tmp> +Foam::fv::skewCorrectedSnGrad::fullGradCorrection +( + const GeometricField& vf +) const +{ + const fvMesh& mesh = this->mesh(); + + auto tssf = tmp>::New + ( + IOobject + ( + "snGradCorr("+vf.name()+')', + vf.instance(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh, + vf.dimensions()*mesh.nonOrthDeltaCoeffs().dimensions() + ); + auto& ssf = tssf.ref(); + + ssf.setOriented(); + ssf = dimensioned(ssf.dimensions(), Zero); + + + typedef typename + outerProduct::cmptType>::type + CmptGradType; + + const labelUList& owner = mesh.owner(); + const labelUList& neighbour = mesh.neighbour(); + + const vectorField& Sf = mesh.Sf().internalField(); + const scalarField& magSf = mesh.magSf().internalField(); + + vectorField nf(Sf/magSf); + + const vectorField& Cf = mesh.Cf().internalField(); + const vectorField& C = mesh.C().internalField(); + + const scalarField& deltaCoeffs = + mesh.deltaCoeffs().internalField(); + + surfaceVectorField kP + ( + IOobject + ( + "kP", + vf.instance(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh, + dimensionedVector(dimLength, Zero) + ); + vectorField& kPI = kP.ref().field(); + + surfaceVectorField kN + ( + IOobject + ( + "kN", + vf.instance(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh, + dimensionedVector(dimLength, Zero) + ); + vectorField& kNI = kN.ref().field(); + + kPI = Cf - vectorField(C, owner); + kPI -= Sf*(Sf&kPI)/sqr(magSf); + + kNI = Cf - vectorField(C, neighbour); + kNI -= Sf*(Sf&kNI)/sqr(magSf); + + forAll(kP.boundaryField(), patchI) + { + if (kP.boundaryField()[patchI].coupled()) + { + kP.boundaryFieldRef()[patchI] = + mesh.boundary()[patchI].Cf() + - mesh.boundary()[patchI].Cn(); + + kP.boundaryFieldRef()[patchI] -= + mesh.boundary()[patchI].Sf() + *( + mesh.boundary()[patchI].Sf() + & kP.boundaryField()[patchI] + ) + /sqr(mesh.boundary()[patchI].magSf()); + + kN.boundaryFieldRef()[patchI] = + mesh.Cf().boundaryField()[patchI] + - ( + mesh.boundary()[patchI].Cn() + + mesh.boundary()[patchI].delta() + ); + + kN.boundaryFieldRef()[patchI] -= + mesh.boundary()[patchI].Sf() + *( + mesh.boundary()[patchI].Sf() + & kN.boundaryField()[patchI] + ) + /sqr(mesh.boundary()[patchI].magSf()); + } + } + + for (direction cmpt = 0; cmpt < pTraits::nComponents; ++cmpt) + { + GeometricField cmptGrad + ( + gradScheme::cmptType>::New + ( + mesh, + mesh.gradScheme("grad(" + vf.name() + ')') + )() + .grad(vf.component(cmpt)) + ); + + const Field& cmptGradI = cmptGrad.internalField(); + + // Skewness and non-rothogonal correction + { + ssf.ref().field().replace + ( + cmpt, + ( + (kNI&Field(cmptGradI, neighbour)) + - (kPI&Field(cmptGradI, owner)) + ) + *deltaCoeffs + ); + } + + forAll(ssf.boundaryField(), patchI) + { + if (ssf.boundaryField()[patchI].coupled()) + { + ssf.boundaryFieldRef()[patchI].replace + ( + cmpt, + ( + ( + kN.boundaryField()[patchI] + & cmptGrad.boundaryField()[patchI] + .patchNeighbourField() + ) + - ( + kP.boundaryField()[patchI] + & cmptGrad.boundaryField()[patchI] + .patchInternalField() + ) + ) + *mesh.deltaCoeffs().boundaryField()[patchI] + ); + } + } + } + + // // construct GeometricField + // tmp> tssf = + // linear::type>(mesh).dotInterpolate + // ( + // mesh.nonOrthCorrectionVectors(), + // gradScheme::New + // ( + // mesh, + // mesh.gradScheme("grad(" + vf.name() + ')') + // )().grad(vf, "grad(" + vf.name() + ')') + // ); + // + // tssf.ref().rename("snGradCorr(" + vf.name() + ')'); + + return tssf; +} + + +template +Foam::tmp> +Foam::fv::skewCorrectedSnGrad::correction +( + const GeometricField& vf +) const +{ + const fvMesh& mesh = this->mesh(); + + auto tssf = tmp>::New + ( + IOobject + ( + "snGradCorr("+vf.name()+')', + vf.instance(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh, + vf.dimensions()*mesh.nonOrthDeltaCoeffs().dimensions() + ); + auto& ssf = tssf.ref(); + ssf.setOriented(); + + for (direction cmpt = 0; cmpt < pTraits::nComponents; ++cmpt) + { + ssf.replace + ( + cmpt, + skewCorrectedSnGrad::cmptType>(mesh) + .fullGradCorrection(vf.component(cmpt)) + ); + } + + return tssf; +} + + +// ************************************************************************* // diff --git a/src/finiteVolume/finiteVolume/snGradSchemes/skewCorrectedSnGrad/skewCorrectedSnGrad.H b/src/finiteVolume/finiteVolume/snGradSchemes/skewCorrectedSnGrad/skewCorrectedSnGrad.H new file mode 100644 index 0000000000..1ab2e696a9 --- /dev/null +++ b/src/finiteVolume/finiteVolume/snGradSchemes/skewCorrectedSnGrad/skewCorrectedSnGrad.H @@ -0,0 +1,162 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2019 Zeljko Tukovic, FSB Zagreb. +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +Class + Foam::fv::skewCorrectedSnGrad + +Group + grpFvSnGradSchemes + +Description + Simple central-difference snGrad scheme with non-orthogonal correction. + +SourceFiles + skewCorrectedSnGrad.C + +\*---------------------------------------------------------------------------*/ + +#ifndef skewCorrectedSnGrad_H +#define skewCorrectedSnGrad_H + +#include "snGradScheme.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace fv +{ + +/*---------------------------------------------------------------------------*\ + Class skewCorrectedSnGrad Declaration +\*---------------------------------------------------------------------------*/ + +template +class skewCorrectedSnGrad +: + public snGradScheme +{ + // Private Member Functions + + //- No copy assignment + void operator=(const skewCorrectedSnGrad&) = delete; + + +public: + + //- Runtime type information + TypeName("skewCorrected"); + + + // Constructors + + //- Construct from mesh + skewCorrectedSnGrad(const fvMesh& mesh) + : + snGradScheme(mesh) + {} + + + //- Construct from mesh and data stream + skewCorrectedSnGrad(const fvMesh& mesh, Istream&) + : + snGradScheme(mesh) + {} + + + //- Destructor + virtual ~skewCorrectedSnGrad() = default; + + + // Member Functions + + //- Return the interpolation weighting factors for the given field + virtual tmp deltaCoeffs + ( + const GeometricField& + ) const + { + return this->mesh().nonOrthDeltaCoeffs(); + } + + //- Return true if this scheme uses an explicit correction + virtual bool corrected() const + { + return true; + } + + //- Return the explicit correction to the skewCorrectedSnGrad + //- for the given field using the gradient of the field + tmp> + fullGradCorrection + ( + const GeometricField& + ) const; + + //- Return the explicit correction to the skewCorrectedSnGrad + //- for the given field using the gradients of the field components + virtual tmp> + correction(const GeometricField&) const; +}; + + +// * * * * * * * * Template Member Function Specialisations * * * * * * * * // + +template<> +tmp skewCorrectedSnGrad::correction +( + const volScalarField& vsf +) const; + + +template<> +tmp skewCorrectedSnGrad::correction +( + const volVectorField& vvf +) const; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace fv + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#ifdef NoRepository + #include "skewCorrectedSnGrad.C" +#endif + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/finiteVolume/finiteVolume/snGradSchemes/skewCorrectedSnGrad/skewCorrectedSnGrads.C b/src/finiteVolume/finiteVolume/snGradSchemes/skewCorrectedSnGrad/skewCorrectedSnGrads.C new file mode 100644 index 0000000000..cba98b1149 --- /dev/null +++ b/src/finiteVolume/finiteVolume/snGradSchemes/skewCorrectedSnGrad/skewCorrectedSnGrads.C @@ -0,0 +1,60 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2019 Zeljko Tukovic, FSB Zagreb. +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +\*---------------------------------------------------------------------------*/ + +#include "skewCorrectedSnGrad.H" +#include "fvMesh.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +makeSnGradScheme(skewCorrectedSnGrad) + + +// * * * * * * * * * * * * Template Specializations * * * * * * * * * * * * // + +template<> +Foam::tmp +Foam::fv::skewCorrectedSnGrad::correction +( + const volScalarField& vsf +) const +{ + return fullGradCorrection(vsf); +} + + +template<> +Foam::tmp +Foam::fv::skewCorrectedSnGrad::correction +( + const volVectorField& vvf +) const +{ + return fullGradCorrection(vvf); +} + + +// ************************************************************************* // diff --git a/src/fvMotionSolver/fvMotionSolvers/velocity/laplacian/velocityLaplacianFvMotionSolver.C b/src/fvMotionSolver/fvMotionSolvers/velocity/laplacian/velocityLaplacianFvMotionSolver.C index a5df5f749a..c6f72e366c 100644 --- a/src/fvMotionSolver/fvMotionSolvers/velocity/laplacian/velocityLaplacianFvMotionSolver.C +++ b/src/fvMotionSolver/fvMotionSolvers/velocity/laplacian/velocityLaplacianFvMotionSolver.C @@ -6,7 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2011-2016 OpenFOAM Foundation - Copyright (C) 2016 OpenCFD Ltd. + Copyright (C) 2016-2019 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -125,22 +125,30 @@ void Foam::velocityLaplacianFvMotionSolver::solve() fv::options& fvOptions(fv::options::New(fvMesh_)); - fvVectorMatrix UEqn + const label nNonOrthCorr ( - fvm::laplacian - ( - dimensionedScalar("viscosity", dimViscosity, 1.0) - *diffusivityPtr_->operator()(), - cellMotionU_, - "laplacian(diffusivity,cellMotionU)" - ) - == - fvOptions(cellMotionU_) + lookupOrDefault