[Hgholami]

Hi
I interest to develop overset version of solids4Foam for ESI openFOAM.
I develop some code for it, but have error in some parts.
The code as below
src/solids4FoamModels/fluidModels/overPimpleFluid/overPimpleFluid.C
overPimpleFluid.C

/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2016 OpenFOAM Foundation
Copyright (C) 2016-2018 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/>.

Application
overPimpleFluid

Description
Transient solver for incompressible flow of Newtonian fluids
on a moving mesh using the PIMPLE (merged PISO-SIMPLE) algorithm.

Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected.

\*---------------------------------------------------------------------------*/
#ifdef OPENFOAMESIORFOUNDATION

#include "overPimpleFluid.H"
#include "addToRunTimeSelectionTable.H"
#include "findRefCell.H"
//#include "fvCFD.H"
//#include "dynamicFvMesh.H"

#include "CorrectPhi.H" //check it
#include "constrainHbyA.H"
#include "constrainPressure.H"

#include "cellCellStencilObject.H"
//#include "zeroGradientFvPatchFields.H"
#include "localMin.H"
#include "interpolationCellPoint.H"
//#include "transform.H"
#include "fvMeshSubset.H"
#include "oversetFvPatchFields.H"
#include "oversetAdjustPhi.H"
//#include "globalOversetAdjustPhi.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

namespace Foam
{

namespace fluidModels
{

// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

defineTypeNameAndDebug(overPimpleFluid, 0);
addToRunTimeSelectionTable(physicsModel, overPimpleFluid, fluid);
addToRunTimeSelectionTable(fluidModel, overPimpleFluid, dictionary);

// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //

overPimpleFluid::overPimpleFluid
(
Time& runTime,
const word& region
)
:
fluidModel(typeName, runTime, region),
laminarTransport_(U(), phi()),
turbulence_
(
incompressible::turbulenceModel::New
(
U(), phi(), laminarTransport_
)
),
rho_
(
IOdictionary
(
IOobject
(
"transportProperties",
runTime.constant(),
mesh(),
IOobject::MUST_READ,
IOobject::NO_WRITE
)
).lookup("rho")
),

fvOptions_(fv::options::New(mesh())), //check

pRefCell_(0),
pRefValue_(0)
{
// osMesh(); //in fluidModel osMesh only available for fe41
setRefCell(p(), pimple().dict(), pRefCell_, pRefValue_);
mesh().setFluxRequired(p().name());

turbulence_->validate();

if (mesh().dynamic())
{
Info<< "Constructing face velocity Uf\n" << endl;

Uf_.set
(
new surfaceVectorField
(
IOobject
(
"Uf",
runTime.timeName(),
mesh(),
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
fvc::interpolate(U())
)
);
}
}

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

tmp<vectorField> overPimpleFluid::patchViscousForce(const label patchID) const
{
tmp<vectorField> tvF
(
new vectorField(mesh().boundary()[patchID].size(), vector::zero)
);

tvF.ref() =
rho_.value()
*(
mesh().boundary()[patchID].nf()
& (-turbulence_->devReff()().boundaryField()[patchID])
);

return tvF;
}


tmp<scalarField> overPimpleFluid::patchPressureForce(const label patchID) const
{
tmp<scalarField> tpF
(
new scalarField(mesh().boundary()[patchID].size(), 0)
);

tpF.ref() = rho_.value()*p().boundaryField()[patchID];

return tpF;
}

tmp<scalarField> overPimpleFluid::faceZoneMuEff //check
(
const label zoneID,
const label patchID
) const
{
scalarField pMuEff =
rho_.value()*turbulence_->nuEff()().boundaryField()[patchID];

tmp<scalarField> tMuEff
(
new scalarField(mesh().faceZones()[zoneID].size(), 0)
);
scalarField& muEff = tMuEff.ref();

const label patchStart =
mesh().boundaryMesh()[patchID].start();

forAll(pMuEff, I)
{
muEff[mesh().faceZones()[zoneID].whichFace(patchStart + I)] =
pMuEff[I];
}

// Parallel data exchange: collect pressure field on all processors
reduce(muEff, sumOp<scalarField>());

return tMuEff;
}

bool overPimpleFluid::evolve()
{
Info<< "Evolving fluid model: " << this->type() << endl;

fvMesh& mesh = fluidModel::mesh(); //check

bool meshChanged = false;
if (fluidModel::fsiMeshUpdate())
{
// The FSI interface is in charge of calling mesh.update()
meshChanged = fluidModel::fsiMeshUpdateChanged();
}
else
{
meshChanged = refCast<dynamicFvMesh>(mesh).update(); //check
reduce(meshChanged, orOp<bool>());
}

if (meshChanged)
{
const Time& runTime = fluidModel::runTime();
# include "volContinuity.H"

// #include "setCellMask.H" //defining cellMask, but not applicable here
// #include "setInterpolatedCells.H"
}

// Make the fluxes relative to the mesh motion
fvc::makeRelative(phi(), U());

// CourantNo

fluidModel::CourantNo();

// --- PIMPLE loop
while (pimple().loop())
{
tmp<fvVectorMatrix> tUEqn
(
fvm::ddt(U())
+ fvm::div(phi(), U())
+ turbulence_->divDevReff(U())
==
fvOptions_(U())
);

fvVectorMatrix& UEqn = tUEqn.ref();

UEqn.relax();

fvOptions_.constrain(UEqn);

if (pimple().momentumPredictor())
{
solve(UEqn == -fvc::grad(p()));
fvOptions_.correct(U());
}

// --- PISO loop

while (pimple().correct())
{
// p().boundaryField().updateCoeffs();

// Option 1: interpolate rAU, do not block out rAU on blocked cells
volScalarField rAU("rAU", 1.0/UEqn.A());
//oversetFvPatchField::oversetInterpolate(rAU);
mesh.interpolate(rAU);

surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
volVectorField H("H", UEqn.H());

volVectorField HbyA("HbyA", U());
HbyA = constrainHbyA(rAU*H, U(), p());
// oversetFvPatchVectorField::oversetInterpolate(U()) ;
if (pimple().nCorrPISO() <= 1)
{
tUEqn.clear();
}

surfaceScalarField phiHbyA = fvc::flux(HbyA);

// if (ddtCorr) default true
// {
surfaceScalarField faceMaskOld
(
localMin<scalar>(mesh).interpolate(cellMask.oldTim e()) //cellMask not defined
);
phiHbyA += rAUf*faceMaskOld*fvc::ddtCorr(U(), Uf_);
// }
// Adjust overset fluxes
// oversetAdjustPhi(phi(), U()); // Fringe flux adjustment
// globalOversetAdjustPhi(phi(), U(), p()); // Global flux adjustment
// WIP
if (p().needReference()) //instead of globalOversetAdjustPhi in fe41
{
fvc::makeRelative(phiHbyA, U());
adjustPhi(phiHbyA, U(), p());
fvc::makeAbsolute(phiHbyA, U());
}


if (adjustFringe)
{
fvc::makeRelative(phiHbyA, U());
oversetAdjustPhi(phiHbyA, U());
fvc::makeAbsolute(phiHbyA, U());
}


// Non-orthogonal pressure corrector loop
while (pimple().correctNonOrthogonal())
{
fvScalarMatrix pEqn
(
fvm::laplacian(rAUf, p())
== fvc::div(phiHbyA)
);

// Adjust non-orthogonal fringe fluxes if necessary
//osMesh().correctNonOrthoFluxes(pEqn, U()); //check it

pEqn.setReference(pRefCell_, pRefValue_);
pEqn.solve
(
mesh.solutionDict().solver
(
p().select(pimple().finalInnerIter())
)
);

if (pimple().finalNonOrthogonalIter())
{
phi() =phiHbyA - pEqn.flux();
}

// Perform overset interpolation (after flux reconstruction)
//oversetFvPatchScalarField::oversetInterpolate(p()) ;
}

fluidModel::ContinuityErrs();

// Explicitly relax pressure for momentum corrector
// except for last corrector
if (!pimple().finalIter())
{
p().relax();
}
volVectorField gradP(fvc::grad(p())); //check if need go to above if

U() = cellMask*(HbyA - rAU*gradP); //cellMask not defined
U().correctBoundaryConditions();
//oversetFvPatchVectorField::oversetInterpolate(U()) ;
fvOptions_.correct(U());
{
Uf_ = fvc::interpolate(U());
surfaceVectorField n(mesh.Sf()/mesh.magSf());
Uf_ += n*(phi/mesh.magSf() - (n & Uf));
}

// Make the fluxes relative to the mesh motion
fvc::makeRelative(phi(), U());
}
surfaceScalarField faceMask
(
localMin<scalar>(mesh).interpolate(cellMask) //still cellMask
);
phi() *= faceMask;

tUEqn.clear();

laminarTransport_.correct();
turbulence_->correct();
}

// Make the fluxes absolut to the mesh motion
fvc::makeAbsolute(phi(), U());

return 0;
}

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

} // End namespace fluidModels
} // End namespace Foam

#endif
//-------------

/* if (changed)
{
#include "setCellMask.H"
#include "setInterpolatedCells.H"

surfaceScalarField faceMaskOld
(
localMin<scalar>(mesh).interpolate(cellMask.oldTim e())
);

// Zero Uf on old faceMask (H-I)
Uf_ *= faceMaskOld;
// Update Uf and phi on new C-I faces
Uf_ += (1-faceMaskOld)*fvc::interpolate(U);
phi() = mesh.Sf() & Uf_;

// Zero phi on current H-I
surfaceScalarField faceMask
(
localMin<scalar>(mesh).interpolate(cellMask)
);
phi() *= faceMask;
}
*/

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

src/solids4FoamModels/fluidModels/overPimpleFluid/overPimpleFluid.H
overPimpleFluid.H

/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright held by original author
\\/ M anipulation |
-------------------------------------------------------------------------------
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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA

Class
overPimpleFluid

Description
Transient solver for incompressible, flow of Newtonian fluids
with dynamic mesh and overset mesh using the PIMPLE
(merged PISO-SIMPLE) algorithm.
Currently this is inconsistent version of the solver.


Author
Hrvoje Jasak, Wikki Ltd. All rights reserved.
Zeljko Tukovic, FSB Zagreb. All rights reserved.
Philip Cardiff, UCD. All rights reserved.

SourceFiles
overPimpleFluid.C

\*---------------------------------------------------------------------------*/

#if OPENFOAMESIORFOUNDATION

#ifndef overPimpleFluid_H
#define overPimpleFluid_H

#include "fluidModel.H"
#include "singlePhaseTransportModel.H"
#include "turbulentTransportModel.H"
#include "fvOptions.H"

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

namespace Foam
{

namespace fluidModels
{

/*---------------------------------------------------------------------------*\
Class overPimpleFluid Declaration
\*---------------------------------------------------------------------------*/

class overPimpleFluid
:
public fluidModel
{
// Private data

//- Transport model
singlePhaseTransportModel laminarTransport_;

//- Turbulence model
autoPtr<incompressible::turbulenceModel> turbulence_;

//- Density
dimensionedScalar rho_;

//- Creates and initialises the velocity field Uf if required
autoPtr<surfaceVectorField> Uf_; //new

//- Create finite volume options
fv::options& fvOptions_; //new

//- Reference pressure cell
label pRefCell_;

//- Reference pressure value
scalar pRefValue_;

// Private Member Functions

//- Disallow default bitwise copy construct
overPimpleFluid(const overPimpleFluid&);

//- Disallow default bitwise assignment
void operator=(const overPimpleFluid&);


public:

//- Runtime type information
TypeName("overPimpleFluid");

// Constructors

//- Construct from components
overPimpleFluid
(
Time& runTime,
const word& region = dynamicFvMesh::defaultRegion
);


// Destructor

virtual ~overPimpleFluid()
{}


// Member Functions

// Access

//- Patch viscous force (N/m2)
virtual tmp<vectorField> patchViscousForce
(
const label patchID
) const;

//- Patch pressure force (N/m2)
virtual tmp<scalarField> patchPressureForce
(
const label patchID
) const;

//- Face zone effective dynamic viscosity
virtual tmp<scalarField> faceZoneMuEff
(
const label zoneID,
const label patchID
) const;

// Edit

//- Evolve the fluid model
virtual bool evolve();
};


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

} // End namespace fluidModels
} // End namespace Foam

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

#endif

#endif // ESIOpenfoam

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

cellMask is not defined,

[bigphil]

Hi Hojatollah,

I suggest you start with an existing fluid model class that does compile and then add small bits of code at a time, where you compile each time. Then you will know which code exactly is causes compilation errors. If the definition of a class or function is missing then you will ned to add this, e.g. via #include header file, and/or by defining the required object.

If you report compilation errors here, then give the full compilation error.

[Hgholami]

Hi Philip
Thank you for your reply.
The pimpleOversetFluid of Fe41, uses

#include "globalOversetAdjustPhi.H"
oversetFvPatchVectorField:versetInterpolate(U()) ;

globalOversetAdjustPhi(phi(), U(), p()); // Global flux adjustment

// Adjust non-orthogonal fringe fluxes if necessary
osMesh().correctNonOrthoFluxes(pEqn, U()); //defined in fluidModel.C

oversetFvPatchScalarField:versetInterpolate(p()) ;

that not defined in ESI openfoam.
So, I pattenting from overPimpleFoam (ESI version). In overPimpleFoam that uses

mesh.interpolate(rAU);
HbyA = constrainHbyA(rAU*H, U(), p());
if (ddtCorr) default true
{
surfaceScalarField faceMaskOld
(
localMin<scalar>(mesh).interpolate(cellMask.oldTim e()) //cellMask not defined
);
phiHbyA += rAUf*faceMaskOld*fvc::ddtCorr(U(), Uf_);
}
.....

The cellMask should declare before, so I add createCellMask in header, but the below error occur.

In file included from fluidModels/overPimpleFluid/overPimpleFluid.C:58:0:
/usr/lib/openfoam/openfoam1912/src/overset/lnInclude/createCellMask.H:36:1: error: ‘Info’ does not name a type; did you mean ‘ynf’?
Info<< "Creating cellMask field to block out hole cells\n" << endl;
^~~~
ynf
/usr/lib/openfoam/openfoam1912/src/overset/lnInclude/createCellMask.H:38:1: error: ‘volScalarField’ does not name a type; did you mean ‘scalarField_H’?
volScalarField cellMask
^~~~~~~~~~~~~~
scalarField_H
In file included from /usr/lib/openfoam/openfoam1912/src/overset/lnInclude/createCellMask.H:53:0,
from fluidModels/overPimpleFluid/overPimpleFluid.C:58:
/usr/lib/openfoam/openfoam1912/src/overset/lnInclude/setCellMask.H:31:1: error: expected unqualified-id before ‘{’ token
{
^

createCellMask.H

Info<< "Creating cellMask field to block out hole cells\n" << endl;

volScalarField cellMask
(
IOobject
(
"cellMask",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh,
dimensionedScalar("one", dimless, 1.0),
zeroGradientFvPatchScalarField::typeName
);

#include "setCellMask.H"

setCellMask.H

{
const cellCellStencilObject& overlap = Stencil::New(mesh);
const labelList& cellTypes = overlap.cellTypes();

cellMask.primitiveFieldRef() = 1.0;
forAll(cellMask, cellI)
{
if (cellTypes[cellI] == cellCellStencil::HOLE)
{
cellMask[cellI] = 0.0;
}
}
cellMask.correctBoundaryConditions();
}

[bigphil]

createCellMask.H seems to be a header file to be included inside the main function of an overset solver, as opposed to a header file with function/class definitions, i.e. createCellMask.H actually creates objects, rather than just defining them. So you cannot include this at the top of a file; the errors you are getting are related to this file being included outside the Foam scope but in any case, you need to copy the code from within that file to be executed within your new fluid model class. Anything that happens before the time loop in the normal OpenFOAM solver should be called in the constructor of the fluid model class.

[Hgholami]

the problem of solver was in createCellMask.H and createInterpolatedCells.H that I extract the files and added to constructor.
Now the constructor modify and the bellow code added to it

cellMask
(
IOobject
(
"cellMask",
runTime.timeName(),
mesh(),
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh(),
dimensionedScalar("one", dimless, 1.0),
zeroGradientFvPatchScalarField::typeName
),
interpolatedCells
(
IOobject
(
"interpolatedCells",
runTime.timeName(),
mesh(),
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh(),
dimensionedScalar("one", dimless, 1.0),
zeroGradientFvPatchScalarField::typeName
),

That compile succesfully.

[bigphil]

Originally Posted by Hgholami

the problem of solver was in createCellMask.H and createInterpolatedCells.H that I extract the files and added to constructor.
Now the constructor modify and the bellow code added to it

That compile succesfully.

Glad to hear it.

If you get it working, feel free to fork solids4foam and make a pull request to add the new code; it sounds like it will be generally useful.

[Hgholami]

Dear Philip
I compile the solver successfully, but the major problem exist in dynamicMotion.
In openfoam1912, the overset patch is different from interface. So the displacement in fluidSolidInterface.H only apply to interface and overset patch is stationary. It seen, we should modify the fluidSolidInterface.H

Originally Posted by bigphil

Glad to hear it.

If you get it working, feel free to fork solids4foam and make a pull request to add the new code; it sounds like it will be generally useful.

[bigphil]

Originally Posted by Hgholami

Dear Philip
I compile the solver successfully, but the major problem exist in dynamicMotion.
In openfoam1912, the overset patch is different from interface. So the displacement in fluidSolidInterface.H only apply to interface and overset patch is stationary. It seen, we should modify the fluidSolidInterface.H

I guess it should be possible to set the overset patch mesh boundary condition as zero-gradient or similar, and this can be done in the case. Have you considered this?

[Hgholami]

Originally Posted by bigphil

I guess it should be possible to set the overset patch mesh boundary condition as zero-gradient or similar, and this can be done in the case. Have you considered this?

I consider zeroGradient, fixedValue, slip, calculated for overset patch in pointMotionU, but the patch still be stationary.
As the moveFluidMesh() class only calculate the interface displacement, I think, this class should modify or new boundary condition for overset mesh derives.

[bigphil]

moveFluidMesh essentially just sets the FSI interface mesh boundary condition and then calls mesh.update() to update the mesh (e.g. solve mesh motion equations). For overset, it should be the same as the overset patch is away from the FSI interface and not directly affect by it.

[Hgholami]

Originally Posted by bigphil

moveFluidMesh essentially just sets the FSI interface mesh boundary condition and then calls mesh.update() to update the mesh (e.g. solve mesh motion equations). For overset, it should be the same as the overset patch is away from the FSI interface and not directly affect by it.

Yes, problem is out of moeFluidMesh() class.
As you say, moveFluidMesh() first moves The FSI interface and then, the "bool meshChanged = fluidMesh().update();" calls mesh.update() and move internal nodes.
In the other hand, in openfoam1912, the movement of overset patch first calculated and then transfer to whole mesh, So available dynamicMesh method for this problem not applicable.
Base on restriction of dynamicMesh in moveFluidMesh(), the available is "fvMeshsolver". So, we need to calculate overset patch movement base on FSI interface movement or transfer velocityLaplacian calculation to overset patch boundary.
The velocityLaplacian that uses in dynamicMeshDict is

void Foam::velocityLaplacianFvMotionSolver::solve()
{
// The points have moved so before interpolation update
// the fvMotionSolver accordingly
movePoints(fvMesh_.points());

diffusivityPtr_->correct();
pointMotionU_.boundaryFieldRef().updateCoeffs();

fv::options& fvOptions(fv::options::New(fvMesh_));

const label nNonOrthCorr
(
getOrDefault<label>("nNonOrthogonalCorrectors", 1)
);

for (label i=0; i<nNonOrthCorr; ++i)
{
fvVectorMatrix UEqn
(
fvm::laplacian
(
dimensionedScalar("viscosity", dimViscosity, 1.0)
* diffusivityPtr_->operator()(),
cellMotionU_,
"laplacian(diffusivity,cellMotionU)"
)
==
fvOptions(cellMotionU_)
);

fvOptions.constrain(UEqn);
UEqn.solveSegregatedOrCoupled(UEqn.solverDict());
fvOptions.correct(cellMotionU_);
}
}

As the code shows, the boundaryField updates before solving Laplacian, so I think, it need to interpolate internal node movement to overset boundary. do you have any idea?
Thanks
Hojatollah

[Hgholami]

For example

for velocityLaplacian (fvMeshSolver) only move internal nodes and the boundary is fixed.
each boundary type that test are fixedValue, slip, zeroGradient, calculated, etc.
how we can corporate the overset boundary to Laplacian equation?