CFD Online Discussion Forums

CFD Online Discussion Forums (http://www.cfd-online.com/Forums/)
-   OpenFOAM Programming & Development (http://www.cfd-online.com/Forums/openfoam-programming-development/)
-   -   Adapting PISO for condition div(fracFluid*U) = 0 (http://www.cfd-online.com/Forums/openfoam-programming-development/133186-adapting-piso-condition-div-fracfluid-u-0-a.html)

Shoji April 11, 2014 09:10

Adapting PISO for condition div(fracFluid*U) = 0
 
2 Attachment(s)
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

Shoji April 15, 2014 07:34

problem solved
 
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


All times are GMT -4. The time now is 08:25.