FreeFOAM The Cross-Platform CFD Toolkit
dsmcFields.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 "dsmcFields.H"
27 #include <finiteVolume/volFields.H>
28 #include <OpenFOAM/dictionary.H>
29 #include <dsmc/dsmcCloud.H>
30 
31 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
32 
33 namespace Foam
34 {
35  defineTypeNameAndDebug(dsmcFields, 0);
36 }
37 
38 
39 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
40 
41 Foam::dsmcFields::dsmcFields
42 (
43  const word& name,
44  const objectRegistry& obr,
45  const dictionary& dict,
46  const bool loadFromFiles
47 )
48 :
49  name_(name),
50  obr_(obr),
51  active_(true)
52 {
53  // Check if the available mesh is an fvMesh, otherwise deactivate
54  if (!isA<fvMesh>(obr_))
55  {
56  active_ = false;
57  WarningIn
58  (
59  "dsmcFields::dsmcFields"
60  "(const objectRegistry&, const dictionary&)"
61  ) << "No fvMesh available, deactivating." << nl
62  << endl;
63  }
64 
65  read(dict);
66 }
67 
68 
69 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
70 
72 {}
73 
74 
75 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
76 
78 {
79  if (active_)
80  {
81 
82  }
83 }
84 
85 
87 {
88  // Do nothing - only valid on write
89 }
90 
91 
93 {
94  // Do nothing - only valid on write
95 }
96 
97 
99 {
100  if (active_)
101  {
102  word rhoNMeanName = "rhoNMean";
103  word rhoMMeanName = "rhoMMean";
104  word momentumMeanName = "momentumMean";
105  word linearKEMeanName = "linearKEMean";
106  word internalEMeanName = "internalEMean";
107  word iDofMeanName = "iDofMean";
108  word fDMeanName = "fDMean";
109 
110  const volScalarField& rhoNMean = obr_.lookupObject<volScalarField>
111  (
112  rhoNMeanName
113  );
114 
115  const volScalarField& rhoMMean = obr_.lookupObject<volScalarField>
116  (
117  rhoMMeanName
118  );
119 
120  const volVectorField& momentumMean = obr_.lookupObject<volVectorField>
121  (
122  momentumMeanName
123  );
124 
125  const volScalarField& linearKEMean = obr_.lookupObject<volScalarField>
126  (
127  linearKEMeanName
128  );
129 
130  const volScalarField& internalEMean = obr_.lookupObject<volScalarField>
131  (
132  internalEMeanName
133  );
134 
135  const volScalarField& iDofMean = obr_.lookupObject<volScalarField>
136  (
137  iDofMeanName
138  );
139 
140  volVectorField fDMean = obr_.lookupObject<volVectorField>
141  (
142  fDMeanName
143  );
144 
145  if (min(mag(rhoNMean)).value() > VSMALL)
146  {
147  Info<< "Calculating dsmcFields." << endl;
148 
149  Info<< " Calculating UMean field." << endl;
151  (
152  IOobject
153  (
154  "UMean",
155  obr_.time().timeName(),
156  obr_,
158  ),
159  momentumMean/rhoMMean
160  );
161 
162  Info<< " Calculating translationalT field." << endl;
163  volScalarField translationalT
164  (
165  IOobject
166  (
167  "translationalT",
168  obr_.time().timeName(),
169  obr_,
171  ),
172  2.0/(3.0*dsmcCloud::kb*rhoNMean)
173  *(linearKEMean - 0.5*rhoMMean*(UMean & UMean))
174  );
175 
176  Info<< " Calculating internalT field." << endl;
177  volScalarField internalT
178  (
179  IOobject
180  (
181  "internalT",
182  obr_.time().timeName(),
183  obr_,
185  ),
186  (2.0/dsmcCloud::kb)*(internalEMean/iDofMean)
187  );
188 
189  Info<< " Calculating overallT field." << endl;
190  volScalarField overallT
191  (
192  IOobject
193  (
194  "overallT",
195  obr_.time().timeName(),
196  obr_,
198  ),
199  2.0/(dsmcCloud::kb*(3.0*rhoNMean + iDofMean))
200  *(linearKEMean - 0.5*rhoMMean*(UMean & UMean) + internalEMean)
201  );
202 
203  Info<< " Calculating pressure field." << endl;
205  (
206  IOobject
207  (
208  "p",
209  obr_.time().timeName(),
210  obr_,
212  ),
213  dsmcCloud::kb*rhoNMean*translationalT
214  );
215 
216  const fvMesh& mesh = fDMean.mesh();
217 
218  forAll(mesh.boundaryMesh(), i)
219  {
220  const polyPatch& patch = mesh.boundaryMesh()[i];
221 
222  if (isA<wallPolyPatch>(patch))
223  {
224  p.boundaryField()[i] =
225  fDMean.boundaryField()[i]
226  & (patch.faceAreas()/mag(patch.faceAreas()));
227  }
228  }
229 
230  Info<< " mag(UMean) max/min : "
231  << max(mag(UMean)).value() << " "
232  << min(mag(UMean)).value() << endl;
233 
234  Info<< " translationalT max/min : "
235  << max(translationalT).value() << " "
236  << min(translationalT).value() << endl;
237 
238  Info<< " internalT max/min : "
239  << max(internalT).value() << " "
240  << min(internalT).value() << endl;
241 
242  Info<< " overallT max/min : "
243  << max(overallT).value() << " "
244  << min(overallT).value() << endl;
245 
246  Info<< " p max/min : "
247  << max(p).value() << " "
248  << min(p).value() << endl;
249 
250  UMean.write();
251 
252  translationalT.write();
253 
254  internalT.write();
255 
256  overallT.write();
257 
258  p.write();
259 
260  Info<< "dsmcFields written." << nl << endl;
261  }
262  else
263  {
264  Info<< "Small value (" << min(mag(rhoNMean))
265  << ") found in rhoNMean field. "
266  << "Not calculating dsmcFields to avoid division by zero."
267  << endl;
268  }
269  }
270 }
271 
272 
273 // ************************ vim: set sw=4 sts=4 et: ************************ //