TrajectoryCollision.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
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
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
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
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
26 #include "TrajectoryCollision.H"
27 
28 // * * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * //
29 
30 template<class CloudType>
32 {
34 }
35 
36 
37 template<class CloudType>
39 (
40  const scalar dt,
41  parcelType& p1,
42  parcelType& p2,
43  scalar& m1,
44  scalar& m2
45 )
46 {
47  bool coalescence = false;
48 
49  const vector& pos1 = p1.position();
50  const vector& pos2 = p2.position();
51 
52  const vector& U1 = p1.U();
53  const vector& U2 = p2.U();
54 
55  vector URel = U1 - U2;
56 
57  vector d = pos2 - pos1;
58  scalar magd = mag(d);
59 
60  scalar vAlign = URel & (d/(magd + ROOTVSMALL));
61 
62  if (vAlign > 0)
63  {
64  const scalar d1 = p1.d();
65  const scalar d2 = p2.d();
66 
67  scalar sumD = d1 + d2;
68 
69  if (vAlign*dt > magd - 0.5*sumD)
70  {
71  scalar magU1 = mag(U1) + ROOTVSMALL;
72  scalar magU2 = mag(U2) + ROOTVSMALL;
73  vector n1 = U1/magU1;
74  vector n2 = U2/magU2;
75 
76  scalar n1n2 = n1 & n2;
77  scalar n1d = n1 & d;
78  scalar n2d = n2 & d;
79 
80  scalar det = 1.0 - sqr(n1n2);
81 
82  scalar alpha = GREAT;
83  scalar beta = GREAT;
84 
85  if (mag(det) > 1.0e-4)
86  {
87  beta = -(n2d - n1n2*n1d)/det;
88  alpha = n1d + n1n2*beta;
89  }
90 
91  alpha /= magU1*dt;
92  beta /= magU2*dt;
93 
94  // is collision possible within this timestep
95  if ((alpha > 0) && (alpha < 1.0) && (beta > 0) && (beta < 1.0))
96  {
97  vector p1c = pos1 + alpha*U1*dt;
98  vector p2c = pos2 + beta*U2*dt;
99 
100  scalar closestDist = mag(p1c - p2c);
101 
102  scalar collProb =
103  pow(0.5*sumD/max(0.5*sumD, closestDist), cSpace_)
104  *exp(-cTime_*mag(alpha - beta));
105 
106  scalar xx = this->owner().rndGen().template sample01<scalar>();
107 
108  // collision occurs
109  if (xx > collProb)
110  {
111  if (d1 > d2)
112  {
113  coalescence = this->collideSorted(dt, p1, p2, m1, m2);
114  }
115  else
116  {
117  coalescence = this->collideSorted(dt, p2, p1, m2, m1);
118  }
119  }
120  }
121  }
122  }
123 
124  return coalescence;
125 }
126 
127 
128 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
129 
130 template<class CloudType>
132 (
133  const dictionary& dict,
134  CloudType& owner
135 )
136 :
137  ORourkeCollision<CloudType>(dict, owner, typeName),
138  cSpace_(readScalar(this->coeffDict().lookup("cSpace"))),
139  cTime_(readScalar(this->coeffDict().lookup("cTime")))
140 {}
141 
142 
143 template<class CloudType>
145 (
147 )
148 :
150  cSpace_(cm.cSpace_),
151  cTime_(cm.cTime_)
152 {}
153 
154 
155 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
156 
157 template<class CloudType>
159 {}
160 
161 
162 // ************************************************************************* //
volVectorField & U1
Definition: createFields.H:18
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
#define readScalar
Definition: doubleScalar.C:38
dimensioned< scalar > mag(const dimensioned< Type > &)
virtual void collide(const scalar dt)
Main collision routine.
scalar cSpace_
Space coefficient.
CloudType::parcelType parcelType
Convenience typedef to the cloud&#39;s parcel type.
volVectorField & U2
Definition: createFields.H:23
dimensionedScalar beta("beta", dimless/dimTemperature, laminarTransport)
dimensionedScalar exp(const dimensionedScalar &ds)
scalar cTime_
Time coefficient.
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:137
Trajectory collision model by N. Nordin, based on O&#39;Rourke&#39;s collision model.
dictionary dict
stressControl lookup("compactNormalStress") >> compactNormalStress
dimensionedScalar det(const dimensionedSphericalTensor &dt)
const dimensionedScalar e
Elementary charge.
Definition: doubleFloat.H:78
const dimensionedScalar alpha
Fine-structure constant: default SI units: [].
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
virtual bool collideParcels(const scalar dt, parcelType &p1, parcelType &p2, scalar &m1, scalar &m2)
Collide parcels and return true if mass has changed.
Collision model by P.J. O&#39;Rourke.
TrajectoryCollision(const dictionary &dict, CloudType &cloud)
Construct from dictionary.
virtual ~TrajectoryCollision()
Destructor.
Templated base class for dsmc cloud.
Definition: DSMCCloud.H:68
dimensionedSymmTensor sqr(const dimensionedVector &dv)