applications/solvers/electromagnetics/mhdFoam/mhdFoam.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2011 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 "OSspecific.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     #include "createFields.H"
  63     #include "initContinuityErrs.H"
  64
  65     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  66
  67     Info<< nl << "Starting time loop" << endl;
  68
  69     while (runTime.loop())
  70     {
  71         #include "readPISOControls.H"
  72         #include "readBPISOControls.H"
  73
  74         Info<< "Time = " << runTime.timeName() << nl << endl;
  75
  76         #include "CourantNo.H"
  77
  78         {
  79             fvVectorMatrix UEqn
  80             (
  81                 fvm::ddt(U)
  82               + fvm::div(phi, U)
  83               - fvc::div(phiB, 2.0*DBU*B)
  84               - fvm::laplacian(nu, U)
  85               + fvc::grad(DBU*magSqr(B))
  86             );
  87
  88             solve(UEqn == -fvc::grad(p));
  89
  90
  91             // --- PISO loop
  92
  93             for (int corr=0; corr<nCorr; corr++)
  94             {
  95                 volScalarField rAU(1.0/UEqn.A());
  96
  97                 U = rAU*UEqn.H();
  98
  99                 phi = (fvc::interpolate(U) & mesh.Sf())
 100                     + fvc::ddtPhiCorr(rAU, U, phi);
 101
 102                 for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
 103                 {
 104                     fvScalarMatrix pEqn
 105                     (
 106                         fvm::laplacian(rAU, p) == fvc::div(phi)
 107                     );
 108
 109                     pEqn.setReference(pRefCell, pRefValue);
 110                     pEqn.solve();
 111
 112                     if (nonOrth == nNonOrthCorr)
 113                     {
 114                         phi -= pEqn.flux();
 115                     }
 116                 }
 117
 118                 #include "continuityErrs.H"
 119
 120                 U -= rAU*fvc::grad(p);
 121                 U.correctBoundaryConditions();
 122             }
 123         }
 124
 125         // --- B-PISO loop
 126
 127         for (int Bcorr=0; Bcorr<nBcorr; Bcorr++)
 128         {
 129             fvVectorMatrix BEqn
 130             (
 131                 fvm::ddt(B)
 132               + fvm::div(phi, B)
 133               - fvc::div(phiB, U)
 134               - fvm::laplacian(DB, B)
 135             );
 136
 137             BEqn.solve();
 138
 139             volScalarField rBA(1.0/BEqn.A());
 140
 141             phiB = (fvc::interpolate(B) & mesh.Sf())
 142                 + fvc::ddtPhiCorr(rBA, B, phiB);
 143
 144             fvScalarMatrix pBEqn
 145             (
 146                 fvm::laplacian(rBA, pB) == fvc::div(phiB)
 147             );
 148             pBEqn.solve();
 149
 150             phiB -= pBEqn.flux();
 151
 152             #include "magneticFieldErr.H"
 153         }
 154
 155         runTime.write();
 156     }
 157
 158     Info<< "End\n" << endl;
 159
 160     return 0;
 161 }
 162
 163
 164 // ************************************************************************* //