FreeFOAM The Cross-Platform CFD Toolkit
MixedDiffuseSpecular.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 "MixedDiffuseSpecular.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  diffuseFraction_(readScalar(this->coeffDict().lookup("diffuseFraction")))
39 {}
40 
41 
42 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
43 
44 template <class CloudType>
46 {}
47 
48 
49 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
50 
51 template <class CloudType>
53 (
54  const wallPolyPatch& wpp,
55  const label faceId,
56  vector& U,
57  scalar& Ei,
58  label typeId
59 )
60 {
61  label wppIndex = wpp.index();
62 
63  label wppLocalFace = wpp.whichFace(faceId);
64 
65  vector nw = wpp.faceAreas()[wppLocalFace];
66 
67  // Normal unit vector
68  nw /= mag(nw);
69 
70  // Normal velocity magnitude
71  scalar U_dot_nw = U & nw;
72 
73  CloudType& cloud(this->owner());
74 
75  Random& rndGen(cloud.rndGen());
76 
77  if (diffuseFraction_ > rndGen.scalar01())
78  {
79  // Diffuse reflection
80 
81  // Wall tangential velocity (flow direction)
82  vector Ut = U - U_dot_nw*nw;
83 
84  while (mag(Ut) < SMALL)
85  {
86  // If the incident velocity is parallel to the face normal, no
87  // tangential direction can be chosen. Add a perturbation to the
88  // incoming velocity and recalculate.
89 
90  U = vector
91  (
92  U.x()*(0.8 + 0.2*rndGen.scalar01()),
93  U.y()*(0.8 + 0.2*rndGen.scalar01()),
94  U.z()*(0.8 + 0.2*rndGen.scalar01())
95  );
96 
97  U_dot_nw = U & nw;
98 
99  Ut = U - U_dot_nw*nw;
100  }
101 
102  // Wall tangential unit vector
103  vector tw1 = Ut/mag(Ut);
104 
105  // Other tangential unit vector
106  vector tw2 = nw^tw1;
107 
108  scalar T = cloud.boundaryT().boundaryField()[wppIndex][wppLocalFace];
109 
110  scalar mass = cloud.constProps(typeId).mass();
111 
112  scalar iDof = cloud.constProps(typeId).internalDegreesOfFreedom();
113 
114  U =
115  sqrt(CloudType::kb*T/mass)
116  *(
117  rndGen.GaussNormal()*tw1
118  + rndGen.GaussNormal()*tw2
119  - sqrt(-2.0*log(max(1 - rndGen.scalar01(), VSMALL)))*nw
120  );
121 
122  U += cloud.boundaryU().boundaryField()[wppIndex][wppLocalFace];
123 
124  Ei = cloud.equipartitionInternalEnergy(T, iDof);
125  }
126  else
127  {
128  // Specular reflection
129 
130  if (U_dot_nw > 0.0)
131  {
132  U -= 2.0*U_dot_nw*nw;
133  }
134  }
135 
136 }
137 
138 
139 // ************************************************************************* //