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/dm4/OpenFOAM/OpenFOAM-dev

Browse Source
This commit is contained in:
andy
2012-02-29 14:23:13 +00:00
parent 9ad38ed168 b5737931e5
commit 3bbe04eae8
20 changed files with 446 additions and 87 deletions
Download Patch File Download Diff File Expand all files Collapse all files

3
applications/solvers/multiphase/bubbleFoam/bubbleFoam.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -90,6 +90,7 @@ int main(int argc, char *argv[])
if (pimple.turbCorr())
{
#include "kEpsilon.H"
nuEffa = sqr(Ct)*nutb + nua;
}
}

15
applications/solvers/multiphase/bubbleFoam/kEpsilon.H
View File

@ -14,16 +14,17 @@ if (turbulence)
// Dissipation equation
fvScalarMatrix epsEqn
(
fvm::ddt(beta, epsilon)
fvm::ddt(epsilon)
+ fvm::div(phib, epsilon)
- fvm::Sp(fvc::div(phib), epsilon)
- fvm::laplacian
(
alphaEps*nuEffb, epsilon,
"laplacian(DepsilonEff,epsilon)"
)
==
C1*beta*G*epsilon/k
- fvm::Sp(C2*beta*epsilon/k, epsilon)
C1*G*epsilon/k
- fvm::Sp(C2*epsilon/k, epsilon)
);
#include "wallDissipation.H"
@ -37,16 +38,17 @@ if (turbulence)
// Turbulent kinetic energy equation
fvScalarMatrix kEqn
(
fvm::ddt(beta, k)
fvm::ddt(k)
+ fvm::div(phib, k)
- fvm::Sp(fvc::div(phib), k)
- fvm::laplacian
(
alphak*nuEffb, k,
"laplacian(DkEff,k)"
)
==
beta*G
- fvm::Sp(beta*epsilon/k, k)
G
- fvm::Sp(epsilon/k, k)
);
kEqn.relax();
kEqn.solve();
@ -59,5 +61,4 @@ if (turbulence)
#include "wallViscosity.H"
}
nuEffa = sqr(Ct)*nutb + nua;
nuEffb = nutb + nub;

1
applications/solvers/multiphase/twoPhaseEulerFoam/Allwclean
View File

@ -6,5 +6,6 @@ wclean libso phaseModel
wclean libso interfacialModels
wclean libso kineticTheoryModels
wclean
wclean MRFtwoPhaseEulerFoam
# ----------------------------------------------------------------- end-of-file

1
applications/solvers/multiphase/twoPhaseEulerFoam/Allwmake
View File

@ -6,5 +6,6 @@ wmake libso phaseModel
wmake libso interfacialModels
wmake libso kineticTheoryModels
wmake
wmake MRFtwoPhaseEulerFoam
# ----------------------------------------------------------------- end-of-file

119
applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/MRFTwoPhaseEulerFoam.C Normal file
View File

@ -0,0 +1,119 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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
twoPhaseEulerFoam
Description
Solver for a system of 2 incompressible fluid phases with one phase
dispersed, e.g. gas bubbles in a liquid or solid particles in a gas.
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "nearWallDist.H"
#include "wallFvPatch.H"
#include "Switch.H"
#include "IFstream.H"
#include "OFstream.H"
#include "dragModel.H"
#include "phaseModel.H"
#include "kineticTheoryModel.H"
#include "pimpleControl.H"
#include "MRFZones.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
#include "setRootCase.H"
#include "createTime.H"
#include "createMesh.H"
#include "readGravitationalAcceleration.H"
#include "createFields.H"
#include "readPPProperties.H"
#include "initContinuityErrs.H"
#include "createMRFZones.H"
#include "readTimeControls.H"
#include "CourantNo.H"
#include "setInitialDeltaT.H"
pimpleControl pimple(mesh);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl;
while (runTime.run())
{
#include "readTwoPhaseEulerFoamControls.H"
#include "CourantNos.H"
#include "setDeltaT.H"
runTime++;
Info<< "Time = " << runTime.timeName() << nl << endl;
// --- Pressure-velocity PIMPLE corrector loop
while (pimple.loop())
{
#include "alphaEqn.H"
#include "liftDragCoeffs.H"
#include "UEqns.H"
// --- Pressure corrector loop
while (pimple.correct())
{
#include "pEqn.H"
if (correctAlpha && !pimple.finalIter())
{
#include "alphaEqn.H"
}
}
#include "DDtU.H"
if (pimple.turbCorr())
{
#include "kEpsilon.H"
}
}
#include "write.H"
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
}
Info<< "End\n" << endl;
return 0;
}
// ************************************************************************* //

3
applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/Make/files Normal file
View File

@ -0,0 +1,3 @@
MRFTwoPhaseEulerFoam.C
EXE = $(FOAM_APPBIN)/MRFTwoPhaseEulerFoam

17
applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/Make/options Normal file
View File

@ -0,0 +1,17 @@
EXE_INC = \
-I.. \
-I../../bubbleFoam \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/transportModels/incompressible/lnInclude \
-IturbulenceModel \
-I../kineticTheoryModels/lnInclude \
-I../interfacialModels/lnInclude \
-I../phaseModel/lnInclude \
EXE_LIBS = \
-lEulerianInterfacialModels \
-lfiniteVolume \
-lmeshTools \
-lincompressibleTransportModels \
-lphaseModel \
-lkineticTheoryModel

99
applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/UEqns.H Normal file
View File

@ -0,0 +1,99 @@
fvVectorMatrix UaEqn(Ua, Ua.dimensions()*dimVol/dimTime);
fvVectorMatrix UbEqn(Ub, Ub.dimensions()*dimVol/dimTime);
{
{
volTensorField gradUaT(T(fvc::grad(Ua)));
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)*nuEffa)*tr(gradUaT) - nuEffa*gradUaT
);
if (kineticTheory.on())
{
Rca -= ((kineticTheory.lambda()/rhoa)*tr(gradUaT))*tensor(I);
}
surfaceScalarField phiRa
(
-fvc::interpolate(nuEffa)*mesh.magSf()*fvc::snGrad(alpha)
/fvc::interpolate(alpha + scalar(0.001))
);
UaEqn =
(
(scalar(1) + Cvm*rhob*beta/rhoa)*
(
fvm::ddt(Ua)
+ fvm::div(phia, Ua, "div(phia,Ua)")
- fvm::Sp(fvc::div(phia), Ua)
)
- fvm::laplacian(nuEffa, Ua)
+ fvc::div(Rca)
+ fvm::div(phiRa, Ua, "div(phia,Ua)")
- fvm::Sp(fvc::div(phiRa), Ua)
+ (fvc::grad(alpha)/(fvc::average(alpha) + scalar(0.001)) & Rca)
==
// g // Buoyancy term transfered to p-equation
- fvm::Sp(beta/rhoa*K, Ua)
//+ beta/rhoa*K*Ub // Explicit drag transfered to p-equation
- beta/rhoa*(liftCoeff - Cvm*rhob*DDtUb)
);
mrfZones.addCoriolis(UaEqn);
UaEqn.relax();
}
{
volTensorField gradUbT(T(fvc::grad(Ub)));
volTensorField Rcb
(
"Rcb",
(((2.0/3.0)*I)*nuEffb)*tr(gradUbT) - nuEffb*gradUbT
);
surfaceScalarField phiRb
(
-fvc::interpolate(nuEffb)*mesh.magSf()*fvc::snGrad(beta)
/fvc::interpolate(beta + scalar(0.001))
);
UbEqn =
(
(scalar(1) + Cvm*rhob*alpha/rhob)*
(
fvm::ddt(Ub)
+ fvm::div(phib, Ub, "div(phib,Ub)")
- fvm::Sp(fvc::div(phib), Ub)
)
- fvm::laplacian(nuEffb, Ub)
+ fvc::div(Rcb)
+ fvm::div(phiRb, Ub, "div(phib,Ub)")
- fvm::Sp(fvc::div(phiRb), Ub)
+ (fvc::grad(beta)/(fvc::average(beta) + scalar(0.001)) & Rcb)
==
// g // Buoyancy term transfered to p-equation
- fvm::Sp(alpha/rhob*K, Ub)
//+ alpha/rhob*K*Ua // Explicit drag transfered to p-equation
+ alpha/rhob*(liftCoeff + Cvm*rhob*DDtUa)
);
mrfZones.addCoriolis(UbEqn);
UbEqn.relax();
}
}

3
applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/createMRFZones.H Normal file
View File

@ -0,0 +1,3 @@
MRFZones mrfZones(mesh);
mrfZones.correctBoundaryVelocity(Ua);
mrfZones.correctBoundaryVelocity(Ub);

121
applications/solvers/multiphase/twoPhaseEulerFoam/MRFtwoPhaseEulerFoam/pEqn.H Normal file
View File

@ -0,0 +1,121 @@
{
surfaceScalarField alphaf(fvc::interpolate(alpha));
surfaceScalarField betaf(scalar(1) - alphaf);
volScalarField rUaA(1.0/UaEqn.A());
volScalarField rUbA(1.0/UbEqn.A());
rUaAf = fvc::interpolate(rUaA);
surfaceScalarField rUbAf(fvc::interpolate(rUbA));
volVectorField HbyAa("HbyAa", Ua);
HbyAa = rUaA*UaEqn.H();
volVectorField HbyAb("HbyAb", Ub);
HbyAb = rUbA*UbEqn.H();
mrfZones.absoluteFlux(phia.oldTime());
mrfZones.absoluteFlux(phia);
mrfZones.absoluteFlux(phib.oldTime());
mrfZones.absoluteFlux(phib);
surfaceScalarField phiDraga
(
fvc::interpolate(beta/rhoa*K*rUaA)*phib + rUaAf*(g & mesh.Sf())
);
if (g0.value() > 0.0)
{
phiDraga -= ppMagf*fvc::snGrad(alpha)*mesh.magSf();
}
if (kineticTheory.on())
{
phiDraga -= rUaAf*fvc::snGrad(kineticTheory.pa()/rhoa)*mesh.magSf();
}
surfaceScalarField phiDragb
(
fvc::interpolate(alpha/rhob*K*rUbA)*phia + rUbAf*(g & mesh.Sf())
);
// Fix for gravity on outlet boundary.
forAll(p.boundaryField(), patchi)
{
if (isA<zeroGradientFvPatchScalarField>(p.boundaryField()[patchi]))
{
phiDraga.boundaryField()[patchi] = 0.0;
phiDragb.boundaryField()[patchi] = 0.0;
}
}
surfaceScalarField phiHbyAa
(
"phiHbyAa",
(fvc::interpolate(HbyAa) & mesh.Sf())
+ fvc::ddtPhiCorr(rUaA, Ua, phia)
+ phiDraga
);
mrfZones.relativeFlux(phiHbyAa);
surfaceScalarField phiHbyAb
(
"phiHbyAb",
(fvc::interpolate(HbyAb) & mesh.Sf())
+ fvc::ddtPhiCorr(rUbA, Ub, phib)
+ phiDragb
);
mrfZones.relativeFlux(phiHbyAb);
surfaceScalarField phiHbyA("phiHbyA", alphaf*phiHbyAa + betaf*phiHbyAb);
surfaceScalarField Dp
(
"Dp",
alphaf*rUaAf/rhoa + betaf*rUbAf/rhob
);
while (pimple.correctNonOrthogonal())
{
fvScalarMatrix pEqn
(
fvm::laplacian(Dp, p) == fvc::div(phiHbyA)
);
pEqn.setReference(pRefCell, pRefValue);
pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
if (pimple.finalNonOrthogonalIter())
{
surfaceScalarField SfGradp(pEqn.flux()/Dp);
phia.boundaryField() ==
(fvc::interpolate(Ua) & mesh.Sf())().boundaryField();
mrfZones.relativeFlux(phia);
phia = phiHbyAa - rUaAf*SfGradp/rhoa;
phib.boundaryField() ==
(fvc::interpolate(Ub) & mesh.Sf())().boundaryField();
mrfZones.relativeFlux(phib);
phib = phiHbyAb - rUbAf*SfGradp/rhob;
phi = alphaf*phia + betaf*phib;
p.relax();
SfGradp = pEqn.flux()/Dp;
Ua = HbyAa + fvc::reconstruct(phiDraga - rUaAf*SfGradp/rhoa);
Ua.correctBoundaryConditions();
Ub = HbyAb + fvc::reconstruct(phiDragb - rUbAf*SfGradp/rhob);
Ub.correctBoundaryConditions();
U = alpha*Ua + beta*Ub;
}
}
}
#include "continuityErrs.H"

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

@ -18,7 +18,7 @@ fvVectorMatrix UbEqn(Ub, Ub.dimensions()*dimVol/dimTime);
volTensorField Rca
(
"Rca",
((2.0/3.0)*I)*(sqr(Ct)*k + nuEffa*tr(gradUaT)) - nuEffa*gradUaT
(((2.0/3.0)*I)*nuEffa)*tr(gradUaT) - nuEffa*gradUaT
);
if (kineticTheory.on())
@ -62,7 +62,7 @@ fvVectorMatrix UbEqn(Ub, Ub.dimensions()*dimVol/dimTime);
volTensorField Rcb
(
"Rcb",
((2.0/3.0)*I)*(k + nuEffb*tr(gradUbT)) - nuEffb*gradUbT
(((2.0/3.0)*I)*nuEffb)*tr(gradUbT) - nuEffb*gradUbT
);
surfaceScalarField phiRb

15
applications/solvers/multiphase/twoPhaseEulerFoam/createFields.H
View File

@ -112,7 +112,9 @@
(
"phi",
runTime.timeName(),
mesh
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
fvc::interpolate(alpha)*phia + fvc::interpolate(beta)*phib
);
@ -194,6 +196,17 @@
}
}
dimensionedScalar residualSlip
(
dimensionedScalar::lookupOrDefault
(
"residualSlip",
interfacialProperties,
0,
dimVelocity
)
);
Info << "dragPhase is " << dragPhase << endl;
kineticTheoryModel kineticTheory
(

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

@ -1,62 +0,0 @@
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/liftDragCoeffs.H
View File

@ -1,5 +1,5 @@
volVectorField Ur(Ua - Ub);
volScalarField magUr(mag(Ur));
volScalarField magUr(mag(Ur) + residualSlip);
volScalarField Ka(draga->K(magUr));
volScalarField K(Ka);

18
applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/initialPointsMethod/autoDensity/autoDensity.C
View File

@ -849,7 +849,10 @@ autoDensity::autoDensity
:
initialPointsMethod(typeName, initialPointsDict, cvMesh),
globalTrialPoints_(0),
minCellSizeLimit_(readScalar(detailsDict().lookup("minCellSizeLimit"))),
minCellSizeLimit_
(
detailsDict().lookupOrDefault<scalar>("minCellSizeLimit", 0.0)
),
minLevels_(readLabel(detailsDict().lookup("minLevels"))),
maxSizeRatio_(readScalar(detailsDict().lookup("maxSizeRatio"))),
volRes_(readLabel(detailsDict().lookup("sampleResolution"))),
@ -899,18 +902,7 @@ List<Vb::Point> autoDensity::initialPoints() const
);
}
// Initialise size of points list.
const scalar volumeBoundBox = Foam::pow3(hierBB.typDim());
const scalar volumeSmallestCell = Foam::pow3(minCellSizeLimit_);
const int initialPointEstimate
= min
(
static_cast<int>(volumeBoundBox/(volumeSmallestCell + SMALL)/10),
1e6
);
DynamicList<Vb::Point> initialPoints(initialPointEstimate);
DynamicList<Vb::Point> initialPoints;
Info<< nl << " " << typeName << endl;

43
src/OpenFOAM/meshes/boundBox/boundBox.C
View File

@ -163,6 +163,49 @@ Foam::tmp<Foam::pointField> Foam::boundBox::points() const
}
Foam::faceList Foam::boundBox::faces()
{
faceList faces(6);
forAll(faces, fI)
{
faces[fI].setSize(4);
}
faces[0][0] = 0;
faces[0][1] = 1;
faces[0][2] = 2;
faces[0][3] = 3;
faces[1][0] = 2;
faces[1][1] = 6;
faces[1][2] = 7;
faces[1][3] = 3;
faces[2][0] = 0;
faces[2][1] = 4;
faces[2][2] = 5;
faces[2][3] = 1;
faces[3][0] = 4;
faces[3][1] = 7;
faces[3][2] = 6;
faces[3][3] = 5;
faces[4][0] = 3;
faces[4][1] = 7;
faces[4][2] = 4;
faces[4][3] = 0;
faces[5][0] = 1;
faces[5][1] = 5;
faces[5][2] = 6;
faces[5][3] = 2;
return faces;
}
void Foam::boundBox::inflate(const scalar s)
{
vector ext = vector::one*s*mag();

4
src/OpenFOAM/meshes/boundBox/boundBox.H
View File

@ -33,6 +33,7 @@ Description
#define boundBox_H
#include "pointField.H"
#include "faceList.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -160,6 +161,9 @@ public:
//- Return corner points in an order corresponding to a 'hex' cell
tmp<pointField> points() const;
//- Return faces with correct point order
static faceList faces();
// Manipulate

1
tutorials/multiphase/twoPhaseEulerFoam/bed/constant/interfacialProperties
View File

@ -21,5 +21,6 @@ dragModelb GidaspowSchillerNaumann;
dragPhase a;
residualSlip 0;
// ************************************************************************* //

1
tutorials/multiphase/twoPhaseEulerFoam/bed2/constant/interfacialProperties
View File

@ -21,5 +21,6 @@ dragModelb GidaspowErgunWenYu;
dragPhase a;
residualSlip 0;
// ************************************************************************* //

1
tutorials/multiphase/twoPhaseEulerFoam/bubbleColumn/constant/interfacialProperties
View File

@ -21,5 +21,6 @@ dragModelb SchillerNaumann;
dragPhase blended;
residualSlip 0;
// ************************************************************************* //