A consistent and conservative scheme for MHD flows with complex boundaries on an unstructured Cartesian adaptive system

Accepted Manuscript A consistent and conservative scheme for MHD flows with complex boundaries on an unstructured cartesian adaptive system Jie Zhang, Ming-Jiu Ni

PII: DOI: Reference:

S0021-9991(13)00531-7 10.1016/j.jcp.2013.08.004 YJCPH 4751

To appear in:

Journal of Computational Physics

Received date: 22 April 2013 Revised date: 20 July 2013 Accepted date: 2 August 2013

Please cite this article in press as: J. Zhang, M.-J. Ni, A consistent and conservative scheme for MHD flows with complex boundaries on an unstructured cartesian adaptive system, Journal of Computational Physics (2013), http://dx.doi.org/10.1016/j.jcp.2013.08.004

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

A Consistent and Conservative Scheme for MHD Flows with Complex Boundaries on an Unstructured Cartesian Adaptive System Jie Zhang, Ming-Jiu Ni* School of Physics, University of Chinese Academy of Sciences, Beijing 100049, China

Abstract The numerical simulation of MHD flows with complex boundaries has been a topic of great interest in the development of a fusion reactor blanket for the difficulty to accurately simulate the Hartmann layers and side layers along arbitrary geometries. An adaptive version of a consistent and conservative scheme has been developed for simulating the magnetohydrodynamics (MHD) flows. Besides, the present study forms the first attempt to apply the cut-cell approach for irregular wall-bounded MHD flows, which is more flexible and conveniently implemented under adaptive mesh refinement(AMR) technique. It employs a Volume-of-Fluid(VOF) approach to represent the fluid-conducting wall interface that makes it possible to solve the fluid-solid coupling magnetic problems, emphasizing at how electric field solver is implemented when conductivity is discontinuous in cut-cell. For the irregular cut-cells, the conservative interpolation technique is applied to calculate the Lorentz force at cell-center. On the other hand, it will be shown how consistent and conservative scheme is implemented on fine/coarse mesh boundaries when using AMR technique. Then, the applied numerical schemes are validated by five test simulations and excellent agreement was obtained for all the cases considered, simultaneously showed good consistency and conservative Email address: [email protected] (Ming-Jiu Ni*)

Preprint submitted to Journal of Computational Physics

September 10, 2013

properties. Keywords: MHD flow, complex boundary, cut-cell, solid-fluid coupling, consistent and conservative scheme 1. Introduction MHD describes the motion of liquids subjected to magnetic fields that the flow state and the current density distribution interact with each other. In a fusion reactor, MHD effect is of great significance that under a strong magnetic field, the conducting fluid will present different characters from ordinary fluid flow[23, 20], such as higher pressure drop, variant velocity distribution pattern, heat transfer enhancement or decline. In more complex cases, liquid flows in domain of intricate shapes are always practical and theoretical interest that the ability to model this complex wall-bounded MHD flows are crucial issues in the design of fusion device. For instance, loopshaped channel flows and sudden-expanding channel flows are both typical flow phenomena in the self-cooled liquid metal fusion blanket[21]. One of the main challenges in numerically simulating such models is handling of the complex boundaries. Traditionally, unstructured grids or boundary-fitted structured grids are always used to represent solid boundaries in MHD numerical simulations while imposition of boundary conditions is greatly simplified[25, 38, 24]. By this way, Ni et al.[25] has developed a consistent and conservative scheme on unstructured meshes, in which they numerically simulated a circular pipe flow under magnetic field. Recently, Vantieghem[38] has also studied the circular pipe flow using unstructured grids and examine how the wall conductivity affects the velocity distribution. In addition, Mutschke et al.[24] have studied the flow past a circular cylinder under magnetic field on exponential polar grids . In the field of hydrodynamics simulation, Cartesian grid methods have been known as an emerging technique attracting more interest for handling complex boundaries due to its convenient grid generation. There are several 2

methods fit into this category[19], such as immersed boundary method[27], ghost-cell method[37] and cut-cell method[4]. In these solution techniques, the governing equations are discretized on a structured Cartesian grid which doesn’t conform to the immersed boundaries. However, the obvious disadvantage of Cartesian grids is that they don’t represent curved and sloping surface well, and special treatment near solid boundaries are needed. During simulating MHD flows, another puzzle of fluid-wall coupling problem may arise based on Cartesian grids. When the solid wall is conductive, current density will flow into the solid wall and the computational domain for electric field is not necessarily limited to the fluid region but may rather extends to the solid wall region. This may not be a problem for unstructured grids or body-fitted grids because the fluid-wall interface conforms to the mesh boundaries, two groups of grids can be applied to solve the hydrodynamics field and the electric field respectively. However, in Cartesian grids, one mesh may contain both fluid and solid, which result in discontinuous conductivity coefficients in one mesh. Furthermore, solid meshes are abandoned when solving the hydrodynamics equations in cut-cell method, but all meshes are needed in obtaining electric field variables. Recently, Grigoriadis et al.[10] have introduce the immersed boundary method to deal with complex boundaries in MHD numerical simulation, in which methodology a discrete forcing field is used to impose the boundary conditions along the immersed surfaces and similarly, a direct discrete forcing field is used on electric current density field to satisfy the charge conservation law. However, the approach is restricted to flows with insulated wall, while no definitive boundary condition can be provided for current density along fluid-wall interface when the solid wall is conductive. The approach adopted here employs a cut cell method to treat complex boundaries, and on uniform grids without solid boundaries, the method reduces to the consistent and conservative scheme developed by Ni et al.[26]. Summarily, the main characteristics of our methodology are as below

3

• Cut cell method is applied on complex boundaries to satisfy the charge conservation law. In order to make this method work more efficiently, an unstructured adaptive octree grids are also adopted. • When conductive wall is considered, we use a VOF method to handle the liquid-solid coupling problem that discontinuous conductivity may exist in one mesh. • Conservative interpolation technique is introduced to keep the current density at the cell center conservative. • A consistent and conservative scheme is applied on fine/coarse mesh boundaries when calculating the electric potential flux and current density flux. The present paper is organized as follows. In section 2 the complete MHD equations are presented. The numerical strategy used is described in section 3, placing a key emphasis on the numerical treatment of current density calculation on fine/coarse mesh faces and how we make the cut-cell method applicable here. The proposed model is validated in section 4 by comparing with analytic and numerical solutions. 2. Governing Equations The set of equations governing the viscous incompressible fluid motion under the influence of an external magnetic field are the mass and momentum conservation equations, ∇·u =0  ρ

∂u + u · ∇u ∂t

(1)

 = −∇p + ∇ · Tv + F m

4

(2)

where ρ is the fluid density, u and p are the velocity vector and kinetic pressure respectively. Tv is the viscous stress tensor given by, Tv = 2μD

(3)

where μ is viscosity and D = 12 (∇u + ∇u T ) is the deformation tensor. F m stands for the volume Lorentz force which can be represented as Fm = J × B

(4)

where J is the induced current density and B is the magnetic field. Generally, two approaches may be used to evaluate the current density[23]. One is through solving a magnetic induction equation and the other is through solving an electric potential equation. In the first approach, the magnetic induction equation provides the coupling between the flow field and the magnetic field and is derived from Ohm’s law and Maxwell’s equation. The second approach is to solve the electric potential equation and calculate the current density using Ohm’s law. Recently, Smolentsev et al. [34] designed a numerical technique for MHD flows based on an induced current density, for which the detailed formula has been obtained. This algorithm overcome the difficulty of boundary condition for the magnetic induction formula, and it is well designed to have good accuracy for MHD flows at high Hartmann number. In this paper, the electrical potential formula is adopted which is proved to be accurate when the magnetic Reynolds number is small. Ohm’s law that defines the current density is given by: J = σE

(5)

where σ is the electric conductivity of the media, E is the produced electric field. For fluid velocity field u in a magnetic field B, Ohm’s law takes the

5

form J = σ(E + u × B)

(6)

in which the electric field can be expressed as E = −∇ϕ +

∂A ∂t

(7)

The first term is the scalar potential, while the second one is the vector potential representing the contribution from the induced magnetic field. Assuming the characteristic length and velocity are respectively L and u0 , if the magnetic Reynolds number Rm = μ∗ σLu (μ∗ is the magnetic permeability) is a small number, the induced magnetic field is much smaller than imposed magnetic intensity. Then Eq.(6) can be written as J = σ(−∇ϕ + u × B)

(8)

with the homogeneous Neumann condition of ϕ and the impenetrable condition of u along the insulated wall, given by ∂ϕ = 0; u · n = 0 ∂n

(9)

The main dimensionless groups governing the flow status are the Reynolds number Re = σLB02 ρu0

ρu0 L μ

and the interaction parameter or Stuart number N =

. In magnetohydrodynamics the interaction parameter gives a measure

for the ratio of the electromagnetic to the inertial force. Another derived dimensionless parameter Hartmann number, expressing the balance of forces established between the electromagnetic force and the viscous force, is given  √ by Ha = ReN = LB0 σμ .

6

3. Numerical Scheme In this paper, an approximate projection method making use of a staggeredin-time discretization is employed to discretize the hydrodynamic part, which is second order accurate. Combining the complete Navier-Stokes equations(2) with the electric field equation(8), a discretised form from time level k to k+1 is presented as below: 1

J k = σ k+ 2 (−∇ϕk + u k × B k )

k+ 12

ρ



1 u∗ − uk + (u · ∇u)k+ 2 Δt



(10)

1

= −∇pk + ∇ · (μk+ 2 (Dk + D∗ ))

(11)

+(J × B)

k

u k+1 = u ∗ −

Δt ρ

k+ 12

∇(pk+1 − pk )

∇ · u k+1 = 0

(12)

(13)

u ∗ , the intermediate velocity, obtained from the predictor step of Eq. (11) does not meet the condition of solenoidal velocity, which is then projected on a divergence-free vector space by Eq.(12). The Poisson equation of pressure difference is given below by combining Eq.(12) with Eq(13)  ∇·

Δt 1

ρk+ 2

 ∇(p

k+1

−p ) k

= ∇ · u∗

(14)

In Eq.(11), the prediction step of calculation u ∗ , the advection term 1

(u · ∇u)k+ 2 is evaluated using a conservative formulation developed from the Bell-Colella-Glaz second order unsplit upwind scheme[3]. On the other hand, an unstructured quad/octree based grid system is adopted for spatial 7

discretisation[8, 28], which is proven to be simple to implement a fully flexible adaptive refinement strategy. The mesh can be adapted at every time-step according to a mesh coarsening or refining criterion. The primary variables of the MHD equations, such as velocity u, pressure p and electric potential ϕ are located in the center of each finite volume. In the above equations, the electric field control equation Eq.(10) seems to be simple as it appears, however, hidden behind is an extra Poisson equation for the electric potential at each computation cycle. In order to speed up the convergence speed, a multi-grid method is adopted that it is high efficiency when using quad/octrees grid arrangement. According to the electric charge conservation law, we must enhance ∇ · Jk = 0

(15)

Combined with Eq.(10), an electric potential Poisson equation is derived at each time step[6] as 1

1

∇ · (σ k+ 2 ∇ϕk ) = ∇ · (σ k+ 2 (u × B)k )

(16)

As the collocated mesh arrangement presented earlier, each discretized volume in domain is called a cell and all the variables are located in the cell center while the current density fluxes are located on cell face. Then a consistent and conservative scheme should be designed to ensure a good solution as it was a true challenge[35, 15] in various mesh arrangement. Recently, Ni et al.[26] have examined the variable location on the consistent and the conservative nature of the collocated algorithm and proposed a numerical scheme for a collocated grid system, which was proven to be consistent and conservative while the simulation results is promising. In the light of the considerations above, we develop the consistent and conservative scheme adapted for the AMR technique that keeps the consistent and conservative properties, which is described in detail in the following 8

subsection By all accounts above, the numerical procedures for a timestep are as follows: • Cell centered electric potential at the time step k are calculated from Eq.(16) by solving a Poisson equation through multilevel approach. • Induced current density is calculated from Eq.(10) after the obtained electric potential is put in. • The Lorentz force is computed at the center of finite volume, then added to the predictor step of Eq.(11), by which the intermediate value u ∗ is obtained using a variant of the multilevel solver. • An approximate projection operator is applied on u ∗ that the pressure difference between time k + 1 and time k is computed by solving the Poisson equation(14) with the multilevel solver • The cell-centered velocity u k+1 is updated using the interpolated value of pressure difference pk+1 − pk by Eq.(12). 3.1. Cut-cell approach for MHD flows with insulated complex boundaries In this paper, the consistent and conservative scheme[26] is extended to an unstructured Cartesian adaptive system. This approach is implemented when applying a cut-cell method for dealing with complex geometries. For the sake of brevity, the flow state without solid boundaries is viewed as a particular case of flows with complex boundaries, and only the solution procedure of electric field is described as well. The solid boundary is first represented by a serial of piecewise linear segments while the cell cut by the boundary is reconstructed through the VOF(volume of fluid) method. When electric field is considered, calculations of the gradient of electric potential and the current density flux across

9

cell faces are crucial issues that should be accurately conducted. Simultaneously, the Hartmann layer and side layer in the vicinity of solid boundaries for MHD flows are relative thin with thickness proportional to 1/Ha and  1/Ha, respectively. Thus precise treatment is necessary near the boundaries while the Cartesian approach is just difficult in applying boundary conditions. 3.1.1. Consistent interpolation operator and gradient operator at cell face to ensure calculated current flux conservative In fact, if we apply a divergence theorem on Eq.(16) in a discrete cell, the discretisation equations yields 

 σ∇ϕ · ndS =

σ(u × B) · ndS

(17)

Then the descretisation form of the equation(17) can be written as, nf 

 σf

f =1

∂ϕ ∂n

 Sf = f

nf 

σf (u × B)f · nf Sf

(18)

f =1

where f denotes the cell face (interface as well if the cell is cut by boundaries), and n is the unit normal of cell face f . nf , denoting number of faces, is 4 in 2-dimensional case and 6 in 3-dimensional case. Sf is the area of face f . Assuming the length of cell edge to be h while the fluid face fraction is sf . Then Sf = sf h2 and a simplified equation is gained as below nf  f =1

 σf sf

∂ϕ ∂n

 = f

nf 

σf sf (u × B)f · nf

(19)

f =1

Usually, in order to keep the second-order accuracy on partial cell face when calculating electric potential gradient, a simple interpolating scheme from the surrounding cell-centers can be applied[40, 14, 7]. However, it’s more difficult within the adaptive framework and the consistent property is 10

not easily retained. Therefore, the full-face-center gradient and interpolation scheme is approximately used here, which results in first-order near solid boundaries. In uniform collocated grid system, it’s usually not difficult to design a consistent form for current density flux and potential gradient calculation on cell faces. But due to the AMR technique, mesh size may vary between neighbours that inconsistency happens when improper interpolation scheme adopted. In practice, three cases must be considered for the construction of the interpolation operator as shown in Fig.1. Take velocity interpolation for example as in the figure: Situation 1 as in Fig.1(a), if the neighbour cell Nd in direction d is at the same level and is a leaf cell, the interpolated value from cell-center is simply uf = (u + ud )/2. Situation 2 as in Fig.1(b), if the neighbour cell Nd with cell-center velocity ud in direction d is at the lower level, a three point interpolating scheme is used to preserve the second-order accuracy on cell face, given by 1 1 uf = u + ud − u∧d 5 5

(20)

where ud is interpolated from ud and its neighbouring cells in directions perpendicular to d. Situation 3 as in Fig.1(c), if the neighbour cell Nd is at the same level but not a leaf cell, an ordinary interpolate scheme of uf = (u + ud )/2 is not yet consistent because the flux of u across cell face calculated from cell C is uf Sf compared to (uf 1 Sf 1 + uf 2 Sf 2 )/2 which is calculated from cell Nd1 and Nd2 , and obviously they are not equal. Therefore we constructed uf as the average value of uf 1 and uf 2 in order to enhance the conservative and consistent property of flux across the cell face. The above approach gives a consistent discretized formula for interpola-

11

tion operators, however, it’s not sufficient for calculating the RHS of Eq.(19) because it’s a cross product of velocity and magnetic field. Consequently, in order to keep (Jnu )f = σf (u × B)f · nf consistent on coarse/fine cell boundaries, we should treat the flux (Jnu )f as a whole. In the same way, with the types of Situation 3, The flux (Jnu )f on coarse face should be calculated as (Jnu )f = ((Jnu )f 1 Sf 1 + (Jnu )f 2 Sf 2 )/(Sf 1 + Sf 2 ). We define an operator f lux(C, f, Jnu ) as the flux of Jnu across the face f of Cell C, and also define interpolate1(f, u) as the interpolation operators of u in Situation 1 and Situation 2 respectively. Then the consistent algorithm can be summarized as, switch(face-type) case(fine-fine) f lux(C, f, Jnu ) = σf (interpolate1(f, u) × interpolate1(f, B)) · nf break case(fine-coarse) f lux(C, f, Jnu ) = σf (interpolate2(f, u) × interpolate2(f, B)) · nf break case(coarse-fine) f lux(C, f, Jnu ) = f lux(Nd1 , fd1 , Jnu ) + f lux(Nd2 , fd2 , Jnu ) break Lastly, to discrete the LHS of Eq.(19), the scheme of electric potential gradient on cell face is constructed, in which the flux calculation of (Jnϕ )f = −σf (∇ϕ)f ·nf on non-uniform mesh boundaries also requires special treatment and the numerical method is similar to the interpolation procedure. Situation 1 as in Fig.1(a), if the neighbour cell Nd in direction d is at the same level and is a leaf cell, the flux of gradient can be given as:  (Jnϕ )f = −σf

∂ϕ ∂n

 = −σf f

12

ϕd − ϕ h

(21)

Situation 2 as in Fig.1(b), if the neighbour cell Nd in direction d is at the lower level with ϕd located at the cell center, the flux can be given as: (Jnϕ )f = −σf

−5ϕ − 3ϕ∧d + 8ϕd 15h

(22)

where ϕd is interpolated from ϕd and its neighbouring cells in directions perpendicular to d. Situation 3 as in Fig.1(c), if the neighbour cell Nd is at the same level but not a leaf cell, the flux (Jnϕ )f across cell face can be calculated from an interpolation as: (Jnϕ )f =

(Jnϕ )f 1 Sf 1 + (Jnϕ )f 2 Sf 2 Sf 1 + Sf 2

(23)

with (Jnϕ )f 1 and (Jnϕ )f 2 calculated from cell Nd1 and Nd2 based on equation (22), which can enhance the conservative and consistent property of flux across the cell face. When the cell is cut by insulated solid boundary, the flux of electric potential gradient (Jnϕ )f and the flux (Jnu )f are both zero under Eq.(9). After the discretezation procedure above is put in Eq.(19), an efficient multilevel method is used to exactly solve the Poisson equation and the electric potential is then obtained in cell-center. After that, according to Eq.(10), current density is calculated on cell face that we should notice the gradient operator used here need to be consistent with the discretization of the gradient operator in the electric potential Poisson equation(19). Moreover, to better enhance the consistency in calculating the current density on cell face, an average scheme is adopted on coarse/fine cell boundaries similar to that used in Situation 3, given by (Jn )f =

(Jn )f 1 Sf 1 + (Jn )f 2 Sf 2 Sf 1 + Sf 2

13

(24)

which is able to tighten the consistent property of current density calculation. The gradient operator, together with the interpolation procedure described above, is second-order accurate in space away from solid boundaries and uses a consistent flux estimation at coarse/fine cell boundaries. 3.1.2. Calculation of the Lorentz force at the cell center We calculate the current density fluxes on the faces of the control finite volume by the consistent and conservative discretized operators developed in section 3.1.1. Since the current density fluxes jx , jy , jz are located on the cell faces normal to the x−, y−, z− directions respectively, interpolation from the cell faces to cell center is needed to calculate the Lorentz force. A simple interpolation operator of J based on area average can be given as, n (ji )c =

f =1 (ji |si |)f nf f =1 |si |f

where i denotes the x−, y−, z− direction components, and |si |f is the dimensional normal area vector of partial face in cut cell. This scheme is adopted in[28] as it can conduct the hydrodynamics flow with good accuracy. However, as depicted in Ni at al.’s paper[25], this interpolation introduces numerical error for the calculation of the Lorentz force on an irregular mesh. Therefore, when cut cell is considered, we adopted the conservative formula-

14

tion as developed in [25], given as 1 Jc = Ωc



1 J dΩ = ∇ · (J (r − r c ))dΩ Ω c Ωc Ωc  1 Jn (r − r c )ds = Ω c sf nf 1  (Jn )f (r f − r c )sf = Ωc f =1 nf 1  = (Jn )f sf df nf Ωc f =1

(25)

here Ωc denotes the effective cell volume and sf is the partial face fraction, and we use a expression of r f − r c = df n f where df denotes the distance from center of partial cell to center of partial face. When the cell is a regular cubic, then Eq.(25) can be reduced as, Jc =

nf  f =1

1 sf df n f = (Jn )f n f Ωc 2 f =1 nf

(Jn )f

(26)

and this is the simple average value of interpolating current density flux from cell face to cell center which also suggests the simple interpolation for current density flux in a complete square cell is conservative and consistent. With the conservative interpolation operator of Eq.(25), we are able to calculate Lorentz force through F m = J c × B c in cells cut by boundaries with good accuracy. 3.2. cut-cell approach for MHD flow with conductive complex boundaries Different from the hydrodynamics flow in which situations the meshes in solid walls are always destroyed and abandoned since no variable needed there. However, in MHD numerical simulation, the existence of conductive wall seriously affected the current density distribution[12], results in dis-

15

tinctive flow characters from the cases with insulated wall. Therefore, the hydrodynamics simulating region is different from that of electric field which also contains the solid wall. Grigoriadis et al.[10] has introduced the immersed boundary method to simulate MHD flows with insulated wall by applying a forcing term during calculating the current density to make J = 0 on fluid-wall interface. But this method can’t be easily extended to MHD flows with conductive wall because no definite boundary condition exists along the rigid wall, and the electric potential and current density are continuous across the interface. In this paper, we develop a conservative approach using the VOF method to deal with cells lying astride the fluid region and solid region, so as to make cut-cell method applicable in MHD flows with conductive wall. And this technique is based on the consistent and conservative scheme covered earlier in section 3.1 for handling insulated wall case. VOF method [11, 29, 31] is originally designed to track the fluid interface, for which the volume fraction is used to represent the fluid in a cell and the interface needs to be reconstructed based on the fraction. Here, PLIC (piecewise linear interface construction)[29] technique is employed to treat the complex wall with good accuracy. 3.2.1. Solving the electric potential Poisson equation with jump conductivity The difficulty within the cut-cell method, evidently, is solving the variable coefficient Poisson equation with jump conductivity in cut-cell. Compared to the pressure Poisson equation in which only the fluid part is considered, the fluid part and solid part should be conjugated in solving electric potential Poisson equation on cut cell. Therefore, the destroyed cells in hydrodynamics calculation should be reconstructed first and all the corresponding variables, such as p, u, are copied to the newly generated grid system. For instance, the fluid and the solid in the cell (i, j) is distributed as Fig.2, the disretized integral form of electric potential Poisson equation of Eq.(17)

16

should be formulated as:

     nf nf   ∂ϕ ∂ϕ Sf 1 + σf 2 Sf 2 = σf 1 (u × B)f · n f Sf 1 (27) σf 1 ∂n f 1 ∂n f 2 f =1 f =1 where the superscripts 1 and 2 indicate fluid and solid respectively, while S1 = ∂ϕ |CF |, S2 = |F B| on face |CB| in Fig.2. When calculating ∂ϕ and ∂n f 1 ∂n f 2 on cell face, Johansen and Colella [14] suggested that a flux on a partial cell is evaluated based on linear interpolation between adjacent fluxes, so ∂ϕ is ∂n f 1 obtained using a quadratic polynomial based on the surrounding nodal values in the fluid part, and alternatively, ∂ϕ is evaluated from the surrounding ∂n f 2 solid nodes. Another second-order accurate method for variable coefficient Poisson equation is presented by Li et al.[13], in which the solid part and fluid part are calculated separately. What’s more, Liu et al.[17] have introduced the Ghost Fluid Method(GFM) to capture the boundary conditions at a contact discontinuity. However, different from the method of Liu et al., the jump condition is not clear along the interface in our problem. Referring to immersed boundary method, LeVeque et al.[16] have developed a secondorder numerical method to preserve the jump condition at the interface too. As mentioned before, it’s complex to develop these interpolation schemes with the adaptive framework and the existing cell-center pressure solvers could not be directly used for the Poisson solution of ϕ if so. Therefore, we choose the uniform gradient of electric potential on mixed cell face and the cut-cell is taken as an integral as well considering that ϕ is continuous across the liquid-solid interface, which results in, nf  f =1

 (σf 1 Sf 1 + σf 2 Sf 2 )

∂ϕ ∂n

 = f

nf 

σf 1 (u × B)f · n f Sf 1

(28)

f =1

If the solid part is assumed to be insulated, which indicates σ2 = 0 almost,

17

then Eq.(28) is simplified as, nf 

 σf 1 Sf 1

f =1

∂ϕ ∂n

 = f

nf 

σf 1 (u × B)f · n f Sf 1

(29)

f =1

which is consistent with the discretized form of Poisson equation of ϕ with insulated wall, as depicted in Eq.(19). Substitution the electric conductivities of fluid and solid, together with the calculated parameters and the interpolation operator introduced earlier, a new Poisson equation of ϕ can be solved using an effective multilevel approach. 3.2.2. Calculation of the Lorentz force at the cell center When cell is cut by boundaries, we calculate Lorentz force merely in fluid part though current density exists in solid domain. Then a consistent and conservative interpolation scheme similar with Eq.(25) is adopted for calculating the current density from cell face to cell center. However, note that the boundary condition of electric potential in Eq.(9) is not effective any more as

∂ϕ ∂n

= 0 on interface, which is |EF | in Fig 2.

It should be noticed that the gradient flux of electric potential through the same piece of the interface |EF | in the fluid part can be different from that in solid part due to the jump condition of electric conductivity across the interface. A simple relationship on interface is given by,  σf

∂ϕ ∂n



 = σ1 f

∂ϕ ∂n



 = σ2 1

∂ϕ ∂n

 (30) 2

For the electric potential is continuous across the interface, here we calculate the

∂ϕ ∂n

on interface using a bilinear interpolation scheme from the

surrounding nodes, whatever lie in solid or lies in fluid. We use the corresponding value of ϕ in neighbour cells to determine the coefficients a, b, c, d

18

of the bilinear function in 2D case ϕ = ax + by + cxy + d We then use these coefficients to approximate

(31) ∂ϕ ∂x

and

∂ϕ ∂y

on |EF | at the

center P of the interface. For example, we have ϕx = a + cyp ϕx = b + cxp

(32)

When evaluating the conductivity on interface, the weighted harmonic mean interpolation scheme is proven to be accurate when treating such coupling problems,written as 1 c 1−c = + σf σ1 σ2

(33)

where c denotes the area ratio of the fluid in mixed-cell and 1 − c is the solid part accordingly. After that, the flux of current density on interface is calculated as  (Jn )f = (Jnϕ )f = −σf

∂ϕ ∂n

 (34) f

where the (Jnu )f part of (Jn )f turns into zero on the interface. Lastly, synthesizing the above-mentioned schemes for conductive wall situation, this paper puts forward the approach for calculating the cell-center current density according to Eq.(25), which scheme is consistent and conservative in irregular finite volume. Moreover, while the wall is insulated which indicates

∂ϕ ∂n

= 0 on interface, the method is consistent with the scheme we

presented in Section 3.1 for flows with insulated complex boundaries

19

4. Numerical Results 4.1. convergence for a simple periodic problem In the former part we have presented a consistent scheme for calculating (Jn )f = (Jnϕ )f + (Jnu )f = −σf ∂ϕ + σf (u × B)f · n f on cell face f when ∂n f solving the electric potential Poisson equation, which is of great importance in keeping the charge conservation law. So the first step is to test the spatial accuracy properties of the method in solving the electric potential, especially when the test is based on non-uniform grids. And also we will compare the results with that when the consistent scheme on coarse/fine cell boundaries in section 3.1.1 is not adopted. In hydrodynamics, there exist various and credible simulations for accuracy test. For MHD simulations, an authentic unsteady problem with exact analytic solution has not been posed and extensively used in previous literatures. Here, we consider a simple periodic problem in a square unit domain which has been studied in hydrodynamics before[18, 1, 28]. The initial conditions are taken as u(x, y) = 1 − 2 cos(2πx) sin(2πy), v(x, y) = 1 + 2 sin(2πx) cos(2πy).

(35)

where the boundary condition for both x− and y− directions are periodic, then the exact solution of the Euler equations is presented as u(x, y, t) = 1 − 2 cos(2π(x − t)) sin(2π(y − t)), v(x, y, t) = 1 + 2 sin(2π(x − t)) cos(2π(y − t)), p(x, y, t) = − cos(4π(x − t)) − cos(4π(y − t)).

(36)

In order to make it useful for MHD simulations, we impose a uniform field Bz perpendicular with the flow plane, and the Lorentz force is not considered

20

in solving the Euler equations. Then the electric potential can be yielded as ∇2 ϕ = ∇ · (u × B) ∂(vBz ) ∂(uBz ) = − ∂x ∂y 1 ⇒ ϕ = − cos(2π(x − t)) cos(2π(y − t)) π

(37)

For our purpose is to determine whether the electric potential calculation procedure is consistent on coarse/fine cell boundaries, we now call the proposed consistent scheme of interpolation algorithm and gradient algorithm in section 3.1.1 as Scheme I. And controversially, if the flux calculation on coarse/fine cell meshes are considered separately in each cell, which means Situation 3 in section 3.1.1 is not considered specially, we name it Scheme II. This test is first conducted on uniform grids with L = 5, 6 and 7 levels of refinement and then one or two additional levels are added in a square domain defined by the points(-0.25,-0.25) and (0.25,0.25). All simulations are performed at the same CFL number of 0.75 and the running time for each is 0.5. When Scheme I is applied in calculating (Jn )f on cell faces, Fig.3 shows the spatial evolution of the L2 and L∞ norms of the error in electric potential ϕ. As the first figure presents, the method preserves a second order accuracy on uniform mesh. Meanwhile, even if the meshes in center of domain are locally refined one or two extra levels, as shown in the second and third figures, the accuracy is still close to second-order. These tests provide the evidence of the Scheme I proposed in section 3.1.1 is consistent and satisfy the charge conservation laws. Subsequently, we use Scheme II to calculate (Jn )f on cell faces, and the results are shown in Fig4. As it clearly depicts, if the mesh is uniform, secondorder convergence for ϕ is still obtained as presented in the first picture. However, in the next two pictures when additional one or two levels are

21

added to the center of domain, the accuracy will decrease rapidly close to first-order. Otherwise, from the comparison between three figures, the error norms are inversely increased when the mesh is refined with extra levels in domain center. This phenomenon is attributed to the inconsistent treatment of (Jn )f , which leads the current density flux not consistent on coarse/fine cell faces. 4.2. MHD flow in a diverging channel To demonstrate the convergence performance of the method in simulating MHD flows with solid boundaries, we select a test case appeared in[2] and the geometrical definition is shown in Fig.5. The diverting channel is constructed in a 4×1 domain and two opposite-facing plates of ybot and ytop cut the domain to restrict the fluid flow between them as detailed depicted as:

ybot

⎧ ⎪ y1 if 0 ≤ x ≤ 1 ⎪ ⎪ ⎨  π  (x − 1) if 1 < x < 3 = y2 + 0.5(y1 − y2 ) 1 + cos ⎪ 2 ⎪ ⎪ ⎩ y2 if 3 ≤ x ≤ 4

(38)

ytop = 1 − ybot A uniform unit inlet velocity is imposed at the inlet while the convective boundary conditions are applied at the outlet. Then we impose an external magnetic field Bz perpendicular to the flow plane, and the Neumann conditions for the electric potential are set at the solid boundary. Since the magnetic field is along z−direction, electric potential is generated along the solid surface. Due to the Ohm’s law of Eq.(5), when the flow recovers after the diverging region, the current density is zero and the velocity field should remain unaffected by the external magnetic field. Correspondingly, analytic solution for the exist velocity and electric potential can be obtained. The simulations are using a CF L = 0.9 and all ran to t = 3.05 which 22

is proven to be sufficient to reach a state condition for the MHD flow. The non-dimensional parameters are respectively Re = 10 and Ha = 10 with the half length of inlet as the characteristic length. Firstly, three simulations are performed on uniform mesh with L = 5, 6 and 7 levels of refinement. Then we examine the variation of the error norms when the cells cut by boundaries are refined with additional level. The solution of x−component velocity along the outlet is given as u = −4.6875(y 2 − 0.16)

(39)

Combined with the Neumann boundary condition, if we solve the electric potential Poisson equation base on the velocity distribution along the outlet, the solution of ϕ is  ϕ = −46.875

 y3 − 0.16y + c 3

(40)

where c is an arbitrary constant and is taken as the average value of the computed electric potential along the outlet. Based on the uniform mesh of L = 6, the numerical solution for velocity and electric potential are presented in Fig.6. A parabolic velocity profile is predicted at the outlet and the distribution of electric potential is also in good agreement with the analytic solution, clearly demonstrate the validity of the cut-cell method applied in wall-bounded MHD flows. The detailed analysis of error norms are presented in Fig.7 for boundaryrefined mesh. The evolutions of L1 , L2 norms of velocity both indicate the preservation of a second-order of accuracy in space while the maximum error L∞ norm of velocity is close to first-order accurate. It’s proven that the numerical solution of velocity with solid boundaries are globally second-order accurate but only first-order convergence is obtained near the solid boundary. However, the evolutions of L1 , L2 and L∞ norms of electrical potential

23

indicate the preservation of a second-order of accuracy in space. 4.3. MHD flow in rectangular channels There are some exact solution for fully developed incompressible laminar flows in ducts with transverse magnetic fields. Shercliff’s solution for rectangular ducts with non-conducting walls and the magnetic field perpendicular to one side[32] and Hunt’s solution for rectangular ducts with two non-conductive side walls and two conductive Hartmann walls[12] have been used to validate the consistent and conservative scheme developed here under AMR framework. Shercliff’s case will be simulated without solid boundaries to test the fundamental elements of the proposed scheme. Hunt’s case will be simulated using the methodology described in Section 3.0 and section 3.1, where the insulated wall and conducting wall will be respectively represented with cut-cell method. 4.3.1. Shercliff ’s Case We simulate Shercliff’s case of a fully developed flow in a square channel with all the walls electrically insulated. We consider the flow travelling in x−direction with the cross section of 2d × 2d in y− and z− directions. A unit magnitude of velocity is imposed at the inlet and periodic boundary conditions are applied aligned along the x−direction, meanwhile a uniform magnetic field is applied along the y−direction. The entrance velocity is selected as the characteristic velocity while the Reynolds number is Re = 25 and the Hartmann number is Ha = 100. The non-uniform meshes are selected to put forward the simulation with fine grids in the side layer and finer grids in the Hartmann layers. In detail, the mesh in the 1.5 times thickness of Hartmann layers is resolved with cells having minimum size 1/(3Ha), and the mesh in the 2.5 times thickness of side layers is resolved with cells having minimum size 1/(5Ha1/2 ). The mesh is shown in Fig.8. With a close look at the grids near a corner as shown in

24

Fig.8(b), we can clearly see that the grids are non-uniform with coarse grids in the core, fine grids in the side layer and finer grids in the Hartmann layer. The consistent and conservation scheme is used to conduct the calculation of the current fluxes on cell faces, which meet the charge conservation law on coarse/fine cell boundaries. The Lorentz force are calculated based on conservative formulation of Eq.(25). A comparison between the exact solution from the Shercliff’s analysis and the numerical solution are presented in Fig.9. The results on adaptive mesh match the analytic results well. The comparison demonstrates a very good computation accuracy. 4.3.2. Hunt’s Case We now consider the fully developed flow in a square channel with Re = 25, Ha = 100 of Hunt’ case. The walls perpendicular to the applied magnetic field are assumed to be conductive and those parallel to the applied magnetic field are electrically insulated. In the numerical result, the presented consistent scheme in section 3.0 are used to deal with the cells cut by insulated wall while the consistent scheme in section 3.1 are used to simulate the cells cut by conducting wall. The Lorentz force is calculated based on the conservative formulation of Eq.(25), especially for the irregular cut-cells. We still adopt an non-uniform grid framework to carry out the simulations as same as in the Shercliff’s case. However, since the velocity distribution is enormously affected by a non-dimensional parameter of conductance ratio, which is given by cw =

σw tw σf L

(41)

where the superscript w and f denotes the solid wall and fluid respectively, this paper simulated the Hunt’s case with symmetrical wall and unsymmetrical wall. Since cw describes the ratio of the electrical conductance of the wall and the fluid material, we first calculate the symmetrical case with cw = 0.4 for both walls perpendicular to the magnetic field and this analytic solution 25

is obtained by Hunt[12]. Recently, Tao & Ni[36] have studied the unsymmetrical cases in Hunt’s case, thus we also consider another two cases: (i) the conductivities are different in the top wall and the bottom wall; ( ii) the thickness of the top wall is different from that of the bottom wall. For symmetrical wall case, the column picture for the current density and the velocity distributions are shown in Fig.10. The pictures illustrate that the resolution of the Hartmann and side layers was very accurately captured and the numerical results don’t suffer from any instabilities. Thus, the treatment of current density on cut-cell face and coarse/fine face boundaries can well keep the conservative properties of charge. Moreover, Fig.11 presents the comparison between the exact solution from Hunt’s analysis and the numerical solution, and the agreement was perfect. Then a calculation model while the top wall’s conductivity is twice the magnitude of the bottom wall is put forward. The velocity profile are presented in Fig.12 which shows the results are in good agreement with Tao & Ni’s analytic solution[36]. The comparison demonstrates the approach settled for MHD flow with conducting wall in section 3.1 has a very good accuracy. Moreover, as the left figure shows, the velocity profile parallel to the magnetic field is not symmetrical any more, the wall with higher conductivity extends the influence of magnetic field on fluid motion. Its existence results in more significant jet effect of velocity, consequently, the velocity decrease more near the top wall Lastly, we perform the simulation while the bottom wall’s thickness is twice the height of the top wall. The comparison between the numerical solution and the analytic solution by Tao & Ni is presented in Fig.13, which shows the two results are in good agreement. Otherwise, due to the different thickness of walls perpendicular to magnetic field, the velocity profile is unsymmetrical as the thicker wall produces large velocity decreasing. The results above prove that the cut-cell approach with VOF method is sufficient simulating MHD flow with complex solid boundaries even when the

26

jump conductivity exists in a cut-cell. 4.4. MHD flow in circular pipes having finite conductivity Concerning the flow through a circular pipe under a uniform transverse magnetic field, several analytic solutions have been obtained with insulated as well as conductive walls[33, 9, 30, 5]. However, numerical difficulties are encountered in analytic methods when Ha number is large due to the infinite series expansions involving modified Bessel functions, and therefore an asymptotic method is developed for circular pipe flow at high Ha numbers. Other than that, under this approximate method, Chang[5] has studied the MHD flow in a circular pipe under the effect of wall’s conductivity and no M type velocity profile appears in his solution. Compared with Chang’s results, a later analytic work by Samad[30] gave the M type velocity profile at the small Ha number and big wall conductivity conditions. In this section, we will simulate the MHD flow in a circular pipe having finite conductivity and finite thickness under uniform transverse magnetic fields and then the results are compared with Samad’s analytic solutions. Let the inner and outer radius of the pipe be, respectively, R and 1.1R. The outer boundary of the pipe is also reconstructed using VOF method. As aforementioned, all the analytic solutions are difficult to evaluate when Ha number gets large that Samad limit his computation to Ha = 18 in his paper, but now we can extend the limitation to Ha = 30 due to the development of computers. So we will compare the numerical results of velocity normal to magnetic field with Samad’s analytic solution at Ha = 18 and Ha = 30 respectively as Fig.14 shows. It’s worthwhile to note that in the figures, c is the conductivity ratio between the wall and the fluid which is defined as c =

σw . σf

and two groups of c are considered in our numerical simulations

as c = 1 and c = 10 are both calculated at each Ha number. As Fig.14 shows, when c is greater, the M type velocity profiles are more obvious and all numerical results suit well with the analytic solutions.

27

Furthermore, since the numerical method is validated at small Ha numbers by comparing with Samad’s analytic solution, we then extend the numerical simulation to Ha = 100 and the velocity distribution along the centroid line vertical to magnetic field are compared with Chang’s asymptotic solution that the comparisons are showed in Fig.15. In the figure, cw is the so-called wall-conductance ratio defined as cw =

σw dw σf R

where dw and R are

respectively the wall thickness and the pipe inner radius. cw = 0.1, 1, 10 are individually selected in our numerical simulations and through the comparison with Chang’s asymptotic solution, there is a significant discrepancy in the results where small zones of overspeed appear near the wall in the numerical solutions, while the profile seems flat in asymptotic approximations. Moreover, the figure shows that the velocity jets appear more apparent when cw increases. This M type velocity distribution is rational and identical with our previous results as well as Samad’s analytic solutions at Ha = 18, 30. Finally, the numerical results at Ha = 100 are put together in Fig.16. 4.5. MHD flow past a circular cylinder The proof of the methodology’s validity is continued by studying the three-dimensional flow around a circular cylinder which has been studied by Mutschke[24] and Grigoriadis[10]. As it is widely recognized that in this flow model, the transition scenario from two-dimensional instability to a treedimensional instability is occurred at Re ≈ 194[39]. When the instability is transformed from 2-D case to 3-D case, then the electric potential is generated and contributes to the current density for an aligned or transverse magnetic field. Base on the above argument, the flow around a circular cylinder is simulated under Re = 200 in order to test the cut-cell method for calculating the electric field. Furthermore, the numerical results are finally compared with what had come before from[24, 10] for a stream-wise applied magnetic field at N = 0, 0.2, 1. We consider the flow of an electrically conducting, incompressible fluid around a circular cylinder and the computational domain is shown in Fig.17. 28

The duct walls and the cylinder are assumed to to be electrically insulating. The cylinder of diameter d is located at a distance of L1 = 8d from the inlet while the outlet at L2 = 24d behind it. In the vertical z−direction, the cylinder is set symmetrically and we extend Lz = 16d in order to minimize the effect of the outer boundary on the development of the wake. A constant value of u is imposed at the inlet of the domain while convective outflow conditions have been set for the outlet. Beyond that, the impervious boundary conditions are applied on the other walls. A homogeneous externally applied magnetic field B is imposed along the mainstream x− direction. To make use of the AMR framework, two different grids have been carried out with and without adaptive refinement technique named as Grid I and Grid II, respectively. The meshes with minimum size of x = 0.03125d around the cylinder are used in both grids which make 32 ×32 cells to resolve the area d × d. Furthermore, the adaptive refinement technique in Grid I is used according to the vortex field and velocity gradient, with maximum size of x = 0.25d and minimum size of x = 0.03125d. In order to compare with the previous data from[24, 10], two Hartmann numbers are considered with Ha = 6.32 and Ha = 14.14, which result in Nx = 0.2 and Nx = 1.0, respectively. During the whole simulations, the major parameters we study and compare are the lift and drag coefficients Cl and Cd , the recirculation pattern length Xr . Additionally, the Strouhal number St is used to describe the frequency of vortex generation when the flow becomes unsteady. All the above non-dimensional parameters are given as Cl =

Fl , 1 ρu2 d 2

Cd =

Fd , 1 ρu2 d 2

St =

fd u

(42)

where d and u are respectively the characteristic length and velocity, selected as the cylindrical diameter and the inflow velocity. The variation of Cd and Cl for both grids are presented in Fig.18. In the 29

Table 1: comparisons for MHD flow around a circular cylinder at Re = 200 Case(grid) Present(I) Present(II) Grigoriadis(coarse grid) Grigoriadis(fine grid) Mutschke Present(I) Present(II) Grigoriadis(coarse grid) Grigoriadis(fine grid) Mutschke Present(I) Present(II) Grigoriadis(coarse grid) Grigoriadis(fine grid) Mutschke

Nx 0 0 0 0 0 0.2 0.2 0.2 0.2 0.2 1.0 1.0 1.0 1.0 1.0

Cd 1.33 1.32 1.29 1.30 1.30 1.10 1.13 1.04 1.08 1.07 1.23 1.20 1.17 1.20 1.19

Cl 0.63 0.57 0.50 0.63 0.60 0.06 0.06 0.07 0.09 0.09 0.0 0.0 0.0 0.0 0.0

Xr −− −− −− −− −− −− −− −− −− −− 3.17 3.18 3.13 3.10 3.11

St 0.197 0.196 0.204 0.204 0.197 0.171 0.167 0.174 0.172 −− −− −− −− −− −−

first scenario as the figure shows, when pure hydrodynamics flow simulations named Nx = 0 are conducted on both grids, the flow developed to a 3-D unstable state as the calculation went on. The flow transition from temporary 2D state to the 3D state is marked by a sudden decrease in the lift and drag forces which has been studied by M¨ uck et al.[22] for flows around the rectangular obstacle. The St number in 2D state (St = 0.210) is slightly higher than that in the following 3D state (St = 0.197) which is also found by M¨ uck. And meanwhile, the slight disturbances of the lift and drag force aptitude represents the 2D to 3D instability transition. A distinguish phenomenon is found that in Grid I, the vibration amplitude of Cd is more intense than that in Grid II, which may be attributed to the dynamic adaptive mesh in domain. In the second scenario, under a small magnetic field with Nx = 0.2, the drag coefficient will drop sightly while a reduction in Cl is noticed to indicate that the flow is more stable compared to hydrodynamics situation. However, it still remains 3D instability marked by the slight disturbance of lift and drag force, this is compatible with the study of Mutschke[24] that transition to be a 2D stable flow was found to take place for interaction parameter values in the range 0.3 < Nx < 0.4. Besides that, we also notice it should 30

take more time to transit from 2D state to 3D. In the last scenario with continuous increasing of magnetic field to Nx = 1, Cd is observed to conversely rise, which is in agreement with previous study[24, 10]. Contrarily, Cl has monotonically continued to drop until zero which represents the flow performance grows more stable due to the Lorentz force. As the figure shows, it reaches the steady-state in an earlier time when tu/d = 30. Table.1 presents a subset of the computed cases and the results in previous study at Re = 200. All non-dimensional parameters on both grids are found in good agreement with previous studies. 5. Conclusions A consistent and conservative scheme for calculating the electric potential and current density is proposed on coarse/fine mesh boundaries under adaptive framework. Thin layers along the boundaries are exactly solved using finer grids compared to coarse mesh used in central region for lower computational cost. The cut-cell approach for hydrodynamics flows has been extended to simulate the MHD flows with complex boundaries. Combined with volume of fraction technique, the cut cell method is used to simulate MHD flows with complex conducting walls, which means the solver region of hydrodynamics is different from solver region of the electric field. When meshes are no longer uniform under adaptive mesh framework, it’s a crucial issue to ensure the current density flux conservative on coarse/fine cell boundaries. We especially consider the fluxes estimation of −σf ∂ϕ ∂n f

and σf (u × B)f · nf on coarse/fine cell boundaries when solving the electric potential Poisson equations and the approach is presented. Furthermore, the methodology dealing with non-uniform meshes in this paper is proven to be consistent and conservative, and simultaneously we give a comparative worse results when no special treatment is considered to evaluate the flux on coarse/fine boundaries. 31

Base on Cartesian grids, numerically simulating MHD flows with conductive wall is more difficult than that with insulated wall, because parameters is not continuous and jump condition may occur in one cut-cell. VOF method is proven to be efficient dealing with this. Several steady and unsteady test cases are carried out to validate the numerical methodology. All the results are in good agreement with the analytic solution and published literature. Based on adaptive mesh refinement framework, Hartmann layers and side layers are exactly resolved while the charge conservation law is satisfied near the cut boundaries. On the whole, as the numerical results showed, the consistent and conservative scheme on an unstructured Cartesian adaptive system combined with cut-cell method and VOF method can efficiently simulate MHD flows in a complex geometry with insulated or conductive walls. 6. Acknowledgment The authors acknowledge the support from the NSFC (Natural Science Foundation of China) under Grants 50936006 & 11125212 and ITER Programme from Ministry of Science and Technology of China under Grant 2013GB114001. References [1] A.S. Almgren, J.B. Bell, P. Colella, L.H. Howell, M.L. Welcome, A conservative adaptive projection method for the variable density incompressible Navier-Stokes equations, Journal of Computational Physics 142 (1998) 1-46. [2] A.S. Almgren, J.B. Bell, P. Colella, T. Marthaler, A Cartesian grid projection method for the incompressible Euler equations in complex geometries, SIAM J. Sci. Comp. 18(1997)1289–1309.

32

[3] J.B. Bell, P. Collela, H.M. Glaz, A second-order projection method for the incompressible Navier-Stokes equations, Journal of Computational Physics 85 (1989) 257-283. [4] D. Calhoun, R.I. LeVeque, Solving the advection-diffusion equation in irregular geometries, Journal of Computational Physics 156 (2000) 1-38. [5] C. Chang, S. Lundren, Duct flow in magnetohydrodynamics. ZAMP 100-114 (1961) [6] P.A. Davidson, An introduction to Magnetohydrodynamics, Combridge Texts in Applied Mathmatics, Combridge University Press, 2001 [7] E.A. Fadlun, R. Verzicco, P. Orlandi, J. Mohd-Yusof, Combined immersed-boundary finite-difference methods for three-dimensional complex flow simulations, Journal of Computational Physics 161 (2000) 35-60. [8] B. Fryxell, K. Olson, P. Ricker, F. X. Timmes, M. Zingale, D. Q. Lamb, P. MacNeice, R. Rosner, J. W. Truran, and H. Tufo, FLASH: an adaptive mesh hydrodynamics code for modeling astrophysical thermonuclear flashes, Astrophysical Journal Suppl. 131 (2000) 273-334. [9] R.R. Gold, Magnetohydrodynamic pipe flow. Part 1. Journal of Fluid Mechanics 13 (1962) 505-512. [10] D.G.E. Grigoriadis, S.C. Kassinos, E.V. Votyakov, Immersed boundary method for the MHD flows of lquid metals, Journal of Computational Physics 228 (2009) 903-920. [11] C.W. Hirt, B.D. Nichols, Volume of fluid (VOF) method for the dynamics of free boundaries, Journal of Computational Physics 39 (1981) 201-225.

33

[12] J.C.R. Hunt, Magnetohydrodynamics flow in rectangular ducts, Journal of Fluid Mechanics 21 (1965) 577-590. [13] H. Ji, F.S. Lien, E. Yee, An efficient second-order accurate cut-cell method for solving the variable coefficient Poisson equation with jump conditions on irregular domains, International Journal of Numericcal Method in Fluids 52 (2006) 723-748. [14] H. Johanson, P. Colella, A cartesian grid embedded boundary method for poisson’s equation on irregular domains, Journal of Computational Physics 147 (1998) 60-85. [15] L. Leboucher, Monotone scheme and boundary conditions for finite volume simulation of magnetohydrodynamics internal flows at high Hartmann number, Journal of Computational Physics 150 (1999) 181-198. [16] R.J. LeVeque, Z. Li, The immersed interface method for elliptic equations with discontinuous coefficients and singular sources, SIAM J. Numer. Anal. 31 (1994) 1019. [17] X.-D. Liu, R.P. Fedkiw, M. Kang, A boundary condition capturing method for Poisson equation on irregular domains, Journal of Computational Physics 160 (2000) 151-178. [18] M.L. Minion, A projection method for locally refined grids, Journal of Computational Physics 127 (1996) 158-178. [19] R. Mittal, G. Iaccarino, Immersed Boundary Methods, Annual Review of Fluid Mechanics 37 (2005) 239-261. [20] R.J. Moreau, Magnetohyedrodynamics. Kluwer Academic Publishers, 1990.

34

[21] N.B. Morley, S. Smolentsev, L. Barleon, I.R. Kirillov, M. Takahashi, Liquid Magnetohydrodynamics-recent progress and future direction for fusion, Fusion Engineering Design 51-52 (2000) 701-703. [22] B. M¨ uck, C. G¨ unther, U. M¨ uhler, L. B¨ uhler Three-dimensional MHD flows in rectangular ducts with internal obstacles, Journal of Fluid Mechanics 418 (2000) 265-295. [23] U. M¨ uller, L. B¨ uhler, Magneto-fluid-dynamics in Channels and Containers, Springer, 2001. [24] G. Mutschke, G. Gerbeth, V. Shatrov, A. Tomboulides, The threedimensional instabilities of the cylinder wake in an external magnetic field:a linear stability analysis, Physics of Fluids 13 (2001) 723-734. [25] M.-J. Ni, R. Munipalli, P. Huang, N.B. Morley, M.A. Abdou, A current conservative scheme for incompressible MHD flows at a low magnetic reynolds numeber. Part 2: on an arbitrary collocated mesh, Journal of Computational Physics 227 (2007) 205-228. [26] M.-J. Ni, R. Munipalli, N.B. Morley, P. Huang, M.A. Abdou, A current conservative scheme for incompressible MHD flows at a low magnetic reynolds number. Part 1: on an rectangular collocated grid system, Journal of Computational Physics 221 (2007) 174-204. [27] C.S. Peskin, Flow patterns around heart valves:a numerical method, Journal of Computational Physics 10 (1972) 252-271. [28] S. Popinet, Gerris: a tree-based adaptive solver for the incompressible Euler equations in complex geometries, Journal of Computational Physics 190 (2003) 572-600. [29] W.J. Rider, D.B. Kothe, Reconstructing volume tracking, Journal of Computational Physics 141 (1998) 112-152. 35

[30] S. Samad, The flow of conducting fluids through circular pipes having finite conductivity under uniform transverse magnetic fields. J. Appl. Mech. 34(1) (1967) 29-36. [31] R. Scardovelli, S. Zaleski, Direct numerical simulation of free-surface and interfacial flow, Annual Review of Fluid Mechanics 31 (1999) 567-603. [32] J.A. Shercliff, Steady motion of conducting fluids in pipes under transverse magnetic fields, Proceedings of Cambridge Philosophical Society 49 (1953) 126-144. [33] J.A. Shercliff, The flow of conducting fluids in circular pipe under transverse magnetic fields, Journal of Fluid Mechanics 1 (1956) 644-666. [34] S. Smolentsev, S. Cuevas, A. Beltr´an, Induced electric current-based formulation in computations of low magnetic Reynolds number magnetohydrodynamic flows, Journal of Computational Physics 229 (2010) 1558-1572. [35] A. Sterl, Numerical simulation of liquid-metal mhd flows in rectangular ducts, Journal of Fluid Mechanics 216(1990). 161-191 [36] Z. Tao, M.-J. Ni, Analytical solutions for MHD flow at a rectangular duct with unsymmetrical walls of arbitrary conductivity, SCIENCE CHINA Physics, Mechanics & Astronomy, submitted for review. [37] Y.-H. Tseng, J.H. Fersiger, A ghost-cell immersed boundary method for flow in complex geometry. Journal of Computational Physics 192 (2003) 593-623 [38] S. Vantieghem, X. Alberts-Chico, B. Knaepen, The velocity profile of laminar MHD flows in circular conducting pipes, Theror. Comput. Fluid. Dyn. 23 (2009) 525-533.

36

[39] H.K. Williamson, Vortex dynamics in the cylinder wake, Journal of Fluid Mechanics 328 (1996) 345-539. [40] T. Ye, R. Mittal, H.S. Udaykumar, W. Shyy, An accurate Cartesian Grid Method for Viscous Incompressible Flows with Complex Immersed Boundaries, Journal of Computational Physics 156 (1999) 209-240.

37

Figure 1: Three cases for face-centered variable interpolation

38

Figure 2: Fluxes at edges of irregular conducting cells

39

(a) a uniform mesh

(b) a non-uniform mesh that the center region are refined with an additional level

40 (c) a non-uniform mesh that the center region are refined with two additional levels Figure 3: Convergence of relative errors in electric potential evolution when the consistent scheme is implemented on coarse/fine mesh boundaries. As it shows, even if one or two additional levels are added in local domain(the second picture and third picture respectively), the convergence rate still keeps second order.

(a) a uniform mesh

(b) a non-uniform mesh that the center region are refined with an additional level

41 (c) a non-uniform mesh that the center region are refined with two additional levels Figure 4: Convergence of relative errors in electric potential evolution when non-consistent scheme is applied on coarse/fine mesh boundaries. As it shows, it’s nearly second order when the mesh is uniform but decrease to first order when one or two additional levels are added in local domain(the second picture and third picture respectively) which indicates this scheme provided on coarse/fine mesh boundaries are not consistent.

Figure 5: The geometry of the diverging channel

42

(a) the velocity distribution at the outlet

(b) the electric potential distribution at the outlet Figure 6: Comparation between the analytic solution given by Eq.(39), Eq.(40) and the numerical results at the outlet of the diverging channel

43

(a) convergence rate of electric potential on non-uniform meshes

(b) convergence rate of velocity on non-uniform meshes Figure 7: Computed order of accuracy in calculating electric potential and velocity on nonuniform grids where cells cut by boundary are refined with one additional level

44

(a) Overall arrangement of computa- (b) Local arrangement of computational grids tional grids Figure 8: The adaptive mesh distribution of the the computational domain. The left picture is for the overall view and the right one is for local view.

45

(a) along a centerline perpendicular to Hartmann walls

(b) along a centerline perpendicular to side walls Figure 9: Results for the simulations using consistent and conservative schemes based on AMR framework. The left is the comparison of analytic solution of Shercliff’s case and the numerical solution along the center-line parallel to the magnetic field. The right is the comparison result along the center-line vertical to the magnetic field.

46

(a) velocity profile

(b) current density streamlines Figure 10: Numerical results for Hunt’s fully developed case of Ha=100 with symmetrical wall.

47

(a) along a centerline perpendicular to Hartmann walls

(b) along a centerline perpendicular to side walls Figure 11: Comparison for Hunt’s symmetrical wall case between the analytic solution and the numerical solution.

48

(a) along a centerline perpendicular to Hartmann walls

(b) along a centerline perpendicular to side walls Figure 12: Results for the simulations describing the Hunt’s variant case when conductivity of top wall is as twice manitude as the bottom wall. And Comparison is given between the analytic solution and the numerical solution.

49

(a) along a centerline perpendicular Hartmann walls

(b) along a centerline perpendicular side walls Figure 13: Results for the simulations describing the Hunt’s variant case when thickness of top wall is half of the bottom wall. And Comparison is given between the analytic solution and the numerical solution.

50

(a) with Ha = 18

(b) with Ha = 30 Figure 14: Comparison of the velocity profile between numerical results and Samad’s analytic solution at Ha = 18(top) and Ha = 30(bottom). Two groups of different conductances of c = 1 and c = 10 are both considered, the M type profile is observed when c increases.

51

(a) with cw = 0.1

(b) with cw = 1

(c) with cw = 10

Figure 15: Comparison of the velocity profile between numerical results (green points) and Chang’s asymptotic solutions (red lines) at Ha = 100, wall-conductance ratio are respectively 0.1, 1, 10. The results match well in the core region but overspeed zone appears near the wall in numerical simulation while the asymptotic approximation profile keeps flat.

Figure 16: Numerical results of different wall-conductances at Ha = 100, which clearly shows when Ha number increases, M type profile is more obvious.

52

Figure 17: Configuration of the study.

53

(a) Nx = 0 based on Grid I

(b) Nx = 0 based on Grid II

(c) Nx = 0.2 based on Grid I

(d) Nx = 0.2 based on Grid II

(e) Nx = 1 based on Grid I

(f) Nx = 1 based on Grid II

Figure 18: The evolution of Cd (red line) and Cl (green line) with different magnetic field. The left three pictures are based on Grid I, while the right three on Grid II.

54