Merge branch 'master' of ssh://noisy/home/noisy3/OpenFOAM/OpenFOAM-dev

src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyDisplacementPointPatchVectorField.C

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

src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.C

@@ -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()),
+    motionState_(sDoFRBM.motionState_),
-    restraints_(sDoFRBM.restraints()),
+    restraints_(sDoFRBM.restraints_),
-    restraintNames_(sDoFRBM.restraintNames()),
+    restraintNames_(sDoFRBM.restraintNames_),
-    constraints_(sDoFRBM.constraints()),
+    constraints_(sDoFRBM.constraints_),
-    constraintNames_(sDoFRBM.constraintNames()),
+    constraintNames_(sDoFRBM.constraintNames_),
-    maxConstraintIterations_(sDoFRBM.maxConstraintIterations()),
+    maxConstraintIterations_(sDoFRBM.maxConstraintIterations_),
-    initialCentreOfMass_(sDoFRBM.initialCentreOfMass()),
+    initialCentreOfMass_(sDoFRBM.initialCentreOfMass_),
-    initialQ_(sDoFRBM.initialQ()),
+    initialQ_(sDoFRBM.initialQ_),
-    momentOfInertia_(sDoFRBM.momentOfInertia()),
+    momentOfInertia_(sDoFRBM.momentOfInertia_),
-    mass_(sDoFRBM.mass()),
+    mass_(sDoFRBM.mass_),
-    report_(sDoFRBM.report())
+    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_)
         {

src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.H

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

src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionIO.C

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

src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/linearAxialAngularSpring/linearAxialAngularSpring.C

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

src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/tabulatedAxialAngularSpring/tabulatedAxialAngularSpring.C

@@ -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;
+        vector oldDir = refQ_ & refDir;
-        newDir = motion.orientation() & refDir;
+
+        vector newDir = motion.orientation() & refDir;
     }
 
     // Removing any axis component from oldDir and newDir and normalising

src/postProcessing/functionObjects/forces/pointPatchFields/derived/uncoupledSixDoFRigidBodyDisplacement/uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField.C

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