rigidBodyDynamics: Added support for restraints and a linear spring with damper

applications/test/rigidBodyDynamics/spring: Test of the linear spring with damper restraint
Damped simple harmonic motion of a weight on a spring is simulated and
the results compared with analytical solution

    Test-spring
    gnuplot spring.gnuplot
    evince spring.eps

This development is sponsored by Carnegie Wave Energy Ltd.

applications/test/rigidBodyDynamics/spring/Make/files

@@ -0,0 +1,3 @@
+spring.C
+
+EXE = $(FOAM_USER_APPBIN)/Test-spring

applications/test/rigidBodyDynamics/spring/Make/options

@@ -0,0 +1,5 @@
+EXE_INC = \
+    -I$(LIB_SRC)/rigidBodyDynamics/lnInclude
+
+EXE_LIBS = \
+    -lrigidBodyDynamics

applications/test/rigidBodyDynamics/spring/spring

@@ -0,0 +1,30 @@
+bodies
+{
+    weight
+    {
+        type            rigidBody;
+        mass            9.6;
+        centreOfMass    (0 0 0);
+        inertia         (0.1052 0 0 0.1052 0 0.1778);
+        parent          root;
+        transform       (1 0 0 0 1 0 0 0 1) (0 0 0);
+        joint
+        {
+            type            Pz;
+        }
+    }
+}
+
+restraints
+{
+    spring
+    {
+        type            linearSpring;
+        body            weight;
+        anchor          (0 0 0.5);
+        refAttachmentPt (0 0 0);
+        stiffness       5000;
+        damping         50;
+        restLength      0.4;
+    }
+}

applications/test/rigidBodyDynamics/spring/spring.C

@@ -0,0 +1,109 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    spring
+
+Description
+    Simple weight and damped-spirng simulation with 1-DoF.
+
+\*---------------------------------------------------------------------------*/
+
+#include "rigidBodyModel.H"
+#include "masslessBody.H"
+#include "sphere.H"
+#include "joints.H"
+#include "IFstream.H"
+#include "OFstream.H"
+
+using namespace Foam;
+using namespace RBD;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    // Create the spring model from dictionary
+    rigidBodyModel spring(dictionary(IFstream("spring")()));
+
+    Info<< spring << endl;
+
+    // Create the joint-space state fields
+    scalarField q(spring.nDoF(), Zero);
+    scalarField w(spring.nw(), Zero);
+    scalarField qDot(spring.nDoF(), Zero);
+    scalarField qDdot(spring.nDoF(), Zero);
+    scalarField tau(spring.nDoF(), Zero);
+    Field<spatialVector> fx(spring.nBodies(), Zero);
+
+    OFstream qFile("qVsTime");
+    OFstream qDotFile("qDotVsTime");
+
+    // Integrate the motion of the spring for 4s using a symplectic method
+    scalar deltaT = 0.002;
+    for (scalar t=0; t<4; t+=deltaT)
+    {
+        qDot += 0.5*deltaT*qDdot;
+        q += deltaT*qDot;
+
+        // Update the body-state prior to the evaluation of the restraints
+        spring.forwardDynamicsCorrection
+        (
+            q,
+            w,
+            qDot,
+            qDdot
+        );
+
+        // Accumulate the restraint forces
+        fx = Zero;
+        spring.applyRestraints(fx);
+
+        // Calculate the body acceleration for the given state
+        // and restraint forces
+        spring.forwardDynamics
+        (
+            q,
+            w,
+            qDot,
+            tau,
+            fx,
+            qDdot
+        );
+
+        // Update the velocity
+        qDot += 0.5*deltaT*qDdot;
+
+        // Write the results for graph generation
+        // using 'gnuplot spring.gnuplot'
+        qFile << t << " " << q[0] << endl;
+        qDotFile << t << " " << qDot[0] << endl;
+    }
+
+    Info<< "\nEnd\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/test/rigidBodyDynamics/spring/spring.gnuplot

@@ -0,0 +1,76 @@
+#------------------------------------------------------------------------------
+# =========                 |
+# \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+#  \\    /   O peration     |
+#   \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
+#    \\/     M anipulation  |
+#-------------------------------------------------------------------------------
+# License
+#     This file is part of OpenFOAM.
+#
+#     OpenFOAM is free software: you can redistribute it and/or modify it
+#     under the terms of the GNU General Public License as published by
+#     the Free Software Foundation, either version 3 of the License, or
+#     (at your option) any later version.
+#
+#     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+#     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+#     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+#     for more details.
+#
+#     You should have received a copy of the GNU General Public License
+#     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+#
+# Script
+#     spring.gnuplot
+#
+# Description
+#     Creates an PostScript graph file of Test-spring results vs
+#     the analytical solution.
+#
+#------------------------------------------------------------------------------
+
+reset
+
+set samples 2000
+
+k = 5000.0
+m = 9.6
+c = 50.0
+a = -0.1
+
+omega = sqrt(k/m)
+zeta = c/(2.0*m*omega)
+
+phi = atan((sqrt(1.0 - zeta**2))/zeta)
+A = a/sin(phi)
+
+pos(A, t, omega, phi, zeta) = A*exp(-zeta*omega*t)*sin(sqrt(1-zeta**2)*omega*t + phi)
+vel(A, t, omega, phi, zeta) = \
+A*exp(-zeta*omega*t)*\
+( \
+  sqrt(1-zeta**2)*omega*cos(sqrt(1-zeta**2)*omega*t + phi) \
+- zeta*omega*sin(sqrt(1-zeta**2)*omega*t + phi) \
+)
+
+set xlabel "Time/[s]"
+set ylabel "Position"
+
+set ytics nomirror
+set y2tics
+
+set yrange [-0.1:0.1]
+set y2range [-2:2]
+
+set xzeroaxis
+
+set terminal postscript eps color enhanced solid
+set output "spring.eps"
+
+plot \
+    "qVsTime" u 1:($2 - 0.1) w l t "Simulation, centre of mass relative to start", \
+    pos(A, x, omega, phi, zeta) w l t "Analytical solution, centre of mass", \
+    "qDotVsTime" u 1:2 w l axes x1y2 t "Simulation, vertical velocity", \
+    vel(A, x, omega, phi, zeta) w l axes x1y2 t "Analytical solution, vertical velocity"
+
+#------------------------------------------------------------------------------