songtianlun/openfoam
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

6-DoF solid body motion: Support PIMPLE iteration and acceleration relaxation

Browse Source
This commit is contained in:
Henry 2013-10-31 10:22:26 +00:00
parent 30db0da817
commit 6df8f705ae
12 changed files with 212 additions and 143 deletions
Show Stats 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
(

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.

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

@ -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);

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

@ -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();

59
src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionI.H
View File

@ -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();

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

@ -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;

6
src/postProcessing/functionObjects/forces/pointPatchFields/derived/uncoupledSixDoFRigidBodyDisplacement/uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField.C
View File

@ -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=
(

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;
}
}