Merge Request #30: Time consistent incompressible solvers. Author: Hrvoje Jasak. Merge: Henrik Rusche
https://sourceforge.net/p/foam-extend/foam-extend-3.2/merge-requests/30/ Conflicts: applications/solvers/incompressible/simpleFoam/pEqn.H applications/solvers/incompressible/simpleSRFFoam/simpleSRFFoam.C tutorials/incompressible/pimpleDyMFoam/movingCone/system/fvSchemes tutorials/incompressible/simpleFoam/motorBike/system/fvSchemes
This commit is contained in:
commit
b858c51ff2
32 changed files with 397 additions and 234 deletions
|
@ -26,6 +26,10 @@ Application
|
|||
|
||||
Description
|
||||
Incompressible LES solver for flow in a channel.
|
||||
Consistent formulation without time-step and relaxation dependence by Jasak
|
||||
|
||||
Author
|
||||
Hrvoje Jasak, Wikki Ltd. All rights reserved
|
||||
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
|
@ -60,38 +64,41 @@ int main(int argc, char *argv[])
|
|||
|
||||
sgsModel->correct();
|
||||
|
||||
fvVectorMatrix UEqn
|
||||
// Convection-diffusion matrix
|
||||
fvVectorMatrix HUEqn
|
||||
(
|
||||
fvm::ddt(U)
|
||||
+ fvm::div(phi, U)
|
||||
fvm::div(phi, U)
|
||||
+ sgsModel->divDevBeff(U)
|
||||
==
|
||||
flowDirection*gradP
|
||||
);
|
||||
|
||||
// Time derivative matrix
|
||||
fvVectorMatrix ddtUEqn(fvm::ddt(U));
|
||||
|
||||
if (momentumPredictor)
|
||||
{
|
||||
solve(UEqn == -fvc::grad(p));
|
||||
solve(ddtUEqn + HUEqn == -fvc::grad(p));
|
||||
}
|
||||
|
||||
// Prepare clean Ap without time derivative contribution
|
||||
// HJ, 26/Oct/2015
|
||||
volScalarField aU = HUEqn.A();
|
||||
|
||||
// --- PISO loop
|
||||
|
||||
volScalarField rUA = 1.0/UEqn.A();
|
||||
|
||||
for (int corr = 0; corr < nCorr; corr++)
|
||||
{
|
||||
U = rUA*UEqn.H();
|
||||
phi = (fvc::interpolate(U) & mesh.Sf())
|
||||
+ fvc::ddtPhiCorr(rUA, U, phi);
|
||||
U = HUEqn.H()/aU;
|
||||
phi = (fvc::interpolate(U) & mesh.Sf());
|
||||
|
||||
adjustPhi(phi, U, p);
|
||||
|
||||
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
|
||||
for (int nonOrth = 0; nonOrth <= nNonOrthCorr; nonOrth++)
|
||||
{
|
||||
fvScalarMatrix pEqn
|
||||
(
|
||||
fvm::laplacian(rUA, p) == fvc::div(phi)
|
||||
fvm::laplacian(1/aU, p) == fvc::div(phi)
|
||||
);
|
||||
|
||||
pEqn.setReference(pRefCell, pRefValue);
|
||||
|
@ -111,13 +118,17 @@ int main(int argc, char *argv[])
|
|||
}
|
||||
}
|
||||
|
||||
#include "continuityErrs.H"
|
||||
# include "continuityErrs.H"
|
||||
|
||||
U -= rUA*fvc::grad(p);
|
||||
// Note: cannot call H(U) here because the velocity is not complete
|
||||
// HJ, 22/Jan/2016
|
||||
U = 1.0/(aU + ddtUEqn.A())*
|
||||
(
|
||||
U*aU - fvc::grad(p) + ddtUEqn.H()
|
||||
);
|
||||
U.correctBoundaryConditions();
|
||||
}
|
||||
|
||||
|
||||
// Correct driving force for a constant mass flow rate
|
||||
|
||||
// Extract the velocity in the flow direction
|
||||
|
@ -127,9 +138,9 @@ int main(int argc, char *argv[])
|
|||
// Calculate the pressure gradient increment needed to
|
||||
// adjust the average flow-rate to the correct value
|
||||
dimensionedScalar gragPplus =
|
||||
(magUbar - magUbarStar)/rUA.weightedAverage(mesh.V());
|
||||
(magUbar - magUbarStar)*aU.weightedAverage(mesh.V());
|
||||
|
||||
U += flowDirection*rUA*gragPplus;
|
||||
U += gragPplus/aU*flowDirection;
|
||||
|
||||
gradP += gragPplus;
|
||||
|
||||
|
|
|
@ -26,6 +26,10 @@ Application
|
|||
|
||||
Description
|
||||
Transient solver for incompressible, laminar flow of Newtonian fluids.
|
||||
Consistent formulation without time-step and relaxation dependence by Jasak
|
||||
|
||||
Author
|
||||
Hrvoje Jasak, Wikki Ltd. All rights reserved
|
||||
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
|
@ -53,32 +57,36 @@ int main(int argc, char *argv[])
|
|||
# include "readPISOControls.H"
|
||||
# include "CourantNo.H"
|
||||
|
||||
fvVectorMatrix UEqn
|
||||
// Convection-diffusion matrix
|
||||
fvVectorMatrix HUEqn
|
||||
(
|
||||
fvm::ddt(U)
|
||||
+ fvm::div(phi, U)
|
||||
fvm::div(phi, U)
|
||||
- fvm::laplacian(nu, U)
|
||||
);
|
||||
|
||||
solve(UEqn == -fvc::grad(p));
|
||||
// Time derivative matrix
|
||||
fvVectorMatrix ddtUEqn(fvm::ddt(U));
|
||||
|
||||
solve(ddtUEqn + HUEqn == -fvc::grad(p));
|
||||
|
||||
// Prepare clean Ap without time derivative contribution
|
||||
// HJ, 26/Oct/2015
|
||||
volScalarField aU = HUEqn.A();
|
||||
|
||||
// --- PISO loop
|
||||
|
||||
for (int corr = 0; corr < nCorr; corr++)
|
||||
{
|
||||
volScalarField rUA = 1.0/UEqn.A();
|
||||
|
||||
U = rUA*UEqn.H();
|
||||
phi = (fvc::interpolate(U) & mesh.Sf())
|
||||
+ fvc::ddtPhiCorr(rUA, U, phi);
|
||||
U = HUEqn.H()/aU;
|
||||
phi = (fvc::interpolate(U) & mesh.Sf());
|
||||
|
||||
adjustPhi(phi, U, p);
|
||||
|
||||
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
|
||||
for (int nonOrth = 0; nonOrth <= nNonOrthCorr; nonOrth++)
|
||||
{
|
||||
fvScalarMatrix pEqn
|
||||
(
|
||||
fvm::laplacian(rUA, p) == fvc::div(phi)
|
||||
fvm::laplacian(1/aU, p) == fvc::div(phi)
|
||||
);
|
||||
|
||||
pEqn.setReference(pRefCell, pRefValue);
|
||||
|
@ -92,7 +100,12 @@ int main(int argc, char *argv[])
|
|||
|
||||
# include "continuityErrs.H"
|
||||
|
||||
U -= rUA*fvc::grad(p);
|
||||
// Note: cannot call H(U) here because the velocity is not complete
|
||||
// HJ, 22/Jan/2016
|
||||
U = 1.0/(aU + ddtUEqn.A())*
|
||||
(
|
||||
U*aU - fvc::grad(p) + ddtUEqn.H()
|
||||
);
|
||||
U.correctBoundaryConditions();
|
||||
}
|
||||
|
||||
|
|
|
@ -26,6 +26,10 @@ Application
|
|||
|
||||
Description
|
||||
Transient solver for incompressible, laminar flow of non-Newtonian fluids.
|
||||
Consistent formulation without time-step and relaxation dependence by Jasak
|
||||
|
||||
Author
|
||||
Hrvoje Jasak, Wikki Ltd. All rights reserved
|
||||
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
|
@ -39,7 +43,7 @@ int main(int argc, char *argv[])
|
|||
# include "setRootCase.H"
|
||||
|
||||
# include "createTime.H"
|
||||
# include "createMeshNoClear.H"
|
||||
# include "createMesh.H"
|
||||
# include "createFields.H"
|
||||
# include "initContinuityErrs.H"
|
||||
|
||||
|
@ -56,32 +60,36 @@ int main(int argc, char *argv[])
|
|||
|
||||
fluid.correct();
|
||||
|
||||
fvVectorMatrix UEqn
|
||||
// Convection-diffusion matrix
|
||||
fvVectorMatrix HUEqn
|
||||
(
|
||||
fvm::ddt(U)
|
||||
+ fvm::div(phi, U)
|
||||
fvm::div(phi, U)
|
||||
- fvm::laplacian(fluid.nu(), U)
|
||||
);
|
||||
|
||||
solve(UEqn == -fvc::grad(p));
|
||||
// Time derivative matrix
|
||||
fvVectorMatrix ddtUEqn(fvm::ddt(U));
|
||||
|
||||
solve(ddtUEqn + HUEqn == -fvc::grad(p));
|
||||
|
||||
// Prepare clean Ap without time derivative contribution
|
||||
// HJ, 26/Oct/2015
|
||||
volScalarField aU = HUEqn.A();
|
||||
|
||||
// --- PISO loop
|
||||
|
||||
for (int corr = 0; corr < nCorr; corr++)
|
||||
{
|
||||
volScalarField rUA = 1.0/UEqn.A();
|
||||
|
||||
U = rUA*UEqn.H();
|
||||
phi = (fvc::interpolate(U) & mesh.Sf())
|
||||
+ fvc::ddtPhiCorr(rUA, U, phi);
|
||||
U = HUEqn.H()/aU;
|
||||
phi = (fvc::interpolate(U) & mesh.Sf());
|
||||
|
||||
adjustPhi(phi, U, p);
|
||||
|
||||
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
|
||||
for (int nonOrth = 0; nonOrth <= nNonOrthCorr; nonOrth++)
|
||||
{
|
||||
fvScalarMatrix pEqn
|
||||
(
|
||||
fvm::laplacian(rUA, p) == fvc::div(phi)
|
||||
fvm::laplacian(1/aU, p) == fvc::div(phi)
|
||||
);
|
||||
|
||||
pEqn.setReference(pRefCell, pRefValue);
|
||||
|
@ -95,7 +103,12 @@ int main(int argc, char *argv[])
|
|||
|
||||
# include "continuityErrs.H"
|
||||
|
||||
U -= rUA*fvc::grad(p);
|
||||
// Note: cannot call H(U) here because the velocity is not complete
|
||||
// HJ, 22/Jan/2016
|
||||
U = 1.0/(aU + ddtUEqn.A())*
|
||||
(
|
||||
U*aU - fvc::grad(p) + ddtUEqn.H()
|
||||
);
|
||||
U.correctBoundaryConditions();
|
||||
}
|
||||
|
||||
|
|
|
@ -1,27 +1,33 @@
|
|||
fvVectorMatrix UEqn
|
||||
// Convection-diffusion matrix
|
||||
fvVectorMatrix HUEqn
|
||||
(
|
||||
fvm::ddt(U)
|
||||
+ fvm::div(phi, U)
|
||||
fvm::div(phi, U)
|
||||
+ turbulence->divDevReff(U)
|
||||
);
|
||||
|
||||
// Time derivative matrix
|
||||
fvVectorMatrix ddtUEqn(fvm::ddt(U));
|
||||
|
||||
// Get under-relaxation factor
|
||||
scalar UUrf = mesh.solutionDict().relaxationFactor(U.name());
|
||||
|
||||
if (oCorr == nOuterCorr - 1)
|
||||
{
|
||||
if (mesh.solutionDict().relax("UFinal"))
|
||||
{
|
||||
UEqn.relax(mesh.solutionDict().relaxationFactor("UFinal"));
|
||||
UUrf = mesh.solutionDict().relaxationFactor("UFinal");
|
||||
}
|
||||
else
|
||||
{
|
||||
UEqn.relax(1);
|
||||
UUrf = 1;
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
UEqn.relax();
|
||||
}
|
||||
|
||||
if (momentumPredictor)
|
||||
{
|
||||
solve(UEqn == -fvc::grad(p));
|
||||
}
|
||||
// Solve momentum predictor
|
||||
solve
|
||||
(
|
||||
ddtUEqn
|
||||
+ relax(HUEqn, UUrf)
|
||||
==
|
||||
- fvc::grad(p)
|
||||
);
|
||||
|
|
|
@ -39,7 +39,7 @@
|
|||
{
|
||||
fvScalarMatrix pcorrEqn
|
||||
(
|
||||
fvm::laplacian(rAU, pcorr) == fvc::div(phi)
|
||||
fvm::laplacian(1/aU, pcorr) == fvc::div(phi)
|
||||
);
|
||||
|
||||
pcorrEqn.setReference(pRefCell, pRefValue);
|
||||
|
|
|
@ -41,18 +41,18 @@
|
|||
incompressible::turbulenceModel::New(U, phi, laminarTransport)
|
||||
);
|
||||
|
||||
Info<< "Reading field rAU if present\n" << endl;
|
||||
volScalarField rAU
|
||||
Info<< "Reading field aU if present\n" << endl;
|
||||
volScalarField aU
|
||||
(
|
||||
IOobject
|
||||
(
|
||||
"rAU",
|
||||
"aU",
|
||||
runTime.timeName(),
|
||||
mesh,
|
||||
IOobject::READ_IF_PRESENT,
|
||||
IOobject::AUTO_WRITE
|
||||
),
|
||||
mesh,
|
||||
runTime.deltaT(),
|
||||
1/runTime.deltaT(),
|
||||
zeroGradientFvPatchScalarField::typeName
|
||||
);
|
||||
|
|
70
applications/solvers/incompressible/pimpleDyMFoam/pEqn.H
Normal file
70
applications/solvers/incompressible/pimpleDyMFoam/pEqn.H
Normal file
|
@ -0,0 +1,70 @@
|
|||
{
|
||||
p.boundaryField().updateCoeffs();
|
||||
|
||||
// Prepare clean Ap without time derivative contribution and
|
||||
// without contribution from under-relaxation
|
||||
// HJ, 26/Oct/2015
|
||||
aU = HUEqn.A();
|
||||
|
||||
// Store velocity under-relaxation point before using U for
|
||||
// the flux precursor
|
||||
U.storePrevIter();
|
||||
|
||||
U = HUEqn.H()/aU;
|
||||
phi = (fvc::interpolate(U) & mesh.Sf());
|
||||
|
||||
// Global flux balance
|
||||
adjustPhi(phi, U, p);
|
||||
|
||||
for (int nonOrth = 0; nonOrth <= nNonOrthCorr; nonOrth++)
|
||||
{
|
||||
fvScalarMatrix pEqn
|
||||
(
|
||||
fvm::laplacian(1/aU, p) == fvc::div(phi)
|
||||
);
|
||||
|
||||
pEqn.setReference(pRefCell, pRefValue);
|
||||
|
||||
if
|
||||
(
|
||||
corr == nCorr - 1
|
||||
&& nonOrth == nNonOrthCorr
|
||||
)
|
||||
{
|
||||
pEqn.solve
|
||||
(
|
||||
mesh.solutionDict().solver(p.name() + "Final")
|
||||
);
|
||||
}
|
||||
else
|
||||
{
|
||||
pEqn.solve(mesh.solutionDict().solver(p.name()));
|
||||
}
|
||||
|
||||
if (nonOrth == nNonOrthCorr)
|
||||
{
|
||||
phi -= pEqn.flux();
|
||||
}
|
||||
}
|
||||
|
||||
// Explicitly relax pressure for momentum corrector
|
||||
if (oCorr != nOuterCorr - 1)
|
||||
{
|
||||
p.relax();
|
||||
}
|
||||
|
||||
// Make the fluxes relative to the mesh motion
|
||||
fvc::makeRelative(phi, U);
|
||||
|
||||
# include "movingMeshContinuityErrs.H"
|
||||
|
||||
U = UUrf*
|
||||
(
|
||||
1.0/(aU + ddtUEqn.A())*
|
||||
(
|
||||
U*aU - fvc::grad(p) + ddtUEqn.H()
|
||||
)
|
||||
)
|
||||
+ (1 - UUrf)*U.prevIter();
|
||||
U.correctBoundaryConditions();
|
||||
}
|
|
@ -30,6 +30,11 @@ Description
|
|||
|
||||
Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected.
|
||||
|
||||
Consistent formulation without time-step and relaxation dependence by Jasak
|
||||
|
||||
Author
|
||||
Hrvoje Jasak, Wikki Ltd. All rights reserved
|
||||
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
#include "fvCFD.H"
|
||||
|
@ -106,61 +111,7 @@ int main(int argc, char *argv[])
|
|||
// --- PISO loop
|
||||
for (int corr = 0; corr < nCorr; corr++)
|
||||
{
|
||||
rAU = 1.0/UEqn.A();
|
||||
|
||||
U = rAU*UEqn.H();
|
||||
phi = (fvc::interpolate(U) & mesh.Sf());
|
||||
// ddtPhiCorr does not work. HJ, 20/Nov/2013
|
||||
|
||||
adjustPhi(phi, U, p);
|
||||
|
||||
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
|
||||
{
|
||||
fvScalarMatrix pEqn
|
||||
(
|
||||
fvm::laplacian(rAU, p) == fvc::div(phi)
|
||||
);
|
||||
|
||||
pEqn.setReference(pRefCell, pRefValue);
|
||||
|
||||
if
|
||||
(
|
||||
// oCorr == nOuterCorr - 1
|
||||
corr == nCorr - 1
|
||||
&& nonOrth == nNonOrthCorr
|
||||
)
|
||||
{
|
||||
pEqn.solve
|
||||
(
|
||||
mesh.solutionDict().solver(p.name() + "Final")
|
||||
);
|
||||
}
|
||||
else
|
||||
{
|
||||
pEqn.solve(mesh.solutionDict().solver(p.name()));
|
||||
}
|
||||
|
||||
if (nonOrth == nNonOrthCorr)
|
||||
{
|
||||
phi -= pEqn.flux();
|
||||
}
|
||||
}
|
||||
|
||||
# include "continuityErrs.H"
|
||||
|
||||
// Explicitly relax pressure for momentum corrector
|
||||
if (oCorr != nOuterCorr - 1)
|
||||
{
|
||||
p.relax();
|
||||
}
|
||||
|
||||
// Make the fluxes relative to the mesh motion
|
||||
fvc::makeRelative(phi, U);
|
||||
|
||||
# include "movingMeshContinuityErrs.H"
|
||||
|
||||
U -= rAU*fvc::grad(p);
|
||||
U.correctBoundaryConditions();
|
||||
# include "pEqn.H"
|
||||
}
|
||||
|
||||
turbulence->correct();
|
||||
|
|
|
@ -1,43 +1,40 @@
|
|||
// Solve the Momentum equation
|
||||
|
||||
tmp<fvVectorMatrix> UEqn
|
||||
// Convection-diffusion matrx
|
||||
fvVectorMatrix HUEqn
|
||||
(
|
||||
fvm::ddt(U)
|
||||
+ fvm::div(phi, U)
|
||||
fvm::div(phi, U)
|
||||
+ turbulence->divDevReff(U)
|
||||
);
|
||||
|
||||
// Time derivative matrix
|
||||
fvVectorMatrix ddtUEqn(fvm::ddt(U));
|
||||
|
||||
// Get under-relaxation factor
|
||||
scalar UUrf = mesh.solutionDict().relaxationFactor(U.name());
|
||||
|
||||
if (oCorr == nOuterCorr - 1)
|
||||
{
|
||||
if (mesh.solutionDict().relax("UFinal"))
|
||||
{
|
||||
UEqn().relax(mesh.solutionDict().relaxationFactor("UFinal"));
|
||||
UUrf = mesh.solutionDict().relaxationFactor("UFinal");
|
||||
}
|
||||
else
|
||||
{
|
||||
UEqn().relax(1);
|
||||
UUrf = 1;
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
UEqn().relax();
|
||||
}
|
||||
|
||||
volScalarField rUA = 1.0/UEqn().A();
|
||||
|
||||
if (momentumPredictor)
|
||||
{
|
||||
if (oCorr == nOuterCorr-1)
|
||||
{
|
||||
solve(UEqn() == -fvc::grad(p), mesh.solutionDict().solver("UFinal"));
|
||||
}
|
||||
else
|
||||
{
|
||||
solve(UEqn() == -fvc::grad(p));
|
||||
}
|
||||
solve
|
||||
(
|
||||
ddtUEqn
|
||||
+ relax(HUEqn, UUrf)
|
||||
==
|
||||
- fvc::grad(p)
|
||||
);
|
||||
}
|
||||
else
|
||||
{
|
||||
U = rUA*(UEqn().H() - fvc::grad(p));
|
||||
U = (ddtUEqn.H() + HUEqn.H() - fvc::grad(p))/(HUEqn.A() + ddtUEqn.A());
|
||||
U.correctBoundaryConditions();
|
||||
}
|
||||
|
|
|
@ -40,3 +40,19 @@ autoPtr<incompressible::turbulenceModel> turbulence
|
|||
(
|
||||
incompressible::turbulenceModel::New(U, phi, laminarTransport)
|
||||
);
|
||||
|
||||
Info<< "Reading field aU if present\n" << endl;
|
||||
volScalarField aU
|
||||
(
|
||||
IOobject
|
||||
(
|
||||
"aU",
|
||||
runTime.timeName(),
|
||||
mesh,
|
||||
IOobject::READ_IF_PRESENT,
|
||||
IOobject::AUTO_WRITE
|
||||
),
|
||||
mesh,
|
||||
1/runTime.deltaT(),
|
||||
zeroGradientFvPatchScalarField::typeName
|
||||
);
|
||||
|
|
|
@ -1,53 +1,67 @@
|
|||
U = rUA*UEqn().H();
|
||||
|
||||
if (nCorr <= 1)
|
||||
{
|
||||
UEqn.clear();
|
||||
}
|
||||
p.boundaryField().updateCoeffs();
|
||||
|
||||
phi = (fvc::interpolate(U) & mesh.Sf())
|
||||
+ fvc::ddtPhiCorr(rUA, U, phi);
|
||||
// Prepare clean Ap without time derivative contribution and
|
||||
// without contribution from under-relaxation
|
||||
// HJ, 26/Oct/2015
|
||||
aU = HUEqn.A();
|
||||
|
||||
adjustPhi(phi, U, p);
|
||||
// Store velocity under-relaxation point before using U for the flux
|
||||
// precursor
|
||||
U.storePrevIter();
|
||||
|
||||
// Non-orthogonal pressure corrector loop
|
||||
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
|
||||
{
|
||||
// Pressure corrector
|
||||
fvScalarMatrix pEqn
|
||||
(
|
||||
fvm::laplacian(rUA, p) == fvc::div(phi)
|
||||
);
|
||||
U = HUEqn.H()/aU;
|
||||
phi = (fvc::interpolate(U) & mesh.Sf());
|
||||
|
||||
pEqn.setReference(pRefCell, pRefValue);
|
||||
// Global flux balance
|
||||
adjustPhi(phi, U, p);
|
||||
|
||||
if
|
||||
(
|
||||
oCorr == nOuterCorr - 1
|
||||
&& corr == nCorr - 1
|
||||
&& nonOrth == nNonOrthCorr
|
||||
)
|
||||
for (int nonOrth = 0; nonOrth <= nNonOrthCorr; nonOrth++)
|
||||
{
|
||||
pEqn.solve(mesh.solutionDict().solver("pFinal"));
|
||||
}
|
||||
else
|
||||
{
|
||||
pEqn.solve();
|
||||
fvScalarMatrix pEqn
|
||||
(
|
||||
fvm::laplacian(1/aU, p) == fvc::div(phi)
|
||||
);
|
||||
|
||||
pEqn.setReference(pRefCell, pRefValue);
|
||||
|
||||
if
|
||||
(
|
||||
corr == nCorr - 1
|
||||
&& nonOrth == nNonOrthCorr
|
||||
)
|
||||
{
|
||||
pEqn.solve
|
||||
(
|
||||
mesh.solutionDict().solver(p.name() + "Final")
|
||||
);
|
||||
}
|
||||
else
|
||||
{
|
||||
pEqn.solve(mesh.solutionDict().solver(p.name()));
|
||||
}
|
||||
|
||||
if (nonOrth == nNonOrthCorr)
|
||||
{
|
||||
phi -= pEqn.flux();
|
||||
}
|
||||
}
|
||||
|
||||
if (nonOrth == nNonOrthCorr)
|
||||
// Explicitly relax pressure for momentum corrector
|
||||
if (oCorr != nOuterCorr - 1)
|
||||
{
|
||||
phi -= pEqn.flux();
|
||||
p.relax();
|
||||
}
|
||||
|
||||
# include "movingMeshContinuityErrs.H"
|
||||
|
||||
U = UUrf*
|
||||
(
|
||||
1.0/(aU + ddtUEqn.A())*
|
||||
(
|
||||
U*aU - fvc::grad(p) + ddtUEqn.H()
|
||||
)
|
||||
)
|
||||
+ (1 - UUrf)*U.prevIter();
|
||||
U.correctBoundaryConditions();
|
||||
}
|
||||
|
||||
#include "continuityErrs.H"
|
||||
|
||||
// Explicitly relax pressure for momentum corrector except for last corrector
|
||||
if (oCorr != nOuterCorr-1)
|
||||
{
|
||||
p.relax();
|
||||
}
|
||||
|
||||
U -= rUA*fvc::grad(p);
|
||||
U.correctBoundaryConditions();
|
||||
|
|
|
@ -30,6 +30,11 @@ Description
|
|||
|
||||
Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected.
|
||||
|
||||
Consistent formulation without time-step and relaxation dependence by Jasak
|
||||
|
||||
Author
|
||||
Hrvoje Jasak, Wikki Ltd. All rights reserved
|
||||
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
#include "fvCFD.H"
|
||||
|
|
|
@ -1,16 +1,20 @@
|
|||
// Solve the momentum equation
|
||||
|
||||
tmp<fvVectorMatrix> UEqn
|
||||
tmp<fvVectorMatrix> HUEqn
|
||||
(
|
||||
fvm::div(phi, U)
|
||||
+ turbulence->divDevReff(U)
|
||||
);
|
||||
|
||||
UEqn().relax();
|
||||
// Get under-relaxation factor
|
||||
const scalar UUrf = mesh.solutionDict().relaxationFactor(U.name());
|
||||
|
||||
// Momentum solution
|
||||
eqnResidual = solve
|
||||
(
|
||||
UEqn() == -fvc::grad(p)
|
||||
relax(HUEqn(), UUrf)
|
||||
==
|
||||
-fvc::grad(p)
|
||||
).initialResidual();
|
||||
|
||||
maxResidual = max(eqnResidual, maxResidual);
|
||||
|
|
|
@ -1,17 +1,28 @@
|
|||
{
|
||||
volScalarField AU = UEqn().A();
|
||||
U = UEqn().H()/AU;
|
||||
UEqn.clear();
|
||||
phi = fvc::interpolate(U) & mesh.Sf();
|
||||
p.boundaryField().updateCoeffs();
|
||||
|
||||
// Prepare clean 1/Ap without contribution from under-relaxation
|
||||
// HJ, 26/Oct/2015
|
||||
volScalarField rUA
|
||||
(
|
||||
"(1|A(U))",
|
||||
1/HUEqn().A()
|
||||
);
|
||||
|
||||
// Store velocity under-relaxation point before using U for
|
||||
// the flux precursor
|
||||
U.storePrevIter();
|
||||
|
||||
U = rUA*HUEqn().H();
|
||||
HUEqn.clear();
|
||||
phi = fvc::interpolate(U) & mesh.Sf();
|
||||
adjustPhi(phi, U, p);
|
||||
|
||||
// Non-orthogonal pressure corrector loop
|
||||
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
|
||||
for (int nonOrth = 0; nonOrth <= nNonOrthCorr; nonOrth++)
|
||||
{
|
||||
fvScalarMatrix pEqn
|
||||
(
|
||||
fvm::laplacian(1.0/AU, p) == fvc::div(phi)
|
||||
fvm::laplacian(rUA, p) == fvc::div(phi)
|
||||
);
|
||||
|
||||
pEqn.setReference(pRefCell, pRefValue);
|
||||
|
@ -39,6 +50,8 @@
|
|||
p.relax();
|
||||
|
||||
// Momentum corrector
|
||||
U -= fvc::grad(p)/AU;
|
||||
// Note: since under-relaxation does not change aU, H/a in U can be
|
||||
// re-used. HJ, 22/Jan/2016
|
||||
U = UUrf*(U - rUA*fvc::grad(p)) + (1 - UUrf)*U.prevIter();
|
||||
U.correctBoundaryConditions();
|
||||
}
|
||||
|
||||
|
|
|
@ -26,6 +26,10 @@ Application
|
|||
|
||||
Description
|
||||
Steady-state solver for incompressible, turbulent flow
|
||||
Consistent formulation without time-step and relaxation dependence by Jasak
|
||||
|
||||
Author
|
||||
Hrvoje Jasak, Wikki Ltd. All rights reserved
|
||||
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
|
|
14
applications/solvers/incompressible/simpleSRFFoam/UEqn.H
Normal file
14
applications/solvers/incompressible/simpleSRFFoam/UEqn.H
Normal file
|
@ -0,0 +1,14 @@
|
|||
// Solve the momentum equation
|
||||
|
||||
tmp<fvVectorMatrix> HUrelEqn
|
||||
(
|
||||
fvm::div(phi, Urel)
|
||||
+ turbulence->divDevReff(Urel)
|
||||
+ SRF->Su()
|
||||
);
|
||||
|
||||
// Get under-relaxation factor
|
||||
const scalar UUrf = mesh.solutionDict().relaxationFactor(Urel.name());
|
||||
|
||||
// Momentum solution
|
||||
solve(relax(HUrelEqn(), UUrf) == -fvc::grad(p));
|
|
@ -71,4 +71,3 @@
|
|||
),
|
||||
Urel + SRF->U()
|
||||
);
|
||||
|
||||
|
|
48
applications/solvers/incompressible/simpleSRFFoam/pEqn.H
Normal file
48
applications/solvers/incompressible/simpleSRFFoam/pEqn.H
Normal file
|
@ -0,0 +1,48 @@
|
|||
p.boundaryField().updateCoeffs();
|
||||
|
||||
// Prepare clean 1/Ap without contribution from under-relaxation
|
||||
// HJ, 26/Oct/2015
|
||||
volScalarField rUA
|
||||
(
|
||||
"(1|A(Urel))",
|
||||
1/HUrelEqn().A()
|
||||
);
|
||||
|
||||
// Store velocity under-relaxation point before using U for the flux
|
||||
// precursor
|
||||
Urel.storePrevIter();
|
||||
|
||||
Urel = rUA*HUrelEqn().H();
|
||||
HUrelEqn.clear();
|
||||
phi = fvc::interpolate(Urel) & mesh.Sf();
|
||||
|
||||
// Global flux balance
|
||||
adjustPhi(phi, Urel, p);
|
||||
|
||||
// Non-orthogonal pressure corrector loop
|
||||
for (int nonOrth = 0; nonOrth <= nNonOrthCorr; nonOrth++)
|
||||
{
|
||||
fvScalarMatrix pEqn
|
||||
(
|
||||
fvm::laplacian(rUA, p) == fvc::div(phi)
|
||||
);
|
||||
|
||||
pEqn.setReference(pRefCell, pRefValue);
|
||||
pEqn.solve();
|
||||
|
||||
if (nonOrth == nNonOrthCorr)
|
||||
{
|
||||
phi -= pEqn.flux();
|
||||
}
|
||||
}
|
||||
|
||||
# include "continuityErrs.H"
|
||||
|
||||
// Explicitly relax pressure for momentum corrector
|
||||
p.relax();
|
||||
|
||||
// Momentum corrector
|
||||
// Note: since under-relaxation does not change aU, H/a in U can be
|
||||
// re-used. HJ, 22/Jan/2016
|
||||
Urel = UUrf*(Urel - rUA*fvc::grad(p)) + (1 - UUrf)*Urel.prevIter();
|
||||
Urel.correctBoundaryConditions();
|
|
@ -25,8 +25,12 @@ Application
|
|||
simpleSRFFoam
|
||||
|
||||
Description
|
||||
Steady-state solver for incompressible, turbulent flow of non-Newtonian
|
||||
Steady-state solver for incompressible, turbulent flow of Newtonian
|
||||
fluids with single rotating frame.
|
||||
Consistent formulation without time-step and relaxation dependence by Jasak
|
||||
|
||||
Author
|
||||
Hrvoje Jasak, Wikki Ltd. All rights reserved
|
||||
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
|
@ -47,8 +51,6 @@ int main(int argc, char *argv[])
|
|||
# include "createFields.H"
|
||||
# include "initContinuityErrs.H"
|
||||
|
||||
//mesh.clearPrimitives();
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
Info<< "\nStarting time loop\n" << endl;
|
||||
|
@ -75,6 +77,7 @@ int main(int argc, char *argv[])
|
|||
|
||||
solve(UrelEqn() == -fvc::grad(p));
|
||||
|
||||
p.boundaryField().updateCoeffs();
|
||||
volScalarField AUrel = UrelEqn().A();
|
||||
Urel = UrelEqn().H()/AUrel;
|
||||
UrelEqn.clear();
|
||||
|
|
|
@ -29,7 +29,7 @@ gradSchemes
|
|||
divSchemes
|
||||
{
|
||||
default none;
|
||||
div(phi,U) Gauss linear;
|
||||
div(phi,U) Gauss linearUpwind Gauss linear;
|
||||
div(phi,k) Gauss limitedLinear 1;
|
||||
div(phi,B) Gauss limitedLinear 1;
|
||||
div(B) Gauss linear;
|
||||
|
|
|
@ -28,7 +28,7 @@ gradSchemes
|
|||
divSchemes
|
||||
{
|
||||
default none;
|
||||
div(phi,U) Gauss limitedLinearV 1;
|
||||
div(phi,U) Gauss linearUpwind leastSquares;
|
||||
}
|
||||
|
||||
laplacianSchemes
|
||||
|
|
|
@ -54,12 +54,6 @@ boundaryField
|
|||
type mixingPlane;
|
||||
}
|
||||
|
||||
top
|
||||
{
|
||||
type fixedValue;
|
||||
value uniform (0 0 0);
|
||||
}
|
||||
|
||||
frontAndBack
|
||||
{
|
||||
type empty;
|
||||
|
|
|
@ -54,12 +54,6 @@ boundaryField
|
|||
type mixingPlane;
|
||||
}
|
||||
|
||||
top
|
||||
{
|
||||
type fixedValue;
|
||||
value uniform (0 0 0);
|
||||
}
|
||||
|
||||
frontAndBack
|
||||
{
|
||||
type empty;
|
||||
|
|
|
@ -54,12 +54,6 @@ boundaryField
|
|||
type mixingPlane;
|
||||
}
|
||||
|
||||
top
|
||||
{
|
||||
type fixedValue;
|
||||
value uniform (0 0 0);
|
||||
}
|
||||
|
||||
frontAndBack
|
||||
{
|
||||
type empty;
|
||||
|
|
|
@ -28,7 +28,7 @@ gradSchemes
|
|||
divSchemes
|
||||
{
|
||||
default none;
|
||||
div(phi,U) Gauss linear;
|
||||
div(phi,U) Gauss linearUpwind leastSquares;
|
||||
}
|
||||
|
||||
laplacianSchemes
|
||||
|
|
|
@ -29,16 +29,16 @@ gradSchemes
|
|||
divSchemes
|
||||
{
|
||||
default none;
|
||||
div(phi,U) Gauss linear;
|
||||
div((nuEff*dev(T(grad(U))))) Gauss linear;
|
||||
div(phi,U) Gauss linearUpwind Gauss linear;
|
||||
div((nuEff*dev(grad(U).T()))) Gauss linear;
|
||||
}
|
||||
|
||||
laplacianSchemes
|
||||
{
|
||||
default none;
|
||||
laplacian(nu,U) Gauss linear corrected;
|
||||
laplacian(rAU,pcorr) Gauss linear corrected;
|
||||
laplacian(rAU,p) Gauss linear corrected;
|
||||
laplacian((1|aU),pcorr) Gauss linear corrected;
|
||||
laplacian((1|aU),p) Gauss linear corrected;
|
||||
laplacian(diffusivity,cellMotionU) Gauss linear uncorrected;
|
||||
laplacian(nuEff,U) Gauss linear uncorrected;
|
||||
}
|
||||
|
|
|
@ -65,4 +65,10 @@ PIMPLE
|
|||
nNonOrthogonalCorrectors 0;
|
||||
}
|
||||
|
||||
relaxationFactors
|
||||
{
|
||||
U 1;
|
||||
UFinal 1;
|
||||
}
|
||||
|
||||
// ************************************************************************* //
|
||||
|
|
|
@ -41,14 +41,15 @@ laplacianSchemes
|
|||
|
||||
laplacian(1,p) Gauss linear limited 0.5;
|
||||
|
||||
laplacian((1|A(U)),p) Gauss linear limited 0.5;
|
||||
laplacian((1|aU),pcorr) Gauss linear limited 0.5;
|
||||
laplacian((1|aU),p) Gauss linear limited 0.5;
|
||||
}
|
||||
|
||||
interpolationSchemes
|
||||
{
|
||||
default linear;
|
||||
interpolate(HbyA) linear;
|
||||
interpolate(1|A) linear;
|
||||
interpolate(1|aU) linear;
|
||||
}
|
||||
|
||||
snGradSchemes
|
||||
|
|
|
@ -43,7 +43,7 @@ laplacianSchemes
|
|||
{
|
||||
default none;
|
||||
laplacian(nuEff,U) Gauss linear corrected;
|
||||
laplacian((1|A(U)),p) Gauss linear corrected;
|
||||
laplacian((1|aU),p) Gauss linear corrected;
|
||||
laplacian(DkEff,k) Gauss linear corrected;
|
||||
laplacian(DepsilonEff,epsilon) Gauss linear corrected;
|
||||
laplacian(DREff,R) Gauss linear corrected;
|
||||
|
|
|
@ -21,13 +21,7 @@ ddtSchemes
|
|||
|
||||
gradSchemes
|
||||
{
|
||||
default Gauss linear;
|
||||
grad(p) Gauss linear;
|
||||
grad(U) Gauss linear;
|
||||
// grad(U) cellLimited Gauss linear 1;
|
||||
|
||||
grad(k) Gauss linear;
|
||||
grad(omega) Gauss linear;
|
||||
default cellLimited leastSquares 1;
|
||||
}
|
||||
|
||||
divSchemes
|
||||
|
@ -36,8 +30,7 @@ divSchemes
|
|||
div(phi,U) Gauss linearUpwindV Gauss linear;
|
||||
div(phi,k) Gauss upwind;
|
||||
div(phi,omega) Gauss upwind;
|
||||
div((nuEff*dev(T(grad(U))))) Gauss linear;
|
||||
div((nuEff*dev(T(grad(U))))) Gauss linear;
|
||||
div((nuEff*dev(grad(U).T()))) Gauss linear;
|
||||
}
|
||||
|
||||
laplacianSchemes
|
||||
|
|
|
@ -57,7 +57,7 @@ solvers
|
|||
|
||||
SIMPLE
|
||||
{
|
||||
nNonOrthogonalCorrectors 0;
|
||||
nNonOrthogonalCorrectors 2;
|
||||
}
|
||||
|
||||
relaxationFactors
|
||||
|
|
|
@ -25,8 +25,8 @@ solvers
|
|||
};
|
||||
Urel
|
||||
{
|
||||
solver PBiCG;
|
||||
preconditioner DILU;
|
||||
solver smoothSolver;
|
||||
smoother GaussSeidel;
|
||||
tolerance 1e-05;
|
||||
relTol 0.1;
|
||||
};
|
||||
|
@ -75,7 +75,7 @@ SIMPLE
|
|||
relaxationFactors
|
||||
{
|
||||
p 0.3;
|
||||
Urel 0.7;
|
||||
Urel 0.5;
|
||||
k 0.7;
|
||||
epsilon 0.7;
|
||||
omega 0.7;
|
||||
|
|
Reference in a new issue