CFD Online Discussion Forums (http://www.cfd-online.com/Forums/)
-   OpenFOAM Programming & Development (http://www.cfd-online.com/Forums/openfoam-programming-development/)

 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

 sharonyue June 30, 2015 12:08

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.

 All times are GMT -4. The time now is 10:59.