setFlow.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 | www.openfoam.com
6 \\/ M anipulation |
7-------------------------------------------------------------------------------
8 Copyright (C) 2017-2020 OpenCFD Ltd.
9-------------------------------------------------------------------------------
10License
11 This file is part of OpenFOAM.
12
13 OpenFOAM is free software: you can redistribute it and/or modify it
14 under the terms of the GNU General Public License as published by
15 the Free Software Foundation, either version 3 of the License, or
16 (at your option) any later version.
17
18 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
19 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
20 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
21 for more details.
22
23 You should have received a copy of the GNU General Public License
24 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
25
26\*---------------------------------------------------------------------------*/
27
28#include "setFlow.H"
29#include "volFields.H"
30#include "surfaceFields.H"
31#include "fvcFlux.H"
33#include "fvcSurfaceIntegrate.H"
34
35// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
36
37namespace Foam
38{
39namespace functionObjects
40{
43}
44}
45
46
47const Foam::Enum
48<
49 Foam::functionObjects::setFlow::modeType
50>
51Foam::functionObjects::setFlow::modeTypeNames
52({
53 { functionObjects::setFlow::modeType::FUNCTION, "function" },
54 { functionObjects::setFlow::modeType::ROTATION, "rotation" },
55 { functionObjects::setFlow::modeType::VORTEX2D, "vortex2D" },
56 { functionObjects::setFlow::modeType::VORTEX3D, "vortex3D" },
57});
58
59
60// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
61
62void Foam::functionObjects::setFlow::setPhi(const volVectorField& U)
63{
64 surfaceScalarField* phiptr =
66
67 if (!phiptr)
68 {
69 return;
70 }
71
72 if (rhoName_ != "none")
73 {
74 const volScalarField* rhoptr =
76
77 if (rhoptr)
78 {
79 const volScalarField& rho = *rhoptr;
80 *phiptr = fvc::flux(rho*U);
81 }
82 else
83 {
85 << "Unable to find rho field'" << rhoName_
86 << "' in the mesh database. Available fields are:"
88 << exit(FatalError);
89 }
90 }
91 else
92 {
93 *phiptr = fvc::flux(U);
94 }
95}
96
97
98// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
99
101(
102 const word& name,
103 const Time& runTime,
104 const dictionary& dict
105)
106:
108 mode_(modeType::FUNCTION),
109 UName_("U"),
110 rhoName_("none"),
111 phiName_("phi"),
112 reverseTime_(VGREAT),
113 scalePtr_(nullptr),
114 origin_(Zero),
115 R_(tensor::I),
116 omegaPtr_(nullptr),
117 velocityPtr_(nullptr)
118{
119 read(dict);
120}
121
122
123// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
124
126{
128 {
129 Info<< name() << ":" << endl;
130
131 modeTypeNames.readEntry("mode", dict, mode_);
132
133 Info<< " operating mode: " << modeTypeNames[mode_] << endl;
134
135 if (dict.readIfPresent("U", UName_))
136 {
137 Info<< " U field name: " << UName_ << endl;
138 }
139
140 if (dict.readIfPresent("rho", rhoName_))
141 {
142 Info<< " rho field name: " << rhoName_ << endl;
143 }
144
145 if (dict.readIfPresent("phi", phiName_))
146 {
147 Info<< " phi field name: " << phiName_ << endl;
148 }
149
150 if (dict.readIfPresent("reverseTime", reverseTime_))
151 {
152 Info<< " reverse flow direction at time: " << reverseTime_
153 << endl;
154 reverseTime_ = mesh_.time().userTimeToTime(reverseTime_);
155 }
156
157 // Scaling is applied across all modes
158 scalePtr_ = Function1<scalar>::New("scale", dict, &mesh_);
159
160 switch (mode_)
161 {
162 case modeType::FUNCTION:
163 {
164 velocityPtr_ = Function1<vector>::New("velocity", dict, &mesh_);
165 break;
166 }
167 case modeType::ROTATION:
168 {
169 omegaPtr_ = Function1<scalar>::New("omega", dict, &mesh_);
170
171 dict.readEntry("origin", origin_);
172 const vector refDir(dict.get<vector>("refDir").normalise());
173 const vector axis(dict.get<vector>("axis").normalise());
174
175 R_ = tensor(refDir, axis, refDir^axis);
176 break;
177 }
178 case modeType::VORTEX2D:
179 case modeType::VORTEX3D:
180 {
181 dict.readEntry("origin", origin_);
182 const vector refDir(dict.get<vector>("refDir").normalise());
183 const vector axis(dict.get<vector>("axis").normalise());
184
185 R_ = tensor(refDir, axis, refDir^axis);
186 break;
187 }
188 }
189
190 Info<< endl;
191
192 return true;
193 }
194
195 return false;
196}
197
198
200{
201 volVectorField* Uptr =
202 mesh_.getObjectPtr<volVectorField>(UName_);
203
204 surfaceScalarField* phiptr =
205 mesh_.getObjectPtr<surfaceScalarField>(phiName_);
206
207 Log << nl << name() << ":" << nl;
208
209 if (!Uptr || !phiptr)
210 {
211 Info<< " Either field " << UName_ << " or " << phiName_
212 << " not found in the mesh database" << nl;
213
214 return true;
215 }
216
217 const scalar t = mesh_.time().timeOutputValue();
218
219 Log << " setting " << UName_ << " and " << phiName_ << nl;
220
221 volVectorField& U = *Uptr;
222 surfaceScalarField& phi = *phiptr;
223
224 switch (mode_)
225 {
226 case modeType::FUNCTION:
227 {
228 const vector Uc = velocityPtr_->value(t);
229 U == dimensionedVector("Uc", dimVelocity, Uc);
230 U.correctBoundaryConditions();
231 setPhi(U);
232
233 break;
234 }
235 case modeType::ROTATION:
236 {
237 const volVectorField& C = mesh_.C();
238 const volVectorField d
239 (
240 typeName + ":d",
241 C - dimensionedVector("origin", dimLength, origin_)
242 );
244 const scalarField z(d.component(vector::Z));
245
246 const scalar omega = omegaPtr_->value(t);
247 vectorField& Uc = U.primitiveFieldRef();
248 Uc.replace(vector::X, -omega*z);
249 Uc.replace(vector::Y, scalar(0));
250 Uc.replace(vector::Z, omega*x);
251
252 volVectorField::Boundary& Ubf = U.boundaryFieldRef();
253 forAll(Ubf, patchi)
254 {
255 vectorField& Uf = Ubf[patchi];
256 if (Uf.size())
257 {
258 const vectorField& Cp = C.boundaryField()[patchi];
259 const vectorField dp(Cp - origin_);
260 const scalarField xp(dp.component(vector::X));
261 const scalarField zp(dp.component(vector::Z));
262 Uf.replace(vector::X, -omega*zp);
263 Uf.replace(vector::Y, scalar(0));
264 Uf.replace(vector::Z, omega*xp);
265 }
266 }
267
268 U = U & R_;
269 U.correctBoundaryConditions();
270 setPhi(U);
271
272 break;
273 }
274 case modeType::VORTEX2D:
275 {
276 const scalar pi = Foam::constant::mathematical::pi;
277
278 const volVectorField& C = mesh_.C();
279
280 const volVectorField d
281 (
282 typeName + ":d",
283 C - dimensionedVector("origin", dimLength, origin_)
284 );
286 const scalarField z(d.component(vector::Z));
287
288 vectorField& Uc = U.primitiveFieldRef();
289 Uc.replace(vector::X, -sin(2*pi*z)*sqr(sin(pi*x)));
290 Uc.replace(vector::Y, scalar(0));
291 Uc.replace(vector::Z, sin(2*pi*x)*sqr(sin(pi*z)));
292
293 U = U & R_;
294 U.correctBoundaryConditions();
295
296 // Calculating phi
297 // Note: R_ rotation not implemented in phi calculation
298 const vectorField Cf(mesh_.Cf().primitiveField() - origin_);
299 const scalarField Xf(Cf.component(vector::X));
300 const scalarField Yf(Cf.component(vector::Y));
301 const scalarField Zf(Cf.component(vector::Z));
302 vectorField Uf(Xf.size());
303 Uf.replace(vector::X, -sin(2*pi*Zf)*sqr(sin(pi*Xf)));
304 Uf.replace(vector::Y, scalar(0));
305 Uf.replace(vector::Z, sin(2*pi*Xf)*sqr(sin(pi*Zf)));
306
307 scalarField& phic = phi.primitiveFieldRef();
308 const vectorField& Sfc = mesh_.Sf().primitiveField();
309 phic = Uf & Sfc;
310
311 surfaceScalarField::Boundary& phibf = phi.boundaryFieldRef();
312 const surfaceVectorField::Boundary& Sfbf =
313 mesh_.Sf().boundaryField();
314 const surfaceVectorField::Boundary& Cfbf =
315 mesh_.Cf().boundaryField();
316
317 forAll(phibf, patchi)
318 {
319 scalarField& phif = phibf[patchi];
320 const vectorField& Sff = Sfbf[patchi];
321 const vectorField& Cff = Cfbf[patchi];
322 const scalarField xf(Cff.component(vector::X));
323 const scalarField yf(Cff.component(vector::Y));
324 const scalarField zf(Cff.component(vector::Z));
325 vectorField Ufb(xf.size());
326 Ufb.replace(vector::X, -sin(2*pi*zf)*sqr(sin(pi*xf)));
327 Ufb.replace(vector::Y, scalar(0));
328 Ufb.replace(vector::Z, sin(2*pi*xf)*sqr(sin(pi*zf)));
329 phif = Ufb & Sff;
330 }
331
332 break;
333 }
334 case modeType::VORTEX3D:
335 {
336 const scalar pi = Foam::constant::mathematical::pi;
337 const volVectorField& C = mesh_.C();
338
339 const volVectorField d
340 (
341 typeName + ":d",
342 C - dimensionedVector("origin", dimLength, origin_)
343 );
346 const scalarField z(d.component(vector::Z));
347
348 vectorField& Uc = U.primitiveFieldRef();
349 Uc.replace(vector::X, 2*sqr(sin(pi*x))*sin(2*pi*y)*sin(2*pi*z));
350 Uc.replace(vector::Y, -sin(2*pi*x)*sqr(sin(pi*y))*sin(2*pi*z));
351 Uc.replace(vector::Z, -sin(2*pi*x)*sin(2*pi*y)*sqr(sin(pi*z)));
352
353 U = U & R_;
354 U.correctBoundaryConditions();
355
356 // Calculating phi
357 // Note: R_ rotation not implemented in phi calculation
358 const vectorField Cf(mesh_.Cf().primitiveField() - origin_);
359 const scalarField Xf(Cf.component(vector::X));
360 const scalarField Yf(Cf.component(vector::Y));
361 const scalarField Zf(Cf.component(vector::Z));
362 vectorField Uf(Xf.size());
363 Uf.replace(0, 2*sqr(sin(pi*Xf))*sin(2*pi*Yf)*sin(2*pi*Zf));
364 Uf.replace(1, -sin(2*pi*Xf)*sqr(sin(pi*Yf))*sin(2*pi*Zf));
365 Uf.replace(2, -sin(2*pi*Xf)*sin(2*pi*Yf)*sqr(sin(pi*Zf)));
366
367 scalarField& phic = phi.primitiveFieldRef();
368 const vectorField& Sfc = mesh_.Sf().primitiveField();
369 phic = Uf & Sfc;
370
371 surfaceScalarField::Boundary& phibf = phi.boundaryFieldRef();
372 const surfaceVectorField::Boundary& Sfbf =
373 mesh_.Sf().boundaryField();
374 const surfaceVectorField::Boundary& Cfbf =
375 mesh_.Cf().boundaryField();
376
377 forAll(phibf, patchi)
378 {
379 scalarField& phif = phibf[patchi];
380 const vectorField& Sff = Sfbf[patchi];
381 const vectorField& Cff = Cfbf[patchi];
382 const scalarField xf(Cff.component(vector::X));
383 const scalarField yf(Cff.component(vector::Y));
384 const scalarField zf(Cff.component(vector::Z));
385 vectorField Uf(xf.size());
386 Uf.replace(0, 2*sqr(sin(pi*xf))*sin(2*pi*yf)*sin(2*pi*zf));
387 Uf.replace(1, -sin(2*pi*xf)*sqr(sin(pi*yf))*sin(2*pi*zf));
388 Uf.replace(2, -sin(2*pi*xf)*sin(2*pi*yf)*sqr(sin(pi*zf)));
389 phif = Uf & Sff;
390 }
391
392 break;
393 }
394 }
395
396 if (t > reverseTime_)
397 {
398 Log << " flow direction: reverse" << nl;
399 U.negate();
400 phi.negate();
401 }
402
403 // Apply scaling
404 const scalar s = scalePtr_->value(t);
405 U *= s;
406 phi *= s;
407
408 U.correctBoundaryConditions();
409
411 Log << " Continuity error: max(mag(sum(phi))) = "
412 << gMax(mag(sumPhi)) << nl << endl;
413
414 return true;
415}
416
417
419{
420 const auto* Uptr = mesh_.findObject<volVectorField>(UName_);
421 if (Uptr)
422 {
423 Uptr->write();
424 }
425
426 const auto* phiptr = mesh_.findObject<surfaceScalarField>(phiName_);
427 if (phiptr)
428 {
429 phiptr->write();
430 }
431
432 return true;
433}
434
435
436// ************************************************************************* //
scalar y
surfaceScalarField phic(mixture.cAlpha() *mag(alphaPhic/mesh.magSf()))
#define Log
Definition: PDRblock.C:35
Macros for easy insertion into run-time selection tables.
#define addToRunTimeSelectionTable(baseType, thisType, argNames)
Add to construction table with typeName as the key.
surfaceScalarField & phi
Graphite solid properties.
Definition: C.H:53
C()
Construct null.
Definition: C.C:43
Enum is a wrapper around a list of names/values that represent particular enumeration (or int) values...
Definition: Enum.H:61
void replace(const direction, const UList< cmptType > &)
Replace a component field of the field.
Definition: Field.C:557
tmp< Field< cmptType > > component(const direction) const
Return a component field of the field.
Definition: Field.C:545
tmp< GeometricField< cmptType, PatchField, GeoMesh > > component(const direction) const
Return a component of the field.
virtual bool read()
Re-read model coefficients if they have changed.
Class to control time during OpenFOAM simulations that is also the top-level objectRegistry.
Definition: Time.H:80
static autoPtr< Time > New()
Construct (dummy) Time - no functionObjects or libraries.
Definition: Time.C:717
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
Vector< Cmpt > & normalise(const scalar tol=ROOTVSMALL)
Inplace normalise the vector by its magnitude.
Definition: VectorI.H:123
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Definition: dictionary.H:126
Abstract base-class for Time/database function objects.
Specialization of Foam::functionObject for an Foam::fvMesh, providing a reference to the Foam::fvMesh...
const fvMesh & mesh_
Reference to the fvMesh.
Computes the magnitude of an input field.
Definition: mag.H:153
Provides options to set the velocity and flux fields as a function of time.
Definition: setFlow.H:275
virtual bool read(const dictionary &dict)
Read the setFlow data.
Definition: setFlow.C:125
virtual bool execute()
Do nothing.
Definition: setFlow.C:199
virtual bool write()
Calculate the setFlow and write.
Definition: setFlow.C:418
refPtr< surfaceScalarField > setPhi()
Return cell face motion fluxes (or null)
wordList names() const
The unsorted names of all objects.
const Type * findObject(const word &name, const bool recursive=false) const
Return const pointer to the object of the given Type.
Type * getObjectPtr(const word &name, const bool recursive=false) const
virtual bool write(const bool valid=true) const
Write using setting from DB.
A class for handling words, derived from Foam::string.
Definition: word.H:68
#define defineTypeNameAndDebug(Type, DebugSwitch)
Define the typeName and debug information.
Definition: className.H:121
U
Definition: pEqn.H:72
const volScalarField & Cp
Definition: EEqn.H:7
engineTime & runTime
autoPtr< surfaceVectorField > Uf
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
Calculate the face-flux of the given field.
Surface integrate surfaceField creating a volField. Surface sum a surfaceField creating a volField.
gmvFile<< "tracers "<< particles.size()<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().x()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().y()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){ word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
constexpr scalar pi(M_PI)
tmp< surfaceScalarField > flux(const volVectorField &vvf)
Return the face-flux field obtained from the given volVectorField.
Definition: fvcFlux.C:34
void surfaceIntegrate(Field< Type > &ivf, const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Namespace for OpenFOAM.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
const dimensionSet dimLength(0, 1, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:52
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition: volFieldsFwd.H:82
dimensionedScalar sin(const dimensionedScalar &ds)
const dimensionSet dimVelocity
messageStream Info
Information stream (stdout output on master, null elsewhere)
static const Identity< scalar > I
Definition: Identity.H:94
GeometricField< scalar, fvsPatchField, surfaceMesh > surfaceScalarField
Tensor< scalar > tensor
Definition: symmTensor.H:61
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
dimensioned< vector > dimensionedVector
Dimensioned vector obtained from generic dimensioned type.
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
error FatalError
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
Definition: exprTraits.C:59
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:130
Type gMax(const FieldField< Field, Type > &f)
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
dictionary dict
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333
Foam::surfaceFields.