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

 Register Blogs Members List Search Today's Posts Mark Forums Read

April 11, 2014, 09:10
Adapting PISO for condition div(fracFluid*U) = 0
#1
New Member

Jean-François Balfroid
Join Date: Oct 2013
Posts: 2
Rep Power: 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
Attached Images
 equations.png (65.4 KB, 18 views) phi.png (7.2 KB, 13 views)

April 15, 2014, 07:34
problem solved
#2
New Member

Jean-François Balfroid
Join Date: Oct 2013
Posts: 2
Rep Power: 0
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

 June 30, 2015, 12:08 #3 Senior Member   Dongyue Li Join Date: Jun 2012 Location: Torino, Italy Posts: 675 Rep Power: 8 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. __________________ Im the translator of OpenFOAM User Guide Chinese Edition! But always newbie on CFD. Never too old to learn CFD. Worship these CFD experts.

 Thread Tools Display Modes Linear Mode

 Posting Rules You may not post new threads You may not post replies You may not post attachments You may not edit your posts BB code is On Smilies are On [IMG] code is On HTML code is OffTrackbacks are On Pingbacks are On Refbacks are On Forum Rules

 Similar Threads Thread Thread Starter Forum Replies Last Post gdeneyer OpenFOAM Programming & Development 1 August 23, 2012 05:19 nishant_hull OpenFOAM Programming & Development 0 August 25, 2009 12:48 Aditya Main CFD Forum 5 April 1, 2006 18:52 Benedikt Flurl Main CFD Forum 0 May 7, 2005 12:02 benedikt flurl Main CFD Forum 2 April 14, 2005 06:54

All times are GMT -4. The time now is 13:42.