| GuilhermeMP |
September 27, 2012 04:40 |
potentialFreeSurfaceFoam with Dynamic Mesh
Hi all,
I'm developing a version of potentialFreeSurfaceFoam with dynamic mesh abilities to be able to solver for 6DOF body motion of floating bodies. The main point of this, instead of using interDymFoam, is to have a cheaper solver for floating body motion, without the need to solve for the interface or atmospheric conditions.
I'm using interDymFoam to base my code, since it has already dynamic mesh and it is able to solve for 6 DOF body motion in the free surface.
I created a correctPhi.H file (non existing in the standard solver) and adapted the potentialFreeSurfaceFoam.C, UEqn.H, pEqn.H files and the Make/options and Make/files files.
The code to the modified solver is below.
I the correctPHi.H file, there's a line that came from interDymFoam:
"dimensionedScalar rAUf("(1|A(U))", dimTime/rho.dimensions(), 1.0);"
This line has the 'rho' variable, but in potentialFreeSurfaceFoam, the pressure is kinematic, so there is no need for 'rho'. How should I change this line?
Since I'm not comfortable with the dynamic mesh programming, I did it in a "blind" way, so I would very much appreciate if someone could check my code and/or give me some hints on how to code the dynamic mesh.
correctPhi.H
Code:
{
#include "continuityErrs.H"
wordList pcorrTypes
(
p_gh.boundaryField().size(),
zeroGradientFvPatchScalarField::typeName
);
forAll (p_gh.boundaryField(), i)
{
if (p_gh.boundaryField()[i].fixesValue())
{
pcorrTypes[i] = fixedValueFvPatchScalarField::typeName;
}
}
volScalarField pcorr
(
IOobject
(
"pcorr",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh,
dimensionedScalar("pcorr", p_gh.dimensions(), 0.0),
pcorrTypes
);
dimensionedScalar rAUf("(1|A(U))", dimTime/rho.dimensions(), 1.0);
adjustPhi(phi, U, pcorr);
fvc::makeAbsolute(phi, U);
while (pimple.correctNonOrthogonal())
{
fvScalarMatrix pcorrEqn
(
fvm::laplacian(rAUf, pcorr) == fvc::div(phi)
);
pcorrEqn.setReference(pRefCell, pRefValue);
pcorrEqn.solve();
if (pimple.finalNonOrthogonalIter())
{
phi -= pcorrEqn.flux();
phiAbs = phi;
fvc::makeRelative(phi, U);
}
}
#include "continuityErrs.H"
}
myPotentialFreeSurfaceFoam.C
Code:
#include "fvCFD.H"
#include "singlePhaseTransportModel.H"
#include "dynamicFvMesh.H"
#include "turbulenceModel.H"
#include "pimpleControl.H"
#include "IObasicSourceList.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
#include "setRootCase.H"
#include "createTime.H"
#include "createDynamicFvMesh.H"
#include "initContinuityErrs.H"
#include "createFields.H"
#include "readTimeControls.H"
pimpleControl pimple(mesh);
surfaceScalarField phiAbs("phiAbs", phi);
fvc::makeAbsolute(phiAbs, U);
#include "correctPhi.H"
#include "CourantNo.H"
#include "setInitialDeltaT.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl;
while (runTime.run())
{
#include "readTimeControls.H"
#include "CourantNo.H"
#include "setDeltaT.H"
runTime++;
Info<< "Time = " << runTime.timeName() << nl << endl;
scalar timeBeforeMeshUpdate = runTime.elapsedCpuTime();
{
// Calculate the relative velocity used to map the relative flux phi
volVectorField Urel("Urel", U);
if (mesh.moving())
{
Urel -= fvc::reconstruct(fvc::meshPhi(U));
}
// Do any mesh changes
mesh.update();
}
/* if (mesh.changing())
{
Info<< "Execution time for mesh.update() = "
<< runTime.elapsedCpuTime() - timeBeforeMeshUpdate
<< " s" << endl;
gh = g & mesh.C();
ghf = g & mesh.Cf();
}
*/
if (mesh.changing() && correctPhi)
{
#include "correctPhi.H"
}
if (mesh.changing() && checkMeshCourantNo)
{
#include "meshCourantNo.H"
}
// --- Pressure-velocity PIMPLE corrector loop
while (pimple.loop())
{
#include "UEqn.H"
// --- Pressure corrector loop
while (pimple.correct())
{
#include "pEqn.H"
}
if (pimple.turbCorr())
{
turbulence->correct();
}
}
runTime.write();
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
}
Info<< "End\n" << endl;
return 0;
}
pEqn.H
Code:
volScalarField rAU(1.0/UEqn().A());
surfaceScalarField rAUf(rAU.name() + 'f', fvc::interpolate(rAU));
U = rAU*(UEqn() == sources(U))().H();
if (pimple.nCorrPISO() <= 1)
{
UEqn.clear();
}
phi = (fvc::interpolate(U) & mesh.Sf())
+ fvc::ddtPhiCorr(rAU, U, phiAbs);
if (p_gh.needReference())
{
fvc::makeRelative(phi, U);
adjustPhi(phi, U, p_gh);
fvc::makeAbsolute(phi, U);
}
// Non-orthogonal pressure corrector loop
while (pimple.correctNonOrthogonal())
{
fvScalarMatrix p_ghEqn
(
fvm::laplacian(rAUf, p_gh) == fvc::div(phi)
);
p_ghEqn.setReference(p_ghRefCell, p_ghRefValue);
p_ghEqn.solve(mesh.solver(p_gh.select(pimple.finalInnerIter())));
if (pimple.finalNonOrthogonalIter())
{
phi -= p_ghEqn.flux();
}
}
#include "continuityErrs.H"
phiAbs = phi;
// Make the fluxes relative to the mesh motion
fvc::makeRelative(phi, U);
// Explicitly relax pressure for momentum corrector
p_gh.relax();
p = p_gh + (g & (mesh.C() + zeta - refLevel));
U -= rAU*fvc::grad(p_gh);
U.correctBoundaryConditions();
sources.correct(U);
UEqn.H
Code:
tmp<fvVectorMatrix> UEqn
(
fvm::ddt(U)
+ fvm::div(phi, U)
+ turbulence->divDevReff(U)
);
UEqn().relax();
sources.constrain(UEqn());
if (pimple.momentumPredictor())
{
solve(UEqn() == -fvc::grad(p_gh) + sources(U));
}
|
