potentialFoam.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) 2011-2016 OpenFOAM Foundation
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 Application
25  potentialFoam
26 
27 Description
28  Potential flow solver which solves for the velocity potential
29  from which the flux-field is obtained and velocity field by reconstructing
30  the flux.
31 
32  This application is particularly useful to generate starting fields for
33  Navier-Stokes codes.
34 
35 \*---------------------------------------------------------------------------*/
36 
37 #include "fvCFD.H"
38 #include "pisoControl.H"
39 
40 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
41 
42 int main(int argc, char *argv[])
43 {
44  argList::addOption
45  (
46  "pName",
47  "pName",
48  "Name of the pressure field"
49  );
50 
51  argList::addBoolOption
52  (
53  "initialiseUBCs",
54  "Initialise U boundary conditions"
55  );
56 
57  argList::addBoolOption
58  (
59  "writePhi",
60  "Write the velocity potential field"
61  );
62 
63  argList::addBoolOption
64  (
65  "writep",
66  "Calculate and write the pressure field"
67  );
68 
69  argList::addBoolOption
70  (
71  "withFunctionObjects",
72  "execute functionObjects"
73  );
74 
75  #include "setRootCase.H"
76  #include "createTime.H"
77  #include "createMesh.H"
78 
79  pisoControl potentialFlow(mesh, "potentialFlow");
80 
81  #include "createFields.H"
82  #include "createMRF.H"
83 
84  // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
85 
86  Info<< nl << "Calculating potential flow" << endl;
87 
88  // Since solver contains no time loop it would never execute
89  // function objects so do it ourselves
90  runTime.functionObjects().start();
91 
92  MRF.makeRelative(phi);
93  adjustPhi(phi, U, p);
94 
95  // Non-orthogonal velocity potential corrector loop
96  while (potentialFlow.correctNonOrthogonal())
97  {
98  fvScalarMatrix PhiEqn
99  (
101  ==
102  fvc::div(phi)
103  );
104 
105  PhiEqn.setReference(PhiRefCell, PhiRefValue);
106  PhiEqn.solve();
107 
108  if (potentialFlow.finalNonOrthogonalIter())
109  {
110  phi -= PhiEqn.flux();
111  }
112  }
113 
114  MRF.makeAbsolute(phi);
115 
116  Info<< "Continuity error = "
117  << mag(fvc::div(phi))().weightedAverage(mesh.V()).value()
118  << endl;
119 
121  U.correctBoundaryConditions();
122 
123  Info<< "Interpolated velocity error = "
124  << (sqrt(sum(sqr((fvc::interpolate(U) & mesh.Sf()) - phi)))
125  /sum(mesh.magSf())).value()
126  << endl;
127 
128  // Write U and phi
129  U.write();
130  phi.write();
131 
132  // Optionally write Phi
133  if (args.optionFound("writePhi"))
134  {
135  Phi.write();
136  }
137 
138  // Calculate the pressure field
139  if (args.optionFound("writep"))
140  {
141  Info<< nl << "Calculating approximate pressure field" << endl;
142 
143  label pRefCell = 0;
144  scalar pRefValue = 0.0;
145  setRefCell
146  (
147  p,
148  potentialFlow.dict(),
149  pRefCell,
150  pRefValue
151  );
152 
153  // Calculate the flow-direction filter tensor
154  volScalarField magSqrU(magSqr(U));
155  volSymmTensorField F(sqr(U)/(magSqrU + SMALL*average(magSqrU)));
156 
157  // Calculate the divergence of the flow-direction filtered div(U*U)
158  // Filtering with the flow-direction generates a more reasonable
159  // pressure distribution in regions of high velocity gradient in the
160  // direction of the flow
161  volScalarField divDivUU
162  (
163  fvc::div
164  (
165  F & fvc::div(phi, U),
166  "div(div(phi,U))"
167  )
168  );
169 
170  // Solve a Poisson equation for the approximate pressure
171  while (potentialFlow.correctNonOrthogonal())
172  {
173  fvScalarMatrix pEqn
174  (
175  fvm::laplacian(p) + divDivUU
176  );
177 
178  pEqn.setReference(pRefCell, pRefValue);
179  pEqn.solve();
180  }
181 
182  p.write();
183  }
184 
185  runTime.functionObjects().end();
186 
187  Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
188  << " ClockTime = " << runTime.elapsedClockTime() << " s"
189  << nl << endl;
190 
191  Info<< "End\n" << endl;
192 
193  return 0;
194 }
195 
196 
197 // ************************************************************************* //
dimensionedScalar sqrt(const dimensionedScalar &ds)
dimensioned< scalar > mag(const dimensioned< Type > &)
const dictionary & potentialFlow(mesh.solutionDict().subDict("potentialFlow"))
dimensioned< scalar > magSqr(const dimensioned< Type > &)
tmp< GeometricField< typename outerProduct< vector, Type >::type, fvPatchField, volMesh >> reconstruct(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
bool adjustPhi(surfaceScalarField &phi, const volVectorField &U, volScalarField &p)
Adjust the balance of fluxes to obey continuity.
Definition: adjustPhi.C:34
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
tmp< GeometricField< Type, fvPatchField, volMesh > > laplacian(const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvcLaplacian.C:45
messageStream Info
tmp< surfaceScalarField > interpolate(const RhoType &rho)
dynamicFvMesh & mesh
const scalar pRefValue
fvMatrix< scalar > fvScalarMatrix
Definition: fvMatricesFwd.H:42
static const char nl
Definition: Ostream.H:260
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
IOMRFZoneList & MRF
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition: volFieldsFwd.H:52
volScalarField & p
Definition: createFields.H:51
bool optionFound(const word &opt) const
Return true if the named option is found.
Definition: argListI.H:108
GeometricField< symmTensor, fvPatchField, volMesh > volSymmTensorField
Definition: volFieldsFwd.H:58
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
surfaceScalarField & phi
void setRefCell(const volScalarField &field, const volScalarField &fieldRef, const dictionary &dict, label &refCelli, scalar &refValue, const bool forceReference=false)
If the field fieldRef needs referencing find the reference cell nearest.
Definition: findRefCell.C:31
Foam::argList args(argc, argv)
tmp< GeometricField< Type, fvPatchField, volMesh > > average(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Area-weighted average a surfaceField creating a volField.
Definition: fvcAverage.C:46
const label pRefCell
const dimensionSet dimless(0, 0, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:47
tmp< GeometricField< Type, fvPatchField, volMesh > > div(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcDiv.C:47
int main(int argc, char *argv[])
Definition: postCalc.C:54
dimensionedSymmTensor sqr(const dimensionedVector &dv)
const dimensionedScalar F
Faraday constant: default SI units: [C/mol].
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
U
Definition: pEqn.H:82