diff --git a/src/sixDoFRigidBodyMotion/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyDisplacementPointPatchVectorField.C b/src/sixDoFRigidBodyMotion/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyDisplacementPointPatchVectorField.C
index dc731f060a39cef715c1c43e1299302974eac315..a22e4e586551af465a9fc12833c7112ca161d613 100644
--- a/src/sixDoFRigidBodyMotion/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyDisplacementPointPatchVectorField.C
+++ b/src/sixDoFRigidBodyMotion/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyDisplacementPointPatchVectorField.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -204,18 +204,14 @@ void sixDoFRigidBodyDisplacementPointPatchVectorField::updateCoeffs()
const pointPatch& ptPatch = this->patch();
// Store the motion state at the beginning of the time-step
+ bool firstIter = false;
if (curTimeIndex_ != t.timeIndex())
{
motion_.newTime();
curTimeIndex_ = t.timeIndex();
+ firstIter = true;
}
- // Patch force data is valid for the current positions, so
- // calculate the forces on the motion object from this data, then
- // update the positions
-
- motion_.updatePosition(t.deltaTValue(), t.deltaT0Value());
-
dictionary forcesDict;
forcesDict.add("type", forces::typeName);
@@ -228,6 +224,11 @@ void sixDoFRigidBodyDisplacementPointPatchVectorField::updateCoeffs()
f.calcForcesMoment();
+ // Patch force data is valid for the current positions, so
+ // calculate the forces on the motion object from this data, then
+ // update the positions
+ motion_.updatePosition(firstIter, t.deltaTValue(), t.deltaT0Value());
+
// Get the forces on the patch faces at the current positions
if (lookupGravity_ == 1)
@@ -243,9 +244,11 @@ void sixDoFRigidBodyDisplacementPointPatchVectorField::updateCoeffs()
motion_.updateAcceleration
(
+ firstIter,
ramp*(f.forceEff() + motion_.mass()*g_),
ramp*(f.momentEff() + motion_.mass()*(motion_.momentArm() ^ g_)),
- t.deltaTValue()
+ t.deltaTValue(),
+ t.deltaT0Value()
);
Field<vector>::operator=
diff --git a/src/sixDoFRigidBodyMotion/pointPatchFields/derived/uncoupledSixDoFRigidBodyDisplacement/uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField.C b/src/sixDoFRigidBodyMotion/pointPatchFields/derived/uncoupledSixDoFRigidBodyDisplacement/uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField.C
index fb832ca6c80aed91d41ecfd95f26b59dba4987f5..03bdfcd77935d87eb114edc2354c58536cd54ff0 100644
--- a/src/sixDoFRigidBodyMotion/pointPatchFields/derived/uncoupledSixDoFRigidBodyDisplacement/uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField.C
+++ b/src/sixDoFRigidBodyMotion/pointPatchFields/derived/uncoupledSixDoFRigidBodyDisplacement/uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -152,13 +152,15 @@ void uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField::updateCoeffs()
const Time& t = mesh.time();
// Store the motion state at the beginning of the time-step
+ bool firstIter = false;
if (curTimeIndex_ != t.timeIndex())
{
motion_.newTime();
curTimeIndex_ = t.timeIndex();
+ firstIter = true;
}
- motion_.updatePosition(t.deltaTValue(), t.deltaT0Value());
+ motion_.updatePosition(firstIter, t.deltaTValue(), t.deltaT0Value());
vector gravity = vector::zero;
@@ -173,9 +175,11 @@ void uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField::updateCoeffs()
// Do not modify the accelerations, except with gravity, where available
motion_.updateAcceleration
(
+ firstIter,
gravity*motion_.mass(),
vector::zero,
- t.deltaTValue()
+ t.deltaTValue(),
+ t.deltaT0Value()
);
Field<vector>::operator=
diff --git a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.C b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.C
index 8155e218dd74e223dcb8ed591003e258f9cff2ea..9bb33ecbf9f677df7662840903c012a779be9b2f 100644
--- a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.C
+++ b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.C
@@ -259,22 +259,38 @@ void Foam::sixDoFRigidBodyMotion::addConstraints
void Foam::sixDoFRigidBodyMotion::updatePosition
(
+ bool firstIter,
scalar deltaT,
scalar deltaT0
)
{
- // First leapfrog velocity adjust and motion part, required before
- // force calculation
-
if (Pstream::master())
{
- v() = tConstraints_ & (v0() + aDamp_*0.5*deltaT0*a());
- pi() = rConstraints_ & (pi0() + aDamp_*0.5*deltaT0*tau());
+ if (firstIter)
+ {
+ // First simplectic step:
+ // Half-step for linear and angular velocities
+ // Update position and orientation
+
+ v() = tConstraints_ & (v0() + aDamp_*0.5*deltaT0*a());
+ pi() = rConstraints_ & (pi0() + aDamp_*0.5*deltaT0*tau());
- // Leapfrog move part
- centreOfRotation() = centreOfRotation0() + deltaT*v();
+ centreOfRotation() = centreOfRotation0() + deltaT*v();
+ }
+ else
+ {
+ // For subsequent iterations use Crank-Nicolson
- // Leapfrog orientation adjustment
+ v() = tConstraints_
+ & (v0() + aDamp_*0.5*deltaT*(a() + motionState0_.a()));
+ pi() = rConstraints_
+ & (pi0() + aDamp_*0.5*deltaT*(tau() + motionState0_.tau()));
+
+ centreOfRotation() =
+ centreOfRotation0() + 0.5*deltaT*(v() + motionState0_.v());
+ }
+
+ // Correct orientation
Tuple2<tensor, vector> Qpi = rotate(Q0(), pi(), deltaT);
Q() = Qpi.first();
pi() = rConstraints_ & Qpi.second();
@@ -286,16 +302,15 @@ void Foam::sixDoFRigidBodyMotion::updatePosition
void Foam::sixDoFRigidBodyMotion::updateAcceleration
(
+ bool firstIter,
const vector& fGlobal,
const vector& tauGlobal,
- scalar deltaT
+ scalar deltaT,
+ scalar deltaT0
)
{
static bool first = false;
- // Second leapfrog velocity adjust part, required after motion and
- // acceleration calculation
-
if (Pstream::master())
{
// Save the previous iteration accelerations for relaxation
@@ -315,9 +330,23 @@ void Foam::sixDoFRigidBodyMotion::updateAcceleration
}
first = false;
- // Correct velocities
- v() += tConstraints_ & aDamp_*0.5*deltaT*a();
- pi() += rConstraints_ & aDamp_*0.5*deltaT*tau();
+ if (firstIter)
+ {
+ // Second simplectic step:
+ // Complete update of linear and angular velocities
+
+ v() += tConstraints_ & aDamp_*0.5*deltaT*a();
+ pi() += rConstraints_ & aDamp_*0.5*deltaT*tau();
+ }
+ else
+ {
+ // For subsequent iterations use Crank-Nicolson
+
+ v() = tConstraints_
+ & (v0() + aDamp_*0.5*deltaT*(a() + motionState0_.a()));
+ pi() = rConstraints_
+ & (pi0() + aDamp_*0.5*deltaT*(tau() + motionState0_.tau()));
+ }
if (report_)
{
@@ -329,6 +358,7 @@ void Foam::sixDoFRigidBodyMotion::updateAcceleration
}
+
void Foam::sixDoFRigidBodyMotion::status() const
{
Info<< "6-DoF rigid body motion" << nl
diff --git a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.H b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.H
index b08c7d011ff70dfb1d19546ee835fc3d976cbd60..70f2e7674f60b37e885eb6125f2878a40857918b 100644
--- a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.H
+++ b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.H
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -42,6 +42,9 @@ Description
url = {http://link.aip.org/link/?JCP/107/5840/1},
doi = {10.1063/1.474310}
+ Changes to Crank-Nicolson integration if subsequent iterations to handle
+ implicit coupling between pressure and motion.
+
Can add restraints (e.g. a spring)
and constraints (e.g. motion may only be on a plane).
@@ -298,18 +301,22 @@ public:
// Update state
- //- First leapfrog velocity adjust and motion part, required
+ //- First symplectic velocity adjust and motion part, required
// before force calculation. Takes old timestep for variable
// timestep cases.
- void updatePosition(scalar deltaT, scalar deltaT0);
+ // Changes to Crank-Nicolson integration for subsequent iterations.
+ void updatePosition(bool firstIter, scalar deltaT, scalar deltaT0);
- //- Second leapfrog velocity adjust part
- // required after motion and force calculation
+ //- Second symplectic velocity adjust part
+ // required after motion and force calculation.
+ // Changes to Crank-Nicolson integration for subsequent iterations.
void updateAcceleration
(
+ bool firstIter,
const vector& fGlobal,
const vector& tauGlobal,
- scalar deltaT
+ scalar deltaT,
+ scalar deltaT0
);
//- Report the status of the motion
diff --git a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionSolver/sixDoFRigidBodyMotionSolver.C b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionSolver/sixDoFRigidBodyMotionSolver.C
index 15bcc0df9316f9c586f235700317d4a2c0f71e9f..39453bb9d7353d95211e663119552b128b05a7ab 100644
--- a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionSolver/sixDoFRigidBodyMotionSolver.C
+++ b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionSolver/sixDoFRigidBodyMotionSolver.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2013-2014 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2013-2015 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -85,6 +85,7 @@ Foam::sixDoFRigidBodyMotionSolver::sixDoFRigidBodyMotionSolver
patchSet_(mesh.boundaryMesh().patchSet(patches_)),
di_(readScalar(coeffDict().lookup("innerDistance"))),
do_(readScalar(coeffDict().lookup("outerDistance"))),
+ test_(coeffDict().lookupOrDefault<Switch>("test", false)),
rhoInf_(1.0),
rhoName_(coeffDict().lookupOrDefault<word>("rhoName", "rho")),
scale_
@@ -179,31 +180,15 @@ void Foam::sixDoFRigidBodyMotionSolver::solve()
<< " points." << exit(FatalError);
}
- // Store the motion state at the beginning of the time-step
+ // Store the motion state at the beginning of the time-stepbool
+ bool firstIter = false;
if (curTimeIndex_ != this->db().time().timeIndex())
{
motion_.newTime();
curTimeIndex_ = this->db().time().timeIndex();
+ firstIter = true;
}
- // Patch force data is valid for the current positions, so
- // calculate the forces on the motion object from this data, then
- // update the positions
-
- motion_.updatePosition(t.deltaTValue(), t.deltaT0Value());
-
- dictionary forcesDict;
-
- forcesDict.add("type", forces::typeName);
- forcesDict.add("patches", patches_);
- forcesDict.add("rhoInf", rhoInf_);
- forcesDict.add("rhoName", rhoName_);
- forcesDict.add("CofR", motion_.centreOfRotation());
-
- forces f("forces", db(), forcesDict);
-
- f.calcForcesMoment();
-
dimensionedVector g("g", dimAcceleration, vector::zero);
if (db().foundObject<uniformDimensionedVectorField>("g"))
@@ -218,12 +203,51 @@ void Foam::sixDoFRigidBodyMotionSolver::solve()
// scalar ramp = min(max((this->db().time().value() - 5)/10, 0), 1);
scalar ramp = 1.0;
- motion_.updateAcceleration
- (
- ramp*(f.forceEff() + motion_.mass()*g.value()),
- ramp*(f.momentEff() + motion_.mass()*(motion_.momentArm() ^ g.value())),
- t.deltaTValue()
- );
+ if (test_)
+ {
+ motion_.updatePosition(firstIter, t.deltaTValue(), t.deltaT0Value());
+
+ motion_.updateAcceleration
+ (
+ firstIter,
+ ramp*(motion_.mass()*g.value()),
+ ramp*(motion_.mass()*(motion_.momentArm() ^ g.value())),
+ t.deltaTValue(),
+ t.deltaT0Value()
+ );
+ }
+ else
+ {
+ dictionary forcesDict;
+
+ forcesDict.add("type", forces::typeName);
+ forcesDict.add("patches", patches_);
+ forcesDict.add("rhoInf", rhoInf_);
+ forcesDict.add("rhoName", rhoName_);
+ forcesDict.add("CofR", motion_.centreOfRotation());
+
+ forces f("forces", db(), forcesDict);
+
+ f.calcForcesMoment();
+
+ // Patch force data is valid for the current positions, so
+ // calculate the forces on the motion object from this data, then
+ // update the positions
+ motion_.updatePosition(firstIter, t.deltaTValue(), t.deltaT0Value());
+
+ motion_.updateAcceleration
+ (
+ firstIter,
+ ramp*(f.forceEff() + motion_.mass()*g.value()),
+ ramp
+ *(
+ f.momentEff()
+ + motion_.mass()*(motion_.momentArm() ^ g.value())
+ ),
+ t.deltaTValue(),
+ t.deltaT0Value()
+ );
+ }
// Update the displacements
pointDisplacement_.internalField() =
diff --git a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionSolver/sixDoFRigidBodyMotionSolver.H b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionSolver/sixDoFRigidBodyMotionSolver.H
index 4f4e1ca31c5a80e14eaee7ab06f00cbfa8ae8685..5bf52d69cfd7805aecba3fdfd0b829570d119368 100644
--- a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionSolver/sixDoFRigidBodyMotionSolver.H
+++ b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionSolver/sixDoFRigidBodyMotionSolver.H
@@ -70,6 +70,10 @@ class sixDoFRigidBodyMotionSolver
//- Outer morphing distance (limit of linear interpolation region)
const scalar do_;
+ //- Switch for test-mode in which only the
+ // gravitational body-force is applied
+ Switch test_;
+
//- Reference density required by the forces object for
// incompressible calculations, required if rhoName == rhoInf
scalar rhoInf_;