surfaceAlignedSBRStressFvMotionSolver.C

Go to the documentation of this file.

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | www.openfoam.com
     \\/     M anipulation  |
-------------------------------------------------------------------------------
    Copyright (C) 2015 OpenFOAM Foundation
    Copyright (C) 2015-2022 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.
 
    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.
 
    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.
 
    You should have received a copy of the GNU General Public License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
\*---------------------------------------------------------------------------*/
 
#include "surfaceAlignedSBRStressFvMotionSolver.H"
#include "addToRunTimeSelectionTable.H"
#include "pointIndexHit.H"
#include "processorPolyPatch.H"
#include "fvmLaplacian.H"
#include "fvcDiv.H"
#include "surfaceInterpolate.H"
#include "unitConversion.H"
#include "motionDiffusivity.H"
#include "fvcSmooth.H"
#include "pointMVCWeight.H"
#include "dimensionedSymmTensor.H"
#include "quaternion.H"
#include "fvOptions.H"
 
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
namespace Foam
{
    defineTypeNameAndDebug(surfaceAlignedSBRStressFvMotionSolver, 0);
 
    addToRunTimeSelectionTable
    (
        motionSolver,
        surfaceAlignedSBRStressFvMotionSolver,
        dictionary
    );
}
 
 
// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
Foam::surfaceAlignedSBRStressFvMotionSolver::
surfaceAlignedSBRStressFvMotionSolver
(
    const polyMesh& mesh,
    const IOdictionary& dict
)
:
    displacementSBRStressFvMotionSolver(mesh, dict),
    surfaceNames_(coeffDict().lookup("surfaces")),
    surfaceMesh_(surfaceNames_.size()),
    cellRot_
    (
        IOobject
        (
            "cellRot",
            fvMesh_.time().timeName(),
            fvMesh_,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        fvMesh_,
        dimensionedVector(dimless, Zero)
    ),
    maxAng_(coeffDict().getOrDefault<scalar>("maxAng", 80)),
    minAng_(coeffDict().getOrDefault<scalar>("minAng", 20)),
    accFactor_(coeffDict().getOrDefault<scalar>("accFactor", 1)),
    smoothFactor_(coeffDict().getOrDefault<scalar>("smoothFactor", 0.9)),
    nNonOrthogonalCorr_(coeffDict().get<label>("nNonOrthogonalCorr")),
    pointDisplacement_(pointDisplacement()),
    sigmaD_
    (
        IOobject
        (
            "sigmaD",
            fvMesh_.time().timeName(),
            fvMesh_,
            IOobject::READ_IF_PRESENT,
            IOobject::NO_WRITE
        ),
        fvMesh_,
        dimensionedSymmTensor(dimless, Zero)
    ),
    minSigmaDiff_(coeffDict().getOrDefault<scalar>("minSigmaDiff", 1e-4))
{
    forAll(surfaceNames_, i)
    {
        surfaceMesh_.set
        (
            i,
            new triSurfaceMesh
            (
                IOobject
                (
                    surfaceNames_[i],
                    mesh.time().constant(),
                    "triSurface",
                    mesh.time(),
                    IOobject::MUST_READ,
                    IOobject::NO_WRITE
                )
            )
        );
    }
}
 
 
// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
Foam::surfaceAlignedSBRStressFvMotionSolver::
~surfaceAlignedSBRStressFvMotionSolver()
{}
 
 
// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 
void Foam::surfaceAlignedSBRStressFvMotionSolver::calculateCellRot()
{
    cellRot_.primitiveFieldRef() = Zero;
    pointDisplacement_.primitiveFieldRef() = Zero;
 
    // Find intersections
    pointField start(fvMesh_.nInternalFaces());
    pointField end(start.size());
 
    const vectorField& Cc = fvMesh_.cellCentres();
    const polyBoundaryMesh& pbm = fvMesh_.boundaryMesh();
 
    for (label faceI = 0; faceI < fvMesh_.nInternalFaces(); faceI++)
    {
        start[faceI] = Cc[fvMesh_.faceOwner()[faceI]];
        end[faceI] = Cc[fvMesh_.faceNeighbour()[faceI]];
    }
 
    DynamicList<label> hitCells;
 
    forAll(surfaceMesh_, surfi)
    {
        List<pointIndexHit> hit(start.size());
        surfaceMesh_[surfi].findLineAny(start, end, hit);
 
        labelField pointsCount(fvMesh_.nPoints(), 1);
 
        const vectorField& nf = surfaceMesh_[surfi].faceNormals();
 
        const vectorField& SfMesh = fvMesh_.faceAreas();
 
        const vectorField nSfMesh(SfMesh/mag(SfMesh));
 
        DynamicList<label> cellsHit;
 
        forAll(hit, facei)
        {
            if (hit[facei].hit())
            {
                label rotCellId(-1);
                const vector& hitPoint = hit[facei].point();
 
                if (fvMesh_.isInternalFace(facei))
                {
                    const point& ownCc = Cc[fvMesh_.faceOwner()[facei]];
                    const point& nbrCc = Cc[fvMesh_.faceNeighbour()[facei]];
 
                    if (hitPoint.distSqr(ownCc) < hitPoint.distSqr(nbrCc))
                    {
                        rotCellId = fvMesh_.faceOwner()[facei];
                    }
                    else
                    {
                        rotCellId = fvMesh_.faceNeighbour()[facei];
                    }
                }
                else
                {
                    label patchi = pbm.whichPatch(facei);
                    if (isA<processorPolyPatch>(pbm[patchi]))
                    {
                        const point& ownCc = Cc[fvMesh_.faceOwner()[facei]];
 
                        const vector cCentreOne = ownCc - hitPoint;
 
                        const vector nbrCc =
                            refCast<const processorPolyPatch>(pbm[patchi])
                                .neighbFaceCellCentres()[facei];
 
                        const vector cCentreTwo = nbrCc - hitPoint;
 
                        if (cCentreOne < cCentreTwo)
                        {
                            rotCellId = fvMesh_.faceOwner()[facei];
                        }
                    }
                }
 
                // Note: faces on boundaries that get hit are not included as
                // the pointDisplacement on boundaries is usually zero for
                // this solver.
 
                // Search for closest direction on faces on mesh
                // and surface normal.
                if (rotCellId != -1)
                {
                    const labelList& cFaces = fvMesh_.cells()[rotCellId];
 
                    scalar cosMax(-GREAT);
                    label faceId(-1);
                    forAll(cFaces, k)
                    {
                        scalar tmp =
                            mag(nf[hit[facei].index()] & nSfMesh[cFaces[k]]);
 
                        if (tmp > cosMax)
                        {
                            cosMax = tmp;
                            faceId = cFaces[k];
                        }
                    }
                    if
                    (
                        faceId != -1
                    &&
                        (
                            ::cos(degToRad(minAng_)) > cosMax
                            || cosMax > ::cos(degToRad(maxAng_))
 
                        )
                    )
                    {
                        cellRot_[rotCellId] =
                            nSfMesh[faceId]^nf[hit[facei].index()];
 
                        const scalar magRot = mag(cellRot_[rotCellId]);
 
                        if (magRot > 0)
                        {
                            const scalar theta = ::asin(magRot);
                            quaternion q(cellRot_[rotCellId]/magRot, theta);
                            const tensor R = q.R();
                            const labelList& cPoints =
                                fvMesh_.cellPoints(rotCellId);
 
                            forAll(cPoints, j)
                            {
                                const label pointId = cPoints[j];
 
                                pointsCount[pointId]++;
 
                                const vector pointPos =
                                    fvMesh_.points()[pointId];
 
                                pointDisplacement_[pointId] +=
                                    (R & (pointPos - hitPoint))
                                  - (pointPos - hitPoint);
                            }
                        }
                    }
                }
            }
        }
 
        vectorField& pd = pointDisplacement_.primitiveFieldRef();
        forAll(pd, pointi)
        {
            vector& point = pd[pointi];
            point /= pointsCount[pointi];
        }
    }
}
 
 
void Foam::surfaceAlignedSBRStressFvMotionSolver::solve()
{
    // The points have moved so before interpolation update
    // the motionSolver accordingly
    this->movePoints(fvMesh_.points());
 
    volVectorField& cellDisp = cellDisplacement();
 
    diffusivity().correct();
 
    // Calculate rotations on surface intersection
    calculateCellRot();
 
    auto tUd = volVectorField::New
    (
        "Ud",
        IOobject::NO_REGISTER,
        fvMesh_,
        dimensionedVector(dimLength, Zero),
        cellMotionBoundaryTypes<vector>
        (
            pointDisplacement().boundaryField()
        )
    );
    auto& Ud = tUd.ref();
 
    const vectorList& C = fvMesh_.C();
    forAll(Ud, i)
    {
        pointMVCWeight pointInter(fvMesh_, C[i], i);
        Ud[i] = pointInter.interpolate(pointDisplacement_);
    }
 
    volScalarField Udx(Ud.component(vector::X));
    fvc::smooth(Udx, smoothFactor_);
 
    volScalarField Udy(Ud.component(vector::Y));
    fvc::smooth(Udy, smoothFactor_);
 
    volScalarField Udz(Ud.component(vector::Z));
    fvc::smooth(Udz, smoothFactor_);
 
    Ud.replace(vector::X, Udx);
    Ud.replace(vector::Y, Udy);
    Ud.replace(vector::Z, Udz);
 
    const volTensorField gradD("gradD", fvc::grad(Ud));
 
    auto tmu = volScalarField::New
    (
        "mu",
        IOobject::NO_REGISTER,
        fvMesh_,
        dimensionedScalar(dimless, Zero)
    );
    auto& mu = tmu.ref();
 
    const scalarList& V = fvMesh_.V();
    mu.primitiveFieldRef() = (1.0/V);
 
    const volScalarField lambda(-mu);
 
    const volSymmTensorField newSigmaD
    (
        mu*twoSymm(gradD) + (lambda*I)*tr(gradD)
    );
 
    const volSymmTensorField magNewSigmaD(sigmaD_ + accFactor_*newSigmaD);
 
    const scalar diffSigmaD =
        gSum(mag(sigmaD_.oldTime().primitiveField()))
     -  gSum(mag(magNewSigmaD.primitiveField()));
 
    if (mag(diffSigmaD) > minSigmaDiff_)
    {
        sigmaD_ = magNewSigmaD;
    }
 
    const dimensionedScalar oneViscosity("viscosity", dimViscosity, 1.0);
 
    const surfaceScalarField Df(oneViscosity*diffusivity().operator()());
 
    pointDisplacement_.boundaryFieldRef().updateCoeffs();
 
    fv::options& fvOptions(fv::options::New(fvMesh_));
 
    const volTensorField gradCd(fvc::grad(cellDisp));
 
    for (int nonOrth=0; nonOrth<=nNonOrthogonalCorr_; nonOrth++)
    {
        fvVectorMatrix DEqn
        (
            fvm::laplacian
            (
                2*Df*fvc::interpolate(mu),
                cellDisp,
                "laplacian(diffusivity,cellDisplacement)"
            )
          + fvc::div
            (
                Df*
                (
                    fvc::interpolate(mu)
                  * (fvMesh_.Sf() & fvc::interpolate(gradCd.T() - gradCd))
                  - fvc::interpolate(lambda)*fvMesh_.Sf()
                  * fvc::interpolate(tr(gradCd))
                )
            )
          ==
            oneViscosity*fvc::div(sigmaD_)
          + fvOptions(cellDisp)
        );
 
        fvOptions.constrain(DEqn);
 
        // Note: solve uncoupled
        DEqn.solveSegregatedOrCoupled();
 
        fvOptions.correct(cellDisp);
    }
}
 
 
// ************************************************************************* //