2.3 Magnetohydrodynamic flow of a liquid

Tutorial path:

In this example we shall investigate the flow of an electrically-conducting liquid through a magnetic field. The problem belongs to the branch of fluid dynamics known as magnetohydrodynamics (MHD), simulated using the mhdFoam solver.

2.3.1 Problem specification

This case is known as the Hartmann problem, chosen as it contains an analytical solution with which mhdFoam can be validated. It is defined as follows:

Solution domain
The domain is 2 dimensional and consists of flow along two parallel plates as shown in Fig. 2.18.
\relax \special {t4ht=

Figure 2.18: Geometry of the Hartmann problem


Governing equations
 
  • Mass continuity for incompressible fluid
    ∇ ∙U =  0
           \relax \special {t4ht=
    (2.16)

  • Momentum equation for incompressible fluid
    ∂U
----+ ∇ ∙(UU  ) + ∇ ∙(2BΓ BUB  ) + ∇ ∙(νU ) + ∇ (Γ BUB ∙∙B ) = − ∇p
∂t
           \relax \special {t4ht=
    (2.17)

    where B  \relax \special {t4ht= is the magnetic flux density,             −1
Γ BU = (2μρ)   \relax \special {t4ht=.

  • Maxwell’s equations
    ∇  × E =  − ∂B--
            ∂t
           \relax \special {t4ht=
    (2.18)

    where E  \relax \special {t4ht= is the electric field strength.

    ∇ ∙B =  0
           \relax \special {t4ht=
    (2.19)

                  ∂D
∇ × H  = J +  ----= J
              ∂t
           \relax \special {t4ht=
    (2.20)

    assuming ∂D  ∕∂t ≪  J  \relax \special {t4ht=. Here, H  \relax \special {t4ht= is the magnetic field strength, J  \relax \special {t4ht= is the current density and D  \relax \special {t4ht= is the electric flux density.

  • Charge continuity
    ∇ ∙J = 0
           \relax \special {t4ht=
    (2.21)

  • Constitutive law
    B  = μH
           \relax \special {t4ht=
    (2.22)

  • Ohm’s law
    J = σ (E + U  × B )
           \relax \special {t4ht=
    (2.23)

  • Combining Equation 2.18, Equation 2.20, Equation 2.23, and taking the curl
    ∂B
----+ ∇ ∙(UB  ) − ∇ ∙(ϕBU ) − ∇ ∙(Γ BB ) = 0
 ∂t
           \relax \special {t4ht=
    (2.24)

Boundary conditions
 
  • inlet is specified the inlet condition with fixed velocity U =  (1,0,0)  \relax \special {t4ht= m/s;
  • outlet is specified as the outlet with with fixed pressure p = 0 Pa  \relax \special {t4ht=;
  • upperWall is specified as a wall where B =  (0, 20,0) T  \relax \special {t4ht=.
  • lowerWall is specified as a wall where B =  (0, 20,0) T  \relax \special {t4ht=.
  • front and back boundaries are specified as empty.
Initial conditions
U  = 0 m ∕s  \relax \special {t4ht=, p = 100 Pa  \relax \special {t4ht=, B =  (0, 20,0) T  \relax \special {t4ht=.
Transport properties
 
  • Kinematic viscosity ν = 1 Pas  \relax \special {t4ht=
  • Density ρ = 1 kgm ∕s  \relax \special {t4ht=
  • Electrical conductivity σ = 1 (Ωm )−1   \relax \special {t4ht=
  • Permeability μ = 1 H ∕m  \relax \special {t4ht=
Solver name
mhdFoam: an incompressible laminar magneto-hydrodynamics code.
Case name
hartmann case located in the $FOAM_TUTORIALS/electromagnetics/mhdFoam directory.

2.3.2 Mesh generation

The geometry is simply modelled with 100 cells in the x  \relax \special {t4ht=-direction and 40 cells in the y  \relax \special {t4ht=-direction; the set of vertices and blocks are given in the mesh description file below:

1/*--------------------------------*- C++ -*----------------------------------*\
2| =========                 |                                                 |
3| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
4|  \\    /   O peration     | Version:  v2006                                 |
5|   \\  /    A nd           | Website:  www.openfoam.com                      |
6|    \\/     M anipulation  |                                                 |
7\*---------------------------------------------------------------------------*/
8FoamFile
9{
10    version     2.0;
11    format      ascii;
12    class       dictionary;
13    object      blockMeshDict;
14}
15// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
16
17scale   1;
18
19vertices
20(
21    (0 -1 0)
22    (20 -1 0)
23    (20 1 0)
24    (0 1 0)
25    (0 -1 0.1)
26    (20 -1 0.1)
27    (20 1 0.1)
28    (0 1 0.1)
29);
30
31blocks
32(
33    hex (0 1 2 3 4 5 6 7) (100 40 1) simpleGrading (1 1 1)
34);
35
36edges
37(
38);
39
40boundary
41(
42    inlet
43    {
44        type patch;
45        faces
46        (
47            (0 4 7 3)
48        );
49    }
50    outlet
51    {
52        type patch;
53        faces
54        (
55            (2 6 5 1)
56        );
57    }
58    lowerWall
59    {
60        type patch;
61        faces
62        (
63            (1 5 4 0)
64        );
65    }
66    upperWall
67    {
68        type patch;
69        faces
70        (
71            (3 7 6 2)
72        );
73    }
74    frontAndBack
75    {
76        type empty;
77        faces
78        (
79            (0 3 2 1)
80            (4 5 6 7)
81        );
82    }
83);
84
85mergePatchPairs
86(
87);
88
89// ************************************************************************* //

2.3.3 Running the case

The user can run the case and view results in ParaView. It is also useful at this stage to run the Ucomponents utility to convert the U  \relax \special {t4ht= vector field into individual scalar components. MHD flow is governed by, amongst other things, the Hartmann number which is a measure of the ratio of electromagnetic body force to viscous force

         ∘ ---
M  = BL    -σ-
           ρν
\relax \special {t4ht=
(2.25)

where L  \relax \special {t4ht= is the characteristic length scale. In this case with By  = 20 T  \relax \special {t4ht=, M  =  20  \relax \special {t4ht= and the electromagnetic body forces dominate the viscous forces. Consequently with the flow fairly steady at t = 2 s  \relax \special {t4ht= the velocity profile is almost planar, viewed at a cross section midway along the domain x =  10 m  \relax \special {t4ht=. The user can plot a graph of the profile of Ux  \relax \special {t4ht= in dxFoam.


\relax \special {t4ht=


Figure 2.19: Velocity profile in the Hartmann problem for By  = 1 T  \relax \special {t4ht= and By  = 20 T  \relax \special {t4ht=.


Now the user should reduce the magnetic flux density B  \relax \special {t4ht= to 1 T  \relax \special {t4ht=and re-run the code and Ucomponents. In this case, M  = 1  \relax \special {t4ht= and the electromagnetic body forces no longer dominate. The velocity profile consequently takes on the parabolic form, characteristic of Poiseuille flow as shown in Figure 2.19. To validate the code the analytical solution for the velocity profile U
  x  \relax \special {t4ht= is superimposed in Figure 2.19, given by:

U  (y)   cosh M  − cosh M (y∕L )
--x--- = -----------------------
Ux (0)         cosh M  − 1
\relax \special {t4ht=
(2.26)

where the characteristic length L  \relax \special {t4ht= is half the width of the domain, i.e. 1 m  \relax \special {t4ht=.