On the solution of fluid flow and heat transfer problem in a 2D channel with backward-facing step

Full Text

On the solution of fluid flow and heat transfer problem in a 2D channel . . . Received: 4th May, 2017 / Revised: 8th June, 2017 / Accepted: 12th June, 2017 / First online: 25th June, 2017 Introduction. Permanent improvement of numerical techniques for modeling fluid dynamics and heat transfer requires adequate methods to assess their effectiveness. For these purposes, consideration of well-known test problems which solutions have all main features of real fluid flows is the best approach. On the one hand, for these problems the huge volume of experimental and theoretical results is saved up that provides their maximum reliability. On the other hand the critical parameters are unambiguously determined for them. The ‘intensification’ of these parameters increases the complexity of the problems. The efficiency evaluation of a new computing technology follows from this: the more critical parameters of the problem are ‘intense’, the more the technology is effective. The well-known example of such problem of hydrodynamics of incompressible viscous liquids is the classical lid-driven cavity problem. Reynolds number is the critical parameter for this problem. The problem of separated flow in a 2D channel with backward-facing step is one more such a problem [1]. The problem is characterized by a simple geometry. Its solution, in general, depends on two parameters: Reynolds number Re and expansion ratio ER, where ER is defined as the ratio of the full height of the channel to the height of its inlet segment. Prandtl number Pr is added when heat transfer is taken into account. Change of these few input parameters allows to receive radically various solutions of the problem-namely, separated flow patterns, which are characterized by the quantity, form, size, and position of the vortices. Furthermore, they can have a very complex structure. The maximum value of Reynolds number (Re = 3000 at ER = 2) was reached in the article [2], and expansion ratio (ER = 4 at Re 300) was reached in the works [3, 4]. Further increase in these parameters causes computational instability which, in particular, is discussed in [2]. E. Erturk recommends to reduce a mesh step for overcoming this instability (see, for example, [5]). But in this case, the systems of the linear algebraic equations (SLAE) with huge number of unknowns are generated, and more effective methods are needed to solve them. Of course, the number of grid nodes (the number of unknowns) can be reduced by reduction of the channel length. But in this case the open outlet border of the channel will approach close to zone of the stream perturbation. And, as is well known, the problem of formulation of conditions at open boundaries is still far from a complete solution. The correction algorithm for solving the incompressible fluid dynamics problems in areas with open borders is proposed in [6]. This algorithm makes it possible to position the open border almost anywhere in the fluid flow with minimal distortions of the problem solution in a narrow border zone. Its application has allowed to obtain solutions of the test problem of this investigation in the case of simultaneous increase in Re and ER up to Re = 3000 and ER = 10. 1. Problem statement and calculation method. The problem of viscous heat-conducting incompressible fluid flow in the plane channel with sudden expansion and the warmed-up bottom wall is considered. It is supposed that temperature differences of the stream are small, and volume forces are absent. As 363 F o m i n A. A., F o m i n a L. N. a result, it is assumed that liquid density, and also coefficients of viscosity and heat conductivity are constants. The scheme of fluid flow structure is presented in Fig. 1. Liquid inflows through the left border B4 , where its movement is a fully-developed flat Poiseuille flow between two parallel plates. It is also assumed that velocity and overheating of the stream practically do not change at the exit from the channel. Other borders of the channel are impermeable walls. The bottom wall B1 is heated to temperature Tw , and walls B2 , B3 and B5 are thermally insulated borders. Flow velocity at the border B4 smoothly increases from zero to the maximum value at initial times. As a result, liquid starts to move. In the final analysis, a stationary solution of the problem is looked for by a relaxation method. Figure 1. (color online) Scheme of flow in the plane channel with sudden expansion and the warmed-up bottom wall with designation of the domain borders The mathematical formulation of the problem in dimensionless form is a system of two-dimensional non-stationary Navier-Stokes equations: ∂u ∂u ∂u ∂p +u +v =- + Re-1 ∂t ∂x ∂y ∂x ∂v ∂v ∂v ∂p +u +v =- + Re-1 ∂t ∂x ∂y ∂y ∂u ∂v + = 0; ∂x ∂y ∂2u ∂2u + 2 , ∂x2 ∂y 2 ∂ v ∂2v + , ∂x2 ∂y 2 and stationary equation of heat balance: u ∂θ ∂θ +v = Pe-1 ∂x ∂y ∂2θ ∂2θ + ; ∂x2 ∂y 2 with initial conditions at t = 0: u = v = θ = 0; and conditions at the borders of the research area (see Fig. 1): B1 (y = 0, lc u = v = 0, θ = 1; ∂θ = 0; B2 (y = hc , 0 x lc ) : u = v = 0, ∂y ∂θ B3 (y = 1, 0 x L) : u = v = 0, = 0; ∂y 364 x L) : On the solution of fluid flow and heat transfer problem in a 2D channel . . . B4 (x = 0, hc y 1) : u = 6f (t)(y - hc )(1 - y)/(1 - hc )2 , v = 0, θ = 0; 0.5 sin 0.5π (2t - 1) + 1 , 0 t 1, t > 1; f (t) = 1; ∂θ = 0; ∂x ∂v ∂ 2 (1/θ) ∂u = 0. = = ∂x ∂x ∂x2 B5 (x = lc , 0 y hc ) : u = v = 0, B6 (x = L, 0 y 1) : Here t is time; x, y are the horizontal and vertical coordinates respectively; u, v are the horizontal and vertical components of the velocity V respectively; p is the pressure; θ = (T - T0 )/(Tw - T0 ) is the relative flow overheating; Pe = Pr · Re is Peclet number; Re = U H/ν is Reynolds number; Pr = ρ0 cp ν/κ is Prandtl number. The dimensional scales are: H is the full height of the channel (at the border B6 ); U is the average velocity of the incoming flow through the border B4 ; ν is the kinematic viscosity; ρ0 is the density; cp is the specific heat at constant pressure; κ is the thermal conductivity; T0 is the flow temperature at the inlet border of the channel. Each stationary solution is calculated in two stages: on the first stage dynamic characteristics of the stream u, v, p, and on the second stage the field of the overheating θ are defined. The non-stationary form of the momentum equations is caused by using a relaxation technique for reaching the stationary solution of the Navier-Stokes equations. On the basis of the received problem solution the distributions of the coefficient of friction Cf , and Nusselt number Nu along the bottom wall of the channel, as well as the coefficient of hydrodynamic resistance to the fluid flow λ are calculated. These coefficients are defined as follows [4, 7-9]. 1. Cf = τw /(1/2ρ0 U 2 ), where τw = ρ0 ν(∂u/∂n)w is the shear stress at the wall. Here (∂/∂n)w is the outward pointing derivative with respect to the wall. The modified coefficient of friction is proposed in [8] which does not depend on Reynolds number Cf∗ = 0.5 Cf Re. In the issue Cf∗ = (∂u/∂n)w . 2. Nu = αH/κ, where α is a local heat transfer coefficient in the formula for the heat flux q between heated surface and fluid stream qw = α(Tw - T0 ). Expression for Nusselt number follows from the definitions of heat flux qw = -κ(∂T /∂n)w and relative flow overheating θ [7] Nu = -(∂θ/∂n)w . 3. The coefficient of hydrodynamic resistance λ relates the pressure drop at the l distance from the entrance of the channel with the velocity head by the formula ∆P = λ lρ0 U 2 /(2H) [9], where ∆P is a difference of average pressure at the entrance to the channel and at l distance from it. Therefore, the formula for calculating λ = λ(x) will be as follows λ= 2 x 1 (p (0, y) - p (x, y)) dy. 0 365 F o m i n A. A., F o m i n a L. N. The technique of splitting of physical factors is applied to solve the of NavierStokes equations at each time step [10]. Herewith two modifications of the technique are used: 1) on the first step of splitting the pressure is taken into account from the previous time layer and implicit difference schemes for the movement equations are used; 2) on the second step of splitting the Neumann problem is formulated for the increment of pressure p , which is equal to a difference of pressure on the current and previous time layers. Details of the technology which allows to construct solution of the problem are presented in [6]. Difference approximation of the initial differential equations is carried out by the control-volume method with fifth-order power-law scheme of second-order accuracy in space and first-order accuracy in time [11]. In all calculations uniform grids are used with an identical step h along both coordinate axes. The resulting SLAE with respect to numerical vectors u = {uij }, v = {vij }, p = {pij }, and θ = {θij } are solved by the implicit iteration line-by-line recurrence method of the second order, accelerated in Krylov subspaces [12]. Dynamic flow characteristics are considered to be found when condition un - un-1 1 + vn - vn-1 ∆t Vn 1 1 < 10-5 is satisfied. Here superscript n is an index of a time level, ∆t is a time step which is defined from the formula ∆t = C min(h, Re h2 ), where h is a step of the spatial grid, C is Courant number. The optimal value of C = 30 has been determined from the results of computational experiments. 2. Algorithm validation. Calculations for different grids of the solution domain at Re = 1000, ER = 2, lc = 0.5, L = 10 showed that for h 1/300 the profiles of both components of velocity u and v are practically identical to each other. Computation errors of the dynamic and heat transfer parameters are given in the table. The relative errors of the parameters have been calculated in relation to their values for the grid step of 1/500. It is clearly seen that the grid step of h = 1/300 is sufficient to calculate these parameters with reasonable accuracy: in this case, in general, the differences of the parameters in neighborhoods of their local extremes are less than 1 %. As a result, all the solutions of the problem were received for this grid step. Also, the increasing of the step length up to lc = 1 has not led to significant differences in the solution [6], so all subsequent calculations have been carried out at the fixed lc = 0.5. Comparisons with the results of other investigators [2, 13-19] are presented in Figs. 2 and 3. The experimental and calculated data are shown in Fig. 2, while only the calculated results are presented in Fig. 3. It is not difficult to see that a good matching of all results take place. At the same time, it should be noted that all of the comparisons with the published data are made in the framework of nondimensionalization which has been accepted in the present work. Therefore, the same parameters and coordinates in the original papers can have, generally speaking, different values. 366 On the solution of fluid flow and heat transfer problem in a 2D channel . . . Maxima of the stream parameters for the set of grid steps at Re = 1000, L = 10, ER = 2, lc = 0.5 Grid step Cf∗ at B1 Cf∗ at B3 λ Nu at B1 1/100 1/200 1/300 1/400 1/500 -7.867 -7.104 -6.934 -6.864 -6.826 13.23 % 11.315 3.92 % 11.279 1.56 % 11.267 0.55 % 11.260 0.0 % 11.257 0.52 % 0.20 % 0.09 % 0.03 % 0.0 % -0.03642 -0.03545 -0.03525 -0.03517 -0.03513 3.54 % 0.91 % 0.33 % 0.11 % 0.0 % 9.316 8.929 8.852 8.824 8.811 5.42 % 1.33 % 0.46 % 0.15 % 0.0 % Figure 2. (color online) Profiles of the horizontal component of velocity u = u(y) in the different sections of the channel at Re = 1000, ER = 2: solid line - present data, • - E. Erturk [2], ◦ - B. F. Armaly et al. [19] Figure 3. (color online) Comparison of the results with literature data: (a) profiles of the vorticity at Re = 800, ER = 2 in the sections: 1 - x - lc = 7, 2 - x - lc = 15; (b) behavior of Nusselt number along the bottom wall: 1 - Re = 300, ER = 1.5, 2 - Re = 600, ER = 1.5, 3 - Re = 200, ER = 2 367 F o m i n A. A., F o m i n a L. N. 3. Calculated results and discussion. Effect of expansion ratio of the stream ER at Re = 500 on flow structure and fluid overheating is presented in Fig. 4. It is visible that increase of ER causes reduction of the perturbance zone length. At the same time, complexity of vortex structures of the flow grows. In particular, the vortex behind the step becomes two-level. Two centers of rotation appear in the vortices. As a result, saddle points are formed between these centers. It is interesting, that the maximum length of the vortices zone behind the step is realized at ER = 2 for the considered expansion ratios. As to the fluid overheating, due to the fact that an increase in ER growths the vertical dimension of the primary vortex behind the step, warm liquid layers approach the roof wall of the channel. As a result, on average, the liquid is more heated-up in all the channel at ER = 10 than at ER = 2. For the reason that in these calculations Re < 600, the results can be regarded as the simulation of the real fluid flow in the channel [3, 19]. In this sense, it would be interesting to obtain experimental confirmation of the existence of a two-level vortex behind the step at ER = 10. From the point of view of computing stability the variant at ER = 10 is the most complex task. Therefore, the influence of Reynolds number on the solution is investigated at the maximum step height. The corresponding results are presented in Fig. 5. These solutions can not already be regarded as a simulation of real fluid flows because Re > 600. However, from the point of view of computing techniques these results are very indicative. Especially in cases of Re 2000, when the right border of the channel cuts the bottom vortex. First, primary vortex behind the step is two-level in all cases of Re. Second, both primary vortex behind the step and secondary vortex at the roof wall are polycentric at Re 1000. Third, the complexity of vortex structure of the stream increases with increase in Re, and the zone of vortices stretches monotonically from the step up to the channel exit. Structure of the overheating field corresponds to structure of the fluid stream (Fig. 5, b). ‘Spots’ of uniform temperature arise in the places of large-scale circulations of the primary vortex. It is also visible that the upper secondary vortex does not allow the flow to get warm in the central part of the channel. Because of large Reynolds numbers (and Peclet numbers, respectively), convective heat transfer considerably surpasses diffusive one. Therefore, in this part of the channel the overheating area is literally pressed to the warm bottom wall. The vertical profiles of the horizontal component of velocity and fluid overheating for various Reynolds numbers at ER = 10 in section x - lc = 5 are presented in Fig. 6. The following numbering of curves is used in this figure and further in Figs. 7 and 8: 1 - Re = 500, 2 - Re = 1000, 3 - Re = 1500, 4 - Re = 2000, 5 - Re = 2500, 6 - Re = 3000. Fragments of the curves with negative values of u(y) detect wall vortices: the roof vortices for Re = 1000 and 1500 and the bottom ones behind the step for 2000 Re 3000. Positive components of the profiles correspond to the kernel of the main stream. It is visible, that the more Re (the less liquid viscosity), the more maximum of u. It is curious, that the curves are almost linear in the central part of the channel. The complex behavior of the u(y) profiles is explained by the interactions of numerous vortices as among themselves, and with the main stream of the fluid flow. Something similar one can observe in [20] where exact solutions of freeconvective fluid flows in a plane layer are considered. However, the driving forces 368 On the solution of fluid flow and heat transfer problem in a 2D channel . . . Figure 4. (color online) Structure of the fluid flow and the liquid overheating at Re = 500 for various ER: (a) streamlines; (b) overheating isolines, levels θ: from 0.95 to 0.05 with step 0.05, and 0.001 369 F o m i n A. A., F o m i n a L. N. Figure 5. (color online) Structure of the fluid flow and the liquid overheating at ER = 10 for various Re: (a) streamlines; (b) overheating isolines, levels θ: from 0.95 to 0.05 with step 0.05, and 0.001 370 On the solution of fluid flow and heat transfer problem in a 2D channel . . . Figure 6. (color online) The cross profiles of the stream parameters in section x - lc = 5 at ER = 10 for various Re: (a) horizontal component of velocity u; (b) fluid overheating θ Figure 7. (color online) Longitudinal profiles of flow parameters along the bottom wall at ER = 10 for various Re: (a) drop pressure (p0 is the pressure at the step basis; p0 has a unique value for each solution of the problem); (b) the modified coefficient of friction Figure 8. (color online) Longitudinal profiles of flow parameters along the bottom wall at ER = 10 for various Re: (a) the coefficient of hydrodynamic resistance; (b) Nusselt number at bottom wall of the channel 371 F o m i n A. A., F o m i n a L. N. are the buoyancy ones in this work in contrast to the present research. This causes more intensive mixing of heat in the field of consideration. As to the profiles of the liquid overheating (Fig. 6, b), existence of the local extrema should be noted here which arise because of vortex structure of the fluid flow. It is clear, that these extrema do not contradict the maximum principle in relation to the 2D overheating field. Also, attention should be paid to the uniform height sections of the curves which correspond to the above-mentioned ‘spots’ of uniform temperature in Fig. 5, b. Behavior of the drop pressure and the modified coefficient of friction along the bottom wall of the channel are shown in Fig. 7. First of all, it is necessary to note the increase in amplitude of the extrema of the profiles with increasing Reynolds number. It is easy to see, that the positions of the local maxima of the pressure correspond exactly to the right edge of the primary vortex behind the step. Moreover, the positions of the local minima correspond approximately to the right edge of the lower-level vortex behind the step. It is interesting, that despite the existence of another large vortex downstream near the bottom wall the pressure in this part of the flow changes a little. Profiles of the modified coefficient of friction Cf∗ behave in the same way. Also, the increase in amplitude with increasing Re takes place. However, the position of the extrema is slightly shifted downstream by 0.2-0.3 units. Therefore, the distinct matching of Cf∗ extremum position to the flow structure is not observed. It is clear, that positive values of the coefficient correspond to direct fluid flow, and negative values - to returned one. Small negative values of Cf∗ at the channel exit speak about the low intensity of the bottom vortex at the right boundary region. Behaviour of the coefficient of hydrodynamic resistance λ and Nusselt number along the channel are presented in Fig. 8. In the narrow entrance section of the channel, λ practically does not depend on longitudinal coordinate and well satisfies approximate ratio λ · ER ≈ 24/Re. This result coincides with the well-known theoretical solution of the problem of viscous fluid flow in a plane channel [9]. It is clear, that the jump of the coefficient of resistance is associated with sudden expansion of the channel. In other words, although the pressure in the transition from the narrow entrance of the channel to its expanded section varies smoothly, however, the cross sectional area of the channel increases abruptly at this time. Therefore, the coefficient of resistance changes abruptly too. The inverse relation of λ and Re is explained by the fact that generation of a stream of less viscous liquid requires smaller drop pressure. Curves of the behavior of Nusselt number say that the number of local maxima of Nu and their locations correlate with the number and positions of the rotation centers inside the primary vortex behind the step. In this sense, they naturally shift to the right with increasing Re because the primary vortex behind the step is stretched to the right. Conclusion. Numerical solutions of the problem of steady flow of heat-conducting viscous incompressible fluid in the plane channel with the backward-facing step and heated bottom wall have been presented in the article for Reynolds number 500 Re 3000, expansion ratio of the channel 1.43 ER 10, and for Prandtl number Pr = 0.71. Reliability of the results has been confirmed by comparison with literature data for corresponding values of the problem parameters. 372 On the solution of fluid flow and heat transfer problem in a 2D channel . . . Patterns of the solutions depending on ER at constant Re = 500 and depending on Re at constant ER = 10 have been given in the work. The last series of solutions (ER = 10) has been considered in detail as the most interesting previously not investigated case. In particular, the behaviors of the coefficients of friction, heat transfer and hydrodynamic resistance along the length of the channel have been analysed. The obtained results may be useful for comparisons in solving problems of this type. It is possible to draw the following conclusions on the basis of the obtained solutions in the considered range of change of the problem parameters: 1. The complex structure of the flow takes place: the vortices can have several centers of rotation and they can also have up to two levels with respect to the walls of the channel. 2. The growth of Re causes the increase in the length of all the vortices, while the dependence of the vortices lengths on ER is not monotonous one. 3. The increase in ER leads to a more uniform heating of the flow along the height of the channel. 4. The bottom vortices promote warming up of the flow and vice versa - roof vortices prevent the liquid from heating up. Competing interests. We have no competing interests. Authors’ contributions and responsibilities. Each author has participated in the article concept development and in the manuscript writing. The authors are absolutely responsible for submitting the final manuscript in print. Each author has approved the final version of manuscript. Funding. The research has not had any sponsorship.

About the authors

Alexander A Fomin

T. F. Gorbachev Kuzbass State Technical University

Email: fomin_aa@mail.ru
Cand. Phys. & Math. Sci.; Senior Researcher; Dept. for Development and International Cooperation 28, Vesennaya str., Kemerovo, 650000, Russian Federation

Liubov N Fomina

Kemerovo State University

Email: lubafomina@mail.ru
Cand. Phys. & Math. Sci., Associate Professor; Associate Professor; Dept. of Information and Computing Technology 6, Krasnaya st., Kemerovo, 650043, Russian Federation

References

Supplementary files