mhdFoam.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2023 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 fictitious magnetic flux pressure pB 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 "argList.H"
52 #include "timeSelector.H"
53 #include "pisoControl.H"
54 #include "pressureReference.H"
55 #include "constrainPressure.H"
56 #include "constrainHbyA.H"
57 
58 #include "fvcDdt.H"
59 #include "fvcGrad.H"
60 #include "fvcFlux.H"
61 
62 #include "fvmDdt.H"
63 #include "fvmDiv.H"
64 #include "fvmLaplacian.H"
65 
66 using namespace Foam;
67 
68 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
69 
70 int main(int argc, char *argv[])
71 {
72  #include "postProcess.H"
73 
74  #include "setRootCase.H"
75  #include "createTime.H"
76  #include "createMesh.H"
77  #include "createControl.H"
78  #include "createFields.H"
79  #include "initContinuityErrs.H"
80 
81  // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
82 
83  Info<< nl << "Starting time loop" << endl;
84 
85  while (runTime.loop())
86  {
87  Info<< "Time = " << runTime.userTimeName() << nl << endl;
88 
89  #include "CourantNo.H"
90 
91  {
93  (
94  fvm::ddt(U)
95  + fvm::div(phi, U)
96  - fvc::div(phiB, 2.0*DBU*B)
97  - fvm::laplacian(nu, U)
98  + fvc::grad(DBU*magSqr(B))
99  );
100 
101  if (piso.momentumPredictor())
102  {
103  solve(UEqn == -fvc::grad(p));
104  }
105 
106 
107  // --- PISO loop
108  while (piso.correct())
109  {
110  volScalarField rAU(1.0/UEqn.A());
111  surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
114  (
115  "phiHbyA",
116  fvc::flux(HbyA)
117  + rAUf*fvc::ddtCorr(U, phi)
118  );
119 
120  // Update the pressure BCs to ensure flux consistency
121  constrainPressure(p, U, phiHbyA, rAUf);
122 
123  while (piso.correctNonOrthogonal())
124  {
125  fvScalarMatrix pEqn
126  (
127  fvm::laplacian(rAUf, p) == fvc::div(phiHbyA)
128  );
129 
130  pEqn.setReference
131  (
134  );
135 
136  pEqn.solve();
137 
139  {
140  phi = phiHbyA - pEqn.flux();
141  }
142  }
143 
144  #include "continuityErrs.H"
145 
146  U = HbyA - rAU*fvc::grad(p);
147  U.correctBoundaryConditions();
148  }
149  }
150 
151  // --- B-PISO loop
152  while (bpiso.correct())
153  {
154  fvVectorMatrix BEqn
155  (
156  fvm::ddt(B)
157  + fvm::div(phi, B)
158  - fvc::div(phiB, U)
159  - fvm::laplacian(DB, B)
160  );
161 
162  BEqn.solve();
163 
164  volScalarField rAB(1.0/BEqn.A());
165  surfaceScalarField rABf("rABf", fvc::interpolate(rAB));
166 
167  phiB = fvc::flux(B);
168 
169  while (bpiso.correctNonOrthogonal())
170  {
171  fvScalarMatrix pBEqn
172  (
173  fvm::laplacian(rABf, pB) == fvc::div(phiB)
174  );
175 
176  pBEqn.solve();
177 
179  {
180  phiB -= pBEqn.flux();
181  }
182  }
183 
184  #include "magneticFieldErr.H"
185  }
186 
187  runTime.write();
188  }
189 
190  Info<< "End\n" << endl;
191 
192  return 0;
193 }
194 
195 
196 // ************************************************************************* //
static const Foam::dimensionedScalar B("B", Foam::dimless, 18.678)
pisoControl bpiso(mesh, "BPISO")
Generic GeometricField class.
bool momentumPredictor() const
Flag to indicate to solve for momentum.
A special matrix type and solver, designed for finite volume solutions of scalar equations....
Definition: fvMatrix.H:118
SolverPerformance< Type > solve(const dictionary &)
Solve segregated or coupled returning the solution statistics.
Definition: fvMatrixSolve.C:58
tmp< volScalarField > A() const
Return the central coefficient.
Definition: fvMatrix.C:808
tmp< VolField< Type > > H() const
Return the H operation source.
Definition: fvMatrix.C:868
tmp< SurfaceField< Type > > flux() const
Return the face-flux field from the matrix.
Definition: fvMatrix.C:964
void setReference(const label celli, const Type &value, const bool forceReference=false)
Set reference level for solution.
Definition: fvMatrix.C:563
bool finalNonOrthogonalIter() const
Flag to indicate the last non-orthogonal iteration.
bool correct(const bool finalIter=false)
Piso loop within outer loop.
Definition: pisoControl.C:79
bool correctNonOrthogonal(const bool finalIter=true)
Non-orthogonal corrector loop.
Definition: pisoControl.C:98
Provides controls for the pressure reference in closed-volume simulations.
scalar refValue() const
Return the pressure reference level.
label refCell() const
Return the cell in which the reference pressure is set.
fvVectorMatrix & UEqn
Definition: UEqn.H:11
Calculates and prints the continuity errors.
surfaceScalarField & phiB
Definition: createPhiB.H:47
pisoControl piso(mesh)
Calculate the first temporal derivative.
Calculate the face-flux of the given field.
Calculate the gradient of the given field.
Calculate the matrix for the first temporal derivative.
Calculate the matrix for the divergence of the given field and flux.
Calculate the matrix for the laplacian of the field.
Declare and initialise the cumulative continuity error.
volScalarField rAU(1.0/UEqn.A())
U
Definition: pEqn.H:72
volVectorField & HbyA
Definition: pEqn.H:13
surfaceScalarField phiHbyA("phiHbyA", fvc::interpolate(rho) *fvc::flux(HbyA))
int main(int argc, char *argv[])
Definition: mhdFoam.C:67
tmp< SurfaceField< typename innerProduct< vector, Type >::type > > flux(const VolField< Type > &vf)
Return the face-flux field obtained from the given volVectorField.
static tmp< SurfaceField< Type > > interpolate(const VolField< Type > &tvf, const surfaceScalarField &faceFlux, Istream &schemeData)
Interpolate field onto faces using scheme given by Istream.
tmp< VolField< typename outerProduct< vector, Type >::type > > grad(const SurfaceField< Type > &ssf)
Definition: fvcGrad.C:46
tmp< VolField< Type > > div(const SurfaceField< Type > &ssf)
Definition: fvcDiv.C:47
tmp< SurfaceField< typename Foam::flux< Type >::type > > ddtCorr(const VolField< Type > &U, const SurfaceField< Type > &Uf)
Definition: fvcDdt.C:196
tmp< fvMatrix< Type > > laplacian(const VolField< Type > &vf, const word &name)
Definition: fvmLaplacian.C:47
tmp< fvMatrix< Type > > div(const surfaceScalarField &flux, const VolField< Type > &vf, const word &name)
Definition: fvmDiv.C:48
tmp< fvMatrix< Type > > ddt(const VolField< Type > &vf)
Definition: fvmDdt.C:46
Namespace for OpenFOAM.
void constrainPressure(volScalarField &p, const RhoType &rho, const volVectorField &U, const surfaceScalarField &phiHbyA, const RAUType &rhorAU, const MRFType &MRF)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:257
messageStream Info
void constrainHbyA(volVectorField &HbyA, const volVectorField &U, const volScalarField &p)
Definition: constrainHbyA.C:34
static const char nl
Definition: Ostream.H:266
dimensioned< scalar > magSqr(const dimensioned< Type > &)
SolverPerformance< Type > solve(fvMatrix< Type > &, const word &)
Solve returning the solution statistics given convergence tolerance.
Execute application functionObjects to post-process existing results.
volScalarField & p