FreeFOAM The Cross-Platform CFD Toolkit
MaxwellianThermal.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) 2009-2010 OpenCFD Ltd.
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 "MaxwellianThermal.H"
27 
28 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
29 
30 template <class CloudType>
32 (
33  const dictionary& dict,
34  CloudType& cloud
35 )
36 :
37  WallInteractionModel<CloudType>(dict, cloud, typeName)
38 {}
39 
40 
41 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
42 
43 template <class CloudType>
45 {}
46 
47 
48 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
49 
50 template <class CloudType>
52 (
53  const wallPolyPatch& wpp,
54  const label faceId,
55  vector& U,
56  scalar& Ei,
57  label typeId
58 )
59 {
60  label wppIndex = wpp.index();
61 
62  label wppLocalFace = wpp.whichFace(faceId);
63 
64  vector nw = wpp.faceAreas()[wppLocalFace];
65 
66  // Normal unit vector
67  nw /= mag(nw);
68 
69  // Normal velocity magnitude
70  scalar magUn = U & nw;
71 
72  // Wall tangential velocity (flow direction)
73  vector Ut = U - magUn*nw;
74 
75  CloudType& cloud(this->owner());
76 
77  Random& rndGen(cloud.rndGen());
78 
79  while (mag(Ut) < SMALL)
80  {
81  // If the incident velocity is parallel to the face normal, no
82  // tangential direction can be chosen. Add a perturbation to the
83  // incoming velocity and recalculate.
84 
85  U = vector
86  (
87  U.x()*(0.8 + 0.2*rndGen.scalar01()),
88  U.y()*(0.8 + 0.2*rndGen.scalar01()),
89  U.z()*(0.8 + 0.2*rndGen.scalar01())
90  );
91 
92  magUn = U & nw;
93 
94  Ut = U - magUn*nw;
95  }
96 
97  // Wall tangential unit vector
98  vector tw1 = Ut/mag(Ut);
99 
100  // Other tangential unit vector
101  vector tw2 = nw^tw1;
102 
103  scalar T = cloud.boundaryT().boundaryField()[wppIndex][wppLocalFace];
104 
105  scalar mass = cloud.constProps(typeId).mass();
106 
107  scalar iDof = cloud.constProps(typeId).internalDegreesOfFreedom();
108 
109  U =
110  sqrt(CloudType::kb*T/mass)
111  *(
112  rndGen.GaussNormal()*tw1
113  + rndGen.GaussNormal()*tw2
114  - sqrt(-2.0*log(max(1 - rndGen.scalar01(), VSMALL)))*nw
115  );
116 
117  U += cloud.boundaryU().boundaryField()[wppIndex][wppLocalFace];
118 
119  Ei = cloud.equipartitionInternalEnergy(T, iDof);
120 }
121 
122 
123 // ************************ vim: set sw=4 sts=4 et: ************************ //