zhyang-dev/OpenFOAM-6
1
0
Fork 0
mirror of https://github.com/OpenFOAM/OpenFOAM-6.git synced 2025-12-08 06:57:46 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

rigidBodyDynamics: Ported tests to general RBD application

Browse Source
This commit is contained in:
Will Bainbridge
2018-05-08 09:08:46 +01:00
parent d501e87ec2
commit 59be3e7113
24 changed files with 307 additions and 618 deletions
Download Patch File Download Diff File Expand all files Collapse all files

0
applications/test/rigidBodyDynamics/rigidBodyDynamics/Make/files → applications/test/rigidBodyDynamics/Make/files
View File

0
applications/test/rigidBodyDynamics/pendulum/Make/options → applications/test/rigidBodyDynamics/Make/options
View File

52
applications/test/rigidBodyDynamics/pendulum Normal file
View File

@ -0,0 +1,52 @@
solver
{
type symplectic;
}
bodies
{
hinge
{
type masslessBody;
parent root;
transform (1 0 0 0 1 0 0 0 1) (0 0 0);
joint
{
type Rz;
}
outline
(
(0 0 0)
(0 -0.95 0)
(0.025 -0.95669873 0)
(0.04330127 -0.975 0)
(0.05 -1 0)
(0.04330127 -1.025 0)
(0.025 -1.04330127 0)
(0 -1.05 0)
(-0.025 -1.04330127 0)
(-0.04330127 -1.025 0)
(-0.05 -1 0)
(-0.04330127 -0.975 0)
(-0.025 -0.95669873 0)
(0 -0.95 0)
);
}
weight
{
type sphere;
mass 1;
radius 0.05;
centreOfMass (0 0 0);
transform (1 0 0 0 1 0 0 0 1) (0 -1 0);
mergeWith hinge;
}
}
g (0 -9.81 0);
q (0.3);
deltaT 0.01;
endTime 4.1;

3
applications/test/rigidBodyDynamics/pendulum/Make/files
View File

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

38
applications/test/rigidBodyDynamics/pendulum/pendulum
View File

@ -1,38 +0,0 @@
bodies0
{
pendulum
{
type rigidBody;
mass 1;
centreOfMass (0 -1 0);
inertia (0.001 0 0 0.001 0 0.001);
parent root;
transform (1 0 0 0 1 0 0 0 1) (0 0 0);
joint
{
type Rz;
}
}
}
bodies
{
hinge
{
type masslessBody;
parent root;
transform (1 0 0 0 1 0 0 0 1) (0 0 0);
joint
{
type Rz;
}
}
weight
{
type sphere;
mass 1;
radius 0.05;
transform (1 0 0 0 1 0 0 0 1) (0 -1 0);
mergeWith hinge;
}
}

124
applications/test/rigidBodyDynamics/pendulum/pendulum.C
View File

@ -1,124 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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
pendulum
Description
Simple swinging pendulum simulation with 1-DoF. The motion is integrated
using a symplectic method for just over 2-periods.
\*---------------------------------------------------------------------------*/
#include "rigidBodyModel.H"
#include "masslessBody.H"
#include "rigidBodyModelState.H"
#include "sphere.H"
#include "joints.H"
#include "IFstream.H"
using namespace Foam;
using namespace RBD;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
/*
bool testMerge = true;
// Create a model for the pendulum
rigidBodyModel pendulum;
// Join a weight to the origin with a centre of mass -1m below the origin
// by a hinge which rotates about the z-axis
if (testMerge)
{
pendulum.join
(
pendulum.bodyID("root"),
Xt(Zero),
joint::New(new joints::Rz()),
autoPtr<rigidBody>(new masslessBody("hinge"))
);
pendulum.merge
(
pendulum.bodyID("hinge"),
Xt(vector(0, -1, 0)),
autoPtr<rigidBody>(new sphere("weight", 1, 0.05))
);
}
else
{
pendulum.join
(
pendulum.bodyID("root"),
Xt(Zero),
joint::New(new joints::Rz()),
rigidBody::New("pendulum", 1, vector(0, -1, 0), 1e-3*I)
);
}
*/
// Create the pendulum model from dictionary
rigidBodyModel pendulum(dictionary(IFstream("pendulum")()));
Info<< pendulum << endl;
// Create the joint-space state fields
rigidBodyModelState pendulumState(pendulum);
scalarField& q = pendulumState.q();
scalarField& qDot = pendulumState.qDot();
scalarField& qDdot = pendulumState.qDdot();
scalarField tau(pendulum.nDoF(), Zero);
// Set the angle of the pendulum to 0.3rad
q[0] = 0.3;
// Set the gravitational acceleration
pendulum.g() = vector(0, -9.81, 0);
// Integrate the motion of the pendulum for 4.1s (~2-periods) using a
// symplectic method
scalar deltaT = 0.01;
for (scalar t=0; t<4.1; t+=deltaT)
{
qDot += 0.5*deltaT*qDdot;
q += deltaT*qDot;
pendulum.forwardDynamics(pendulumState, tau, Field<spatialVector>());
qDot += 0.5*deltaT*qDdot;
Info<< "Time " << t << "s, angle = " << q[0] << "rad" << endl;
}
Info<< "\nEnd\n" << endl;
return 0;
}
// ************************************************************************* //

42
applications/test/rigidBodyDynamics/pendulumAndSpring/pendulumAndSpring → applications/test/rigidBodyDynamics/pendulumAndSpring
View File

@ -3,9 +3,6 @@ solver
type Newmark;
}
// It is necessary to iterate for the Newmark solver
nIter 2;
bodies
{
M
@ -14,6 +11,7 @@ bodies
parent root;
mass 6e4;
radius 0.01;
centreOfMass (0 0 0);
transform (1 0 0 0 1 0 0 0 1) (0 0 0);
joint
{
@ -29,6 +27,21 @@ bodies
}
);
}
outline
(
(0 0.01 0)
(0.005 0.00866025 0)
(0.00866025 0.005 0)
(0.01 0 0)
(0.00866025 -0.005 0)
(0.005 -0.00866025 0)
(0 -0.01 0)
(-0.00866025 -0.005 0)
(-0.01 0 0)
(-0.00866025 0.005 0)
(-0.005 0.00866025 0)
(0 0.01 0)
);
}
m
@ -37,8 +50,24 @@ bodies
parent M;
mass 6e3;
radius 0.01;
centreOfMass (0 0 0);
transform (1 0 0 0 1 0 0 0 1) (0 -5 0);
mergeWith M;
outline
(
(0 0.01 0)
(0.005 0.00866025 0)
(0.00866025 0.005 0)
(0.01 0 0)
(0.00866025 -0.005 0)
(0.005 -0.00866025 0)
(0 -0.01 0)
(-0.00866025 -0.005 0)
(-0.01 0 0)
(-0.00866025 0.005 0)
(-0.005 0.00866025 0)
(0 0.01 0)
);
}
}
@ -62,3 +91,10 @@ g (0 -9.81 0);
// Set the initial offset of the pendulum to 1.5m
// and the angle of the pendulum to 30deg
q (1.5 0.5235987);
deltaT 0.01;
// It is necessary to iterate for the Newmark solver
nIter 2;
endTime 20;

6
applications/test/rigidBodyDynamics/pendulumAndSpring/pendulumAndSpring.gnuplot → applications/test/rigidBodyDynamics/pendulumAndSpring.gnuplot
View File

@ -2,7 +2,7 @@
# ========= |
# \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
# \\ / O peration |
# \\ / A nd | Copyright (C) 2016 OpenFOAM Foundation
# \\ / A nd | Copyright (C) 2016-2018 OpenFOAM Foundation
# \\/ M anipulation |
#-------------------------------------------------------------------------------
# License
@ -47,7 +47,7 @@ set terminal postscript eps color enhanced solid
set output "pendulumAndSpring.eps"
plot \
"xVsTime" u 1:2 w l t "x", \
"omegaVsTime" u 1:(57.29578*$2) w l axes x1y2 t "omega"
"pendulumAndSpring.dat" u 1:2 w l t "x", \
"pendulumAndSpring.dat" u 1:(57.29578*$2) w l axes x1y2 t "omega"
#------------------------------------------------------------------------------

3
applications/test/rigidBodyDynamics/pendulumAndSpring/Make/files
View File

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

5
applications/test/rigidBodyDynamics/pendulumAndSpring/Make/options
View File

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

92
applications/test/rigidBodyDynamics/pendulumAndSpring/pendulumAndSpring.C
View File

@ -1,92 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2016-2017 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
pendulumAndSpring
Description
Simple 2-DoF pendulum and spring simulation.
\*---------------------------------------------------------------------------*/
#include "rigidBodyMotion.H"
#include "masslessBody.H"
#include "sphere.H"
#include "joints.H"
#include "rigidBodyRestraint.H"
#include "rigidBodyModelState.H"
#include "IFstream.H"
#include "OFstream.H"
#include "constants.H"
using namespace Foam;
using namespace RBD;
using namespace Foam::constant::mathematical;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
dictionary pendulumAndSpringDict(IFstream("pendulumAndSpring")());
// Create the pendulumAndSpring model from dictionary
rigidBodyMotion pendulumAndSpring(pendulumAndSpringDict);
label nIter(readLabel(pendulumAndSpringDict.lookup("nIter")));
Info<< pendulumAndSpring << endl;
Info<< "// Joint state " << endl;
pendulumAndSpring.state().write(Info);
// Create the joint-space force field
scalarField tau(pendulumAndSpring.nDoF(), Zero);
// Create the external body force field
Field<spatialVector> fx(pendulumAndSpring.nBodies(), Zero);
OFstream xFile("xVsTime");
OFstream omegaFile("omegaVsTime");
// Integrate the motion of the pendulumAndSpring for 100s
scalar deltaT = 0.01;
for (scalar t=0; t<20; t+=deltaT)
{
pendulumAndSpring.newTime();
for (label i=0; i<nIter; i++)
{
pendulumAndSpring.solve(t + deltaT, deltaT, tau, fx);
}
// Write the results for graph generation
xFile << t << " " << pendulumAndSpring.state().q()[0] << endl;
omegaFile << t << " " << pendulumAndSpring.state().q()[1] << endl;
}
Info<< "\nEnd\n" << endl;
return 0;
}
// ************************************************************************* //

186
applications/test/rigidBodyDynamics/rigidBodyDynamics.C Normal file
View File

@ -0,0 +1,186 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2018 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
RBD
Description
Does an RBD simulation and outputs the result as a gnuplot animation
\*---------------------------------------------------------------------------*/
#include "rigidBodyMotion.H"
#include "IFstream.H"
#include "OFstream.H"
#include "boundBox.H"
#include "argList.H"
using namespace Foam;
using namespace RBD;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
// Create the input dictionary
argList::validArgs.append("dictionary");
argList args(argc, argv);
const word dictName(args[1]);
Info<< "Reading " << dictName << nl << endl;
const dictionary dict = IFstream(dictName)();
// Read the model, time controls and plot outlines from the dictionary
rigidBodyMotion motion(dict);
const scalar deltaT(readScalar(dict.lookup("deltaT")));
const label nIter(dict.lookupOrDefault<label>("nIter", 1));
const scalar endTime(readScalar(dict.lookup("endTime")));
const dictionary& bodiesDict = dict.subDict("bodies");
List<vectorField> outlines;
labelList outlineBodyIDs;
forAll(motion.bodies(), bodyID)
{
const word& bodyName = motion.bodies()[bodyID].name();
if (bodiesDict.isDict(bodyName))
{
const dictionary& bodyDict = bodiesDict.subDict(bodyName);
if (bodyDict.found("outline"))
{
outlines.append(vectorField(bodyDict.lookup("outline")));
outlineBodyIDs.append(bodyID);
}
}
}
// Set up motion fields
const label nDoF = motion.nDoF();
Info << nDoF << " degrees of freedom" << endl;
scalarField tau(nDoF, Zero);
Field<spatialVector> fx(motion.nBodies(), Zero);
// Initialise output files and the bound box
OFstream dataFile(dictName + ".dat");
OFstream animationFile(dictName + ".animate");
boundBox animationBox;
dataFile<< "# time(s)";
forAll(motion.state().q(), stateI)
{
dataFile << " q_" << stateI;
}
forAll(motion.state().q(), stateI)
{
dataFile << " qDot_" << stateI;
}
forAll(motion.state().q(), stateI)
{
dataFile << " qDdot_" << stateI;
}
dataFile<< endl;
animationFile
<< "#!/usr/bin/env gnuplot" << endl << endl
<< "$data << end" << endl;
// Run the RBD simulation
for (scalar t = 0; t <= endTime + 0.5*deltaT; t += deltaT)
{
Info().stdStream() << "\33[2KTime = " << t << '\r' << std::flush;
motion.newTime();
for (label i = 0; i < nIter; ++ i)
{
motion.solve(t + deltaT, deltaT, tau, fx);
}
// Write the current state
dataFile<< t;
forAll(motion.state().q(), stateI)
{
dataFile << ' ' << motion.state().q()[stateI];
}
forAll(motion.state().q(), stateI)
{
dataFile << ' ' << motion.state().qDot()[stateI];
}
forAll(motion.state().q(), stateI)
{
dataFile << ' ' << motion.state().qDdot()[stateI];
}
dataFile<< endl;
// Write the bodies' outlines at the current transformation
forAll(outlines, outlineI)
{
const label bodyID = outlineBodyIDs[outlineI];
const spatialTransform& invX0 = motion.X0(bodyID).inv();
forAll(outlines[outlineI], i)
{
const vector p = invX0.transformPoint(outlines[outlineI][i]);
animationFile<< t << ' ' << p.x() << ' ' << p.y() << endl;
animationBox.min() = min(animationBox.min(), p);
animationBox.max() = max(animationBox.max(), p);
}
animationFile<< endl;
}
}
Info << endl;
// Expand the bound box
{
const vector c = animationBox.midpoint(), d = animationBox.span();
if (d.x() < d.y()/2)
{
animationBox.min().x() = c.x() - d.y()/4;
animationBox.max().x() = c.x() + d.y()/4;
}
else if (d.y() < d.x()/2)
{
animationBox.min().y() = c.y() - d.x()/4;
animationBox.max().y() = c.y() + d.x()/4;
}
}
// Write the animation commands
animationFile
<< "end" << endl << endl
<< "set size ratio -1" << endl
<< "do for [i=1:" << label(endTime/deltaT - 0.5) << "] {" << endl
<< " set title sprintf('\%g s', i*" << deltaT << ')' << endl
<< " plot ["
<< animationBox.min().x() << ':' << animationBox.max().x() << "]["
<< animationBox.min().y() << ':' << animationBox.max().y() << ']';
forAll(outlines, outlineI)
{
const label bodyID = outlineBodyIDs[outlineI];
animationFile<< (outlineI ? "," : "") << " \\" << endl
<< " '$data' us 2:3 every :::"
<< outlines.size() << "*i+" << outlineI << "::"
<< outlines.size() << "*i+" << outlineI << " w l lw 2 "
<< "t '" << motion.bodies()[bodyID].name() << '\'';
}
animationFile<< endl << " pause " << deltaT << endl << "}" << endl;
return 0;
}
// ************************************************************************* //

5
applications/test/rigidBodyDynamics/rigidBodyDynamics/Make/options
View File

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

123
applications/test/rigidBodyDynamics/rigidBodyDynamics/rigidBodyDynamics.C
View File

@ -1,123 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2018 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
RBD
Description
Does an RBD simulation and outputs the result as a gnuplot animation
\*---------------------------------------------------------------------------*/
#include "rigidBodyMotion.H"
#include "IFstream.H"
#include "OFstream.H"
#include "boundBox.H"
#include "argList.H"
using namespace Foam;
using namespace RBD;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
// Create the input dictionary
argList::validArgs.append("dictionary");
argList args(argc, argv);
const word dictName(args[1]);
Info<< "Reading " << dictName << nl << endl;
const dictionary dict = IFstream(dictName)();
// Read the model, time controls and plot outlines from the dictionary
rigidBodyMotion motion(dict);
const scalar deltaT(readScalar(dict.lookup("deltaT")));
const scalar endTime(readScalar(dict.lookup("endTime")));
const dictionary& bodiesDict = dict.subDict("bodies");
List<vectorField> bodiesOutline(bodiesDict.size());
const label i0 = motion.bodies().size() - bodiesOutline.size();
forAll(bodiesOutline, i)
{
const dictionary& bodyDict =
bodiesDict.subDict(motion.bodies()[i0 + i].name());
bodiesOutline[i] = vectorField(bodyDict.lookup("outline"));
}
// Set up motion fields
scalarField tau(motion.nDoF(), Zero);
Field<spatialVector> fx(motion.nBodies(), Zero);
Info << motion.nDoF() << " degrees of freedom" << endl;
// Initialise plot bound box and file
boundBox plotBox;
OFstream plotFile(dictName + ".gnuplot");
plotFile<< "$data << end" << endl;
// Run the RBD simulation
for (scalar t = 0; t <= endTime + 0.5*deltaT; t += deltaT)
{
Info().stdStream() << "\33[2KTime = " << t << '\r' << std::flush;
motion.newTime();
motion.solve(t + deltaT, deltaT, tau, fx);
// Write the bodies' outlines at the current transformation
forAll(bodiesOutline, i)
{
const spatialTransform& invX0 = motion.X0(i0 + i).inv();
forAll(bodiesOutline[i], j)
{
const vector p = invX0.transformPoint(bodiesOutline[i][j]);
plotFile<< t << ' ' << p.x() << ' ' << p.y() << endl;
plotBox.min() = min(plotBox.min(), p);
plotBox.max() = max(plotBox.max(), p);
}
plotFile<< endl;
}
}
Info << endl;
// Write the plot commands
plotFile
<< "end" << endl << endl
<< "set size ratio -1" << endl
<< "do for [i=1:" << label(endTime/deltaT - 0.5) << "] {" << endl
<< " set title sprintf('\%g s', i*" << deltaT << ')' << endl
<< " plot "
<< "[" << plotBox.min().x() << ':' << plotBox.max().x() << ']'
<< "[" << plotBox.min().y() << ':' << plotBox.max().y() << ']';
forAll(bodiesOutline, i)
{
plotFile<< (i ? "," : "") << " \\" << endl
<< " '$data' us 2:3 every :::"
<< bodiesOutline.size() << "*i+" << i << "::"
<< bodiesOutline.size() << "*i+" << i << " w lp lw 2 pt 7 "
<< "t '" << motion.bodies()[i0 + i].name() << '\'';
}
plotFile<< endl << " pause " << deltaT << endl << "}" << endl;
return 0;
}
// ************************************************************************* //

15
applications/test/rigidBodyDynamics/sphericalJoint/sphericalJoint → applications/test/rigidBodyDynamics/sphericalJoint
View File

@ -3,9 +3,6 @@ solver
type Newmark;
}
// It is necessary to iterate for the Newmark solver
nIter 2;
bodies
{
pendulum
@ -20,5 +17,17 @@ bodies
{
type Rs;
}
outline ((0 0 0) (0 -2 0));
}
}
g (0 -9.81 0);
q (0 0 0.149438);
deltaT 0.01;
// It is necessary to iterate for the Newmark solver
nIter 2;
endTime 4.1;

3
applications/test/rigidBodyDynamics/sphericalJoint/Make/files
View File

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

5
applications/test/rigidBodyDynamics/sphericalJoint/Make/options
View File

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

103
applications/test/rigidBodyDynamics/sphericalJoint/sphericalJoint.C
View File

@ -1,103 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2016-2017 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
sphericalJoint
Description
Simple spherical-joint pendulum.
\*---------------------------------------------------------------------------*/
#include "rigidBodyMotion.H"
#include "masslessBody.H"
#include "sphere.H"
#include "joints.H"
#include "rigidBodyRestraint.H"
#include "rigidBodyModelState.H"
#include "IFstream.H"
#include "OFstream.H"
using namespace Foam;
using namespace RBD;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
dictionary sphericalJointDict(IFstream("sphericalJoint")());
// Create the sphericalJoint model from dictionary
rigidBodyMotion sphericalJoint(sphericalJointDict);
label nIter(readLabel(sphericalJointDict.lookup("nIter")));
Info<< sphericalJoint << endl;
// Create the joint-space force field
scalarField tau(sphericalJoint.nDoF(), Zero);
// Create the external body force field
Field<spatialVector> fx(sphericalJoint.nBodies(), Zero);
// Set the angle of the pendulum to 0.3rad
sphericalJoint.joints()[1].unitQuaternion
(
quaternion(quaternion::ZYX, vector(0.3, 0, 0)),
sphericalJoint.state().q()
);
// Set the gravitational acceleration
sphericalJoint.g() = vector(0, -9.81, 0);
OFstream omegaFile("omegaVsTime");
// Integrate the motion of the sphericalJoint for 4.1s
scalar deltaT = 0.01;
for (scalar t=0; t<4.1; t+=deltaT)
{
sphericalJoint.newTime();
for (label i=0; i<nIter; i++)
{
sphericalJoint.solve(t + deltaT, deltaT, tau, fx);
}
// Write the results for graph generation
// using 'gnuplot sphericalJoint.gnuplot'
omegaFile
<< t << " "
<< sphericalJoint.joints()[1].unitQuaternion
(
sphericalJoint.state().q()
).eulerAngles(quaternion::ZYX).x()
<< endl;
}
Info<< "\nEnd\n" << endl;
return 0;
}
// ************************************************************************* //

12
applications/test/rigidBodyDynamics/spring/spring → applications/test/rigidBodyDynamics/spring
View File

@ -3,9 +3,6 @@ solver
type Newmark;
}
// It is necessary to iterate for the Newmark solver
nIter 2;
bodies
{
weight
@ -36,3 +33,12 @@ restraints
restLength 0.4;
}
}
g (0 0 0);
deltaT 0.002;
// It is necessary to iterate for the Newmark solver
nIter 2;
endTime 4;

11
applications/test/rigidBodyDynamics/spring/spring.gnuplot → applications/test/rigidBodyDynamics/spring.gnuplot
View File

@ -2,7 +2,7 @@
# ========= |
# \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
# \\ / O peration |
# \\ / A nd | Copyright (C) 2016 OpenFOAM Foundation
# \\ / A nd | Copyright (C) 2016-2018 OpenFOAM Foundation
# \\/ M anipulation |
#-------------------------------------------------------------------------------
# License
@ -53,8 +53,9 @@ A*exp(-zeta*omega*t)*\
- zeta*omega*sin(sqrt(1-zeta**2)*omega*t + phi) \
)
set xlabel "Time/[s]"
set ylabel "Position"
set xlabel "Time (s)"
set ylabel "Position (m)"
set y2label "Velocity (m/s)"
set ytics nomirror
set y2tics
@ -68,9 +69,9 @@ 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", \
"spring.dat" 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", \
"spring.dat" u 1:3 w l axes x1y2 t "Simulation, vertical velocity", \
vel(A, x, omega, phi, zeta) w l axes x1y2 t "Analytical solution, vertical velocity"
#------------------------------------------------------------------------------

3
applications/test/rigidBodyDynamics/spring/Make/files
View File

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

5
applications/test/rigidBodyDynamics/spring/Make/options
View File

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

89
applications/test/rigidBodyDynamics/spring/spring.C
View File

@ -1,89 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2016-2017 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-spring simulation with 1-DoF.
\*---------------------------------------------------------------------------*/
#include "rigidBodyMotion.H"
#include "masslessBody.H"
#include "sphere.H"
#include "joints.H"
#include "rigidBodyRestraint.H"
#include "rigidBodyModelState.H"
#include "IFstream.H"
#include "OFstream.H"
using namespace Foam;
using namespace RBD;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
dictionary springDict(IFstream("spring")());
// Create the spring model from dictionary
rigidBodyMotion spring(springDict);
label nIter(readLabel(springDict.lookup("nIter")));
Info<< spring << endl;
// Create the joint-space force field
scalarField tau(spring.nDoF(), Zero);
// Create the external body force field
Field<spatialVector> fx(spring.nBodies(), Zero);
OFstream qFile("qVsTime");
OFstream qDotFile("qDotVsTime");
// Integrate the motion of the spring for 4s
scalar deltaT = 0.002;
for (scalar t=0; t<4; t+=deltaT)
{
spring.newTime();
for (label i=0; i<nIter; i++)
{
spring.solve(t + deltaT, deltaT, tau, fx);
}
// Write the results for graph generation
// using 'gnuplot spring.gnuplot'
qFile << t << " " << spring.state().q()[0] << endl;
qDotFile << t << " " << spring.state().qDot()[0] << endl;
}
Info<< "\nEnd\n" << endl;
return 0;
}
// ************************************************************************* //

0
applications/test/rigidBodyDynamics/rigidBodyDynamics/triplePendulum → applications/test/rigidBodyDynamics/triplePendulum
View File