Plane Thermoelastic Waves in Infinite Half-Space Caused FACTA UNIVERSITATIS Series: Mechanical Engineering Vol. 14, N o 2, 2016, pp. 179 - 197 Original scientific paper INSTABILITY OF THE RAYLEIGH-BENARD CONVECTION FOR INCLINED LOWER WALL WITH TEMPERATURE VARIATION UDC 536.2 Sadoon Ayed 1 , Gradimir Ilić 2 , Predrag Živković 2 , Mića Vukić 2 , Mladen Tomić 3 1 University of Technology, Department of Mechanical Engineering, Baghdad, Iraq 2 University of Niš, Faculty of Mechanical Engineering, Niš, Serbia 3 College of Applied Technical Sciences, Niš, Serbia Abstract. This paper deals with an analysis of a two-dimensional viscous fluid flow between the two parallel plates inclined with respect to the horizontal plane, where the lower plate is heated and the upper one is cooled. The temperature difference between the plates is gradually increased during a certain time period after which it is temporarily constant. The temperature distribution on the lower plate is not constant in x-direction, there is a longitudinal sinusoidal temperature variation imposed on the mean temperature. We have investigated the wave number and amplitude influence of this variation on the subcritical stability and the onset of the Rayleigh-Bénard convective cells, by direct numerical simulation of 2D Navier-Stokes and energy equation. Key Words: DNS, Rayleigh Number, Boussinesq Approximation, Inclined Walls, Non- linear Stability Analysis 1. INTRODUCTION The Rayleigh-Bénard flow is a type of natural convection. It is one of the classical problems in fluid mechanics where a fluid is typically bounded by bottom and top walls which are heated and cooled, respectively (Fig. (1)). In this research, the fluid is subjected to the spatial temperature modulation at the lower wall and a small inclination of the both walls with respect to the horizontal plane. The reason for this convection is the temperature gradient in the vertical direction, which causes instable density stratification and, consequently, the fluid motion. Received September 29, 2015 / Accepted April 12, 2016 Corresponding author: Predrag Živković University of Niš, Faculty of Mechanical Engineering, Aleksandra Medvedeva 14, Niš, Serbia E-mail: pzivkovic@masfak.ni.ac.rs 180 S. AYED, G. ILIĆ, P. ŽIVKOVIĆ, M. VUKIĆ, M. TOMIĆ Fig. 1 Formation of Rayleigh-Bénard cells Rayleigh gave the first solution of the problem of the conducting state stability when the fluid in a gravitational field is bounded from above and below with constant but unequal temperatures. He obtained the critical value of a dimensionless parameter at which the flow starts. This parameter is named after him: the Rayleigh number. Another parameter, which measures the relative strength of the non-linearity in the momentum equations, is the Prandtl number. They are defined in the following way: 3 1 2 ( ) a g T T l R a     ; rP a   (1) where g is gravitational acceleration, β is thermal expansion coefficient, T1 is temperature of lower plate and T2 is the temperature of upper plate, l = 2H is the distance between the plates, υ is kinematic viscosity and a is thermal diffusivity. This is a non-dimensional parameter for the measure of the buoyancy/diffusive forces ratio. In the above definition of the Rayleigh number, the fluid properties are calculated at mean temperature Tm = (T1+T2)/2 as the best reference temperature. In our calculation T1 is a temporal and x-direction variable function as, according to Eq. (1), is the Ra number. The value of the critical Rayleigh number according to the linear stability theory is Rac = 1708 at wave number qc = 3.117. Beyond this value the fluid starts moving thus forming counter-rotating two-dimensional cells, the cross-section of which is almost square. The cellular flow becomes considerably more complicated as the Ra increases. The two-dimensional cells break up in three-dimensional cells, which appear hexagonal in shape when viewed from above. With larger Ra numbers, the cells multiply, become oscillatory and, finally, turbulent. The flow becomes turbulent at Ra = 10 4 and Pr = 1, for water Pr = 7 there is Ra = 10 5 , and for higher Pr numbers there is Ra = 10 6 -10 7 . Most of the above is valid for small temperature differences between the plates when the so-called Oberbeck-Boussinesq approximation is valid. In the Oberbeck-Boussinesq approximation all the fluid properties are considered constant except density, which is assumed to be a linear function of temperature: 2 2 [1 ( )]T T     (2) where ρ2 is fluid density at upper plate, T is fluid temperature. One of the first papers that analyzed the Rayleigh-Bénard flow with variable viscosity is that of Tipelkirch et al. [1]. Thereafter, many works have appeared treating the Rayleigh-Bénard flow with temperature dependant viscosity. Some of the authors have treated basic parameters [2,3], variable fluid viscosity [4,5], different enclosures [6,7,8], two-phase flows [9,10], etc. Some modern approaches to numerical solution can be found in the recent papers [8,10,11]. Instability of Rayleigh-Bénard Convection for Inclined Lower Wall with Temperature Variation 181 Probably the first work where the Rayleigh-Bénard flow was treated through assuming all the fluid properties as a function of temperature was that of Paolucci et al. [12]. The working fluid was air, its dynamic viscosity and thermal conductivity following the Sutherland law. The work of Fröhlich et all [13] and Severin and Herwig et al. [14] followed, which is the most recent work on the Rayleigh-Bénard convection considering all the fluid properties variable and angle-inclined [15-18]. Severin and Herwig calculated the critical Reynolds number using the method of asymptotic expansions and the results were valid only for small heat transfer rates. In this paper we consider a fluid flow for the Rayleigh number below the critical value against the wave number very close to the critical one. The numerical simulation of 2D Navier-Stokes equation in vorticity-stream function form is carried out for temperature- dependant thermophysical properties. In those cases there is spatial temperature modulation at the lower wall with small angle γ, where δm is amplitude and qm is wave number of the temperature modulation around some average value at the lower plate. 2. MATHEMATICAL MODEL At the Rayleigh-Bénard convection we can use three different approaches: Boussinesq approximation, low Mach number approximation, and compressible Navier-Stokes equations. The incompressible Navier-Stokes equations for Bousinesq approximation, continuity and energy equation read: 2 2 1 ( ) [1 ( )] 0 ( ) v v v p v T T g t v T v T T t                         (3) where kwjviuv   is velocity vector, p pressure, a thermal diffusivity and υ is kinematic viscosity. For the case of the fluid flow between two parallel plates inclined at an angle γ, relative to the horizontal plane, force g  may be expressed as: sin cosg g i g j   At this point in the research, considering the limited computational resources, the problem will be considered only in 2D domain. The boundary conditions are: ( , , ) 0, ( , , ) 0u x y H t u x y H t     (4) ( , , ) 0, ( , , ) 0v x y H t v x y H t     (5) 2 1 ( , , ) , ( , , )T x y H t T T x y H t T     (6) For the plane flow in the vorticity stream-function formulation we have: 2 ( ) 1 ( )( ) ( ) [1 ( )]( sin cos ) v v v p v T T g i g j t                      (7) 182 S. AYED, G. ILIĆ, P. ŽIVKOVIĆ, M. VUKIĆ, M. TOMIĆ Having in mind that per definition p and  υ, and that for plane flow we have kjviuv  0 , and x for any physical value , we obtain the following vorticity equation: 2 2( ) sin ( ) cos ( ) d d v T T T T gk t dy dx                      (8) Since 0 0x y z zi j k i j k          in two dimensional fluid flow, our equation after projection onto z direction becomes: 2 2( ) sin ( ) cos ( ) d d v T T T T g t dy dx                      (9) where we have denoted   z. The problem can be made dimensionless by using the following characteristic scales: L = H for length, T = H2/a for time and V = a/H for velocity. After rendering those equations dimensionless, we change the variables introduced in order to transform the physical domain into the computational one which is adopted to the Fourier-Chebyshev approximation, namely   x *  2 and   y   . The flow variables are reduced by scale factor to obtain the following coordinate transformation: * * * * 2 * * * * * * * * 2 2 * * * * 2 , , 1 1 1 1 , , x Lx Hx y Ly Hy a H v Vv v t Tt t H a L H L H V a a VL L HH                             (10) Substituting the above coordinate transformation in the vorticity equation, we have: * * * * * * * 2 2 2* * 1 1 ( ) sin ( ) cos ( ) V V VL Vv L L LTt L g d d T T T T L dy dx                                        (11) or in the following forms: * 2 * * * * 2 2* 2 3 * * ( ) sin ( ) cos ( ) V V V g d d v T T T T LT Lt L L dy dx                      (12) 2 * 2 * * * * 2 24 * 4 4 * * ( ) sin ( ) cos ( ) a a a g d d v T T T T HH t H H dy dx                      (13) Instability of Rayleigh-Bénard Convection for Inclined Lower Wall with Temperature Variation 183 If we multiply this equation with H 4 /a 2 the above equation yields: 3* * * * * 1 2 2 1 2 * * * 1 2 1 2 ( ) ( ) ( ) ( ) sin cos ( ) ( ) H g T T T T T Td d v a a T T T Tt dy dx                      (14) 3* * * * * 1 2 2 1 2 * * * 1 2 1 2 ( ) ( ) ( ) ( ) sin cos ( ) ( ) H g T T T T T Td d v a a a T T T Tt dy dx                       (15) Non-dimensional temperature is defined with: * 2 1 2 T T T T     (16) Rayleigh number by Eq. (1) and Prandtl number as ratio of kinematic viscosity to thermal diffusivity: Pr a   (17) Now, after substitution of Eqs. (17) and (16) in (15) we have: * * * * * * * * * ( ) Pr Pr sin cosv Ra t dy dx                       (18) Instead of continuity equation the definition of vorticity is used here:    (19) After scalar product with k we obtain: * * * 2 1 0 V v VL L L         (20) and after multiplication with L/V we obtain the above equation in dimensionless form: * * * 0    (21) The energy equation can be transformed in the following way: 2 2 2 ( ) ( )( ) ( ) T T v T T a T T t          (22) and dividing the equation by T1-T2, it yields: 2 2 2 1 2 1 2 1 2 ( ) ( ) ( ) ( ) ( ) ( ) ( ) T T T T T T v a t T T T T T T            (23) or: * * * ( )v a t          (24) 184 S. AYED, G. ILIĆ, P. ŽIVKOVIĆ, M. VUKIĆ, M. TOMIĆ Substituting the expressions for coordinate transformations: * * * * * 2 1 1 1 Vv a T t L L                (25) or: * * * * * 2 2 1 1a a v a t H HH H                (26) and after dividing by a/H 2 we obtain the energy equation in non-dimensional form: * * * * * * * ( )v t          (27) The nondimensional Oberbeck-Bousinesq system of equations to be solved (18), (21) and (27): * * * * * * * * * * * * * * * * * * * ( ) Pr Pr 0 ( ) v Ra t x v t                              The boundary conditions are the following: * * * * * * * * * * * * * * 0, , 0, , 0, 1 1, , 0, , 0, 1 {( , ) R| 0 2 1 1} at y y at y y D x y x y                             (28) And the initial conditions will be defined later. We have two boundary conditions for stream function and none for vorticity. The problem is solved by influence matrix method [19, 20], so that the first two equations can be solved simultaneously, and the numerical method is described in detail [21]. We describe here only the numerical method used for energy equation. 3. SOLUTION PROCEDURE The above set of equations should be solved simultaneously in time, so we describe the procedure only for the energy equation. For the solution of the described problem we use the Fourier-Chebyshev pseudo-spectral method, with the Fourier-Galerkin method for approximation in homogeneous x-direction and the Chebyshev collocation method for non-homogenous y-direction. Eq. (16) in developed form is: Instability of Rayleigh-Bénard Convection for Inclined Lower Wall with Temperature Variation 185 * * * 2 * * * * * * * 2* u v S t x y x                (29) Where S * is source term, which, in our case, is zero. The Fourier approximation in streamwise direction can be expressed by trigonometric polynomials in exponential form: * * * * * * * ( , , ) ( , ) K ikx k k K x y t y t e     (30) * * * * * * * ( , , ) ( , ) K ikx k k K S x y t S y t e    (31) * * * * * * * ( , , ) ( , ) K ikx k k K u x y t u y t e    (32) * * * * * * * ( , , ) ( , ) K ikx k k K v x y t v y t e    (33) * * * * * * * * * * * * * 1,* * ( , , ) ( , ) ( , ) K K ikx ikx k k K k Kk u x y t u y t e B y t e x x                (34) * * * * * * * * * * * * * 2,* * ( , , ) ( , ) ( , ) K K ikx ikx k k K k Kk v x y t v y t e B y t e x y                (35) Here * * ( , )y t designate Fourier coefficients for the exponential form of trigonometric polynomials in Eqs. (30)-(35) for any physical value  and 1i . In the following, we will always refer to dimensionless quantities; consequently * are suppressed from superscript, and after substitution of Eqs. (30), (31), (34) and (35) in (29), we obtain: * 2 2 2 2 ( , ) ( , ) ( , ) ( , ) ( , ) ( , ) K K K ikx ikx ikx k k K k K k Kk k K K K ikx ikx ikx k k k k K k K k K y t e u y t e v y t e t x y y t e y t e S y t e x y                                        (36) After differentiation: * 1, 2, 22 2 2 2 ( , ) ( , ) ( , ) ( ) ( , )( , ) ( , ) ( , ) K K K ikx ikx ikxk k k k K k K k K K K K ikx ilx ikx ikxk k k k K k K k K y t e B y t e B y t e t ik y t e e y t e S y t e y y                          (37) If multiply both the sides with the same weight function as the basis functions, we have: 186 S. AYED, G. ILIĆ, P. ŽIVKOVIĆ, M. VUKIĆ, M. TOMIĆ 1, 2, 2 2 2 ( , )( , ) ( , )( , ) ( , )( , ) ( , )( , ) ( , )( , ) ( , )( , ) ,..., 2 2 K K K ikx ilx ikx ilx ikx ilxk k k k K k K k K K K K ikx ilx ikx ilx ikx ilxk k k k K k K k K x x y t e e B y t e e B y t e e t k y t e e y t e e S y t e e t N N l k                            (38) If we apply the orthogonality condition, so that inner product of exponential functions yields: 1, 2, 2 0 ( , )( , ) ( , )( , ) ( , ) 0 for ( , ) 2 for K K K ikx ilx ikx ilxk k k k K k K k K ikx ilx ikx ilx y t e e B y t e e B y t t k l e e e e k l                    (39) Our system of equations is reduced to: 2 2 1, 2, 2 ( , ) ( , ) ( , ) ( , ) ( , ) ( , ) ,..., 2 2 k k k k k k x x y t B y t B y t k y t y t S y t t t N N l k                (40) For approximation in y-direction we use the Chebyshev-collocation method in the following way: 2 (2) 1, 2, , 0 ( , ) ( , ) ( , ) ( , ) ( , ) ( , ) ,..., , cos , 1,..., 1 2 2 yN k j k j k j k j j i j k j i jx x j y y y t B y t B y t k y t d y t S y t t N N k y j N N                  (41) where )2( ,ij d are elements of the Chebyshev differentiation matrix. The boundary points j = 0 and j = Ny are not included in the above system of equations. If we designate: 1, 2, ( , ) ( , ) ( , ) k j k j k j B y t B y t B y t  (42) For discretisation in time, the second order finite difference method in generalized form yields: 1 1 1 1 2 1 1 2 1 2 ( 2) 1 1 , 0 (1 ) ( ) 2 ( ) (1 ) ( ) 2 [ ( ) ( ) (1 ) ( )] [ ( ) (1 ) ( )] [ ( ) (1 ) ( )] ( ) (1 ) ( ) ,..., 2 y n n n k j k j k j n n n n n k j k j k j k j k j N n n n n j i k l k l k l k l i x y y y t B y B y B y k y y d y y S y S y N N k                                                , cos , 1,..., 1, , 0,1,..., 2 jx j y n t y y j N t n t n N N        (43) Instability of Rayleigh-Bénard Convection for Inclined Lower Wall with Temperature Variation 187 With these values the above system of equations becomes: 1 2 1 2 1 2 2, 0, 0, 1, 0, . . , 1i e                (44) We use the semi-implicit Adams-Bashworth backward differentiation method for which we have: 1 2 1 ( 2) 1 , 0 1 1 1 3 ( ) ( ) ( ) 2 4 ( ) ( ) ( ) 2[ ( ) ( )] 2 ,..., , cos , 1,..., 1, , 0,1,..., 2 2 y n N k j n n k j j i k j i n n k j k jn n n k j k j k j jx x j y n t y y k y d y t y y S y B y B y t N N k y j N t n t n N N                                 (45) The values on the right hand side are known, and the values on left hand side are unknown and have to be calculated for each time step. Values 1n k   have to be calculated for each wave number k = 0,…,Nx/2, and each collocation point yj = cos(j / Ny), j = 1,…,Ny. The generalized form of a set of Eqs. (43) can be solved for γ1=0, γ2=0, in the following way: 2 ( 2) 1 2 ( 2) , , 0 0 1 1 1 (1 ) 2 ( ) (1 ) (1 ) ( ) 2 2 (1 ) ( ) ( ) (1 ) ( ) [2 ( ) ( )] 2 ,..., , cos , 1,..., 1, 2 2 y yN N n n j i k j j i k j i i n n n n n k j k j k j k j k j jx x j y n y k d y k d y t t y S y S y y B y t N N k y j N t N                                                        , 0,1,..., t n t n N  (46) If introduce the following matrices: ( 2) 1 , 4 , 1,..., 1 1,..., 1 [ ] , [ ] 0,..., 0,..., x x i j i j y y i N i N I D d j N j N          (47) and scalars: 2 2 1 0 1 (1 ) 1 , (1 ), 2 2 k k t t t                       (48) The above system of equation is reduced to the following form: 1 1 1 1 4 0 1 4 1 1 1 1 1 1 [ ] [ (1 ) ] (1 ) 2 ( ) ( ) ,..., , , 0,1,..., 2 2 n n n k k k n n n n k k k j k j x x n t I D I D I I S I S B y B y N N k t n t n N                               (49) 188 S. AYED, G. ILIĆ, P. ŽIVKOVIĆ, M. VUKIĆ, M. TOMIĆ where: 0 0 0 [ ( ),..., ( )] , [ ( ),..., ( )] , [ ( ),..., ( )] n n n T n n n T n n n T k k k N k k k N k k k N y y S S y S y B B y B y     (50) If we introduce the following notation: 1 1 1 4 0 0 1 4 1 1 1 [ ], [ (1 ) ],B I D B I D B I               (51) We finally obtain these matrix equations: 1 1 1 1 1 0 1 1 1 (1 ) 2 0,..., 2 , , 0,1,..., n n n n n n n k k k k k k k x n t B B B I S I S B B k N t n t n N                      (52) This system of equations should be supplemented by boundary conditions. The horizontal walls are assumed to be isothermal with non-dimensional temperatures θhot = 1 and θcold = 0 at the bottom and top walls: 2 2 2 1 2 1 2 1 1 2 ( ) ( 1) , ( 1) 0 ( ) ( ) ( 1) , ( 1) 1 ( ) cold hot T T T y T y T T T T T y T y T T                   (53) The generalized Robin boundary conditions have the form: ( 1) ( 1) ( 1) ( 1) y y g y y y g y                             (54) After the implementation of the Chebyshev collocation method we have: (1) 0 , 0 (1) , 0 ( ) ( ) ( ) ( ) cos , 1,..., N N N j N j j N N N N j N j j j y d y g y d y g j y j N N                          (55) In these equations N = Ny. Boundary conditions on the upper and lower plates can be represented in this way: , , ( ,1, ) ( , ) ( , 1, ) ( , ) K ikx k k K K ikx k k K g x t g y t e g x t g y t e            (56) Instability of Rayleigh-Bénard Convection for Inclined Lower Wall with Temperature Variation 189 If we multiply by e ilx and implement orthogonality condition (39) we get: (1) 0 0, , 0 (1) , , 0 ( , ) ( , ) ( ) ( , ) ( , ) ( ) cos , 1,..., N k j k j k j N k N N j k j k j j y t d y t g t y t d y t g t j y j N N                          (57) In our case we have: 1 1(1) (1) (1) (1) 0 0,0 0,1 0,2 0, , (1) (1) (1) (1) ,,0 ,1 ,2 , ( ) ( ) n nk N k kN N N N N k N y d d d d g gd d d d y                                             (58) (1) (1) (1) (1) 0,0 0,1 0,2 0, 6 (1) (1) (1) (1) ,0 ,1 ,2 , 1 ,1 , 1 1 6 , , 0,..., , , 0,1,..., 2 N N N N N N n kn k k n n x k k n t d d d d D d d d d g G g N D G k t n t n N                                              (59) System of equations (52) and (59) can be solved by direct methods: 1 1 1 1 1 1 1 1 1 1 1 6 (1 ) 2 , 0,..., , , 0,1,..., 2 n n n n n n n k k k k k k k n n x k k n t B I S I S B B B H N D G k t n t n N                          (60) If we introduce the following notation:                 1 1 1 6 1 , n k n kn k G H F D B A   (61) The system is reduced to: 1 1 0,..., , 0,1,..., 2 n n k k x t A F N k n N       (62) This system of equations should be solved in a given time step prior the solution of momentum equation in the vorticity-stream function form Eqs. (16) and (19), and the calculated values of θn+1 should be used in Eq. (16). The system of Eqs. (16), (19) and (25) is solved for each time step n = 0,1,…,Nt for all wave numbers k = 0,1,…,Nx/2. 190 S. AYED, G. ILIĆ, P. ŽIVKOVIĆ, M. VUKIĆ, M. TOMIĆ 4. INITIAL AND BOUNDARY CONDITIONS We consider the initial condition for our simulation for the case of the Rayleigh-Bénard convection: * ( , , 0) 0, , ( , , 0) 0, , ( , , 0) 0 ( , ) {( , ) R| 0 2 , 1 1} x y t x y x y x y D x y x y                  (63) The temperature at the lower plate is raised gradually according to the following law: ( , 1, ) (1 sin cos ) sin 0, 0 / 2 ( , 1, ) (1 sin cos ), / 2 m m m m m m m m x y t q x q x t t x y t q x q x t                          (64) The temperature at the lower wall is not constant in x-direction, it depends on qm- wave number modulation, amplitude δm and frequency ω. Rayleigh number Ra measures average temperature gradient, while the additional spatial modulation is characterized by small amplitude δm and wave number qm. In the absence of forcing (δm = 0), convection cells with wave number qc, bifurcate only for Ra above critical Rayleigh number Rac. The onset of convection is characterized by parabolic linear stability - neutral curve in parameter space with Ra on abscise and q on ordinate axis. The neutral curve has its minimum at critical wave qc = 3.117 and Rayleigh number Rac = 1707.8. In contrast, for δm = 0, convection is unavoidable for any finite Ra in the simplest case in the form of “forced cells” with wave vector qm. The main goal of the present work is to provide direct numerical simulation of cells and their stability in presence of forcing with small amplitude δm ≈ 0 (0.01) and in ratio qm / qc = 1.2. The stability of the forced cells strongly depends on ratio qm/qc. It is very important to emphasize that our Ra number varies temporally and spatially, since temperature at the lower wall is described by the Eq. (64), while temperature difference T1-T2 in Eq. (1) varies with time t and x- coordinate. Our idea is to show the evolution of stream function, vorticity, velocity and temperature field in the transition period 0  t  , for ω = 1. The values we have chosen are subcritical values Ra=1000, qm=3.7 according to linear stability analysis, for Pr = 7, Δt = /300, δm = 0.01, number of Fourier modes K = 96, number of nodes Nx = 192, number of Chebyshev collocation points Ny = 192. 5. RESULTS OF THE NUMERICAL SIMULATION In a thin layer of the fluid heated from below, convection occurs as a steady pattern of two dimensional cells. The two-dimensional convection cells and stability properties are investigated in detail in [22] and [23]. For the heated layer corresponding to Ra > Rac the stable roll patterns occur only within a band of wave number centered approximately about critical wave number qc. Within the stable band, the cells realized do not necessarily have a preferred length scale. Indeed, their wavelengths appear to be dictated by the initial conditions used to select the cells and by the manner that the basic state temperature is prescribed spatially and temporally. The band is bounded on both sides by instabilities that pertain to changing the wavelength of the cells but without changing the form. As the induced cells acquire wave length too large or too small, instability will occur to shift their length scale back to a value close to the critical one. As the value of Ra increases, the cells at some point will become unstable and the convection structure will converge to a pattern with a more complex spatial and temporal structure. Instability of Rayleigh-Bénard Convection for Inclined Lower Wall with Temperature Variation 191 In this section we present the results of numerical simulation obtained by our 2D pseudo-spectral method developed in MATLAB code for the Navier-Stokes equation in stream function-vorticity formulation and energy equation, with the initial and boundary conditions described above and the walls inclined for 1 angular degree. In Fig. (2) we have shown the vorticity distribution for 4 instants of non-dimensional time t = /2, , 3/2 and 2. Since we have chosen qm = 3.7, such number of the pair of cells can be seen at t = /2 at the lower wall. This is the moment when temperature on the lower plate reaches its final value with its modulation in x-direction. The distribution in vorticity is not noticeably perturbed. The maximal and minimal values of vorticity are between -11 and 3. In the next instant of time t =  we can see an increase in the vorticity intensity, where it attains the values in the range from -30 to 12, but the distribution is already obviously perturbed in the direction of lifted walls. In the next instant of time t = 3/2 we can see even more vorticity distribution deformation, and the jump in the range of the values from -52 to 33. At t = 2, the convective terms in energy equation start to transport vorticity in upwards left direction and the external vorticity values increase further -75  ω  85, although the temperature at the lower wall is temporarily constant now, but with modulation in x-direction. The results of numerical simulation for stream function are given in Fig. (3). We can see that the range of values is increased not only for time period when temperature at the lower wall increases but also when the temperature attains its constant value. The spatial stream function distribution is slightly changed in time period 0  t  /2, but afterwards /2  t   convective terms attain values that significantly change the stream function’s intensity and distribution. The distribution pattern is more perturbed with the passage of time. The results of numerical simulation for velocity in x-direction are given in Fig. (4). The u-velocity keeps its form in time period 0  t  /2, but afterwards /2  t   buoyancy effect becomes significant, and convective terms attain values that significantly change the intensity and distribution of u-velocity. In time period /2  t   distribution pattern has been changing and becomes partly disordered. The range of values is increasing between -12 and 14. The fluid heating through the lower wall and fluid cooling through the upper wall are not in balance so that the increase in kinetic and thermal energy in the fluid flow is remarkable. The range of values in t = 3/2 is -22  u 25, and in t = 2 is-47  u  31, and the position of the extremely values is pushed from the lower wall toward the middle of the channel and to the direction of the lifted wall. The lower wall with its temperature modulation is the source of momentum in the initial stages and later the external values are shifted to the middle of the channel. The values near the lower wall are constantly increased. The lower wall serves as a source of momentum in x-direction and feeds the momentum in the middle of channel. In Fig. (5) we have shown the v-velocity evolution in period of time 0  t  2. We can see that the values of v-velocity constantly increase and that the upflow and downflow for t  /2 at the center of the channel are disturbed. In this case both buoyancy and shear have a destabilizing effect on the flow pattern and traveling cells are possible flow structures at t  /2. 192 S. AYED, G. ILIĆ, P. ŽIVKOVIĆ, M. VUKIĆ, M. TOMIĆ Fig. 2 Non-dimensional vorticity evolution in the forced Rayleigh-Bénard convection Instability of Rayleigh-Bénard Convection for Inclined Lower Wall with Temperature Variation 193 Fig. 3 Non-dimensional stream function evolution in the forced Rayleigh-Bénard convection 194 S. AYED, G. ILIĆ, P. ŽIVKOVIĆ, M. VUKIĆ, M. TOMIĆ Fig. 4 Non-dimensional u-velocity evolution in the forced Rayleigh-Bénard convection Instability of Rayleigh-Bénard Convection for Inclined Lower Wall with Temperature Variation 195 Fig. 5 Non-dimensional v-velocity evolution in the forced Rayleigh-Bénard convection 196 S. AYED, G. ILIĆ, P. ŽIVKOVIĆ, M. VUKIĆ, M. TOMIĆ We can see that the convection is unavoidable for Rayleigh number Ra = 1000, for δ = 0.01 (forced RBC flow) for each instant of time in the form of forced cells with wave number qm = 3.7 as can be seen in Fig. 5. These forced cells develop instability against resonant modes for Ra > Rac=1708 in wide range of qm/qc; in our case it is qm/qc = 3.7/3.12. Only for qm in vicinity of qc the forced cells remain stable up to fairly large Ra>Rac. In the case of inclined walls this cell has lost its stability at t  /2, in spite of the fact that our qm = 3.7 is in the vicinity of the critical wave number, and Rayleigh number Ra = 1000 is significantly below critical value Rac = 1708. In our simulation the instability is reached at lower values of Ra at the value of wave number which is close to the critical one. We have Ra < Rac but a periodic roll solution exists, since we have δm ≠ 0. For δm = 0 the periodic roll solution can exist only if Ra > Rac. Exploring the stability regime of cells is a demanding task, and even more difficult is the pattern selection, i.e. understating which is spontaneously chosen by system dynamics. The compression of the cells in the interior, which accompanies the enhanced cells curvature, causes the wave number in the cell center to exceed the instability boundary and leads to the formation of the dislocation pairs. The defects then travel toward the wall by climbing in the direction opposite to that of gravity acceleration. The result of this process is reduction of qm to the values less than qc and thus to the re- stabilization of the pattern. However, the domain walls emit new cells, which gradually re-compress the ones in the center and thus lead to persistent time dependence. 6. CONCLUSIONS In this paper, DNS of the Rayleigh-Bénard convection (RBC) is performed. In addition to the applied temperature gradient, the parallel walls are inclined and a sinusoidal temperature modulation with amplitude δm and wave number qm is introduced at the lower plate. While in the unforced RBC the heat conduction state becomes unstable at critical Rayleigh number Rac against convection cells with wave number qc, the forced cell solutions with wave number qm exist at any given Ra. Various destabilization mechanisms acting on the forced cells depend sensitively on the ratio qm/qc; here the results of simulation when this ratio is 1.2 are presented. Additional spatial forcing appears at upper plate and in the horizontal direction as the walls are inclined. This opens up new possibilities, in particular when the system is exposed to the presence of non- equal forcing wave numbers imposed at the two confining plates. It should be also mentioned that even in the absence of an applied uniform temperature gradient, a pure temperature modulation leads to periodic convection patterns. Acknowledgements: The authors would like to thank Prof. Dr. Miloš Jovanović for his invaluable insights into the problems of the Rayleigh-Bénard convection. The paper is part of the work within the project with the Government of the Republic of Serbia TR33036. Instability of Rayleigh-Bénard Convection for Inclined Lower Wall with Temperature Variation 197 REFERENCES 1. Tippelskirch, H., 1956, Über Konvektions-zellen, insbesondere in flüssigemSchwefel, Beitr. Phys. Atmos., 29, pp. 37-54. 2. Segel, L.A., Stuart, J.T., 1962, On the question of the preferred mode in cellular thermal convection, Journal of Fluid Mechanics, 13(2), pp. 289-306. 3. Busse, F.H., 1967, The stability of finite amplitude cellular convection and its relation to an extremum principle, Journal of Fluid Mechanics, 30(4), pp. 625-649. 4. Stengel, K. C., Oliver, D. S. and Booker, J. R., 1982, Onset of convection in a variable-viscosity fluid, Journal of Fluid Mechanics, 120, pp. 411-431. 5. Busse, F.H., Frick, H., 1985, Square-pattern convection in fluids with strongly temperature-dependent viscosity, Journal of Fluid Mechanics, 150, pp. 451-465. 6. Subha, S., 2012, Numerical stimulation of Rayleigh Bernard convection in wavy enclosures, International Journal of Instrumentation, Control and Automation (IJICA), 1(3-4), pp. 59-62. 7. Tracy, N.I., Crunkleton, D.W., 2012, Oscillatory natural convection in trapezoidal enclosures, International Journal of Heat and Mass Transfer, 55, pp 4498–4510. 8. Yigit, S., Poole, R. J., Chakraborty, N., 2015, Effects of aspect ratio on laminar Rayleigh–Bénard convection of power-law fluids in rectangular enclosures: A numerical investigation, International Journal of Heat and Mass Transfer 91, pp.1292–1307 9. Park, H.M., 2015, Rayleigh-Bénard convection of nanofluids based on the pseudo-single-phase continuum model, International Journal of Thermal Sciences, 90, pp. 267-278. 10. Matveev, L.V., 2016, Impurity transport in developed Rayleigh–Bénard convection, International Journal of Heat and Mass Transfer 95 pp. 15–21. 11. Zimmermann, C., Groll, R., 2015, Computational investigation of thermal boundary layers in a turbulent Rayleigh–Bénard problem, International Journal of Heat and Fluid Flow, 54, pp. 276–291. 12. Paolucci, S., Chenoweth, D.R., 1987, Departures from the Boussinesq approximation in the laminar Bénard convection, Physics of Fluids, 30(5), pp. 1561-1564. 13. Fröhlich, J., Laure, P., Peyret, R., 1992, Large departures from Boussinesq approximation in the Rayleigh-Bénard problem, Physics of Fluids A, 4(7), pp. 1355-1372. 14. Severin, J., Herwig, H., 2001, Onset of convection in the Rayleigh- Bénard flow under non-Boussinesq conditions: an asymptotic approach, Forschung im Ingenieurwesen, 66(4), pp. 185-191. 15. Ayed, S., Jovanović, M., Tomić, M., Ilić, G., Živković, P., Vukić, M., Dobrnjac, M., 2015, Instability Of Rayleigh-Bénard Convection Affected by Inclined Wall Temperature Variation, Proc. Int. Conf. DEMI, Banjaluka, pp. 397-402. 16. Sharma, A.K., Velusamy, K., Balaji, C., 2008, Interaction of turbulent natural convection and surface thermal radiation in inclined square enclosure, Heat and mass transfer, 44(10), pp. 1153-1170. 17. Crunkleton, D.W., Anderson, T.J., 2006, Numerical study of flow and thermal convection, International communications in heat and mass transfer, 33(1), pp. 24-29. 18. Beya, B.B., Lili, T., 2009, Transient Natural Convection in 3D Tilted Heated from Two Opposite Sites, International communications in heat and mass transfer, 36(6), pp. 604-613. 19. Kleiser, L., Schumann U., 1980, Treatment of incompressibility and boundary conditions in 3D numerical spectral simulation of plane channel flows, Proc. Third GAMM Conference Numerical Methods in Fluid Dynamics, Vieweg, Braunschweig, pp. 165-173. 20. Pulicani, J.P., 1988, A spectral multi-domain method for the solution of 1-D Helmholtz and Stokes-type equations, Computers and Fluids, 16(2), pp. 207-215. 21. Jovanović, M., 2009, Simulation of temporal hydrodynamic instability in plane channel flow, Proc. 2 nd Int. Congress of Serbian Soc. of Mechnics - IConSSM, Palić, pp. B09-B23. 22. Clever, R.M., Busse F.H., 1974, Transition to time dependant convection, Journal of Fluid Mechanics, 65(4), pp. 625-645. 23. Busse, F.H., Clever R.M. 1979, Instability of convection rolls in a fluid of moderate Prandtl number, Journal of Fluid Mechanics, 91(2), pp. 319-335.