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 | www.openfoam.com
6 \\/ M anipulation |
7-------------------------------------------------------------------------------
8 Copyright (C) 2011-2016 OpenFOAM Foundation
9 Copyright (C) 2019 OpenCFD Ltd.
10-------------------------------------------------------------------------------
11License
12 This file is part of OpenFOAM.
13
14 OpenFOAM is free software: you can redistribute it and/or modify it
15 under the terms of the GNU General Public License as published by
16 the Free Software Foundation, either version 3 of the License, or
17 (at your option) any later version.
18
19 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22 for more details.
23
24 You should have received a copy of the GNU General Public License
25 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26
27Application
28 potentialFoam
29
30Group
31 grpBasicSolvers
32
33Description
34 Potential flow solver which solves for the velocity potential, to
35 calculate the flux-field, from which the velocity field is obtained by
36 reconstructing the flux.
37
38 \heading Solver details
39 The potential flow solution is typically employed to generate initial fields
40 for full Navier-Stokes codes. The flow is evolved using the equation:
41
42 \f[
43 \laplacian \Phi = \div(\vec{U})
44 \f]
45
46 Where:
47 \vartable
48 \Phi | Velocity potential [m2/s]
49 \vec{U} | Velocity [m/s]
50 \endvartable
51
52 The corresponding pressure field could be calculated from the divergence
53 of the Euler equation:
54
55 \f[
56 \laplacian p + \div(\div(\vec{U}\otimes\vec{U})) = 0
57 \f]
58
59 but this generates excessive pressure variation in regions of large
60 velocity gradient normal to the flow direction. A better option is to
61 calculate the pressure field corresponding to velocity variation along the
62 stream-lines:
63
64 \f[
65 \laplacian p + \div(\vec{F}\cdot\div(\vec{U}\otimes\vec{U})) = 0
66 \f]
67 where the flow direction tensor \f$\vec{F}\f$ is obtained from
68 \f[
69 \vec{F} = \hat{\vec{U}}\otimes\hat{\vec{U}}
70 \f]
71
72 \heading Required fields
73 \plaintable
74 U | Velocity [m/s]
75 \endplaintable
76
77 \heading Optional fields
78 \plaintable
79 p | Kinematic pressure [m2/s2]
80 Phi | Velocity potential [m2/s]
81 | Generated from p (if present) or U if not present
82 \endplaintable
83
84 \heading Options
85 \plaintable
86 -writep | write the Euler pressure
87 -writephi | Write the final volumetric flux
88 -writePhi | Write the final velocity potential
89 -initialiseUBCs | Update the velocity boundaries before solving for Phi
90 \endplaintable
91
92\*---------------------------------------------------------------------------*/
93
94#include "fvCFD.H"
95#include "pisoControl.H"
96
97// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
98
99int main(int argc, char *argv[])
100{
101 argList::addNote
102 (
103 "Potential flow solver which solves for the velocity potential"
104 );
105
106 argList::addOption
107 (
108 "pName",
109 "pName",
110 "Name of the pressure field"
111 );
112
113 argList::addBoolOption
114 (
115 "initialiseUBCs",
116 "Initialise U boundary conditions"
117 );
118
119 argList::addBoolOption
120 (
121 "writephi",
122 "Write the final volumetric flux field"
123 );
124
125 argList::addBoolOption
126 (
127 "writePhi",
128 "Write the final velocity potential field"
129 );
130
131 argList::addBoolOption
132 (
133 "writep",
134 "Calculate and write the Euler pressure field"
135 );
136
137 argList::addBoolOption
138 (
139 "withFunctionObjects",
140 "Execute functionObjects"
141 );
142
143 #include "addRegionOption.H"
144 #include "addCheckCaseOptions.H"
145 #include "setRootCaseLists.H"
146 #include "createTime.H"
147 #include "createMesh.H"
148
149 pisoControl potentialFlow(mesh, "potentialFlow");
150
151 #include "createFields.H"
152
153 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
154
155 Info<< nl << "Calculating potential flow" << endl;
156
157 // Since solver contains no time loop it would never execute
158 // function objects so do it ourselves
159 runTime.functionObjects().start();
160
161 MRF.makeRelative(phi);
162 adjustPhi(phi, U, p);
163
164 // Non-orthogonal velocity potential corrector loop
165 while (potentialFlow.correctNonOrthogonal())
166 {
167 fvScalarMatrix PhiEqn
168 (
169 fvm::laplacian(dimensionedScalar("1", dimless, 1), Phi)
170 ==
171 fvc::div(phi)
172 );
173
174 PhiEqn.setReference(PhiRefCell, PhiRefValue);
175 PhiEqn.solve();
176
177 if (potentialFlow.finalNonOrthogonalIter())
178 {
179 phi -= PhiEqn.flux();
180 }
181 }
182
183 MRF.makeAbsolute(phi);
184
185 Info<< "Continuity error = "
186 << mag(fvc::div(phi))().weightedAverage(mesh.V()).value()
187 << endl;
188
189 U = fvc::reconstruct(phi);
190 U.correctBoundaryConditions();
191
192 Info<< "Interpolated velocity error = "
193 << (sqrt(sum(sqr(fvc::flux(U) - phi)))/sum(mesh.magSf())).value()
194 << endl;
195
196 // Write U
197 U.write();
198
199 // Optionally write the volumetric flux, phi
200 if (args.found("writephi"))
201 {
202 phi.write();
203 }
204
205 // Optionally write velocity potential, Phi
206 if (args.found("writePhi"))
207 {
208 Phi.write();
209 }
210
211 // Calculate the pressure field from the Euler equation
212 if (args.found("writep"))
213 {
214 Info<< nl << "Calculating approximate pressure field" << endl;
215
216 label pRefCell = 0;
217 scalar pRefValue = 0.0;
219 (
220 p,
221 potentialFlow.dict(),
222 pRefCell,
224 );
225
226 // Calculate the flow-direction filter tensor
227 volScalarField magSqrU(magSqr(U));
228 volSymmTensorField F(sqr(U)/(magSqrU + SMALL*average(magSqrU)));
229
230 // Calculate the divergence of the flow-direction filtered div(U*U)
231 // Filtering with the flow-direction generates a more reasonable
232 // pressure distribution in regions of high velocity gradient in the
233 // direction of the flow
234 volScalarField divDivUU
235 (
236 fvc::div
237 (
238 F & fvc::div(phi, U),
239 "div(div(phi,U))"
240 )
241 );
242
243 // Solve a Poisson equation for the approximate pressure
244 while (potentialFlow.correctNonOrthogonal())
245 {
246 fvScalarMatrix pEqn
247 (
248 fvm::laplacian(p) + divDivUU
249 );
250
251 pEqn.setReference(pRefCell, pRefValue);
252 pEqn.solve();
253 }
254
255 p.write();
256 }
257
258 runTime.functionObjects().end();
259
260 runTime.printExecutionTime(Info);
261
262 Info<< "End\n" << endl;
263
264 return 0;
265}
266
267
268// ************************************************************************* //
Required Classes.
surfaceScalarField & phi
const scalar pRefValue
const label pRefCell
IOMRFZoneList & MRF
bool found(const word &optName) const
Return true if the named option is found.
Definition: argListI.H:178
U
Definition: pEqn.H:72
volScalarField & p
dynamicFvMesh & mesh
engineTime & runTime
Required Variables.
adjustPhi(phiHbyA, U, p_rgh)
volVectorField F(fluid.F())
tmp< GeometricField< Type, faPatchField, areaMesh > > average(const GeometricField< Type, faePatchField, edgeMesh > &ssf)
Area-weighted average a edgeField creating a areaField.
Definition: facAverage.C:47
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
fvMatrix< scalar > fvScalarMatrix
Definition: fvMatricesFwd.H:45
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition: volFieldsFwd.H:82
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
messageStream Info
Information stream (stdout output on master, null elsewhere)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
dimensionedScalar sqrt(const dimensionedScalar &ds)
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
GeometricField< symmTensor, fvPatchField, volMesh > volSymmTensorField
Definition: volFieldsFwd.H:86
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
Foam::argList args(argc, argv)
const dictionary & potentialFlow(mesh.solutionDict().subDict("potentialFlow"))
setRefCell(p, pimple.dict(), pRefCell, pRefValue)