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