CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM Programming & Development

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

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

Like Tree1Likes
  • 1 Post By Shoji

Reply
 
LinkBack Thread Tools Display Modes
Old   April 11, 2014, 09:10
Default Adapting PISO for condition div(fracFluid*U) = 0
  #1
New Member
 
Jean-François Balfroid
Join Date: Oct 2013
Posts: 2
Rep Power: 0
Shoji is on a distinguished road
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
File Type: png equations.png (65.4 KB, 18 views)
File Type: png phi.png (7.2 KB, 13 views)
Shoji is offline   Reply With Quote

Old   April 15, 2014, 07:34
Default problem solved
  #2
New Member
 
Jean-François Balfroid
Join Date: Oct 2013
Posts: 2
Rep Power: 0
Shoji is on a distinguished road
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
jherb likes this.
Shoji is offline   Reply With Quote

Old   June 30, 2015, 12:08
Default
  #3
Senior Member
 
Dongyue Li
Join Date: Jun 2012
Location: Torino, Italy
Posts: 675
Rep Power: 8
sharonyue is on a distinguished road
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.
sharonyue is offline   Reply With Quote

Reply

Thread Tools
Display Modes

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 Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Non-linearity Pressure Equation -- PISO algorithm gdeneyer OpenFOAM Programming & Development 1 August 23, 2012 05:19
Development of a Low mach PISO solver nishant_hull OpenFOAM Programming & Development 0 August 25, 2009 12:48
SIMPLE and PISO Aditya Main CFD Forum 5 April 1, 2006 18:52
PISO - new questions Benedikt Flurl Main CFD Forum 0 May 7, 2005 12:02
PISO vs. SIMPLE benedikt flurl Main CFD Forum 2 April 14, 2005 06:54


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