Home
Downloads
Documentation
Installation
User Guide
man-pages
API Documentation
README
Release Notes
Changes
License
Support
SourceForge Project
Main Page
Related Pages
Namespaces
Classes
Files
File List
File Members
All
Classes
Namespaces
Files
Functions
Variables
Typedefs
Enumerations
Enumerator
Friends
Macros
Pages
applications
solvers
electromagnetics
mhdFoam
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 | Copyright (C) 1991-2010 OpenCFD Ltd.
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
mhdFoam
26
27
Description
28
Solver for magnetohydrodynamics (MHD): incompressible, laminar flow of a
29
conducting fluid under the influence of a magnetic field.
30
31
An applied magnetic field H acts as a driving force,
32
at present boundary conditions cannot be set via the
33
electric field E or current density J. The fluid viscosity nu,
34
conductivity sigma and permeability mu are read in as uniform
35
constants.
36
37
A fictitous magnetic flux pressure pH is introduced in order to
38
compensate for discretisation errors and create a magnetic face flux
39
field which is divergence free as required by Maxwell's equations.
40
41
However, in this formulation discretisation error prevents the normal
42
stresses in UB from cancelling with those from BU, but it is unknown
43
whether this is a serious error. A correction could be introduced
44
whereby the normal stresses in the discretised BU term are replaced
45
by those from the UB term, but this would violate the boundedness
46
constraint presently observed in the present numerics which
47
guarantees div(U) and div(H) are zero.
48
49
Usage
50
- mhdFoam [OPTION]
51
52
@param -case <dir> \n
53
Specify the case directory
54
55
@param -parallel \n
56
Run the case in parallel
57
58
@param -help \n
59
Display short usage message
60
61
@param -doc \n
62
Display Doxygen documentation page
63
64
@param -srcDoc \n
65
Display source code
66
67
\*---------------------------------------------------------------------------*/
68
69
#include <
finiteVolume/fvCFD.H
>
70
#include <
OpenFOAM/OSspecific.H
>
71
72
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
73
74
int
main
(
int
argc,
char
*argv[])
75
{
76
#include <
OpenFOAM/setRootCase.H
>
77
78
#include <
OpenFOAM/createTime.H
>
79
#include <
OpenFOAM/createMesh.H
>
80
#include "
createFields.H
"
81
#include <
finiteVolume/initContinuityErrs.H
>
82
83
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
84
85
Info
<<
nl
<<
"Starting time loop"
<<
endl
;
86
87
while
(runTime.loop())
88
{
89
#include <
finiteVolume/readPISOControls.H
>
90
#include "
readBPISOControls.H
"
91
92
Info
<<
"Time = "
<< runTime.timeName() <<
nl
<<
endl
;
93
94
#include <
finiteVolume/CourantNo.H
>
95
96
{
97
fvVectorMatrix
UEqn
98
(
99
fvm::ddt
(
U
)
100
+
fvm::div
(
phi
,
U
)
101
-
fvc::div
(phiB, 2.0*DBU*B)
102
-
fvm::laplacian
(
nu
,
U
)
103
+
fvc::grad
(DBU*
magSqr
(B))
104
);
105
106
solve
(
UEqn
== -
fvc::grad
(
p
));
107
108
109
// --- PISO loop
110
111
for
(
int
corr=0; corr<
nCorr
; corr++)
112
{
113
volScalarField
rUA
= 1.0/
UEqn
.
A
();
114
115
U
= rUA*
UEqn
.
H
();
116
117
phi
= (
fvc::interpolate
(
U
) &
mesh
.
Sf
())
118
+
fvc::ddtPhiCorr
(rUA,
U
,
phi
);
119
120
for
(
int
nonOrth=0; nonOrth<=
nNonOrthCorr
; nonOrth++)
121
{
122
fvScalarMatrix
pEqn
123
(
124
fvm::laplacian
(rUA,
p
) ==
fvc::div
(
phi
)
125
);
126
127
pEqn.setReference(
pRefCell
,
pRefValue
);
128
pEqn.solve();
129
130
if
(nonOrth ==
nNonOrthCorr
)
131
{
132
phi
-= pEqn.flux();
133
}
134
}
135
136
#include <
finiteVolume/continuityErrs.H
>
137
138
U
-= rUA*
fvc::grad
(
p
);
139
U
.correctBoundaryConditions();
140
}
141
}
142
143
// --- B-PISO loop
144
145
for
(
int
Bcorr=0; Bcorr<
nBcorr
; Bcorr++)
146
{
147
fvVectorMatrix
BEqn
148
(
149
fvm::ddt
(B)
150
+
fvm::div
(
phi
, B)
151
-
fvc::div
(phiB,
U
)
152
-
fvm::laplacian
(DB, B)
153
);
154
155
BEqn.solve();
156
157
volScalarField
rBA = 1.0/BEqn.A();
158
159
phiB = (
fvc::interpolate
(B) &
mesh
.
Sf
())
160
+
fvc::ddtPhiCorr
(rBA, B, phiB);
161
162
fvScalarMatrix
pBEqn
163
(
164
fvm::laplacian
(rBA, pB) ==
fvc::div
(phiB)
165
);
166
pBEqn.solve();
167
168
phiB -= pBEqn.flux();
169
170
#include "
magneticFieldErr.H
"
171
}
172
173
runTime.write();
174
}
175
176
Info
<<
"End\n"
<<
endl
;
177
178
return
0;
179
}
180
181
182
// ************************ vim: set sw=4 sts=4 et: ************************ //