TrajectoryCollision.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2020 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 (
33  typename CloudType::parcelType::trackingData& td,
34  const scalar dt
35 )
36 {
38 }
39 
40 
41 template<class CloudType>
43 (
44  const scalar dt,
45  parcelType& p1,
46  parcelType& p2,
47  scalar& m1,
48  scalar& m2
49 )
50 {
51  bool coalescence = false;
52 
53  const vector& pos1 = p1.position();
54  const vector& pos2 = p2.position();
55 
56  const vector& U1 = p1.U();
57  const vector& U2 = p2.U();
58 
59  vector URel = U1 - U2;
60 
61  vector d = pos2 - pos1;
62  scalar magd = mag(d);
63 
64  scalar vAlign = URel & (d/(magd + rootVSmall));
65 
66  if (vAlign > 0)
67  {
68  const scalar d1 = p1.d();
69  const scalar d2 = p2.d();
70 
71  scalar sumD = d1 + d2;
72 
73  if (vAlign*dt > magd - 0.5*sumD)
74  {
75  scalar magU1 = mag(U1) + rootVSmall;
76  scalar magU2 = mag(U2) + rootVSmall;
77  vector n1 = U1/magU1;
78  vector n2 = U2/magU2;
79 
80  scalar n1n2 = n1 & n2;
81  scalar n1d = n1 & d;
82  scalar n2d = n2 & d;
83 
84  scalar det = 1.0 - sqr(n1n2);
85 
86  scalar alpha = great;
87  scalar beta = great;
88 
89  if (mag(det) > 1.0e-4)
90  {
91  beta = -(n2d - n1n2*n1d)/det;
92  alpha = n1d + n1n2*beta;
93  }
94 
95  alpha /= magU1*dt;
96  beta /= magU2*dt;
97 
98  // is collision possible within this timestep
99  if ((alpha > 0) && (alpha < 1.0) && (beta > 0) && (beta < 1.0))
100  {
101  vector p1c = pos1 + alpha*U1*dt;
102  vector p2c = pos2 + beta*U2*dt;
103 
104  scalar closestDist = mag(p1c - p2c);
105 
106  scalar collProb =
107  pow(0.5*sumD/max(0.5*sumD, closestDist), cSpace_)
108  *exp(-cTime_*mag(alpha - beta));
109 
110  scalar xx = this->owner().rndGen().template sample01<scalar>();
111 
112  // collision occurs
113  if (xx > collProb)
114  {
115  if (d1 > d2)
116  {
117  coalescence = this->collideSorted(dt, p1, p2, m1, m2);
118  }
119  else
120  {
121  coalescence = this->collideSorted(dt, p2, p1, m2, m1);
122  }
123  }
124  }
125  }
126  }
127 
128  return coalescence;
129 }
130 
131 
132 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
133 
134 template<class CloudType>
136 (
137  const dictionary& dict,
138  CloudType& owner
139 )
140 :
141  ORourkeCollision<CloudType>(dict, owner, typeName),
142  cSpace_(this->coeffDict().template lookup<scalar>("cSpace")),
143  cTime_(this->coeffDict().template lookup<scalar>("cTime"))
144 {}
145 
146 
147 template<class CloudType>
149 (
151 )
152 :
154  cSpace_(cm.cSpace_),
155  cTime_(cm.cTime_)
156 {}
157 
158 
159 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
160 
161 template<class CloudType>
163 {}
164 
165 
166 // ************************************************************************* //
virtual ~TrajectoryCollision()
Destructor.
dictionary dict
scalar cSpace_
Space coefficient.
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:156
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
dimensionedSymmTensor sqr(const dimensionedVector &dv)
volScalarField alpha(IOobject("alpha", runTime.timeName(), mesh, IOobject::READ_IF_PRESENT, IOobject::AUTO_WRITE), lambda *max(Ua &U, zeroSensitivity))
dimensionedScalar det(const dimensionedSphericalTensor &dt)
dimensionedScalar exp(const dimensionedScalar &ds)
virtual bool collideParcels(const scalar dt, parcelType &p1, parcelType &p2, scalar &m1, scalar &m2)
Collide parcels and return true if mass has changed.
CloudType::parcelType parcelType
Convenience typedef to the cloud&#39;s parcel type.
Trajectory collision model by N. Nordin, based on O&#39;Rourke&#39;s collision model.
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
Collision model by P.J. O&#39;Rourke.
virtual void collide(typename CloudType::parcelType::trackingData &td, const scalar dt)
Main collision routine.
dimensioned< scalar > mag(const dimensioned< Type > &)
scalar cTime_
Time coefficient.
Templated base class for dsmc cloud.
Definition: DSMCCloud.H:75
const dimensionedScalar e
Elementary charge.
Definition: doubleScalar.H:105
TrajectoryCollision(const dictionary &dict, CloudType &cloud)
Construct from dictionary.