1 /*---------------------------------------------------------------------------*\
3 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
5 \\ / A nd | Copyright held by original author
7 -------------------------------------------------------------------------------
9 This file is part of OpenFOAM.
11 OpenFOAM is free software; you can redistribute it and/or modify it
12 under the terms of the GNU General Public License as published by the
13 Free Software Foundation; either version 2 of the License, or (at your
14 option) any later version.
16 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
21 You should have received a copy of the GNU General Public License
22 along with OpenFOAM; if not, write to the Free Software Foundation,
23 Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
25 \*---------------------------------------------------------------------------*/
27 #include "harmonicOscillation.H"
28 #include "addToRunTimeSelectionTable.H"
29 #include "mathematicalConstants.H"
31 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
35 namespace solidBodyMotionFunctions
38 defineTypeNameAndDebug(harmonicOscillation, 0);
39 addToRunTimeSelectionTable
41 solidBodyMotionFunction,
49 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
52 Foam::solidBodyMotionFunctions::harmonicOscillation::calcTransformation
57 const vector translation =
63 sin(transAngularFreq_.x()*t + transPhaseAngles_.x()),
64 sin(transAngularFreq_.y()*t + transPhaseAngles_.y()),
65 sin(transAngularFreq_.z()*t + transPhaseAngles_.z())
70 const vector rotation =
76 sin(rotAngularFreq_.x()*t + rotPhaseAngles_.x()),
77 sin(rotAngularFreq_.y()*t + rotPhaseAngles_.y()),
78 sin(rotAngularFreq_.z()*t + rotPhaseAngles_.z())
82 const quaternion R(rotation.x(), rotation.y(), rotation.z());
85 septernion(origin_ + translation)*R*septernion(-origin_)
91 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
93 Foam::solidBodyMotionFunctions::harmonicOscillation::harmonicOscillation
95 const dictionary& SBMFCoeffs,
99 solidBodyMotionFunction(SBMFCoeffs, runTime)
104 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
107 Foam::solidBodyMotionFunctions::harmonicOscillation::transformation() const
109 const scalar t = time_.value();
111 const septernion TR = calcTransformation(t);
113 Info<< "solidBodyMotionFunctions::harmonicOscillation::transformation(): "
114 << "Time = " << t << " transformation: " << TR << endl;
121 Foam::solidBodyMotionFunctions::harmonicOscillation::velocity() const
123 const scalar t = time_.value();
124 const scalar dt = time_.deltaT().value();
126 const septernion velocity
128 (calcTransformation(t).t() - calcTransformation(t - dt).t())/dt,
129 (calcTransformation(t).r()/calcTransformation(t - dt).r())/dt
136 bool Foam::solidBodyMotionFunctions::harmonicOscillation::read
138 const dictionary& SBMFCoeffs
141 solidBodyMotionFunction::read(SBMFCoeffs);
143 SBMFCoeffs_.lookup("origin") >> origin_;
144 SBMFCoeffs_.lookup("translationalAmplitudes") >> transAmplitudes_;
145 SBMFCoeffs_.lookup("translationalAngularFrequencies") >> transAngularFreq_;
146 SBMFCoeffs_.lookup("translationalPhaseAngles") >> transPhaseAngles_;
147 SBMFCoeffs_.lookup("rotationalAmplitudes") >> rotAmplitudes_;
148 SBMFCoeffs_.lookup("rotationalAngularFrequencies") >> rotAngularFreq_;
149 SBMFCoeffs_.lookup("rotationalPhaseAngles") >> rotPhaseAngles_;
150 SBMFCoeffs_.lookup("inDegrees") >> inDegrees_;
152 // Sanity check for negative frequencies
155 transAngularFreq_.x() < -SMALL
156 || transAngularFreq_.y() < -SMALL
157 || transAngularFreq_.z() < -SMALL
158 || rotAngularFreq_.x() < -SMALL
159 || rotAngularFreq_.y() < -SMALL
160 || rotAngularFreq_.z() < -SMALL
165 "harmonicOscillation::read"
167 " const fvMesh& mesh\n"
168 " const dictionary& dict\n"
170 ) << "Negative angular frequency detected."
171 << abort(FatalError);
174 // Convert to radians if necessary, printing data
177 const scalar piBy180 = mathematicalConstant::pi/180.0;
179 transPhaseAngles_ *= piBy180;
180 rotAmplitudes_ *= piBy180;
181 rotPhaseAngles_ *= piBy180;
183 Info<< "Data in degrees converted:" << nl
184 << "translationalPhaseAngles = " << transPhaseAngles_ << nl
185 << "rotationalAmplitudes = " << rotAmplitudes_ << nl
186 << "rotationalPhaseAngles = " << rotPhaseAngles_ << endl;
189 // Count prescribed rotations: only two harmonic rotations can be prescribed
193 mag(rotAmplitudes_.x()) > SMALL ? ++nRot : /* null expression */ 0 ;
194 mag(rotAmplitudes_.y()) > SMALL ? ++nRot : /* null expression */ 0 ;
195 mag(rotAmplitudes_.z()) > SMALL ? ++nRot : /* null expression */ 0 ;
197 // Check if more than two rotations prescribed
202 "harmonicOscillation::read"
204 " const fvMesh& mesh\n"
205 " const dictionary& dict\n"
207 ) << "Only two harmonic rotations can be prescribed. Prescribing "
208 << "more than two yields unphysical results. "
209 << abort(FatalError);
216 // ************************************************************************* //