diff --git a/src/OpenFOAM/meshes/primitiveShapes/plane/plane.C b/src/OpenFOAM/meshes/primitiveShapes/plane/plane.C index 2e90e4d62c..6f711cc603 100644 --- a/src/OpenFOAM/meshes/primitiveShapes/plane/plane.C +++ b/src/OpenFOAM/meshes/primitiveShapes/plane/plane.C @@ -128,7 +128,20 @@ Foam::plane::plane(const vector& normalVector) : unitVector_(normalVector), basePoint_(vector::zero) -{} +{ + scalar magUnitVector(mag(unitVector_)); + + if (magUnitVector > VSMALL) + { + unitVector_ /= magUnitVector; + } + else + { + FatalErrorIn("plane::plane(const point&, const vector&)") + << "plane normal has zero length" + << abort(FatalError); + } +} // Construct from point and normal vector @@ -146,8 +159,8 @@ Foam::plane::plane(const point& basePoint, const vector& normalVector) else { FatalErrorIn("plane::plane(const point&, const vector&)") - << "plane normal has got zero length" - << abort(FatalError); + << "plane normal has zero length" + << abort(FatalError); } } @@ -217,8 +230,8 @@ Foam::plane::plane(const dictionary& dict) "plane::plane(const dictionary&)", dict ) - << "Invalid plane type: " << planeType - << abort(FatalIOError); + << "Invalid plane type: " << planeType + << abort(FatalIOError); } } @@ -238,7 +251,7 @@ Foam::plane::plane(Istream& is) else { FatalErrorIn("plane::plane(Istream& is)") - << "plane normal has got zero length" + << "plane normal has zero length" << abort(FatalError); } } diff --git a/src/meshTools/momentOfInertia/momentOfInertia.C b/src/meshTools/momentOfInertia/momentOfInertia.C new file mode 100644 index 0000000000..5207850172 --- /dev/null +++ b/src/meshTools/momentOfInertia/momentOfInertia.C @@ -0,0 +1,403 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2007-2009 OpenCFD Ltd. + \\/ 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 2 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, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +Class + momentOfInertia + +Description + Reimplementation of volInt.c by Brian Mirtich. + * mirtich@cs.berkeley.edu * + * http://www.cs.berkeley.edu/~mirtich * + +------------------------------------------------------------------------------- +*/ + +#include "momentOfInertia.H" +//#include "pyramidPointFaceRef.H" + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +//Foam::tensor Foam::momentOfInertia +//( +// const pointField& points, +// const faceList& faces, +// const cell& cFaces, +// const point& cc +//) +//{ +// tensor t(tensor::zero); +// +// forAll(cFaces, i) +// { +// const face& f = faces[cFaces[i]]; +// +// scalar pyrVol = pyramidPointFaceRef(f, cc).mag(points); +// +// vector pyrCentre = pyramidPointFaceRef(f, cc).centre(points); +// +// vector d = pyrCentre - cc; +// +// t.xx() += pyrVol*(sqr(d.y()) + sqr(d.z())); +// t.yy() += pyrVol*(sqr(d.x()) + sqr(d.z())); +// t.zz() += pyrVol*(sqr(d.x()) + sqr(d.y())); +// +// t.xy() -= pyrVol*d.x()*d.y(); +// t.xz() -= pyrVol*d.x()*d.z(); +// t.yz() -= pyrVol*d.y()*d.z(); +// } +// +// // Symmetric +// t.yx() = t.xy(); +// t.zx() = t.xz(); +// t.zy() = t.yz(); +// +// return t; +//} + + +#define sqr(x) ((x)*(x)) +#define pow3(x) ((x)*(x)*(x)) + +// compute various integrations over projection of face +void Foam::compProjectionIntegrals +( + const pointField& points, + const face& f, + const direction A, + const direction B, + + scalar& P1, + scalar& Pa, + scalar& Pb, + scalar& Paa, + scalar& Pab, + scalar& Pbb, + scalar& Paaa, + scalar& Paab, + scalar& Pabb, + scalar& Pbbb +) +{ + P1 = Pa = Pb = Paa = Pab = Pbb = Paaa = Paab = Pabb = Pbbb = 0.0; + + forAll(f, i) + { + scalar a0 = points[f[i]][A]; + scalar b0 = points[f[i]][B]; + scalar a1 = points[f[(i+1) % f.size()]][A]; + scalar b1 = points[f[(i+1) % f.size()]][B]; + scalar da = a1 - a0; + scalar db = b1 - b0; + + scalar a0_2 = a0 * a0; + scalar a0_3 = a0_2 * a0; + scalar a0_4 = a0_3 * a0; + + scalar b0_2 = b0 * b0; + scalar b0_3 = b0_2 * b0; + scalar b0_4 = b0_3 * b0; + + scalar a1_2 = a1 * a1; + scalar a1_3 = a1_2 * a1; + + scalar b1_2 = b1 * b1; + scalar b1_3 = b1_2 * b1; + + scalar C1 = a1 + a0; + + scalar Ca = a1*C1 + a0_2; + scalar Caa = a1*Ca + a0_3; + scalar Caaa = a1*Caa + a0_4; + + scalar Cb = b1*(b1 + b0) + b0_2; + scalar Cbb = b1*Cb + b0_3; + scalar Cbbb = b1*Cbb + b0_4; + + scalar Cab = 3*a1_2 + 2*a1*a0 + a0_2; + scalar Kab = a1_2 + 2*a1*a0 + 3*a0_2; + + scalar Caab = a0*Cab + 4*a1_3; + scalar Kaab = a1*Kab + 4*a0_3; + + scalar Cabb = 4*b1_3 + 3*b1_2*b0 + 2*b1*b0_2 + b0_3; + scalar Kabb = b1_3 + 2*b1_2*b0 + 3*b1*b0_2 + 4*b0_3; + + P1 += db*C1; + Pa += db*Ca; + Paa += db*Caa; + Paaa += db*Caaa; + Pb += da*Cb; + Pbb += da*Cbb; + Pbbb += da*Cbbb; + Pab += db*(b1*Cab + b0*Kab); + Paab += db*(b1*Caab + b0*Kaab); + Pabb += da*(a1*Cabb + a0*Kabb); + } + + P1 /= 2.0; + Pa /= 6.0; + Paa /= 12.0; + Paaa /= 20.0; + Pb /= -6.0; + Pbb /= -12.0; + Pbbb /= -20.0; + Pab /= 24.0; + Paab /= 60.0; + Pabb /= -60.0; +} + + +void Foam::compFaceIntegrals +( + const pointField& points, + const face& f, + const vector& n, + const scalar w, + const direction A, + const direction B, + const direction C, + + scalar& Fa, + scalar& Fb, + scalar& Fc, + scalar& Faa, + scalar& Fbb, + scalar& Fcc, + scalar& Faaa, + scalar& Fbbb, + scalar& Fccc, + scalar& Faab, + scalar& Fbbc, + scalar& Fcca +) +{ + scalar P1, Pa, Pb, Paa, Pab, Pbb, Paaa, Paab, Pabb, Pbbb; + + compProjectionIntegrals + ( + points, + f, + A, + B, + + P1, + Pa, + Pb, + Paa, + Pab, + Pbb, + Paaa, + Paab, + Pabb, + Pbbb + ); + + scalar k1 = 1 / n[C]; + scalar k2 = k1 * k1; + scalar k3 = k2 * k1; + scalar k4 = k3 * k1; + + Fa = k1 * Pa; + Fb = k1 * Pb; + Fc = -k2 * (n[A]*Pa + n[B]*Pb + w*P1); + + Faa = k1 * Paa; + Fbb = k1 * Pbb; + Fcc = k3 * (sqr(n[A])*Paa + 2*n[A]*n[B]*Pab + sqr(n[B])*Pbb + + w*(2*(n[A]*Pa + n[B]*Pb) + w*P1)); + + Faaa = k1 * Paaa; + Fbbb = k1 * Pbbb; + Fccc = -k4 * (pow3(n[A])*Paaa + 3*sqr(n[A])*n[B]*Paab + + 3*n[A]*sqr(n[B])*Pabb + pow3(n[B])*Pbbb + + 3*w*(sqr(n[A])*Paa + 2*n[A]*n[B]*Pab + sqr(n[B])*Pbb) + + w*w*(3*(n[A]*Pa + n[B]*Pb) + w*P1)); + + Faab = k1 * Paab; + Fbbc = -k2 * (n[A]*Pabb + n[B]*Pbbb + w*Pbb); + Fcca = k3 * (sqr(n[A])*Paaa + 2*n[A]*n[B]*Paab + sqr(n[B])*Pabb + + w*(2*(n[A]*Paa + n[B]*Pab) + w*Pa)); +} + + +void Foam::compVolumeIntegrals +( + const pointField& points, + const faceList& faces, + const cell& cFaces, + const vectorField& fNorm, + const scalarField& fW, + + scalar& T0, + vector& T1, + vector& T2, + vector& TP +) +{ + T0 = 0; + T1 = vector::zero; + T2 = vector::zero; + TP = vector::zero; + + forAll(cFaces, i) + { + const vector& n = fNorm[i]; + + scalar nx = mag(n[0]); + scalar ny = mag(n[1]); + scalar nz = mag(n[2]); + + direction A, B, C; + + if (nx > ny && nx > nz) + { + C = 0; + } + else + { + C = (ny > nz) ? 1 : 2; + } + + A = (C + 1) % 3; + B = (A + 1) % 3; + + scalar Fa, Fb, Fc, Faa, Fbb, Fcc, Faaa, Fbbb, Fccc, Faab, Fbbc, Fcca; + compFaceIntegrals + ( + points, + faces[cFaces[i]], + n, + fW[i], + A, + B, + C, + + Fa, + Fb, + Fc, + Faa, + Fbb, + Fcc, + Faaa, + Fbbb, + Fccc, + Faab, + Fbbc, + Fcca + ); + + T0 += n[0] * ((A == 0) ? Fa : ((B == 0) ? Fb : Fc)); + + T1[A] += n[A] * Faa; + T1[B] += n[B] * Fbb; + T1[C] += n[C] * Fcc; + + T2[A] += n[A] * Faaa; + T2[B] += n[B] * Fbbb; + T2[C] += n[C] * Fccc; + + TP[A] += n[A] * Faab; + TP[B] += n[B] * Fbbc; + TP[C] += n[C] * Fcca; + } + + T1 /= 2; + T2 /= 3; + TP /= 2; +} + + +// Calculate +// - r: centre of mass +// - J: inertia around origin (point 0,0,0) +void Foam::momentOfIntertia +( + const pointField& points, + const faceList& faces, + const cell& cFaces, + point& r, + tensor& J +) +{ + // Face normals + vectorField fNorm(cFaces.size()); + scalarField fW(cFaces.size()); + + forAll(cFaces, i) + { + label faceI = cFaces[i]; + + const face& f = faces[faceI]; + + fNorm[i] = f.normal(points); + fNorm[i] /= mag(fNorm[i]) + VSMALL; + + fW[i] = - (fNorm[i] & points[f[0]]); + } + + + scalar T0; + vector T1, T2, TP; + + compVolumeIntegrals + ( + points, + faces, + cFaces, + fNorm, + fW, + + T0, + T1, + T2, + TP + ); + + const scalar density = 1.0; /* assume unit density */ + + scalar mass = density * T0; + + /* compute center of mass */ + r = T1 / T0; + + /* compute inertia tensor */ + J.xx() = density * (T2[1] + T2[2]); + J.yy() = density * (T2[2] + T2[0]); + J.zz() = density * (T2[0] + T2[1]); + J.xy() = J.yx() = - density * TP[0]; + J.yz() = J.zy() = - density * TP[1]; + J.zx() = J.xz() = - density * TP[2]; + + ///* translate inertia tensor to center of mass */ + //J[XX] -= mass * (r[1]*r[1] + r[2]*r[2]); + //J[YY] -= mass * (r[2]*r[2] + r[0]*r[0]); + //J[ZZ] -= mass * (r[0]*r[0] + r[1]*r[1]); + //J[XY] = J[YX] += mass * r[0] * r[1]; + //J[YZ] = J[ZY] += mass * r[1] * r[2]; + //J[ZX] = J[XZ] += mass * r[2] * r[0]; +} + + + +// ************************************************************************* // diff --git a/src/meshTools/momentOfInertia/momentOfInertia.H b/src/meshTools/momentOfInertia/momentOfInertia.H new file mode 100644 index 0000000000..66a8d357ad --- /dev/null +++ b/src/meshTools/momentOfInertia/momentOfInertia.H @@ -0,0 +1,75 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2007-2009 OpenCFD Ltd. + \\/ 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 2 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, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +Class + momentOfInertia + +Description + +SourceFiles + momentOfInertia.H + +\*---------------------------------------------------------------------------*/ + +#ifndef momentOfInertia_H +#define momentOfInertia_H + +#include "tensor.H" +#include "primitiveMesh.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +////- Moment of inertia around cell centre for single cell. +//tensor momentOfInertia +//( +// const pointField&, +// const faceList&, +// const cell&, +// const point& cc +//); + +// Calculate +// - centre of mass +// - inertia tensor around (0,0,0) +void momentOfIntertia +( + const pointField&, + const faceList&, + const cell&, + point& r, + tensor& Jorigin +); + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/meshTools/momentOfInertia/volInt.ps.gz b/src/meshTools/momentOfInertia/volInt.ps.gz new file mode 100644 index 0000000000..19fb295421 Binary files /dev/null and b/src/meshTools/momentOfInertia/volInt.ps.gz differ diff --git a/src/meshTools/momentOfInertia/volumeIntegration/README b/src/meshTools/momentOfInertia/volumeIntegration/README new file mode 100644 index 0000000000..dadca677b7 --- /dev/null +++ b/src/meshTools/momentOfInertia/volumeIntegration/README @@ -0,0 +1,146 @@ +1. OVERVIEW + + This code accompanies the paper: + + Brian Mirtich, "Fast and Accurate Computation of + Polyhedral Mass Properties," journal of graphics + tools, volume 1, number 2, 1996. + + It computes the ten volume integrals needed for + determining the center of mass, moments of + inertia, and products of inertia for a uniform + density polyhedron. From this information, a + body frame can be computed. + + To compile the program, use an ANSI compiler, and + type something like + + % cc volInt.c -O2 -lm -o volInt + + + Revision history + + 26 Jan 1996 Program creation. + + 3 Aug 1996 Corrected bug arising when polyhedron density + is not 1.0. Changes confined to function main(). + Thanks to Zoran Popovic for catching this one. + + + +2. POLYHEDRON GEOMETRY FILES + + The program reads a data file specified on the + command line. This data file describes the + geometry of a polyhedron, and has the following + format: + + N + + x_0 y_0 z_0 + x_1 y_1 z_1 + . + . + . + x_{N-1} y_{N-1} z_{N-1} + + M + + k1 v_{1,1} v_{1,2} ... v_{1,k1} + k2 v_{2,1} v_{2,2} ... v_{2,k2} + . + . + . + kM v_{M,1} v_{M,2} ... v_{M,kM} + + + where: + N number of vertices on polyhedron + x_i y_i z_i x, y, and z coordinates of ith vertex + M number of faces on polyhedron + ki number of vertices on ith face + v_{i,j} jth vertex on ith face + + In English: + + First the number of vertices are specified. Next + the vertices are defined by listing the 3 + coordinates of each one. Next the number of faces + are specified. Finally, the faces themselves are + specified. A face is specified by first giving + the number of vertices around the polygonal face, + followed by the (integer) indices of these + vertices. When specifying indices, note that + they must be given in counter-clockwise order + (when looking at the face from outside the + polyhedron), and the vertices are indexed from 0 + to N-1 for a polyhedron with N faces. + + White space can be inserted (or not) as desired. + Three example polyhedron geometry files are included: + + cube A cube, 20 units on a side, centered at + the origin and aligned with the coordinate axes. + + tetra A tetrahedron formed by taking the convex + hull of the origin, and the points (5,0,0), + (0,4,0), and (0,0,3). + + icosa An icosahedron with vertices lying on the unit + sphere, centered at the origin. + + + +3. RUNNING THE PROGRAM + + Simply type, + + % volInt + + The program will read in the geometry of the + polyhedron, and the print out the ten volume + integrals. + + The program also computes some of the mass + properties which may be inferred from the volume + integrals. A density of 1.0 is assumed, although + this may be changed in function main(). The + center of mass is computed as well as the inertia + tensor relative to a frame with origin at the + center of mass. Note, however, that this matrix + may not be diagonal. If not, a diagonalization + routine must be performed, to determine the + orientation of the true body frame relative to + the original model frame. The Jacobi method + works quite well (see Numerical Recipes in C by + Press, et. al.). + + + +4. DISCLAIMERS + + 1. The volume integration code has been written + to match the development and algorithms presented + in the paper, and not with maximum optimization + in mind. While inherently very efficient, a few + more cycles can be squeezed out of the algorithm. + This is left as an exercise. :) + + 2. Don't like global variables? The three + procedures which evaluate the volume integrals + can be combined into a single procedure with two + nested loops. In addition to providing some + speedup, all of the global variables can then be + made local. + + 3. The polyhedron data structure used by the + program is admittedly lame; much better schemes + are possible. The idea here is just to give the + basic integral evaluation code, which will have + to be adjusted for other polyhedron data + structures. + + 4. There is no error checking for the input + files. Be careful. Note the hard limits + #defined for the number of vertices, number of + faces, and number of vertices per faces. diff --git a/src/meshTools/momentOfInertia/volumeIntegration/cube b/src/meshTools/momentOfInertia/volumeIntegration/cube new file mode 100644 index 0000000000..07ab26ad16 --- /dev/null +++ b/src/meshTools/momentOfInertia/volumeIntegration/cube @@ -0,0 +1,21 @@ +8 + +-10 -10 -10 ++10 -10 -10 ++10 +10 -10 +-10 +10 -10 +-10 -10 +10 ++10 -10 +10 ++10 +10 +10 +-10 +10 +10 + +6 + +4 0 3 2 1 +4 4 5 6 7 +4 0 1 5 4 +4 6 2 3 7 +4 1 2 6 5 +4 0 4 7 3 + + diff --git a/src/meshTools/momentOfInertia/volumeIntegration/icosa b/src/meshTools/momentOfInertia/volumeIntegration/icosa new file mode 100644 index 0000000000..16339ece89 --- /dev/null +++ b/src/meshTools/momentOfInertia/volumeIntegration/icosa @@ -0,0 +1,38 @@ +12 + ++0.00000000 +0.00000000 +1.00000000 ++0.00000000 +0.00000000 -1.00000000 ++0.89442719 +0.00000000 +0.44721360 ++0.27639320 +0.85065081 +0.44721360 +-0.72360680 +0.52573111 +0.44721360 +-0.72360680 -0.52573111 +0.44721360 ++0.27639320 -0.85065081 +0.44721360 ++0.72360680 +0.52573111 -0.44721360 +-0.27639320 +0.85065081 -0.44721360 +-0.89442719 +0.00000000 -0.44721360 +-0.27639320 -0.85065081 -0.44721360 ++0.72360680 -0.52573111 -0.44721360 + +20 + +3 6 11 2 +3 3 2 7 +3 7 2 11 +3 0 2 3 +3 0 6 2 +3 10 1 11 +3 1 7 11 +3 10 11 6 +3 8 7 1 +3 8 3 7 +3 5 10 6 +3 5 6 0 +3 4 3 8 +3 4 0 3 +3 9 8 1 +3 9 1 10 +3 4 5 0 +3 9 10 5 +3 9 5 4 +3 9 4 8 + diff --git a/src/meshTools/momentOfInertia/volumeIntegration/tetra b/src/meshTools/momentOfInertia/volumeIntegration/tetra new file mode 100644 index 0000000000..159e5a2bf8 --- /dev/null +++ b/src/meshTools/momentOfInertia/volumeIntegration/tetra @@ -0,0 +1,16 @@ +4 + +0 0 0 +5 0 0 +0 4 0 +0 0 3 + +4 + +3 0 3 2 +3 3 0 1 +3 2 1 0 +3 1 2 3 + + + diff --git a/src/meshTools/momentOfInertia/volumeIntegration/volInt.c b/src/meshTools/momentOfInertia/volumeIntegration/volInt.c new file mode 100644 index 0000000000..7601655488 --- /dev/null +++ b/src/meshTools/momentOfInertia/volumeIntegration/volInt.c @@ -0,0 +1,380 @@ + + + /******************************************************* + * * + * volInt.c * + * * + * This code computes volume integrals needed for * + * determining mass properties of polyhedral bodies. * + * * + * For more information, see the accompanying README * + * file, and the paper * + * * + * Brian Mirtich, "Fast and Accurate Computation of * + * Polyhedral Mass Properties," journal of graphics * + * tools, volume 1, number 1, 1996. * + * * + * This source code is public domain, and may be used * + * in any way, shape or form, free of charge. * + * * + * Copyright 1995 by Brian Mirtich * + * * + * mirtich@cs.berkeley.edu * + * http://www.cs.berkeley.edu/~mirtich * + * * + *******************************************************/ + +/* + Revision history + + 26 Jan 1996 Program creation. + + 3 Aug 1996 Corrected bug arising when polyhedron density + is not 1.0. Changes confined to function main(). + Thanks to Zoran Popovic for catching this one. + + 27 May 1997 Corrected sign error in translation of inertia + product terms to center of mass frame. Changes + confined to function main(). Thanks to + Chris Hecker. +*/ + + + +#include +#include +#include + +/* + ============================================================================ + constants + ============================================================================ +*/ + +#define MAX_VERTS 100 /* maximum number of polyhedral vertices */ +#define MAX_FACES 100 /* maximum number of polyhedral faces */ +#define MAX_POLYGON_SZ 10 /* maximum number of verts per polygonal face */ + +#define X 0 +#define Y 1 +#define Z 2 + +/* + ============================================================================ + macros + ============================================================================ +*/ + +#define SQR(x) ((x)*(x)) +#define CUBE(x) ((x)*(x)*(x)) + +/* + ============================================================================ + data structures + ============================================================================ +*/ + +typedef struct { + int numVerts; + double norm[3]; + double w; + int verts[MAX_POLYGON_SZ]; + struct polyhedron *poly; +} FACE; + +typedef struct polyhedron { + int numVerts, numFaces; + double verts[MAX_VERTS][3]; + FACE faces[MAX_FACES]; +} POLYHEDRON; + + +/* + ============================================================================ + globals + ============================================================================ +*/ + +static int A; /* alpha */ +static int B; /* beta */ +static int C; /* gamma */ + +/* projection integrals */ +static double P1, Pa, Pb, Paa, Pab, Pbb, Paaa, Paab, Pabb, Pbbb; + +/* face integrals */ +static double Fa, Fb, Fc, Faa, Fbb, Fcc, Faaa, Fbbb, Fccc, Faab, Fbbc, Fcca; + +/* volume integrals */ +static double T0, T1[3], T2[3], TP[3]; + + +/* + ============================================================================ + read in a polyhedron + ============================================================================ +*/ + +void readPolyhedron(char *name, POLYHEDRON *p) +{ + FILE *fp; + char line[200], *c; + int i, j, n; + double dx1, dy1, dz1, dx2, dy2, dz2, nx, ny, nz, len; + FACE *f; + + + if (!(fp = fopen(name, "r"))) { + printf("i/o error\n"); + exit(1); + } + + fscanf(fp, "%d", &p->numVerts); + printf("Reading in %d vertices\n", p->numVerts); + for (i = 0; i < p->numVerts; i++) + fscanf(fp, "%lf %lf %lf", + &p->verts[i][X], &p->verts[i][Y], &p->verts[i][Z]); + + fscanf(fp, "%d", &p->numFaces); + printf("Reading in %d faces\n", p->numFaces); + for (i = 0; i < p->numFaces; i++) { + f = &p->faces[i]; + f->poly = p; + fscanf(fp, "%d", &f->numVerts); + for (j = 0; j < f->numVerts; j++) fscanf(fp, "%d", &f->verts[j]); + + /* compute face normal and offset w from first 3 vertices */ + dx1 = p->verts[f->verts[1]][X] - p->verts[f->verts[0]][X]; + dy1 = p->verts[f->verts[1]][Y] - p->verts[f->verts[0]][Y]; + dz1 = p->verts[f->verts[1]][Z] - p->verts[f->verts[0]][Z]; + dx2 = p->verts[f->verts[2]][X] - p->verts[f->verts[1]][X]; + dy2 = p->verts[f->verts[2]][Y] - p->verts[f->verts[1]][Y]; + dz2 = p->verts[f->verts[2]][Z] - p->verts[f->verts[1]][Z]; + nx = dy1 * dz2 - dy2 * dz1; + ny = dz1 * dx2 - dz2 * dx1; + nz = dx1 * dy2 - dx2 * dy1; + len = sqrt(nx * nx + ny * ny + nz * nz); + f->norm[X] = nx / len; + f->norm[Y] = ny / len; + f->norm[Z] = nz / len; + f->w = - f->norm[X] * p->verts[f->verts[0]][X] + - f->norm[Y] * p->verts[f->verts[0]][Y] + - f->norm[Z] * p->verts[f->verts[0]][Z]; + + } + + fclose(fp); + +} + +/* + ============================================================================ + compute mass properties + ============================================================================ +*/ + + +/* compute various integrations over projection of face */ +void compProjectionIntegrals(FACE *f) +{ + double a0, a1, da; + double b0, b1, db; + double a0_2, a0_3, a0_4, b0_2, b0_3, b0_4; + double a1_2, a1_3, b1_2, b1_3; + double C1, Ca, Caa, Caaa, Cb, Cbb, Cbbb; + double Cab, Kab, Caab, Kaab, Cabb, Kabb; + int i; + + P1 = Pa = Pb = Paa = Pab = Pbb = Paaa = Paab = Pabb = Pbbb = 0.0; + + for (i = 0; i < f->numVerts; i++) { + a0 = f->poly->verts[f->verts[i]][A]; + b0 = f->poly->verts[f->verts[i]][B]; + a1 = f->poly->verts[f->verts[(i+1) % f->numVerts]][A]; + b1 = f->poly->verts[f->verts[(i+1) % f->numVerts]][B]; + da = a1 - a0; + db = b1 - b0; + a0_2 = a0 * a0; a0_3 = a0_2 * a0; a0_4 = a0_3 * a0; + b0_2 = b0 * b0; b0_3 = b0_2 * b0; b0_4 = b0_3 * b0; + a1_2 = a1 * a1; a1_3 = a1_2 * a1; + b1_2 = b1 * b1; b1_3 = b1_2 * b1; + + C1 = a1 + a0; + Ca = a1*C1 + a0_2; Caa = a1*Ca + a0_3; Caaa = a1*Caa + a0_4; + Cb = b1*(b1 + b0) + b0_2; Cbb = b1*Cb + b0_3; Cbbb = b1*Cbb + b0_4; + Cab = 3*a1_2 + 2*a1*a0 + a0_2; Kab = a1_2 + 2*a1*a0 + 3*a0_2; + Caab = a0*Cab + 4*a1_3; Kaab = a1*Kab + 4*a0_3; + Cabb = 4*b1_3 + 3*b1_2*b0 + 2*b1*b0_2 + b0_3; + Kabb = b1_3 + 2*b1_2*b0 + 3*b1*b0_2 + 4*b0_3; + + P1 += db*C1; + Pa += db*Ca; + Paa += db*Caa; + Paaa += db*Caaa; + Pb += da*Cb; + Pbb += da*Cbb; + Pbbb += da*Cbbb; + Pab += db*(b1*Cab + b0*Kab); + Paab += db*(b1*Caab + b0*Kaab); + Pabb += da*(a1*Cabb + a0*Kabb); + } + + P1 /= 2.0; + Pa /= 6.0; + Paa /= 12.0; + Paaa /= 20.0; + Pb /= -6.0; + Pbb /= -12.0; + Pbbb /= -20.0; + Pab /= 24.0; + Paab /= 60.0; + Pabb /= -60.0; +} + +compFaceIntegrals(FACE *f) +{ + double *n, w; + double k1, k2, k3, k4; + + compProjectionIntegrals(f); + + w = f->w; + n = f->norm; + k1 = 1 / n[C]; k2 = k1 * k1; k3 = k2 * k1; k4 = k3 * k1; + + Fa = k1 * Pa; + Fb = k1 * Pb; + Fc = -k2 * (n[A]*Pa + n[B]*Pb + w*P1); + + Faa = k1 * Paa; + Fbb = k1 * Pbb; + Fcc = k3 * (SQR(n[A])*Paa + 2*n[A]*n[B]*Pab + SQR(n[B])*Pbb + + w*(2*(n[A]*Pa + n[B]*Pb) + w*P1)); + + Faaa = k1 * Paaa; + Fbbb = k1 * Pbbb; + Fccc = -k4 * (CUBE(n[A])*Paaa + 3*SQR(n[A])*n[B]*Paab + + 3*n[A]*SQR(n[B])*Pabb + CUBE(n[B])*Pbbb + + 3*w*(SQR(n[A])*Paa + 2*n[A]*n[B]*Pab + SQR(n[B])*Pbb) + + w*w*(3*(n[A]*Pa + n[B]*Pb) + w*P1)); + + Faab = k1 * Paab; + Fbbc = -k2 * (n[A]*Pabb + n[B]*Pbbb + w*Pbb); + Fcca = k3 * (SQR(n[A])*Paaa + 2*n[A]*n[B]*Paab + SQR(n[B])*Pabb + + w*(2*(n[A]*Paa + n[B]*Pab) + w*Pa)); +} + +void compVolumeIntegrals(POLYHEDRON *p) +{ + FACE *f; + double nx, ny, nz; + int i; + + T0 = T1[X] = T1[Y] = T1[Z] + = T2[X] = T2[Y] = T2[Z] + = TP[X] = TP[Y] = TP[Z] = 0; + + for (i = 0; i < p->numFaces; i++) { + + f = &p->faces[i]; + + nx = fabs(f->norm[X]); + ny = fabs(f->norm[Y]); + nz = fabs(f->norm[Z]); + if (nx > ny && nx > nz) C = X; + else C = (ny > nz) ? Y : Z; + A = (C + 1) % 3; + B = (A + 1) % 3; + + compFaceIntegrals(f); + + T0 += f->norm[X] * ((A == X) ? Fa : ((B == X) ? Fb : Fc)); + + T1[A] += f->norm[A] * Faa; + T1[B] += f->norm[B] * Fbb; + T1[C] += f->norm[C] * Fcc; + T2[A] += f->norm[A] * Faaa; + T2[B] += f->norm[B] * Fbbb; + T2[C] += f->norm[C] * Fccc; + TP[A] += f->norm[A] * Faab; + TP[B] += f->norm[B] * Fbbc; + TP[C] += f->norm[C] * Fcca; + } + + T1[X] /= 2; T1[Y] /= 2; T1[Z] /= 2; + T2[X] /= 3; T2[Y] /= 3; T2[Z] /= 3; + TP[X] /= 2; TP[Y] /= 2; TP[Z] /= 2; +} + + +/* + ============================================================================ + main + ============================================================================ +*/ + + +int main(int argc, char *argv[]) +{ + POLYHEDRON p; + double density, mass; + double r[3]; /* center of mass */ + double J[3][3]; /* inertia tensor */ + + if (argc != 2) { + printf("usage: %s \n", argv[0]); + exit(0); + } + + readPolyhedron(argv[1], &p); + + compVolumeIntegrals(&p); + + + printf("\nT1 = %+20.6f\n\n", T0); + + printf("Tx = %+20.6f\n", T1[X]); + printf("Ty = %+20.6f\n", T1[Y]); + printf("Tz = %+20.6f\n\n", T1[Z]); + + printf("Txx = %+20.6f\n", T2[X]); + printf("Tyy = %+20.6f\n", T2[Y]); + printf("Tzz = %+20.6f\n\n", T2[Z]); + + printf("Txy = %+20.6f\n", TP[X]); + printf("Tyz = %+20.6f\n", TP[Y]); + printf("Tzx = %+20.6f\n\n", TP[Z]); + + density = 1.0; /* assume unit density */ + + mass = density * T0; + + /* compute center of mass */ + r[X] = T1[X] / T0; + r[Y] = T1[Y] / T0; + r[Z] = T1[Z] / T0; + + /* compute inertia tensor */ + J[X][X] = density * (T2[Y] + T2[Z]); + J[Y][Y] = density * (T2[Z] + T2[X]); + J[Z][Z] = density * (T2[X] + T2[Y]); + J[X][Y] = J[Y][X] = - density * TP[X]; + J[Y][Z] = J[Z][Y] = - density * TP[Y]; + J[Z][X] = J[X][Z] = - density * TP[Z]; + + /* translate inertia tensor to center of mass */ + J[X][X] -= mass * (r[Y]*r[Y] + r[Z]*r[Z]); + J[Y][Y] -= mass * (r[Z]*r[Z] + r[X]*r[X]); + J[Z][Z] -= mass * (r[X]*r[X] + r[Y]*r[Y]); + J[X][Y] = J[Y][X] += mass * r[X] * r[Y]; + J[Y][Z] = J[Z][Y] += mass * r[Y] * r[Z]; + J[Z][X] = J[X][Z] += mass * r[Z] * r[X]; + + printf("center of mass: (%+12.6f,%+12.6f,%+12.6f)\n\n", r[X], r[Y], r[Z]); + + printf("inertia tensor with origin at c.o.m. :\n"); + printf("%+15.6f %+15.6f %+15.6f\n", J[X][X], J[X][Y], J[X][Z]); + printf("%+15.6f %+15.6f %+15.6f\n", J[Y][X], J[Y][Y], J[Y][Z]); + printf("%+15.6f %+15.6f %+15.6f\n\n", J[Z][X], J[Z][Y], J[Z][Z]); + +} diff --git a/src/postProcessing/functionObjects/forces/Make/files b/src/postProcessing/functionObjects/forces/Make/files index 511568cab1..6538e133ef 100644 --- a/src/postProcessing/functionObjects/forces/Make/files +++ b/src/postProcessing/functionObjects/forces/Make/files @@ -4,10 +4,34 @@ forces/forcesFunctionObject.C forceCoeffs/forceCoeffs.C forceCoeffs/forceCoeffsFunctionObject.C +sDoFRBM = pointPatchFields/derived/sixDoFRigidBodyMotion + +$(sDoFRBM)/sixDoFRigidBodyMotion.C +$(sDoFRBM)/sixDoFRigidBodyMotionIO.C +$(sDoFRBM)/sixDoFRigidBodyMotionState.C +$(sDoFRBM)/sixDoFRigidBodyMotionStateIO.C + +sDoFRBMR = $(sDoFRBM)/sixDoFRigidBodyMotionRestraint + +$(sDoFRBMR)/sixDoFRigidBodyMotionRestraint/sixDoFRigidBodyMotionRestraint.C +$(sDoFRBMR)/sixDoFRigidBodyMotionRestraint/newSixDoFRigidBodyMotionRestraint.C +$(sDoFRBMR)/linearAxialAngularSpring/linearAxialAngularSpring.C +$(sDoFRBMR)/linearSpring/linearSpring.C +$(sDoFRBMR)/sphericalAngularSpring/sphericalAngularSpring.C +$(sDoFRBMR)/tabulatedAxialAngularSpring/tabulatedAxialAngularSpring.C + +sDoFRBMC = $(sDoFRBM)/sixDoFRigidBodyMotionConstraint + +$(sDoFRBMC)/sixDoFRigidBodyMotionConstraint/sixDoFRigidBodyMotionConstraint.C +$(sDoFRBMC)/sixDoFRigidBodyMotionConstraint/newSixDoFRigidBodyMotionConstraint.C +$(sDoFRBMC)/fixedAxis/fixedAxis.C +$(sDoFRBMC)/fixedLine/fixedLine.C +$(sDoFRBMC)/fixedOrientation/fixedOrientation.C +$(sDoFRBMC)/fixedPlane/fixedPlane.C +$(sDoFRBMC)/fixedPoint/fixedPoint.C + pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyDisplacementPointPatchVectorField.C -pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.C -pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionIO.C -pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionState.C -pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionStateIO.C +pointPatchFields/derived/uncoupledSixDoFRigidBodyDisplacement/uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField.C +pointPatchFields/derived/constrainedSixDoFRigidBodyDisplacement/constrainedSixDoFRigidBodyDisplacementPointPatchVectorField.C LIB = $(FOAM_LIBBIN)/libforces diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/constrainedSixDoFRigidBodyDisplacement/constrainedSixDoFRigidBodyDisplacementPointPatchVectorField.C b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/constrainedSixDoFRigidBodyDisplacement/constrainedSixDoFRigidBodyDisplacementPointPatchVectorField.C new file mode 100644 index 0000000000..01e47944f4 --- /dev/null +++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/constrainedSixDoFRigidBodyDisplacement/constrainedSixDoFRigidBodyDisplacementPointPatchVectorField.C @@ -0,0 +1,206 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd. + \\/ 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 2 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, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +\*---------------------------------------------------------------------------*/ + +#include "constrainedSixDoFRigidBodyDisplacementPointPatchVectorField.H" +#include "pointPatchFields.H" +#include "addToRunTimeSelectionTable.H" +#include "Time.H" +#include "fvMesh.H" +#include "volFields.H" +#include "uniformDimensionedFields.H" +#include "forces.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +constrainedSixDoFRigidBodyDisplacementPointPatchVectorField:: +constrainedSixDoFRigidBodyDisplacementPointPatchVectorField +( + const pointPatch& p, + const DimensionedField& iF +) +: + fixedValuePointPatchField(p, iF), + motion_(), + p0_(p.localPoints()), + rhoInf_(1.0) +{} + + +constrainedSixDoFRigidBodyDisplacementPointPatchVectorField:: +constrainedSixDoFRigidBodyDisplacementPointPatchVectorField +( + const pointPatch& p, + const DimensionedField& iF, + const dictionary& dict +) +: + fixedValuePointPatchField(p, iF, dict), + motion_(dict), + rhoInf_(readScalar(dict.lookup("rhoInf"))) +{ + if (!dict.found("value")) + { + updateCoeffs(); + } + + if (dict.found("p0")) + { + p0_ = vectorField("p0", dict , p.size()); + } + else + { + p0_ = p.localPoints(); + } +} + + +constrainedSixDoFRigidBodyDisplacementPointPatchVectorField:: +constrainedSixDoFRigidBodyDisplacementPointPatchVectorField +( + const constrainedSixDoFRigidBodyDisplacementPointPatchVectorField& ptf, + const pointPatch& p, + const DimensionedField& iF, + const pointPatchFieldMapper& mapper +) +: + fixedValuePointPatchField(ptf, p, iF, mapper), + motion_(ptf.motion_), + p0_(ptf.p0_, mapper), + rhoInf_(ptf.rhoInf_) +{} + + +constrainedSixDoFRigidBodyDisplacementPointPatchVectorField:: +constrainedSixDoFRigidBodyDisplacementPointPatchVectorField +( + const constrainedSixDoFRigidBodyDisplacementPointPatchVectorField& ptf, + const DimensionedField& iF +) +: + fixedValuePointPatchField(ptf, iF), + motion_(ptf.motion_), + p0_(ptf.p0_), + rhoInf_(ptf.rhoInf_) +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +void constrainedSixDoFRigidBodyDisplacementPointPatchVectorField::updateCoeffs() +{ + if (this->updated()) + { + return; + } + + const polyMesh& mesh = this->dimensionedInternalField().mesh()(); + const Time& t = mesh.time(); + const pointPatch& ptPatch = this->patch(); + + // 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()); + + dictionary forcesDict; + + forcesDict.add("patches", wordList(1, ptPatch.name())); + forcesDict.add("rhoInf", rhoInf_); + forcesDict.add("CofR", motion_.centreOfMass()); + + forces f("forces", db(), forcesDict); + + forces::forcesMoments fm = f.calcForcesMoment(); + + // Get the forces on the patch faces at the current positions + + vector gravity = vector::zero; + + if (db().foundObject("g")) + { + uniformDimensionedVectorField g = + db().lookupObject("g"); + + gravity = g.value(); + } + + vector rotationAxis(0, 1, 0); + + vector torque + ( + ( + (fm.second().first() + fm.second().second()) + & rotationAxis + ) + *rotationAxis + ); + + motion_.updateForce + ( + vector::zero, // Force no centre of mass motion + torque, // Only rotation allowed around the unit rotationAxis + t.deltaTValue() + ); + + Field::operator=(motion_.currentPosition(p0_) - p0_); + + fixedValuePointPatchField::updateCoeffs(); +} + + +void constrainedSixDoFRigidBodyDisplacementPointPatchVectorField::write +( + Ostream& os +) const +{ + pointPatchField::write(os); + motion_.write(os); + os.writeKeyword("rhoInf") + << rhoInf_ << token::END_STATEMENT << nl; + p0_.writeEntry("p0", os); + writeEntry("value", os); +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +makePointPatchTypeField +( + pointPatchVectorField, + constrainedSixDoFRigidBodyDisplacementPointPatchVectorField +); + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// ************************************************************************* // diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/constrainedSixDoFRigidBodyDisplacement/constrainedSixDoFRigidBodyDisplacementPointPatchVectorField.H b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/constrainedSixDoFRigidBodyDisplacement/constrainedSixDoFRigidBodyDisplacementPointPatchVectorField.H new file mode 100644 index 0000000000..a17b9aa20c --- /dev/null +++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/constrainedSixDoFRigidBodyDisplacement/constrainedSixDoFRigidBodyDisplacementPointPatchVectorField.H @@ -0,0 +1,157 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd. + \\/ 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 2 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, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +Class + Foam::constrainedSixDoFRigidBodyDisplacementPointPatchVectorField + +Description + Foam::constrainedSixDoFRigidBodyDisplacementPointPatchVectorField + +SourceFiles + constrainedSixDoFRigidBodyDisplacementPointPatchVectorField.C + +\*---------------------------------------------------------------------------*/ + +#ifndef constrainedSixDoFRigidBodyDisplacementPointPatchVectorField_H +#define constrainedSixDoFRigidBodyDisplacementPointPatchVectorField_H + +#include "fixedValuePointPatchField.H" +#include "sixDoFRigidBodyMotion.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class constrainedSixDoFRigidBodyDisplacementPointPatchVectorField Declaration +\*---------------------------------------------------------------------------*/ + +class constrainedSixDoFRigidBodyDisplacementPointPatchVectorField +: + public fixedValuePointPatchField +{ + // Private data + + //- Six dof motion object + sixDoFRigidBodyMotion motion_; + + //- Reference positions of points on the patch + pointField p0_; + + //- Reference density required by the forces object for + // incompressible calculations + scalar rhoInf_; + + +public: + + //- Runtime type information + TypeName("constrainedSixDoFRigidBodyDisplacement"); + + + // Constructors + + //- Construct from patch and internal field + constrainedSixDoFRigidBodyDisplacementPointPatchVectorField + ( + const pointPatch&, + const DimensionedField& + ); + + //- Construct from patch, internal field and dictionary + constrainedSixDoFRigidBodyDisplacementPointPatchVectorField + ( + const pointPatch&, + const DimensionedField&, + const dictionary& + ); + + //- Construct by mapping given patchField onto a new patch + constrainedSixDoFRigidBodyDisplacementPointPatchVectorField + ( + const constrainedSixDoFRigidBodyDisplacementPointPatchVectorField&, + const pointPatch&, + const DimensionedField&, + const pointPatchFieldMapper& + ); + + //- Construct and return a clone + virtual autoPtr > clone() const + { + return autoPtr > + ( + new constrainedSixDoFRigidBodyDisplacementPointPatchVectorField + ( + *this + ) + ); + } + + //- Construct as copy setting internal field reference + constrainedSixDoFRigidBodyDisplacementPointPatchVectorField + ( + const constrainedSixDoFRigidBodyDisplacementPointPatchVectorField&, + const DimensionedField& + ); + + //- Construct and return a clone setting internal field reference + virtual autoPtr > clone + ( + const DimensionedField& iF + ) const + { + return autoPtr > + ( + new constrainedSixDoFRigidBodyDisplacementPointPatchVectorField + ( + *this, + iF + ) + ); + } + + + // Member functions + + // Evaluation functions + + //- Update the coefficients associated with the patch field + virtual void updateCoeffs(); + + + //- Write + virtual void write(Ostream&) const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyDisplacementPointPatchVectorField.C b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyDisplacementPointPatchVectorField.C index 8fb0408e67..dce0389466 100644 --- a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyDisplacementPointPatchVectorField.C +++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyDisplacementPointPatchVectorField.C @@ -160,7 +160,27 @@ void sixDoFRigidBodyDisplacementPointPatchVectorField::updateCoeffs() t.deltaTValue() ); - Field::operator=(motion_.generatePositions(p0_) - p0_); + // ---------------------------------------- + // vector rotationAxis(0, 1, 0); + + // vector torque + // ( + // ( + // (fm.second().first() + fm.second().second()) + // & rotationAxis + // ) + // *rotationAxis + // ); + + // motion_.updateForce + // ( + // vector::zero, // Force no centre of mass motion + // torque, // Only rotation allowed around the unit rotationAxis + // t.deltaTValue() + // ); + // ---------------------------------------- + + Field::operator=(motion_.currentPosition(p0_) - p0_); fixedValuePointPatchField::updateCoeffs(); } diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.C b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.C deleted file mode 100644 index dae1ba9f3d..0000000000 --- a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.C +++ /dev/null @@ -1,210 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | OpenFOAM: The Open Source CFD Toolbox - \\ / O peration | - \\ / A nd | Copyright (C) 2009-2009 OpenCFD Ltd. - \\/ 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 2 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, write to the Free Software Foundation, - Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA - -\*---------------------------------------------------------------------------*/ - -#include "sixDoFRigidBodyMotion.H" - -// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // - -Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion() -: - motionState_(), - refCentreOfMass_(vector::zero), - momentOfInertia_(diagTensor::one*VSMALL), - mass_(VSMALL) -{} - - -Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion -( - const point& centreOfMass, - const tensor& Q, - const vector& v, - const vector& a, - const vector& pi, - const vector& tau, - scalar mass, - const point& refCentreOfMass, - const diagTensor& momentOfInertia -) -: - motionState_ - ( - centreOfMass, - Q, - v, - a, - pi, - tau - ), - refCentreOfMass_(refCentreOfMass), - momentOfInertia_(momentOfInertia), - mass_(mass) -{} - - -Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion(const dictionary& dict) -: - motionState_(dict), - refCentreOfMass_(dict.lookupOrDefault("refCentreOfMass", centreOfMass())), - momentOfInertia_(dict.lookup("momentOfInertia")), - mass_(readScalar(dict.lookup("mass"))) -{} - - -Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion -( - const sixDoFRigidBodyMotion& sDoFRBM -) -: - motionState_(sDoFRBM.motionState()), - refCentreOfMass_(sDoFRBM.refCentreOfMass()), - momentOfInertia_(sDoFRBM.momentOfInertia()), - mass_(sDoFRBM.mass()) -{} - - -// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // - -Foam::sixDoFRigidBodyMotion::~sixDoFRigidBodyMotion() -{} - - -// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * // - -void Foam::sixDoFRigidBodyMotion::updatePosition -( - scalar deltaT -) -{ - // First leapfrog velocity adjust and motion part, required before - // force calculation - - if (Pstream::master()) - { - v() += 0.5*deltaT*a(); - - pi() += 0.5*deltaT*tau(); - - // Leapfrog move part - centreOfMass() += deltaT*v(); - - // Leapfrog orientation adjustment - - tensor R; - - R = rotationTensorX(0.5*deltaT*pi().x()/momentOfInertia_.xx()); - pi() = pi() & R; - Q() = Q() & R; - - R = rotationTensorY(0.5*deltaT*pi().y()/momentOfInertia_.yy()); - pi() = pi() & R; - Q() = Q() & R; - - R = rotationTensorZ(deltaT*pi().z()/momentOfInertia_.zz()); - pi() = pi() & R; - Q() = Q() & R; - - R = rotationTensorY(0.5*deltaT*pi().y()/momentOfInertia_.yy()); - pi() = pi() & R; - Q() = Q() & R; - - R = rotationTensorX(0.5*deltaT*pi().x()/momentOfInertia_.xx()); - pi() = pi() & R; - Q() = Q() & R; - - } - - Pstream::scatter(motionState_); -} - - -void Foam::sixDoFRigidBodyMotion::updateForce -( - const vector& fGlobal, - const vector& tauGlobal, - scalar deltaT -) -{ - // Second leapfrog velocity adjust part, required after motion and - // force calculation part - - if (Pstream::master()) - { - a() = fGlobal/mass_; - - tau() = (Q().T() & tauGlobal); - - v() += 0.5*deltaT*a(); - - pi() += 0.5*deltaT*tau(); - } - - Pstream::scatter(motionState_); -} - - -void Foam::sixDoFRigidBodyMotion::updateForce -( - const pointField& positions, - const vectorField& forces, - scalar deltaT -) -{ - // Second leapfrog velocity adjust part, required after motion and - // force calculation part - - if (Pstream::master()) - { - a() = vector::zero; - - tau() = vector::zero; - - forAll(positions, i) - { - const vector& f = forces[i]; - - a() += f/mass_; - - tau() += (positions[i] ^ (Q().T() & f)); - } - - v() += 0.5*deltaT*a(); - - pi() += 0.5*deltaT*tau(); - } - - Pstream::scatter(motionState_); -} - - -Foam::tmp -Foam::sixDoFRigidBodyMotion::generatePositions(const pointField& pts) const -{ - return (centreOfMass() + (Q() & (pts - refCentreOfMass_))); -} - - -// ************************************************************************* // diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.C b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.C new file mode 100644 index 0000000000..1088ab1506 --- /dev/null +++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.C @@ -0,0 +1,447 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2009-2009 OpenCFD Ltd. + \\/ 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 2 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, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +\*---------------------------------------------------------------------------*/ + +#include "sixDoFRigidBodyMotion.H" + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +void Foam::sixDoFRigidBodyMotion::applyRestraints() +{ + if (Pstream::master()) + { + forAll(restraints_, rI) + { + // restraint position + point rP = vector::zero; + + // restraint force + vector rF = vector::zero; + + // restraint moment + vector rM = vector::zero; + + restraints_[rI].restrain(*this, rP, rF, rM); + + if (report_) + { + Info<< "Restraint " << restraints_[rI].name() << ": " + << "force " << rF << " moment " << rM + << endl; + } + + a() += rF/mass_; + + // Moments are returned in global axes, transforming to + // body local to add to torque. + tau() += Q().T() & (rM + ((rP - centreOfMass()) ^ rF)); + } + } +} + + +void Foam::sixDoFRigidBodyMotion::applyConstraints(scalar deltaT) +{ + if (Pstream::master()) + { + label iter = 0; + + bool allConverged = true; + + // constraint force accumulator + vector cFA = vector::zero; + + // constraint moment accumulator + vector cMA = vector::zero; + + do + { + allConverged = true; + + forAll(constraints_, cI) + { + // Info<< nl << "Constraint " << cI << endl; + + // constraint position + point cP = vector::zero; + + // constraint force + vector cF = vector::zero; + + // constraint moment + vector cM = vector::zero; + + bool constraintConverged = constraints_[cI].constrain + ( + *this, + cFA, + cMA, + deltaT, + cP, + cF, + cM + ); + + allConverged = allConverged && constraintConverged; + + // Accumulate constraint force + cFA += cF; + + // Accumulate constraint moment + cMA += cM + ((cP - centreOfMass()) ^ cF); + } + + } while(!allConverged && ++iter < maxConstraintIters_); + + if (iter >= maxConstraintIters_) + { + FatalErrorIn + ( + "Foam::sixDoFRigidBodyMotion::applyConstraints" + "(scalar deltaT)" + ) + << nl << "Maximum number of sixDoFRigidBodyMotion constraint " + << "iterations (" << maxConstraintIters_ << ") exceeded." << nl + << exit(FatalError); + } + else if (report_) + { + Info<< "sixDoFRigidBodyMotion constraints converged in " + << iter << " iterations" + << nl << "Constraint force: " << cFA << nl + << "Constraint moment: " << cMA + << endl; + } + + // Add the constraint forces and moments to the motion state variables + a() += cFA/mass_; + + // The moment of constraint forces has already been added + // during accumulation. Moments are returned in global axes, + // transforming to body local + tau() += Q().T() & cMA; + } +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion() +: + motionState_(), + restraints_(), + constraints_(), + maxConstraintIters_(0), + refCentreOfMass_(vector::zero), + momentOfInertia_(diagTensor::one*VSMALL), + mass_(VSMALL), + report_(false) +{} + + +Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion +( + const point& centreOfMass, + const tensor& Q, + const vector& v, + const vector& a, + const vector& pi, + const vector& tau, + scalar mass, + const point& refCentreOfMass, + const diagTensor& momentOfInertia, + bool report +) +: + motionState_ + ( + centreOfMass, + Q, + v, + a, + pi, + tau + ), + restraints_(), + constraints_(), + maxConstraintIters_(0), + refCentreOfMass_(refCentreOfMass), + momentOfInertia_(momentOfInertia), + mass_(mass), + report_(report) +{} + + +Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion(const dictionary& dict) +: + motionState_(dict), + restraints_(), + constraints_(), + maxConstraintIters_(0), + refCentreOfMass_(dict.lookupOrDefault("refCentreOfMass", centreOfMass())), + momentOfInertia_(dict.lookup("momentOfInertia")), + mass_(readScalar(dict.lookup("mass"))), + report_(dict.lookupOrDefault("report", false)) +{ + addRestraints(dict); + + addConstraints(dict); +} + + +Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion +( + const sixDoFRigidBodyMotion& sDoFRBM +) +: + motionState_(sDoFRBM.motionState()), + restraints_(sDoFRBM.restraints()), + constraints_(sDoFRBM.constraints()), + maxConstraintIters_(sDoFRBM.maxConstraintIters()), + refCentreOfMass_(sDoFRBM.refCentreOfMass()), + momentOfInertia_(sDoFRBM.momentOfInertia()), + mass_(sDoFRBM.mass()), + report_(sDoFRBM.report()) +{} + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +Foam::sixDoFRigidBodyMotion::~sixDoFRigidBodyMotion() +{} + + +// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * // + +void Foam::sixDoFRigidBodyMotion::addRestraints +( + const dictionary& dict +) +{ + if (dict.found("restraints")) + { + const dictionary& restraintDict = dict.subDict("restraints"); + + label i = 0; + + restraints_.setSize(restraintDict.size()); + + forAllConstIter(IDLList, restraintDict, iter) + { + if (iter().isDict()) + { + Info<< "Adding restraint: " << iter().keyword() << endl; + + restraints_.set + ( + i, + sixDoFRigidBodyMotionRestraint::New(iter().dict()) + ); + } + + i++; + } + + restraints_.setSize(i); + } +} + + +void Foam::sixDoFRigidBodyMotion::addConstraints +( + const dictionary& dict +) +{ + if (dict.found("constraints")) + { + const dictionary& constraintDict = dict.subDict("constraints"); + + label i = 0; + + constraints_.setSize(constraintDict.size()); + + forAllConstIter(IDLList, constraintDict, iter) + { + if (iter().isDict()) + { + Info<< "Adding constraint: " << iter().keyword() << endl; + + constraints_.set + ( + i++, + sixDoFRigidBodyMotionConstraint::New(iter().dict()) + ); + } + } + + constraints_.setSize(i); + + if (constraints_.size()) + { + maxConstraintIters_ = readLabel + ( + constraintDict.lookup("maxIterations") + ); + } + } +} + + +void Foam::sixDoFRigidBodyMotion::updatePosition +( + scalar deltaT +) +{ + // First leapfrog velocity adjust and motion part, required before + // force calculation + + if (Pstream::master()) + { + v() += 0.5*deltaT*a(); + + pi() += 0.5*deltaT*tau(); + + // Leapfrog move part + centreOfMass() += deltaT*v(); + + // Leapfrog orientation adjustment + + rotate(Q(), pi(), deltaT); + } + + Pstream::scatter(motionState_); +} + + +void Foam::sixDoFRigidBodyMotion::updateForce +( + const vector& fGlobal, + const vector& tauGlobal, + scalar deltaT +) +{ + // Second leapfrog velocity adjust part, required after motion and + // force calculation + + if (Pstream::master()) + { + a() = fGlobal/mass_; + + tau() = (Q().T() & tauGlobal); + + applyRestraints(); + + applyConstraints(deltaT); + + v() += 0.5*deltaT*a(); + + pi() += 0.5*deltaT*tau(); + + if(report_) + { + status(); + } + } + + Pstream::scatter(motionState_); +} + + +void Foam::sixDoFRigidBodyMotion::updateForce +( + const pointField& positions, + const vectorField& forces, + scalar deltaT +) +{ + vector a = vector::zero; + + vector tau = vector::zero; + + if (Pstream::master()) + { + forAll(positions, i) + { + const vector& f = forces[i]; + + a += f/mass_; + + tau += Q().T() & ((positions[i] - centreOfMass()) ^ f); + } + } + + updateForce(a, tau, deltaT); +} + + +Foam::point Foam::sixDoFRigidBodyMotion::predictedPosition +( + const point& pt, + const vector& deltaForce, + const vector& deltaMoment, + scalar deltaT +) 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); + + return (centreOfMassTemp + (QTemp & (pt - refCentreOfMass_))); +} + + +Foam::vector Foam::sixDoFRigidBodyMotion::predictedOrientation +( + const vector& v, + const vector& deltaMoment, + scalar deltaT +) const +{ + vector piTemp = pi() + deltaT*(tau() + (Q().T() & deltaMoment)); + + tensor QTemp = Q(); + + rotate(QTemp, piTemp, deltaT); + + return (QTemp & v); +} + + +void Foam::sixDoFRigidBodyMotion::status() const +{ + Info<< "Centre of mass: " << centreOfMass() << nl + << "Linear velocity: " << v() << nl + << "Angular velocity: " << omega() + << endl; +} + + +// ************************************************************************* // diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.H b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.H similarity index 60% rename from src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.H rename to src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.H index 16ea1d42b2..ecde4ac6d7 100644 --- a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.H +++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.H @@ -26,10 +26,10 @@ Class Foam::sixDoFRigidBodyMotion Description - Six degree of freedom motion for a rigid body. Angular momentum stored in - body fixed reference frame. Reference orientation of the body must align - with the cartesian axes such that the Inertia tensor is in principle - component form. + Six degree of freedom motion for a rigid body. Angular momentum + stored in body fixed reference frame. Reference orientation of + the body (where Q = I) must align with the cartesian axes such + that the Inertia tensor is in principle component form. Symplectic motion as per: @@ -43,6 +43,9 @@ Description url = {http://link.aip.org/link/?JCP/107/5840/1}, doi = {10.1063/1.474310} + Can add restraints (i.e. a spring) and constraints (i.e. motion + may only be on a plane). + SourceFiles sixDoFRigidBodyMotionI.H sixDoFRigidBodyMotion.C @@ -55,6 +58,9 @@ SourceFiles #include "sixDoFRigidBodyMotionState.H" #include "pointField.H" +#include "sixDoFRigidBodyMotionRestraint.H" +#include "sixDoFRigidBodyMotionConstraint.H" + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -79,9 +85,19 @@ class sixDoFRigidBodyMotion { // Private data - // state data object + //- Motion state data object sixDoFRigidBodyMotionState motionState_; + //- Restraints on the motion + PtrList restraints_; + + //- Constaints on the motion + PtrList constraints_; + + //- Maximum number of iterations allowed to attempt to obey + // constraints + label maxConstraintIters_; + //- Centre of mass of reference state point refCentreOfMass_; @@ -91,6 +107,9 @@ class sixDoFRigidBodyMotion //- Mass of the body scalar mass_; + //- Switch to turn reporting of motion data on and off + Switch report_; + // Private Member Functions @@ -106,6 +125,74 @@ 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 the restraints to the object + void applyRestraints(); + + //- 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 + + // Access + + //- Return access to the motion state + inline const sixDoFRigidBodyMotionState& motionState() const; + + //- Return access to the restraints + inline const PtrList& + restraints() const; + + //- Return access to the constraints + inline const PtrList& + constraints() const; + + //- Return access to the maximum allowed number of + // constraint iterations + inline label maxConstraintIters() const; + + //- Return access to the centre of mass + inline const point& refCentreOfMass() const; + + //- Return access to the orientation + inline const tensor& Q() const; + + //- Return access to velocity + inline const vector& v() const; + + //- Return access to acceleration + inline const vector& a() const; + + //- Return access to angular momentum + inline const vector& pi() const; + + //- Return access to torque + inline const vector& tau() const; + + + // Edit + + //- Return access to the centre of mass + inline point& refCentreOfMass(); + + //- Return non-const access to the orientation + inline tensor& Q(); + + //- Return non-const access to vector + inline vector& v(); + + //- Return non-const access to acceleration + inline vector& a(); + + //- Return non-const access to angular momentum + inline vector& pi(); + + //- Return non-const access to torque + inline vector& tau(); + public: @@ -125,7 +212,8 @@ public: const vector& tau, scalar mass, const point& refCentreOfMass, - const diagTensor& momentOfInertia + const diagTensor& momentOfInertia, + bool report = false ); //- Construct from dictionary @@ -141,11 +229,23 @@ public: // Member Functions + //- Add restraints to the motion, public to allow external + // addition of restraints after construction + void addRestraints(const dictionary& dict); + + //- Add restraints to the motion, public to allow external + // addition of restraints after construction + void addConstraints(const dictionary& dict); + + //- First leapfrog velocity adjust and motion part, required + // before force calculation void updatePosition ( scalar deltaT ); + //- Second leapfrog velocity adjust part, required after motion and + // force calculation void updateForce ( const vector& fGlobal, @@ -153,6 +253,8 @@ public: scalar deltaT ); + //- Global forces supplied at locations, calculating net force + // and moment void updateForce ( const pointField& positions, @@ -160,39 +262,68 @@ public: scalar deltaT ); - tmp generatePositions(const pointField& pts) const; + //- Transform the given reference state pointField by the + // current motion state + inline tmp currentPosition(const pointField& refPts) const; + + //- Transform the given reference state point by the current + // motion state + inline point currentPosition(const point& refPt) const; + + //- Transform the given reference state direction by the current + // motion state + inline vector currentOrientation(const vector& refDir) const; + + //- Access the orientation tensor, Q. + // globalVector = Q & bodyLocalVector + // bodyLocalVector = Q.T() & globalVector + inline const tensor& currentOrientation() const; + + //- Predict the position of the supplied point after deltaT + // given the current motion state and the additional supplied + // force and moment + point predictedPosition + ( + const point& pt, + const vector& deltaForce, + const vector& deltaMoment, + scalar deltaT + ) const; + + //- Predict the orientation of the supplied vector after deltaT + // given the current motion state and the additional supplied + // moment + vector predictedOrientation + ( + const vector& v, + const vector& deltaMoment, + scalar deltaT + ) const; + + //- Return the angular velocity in the global frame + inline vector omega() const; + + //- Return the velocity of a position given by the current + // motion state + inline point currentVelocity(const point& pt) const; + + //- Report the status of the motion + void status() const; + // Access - //- Return access to the motion state - inline const sixDoFRigidBodyMotionState& motionState() const; - - //- Return access to the centre of mass + //- Return const access to the centre of mass inline const point& centreOfMass() const; - //- Return access to the centre of mass - inline const point& refCentreOfMass() const; - //- Return access to the inertia tensor inline const diagTensor& momentOfInertia() const; - //- Return access to the mass + //- Return const access to the mass inline scalar mass() const; - //- Return access to the orientation - inline const tensor& Q() const; - - //- Return access to velocity - inline const vector& v() const; - - //- Return access to acceleration - inline const vector& a() const; - - //- Return access to angular momentum - inline const vector& pi() const; - - //- Return access to torque - inline const vector& tau() const; + //- Return the report Switch + inline bool report() const; // Edit @@ -200,30 +331,12 @@ public: //- Return non-const access to the centre of mass inline point& centreOfMass(); - //- Return access to the centre of mass - inline point& refCentreOfMass(); - //- Return non-const access to the inertia tensor inline diagTensor& momentOfInertia(); //- Return non-const access to the mass inline scalar& mass(); - //- Return non-const access to the orientation - inline tensor& Q(); - - //- Return non-const access to vector - inline vector& v(); - - //- Return non-const access to acceleration - inline vector& a(); - - //- Return non-const access to angular momentum - inline vector& pi(); - - //- Return non-const access to torque - inline vector& tau(); - //- Write void write(Ostream&) const; diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedAxis/fixedAxis.C b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedAxis/fixedAxis.C new file mode 100644 index 0000000000..996b76b06c --- /dev/null +++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedAxis/fixedAxis.C @@ -0,0 +1,157 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd. + \\/ 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 2 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, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +\*---------------------------------------------------------------------------*/ + +#include "fixedAxis.H" +#include "addToRunTimeSelectionTable.H" +#include "sixDoFRigidBodyMotion.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ +namespace sixDoFRigidBodyMotionConstraints +{ + defineTypeNameAndDebug(fixedAxis, 0); + addToRunTimeSelectionTable + ( + sixDoFRigidBodyMotionConstraint, + fixedAxis, + dictionary + ); +}; +}; + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::sixDoFRigidBodyMotionConstraints::fixedAxis::fixedAxis +( + const dictionary& sDoFRBMCDict +) +: + sixDoFRigidBodyMotionConstraint(sDoFRBMCDict), + fixedAxis_() +{ + read(sDoFRBMCDict); +} + + +// * * * * * * * * * * * * * * * * Destructors * * * * * * * * * * * * * * * // + +Foam::sixDoFRigidBodyMotionConstraints::fixedAxis::~fixedAxis() +{} + + +// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * // + +bool Foam::sixDoFRigidBodyMotionConstraints::fixedAxis::constrain +( + const sixDoFRigidBodyMotion& motion, + const vector& existingConstraintForce, + const vector& existingConstraintMoment, + scalar deltaT, + vector& constraintPosition, + vector& constraintForceIncrement, + vector& constraintMomentIncrement +) const +{ + constraintMomentIncrement = vector::zero; + + vector predictedDir = motion.predictedOrientation + ( + fixedAxis_, + existingConstraintMoment, + deltaT + ); + + scalar theta = acos(predictedDir & fixedAxis_); + + vector rotationAxis = fixedAxis_ ^ predictedDir; + + scalar magRotationAxis = mag(rotationAxis); + + if (magRotationAxis > VSMALL) + { + rotationAxis /= magRotationAxis; + + const tensor& Q = motion.currentOrientation(); + + // Transform rotationAxis to body local system + rotationAxis = Q.T() & rotationAxis; + + constraintMomentIncrement = + -relaxationFactor_ + *(motion.momentOfInertia() & rotationAxis) + *theta/sqr(deltaT); + + // Transform moment increment to global system + constraintMomentIncrement = Q & constraintMomentIncrement; + + // Remove any moment that is around the fixedAxis_ + constraintMomentIncrement -= + (constraintMomentIncrement & fixedAxis_)*fixedAxis_; + } + + constraintPosition = motion.centreOfMass(); + + constraintForceIncrement = vector::zero; + + return (mag(theta) < tolerance_); +} + + +bool Foam::sixDoFRigidBodyMotionConstraints::fixedAxis::read +( + const dictionary& sDoFRBMCDict +) +{ + sixDoFRigidBodyMotionConstraint::read(sDoFRBMCDict); + + sDoFRBMCCoeffs_.lookup("axis") >> fixedAxis_; + + scalar magFixedAxis(mag(fixedAxis_)); + + if (magFixedAxis > VSMALL) + { + fixedAxis_ /= magFixedAxis; + } + else + { + FatalErrorIn + ( + "Foam::sixDoFRigidBodyMotionConstraints::fixedAxis::read" + "(" + "const dictionary& sDoFRBMCDict" + ")" + ) + << "axis has zero length" + << abort(FatalError); + } + + return true; +} + +// ************************************************************************* // diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedAxis/fixedAxis.H b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedAxis/fixedAxis.H new file mode 100644 index 0000000000..5ba3dcb60d --- /dev/null +++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedAxis/fixedAxis.H @@ -0,0 +1,126 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd. + \\/ 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 2 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, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +Class + Foam::sixDoFRigidBodyMotionConstraints::fixedAxis + +Description + sixDoFRigidBodyMotionConstraint. Axis of body fixed global + space. + +SourceFiles + fixedAxis.C + +\*---------------------------------------------------------------------------*/ + +#ifndef fixedAxis_H +#define fixedAxis_H + +#include "sixDoFRigidBodyMotionConstraint.H" +#include "point.H" +#include "tensor.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +namespace sixDoFRigidBodyMotionConstraints +{ + +/*---------------------------------------------------------------------------*\ + Class fixedAxis Declaration +\*---------------------------------------------------------------------------*/ + +class fixedAxis +: + public sixDoFRigidBodyMotionConstraint +{ + + // Private data + + //- Reference axis in global space + vector fixedAxis_; + + +public: + + //- Runtime type information + TypeName("fixedAxis"); + + + // Constructors + + //- Construct from components + fixedAxis + ( + const dictionary& sDoFRBMCDict + ); + + //- Construct and return a clone + virtual autoPtr clone() const + { + return autoPtr + ( + new fixedAxis(*this) + ); + } + + + // Destructor + + virtual ~fixedAxis(); + + + // Member Functions + + //- Calculate the constraint position, force and moment. + // Global reference frame vectors. Returns boolean stating + // whether the constraint been converged to tolerance. + virtual bool constrain + ( + const sixDoFRigidBodyMotion& motion, + const vector& existingConstraintForce, + const vector& existingConstraintMoment, + scalar deltaT, + vector& constraintPosition, + vector& constraintForceIncrement, + vector& constraintMomentIncrement + ) const; + + //- Update properties from given dictionary + virtual bool read(const dictionary& sDoFRBMCCoeff); +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace solidBodyMotionFunctions +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedLine/fixedLine.C b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedLine/fixedLine.C new file mode 100644 index 0000000000..b5cee9ecc2 --- /dev/null +++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedLine/fixedLine.C @@ -0,0 +1,145 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd. + \\/ 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 2 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, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +\*---------------------------------------------------------------------------*/ + +#include "fixedLine.H" +#include "addToRunTimeSelectionTable.H" +#include "sixDoFRigidBodyMotion.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ +namespace sixDoFRigidBodyMotionConstraints +{ + defineTypeNameAndDebug(fixedLine, 0); + addToRunTimeSelectionTable + ( + sixDoFRigidBodyMotionConstraint, + fixedLine, + dictionary + ); +}; +}; + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::sixDoFRigidBodyMotionConstraints::fixedLine::fixedLine +( + const dictionary& sDoFRBMCDict +) +: + sixDoFRigidBodyMotionConstraint(sDoFRBMCDict), + refPt_(), + dir_() +{ + read(sDoFRBMCDict); +} + + +// * * * * * * * * * * * * * * * * Destructors * * * * * * * * * * * * * * * // + +Foam::sixDoFRigidBodyMotionConstraints::fixedLine::~fixedLine() +{} + + +// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * // + +bool Foam::sixDoFRigidBodyMotionConstraints::fixedLine::constrain +( + const sixDoFRigidBodyMotion& motion, + const vector& existingConstraintForce, + const vector& existingConstraintMoment, + scalar deltaT, + vector& constraintPosition, + vector& constraintForceIncrement, + vector& constraintMomentIncrement +) const +{ + point predictedPosition = motion.predictedPosition + ( + refPt_, + existingConstraintForce, + existingConstraintMoment, + deltaT + ); + + constraintPosition = motion.currentPosition(refPt_); + + // Info<< "current position " << constraintPosition << nl + // << "next predictedPosition " << predictedPosition + // << endl; + + // Vector from reference point to predicted point + vector rC = predictedPosition - refPt_; + + vector error = rC - ((rC) & dir_)*dir_; + + // Info<< "error " << error << endl; + + constraintForceIncrement = + -relaxationFactor_*error*motion.mass()/sqr(deltaT); + + constraintMomentIncrement = vector::zero; + + return (mag(error) < tolerance_); +} + + +bool Foam::sixDoFRigidBodyMotionConstraints::fixedLine::read +( + const dictionary& sDoFRBMCDict +) +{ + sixDoFRigidBodyMotionConstraint::read(sDoFRBMCDict); + + sDoFRBMCCoeffs_.lookup("refPoint") >> refPt_; + + sDoFRBMCCoeffs_.lookup("direction") >> dir_; + + scalar magDir(mag(dir_)); + + if (magDir > VSMALL) + { + dir_ /= magDir; + } + else + { + FatalErrorIn + ( + "Foam::sixDoFRigidBodyMotionConstraints::fixedLine::read" + "(" + "const dictionary& sDoFRBMCDict" + ")" + ) + << "line direction has zero length" + << abort(FatalError); + } + + return true; +} + +// ************************************************************************* // diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedLine/fixedLine.H b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedLine/fixedLine.H new file mode 100644 index 0000000000..01cc8ff52f --- /dev/null +++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedLine/fixedLine.H @@ -0,0 +1,127 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd. + \\/ 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 2 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, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +Class + Foam::sixDoFRigidBodyMotionConstraints::fixedLine + +Description + sixDoFRigidBodyMotionConstraint. Reference point may only move + along a line. + +SourceFiles + fixedLine.C + +\*---------------------------------------------------------------------------*/ + +#ifndef fixedLine_H +#define fixedLine_H + +#include "sixDoFRigidBodyMotionConstraint.H" +#include "point.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +namespace sixDoFRigidBodyMotionConstraints +{ + +/*---------------------------------------------------------------------------*\ + Class fixedLine Declaration +\*---------------------------------------------------------------------------*/ + +class fixedLine +: + public sixDoFRigidBodyMotionConstraint +{ + // Private data + + //- Reference point for the constraining line + point refPt_; + + //- Direction of the constraining line + vector dir_; + + +public: + + //- Runtime type information + TypeName("fixedLine"); + + + // Constructors + + //- Construct from components + fixedLine + ( + const dictionary& sDoFRBMCDict + ); + + //- Construct and return a clone + virtual autoPtr clone() const + { + return autoPtr + ( + new fixedLine(*this) + ); + } + + + // Destructor + + virtual ~fixedLine(); + + + // Member Functions + + //- Calculate the constraint position, force and moment. + // Global reference frame vectors. Returns boolean stating + // whether the constraint been converged to tolerance. + virtual bool constrain + ( + const sixDoFRigidBodyMotion& motion, + const vector& existingConstraintForce, + const vector& existingConstraintMoment, + scalar deltaT, + vector& constraintPosition, + vector& constraintForceIncrement, + vector& constraintMomentIncrement + ) const; + + //- Update properties from given dictionary + virtual bool read(const dictionary& sDoFRBMCCoeff); +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace solidBodyMotionFunctions +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedOrientation/fixedOrientation.C b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedOrientation/fixedOrientation.C new file mode 100644 index 0000000000..03ba862b71 --- /dev/null +++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedOrientation/fixedOrientation.C @@ -0,0 +1,180 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd. + \\/ 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 2 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, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +\*---------------------------------------------------------------------------*/ + +#include "fixedOrientation.H" +#include "addToRunTimeSelectionTable.H" +#include "sixDoFRigidBodyMotion.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ +namespace sixDoFRigidBodyMotionConstraints +{ + defineTypeNameAndDebug(fixedOrientation, 0); + addToRunTimeSelectionTable + ( + sixDoFRigidBodyMotionConstraint, + fixedOrientation, + dictionary + ); +}; +}; + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::sixDoFRigidBodyMotionConstraints::fixedOrientation::fixedOrientation +( + const dictionary& sDoFRBMCDict +) +: + sixDoFRigidBodyMotionConstraint(sDoFRBMCDict), + fixedOrientation_() +{ + read(sDoFRBMCDict); +} + + +// * * * * * * * * * * * * * * * * Destructors * * * * * * * * * * * * * * * // + +Foam::sixDoFRigidBodyMotionConstraints::fixedOrientation::~fixedOrientation() +{} + + +// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * // + +bool Foam::sixDoFRigidBodyMotionConstraints::fixedOrientation::constrain +( + const sixDoFRigidBodyMotion& motion, + const vector& existingConstraintForce, + const vector& existingConstraintMoment, + scalar deltaT, + vector& constraintPosition, + vector& constraintForceIncrement, + vector& constraintMomentIncrement +) const +{ + constraintMomentIncrement = vector::zero; + + scalar maxTheta = -SMALL; + + for (direction cmpt=0; cmpt VSMALL) + { + predictedDir /= magPredictedDir; + + theta = acos(predictedDir & refDir); + + // Recalculating axis to give correct sign to angle + axis = (refDir ^ predictedDir); + + scalar magAxis = mag(axis); + + if (magAxis > VSMALL) + { + axis /= magAxis; + } + else + { + axis = vector::zero; + } + } + + if (theta > maxTheta) + { + maxTheta = theta; + } + + constraintMomentIncrement += + -relaxationFactor_ + *theta*axis + *motion.momentOfInertia()[cmpt]/sqr(deltaT); + } + + constraintPosition = motion.centreOfMass(); + + constraintForceIncrement = vector::zero; + + return (mag(maxTheta) < tolerance_); +} + + +bool Foam::sixDoFRigidBodyMotionConstraints::fixedOrientation::read +( + const dictionary& sDoFRBMCDict +) +{ + sixDoFRigidBodyMotionConstraint::read(sDoFRBMCDict); + + fixedOrientation_ = + sDoFRBMCCoeffs_.lookupOrDefault("fixedOrientation", I); + + if (mag(mag(fixedOrientation_) - sqrt(3.0)) > 1e-9) + { + FatalErrorIn + ( + "Foam::sixDoFRigidBodyMotionConstraints::fixedOrientation::read" + "(" + "const dictionary& sDoFRBMCDict" + ")" + ) + << "fixedOrientation " << fixedOrientation_ + << " is not a rotation tensor." + << exit(FatalError); + } + + return true; +} + +// ************************************************************************* // diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedOrientation/fixedOrientation.H b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedOrientation/fixedOrientation.H new file mode 100644 index 0000000000..0d46cddea9 --- /dev/null +++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedOrientation/fixedOrientation.H @@ -0,0 +1,127 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd. + \\/ 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 2 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, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +Class + Foam::sixDoFRigidBodyMotionConstraints::fixedOrientation + +Description + sixDoFRigidBodyMotionConstraint. Orientation of body fixed global + space. Only valid where the predicted deviation from alignment is + < 90 degrees. + +SourceFiles + fixedOrientation.C + +\*---------------------------------------------------------------------------*/ + +#ifndef fixedOrientation_H +#define fixedOrientation_H + +#include "sixDoFRigidBodyMotionConstraint.H" +#include "point.H" +#include "tensor.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +namespace sixDoFRigidBodyMotionConstraints +{ + +/*---------------------------------------------------------------------------*\ + Class fixedOrientation Declaration +\*---------------------------------------------------------------------------*/ + +class fixedOrientation +: + public sixDoFRigidBodyMotionConstraint +{ + + // Private data + + //- Reference orientation where there is no moment + tensor fixedOrientation_; + + +public: + + //- Runtime type information + TypeName("fixedOrientation"); + + + // Constructors + + //- Construct from components + fixedOrientation + ( + const dictionary& sDoFRBMCDict + ); + + //- Construct and return a clone + virtual autoPtr clone() const + { + return autoPtr + ( + new fixedOrientation(*this) + ); + } + + + // Destructor + + virtual ~fixedOrientation(); + + + // Member Functions + + //- Calculate the constraint position, force and moment. + // Global reference frame vectors. Returns boolean stating + // whether the constraint been converged to tolerance. + virtual bool constrain + ( + const sixDoFRigidBodyMotion& motion, + const vector& existingConstraintForce, + const vector& existingConstraintMoment, + scalar deltaT, + vector& constraintPosition, + vector& constraintForceIncrement, + vector& constraintMomentIncrement + ) const; + + //- Update properties from given dictionary + virtual bool read(const dictionary& sDoFRBMCCoeff); +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace solidBodyMotionFunctions +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedPlane/fixedPlane.C b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedPlane/fixedPlane.C new file mode 100644 index 0000000000..fad22e7c13 --- /dev/null +++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedPlane/fixedPlane.C @@ -0,0 +1,128 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd. + \\/ 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 2 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, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +\*---------------------------------------------------------------------------*/ + +#include "fixedPlane.H" +#include "addToRunTimeSelectionTable.H" +#include "sixDoFRigidBodyMotion.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ +namespace sixDoFRigidBodyMotionConstraints +{ + defineTypeNameAndDebug(fixedPlane, 0); + addToRunTimeSelectionTable + ( + sixDoFRigidBodyMotionConstraint, + fixedPlane, + dictionary + ); +}; +}; + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::sixDoFRigidBodyMotionConstraints::fixedPlane::fixedPlane +( + const dictionary& sDoFRBMCDict +) +: + sixDoFRigidBodyMotionConstraint(sDoFRBMCDict), + fixedPlane_(vector::one) +{ + read(sDoFRBMCDict); +} + + +// * * * * * * * * * * * * * * * * Destructors * * * * * * * * * * * * * * * // + +Foam::sixDoFRigidBodyMotionConstraints::fixedPlane::~fixedPlane() +{} + + +// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * // + +bool Foam::sixDoFRigidBodyMotionConstraints::fixedPlane::constrain +( + const sixDoFRigidBodyMotion& motion, + const vector& existingConstraintForce, + const vector& existingConstraintMoment, + scalar deltaT, + vector& constraintPosition, + vector& constraintForceIncrement, + vector& constraintMomentIncrement +) const +{ + const point& refPt = fixedPlane_.refPoint(); + + const vector& n = fixedPlane_.normal(); + + point predictedPosition = motion.predictedPosition + ( + refPt, + existingConstraintForce, + existingConstraintMoment, + deltaT + ); + + constraintPosition = motion.currentPosition(refPt); + + // Info<< "current position " << constraintPosition << nl + // << "next predictedPosition " << predictedPosition + // << endl; + + vector error = ((predictedPosition - refPt) & n)*n; + + // Info<< "error " << error << endl; + + constraintForceIncrement = + -relaxationFactor_*error*motion.mass()/sqr(deltaT); + + constraintMomentIncrement = vector::zero; + + return (mag(error) < tolerance_); +} + + +bool Foam::sixDoFRigidBodyMotionConstraints::fixedPlane::read +( + const dictionary& sDoFRBMCDict +) +{ + sixDoFRigidBodyMotionConstraint::read(sDoFRBMCDict); + + point refPt = sDoFRBMCCoeffs_.lookup("refPoint"); + + vector normal = sDoFRBMCCoeffs_.lookup("normal"); + + fixedPlane_ = plane(refPt, normal); + + return true; +} + +// ************************************************************************* // diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedPlane/fixedPlane.H b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedPlane/fixedPlane.H new file mode 100644 index 0000000000..45f35dbbb4 --- /dev/null +++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedPlane/fixedPlane.H @@ -0,0 +1,125 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd. + \\/ 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 2 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, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +Class + Foam::sixDoFRigidBodyMotionConstraints::fixedPlane + +Description + sixDoFRigidBodyMotionConstraint. Reference point may only move + along a plane. + +SourceFiles + fixedPlane.C + +\*---------------------------------------------------------------------------*/ + +#ifndef fixedPlane_H +#define fixedPlane_H + +#include "sixDoFRigidBodyMotionConstraint.H" +#include "plane.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +namespace sixDoFRigidBodyMotionConstraints +{ + +/*---------------------------------------------------------------------------*\ + Class fixedPlane Declaration +\*---------------------------------------------------------------------------*/ + +class fixedPlane +: + public sixDoFRigidBodyMotionConstraint +{ + // Private data + + //- Plane which the transformed reference point is constrained + // to move along + plane fixedPlane_; + + +public: + + //- Runtime type information + TypeName("fixedPlane"); + + + // Constructors + + //- Construct from components + fixedPlane + ( + const dictionary& sDoFRBMCDict + ); + + //- Construct and return a clone + virtual autoPtr clone() const + { + return autoPtr + ( + new fixedPlane(*this) + ); + } + + + // Destructor + + virtual ~fixedPlane(); + + + // Member Functions + + //- Calculate the constraint position, force and moment. + // Global reference frame vectors. Returns boolean stating + // whether the constraint been converged to tolerance. + virtual bool constrain + ( + const sixDoFRigidBodyMotion& motion, + const vector& existingConstraintForce, + const vector& existingConstraintMoment, + scalar deltaT, + vector& constraintPosition, + vector& constraintForceIncrement, + vector& constraintMomentIncrement + ) const; + + //- Update properties from given dictionary + virtual bool read(const dictionary& sDoFRBMCCoeff); +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace solidBodyMotionFunctions +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedPoint/fixedPoint.C b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedPoint/fixedPoint.C new file mode 100644 index 0000000000..cfe8b5421d --- /dev/null +++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedPoint/fixedPoint.C @@ -0,0 +1,133 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd. + \\/ 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 2 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, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +\*---------------------------------------------------------------------------*/ + +#include "fixedPoint.H" +#include "addToRunTimeSelectionTable.H" +#include "sixDoFRigidBodyMotion.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ +namespace sixDoFRigidBodyMotionConstraints +{ + defineTypeNameAndDebug(fixedPoint, 0); + addToRunTimeSelectionTable + ( + sixDoFRigidBodyMotionConstraint, + fixedPoint, + dictionary + ); +}; +}; + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::sixDoFRigidBodyMotionConstraints::fixedPoint::fixedPoint +( + const dictionary& sDoFRBMCDict +) +: + sixDoFRigidBodyMotionConstraint(sDoFRBMCDict), + fixedPoint_() +{ + read(sDoFRBMCDict); +} + + +// * * * * * * * * * * * * * * * * Destructors * * * * * * * * * * * * * * * // + +Foam::sixDoFRigidBodyMotionConstraints::fixedPoint::~fixedPoint() +{} + + +// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * // + +bool Foam::sixDoFRigidBodyMotionConstraints::fixedPoint::constrain +( + const sixDoFRigidBodyMotion& motion, + const vector& existingConstraintForce, + const vector& existingConstraintMoment, + scalar deltaT, + vector& constraintPosition, + vector& constraintForceIncrement, + vector& constraintMomentIncrement +) const +{ + point predictedPosition = motion.predictedPosition + ( + fixedPoint_, + existingConstraintForce, + existingConstraintMoment, + deltaT + ); + + constraintPosition = motion.currentPosition(fixedPoint_); + + // Info<< "current position " << constraintPosition << nl + // << "next predictedPosition " << predictedPosition + // << endl; + + vector error = predictedPosition - fixedPoint_; + + // Info<< "error " << error << endl; + + // Correction force derived from Lagrange multiplier: + // G = -lambda*grad(sigma) + // where + // sigma = mag(error) = 0 + // so + // grad(sigma) = error/mag(error) + // Solving for lambda using the SHAKE methodology gives + // lambda = mass*mag(error)/sqr(deltaT) + // This is only strictly applicable (i.e. will converge in one + // iteration) to constraints at the centre of mass. Everything + // else will need to iterate, and may need under-relaxed to be + // stable. + + constraintForceIncrement = + -relaxationFactor_*error*motion.mass()/sqr(deltaT); + + constraintMomentIncrement = vector::zero; + + return (mag(error) < tolerance_); +} + + +bool Foam::sixDoFRigidBodyMotionConstraints::fixedPoint::read +( + const dictionary& sDoFRBMCDict +) +{ + sixDoFRigidBodyMotionConstraint::read(sDoFRBMCDict); + + sDoFRBMCCoeffs_.lookup("fixedPoint") >> fixedPoint_; + + return true; +} + +// ************************************************************************* // diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedPoint/fixedPoint.H b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedPoint/fixedPoint.H new file mode 100644 index 0000000000..933ba2db09 --- /dev/null +++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedPoint/fixedPoint.H @@ -0,0 +1,126 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd. + \\/ 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 2 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, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +Class + Foam::sixDoFRigidBodyMotionConstraints::fixedPoint + +Description + sixDoFRigidBodyMotionConstraint. Point fixed in space. + +SourceFiles + fixedPoint.C + +\*---------------------------------------------------------------------------*/ + +#ifndef fixedPoint_H +#define fixedPoint_H + +#include "sixDoFRigidBodyMotionConstraint.H" +#include "point.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +namespace sixDoFRigidBodyMotionConstraints +{ + +/*---------------------------------------------------------------------------*\ + Class fixedPoint Declaration +\*---------------------------------------------------------------------------*/ + +class fixedPoint +: + public sixDoFRigidBodyMotionConstraint +{ + // Private data + + //- Point which is rigidly attached to the body and constrains + // it so that this point remains fixed. This serves as the + // reference point for displacements as well as the target + // position. + point fixedPoint_; + + +public: + + //- Runtime type information + TypeName("fixedPoint"); + + + // Constructors + + //- Construct from components + fixedPoint + ( + const dictionary& sDoFRBMCDict + ); + + //- Construct and return a clone + virtual autoPtr clone() const + { + return autoPtr + ( + new fixedPoint(*this) + ); + } + + + // Destructor + + virtual ~fixedPoint(); + + + // Member Functions + + //- Calculate the constraint position, force and moment. + // Global reference frame vectors. Returns boolean stating + // whether the constraint been converged to tolerance. + virtual bool constrain + ( + const sixDoFRigidBodyMotion& motion, + const vector& existingConstraintForce, + const vector& existingConstraintMoment, + scalar deltaT, + vector& constraintPosition, + vector& constraintForceIncrement, + vector& constraintMomentIncrement + ) const; + + //- Update properties from given dictionary + virtual bool read(const dictionary& sDoFRBMCCoeff); +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace solidBodyMotionFunctions +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/sixDoFRigidBodyMotionConstraint/newSixDoFRigidBodyMotionConstraint.C b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/sixDoFRigidBodyMotionConstraint/newSixDoFRigidBodyMotionConstraint.C new file mode 100644 index 0000000000..890a83b2c0 --- /dev/null +++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/sixDoFRigidBodyMotionConstraint/newSixDoFRigidBodyMotionConstraint.C @@ -0,0 +1,65 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd. + \\/ 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 2 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, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +\*---------------------------------------------------------------------------*/ + +#include "sixDoFRigidBodyMotionConstraint.H" + +// * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * * // + +Foam::autoPtr +Foam::sixDoFRigidBodyMotionConstraint::New(const dictionary& sDoFRBMCDict) +{ + word sixDoFRigidBodyMotionConstraintTypeName = + sDoFRBMCDict.lookup("sixDoFRigidBodyMotionConstraint"); + + Info<< "Selecting sixDoFRigidBodyMotionConstraint function " + << sixDoFRigidBodyMotionConstraintTypeName << endl; + + dictionaryConstructorTable::iterator cstrIter = + dictionaryConstructorTablePtr_->find + ( + sixDoFRigidBodyMotionConstraintTypeName + ); + + if (cstrIter == dictionaryConstructorTablePtr_->end()) + { + FatalErrorIn + ( + "sixDoFRigidBodyMotionConstraint::New" + "(" + "const dictionary& sDoFRBMCDict" + ")" + ) << "Unknown sixDoFRigidBodyMotionConstraint type " + << sixDoFRigidBodyMotionConstraintTypeName << endl << endl + << "Valid sixDoFRigidBodyMotionConstraints are : " << endl + << dictionaryConstructorTablePtr_->sortedToc() + << exit(FatalError); + } + + return autoPtr(cstrIter()(sDoFRBMCDict)); +} + + +// ************************************************************************* // diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/sixDoFRigidBodyMotionConstraint/sixDoFRigidBodyMotionConstraint.C b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/sixDoFRigidBodyMotionConstraint/sixDoFRigidBodyMotionConstraint.C new file mode 100644 index 0000000000..082bbb34af --- /dev/null +++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/sixDoFRigidBodyMotionConstraint/sixDoFRigidBodyMotionConstraint.C @@ -0,0 +1,86 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd. + \\/ 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 2 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, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +\*---------------------------------------------------------------------------*/ + +#include "sixDoFRigidBodyMotionConstraint.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +defineTypeNameAndDebug(Foam::sixDoFRigidBodyMotionConstraint, 0); + +defineRunTimeSelectionTable(Foam::sixDoFRigidBodyMotionConstraint, dictionary); + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::sixDoFRigidBodyMotionConstraint::sixDoFRigidBodyMotionConstraint +( + const dictionary& sDoFRBMCDict +) +: + name_(fileName(sDoFRBMCDict.name().name()).components(token::COLON).last()), + sDoFRBMCCoeffs_ + ( + sDoFRBMCDict.subDict + ( + word(sDoFRBMCDict.lookup("sixDoFRigidBodyMotionConstraint")) + + "Coeffs" + ) + ), + tolerance_(readScalar(sDoFRBMCDict.lookup("tolerance"))), + relaxationFactor_ + ( + sDoFRBMCDict.lookupOrDefault("relaxationFactor", 1) + ) +{} + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +Foam::sixDoFRigidBodyMotionConstraint::~sixDoFRigidBodyMotionConstraint() +{} + + +// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * // + +bool Foam::sixDoFRigidBodyMotionConstraint::read +( + const dictionary& sDoFRBMCDict +) +{ + tolerance_ = (readScalar(sDoFRBMCDict.lookup("tolerance"))); + + relaxationFactor_ = sDoFRBMCDict.lookupOrDefault + ( + "relaxationFactor", + 1 + ); + + sDoFRBMCCoeffs_ = sDoFRBMCDict.subDict(type() + "Coeffs"); + + return true; +} + + +// ************************************************************************* // diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/sixDoFRigidBodyMotionConstraint/sixDoFRigidBodyMotionConstraint.H b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/sixDoFRigidBodyMotionConstraint/sixDoFRigidBodyMotionConstraint.H new file mode 100644 index 0000000000..22e5770202 --- /dev/null +++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/sixDoFRigidBodyMotionConstraint/sixDoFRigidBodyMotionConstraint.H @@ -0,0 +1,186 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd. + \\/ 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 2 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, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +Namespace + Foam::sixDoFRigidBodyMotionConstraints + +Description + Namespace for six DoF motion constraints + + +Class + Foam::sixDoFRigidBodyMotionConstraint + +Description + Base class for defining constraints for sixDoF motions + +SourceFiles + sixDoFRigidBodyMotionConstraint.C + newDynamicFvMesh.C + +\*---------------------------------------------------------------------------*/ + +#ifndef sixDoFRigidBodyMotionConstraint_H +#define sixDoFRigidBodyMotionConstraint_H + +#include "Time.H" +#include "dictionary.H" +#include "autoPtr.H" +#include "vector.H" +#include "runTimeSelectionTables.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// Forward declaration of classes +class sixDoFRigidBodyMotion; + + +/*---------------------------------------------------------------------------*\ + Class sixDoFRigidBodyMotionConstraint Declaration +\*---------------------------------------------------------------------------*/ + +class sixDoFRigidBodyMotionConstraint +{ + +protected: + + // Protected data + + //- Name of the constraint in dictionary + word name_; + + //- Constraint model specific coefficient dictionary + dictionary sDoFRBMCCoeffs_; + + //- Solution tolerance. Meaning depends on model, usually an + // absolute distance or angle. + scalar tolerance_; + + //- Relaxation factor for solution, default to one + scalar relaxationFactor_; + + +public: + + //- Runtime type information + TypeName("sixDoFRigidBodyMotionConstraint"); + + + // Declare run-time constructor selection table + + declareRunTimeSelectionTable + ( + autoPtr, + sixDoFRigidBodyMotionConstraint, + dictionary, + (const dictionary& sDoFRBMCDict), + (sDoFRBMCDict) + ); + + + // Constructors + + //- Construct from the sDoFRBMCDict dictionary and Time + sixDoFRigidBodyMotionConstraint + ( + const dictionary& sDoFRBMCDict + ); + + //- Construct and return a clone + virtual autoPtr clone() const = 0; + + + // Selectors + + //- Select constructed from the sDoFRBMCDict dictionary and Time + static autoPtr New + ( + const dictionary& sDoFRBMCDict + ); + + + // Destructor + + virtual ~sixDoFRigidBodyMotionConstraint(); + + + // Member Functions + + //- Calculate the constraint position, force and moment. + // Global reference frame vectors. Returns boolean stating + // whether the constraint been converged to tolerance. + virtual bool constrain + ( + const sixDoFRigidBodyMotion& motion, + const vector& existingConstraintForce, + const vector& existingConstraintMoment, + scalar deltaT, + vector& constraintPosition, + vector& constraintForceIncrement, + vector& constraintMomentIncrement + ) const = 0; + + //- Update properties from given dictionary + virtual bool read(const dictionary& sDoFRBMCDict) = 0; + + // Access + + //- Return access to the name of the restraint + inline const word& name() const + { + return name_; + } + + // Return access to sDoFRBMCCoeffs + inline const dictionary& coeffDict() const + { + return sDoFRBMCCoeffs_; + } + + //- Return access to the tolerance + inline scalar tolerance() const + { + return tolerance_; + } + + //- Return access to the relaxationFactor + inline scalar relaxationFactor() const + { + return relaxationFactor_; + } +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionI.H b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionI.H similarity index 66% rename from src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionI.H rename to src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionI.H index aeb75bdea4..f112ef42ef 100644 --- a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionI.H +++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionI.H @@ -24,8 +24,6 @@ License \*---------------------------------------------------------------------------*/ -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // inline Foam::tensor @@ -64,6 +62,37 @@ Foam::sixDoFRigidBodyMotion::rotationTensorZ(scalar phi) const } +inline void Foam::sixDoFRigidBodyMotion::rotate +( + tensor& Q, + vector& pi, + scalar deltaT +) const +{ + tensor R; + + R = rotationTensorX(0.5*deltaT*pi.x()/momentOfInertia_.xx()); + pi = pi & R; + Q = Q & R; + + R = rotationTensorY(0.5*deltaT*pi.y()/momentOfInertia_.yy()); + pi = pi & R; + Q = Q & R; + + R = rotationTensorZ(deltaT*pi.z()/momentOfInertia_.zz()); + pi = pi & R; + Q = Q & R; + + R = rotationTensorY(0.5*deltaT*pi.y()/momentOfInertia_.yy()); + pi = pi & R; + Q = Q & R; + + R = rotationTensorX(0.5*deltaT*pi.x()/momentOfInertia_.xx()); + pi = pi & R; + Q = Q & R; +} + + inline const Foam::sixDoFRigidBodyMotionState& Foam::sixDoFRigidBodyMotion::motionState() const { @@ -71,9 +100,23 @@ Foam::sixDoFRigidBodyMotion::motionState() const } -inline const Foam::point& Foam::sixDoFRigidBodyMotion::centreOfMass() const +inline const Foam::PtrList& +Foam::sixDoFRigidBodyMotion::restraints() const { - return motionState_.centreOfMass(); + return restraints_; +} + + +inline const Foam::PtrList& +Foam::sixDoFRigidBodyMotion::constraints() const +{ + return constraints_; +} + + +inline Foam::label Foam::sixDoFRigidBodyMotion::maxConstraintIters() const +{ + return maxConstraintIters_; } @@ -83,19 +126,6 @@ inline const Foam::point& Foam::sixDoFRigidBodyMotion::refCentreOfMass() const } -inline const Foam::diagTensor& -Foam::sixDoFRigidBodyMotion::momentOfInertia() const -{ - return momentOfInertia_; -} - - -inline Foam::scalar Foam::sixDoFRigidBodyMotion::mass() const -{ - return mass_; -} - - inline const Foam::tensor& Foam::sixDoFRigidBodyMotion::Q() const { return motionState_.Q(); @@ -126,30 +156,12 @@ inline const Foam::vector& Foam::sixDoFRigidBodyMotion::tau() const } -inline Foam::point& Foam::sixDoFRigidBodyMotion::centreOfMass() -{ - return motionState_.centreOfMass(); -} - - inline Foam::point& Foam::sixDoFRigidBodyMotion::refCentreOfMass() { return refCentreOfMass_; } -inline Foam::diagTensor& Foam::sixDoFRigidBodyMotion::momentOfInertia() -{ - return momentOfInertia_; -} - - -inline Foam::scalar& Foam::sixDoFRigidBodyMotion::mass() -{ - return mass_; -} - - inline Foam::tensor& Foam::sixDoFRigidBodyMotion::Q() { return motionState_.Q(); @@ -180,4 +192,95 @@ inline Foam::vector& Foam::sixDoFRigidBodyMotion::tau() } +// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * // + +inline Foam::tmp +Foam::sixDoFRigidBodyMotion::currentPosition(const pointField& refPts) const +{ + return (centreOfMass() + (Q() & (refPts - refCentreOfMass_))); +} + + +inline Foam::point Foam::sixDoFRigidBodyMotion::currentPosition +( + const point& refPt +) const +{ + return (centreOfMass() + (Q() & (refPt - refCentreOfMass_))); +} + + +inline Foam::vector Foam::sixDoFRigidBodyMotion::currentOrientation +( + const vector& refDir +) const +{ + return (Q() & refDir); +} + + +inline const Foam::tensor& +Foam::sixDoFRigidBodyMotion::currentOrientation() const +{ + return Q(); +} + + +inline Foam::vector Foam::sixDoFRigidBodyMotion::omega() const +{ + return Q() & (inv(momentOfInertia_) & pi()); +} + + +inline Foam::point Foam::sixDoFRigidBodyMotion::currentVelocity +( + const point& pt +) const +{ + return (omega() ^ (pt - centreOfMass())) + v(); +} + + +inline const Foam::point& Foam::sixDoFRigidBodyMotion::centreOfMass() const +{ + return motionState_.centreOfMass(); +} + + +inline const Foam::diagTensor& +Foam::sixDoFRigidBodyMotion::momentOfInertia() const +{ + return momentOfInertia_; +} + + +inline Foam::scalar Foam::sixDoFRigidBodyMotion::mass() const +{ + return mass_; +} + + +inline bool Foam::sixDoFRigidBodyMotion::report() const +{ + return report_; +} + + +inline Foam::point& Foam::sixDoFRigidBodyMotion::centreOfMass() +{ + return motionState_.centreOfMass(); +} + + +inline Foam::diagTensor& Foam::sixDoFRigidBodyMotion::momentOfInertia() +{ + return momentOfInertia_; +} + + +inline Foam::scalar& Foam::sixDoFRigidBodyMotion::mass() +{ + return mass_; +} + // ************************************************************************* // diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionIO.C b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionIO.C similarity index 61% rename from src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionIO.C rename to src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionIO.C index 005cbbe6aa..55245c3bfe 100644 --- a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionIO.C +++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionIO.C @@ -39,6 +39,72 @@ void Foam::sixDoFRigidBodyMotion::write(Ostream& os) const << momentOfInertia_ << token::END_STATEMENT << nl; os.writeKeyword("mass") << mass_ << token::END_STATEMENT << nl; + os.writeKeyword("report") + << report_ << token::END_STATEMENT << nl; + + if (restraints_.size()) + { + dictionary restraintsDict; + + forAll(restraints_, rI) + { + word restraintType = restraints_[rI].type(); + + dictionary restraintDict; + + restraintDict.add("sixDoFRigidBodyMotionRestraint", restraintType); + + restraintDict.add + ( + word(restraintType + "Coeffs"), restraints_[rI].coeffDict() + ); + + restraintsDict.add(restraints_[rI].name(), restraintDict); + } + + os.writeKeyword("restraints") << restraintsDict; + } + + if (constraints_.size()) + { + dictionary constraintsDict; + + constraintsDict.add("maxIterations", maxConstraintIters_); + + forAll(constraints_, rI) + { + word constraintType = constraints_[rI].type(); + + dictionary constraintDict; + + constraintDict.add + ( + "sixDoFRigidBodyMotionConstraint", + constraintType + ); + + constraintDict.add + ( + "tolerance", + constraints_[rI].tolerance() + ); + + constraintDict.add + ( + "relaxationFactor", + constraints_[rI].relaxationFactor() + ); + + constraintDict.add + ( + word(constraintType + "Coeffs"), constraints_[rI].coeffDict() + ); + + constraintsDict.add(constraints_[rI].name(), constraintDict); + } + + os.writeKeyword("constraints") << constraintsDict; + } } @@ -71,7 +137,7 @@ Foam::Ostream& Foam::operator<< os << sDoFRBM.motionState() << token::SPACE << sDoFRBM.refCentreOfMass() << token::SPACE << sDoFRBM.momentOfInertia() - << token::SPACE << sDoFRBM.mass() ; + << token::SPACE << sDoFRBM.mass(); // Check state of Ostream os.check diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/linearAxialAngularSpring/linearAxialAngularSpring.C b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/linearAxialAngularSpring/linearAxialAngularSpring.C new file mode 100644 index 0000000000..bed025bae6 --- /dev/null +++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/linearAxialAngularSpring/linearAxialAngularSpring.C @@ -0,0 +1,193 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd. + \\/ 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 2 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, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +\*---------------------------------------------------------------------------*/ + +#include "linearAxialAngularSpring.H" +#include "addToRunTimeSelectionTable.H" +#include "sixDoFRigidBodyMotion.H" +#include "transform.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ +namespace sixDoFRigidBodyMotionRestraints +{ + defineTypeNameAndDebug(linearAxialAngularSpring, 0); + addToRunTimeSelectionTable + ( + sixDoFRigidBodyMotionRestraint, + linearAxialAngularSpring, + dictionary + ); +}; +}; + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::sixDoFRigidBodyMotionRestraints::linearAxialAngularSpring:: +linearAxialAngularSpring +( + const dictionary& sDoFRBMRDict +) +: + sixDoFRigidBodyMotionRestraint(sDoFRBMRDict), + refQ_(), + axis_(), + stiffness_(), + damping_() +{ + read(sDoFRBMRDict); +} + + +// * * * * * * * * * * * * * * * * Destructors * * * * * * * * * * * * * * * // + +Foam::sixDoFRigidBodyMotionRestraints::linearAxialAngularSpring:: +~linearAxialAngularSpring() +{} + + +// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * // + +void +Foam::sixDoFRigidBodyMotionRestraints::linearAxialAngularSpring::restrain +( + const sixDoFRigidBodyMotion& motion, + vector& restraintPosition, + vector& restraintForce, + vector& restraintMoment +) const +{ + vector refDir = rotationTensor(vector(1, 0 ,0), axis_) & vector(0, 1, 0); + + vector oldDir = refQ_ & refDir; + + vector newDir = motion.currentOrientation(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); + + vector oldDir = refQ_ & refDir; + + vector newDir = motion.currentOrientation(refDir); + } + + // Removing any axis component from oldDir and newDir and normalising + oldDir -= (axis_ & oldDir)*axis_; + oldDir /= mag(oldDir); + + newDir -= (axis_ & newDir)*axis_; + newDir /= mag(newDir); + + scalar theta = mag(acos(oldDir & newDir)); + + // Temporary axis with sign information. + vector a = (oldDir ^ newDir); + + // Remove any component that is not along axis that may creep in + a = (a & axis_)*axis_; + + scalar magA = mag(a); + + if (magA > VSMALL) + { + a /= magA; + } + else + { + a = vector::zero; + } + + // Damping of along axis angular velocity only + restraintMoment = -stiffness_*theta*a - damping_*(motion.omega() & a)*a; + + restraintForce = vector::zero; + + // Not needed to be altered as restraintForce is zero, but set to + // centreOfMass to be sure of no spurious moment + restraintPosition = motion.centreOfMass(); +} + + +bool Foam::sixDoFRigidBodyMotionRestraints::linearAxialAngularSpring::read +( + const dictionary& sDoFRBMRDict +) +{ + sixDoFRigidBodyMotionRestraint::read(sDoFRBMRDict); + + refQ_ = sDoFRBMRCoeffs_.lookupOrDefault("referenceOrientation", I); + + if (mag(mag(refQ_) - sqrt(3.0)) > 1e-9) + { + FatalErrorIn + ( + "Foam::sixDoFRigidBodyMotionRestraints::linearAxialAngularSpring::" + "read" + "(" + "const dictionary& sDoFRBMRDict" + ")" + ) + << "referenceOrientation " << refQ_ << " is not a rotation tensor. " + << "mag(referenceOrientation) - sqrt(3) = " + << mag(refQ_) - sqrt(3.0) << nl + << exit(FatalError); + } + + axis_ = sDoFRBMRCoeffs_.lookup("axis"); + + scalar magAxis(mag(axis_)); + + if (magAxis > VSMALL) + { + axis_ /= magAxis; + } + else + { + FatalErrorIn + ( + "Foam::sixDoFRigidBodyMotionRestraints::linearAxialAngularSpring::" + "read" + "(" + "const dictionary& sDoFRBMCDict" + ")" + ) + << "axis has zero length" + << abort(FatalError); + } + + sDoFRBMRCoeffs_.lookup("stiffness") >> stiffness_; + + sDoFRBMRCoeffs_.lookup("damping") >> damping_; + + return true; +} + +// ************************************************************************* // diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/linearAxialAngularSpring/linearAxialAngularSpring.H b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/linearAxialAngularSpring/linearAxialAngularSpring.H new file mode 100644 index 0000000000..780598afb2 --- /dev/null +++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/linearAxialAngularSpring/linearAxialAngularSpring.H @@ -0,0 +1,129 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd. + \\/ 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 2 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, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +Class + Foam::sixDoFRigidBodyMotionRestraints::linearAxialAngularSpring + +Description + sixDoFRigidBodyMotionRestraints model. Linear axial angular spring. + +SourceFiles + linearAxialAngularSpring.C + +\*---------------------------------------------------------------------------*/ + +#ifndef linearAxialAngularSpring_H +#define linearAxialAngularSpring_H + +#include "sixDoFRigidBodyMotionRestraint.H" +#include "point.H" +#include "tensor.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +namespace sixDoFRigidBodyMotionRestraints +{ + +/*---------------------------------------------------------------------------*\ + Class linearAxialAngularSpring Declaration +\*---------------------------------------------------------------------------*/ + +class linearAxialAngularSpring +: + public sixDoFRigidBodyMotionRestraint +{ + // Private data + + //- Reference orientation where there is no moment + tensor refQ_; + + //- Global unit axis around which the motion is sprung + vector axis_; + + //- Spring stiffness coefficient (Nm/rad) + scalar stiffness_; + + //- Damping coefficient (Nms/rad) + scalar damping_; + + +public: + + //- Runtime type information + TypeName("linearAxialAngularSpring"); + + + // Constructors + + //- Construct from components + linearAxialAngularSpring + ( + const dictionary& sDoFRBMRDict + ); + + //- Construct and return a clone + virtual autoPtr clone() const + { + return autoPtr + ( + new linearAxialAngularSpring(*this) + ); + } + + + // Destructor + + virtual ~linearAxialAngularSpring(); + + + // Member Functions + + //- Calculate the restraint position, force and moment. + // Global reference frame vectors. + virtual void restrain + ( + const sixDoFRigidBodyMotion& motion, + vector& restraintPosition, + vector& restraintForce, + vector& restraintMoment + ) const; + + //- Update properties from given dictionary + virtual bool read(const dictionary& sDoFRBMRCoeff); +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace solidBodyMotionFunctions +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/linearSpring/linearSpring.C b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/linearSpring/linearSpring.C new file mode 100644 index 0000000000..966e911c6e --- /dev/null +++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/linearSpring/linearSpring.C @@ -0,0 +1,119 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd. + \\/ 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 2 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, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +\*---------------------------------------------------------------------------*/ + +#include "linearSpring.H" +#include "addToRunTimeSelectionTable.H" +#include "sixDoFRigidBodyMotion.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ +namespace sixDoFRigidBodyMotionRestraints +{ + defineTypeNameAndDebug(linearSpring, 0); + addToRunTimeSelectionTable + ( + sixDoFRigidBodyMotionRestraint, + linearSpring, + dictionary + ); +}; +}; + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::sixDoFRigidBodyMotionRestraints::linearSpring::linearSpring +( + const dictionary& sDoFRBMRDict +) +: + sixDoFRigidBodyMotionRestraint(sDoFRBMRDict), + anchor_(), + refAttachmentPt_(), + stiffness_(), + damping_(), + restLength_() +{ + read(sDoFRBMRDict); +} + + +// * * * * * * * * * * * * * * * * Destructors * * * * * * * * * * * * * * * // + +Foam::sixDoFRigidBodyMotionRestraints::linearSpring::~linearSpring() +{} + + +// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * // + +void Foam::sixDoFRigidBodyMotionRestraints::linearSpring::restrain +( + const sixDoFRigidBodyMotion& motion, + vector& restraintPosition, + vector& restraintForce, + vector& restraintMoment +) const +{ + restraintPosition = motion.currentPosition(refAttachmentPt_); + + vector r = restraintPosition - anchor_; + + scalar magR = mag(r); + + // r is now the r unit vector + r /= magR; + + vector v = motion.currentVelocity(restraintPosition); + + restraintForce = -stiffness_*(magR - restLength_)*r - damping_*(r & v)*r; + + restraintMoment = vector::zero; +} + + +bool Foam::sixDoFRigidBodyMotionRestraints::linearSpring::read +( + const dictionary& sDoFRBMRDict +) +{ + sixDoFRigidBodyMotionRestraint::read(sDoFRBMRDict); + + sDoFRBMRCoeffs_.lookup("anchor") >> anchor_; + + sDoFRBMRCoeffs_.lookup("refAttachmentPt") >> refAttachmentPt_; + + sDoFRBMRCoeffs_.lookup("stiffness") >> stiffness_; + + sDoFRBMRCoeffs_.lookup("damping") >> damping_; + + sDoFRBMRCoeffs_.lookup("restLength") >> restLength_; + + return true; +} + +// ************************************************************************* // diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/linearSpring/linearSpring.H b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/linearSpring/linearSpring.H new file mode 100644 index 0000000000..c5dc9bc074 --- /dev/null +++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/linearSpring/linearSpring.H @@ -0,0 +1,132 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd. + \\/ 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 2 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, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +Class + Foam::sixDoFRigidBodyMotionRestraints::linearSpring + +Description + sixDoFRigidBodyMotionRestraints model. Linear spring. + +SourceFiles + linearSpring.C + +\*---------------------------------------------------------------------------*/ + +#ifndef linearSpring_H +#define linearSpring_H + +#include "sixDoFRigidBodyMotionRestraint.H" +#include "point.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +namespace sixDoFRigidBodyMotionRestraints +{ + +/*---------------------------------------------------------------------------*\ + Class linearSpring Declaration +\*---------------------------------------------------------------------------*/ + +class linearSpring +: + public sixDoFRigidBodyMotionRestraint +{ + // Private data + + //- Anchor point, where the spring is attached to an immovable + // object + point anchor_; + + //- Reference point of attachment to the solid body + point refAttachmentPt_; + + //- Spring stiffness coefficient (N/m) + scalar stiffness_; + + //- Damping coefficient (Ns/m) + scalar damping_; + + //- Rest length - length of spring when no forces are applied to it + scalar restLength_; + + +public: + + //- Runtime type information + TypeName("linearSpring"); + + + // Constructors + + //- Construct from components + linearSpring + ( + const dictionary& sDoFRBMRDict + ); + + //- Construct and return a clone + virtual autoPtr clone() const + { + return autoPtr + ( + new linearSpring(*this) + ); + } + + + // Destructor + + virtual ~linearSpring(); + + + // Member Functions + + //- Calculate the restraint position, force and moment. + // Global reference frame vectors. + virtual void restrain + ( + const sixDoFRigidBodyMotion& motion, + vector& restraintPosition, + vector& restraintForce, + vector& restraintMoment + ) const; + + //- Update properties from given dictionary + virtual bool read(const dictionary& sDoFRBMRCoeff); +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace solidBodyMotionFunctions +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/sixDoFRigidBodyMotionRestraint/newSixDoFRigidBodyMotionRestraint.C b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/sixDoFRigidBodyMotionRestraint/newSixDoFRigidBodyMotionRestraint.C new file mode 100644 index 0000000000..7fee5eb84e --- /dev/null +++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/sixDoFRigidBodyMotionRestraint/newSixDoFRigidBodyMotionRestraint.C @@ -0,0 +1,65 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd. + \\/ 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 2 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, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +\*---------------------------------------------------------------------------*/ + +#include "sixDoFRigidBodyMotionRestraint.H" + +// * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * * // + +Foam::autoPtr +Foam::sixDoFRigidBodyMotionRestraint::New(const dictionary& sDoFRBMRDict) +{ + word sixDoFRigidBodyMotionRestraintTypeName = + sDoFRBMRDict.lookup("sixDoFRigidBodyMotionRestraint"); + + Info<< "Selecting sixDoFRigidBodyMotionRestraint function " + << sixDoFRigidBodyMotionRestraintTypeName << endl; + + dictionaryConstructorTable::iterator cstrIter = + dictionaryConstructorTablePtr_->find + ( + sixDoFRigidBodyMotionRestraintTypeName + ); + + if (cstrIter == dictionaryConstructorTablePtr_->end()) + { + FatalErrorIn + ( + "sixDoFRigidBodyMotionRestraint::New" + "(" + "const dictionary& sDoFRBMRDict" + ")" + ) << "Unknown sixDoFRigidBodyMotionRestraint type " + << sixDoFRigidBodyMotionRestraintTypeName << endl << endl + << "Valid sixDoFRigidBodyMotionRestraints are : " << endl + << dictionaryConstructorTablePtr_->sortedToc() + << exit(FatalError); + } + + return autoPtr(cstrIter()(sDoFRBMRDict)); +} + + +// ************************************************************************* // diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/sixDoFRigidBodyMotionRestraint/sixDoFRigidBodyMotionRestraint.C b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/sixDoFRigidBodyMotionRestraint/sixDoFRigidBodyMotionRestraint.C new file mode 100644 index 0000000000..7cfe723c01 --- /dev/null +++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/sixDoFRigidBodyMotionRestraint/sixDoFRigidBodyMotionRestraint.C @@ -0,0 +1,73 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd. + \\/ 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 2 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, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +\*---------------------------------------------------------------------------*/ + +#include "sixDoFRigidBodyMotionRestraint.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +defineTypeNameAndDebug(Foam::sixDoFRigidBodyMotionRestraint, 0); + +defineRunTimeSelectionTable(Foam::sixDoFRigidBodyMotionRestraint, dictionary); + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::sixDoFRigidBodyMotionRestraint::sixDoFRigidBodyMotionRestraint +( + const dictionary& sDoFRBMRDict +) +: + name_(fileName(sDoFRBMRDict.name().name()).components(token::COLON).last()), + sDoFRBMRCoeffs_ + ( + sDoFRBMRDict.subDict + ( + word(sDoFRBMRDict.lookup("sixDoFRigidBodyMotionRestraint")) + + "Coeffs" + ) + ) +{} + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +Foam::sixDoFRigidBodyMotionRestraint::~sixDoFRigidBodyMotionRestraint() +{} + + +// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * // + +bool Foam::sixDoFRigidBodyMotionRestraint::read +( + const dictionary& sDoFRBMRDict +) +{ + sDoFRBMRCoeffs_ = sDoFRBMRDict.subDict(type() + "Coeffs"); + + return true; +} + + +// ************************************************************************* // diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/sixDoFRigidBodyMotionRestraint/sixDoFRigidBodyMotionRestraint.H b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/sixDoFRigidBodyMotionRestraint/sixDoFRigidBodyMotionRestraint.H new file mode 100644 index 0000000000..1bdabb57e3 --- /dev/null +++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/sixDoFRigidBodyMotionRestraint/sixDoFRigidBodyMotionRestraint.H @@ -0,0 +1,163 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd. + \\/ 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 2 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, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +Namespace + Foam::sixDoFRigidBodyMotionRestraints + +Description + Namespace for six DoF motion restraints + + +Class + Foam::sixDoFRigidBodyMotionRestraint + +Description + Base class for defining restraints for sixDoF motions + +SourceFiles + sixDoFRigidBodyMotionRestraint.C + newDynamicFvMesh.C + +\*---------------------------------------------------------------------------*/ + +#ifndef sixDoFRigidBodyMotionRestraint_H +#define sixDoFRigidBodyMotionRestraint_H + +#include "Time.H" +#include "dictionary.H" +#include "autoPtr.H" +#include "vector.H" +#include "runTimeSelectionTables.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// Forward declaration of classes +class sixDoFRigidBodyMotion; + + +/*---------------------------------------------------------------------------*\ + Class sixDoFRigidBodyMotionRestraint Declaration +\*---------------------------------------------------------------------------*/ + +class sixDoFRigidBodyMotionRestraint +{ + +protected: + + // Protected data + + //- Name of the constraint in dictionary + word name_; + + //- Restraint model specific coefficient dictionary + dictionary sDoFRBMRCoeffs_; + + +public: + + //- Runtime type information + TypeName("sixDoFRigidBodyMotionRestraint"); + + + // Declare run-time constructor selection table + + declareRunTimeSelectionTable + ( + autoPtr, + sixDoFRigidBodyMotionRestraint, + dictionary, + (const dictionary& sDoFRBMRDict), + (sDoFRBMRDict) + ); + + + // Constructors + + //- Construct from the sDoFRBMRDict dictionary and Time + sixDoFRigidBodyMotionRestraint + ( + const dictionary& sDoFRBMRDict + ); + + //- Construct and return a clone + virtual autoPtr clone() const = 0; + + + // Selectors + + //- Select constructed from the sDoFRBMRDict dictionary and Time + static autoPtr New + ( + const dictionary& sDoFRBMRDict + ); + + + // Destructor + + virtual ~sixDoFRigidBodyMotionRestraint(); + + + // Member Functions + + //- Calculate the restraint position, force and moment. + // Global reference frame vectors. + virtual void restrain + ( + const sixDoFRigidBodyMotion& motion, + vector& restraintPosition, + vector& restraintForce, + vector& restraintMoment + ) const = 0; + + //- Update properties from given dictionary + virtual bool read(const dictionary& sDoFRBMRDict) = 0; + + // Access + + //- Return access to the name of the restraint + inline const word& name() const + { + return name_; + } + + // Return access to sDoFRBMRCoeffs + inline const dictionary& coeffDict() const + { + return sDoFRBMRCoeffs_; + } +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/sphericalAngularSpring/sphericalAngularSpring.C b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/sphericalAngularSpring/sphericalAngularSpring.C new file mode 100644 index 0000000000..af11dcb421 --- /dev/null +++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/sphericalAngularSpring/sphericalAngularSpring.C @@ -0,0 +1,148 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd. + \\/ 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 2 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, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +\*---------------------------------------------------------------------------*/ + +#include "sphericalAngularSpring.H" +#include "addToRunTimeSelectionTable.H" +#include "sixDoFRigidBodyMotion.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ +namespace sixDoFRigidBodyMotionRestraints +{ + defineTypeNameAndDebug(sphericalAngularSpring, 0); + addToRunTimeSelectionTable + ( + sixDoFRigidBodyMotionRestraint, + sphericalAngularSpring, + dictionary + ); +}; +}; + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::sixDoFRigidBodyMotionRestraints::sphericalAngularSpring:: +sphericalAngularSpring +( + const dictionary& sDoFRBMRDict +) +: + sixDoFRigidBodyMotionRestraint(sDoFRBMRDict), + refQ_(), + stiffness_(), + damping_() +{ + read(sDoFRBMRDict); +} + + +// * * * * * * * * * * * * * * * * Destructors * * * * * * * * * * * * * * * // + +Foam::sixDoFRigidBodyMotionRestraints::sphericalAngularSpring:: +~sphericalAngularSpring() +{} + + +// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * // + +void +Foam::sixDoFRigidBodyMotionRestraints::sphericalAngularSpring::restrain +( + const sixDoFRigidBodyMotion& motion, + vector& restraintPosition, + vector& restraintForce, + vector& restraintMoment +) const +{ + restraintMoment = vector::zero; + + for (direction cmpt=0; cmpt("referenceOrientation", I); + + if (mag(mag(refQ_) - sqrt(3.0)) > 1e-9) + { + FatalErrorIn + ( + "Foam::sixDoFRigidBodyMotionRestraints::sphericalAngularSpring::" + "read" + "(" + "const dictionary& sDoFRBMRDict" + ")" + ) + << "referenceOrientation " << refQ_ << " is not a rotation tensor. " + << "mag(referenceOrientation) - sqrt(3) = " + << mag(refQ_) - sqrt(3.0) << nl + << exit(FatalError); + } + + sDoFRBMRCoeffs_.lookup("stiffness") >> stiffness_; + + sDoFRBMRCoeffs_.lookup("damping") >> damping_; + + return true; +} + +// ************************************************************************* // diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/sphericalAngularSpring/sphericalAngularSpring.H b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/sphericalAngularSpring/sphericalAngularSpring.H new file mode 100644 index 0000000000..6ee6636b29 --- /dev/null +++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/sphericalAngularSpring/sphericalAngularSpring.H @@ -0,0 +1,126 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd. + \\/ 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 2 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, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +Class + Foam::sixDoFRigidBodyMotionRestraints::sphericalAngularSpring + +Description + sixDoFRigidBodyMotionRestraints model. Spherical angular spring. + +SourceFiles + sphericalAngularSpring.C + +\*---------------------------------------------------------------------------*/ + +#ifndef sphericalAngularSpring_H +#define sphericalAngularSpring_H + +#include "sixDoFRigidBodyMotionRestraint.H" +#include "point.H" +#include "tensor.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +namespace sixDoFRigidBodyMotionRestraints +{ + +/*---------------------------------------------------------------------------*\ + Class sphericalAngularSpring Declaration +\*---------------------------------------------------------------------------*/ + +class sphericalAngularSpring +: + public sixDoFRigidBodyMotionRestraint +{ + // Private data + + //- Reference orientation where there is no moment + tensor refQ_; + + //- Spring stiffness coefficient (Nm/rad) + scalar stiffness_; + + //- Damping coefficient (Nms/rad) + scalar damping_; + + +public: + + //- Runtime type information + TypeName("sphericalAngularSpring"); + + + // Constructors + + //- Construct from components + sphericalAngularSpring + ( + const dictionary& sDoFRBMRDict + ); + + //- Construct and return a clone + virtual autoPtr clone() const + { + return autoPtr + ( + new sphericalAngularSpring(*this) + ); + } + + + // Destructor + + virtual ~sphericalAngularSpring(); + + + // Member Functions + + //- Calculate the restraint position, force and moment. + // Global reference frame vectors. + virtual void restrain + ( + const sixDoFRigidBodyMotion& motion, + vector& restraintPosition, + vector& restraintForce, + vector& restraintMoment + ) const; + + //- Update properties from given dictionary + virtual bool read(const dictionary& sDoFRBMRCoeff); +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace solidBodyMotionFunctions +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/tabulatedAxialAngularSpring/tabulatedAxialAngularSpring.C b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/tabulatedAxialAngularSpring/tabulatedAxialAngularSpring.C new file mode 100644 index 0000000000..c2e164ecd8 --- /dev/null +++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/tabulatedAxialAngularSpring/tabulatedAxialAngularSpring.C @@ -0,0 +1,230 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd. + \\/ 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 2 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, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +\*---------------------------------------------------------------------------*/ + +#include "tabulatedAxialAngularSpring.H" +#include "addToRunTimeSelectionTable.H" +#include "sixDoFRigidBodyMotion.H" +#include "transform.H" +#include "unitConversion.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ +namespace sixDoFRigidBodyMotionRestraints +{ + defineTypeNameAndDebug(tabulatedAxialAngularSpring, 0); + addToRunTimeSelectionTable + ( + sixDoFRigidBodyMotionRestraint, + tabulatedAxialAngularSpring, + dictionary + ); +}; +}; + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::sixDoFRigidBodyMotionRestraints::tabulatedAxialAngularSpring:: +tabulatedAxialAngularSpring +( + const dictionary& sDoFRBMRDict +) +: + sixDoFRigidBodyMotionRestraint(sDoFRBMRDict), + refQ_(), + axis_(), + stiffness_(), + convertToDegrees_(), + damping_() +{ + read(sDoFRBMRDict); +} + + +// * * * * * * * * * * * * * * * * Destructors * * * * * * * * * * * * * * * // + +Foam::sixDoFRigidBodyMotionRestraints::tabulatedAxialAngularSpring:: +~tabulatedAxialAngularSpring() +{} + + +// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * // + +void +Foam::sixDoFRigidBodyMotionRestraints::tabulatedAxialAngularSpring::restrain +( + const sixDoFRigidBodyMotion& motion, + vector& restraintPosition, + vector& restraintForce, + vector& restraintMoment +) const +{ + vector refDir = rotationTensor(vector(1, 0 ,0), axis_) & vector(0, 1, 0); + + vector oldDir = refQ_ & refDir; + + vector newDir = motion.currentOrientation(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); + + vector oldDir = refQ_ & refDir; + + vector newDir = motion.currentOrientation(refDir); + } + + // Removing any axis component from oldDir and newDir and normalising + oldDir -= (axis_ & oldDir)*axis_; + oldDir /= mag(oldDir); + + newDir -= (axis_ & newDir)*axis_; + newDir /= mag(newDir); + + scalar theta = mag(acos(oldDir & newDir)); + + // Temporary axis with sign information. + vector a = (oldDir ^ newDir); + + // Remove any component that is not along axis that may creep in + a = (a & axis_)*axis_; + + scalar magA = mag(a); + + if (magA > VSMALL) + { + a /= magA; + } + else + { + a = vector::zero; + } + + scalar stiffness; + + if (convertToDegrees_) + { + stiffness = stiffness_(radToDeg(theta)); + } + else + { + stiffness = stiffness_(theta); + } + + // Damping of along axis angular velocity only + restraintMoment = -stiffness*theta*a - damping_*(motion.omega() & a)*a; + + restraintForce = vector::zero; + + // Not needed to be altered as restraintForce is zero, but set to + // centreOfMass to be sure of no spurious moment + restraintPosition = motion.centreOfMass(); +} + + +bool Foam::sixDoFRigidBodyMotionRestraints::tabulatedAxialAngularSpring::read +( + const dictionary& sDoFRBMRDict +) +{ + sixDoFRigidBodyMotionRestraint::read(sDoFRBMRDict); + + refQ_ = sDoFRBMRCoeffs_.lookupOrDefault("referenceOrientation", I); + + if (mag(mag(refQ_) - sqrt(3.0)) > 1e-9) + { + FatalErrorIn + ( + "Foam::sixDoFRigidBodyMotionRestraints::" + "tabulatedAxialAngularSpring::read" + "(" + "const dictionary& sDoFRBMRDict" + ")" + ) + << "referenceOrientation " << refQ_ << " is not a rotation tensor. " + << "mag(referenceOrientation) - sqrt(3) = " + << mag(refQ_) - sqrt(3.0) << nl + << exit(FatalError); + } + + axis_ = sDoFRBMRCoeffs_.lookup("axis"); + + scalar magAxis(mag(axis_)); + + if (magAxis > VSMALL) + { + axis_ /= magAxis; + } + else + { + FatalErrorIn + ( + "Foam::sixDoFRigidBodyMotionRestraints::" + "tabulatedAxialAngularSpring::read" + "(" + "const dictionary& sDoFRBMCDict" + ")" + ) + << "axis has zero length" + << abort(FatalError); + } + + stiffness_ = interpolationTable(sDoFRBMRCoeffs_); + + word angleFormat = sDoFRBMRCoeffs_.lookup("angleFormat"); + + if (angleFormat == "degrees" || angleFormat == "degree") + { + convertToDegrees_ = true; + } + else if (angleFormat == "radians" || angleFormat == "radian") + { + convertToDegrees_ = false; + } + else + { + FatalErrorIn + ( + "Foam::sixDoFRigidBodyMotionRestraints::" + "tabulatedAxialAngularSpring::read" + "(" + "const dictionary& sDoFRBMCDict" + ")" + ) + << "angleFormat must be degree, degrees, radian or radians" + << abort(FatalError); + } + + sDoFRBMRCoeffs_.lookup("damping") >> damping_; + + return true; +} + +// ************************************************************************* // diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/tabulatedAxialAngularSpring/tabulatedAxialAngularSpring.H b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/tabulatedAxialAngularSpring/tabulatedAxialAngularSpring.H new file mode 100644 index 0000000000..61c2c64197 --- /dev/null +++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/tabulatedAxialAngularSpring/tabulatedAxialAngularSpring.H @@ -0,0 +1,137 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd. + \\/ 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 2 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, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +Class + Foam::sixDoFRigidBodyMotionRestraints::tabulatedAxialAngularSpring + +Description + sixDoFRigidBodyMotionRestraints model. Axial angular spring with stiffness + values drawn from an interpolation table. Linear damping. + +SourceFiles + tabulatedAxialAngularSpring.C + +\*---------------------------------------------------------------------------*/ + +#ifndef tabulatedAxialAngularSpring_H +#define tabulatedAxialAngularSpring_H + +#include "sixDoFRigidBodyMotionRestraint.H" +#include "point.H" +#include "tensor.H" +#include "interpolationTable.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +namespace sixDoFRigidBodyMotionRestraints +{ + +/*---------------------------------------------------------------------------*\ + Class tabulatedAxialAngularSpring Declaration +\*---------------------------------------------------------------------------*/ + +class tabulatedAxialAngularSpring +: + public sixDoFRigidBodyMotionRestraint +{ + // Private data + + //- Reference orientation where there is no moment + tensor refQ_; + + //- Global unit axis around which the motion is sprung + vector axis_; + + //- Spring stiffness coefficient interpolation table (Nm/rad + // or Nm/deg, depending on angleFormat) + interpolationTable stiffness_; + + //- Boolean stating whether the angle around the axis needs to + // be converted to degrees before asking the + // interpolationTable for a value + bool convertToDegrees_; + + //- Damping coefficient (Nms/rad) + scalar damping_; + + +public: + + //- Runtime type information + TypeName("tabulatedAxialAngularSpring"); + + + // Constructors + + //- Construct from components + tabulatedAxialAngularSpring + ( + const dictionary& sDoFRBMRDict + ); + + //- Construct and return a clone + virtual autoPtr clone() const + { + return autoPtr + ( + new tabulatedAxialAngularSpring(*this) + ); + } + + + // Destructor + + virtual ~tabulatedAxialAngularSpring(); + + + // Member Functions + + //- Calculate the restraint position, force and moment. + // Global reference frame vectors. + virtual void restrain + ( + const sixDoFRigidBodyMotion& motion, + vector& restraintPosition, + vector& restraintForce, + vector& restraintMoment + ) const; + + //- Update properties from given dictionary + virtual bool read(const dictionary& sDoFRBMRCoeff); +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace solidBodyMotionFunctions +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionState.C b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionState.C similarity index 100% rename from src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionState.C rename to src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionState.C diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionState.H b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionState.H similarity index 100% rename from src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionState.H rename to src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionState.H diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionStateI.H b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionStateI.H similarity index 100% rename from src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionStateI.H rename to src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionStateI.H diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionStateIO.C b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionStateIO.C similarity index 100% rename from src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionStateIO.C rename to src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionStateIO.C diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/uncoupledSixDoFRigidBodyDisplacement/uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField.C b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/uncoupledSixDoFRigidBodyDisplacement/uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField.C new file mode 100644 index 0000000000..d35eefd671 --- /dev/null +++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/uncoupledSixDoFRigidBodyDisplacement/uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField.C @@ -0,0 +1,173 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd. + \\/ 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 2 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, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +\*---------------------------------------------------------------------------*/ + +#include "uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField.H" +#include "pointPatchFields.H" +#include "addToRunTimeSelectionTable.H" +#include "Time.H" +#include "fvMesh.H" +#include "volFields.H" +#include "uniformDimensionedFields.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField:: +uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField +( + const pointPatch& p, + const DimensionedField& iF +) +: + fixedValuePointPatchField(p, iF), + motion_(), + p0_(p.localPoints()), + rhoInf_(1.0) +{} + + +uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField:: +uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField +( + const pointPatch& p, + const DimensionedField& iF, + const dictionary& dict +) +: + fixedValuePointPatchField(p, iF, dict), + motion_(dict), + rhoInf_(readScalar(dict.lookup("rhoInf"))) +{ + if (!dict.found("value")) + { + updateCoeffs(); + } + + if (dict.found("p0")) + { + p0_ = vectorField("p0", dict , p.size()); + } + else + { + p0_ = p.localPoints(); + } +} + + +uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField:: +uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField +( + const uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField& ptf, + const pointPatch& p, + const DimensionedField& iF, + const pointPatchFieldMapper& mapper +) +: + fixedValuePointPatchField(ptf, p, iF, mapper), + motion_(ptf.motion_), + p0_(ptf.p0_, mapper), + rhoInf_(ptf.rhoInf_) +{} + + +uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField:: +uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField +( + const uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField& ptf, + const DimensionedField& iF +) +: + fixedValuePointPatchField(ptf, iF), + motion_(ptf.motion_), + p0_(ptf.p0_), + rhoInf_(ptf.rhoInf_) +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +void uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField::updateCoeffs() +{ + if (this->updated()) + { + return; + } + + const polyMesh& mesh = this->dimensionedInternalField().mesh()(); + const Time& t = mesh.time(); + + motion_.updatePosition(t.deltaTValue()); + + vector gravity = vector::zero; + + if (db().foundObject("g")) + { + uniformDimensionedVectorField g = + db().lookupObject("g"); + + gravity = g.value(); + } + + // Do not modify the accelerations, except with gravity, where available + motion_.updateForce(gravity*motion_.mass(), vector::zero, t.deltaTValue()); + + Field::operator=(motion_.currentPosition(p0_) - p0_); + + fixedValuePointPatchField::updateCoeffs(); +} + + +void uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField::write +( + Ostream& os +) const +{ + pointPatchField::write(os); + motion_.write(os); + os.writeKeyword("rhoInf") + << rhoInf_ << token::END_STATEMENT << nl; + p0_.writeEntry("p0", os); + writeEntry("value", os); +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +makePointPatchTypeField +( + pointPatchVectorField, + uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField +); + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// ************************************************************************* // diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/uncoupledSixDoFRigidBodyDisplacement/uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField.H b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/uncoupledSixDoFRigidBodyDisplacement/uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField.H new file mode 100644 index 0000000000..fb4177d104 --- /dev/null +++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/uncoupledSixDoFRigidBodyDisplacement/uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField.H @@ -0,0 +1,158 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd. + \\/ 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 2 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, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +Class + Foam::uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField + +Description + Foam::uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField + +SourceFiles + uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField.C + +\*---------------------------------------------------------------------------*/ + +#ifndef uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField_H +#define uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField_H + +#include "fixedValuePointPatchField.H" +#include "sixDoFRigidBodyMotion.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField Declaration +\*---------------------------------------------------------------------------*/ + +class uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField +: + public fixedValuePointPatchField +{ + // Private data + + //- Six dof motion object + sixDoFRigidBodyMotion motion_; + + //- Reference positions of points on the patch + pointField p0_; + + //- Reference density required by the forces object for + // incompressible calculations. Retained here to give + // dictionary compatibility with other sixDoF patches. + scalar rhoInf_; + + +public: + + //- Runtime type information + TypeName("uncoupledSixDoFRigidBodyDisplacement"); + + + // Constructors + + //- Construct from patch and internal field + uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField + ( + const pointPatch&, + const DimensionedField& + ); + + //- Construct from patch, internal field and dictionary + uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField + ( + const pointPatch&, + const DimensionedField&, + const dictionary& + ); + + //- Construct by mapping given patchField onto a new patch + uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField + ( + const uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField&, + const pointPatch&, + const DimensionedField&, + const pointPatchFieldMapper& + ); + + //- Construct and return a clone + virtual autoPtr > clone() const + { + return autoPtr > + ( + new uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField + ( + *this + ) + ); + } + + //- Construct as copy setting internal field reference + uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField + ( + const uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField&, + const DimensionedField& + ); + + //- Construct and return a clone setting internal field reference + virtual autoPtr > clone + ( + const DimensionedField& iF + ) const + { + return autoPtr > + ( + new uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField + ( + *this, + iF + ) + ); + } + + + // Member functions + + // Evaluation functions + + //- Update the coefficients associated with the patch field + virtual void updateCoeffs(); + + + //- Write + virtual void write(Ostream&) const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* //