Chebyshev Spectral Collocation Method Approximation to Thermally Coupled MHD Equations

In this study, a Chebyshev spectral collocation method (CSCM) approximation is proposed for solving the full magnetohydrodynamics (MHD) equations coupled with energy equation. The MHD flow is two-dimensional, unsteady, laminar and incompressible, and the heat transfer is considered using the Boussinesq approximation for thermal coupling. The flow takes place in a square cavity which is subjected to a vertically applied external magnetic field, and the presence of the induced magnetic field is also taken into account due to the electrical conductivity of the fluid. The governing equations given in terms of stream function, vorticity, temperature, magnetic stream function, and current density, are solved iteratively using CSCM for the spatial discretisation, and an unconditionally stable backward difference scheme for the time integration. The induced magnetic field is obtained by means of its relation to the magnetic stream function. The behaviours of the flow and the heat transfer are investigated for varying values of Reynolds ($Re$), magnetic Reynolds ($Rem$), Rayleigh ($Ra$) and Hartmann ($Ha$) numbers.


Introduction
Magnetohydrodynamic (MHD) flow and heat transfer of electrically conducting fluids in enclosures has been the subject of a great number of investigations due to the wide range of significant applications in various fields such as nuclear reactors, MHD generators, and metallurgical industries.The physical models that describe MHD flow and heat transfer consist of coupling the hydrodynamic problem by means of the convective term in the magnetic field equation, coupling the magnetic field problem in Navier-Stokes equations with Lorentz's force, and the thermal coupling through Boussinesq's assumption.Extensive research is ongoing in the direction of designing numerical techniques applicable to the MHD flow and heat transfer models, the reason being that analytical solutions apply only under very restrictive conditions.In majority of the studies, the aim is to determine the influence of the characteristic problem parameters such as Reynolds, Rayleigh, and Hartmann numbers on the flow and heat transfer.In particular, the influence of an externally applied magnetic field on the flow and heat transfer is of primary importance.Many different numerical methods have been introduced for approximating solutions to MHD flow and heat transfer problems.
Ece and Büyük [1], presented differential quadrature solutions of natural convection flow under a magnetic field in an inclined rectangular enclosure heated or cooled on adjacent walls, and with isothermal vertical or adiabatic horizontal walls.The radial basis function formulation has been used by Colaço et al., in [2], to approximate a steady MHD problem with heat transfer in a square cavity whose horizontal walls are kept at different and constant temperatures.In addition, a meshless method based on local radial basis function collocation method and fractional step method for pressure-velocity coupling, is proposed by Mramor et al. for solving MHD convection in a square cavity [3].Oztop et al. [4], have analysed a mixed convection flow in the presence of a magnetic field in a top sided lid driven cavity heated by a corner heater using finite volume method in combination with SIMPLE algorithm.Another finite volume application is the work of Al-Salem et al. [5], where the MHD mixed convection in a non-isothermally heated square enclosure is considered.The finite element method has been implemented for the approximation to steady, laminar, natural convection flow in inclined enclosures in the presence of an oblique magnetic field in the work of Türk and Tezer-Sezgin [6].
The majority of the above listed studies have assumed a very low magnetic Reynolds number that represents the ratio of advection to diffusion in the magnetic field, and have simplified the equation of magnetic induction.This assumption reduces the number of equations in the system to be solved and hence significantly lowers the computational cost.On the other hand, Sarris et al. [7] have examined the MHD natural convection model and investigated the limitations in predicting the flow and heat transfer characteristics by the hypothesis that Lorentz's force is suppressed into a damping term resisting the motion of the fluid.Their results set forth that for large values of Hartmann number, the magnetic induction equations should be taken into consideration in the mathematical model.There are several studies, implementing different numerical methods, where the full model for liquid metal flows is considered, and the existence of external and internal magnetic fields is taken into account.Şentürk et al. proposed a two-dimensional scheme in [8] to simulate the natural convection with internal heat generation and absorption in addition to non-linear time dependent evolution of heated and magnetized liquid metals exposed to external fields.The algorithm consists of a Lax-Wendroftype matrix distribution together with a dual time-stepping scheme employing third-order Runge-Kutta steps.The solution of the incompressible MHD flow equations using a dual reciprocity boundary element method formulation is considered by Bozkaya and Tezer-Sezgin in [9].The full MHD model applicable to three-dimensional problems with thermal coupling have been considered numerically by a stabilized finite element method in the work of Codina and Hernández [10].Pekmen and Tezer-Sezgin, have proposed a dual reciprocity boundary element solution for full MHD equations in a lid driven square cavity in [11].The full MHD equations in an alternative formulation has been considered by Sivakumar et al. [12] to study the influence of an induced magnetic field on the forced convective heat transfer from an isothermal sphere in the presence of an external magnetic field.The equations are expressed in stream function, vorticity and magnetic stream function, and are solved using finite difference methods combined with multi-grid techniques.In a recent study [13], Selimli and Recebli investigated the cooling period of a liquid metal flowing under the effect of an externally applied magnetic field with the use of a commercial software, AN-SYS.The Chebyshev spectral collocation method (CSCM) is a well developed and widely used technique based on expanding the unknown fields in tensor product of Chebyshev polynomials.The polynomials are differentiated analytically on the abscissae of the extreme points of the Chebyshev polynomials with the construction of differentiation matrices, where higher order derivatives can easily be obtained by multiplying these matrices.Therefore, a high order accuracy is achieved in combination with computational efficiency.The method is also being successfully implemented in the models of computational fluid dynamics.An application of CSCM enriched with a multi-domain technique to the approximation of unsteady natural convection heat transfer (without the magnetic field effect) in an enclosure with a square body has been presented in [14].The objective of the current study is to present a numerical approach based on CSCM for the approximation of thermally coupled MHD problems.The unsteady regularised driven cavity flow is considered where the fluid is incompressible and assumed to be electrically conducting.The external magnetic field is applied in the vertical direction.The governing equations are given in stream function, vorticity, temperature, magnetic stream function, and current density.A novel iterative scheme is designed based on CSCM for the spatial discretisation in combination with an implicit backward finite difference for the time integration.The mechanism additionally provides a means of approximating the unknown vorticity and current density boundary conditions from the velocity and magnetic field, respectively, with the use of differentiation matrices.In essence, this study extends the CSCM applicability to full MHD equations with heat transfer and magnetic induction taking place on a regularised lid driven cavity, for the first time, to the best of the author's knowledge.The effects of the characteristic flow parameters on the MHD flow and heat transfer are investigated, and the results obtained from the simulations are presented in terms of the contours of the unknowns.

Mathematical Model
The two-dimensional, transient, laminar, incompressible MHD flow and heat transfer taking place in a square cavity is considered.The flow is subjected to an externally applied magnetic field of a constant intensity B 0 in the positive vertical direction.The fluid is taken to be electrically conducting and hence, the induced magnetic field is also included in the model.It is assumed that the Joule heating and Hall effects, as well as the influence of the viscous dissipation, displacement and convection currents are negligible.Under these assumptions, the MHD equations can be obtained by coupling the Navier-Stokes equations with Maxwell's equations by means of Ohm's law.If the thermal coupling is modelled through Boussinesq's assumption, the governing equations are given as [15,16] ∂ ū (1) In these equations, the unknowns are; ū and v, the respective x− and ȳ− components of the velocity of the fluid, p, the pressure, Bx and By , the induced magnetic field components, and T , the temperature of the fluid.T c represents the reference temperature (of the cold boundary), ν is the kinematic viscosity, ρ is the reference density of the fluid, µ m is the magnetic permeability, β is the thermal expansion coefficient, σ is the electrical conductivity, and α is the thermal diffusivity of the fluid.The direct solution strategies applied to system (1) have traditional difficulties due to existence of the pressure terms.Therefore, alternative formulations have been introduced to represent these equations and have been very successfully adopted by many researchers.For the two-dimensional flows, the most consistently used one is the well established stream function-vorticity formulation.The pressure terms in the momentum equations in system (1) can be eliminated by subtracting the derivative of the second equation of (1) with respect to ȳ from the derivative of the third equation of (1) with respect to x.If the stream function ψ is defined so that ū = ∂ ψ/∂ ȳ and v = −∂ ψ/∂ x, then the continuity equation is automatically satisfied.Now, the only non-zero component of the vorticity field is introduced as w = ∂ v/∂ x − ∂ ū/∂ ȳ, yielding the relation between ψ and w as w = −∆ ψ.The divergence-free condition for the magnetic field in (1), that is, ∂ Bx /∂ x + ∂ By /∂ ȳ = 0, could be considered as redundant.In numerical point of view, the system can be solved without the use of this equation.In the continuous level, this condition is implied in combination of Ampere's and Ohm's laws to reduce the Maxwell equations.In other words, the divergence-free condition of the magnetic field is automatically satisfied through the domain of the problem provided that the initial field is solenoidal.However, at the discrete level, this condition may not be satisfied, and it is well known that violation of it in the numeri-cal computations, results in a non-physical force (see, e.g.[16][17][18]).The previous argument sets another restriction on the applicability of direct solution strategies.There are several remedies successfully applied by many researchers (see [10,16,19] and the references therein).A natural approach in analogy to the stream function vorticity formulation, is to introduce the magnetic stream function Ā, relating to the magnetic field components as Bx = ∂ Ā/∂ ȳ and By = −∂ Ā/∂ x, and the current density j defined by j = ∂ By /∂ x − ∂ Bx /∂ ȳ (see e.g.[16] for the details).Further, the following dimensionless variables are defined where is a characteristic length, u 0 is a characteristic velocity of the fluid, δ T is the temperature difference between hot and cold walls.Consequently, the governing equations given in system (1), can be written in nondimensional form in terms of ψ, w, T , A, and j as In the equations above, dimensionless parameters, namely, Reynolds number (Re), Prandtl number (Pr), magnetic Reynolds number (Rem), Rayleigh number (Ra), and Hartmann number (Ha), have been introduced as follows where µ is the dynamic viscosity.
Homogeneous initial conditions all holding in the spatial domain at time t = 0 accompany the equations in the given model.The boundary conditions are imposed as follows.
The velocity at the top wall of the cavity is given by ( û, 0) for a prescribed û, whereas the (0,0) condition is imposed on the other walls.Accordingly, the value of the stream function on boundary is known, as the walls are natural streamlines, hence the condition ψ = 0 is considered.In the given configuration, the induced magnetic field boundary conditions are B x = 0 and B y = 1, accordingly the condition A = −x is imposed on all of the walls.The vertical walls are considered to be adiabatic, while the horizontal walls are kept at constant temperatures which are given as T h = 0.5 at the upper wall and T c = −0.5 at the lower wall (see Figure 1).The vorticity and current density boundary conditions are not known, and they are calculated numerically with the use of stream function values for the former, and magnetic stream function values for the latter as described in the next section.
The problem configuration and the boundary conditions.

CSCM Application and the Time Integration
System (2) consists of coupled equations which are complicated not only by the non-linearity but also by the lack of the vorticity and current density boundary conditions.Therefore, it is compulsory to devise a solution methodology that is capable of handling these complications.The approach implemented in the present study is an iterative procedure based on space and time discretisation of the equations in the given system as described in the sequel.The spatial discretisation of the equations in ( 2) is based on requiring the numerical approximation of each unknown to be exactly satisfied on the abscissae of the extreme points of the Chebyshev polynomials.Therefore, the technique is referred as Chebyshev spectral collocation method.In this method, each function spans the whole domain under consideration and thus, the derivatives of the function depend on the entire discretisation.A function Φ(x) defined in [−1, 1] is interpolated by the polynomial Φ(x) of degree at most N, of the form [20][21][22] with Φ(x j ) = Φ(x j ), and C j (x) is a Cardinal function (or Lagrange basis) of degree N defined using the Chebyshev polynomials of the first kind (T n (x) = cos(n arccos x), n = 0, 1, . . ., N) by where c 0 = c N = 2, and c j = 1, for 1, . . ., N − 1.
They possess the desired property of being clustered through the end points of the interval, consequently in a two-dimensional domain, having a concentration of grid lines near the boundaries (see Figure 2).The n−th derivative of Φ(x) is then approximated by The first derivative at these nodes satisfy C (1) j (x i ) = d i j where (−1) i+ j x i − x j , i = j, i, j = 0, . . ., N, Now, the discrete values of the first derivative of the function Φ can be obtained as This equation can be written in a matrix-vector form as N Φ where D N = [d i j ] is called the first order Chebyshev spectral differentiation matrix which is of size (N + 1) × (N + 1).In order to minimise the round-off errors for the calculation of the first derivatives, the diagonal entries d ii are computed as [22] The approximation to the n−th order derivative of the function Φ(x) now becomes where D (n) = [D (1) ] n , that is, the n-times matrix multiplication of D (1) .In general, D (n) is referred as the n-th order Chebyshev spectral differentiation matrix.The use of matrix multiplication for higher order derivatives and the use of Equation ( 8) for obtaining the diagonal entries, lead to a significantly greater accuracy in the computation of second and higher order derivatives for a wide range of functions.The Chebyshev spectral differentiation matrix for functions defined on an arbitrary interval In an analogous way, the first order differentiation matrix in the ydirection, denoted by i j ], is computed.The second order Chebyshev differentiation matrix in this direction, E (2) , is calculated by simply squaring E (1) .For the time integration, the unconditionally stable backward difference scheme is implemented which is defined by ∂ φ /∂t s+1 = (φ s+1 − φ s )/δt, s and δt being the time level and the time step, respectively.Having constructed all the necessary ingredients, one can obtain the fully space-and-time discretised equations by substituting the respective approximations ψ, w, T , A and j, and to ψ, w, T , A, and j, into the equations of system (2).The approximation vectors are of order (N + 1) 2 and computed in the following pattern φ = [φ (x 0 , y 0 ), . . ., φ (x N , y 0 ), . . ., φ (x 0 , y N ), . . ., φ (x N , y N )] T .
Let D( φ ) denote the diagonal matrix with the entries of a vector φ on its diagonal, and introduce the Kronecker product of an m 1 × n 1 matrix P and an m 2 × n 2 matrix Q that is defined as the m 1 m 2 × n 1 n 2 matrix P ⊗ Q, which is the m 1 × n 1 block matrix whose i j-th block is the m 2 × n 2 matrix p i j Q.Then, the CSCM-and-time discretised form of the equations in (2) can be written as In these equations, the (N + 1) 2 × (N + 1) 2 matrix K is defined as where D(i) = I ⊗ D (i) and Ê(i) = E (i) ⊗ I, D (i) and E (i) , i = 1, 2, being the order, are the Chebyshev differentiation matrices in x− and y− directions, respectively.I is the identity matrix of order (N + 1) 2 .Moreover, P( φ , φ) denotes the vector formed by multiplication of the approximations to the first partial derivatives of its argument vectors, and is defined as In (10), ι is the vector of order (N + 1) 2 whose all entries are 1.

The iterative solution procedure
As already noted, there are critical difficulties in characterising a flow using the stream function vorticity formulation.One of the causes of these difficulties is the existence of an infinite sequence of eddies in the corners where no-slip boundary conditions are imposed.In such instances, the second order derivatives of the vorticity are unbounded in the vicinity of the corners.The derivation of boundary conditions for the vorticity is another cause when the specified boundary conditions are discontinuous at the corners (see, e.g., [23]).These irregularities significantly affect the high accuracy of spectral methods as the contamination created by the corner singularities is extended to the whole solution due to the global nature possessed by the methods.The primary interest of the present study is the investigation of characteristic flow parameters on the MHD flow and heat transfer with the use of CSCM, therefore, the regularised driven cavity flow is considered.This is a setting widely used in incompressible flow models where the boundary conditions are smoothed in such a way that the velocity and its derivatives vanish at the corners (see, [24], [25], [26], and the references therein).
Remark 2.1.The cases of more practical interests on account of fluid flow mechanism such as the singular lid driven cavity flow, could be handled by several methodologies alternatively, (see [23] and [27] for a detailed description).
The essential feature of the approach followed in this study is its iterative nature derived to handle the non-linearity of the equations, additionally providing a means of approximating the vorticity and current density boundary conditions based on differentiation matrices.In accordance with this, the equations in system (10) are solved in the written order starting with the given initial conditions.First, the stream function equation is solved.Next, the velocity components are updated by the relation and, the vorticity boundary values are calculated as where l denotes the l-th boundary node.Next, the equation for the magnetic stream function is solved.A similar approach as in the case of vorticity is followed for the computation of the current density boundary values.The magnetic field components are obtained at the time level s + 1, and, these value are used to calculate the current density boundary values as Having obtained its boundary values, the current density equation is solved.This is followed by the solution of the energy equation.The vorticity equation is solved at the last step of each cycle with the imposition of the boundary values approximated by using the updated velocity approximations at the time level s + 1.This iterative procedure continues until a preassigned convergence tolerance between two successive iterations is reached for all the unknowns on the whole problem region.

Results
The numerical simulations for the MHD flow and heat transfer in a square cavity under the influence of externally applied magnetic field are presented in this section.The solution procedure described in Section 2.2 has been implemented.The Prandtl number is taken as 0.1 throughout the study, and the influence of different combinations of Re, Ra, Ha, and Rem on the flow and heat transfer is investigated.The time integration has been carried out by using a constant time step δt = 0.25.The stopping criteria of the iterative scheme has been set to be 10 −5 for all the unknowns, and the solutions in regard to this criteria are referred as the steady-state solutions.The present work is mainly concerned with the assessment of the effectiveness of CSCM applied to MHD flow and heat transfer.Thus, the regularised driven cavity flow in which the horizontal component of the velocity at the upper wall is specified by û = 4x 2 (1 − x 2 ), has been considered in all the numerical simulations presented in the sequel (see Section 2.2).First, the mesh dependency of the solution is investigated in terms of the maximum absolute values of the stream function, vorticity, and the current density when Ra = 10 5 , Ha = 100, Re = 1000, and Rem = 100.The results are illustrated in Figure 3.It is observed that the numerical solutions do not show a significant difference for N values that are larger than 50, and hence, the case N = 50 is considered in all the simulations.Next, the regularised lid-driven cavity flow studied in [9,16] is considered in order to validate the computational mechanism designed in the present study.In this setting, the flow takes place in a square duct subjected to a transverse magnetic field in the absence of heat sources or sinks.The magnetic field distribution for Rem = 0.1 and Rem = 100, where Ha = 10 and Re = 100, are depicted in Figure 4.The streamlines of the magnetic field can be compared to the ones given in [9] and [16] revealing a good agreement.Further, to allow a quantitative comparison, the results in terms of the primary vortex and the magnetic field intensity for the case of higher magnetic Reynolds number, are tabulated in Table 1.These results agree reasonably well with the ones given in [16].
The results presented below, concern cases in which the iterative solution procedure applied to thermally coupled MHD equations given in (10).The numerical results in terms of the contours of the unknowns, namely, stream function, vorticity, temperature, magnetic stream function, and current density, are presented.The time evolution of the flow and heat transfer is depicted in Figure 5, for the values of Ra = 10 4 , Re = 400, Ha = 50, and Rem = 100.It is observed that the vortex in the flow formed along the upper wall due to the moving boundary, is transferred through the central region of the cavity with the time advance.The vorticity contours evolve accordingly with an increase in magnitude, and a formation of boundary layers in the upper-right corner is observed.The linear initial profile in the isotherms is distorted through the steady-state at which there exists a small curvature indicating a slight influence of convection due to a moderate value of Ra.The initial linear distribution of the magnetic streamlines following the external magnetic field is reformed as time advances, and a circulation near the right wall of the cavity is observed due to the convection effect.Following the flow behaviour, the current density contours form a vortex in the upper corner which is spreading through the central region as time advances, with boundary layers near the upper wall.The flow and heat transfer as well as the magnetic field and current density distributions inside the enclosure settle down at approximately t = 2136, where the steady-state is reached according to the given time-convergence tolerance.
In the rest of the study, the steady-state results are presented for several combinations of characteristic parameters, and hence, the time asymptotic behaviours of the solutions are investigated.of the similarity in their pattern, the magnitude of the current density differs significantly, rising with an increase in Rem.The simulation is repeated for Ha = 100, and the corresponding results are plotted in Figure 8.A comparison with Figure 7 shows that the effect of increasing Ha is prominent when Rem = 100 for the corresponding flow parameters.The retarding effect of the externally applied magnetic field is responsible for the circulation of the magnetic streamlines near the upper-right corner of the cavity.
In the same way as the lower Ha case, the distribution of the current density has a similar profile with increasing magnitude as Rem increases.
Figure 9 depicts the variation in the streamlines with respect to the Rayleigh and Hartmann numbers, where Re = 1000 and Rem = 1.As Rayleigh number increases, the flow is separated following the formation of a secondary eddy in the opposite direction.This phenomenon occurs at Ra = 10 5 when Ha = 25, however, it is observed for the higher value of Ra when Ha = 100, since higher intensity of the magnetic field reduces the effect of the Rayleigh number.As Hartmann number is increased from 25 to 100, the flow intensity caused by the motion of the upper lid moves towards the centre of the cavity.
The effect of the Rayleigh and Hartmann numbers on the heat transfer inside the enclosure when Re = 1000 and Rem = 1, is demonstrated by considering the isotherms presented in Figure 10.As Ra increases, for a fixed of Hartmann number, the dominance of the buoyancy force results in a convective flow.The higher value of Ha increases the resistance to the fluid motion, resulting in an enhancement of conductive heat transfer.

Discussion and Conclusion
In this paper, a numerical approximation for solving the thermally coupled full MHD equations is presented.The formulation is based on the Chebyshev spectral collocation method for the spatial derivatives, and a backward difference scheme for temporal derivatives.The essential feature of the solution procedure is its iterative nature that allows handling the non-linearity and providing a means of approximating the unknown vorticity and current density boundary conditions.The effect of the variation in several characteristic flow parameters on the flow and heat transfer is studied.The analysis of the results obtained shows that the flow field is influenced appreciably by the variation of the characteristic parameters tested.As Reynolds number increases, the flow concentrates through the upper and right walls with an increase in the intensity of the boundary layers.The existence of strong gradients clustered through the top-left and bottom-right regions indicates the dominance of convection over conduction for larger Re values.
The magnetic field distribution in the cavity is not affected much by increasing the magnetic Reynolds number to moderate values.However, for large Rem, such as 100, a strong circulation takes place near the upper-right region of the cavity, due to the dominance of convection when the Hartmann number is fixed.The distribution of the current density has a similar form for all tested magnetic Reynolds number values, even though the magnitude increases in   The increase in Rayleigh number causes a separation in the flow and a formation of a secondary vortex in the opposite direction.This behaviour is more pronounced for lower Hartmann number values, as an increase in Hartman number that amounts to a higher intensity of the magnetic field, reduces this effect.The convection dominated flow due to the buoyancy force can be observed for high Rayleigh number values when Hartmann number is fixed.On the other hand, an increase in Hartmann number, causes a rise in the resistance to the fluid motion where the conductive heat transfer is enhanced.

Figure 2 .
Figure 2. A sample node distribution for N = 16.
[a, b] can be constructed by a linear transformation η = (b − a)/2x + (a + b)/2, which maps the standard interval [−1, 1] to any finite interval [a, b].A distribution of the so-called Chebyshev-Gauss-Lobatto points used as collocation points in the problem domain where N = 16 is illustrated in Figure 2.

Table 1 .
The characteristics of the primary vortex and the magnetic field intensity.) max(B x ) min(B y ) max(B y ) Rem.When Hartmann number is increased, the effects of the magnetic Reynolds number on the magnetic field and current density are reduced, especially for the case of the highest Rem value tested. x