- Email: [email protected]
Rayleigh–Benard convection of viscoelastic fluid
Applied Mathematics and Computation 136 (2003) 251–267 www.elsevier.com/locate/amc
Rayleigh–Benard convection of viscoelastic fluid H. Demir Department of Mathematics, Arts and Science Faculty, Ondokuz Mayıs University, 55139 Kurupelit, Samsun, Turkey
Abstract We consider two-dimensional unsteady Rayleigh–Benard convective motion of a viscoelastic fluid in a square cavity. The governing vorticity and energy transport equations are discretised by using finite difference approximations and solved numerically by either simple explicit or ADI methods, respectively. The two-dimensional convective motion is generated by buoyancy forces on the fluid in a square cavity, when the horizontal walls are maintained at two constant temperatures and the vertical walls are perfectly insulated. In this spirit we consider Rayleigh–Benard convection heat transfer in a square cavity shear-thinning and elastic effects (the first normal stress difference), as embodied in the dimensionless shear-thinning and elastic numbers, have an influence in shaping the flow field and determining the heat transfer characteristics with respect to the Rayleigh numbers. Their combined effect acts to increase and decrease the heat transfer as represented by the local Nusselt number with the results that overall heat transfer for the viscoinelastic and viscoelastic cases and all results in this paper have been documented first time for the viscoelastic fluid. Ó 2002 Elsevier Science Inc. All rights reserved. Keywords: Criminale–Erickson–Filbey (CEF) model; Thermal convection; Rayleigh–Benard convection; Weissenberg number
1. Introduction Two-dimensional natural convection heat transfer of a non-linear fluid in enclosures of rectangular cross-section is encountered in many practical E-mail address: [email protected] (H. Demir). 0096-3003/02/$ - see front matter Ó 2002 Elsevier Science Inc. All rights reserved. PII: S 0 0 9 6 - 3 0 0 3 ( 0 2 ) 0 0 0 3 6 - X
252
H. Demir / Appl. Math. Comput. 136 (2003) 251–267
situations due to its wide application in various engineering devices, which involve heat transfer across double glazing windows, solar energy, the electrical and nuclear industries, sterilisation of foods, building insulation. Such flows are also of interest in geophysics and can be applied in the circulation of the atmosphere and of magma in the Earth’s upper mantle. Many fluids exhibit non-Newtonian behaviour in industrial applications, but there are very few studies reported for natural convection based on non-Newtonian fluids in the literature for both viscoinelastic or viscoelastic fluids. Therefore those concerning the natural convection of viscoelastic fluids are almost non-existent. Review of the literature for two-dimensional free convection in rectangular enclosures can be found in [1–5]. In this present work, we examine two-dimensional Rayleigh–Benard convection for viscoelastic fluid namely CEF fluid in a horizontal enclosure of cross-section heated from below is studied when two vertical sides are insulated with shear-free boundary conditions. Rayleigh–Benard convection is one of the best known problems of fluid mechanics and has been subject of theoretical research and this problem has been investigated by many researchers with various aims due to its practical applicability. In this paper we present numerical results for Rayleigh–Benard convection with heating from below. In general, the use of this type of model would give realistic results in motions with small strain rates. In free convection, although the strains may be large, the strain rates are relatively small. In most buoyancy driven motions is also slow due to moderate temperature gradient. Therefore, CEF model is ideally suited for the study of this class of motions. This provides the motivation for the present work where the effect of the Weissenberg number on the flow field is explored numerically first time. For calibration of the code, firstly two-dimensional plane natural convection of the Newtonian fluid in a square cavity is solved. Excellent agreement is found with the work of De Vahl Davis [6], Torrance and Turcotte [7] and Mc Kensie [8]. Then for further calibration of the code, two-dimensional plane natural of a viscous pseudoplastic flow in a square cavity is studied. The recirculating flows are depicted properly in detail along with the temperature distributions which have not been given before in detail. Following this, to test the capabilities of the code Weissenberg effect on the flow field for CEF fluid is studied. Finally, the dependence of the Nusselt number on the characteristics is investigated.
2. Mathematical formulation We consider the flow regime in a two-dimensional square cavity as discussed in [4]. The two end walls are kept at different temperatures T1 > T0 and the remaining walls are insulated or conducted with Biot boundary conditions. The square enclosure considered is assumed to be very long in the third dimension.
H. Demir / Appl. Math. Comput. 136 (2003) 251–267
253
Then the governing differential equations subject to Boussinesq approximation are, DT ¼ kr2 T þ U; r u ¼ 0; qCv ð2:1Þ Dt q
Du ¼ r / þ qgaey ðT T0 Þ þ r S ; Dt
ð2:2Þ
where S is the extra stress and / is the reduced pressure field, / ¼ p p0 , p0 is the static pressure. q, a and k are the density, the coefficient of thermal expansivity and thermal conductivity, respectively, all evaluated at some average temperature. The operator D=Dt is the material time derivative o=ot þ u r. For the CEF model [4], the extra stress can be written as r
Sik ¼ 2gs ðc_ Þdik þ 4vðc_ Þdij djk 2fðc_ Þd ik ; where dik is the first rate-of-strain tensor given by 1 oui ouk dik ¼ þ ; 2 oxk oxi
Fig. 1. Streamline contours for viscous pseudoplastic fluid for We ¼ 0, Ra ¼ 104 .
ð2:3Þ
ð2:4Þ
254
H. Demir / Appl. Math. Comput. 136 (2003) 251–267
where gs ðc_ Þ is known as solvent viscosity and is a function of the second rateof-strain invariant c_ 2 ¼ 2dik dik . For our case we find that 2 2 ou ou ov _c2 ¼ 4 þ : ð2:5Þ þ ox oy ox In general the viscosity function gs ðc_ Þ should be continuous, well defined and able to attain a finite zero shear-rate viscosity. Consequently the viscosities are calculated using the Cross model g g1 1 ¼ ; g0 g1 1 þ ðkc: Þð1nÞ
ð2:6Þ
where g0 and g1 refer to the asymptotic values of viscosity at very low and very high shear rates, respectively, k is a constant parameter with the dimension of time and n is a dimensionless constant. Here, vðc_ Þ and 1ðc_ Þ generally are defined as 4vðc_ Þ ¼ N2 ðc_ Þ;
21ðc_ Þ ¼ N1 ðc_ Þ;
Fig. 2. Temperature contours for viscous pseudoplastic fluid for We ¼ 0, Ra ¼ 104 .
ð2:7Þ
H. Demir / Appl. Math. Comput. 136 (2003) 251–267
255
Fig. 3. Streamline contours for CEF pseudoplastic fluid for We ¼ 0:001, Ra ¼ 104 .
where N1 ðc_ Þ and N2 ðc_ Þ are material functions known as the primary and secondary normal stress coefficients, respectively. However, in this work, we use the following forms of this: N1 ¼ 2k1
ðg0 g1 Þ ; 1n 1 þ ðkc_ Þ
ð2:8Þ
N2 ¼ 0; and upper convected derivative in Eq. (2.3) is defined as r
d ik ¼
D dik Ldik dik LT ; Dt
ð2:9Þ
where L ¼ rVi
and
T
LT ¼ ðrVi Þ :
ð2:10Þ
Now, substituting Sik ¼ 2gs ðc_ Þdik þ rik in Eqs. (2.2) and (2.3), the stress constitutive equation and stress equation of motion in the form of elastic–viscous
256
H. Demir / Appl. Math. Comput. 136 (2003) 251–267
pseudoplastic–split-stress (EVp SS) which is similar to the elastic viscous split stress (EVSS) form of Rajagopalan et al. [10] can be written as Sik ¼ 2gs ðc_ Þdik þ rik ;
ð2:11Þ
where r
rik ¼ N1 ðc_ Þd ik ; ð2:12Þ ou ou ou o/ o o ou ov ou þ 2gs c g c þu þv ¼ þ þ q ot ox oy ox ox ox oy s oy ox o o þ rxx þ rxy ; ð2:13Þ ox oy ov ov ov o/ o ou ov þ qgaðT T0 Þ gs c_ q þu þv ¼ þ ot ox oy oy ox oy ox o ov o o 2gs ðc_ Þ ð2:14Þ þ þ rxy þ ryy : oy oy ox oy
Fig. 4. Temperature contours for CEF pseudoplastic fluid for We ¼ 0:001, Ra ¼ 104 .
H. Demir / Appl. Math. Comput. 136 (2003) 251–267
257
For a 2-D planar case, if the velocities are defined as u¼
ow ; oy
v¼
ow ; ox
x¼
ov ou ; ox oy
ð2:15Þ
and non-dimensional variables are introduced as x X ¼ ; L g¼
g ; g0
y Y ¼ ; L W¼
s¼
W g0 ; q
tg0 ; qL
U¼
ug0 ; qL
h¼
ðT T0 Þ ; ðT1 T0 Þ
V ¼
X¼
vg0 ; qL
xg0 ; qL2
ð2:16Þ
ð2:17Þ
then, taking the curl of the equation of motion which eliminates the unknown reduced pressure gradient, the resulting differences of the local velocity gradients can be abbreviated by the vorticity function x. Thus, the partial differential equation for stream function, vorticity, stresses and energy equation can be written as X ¼ r2 W;
Fig. 5. Streamline contours for viscous pseudoplastic fluid for We ¼ 0, Ra ¼ 105 .
ð2:18Þ
258
H. Demir / Appl. Math. Comput. 136 (2003) 251–267
2 oX 1 oX oX oh 1 ¼ þ F; H g ;X g U þV þ Gr os gRe oX oY oY Re where 2
H ðg ; XÞ ¼
o oX
oX g oX 2
o þ oY
oX g oY 2
;
1 F ¼ M ðXÞMðgÞ þ LðXÞLðgÞ Mðrxx ryy Þ Lðrxy Þ; 2 and the operators are defined as 2 o2 o ðÞ o2 ðÞ MðÞ ¼ 2 ; LðÞ ¼ ; oY 2 oX 2 oX oY oh 1 ¼ r2 h os Pr Re
ð2:19Þ
o o Ec ðU hÞ þ ðV hÞ þ U; oX oY Re
ð2:20Þ
ð2:21Þ
ð2:22Þ
ð2:23Þ
Fig. 6. Temperature contours for viscous pseudoplastic fluid for We ¼ 0, Ra ¼ 105 .
H. Demir / Appl. Math. Comput. 136 (2003) 251–267
259
"
2 2 # o oU o2 U o2 U oU oU oV oU 2 rxx ¼ 2WeN þU 2 þV ; os oX oX oX oY oX oY oX oY o 1 oU oV U o2 U o2 V V o2 U o2 V rxy ¼ 2WeN þ þ þ 2 þ þ os 2 oY oX 2 oX oY oX 2 oY 2 oX oY oV oU oV oU ; oY oY oX oX " 2 2 # o oV o2 V o2 V oV oU oV oU þV 2 2 þU : rxx ¼ 2WeN os oY oX oY oY oY oX oX oY ð2:24Þ On taking the non-dimensionalised form of the governing equations we have several non-dimensional parameters namely the Grashof number (Gr), the Prandtl number (Pr), the Eckert number (Ec) and Weissenberg number ðWeÞ. The Grashof number is a measure of the ratio of buoyancy force to viscous
Fig. 7. Streamline contours for CEF pseudoplastic fluid for We ¼ 0:001, Ra ¼ 105 .
260
H. Demir / Appl. Math. Comput. 136 (2003) 251–267
Fig. 8. Temperature contours for CEF pseudoplastic fluid for We ¼ 0:001, Ra ¼ 105 .
effects within the flow and defined as Gr ¼ gq2 bL3 ðT1 T0 Þ=g20 while the Prandtl number is a measure the ratio of the fluid’s kinematic viscosity to the thermal conductivity constant k and defined as Pr ¼ Cv g=k. The Eckert number is a measure of the viscous heating or the viscous dissipation part of the energy equation and defined as Ec ¼ U 2 =Cv ðT1 T0 Þ. The product of Pr and Gr gives the Rayleigh number Ra. The wavenumber is the measure of the elasticity of the fluid and is defined as k1 u=L and N is a dimensionless form of 2ðg0 g1 Þ= 1n ð1 þ ðkc_ ÞÞ . The dimensionless boundary and initial conditions that correspond to our problem are s60
whole space h ¼ X ¼ W ¼ U ¼ V ¼ 0
s>0
U ¼ V ¼ W ¼ X ¼ 0; rxy ¼
Weðo2 U =osoY Þ 1 þ ðkjoU =oY jÞ
1n
rxx ¼ ;
2Weðo2 U =oY 2 Þ 1 þ ðkjoU =oY jÞ1n
ryy ¼ 0
;
at X ¼ 0 and X ¼ 1:
H. Demir / Appl. Math. Comput. 136 (2003) 251–267
261
Fig. 9. Streamline contours for viscous pseudoplastic fluid for We ¼ 0, Ra ¼ 106 .
s > 0 U ¼ V ¼ W ¼ X ¼ 0; ryy ¼
2Weðo2 V =oX 2 Þ 1 þ ðkjoV =oX jÞ
1n
rxx ¼ 0;
rxy ¼
Weðo2 V =osoY Þ 1 þ ðkjoV =oX jÞ
1n
;
at Y ¼ 0 and Y ¼ 1:
Also we have h ¼ 0 at Y ¼ 1, h ¼ 1 at Y ¼ 0 and oh=oX ¼ 0 at X ¼ 1 and oh=oX ¼ 0 at X ¼ 0. Here T0 is used as the temperature of the cold wall. The local Nusselt numbers can be calculated from the temperature distribution as follows:
oh
: ð2:25Þ Nu ¼ oY wall
3. Method of solution and grid independence The modified equations were solved on a square mesh by a finite difference method; second-order central differences were used for all space variables and a fully implicit scheme for evaluating the time derivatives as described [11]. In the
262
H. Demir / Appl. Math. Comput. 136 (2003) 251–267
Fig. 10. Temperature contours for viscous pseudoplastic fluid for We ¼ 0, Ra ¼ 106 .
present work, the mixed equations (2.1) and (2.2) are solved by the alternating direction implicit (ADI) method developed by Peacemann–Rachford [12–15] and tridiagonal matrix algorithm for solving the discretised equations. The elliptic stream function is solved iteratively using the successive-over-relaxation (SOR) procedure. ADI method was used for this investigation because of earlier experience with this method in problems of buoyancy induced flow and heat transfer (e.g. [6]). Richardson’s extrapolation has been used; this leads to the high accuracy benchmark solution as described in [9]. Grid and time independence are carefully tested as shown in Fig. 1. We also tested our computer code against the earlier published numerical result of [6–8] for accuracy and excellent agreement was found.
4. Results and discussion In order to compare our original results, we present the streamline structure and temperature distribution for viscous pseudoplastic and CEF fluids together
H. Demir / Appl. Math. Comput. 136 (2003) 251–267
263
for the problem under consideration. In the results presented, Rayleigh, Eccert and Weissenberg numbers are set at 104 6 Ra 6 106 , 0 and 0 6 We 6 0:01, respectively. Despite the evident progress over the last few years, obtaining accurate numerical solutions (or worse, any solution at all) at high values of the Weissenberg number (We) remains a challenge. The high Weissenberg number problem, i.e. the divergence of conventional iterative schemes beyond some critical value of the Weissenberg number, has been reported in virtually all published works [12]. There are several reasons for this as discussed in [12]. Despite our work limited to the small values of the Weissenberg number, we found interesting structural changes in both streamline and temperature profiles as Weissenberg number increased from zero. For a viscous fluid, one cell filled whole cavity at Ra ¼ 104 which can be seen in Fig. 1 was found and it is symmetric both horizontally and vertically. The corresponding temperature profiles are shown in Fig. 2. In Fig. 2, warm fluid rises near the bottom wall and cold fluid descends along the top. As Weissenberg number increases, in this case, the CEF fluid flow results show that the fluid exhibits similar behaviour when We number is used, we have one vortex
Fig. 11. Streamline contours for CEF pseudoplastic fluid for We ¼ 0:001, Ra ¼ 106 .
264
H. Demir / Appl. Math. Comput. 136 (2003) 251–267
with the direction of rotation as shown in Fig. 3. The corresponding temperature profiles are shown in Fig. 4, where warm fluid rises near the bottom and left walls and cold fluid descends along the top and right walls. The last one of them indicates that the left of cavity is always warmer than right and warm fluid rises near the left wall. We are unable to get stable solution for We P 0:01, therefore we indicate this value of the We number as Wecrit: and it appeared to be impossible (as far as we know) to have stable solution beyond that value. This conclusion also indicates that the critical Weissenberg number is not artefact of our method. We have one cell vortex at Ra ¼ 105 in Fig. 5 for viscous fluid and it shows that the steady solution gives one large vortex which fills almost the whole cavity. The corresponding temperature profiles are given in Fig. 6. In Fig. 6 the temperature throughout most of the cavity is uniform obtaining a value which is average of the top and bottom walls. As Weissenberg number increases, in this case, the CEF fluid flow results show that the flow emerged as two unequal cells and they were not symmetric and one vortex dominant to other vortex for viscous fluid as shown in Fig. 7. The corresponding temperature profile is
Fig. 12. Temperature contours for CEF pseudoplastic fluid for We ¼ 0:001, Ra ¼ 106 .
H. Demir / Appl. Math. Comput. 136 (2003) 251–267
265
shown that it is almost symmetric vertically and warm fluid rises near the bottom while colder fluid located near the top wall as shown in Fig. 8. In this case, the CEF fluid flow results show that the fluid exhibits different behaviours when We number is used when Ra number increases to 105 . We are unable to get stable solution for We P 0:01, therefore we indicate this value of the We number as Wecrit: . We have observed one cell structure at Ra ¼ 106 for a viscous fluid and one vortex only and it is more compact as shown in Fig. 9. Fluid flow results show that the fluid exhibits similar behaviour. The corresponding temperature profiles are given in Fig. 10, and the temperature throughout most of the cavity is uniform obtaining a value which is average of the top and bottom walls. As We
Fig. 13. The effect of the Nusselt number on the heat transfer for viscous pseudoplastic fluid for We ¼ 0.
Fig. 14. The effect of the Nusselt number on the heat transfer for CEF pseudoplastic fluid for We ¼ 0:001.
266
H. Demir / Appl. Math. Comput. 136 (2003) 251–267
number is used there is not much difference between the CEF and viscous pseudoplastic fluids. The corresponding streamline and temperature profiles are shown in Figs. 11 and 12 respectively. The effect of the Nusselt number on the heat transfer is given in Fig. 14 for the CEF fluid and comparison with the viscous pseudoplastic fluid cases is shown in Fig. 13. There is no relationship between them. We believe that experimental results are necessary to clear this phenomena.
5. Conclusion In this paper, we have attempted to determine the effect of elasticity and shear rate dependent characteristics of the liquid to the non-Newtonian behaviour. The two-dimensional unsteady convective motion was generated by buoyancy forces on the fluids. Investigation was made numerically by using a second order accurate scheme for viscoelastic namely CEF model when N1 is in the form of Eq. (2.8) and finite difference method was used for the discretisation of time and space respectively. A range of Ra and We up to 106 and 0.001 has been investigated respectively for both vertical walls are perfectly insulated. The results show that overall heat transfer for the viscoinelastic and viscoelastic CEF fluid cases is quite close at Ra ¼ 106 . Beyond the critical We number, we have observed that instability occurs in the flow and the system becomes unstable and loses its equilibrium. These numerical simulations have enlarged our physical understanding of two-dimensional thermal convection between shearfree and insulated boundaries and it is not obvious that extensions to higher Rayleigh numbers will introduce any new phenomena. There are some points which still need to be clarified. For example, we believe that more works need to be done in this subject such as pseudospectral finite difference (PSFD) in order to explain bifurcation phenomena which have been noted in [13] and as is indicated [3].
References [1] A.V. Shoney, Natural convection heat transfer to viscoelastic fluid, in: Cheremisinoff (Ed.), Encyclopedia of Fluid Mechanics, Gulf, Houston, 1988. [2] B. Gebhard, Y. Jaluria, R.L. Mahajan, B. Sammakia, Buoyancy-Induced Flows and Transport, Hemisphere, Washington, DC, 1988. [3] H. Demir, F.T. Akyildiz, Unsteady thermal convection of a non-Newtonian fluid, Int. J. Eng. Sci. 38 (2000) 1923–1938. [4] H. Demir, Ph.D. Thesis, The stability properties of some rheological flows, school of accounting and mathematics, Division of Mathematical and Computing, The University of Glamorgan, Wales, UK, 1996. [5] H. Demir, Thermal convection of viscoelastic fluid with biot boundary conduction, Math. Comput. Simulation 56 (2001) 277–296.
H. Demir / Appl. Math. Comput. 136 (2003) 251–267
267
[6] G. De Vahl Davis, Natural convection of air in a square cavity: a bench mark numerical solution, Int. J. Numer. Meth. Fluids 3 (1983) 249. [7] K.G. Torrance, D.L. Turcotte, Thermal convection with large viscosity variations, J. Fluid Mech. 47 (1) (1971) 113–125. [8] D.P. McKensie et al., Numerical experiments on convection in the earth’s mantle, J. Fluid Mech. 62 (1974) 465–538. [9] R.B. Bird, W.E. Steward, E.D. Lightfoot, Transport Phenomena, Wiley, New York, 1960. [10] D. Rajagopalan, R.C. Amstrang, R.A. Brown, J. Non-Newtonian Fluid Mech. 36 (1990) 159. [11] J. Douglas Jr., D.W. Peacemant, Numerical solution of two dimensional heat flow problems, AIChE J. 12 (1966) 161. [12] L. Charles Tucker III, Fundamentals of Computer Modelling for Polymer Processing, Hanser, Munich, 1989. [13] D.A. Siginer, A. Valenzuaela-Rendon, Natural convection of viscoelastic fluid, in: Numerical Methods for Non-Newtonian Fluid Dynamics, FED-vol. 174, ASME, New York, 1994. [14] K.Y. Mortom, D.F. Mayers, Numerical Solution of Partial Differential Equations, Cambridge University Press, Cambridge, MA, 1994. [15] R.E. Scraton, Further Numerical Methods in Basic, Arnold, Paris, 1988.
Rayleigh–Benard convection of viscoelastic fluid H. Demir Department of Mathematics, Arts and Science Faculty, Ondokuz Mayıs University, 55139 Kurupelit, Samsun, Turkey
Abstract We consider two-dimensional unsteady Rayleigh–Benard convective motion of a viscoelastic fluid in a square cavity. The governing vorticity and energy transport equations are discretised by using finite difference approximations and solved numerically by either simple explicit or ADI methods, respectively. The two-dimensional convective motion is generated by buoyancy forces on the fluid in a square cavity, when the horizontal walls are maintained at two constant temperatures and the vertical walls are perfectly insulated. In this spirit we consider Rayleigh–Benard convection heat transfer in a square cavity shear-thinning and elastic effects (the first normal stress difference), as embodied in the dimensionless shear-thinning and elastic numbers, have an influence in shaping the flow field and determining the heat transfer characteristics with respect to the Rayleigh numbers. Their combined effect acts to increase and decrease the heat transfer as represented by the local Nusselt number with the results that overall heat transfer for the viscoinelastic and viscoelastic cases and all results in this paper have been documented first time for the viscoelastic fluid. Ó 2002 Elsevier Science Inc. All rights reserved. Keywords: Criminale–Erickson–Filbey (CEF) model; Thermal convection; Rayleigh–Benard convection; Weissenberg number
1. Introduction Two-dimensional natural convection heat transfer of a non-linear fluid in enclosures of rectangular cross-section is encountered in many practical E-mail address: [email protected] (H. Demir). 0096-3003/02/$ - see front matter Ó 2002 Elsevier Science Inc. All rights reserved. PII: S 0 0 9 6 - 3 0 0 3 ( 0 2 ) 0 0 0 3 6 - X
252
H. Demir / Appl. Math. Comput. 136 (2003) 251–267
situations due to its wide application in various engineering devices, which involve heat transfer across double glazing windows, solar energy, the electrical and nuclear industries, sterilisation of foods, building insulation. Such flows are also of interest in geophysics and can be applied in the circulation of the atmosphere and of magma in the Earth’s upper mantle. Many fluids exhibit non-Newtonian behaviour in industrial applications, but there are very few studies reported for natural convection based on non-Newtonian fluids in the literature for both viscoinelastic or viscoelastic fluids. Therefore those concerning the natural convection of viscoelastic fluids are almost non-existent. Review of the literature for two-dimensional free convection in rectangular enclosures can be found in [1–5]. In this present work, we examine two-dimensional Rayleigh–Benard convection for viscoelastic fluid namely CEF fluid in a horizontal enclosure of cross-section heated from below is studied when two vertical sides are insulated with shear-free boundary conditions. Rayleigh–Benard convection is one of the best known problems of fluid mechanics and has been subject of theoretical research and this problem has been investigated by many researchers with various aims due to its practical applicability. In this paper we present numerical results for Rayleigh–Benard convection with heating from below. In general, the use of this type of model would give realistic results in motions with small strain rates. In free convection, although the strains may be large, the strain rates are relatively small. In most buoyancy driven motions is also slow due to moderate temperature gradient. Therefore, CEF model is ideally suited for the study of this class of motions. This provides the motivation for the present work where the effect of the Weissenberg number on the flow field is explored numerically first time. For calibration of the code, firstly two-dimensional plane natural convection of the Newtonian fluid in a square cavity is solved. Excellent agreement is found with the work of De Vahl Davis [6], Torrance and Turcotte [7] and Mc Kensie [8]. Then for further calibration of the code, two-dimensional plane natural of a viscous pseudoplastic flow in a square cavity is studied. The recirculating flows are depicted properly in detail along with the temperature distributions which have not been given before in detail. Following this, to test the capabilities of the code Weissenberg effect on the flow field for CEF fluid is studied. Finally, the dependence of the Nusselt number on the characteristics is investigated.
2. Mathematical formulation We consider the flow regime in a two-dimensional square cavity as discussed in [4]. The two end walls are kept at different temperatures T1 > T0 and the remaining walls are insulated or conducted with Biot boundary conditions. The square enclosure considered is assumed to be very long in the third dimension.
H. Demir / Appl. Math. Comput. 136 (2003) 251–267
253
Then the governing differential equations subject to Boussinesq approximation are, DT ¼ kr2 T þ U; r u ¼ 0; qCv ð2:1Þ Dt q
Du ¼ r / þ qgaey ðT T0 Þ þ r S ; Dt
ð2:2Þ
where S is the extra stress and / is the reduced pressure field, / ¼ p p0 , p0 is the static pressure. q, a and k are the density, the coefficient of thermal expansivity and thermal conductivity, respectively, all evaluated at some average temperature. The operator D=Dt is the material time derivative o=ot þ u r. For the CEF model [4], the extra stress can be written as r
Sik ¼ 2gs ðc_ Þdik þ 4vðc_ Þdij djk 2fðc_ Þd ik ; where dik is the first rate-of-strain tensor given by 1 oui ouk dik ¼ þ ; 2 oxk oxi
Fig. 1. Streamline contours for viscous pseudoplastic fluid for We ¼ 0, Ra ¼ 104 .
ð2:3Þ
ð2:4Þ
254
H. Demir / Appl. Math. Comput. 136 (2003) 251–267
where gs ðc_ Þ is known as solvent viscosity and is a function of the second rateof-strain invariant c_ 2 ¼ 2dik dik . For our case we find that 2 2 ou ou ov _c2 ¼ 4 þ : ð2:5Þ þ ox oy ox In general the viscosity function gs ðc_ Þ should be continuous, well defined and able to attain a finite zero shear-rate viscosity. Consequently the viscosities are calculated using the Cross model g g1 1 ¼ ; g0 g1 1 þ ðkc: Þð1nÞ
ð2:6Þ
where g0 and g1 refer to the asymptotic values of viscosity at very low and very high shear rates, respectively, k is a constant parameter with the dimension of time and n is a dimensionless constant. Here, vðc_ Þ and 1ðc_ Þ generally are defined as 4vðc_ Þ ¼ N2 ðc_ Þ;
21ðc_ Þ ¼ N1 ðc_ Þ;
Fig. 2. Temperature contours for viscous pseudoplastic fluid for We ¼ 0, Ra ¼ 104 .
ð2:7Þ
H. Demir / Appl. Math. Comput. 136 (2003) 251–267
255
Fig. 3. Streamline contours for CEF pseudoplastic fluid for We ¼ 0:001, Ra ¼ 104 .
where N1 ðc_ Þ and N2 ðc_ Þ are material functions known as the primary and secondary normal stress coefficients, respectively. However, in this work, we use the following forms of this: N1 ¼ 2k1
ðg0 g1 Þ ; 1n 1 þ ðkc_ Þ
ð2:8Þ
N2 ¼ 0; and upper convected derivative in Eq. (2.3) is defined as r
d ik ¼
D dik Ldik dik LT ; Dt
ð2:9Þ
where L ¼ rVi
and
T
LT ¼ ðrVi Þ :
ð2:10Þ
Now, substituting Sik ¼ 2gs ðc_ Þdik þ rik in Eqs. (2.2) and (2.3), the stress constitutive equation and stress equation of motion in the form of elastic–viscous
256
H. Demir / Appl. Math. Comput. 136 (2003) 251–267
pseudoplastic–split-stress (EVp SS) which is similar to the elastic viscous split stress (EVSS) form of Rajagopalan et al. [10] can be written as Sik ¼ 2gs ðc_ Þdik þ rik ;
ð2:11Þ
where r
rik ¼ N1 ðc_ Þd ik ; ð2:12Þ ou ou ou o/ o o ou ov ou þ 2gs c g c þu þv ¼ þ þ q ot ox oy ox ox ox oy s oy ox o o þ rxx þ rxy ; ð2:13Þ ox oy ov ov ov o/ o ou ov þ qgaðT T0 Þ gs c_ q þu þv ¼ þ ot ox oy oy ox oy ox o ov o o 2gs ðc_ Þ ð2:14Þ þ þ rxy þ ryy : oy oy ox oy
Fig. 4. Temperature contours for CEF pseudoplastic fluid for We ¼ 0:001, Ra ¼ 104 .
H. Demir / Appl. Math. Comput. 136 (2003) 251–267
257
For a 2-D planar case, if the velocities are defined as u¼
ow ; oy
v¼
ow ; ox
x¼
ov ou ; ox oy
ð2:15Þ
and non-dimensional variables are introduced as x X ¼ ; L g¼
g ; g0
y Y ¼ ; L W¼
s¼
W g0 ; q
tg0 ; qL
U¼
ug0 ; qL
h¼
ðT T0 Þ ; ðT1 T0 Þ
V ¼
X¼
vg0 ; qL
xg0 ; qL2
ð2:16Þ
ð2:17Þ
then, taking the curl of the equation of motion which eliminates the unknown reduced pressure gradient, the resulting differences of the local velocity gradients can be abbreviated by the vorticity function x. Thus, the partial differential equation for stream function, vorticity, stresses and energy equation can be written as X ¼ r2 W;
Fig. 5. Streamline contours for viscous pseudoplastic fluid for We ¼ 0, Ra ¼ 105 .
ð2:18Þ
258
H. Demir / Appl. Math. Comput. 136 (2003) 251–267
2 oX 1 oX oX oh 1 ¼ þ F; H g ;X g U þV þ Gr os gRe oX oY oY Re where 2
H ðg ; XÞ ¼
o oX
oX g oX 2
o þ oY
oX g oY 2
;
1 F ¼ M ðXÞMðgÞ þ LðXÞLðgÞ Mðrxx ryy Þ Lðrxy Þ; 2 and the operators are defined as 2 o2 o ðÞ o2 ðÞ MðÞ ¼ 2 ; LðÞ ¼ ; oY 2 oX 2 oX oY oh 1 ¼ r2 h os Pr Re
ð2:19Þ
o o Ec ðU hÞ þ ðV hÞ þ U; oX oY Re
ð2:20Þ
ð2:21Þ
ð2:22Þ
ð2:23Þ
Fig. 6. Temperature contours for viscous pseudoplastic fluid for We ¼ 0, Ra ¼ 105 .
H. Demir / Appl. Math. Comput. 136 (2003) 251–267
259
"
2 2 # o oU o2 U o2 U oU oU oV oU 2 rxx ¼ 2WeN þU 2 þV ; os oX oX oX oY oX oY oX oY o 1 oU oV U o2 U o2 V V o2 U o2 V rxy ¼ 2WeN þ þ þ 2 þ þ os 2 oY oX 2 oX oY oX 2 oY 2 oX oY oV oU oV oU ; oY oY oX oX " 2 2 # o oV o2 V o2 V oV oU oV oU þV 2 2 þU : rxx ¼ 2WeN os oY oX oY oY oY oX oX oY ð2:24Þ On taking the non-dimensionalised form of the governing equations we have several non-dimensional parameters namely the Grashof number (Gr), the Prandtl number (Pr), the Eckert number (Ec) and Weissenberg number ðWeÞ. The Grashof number is a measure of the ratio of buoyancy force to viscous
Fig. 7. Streamline contours for CEF pseudoplastic fluid for We ¼ 0:001, Ra ¼ 105 .
260
H. Demir / Appl. Math. Comput. 136 (2003) 251–267
Fig. 8. Temperature contours for CEF pseudoplastic fluid for We ¼ 0:001, Ra ¼ 105 .
effects within the flow and defined as Gr ¼ gq2 bL3 ðT1 T0 Þ=g20 while the Prandtl number is a measure the ratio of the fluid’s kinematic viscosity to the thermal conductivity constant k and defined as Pr ¼ Cv g=k. The Eckert number is a measure of the viscous heating or the viscous dissipation part of the energy equation and defined as Ec ¼ U 2 =Cv ðT1 T0 Þ. The product of Pr and Gr gives the Rayleigh number Ra. The wavenumber is the measure of the elasticity of the fluid and is defined as k1 u=L and N is a dimensionless form of 2ðg0 g1 Þ= 1n ð1 þ ðkc_ ÞÞ . The dimensionless boundary and initial conditions that correspond to our problem are s60
whole space h ¼ X ¼ W ¼ U ¼ V ¼ 0
s>0
U ¼ V ¼ W ¼ X ¼ 0; rxy ¼
Weðo2 U =osoY Þ 1 þ ðkjoU =oY jÞ
1n
rxx ¼ ;
2Weðo2 U =oY 2 Þ 1 þ ðkjoU =oY jÞ1n
ryy ¼ 0
;
at X ¼ 0 and X ¼ 1:
H. Demir / Appl. Math. Comput. 136 (2003) 251–267
261
Fig. 9. Streamline contours for viscous pseudoplastic fluid for We ¼ 0, Ra ¼ 106 .
s > 0 U ¼ V ¼ W ¼ X ¼ 0; ryy ¼
2Weðo2 V =oX 2 Þ 1 þ ðkjoV =oX jÞ
1n
rxx ¼ 0;
rxy ¼
Weðo2 V =osoY Þ 1 þ ðkjoV =oX jÞ
1n
;
at Y ¼ 0 and Y ¼ 1:
Also we have h ¼ 0 at Y ¼ 1, h ¼ 1 at Y ¼ 0 and oh=oX ¼ 0 at X ¼ 1 and oh=oX ¼ 0 at X ¼ 0. Here T0 is used as the temperature of the cold wall. The local Nusselt numbers can be calculated from the temperature distribution as follows:
oh
: ð2:25Þ Nu ¼ oY wall
3. Method of solution and grid independence The modified equations were solved on a square mesh by a finite difference method; second-order central differences were used for all space variables and a fully implicit scheme for evaluating the time derivatives as described [11]. In the
262
H. Demir / Appl. Math. Comput. 136 (2003) 251–267
Fig. 10. Temperature contours for viscous pseudoplastic fluid for We ¼ 0, Ra ¼ 106 .
present work, the mixed equations (2.1) and (2.2) are solved by the alternating direction implicit (ADI) method developed by Peacemann–Rachford [12–15] and tridiagonal matrix algorithm for solving the discretised equations. The elliptic stream function is solved iteratively using the successive-over-relaxation (SOR) procedure. ADI method was used for this investigation because of earlier experience with this method in problems of buoyancy induced flow and heat transfer (e.g. [6]). Richardson’s extrapolation has been used; this leads to the high accuracy benchmark solution as described in [9]. Grid and time independence are carefully tested as shown in Fig. 1. We also tested our computer code against the earlier published numerical result of [6–8] for accuracy and excellent agreement was found.
4. Results and discussion In order to compare our original results, we present the streamline structure and temperature distribution for viscous pseudoplastic and CEF fluids together
H. Demir / Appl. Math. Comput. 136 (2003) 251–267
263
for the problem under consideration. In the results presented, Rayleigh, Eccert and Weissenberg numbers are set at 104 6 Ra 6 106 , 0 and 0 6 We 6 0:01, respectively. Despite the evident progress over the last few years, obtaining accurate numerical solutions (or worse, any solution at all) at high values of the Weissenberg number (We) remains a challenge. The high Weissenberg number problem, i.e. the divergence of conventional iterative schemes beyond some critical value of the Weissenberg number, has been reported in virtually all published works [12]. There are several reasons for this as discussed in [12]. Despite our work limited to the small values of the Weissenberg number, we found interesting structural changes in both streamline and temperature profiles as Weissenberg number increased from zero. For a viscous fluid, one cell filled whole cavity at Ra ¼ 104 which can be seen in Fig. 1 was found and it is symmetric both horizontally and vertically. The corresponding temperature profiles are shown in Fig. 2. In Fig. 2, warm fluid rises near the bottom wall and cold fluid descends along the top. As Weissenberg number increases, in this case, the CEF fluid flow results show that the fluid exhibits similar behaviour when We number is used, we have one vortex
Fig. 11. Streamline contours for CEF pseudoplastic fluid for We ¼ 0:001, Ra ¼ 106 .
264
H. Demir / Appl. Math. Comput. 136 (2003) 251–267
with the direction of rotation as shown in Fig. 3. The corresponding temperature profiles are shown in Fig. 4, where warm fluid rises near the bottom and left walls and cold fluid descends along the top and right walls. The last one of them indicates that the left of cavity is always warmer than right and warm fluid rises near the left wall. We are unable to get stable solution for We P 0:01, therefore we indicate this value of the We number as Wecrit: and it appeared to be impossible (as far as we know) to have stable solution beyond that value. This conclusion also indicates that the critical Weissenberg number is not artefact of our method. We have one cell vortex at Ra ¼ 105 in Fig. 5 for viscous fluid and it shows that the steady solution gives one large vortex which fills almost the whole cavity. The corresponding temperature profiles are given in Fig. 6. In Fig. 6 the temperature throughout most of the cavity is uniform obtaining a value which is average of the top and bottom walls. As Weissenberg number increases, in this case, the CEF fluid flow results show that the flow emerged as two unequal cells and they were not symmetric and one vortex dominant to other vortex for viscous fluid as shown in Fig. 7. The corresponding temperature profile is
Fig. 12. Temperature contours for CEF pseudoplastic fluid for We ¼ 0:001, Ra ¼ 106 .
H. Demir / Appl. Math. Comput. 136 (2003) 251–267
265
shown that it is almost symmetric vertically and warm fluid rises near the bottom while colder fluid located near the top wall as shown in Fig. 8. In this case, the CEF fluid flow results show that the fluid exhibits different behaviours when We number is used when Ra number increases to 105 . We are unable to get stable solution for We P 0:01, therefore we indicate this value of the We number as Wecrit: . We have observed one cell structure at Ra ¼ 106 for a viscous fluid and one vortex only and it is more compact as shown in Fig. 9. Fluid flow results show that the fluid exhibits similar behaviour. The corresponding temperature profiles are given in Fig. 10, and the temperature throughout most of the cavity is uniform obtaining a value which is average of the top and bottom walls. As We
Fig. 13. The effect of the Nusselt number on the heat transfer for viscous pseudoplastic fluid for We ¼ 0.
Fig. 14. The effect of the Nusselt number on the heat transfer for CEF pseudoplastic fluid for We ¼ 0:001.
266
H. Demir / Appl. Math. Comput. 136 (2003) 251–267
number is used there is not much difference between the CEF and viscous pseudoplastic fluids. The corresponding streamline and temperature profiles are shown in Figs. 11 and 12 respectively. The effect of the Nusselt number on the heat transfer is given in Fig. 14 for the CEF fluid and comparison with the viscous pseudoplastic fluid cases is shown in Fig. 13. There is no relationship between them. We believe that experimental results are necessary to clear this phenomena.
5. Conclusion In this paper, we have attempted to determine the effect of elasticity and shear rate dependent characteristics of the liquid to the non-Newtonian behaviour. The two-dimensional unsteady convective motion was generated by buoyancy forces on the fluids. Investigation was made numerically by using a second order accurate scheme for viscoelastic namely CEF model when N1 is in the form of Eq. (2.8) and finite difference method was used for the discretisation of time and space respectively. A range of Ra and We up to 106 and 0.001 has been investigated respectively for both vertical walls are perfectly insulated. The results show that overall heat transfer for the viscoinelastic and viscoelastic CEF fluid cases is quite close at Ra ¼ 106 . Beyond the critical We number, we have observed that instability occurs in the flow and the system becomes unstable and loses its equilibrium. These numerical simulations have enlarged our physical understanding of two-dimensional thermal convection between shearfree and insulated boundaries and it is not obvious that extensions to higher Rayleigh numbers will introduce any new phenomena. There are some points which still need to be clarified. For example, we believe that more works need to be done in this subject such as pseudospectral finite difference (PSFD) in order to explain bifurcation phenomena which have been noted in [13] and as is indicated [3].
References [1] A.V. Shoney, Natural convection heat transfer to viscoelastic fluid, in: Cheremisinoff (Ed.), Encyclopedia of Fluid Mechanics, Gulf, Houston, 1988. [2] B. Gebhard, Y. Jaluria, R.L. Mahajan, B. Sammakia, Buoyancy-Induced Flows and Transport, Hemisphere, Washington, DC, 1988. [3] H. Demir, F.T. Akyildiz, Unsteady thermal convection of a non-Newtonian fluid, Int. J. Eng. Sci. 38 (2000) 1923–1938. [4] H. Demir, Ph.D. Thesis, The stability properties of some rheological flows, school of accounting and mathematics, Division of Mathematical and Computing, The University of Glamorgan, Wales, UK, 1996. [5] H. Demir, Thermal convection of viscoelastic fluid with biot boundary conduction, Math. Comput. Simulation 56 (2001) 277–296.
H. Demir / Appl. Math. Comput. 136 (2003) 251–267
267
[6] G. De Vahl Davis, Natural convection of air in a square cavity: a bench mark numerical solution, Int. J. Numer. Meth. Fluids 3 (1983) 249. [7] K.G. Torrance, D.L. Turcotte, Thermal convection with large viscosity variations, J. Fluid Mech. 47 (1) (1971) 113–125. [8] D.P. McKensie et al., Numerical experiments on convection in the earth’s mantle, J. Fluid Mech. 62 (1974) 465–538. [9] R.B. Bird, W.E. Steward, E.D. Lightfoot, Transport Phenomena, Wiley, New York, 1960. [10] D. Rajagopalan, R.C. Amstrang, R.A. Brown, J. Non-Newtonian Fluid Mech. 36 (1990) 159. [11] J. Douglas Jr., D.W. Peacemant, Numerical solution of two dimensional heat flow problems, AIChE J. 12 (1966) 161. [12] L. Charles Tucker III, Fundamentals of Computer Modelling for Polymer Processing, Hanser, Munich, 1989. [13] D.A. Siginer, A. Valenzuaela-Rendon, Natural convection of viscoelastic fluid, in: Numerical Methods for Non-Newtonian Fluid Dynamics, FED-vol. 174, ASME, New York, 1994. [14] K.Y. Mortom, D.F. Mayers, Numerical Solution of Partial Differential Equations, Cambridge University Press, Cambridge, MA, 1994. [15] R.E. Scraton, Further Numerical Methods in Basic, Arnold, Paris, 1988.