overPotentialFoam.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 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
26Application
27 overPotentialFoam
28
29Group
30 grpBasicSolvers
31
32Description
33 Potential flow solver which solves for the velocity potential, to
34 calculate the flux-field, from which the velocity field is obtained by
35 reconstructing the flux.
36
37 \heading Solver details
38 The potential flow solution is typically employed to generate initial fields
39 for full Navier-Stokes codes. The flow is evolved using the equation:
40
41 \f[
42 \laplacian \Phi = \div(\vec{U})
43 \f]
44
45 Where:
46 \vartable
47 \Phi | Velocity potential [m2/s]
48 \vec{U} | Velocity [m/s]
49 \endvartable
50
51 The corresponding pressure field could be calculated from the divergence
52 of the Euler equation:
53
54 \f[
55 \laplacian p + \div(\div(\vec{U}\otimes\vec{U})) = 0
56 \f]
57
58 but this generates excessive pressure variation in regions of large
59 velocity gradient normal to the flow direction. A better option is to
60 calculate the pressure field corresponding to velocity variation along the
61 stream-lines:
62
63 \f[
64 \laplacian p + \div(\vec{F}\cdot\div(\vec{U}\otimes\vec{U})) = 0
65 \f]
66 where the flow direction tensor \f$\vec{F}\f$ is obtained from
67 \f[
68 \vec{F} = \hat{\vec{U}}\otimes\hat{\vec{U}}
69 \f]
70
71 \heading Required fields
72 \plaintable
73 U | Velocity [m/s]
74 \endplaintable
75
76 \heading Optional fields
77 \plaintable
78 p | Kinematic pressure [m2/s2]
79 Phi | Velocity potential [m2/s]
80 | Generated from p (if present) or U if not present
81 \endplaintable
82
83 \heading Options
84 \plaintable
85 -writep | write the Euler pressure
86 -writePhi | Write the final velocity potential
87 -initialiseUBCs | Update the velocity boundaries before solving for Phi
88 \endplaintable
89
90\*---------------------------------------------------------------------------*/
91
92#include "fvCFD.H"
93#include "pisoControl.H"
94#include "dynamicFvMesh.H"
96#include "localMin.H"
97
98// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
99
100int main(int argc, char *argv[])
101{
102 argList::addNote
103 (
104 "Overset potential flow solver which solves for the velocity potential"
105 );
106
107 argList::addOption
108 (
109 "pName",
110 "pName",
111 "Name of the pressure field"
112 );
113
114 argList::addBoolOption
115 (
116 "initialiseUBCs",
117 "Initialise U boundary conditions"
118 );
119
120 argList::addBoolOption
121 (
122 "writePhi",
123 "Write the final velocity potential field"
124 );
125
126 argList::addBoolOption
127 (
128 "writep",
129 "Calculate and write the Euler pressure field"
130 );
131
132 argList::addBoolOption
133 (
134 "withFunctionObjects",
135 "Execute functionObjects"
136 );
137
138 #include "setRootCaseLists.H"
139 #include "createTime.H"
141
142 pisoControl potentialFlow(mesh, "potentialFlow");
143
144 #include "createFields.H"
145
146 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
147
148 Info<< nl << "Calculating potential flow" << endl;
149
150 mesh.update();
151
152 surfaceScalarField faceMask(localMin<scalar>(mesh).interpolate(cellMask));
153
154 // Since solver contains no time loop it would never execute
155 // function objects so do it ourselves
156 runTime.functionObjects().start();
157
158 MRF.makeRelative(phi);
159 adjustPhi(phi, U, p);
160
161 // Non-orthogonal velocity potential corrector loop
162 while (potentialFlow.correct())
163 {
164 phi = fvc::flux(U);
165
166 while (potentialFlow.correctNonOrthogonal())
167 {
168 fvScalarMatrix PhiEqn
169 (
170 fvm::laplacian(faceMask, Phi)
171 ==
172 fvc::div(phi)
173 );
174
175 PhiEqn.setReference(PhiRefCell, PhiRefValue);
176 PhiEqn.solve();
177
178 if (potentialFlow.finalNonOrthogonalIter())
179 {
180 phi -= PhiEqn.flux();
181 }
182 }
183
184 MRF.makeAbsolute(phi);
185
186 Info<< "Continuity error = "
187 << mag(fvc::div(phi))().weightedAverage(mesh.V()).value()
188 << endl;
189
190 U = fvc::reconstruct(phi);
191 U.correctBoundaryConditions();
192
193 Info<< "Interpolated velocity error = "
194 << (sqrt(sum(sqr(fvc::flux(U) - phi)))/sum(mesh.magSf())).value()
195 << endl;
196 }
197
198 // Write U and phi
199 U.write();
200 phi.write();
201
202 // Optionally write Phi
203 if (args.found("writePhi"))
204 {
205 Phi.write();
206 }
207
208 // Calculate the pressure field from the Euler equation
209 if (args.found("writep"))
210 {
211 Info<< nl << "Calculating approximate pressure field" << endl;
212
213 label pRefCell = 0;
214 scalar pRefValue = 0.0;
216 (
217 p,
218 potentialFlow.dict(),
219 pRefCell,
221 );
222
223 // Calculate the flow-direction filter tensor
224 volScalarField magSqrU(magSqr(U));
225 volSymmTensorField F(sqr(U)/(magSqrU + SMALL*average(magSqrU)));
226
227 // Calculate the divergence of the flow-direction filtered div(U*U)
228 // Filtering with the flow-direction generates a more reasonable
229 // pressure distribution in regions of high velocity gradient in the
230 // direction of the flow
231 volScalarField divDivUU
232 (
233 fvc::div
234 (
235 F & fvc::div(phi, U),
236 "div(div(phi,U))"
237 )
238 );
239
240 // Solve a Poisson equation for the approximate pressure
241 while (potentialFlow.correctNonOrthogonal())
242 {
243 fvScalarMatrix pEqn
244 (
245 fvm::laplacian(p) + divDivUU
246 );
247
248 pEqn.setReference(pRefCell, pRefValue);
249 pEqn.solve();
250 }
251
252 p.write();
253 }
254
255 runTime.functionObjects().end();
256
257 runTime.printExecutionTime(Info);
258
259 Info<< "End\n" << endl;
260
261 return 0;
262}
263
264
265// ************************************************************************* //
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
surfaceScalarField faceMask(localMin< scalar >(mesh).interpolate(cellMask))
mesh interpolate(rAU)
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
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)
GeometricField< scalar, fvsPatchField, surfaceMesh > surfaceScalarField
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)