A spectral/hp element MHD solver
Alexander V. Proskurin Anatoly M. Sagalakov
physics.comp-ph


Analysis of Deformation Fields in Spatio-temporal CBCT images of lungs for radiotherapy patients
    Sriram Prasath
    Received: date / Accepted: date
================================================================================================





A new MHD solver, based on the  spectral/hp element framework, is presented in this paper. The velocity and electric potential quasi-static MHD model is used. The Hartmann flow in plane channel and its stability, the Hartmann flow in rectangular duct, and the stability of Hunt's flow are explored as examples.  Exponential convergence is achieved and the resulting numerical values were found to have an accuracy up to 10^-12 for the state flows compared to an exact solution, and 10^-5 for the stability eigenvalues compared to independent numerical results.






§ INTRODUCTION


It is well known that high-order methods have good computational properties, fast convergence, small errors, and the most compact data representation. For many problems in hydrodynamics, high-order methods are necessary. Such problems include the time-dependent simulation of transient flow regimes and the investigation of hydrodynamic stability. Of course, turbulent flows can be investigated using low accuracy schemes (achieving low precision results), but in most cases problems in channels of hydrodynamic stability  require the use of spectral methods. The classical Orr-Sommerfield equation has a small parameter 1/Re at the highest derivative which causes rapidly oscillating solutions. The first numerical calculation of eigenvalues for this equation <cit.> used a high-order finite difference scheme. Later, Orszag <cit.>  achieved more accurate results and explained the convenience of using high-order methods for problems of flow stability. A recent review of flow stability in complex geometries and the advantages and disadvantages of high-order methods can be found in <cit.>.



§ PROBLEM FORMULATION


Consider a flow of incompressible, electrically conducting liquid in the presence of an imposed magnetic field. We suppose that Re_m ≪ Re. In this case a magnetic field generated by the fluid movement does not affect the imposed magnetic field. This is correct for most engineering applications <cit.>.  It is now possible to write the Navier-Stokes equation in the form

    ∂v/∂t+( v∇)v   = -1/ρ∇ p + νΔv + F(v,H),
    
    div v   = 0,

where v is the velocity, p is the pressure, ν is the viscosity, ρ is the density, F is the magnetic force, and H is the magnitude of the imposed magnetic field.

Ohm's law is:

    j = σ( -∇φ+v×H),

where j is the density of electirc current, φ is the electric potential, and σ is the conductivity. Using the law of conservation of electric charge (div j = 0), it is possible to derive the equation for electric potential as:

    Δφ = ∇(v×H).


The system (<ref>) can be written in the form: 

    ∂v/∂t+( v∇)v   = -∇ p + 1/ReΔv + St ( -∇φ+v×H)×H ,
       div v = 0,
       Δφ = ∇(v×H),

where Re = L_0V_0/ν is the Reynolds number, St=σ H_0^2 L_0/ρ V_0 is the magnetic interaction parameter (Stuart number), and L_0, V_0, and H_0 represent the scales of length, velocity and magnitude of the imposed magnetic field, respectively. The system (<ref>) is also known as the MHD system in quasi-static approximation in electric potential form. This system is widely used in theoretical investigations and accurately approximate many cases of liquid metal flows (for example, see the appropriate discussion and reference in <cit.>).

The boundary condition for velocity at walls have the form:

    v = 0,

and the boundary condition for the electric potential at perfectly conducting walls is:

    φ = const.

The boundary condition for insulating walls is:

    ∂φ/∂n = 0.




§ NUMERICAL TECHNIQUE OVERVIEW


Our new MHD solver has been developed on the basis of an open source spectral/hp element framework  <cit.>.
The incompressible Navier-Stokes solver () from the framework has been taken as the source for the MHD solver.

 uses the velocity correction scheme as described in <cit.>, assuming the time grid t_0, t_1, …,t_n-1,t_n,t_n+1. Using a first-order difference scheme, it is possible to define intermediate velocity ṽ by the equation:

    ṽ-v_n/δ t = -(v_n ∇)v_n + F_n+St · v_n×( -∇φ_n+v_n×H),

where F is a force acting of a fluid. At this stage, the  software allows force to be imposed which acts on the flow. The electrical potential is found by solving the equation Δφ_n = ∇(v_n×H) and, after this, the magnetic force in (<ref>) is calculated.
We define the second intermediate velocity as:

    ṽ̂-ṽ/δ t = - ∇ p_n+1.

The Poisson equation

    Δ p_n+1 = ∇(ṽ/δ t)

is immediately derived using div ṽ̂ = 0. Thus, at this stage, the divergence-free condition is approximately satisfied. The boundary conditions for pressure are discussed in <cit.>. The last step of the velocity correction procedure is the equation:

    (Δ - Re/δ t)v_n+1 = - Re/δ tṽ + Re ∇ p_n+1,

wich allows us to find the next time-step velocity v_n+1.  can use first, second and third order schemes.




§ THE HARTMANN FLOW IN A PLANE CHANNEL


Figure <ref> illustrates a flow in a plane channel. Two parallel infinite planes are installed at points y = ± 1. The liquid between the planes flows under a constant pressure gradient ∇ p in a x direction. The magnetic field H is perpendicular to the planes. This is the Hartmann flow in the plane channel. According to <cit.> the solution of (<ref>) is 

    u(y)/u(0)=cosh (M  ) - cosh (M y  )/cosh (M  ) - 1

where M = √(St · Re), u(0) is a centreline velocity. Velocity graphs are shown in Figure <ref> for M=10,100,10000.

Figure <ref> shows a 2D mesh for numeric calculations of the flow. In Figure <ref> one can see large gradients of velocity near the walls in cases of large M and we should take special attention to this part of the flow. It is possible to increase accuracy by mesh concentration near the walls where there are large velocity gradients. This is the h-type solution refinement. The high-order method can increase accuracy by increasing the polynomials' order of an approximation, this is p-refinement. For the flow calculations we will combine both methods by using the order of approximation p from 5 to 25 and mesh condensation near the walls with a coefficient β (β=1 for a uniform grid).  

When it is supposed that the flow is two-dimensional, all functions are independent from z-coordinate and v_z=0. This proposition leads to the equation Δφ = 0,  this means that boundary conditions for φ are not required.
The boundary conditions for velocity and pressure are  

    v   = 0,  ∂ p/∂ n = 0  at walls,
    ∂v/∂ n   = 0, ∂ p/∂ n = 0  at inflow and outflow.


In Table <ref>, maximum deviations from the exact solution are presented at M=10∼ 10^4 for different orders of polynomial approximation p. The state flow (<ref>) is calculated as a time-dependent flow with zero initial conditions. Table <ref> includes the running time of the solver on an AMD Ryzen Threadripper 1920X machine with 12 threads. 



§ FLOW IN RECTANGULAR DUCT

Consider a steady flow in a rectangular duct. A sketch of this flow is shown in Figure <ref>. The rectangle in Figure <ref> is the cross section of the channel and a uniform magnetic field is applied vertically. Liquid flows under a constant pressure gradient, perpendicular to the plane of the diagram. This flow was investigated analytically in <cit.>.

For flow computations in the rectangular duct, the authors use a mesh shown in Figure <ref>.  software allows us to set up 3D problems where the homogeneous direction is z, the number of Fourier modes is 2(the minimal possible value).

The velocity convergence at points (0.95,0.0) for M=10^3 and (0.98,0.0) for M=10^4 is presented in Table <ref> for the case of perfectly electro-conducting walls. The table includes the mesh concentration coefficient β and velocity values for different p from 5 to 25, the time of calculation and memory usage, for an AMD Ryzen Threadripper 1920X machine with 12 threads. The velocity graph for M=10^3, M=10^4 is shown in Figure <ref>.




§ THE STABILITY PROBLEM

For a stability analysis let us decompose velocity, pressure and electric potential to form

    v   = U+v,
    φ   = φ_0+φ,
    
    p    = p_0 +p,

where U, φ_0, and p_0 are values of a steady-state flow and v, φ and p are small disturbances. The system (<ref>) becomes linear:

    ∂v/∂t+( U∇)v+( v∇)U   = -∇ p + 1/ReΔv + St ( -∇φ+v×H)×H ,
       div v = 0,
       Δφ = ∇(v×H).


From equations (<ref>) a linear operator A can be constructed:

    v(x,y,z,T) = A(T)v(x,y,z,0) = λ(T) v(x,y,z,0),

where T is a time interval.
The linear operator A(T) is constructed numerically by the splitting procedure in the same way as in the case of  nonlinear equations (<ref>). In order to find an eigenvalue λ(T), it is convenient to construct a Krylov subspace

    K_n(A,v_0) = span{v_0, A(T) v_0,  A(T)^2 v_0, …,  A(T)^n-1v_0 },

where A(T)^i v_0 is obtained by direct calculation v_1 = A(T)v_0, v_2 = A(T)v_1, …. Further eigenvalue calculations are carried out by standard numerical algebraic techniques, such as the Arnoldi method. The eigenvalues are obtained by  in the form:

    λ(T) = m · e^θ i,

and if m>1 then the flow is unstable. The time-independent growth is σ = ln(m)/T and the time-independent frequency is ω=θ/T.




§ STABILITY OF THE HARTMANN FLOW

In Section <ref>, the Hartmann flow was considered. Now, we will explore the stability of this flow. We take a Hartmann 2D flow disturbance (<ref>) in the form:

    {v,φ,p } = q(x,y) e^iα(x-C t),

where q(x,y)  is the amplitude of disturbance, α=2π/γ is the wave-vector, γ is the wavelength, C=X+iY is the phase velocity of disturbance, α X=ω is the frequency, α Y = σ is the growth of disturbance. When σ≤ 0, it means that the flow is stable. 

The disturbance form (<ref>) is widely used in hydrodynamics stability analysis and leads to the eigenvalue problem, equivalent to (<ref>). As reference data, we take Takashima critical values <cit.> for 2D disturbances.  
In Table  <ref>, growths and frequencies of the disturbances are given for several cases. The values of Re and M 
are taken from the article <cit.>, and these are critical values of Hartmann flow. The computational grid is shown in Figure <ref>; nx and ny are the number of cells in the horizontal and vertical directions. The length L of the grid is set up by using a critical wave-vector  α_c = 2π/L. Boundary conditions at inlet and outlet are periodical. In Table <ref> complete coincidence is observed with the reference data from <cit.> and convergence of the eigenvalues is achieved up to 10^-7. Additionally, the table includes the time and memory usage (for the AMD Phenom FX-8150 processor with 8 threads).



§ STABILITY OF HUNT'S FLOW

Consider a steady flow in the rectangular duct (Figure <ref>), where horizontal walls are perfectly electrically conducting and vertical walls are perfectly electrically insulating. The flow was investigated in <cit.> and it is known as the Hunt's flow. A mesh for base flow and stability calculations is shown at Figure <ref>. Figure <ref> presents a graph of velocity over a line y=0 at M=10, 10^2 10^3. In Table <ref> our calculated eigenvalues are compared with reference values from article <cit.>; the time of calculation and the memory usage are presented. It is possible to see numerical convergence by increasing the order p and match with the reference values from <cit.> up to 10^-5, excluding the case M=10^3.



§ CONCLUSION

This article presents the spectral/hp element solver for MHD problems based on the  framework. The solver also makes it possible to investigate the stability of such flows. In order to demonstrate the solver's capacity, several examples were considered: the Hartmann flow in a plane channel and its stability, the Hartmann flow in a rectangular duct, and the stability of 
Hunt's flow. For the flows, it is easy to find steady-state solutions analytically, and these results were used as the reference test solutions. It was found that the margin of error decreases exponentially with an increasing degree of approximating polynomials, an accuracy 10^-12 can be achieved. To estimate the costs of computer time and memory, these data were listed in the tables for several cases. The computational costs for the stability calculations are large. The first reason for this is the fact that a non-stationary algorithm was used, which allowed use of the non-stationary solver with small adaptations. To obtain eigenvalues with high accuracy, we should set a small time step. The second reason is that we are considering the test examples in 2D for the Hartmann flow and 3D for flows in duct. Usually, these problems can be reduced to the more simple cases described in <cit.>, which can be investigated with much lower computational costs. In this article we demonstrated the accuracy of the method using the well investigated examples. 
In general, our numerical technique is intended for complex geometry flows where such simplifications are not possible. 


mhd










































