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-------------------------------------------------------------------------------
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 mhdFoam
28
29Group
30 grpElectroMagneticsSolvers
31
32Description
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
61int 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 {
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));
113 (
114 "phiHbyA",
115 fvc::flux(HbyA)
116 + rAUf*fvc::ddtCorr(U, phi)
117 );
118
119 // Update the pressure BCs to ensure flux consistency
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// ************************************************************************* //
static const Foam::dimensionedScalar B("", Foam::dimless, 18.678)
Required Classes.
pisoControl bpiso(mesh, "BPISO")
surfaceScalarField & phi
const scalar pRefValue
const label pRefCell
U
Definition: pEqn.H:72
volScalarField & p
fvVectorMatrix & UEqn
Definition: UEqn.H:13
phiHbyA
Definition: pcEqn.H:73
HbyA
Definition: pcEqn.H:74
dynamicFvMesh & mesh
engineTime & runTime
Required Variables.
surfaceScalarField & phiB
Definition: createPhiB.H:47
pisoControl piso(mesh)
surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU))
tmp< volScalarField > rAU
Definition: initCorrectPhi.H:1
GeometricField< vector, fvPatchField, volMesh > volVectorField
Definition: volFieldsFwd.H:83
fvMatrix< scalar > fvScalarMatrix
Definition: fvMatricesFwd.H:45
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition: volFieldsFwd.H:82
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
fvMatrix< vector > fvVectorMatrix
Definition: fvMatricesFwd.H:46
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
Execute application functionObjects to post-process existing results.
volScalarField & nu
CEqn solve()