applications/solvers/electromagnetics/mhdFoam/mhdFoam.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright held by original author
   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 the
  13     Free Software Foundation; either version 2 of the License, or (at your
  14     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, write to the Free Software Foundation,
  23     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
  24
  25 Application
  26     mhdFoam
  27
  28 Description
  29     Solver for magnetohydrodynamics (MHD): incompressible, laminar flow of a
  30     conducting fluid under the influence of a magnetic field.
  31
  32     An applied magnetic field H acts as a driving force,
  33     at present boundary conditions cannot be set via the
  34     electric field E or current density J. The fluid viscosity nu,
  35     conductivity sigma and permeability mu are read in as uniform
  36     constants.
  37
  38     A fictitous magnetic flux pressure pH is introduced in order to
  39     compensate for discretisation errors and create a magnetic face flux
  40     field which is divergence free as required by Maxwell's equations.
  41
  42     However, in this formulation discretisation error prevents the normal
  43     stresses in UB from cancelling with those from BU, but it is unknown
  44     whether this is a serious error.  A correction could be introduced
  45     whereby the normal stresses in the discretised BU term are replaced
  46     by those from the UB term, but this would violate the boundedness
  47     constraint presently observed in the present numerics which
  48     guarantees div(U) and div(H) are zero.
  49
  50 \*---------------------------------------------------------------------------*/
  51
  52 #include "fvCFD.H"
  53 #include "OSspecific.H"
  54
  55 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  56
  57 int main(int argc, char *argv[])
  58 {
  59 #   include "setRootCase.H"
  60
  61 #   include "createTime.H"
  62 #   include "createMesh.H"
  63 #   include "createFields.H"
  64 #   include "initContinuityErrs.H"
  65
  66
  67 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  68
  69     Info<< nl << "Starting time loop" << endl;
  70
  71     while (runTime.loop())
  72     {
  73 #       include "readPISOControls.H"
  74 #       include "readBPISOControls.H"
  75
  76         Info<< "Time = " << runTime.timeName() << nl << endl;
  77
  78 #       include "CourantNo.H"
  79
  80         {
  81             fvVectorMatrix UEqn
  82             (
  83                 fvm::ddt(U)
  84               + fvm::div(phi, U)
  85               - fvc::div(phiB, 2.0*DBU*B)
  86               - fvm::laplacian(nu, U)
  87               + fvc::grad(DBU*magSqr(B))
  88             );
  89
  90             solve(UEqn == -fvc::grad(p));
  91
  92
  93             // --- PISO loop
  94
  95             for (int corr = 0; corr < nCorr; corr++)
  96             {
  97                 volScalarField rUA = 1.0/UEqn.A();
  98
  99                 U = rUA*UEqn.H();
 100
 101                 phi = (fvc::interpolate(U) & mesh.Sf())
 102                     + fvc::ddtPhiCorr(rUA, U, phi);
 103
 104                 for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
 105                 {
 106                     fvScalarMatrix pEqn
 107                     (
 108                         fvm::laplacian(rUA, p) == fvc::div(phi)
 109                     );
 110
 111                     pEqn.setReference(pRefCell, pRefValue);
 112                     pEqn.solve();
 113
 114                     if (nonOrth == nNonOrthCorr)
 115                     {
 116                         phi -= pEqn.flux();
 117                     }
 118                 }
 119
 120 #               include "continuityErrs.H"
 121
 122                 U -= rUA*fvc::grad(p);
 123                 U.correctBoundaryConditions();
 124             }
 125         }
 126
 127         // --- B-PISO loop
 128
 129         for (int Bcorr=0; Bcorr<nBcorr; Bcorr++)
 130         {
 131             fvVectorMatrix BEqn
 132             (
 133                 fvm::ddt(B)
 134               + fvm::div(phi, B)
 135               - fvc::div(phiB, U)
 136               - fvm::laplacian(DB, B)
 137             );
 138
 139             BEqn.solve();
 140
 141             volScalarField rBA = 1.0/BEqn.A();
 142
 143             phiB = (fvc::interpolate(B) & mesh.Sf())
 144                 + fvc::ddtPhiCorr(rBA, B, phiB);
 145
 146             fvScalarMatrix pBEqn
 147             (
 148                 fvm::laplacian(rBA, pB) == fvc::div(phiB)
 149             );
 150             pBEqn.solve();
 151
 152             phiB -= pBEqn.flux();
 153
 154 #           include "magneticFieldErr.H"
 155         }
 156
 157         runTime.write();
 158     }
 159
 160     Info<< "End\n" << endl;
 161
 162     return 0;
 163 }
 164
 165
 166 // ************************************************************************* //