mirror of
https://develop.openfoam.com/Development/openfoam.git
synced 2025-11-28 03:28:01 +00:00
Retire settlingFoam: replaced by driftFluxFoam
This commit is contained in:
Henry
@ -1,3 +0,0 @@
|
||||
settlingFoam.C
|
||||
|
||||
EXE = $(FOAM_APPBIN)/settlingFoam
|
||||
@ -1,5 +0,0 @@
|
||||
EXE_INC = \
|
||||
-I$(LIB_SRC)/finiteVolume/lnInclude
|
||||
|
||||
EXE_LIBS = \
|
||||
-lfiniteVolume
|
||||
@ -1,32 +0,0 @@
|
||||
// Solve the Momentum equation
|
||||
|
||||
fvVectorMatrix UEqn
|
||||
(
|
||||
fvm::ddt(rho, U)
|
||||
+ fvm::div(phi, U)
|
||||
+ fvc::div
|
||||
(
|
||||
(Alpha/(scalar(1.001) - Alpha))*(sqr(rhoc)/rho)*Vdj*Vdj,
|
||||
"div(phiVdj,Vdj)"
|
||||
)
|
||||
- fvm::laplacian(muEff, U, "laplacian(muEff,U)")
|
||||
- (fvc::grad(U) & fvc::grad(muEff))
|
||||
);
|
||||
|
||||
UEqn.relax();
|
||||
|
||||
if (pimple.momentumPredictor())
|
||||
{
|
||||
solve
|
||||
(
|
||||
UEqn
|
||||
==
|
||||
fvc::reconstruct
|
||||
(
|
||||
(
|
||||
- ghf*fvc::snGrad(rho)
|
||||
- fvc::snGrad(p_rgh)
|
||||
)*mesh.magSf()
|
||||
)
|
||||
);
|
||||
}
|
||||
@ -1,34 +0,0 @@
|
||||
{
|
||||
surfaceScalarField phiAlpha
|
||||
(
|
||||
IOobject
|
||||
(
|
||||
"phiAlpha",
|
||||
runTime.timeName(),
|
||||
mesh
|
||||
),
|
||||
phi + rhoc*(mesh.Sf() & fvc::interpolate(Vdj))
|
||||
);
|
||||
|
||||
fvScalarMatrix AlphaEqn
|
||||
(
|
||||
fvm::ddt(rho, Alpha)
|
||||
+ fvm::div(phiAlpha, Alpha)
|
||||
- fvm::laplacian(mut, Alpha)
|
||||
);
|
||||
|
||||
AlphaEqn.relax();
|
||||
AlphaEqn.solve();
|
||||
|
||||
Info<< "Solid phase fraction = "
|
||||
<< Alpha.weightedAverage(mesh.V()).value()
|
||||
<< " Min(Alpha) = " << min(Alpha).value()
|
||||
<< " Max(Alpha) = " << max(Alpha).value()
|
||||
<< endl;
|
||||
|
||||
Alpha.min(1.0);
|
||||
Alpha.max(0.0);
|
||||
|
||||
rho == rhoc/(scalar(1) + (rhoc/rhod - 1.0)*Alpha);
|
||||
alpha == rho*Alpha/rhod;
|
||||
}
|
||||
@ -1,20 +0,0 @@
|
||||
if (VdjModel == "general")
|
||||
{
|
||||
Vdj = V0*
|
||||
(
|
||||
exp(-a*max(alpha - alphaMin, scalar(0)))
|
||||
- exp(-a1*max(alpha - alphaMin, scalar(0)))
|
||||
);
|
||||
}
|
||||
else if (VdjModel == "simple")
|
||||
{
|
||||
Vdj = V0*pow(10.0, -a*alpha);
|
||||
}
|
||||
else
|
||||
{
|
||||
FatalErrorIn(args.executable())
|
||||
<< "Unknown VdjModel : " << VdjModel
|
||||
<< abort(FatalError);
|
||||
}
|
||||
|
||||
Vdj.correctBoundaryConditions();
|
||||
@ -1,21 +0,0 @@
|
||||
scalar sumLocalContErr =
|
||||
runTime.deltaTValue()*
|
||||
mag
|
||||
(
|
||||
fvc::ddt(rho)
|
||||
+ fvc::div(phi)
|
||||
)().weightedAverage(rho*mesh.V()).value();
|
||||
|
||||
scalar globalContErr =
|
||||
runTime.deltaTValue()*
|
||||
(
|
||||
fvc::ddt(rho)
|
||||
+ fvc::div(phi)
|
||||
)().weightedAverage(rho*mesh.V()).value();
|
||||
|
||||
cumulativeContErr += globalContErr;
|
||||
|
||||
Info<< "time step continuity errors : sum local = " << sumLocalContErr
|
||||
<< ", global = " << globalContErr
|
||||
<< ", cumulative = " << cumulativeContErr
|
||||
<< endl;
|
||||
@ -1,39 +0,0 @@
|
||||
{
|
||||
mul = muc +
|
||||
plasticViscosity
|
||||
(
|
||||
plasticViscosityCoeff,
|
||||
plasticViscosityExponent,
|
||||
alpha
|
||||
);
|
||||
|
||||
if (BinghamPlastic)
|
||||
{
|
||||
volScalarField tauy = yieldStress
|
||||
(
|
||||
yieldStressCoeff,
|
||||
yieldStressExponent,
|
||||
yieldStressOffset,
|
||||
alpha
|
||||
);
|
||||
|
||||
mul =
|
||||
tauy/
|
||||
(
|
||||
mag(fvc::grad(U))
|
||||
+ 1.0e-4*
|
||||
(
|
||||
tauy
|
||||
+ dimensionedScalar
|
||||
(
|
||||
"deltaTauy",
|
||||
tauy.dimensions(),
|
||||
1.0e-15
|
||||
)
|
||||
)/mul
|
||||
)
|
||||
+ mul;
|
||||
}
|
||||
|
||||
mul = min(mul, muMax);
|
||||
}
|
||||
@ -1,383 +0,0 @@
|
||||
Info<< "Reading field p_rgh\n" << endl;
|
||||
volScalarField p_rgh
|
||||
(
|
||||
IOobject
|
||||
(
|
||||
"p_rgh",
|
||||
runTime.timeName(),
|
||||
mesh,
|
||||
IOobject::MUST_READ,
|
||||
IOobject::AUTO_WRITE
|
||||
),
|
||||
mesh
|
||||
);
|
||||
|
||||
Info<< "Reading field Alpha\n" << endl;
|
||||
volScalarField Alpha
|
||||
(
|
||||
IOobject
|
||||
(
|
||||
"Alpha",
|
||||
runTime.timeName(),
|
||||
mesh,
|
||||
IOobject::MUST_READ,
|
||||
IOobject::AUTO_WRITE
|
||||
),
|
||||
mesh
|
||||
);
|
||||
|
||||
Info<< "Reading field U\n" << endl;
|
||||
volVectorField U
|
||||
(
|
||||
IOobject
|
||||
(
|
||||
"U",
|
||||
runTime.timeName(),
|
||||
mesh,
|
||||
IOobject::MUST_READ,
|
||||
IOobject::AUTO_WRITE
|
||||
),
|
||||
mesh
|
||||
);
|
||||
|
||||
|
||||
Info<< "Reading transportProperties\n" << endl;
|
||||
|
||||
IOdictionary transportProperties
|
||||
(
|
||||
IOobject
|
||||
(
|
||||
"transportProperties",
|
||||
runTime.constant(),
|
||||
mesh,
|
||||
IOobject::MUST_READ_IF_MODIFIED,
|
||||
IOobject::NO_WRITE
|
||||
)
|
||||
);
|
||||
|
||||
|
||||
dimensionedScalar rhoc(transportProperties.lookup("rhoc"));
|
||||
|
||||
dimensionedScalar rhod(transportProperties.lookup("rhod"));
|
||||
|
||||
dimensionedScalar muc(transportProperties.lookup("muc"));
|
||||
dimensionedScalar muMax(transportProperties.lookup("muMax"));
|
||||
|
||||
dimensionedScalar plasticViscosityCoeff
|
||||
(
|
||||
transportProperties.lookup("plasticViscosityCoeff")
|
||||
);
|
||||
|
||||
dimensionedScalar plasticViscosityExponent
|
||||
(
|
||||
transportProperties.lookup("plasticViscosityExponent")
|
||||
);
|
||||
|
||||
dimensionedScalar yieldStressCoeff
|
||||
(
|
||||
transportProperties.lookup("yieldStressCoeff")
|
||||
);
|
||||
|
||||
dimensionedScalar yieldStressExponent
|
||||
(
|
||||
transportProperties.lookup("yieldStressExponent")
|
||||
);
|
||||
|
||||
dimensionedScalar yieldStressOffset
|
||||
(
|
||||
transportProperties.lookup("yieldStressOffset")
|
||||
);
|
||||
|
||||
Switch BinghamPlastic(transportProperties.lookup("BinghamPlastic"));
|
||||
|
||||
volScalarField rho
|
||||
(
|
||||
IOobject
|
||||
(
|
||||
"rho",
|
||||
runTime.timeName(),
|
||||
mesh,
|
||||
IOobject::NO_READ,
|
||||
IOobject::NO_WRITE
|
||||
),
|
||||
rhoc/(scalar(1) + (rhoc/rhod - 1.0)*Alpha)
|
||||
);
|
||||
|
||||
volScalarField alpha
|
||||
(
|
||||
IOobject
|
||||
(
|
||||
"alpha",
|
||||
runTime.timeName(),
|
||||
mesh,
|
||||
IOobject::NO_READ,
|
||||
IOobject::AUTO_WRITE
|
||||
),
|
||||
rho*Alpha/rhod
|
||||
);
|
||||
|
||||
#include "compressibleCreatePhi.H"
|
||||
|
||||
|
||||
Info<< "Calculating field mul\n" << endl;
|
||||
volScalarField mul
|
||||
(
|
||||
IOobject
|
||||
(
|
||||
"mul",
|
||||
runTime.timeName(),
|
||||
mesh,
|
||||
IOobject::NO_READ,
|
||||
IOobject::AUTO_WRITE
|
||||
),
|
||||
muc
|
||||
+ plasticViscosity
|
||||
(
|
||||
plasticViscosityCoeff,
|
||||
plasticViscosityExponent,
|
||||
Alpha
|
||||
)
|
||||
);
|
||||
|
||||
|
||||
Info<< "Initialising field Vdj\n" << endl;
|
||||
volVectorField Vdj
|
||||
(
|
||||
IOobject
|
||||
(
|
||||
"Vdj",
|
||||
runTime.timeName(),
|
||||
mesh,
|
||||
IOobject::NO_READ,
|
||||
IOobject::AUTO_WRITE
|
||||
),
|
||||
mesh,
|
||||
dimensionedVector("0.0", U.dimensions(), vector::zero),
|
||||
U.boundaryField().types()
|
||||
);
|
||||
|
||||
|
||||
Info<< "Selecting Drift-Flux model " << endl;
|
||||
|
||||
const word VdjModel(transportProperties.lookup("VdjModel"));
|
||||
|
||||
Info<< tab << VdjModel << " selected\n" << endl;
|
||||
|
||||
const dictionary& VdjModelCoeffs
|
||||
(
|
||||
transportProperties.subDict(VdjModel + "Coeffs")
|
||||
);
|
||||
|
||||
dimensionedVector V0(VdjModelCoeffs.lookup("V0"));
|
||||
|
||||
dimensionedScalar a(VdjModelCoeffs.lookup("a"));
|
||||
|
||||
dimensionedScalar a1(VdjModelCoeffs.lookup("a1"));
|
||||
|
||||
dimensionedScalar alphaMin(VdjModelCoeffs.lookup("alphaMin"));
|
||||
|
||||
|
||||
IOdictionary RASProperties
|
||||
(
|
||||
IOobject
|
||||
(
|
||||
"RASProperties",
|
||||
runTime.constant(),
|
||||
mesh,
|
||||
IOobject::MUST_READ_IF_MODIFIED,
|
||||
IOobject::NO_WRITE
|
||||
)
|
||||
);
|
||||
|
||||
|
||||
Switch turbulence(RASProperties.lookup("turbulence"));
|
||||
|
||||
dictionary kEpsilonDict(RASProperties.subDictPtr("kEpsilonCoeffs"));
|
||||
|
||||
dimensionedScalar Cmu
|
||||
(
|
||||
dimensionedScalar::lookupOrAddToDict
|
||||
(
|
||||
"Cmu",
|
||||
kEpsilonDict,
|
||||
0.09
|
||||
)
|
||||
);
|
||||
|
||||
dimensionedScalar C1
|
||||
(
|
||||
dimensionedScalar::lookupOrAddToDict
|
||||
(
|
||||
"C1",
|
||||
kEpsilonDict,
|
||||
1.44
|
||||
)
|
||||
);
|
||||
|
||||
dimensionedScalar C2
|
||||
(
|
||||
dimensionedScalar::lookupOrAddToDict
|
||||
(
|
||||
"C2",
|
||||
kEpsilonDict,
|
||||
1.92
|
||||
)
|
||||
);
|
||||
|
||||
dimensionedScalar C3
|
||||
(
|
||||
dimensionedScalar::lookupOrAddToDict
|
||||
(
|
||||
"C3",
|
||||
kEpsilonDict,
|
||||
0.85
|
||||
)
|
||||
);
|
||||
|
||||
dimensionedScalar sigmak
|
||||
(
|
||||
dimensionedScalar::lookupOrAddToDict
|
||||
(
|
||||
"sigmak",
|
||||
kEpsilonDict,
|
||||
1.0
|
||||
)
|
||||
);
|
||||
|
||||
dimensionedScalar sigmaEps
|
||||
(
|
||||
dimensionedScalar::lookupOrAddToDict
|
||||
(
|
||||
"sigmaEps",
|
||||
kEpsilonDict,
|
||||
1.3
|
||||
)
|
||||
);
|
||||
|
||||
dictionary wallFunctionDict(RASProperties.subDictPtr("wallFunctionCoeffs"));
|
||||
|
||||
dimensionedScalar kappa
|
||||
(
|
||||
dimensionedScalar::lookupOrAddToDict
|
||||
(
|
||||
"kappa",
|
||||
wallFunctionDict,
|
||||
0.41
|
||||
)
|
||||
);
|
||||
|
||||
dimensionedScalar E
|
||||
(
|
||||
dimensionedScalar::lookupOrAddToDict
|
||||
(
|
||||
"E",
|
||||
wallFunctionDict,
|
||||
9.8
|
||||
)
|
||||
);
|
||||
|
||||
if (RASProperties.lookupOrDefault("printCoeffs", false))
|
||||
{
|
||||
Info<< "kEpsilonCoeffs" << kEpsilonDict << nl
|
||||
<< "wallFunctionCoeffs" << wallFunctionDict << endl;
|
||||
}
|
||||
|
||||
|
||||
nearWallDist y(mesh);
|
||||
|
||||
Info<< "Reading field k\n" << endl;
|
||||
volScalarField k
|
||||
(
|
||||
IOobject
|
||||
(
|
||||
"k",
|
||||
runTime.timeName(),
|
||||
mesh,
|
||||
IOobject::MUST_READ,
|
||||
IOobject::AUTO_WRITE
|
||||
),
|
||||
mesh
|
||||
);
|
||||
|
||||
Info<< "Reading field epsilon\n" << endl;
|
||||
volScalarField epsilon
|
||||
(
|
||||
IOobject
|
||||
(
|
||||
"epsilon",
|
||||
runTime.timeName(),
|
||||
mesh,
|
||||
IOobject::MUST_READ,
|
||||
IOobject::AUTO_WRITE
|
||||
),
|
||||
mesh
|
||||
);
|
||||
|
||||
Info<< "Calculating field mut\n" << endl;
|
||||
volScalarField mut
|
||||
(
|
||||
IOobject
|
||||
(
|
||||
"mut",
|
||||
runTime.timeName(),
|
||||
mesh,
|
||||
IOobject::NO_READ,
|
||||
IOobject::AUTO_WRITE
|
||||
),
|
||||
Cmu*rho*sqr(k)/epsilon
|
||||
);
|
||||
|
||||
|
||||
Info<< "Calculating field muEff\n" << endl;
|
||||
volScalarField muEff
|
||||
(
|
||||
IOobject
|
||||
(
|
||||
"muEff",
|
||||
runTime.timeName(),
|
||||
mesh,
|
||||
IOobject::NO_READ,
|
||||
IOobject::AUTO_WRITE
|
||||
),
|
||||
mut + mul
|
||||
);
|
||||
|
||||
|
||||
Info<< "Calculating field (g.h)f\n" << endl;
|
||||
volScalarField gh("gh", g & mesh.C());
|
||||
surfaceScalarField ghf("gh", g & mesh.Cf());
|
||||
|
||||
volScalarField p
|
||||
(
|
||||
IOobject
|
||||
(
|
||||
"p",
|
||||
runTime.timeName(),
|
||||
mesh,
|
||||
IOobject::NO_READ,
|
||||
IOobject::AUTO_WRITE
|
||||
),
|
||||
p_rgh + rho*gh
|
||||
);
|
||||
|
||||
label pRefCell = 0;
|
||||
scalar pRefValue = 0.0;
|
||||
setRefCell
|
||||
(
|
||||
p,
|
||||
p_rgh,
|
||||
mesh.solutionDict().subDict("PIMPLE"),
|
||||
pRefCell,
|
||||
pRefValue
|
||||
);
|
||||
|
||||
if (p_rgh.needReference())
|
||||
{
|
||||
p += dimensionedScalar
|
||||
(
|
||||
"p",
|
||||
p.dimensions(),
|
||||
pRefValue - getRefCellValue(p, pRefCell)
|
||||
);
|
||||
p_rgh = p - rho*gh;
|
||||
}
|
||||
@ -1,78 +0,0 @@
|
||||
if (turbulence)
|
||||
{
|
||||
if (mesh.changing())
|
||||
{
|
||||
y.correct();
|
||||
}
|
||||
|
||||
dimensionedScalar k0("k0", k.dimensions(), 0);
|
||||
dimensionedScalar kMin("kMin", k.dimensions(), SMALL);
|
||||
dimensionedScalar epsilon0("epsilon0", epsilon.dimensions(), 0);
|
||||
dimensionedScalar epsilonMin("epsilonMin", epsilon.dimensions(), SMALL);
|
||||
|
||||
volScalarField divU(fvc::div(phi/fvc::interpolate(rho)));
|
||||
|
||||
tmp<volTensorField> tgradU = fvc::grad(U);
|
||||
volScalarField G(2*mut*(tgradU() && dev(symm(tgradU()))));
|
||||
tgradU.clear();
|
||||
|
||||
volScalarField Gcoef
|
||||
(
|
||||
Cmu*k/sigmak*(g & fvc::grad(rho))/(epsilon + epsilonMin)
|
||||
);
|
||||
|
||||
#include "wallFunctions.H"
|
||||
|
||||
// Dissipation equation
|
||||
fvScalarMatrix epsEqn
|
||||
(
|
||||
fvm::ddt(rho, epsilon)
|
||||
+ fvm::div(phi, epsilon)
|
||||
- fvm::laplacian
|
||||
(
|
||||
mut/sigmaEps + muc, epsilon,
|
||||
"laplacian(DepsilonEff,epsilon)"
|
||||
)
|
||||
==
|
||||
C1*G*epsilon/(k + kMin)
|
||||
- fvm::SuSp(C1*(1.0 - C3)*Gcoef + (2.0/3.0*C1)*rho*divU, epsilon)
|
||||
- fvm::Sp(C2*rho*epsilon/(k + kMin), epsilon)
|
||||
);
|
||||
|
||||
#include "wallDissipation.H"
|
||||
|
||||
epsEqn.relax();
|
||||
epsEqn.solve();
|
||||
|
||||
bound(epsilon, epsilon0);
|
||||
|
||||
|
||||
// Turbulent kinetic energy equation
|
||||
fvScalarMatrix kEqn
|
||||
(
|
||||
fvm::ddt(rho, k)
|
||||
+ fvm::div(phi, k)
|
||||
- fvm::laplacian
|
||||
(
|
||||
mut/sigmak + muc, k,
|
||||
"laplacian(DkEff,k)"
|
||||
)
|
||||
==
|
||||
G
|
||||
- fvm::SuSp(Gcoef + 2.0/3.0*rho*divU, k)
|
||||
- fvm::Sp(rho*epsilon/(k + kMin), k)
|
||||
);
|
||||
|
||||
kEqn.relax();
|
||||
kEqn.solve();
|
||||
|
||||
bound(k, k0);
|
||||
|
||||
|
||||
//- Re-calculate viscosity
|
||||
mut = rho*Cmu*sqr(k)/(epsilon + epsilonMin);
|
||||
|
||||
#include "wallViscosity.H"
|
||||
}
|
||||
|
||||
muEff = mut + mul;
|
||||
@ -1,70 +0,0 @@
|
||||
{
|
||||
volScalarField rAU("rAU", 1.0/UEqn.A());
|
||||
surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
|
||||
|
||||
volVectorField HbyA("HbyA", U);
|
||||
HbyA = rAU*UEqn.H();
|
||||
|
||||
surfaceScalarField phiHbyA
|
||||
(
|
||||
"phiHbyA",
|
||||
(
|
||||
(fvc::interpolate(rho*HbyA) & mesh.Sf())
|
||||
+ rhorAUf*fvc::ddtCorr(rho, U, phi)
|
||||
)
|
||||
);
|
||||
|
||||
surfaceScalarField phig
|
||||
(
|
||||
- ghf*fvc::snGrad(rho)*rhorAUf*mesh.magSf()
|
||||
);
|
||||
|
||||
phiHbyA += phig;
|
||||
|
||||
// Update the fixedFluxPressure BCs to ensure flux consistency
|
||||
setSnGrad<fixedFluxPressureFvPatchScalarField>
|
||||
(
|
||||
p_rgh.boundaryField(),
|
||||
(
|
||||
phiHbyA.boundaryField()
|
||||
- (mesh.Sf().boundaryField() & U.boundaryField())
|
||||
*rho.boundaryField()
|
||||
)/(mesh.magSf().boundaryField()*rhorAUf.boundaryField())
|
||||
);
|
||||
|
||||
while (pimple.correctNonOrthogonal())
|
||||
{
|
||||
fvScalarMatrix p_rghEqn
|
||||
(
|
||||
fvm::laplacian(rhorAUf, p_rgh) == fvc::ddt(rho) + fvc::div(phiHbyA)
|
||||
);
|
||||
|
||||
p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
|
||||
|
||||
p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
|
||||
|
||||
if (pimple.finalNonOrthogonalIter())
|
||||
{
|
||||
phi = phiHbyA - p_rghEqn.flux();
|
||||
|
||||
U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rhorAUf);
|
||||
U.correctBoundaryConditions();
|
||||
}
|
||||
}
|
||||
|
||||
p == p_rgh + rho*gh;
|
||||
|
||||
if (p_rgh.needReference())
|
||||
{
|
||||
p += dimensionedScalar
|
||||
(
|
||||
"p",
|
||||
p.dimensions(),
|
||||
pRefValue - getRefCellValue(p, pRefCell)
|
||||
);
|
||||
p_rgh = p - rho*gh;
|
||||
}
|
||||
|
||||
#include "rhoEqn.H"
|
||||
#include "compressibleContinuityErrs.H"
|
||||
}
|
||||
@ -1,17 +0,0 @@
|
||||
volScalarField plasticViscosity
|
||||
(
|
||||
const dimensionedScalar& plasticViscosityCoeff,
|
||||
const dimensionedScalar& plasticViscosityExponent,
|
||||
const volScalarField& alpha
|
||||
)
|
||||
{
|
||||
tmp<volScalarField> tfld
|
||||
(
|
||||
plasticViscosityCoeff*
|
||||
(
|
||||
pow(10.0, plasticViscosityExponent*alpha + SMALL) - scalar(1)
|
||||
)
|
||||
);
|
||||
|
||||
return tfld();
|
||||
}
|
||||
@ -1,109 +0,0 @@
|
||||
/*---------------------------------------------------------------------------*\
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration |
|
||||
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
License
|
||||
This file is part of OpenFOAM.
|
||||
|
||||
OpenFOAM is free software: you can redistribute it and/or modify it
|
||||
under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation, either version 3 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
|
||||
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
|
||||
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
|
||||
for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License
|
||||
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
|
||||
|
||||
Application
|
||||
settlingFoam
|
||||
|
||||
Description
|
||||
Solver for 2 incompressible fluids for simulating the settling of the
|
||||
dispersed phase.
|
||||
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
#include "fvCFD.H"
|
||||
#include "nearWallDist.H"
|
||||
#include "wallFvPatch.H"
|
||||
#include "bound.H"
|
||||
#include "Switch.H"
|
||||
#include "plasticViscosity.H"
|
||||
#include "yieldStress.H"
|
||||
#include "pimpleControl.H"
|
||||
#include "fixedFluxPressureFvPatchScalarField.H"
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
int main(int argc, char *argv[])
|
||||
{
|
||||
#include "setRootCase.H"
|
||||
|
||||
#include "createTime.H"
|
||||
#include "createMesh.H"
|
||||
#include "readGravitationalAcceleration.H"
|
||||
#include "createFields.H"
|
||||
#include "initContinuityErrs.H"
|
||||
|
||||
pimpleControl pimple(mesh);
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
Info<< "\nStarting time loop\n" << endl;
|
||||
|
||||
while (runTime.run())
|
||||
{
|
||||
#include "readTimeControls.H"
|
||||
#include "compressibleCourantNo.H"
|
||||
#include "setDeltaT.H"
|
||||
|
||||
runTime++;
|
||||
|
||||
Info<< "Time = " << runTime.timeName() << nl << endl;
|
||||
|
||||
#include "rhoEqn.H"
|
||||
|
||||
// --- Pressure-velocity PIMPLE corrector loop
|
||||
while (pimple.loop())
|
||||
{
|
||||
#include "calcVdj.H"
|
||||
|
||||
#include "UEqn.H"
|
||||
|
||||
#include "alphaEqn.H"
|
||||
|
||||
#include "correctViscosity.H"
|
||||
|
||||
// --- Pressure corrector loop
|
||||
while (pimple.correct())
|
||||
{
|
||||
#include "pEqn.H"
|
||||
}
|
||||
|
||||
if (pimple.turbCorr())
|
||||
{
|
||||
#include "kEpsilon.H"
|
||||
}
|
||||
}
|
||||
|
||||
runTime.write();
|
||||
|
||||
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
|
||||
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
|
||||
<< nl << endl;
|
||||
}
|
||||
|
||||
Info<< "End\n" << endl;
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
||||
|
||||
// ************************************************************************* //
|
||||
@ -1,50 +0,0 @@
|
||||
/*---------------------------------------------------------------------------*\
|
||||
========= |
|
||||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||||
\\ / O peration |
|
||||
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
License
|
||||
This file is part of OpenFOAM.
|
||||
|
||||
OpenFOAM is free software: you can redistribute it and/or modify it
|
||||
under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation, either version 3 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
|
||||
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
|
||||
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
|
||||
for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License
|
||||
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
|
||||
|
||||
Global
|
||||
wallDissipation
|
||||
|
||||
Description
|
||||
Set wall dissipation in the epsilon matrix
|
||||
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
{
|
||||
const fvPatchList& patches = mesh.boundary();
|
||||
|
||||
forAll(patches, patchi)
|
||||
{
|
||||
const fvPatch& p = patches[patchi];
|
||||
|
||||
if (isA<wallFvPatch>(p))
|
||||
{
|
||||
epsEqn.setValues
|
||||
(
|
||||
p.faceCells(),
|
||||
epsilon.boundaryField()[patchi].patchInternalField()
|
||||
);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// ************************************************************************* //
|
||||
@ -1,85 +0,0 @@
|
||||
{
|
||||
labelList cellBoundaryFaceCount(epsilon.size(), 0);
|
||||
|
||||
const scalar Cmu25 = ::pow(Cmu.value(), 0.25);
|
||||
const scalar Cmu75 = ::pow(Cmu.value(), 0.75);
|
||||
const scalar kappa_ = kappa.value();
|
||||
const scalar muc_ = muc.value();
|
||||
|
||||
const fvPatchList& patches = mesh.boundary();
|
||||
|
||||
//- Initialise the near-wall P field to zero
|
||||
forAll(patches, patchi)
|
||||
{
|
||||
const fvPatch& curPatch = patches[patchi];
|
||||
|
||||
if (isA<wallFvPatch>(curPatch))
|
||||
{
|
||||
forAll(curPatch, facei)
|
||||
{
|
||||
label faceCelli = curPatch.faceCells()[facei];
|
||||
|
||||
epsilon[faceCelli] = 0.0;
|
||||
G[faceCelli] = 0.0;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
//- Accumulate the wall face contributions to epsilon and G
|
||||
// Increment cellBoundaryFaceCount for each face for averaging
|
||||
forAll(patches, patchi)
|
||||
{
|
||||
const fvPatch& curPatch = patches[patchi];
|
||||
|
||||
if (isA<wallFvPatch>(curPatch))
|
||||
{
|
||||
const scalarField& mutw = mut.boundaryField()[patchi];
|
||||
|
||||
scalarField magFaceGradU
|
||||
(
|
||||
mag(U.boundaryField()[patchi].snGrad())
|
||||
);
|
||||
|
||||
forAll(curPatch, facei)
|
||||
{
|
||||
label faceCelli = curPatch.faceCells()[facei];
|
||||
|
||||
// For corner cells (with two boundary or more faces),
|
||||
// epsilon and G in the near-wall cell are calculated
|
||||
// as an average
|
||||
|
||||
cellBoundaryFaceCount[faceCelli]++;
|
||||
|
||||
epsilon[faceCelli] +=
|
||||
Cmu75*::pow(k[faceCelli], 1.5)
|
||||
/(kappa_*y[patchi][facei]);
|
||||
|
||||
G[faceCelli] +=
|
||||
(mutw[facei] + muc_)
|
||||
*magFaceGradU[facei]
|
||||
*Cmu25*::sqrt(k[faceCelli])
|
||||
/(kappa_*y[patchi][facei]);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// perform the averaging
|
||||
|
||||
forAll(patches, patchi)
|
||||
{
|
||||
const fvPatch& curPatch = patches[patchi];
|
||||
|
||||
if (isA<wallFvPatch>(curPatch))
|
||||
{
|
||||
forAll(curPatch, facei)
|
||||
{
|
||||
label faceCelli = curPatch.faceCells()[facei];
|
||||
|
||||
epsilon[faceCelli] /= cellBoundaryFaceCount[faceCelli];
|
||||
G[faceCelli] /= cellBoundaryFaceCount[faceCelli];
|
||||
cellBoundaryFaceCount[faceCelli] = 1;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
@ -1,38 +0,0 @@
|
||||
{
|
||||
const scalar Cmu25 = ::pow(Cmu.value(), 0.25);
|
||||
const scalar kappa_ = kappa.value();
|
||||
const scalar E_ = E.value();
|
||||
const scalar muc_ = muc.value();
|
||||
const scalar nuc_ = muc_/rhoc.value();
|
||||
|
||||
const fvPatchList& patches = mesh.boundary();
|
||||
|
||||
forAll(patches, patchi)
|
||||
{
|
||||
const fvPatch& curPatch = patches[patchi];
|
||||
|
||||
if (isA<wallFvPatch>(curPatch))
|
||||
{
|
||||
scalarField& mutw = mut.boundaryField()[patchi];
|
||||
|
||||
forAll(curPatch, facei)
|
||||
{
|
||||
label faceCelli = curPatch.faceCells()[facei];
|
||||
|
||||
scalar yPlus =
|
||||
Cmu25*y[patchi][facei]*::sqrt(k[faceCelli])
|
||||
/nuc_;
|
||||
|
||||
if (yPlus > 11.6)
|
||||
{
|
||||
mutw[facei] =
|
||||
muc_*(yPlus*kappa_/::log(E_*yPlus) - 1);
|
||||
}
|
||||
else
|
||||
{
|
||||
mutw[facei] = 0.0;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
@ -1,19 +0,0 @@
|
||||
volScalarField yieldStress
|
||||
(
|
||||
const dimensionedScalar& yieldStressCoeff,
|
||||
const dimensionedScalar& yieldStressExponent,
|
||||
const dimensionedScalar& yieldStressOffset,
|
||||
const volScalarField& alpha
|
||||
)
|
||||
{
|
||||
tmp<volScalarField> tfld
|
||||
(
|
||||
yieldStressCoeff*
|
||||
(
|
||||
pow(10.0, yieldStressExponent*(alpha + yieldStressOffset))
|
||||
- pow(10.0, yieldStressExponent*yieldStressOffset)
|
||||
)
|
||||
);
|
||||
|
||||
return tfld();
|
||||
}
|
||||
Reference in New Issue
Block a user