FreeFOAM The Cross-Platform CFD Toolkit
reflectParcel.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) 1991-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 "reflectParcel.H"
28 #include <OpenFOAM/wallPolyPatch.H>
29 
30 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
31 
32 namespace Foam
33 {
34 
35 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
36 
38 
40 (
41  wallModel,
44 );
45 
46 
47 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
48 
49 // Construct from components
51 (
52  const dictionary& dict,
53  const volVectorField& U,
54  spray& sm
55 )
56 :
57  wallModel(dict, U, sm),
58  U_(U),
59  coeffsDict_(dict.subDict(typeName + "Coeffs")),
60  elasticity_(readScalar(coeffsDict_.lookup("elasticity")))
61 {}
62 
63 
64 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
65 
67 {}
68 
69 
70 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
71 
72 // Return 'keepParcel'
74 (
75  parcel& p,
76  const label globalFacei
77 ) const
78 {
79  label patchi = p.patch(globalFacei);
80  label facei = p.patchFace(patchi, globalFacei);
81 
82  const polyMesh& mesh = spray_.mesh();
83 
84  if (isA<wallPolyPatch>(mesh_.boundaryMesh()[patchi]))
85  {
86  // wallNormal defined to point outwards of domain
87  vector Sf = mesh_.Sf().boundaryField()[patchi][facei];
88  Sf /= mag(Sf);
89 
90  if (!mesh.moving())
91  {
92  // static mesh
93  scalar Un = p.U() & Sf;
94 
95  if (Un > 0)
96  {
97  p.U() -= (1.0 + elasticity_)*Un*Sf;
98  }
99 
100  }
101  else
102  {
103  // moving mesh
104  vector Ub1 = U_.boundaryField()[patchi][facei];
105  vector Ub0 = U_.oldTime().boundaryField()[patchi][facei];
106 
107  scalar dt = spray_.runTime().deltaT().value();
108  const vectorField& oldPoints = mesh.oldPoints();
109 
110  const vector& Cf1 = mesh.faceCentres()[globalFacei];
111 
112  vector Cf0 = mesh.faces()[globalFacei].centre(oldPoints);
113  vector Cf = Cf0 + p.stepFraction()*(Cf1 - Cf0);
114  vector Sf0 = mesh.faces()[globalFacei].normal(oldPoints);
115 
116  // for layer addition Sf0 = vector::zero and we use Sf
117  if (mag(Sf0) > SMALL)
118  {
119  Sf0 /= mag(Sf0);
120  }
121  else
122  {
123  Sf0 = Sf;
124  }
125 
126  scalar magSfDiff = mag(Sf - Sf0);
127 
128  vector Ub = Ub0 + p.stepFraction()*(Ub1 - Ub0);
129 
130  if (magSfDiff > SMALL)
131  {
132  // rotation + translation
133  vector Sfp = Sf0 + p.stepFraction()*(Sf - Sf0);
134 
135  vector omega = Sf0 ^ Sf;
136  scalar magOmega = mag(omega);
137  omega /= magOmega+SMALL;
138 
139  scalar phiVel = ::asin(magOmega)/dt;
140 
141  scalar dist = (p.position() - Cf) & Sfp;
142  vector pos = p.position() - dist*Sfp;
143  vector vrot = phiVel*(omega ^ (pos - Cf));
144 
145  vector v = Ub + vrot;
146 
147  scalar Un = ((p.U() - v) & Sfp);
148 
149  if (Un > 0.0)
150  {
151  p.U() -= (1.0 + elasticity_)*Un*Sfp;
152  }
153  }
154  else
155  {
156  // translation
157  vector Ur = p.U() - Ub;
158  scalar Urn = Ur & Sf;
159  /*
160  if (mag(Ub-v) > SMALL)
161  {
162  Info << "reflectParcel:: v = " << v
163  << ", Ub = " << Ub
164  << ", facei = " << facei
165  << ", patchi = " << patchi
166  << ", globalFacei = " << globalFacei
167  << ", name = " << mesh_.boundaryMesh()[patchi].name()
168  << endl;
169  }
170  */
171  if (Urn > 0.0)
172  {
173  p.U() -= (1.0 + elasticity_)*Urn*Sf;
174  }
175  }
176  }
177 
178  }
179  else
180  {
181  FatalError
182  << "bool reflectParcel::wallTreatment(parcel& parcel) const "
183  << " parcel has hit a boundary "
184  << mesh_.boundary()[patchi].type()
185  << " which not yet has been implemented."
186  << abort(FatalError);
187  }
188  return true;
189 }
190 
191 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
192 
193 } // End namespace Foam
194 
195 // ************************ vim: set sw=4 sts=4 et: ************************ //