Go to the documentation of this file.00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026 #include <molecule/moleculeCloud.H>
00027 #include "molecule.H"
00028 #include <OpenFOAM/Random.H>
00029 #include <OpenFOAM/Time.H>
00030
00031
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
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
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
00095
00096
00097 v_ += 0.5*deltaT*a_;
00098
00099 pi_ += 0.5*deltaT*tau_;
00100 }
00101 else if (td.part() == 1)
00102 {
00103
00104
00105 scalar tEnd = (1.0 - stepFraction())*deltaT;
00106 scalar dtMax = tEnd;
00107
00108 while (td.keepParticle && !td.switchProcessor && tEnd > ROOTVSMALL)
00109 {
00110
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
00122
00123
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
00163
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
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