International Journal of Analysis and Applications Volume 19, Number 3 (2021), 440-454 URL: https://doi.org/10.28924/2291-8639 DOI: 10.28924/2291-8639-19-2021-440 NUMERICAL STUDY OF RAYLEIGH-BENARD PROBLEM UNDER THE EFFECT OF MAGNETIC FIELD ABDELFATAH ABASHER1,∗, MOHAMMED GUBARA2 , SULIMAN SHEEN2, IBRAHIM BASHIR2 1Mathematics Department, Faculty of Science, Jazan University, Jazan, KSA 2Mathematics Department, Faculty of Mathematical Science and Statistic, Al-Neelain University, Khartoum, Sudan ∗Corresponding author: amoaf84@gmail.com Abstract. In this paper, a linear stability analysis is studied for Rayliegh-Benard problem with the effect of magnetic field, a perturbation equations is solved numerically by using spectral Chebyshev tau method, the boundaries are considered both are free, both are rigid, the lower is free and the upper is rigid, the results were illustrated graphically and compared with previous studies. 1. Introduction The aim of this work is to solve the perturbation equations that represent Rayleigh-Benard problem under the effect of magnetic field ( [1], Page (160, 161)) as follows (1.1) ∂u ∂t = −∇ ( δp ρ0 + µ H � h 4πρ0 ) + gαθk + v∇2u + µ 4πρ0 (H �∇) h, (1.2) ∂h ∂t = (H �∇) u + η∇2h, Received February 24th, 2021; accepted March 26th, 2021; published April 28th, 2021. 2010 Mathematics Subject Classification. 00A69. Key words and phrases. Rayleigh-Benard problem; linear stability analysis; Chebyshev Tau method. ©2021 Authors retain the copyrights of their papers, and all open access articles are distributed under the terms of the Creative Commons Attribution License. 440 https://doi.org/10.28924/2291-8639 https://doi.org/10.28924/2291-8639-19-2021-440 Int. J. Anal. Appl. 19 (3) (2021) 441 (1.3) ∇ � h = 0, (1.4) ∂θ ∂t = βw + κ∇2θ, (1.5) ∇ � u = 0, where, u, δp, θ and h are the perturbation in the velocity, pressure, temperature and magnetic field respec- tively, ρ0 is the density at mean temperature T0, v is the kinematic viscosity, κ is the thermal diffusivity, α the coefficient of volume expansion, µ is the magnetic permeability, η is the resistivity and β is the temperature gradient defined as (1.6) β = T0 −Td d , The problem is studied analytically by various authors ( [1]- [4], [12]- [14]). Also many numerical methods were used, one of these method is Galerkin method which discussed in [14]. In this paper we used the numerical method namely Spectral Chebyshev tau method to determine the conditions of instability for various cases of the boundary conditions. The paper outlined as follows. In section 1 we formulate the governing mathematical perturbation equations. In section 2 a linear stability analysis for the perturbation equations. In section 3 we describe the method of solution Spectral Chebyshev Tau method for the three cases of the boundary conditions. Finally, we present our numerical results in which are computed using MATLAB. 2. Linear Stability analysis By taking the curl operator of equation (1.1), we get (2.1) ∂ω ∂t = gα ( ∂θ ∂y i− ∂θ ∂x j ) + v∇2ω + µ 4πρ0 (H �∇) v, where ω = ∇× u is the vorticity vector and v = ∇× h is the current density induced by the perturbation. Taking the curl operator of (2.1) again, we get (2.2) ∂∇2u ∂t = −gα ( ∂2θ ∂z∂x i + ∂2θ ∂z∂y j − ( ∂2θ ∂x2 + ∂2θ ∂y2 ) k ) + v∇4u + µ 4πρ0 (H �∇)∇2h. Take the curl of equation (1.2), we get (2.3) ∂v ∂t = (H �∇) ω + η∇2v. Int. J. Anal. Appl. 19 (3) (2021) 442 Now by equating the z-component of equation (1.2), (2.1), (2.2) and (2.3) respectively, we get (2.4) ∂hz ∂t = η∇2hz + (H �∇) w, (2.5) ∂ζ ∂t = v∇2ζ + µ 4πρ0 (H �∇) ξ, (2.6) ∂∇2w ∂t = gα ( ∂2θ ∂x2 + ∂2θ ∂y2 ) + v∇4w + µ 4πρ0 (H �∇) hz, (2.7) ∂ξ ∂t = η∇2ξ + (H �∇) ζ, where w, ζ and ξ and hz are the z-components of the velocity, the vorticity, the current density and the magnetic field respectively. When the direction of the magnetic field coincides with the vertical direction H = (0, 0,H), the required perturbation equations become (2.8) ∂ ∂t ∇2w = gα ( ∂2θ ∂x2 + ∂2θ ∂y2 ) + v∇4w + µH 4πρ0 ∂ ∂z ∇2hz, (2.9) ∂hz ∂t = η∇2hz + H ∂w ∂z , (2.10) ∂ξ ∂t = η∇2ξ + H ∂ζ ∂z , (2.11) ∂ζ ∂t = v∇2ζ + µH 4πρ0 ∂ξ ∂z . In addition to equation (1.4). The boundary conditions of θ, w and ζ for the three cases free-free, rigid-rigid and rigid-free are (2.12)   θ = 0,w = 0 at z = 0 and z = d, dζ dz = 0, d 2w dz2 = 0 at free boundary, ζ = 0, dw dz = 0 at rigid boundary. Int. J. Anal. Appl. 19 (3) (2021) 443 3. Normal Mode Analysis Let w,θ,ζ,hz and ξ be defined as (3.1)   zθ ζ w ξ hz   =   Θ(z) Z(z) W(z) X(z) K(z)   exp (i(kxx + kyy) + γt), where k = √ k2x + k 2 y is the wave number and γ is a constant. Substitute (3.1) into the system (1.4), (2.8)-(2.11), the system become (3.2) γΘ = βW + κ ( d2 dz2 −k2 ) Θ, (3.3) γK = η ( d2 dz2 −k2 ) K + H dW dz , (3.4) γ ( d2 dz2 −k2 ) W = −gαk2Θ + v ( d2 dz2 −k2 )2 W + µH 4πρ0 d dz ( d2 dz2 −k2 ) K, (3.5) γX = η ( d2 dz2 −k2 ) X + H dZ dz , (3.6) γZ = v ( d2 dz2 −k2 ) Z + µH 4πρ0 dX dz . And the boundary conditions (2.12) become (3.7)   Θ = 0,W = 0 at z = 0 and z = d, dZ dz = 0, d 2W dz2 = 0 at free boundary, Z = 0, dW dz = 0 at rigid boundary. Define the following non-dimensional variables, (3.8) a = kd, σ = γd2 v , z∗ = z d , P1 = v κ , and P2 = v η , the operators d dz = 1 d d dz∗ and d 2 dz2 = 1 d2 d2 dz∗2 , assume D = d dz∗ , then by substituting (3.8) into the system (3.2)-(3.6), we get (3.9) ( D2 −a2 −P1σ ) Θ = − ( βd2 κ ) W, Int. J. Anal. Appl. 19 (3) (2021) 444 (3.10) ( D2 −a2 −P2σ ) K = − ( Hd η ) DW, (3.11) ( D2 −a2 )( D2 −a2 −σ ) W + ( µHd 4πρ0v ) D ( D2 −a2 ) K = ( gαd2 v ) a2Θ, (3.12) ( D2 −a2 −P2σ ) X = − ( Hd η ) DZ, (3.13) ( D2 −a2 −σ ) Z = − ( µHd 4πρ0v ) DX, At the marginal state (σ = 0), then equations (3.9)-(3.13) become (3.14) ( D2 −a2 ) Θ = − ( βd2 κ ) W, (3.15) ( D2 −a2 ) K = − ( Hd η ) DW, (3.16) ( D2 −a2 )2 W + ( µHd 4πρ0v ) D ( D2 −a2 ) K = ( gαd2 v ) a2Θ (3.17) ( D2 −a2 ) X = − ( Hd η ) DZ, (3.18) ( D2 −a2 ) Z = − ( µHd 4πρ0v ) DX, and the boundary conditions (3.7) become (3.19)   Θ = 0, W = 0 at z = 0 and z = 1, DZ = 0, D2W = 0 at free boundary, Z = 0, DW = 0 at rigid boundary. By Taking the operator D for equation (3.15) and substituting in (3.16), we get (3.20) ( D2 −a2 )2 W −QD2W = ( gαd2 v ) a2Θ, Int. J. Anal. Appl. 19 (3) (2021) 445 where, (3.21) Q = µH2d2 4πρ0vη , is Chandrasekhar number [1]. Taking the operator (D2 −a2) for equation (3.20), and using equation (3.14), we get (3.22) (D2 −a2) [( D2 −a2 )2 −QD2 ] W = −Ra2W, where, (3.23) R = gαβd4 κv , is Rayleigh number. By Taking the operator (D2 −a2) for equation (3.18) and using (3.17), we get (3.24) [( D2 −a2 )2 −QD2 ] Z = 0, we must seek the solution of (3.22) and (3.24) subject to the boundary conditions (3.19). 4. Spectral Chebyshev tau Method Spectral Chebyshev tau method is a numerical method to solve the differential equations and the eigen values problems, (see [6], [7], [8] and [10]). To solve equations (3.22) and (3.24) subject to the boundary conditions (3.19), first convert the domain to Chebyshev polynomials domain [−1, 1], use the relation x = 2z − 1, if z ∈ [0, 1] implies x ∈ [−1, 1] and the derivative d dz = 2 d dx , d 2 dz2 = 4 d 2 dx2 , then (3.22), (3.24) and (3.19) become (4.1) (4D2 −a2) [( 4D2 −a2 )2 − 4QD2 ] W = −Ra2W, (4.2) [( 4D2 −a2 )2 − 4QD2 ] Z = 0, (4.3)   W = 0 at x = −1 and x = 1, DZ = 0, D2W = 0 at free boundary, Z = 0, DW = 0 at rigid boundary, where D = d dx . Let S = [( 4D2 −a2 )2 − 4QD2]W then we can write (4.1) - (4.3) as (4.4) [( 4D2 −a2 )2 − 4QD2 ] W −S = 0, Int. J. Anal. Appl. 19 (3) (2021) 446 (4.5) ( 4D2 −a2 ) S = −Ra2W, (4.6) [( 4D2 −a2 )2 − 4QD2 ] Z = 0, (4.7)   W = 0,S = 0 at x = −1 and x = 1, DZ = 0, D2W = 0 at free boundary, Z = 0, DW = 0 at rigid boundary. Now, expand W,S and Z as Chebyshev polynomials W = N∑ n=0 wnTn(x) = [ w0 · · · wN ]   T0 ... TN   = Wφ,(4.8) S = N∑ n=0 snTn(x) = [ s0 · · · sN ]   T0 ... TN   = Sφ,(4.9) Z = N∑ n=0 znTn(x) = [ z0 · · · zN ]   T0 ... TN   = Zφ,(4.10) where W, S and Z are row vectors represent the coefficients of W,S and Z respectively, and φ is the vector of chebyshev polynomials T0 up to TN . Furthermore we can expand the derivatives DW, D 2W and D4W as chebyshev polynomials as (4.11) DW = N∑ n=0 w(1)n Tn(x) = WDφ, (4.12) D2W = N∑ n=0 w(2)n Tn(x) = WD 2φ, (4.13) D4W = N∑ n=0 w(4)n Tn(x) = WD 4φ, Int. J. Anal. Appl. 19 (3) (2021) 447 similarly for the derivatives of S and Z (4.14) D2S = N∑ n=0 s(2)n Tn(x) = SD 2φ, (4.15) D2Z = N∑ n=0 z(2)n Tn(x) = ZD 2φ, Where D and D2 are (N + 1)×(N + 1) matrices represent the coefficients of the first and second chebyshev derivatives, (see [6] and [10] for more details of formulations of D and D2). D =   0 0 0 0 0 · · · 0 1 0 0 0 0 · · · 0 0 4 0 0 0 · · · 0 3 0 6 0 0 · · · 0 0 8 0 8 0 · · · 0 ... ... ... ... ... . . . 0 N 0 2N 0 2N · · · 0   ,(4.16) D2 =   0 0 0 0 0 0 · · · 0 0 0 0 0 0 0 · · · 0 4 0 0 0 0 0 · · · 0 0 24 0 0 0 0 · · · 0 32 0 48 0 0 0 · · · 0 0 120 0 80 0 0 · · · 0 ... ... ... ... ... ... . . . 0 0 N(N2 − 1) 0 N(N2 − 9) 0 N(N2 − 25) · · · 0   ,(4.17) and D4 = (D)4 = D2D2. By substitute (4.8)–(4.15) into equations (4.4)-(4.6) we get, (4.18) W [ 16D4 − (8a2 + 4Q)D2 + a4I ] φ− Sφ = 0, (4.19) S [ 4D2 −a2I ] φ = −Ra2Wφ, (4.20) Z [ 16D4 − (8a2 + 4Q)D2 + a4I ] φ = 0, Int. J. Anal. Appl. 19 (3) (2021) 448 where I is the identity matrix of order (N + 1). By taking the inner product with Tn,n = 0, 1, ..,N, for each equation in the system (4.18) – (4.20) and using the property of orthogonality of Chebyshev polynomials, we obtain 3(N + 1) equations as follows (4.21) [ 16w(4)n − ( 8a2 + 4Q ) w(2)n + a 4wn ] −sn = 0, n = 0, 1, . . . ,N (4.22) 4s(2)n −a 2sn = −Ra2wn, n = 0, 1, . . . ,N (4.23) 16z(4)n − (8a 2 + 4Q)z(2)n + a 4zn = 0. n = 0, 1, . . . ,N Rewriting these equations in matrices form as (4.24) W [ 16D4 − (8a2 + 4Q)D2 + a4I ] − S = 0, (4.25) S [ 4D2 −a2I ] = −Ra2W, (4.26) Z [ 16D4 − (8a2 + 4Q)D2 + a4I ] = 0, or, (4.27) AX = RBX, where A and B are square matrices of order 3(N + 1) given as A =   L1 −I 0 0 L2 0 0 0 L1   , B =   0 0 0 −a2I 0 0 0 0 0   , and X =   WT ST ZT  (4.28) where L1 = [ 16D4 − (8a2 + 4Q)D2 + a4I ]T and L2 = [ 4D2 −a2I ]T . Now returning to the boundary conditions (4.7), the p-derivative of Chebyshev polynomials at x = ±1 is given by the formula (4.29) dpTn dxp ∣∣∣∣ x=±1 = (±1)n+p p−1∏ k=0 n2 −k2 2k + 1 , For Free-Free boundaries, We have eight boundary conditions as W = 0, D2W = 0, S = 0 and DZ = 0 for x = ±1. Int. J. Anal. Appl. 19 (3) (2021) 449 BC1 : W(−1) = 0 ⇒ N∑ n=0 wnTn(−1) = 0 ⇒ N∑ n=0 (−1)n wn = 0, BC2 : W(1) = 0 ⇒ N∑ n=0 wnTn(1) = 0 ⇒ N∑ n=0 wn = 0, BC3 : D2W(−1) = 0 ⇒ N∑ n=0 wnT ′′ n (−1) = 0 ⇒ N∑ n=0 (−1)n+2 n2(n2 − 1) 3 wn = 0, BC4 : D2W(1) = 0 ⇒ N∑ n=0 wnT ′′ n (1) = 0 ⇒ N∑ n=0 n2(n2 − 1) 3 wn = 0, BC5 : S(−1) = 0 ⇒ N∑ n=0 snTn(−1) = 0 ⇒ N∑ n=0 (−1)n sn = 0, BC6 : S(1) = 0 ⇒ N∑ n=0 snTn(1) = 0 ⇒ N∑ n=0 sn = 0, BC7 : DZ(−1) = 0 ⇒ N∑ n=0 znT ′ n(−1) = 0 ⇒ N∑ n=0 (−1)n+1 n2zn = 0, BC8 : DZ(1) = 0 ⇒ N∑ n=0 znT ′ n(1) = 0 ⇒ N∑ n=0 n2zn = 0, similarly, for rigid-rigid boundaries, W = 0, DW = 0, S = 0 and Z = 0 for x = ±1. and for rigid-free boundaries, W = 0, DW = 0, S = 0 and Z = 0 for x = −1 rigid. W = 0, D2W = 0, S = 0 and DZ = 0 for x = 1 free. For each case of the boundary conditions, insert BC1 up to BC4 into the rows (N −2)th up to (N + 1)th of the first column in the matrix A in (4.27), BC5 and BC6 into the rows (2N + 1)th and (2N + 2)th of the second column in A, BC7 and BC8 into the rows (3N + 2)th and (3N + 3)th of the third column in A. The Int. J. Anal. Appl. 19 (3) (2021) 450 corresponding rows in the matrix B are zeros, then we can write the system (4.27) as follows  L1 −I 0 BC1 0 . . . 0 0 . . . 0 BC2 0 . . . 0 0 . . . 0 BC3 0 . . . 0 0 . . . 0 BC4 0 . . . 0 0 . . . 0 0 L2 0 0 . . . 0 BC5 0 . . . 0 0 . . . 0 BC6 0 . . . 0 0 0 L1 0 . . . 0 0 . . . 0 BC7 0 . . . 0 0 . . . 0 BC8     w0 ... wn s0 ... sn z0 ... zn   = R   0 0 0 0 . . . 0 0 . . . 0 0 . . . 0 0 . . . 0 0 . . . 0 0 . . . 0 0 . . . 0 0 . . . 0 0 . . . 0 0 . . . 0 0 . . . 0 0 . . . 0 −a2I 0 0 0 . . . 0 0 . . . 0 0 . . . 0 0 . . . 0 0 . . . 0 0 . . . 0 0 0 0 0 . . . 0 0 . . . 0 0 . . . 0 0 . . . 0 0 . . . 0 0 . . . 0     w0 ... wn s0 ... sn z0 ... zn   (4.30) where, L1, L2, 0 and I are the sub matrices of L1, L2, 0 and I respectively. Now we have generalized eigen value (4.30), we can find the minimum eigen values R(a), then the critical Rayleigh number values Rc and the corresponding wave number ac for various values of Q using MATLAB software, the results are illustrated in the next section. 5. Results and Conclusion • Spectral Chebyshev Tau method give results in full agreement with the results that obtained by the analytical solution given by Chandrasekhar [1]. • For the three cases of the boundaries free-free, rigid-rigid and rigid-free the critical Rayleigh number determine the stability. if R < Rc the motion is stable and no convection, when R = Rc it is stationary convection, and unstable motion when R > Rc. • It is also observed the effect of the magnetic field, as Q increases, the value of the critical Rayleigh number and the critical wave number also increases. Int. J. Anal. Appl. 19 (3) (2021) 451 Q Chandrasekhar [1] Present study ac Rc ac Rc 0 2.233 657.511 2.22 657.51 10 2.590 923.070 2.59 923.07 50 3.270 1762.04 3.27 1762.04 100 3.702 2653.71 3.7 2653.71 200 4.210 4258.49 4.21 4258.49 500 4.998 8578.88 5.00 8578.89 1000 5.684 15207.0 5.68 15207 2000 6.453 27699.9 6.46 27699.79 5000 7.585 63135.9 7.61 63135.48 10000 8.588 119832 8.61 119831.71 20000 9.706 230038 9.72 230038.48 40000 10.95 445507 10.96 445506.86 Table 1. The relation between between Rc and ac when the boundaries are free-free and Q = 0, 10, 50, . . . , 40000. wave number a 0 5 10 15 R a y li e g h n u m b e r R ×10 4 0 1 2 3 4 5 6 7 8 9 10 Free Free boundaries Figure 1. The relation between between Rc and ac when the boundaries are free-free and Q = 0, 10, 50, . . . , 40000. Int. J. Anal. Appl. 19 (3) (2021) 452 Q Chandrasekhar [1] Present study ac Rc ac Rc 0 3.13 1707.8 3.12 1707.77 10 3.25 1945.9 3.27 1945.75 50 3.68 2803.1 3.68 2802.01 100 4.00 3757.4 4.01 3757.23 200 4.45 5488.6 4.45 5488.54 500 5.16 10110 5.16 10109.78 1000 5.80 17103 5.81 17102.85 2000 6.55 30125 6.56 30124.81 4000 7.40 54697 7.39 54700.75 6000 7.94 78391 7.93 78441.25 8000 8.34 101606 8.30 101930.17 10000 8.66 124509 8.52 125856.50 Table 2. The relation between between Rc and ac when the boundaries are rigid-rigid and Q = 0, 10, 50, . . . , 10000. wave number a 0 5 10 15 R a y li e g h n u m b e r R ×10 4 0 1 2 3 4 5 6 7 8 9 10 Rigid Rigid boundaries Figure 2. The relation between between Rc and ac when the boundaries are rigid-rigid and Q = 0, 10, . . . , 10000. Int. J. Anal. Appl. 19 (3) (2021) 453 Q Chandrasekhar [1] Present study ac Rc ac Rc 0 2.68 1100.75 2.7 1100.81 2.5 2.75 1167.2 2.71 1167.45 12.5 2.97 1415.5 3.01 1415.81 25 3.17 1699.4 3.21 1699.79 50 3.45 2217.6 3.51 2217.98 125 4.00 3586.1 4.01 3585.85 250 4.50 5613.3 4.51 5612.93 500 5.10 9304.5 5.11 9303.95 1000 5.75 16119 5.71 16118.68 1500 6.20 22592 6.21 22591.30 2000 6.50 28879 6.51 28877.86 2500 6.75 35044 6.81 35043.57 5000 7.65 64847 7.61 64836.70 10000 8.65 122140 8.71 121512.72 Table 3. The relation between between Rc and ac when the boundaries are rigid-free and Q = 0, 2.5, 12.5, . . . , 10000. wave number a 0 5 10 15 R a y li e g h n u m b e r R ×10 4 0 1 2 3 4 5 6 7 8 9 10 Rigid Free boundaries Figure 3. The relation between between Rc and ac when the boundaries are rigid-free and Q = 0, 2.5, 12.5. . . , 10000. Int. J. Anal. Appl. 19 (3) (2021) 454 Conflicts of Interest: The author(s) declare that there are no conflicts of interest regarding the publication of this paper. References [1] S. Chandrasekhar, Hydrodynamic and hydromagnetic stability, Oxford University Press, Oxford, 1961. [2] S. Chandrasekhar, The instability of a layer of fluid heated below and subject to Coriolis forces, Proc. R. Soc. Lond. A. 217 (1953), 306–327. [3] S. Chandrasekhar, The instability of a layer of fluid heated below and subject to the simultaneous action of a magnetic field and rotation, Proc. R. Soc. Lond. A. 225 (1954), 173–184. [4] S. Chandrasekhar, D.D. Elbert, The instability of a layer of fluid heated below and subject to Coriolis forces. II, Proc. R. Soc. Lond. A. 231 (1955), 198–210. [5] A.J. Chorin, A numerical method for solving incompressible viscous flow problems, J. Comput. Phys. 2 (1967), 12–26. [6] J.J. Dongarra, B. Straughan, D.W. Walker, Chebyshev tau-QZ algorithm methods for calculating spectra of hydrodynamic stability problems, Appl. Numer. Math. 22 (1996), 399–434. [7] L. Fox, Chebyshev Methods for Ordinary Differential Equations, Computer J. 4 (1962), 318–331. [8] M.A. Hernandez, Chebyshev’s approximation algorithms and applications, Comput. Math. Appl. 41 (2001), 433–445. [9] A.A. Hill, Convection in porous media and Legendre, Chebyshev Galerkin methods, Ph.D. Thesis, Durham University, 2005. [10] D. Johnson, Chebyshev Polynomials in the Spectral Tau Method and Applications to Eigenvalue Problems, NASA Con- tractor Report 198451, 1996. [11] S. Talukdar, A study of some steady and unsteady MHD flow problems with heat and mass transfer, Ph.D. Thesis, Gauhati University, 2013. [12] D.Y. Tzou, Instability of Nanofluids in Natural Convection, J. Heat Transfer, 130 (2008), 072401. [13] B.-F. Wang, D.-J. Ma, Z.-W. Guo, D.-J. Sun, Linear instability analysis of Rayleigh–Bénard convection in a cylinder with traveling magnetic field, J. Crystal Growth, 400 (2014), 49–53. [14] D. Yadav, J. Wang, R. Bhargava, J. Lee, H.H. Cho, Numerical investigation of the effect of magnetic field on the onset of nanofluid convection, Appl. Therm. Eng. 103 (2016), 1441–1449. 1. Introduction 2. Linear Stability analysis 3. Normal Mode Analysis 4. Spectral Chebyshev tau Method 5. Results and Conclusion References