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) 2011-2015 OpenFOAM Foundation
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 \*---------------------------------------------------------------------------*/
50 
51 #include "fvCFD.H"
52 #include "pisoControl.H"
53 
54 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
55 
56 int main(int argc, char *argv[])
57 {
58  #include "setRootCase.H"
59 
60  #include "createTime.H"
61  #include "createMesh.H"
62 
63  pisoControl piso(mesh);
64  pisoControl bpiso(mesh, "BPISO");
65 
66  #include "createFields.H"
67  #include "initContinuityErrs.H"
68 
69  // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
70 
71  Info<< nl << "Starting time loop" << endl;
72 
73  while (runTime.loop())
74  {
75  Info<< "Time = " << runTime.timeName() << nl << endl;
76 
77  #include "CourantNo.H"
78 
79  {
81  (
82  fvm::ddt(U)
83  + fvm::div(phi, U)
84  - fvc::div(phiB, 2.0*DBU*B)
85  - fvm::laplacian(nu, U)
86  + fvc::grad(DBU*magSqr(B))
87  );
88 
89  if (piso.momentumPredictor())
90  {
91  solve(UEqn == -fvc::grad(p));
92  }
93 
94 
95  // --- PISO loop
96  while (piso.correct())
97  {
98  volScalarField rAU(1.0/UEqn.A());
100 
101  volVectorField HbyA("HbyA", U);
102  HbyA = rAU*UEqn.H();
103 
105  (
106  "phiHbyA",
107  (fvc::interpolate(HbyA) & mesh.Sf())
108  + rAUf*fvc::ddtCorr(U, phi)
109  );
110 
111  while (piso.correctNonOrthogonal())
112  {
113  fvScalarMatrix pEqn
114  (
116  );
117 
118  pEqn.setReference(pRefCell, pRefValue);
119  pEqn.solve(mesh.solver(p.select(piso.finalInnerIter())));
120 
121  if (piso.finalNonOrthogonalIter())
122  {
123  phi = phiHbyA - pEqn.flux();
124  }
125  }
126 
127  #include "continuityErrs.H"
128 
129  U = HbyA - rAU*fvc::grad(p);
130  U.correctBoundaryConditions();
131  }
132  }
133 
134  // --- B-PISO loop
135  while (bpiso.correct())
136  {
137  fvVectorMatrix BEqn
138  (
139  fvm::ddt(B)
140  + fvm::div(phi, B)
141  - fvc::div(phiB, U)
142  - fvm::laplacian(DB, B)
143  );
144 
145  BEqn.solve();
146 
147  volScalarField rAB(1.0/BEqn.A());
148  surfaceScalarField rABf("rABf", fvc::interpolate(rAB));
149 
150  phiB = (fvc::interpolate(B) & mesh.Sf())
151  + rABf*fvc::ddtCorr(B, phiB);
152 
153  while (bpiso.correctNonOrthogonal())
154  {
155  fvScalarMatrix pBEqn
156  (
157  fvm::laplacian(rABf, pB) == fvc::div(phiB)
158  );
159 
160  pBEqn.solve(mesh.solver(pB.select(bpiso.finalInnerIter())));
161 
162  if (bpiso.finalNonOrthogonalIter())
163  {
164  phiB -= pBEqn.flux();
165  }
166  }
167 
168  #include "magneticFieldErr.H"
169  }
170 
171  runTime.write();
172  }
173 
174  Info<< "End\n" << endl;
175 
176  return 0;
177 }
178 
179 
180 // ************************************************************************* //
rhoEqn solve()
surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU))
dimensioned< scalar > magSqr(const dimensioned< Type > &)
phiHbyA
Definition: pcEqn.H:74
tmp< GeometricField< typename flux< Type >::type, fvsPatchField, surfaceMesh > > ddtCorr(const GeometricField< Type, fvPatchField, volMesh > &U, const GeometricField< Type, fvsPatchField, surfaceMesh > &Uf)
Definition: fvcDdt.C:155
tmp< GeometricField< Type, fvPatchField, volMesh > > laplacian(const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvcLaplacian.C:45
messageStream Info
tmp< surfaceScalarField > interpolate(const RhoType &rho)
tmp< GeometricField< Type, fvPatchField, volMesh > > ddt(const dimensioned< Type > dt, const fvMesh &mesh)
Definition: fvcDdt.C:45
dynamicFvMesh & mesh
fvMatrix< vector > fvVectorMatrix
Definition: fvMatricesFwd.H:45
const scalar pRefValue
fvMatrix< scalar > fvScalarMatrix
Definition: fvMatricesFwd.H:42
static const char nl
Definition: Ostream.H:260
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
fvVectorMatrix UEqn(fvm::ddt(rho, U)+fvm::div(phi, U)+turbulence->divDevRhoReff(U)==fvOptions(rho, U))
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition: volFieldsFwd.H:52
volScalarField & p
Definition: createFields.H:51
GeometricField< scalar, fvsPatchField, surfaceMesh > surfaceScalarField
surfaceScalarField & phi
GeometricField< vector, fvPatchField, volMesh > volVectorField
Definition: volFieldsFwd.H:55
volScalarField rAU(1.0/UEqn.A())
const label pRefCell
volScalarField & nu
tmp< GeometricField< Type, fvPatchField, volMesh > > div(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcDiv.C:47
tmp< GeometricField< typename outerProduct< vector, Type >::type, fvPatchField, volMesh >> grad(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcGrad.C:52
int main(int argc, char *argv[])
Definition: postCalc.C:54
surfaceScalarField & phiB
Definition: createPhiB.H:47
U
Definition: pEqn.H:82
HbyA
Definition: pEqn.H:7