Merge branch 'cvm' of /home/noisy3/OpenFOAM/OpenFOAM-dev into cvm

This commit is contained in:
mattijs
2011-07-06 15:28:25 +01:00
20 changed files with 231 additions and 116 deletions

View File

@ -1,4 +1,4 @@
volScalarField rAUrel = 1.0/UrelEqn().A();
volScalarField rAUrel(1.0/UrelEqn().A());
Urel = rAUrel*UrelEqn().H();
if (pimple.nCorr() <= 1)

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -73,34 +73,23 @@ Foam::tmp<Foam::volScalarField> Foam::GidaspowErgunWenYu::K
volScalarField bp(pow(beta, -2.65));
volScalarField Re(max(Ur*phasea_.d()/phaseb_.nu(), scalar(1.0e-3)));
volScalarField Cds(24.0*(1.0 + 0.15*pow(Re, 0.687))/Re);
forAll(Re, celli)
{
if (Re[celli] > 1000.0)
{
Cds[celli] = 0.44;
}
}
volScalarField Cds
(
neg(Re - 1000)*(24.0*(1.0 + 0.15*pow(Re, 0.687))/Re)
+ pos(Re - 1000)*0.44
);
// Wen and Yu (1966)
tmp<volScalarField> tKWenYu = 0.75*Cds*phaseb_.rho()*Ur*bp/phasea_.d();
volScalarField& KWenYu = tKWenYu();
// Ergun
forAll (beta, cellj)
{
if (beta[cellj] <= 0.8)
{
KWenYu[cellj] =
150.0*alpha_[cellj]*phaseb_.nu().value()*phaseb_.rho().value()
/sqr(beta[cellj]*phasea_.d().value())
+ 1.75*phaseb_.rho().value()*Ur[cellj]
/(beta[cellj]*phasea_.d().value());
}
}
return tKWenYu;
return
(
pos(beta - 0.8)
*(0.75*Cds*phaseb_.rho()*Ur*bp/phasea_.d())
+ neg(beta - 0.8)
*(
150.0*alpha_*phaseb_.nu()*phaseb_.rho()/(sqr(beta*phasea_.d()))
+ 1.75*phaseb_.rho()*Ur/(beta*phasea_.d())
)
);
}

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -72,15 +72,12 @@ Foam::tmp<Foam::volScalarField> Foam::GidaspowSchillerNaumann::K
volScalarField bp(pow(beta, -2.65));
volScalarField Re(max(beta*Ur*phasea_.d()/phaseb_.nu(), scalar(1.0e-3)));
volScalarField Cds(24.0*(scalar(1) + 0.15*pow(Re, 0.687))/Re);
forAll(Re, celli)
{
if (Re[celli] > 1000.0)
{
Cds[celli] = 0.44;
}
}
volScalarField Cds
(
neg(Re - 1000)*(24.0*(1.0 + 0.15*pow(Re, 0.687))/Re)
+ pos(Re - 1000)*0.44
);
return 0.75*Cds*phaseb_.rho()*Ur*bp/phasea_.d();
}

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -69,15 +69,12 @@ Foam::tmp<Foam::volScalarField> Foam::SchillerNaumann::K
) const
{
volScalarField Re(max(Ur*phasea_.d()/phaseb_.nu(), scalar(1.0e-3)));
volScalarField Cds(24.0*(scalar(1) + 0.15*pow(Re, 0.687))/Re);
forAll(Re, celli)
{
if (Re[celli] > 1000.0)
{
Cds[celli] = 0.44;
}
}
volScalarField Cds
(
neg(Re - 1000)*(24.0*(1.0 + 0.15*pow(Re, 0.687))/Re)
+ pos(Re - 1000)*0.44
);
return 0.75*Cds*phaseb_.rho()*Ur/phasea_.d();
}

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -70,15 +70,11 @@ Foam::tmp<Foam::volScalarField> Foam::SyamlalOBrien::K
{
volScalarField beta(max(scalar(1) - alpha_, scalar(1.0e-6)));
volScalarField A(pow(beta, 4.14));
volScalarField B(0.8*pow(beta, 1.28));
forAll (beta, celli)
{
if (beta[celli] > 0.85)
{
B[celli] = pow(beta[celli], 2.65);
}
}
volScalarField B
(
neg(beta - 0.85)*(0.8*pow(beta, 1.28))
+ pos(beta - 0.85)*(pow(beta, 2.65))
);
volScalarField Re(max(Ur*phasea_.d()/phaseb_.nu(), scalar(1.0e-3)));

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -72,15 +72,11 @@ Foam::tmp<Foam::volScalarField> Foam::WenYu::K
volScalarField bp(pow(beta, -2.65));
volScalarField Re(max(Ur*phasea_.d()/phaseb_.nu(), scalar(1.0e-3)));
volScalarField Cds(24.0*(scalar(1) + 0.15*pow(Re, 0.687))/Re);
forAll(Re, celli)
{
if (Re[celli] > 1000.0)
{
Cds[celli] = 0.44;
}
}
volScalarField Cds
(
neg(Re - 1000)*(24.0*(1.0 + 0.15*pow(Re, 0.687))/Re)
+ pos(Re - 1000)*0.44
);
return 0.75*Cds*phaseb_.rho()*Ur*bp/phasea_.d();
}

View File

@ -150,7 +150,7 @@ void drawHitProblem
}
scalarField curvature(const triSurface& surf)
scalarField calcCurvature(const triSurface& surf)
{
scalarField k(surf.points().size(), 0);
@ -400,6 +400,16 @@ int main(int argc, char *argv[])
"writeObj",
"write extendedFeatureEdgeMesh obj files"
);
argList::addBoolOption
(
"writeVTK",
"write extendedFeatureEdgeMesh vtk files"
);
argList::addBoolOption
(
"calcCurvature",
"calculate curvature and closeness fields"
);
argList::addOption
(
"closeness",
@ -427,7 +437,13 @@ int main(int argc, char *argv[])
# include "setRootCase.H"
# include "createTime.H"
if (env("FOAM_SIGFPE"))
bool writeVTK = args.optionFound("writeVTK");
bool writeObj = args.optionFound("writeObj");
bool curvature = args.optionFound("curvature");
if (curvature && env("FOAM_SIGFPE"))
{
WarningIn(args.executable())
<< "Detected floating point exception trapping (FOAM_SIGFPE)."
@ -441,10 +457,6 @@ int main(int argc, char *argv[])
Info<< "Feature line extraction is only valid on closed manifold surfaces."
<< endl;
bool writeVTK = args.optionFound("writeVTK");
bool writeObj = args.optionFound("writeObj");
const fileName surfFileName = args[1];
const fileName outFileName = args[2];
@ -674,6 +686,13 @@ int main(int argc, char *argv[])
surf
);
if (!curvature)
{
Info<< "End\n" << endl;
return 0;
}
// Find close features
// // Dummy trim operation to mark features
@ -741,7 +760,8 @@ int main(int argc, char *argv[])
// )
// );
// Examine curvature, feature proximity and internal and external closeness.
Info<< "Examine curvature, feature proximity and internal and "
<< "external closeness." << endl;
// Internal and external closeness
@ -935,7 +955,7 @@ int main(int argc, char *argv[])
externalClosenessField.write();
scalarField k = curvature(surf);
scalarField k = calcCurvature(surf);
// Modify the curvature values on feature edges and points to be zero.

View File

@ -59,6 +59,10 @@ Description
#ifdef USE_RANDOM
# include <climits>
# if INT_MAX != 2147483647
# error "INT_MAX != 2147483647"
# error "The random number generator may not work!"
# endif
#endif
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -1250,7 +1254,7 @@ Foam::label Foam::osRandomInteger()
Foam::scalar Foam::osRandomDouble()
{
#ifdef USE_RANDOM
return (scalar)random();
return (scalar)random()/INT_MAX;
#else
return drand48();
#endif

View File

@ -1,6 +1,6 @@
global/global.Cver
global/constants/constants.C
global/constants/dimensionedConstants.C
/* global/constants/constants.C in global.Cver */
/* global/constants/dimensionedConstants.C in global.Cver */
global/argList/argList.C
global/clock/clock.C
@ -136,7 +136,7 @@ $(StringStreams)/StringStreamsPrint.C
Pstreams = $(Streams)/Pstreams
$(Pstreams)/UIPstream.C
$(Pstreams)/IPstream.C
$(Pstreams)/UPstream.C
/* $(Pstreams)/UPstream.C in global.Cver */
$(Pstreams)/UPstreamCommsStruct.C
$(Pstreams)/Pstream.C
$(Pstreams)/UOPstream.C
@ -181,7 +181,7 @@ $(IOobject)/IOobjectReadHeader.C
$(IOobject)/IOobjectWriteHeader.C
regIOobject = db/regIOobject
$(regIOobject)/regIOobject.C
/* $(regIOobject)/regIOobject.C in global.Cver */
$(regIOobject)/regIOobjectRead.C
$(regIOobject)/regIOobjectWrite.C

View File

@ -197,8 +197,22 @@ Foam::functionEntries::codeStream::getFunction
}
}
// all processes must wait for compile to finish
// all processes must wait for compile to finish - except if this
// file is only read on the master
bool masterOnly =
(
regIOobject::fileModificationChecking
== regIOobject::timeStampMaster
)
|| (
regIOobject::fileModificationChecking
== regIOobject::inotifyMaster
);
if (!masterOnly)
{
reduce(create, orOp<bool>());
}
if (isA<IOdictionary>(topDict(parentDict)))
{

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -65,6 +65,22 @@ bool Foam::JobInfo::constructed(false);
#include "debug.C"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// Read file modification checking switches
#include "regIOobject.C"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// Read parallel communication switches
#include "UPstream.C"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// Read constants
#include "constants.C"
#include "dimensionedConstants.C"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// Read and set cell models

View File

@ -83,10 +83,8 @@ class fanPressureFvPatchScalarField
:
public totalPressureFvPatchScalarField
{
// Private data
//- Tabulated fan curve
interpolationTable<scalar> fanCurve_;
public:
//- Fan flow direction
enum fanFlowDirection
@ -97,6 +95,13 @@ class fanPressureFvPatchScalarField
static const NamedEnum<fanFlowDirection, 2> fanFlowDirectionNames_;
private:
// Private data
//- Tabulated fan curve
interpolationTable<scalar> fanCurve_;
//- Direction of flow through the fan relative to patch
fanFlowDirection direction_;

View File

@ -1318,7 +1318,7 @@ bool Foam::conformalVoronoiMesh::distributeBackground()
}
Info<< " Total number of vertices before redistribution "
<< globalNRealVertices
<< returnReduce(label(number_of_vertices()), sumOp<label>())
<< endl;
// Pout<< " Real vertices before distribution "
@ -1446,9 +1446,8 @@ bool Foam::conformalVoronoiMesh::distributeBackground()
// processor each point should be on, so give distribute = false
insertPoints(pointsToInsert, indices, types, false);
// Adjust by 8 because of insertion of bounding points by reset
Info<< " Total number of vertices after redistribution "
<< returnReduce(label(number_of_vertices() - 8), sumOp<label>())
<< returnReduce(label(number_of_vertices()), sumOp<label>())
<< endl;
// Pout<< " Real vertices after distribution "

View File

@ -203,7 +203,7 @@ void sixDoFRigidBodyDisplacementPointPatchVectorField::updateCoeffs()
// calculate the forces on the motion object from this data, then
// update the positions
motion_.updatePosition(t.deltaTValue());
motion_.updatePosition(t.deltaTValue(), t.deltaT0Value());
dictionary forcesDict;

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2009-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2009-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -174,6 +174,8 @@ Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion()
initialQ_(I),
momentOfInertia_(diagTensor::one*VSMALL),
mass_(VSMALL),
cDamp_(0.0),
aLim_(VGREAT),
report_(false)
{}
@ -190,6 +192,8 @@ Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion
const point& initialCentreOfMass,
const tensor& initialQ,
const diagTensor& momentOfInertia,
scalar cDamp,
scalar aLim,
bool report
)
:
@ -211,6 +215,8 @@ Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion
initialQ_(initialQ),
momentOfInertia_(momentOfInertia),
mass_(mass),
cDamp_(cDamp),
aLim_(aLim),
report_(report)
{}
@ -233,6 +239,8 @@ Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion(const dictionary& dict)
),
momentOfInertia_(dict.lookup("momentOfInertia")),
mass_(readScalar(dict.lookup("mass"))),
cDamp_(dict.lookupOrDefault<scalar>("accelerationDampingCoeff", 0.0)),
aLim_(dict.lookupOrDefault<scalar>("accelerationLimit", VGREAT)),
report_(dict.lookupOrDefault<Switch>("report", false))
{
addRestraints(dict);
@ -246,17 +254,19 @@ Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion
const sixDoFRigidBodyMotion& sDoFRBM
)
:
motionState_(sDoFRBM.motionState()),
restraints_(sDoFRBM.restraints()),
restraintNames_(sDoFRBM.restraintNames()),
constraints_(sDoFRBM.constraints()),
constraintNames_(sDoFRBM.constraintNames()),
maxConstraintIterations_(sDoFRBM.maxConstraintIterations()),
initialCentreOfMass_(sDoFRBM.initialCentreOfMass()),
initialQ_(sDoFRBM.initialQ()),
momentOfInertia_(sDoFRBM.momentOfInertia()),
mass_(sDoFRBM.mass()),
report_(sDoFRBM.report())
motionState_(sDoFRBM.motionState_),
restraints_(sDoFRBM.restraints_),
restraintNames_(sDoFRBM.restraintNames_),
constraints_(sDoFRBM.constraints_),
constraintNames_(sDoFRBM.constraintNames_),
maxConstraintIterations_(sDoFRBM.maxConstraintIterations_),
initialCentreOfMass_(sDoFRBM.initialCentreOfMass_),
initialQ_(sDoFRBM.initialQ_),
momentOfInertia_(sDoFRBM.momentOfInertia_),
mass_(sDoFRBM.mass_),
cDamp_(sDoFRBM.cDamp_),
aLim_(sDoFRBM.aLim_),
report_(sDoFRBM.report_)
{}
@ -358,7 +368,8 @@ void Foam::sixDoFRigidBodyMotion::addConstraints
void Foam::sixDoFRigidBodyMotion::updatePosition
(
scalar deltaT
scalar deltaT,
scalar deltaT0
)
{
// First leapfrog velocity adjust and motion part, required before
@ -366,9 +377,32 @@ void Foam::sixDoFRigidBodyMotion::updatePosition
if (Pstream::master())
{
v() += 0.5*deltaT*a();
vector aClip = a();
scalar aMag = mag(aClip);
pi() += 0.5*deltaT*tau();
if (aMag > SMALL)
{
aClip /= aMag;
}
if (aMag > aLim_)
{
WarningIn
(
"void Foam::sixDoFRigidBodyMotion::updatePosition"
"("
"scalar deltaT, "
"scalar deltaT0"
")"
)
<< "Limited acceleration " << a()
<< " to " << aClip*aLim_
<< endl;
}
v() += 0.5*(1 - cDamp_)*deltaT0*aClip*min(aMag, aLim_);
pi() += 0.5*(1 - cDamp_)*deltaT0*tau();
// Leapfrog move part
centreOfMass() += deltaT*v();
@ -401,9 +435,33 @@ void Foam::sixDoFRigidBodyMotion::updateForce
applyConstraints(deltaT);
v() += 0.5*deltaT*a();
vector aClip = a();
scalar aMag = mag(aClip);
pi() += 0.5*deltaT*tau();
if (aMag > SMALL)
{
aClip /= aMag;
}
if (aMag > aLim_)
{
WarningIn
(
"void Foam::sixDoFRigidBodyMotion::updateForce"
"("
"const vector& fGlobal, "
"const vector& tauGlobal, "
"scalar deltaT"
")"
)
<< "Limited acceleration " << a()
<< " to " << aClip*aLim_
<< endl;
}
v() += 0.5*(1 - cDamp_)*deltaT*aClip*min(aMag, aLim_);
pi() += 0.5*(1 - cDamp_)*deltaT*tau();
if(report_)
{

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2009-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2009-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -116,6 +116,15 @@ class sixDoFRigidBodyMotion
//- Mass of the body
scalar mass_;
//- Acceleration damping coefficient. Modify applied acceleration:
// v1 = v0 + a*dt - cDamp*a*dt
// = v0 + dt*f*(1 - cDamp)/m
// Increases effective mass by 1/(1 - cDamp).
scalar cDamp_;
//- Acceleration magnitude limit - clips large accelerations
scalar aLim_;
//- Switch to turn reporting of motion data on and off
Switch report_;
@ -235,6 +244,8 @@ public:
const point& initialCentreOfMass,
const tensor& initialQ,
const diagTensor& momentOfInertia,
scalar cDamp = 0.0,
scalar aLim = VGREAT,
bool report = false
);
@ -260,10 +271,12 @@ public:
void addConstraints(const dictionary& dict);
//- First leapfrog velocity adjust and motion part, required
// before force calculation
// before force calculation. Takes old timestep for variable
// timestep cases.
void updatePosition
(
scalar deltaT
scalar deltaT,
scalar deltaT0
);
//- Second leapfrog velocity adjust part, required after motion and

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2009-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2009-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -40,6 +40,10 @@ void Foam::sixDoFRigidBodyMotion::write(Ostream& os) const
<< momentOfInertia_ << token::END_STATEMENT << nl;
os.writeKeyword("mass")
<< mass_ << token::END_STATEMENT << nl;
os.writeKeyword("accelerationDampingCoeff")
<< cDamp_ << token::END_STATEMENT << nl;
os.writeKeyword("accelerationLimit")
<< aLim_ << token::END_STATEMENT << nl;
os.writeKeyword("report")
<< report_ << token::END_STATEMENT << nl;

View File

@ -83,15 +83,20 @@ Foam::sixDoFRigidBodyMotionRestraints::linearAxialAngularSpring::restrain
) const
{
vector refDir = rotationTensor(vector(1, 0 ,0), axis_) & vector(0, 1, 0);
vector oldDir = refQ_ & refDir;
vector newDir = motion.orientation() & refDir;
if (mag(oldDir & axis_) > 0.95 || mag(newDir & axis_) > 0.95)
{
// Directions getting close to the axis, change reference
refDir = rotationTensor(vector(1, 0 ,0), axis_) & vector(0, 0, 1);
oldDir = refQ_ & refDir;
newDir = motion.orientation() & refDir;
vector oldDir = refQ_ & refDir;
vector newDir = motion.orientation() & refDir;
}
// Removing any axis component from oldDir and newDir and normalising

View File

@ -87,6 +87,7 @@ Foam::sixDoFRigidBodyMotionRestraints::tabulatedAxialAngularSpring::restrain
vector refDir = rotationTensor(vector(1, 0 ,0), axis_) & vector(0, 1, 0);
vector oldDir = refQ_ & refDir;
vector newDir = motion.orientation() & refDir;
if (mag(oldDir & axis_) > 0.95 || mag(newDir & axis_) > 0.95)
@ -95,8 +96,9 @@ Foam::sixDoFRigidBodyMotionRestraints::tabulatedAxialAngularSpring::restrain
refDir = rotationTensor(vector(1, 0 ,0), axis_) & vector(0, 0, 1);
oldDir = refQ_ & refDir;
newDir = motion.orientation() & refDir;
vector oldDir = refQ_ & refDir;
vector newDir = motion.orientation() & refDir;
}
// Removing any axis component from oldDir and newDir and normalising

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd.
\\ / A nd | Copyright (C) 2009-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -147,7 +147,7 @@ void uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField::updateCoeffs()
const polyMesh& mesh = this->dimensionedInternalField().mesh()();
const Time& t = mesh.time();
motion_.updatePosition(t.deltaTValue());
motion_.updatePosition(t.deltaTValue(), t.deltaT0Value());
vector gravity = vector::zero;