Commit 079eb373 authored by sergio's avatar sergio
Browse files

ENH: Updating overPimpleDyMFoam solver

parent 4637e458
......@@ -39,19 +39,30 @@ mesh.setFluxRequired(p.name());
// Add solver-specific interpolations
{
dictionary oversetDict;
oversetDict.add("U", true);
oversetDict.add("p", true);
oversetDict.add("HbyA", true);
oversetDict.add("grad(p)", true);
dictionary suppressDict;
{
const wordHashSet& nonInt = Stencil::New(mesh).nonInterpolatedFields();
for (auto fldName : nonInt)
{
suppressDict.add(fldName, true);
}
suppressDict.add("HbyA", true);
suppressDict.add("grad(p)", true);
suppressDict.add("surfaceIntegrate(phi)", true);
suppressDict.add("surfaceIntegrate(phiHbyA)", true);
suppressDict.add("cellMask", true);
suppressDict.add("cellDisplacement", true);
suppressDict.add("interpolatedCells", true);
suppressDict.add("cellInterpolationWeight", true);
}
const_cast<dictionary&>
(
mesh.schemesDict()
).add
(
"oversetInterpolationRequired",
oversetDict,
"oversetInterpolationSuppressed",
suppressDict,
true
);
}
......@@ -68,3 +79,17 @@ autoPtr<incompressible::turbulenceModel> turbulence
(
incompressible::turbulenceModel::New(U, phi, laminarTransport)
);
surfaceScalarField phiHbyA
(
IOobject
(
"phiHbyA",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
fvc::flux(U)
);
......@@ -3,7 +3,7 @@
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
\\/ M anipulation | Copyright (C) 2016-2017 OpenCFD Ltd.
\\/ M anipulation | Copyright (C) 2016-2018 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
......@@ -94,13 +94,24 @@ int main(int argc, char *argv[])
{
#include "setCellMask.H"
#include "setInterpolatedCells.H"
surfaceScalarField faceMaskOld
(
localMin<scalar>(mesh).interpolate(cellMask.oldTime())
);
// 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;
}
// Calculate absolute flux from the mapped surface velocity
phi = mesh.Sf() & Uf;
if (mesh.changing() && correctPhi)
{
// Calculate absolute flux from the mapped surface velocity
#include "correctPhi.H"
}
......
volScalarField rAU(1.0/UEqn.A());
// Option 1: interpolate rAU, do not block out rAU on blocked cells
//mesh.interpolate(rAU, false);
//surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
volScalarField rAU("rAU", 1.0/UEqn.A());
// Option 2: do not interpolate rAU but block out rAU
//surfaceScalarField rAUf("rAUf", fvc::interpolate(blockedCells*rAU));
......@@ -15,14 +12,12 @@ volScalarField rAU(1.0/UEqn.A());
// H I C C C C
// H
//
surfaceScalarField faceMask(localMin<scalar>(mesh).interpolate(cellMask));
surfaceScalarField rAUf("rAUf", faceMask*fvc::interpolate(rAU));
surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
volVectorField H("H", UEqn.H());
volVectorField HbyA("HbyA", U);
HbyA = constrainHbyA(cellMask*rAU*UEqn.H(), U, p);
HbyA = constrainHbyA(rAU*H, U, p);
//mesh.interpolate(HbyA);
if (massFluxInterpolation)
{
#include "interpolatedFaces.H"
......@@ -33,21 +28,28 @@ if (pimple.nCorrPISO() <= 1)
tUEqn.clear();
}
surfaceScalarField phiHbyA("phiHbyA", fvc::flux(HbyA));
phiHbyA = fvc::flux(HbyA);
if (ddtCorr)
{
phiHbyA += rAUf*fvc::ddtCorr(U, Uf);
surfaceScalarField faceMaskOld
(
localMin<scalar>(mesh).interpolate(cellMask.oldTime())
);
phiHbyA += rAUf*faceMaskOld*fvc::ddtCorr(U, Uf);
}
MRF.makeRelative(phiHbyA);
// WIP
/*
if (p.needReference())
{
fvc::makeRelative(phiHbyA, U);
adjustPhi(phiHbyA, U, p);
fvc::makeAbsolute(phiHbyA, U);
}
*/
if (adjustFringe)
{
......@@ -70,21 +72,20 @@ while (pimple.correctNonOrthogonal())
if (pimple.finalNonOrthogonalIter())
{
phi = phiHbyA - pEqn.flux();
// option 2:
// rAUf*fvc::snGrad(p)*mesh.magSf();
}
}
#include "continuityErrs.H"
// Explicitly relax pressure for momentum corrector
p.relax();
volVectorField gradP(fvc::grad(p));
//mesh.interpolate(gradP);
// Option 1: leave velocity intact on blocked out cells
//U = HbyA - rAU*gradP;
// Option 2: zero out velocity on blocked out cells
U = (HbyA - rAU*cellMask*gradP);
U = HbyA - rAU*cellMask*gradP;
U.correctBoundaryConditions();
fvOptions.correct(U);
......
......@@ -335,9 +335,10 @@ void Foam::dynamicOversetFvMesh::addInterpolation(fvMatrix<Type>& m) const
lower[facei] = 0.0;
}
// For safety we make zero the HOLES
const scalar normalisation = V()[celli];
diag[celli] = normalisation;
source[celli] = normalisation*m.psi()[celli];
source[celli] = pTraits<Type>::zero;//normalisation*m.psi()[celli];
}
}
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment