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:
mattijs
2013-10-31 12:57:33 +00:00
parent ed99c9d05d 7ecf8fa307
commit ced5391952
28 changed files with 277 additions and 199 deletions
Download Patch File Download Diff File Expand all files Collapse all files

3
applications/solvers/multiphase/interFoam/alphaEqn.H
View File

@ -6,7 +6,8 @@
phic = min(interface.cAlpha()*phic, max(phic));
surfaceScalarField phir(phic*interface.nHatf());
if (pimple.firstIter() && MULESCorr)
//***HGW if (pimple.firstIter() && MULESCorr)
if (MULESCorr)
{
fvScalarMatrix alpha1Eqn
(

9
applications/test/ODE/Test-ODE.C
View File

@ -27,7 +27,7 @@ Description
#include "argList.H"
#include "IOmanip.H"
#include "ODE.H"
#include "ODESystem.H"
#include "ODESolver.H"
#include "RK.H"
@ -37,7 +37,7 @@ using namespace Foam;
class testODE
:
public ODE
public ODESystem
{
public:
@ -107,9 +107,13 @@ int main(int argc, char *argv[])
argList::validArgs.append("ODESolver");
argList args(argc, argv);
// Create the ODE system
testODE ode;
// Create the selected ODE system solver
autoPtr<ODESolver> odeSolver = ODESolver::New(args[1], ode);
// Initialise the ODE system fields
scalar xStart = 1.0;
scalarField yStart(ode.nEqns());
yStart[0] = ::Foam::j0(xStart);
@ -117,6 +121,7 @@ int main(int argc, char *argv[])
yStart[2] = ::Foam::jn(2, xStart);
yStart[3] = ::Foam::jn(3, xStart);
// Print the evolution of the solution and the time-step
scalarField dyStart(ode.nEqns());
ode.derivatives(xStart, yStart, dyStart);

8
src/ODE/ODESolvers/KRR4/KRR4.C
View File

@ -55,7 +55,7 @@ const scalar
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::KRR4::KRR4(const ODE& ode)
Foam::KRR4::KRR4(const ODESystem& ode)
:
ODESolver(ode),
yTemp_(n_, 0.0),
@ -76,7 +76,7 @@ Foam::KRR4::KRR4(const ODE& ode)
void Foam::KRR4::solve
(
const ODE& ode,
const ODESystem& ode,
scalar& x,
scalarField& y,
scalarField& dydx,
@ -168,7 +168,7 @@ void Foam::KRR4::solve
(
"void Foam::KRR4::solve"
"("
"const ODE&, "
"const ODESystem&, "
"scalar&, "
"scalarField&, "
"scalarField&, "
@ -206,7 +206,7 @@ void Foam::KRR4::solve
(
"void Foam::KRR4::solve"
"("
"const ODE&, "
"const ODESystem&, "
"scalar&, "
"scalarField&, "
"scalarField&, "

4
src/ODE/ODESolvers/KRR4/KRR4.H
View File

@ -88,14 +88,14 @@ public:
// Constructors
//- Construct from ODE
KRR4(const ODE& ode);
KRR4(const ODESystem& ode);
// Member Functions
void solve
(
const ODE& ode,
const ODESystem& ode,
scalar& x,
scalarField& y,
scalarField& dydx,

6
src/ODE/ODESolvers/ODESolver/ODESolver.C
View File

@ -36,7 +36,7 @@ namespace Foam
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::ODESolver::ODESolver(const ODE& ode)
Foam::ODESolver::ODESolver(const ODESystem& ode)
:
n_(ode.nEqns()),
yScale_(n_),
@ -48,7 +48,7 @@ Foam::ODESolver::ODESolver(const ODE& ode)
void Foam::ODESolver::solve
(
const ODE& ode,
const ODESystem& ode,
const scalar xStart,
const scalar xEnd,
scalarField& y,
@ -102,7 +102,7 @@ void Foam::ODESolver::solve
FatalErrorIn
(
"ODESolver::solve"
"(const ODE& ode, const scalar xStart, const scalar xEnd,"
"(const ODESystem& ode, const scalar xStart, const scalar xEnd,"
"scalarField& yStart, const scalar eps, scalar& hEst) const"
) << "Too many integration steps"
<< exit(FatalError);

16
src/ODE/ODESolvers/ODESolver/ODESolver.H
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-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -25,7 +25,7 @@ Class
Foam::ODESolver
Description
Selection for ODE solver
Abstract base-class for ODE system solvers
SourceFiles
ODESolver.C
@ -35,7 +35,7 @@ SourceFiles
#ifndef ODESolver_H
#define ODESolver_H
#include "ODE.H"
#include "ODESystem.H"
#include "typeInfo.H"
#include "autoPtr.H"
@ -82,7 +82,7 @@ public:
autoPtr,
ODESolver,
ODE,
(const ODE& ode),
(const ODESystem& ode),
(ode)
);
@ -90,7 +90,7 @@ public:
// Constructors
//- Construct for given ODE
ODESolver(const ODE& ode);
ODESolver(const ODESystem& ode);
// Selectors
@ -99,7 +99,7 @@ public:
static autoPtr<ODESolver> New
(
const word& ODESolverTypeName,
const ODE& ode
const ODESystem& ode
);
@ -112,7 +112,7 @@ public:
virtual void solve
(
const ODE& ode,
const ODESystem& ode,
scalar& x,
scalarField& y,
scalarField& dydx,
@ -126,7 +126,7 @@ public:
virtual void solve
(
const ODE& ode,
const ODESystem& ode,
const scalar xStart,
const scalar xEnd,
scalarField& y,

7
src/ODE/ODESolvers/ODESolver/ODESolverNew.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-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -30,7 +30,7 @@ License
Foam::autoPtr<Foam::ODESolver> Foam::ODESolver::New
(
const Foam::word& ODESolverTypeName,
const Foam::ODE& ode
const Foam::ODESystem& ode
)
{
Info<< "Selecting ODE solver " << ODESolverTypeName << endl;
@ -42,7 +42,8 @@ Foam::autoPtr<Foam::ODESolver> Foam::ODESolver::New
{
FatalErrorIn
(
"ODESolver::New(const word& ODESolverTypeName, const ODE& ode)"
"ODESolver::New"
"(const word& ODESolverTypeName, const ODESystem& ode)"
) << "Unknown ODESolver type "
<< ODESolverTypeName << nl << nl
<< "Valid ODESolvers are : " << endl

6
src/ODE/ODESolvers/RK/RK.C
View File

@ -54,7 +54,7 @@ const scalar
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::RK::RK(const ODE& ode)
Foam::RK::RK(const ODESystem& ode)
:
ODESolver(ode),
yTemp_(n_, 0.0),
@ -72,7 +72,7 @@ Foam::RK::RK(const ODE& ode)
void Foam::RK::solve
(
const ODE& ode,
const ODESystem& ode,
const scalar x,
const scalarField& y,
const scalarField& dydx,
@ -142,7 +142,7 @@ void Foam::RK::solve
void Foam::RK::solve
(
const ODE& ode,
const ODESystem& ode,
scalar& x,
scalarField& y,
scalarField& dydx,

6
src/ODE/ODESolvers/RK/RK.H
View File

@ -81,14 +81,14 @@ public:
// Constructors
//- Construct from ODE
RK(const ODE& ode);
RK(const ODESystem& ode);
// Member Functions
void solve
(
const ODE& ode,
const ODESystem& ode,
const scalar x,
const scalarField& y,
const scalarField& dydx,
@ -100,7 +100,7 @@ public:
void solve
(
const ODE& ode,
const ODESystem& ode,
scalar& x,
scalarField& y,
scalarField& dydx,

4
src/ODE/ODESolvers/SIBS/SIBS.C
View File

@ -47,7 +47,7 @@ namespace Foam
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::SIBS::SIBS(const ODE& ode)
Foam::SIBS::SIBS(const ODESystem& ode)
:
ODESolver(ode),
a_(iMaxX_, 0.0),
@ -70,7 +70,7 @@ Foam::SIBS::SIBS(const ODE& ode)
void Foam::SIBS::solve
(
const ODE& ode,
const ODESystem& ode,
scalar& x,
scalarField& y,
scalarField& dydx,

6
src/ODE/ODESolvers/SIBS/SIBS.H
View File

@ -78,7 +78,7 @@ class SIBS
void SIMPR
(
const ODE& ode,
const ODESystem& ode,
const scalar xStart,
const scalarField& y,
const scalarField& dydx,
@ -110,14 +110,14 @@ public:
// Constructors
//- Construct from ODE
SIBS(const ODE& ode);
SIBS(const ODESystem& ode);
// Member Functions
void solve
(
const ODE& ode,
const ODESystem& ode,
scalar& x,
scalarField& y,
scalarField& dydx,

2
src/ODE/ODESolvers/SIBS/SIMPR.C
View File

@ -29,7 +29,7 @@ License
void Foam::SIBS::SIMPR
(
const ODE& ode,
const ODESystem& ode,
const scalar xStart,
const scalarField& y,
const scalarField& dydx,

25
src/ODE/ODE/ODE.H → src/ODE/ODESystem/ODESystem.H
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-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -22,15 +22,15 @@ License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::ODE
Foam::ODESystem
Description
Abstract base class for the ODE solvers.
Abstract base class for the systems of ordinary differential equations.
\*---------------------------------------------------------------------------*/
#ifndef ODE_H
#define ODE_H
#ifndef ODESystem_H
#define ODESystem_H
#include "scalarField.H"
#include "scalarMatrices.H"
@ -41,30 +41,32 @@ namespace Foam
{
/*---------------------------------------------------------------------------*\
Class ODE Declaration
Class ODESystem Declaration
\*---------------------------------------------------------------------------*/
class ODE
class ODESystem
{
public:
// Constructors
//- Construct null
ODE()
//- Construct null8
ODESystem()
{}
//- Destructor
virtual ~ODE()
virtual ~ODESystem()
{}
// Member Functions
//- Return the number of equations in the system
virtual label nEqns() const = 0;
//- Calculate the derivatives in dydx
virtual void derivatives
(
const scalar x,
@ -72,7 +74,8 @@ public:
scalarField& dydx
) const = 0;
//- Calculate the Jacobian of the system
// Need by the stiff-system solvers
virtual void jacobian
(
const scalar x,

2
src/OpenFOAM/primitives/Tuple2/Tuple2.H
View File

@ -25,7 +25,7 @@ Class
Foam::Tuple2
Description
A 2-tuple.
A 2-tuple for storing two objects of different types.
SeeAlso
Foam::Pair for storing two objects of identical types.

28
src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyDisplacementPointPatchVectorField.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -53,7 +53,7 @@ sixDoFRigidBodyDisplacementPointPatchVectorField
rhoName_("rho"),
lookupGravity_(-1),
g_(vector::zero),
relaxationFactor_(1)
curTimeIndex_(-1)
{}
@ -71,7 +71,7 @@ sixDoFRigidBodyDisplacementPointPatchVectorField
rhoName_(dict.lookupOrDefault<word>("rhoName", "rho")),
lookupGravity_(-1),
g_(vector::zero),
relaxationFactor_(dict.lookupOrDefault<scalar>("relaxationFactor", 1))
curTimeIndex_(-1)
{
if (rhoName_ == "rhoInf")
{
@ -115,7 +115,7 @@ sixDoFRigidBodyDisplacementPointPatchVectorField
rhoName_(ptf.rhoName_),
lookupGravity_(ptf.lookupGravity_),
g_(ptf.g_),
relaxationFactor_(ptf.relaxationFactor_)
curTimeIndex_(-1)
{}
@ -133,7 +133,7 @@ sixDoFRigidBodyDisplacementPointPatchVectorField
rhoName_(ptf.rhoName_),
lookupGravity_(ptf.lookupGravity_),
g_(ptf.g_),
relaxationFactor_(ptf.relaxationFactor_)
curTimeIndex_(-1)
{}
@ -203,6 +203,13 @@ void sixDoFRigidBodyDisplacementPointPatchVectorField::updateCoeffs()
const Time& t = mesh.time();
const pointPatch& ptPatch = this->patch();
// Store the motion state at the beginning of the time-step
if (curTimeIndex_ != this->db().time().timeIndex())
{
motion_.newTime();
curTimeIndex_ = this->db().time().timeIndex();
}
// Patch force data is valid for the current positions, so
// calculate the forces on the motion object from this data, then
// update the positions
@ -231,7 +238,7 @@ void sixDoFRigidBodyDisplacementPointPatchVectorField::updateCoeffs()
g_ = g.value();
}
motion_.updateForce
motion_.updateAcceleration
(
f.forceEff() + g_*motion_.mass(),
f.momentEff(),
@ -240,11 +247,7 @@ void sixDoFRigidBodyDisplacementPointPatchVectorField::updateCoeffs()
Field<vector>::operator=
(
relaxationFactor_
*(
motion_.currentPosition(initialPoints_)
- initialPoints_
)
motion_.currentPosition(initialPoints_) - initialPoints_
);
fixedValuePointPatchField<vector>::updateCoeffs();
@ -267,9 +270,6 @@ void sixDoFRigidBodyDisplacementPointPatchVectorField::write(Ostream& os) const
os.writeKeyword("g") << g_ << token::END_STATEMENT << nl;
}
os.writeKeyword("relaxationFactor")
<< relaxationFactor_ << token::END_STATEMENT << nl;
motion_.write(os);
initialPoints_.writeEntry("initialPoints", os);

6
src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyDisplacementPointPatchVectorField.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -84,8 +84,8 @@ class sixDoFRigidBodyDisplacementPointPatchVectorField
//- Gravity vector to store when not available from the db
vector g_;
//- Optional under-relaxation factor for the motion
scalar relaxationFactor_;
//- Current time index (used for updating)
label curTimeIndex_;
public:

147
src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.C
View File

@ -164,6 +164,7 @@ void Foam::sixDoFRigidBodyMotion::applyConstraints(scalar deltaT)
Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion()
:
motionState_(),
motionState0_(),
restraints_(),
restraintNames_(),
constraints_(),
@ -173,8 +174,7 @@ Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion()
initialQ_(I),
momentOfInertia_(diagTensor::one*VSMALL),
mass_(VSMALL),
cDamp_(0.0),
aLim_(VGREAT),
aRelax_(1.0),
report_(false)
{}
@ -191,8 +191,7 @@ Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion
const point& initialCentreOfMass,
const tensor& initialQ,
const diagTensor& momentOfInertia,
scalar cDamp,
scalar aLim,
scalar aRelax,
bool report
)
:
@ -205,6 +204,7 @@ Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion
pi,
tau
),
motionState0_(motionState_),
restraints_(),
restraintNames_(),
constraints_(),
@ -214,8 +214,7 @@ Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion
initialQ_(initialQ),
momentOfInertia_(momentOfInertia),
mass_(mass),
cDamp_(cDamp),
aLim_(aLim),
aRelax_(aRelax),
report_(report)
{}
@ -223,6 +222,7 @@ Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion
Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion(const dictionary& dict)
:
motionState_(dict),
motionState0_(motionState_),
restraints_(),
restraintNames_(),
constraints_(),
@ -238,8 +238,7 @@ 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)),
aRelax_(dict.lookupOrDefault<scalar>("accelerationRelaxation", 1.0)),
report_(dict.lookupOrDefault<Switch>("report", false))
{
addRestraints(dict);
@ -254,6 +253,7 @@ Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion
)
:
motionState_(sDoFRBM.motionState_),
motionState0_(sDoFRBM.motionState0_),
restraints_(sDoFRBM.restraints_),
restraintNames_(sDoFRBM.restraintNames_),
constraints_(sDoFRBM.constraints_),
@ -263,8 +263,7 @@ Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion
initialQ_(sDoFRBM.initialQ_),
momentOfInertia_(sDoFRBM.momentOfInertia_),
mass_(sDoFRBM.mass_),
cDamp_(sDoFRBM.cDamp_),
aLim_(sDoFRBM.aLim_),
aRelax_(sDoFRBM.aRelax_),
report_(sDoFRBM.report_)
{}
@ -376,91 +375,59 @@ void Foam::sixDoFRigidBodyMotion::updatePosition
if (Pstream::master())
{
vector aClip = a();
scalar aMag = mag(aClip);
if (aMag > SMALL)
{
aClip /= aMag;
}
if (aMag > aLim_)
{
WarningIn
(
"void Foam::sixDoFRigidBodyMotion::updatePosition"
"("
"scalar, "
"scalar"
")"
)
<< "Limited acceleration " << a()
<< " to " << aClip*aLim_
<< endl;
}
v() += 0.5*(1 - cDamp_)*deltaT0*aClip*min(aMag, aLim_);
pi() += 0.5*(1 - cDamp_)*deltaT0*tau();
v() = v0() + 0.5*deltaT0*a();
pi() = pi0() + 0.5*deltaT0*tau();
// Leapfrog move part
centreOfMass() += deltaT*v();
centreOfMass() = centreOfMass0() + deltaT*v();
// Leapfrog orientation adjustment
rotate(Q(), pi(), deltaT);
Tuple2<tensor, vector> Qpi = rotate(Q0(), pi(), deltaT);
Q() = Qpi.first();
pi() = Qpi.second();
}
Pstream::scatter(motionState_);
}
void Foam::sixDoFRigidBodyMotion::updateForce
void Foam::sixDoFRigidBodyMotion::updateAcceleration
(
const vector& fGlobal,
const vector& tauGlobal,
scalar deltaT
)
{
static bool first = true;
// Second leapfrog velocity adjust part, required after motion and
// force calculation
// acceleration calculation
if (Pstream::master())
{
// Save the previous iteration accelerations for relaxation
vector aPrevIter = a();
vector tauPrevIter = tau();
// Calculate new accelerations
a() = fGlobal/mass_;
tau() = (Q().T() & tauGlobal);
applyRestraints();
// Relax accelerations on all but first iteration
if (!first)
{
a() = aRelax_*a() + (1 - aRelax_)*aPrevIter;
tau() = aRelax_*tau() + (1 - aRelax_)*tauPrevIter;
}
first = false;
// Apply constraints after relaxation
applyConstraints(deltaT);
vector aClip = a();
scalar aMag = mag(aClip);
if (aMag > SMALL)
{
aClip /= aMag;
}
if (aMag > aLim_)
{
WarningIn
(
"void Foam::sixDoFRigidBodyMotion::updateForce"
"("
"const vector&, "
"const vector&, "
"scalar"
")"
)
<< "Limited acceleration " << a()
<< " to " << aClip*aLim_
<< endl;
}
v() += 0.5*(1 - cDamp_)*deltaT*aClip*min(aMag, aLim_);
pi() += 0.5*(1 - cDamp_)*deltaT*tau();
// Correct velocities
v() += 0.5*deltaT*a();
pi() += 0.5*deltaT*tau();
if (report_)
{
@ -472,7 +439,27 @@ void Foam::sixDoFRigidBodyMotion::updateForce
}
void Foam::sixDoFRigidBodyMotion::updateForce
void Foam::sixDoFRigidBodyMotion::updateVelocity(scalar deltaT)
{
// Second leapfrog velocity adjust part, required after motion and
// acceleration calculation
if (Pstream::master())
{
v() += 0.5*deltaT*a();
pi() += 0.5*deltaT*tau();
if (report_)
{
status();
}
}
Pstream::scatter(motionState_);
}
void Foam::sixDoFRigidBodyMotion::updateAcceleration
(
const pointField& positions,
const vectorField& forces,
@ -493,7 +480,7 @@ void Foam::sixDoFRigidBodyMotion::updateForce
}
}
updateForce(fGlobal, tauGlobal, deltaT);
updateAcceleration(fGlobal, tauGlobal, deltaT);
}
@ -506,19 +493,14 @@ Foam::point Foam::sixDoFRigidBodyMotion::predictedPosition
) const
{
vector vTemp = v() + deltaT*(a() + deltaForce/mass_);
vector piTemp = pi() + deltaT*(tau() + (Q().T() & deltaMoment));
point centreOfMassTemp = centreOfMass() + deltaT*vTemp;
tensor QTemp = Q();
rotate(QTemp, piTemp, deltaT);
point centreOfMassTemp = centreOfMass0() + deltaT*vTemp;
Tuple2<tensor, vector> QpiTemp = rotate(Q0(), piTemp, deltaT);
return
(
centreOfMassTemp
+ (QTemp & initialQ_.T() & (pInitial - initialCentreOfMass_))
+ (QpiTemp.first() & initialQ_.T() & (pInitial - initialCentreOfMass_))
);
}
@ -531,12 +513,9 @@ Foam::vector Foam::sixDoFRigidBodyMotion::predictedOrientation
) const
{
vector piTemp = pi() + deltaT*(tau() + (Q().T() & deltaMoment));
Tuple2<tensor, vector> QpiTemp = rotate(Q0(), piTemp, deltaT);
tensor QTemp = Q();
rotate(QTemp, piTemp, deltaT);
return (QTemp & initialQ_.T() & vInitial);
return (QpiTemp.first() & initialQ_.T() & vInitial);
}

64
src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.H
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-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -59,7 +59,7 @@ SourceFiles
#include "pointField.H"
#include "sixDoFRigidBodyMotionRestraint.H"
#include "sixDoFRigidBodyMotionConstraint.H"
#include "Tuple2.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -87,6 +87,9 @@ class sixDoFRigidBodyMotion
//- Motion state data object
sixDoFRigidBodyMotionState motionState_;
//- Motion state data object for previous time-step
sixDoFRigidBodyMotionState motionState0_;
//- Restraints on the motion
PtrList<sixDoFRigidBodyMotionRestraint> restraints_;
@ -116,14 +119,8 @@ 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_;
//- Acceleration relaxation coefficient
scalar aRelax_;
//- Switch to turn reporting of motion data on and off
Switch report_;
@ -143,8 +140,14 @@ class sixDoFRigidBodyMotion
// frame z-axis by the given angle
inline tensor rotationTensorZ(scalar deltaT) const;
//- Apply rotation tensors to Q for the given torque (pi) and deltaT
inline void rotate(tensor& Q, vector& pi, scalar deltaT) const;
//- Apply rotation tensors to Q0 for the given torque (pi) and deltaT
// and return the rotated Q and pi as a tuple
inline Tuple2<tensor, vector> rotate
(
const tensor& Q0,
const vector& pi,
const scalar deltaT
) const;
//- Apply the restraints to the object
void applyRestraints();
@ -152,6 +155,7 @@ class sixDoFRigidBodyMotion
//- Apply the constraints to the object
void applyConstraints(scalar deltaT);
// Access functions retained as private because of the risk of
// confusion over what is a body local frame vector and what is global
@ -199,6 +203,21 @@ class sixDoFRigidBodyMotion
//- Return access to torque
inline const vector& tau() const;
//- Return access to the orientation at previous time-step
inline const tensor& Q0() const;
//- Return access to velocity at previous time-step
inline const vector& v0() const;
//- Return access to acceleration at previous time-step
inline const vector& a0() const;
//- Return access to angular momentum at previous time-step
inline const vector& pi0() const;
//- Return access to torque at previous time-step
inline const vector& tau0() const;
// Edit
@ -244,8 +263,7 @@ public:
const point& initialCentreOfMass,
const tensor& initialQ,
const diagTensor& momentOfInertia,
scalar cDamp = 0.0,
scalar aLim = VGREAT,
scalar aRelax = 1.0,
bool report = false
);
@ -279,18 +297,22 @@ public:
scalar deltaT0
);
//- Second leapfrog velocity adjust part, required after motion and
// force calculation
void updateForce
//- Second leapfrog velocity adjust part
// required after motion and force calculation
void updateAcceleration
(
const vector& fGlobal,
const vector& tauGlobal,
scalar deltaT
);
//- Second leapfrog velocity adjust part
// required after motion and force calculation
void updateVelocity(scalar deltaT);
//- Global forces supplied at locations, calculating net force
// and moment
void updateForce
void updateAcceleration
(
const pointField& positions,
const vectorField& forces,
@ -354,6 +376,9 @@ public:
//- Return const access to the centre of mass
inline const point& centreOfMass() const;
//- Return const access to the centre of mass at previous time-step
inline const point& centreOfMass0() const;
//- Return access to the inertia tensor
inline const diagTensor& momentOfInertia() const;
@ -366,6 +391,9 @@ public:
// Edit
//- Store the motion state at the beginning of the time-step
inline void newTime();
//- Return non-const access to the centre of mass
inline point& centreOfMass();

61
src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionI.H
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-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -61,16 +61,19 @@ Foam::sixDoFRigidBodyMotion::rotationTensorZ(scalar phi) const
}
inline void Foam::sixDoFRigidBodyMotion::rotate
inline Foam::Tuple2<Foam::tensor, Foam::vector>
Foam::sixDoFRigidBodyMotion::rotate
(
tensor& Q,
vector& pi,
scalar deltaT
const tensor& Q0,
const vector& pi0,
const scalar deltaT
) const
{
tensor R;
Tuple2<tensor, vector> Qpi(Q0, pi0);
tensor& Q = Qpi.first();
vector& pi = Qpi.second();
R = rotationTensorX(0.5*deltaT*pi.x()/momentOfInertia_.xx());
tensor R = rotationTensorX(0.5*deltaT*pi.x()/momentOfInertia_.xx());
pi = pi & R;
Q = Q & R;
@ -89,6 +92,8 @@ inline void Foam::sixDoFRigidBodyMotion::rotate
R = rotationTensorX(0.5*deltaT*pi.x()/momentOfInertia_.xx());
pi = pi & R;
Q = Q & R;
return Qpi;
}
@ -176,6 +181,36 @@ inline const Foam::vector& Foam::sixDoFRigidBodyMotion::tau() const
}
inline const Foam::tensor& Foam::sixDoFRigidBodyMotion::Q0() const
{
return motionState0_.Q();
}
inline const Foam::vector& Foam::sixDoFRigidBodyMotion::v0() const
{
return motionState0_.v();
}
inline const Foam::vector& Foam::sixDoFRigidBodyMotion::a0() const
{
return motionState0_.a();
}
inline const Foam::vector& Foam::sixDoFRigidBodyMotion::pi0() const
{
return motionState0_.pi();
}
inline const Foam::vector& Foam::sixDoFRigidBodyMotion::tau0() const
{
return motionState0_.tau();
}
inline Foam::point& Foam::sixDoFRigidBodyMotion::initialCentreOfMass()
{
return initialCentreOfMass_;
@ -281,6 +316,12 @@ inline const Foam::point& Foam::sixDoFRigidBodyMotion::centreOfMass() const
}
inline const Foam::point& Foam::sixDoFRigidBodyMotion::centreOfMass0() const
{
return motionState0_.centreOfMass();
}
inline const Foam::diagTensor&
Foam::sixDoFRigidBodyMotion::momentOfInertia() const
{
@ -300,6 +341,12 @@ inline bool Foam::sixDoFRigidBodyMotion::report() const
}
inline void Foam::sixDoFRigidBodyMotion::newTime()
{
motionState0_ = motionState_;
}
inline Foam::point& Foam::sixDoFRigidBodyMotion::centreOfMass()
{
return motionState_.centreOfMass();

8
src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionIO.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-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -40,10 +40,8 @@ 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("accelerationRelaxation")
<< aRelax_ << token::END_STATEMENT << nl;
os.writeKeyword("report")
<< report_ << token::END_STATEMENT << nl;

8
src/postProcessing/functionObjects/forces/pointPatchFields/derived/uncoupledSixDoFRigidBodyDisplacement/uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField.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-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -160,7 +160,11 @@ void uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField::updateCoeffs()
}
// Do not modify the accelerations, except with gravity, where available
motion_.updateForce(gravity*motion_.mass(), vector::zero, t.deltaTValue());
motion_.updateAcceleration
(
gravity*motion_.mass(),
vector::zero, t.deltaTValue()
);
Field<vector>::operator=
(

2
src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.C
View File

@ -36,7 +36,7 @@ Foam::chemistryModel<CompType, ThermoType>::chemistryModel
)
:
CompType(mesh),
ODE(),
ODESystem(),
Y_(this->thermo().composition().Y()),
reactions_
(

4
src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.H
View File

@ -39,7 +39,7 @@ SourceFiles
#define chemistryModel_H
#include "Reaction.H"
#include "ODE.H"
#include "ODESystem.H"
#include "volFieldsFwd.H"
#include "simpleMatrix.H"
#include "DimensionedField.H"
@ -60,7 +60,7 @@ template<class CompType, class ThermoType>
class chemistryModel
:
public CompType,
public ODE
public ODESystem
{
// Private Member Functions

2
src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModel.C
View File

@ -37,7 +37,7 @@ solidChemistryModel
)
:
CompType(mesh),
ODE(),
ODESystem(),
Ys_(this->solidThermo().composition().Y()),
reactions_
(

4
src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModel.H
View File

@ -40,7 +40,7 @@ SourceFiles
#define solidChemistryModel_H
#include "Reaction.H"
#include "ODE.H"
#include "ODESystem.H"
#include "volFieldsFwd.H"
#include "DimensionedField.H"
#include "simpleMatrix.H"
@ -61,7 +61,7 @@ template<class CompType, class SolidThermo>
class solidChemistryModel
:
public CompType,
public ODE
public ODESystem
{
// Private Member Functions

10
tutorials/multiphase/interDyMFoam/ras/floatingObject/system/controlDict
View File

@ -17,7 +17,7 @@ FoamFile
application interDyMFoam;
startFrom latestTime;
startFrom startTime;
startTime 0;
@ -29,7 +29,7 @@ deltaT 0.01;
writeControl adjustableRunTime;
writeInterval 0.1;
writeInterval 0.2;
purgeWrite 0;
@ -47,9 +47,9 @@ runTimeModifiable yes;
adjustTimeStep yes;
maxCo 0.5;
maxAlphaCo 0.5;
maxDeltaT 0.01;
maxCo 5;
maxAlphaCo 5;
maxDeltaT 1;
libs
(

4
tutorials/multiphase/interDyMFoam/ras/floatingObject/system/fvSchemes
View File

@ -27,7 +27,7 @@ gradSchemes
divSchemes
{
div(rho*phi,U) Gauss vanLeerV;
div(rhoPhi,U) Gauss vanLeerV;
div(phi,alpha) Gauss vanLeer;
div(phirb,alpha) Gauss interfaceCompression;
div(phi,k) Gauss upwind;
@ -55,7 +55,7 @@ fluxRequired
default no;
p_rgh;
pcorr;
alpha;
alpha.water;
}

24
tutorials/multiphase/interDyMFoam/ras/floatingObject/system/fvSolution
View File

@ -17,7 +17,7 @@ FoamFile
solvers
{
cellDisplacement
"cellDisplacement.*"
{
solver GAMG;
tolerance 1e-5;
@ -29,14 +29,24 @@ solvers
mergeLevels 1;
}
alpha.water
"alpha.water.*"
{
nAlphaCorr 1;
nAlphaSubCycles 3;
nAlphaSubCycles 1;
cAlpha 1;
alphaOuterCorrectors yes;
MULESCorr yes;
nLimiterIter 5;
solver PBiCG;
preconditioner DILU;
tolerance 1e-8;
relTol 0;
}
pcorr
"pcorr.*"
{
solver PCG;
preconditioner
@ -120,9 +130,11 @@ solvers
PIMPLE
{
momentumPredictor no;
nCorrectors 2;
nOuterCorrectors 6;
nCorrectors 1;
nNonOrthogonalCorrectors 0;
correctPhi yes;
moveMeshOuterCorrectors yes;
}
relaxationFactors
@ -132,7 +144,7 @@ relaxationFactors
}
equations
{
"(U|k|epsilon|omega|R|nuTilda).*" 1;
".*" 1;
}
}