CFD Online Discussion Forums - Adapting PISO for condition div(fracFluid*U) = 0

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)

- OpenFOAM Programming & Development (https://www.cfd-online.com/Forums/openfoam-programming-development/)

- - Adapting PISO for condition div(fracFluid*U) = 0 (https://www.cfd-online.com/Forums/openfoam-programming-development/133186-adapting-piso-condition-div-fracfluid-u-0-a.html)

Adapting PISO for condition div(fracFluid*U) = 0

Hello everyone,

I'm studying the heat transfert and water flow in a nuclear pool.
The model i use for nuclear waste is porosity. So, i have introduced a variable "fracFluid" that represent the fraction of fluid.

So, i have 2 regions in my pool : the porous medium with phi = 0.57, and the water above with phi = 1.
Then, my continuity equation is not anymore div(U) = 0, but div(fracFluid*U) =0.

I've done some modification on my PISO loop, but it doesn't work.
Do you have any ideas why it's not good ?

here is my code :

Quote:

{
volScalarField rAU("rAU", 1.0/UEqn.A()); // (1/ap)
surfaceScalarField rAUf("Dp", fvc::interpolate(rAU)); // (1/ap)f

volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn.H(); // H(U)/ap

surfaceScalarField phig(-rAUf*ghf*fvc::snGrad(rhok)*mesh.magSf()); // (rho_k*g/ap)

surfaceScalarField phiHbyA
(
"phiHbyA",
(fvc::interpolate(HbyA) & mesh.Sf())
+ fvc::ddtPhiCorr(rAU, U, phi)
+ phig
); // (H(u)/ap)f - (rho_k*g/ap)f ddtphicorr ????

surfaceScalarField fracFluidrAUf("fracFluidrAUf", fvc::interpolate(fracFluid)*rAUf);
surfaceScalarField fracFluidphiHbyA("fracFluidphiHbyA", fvc::interpolate(fracFluid)*phiHbyA);

while (pimple.correctNonOrthogonal())
{
fvScalarMatrix p_rghEqn
(
fvm::laplacian(fracFluidrAUf, p_rgh) == fvc::div(fracFluidphiHbyA)
);

p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));

p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.fin alInnerIter())));

if (pimple.finalNonOrthogonalIter())
{
// Calculate the conservative fluxes
phi = phiHbyA - p_rghEqn.flux();

// Explicitly relax pressure for momentum corrector
p_rgh.relax();

// Correct the momentum source with the pressure gradient flux
// calculated from the relaxed pressure
U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
U.correctBoundaryConditions();
fvOptions.correct(U);
}
}

#include "continuityErrs.H"

p = p_rgh + rhok*gh;

if (p_rgh.needReference())
{
p += dimensionedScalar
(
"p",
p.dimensions(),
pRefValue - getRefCellValue(p, pRefCell)
);
p_rgh = p - rhok*gh;
}
}

my equations and the geometry in attachment.

Thanks a lot :-)

Shoji

Hi everyone !

I solved my problem. I post the solution since it may be usefull for some of you...

The trick was to make interpolation(fracFluid*rAU) instead of interpolation(fracFluid)*interpolation(rAu) and so on for all variables.

here is the new equation of p_rgh :

Quote:

{
volScalarField rAU("rAU", 1.0/UEqn.A());

// New variable to put in ddtPhiCorr
volScalarField fracFluidrAU("fracFluidrAU", rAU*fracFluid);

surfaceScalarField rAUf("Dp", fvc::interpolate(rAU));

// New variable to put in phig and laplacian of p_rghEqn
surfaceScalarField fracFluidrAUf("fracFluidrAUf", fvc::interpolate(rAU*fracFluid));

volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn.H();

surfaceScalarField phig(-fracFluidrAUf*ghf*fvc::snGrad(rhok)*mesh.magSf());

surfaceScalarField phiHbyA
(
"phiHbyA",
(fvc::interpolate(HbyA*fracFluid) & mesh.Sf()) // interpolation of (HbyA*fracFluid)
+ fvc::ddtPhiCorr(fracFluidrAU, U, phi)
+ phig
);

while (pimple.correctNonOrthogonal())
{
fvScalarMatrix p_rghEqn
(
fvm::laplacian(fracFluidrAUf, p_rgh) == fvc::div(phiHbyA)
);

p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));

p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.fin alInnerIter())));

if (pimple.finalNonOrthogonalIter())
{
// Calculate the conservative fluxes
phi = phiHbyA - p_rghEqn.flux();

// Explicitly relax pressure for momentum corrector
p_rgh.relax();

// Correct the momentum source with the pressure gradient flux
// calculated from the relaxed pressure
U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
U.correctBoundaryConditions();
fvOptions.correct(U);
}
}

#include "continuityErrs.H"

p = p_rgh + rhok*gh;

if (p_rgh.needReference())
{
p += dimensionedScalar
(
"p",
p.dimensions(),
pRefValue - getRefCellValue(p, pRefCell)
);
p_rgh = p - rhok*gh;
}
}

Best regards,

Shoji

Hell Jean,

Why is it U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf)
instead of U = HbyA - rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);

Thanks.