FreeFOAM The Cross-Platform CFD Toolkit
Hosted by SourceForge:
Get FreeFOAM at SourceForge.net.
            Fast, secure and Free Open Source software downloads

molecule.C

Go to the documentation of this file.
00001 /*---------------------------------------------------------------------------*\
00002   =========                 |
00003   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
00004    \\    /   O peration     |
00005     \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
00006      \\/     M anipulation  |
00007 -------------------------------------------------------------------------------
00008 License
00009     This file is part of OpenFOAM.
00010 
00011     OpenFOAM is free software: you can redistribute it and/or modify it
00012     under the terms of the GNU General Public License as published by
00013     the Free Software Foundation, either version 3 of the License, or
00014     (at your option) any later version.
00015 
00016     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
00017     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00018     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
00019     for more details.
00020 
00021     You should have received a copy of the GNU General Public License
00022     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
00023 
00024 \*----------------------------------------------------------------------------*/
00025 
00026 #include <molecule/moleculeCloud.H>
00027 #include "molecule.H"
00028 #include <OpenFOAM/Random.H>
00029 #include <OpenFOAM/Time.H>
00030 
00031 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
00032 
00033 Foam::tensor Foam::molecule::rotationTensorX(scalar phi) const
00034 {
00035     return tensor
00036     (
00037         1, 0, 0,
00038         0, Foam::cos(phi), -Foam::sin(phi),
00039         0, Foam::sin(phi), Foam::cos(phi)
00040     );
00041 }
00042 
00043 
00044 Foam::tensor Foam::molecule::rotationTensorY(scalar phi) const
00045 {
00046     return tensor
00047     (
00048         Foam::cos(phi), 0, Foam::sin(phi),
00049         0, 1, 0,
00050         -Foam::sin(phi), 0, Foam::cos(phi)
00051     );
00052 }
00053 
00054 
00055 Foam::tensor Foam::molecule::rotationTensorZ(scalar phi) const
00056 {
00057     return tensor
00058     (
00059         Foam::cos(phi), -Foam::sin(phi), 0,
00060         Foam::sin(phi), Foam::cos(phi), 0,
00061         0, 0, 1
00062     );
00063 }
00064 
00065 
00066 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
00067 
00068 Foam::molecule::trackData::trackData
00069 (
00070     moleculeCloud& molCloud,
00071     label part
00072 )
00073 :
00074     Particle<molecule>::trackData(molCloud),
00075     molCloud_(molCloud),
00076     part_(part)
00077 {}
00078 
00079 
00080 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
00081 
00082 
00083 bool Foam::molecule::move(molecule::trackData& td)
00084 {
00085     td.switchProcessor = false;
00086     td.keepParticle = true;
00087 
00088     const constantProperties& constProps(td.molCloud().constProps(id_));
00089 
00090     scalar deltaT = cloud().pMesh().time().deltaT().value();
00091 
00092     if (td.part() == 0)
00093     {
00094         // First leapfrog velocity adjust part, required before tracking+force
00095         // part
00096 
00097         v_ += 0.5*deltaT*a_;
00098 
00099         pi_ += 0.5*deltaT*tau_;
00100     }
00101     else if (td.part() == 1)
00102     {
00103         // Leapfrog tracking part
00104 
00105         scalar tEnd = (1.0 - stepFraction())*deltaT;
00106         scalar dtMax = tEnd;
00107 
00108         while (td.keepParticle && !td.switchProcessor && tEnd > ROOTVSMALL)
00109         {
00110             // set the lagrangian time-step
00111             scalar dt = min(dtMax, tEnd);
00112 
00113             dt *= trackToFace(position() + dt*v_, td);
00114 
00115             tEnd -= dt;
00116             stepFraction() = 1.0 - tEnd/deltaT;
00117         }
00118     }
00119     else if (td.part() == 2)
00120     {
00121         // Leapfrog orientation adjustment, carried out before force calculation
00122         // but after tracking stage, i.e. rotation carried once linear motion
00123         // complete.
00124 
00125         if (!constProps.pointMolecule())
00126         {
00127             const diagTensor& momentOfInertia(constProps.momentOfInertia());
00128 
00129             tensor R;
00130 
00131             if (!constProps.linearMolecule())
00132             {
00133                 R = rotationTensorX(0.5*deltaT*pi_.x()/momentOfInertia.xx());
00134                 pi_ = pi_ & R;
00135                 Q_ = Q_ & R;
00136             }
00137 
00138             R = rotationTensorY(0.5*deltaT*pi_.y()/momentOfInertia.yy());
00139             pi_ = pi_ & R;
00140             Q_ = Q_ & R;
00141 
00142             R = rotationTensorZ(deltaT*pi_.z()/momentOfInertia.zz());
00143             pi_ = pi_ & R;
00144             Q_ = Q_ & R;
00145 
00146             R = rotationTensorY(0.5*deltaT*pi_.y()/momentOfInertia.yy());
00147             pi_ = pi_ & R;
00148             Q_ = Q_ & R;
00149 
00150             if (!constProps.linearMolecule())
00151             {
00152                 R = rotationTensorX(0.5*deltaT*pi_.x()/momentOfInertia.xx());
00153                 pi_ = pi_ & R;
00154                 Q_ = Q_ & R;
00155             }
00156         }
00157 
00158         setSitePositions(constProps);
00159     }
00160     else if (td.part() == 3)
00161     {
00162         // Second leapfrog velocity adjust part, required after tracking+force
00163         // part
00164 
00165         scalar m = constProps.mass();
00166 
00167         a_ = vector::zero;
00168 
00169         tau_ = vector::zero;
00170 
00171         forAll(siteForces_, s)
00172         {
00173             const vector& f = siteForces_[s];
00174 
00175             a_ += f/m;
00176 
00177             tau_ += (constProps.siteReferencePositions()[s] ^ (Q_.T() & f));
00178         }
00179 
00180         v_ += 0.5*deltaT*a_;
00181 
00182         pi_ += 0.5*deltaT*tau_;
00183 
00184         if (constProps.pointMolecule())
00185         {
00186             tau_ = vector::zero;
00187 
00188             pi_ = vector::zero;
00189         }
00190 
00191         if (constProps.linearMolecule())
00192         {
00193             tau_.x() = 0.0;
00194 
00195             pi_.x() = 0.0;
00196         }
00197     }
00198     else
00199     {
00200         FatalErrorIn("molecule::move(molecule::trackData& td)") << nl
00201             << td.part()
00202             << " is an invalid part of the integration method."
00203             << abort(FatalError);
00204     }
00205 
00206     return td.keepParticle;
00207 }
00208 
00209 
00210 void Foam::molecule::transformProperties(const tensor& T)
00211 {
00212     Q_ = T & Q_;
00213 
00214     sitePositions_ = position_ + (T & (sitePositions_ - position_));
00215 
00216     siteForces_ = T & siteForces_;
00217 }
00218 
00219 
00220 void Foam::molecule::transformProperties(const vector& separation)
00221 {
00222     if (special_ == SPECIAL_TETHERED)
00223     {
00224         specialPosition_ += separation;
00225     }
00226 }
00227 
00228 
00229 void Foam::molecule::setSitePositions(const constantProperties& constProps)
00230 {
00231     sitePositions_ = position_ + (Q_ & constProps.siteReferencePositions());
00232 }
00233 
00234 
00235 void Foam::molecule::setSiteSizes(label size)
00236 {
00237     sitePositions_.setSize(size);
00238 
00239     siteForces_.setSize(size);
00240 }
00241 
00242 
00243 bool Foam::molecule::hitPatch
00244 (
00245     const polyPatch&,
00246     molecule::trackData&,
00247     const label
00248 )
00249 {
00250     return false;
00251 }
00252 
00253 
00254 bool Foam::molecule::hitPatch
00255 (
00256     const polyPatch&,
00257     int&,
00258     const label
00259 )
00260 {
00261     return false;
00262 }
00263 
00264 
00265 void Foam::molecule::hitProcessorPatch
00266 (
00267     const processorPolyPatch&,
00268     molecule::trackData& td
00269 )
00270 {
00271     td.switchProcessor = true;
00272 }
00273 
00274 
00275 void Foam::molecule::hitProcessorPatch
00276 (
00277     const processorPolyPatch&,
00278     int&
00279 )
00280 {}
00281 
00282 
00283 void Foam::molecule::hitWallPatch
00284 (
00285     const wallPolyPatch& wpp,
00286     molecule::trackData& td
00287 )
00288 {
00289     vector nw = wpp.faceAreas()[wpp.whichFace(face())];
00290     nw /= mag(nw);
00291 
00292     scalar vn = v_ & nw;
00293 
00294     // Specular reflection
00295     if (vn > 0)
00296     {
00297         v_ -= 2*vn*nw;
00298     }
00299 }
00300 
00301 
00302 void Foam::molecule::hitWallPatch
00303 (
00304     const wallPolyPatch&,
00305     int&
00306 )
00307 {}
00308 
00309 
00310 void Foam::molecule::hitPatch
00311 (
00312     const polyPatch&,
00313     molecule::trackData& td
00314 )
00315 {
00316     td.keepParticle = false;
00317 }
00318 
00319 
00320 void Foam::molecule::hitPatch
00321 (
00322     const polyPatch&,
00323     int&
00324 )
00325 {}
00326 
00327 
00328 // ************************ vim: set sw=4 sts=4 et: ************************ //
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines