zhyang-dev/openfoam
1
0
Fork 0
mirror of https://develop.openfoam.com/Development/openfoam.git synced 2025-11-28 03:28:01 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

Merge branch 'master' of /home/noisy3/OpenFOAM/OpenFOAM-dev

Browse Source
This commit is contained in:
mattijs
2010-09-07 14:11:45 +01:00
parent 97dff2650e 08a9c9da8f
commit d400262c70
21 changed files with 659 additions and 108 deletions
Download Patch File Download Diff File Expand all files Collapse all files

12
applications/solvers/multiphase/twoPhaseEulerFoam/UEqns.H
View File

@ -4,12 +4,24 @@ fvVectorMatrix UbEqn(Ub, Ub.dimensions()*dimVol/dimTime);
{
{
volTensorField gradUaT = fvc::grad(Ua)().T();
if (kineticTheory.on())
{
kineticTheory.solve(gradUaT);
nuEffa = kineticTheory.mua()/rhoa;
}
else // If not using kinetic theory is using Ct model
{
nuEffa = sqr(Ct)*nutb + nua;
}
volTensorField Rca
(
"Rca",
((2.0/3.0)*I)*(sqr(Ct)*k + nuEffa*tr(gradUaT)) - nuEffa*gradUaT
);
if (kineticTheory.on())
{
Rca -= ((kineticTheory.lambda()/rhoa)*tr(gradUaT))*tensor(I);

62
applications/solvers/multiphase/twoPhaseEulerFoam/kEpsilon.H Normal file
View File

@ -0,0 +1,62 @@
if(turbulence)
{
if (mesh.changing())
{
y.correct();
}
tmp<volTensorField> tgradUb = fvc::grad(Ub);
volScalarField G = 2*nutb*(tgradUb() && dev(symm(tgradUb())));
tgradUb.clear();
#include "wallFunctions.H"
// Dissipation equation
fvScalarMatrix epsEqn
(
fvm::ddt(beta, epsilon)
+ fvm::div(phib, epsilon)
- fvm::laplacian
(
alphaEps*nuEffb, epsilon,
"laplacian(DepsilonEff,epsilon)"
)
==
C1*beta*G*epsilon/k
- fvm::Sp(C2*beta*epsilon/k, epsilon)
);
#include "wallDissipation.H"
epsEqn.relax();
epsEqn.solve();
epsilon.max(dimensionedScalar("zero", epsilon.dimensions(), 1.0e-15));
// Turbulent kinetic energy equation
fvScalarMatrix kEqn
(
fvm::ddt(beta, k)
+ fvm::div(phib, k)
- fvm::laplacian
(
alphak*nuEffb, k,
"laplacian(DkEff,k)"
)
==
beta*G
- fvm::Sp(beta*epsilon/k, k)
);
kEqn.relax();
kEqn.solve();
k.max(dimensionedScalar("zero", k.dimensions(), 1.0e-8));
//- Re-calculate turbulence viscosity
nutb = Cmu*sqr(k)/epsilon;
#include "wallViscosity.H"
}
nuEffb = nutb + nub;

2
applications/solvers/multiphase/twoPhaseEulerFoam/kineticTheoryModels/frictionalStressModel/JohnsonJackson/JohnsonJacksonFrictionalStress.C
View File

@ -102,7 +102,7 @@ Foam::tmp<Foam::volScalarField> Foam::JohnsonJacksonFrictionalStress::muf
const volScalarField& alpha,
const dimensionedScalar& alphaMax,
const volScalarField& pf,
const volTensorField& D,
const volSymmTensorField& D,
const dimensionedScalar& phi
) const
{

2
applications/solvers/multiphase/twoPhaseEulerFoam/kineticTheoryModels/frictionalStressModel/JohnsonJackson/JohnsonJacksonFrictionalStress.H
View File

@ -93,7 +93,7 @@ public:
const volScalarField& alpha,
const dimensionedScalar& alphaMax,
const volScalarField& pf,
const volTensorField& D,
const volSymmTensorField& D,
const dimensionedScalar& phi
) const;
};

6
applications/solvers/multiphase/twoPhaseEulerFoam/kineticTheoryModels/frictionalStressModel/Schaeffer/SchaefferFrictionalStress.C
View File

@ -99,7 +99,7 @@ Foam::tmp<Foam::volScalarField> Foam::SchaefferFrictionalStress::muf
const volScalarField& alpha,
const dimensionedScalar& alphaMax,
const volScalarField& pf,
const volTensorField& D,
const volSymmTensorField& D,
const dimensionedScalar& phi
) const
{
@ -124,9 +124,9 @@ Foam::tmp<Foam::volScalarField> Foam::SchaefferFrictionalStress::muf
volScalarField& muff = tmuf();
forAll(D, celli)
forAll (D, celli)
{
if (alpha[celli] > alphaMax.value()-5e-2)
if (alpha[celli] > alphaMax.value() - 5e-2)
{
muff[celli] =
0.5*pf[celli]*sin(phi.value())

2
applications/solvers/multiphase/twoPhaseEulerFoam/kineticTheoryModels/frictionalStressModel/Schaeffer/SchaefferFrictionalStress.H
View File

@ -93,7 +93,7 @@ public:
const volScalarField& alpha,
const dimensionedScalar& alphaMax,
const volScalarField& pf,
const volTensorField& D,
const volSymmTensorField& D,
const dimensionedScalar& phi
) const;
};

4
applications/solvers/multiphase/twoPhaseEulerFoam/kineticTheoryModels/frictionalStressModel/frictionalStressModel/frictionalStressModel.H
View File

@ -48,7 +48,7 @@ namespace Foam
class frictionalStressModel
{
// Private Member Functions
// Private member functions
//- Disallow default bitwise copy construct
frictionalStressModel(const frictionalStressModel&);
@ -127,7 +127,7 @@ public:
const volScalarField& alpha,
const dimensionedScalar& alphaMax,
const volScalarField& pf,
const volTensorField& D,
const volSymmTensorField& D,
const dimensionedScalar& phi
) const = 0;
};

83
applications/solvers/multiphase/twoPhaseEulerFoam/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.C
View File

@ -56,12 +56,13 @@ Foam::kineticTheoryModel::kineticTheoryModel
"kineticTheoryProperties",
Ua_.time().constant(),
Ua_.mesh(),
IOobject::MUST_READ_IF_MODIFIED,
IOobject::MUST_READ,
IOobject::NO_WRITE
)
),
kineticTheory_(kineticTheoryProperties_.lookup("kineticTheory")),
equilibrium_(kineticTheoryProperties_.lookup("equilibrium")),
viscosityModel_
(
kineticTheoryModels::viscosityModel::New
@ -192,24 +193,19 @@ Foam::kineticTheoryModel::~kineticTheoryModel()
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::kineticTheoryModel::solve()
void Foam::kineticTheoryModel::solve(const volTensorField& gradUat)
{
if (!kineticTheory_)
{
return;
}
word scheme("div(phi,Theta)");
volScalarField alpha = alpha_;
alpha.max(1.0e-6);
const scalar sqrtPi = sqrt(constant::mathematical::pi);
surfaceScalarField phi = 1.5*rhoa_*phia_*fvc::interpolate(alpha_);
volTensorField dU = fvc::grad(Ua_);
volTensorField dUT = dU.T();
volTensorField D = 0.5*(dU + dUT);
volTensorField dU = gradUat.T();//fvc::grad(Ua_);
volSymmTensorField D = symm(dU);
// NB, drag = K*alpha*beta,
// (the alpha and beta has been extracted from the drag function for
@ -220,45 +216,52 @@ void Foam::kineticTheoryModel::solve()
// Calculating the radial distribution function (solid volume fraction is
// limited close to the packing limit, but this needs improvements)
// The solution is higly unstable close to the packing limit.
gs0_ = radialModel_->g0(min(alpha, alphaMax_-1.0e-2), alphaMax_);
gs0_ = radialModel_->g0
(
min(max(alpha_, 1e-6), alphaMax_ - 0.01),
alphaMax_
);
// particle pressure - coefficient in front of Theta (Eq. 3.22, p. 45)
volScalarField PsCoeff =
granularPressureModel_->granularPressureCoeff(alpha_,gs0_,rhoa_,e_ );
volScalarField PsCoeff = granularPressureModel_->granularPressureCoeff
(
alpha_,
gs0_,
rhoa_,
e_
);
// 'thermal' conductivity (Table 3.3, p. 49)
kappa_ = conductivityModel_->kappa(alpha_, Theta_, gs0_, rhoa_, da_, e_);
// particle viscosity (Table 3.2, p.47)
mua_ = viscosityModel_->mua(alpha_, Theta_, gs0_, rhoa_, da_, e_);
dimensionedScalar Tsmall
(
"small",
dimensionSet(0,2,-2,0,0,0,0),
dimensionSet(0 , 2 ,-2 ,0 , 0, 0, 0),
1.0e-6
);
dimensionedScalar TsmallSqrt = sqrt(Tsmall);
volScalarField ThetaSqrt = sqrt(Theta_);
// 'thermal' conductivity (Table 3.3, p. 49)
kappa_ = conductivityModel_->kappa(alpha, Theta_, gs0_, rhoa_, da_, e_);
// particle viscosity (Table 3.2, p.47)
mua_ = viscosityModel_->mua(alpha, Theta_, gs0_, rhoa_, da_, e_);
// dissipation (Eq. 3.24, p.50)
volScalarField gammaCoeff =
12.0*(1.0 - e_*e_)*sqr(alpha)*rhoa_*gs0_*(1.0/da_)
*ThetaSqrt/sqrtPi;
12.0*(1.0 - sqr(e_))*sqr(alpha_)*rhoa_*gs0_*(1.0/da_)*ThetaSqrt/sqrtPi;
// Eq. 3.25, p. 50 Js = J1 - J2
volScalarField J1 = 3.0*betaPrim;
volScalarField J2 =
0.25*sqr(betaPrim)*da_*sqr(Ur)
/(alpha*rhoa_*sqrtPi*(ThetaSqrt + TsmallSqrt));
/(max(alpha_, 1e-6)*rhoa_*sqrtPi*(ThetaSqrt + TsmallSqrt));
// bulk viscosity p. 45 (Lun et al. 1984).
lambda_ = (4.0/3.0)*sqr(alpha_)*rhoa_*da_*gs0_*(1.0+e_)*ThetaSqrt/sqrtPi;
// stress tensor, Definitions, Table 3.1, p. 43
volTensorField tau = 2.0*mua_*D + (lambda_ - (2.0/3.0)*mua_)*tr(D)*I;
volSymmTensorField tau = 2.0*mua_*D + (lambda_ - (2.0/3.0)*mua_)*tr(D)*I;
if (!equilibrium_)
{
@ -268,8 +271,8 @@ void Foam::kineticTheoryModel::solve()
// wrong sign infront of laplacian
fvScalarMatrix ThetaEqn
(
fvm::ddt(1.5*alpha*rhoa_, Theta_)
+ fvm::div(phi, Theta_, scheme)
fvm::ddt(1.5*alpha_*rhoa_, Theta_)
+ fvm::div(phi, Theta_, "div(phi,Theta)")
==
fvm::SuSp(-((PsCoeff*I) && dU), Theta_)
+ (tau && dU)
@ -290,33 +293,31 @@ void Foam::kineticTheoryModel::solve()
volScalarField K3 = 0.5*da_*rhoa_*
(
(sqrtPi/(3.0*(3.0-e_)))
*(1.0 + 0.4*(1.0 + e_)*(3.0*e_ - 1.0)*alpha*gs0_)
+ 1.6*alpha*gs0_*(1.0 + e_)/sqrtPi
*(1.0 + 0.4*(1.0 + e_)*(3.0*e_ - 1.0)*alpha_*gs0_)
+1.6*alpha_*gs0_*(1.0 + e_)/sqrtPi
);
volScalarField K2 =
4.0*da_*rhoa_*(1.0 + e_)*alpha*gs0_/(3.0*sqrtPi) - 2.0*K3/3.0;
4.0*da_*rhoa_*(1.0 + e_)*alpha_*gs0_/(3.0*sqrtPi) - 2.0*K3/3.0;
volScalarField K4 = 12.0*(1.0 - e_*e_)*rhoa_*gs0_/(da_*sqrtPi);
volScalarField K4 = 12.0*(1.0 - sqr(e_))*rhoa_*gs0_/(da_*sqrtPi);
volScalarField trD = tr(D);
volTensorField D2 = D & D;
volScalarField tr2D = trD*trD;
volScalarField trD2 = tr(D2);
volScalarField tr2D = sqr(trD);
volScalarField trD2 = tr(D & D);
volScalarField t1 = K1*alpha + rhoa_;
volScalarField t1 = K1*alpha_ + rhoa_;
volScalarField l1 = -t1*trD;
volScalarField l2 = sqr(t1)*tr2D;
volScalarField l3 = 4.0*K4*alpha*(2.0*K3*trD2 + K2*tr2D);
volScalarField l3 = 4.0*K4*max(alpha_, 1e-6)*(2.0*K3*trD2 + K2*tr2D);
Theta_ = sqr((l1 + sqrt(l2 + l3))/(2.0*(alpha + 1.0e-4)*K4));
Theta_ = sqr((l1 + sqrt(l2 + l3))/(2.0*(alpha_ + 1.0e-4)*K4));
}
Theta_.max(1.0e-15);
Theta_.min(1.0e+3);
volScalarField pf =
frictionalStressModel_->frictionalPressure
volScalarField pf = frictionalStressModel_->frictionalPressure
(
alpha_,
alphaMinFriction_,
@ -344,13 +345,11 @@ void Foam::kineticTheoryModel::solve()
phi_
);
// add frictional stress for alpha > alphaMinFriction
mua_ = viscosityModel_->mua(alpha, Theta_, gs0_, rhoa_, da_, e_) + muf;
// add frictional stress
mua_ += muf;
mua_.min(1.0e+2);
mua_.max(0.0);
lambda_ = (4.0/3.0)*sqr(alpha_)*rhoa_*da_*gs0_*(1.0 + e_)*ThetaSqrt/sqrtPi;
Info<< "kinTheory: max(Theta) = " << max(Theta_).value() << endl;
volScalarField ktn = mua_/rhoa_;

2
applications/solvers/multiphase/twoPhaseEulerFoam/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.H
View File

@ -156,7 +156,7 @@ public:
// Member Functions
void solve();
void solve(const volTensorField& gradUat);
bool on() const
{

8
applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H
View File

@ -5,6 +5,9 @@
volScalarField rUaA = 1.0/UaEqn.A();
volScalarField rUbA = 1.0/UbEqn.A();
phia == (fvc::interpolate(Ua) & mesh.Sf());
phib == (fvc::interpolate(Ub) & mesh.Sf());
rUaAf = fvc::interpolate(rUaA);
surfaceScalarField rUbAf = fvc::interpolate(rUbA);
@ -47,11 +50,10 @@
surfaceScalarField Dp
(
"(rho*(1|A(U)))",
alphaf*rUaAf/rhoa + betaf*rUbAf/rhob
"(rho*(1|A(U)))", alphaf*rUaAf/rhoa + betaf*rUbAf/rhob
);
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
for(int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
fvScalarMatrix pEqn
(

5
applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C
View File

@ -92,11 +92,6 @@ int main(int argc, char *argv[])
#include "kEpsilon.H"
if (kineticTheory.on())
{
kineticTheory.solve();
nuEffa += kineticTheory.mua()/rhoa;
}
#include "write.H"
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"

4
applications/utilities/postProcessing/graphics/PV3Readers/Allwmake
View File

@ -7,9 +7,9 @@ then
case "$ParaView_VERSION" in
3* | git)
if [ ! -d "${PV_PLUGIN_PATH}" ]
if [ ! -n "${PV_PLUGIN_PATH}" ]
then
echo "$0 : PV_PLUGIN_PATH not a valid directory."
echo "$0 : PV_PLUGIN_PATH not a valid."
exit 1
fi