OpenFOAM: API Guide
v2112
The open source CFD toolbox
mhdFoam.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-2018 OpenFOAM Foundation
9
-------------------------------------------------------------------------------
10
License
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
Application
27
mhdFoam
28
29
Group
30
grpElectroMagneticsSolvers
31
32
Description
33
Solver for magnetohydrodynamics (MHD): incompressible, laminar flow of a
34
conducting fluid under the influence of a magnetic field.
35
36
An applied magnetic field H acts as a driving force,
37
at present boundary conditions cannot be set via the
38
electric field E or current density J. The fluid viscosity nu,
39
conductivity sigma and permeability mu are read in as uniform
40
constants.
41
42
A fictitous magnetic flux pressure pH is introduced in order to
43
compensate for discretisation errors and create a magnetic face flux
44
field which is divergence free as required by Maxwell's equations.
45
46
However, in this formulation discretisation error prevents the normal
47
stresses in UB from cancelling with those from BU, but it is unknown
48
whether this is a serious error. A correction could be introduced
49
whereby the normal stresses in the discretised BU term are replaced
50
by those from the UB term, but this would violate the boundedness
51
constraint presently observed in the present numerics which
52
guarantees div(U) and div(H) are zero.
53
54
\*---------------------------------------------------------------------------*/
55
56
#include "
fvCFD.H
"
57
#include "
pisoControl.H
"
58
59
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
60
61
int
main(
int
argc,
char
*argv[])
62
{
63
argList::addNote
64
(
65
"Solver for magnetohydrodynamics (MHD):"
66
" incompressible, laminar flow of a conducting fluid"
67
" under the influence of a magnetic field."
68
);
69
70
#include "
postProcess.H
"
71
72
#include "
addCheckCaseOptions.H
"
73
#include "
setRootCaseLists.H
"
74
#include "
createTime.H
"
75
#include "
createMesh.H
"
76
#include "createControl.H"
77
#include "createFields.H"
78
#include "initContinuityErrs.H"
79
80
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
81
82
Info
<<
nl
<<
"Starting time loop"
<<
endl
;
83
84
while
(
runTime
.loop())
85
{
86
Info
<<
"Time = "
<<
runTime
.timeName() <<
nl
<<
endl
;
87
88
#include "CourantNo.H"
89
90
{
91
fvVectorMatrix
UEqn
92
(
93
fvm::ddt(
U
)
94
+ fvm::div(
phi
,
U
)
95
- fvc::div(
phiB
, 2.0*DBU*
B
)
96
- fvm::laplacian(
nu
,
U
)
97
+ fvc::grad(DBU*
magSqr
(
B
))
98
);
99
100
if
(
piso
.momentumPredictor())
101
{
102
solve
(
UEqn
== -fvc::grad(
p
));
103
}
104
105
106
// --- PISO loop
107
while
(
piso
.correct())
108
{
109
volScalarField
rAU
(1.0/
UEqn
.A());
110
surfaceScalarField
rAUf
(
"rAUf"
, fvc::interpolate(
rAU
));
111
volVectorField
HbyA
(
constrainHbyA
(
rAU
*
UEqn
.H(),
U
,
p
));
112
surfaceScalarField
phiHbyA
113
(
114
"phiHbyA"
,
115
fvc::flux(
HbyA
)
116
+
rAUf
*fvc::ddtCorr(
U
,
phi
)
117
);
118
119
// Update the pressure BCs to ensure flux consistency
120
constrainPressure
(
p
,
U
,
phiHbyA
,
rAUf
);
121
122
while
(
piso
.correctNonOrthogonal())
123
{
124
fvScalarMatrix
pEqn
125
(
126
fvm::laplacian(
rAUf
,
p
) == fvc::div(
phiHbyA
)
127
);
128
129
pEqn.setReference(
pRefCell
,
pRefValue
);
130
pEqn.solve(
mesh
.solver(
p
.select(
piso
.finalInnerIter())));
131
132
if
(
piso
.finalNonOrthogonalIter())
133
{
134
phi
=
phiHbyA
- pEqn.flux();
135
}
136
}
137
138
#include "continuityErrs.H"
139
140
U
=
HbyA
-
rAU
*fvc::grad(
p
);
141
U
.correctBoundaryConditions();
142
}
143
}
144
145
// --- B-PISO loop
146
while
(
bpiso
.correct())
147
{
148
fvVectorMatrix
BEqn
149
(
150
fvm::ddt(
B
)
151
+ fvm::div(
phi
,
B
)
152
- fvc::div(
phiB
,
U
)
153
- fvm::laplacian(DB,
B
)
154
);
155
156
BEqn.solve();
157
158
volScalarField
rAB(1.0/BEqn.A());
159
surfaceScalarField
rABf(
"rABf"
, fvc::interpolate(rAB));
160
161
phiB
= fvc::flux(
B
);
162
163
while
(
bpiso
.correctNonOrthogonal())
164
{
165
fvScalarMatrix
pBEqn
166
(
167
fvm::laplacian(rABf, pB) == fvc::div(
phiB
)
168
);
169
170
pBEqn.solve(
mesh
.solver(pB.select(
bpiso
.finalInnerIter())));
171
172
if
(
bpiso
.finalNonOrthogonalIter())
173
{
174
phiB
-= pBEqn.flux();
175
}
176
}
177
178
#include "
magneticFieldErr.H
"
179
}
180
181
runTime
.write();
182
}
183
184
Info
<<
"End\n"
<<
endl
;
185
186
return
0;
187
}
188
189
190
// ************************************************************************* //
B
static const Foam::dimensionedScalar B("", Foam::dimless, 18.678)
addCheckCaseOptions.H
Required Classes.
bpiso
pisoControl bpiso(mesh, "BPISO")
phi
surfaceScalarField & phi
Definition:
setRegionFluidFields.H:8
pRefValue
const scalar pRefValue
Definition:
setRegionFluidFields.H:37
pRefCell
const label pRefCell
Definition:
setRegionFluidFields.H:36
constrainHbyA
constrainPressure
U
U
Definition:
pEqn.H:72
p
volScalarField & p
Definition:
createFieldRefs.H:8
UEqn
fvVectorMatrix & UEqn
Definition:
UEqn.H:13
phiHbyA
phiHbyA
Definition:
pcEqn.H:73
HbyA
HbyA
Definition:
pcEqn.H:74
mesh
dynamicFvMesh & mesh
Definition:
createDynamicFvMesh.H:6
runTime
engineTime & runTime
Definition:
createEngineTime.H:13
createMesh.H
Required Variables.
phiB
surfaceScalarField & phiB
Definition:
createPhiB.H:47
piso
pisoControl piso(mesh)
createTime.H
fvCFD.H
rAUf
surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU))
rAU
tmp< volScalarField > rAU
Definition:
initCorrectPhi.H:1
magneticFieldErr.H
Foam::volVectorField
GeometricField< vector, fvPatchField, volMesh > volVectorField
Definition:
volFieldsFwd.H:83
Foam::fvScalarMatrix
fvMatrix< scalar > fvScalarMatrix
Definition:
fvMatricesFwd.H:45
Foam::volScalarField
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition:
volFieldsFwd.H:82
Foam::Info
messageStream Info
Information stream (stdout output on master, null elsewhere)
Foam::surfaceScalarField
GeometricField< scalar, fvsPatchField, surfaceMesh > surfaceScalarField
Definition:
surfaceFieldsFwd.H:66
Foam::endl
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition:
Ostream.H:372
Foam::fvVectorMatrix
fvMatrix< vector > fvVectorMatrix
Definition:
fvMatricesFwd.H:46
Foam::magSqr
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
Foam::nl
constexpr char nl
The newline '\n' character (0x0a)
Definition:
Ostream.H:53
pisoControl.H
postProcess.H
Execute application functionObjects to post-process existing results.
nu
volScalarField & nu
Definition:
readMechanicalProperties.H:176
setRootCaseLists.H
solve
CEqn solve()
applications
solvers
electromagnetics
mhdFoam
mhdFoam.C
Generated by
1.9.5
OPENFOAM® is a registered
trademark
of OpenCFD Ltd.