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-2016 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 "postProcess.H"
59 
60  #include "setRootCase.H"
61  #include "createTime.H"
62  #include "createMesh.H"
63  #include "createControl.H"
64  #include "createFields.H"
65  #include "initContinuityErrs.H"
66 
67  // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
68 
69  Info<< nl << "Starting time loop" << endl;
70 
71  while (runTime.loop())
72  {
73  Info<< "Time = " << runTime.timeName() << nl << endl;
74 
75  #include "CourantNo.H"
76 
77  {
79  (
80  fvm::ddt(U)
81  + fvm::div(phi, U)
82  - fvc::div(phiB, 2.0*DBU*B)
83  - fvm::laplacian(nu, U)
84  + fvc::grad(DBU*magSqr(B))
85  );
86 
87  if (piso.momentumPredictor())
88  {
89  solve(UEqn == -fvc::grad(p));
90  }
91 
92 
93  // --- PISO loop
94  while (piso.correct())
95  {
96  volScalarField rAU(1.0/UEqn.A());
100  (
101  "phiHbyA",
102  fvc::flux(HbyA)
103  + rAUf*fvc::ddtCorr(U, phi)
104  );
105 
106  // Update the pressure BCs to ensure flux consistency
108 
109  while (piso.correctNonOrthogonal())
110  {
111  fvScalarMatrix pEqn
112  (
114  );
115 
116  pEqn.setReference(pRefCell, pRefValue);
117  pEqn.solve(mesh.solver(p.select(piso.finalInnerIter())));
118 
119  if (piso.finalNonOrthogonalIter())
120  {
121  phi = phiHbyA - pEqn.flux();
122  }
123  }
124 
125  #include "continuityErrs.H"
126 
127  U = HbyA - rAU*fvc::grad(p);
128  U.correctBoundaryConditions();
129  }
130  }
131 
132  // --- B-PISO loop
133  while (bpiso.correct())
134  {
135  fvVectorMatrix BEqn
136  (
137  fvm::ddt(B)
138  + fvm::div(phi, B)
139  - fvc::div(phiB, U)
140  - fvm::laplacian(DB, B)
141  );
142 
143  BEqn.solve();
144 
145  volScalarField rAB(1.0/BEqn.A());
146  surfaceScalarField rABf("rABf", fvc::interpolate(rAB));
147 
148  phiB = fvc::flux(B) + rABf*fvc::ddtCorr(B, phiB);
149 
150  while (bpiso.correctNonOrthogonal())
151  {
152  fvScalarMatrix pBEqn
153  (
154  fvm::laplacian(rABf, pB) == fvc::div(phiB)
155  );
156 
157  pBEqn.solve(mesh.solver(pB.select(bpiso.finalInnerIter())));
158 
159  if (bpiso.finalNonOrthogonalIter())
160  {
161  phiB -= pBEqn.flux();
162  }
163  }
164 
165  #include "magneticFieldErr.H"
166  }
167 
168  runTime.write();
169  }
170 
171  Info<< "End\n" << endl;
172 
173  return 0;
174 }
175 
176 
177 // ************************************************************************* //
surfaceScalarField & phi
tmp< GeometricField< typename outerProduct< vector, Type >::type, fvPatchField, volMesh >> grad(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcGrad.C:52
fvMatrix< scalar > fvScalarMatrix
Definition: fvMatricesFwd.H:42
U
Definition: pEqn.H:83
tmp< GeometricField< Type, fvPatchField, volMesh > > div(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcDiv.C:47
void constrainPressure(volScalarField &p, const RhoType &rho, const volVectorField &U, const surfaceScalarField &phiHbyA, const RAUType &rhorAU, const MRFType &MRF)
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
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:253
pisoControl piso(mesh)
tmp< surfaceScalarField > interpolate(const RhoType &rho)
GeometricField< vector, fvPatchField, volMesh > volVectorField
Definition: volFieldsFwd.H:55
tmp< GeometricField< Type, fvPatchField, volMesh > > laplacian(const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvcLaplacian.C:45
HbyA
Definition: pcEqn.H:74
pisoControl bpiso(mesh,"BPISO")
phiHbyA
Definition: pcEqn.H:73
tmp< GeometricField< Type, fvPatchField, volMesh > > ddt(const dimensioned< Type > dt, const fvMesh &mesh)
Definition: fvcDdt.C:45
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition: volFieldsFwd.H:52
dynamicFvMesh & mesh
rhoEqn solve()
tmp< volVectorField > constrainHbyA(const tmp< volVectorField > &tHbyA, const volVectorField &U, const volScalarField &p)
Definition: constrainHbyA.C:33
dimensioned< scalar > magSqr(const dimensioned< Type > &)
static const char nl
Definition: Ostream.H:262
surfaceScalarField & phiB
Definition: createPhiB.H:47
const scalar pRefValue
const label pRefCell
surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU))
fvVectorMatrix & UEqn
Definition: UEqn.H:13
fvMatrix< vector > fvVectorMatrix
Definition: fvMatricesFwd.H:45
messageStream Info
volScalarField & p
tmp< surfaceScalarField > flux(const volVectorField &vvf)
Return the face-flux field obtained from the given volVectorField.
Definition: fvcFlux.C:32
GeometricField< scalar, fvsPatchField, surfaceMesh > surfaceScalarField
Execute application functionObjects to post-process existing results.
volScalarField & nu
volScalarField rAU(1.0/UEqn.A())