FreeFOAM The Cross-Platform CFD Toolkit
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"
82 
83  // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
84 
85  Info<< nl << "Starting time loop" << endl;
86 
87  while (runTime.loop())
88  {
90  #include "readBPISOControls.H"
91 
92  Info<< "Time = " << runTime.timeName() << nl << endl;
93 
94  #include <finiteVolume/CourantNo.H>
95 
96  {
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 
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: ************************ //