In this paper, a localized integrated radial basis function-based differential quadrature (iRBF-DQ) method to solve steady, magnetohydrodynamic (MHD) duct flows is presented. Local iRBF-DQ method is a truly meshless and efficient method, which discretizes any derivative at a knot by weighted linear sum of functional values at its nearby nodes. The integrated RBF based approaches are more stable than conventional meshless procedure ones. The high Hartmann numbers (high magnetic field) MHD problems apply to the design of cooling systems with liquid metals for a thermal nuclear fusion blanket. We present results for Hartmann number up to${10}^5$ with fully insulating, partially insulating, and partially conducting walls, with the shapes of rectangular, circular or arbitrary cross sections. The results show that by using local iRBF-DQ method we can obtain accurate and stable approximations for the velocity and induced magnetic field at the range of Hartmann numbers of ${10}^2\le \text{M}\le {10}^5\mathrm{.}$ aa

Keywords: meshless, iRBF collocation method, local iRBF-DQ method, MHD, Hartmann number

MHD, magnetohydrodynamic; FDM, finite difference method; AFEM, analytical finite element method; RBFs, radial basis functions; PDEs, partial differential equations

Magnetohydrodynamics (MHD) is the study of the interaction between moving, conducting fluids and magnetic fields. In this paper, the discussion is restricted to incompressible, viscous and electrically conducting fluids through pipes and ducts under an imposed uniform, oblique magnetic field. The MHD issues have become increasingly important because of the practical engineering applications such as the design of cooling systems with liquid metals for a thermal nuclear fusion blanket. Blankets which rely only on the heat transverse ability of liquid metal are known as self-cooled liquid metal blankets. There are numbers of amount of data and papers concerning with the physical phenomena of liquid metal flows in the absence of a magnetic field. However, in a fusion reactor environment such as the first wall and the blankets, a strong magnetic field (high Hartmann numbers $\text{M}\mathrm{<mo><</mo>1000}$ Molokov S & Buhler L^1,2 is necessary and confines the liquid-metal (electrically conducting) flows. The hydraulic and thermal behaviors of liquid metal flows in the existence of a magnetic field are quite diff erent fr om the behaviors of liquid metal flows without a magnetic field.^3,4 Furthermore, the exact solutions of the MHD problems can only be obtained for some special cases.^5,6 Therefore, it is important to explore a stable, accurate, and eff ective numerical method to obtain the approximate solutions of MHD problems.

Currently, researchers have investigated the two-dimensional MHD problem using several numerical methods. Singh & Lal^7,8 solved MHD flows through pipes of triangular cross-section for small values of Hartmann number using finite difference method (FDM). Since it is hard to fit the arbitrary cross-section ofthe channel, the finite element method (FEM) for Hartmann number less than 10 was presented in.^9-11 Tezer-Sezgin & Koksal¹² used standard FEM with linear and quadratic elements for Hartmann numbers up to 100. Aft er that, Demendy and Nagy obtained the numerical solutions for Hartmann number $\text{M}\mathrm{<mo><</mo>1000}$ using analytical finite element method (AFEM).¹³ Furthermore, the boundary element method (BEM),¹⁴ and the fundamental solution method¹⁵ also have been used to obtain approximate solutions for MHD flow problems. To our best knowledge, the numerical results for MHD flow problems are mostly restricted to the moderate Hartmann numbers$\text{M}\mathrm{<mo><</mo>1000}$.

However, Nesliturk & Tezer-Sezgin¹⁶ have obtained the results for Hartmann numbers up to ${10}^5$ is only for one case, a rectangular duct with partly conducting wall, which is computationally very expensive, due to the need of large memory and time consuming in building finer meshes. The difficulty for high Hartmann numbers is analogous to the advectiondiff usion equation. When advection process dominates diffusion, the physical values near the boundary would suddenly decrease to zero to satisfy boundary conditions and form very thin boundary layers. Radial basis functions (RBFs) have been originally developed for scattered data approximation, especially for higher dimensional problems. In recent years, RBFs have been further applied to solve numerical partial diff erential equations (PDEs). This new development is very useful due to the fact that this method is a truly mesh or grid free technique. Kansa¹⁷ initially published a series of papers in this area, and then Franke & Schaback^18,19 were able to give a convergence proof and error bounds of numerical approach. In a comparative study, Franke²⁰ surveyed an extensive number of techniques for interpolation/approximation. He concluded that overall the Multiquadrics (MQ) and Thin Plate Splines (TPS) are the most favorable techniques for scattered data approximation. Nevertheless, the accuracy of MQ depends on a shape parameter which is still an open research topic. Hence most applications of MQ used experimental tuning parameters or expensive optimization techniques to evaluate the optimum shape parameter.²¹ In this work, the integrated TPS is used as it combines good accuracy without the additional burden of computing a shape parameter.

Furthermore, integrated TPS is based on sound mathematical theory Duchon J²² whereas for $\text{MQ}$ , although it works well, its construction theory is yet to be established.²³ Furthermore, the numerical schemes based on RBFs to solve PDEs proposed up to date have one more feature concerned, which are based on function approximation. That is to say, these methods directly substitute the expression of function approximation by RBFs to a PDE, and aft er that we replace the dependent variables into the coefficients of function approximation. The process is very complex, especially for nonlinear problems. Later, Wu & Shu²⁴ proposed a new approach using the so called RBF-DQ method, which focuses on the derivative approximation through the diff erential quadrature (DQ) methodology. The basic idea of the DQ method is that the derivatives of unknown function can be approximated as the function values at a group of nodes. Unfortunately, full exploitation of the feature of RBFs-based methods causes the ill-conditioned coefficient matrix as the rank increases. The RBF collocation method and RBF-DQ method also suff ered this kind of problem. To circumvent this difficulty, Kansa²⁵ suggested substituting the global solvers into block partitioning and LU decomposition schemes for huge simulation problems. Shu & Yeo²⁶ further proposed the local RBF-DQ method using local supporting nodes. Recently, Tran-Cong and his collaborators,^27-31 further proposed the indirect radial basis function networks (IRBFN). They found that the IRBFN are more stable than the conventional RBF method. In the IRBFN approach, the RBF approximation is applied to a specific derivative of the solution function, and thus the solution function can be derived by integrating the derivative with the RBF expression. Consequently, the expressions of IRBFN are derived from integration of conventional RBFs. Mai-Duy & Tran-Cong^27,28 proved that this procedure can have better approximation for the targeted derivative. Later, based on the idea of IRBFN, Shu & Wu³² proposed the local iRBF-DQ method to improve its performance.

The present paper is an application of techniques mentioned above, namely the iRBF collocation method and the local iRBF-DQ method, to the MHD duct flows. And its aim is to establish an accurate numerical approximation to the solutions of benchmark problems of MHD duct flows for the range of Hartmann number${10}^2\le \text{M}\le {10}^5$ . In general, the important range of Hartmann number in engineering practice is ${10}^2\le \text{M}\le {10}^5$ . Since these two meshless methods are based on RBFs, they diff er fr om traditionally mesh-dependent numerical methods such as FEM with the need of meshing.

MHD equations derivation
The discussed geometry of cross section of blankets is restricted to the directions of imposed magnetic field lines which are orthogonal to the direction of flow velocity. We now consider a steady, laminar, fully developed flow with incompressible, viscous and finite conducting fluid in a duct which is infinitely long in $\text{z}$ -direction. Hence, we can assume this is a two-dimensional duct flow along $\text{x}$ and $\text{y}$  directions.6

There is an imposed uniform, magnetic field $H_0$  outside the field forming an angle $\alpha$  with the $\text{y}$ -axis. The fluid velocity $u$  is parallel to $\text{z}$ -axis, the half channel width is $a$ , and the origin is located at the center ofthe channel (Figure 1).

The Maxwell’s equations in steady state can be expressed as follows:

$\nabla \times \text{E}\mathrm{<mo>=</mo>0}$,(1)

$\nabla \times \text{H}\mathrm{<mo>=</mo>}\text{j}$ , (2)

$$\text{B}\mathrm{<mo>=</mo>}\mu \text{H}$$. (3)

$\text{j}\mathrm{<mo>=</mo>}\sigma \mathrm{<mo\; stretchy="false">(</mo>}\text{E}+\text{v}\times \text{B}\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}\sigma \mathrm{<mo\; stretchy="false">(</mo>}\text{E}+\mu \text{v}\times \text{H}\mathrm{<mo\; stretchy="false">)</mo>}$. (4) from (1) and (2), we obtain

$\nabla \times \text{j}\mathrm{<mo>=</mo>}\nabla \times \nabla \times \text{H}\mathrm{<mo>=</mo>}\sigma \mu \nabla \times \mathrm{<mo\; stretchy="false">(</mo>}\text{v}\times \text{H}\mathrm{<mo\; stretchy="false">)</mo>}$. (5 )

Through vector analysis, (5) can be written as follows

$\nabla \mathrm{<mo\; stretchy="false">(</mo>}\nabla \cdot \text{H}\mathrm{<mo\; stretchy="false">)</mo>}-{\nabla }^2\text{H}\mathrm{<mo>=</mo>}\mu \sigma \mathrm{<mo\; stretchy="false">[</mo><mo\; stretchy="false">(</mo>}\text{H}\cdot \nabla \mathrm{<mo\; stretchy="false">)</mo>}\text{v}-\mathrm{<mo\; stretchy="false">(</mo>}\text{v}\cdot \nabla \mathrm{<mo\; stretchy="false">)</mo>}\text{H}+\mathrm{<mo\; stretchy="false">(</mo>}\nabla \cdot \text{H}\mathrm{<mo\; stretchy="false">)</mo>}\text{v}-\mathrm{<mo\; stretchy="false">(</mo>}\nabla \cdot \text{v}\mathrm{<mo\; stretchy="false">)</mo>}\text{H}\mathrm{<mo\; stretchy="false">]</mo>.}$  (6)

Adding Lorentz force (see Appendix A for more details) to the Navier-Stokes equations in a steady state problem, we have the following form:

$\mu \text{j}\times \text{H}-\nabla p+\eta {\nabla }^2\text{v}\mathrm{<mo>=</mo>}\rho \mathrm{<mo\; stretchy="false">(</mo>}\text{v}\cdot \nabla \mathrm{<mo\; stretchy="false">)</mo>}\text{v}$ . (8)

From (2) and using a vector identity, we have

$\text{j}\times \text{H}\mathrm{<mo>=</mo><mo\; stretchy="false">(</mo>}\text{H}\cdot \nabla \mathrm{<mo\; stretchy="false">)</mo>}\text{H}-\frac{1}{2}\nabla \mathrm{<mo>|</mo>}\text{H}{\mathrm{<mo>|</mo>}}^2$ (9)

Substituting (8) into ($\mathrm{7<mo\; stretchy="false">)</mo>}$  , we will obtain the modified Navier-Stokes equations as follows:

$\eta {\nabla }^2\text{v}+\mu \mathrm{<mo\; stretchy="false">(</mo>}\text{H}\cdot \nabla \mathrm{<mo\; stretchy="false">)</mo>}\text{H}\mathrm{<mo>=</mo>}\rho \mathrm{<mo\; stretchy="false">(</mo>}\text{v}\cdot \nabla \mathrm{<mo\; stretchy="false">)</mo>}\text{v}+\nabla \mathrm{<mo\; stretchy="false">(</mo>}p+\frac{\mu }{2}\mathrm{<mo>|</mo>}\text{H}{\mathrm{<mo>|</mo>}}^2\mathrm{<mo\; stretchy="false">)</mo>}$ . (10)

Notice that all physical values except pressure gradient are unchanged in the $\text{z}$ -direction. Thus the fluid velocity $\text{v}$ , the magnetic field$\text{H}$ , and the pressure gradient $\nabla p$  can be expressed in the following form:

$\text{v}\mathrm{<mo>=</mo><mo\; stretchy="false">(</mo>0,0,}{\text{v}}_{\text{z}}\mathrm{<mo\; stretchy="false">)</mo>,}\text{H}\mathrm{<mo>=</mo><mo\; stretchy="false">(</mo>}{H}_x\mathrm{,}{H}_y\mathrm{,}\frac{B_Z}{\mu }\mathrm{<mo\; stretchy="false">)</mo>,}\nabla p\mathrm{<mo>=</mo>}-k\eta$ , (11)

$\{\begin{array}{l}{H}_X\mathrm{<mo>=</mo>}{H}_{0}sin\alpha \mathrm{,}\hfill \\ H_y\mathrm{<mo>=</mo>}{H}_{0}cos\alpha \mathrm{,}\hfill \end{array}$(12)

$\frac{1}{\mu }{\nabla }^2{B}_{\text{z}\mathrm{,}\to }+\mu \sigma H_X\frac{\partial {\text{v}}_{\text{z}}}{\partial x}+\mu \sigma H_y\frac{\partial {\text{v}}_{\text{z}}}{\partial y}\mathrm{<mo>=</mo>0}$(13)

$\eta {\nabla }^2{\text{v}}_{\text{z}}+H_X\frac{\partial B_{-}-}{\partial x}+\mu H_y\frac{\partial B_{-}-}{\partial y}+k\eta \mathrm{<mo>=</mo>0.}$ (14)

It follows that we can substitute some non-dimensional variables [33]

${\text{v}}_{\text{z}}\mathrm{<mo>=</mo>}\frac{H_0}{\sqrt{\sigma \eta }}{\text{v}}_{\text{z}}^{\mathrm{<mo>*</mo>}}\mathrm{,}{B}_{\underset{\_}{\text{z}}}\mathrm{<mo>=</mo>}\mu H_0{B}_{\underset{\_}{\text{z}}}^{\mathrm{<mo>*</mo>}}\mathrm{,<mo\; stretchy="false">(</mo>}x\mathrm{,}y\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}L\mathrm{<mo\; stretchy="false">(</mo>}{x}^{\mathrm{<mo>*</mo>}}\mathrm{,}{y}^{\mathrm{<mo>*</mo>}}\mathrm{<mo\; stretchy="false">)</mo>,}k\mathrm{<mo>=</mo>}\frac{H_0}{L^2\sqrt{\sigma \eta }}{k}^{\mathrm{<mo>*</mo>}}$ (15)

${\nabla }^2{B}_{-}^{\mathrm{<mo>*</mo>}}-+\mu H_{0}a\sqrt{\frac{\sigma }{\eta }}sin\theta \frac{\partial {\text{v}}_{\text{z}}^{\mathrm{<mo>*</mo>}}}{\partial x}+\mu H_{0}a\sqrt{\frac{\sigma }{\eta }}cos\theta \frac{{\partial }_{{\text{V}}_{\text{z}}^{\mathrm{<mo>*</mo>}}}}{\partial y}\mathrm{<mo>=</mo>0}$ (16)

${\nabla }^2{\text{v}}_{\text{z}}^{\mathrm{<mo>*</mo>}}+\mu H_{0}a\sqrt{\frac{\sigma }{\eta }}sin\alpha \frac{\partial B_{--}^{\mathrm{<mo>*</mo>}}}{\partial x}+\mu H_{0}a\sqrt{\frac{\sigma }{\eta }}cos\alpha \frac{\partial B_{--}^{\mathrm{<mo>*</mo>}}}{\partial y}+k^{\mathrm{<mo>*</mo>}}\mathrm{<mo>=</mo>0}$ (17)

It is noted that there is a convenient independent dimensionless number, which is called Hartmann number,$M$; and is defined by as follows:

$$M\mathrm{<mo>=</mo>}\mu H_{0}a\sqrt{\sigma }\mathrm{<mo>/</mo>}\eta$$ ,      (18)

$$\{\begin{array}{l}{M}_X\mathrm{<mo>=</mo>}Msin\alpha \mathrm{,}\hfill \\ M_y\mathrm{<mo>=</mo>}Mcos\alpha \mathrm{,}\hfill \end{array}$$ (19)

$${\nabla }^2{B}_Z^{\mathrm{<mo>*</mo>}}+M_x\frac{\partial {\text{v}}_{\text{z}}^{\mathrm{<mo>*</mo>}}}{\partial x}+M_y\frac{\partial {\text{v}}_{\text{z}}^{\mathrm{<mo>*</mo>}}}{\partial y}\mathrm{<mo>=</mo>0}$$ , (20)

$${\nabla }_{{\text{V}}_{\text{z}}^{\mathrm{<mo>*</mo>}}}^2+M_x\frac{\partial B_Z^{\mathrm{<mo>*</mo>}}}{\partial x}+M_y\frac{\partial B_Z^{\mathrm{<mo>*</mo>}}}{\partial y}+k^{\mathrm{<mo>*</mo>}}\mathrm{<mo>=</mo>0}$$ (21)

where $k^{\mathrm{<mo>*</mo>}}\mathrm{<mo>=</mo>1}$  is assumed.

To facilitate easy reading, we omit the symbol   and make   and   replace  and $B_{\text{z}}^{\mathrm{<mo>*</mo>}}$ respectively in the following sections.

Boundary conditions

We consider the following simplified dimensionless coupled diff erential equations

$\Omega \mathrm{,}$, in $\Omega \mathrm{,}$ (22)

$\Omega \mathrm{,}$ , in$\Omega \mathrm{,}$,   (23)

$\{\begin{array}{ll}u\mathrm{<mo>=</mo>0,}\hfill & \text{on}\partial \Omega \mathrm{,}\hfill \\ B\mathrm{<mo>=</mo>0,}\hfill & \text{on}{\Gamma }_1\mathrm{,}\hfill \\ \frac{\partial B}{\partial n}\mathrm{<mo>=</mo>0,}\hfill & \text{on}{\Gamma }_2\mathrm{,}\hfill \end{array}$

Figure 1Duct problem setting.

Figure 2A 2D duct flow.

Decoupled governing equations

The MHD equations are coupled convection-diff usion equations. In general, solving coupled equations will consume more computational time and memory space than solving decoupled equations. Considering a rectangular duct with insulating walls, we can decouple the MHD equations.

$U_1\mathrm{<mo>=</mo>}u+B$ ; $U_2\mathrm{<mo>=</mo>}u-B$,    (24)

$${\nabla }^2{U}_1+M_x\frac{\partial U_1}{\partial x}+M_y\frac{\partial U_1}{\partial y}\mathrm{<mo>=</mo>}-1$$ , in$\Omega \mathrm{,}$,   (25)

<${\nabla }^2{U}_2-M_x\frac{\partial U_2}{\partial x}-M_y\frac{\partial U_2}{\partial y}\mathrm{<mo>=</mo>}-1$ , in$$\Omega \mathrm{,}$$(26)

$U_1\mathrm{<mo>=</mo>}{U}_2\mathrm{<mo>=</mo>0}$ on$\Gamma$ >.           (27)

Radial basis functions
A m-th order integrated TPS is defined as

$\phi \mathrm{<mo\; stretchy="false">(</mo>}{r}_j\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{\nabla }^2\Phi \mathrm{<mo\; stretchy="false">(</mo>}{r}_j\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{r}_j^{2m}log\mathrm{<mo\; stretchy="false">(</mo>}{r}_j\mathrm{<mo\; stretchy="false">)</mo>}$,      (28)

where $r_j\mathrm{<mo>=</mo>}P\text{x}-{\text{x}}_{j}P$ is the Euclidean norm. Since $\varphi$ given by (28) is $$C^{2m-1}$$ , a higher-order integrated Thin Plate Splines must be used.

Local $i\sqrt{BF}-DQ$   method

The basic idea of iRBF-DQ method was originated fr om the traditional DQ method. The basic approach of DQ method implies that the partial derivative of a continuous function at a reference node can be approximated by a weighted linear sum of function values at all discrete points within its support. Thus, DQ approximation of m-th order derivate with respect to$x$  -direction at a reference node ${\text{x}}_i$ can be expressed as

$\frac{{\partial }^{p}f}{\partial x^p}\text{x}\mathrm{<mo>=</mo>}{\text{x}}_j\mathrm{<mo>=</mo>}{\displaystyle {\sum }_{j\mathrm{<mo>=</mo>1}}}{\mathrm{<mo\; stretchy="false">(</mo>}wx\mathrm{<mo\; stretchy="false">)</mo>}}_{ij}^{\mathrm{<mo\; stretchy="false">(</mo>}p\mathrm{<mo\; stretchy="false">)</mo>}}f\mathrm{<mo\; stretchy="false">(</mo>}{\text{x}}_j\mathrm{<mo\; stretchy="false">)</mo>}N\mathrm{,}j\mathrm{<mo>=</mo>1,2,}\dots \mathrm{,}N$ ,     (29)

${\frac{{\partial }^{p}f}{\partial x^p}|}_{\text{x}\mathrm{<mo>=</mo>}{\text{x}}_j}\mathrm{<mo>=</mo>}\sum _{k\mathrm{<mo>=</mo>1}}^{N_l}{\mathrm{<mo\; stretchy="false">(</mo>}wx\mathrm{<mo\; stretchy="false">)</mo>}}_{jk}^{\mathrm{<mo\; stretchy="false">(</mo>}p\mathrm{<mo\; stretchy="false">)</mo>}}f\mathrm{<mo\; stretchy="false">(</mo>}{r}_{jk}\mathrm{<mo\; stretchy="false">)</mo>,}j\mathrm{<mo>=</mo>1,2,}\dots \mathrm{,}{N}_l\mathrm{,}k\mathrm{<mo>=</mo>1,2,}\dots \mathrm{,}{N}_l$     (30)

For the local iRBF-DQ method, $m\mathrm{<mo>=</mo>}l$ is chosen (i.e. integrated first-order TPS) as the basis functions. Integrated TPS can be obtained as follows

$\Phi \mathrm{<mo\; stretchy="false">(</mo>}{r}_k\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}\frac{{\mathrm{<mo\; stretchy="false">(</mo>}{r}_k\mathrm{<mo\; stretchy="false">)</mo>}}^4\mathrm{<mo\; stretchy="false">(</mo>2}log{r}_k-\mathrm{1<mo\; stretchy="false">)</mo>}}{32}$(31)

$\frac{\partial {\Phi }_k}{\partial x}\mathrm{<mo\; stretchy="false">(</mo>}{\text{x}}_j\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{\displaystyle {\sum }_{k\mathrm{<mo>=</mo>1}}^{N_l}}{\mathrm{<mo\; stretchy="false">(</mo>}wx\mathrm{<mo\; stretchy="false">)</mo>}}_{jk}^{\mathrm{<mo\; stretchy="false">(</mo>1<mo\; stretchy="false">)</mo>}}{\Phi }_{jk}\mathrm{<mo\; stretchy="false">(</mo>}{r}_{jk}\mathrm{<mo\; stretchy="false">)</mo>,}j\mathrm{<mo>=</mo>1,2,}\dots \mathrm{,}{N}_T\mathrm{,}k\mathrm{<mo>=</mo>1,2,}\dots \mathrm{,}{N}_T$       (32)

$\left[\begin{array}{l}\frac{\partial {\Phi }_1}{\partial x}\\ \frac{\partial {\Phi }_2}{\begin{array}{l}\partial x\\ \frac{\partial {\Phi }_{N_l}}{\partial x}\end{array}}\end{array}\right]=\left[\begin{array}{l}{\Phi }_{11}{\Phi }_{12}\cdots {\Phi }_{1N_l}\\ {\Phi }_{21}{\Phi }_{22\cdots }{\Phi }_{2N_l}\\ ⋮⋮\ddots ⋮\\ {\Phi }_{N_{l}1}{\Phi }_{N_{l}2}\cdots {\Phi }_{N_l{N}_l}\\ \end{array}\right]\left\{\begin{array}{l}{\mathrm{<mo\; stretchy="false">(</mo>}wx\mathrm{<mo\; stretchy="false">)</mo>}}_{\text{l}}^{\mathrm{<mo\; stretchy="false">(</mo>}\text{l}\mathrm{<mo\; stretchy="false">)</mo>}}\hfill \\ {\mathrm{<mo\; stretchy="false">(</mo>}wx\mathrm{<mo\; stretchy="false">)</mo>}}_2^{\mathrm{<mo\; stretchy="false">(</mo>}\text{l}\mathrm{<mo\; stretchy="false">)</mo>}}\hfill \\ ⋮\hfill \\ {\mathrm{<mo\; stretchy="false">(</mo>}wx\mathrm{<mo\; stretchy="false">)</mo>}}_{N_1}^{\mathrm{<mo\; stretchy="false">(</mo>}\text{l}\mathrm{<mo\; stretchy="false">)</mo>}}\hfill \end{array}\right\}$ (33)

$\{\begin{array}{ll}{\displaystyle \sum _{i\mathrm{<mo>=</mo>1}}^N}{\alpha }_{ik}^{U_1}{U}_1{}_{\left(j\right)}\mathrm{<mo>=</mo>}-\mathrm{1,}i\mathrm{<mo>=</mo>1,2,}\dots \mathrm{,}N\mathrm{,}j\mathrm{<mo>=</mo>1,2,}\dots \mathrm{,}{N}_T\mathrm{,}k\mathrm{<mo>=</mo>1,2,}\dots \mathrm{,}{N}_T\hfill & \text{in}\Omega \mathrm{,}\hfill \\ {\displaystyle \sum _{i\mathrm{<mo>=</mo>1}}^N}{\beta }_{ik}^{U_2}{U}_{2_{\left(j\right)}}\mathrm{<mo>=</mo>}-\mathrm{1,}i\mathrm{<mo>=</mo>1,2,}\dots \mathrm{,}N\mathrm{,}j\mathrm{<mo>=</mo>1,2,}\dots \mathrm{,}{N}_T\mathrm{,}k\mathrm{<mo>=</mo>1,2,}\dots \mathrm{,}{N}_T\hfill & \text{in}\Omega \mathrm{,}\hfill \end{array}$  (34)

$U_1\mathrm{<mo>=</mo>}{U}_2\mathrm{<mo>=</mo>0}$ , on$\Gamma$ , (35)

$\{\begin{array}{l}{a}_{ik}^{U_1}\mathrm{<mo>=</mo>}{\displaystyle \sum _j^N}\mathrm{<mo\; stretchy="false">(</mo><mo\; stretchy="false">(</mo>}wx{\mathrm{<mo\; stretchy="false">)</mo>}}_{jk}^{\mathrm{<mo\; stretchy="false">(</mo>2<mo\; stretchy="false">)</mo>}}+{\mathrm{<mo\; stretchy="false">(</mo>}wy\mathrm{<mo\; stretchy="false">)</mo>}}_j^{\mathrm{<mo\; stretchy="false">(</mo>2<mo\; stretchy="false">)</mo>}}k+M_X{\mathrm{<mo\; stretchy="false">(</mo>}wx\mathrm{<mo\; stretchy="false">)</mo>}}_{jk}^{\mathrm{<mo\; stretchy="false">(</mo>1<mo\; stretchy="false">)</mo>}}+M_y{\mathrm{<mo\; stretchy="false">(</mo>}wy\mathrm{<mo\; stretchy="false">)</mo>}}_j^{\mathrm{<mo\; stretchy="false">(</mo>1<mo\; stretchy="false">)</mo>}}k\mathrm{<mo\; stretchy="false">)</mo>}\mathrm{,}k\mathrm{<mo>=</mo>1,2,}\dots \mathrm{,}{N}_l\hfill \\ {\beta }_{ik}^{U_2}\mathrm{<mo>=</mo>}{\displaystyle \sum _{i\mathrm{<mo>=</mo>1}}^{N_L}}\mathrm{<mo\; stretchy="false">(</mo><mo\; stretchy="false">(</mo>}wx{\mathrm{<mo\; stretchy="false">)</mo>}}_{jk}^{\mathrm{<mo\; stretchy="false">(</mo>2<mo\; stretchy="false">)</mo>}}+{\mathrm{<mo\; stretchy="false">(</mo>}wy\mathrm{<mo\; stretchy="false">)</mo>}}_j^{\mathrm{<mo\; stretchy="false">(</mo>2<mo\; stretchy="false">)</mo>}}k-M_X{\mathrm{<mo\; stretchy="false">(</mo>}wx\mathrm{<mo\; stretchy="false">)</mo>}}_{jk}^{\mathrm{<mo\; stretchy="false">(</mo>1<mo\; stretchy="false">)</mo>}}-M_y{\mathrm{<mo\; stretchy="false">(</mo>}wy\mathrm{<mo\; stretchy="false">)</mo>}}_j^{\mathrm{<mo\; stretchy="false">(</mo>1<mo\; stretchy="false">)</mo>}}k\mathrm{<mo\; stretchy="false">)</mo>}\mathrm{,}k\mathrm{<mo>=</mo>1,2,}\dots \mathrm{,}{N}_l\hfill \end{array}$  (36)

$\left\{\begin{array}{l}{U}_{\mathrm{1<mo\; stretchy="false">(</mo>1<mo\; stretchy="false">)</mo>}}\hfill \\ U_{\mathrm{1<mo\; stretchy="false">(</mo>2<mo\; stretchy="false">)</mo>}}\hfill \\ ⋮\hfill \\ U_{\mathrm{1<mo\; stretchy="false">(</mo>}B\text{C}\mathrm{<mo\; stretchy="false">)</mo>}}\hfill \\ U_{\mathrm{1<mo\; stretchy="false">(</mo>}N\mathrm{<mo\; stretchy="false">)</mo>}}\hfill \end{array}\right\}\mathrm{<mo>=</mo>}\left[\begin{array}{l}{\alpha }_{1,1}{\alpha }_{1,2}\cdots {\alpha }_{1,N_l}00\\ {\alpha }_{1,2{\alpha }_{2,2\cdots }{\alpha }_{2,N_l}00}\\ ⋮⋮\ddots \ddots ⋮⋮\\ 001000\\ 00{\alpha }_{N,1{\alpha }_{N,2\cdots }}{\alpha }_{N,N_l}\end{array}\right]\left\{\begin{array}{l}-1\hfill \\ -1\hfill \\ ⋮\hfill \\ 0\hfill \\ -1\hfill \end{array}\right\}$ (37)

$\left\{\begin{array}{l}{U}_{\mathrm{2<mo\; stretchy="false">(</mo>1<mo\; stretchy="false">)</mo>}}\hfill \\ U_{\mathrm{2<mo\; stretchy="false">(</mo>2<mo\; stretchy="false">)</mo>}}\hfill \\ ⋮\hfill \\ U_{\mathrm{2<mo\; stretchy="false">(</mo>}B\text{C}\mathrm{<mo\; stretchy="false">)</mo>}}\hfill \\ U_{\mathrm{2<mo\; stretchy="false">(</mo>}N\mathrm{<mo\; stretchy="false">)</mo>}}\hfill \end{array}\right\}\mathrm{<mo>=</mo>}\left\{\begin{array}{llllllll}{\beta }_{\mathrm{1,1}}\hfill & {\beta }_{\mathrm{1,2}}\hfill & \cdots \hfill & {\beta }_{\mathrm{1,}{N}_l}\hfill & 0\hfill & \hfill & \hfill & 0\hfill \\ {\beta }_{\mathrm{2,1}}\hfill & {\beta }_{\mathrm{2,2}}\hfill & \cdots \hfill & \cdots \hfill & {\beta }_{\mathrm{2,}{N}_l}0\hfill & \hfill & \hfill & 0\hfill \\ ⋮\hfill & ⋮\hfill & \hfill & \ddots \hfill & ⋮\hfill & \ddots \hfill & \hfill & ⋮\hfill \\ 0\hfill & 0\hfill & \hfill & \hfill & 10\hfill & 0\hfill & \ddots \hfill & 0\hfill \\ 0\hfill & 0\hfill & \hfill & \hfill & {\beta }_{N\mathrm{,,1}}\hfill & {\beta }_{N\mathrm{,,2}}\hfill & \cdots \hfill & {\beta }_{N\mathrm{,,}{N}_l}\hfill \end{array}\right\}\left\{\begin{array}{l}-1\hfill \\ -1\hfill \\ ⋮\hfill \\ 0\hfill \\ -1\hfill \end{array}\right\}$ (38)

$$u\mathrm{<mo>=</mo>}\frac{U_1+U_2}{2}\mathrm{<mo>;</mo>}B\mathrm{<mo>=</mo>}\frac{U_1-U_2}{2}$$ (39)

$i\sqrt{BF}$   collocation method

$U_1\mathrm{<mo\; stretchy="false">(</mo>}{\text{x}}_j\mathrm{<mo\; stretchy="false">)</mo>}\equiv {\displaystyle {\sum }_{j\mathrm{<mo>=</mo>1}}^N}{\alpha }_j^{U_1}\Phi \mathrm{<mo\; stretchy="false">(</mo>}{r}_{ij}\mathrm{<mo\; stretchy="false">)</mo>}$ , for ${\text{x}}_i\in \Omega$ (40)

$U_2\mathrm{<mo\; stretchy="false">(</mo>}{\text{x}}_i\mathrm{<mo\; stretchy="false">)</mo>}\equiv {\displaystyle {\sum }_{j\mathrm{<mo>=</mo>1}}^N}{\alpha }_j^{U_2}\Phi \mathrm{<mo\; stretchy="false">(</mo>}{r}_{ij}\mathrm{<mo\; stretchy="false">)</mo>}$ , for <${\text{x}}_i\in \Omega$ (41)

$\{\begin{array}{l}{L}_1{U}_1\mathrm{<mo>=</mo>}{f}_1\hfill \\ L_2{U}_2\mathrm{<mo>=</mo>}{f}_2\mathrm{<mo>\text{'}</mo>}\hfill \end{array}$ (42)

$L_1\mathrm{<mo>=</mo><mo\; stretchy="false">[</mo>}{\nabla }^2+M_x\frac{\partial }{\partial x}+M_y\frac{\partial }{\partial y}\mathrm{<mo\; stretchy="false">]</mo>}$   (43)

$L_2\mathrm{<mo>=</mo><mo\; stretchy="false">[</mo>}{\nabla }^2-M_x\frac{\partial }{\partial x}-M_y\frac{\partial }{\partial y}\mathrm{<mo\; stretchy="false">]</mo>}$ (44) and

$f_1\mathrm{<mo>=</mo><mo\; stretchy="false">[</mo>}-\mathrm{1<mo\; stretchy="false">]</mo>,}{f}_2\mathrm{<mo>=</mo><mo\; stretchy="false">[</mo>}-\mathrm{1<mo\; stretchy="false">]</mo>}$ (45)

${\displaystyle {\sum }_{j\mathrm{<mo>=</mo>1}}^N}{\alpha }_j^{U_1}\mathrm{<mo\; stretchy="false">[</mo>}{\nabla }^2\Phi \mathrm{<mo\; stretchy="false">(</mo>}{r}_j\mathrm{<mo\; stretchy="false">)</mo>}+M_x\frac{\partial \Phi \mathrm{<mo\; stretchy="false">(</mo>}{r}_j\mathrm{<mo\; stretchy="false">)</mo>}}{\partial x}+M_y\frac{\partial \Phi \mathrm{<mo\; stretchy="false">(</mo>}{r}_j\mathrm{<mo\; stretchy="false">)</mo>}}{\partial y}\mathrm{<mo\; stretchy="false">]</mo><mo>=</mo>}-\mathrm{1,}\text{in}\Omega \mathrm{,}$ (46)

${\displaystyle {\sum }_j^N}{\alpha }_j^{U_2}\mathrm{<mo\; stretchy="false">[</mo>}{\nabla }^2\Phi \mathrm{<mo\; stretchy="false">(</mo>}{r}_j\mathrm{<mo\; stretchy="false">)</mo>}-M_x\frac{\partial \Phi \mathrm{<mo\; stretchy="false">(</mo>}{r}_j\mathrm{.<mo\; stretchy="false">)</mo>}}{\partial x}-M_y\frac{\partial \Phi \mathrm{<mo\; stretchy="false">(</mo>}{r}_j\mathrm{.<mo\; stretchy="false">)</mo>}}{\partial y}\mathrm{<mo\; stretchy="false">]</mo><mo>=</mo>}-1$ , in$\Omega$ ,  (47)

${\displaystyle {\sum }_j^N}{\alpha }_J^{U_1}\Phi \mathrm{<mo\; stretchy="false">(</mo>}{r}_j\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{\displaystyle {\sum }_j^N}{\alpha }_J^{U_2}\Phi \mathrm{<mo\; stretchy="false">(</mo>}{r}_j\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>0}$ on$\Gamma$ ,                   (48)

${\nabla }^2\Phi \mathrm{<mo\; stretchy="false">(</mo>}{r}_j\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}\varphi \mathrm{<mo\; stretchy="false">(</mo>}{r}_j\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{r}_J^{4}log\mathrm{<mo\; stretchy="false">(</mo>}{r}_j\mathrm{<mo\; stretchy="false">)</mo>}$ (49)

$\Phi \mathrm{<mo\; stretchy="false">(</mo>}{r}_j\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}\frac{r_J^6\mathrm{<mo\; stretchy="false">(</mo>3}log{r}_J\mathrm{.}-\mathrm{1<mo\; stretchy="false">)</mo>}}{108}$ (50)

$\frac{\partial \Phi \mathrm{<mo\; stretchy="false">(</mo>}{r}_j\mathrm{<mo\; stretchy="false">)</mo>}}{\partial x}\mathrm{<mo>=</mo>}\frac{x{r}_j^4\mathrm{<mo\; stretchy="false">(</mo>6}log{r}_j-\mathrm{1<mo\; stretchy="false">)</mo>}}{36}$       (51)

$\frac{\partial \Phi \mathrm{<mo\; stretchy="false">(</mo>}{r}_j\mathrm{<mo\; stretchy="false">)</mo>}}{\partial y}\mathrm{<mo>=</mo>}\frac{y{r}_j^4\mathrm{<mo\; stretchy="false">(</mo>6}log{r}_j-\mathrm{1<mo\; stretchy="false">)</mo>}}{36}$ (52)

$\left\{\begin{array}{lll}{L}_1\mathrm{<mo\; stretchy="false">(</mo>}{\Phi }_{11}\mathrm{<mo\; stretchy="false">)</mo>}\hfill & L_1\mathrm{<mo\; stretchy="false">(</mo>}{\Phi }_{12}\mathrm{<mo\; stretchy="false">)</mo>}\hfill & L_1\mathrm{<mo\; stretchy="false">(</mo>}{\Phi }_{1N}\mathrm{<mo\; stretchy="false">)</mo>}\hfill \\ L_1\mathrm{<mo\; stretchy="false">(</mo>}{\Phi }_{21}\mathrm{<mo\; stretchy="false">)</mo>}\hfill & L_1\mathrm{<mo\; stretchy="false">(</mo>}{\Phi }_{22}\mathrm{<mo\; stretchy="false">)</mo>}\hfill & L_1\mathrm{<mo\; stretchy="false">(</mo>}{\Phi }_{2N}\mathrm{<mo\; stretchy="false">)</mo>}\hfill \\ ⋮\hfill & ⋮\hfill & ⋮\hfill \\ {\Phi }_{BC.1}\hfill & {\Phi }_{BC.2}\hfill & {\Phi }_{BCN}\hfill \\ L_1\mathrm{<mo\; stretchy="false">(</mo>}{\Phi }_{N1}\mathrm{<mo\; stretchy="false">)</mo>}\hfill & L_1\mathrm{<mo\; stretchy="false">(</mo>}{\Phi }_{N2}\mathrm{<mo\; stretchy="false">)</mo>}\hfill & L_1\mathrm{<mo\; stretchy="false">(</mo>}{\Phi }_{NN}\mathrm{<mo\; stretchy="false">)</mo>}\hfill \end{array}\right\}\left\{\begin{array}{l}{\alpha }_1^{U_1}\hfill \\ {\alpha }_2^{U_1}\hfill \\ ⋮\hfill \\ {\alpha }_N^{U_1}\hfill \end{array}\right\}\mathrm{<mo>=</mo>}\left\{\begin{array}{l}-1\hfill \\ -1\hfill \\ ⋮\hfill \\ 0\hfill \\ -1\hfill \end{array}\right\}$ (59)

$\left\{\begin{array}{lll}{L}_2\mathrm{<mo\; stretchy="false">(</mo>}{\Phi }_{11}\mathrm{<mo\; stretchy="false">)</mo>}\hfill & L_2\mathrm{<mo\; stretchy="false">(</mo>}{\Phi }_{12}\mathrm{<mo\; stretchy="false">)</mo>}\hfill & L_2\mathrm{<mo\; stretchy="false">(</mo>}{\Phi }_{1N}\mathrm{<mo\; stretchy="false">)</mo>}\hfill \\ L_2\mathrm{<mo\; stretchy="false">(</mo>}{\Phi }_{21}\mathrm{<mo\; stretchy="false">)</mo>}\hfill & L_2\mathrm{<mo\; stretchy="false">(</mo>}{\Phi }_{22}\mathrm{<mo\; stretchy="false">)</mo>}\hfill & L_2\mathrm{<mo\; stretchy="false">(</mo>}{\Phi }_{2N}\mathrm{<mo\; stretchy="false">)</mo>}\hfill \\ ⋮\hfill & ⋮\hfill & ⋮\hfill \\ {\Phi }_{BC.1}\hfill & {\Phi }_{BC.2}\hfill & {\Phi }_{BCN}\hfill \\ L_2\mathrm{<mo\; stretchy="false">(</mo>}{\Phi }_{N1}\mathrm{<mo\; stretchy="false">)</mo>}\hfill & L_2\mathrm{<mo\; stretchy="false">(</mo>}{\Phi }_{N2}\mathrm{<mo\; stretchy="false">)</mo>}\hfill & L_2\mathrm{<mo\; stretchy="false">(</mo>}{\Phi }_{NN}\mathrm{<mo\; stretchy="false">)</mo>}\hfill \end{array}\right\}\left\{\begin{array}{l}{\alpha }_1^{U_2}\hfill \\ {\alpha }_2^{U_2}\hfill \\ ⋮\hfill \\ {\alpha }_N^{U_2}\hfill \end{array}\right\}\mathrm{<mo>=</mo>}\left\{\begin{array}{l}-1\hfill \\ -1\hfill \\ ⋮\hfill \\ 0\hfill \\ -1\hfill \end{array}\right\}$ (60)

$\left\{\begin{array}{l}{\alpha }_1^{U_1}\hfill \\ {\alpha }_2^{U_2}\hfill \\ ⋮\hfill \\ {\alpha }_N^{U_1}\hfill \end{array}\right\}\mathrm{<mo>=</mo>}\left\{\begin{array}{lll}{L}_1\mathrm{<mo\; stretchy="false">(</mo>}{\Phi }_{11}\mathrm{<mo\; stretchy="false">)</mo>}\hfill & L_1\mathrm{<mo\; stretchy="false">(</mo>}{\Phi }_{12}\mathrm{<mo\; stretchy="false">)</mo>}\hfill & L_1\mathrm{<mo\; stretchy="false">(</mo>}{\Phi }_{1N}\mathrm{<mo\; stretchy="false">)</mo>}\hfill \\ L_1\mathrm{<mo\; stretchy="false">(</mo>}{\Phi }_{21}\mathrm{<mo\; stretchy="false">)</mo>}\hfill & L_1\mathrm{<mo\; stretchy="false">(</mo>}{\Phi }_{22}\mathrm{<mo\; stretchy="false">)</mo>}\hfill & L_1\mathrm{<mo\; stretchy="false">(</mo>}{\Phi }_{2N}\mathrm{<mo\; stretchy="false">)</mo>}\hfill \\ ⋮\hfill & ⋮\hfill & ⋮\hfill \\ {\Phi }_{BC.1}\hfill & {\Phi }_{BC.2}\hfill & {\Phi }_{BCN}\hfill \\ L_1\mathrm{<mo\; stretchy="false">(</mo>}{\Phi }_{N1}\mathrm{<mo\; stretchy="false">)</mo>}\hfill & L_1\mathrm{<mo\; stretchy="false">(</mo>}{\Phi }_{N2}\mathrm{<mo\; stretchy="false">)</mo>}\hfill & L_1\mathrm{<mo\; stretchy="false">(</mo>}{\Phi }_{NN}\mathrm{<mo\; stretchy="false">)</mo>}\hfill \end{array}\right\}\left\{\begin{array}{l}-1\hfill \\ -1\hfill \\ ⋮\hfill \\ 0\hfill \\ -1\hfill \end{array}\right\}$ (61)

$\left\{\begin{array}{l}{\alpha }_1^{U_2}\hfill \\ {\alpha }_2^{U_2}\hfill \\ ⋮\hfill \\ {\alpha }_N^{U_2}\hfill \end{array}\right\}\mathrm{<mo>=</mo>}\left\{\begin{array}{lll}{L}_2\mathrm{<mo\; stretchy="false">(</mo>}{\Phi }_{11}\mathrm{<mo\; stretchy="false">)</mo>}\hfill & L_2\mathrm{<mo\; stretchy="false">(</mo>}{\Phi }_{12}\mathrm{<mo\; stretchy="false">)</mo>}\hfill & L_2\mathrm{<mo\; stretchy="false">(</mo>}{\Phi }_{1N}\mathrm{<mo\; stretchy="false">)</mo>}\hfill \\ L_2\mathrm{<mo\; stretchy="false">(</mo>}{\Phi }_{21}\mathrm{<mo\; stretchy="false">)</mo>}\hfill & L_2\mathrm{<mo\; stretchy="false">(</mo>}{\Phi }_{22}\mathrm{<mo\; stretchy="false">)</mo>}\hfill & L_2\mathrm{<mo\; stretchy="false">(</mo>}{\Phi }_{2N}\mathrm{<mo\; stretchy="false">)</mo>}\hfill \\ ⋮\hfill & ⋮\hfill & ⋮\hfill \\ {\Phi }_{BC.1}\hfill & {\Phi }_{BC.2}\hfill & {\Phi }_{BCN}\hfill \\ L_2\mathrm{<mo\; stretchy="false">(</mo>}{\Phi }_{N1}\mathrm{<mo\; stretchy="false">)</mo>}\hfill & L_2\mathrm{<mo\; stretchy="false">(</mo>}{\Phi }_{N2}\mathrm{<mo\; stretchy="false">)</mo>}\hfill & L_2\mathrm{<mo\; stretchy="false">(</mo>}{\Phi }_{NN}\mathrm{<mo\; stretchy="false">)</mo>}\hfill \end{array}\right\}\left\{\begin{array}{l}-1\hfill \\ -1\hfill \\ ⋮\hfill \\ 0\hfill \\ -1\hfill \end{array}\right\}$        (62)

$U_1\mathrm{<mo\; stretchy="false">(</mo>}{\text{x}}_i\mathrm{<mo\; stretchy="false">)</mo>}\equiv {\displaystyle {\sum }_{j\mathrm{<mo>=</mo>1}}^N}{\alpha }_j^{U_1}\Phi \mathrm{<mo\; stretchy="false">(</mo>}{r}_{ij}\mathrm{<mo\; stretchy="false">)</mo>}$ , for ${\text{x}}_i\in \Omega$ (63)

$U_2\mathrm{<mo\; stretchy="false">(</mo>}{\text{x}}_i\mathrm{<mo\; stretchy="false">)</mo>}\equiv {\displaystyle {\sum }_{j\mathrm{<mo>=</mo>1}}^N}{\alpha }_j^{U_2}\Phi \mathrm{<mo\; stretchy="false">(</mo>}{r}_{ij}\mathrm{<mo\; stretchy="false">)</mo>}$ , for ${\text{x}}_i\in \Omega$ (64)

$u\mathrm{<mo>=</mo>}\frac{U_1+U_2}{2}\mathrm{<mo>;</mo>}B\mathrm{<mo>=</mo>}\frac{U_1-U_2}{2}$ (65)

In this section we show the results of the MHD duct flows using the numerical methods introduced in the previous section. It is obvious that RBF collocation method and local RBF-DQ method do not need to mesh the complicated physical problems. In the following subsections, under diff erent influence of magnetic field and boundary conditions, the numerical results using small $\mathrm{<mo\; stretchy="false">(</mo><mo><</mo>100<mo\; stretchy="false">)</mo>}$  , moderate $\mathrm{<mo\; stretchy="false">(</mo><mo>></mo>100}$  and $\mathrm{<mo><</mo>1000<mo\; stretchy="false">)</mo>}$  , high (10000), and even to ${10}^5$  Hartmann numbers are obtained.

Shercliff $s$   problem:  imposed transverse magnetic field outside a square channel with insulating walls
The MHD flow is in along duct channel with a cross section$\mathrm{<mo\; stretchy="false">[</mo>}-\mathrm{1,1<mo\; stretchy="false">]</mo>}\times \mathrm{<mo\; stretchy="false">[</mo>}-\mathrm{1,1<mo\; stretchy="false">]</mo>}$ , and the angle $\alpha$  between $\text{y}$ -axis and imposed magnetic field is $\pi \mathrm{<mo>/</mo>2}$ (Figure 3). Shercliff [5] proposed this example where the walls ofthe channel are insulated $\mathrm{<mo\; stretchy="false">(</mo>}B\mathrm{<mo>=</mo>0<mo\; stretchy="false">)</mo>}$  . The velocity is zero on the solid walls $\mathrm{<mo\; stretchy="false">(</mo>}u\mathrm{<mo>=</mo>0<mo\; stretchy="false">)</mo>}$  . In Tables 1&2, the approximate solutions of iRBF collocation method and local iRBF-DQ method for velocity at Hartmann number 100 on $\text{x}\mathrm{<mo>=</mo>0,}\text{y=}-1\sim 1$  are compared with uniform triangular element FEM and analytical solution, respectively. There are 10201 global nodes and 400 supporting nodes for local iRBF-DQ method, and 20000 elements for FEM. In Tables 1 and 2, we use 1681 and 2601 nodes for iRBF collocation method independently. The definition of relative error (i.e. R.E.) is $E_{\text{R}}\mathrm{<mo>=</mo><mo>|</mo>}\frac{e_{\text{experimenta}1}-e_{\text{exact}}}{e_{\text{exact}}}\mathrm{<mo>|</mo>}$ . We can find that the accuracy close to the walls becomes worse. In Figure 4, the relative errors at $\text{M}\mathrm{<mo>=</mo>100}$  of uniform FEM can only reach to at least about 10th power minus 1, but iRBF collocation method and local iRBF-DQ method can up to at least about 10th power minus 2 and 10th power minus 3, respectively. On the other hand, the relative errors at $\text{M}\mathrm{<mo>=</mo>500}$  of uniform FEM and iRBF collocation method can only reach to at least about 10th power minus 1, but local iRBF-DQ method can get up to at least about 10th power minus 2. In these tables, the exact solutions of velocity and magnetic induction are presented in,⁵.

In order to apply the approximation data to real blanket channels (high magneticfield), the high Hartmann numbers ${\mathrm{<mo\; stretchy="false">(</mo>10}}^3\le M\le {10}^5\mathrm{<mo\; stretchy="false">)</mo>}$  MHD problems must be solved. In Figures 5&6, the approximate solutions of velocity and magnetic induction using iRBF collocation method are compared with the approximate solutions of local iRBF-DQ method for Hartmann number $M{\mathrm{<mo>=</mo>10}}^3$  and ${10}^5$ . In Figure 7, the approximate solutions of velocity and magnetic induction using the local iRBF-DQ method at $M{\mathrm{<mo>=</mo>10}}^5$  is also presented independently. There are some oscillations near the walls of contour lines by iRBF collocation method as shown in Figures 5&6. However, the contour lines by local iRBF-DQ method are much smoother than by iRBF collocation method. The iRBF collocation method cannot obtain good numerical results of$B$ and $u$  by the reason of ill-conditional matrix. However, local iRBF-DQ method can still obtain accurate numerical results of velocity and induced magnetic field, especially when Hartmann number is up to105. It can be seen that the boundary layers become thinner as Hartmann number increases (Figures 5-8). We also observe that the Hartmann number and the velocity and magnetic induction are varied inversely. It means that the values of velocity and magnetic induction approximate to 1/10 of their original values when Hartmann number becomes 10 times larger (see more detail in Appendix A).

Drectangular duct flow with insulating walls, under the influence of an oblique magnetic field

In this subsection, the numerical experiments are the same as the one in section 4.1 except the angle $\alpha$ between $\text{y}$ -axis and imposed magnetic field changes fr om $\alpha \mathrm{<mo>=</mo>}\pi \mathrm{<mo>/</mo>2}$  to <$\alpha \mathrm{<mo>=</mo>}\pi \mathrm{<mo>/</mo>3}$ and $\pi \mathrm{<mo>/</mo>4}$ . In order to understand how diff erent angles between imposed magnetic field and $\text{y}$ -axis will aff ect the velocity and magnetic induction, the approximate solutions of $u$  and $B$  by iRBF collocation method are compared with the approximate solutions of$u$  and$B$  by local iRBF-DQ method for $\alpha \mathrm{<mo>=</mo>}\pi \mathrm{<mo>/</mo>3}$  and$\pi \mathrm{<mo>/</mo>4}$ in section 4.2. In Figures 9&10, we approximate a medium high Hartmann number $M{\mathrm{<mo>=</mo>10}}^3$  and${10}^5$ , and the patterns show that the contour lines by iRBF collocation method are not smooth enough. There are some oscillations near the walls of velocity contour lines. We also show the capability to approximate the velocity and induced magnetic field by local iRBF-DQ method at $M{\mathrm{<mo>=</mo>10}}^5$  in Figure 11. The results show that as Hartmann number increases the boundary layers formation close to the walls for both velocity and magnetic induction is well observed. Velocity becomes stagnant at the center of the channel. We can also find an inverse ratio relationship between Hartmann numbers and the values of velocity and magnetic induction. It means that the values of velocity and magnetic induction approximate to $\frac{1}{10}$ when Hartmann number becomes 10 times larger. local nodes  In Figsures 12&13, we approximate a medium high Hartmann number $\text{M}{\mathrm{<mo>=</mo>10}}^3$  and <${10}^4$  and the patterns show that the contour lines by iRBF collocation method are more smoother than the contour lines for$\alpha \mathrm{<mo>=</mo>}\pi \mathrm{<mo>/</mo>3}$ . This is because that the symmetrical feature of $\alpha \mathrm{<mo>=</mo>}\pi \mathrm{<mo>/</mo>4}$  causes the velocity contour lines parallel to$x\mathrm{<mo>=</mo>}y$ . There are still some oscillations near the walls of velocity contour lines. We also show the capability to approximate the velocity and induced magnetic field by local iRBF-DQ method at $M{\mathrm{<mo>=</mo>10}}^5$  in Figure 14. The results show that as Hartmann number increases the boundary layers formation close to the walls for both velocity and magnetic induction is very obviously observed. The velocity becomes a little stagnant at the center of the channel. When external magnetic field applies obliquely, the boundary layers are concentrated near the corners in the direction of the field for both solutions of velocity and induced magnetic field. These are the well-known behaviors of the magnetohydrodynamic flows (see more details in Appendix A). We can also find an inverse ratio relationship between Hartmann number and the values of velocity and magnetic induction. It means that the values of velocity and magnetic induction approximate to $\frac{1}{10}$  when Hartmann number becomes 10 times larger. In case 4.1 and case 4.2, we show the capability of our local iRBF-DQ method scheme to establish numerical approximations to the solution of MHD duct flows for different angles between imposed magnetic field and $\text{y}$ -axis even at very high Hartmann numbers

Imposed transverse magnetic field outside a square channel with partly insulating walls andpartly conducting walls

In Figure 16&17, we compare the approximate solutions of iRBF collocation method with the approximate solutions of local iRBF-DQ method at high Hartmann numbers $\text{M}{\mathrm{<mo>=</mo>10}}^3$  and ${10}^4$ . The patterns show that the contour lines by iRBF collocation method are almost the same with the contour lines by local iRBF-DQ method. However the approximate solutions near the walls of iRBF collocation method are worse than the approximate solution of local iRBF-DQ method. The contour lines near the edge by iRBF collocation method are not smooth enough. These figures show that there are more variations near the walls of velocity contour lines by iRBF collocation method as Hartmann number increases. However, the results for both $u$  and $B$ by local iRBF-DQ method are still accurate. We can find that the boundary layers become thinner as Hartmann number increases. In this case, the boundary layers are not only near to the walls but also inside the domain. Parabolic boundary layers are formed in the direction of imposed transverse magnetic field close to the discontinuity points $y\mathrm{<mo>=</mo>}-l$ and $y\mathrm{<mo>=</mo>}+l\mathrm{,}$ and thus cause three stagnant regions for the velocity (Figures 16&17). We can also find the dominant characteristics ofMHD flows at very high Hartmann numbers. Flow in the channel is totally separately into two rectangular flow regions with the boundary layers leaving the heart regions stagnant. Therefore our local iRBF-DQ scheme has the capability to deal with this kind boundary condition of MHD problems for high Hartmann numbers.

2-D circular duct flow with insulating walls, under the influence of an imposed transverse magnetic field

Figure 3Shercliff’s problem.

Figure 4 The comparison of relative errors of velocity.

Figure 7 The velocity and magnetic induction oflocal iRBF-DQ method at $\text{M}{\mathrm{<mo>=</mo>10}}^5$ ge002_0009.gif" />
local iRBF-DQ method (401 401); local nodes(9 9) .

Figure 8 A rectangular flow with an oblique magnetic field.

Figure 9 The velocity and magnetic induction at $\text{M}\mathrm{<mo>=</mo>1000<mo\; stretchy="false">(</mo>}\alpha \mathrm{<mo>=</mo>}\chi \mathrm{<mo>/</mo>3<mo\; stretchy="false">)</mo>}$  : (a) iRBF collocation
method (61 61) and (b) local iRBF-DQ method (101 101); local nodes (20 20) 0.8.

Figure 10 The velocity and magnetic induction at $\text{M}\mathrm{<mo>=</mo>10000<mo\; stretchy="false">(</mo>}\alpha \mathrm{<mo>=</mo>}\chi \mathrm{<mo>/</mo>3<mo\; stretchy="false">)</mo>}$  : (a) iRBF collocation method (81 81) and (b) local iRBF-DQ method (191 191); local nodes${\mathrm{<mo\; stretchy="false">(</mo>9}}^{\mathrm{<mo>*</mo>}}\mathrm{9<mo\; stretchy="false">)</mo>}$

Figure 11 The velocity and magnetic induction oflocal iRBF-DQ method at, local iRBF-DQ method (401 401); local nodes .

Figure 12 The velocity and magnetic induction at $\text{M}\mathrm{<mo>=</mo>1000<mo\; stretchy="false">(</mo>}\alpha \mathrm{<mo>=</mo>}\chi \mathrm{<mo>/</mo>4<mo\; stretchy="false">)</mo>}$  :(a)iRBF collocation method (61 61) and (b) local iRBF-DQ method (101 101); local nodes (2020).

Figure 13 The velocity and magnetic induction at $\text{M}\mathrm{<mo>=</mo>10000<mo\; stretchy="false">(</mo>}\alpha \mathrm{<mo>=</mo>}\chi \mathrm{<mo>/</mo>4<mo\; stretchy="false">)</mo>}$ :(a)iRBF collocation method (81 81) and (b) local iRBF-DQ method (184 184); local nodes27

Figure 14 The velocity and magnetic induction oflocal iRBF-DQ method at , local iRBF-DQ method (391 391);

Figure 15 A rectangular flow with partly conducting walls.

Figure 16 The velocity and magnetic induction at $\text{M}\mathrm{<mo>=</mo>1000}$ :(a)iRBF collocation method(60 60) and (b) local iRBF-DQ method (101 101); local nodes (20 20).

Figure 17 The velocity and magnetic induction at $\text{M}\mathrm{<mo>=</mo>10000}$ :(a)iRBF collocation method(90 90) and (b) local iRBF-DQ method (299 299); local nodes .

Figure 18 A circular flow with a transverse magnetic field.

Figure 19 The velocity and magnetic induction oflocal iRBF-DQ method at $\text{M}\mathrm{<mo>=</mo>100}$ local iRBF-DQ method (71082) ; local nodes(138).

Figure 20The velocity and magnetic induction oflocal iRBF-DQ method at $\text{M}\mathrm{<mo>=</mo>500}$ local iRBF-DQ method (71082) ; local nodes(138).

Figure 21The velocity and magnetic induction oflocal iRBF-DQ method at $\text{M}\mathrm{<mo>=</mo>1000}$ local iRBF-DQ method (102263) ; local nodes (200).

Figure 22The velocity and magnetic induction oflocal iRBF-DQ method at $\text{M}\mathrm{<mo>=</mo>10000}$ local iRBF-DQ method (102263); local nodes(200).

Figure 23Formation of MHD Lorentz force.

Figure 24The velocity contour lines when x = 0, y = -1 ~ 1 at M = 100 ~ 105.

Figure 25 The diagram of Hartmann effect and Hartmann layer.

In this study, we have considered the application of local iRBF-DQ method to the MHD flow problems in a straight channel of uniform or arbitrary cross-sections. The main work of this investigation is an accurate approximation and easy implementation of the two-dimensional MHD flow problems. The results indicate that the local iRBF-DQ method is stable at moderate and high values of Hartmann number as comparing to the iRBF collocation method and FEM method. The localization feature reduces the computational time and solves the ill-condition problem of conventional RBFs methods. As the results shown we can conclude the present meshless local iRBF-DQ method can stably, accurately and quickly obtain the weighting coefficients of derivative approximation with arbitrary cross-sections even at very high Hartmann numbers in MHD flow problems.

None.

The authors declare that there is no conflict of interest.

1. Molokov S, Reed CB. Liquid metal magnetohydrodynamic flows in circular ducts at intermediate hartmann numbers and interaction parameters. Magnetohydrodynamics. 2003;39(4):539–546.
2. Buhler L. Liquid metal magnetohydrodynamics for fusion blankets. In: Molokov S, Moreau R, editors. Magnetohydrodynamics, Historical Evolution and Trends. 2007.
3. Branover H. Manetohydrodynamic Flow in Ducts. New York: Wiley and Sons; 1984.
4. Garner RA, Lykoudis PS. Magneto-fluid-mechanic pipe flow in a transverse magnetic field, Part 2: Heat transfer. J Fluid Mech. 1997;48(1):29–141.
5. Dragos L. Magneto-Fluid dynamics. Abacus Press, UK, 1975.
6. Shercliff JA. Steady motion of conducting fluids in pipes under transverse magnetic fields. Proc Cambridge Philos Soc. 1953;49:136–144.
7. Singh B, Lal J. MHD axial flow in a triangular pipe under transverse magnetic field. Indian J Pure Appl Math 1978;9:01–115.
8. Singh B, Lal J. MHD axial flow in a triangular pipe under transverse magnetic field parallel to a side of the triangle. IndJ Techno.1979;117(5):184–189.
9. Singh B, Lal J. FEM in MHD channel flow problems. Int J Numer Methods Eng. 1982;18(7):1104–1111.
10. Singh B, Lal J. FEM for MHD channel flow with arbitrary wall conductivity. J Math Phys Sci. 1984;18:501–516.
11. Singh B, Lal J. FEM for unsteady MHD flow through pipes with arbitrary wall conductivity. Int J Numer Methods Fluids. 1984;4(3):291–302.
12. Tezer Sezgin CM. Koksal S, Finite element method for solving MHD flow in a rectangular duct. Int J Numer Methods Eng. 1989;28(2):445–459.
13. Demendy Z, Nagy T. A new algorithm for solution of equations of MHD channel flows at moderate Hartmann numbers. Acta Mech 1997;123(4):135–149.
14. Tezer–Sezgin CM, Han Aydn S, Solution of magnetohydrodynamic flow problems using the boundary element method. Eng Anal Bound Elem 2006;30(5): 411–418.
15. Bozkaya C. Tezer Sezgin CM, Fundamental solution for coupled magnetohydrodynamic flow equations. J Comput Appl Math. 2007;203(1):125–144.
16. Nesliturk AI, Tezer Sezgin CM, The finite element method for MHD flow at high Hartmann numbers. Comput Methods Appl Mech Eng. 2005;194(9–11):1201–1224.
17. Kansa EJ. Multiquadrics: a scattered data approximation scheme with applications to computational fluid-dynamics. Comp Math Appl. 1990;19(8–9):147–161.
18. Franke C, Schaback R, Solving partial diff erential equations by collocation using radial basis functions. Appl Math Comput 1998;93(1):73–82.
19. Franke C, Schaback R. Convergence order estimates of mesh less collocation methods using radial basis functions. Adv Comput Math. 1998;8(4):381–399.
20. Franke C, Scattered data interpolation: test of some methods. Math Comput. 1982;38:181–200.
21. Carlson RE, Foley T. The parameter R in multiquadric interpolation. Comp Math App. 1991;21(9):29–42.
22. Duchon J. Splines minimizing rotation invariant semi–norms in Sobolev spaces. Constructive Theory of Functions of Several Variables. 1977;38:85–100.
23. Hardy RL. Theory and application of the multiquadric–biharmonic method. Comp Math Appl 1990;19(8–9):163–208.
24. Wu YL, Shu C. Development of RBF-DQ method for derivative approximation and its application to simulate natural convection in concentric annuli. Comput Mech. 2002;29(6):477–485.
25. Kansa EJ, Hon YC. Circumventing the ill-conditioning problem with multiquadric radial basis functions: applications to elliptic partial diff erential equations. Comput Math App. 2000;39(7–8):123–137.
26. Shu C, Yeo KS. Local radial basis diff erential quadrature method and its application to solve two–dimensional incompressible Navier–Strokes equations. Comput Meth Appl Mech Eng. 2003;192(7–8):941–954.  
27. Mai Duy N, Tran Con T. Approximation of function and its derivatives using radial basis function networks. Appl Math Modelling. 2003;27(3):97–220.
28. Mai Duy N, Tran Con T. Numerical solution of diff erential equations using multiquadric radial basis function networks. Neural Networks 2001;14(2):185–199.
29. Ngo Cong D, Mai–Duy N, Karunasena K, et al. A numerical procedure based on –IRBFN and local MLS–ID–IRBFN methods for fluid–structure interaction analysis. CMES: Comput Modeling in Eng Sci. 2012;83(5):459–498.
30. Thai Quang N, Le Cao K, Mai Duy N, et al. A numerical scheme based on compact integrated–RBFs and Adams Bash forth/Crank Nicolson algorithms for diffusion and unsteady fluid flow problems. Eng Anal Bound Elem. 2013;37(12):1657–1667.
31. Ngo Cong D, Tien CM, Nguyen Ky T, et al. A generalised finite difference scheme based on compact integrated radial basis function for flow in heterogeneous soils. Int J Numer Meth Fluid. 2017;85(7):404–429.
32. Shu C, Wu YL. Integrated radial basis functions–based diff erential quadrature method and its performance. Int J Numer Meth Fluid. 2007;53(6):969–84.
33. Barrett KE. Duct flow with a transverse magnetic field at high Hartmann numbers. Int J Numer Meth Eng. 2001;50(8):1893–1906.