exaFOAM/foam-extend4.1-coherent-io
Archived
1
1
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity

Consistency update: channelFoam

Browse source
This commit is contained in:
Vuko Vukcevic 2017-01-04 14:37:08 +01:00
parent 5d7793cf06
commit cd4c08fedc
1 changed files with 24 additions and 19 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

43
applications/solvers/incompressible/channelFoam/channelFoam.C
View file

@ -26,7 +26,8 @@ Application
Description Description
Incompressible LES solver for flow in a channel. Incompressible LES solver for flow in a channel.
Consistent formulation without time-step and relaxation dependence by Jasak Consistent formulation without time-step and relaxation dependence by
Jasak and Tukovic.
Author Author
Hrvoje Jasak, Wikki Ltd. All rights reserved Hrvoje Jasak, Wikki Ltd. All rights reserved
@ -66,6 +67,9 @@ int main(int argc, char *argv[])
sgsModel->correct(); sgsModel->correct();
// Time derivative matrix
fvVectorMatrix ddtUEqn(fvm::ddt(U));
// Convection-diffusion matrix // Convection-diffusion matrix
fvVectorMatrix HUEqn fvVectorMatrix HUEqn
( (
@ -75,24 +79,23 @@ int main(int argc, char *argv[])
flowDirection*gradP flowDirection*gradP
); );
// Time derivative matrix
fvVectorMatrix ddtUEqn(fvm::ddt(U));
if (piso.momentumPredictor()) if (piso.momentumPredictor())
{ {
solve(ddtUEqn + HUEqn == -fvc::grad(p)); solve(ddtUEqn + HUEqn == -fvc::grad(p));
} }
// Prepare clean Ap without time derivative contribution // Prepare clean 1/a_p without time derivative contribution
// HJ, 26/Oct/2015 volScalarField rAU = 1.0/HUEqn.A();
volScalarField aU = HUEqn.A();
// --- PISO loop // --- PISO loop
while (piso.correct()) while (piso.correct())
{ {
U = HUEqn.H()/aU; // Calculate U from convection-diffusion matrix
phi = (fvc::interpolate(U) & mesh.Sf()); U = rAU*HUEqn.H();
// Consistently calculate flux
piso.calcTransientConsistentFlux(phi, U, rAU, ddtUEqn);
adjustPhi(phi, U, p); adjustPhi(phi, U, p);
@ -100,7 +103,14 @@ int main(int argc, char *argv[])
{ {
fvScalarMatrix pEqn fvScalarMatrix pEqn
( (
fvm::laplacian(1/aU, p) == fvc::div(phi) fvm::laplacian
(
fvc::interpolate(rAU)/piso.aCoeff(),
p,
"laplacian(rAU," + p.name() + ')'
)
==
fvc::div(phi)
); );
pEqn.setReference(pRefCell, pRefValue); pEqn.setReference(pRefCell, pRefValue);
@ -117,13 +127,8 @@ int main(int argc, char *argv[])
# include "continuityErrs.H" # include "continuityErrs.H"
// Note: cannot call H(U) here because the velocity is not complete // Consistently reconstruct velocity after pressure equation
// HJ, 22/Jan/2016 piso.reconstructTransientVelocity(U, ddtUEqn, rAU, p, phi);
U = 1.0/(aU + ddtUEqn.A())*
(
U*aU - fvc::grad(p) + ddtUEqn.H()
);
U.correctBoundaryConditions();
} }
// Correct driving force for a constant mass flow rate // Correct driving force for a constant mass flow rate
@ -135,9 +140,9 @@ int main(int argc, char *argv[])
// Calculate the pressure gradient increment needed to // Calculate the pressure gradient increment needed to
// adjust the average flow-rate to the correct value // adjust the average flow-rate to the correct value
dimensionedScalar gragPplus = dimensionedScalar gragPplus =
(magUbar - magUbarStar)*aU.weightedAverage(mesh.V()); (magUbar - magUbarStar)/rAU.weightedAverage(mesh.V());
U += gragPplus/aU*flowDirection; U += gragPplus*rAU*flowDirection;
gradP += gragPplus; gradP += gragPplus;