• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    Laminar natural convection characteristics in an enclosure with heated hexagonal block for non-Newtonian power law fluids

    2017-05-28 10:22:50KrunalGangawaneManikandan

    Krunal M.Gangawane*,B.Manikandan

    Department of Chemical Engineering,College of Engineering Studies,University of Petroleum and Energy Studies,Dehradun 248007,Uttarakhand,India

    1.Introduction

    The analysis of heat transfer processes along with mechanics of non-Newtonian fluids is of great significance due to its wide range of applications in nature.Many fluids which quite frequently observed in various industries exhibit non-Newtonian behavior,such as,shearthinning(Pseudoplastic),visco-elastic,shear-thickening(dilatant),etc.under suitable concentrations and/or flow conditions.Such materials include oil–water emulsions,froths and foams,sewage sludge,butter,gas liquid dispersions,biological fluids,paints,foodstuffs and many more[1,2].Moreover,it has also been acknowledged that,in day to day activities, fluids with non-Newtonian behavior are quite frequently encountered than Newtonian one.Non-Newtonian fluids differ from Newtonian one due to its non-linear dependence of shear stress with shear rate.As the dynamic viscosity varies with shear stress,its flow characteristics and thereby heat transfer become more complex[1,2].In order to delineate the non-Newtonian fluid behavior,several mathematical models are present(such as,power-law model,Carrieu–Yesuda model,etc.).Out of all these models,the power law(or Ostwaald de Waele)model has gained much popularity due to simplicity along with adequate representation of fluids with broad range of shear rates[3].

    Owing to its fundamental and pragmatic significance,natural convection in square/rectangular enclosures has been extensively studied and it is considered as one of the largely invaded area of heat transfer.The reason for overwhelming popularity of heat transfer in enclosures(closed as well open)con figurations is due to simple structure with ability to explore broad varieties of fluid flow and heat transfer fundamentals(boundary layer,vortex size and location,circulation of fluid,etc.).Typical examples which can be idealized as enclosures are,for instance,cooling systems of electronic gadgets,high performance building insulation,multi-shed structures,furnace,food processing(heating of various food stuffs like beans,carrot and potato chips,etc.),lubrication technologies, fluidized bed drying of fibrous substances,solar heat collectors,drying,etc.[4–11].

    On the other hand,detached obstacle of varying shapes(circular,square,triangular,etc.)is used for controlling fluid flow due to convection[12]in enclosure,which is needed for applications where fluid flow due to convection should be restricted or bifurcated(crystal growth,solidification,etc.).Though many studies have revealed the natural/forced/mixed convection characteristics in enclosure with adiabatic/isothermal obstacle of circular/square shape[12–20],no information,as much known to authors,is available for convection characteristics from hexagonal block.Moreover,the natural/forced or combined convection from heated objects in enclosures have practical findings such as,heat exchangers[21].Further,it is also known that the square block is more bluff body than circular one.The basic difference in fluid flow structure occurs due to the presence of cornered edges of block[22].Therefore,the presence of six,slightly inclined edges in hexagonal block can significantly affect momentum and heat transfer characteristics as compared to block of other shape.

    Additionally,in heat transfer studies,we come across various thermal boundary conditions,for example,constant temperature,heat flux,linear,sinusoidal,etc.These thermal conditions have remarkable in fluence on fluid flow and heat transfer characteristics in enclosures(due to change in boundary layer structure).As compared to natural or mixed convection characteristics from heated obstacles in cavity/enclosures,large amount of studies pertain to constant temperature condition.Very modest is known about the effect of other important thermal conditions(linear,heat flux,etc.)on the natural convection characteristics in enclosure.Therefore,in present work,difference in flow pattern and heat transfer in enclosure is analyzed for heated hexagonal block with either of two limiting thermal conditions of constant wall temperature(CWT)and uniform heat flux(UHF).It constitutes the main objective of the work,i.e.,exploring momentum and thermal transfer characteristics in an enclosure with centrally placed heated hexagonal block.But before presenting the detailed problem formulation and new simulation results,it is practicable to present the selected previous literature pertaining to natural/mixed convection in enclosures with/without heated body/block for Newtonian as well non-Newtonian fluids.

    2.Previous Work

    Literature covering heat transfer aspects in enclosures with heated blocks is large for Newtonian fluids[12–20].Many studies have already explored the influence of heated/adiabatic circular[21–33],sinusoidal[34],elliptical[35],triangular[36–39],square[12,40]blocks on natural/mixed convection characteristics for Newtonian fluids.However,limited information is available on natural convection characteristics for non-Newtonian fluids in enclosure with/without detached bodies of varying shapes and thermal conditions.For instance,Hussain and Hussein[25]investigated detailed convection characteristics due to buoyancy in cavity containing circular cylinder with uniform heat flux boundary condition.Lamsaadiet al.[41]presented the natural convection for non-Newtonian power law fluids in rectangular shallow cavity by numerical as well as analytical approach.They developed theoretical solution on the basis of non-Newtonian fluids.Cheng[42]reported influence of mixed thermal boundary condition(heat flux and temperature)on natural convection characteristics in a porous media from a vertical cone for power law fluids.Lower values of surface heat flux and temperature are obtained for shear thickening fluids.Few studies also have delineated the Bingham fluid flow features due to natural/mixed convection[43–46].In addition to that,Turanet al.[47]studied on natural convection analysis in differentially heated cavity for non-Newtonian power law fluids for range of governing parameters(Rayleigh number,Prandtl number,power law index,etc.).The averageNuapproaches unity for shear thickening fluids.Whereas,Khalifaet al.[48]chose horizontal porous cavity as computational model for analyzing non-Newtonian power law behavior for range of governing parameters.They predicted criteria for sub and supercritical onset of fluid flow.Sasmal and Chhabra[49]explored momentum and heat transfer characteristics due to natural convection from heated square cylinder.For numerical experimentation,they used range of parameters,such as,Grashof number(10≤Gr≤105),Prandtl number(0.72≤Pr≤100)and powerlaw index(0.3≤n≤1.8).Forotherwise similar conditions,shear thinning fluids showed higher heat transfer rate.Similar study is conducted by Chandra and Chhabra[50]for another shape of obstacle/block,i.e.,semi-circular cylinder.Subsequently,Matinet al.[51]studied natural convection for non-Newtonian power law fluids between two-square shaped eccentric duct annuli.There numerical study was focused on effect of eccentricity for range of governing parameters.They reported negligible effect of Prandtl number on averagedNuvariation.

    On the other hand,numerical study of double diffusive natural convection for Newtonian fluid in cavity with heated square body is reported by Nazariet al.[52]by using kinetic,lattice Boltzmann method.The results of the study shown increase in averageNuwith Rayleigh and Sherwood numbers.Selimefendigil and Oztop[53]shown comparison of effect of shape of obstacle placed inside cavity on flow structure due to magneto-hydrodynamics(MHD)natural convection.They investigated the effect of different shapes(circular,square and diamond)on convection characteristics in cavity for nano fluids.They reported reduction of averageNuvalues due to the presence of obstacle by 21.35%(circular),32.85%(diamond)and 34.64%(square)as compared to cavity without block.Recently,Ren and Chan[54,55]elucidated the influence of array of square blocks on natural convection characteristics by using lattice Boltzmann method and CUDA platform.It is observed that GPU can enhance the computation speed by a factor up to 20 as compared to the non-parallel CPU code.

    Therefore,from review of available literature,it can be said that immense knowledge is available for various aspects of hydrodynamics and heat transfer behavior in cavity with heated block for Newtonian fluids[23–40].Much less is known for non-Newtonian fluid flow characteristics in cavity with detached heated body for different thermal conditions.Few studies have delineated natural convection characteristics for non-Newtonian fluids in cavity with square/circular body[50,51].No information,as much know to authors,is available on natural convection characteristics from heated hexagonal block detached in enclosure for different thermal conditions and power law fluids.The present study aims to ful fill the discrepancy found in literature.In particular,the natural convection heat transfer characteristics from heated hexagonal block have been explored for range of Grashof number(103≤Gr≤106),Prandtl number(1≤Pr≤100)for power law index(0.5≤n≤1.5)for two limiting thermal conditions of constant wall temperature(CWT)and uniform heat flux(UHF).

    3.Problem Statement and Mathematical Formulation

    Fig.1.Schematics of physical as well as computational domain for problem under consideration along with mesh structure.

    The schematic of geometry and coordinate system is depicted in Fig.1.Geometry of problem consists of square enclosure containing centrally placed hexagonal block.The cavity contains fluid which is assumed to be steady,incompressible and obeys non-Newtonian power law.Hexagonal block is maintained either at constant wall temperature,Tw,(CWT)or uniform heat flux,q,(UHF).Vertical walls of cavity are exposed to ambient temperature(Tc<Tw),while horizontal walls are maintained at adiabatic condition.All walls of cavity are maintained at no-slip condition.Further,the thermo-physical properties of the working fluid are assumed to be independent of the temperature and the viscous dissipation effects as well as radiation heat transfer are assumed to be neglected.The extent of the variation of density with temperature is expressed by using Boussinesq approximation.Hence,the mass,momentum and energy equations for natural convection in non-dimensional form are expressed as follow[50].

    τijrepresents shear stress tensor.For fluids obeying non-Newtonian power law(or Ostwald-De Waele)model,it is expressed as,

    where,Dijrepresents the rate of deformation tensor for the two dimensional Cartesian coordinate and η is apparent viscosity which is also expressed as,

    Above expression contains,mandnwhich are consistency index and power law index,respectively.Forn=1,Eq.(5)simplifies to Newtonian fluids.The general energy equation for two-dimensional flow is expressed as follows:

    Non-dimensional equations for natural convection are obtained by following mathematical rearrangements of flow parameters,

    where,arched(*)variables are dimensional one.The consistent boundary conditions for the problem under consideration(Fig.1)are expressed in non-dimensional form as follows:

    ?West(0,y)and East(1,y)walls:These have no-slip walls and are exposed to ambient isothermal temperature(Tc).

    ?Hexagonal block:The centrally placed block is stationary and maintained either at CWT or UHF conditions.

    The dimensionless groups(for non-Newtonian fluids)appearing in above equations are represented as follows[50]:

    The numerical simulation of governing field equations(Eqs.(1)–(6))in conjunction with the above noted boundary conditions(Eqs.(8)–(11))yields the primitive variable fields,such as,velocity,pressure,temperature,etc.These parameters are further utilized and processed to obtain the engineering parameter,such as,streamlines,vorticity,Nusselt number,etc.These physical characteristics can be obtained as follows:

    ?Stream function(ψ):It is estimated from velocity field as,

    ?Vorticity(?)or tendency to rotate at point is expressed as,

    ?Nusselt number:The local Nusselt number on the surface of the heated

    hexagonal block is estimated by the following expressions:

    where,sis normal drawn to the edge of hexagonal block.The obtained local values have been further averaged over the surface of a block to estimate the surface averaged Nusseltnumber.Itis expressed as follows[57]:

    Itis known from dimensional analysis that the average Nusselt number is a function of the Grashof number,Prandtl number,and type of thermal boundary condition for the problem studied herein[57].The selection of simulation parameters(i.e.,grid,domain size,numerical tool,discretization type,etc.)has been discussed in subsequent section.

    4.Procedure for Numerical Simulation

    Governing partial differential equations(mass,momentum and energy)are discretized by using finite volume method(FVM).The solution of these equations has been obtained by using the ANSYS FLUENT(version 14)commercial CFD solver.The unstructured ‘triangular’cells of non-uniform grid spacing have been generated as shown in Fig.1.Furthermore,well-known semi-implicit method for the pressure linked equations(SIMPLE)scheme is used to solve the pressure–velocity decoupling.The set of algebraic equations are further solved by using the Gauss–Siedel(G–S)point-by-point iterative scheme.The momentum and energy terms are discretized by means of third order accurate QUICK(Quadratic Upwind Interpolation for Convection Kinetics).It uses 3-point to point upstream weighted quadratic interpolation for calculation of cell face values.The details of QUICK scheme can be found in[56].Iterative process is terminated for convergence criteria of 10?7based on the normalized residuals for each field equation.For general variable of ζ,the generalized convergence equation at node‘M’is expressed as follows[32]:

    where,M,N,ΓNBandBare node,neighboring node,in fluence coefficients for the neighboring nodes and constant part,respectively.Moreover,under relaxation of variables is performed in order to avoid divergence.These factors are maintained at 0.7,0.3 and 0.7 for velocity,pressure and temperature,respectively.In numerical simulation studies,proper choice of grid and domain sizes is a very important step as the accuracy of results is dependenton these parameters.Therefore,prudent choice of these parameters is a must.Next section examines grid independence study as well as validation of present numerical procedure.

    4.1.Choice of numerical parameters and validation

    The dependability as well as accuracy of the numerical simulation procedure is naturally dependent upon a right choice of optimal parameters such as sizes of the computational domain and computational grid.In this work,the size of computational domain is itself delineated by the problem;therefore,a thorough grid independence study has been carried out using five non-uniform of triangular grids(G1,G2,G3,G4 and G5).The details of grid size showing number of nodes and elements are given in Table 1.The influence of different grid sizes on average Nusselt number values is shown in Fig.2.It is observed that minimum and maximum deviations between G1 and G5 is 22%(Gr=103,n=0.5)and 13%(Gr=103,n=1.5),respectively.The average difference between G1 and G5 is estimated to be 2.4%.On the other hand,if deviation between G3 and G5 is concerned,maximum and minimum differences reduce to 5.8%(Gr=103,n=1.5)and 3.2%(Gr=103,n=0.5),respectively.It can also be clearly seen(Fig.2)that the effect of change in grid size after G3 is insignificant with enormous increase in the computational time as well as memory.Further,in changing grid size after G3,the computational time increases by almost two–three times.

    Thus,keeping in mind the factors such as,accuracy,computational time and space,triangular grid G3(nodes:18200;elements:35776)are found to be optimum and have been used for obtaining results.

    Table 1Different grid structures(number of nodes and elements)utilized for grid independence study

    Table 2Validation of numerical procedure based on average Nusselt number of heated wall of differentially heated cavity at Pr=0.71 and n=1

    Fig.2.Variation of average Nu for various grid structures at Gr=103;106 power law index of n=0.5,1.5 and Pr=100 for CWT and UHF thermal conditions.

    In order to ascertain the accuracy of present work,additional simulations are performed for problem of square cavity with centrally placed heated square block[51]for non-Newtonian power law fluids based on average Nusselt number of block.Table 3 shows comparison of simulatedNuvalues with that of[51]for Rayleigh number(104,105);power law index(n=0.6 and 1.4)forPr=10.The analysis of Table 3 reveals that maximum and minimum deviations are about 4.3%and 3.2%,forRa=104,n=1.4 andRa=105,n=0.6,respectively.An average error is estimated to be 0.45%.Such discrepancies in the numerical results are quite common in modeling studies.The factors responsible for such deviations are numerical method,grid size and structure,convergence criteria,etc.Hence,minding the above mentioned inadvertent factors affecting the accuracy of numerical results,the validation presented in Tables 2,3 ascertains the confidence in the accuracy and reliability of the present numerical simulation procedure.The results presented herein are,therefore,believed to be accurate and reliable within ±(3%–4%).Having gained the assurance in the accuracy for the present numerical simulation procedure,the subsequent sections presentthe influence of flow governing parameters(i.e.,Prandtl number,Grashof number,power law index,type of boundary condition,etc.)on the detailed natural convection heat transfer square cavity containing centrally placed heated hexagonal block in terms of the streamline and isotherm patterns and average Nusselt numbers.

    Table 3Comparison of average Nu values between present results with that of Martin et al.'s[51]based on problem of cavity with centrally placed heated square block for different Ra and power law index values at Pr=10

    5.Results and Discussion

    In this numerical work,the characteristics of natural convection in square cavity with centrally placed hexagonal block have been explored.The dependency of Nusselt number(dimensionless heat transfer coefficient)on governing parameters,such as,Grashof number(Gr),Prandtl number(Pr)and power law index(n)has been elucidated.In particular,the influence of range of various governing parameters that have been used for numerical experimentations is given as follows:

    Thermal conditions:constant wall temperature and uniform heat flux(or heat source).

    Extensive results are obtained for above mentioned flow parameters revealing the local and global flow and heat transfer characteristics(such as streamline,isotherm contours,average Nusselt number and Colburn factor of heat transfer(jnH)etc.)are presented and discussed herein the ensuing section.Further,for contour plots,the numbers of contour lines are kept constant at 15.

    Fig.3.Variation of streamlines with Grashof(Gr=I-103;II-104;III-105;IV-106)and Prandtl numbers(Pr=a-1;b-50 and c-100)for power law index of n=0.5 for constant wall temperature(CWT)condition.

    5.1.Streamline and isotherm patterns

    It is widely acknowledged that the flow and heat transfer from heated body/wall is affected mainly due to thickness of boundary layers,i.e.,momentum as well as thermal boundary layers.In natural convection studies,the parameters which in fluence these boundary layers are mainly,Grashof number,Prandtl number and power law index(in case of non-Newtonian fluids)[7].Therefore,these dependencies have been explored in ensuing paragraphs.

    Figs.3–8 depict the streamline structures for range of Grashof numbers(I-103;II-104;III-105;IV-106)corresponding to laminar flow and Prandtl numbers(a-1;b-50;c-100)for thermal conditions of CWT(Figs.3-5)followed by UHF(Figs.6-8)for power law index in increasing order,i.e.,n=0.5,1 and 1.5.It is obvious from these figures that the variations are largely due to Prandtl number followed by Grashof number and then power law index.In general, fluid near the vicinity of heated block experiences decrease in density due to rise in temperature.These lighter fluid elements then move in upward direction and strike the top wall(which is adiabatic).Fluid elements,then,start approaching cold vertical wall with gradual decrease in the temperature.This again accelerates the density of fluid,making it heavier.The fluid then settles at the bottom wall.These settled cold fluid elements then again approach the heated hexagonal body with gradual augmentation in temperature,thus completing the fluid circulation/rotation.As the fluid circulation remains mostly close the active walls(i.e.,from surface of block to front vertical wall),the core region of circulation remains nearly stagnant.This quasi-motionless region formed in the core of circulation can be named as convection cell or vortex.Furthermore,the flow is observed to be symmetric about the vertical axis,on either sides of heated block.Such behavior is quite common in the study of cavity with heated block.Such symmetric/bifurcated flow around vertical block are reported for other shapes(circular,square,etc.)also[12,25,32,55].Moreover,Grashof number is associated with the strength of buoyancy driven flow.The influence of increase in Grashof number(I-103;II-104;III-105;IV-106)can be clearly seen on the structure of streamlines(Figs.3-8).The significant variation in the shape of convection cell is observed withGrfor otherwise identical conditions.The increase inGrpromotes rapid circulation of fluid between active walls(i.e.,hot and cold surfaces).This rapid circulation stretches the convection cell and pulls it vertically upward,clearly indicating rise in natural convection phenomenon.Prandtl number shows significant effect on overall flow structure in cavity.The rise inPrcauses shrinking of the convection cells and also shifting towards the upper half of cavity.The dominant viscous forces slightly disturb the symmetric of flow due to thickened hydrodynamic boundary layer,which controls the flow circulation near heated hexagonal block.In addition to this,formationof plume can be noted forPr≥ 50,(Figs 3–5,b and c).It is due to the fact that higher Prandtl number fluids have dominant momentum/hydrodynamic boundary layer than that of thermal one.Simply saying,such fluids are dominated by viscous forces,thus impedes the thermal penetration in fluid.The formation of plume is reported for higherPrfluids,due to striking of fluid elements completing half circulation over the upper region of hexagonal body(as fluid circulation is mainly between part of hexagonal block facing vertical wall).The size of plume increases withGrdue to prominent circulation fluids for higher Prandtl number cases.The influence of power law index of the plume structure is observed to be nearly insignificant.For lowPrfluid(Pr=1),due to higher fluid movement and lower viscous effects,the formation of plume is not visible forGr≤104.

    Fig.4.Variation of streamlines with Grashof(Gr=I-103;II-104;III-105;IV-106)and Prandtl numbers(Pr=a-1;b-50 and c-100)for power law index of n=1.0 for constant wall temperature(CWT)condition.

    Fig.5.Variation of streamlines with Grashof(Gr=I-103;II-104;III-105;IV-106)and Prandtl numbers(Pr=a-1;b-50 and c-100)for power law index of n=1.5 for constant wall temperature(CWT)condition.

    Fig.6.Variation of streamlines with Grashof(Gr=I-103;II-104;III-105;IV-106)and Prandtl numbers(Pr=a-1;b-50 and c-100)for power law index of n=0.5 for uniform heat flux(UHF)condition.

    Fig.7.Variation of streamlines with Grashof(Gr=I-103;II-104;III-105;IV-106)and Prandtl numbers(Pr=a-1;b-50 and c-100)forpowerlaw index of n=1.0 for uniform heat flux(UHF)condition.

    Fig.8.Variation of streamlines with Grashof(Gr=I-103;II-104;III-105;IV-106)and Prandtl numbers(Pr=a-1;b-50 and c-100)for power law index of n=1.5 for uniform heat flux(UHF)condition.

    Streamline distributions for hexagonal block with uniform heat flux(UHF)case are shown in Figs.6–8,for power law indexes ofn=0.5,1.0 and 1.5 respectively.A remarkable difference is visible in fluid structure for both thermal conditions.For UHF condition, flow structure is observed to be vertically reverse of that of CWT case.Due to the presence of heat source and gravity effect,the flow circulation originates from the bottom.The convection cells,which are present in upper half of cavity for CWT case,now shift in lower half of cavity.One more difference can be noted that, fluid circulation for UHF case is much predominant than CWT for otherwise same parameters.For higher buoyancy forces,flow becomes chaotic with enlargement of convection cell.Thus, flow circulation becomes rapid and it limits close to the active walls of cavity.Similar way,formation of plume is also visible for UHF case even for low Prandtl fluids(Pr=1).As the direction of flow is reverse of CWT,the vertical plume formation takes place in downward direction only.The increase inGr,causes stretching of plume towards lower wall with reduction in its width.Another significant change on the streamline structure can be noted for higherPr(Pr≥50)andGr(Gr≥104),i.e., flow becomes parallel to the horizontal axis with elongation(horizontally)and shifting(towards the lower wall)of convection cells for a given power law index.Such behavior is observed due to accelerated flow due to high intensity buoyancy force.Thus from above illustration,it can be clearly stated that the flow structure is predominantly affected due to type of boundary condition,followed by Prandtl and Grashof number for given power law index.

    Fig.9.Variation of isotherms with Grashof(Gr=I-103;II-104;III-105;IV-106)and Prandtl numbers(Pr=a-1;b-50 and c-100)for power law index of n=0.5 for constant wall temperature(CWT)condition.

    Fig.10.Variation of isotherms with Grashof(Gr=I-103;II-104;III-105;IV-106)and Prandtl numbers(Pr=a-1;b-50 and c-100)for power law index of n=1.0 for constant wall temperature(CWT)condition.

    Fig.11.Variation of isotherms with Grashof(Gr=I-103;II-104;III-105;IV-106)and Prandtl numbers(Pr=a-1;b-50 and c-100)for power law index of n=1.5 for constant wall temperature(CWT)condition.

    Fig.12.Variation ofisotherms with Grashof(Gr=I-103;II-104;III-105;IV-106)and Prandtlnumbers(Pr=a-1;b-50 and c-100)forpower law index of n=0.5 for uniform heat flux(UHF)condition.

    Fig.13.Variation ofisotherms with Grashof(Gr=I-103;II-104;III-105;IV-106)and Prandtlnumbers(Pr=a-1;b-50 and c-100)forpower law index of n=1.0 for uniform heat flux(UHF)condition.

    On the other hand,Figs.9–14 show temperature pro files in cavity with heated hexagonal block.Fig.9 delineates isotherm contours for shear thinning fluids(n=0.5)and other similar condition considered herein.Atlow Grashofnumber(Gr=103),due to conduction dominant heat transfer,isotherm lines are nearly parallel to vertical cold walls.The augmentation inGr,causes rise in flow circulation,thus making isotherm lines parallel to horizontal walls.Also,crowding of isotherms around hexagonal block was noticed.It can also be seen for higher Prandtl number cases(Pr=a-1;b-50 and c-100)that the thickness of isotherm crowding decreases.It is due to the thinning of boundary layer,resulting in higher rate of heat transfer from block to surrounding fluid.Prandtl number increase also restricts the growth of thermal plume.Though,Prandtl number enhancement decreases the thermal penetration capacity of fluid,it enhances the flow circulation with increase in heat transfer rate.Further,shear thinning fluids show higher thermal gradient over the heated body,thus yielding higher heat transfer rate.The change of fluids from shear thinning(Fig.9)to shear thickening(Fig.10),the heat transfer rate decreases.

    Figs.12–14 describe the isotherm distribution in cavity with heated block for uniform heat flux(UHF)condition.As per previous discussion of streamlines(Figs.6–8),similar isotherm pattern are visible,i.e.,vertical reverse of CWT condition.The origin and formation of thermal plume takes place in vertically downward direction for UHF and for same range of parameters.

    Also,it can be marked that higher crowding of isotherms near block with lower boundary layer as compared to CWT.The dependency of isotherm as well as streamline found to have complex relationship with range of pertinent parameters(Grashof number,Prandtl number,power law index)as well as type of thermal condition.The influence of power law index on flow and isotherm patterns is explored in subsequent section.

    5.2.Effect of power law index

    The present work has explored the influence of non-Newtonian power law fluids with power law index values,n=0.5,1 and 1.5,thus covering Pseudoplastic(shear thinning),Newtonian and dilatant fluids(shear thickening),respectively.The influence of power law index on streamline patterns can be seen from Figs.3–5(CWT)and Figs.6–8(UHF)with increasing order ofn.It can be seen from Figs.3–8,power law index slightly modify flow pattern,which is visible from the change in recirculation/convection cell zone.The size of convection cell shrinks with increase in power law index indicating reduction in heat transfer rate.Another significant change can be noticed on size of plume formed above the top surface of hexagonal block.Similar effects(as that of streamline patterns)can be noticed on isotherm patterns(Figs.9–14).The size of thermal plume decreases withn.While less crowding of isotherm patterns near hexagonal block is observed,indicating low temperature gradient,thus,diminishing heat transfer rate due to convection.All in all,it can be clearly concluded from above mentioned discussions that heat transfer rate in cavity is less significantly affected by power law index than Grashof and Prandtl number.The heat transfer characteristics in cavity are delineated in subsequent section.

    Fig.14.Variation of isotherms with Grashof(Gr=I-103;II-104;III-105;IV-106)and Prandtl numbers(Pr=a-1;b-50 and c-100)for power law index of n=1.5 for uniform heat flux(UHF)condition.

    5.3.Heat transfer characteristics

    In this section,dependence of average Nusselt number(Eq.(14))as well as Colburn factor on natural convection heat transfer(jnH)is explored for range of governing parameters considered herein.Fig.15 represents the variation of average Nusselt number over the surface of hexagonal block for the considered range of conditions.Plot shows linear increase of average Nu with Grashof numbers,which is due to increased flow due to buoyancy for given power law index.The increase in power law index,i.e.,from shear thinning to shear thickening fluid,the density of isotherm clustering reduces,thereby reduces the value of average Nusselt number.Also,for shear thinning fluids,the remarkable deviation in the averageNuvalues is observed for both thermal conditions.Generally,hexagonal block with heat source or UHF condition yield higher Nu values than that of CWT condition especially for shear shinning fluids.The change of fluid from Newtonian(n=1)and thereby dilatant(n=1.5),and the discrepancy inNuvalues for two boundary conditions nearly vanished.Thus it can be proclaimed that the impact of boundary condition over the hexagonal block diminishes for shear-thinning(n=0.5),followed by Newtonian(n=1)and shear-thickening(n=1.5).

    In heat transfer studies,the rate of heat transfer is expressed by using average Nusselt number,which is principally a function of Reynolds/Grashof number and Prandtl numbers.All these three parameters(Nu,Re/Gr,Pr)can be combined using a single parameter,called Colburn factor for heat transfer(jH).For pure forced convection case,it can be expressed as follows:

    Fig.16 depicts the variation ofjnHwith Grashof number for different Prandtl numbers and power law indexes.It is known from our previous discussions,higherNuvalues are obtained for UHF condition than that of CWT.This effect is seen in variation ofjnHalso.For given value of power law index,higherjnHvalues are obtained for UHF case.Moreover,increasing power law index,jnHare predominantly affected at higherGrvalues,i.e.,Gr≥ 104.The higher effect of‘n’on Colburn factor pattern at higherGrnumbers is due to increased buoyancy effect,responsible of density driven flow.Thus,jnHhas complex dependence on physical parameters as well as type of boundary condition.

    In order to delineate the effect of shape of heated block placed centrally in cavity with otherwise similar boundary as well as parametric conditions,additional computations are performed for standard shapes of circular and square for similar grid structure and size.Comparison is shown in Fig.17 for Newtonian fluid(n=1),Pr=1100 and both thermal conditions.It is evident from observation of Fig.17 that higherNuvalues are obtained for hexagonal shape,followed by,square and circular.For CWT case,increase inNuwithGr,is more significant than UHF.Also,as per earlier discussions,it is clear that higherNuvalues are obtained for UHF than CWT thermal condition.The variationNuwithGrandPrfor different shapes is found to be linear/proportional,whereas,some complex dependency betweenNuwithGrandPrfor different shapes is observed for UHF case.Nushows very little variation withGr.HigherNuvalues are indicated by square block followed by hexagonal and circular blocks.

    Fig.15.Variation of Nu for different Grashof and Prandtl numbers at given power law index(n=0.5,1 or 1.5)for thermal conditions of CWT(lines)and UHF(symbols).

    Fig.16.Dependence of Colburn factor of natural convection(j nH)with Grashof and Prandtl number for different power law indexes(n=a-0.5;b-1 and c-1.5).

    5.4.Correlations for Nusselt number and Colburn factor of heat transfer

    In heat transfer studies,predictive correlations expressing the functional dependence of the experimental/numerical results are highly valued due to their possible use in the process engineering design and scientific developments.In the present numerical study,the functional dependence of the average Nusselt number(Nu)of heated hexagonal block and Colburn factor for heat transfer(jnH)on the flow governing parameters,namely,Grashof number(103≤Gr≤106),Prandtl number(1≤Pr≤100 and power law index(n=0.5,1 and 1.5),for thermal conditions of CWT and UHF is expressed as follows in general form(basic form for heat transfer correlations):

    Fig.17.Effect of shape of heated block(hexagonal,circular and square)on average Nu for range of Pr and Gr for n=1(Newtonian).

    Empirical correlations are developed by using in-house developed MATLAB code for least square curve fitting method.Table 4 presents the developed empirical correlations along with respectiveR2values over the ranges of the parameters considered herein.It can be clearly seen that thermal condition has remarkable influence on exponent of the concerned dimensionless parameters(i.e.,Gr,Prandn).It is observed from analysis of Table 4 that for both thermal conditions,Nusselt number shows linear variation withGrandPr;while it varies in inverse proportion with power law index(n).On the other hand,remarkable change is observed between correlations ofjnHfor both thermal conditions.In particular,exponent of‘n’shows opposite behavior for CWT and UHF.Fig.18 depicts the comparison between the numerical and predicted(by correlations given in Table 4)values of the average Nusselt number.It is clearly visible that the proposed predictive closure correlations predict the average Nusselt number(Nu)as well as Colburn factor(jnH)values within an acceptable level of deviations(±10%)from the computed values.

    Table 4Empirical correlations expressing functional dependence of Nu and j nH with Grashof(Gr),Prandtl(Pr)numbers and power law index(n)for CWT and UHF thermal conditions

    6.Concluding Remarks

    The numerical investigation of two dimensional,steady,laminar natural convection for power law fluids in square enclosure with built-in hexagonal heated block for two limiting thermal conditions of constant wall temperature(CWT)and uniform heat flux(UHF)is carried out by using finite volume method with SIMPLE algorithm and QUICK scheme for discretization of convection terms.Natural convection characteristics are analyzed for wide range of pertinent flow governing parameters such as,Grashof number(103≤Gr≤106)corresponding to laminar flow,Prandtl number(1≤Pr≤100)and power law index(0.5≤n≤1.5).For numerical results,optimum grid size is chosen by carrying out systematic grid independence study.The numerical validation of present work with literature has shown excellent agreement.For considered parameter range,hexagonal block with uniform heat flux(UHF)showed higher heat transfer rate.The dependence of average Nusselt number is more significantly affected by Prandtl number,followed by Grashof number,type of thermal boundary condition and lastly power law index.Further,rate of heat transfer diminishes from shear-thinning to Newtonian followed by shear thickening.The Colburn factor for natural convection heat transfer(jnH)is also determined.The UHF boundary condition shows higher values ofjnHthan CWT condition,proclaiming higher heat transfer rate for UHF case.In the end,numerical results are summarized by means of closure relationship revealing dependency of Nusselt number and Colburn factor with considered flow governing parameters.

    Fig.18.Parity plot between simulated and predicted(Eq.(20)and Table 4)values of average Nusselt number.

    Nomenclature

    ARaspect ratio

    bedge of equilateral hexagonal block,m

    Cpheat capacity,J·kg?1·K?1

    CWT constant wall temperature,K

    DaDarcy number

    GrGrashof number

    Hheight of cavity,m

    hheat transfer coefficient,W·m?2·K?1

    kthermal conductivity,W·m?1·K?1

    Llength of cavity,m

    MaMach number

    Nulocal Nusselt number

    Nuaverage Nusselt number

    Ppressure,N·m?2

    PePeclet number

    PrPrandtl number

    qwheat flux,W·m?2

    RaRayleigh number

    ReReynolds number

    RiRichardson number

    Ttemperature

    Tcambient temperature

    Trefreference temperature,K

    Twsurface temperature of triangular block

    T*dimensional temperature,K

    ΔTtemperature difference,K

    UHF uniform heat flux,W·m?2

    ux,uy,velocity components,m2·s?1

    Vcharacteristics velocity,m2·s?1

    x,yco-ordinates,m

    α thermal diffusivity,m2·s?1

    β coefficient of thermal expansion,K?1

    μ dynamic viscosity,N·s·m?2

    ν kinematic viscosity,m2·s?1

    ρ density,kg·m?3

    ρ average density determined atTref,K

    ? vorticity,m2·s?1

    Ψ stream function,m2·s?1

    Ω angular rotational speed,r·min?1

    [1]H.Xu,S.J.Liao,Laminar flow and heat transfer in the boundary-layer of non-Newtonian fluids over a stretching flat sheet,Comput.Math.Appl.57(2009)1425–1431.

    [2]J.Su,J.Ouyang,X.Wang,B.Yang,W.Zhou,Lattice Boltzmann method for the simulation of viscoelastic fluid flows over a large range of Weissenberg numbers,J.Non-Newtonian Fluid Mech.194(2013)42–59.

    [3]R.P.Chhabra,Bubbles,Drops,and Particles in Non-Newtonian Fluids,CRC Press,Boca Raton,FL,2007.

    [4]S.S.Mendu,P.K.Das,Flow of power-law fluids in a cavity driven by the motion of two facing lids-a simulation by lattice Boltzmann method,J.Non-Newtonian Fluid Mech.175-176(2012)10–24.

    [5]K.M.Gangawane,R.P.Bharti,S.Kumar,Lattice Boltzmann analysis of natural convection in a partially heated open ended enclosure for different fluids,J.Taiwan Inst.Chem.Eng.49(2015)27–39.

    [6]K.M.Gangawane,R.P.Bharti,S.Kumar,Lattice Boltzmann analysis of effect of heating location and Rayleigh number on natural convection in partially heated open ended cavity,Korean J.Chem.Eng.32(8)(2015)1498–1514.

    [7]K.M.Gangawane,R.P.Bharti,S.Kumar,Two-dimensional lattice Boltzmann simulation ofnatural convection in differentially heated square cavity:effect of Prandtl and Rayleigh numbers,Can.J.Chem.Eng.93(2015)766–780.

    [8]K.M.Gangawane,R.P.Bharti,S.Kumar,Effects of heating location and size on natural convection in partially heated open ended enclosure by using lattice Boltzmann method,Heat Transfer Eng.36(7)(2016)507–522.

    [9]J.R.Koseff,A.K.Prasad,The lid driven cavity flow:a synthesis of quantitative and qualitative observations,ASME J.Fluids Eng.106(1984)390–398.

    [10]M.Morzinski,C.O.Popiel,Linear heat transfer in a two-dimensional cavity covered by a moving wall,Numer.Heat Transfer A12(1988)265–273.

    [11]H.F.Oztop,Mixed convection in two-sided lid-driven differentially heated square cavity,Int.J.Heat Mass Transf.47(2004)1761–1769.

    [12]M.Jami,A.Mezrhab,M.Bouzidi,P.Lallemand,Lattice Boltzmann method applied to the laminar natural convection in an enclosure with heat cylinder conducting body,Int.J.Therm.Sci.46(2007)38–47.

    [13]D.Angeli,P.Levoni,G.S.Barozzi,Numerical predictions for stable buoyant regimes within a square cavity containing a heated horizontal cylinder,Int.J.Heat Mass Transf.51(2008)553–565.

    [14]G.Cesini,M.Paroncini,G.Cortella,M.Manzan,Natural convection from a horizontal cylinder in a rectangular cavity,Int.J.Heat Mass Transf.42(1999)1801–1811.

    [15]S.F.Dong,Y.T.Li,Conjugate of natural convection and conduction in a complicated enclosure,Int.J.Heat Mass Transf.24(2004)2233–2239.

    [16]J.M.House,C.Beckermann,T.F.Smith,Effect of a centered conducting body on natural convection heat transfer in an enclosure,Numer.Heat Transfer A18(1990)213–225.

    [17]J.Y.Oh,M.Y.Ha,Y.S.Kim,Numerical study of heat transfer and flow of natural convection in an enclosure with a heat-generating conducting body,Numer.HeatTransfer A31(1997)289–304.

    [18]M.Y.Ha,M.J.Jung,Y.S.Kim,A numerical study on transient heat transfer and fluid flow of natural convection in an enclosure with a heat-generating conducting body,Numer.Heat Transfer A35(1999)415–434.

    [19]M.Y.Ha,M.J.Jung,A numerical study on three-dimensional conjugate heat transfer of natural convection and conduction in a differentially heated cubic enclosure with a heat-generating cubic conducting body,Int.J.Heat Mass Transf.43(2000)4229–4248.

    [20]A.Mezrhab,H.Bouali,C.Abid,Modelling of combined radiative and convective heat transfer in an enclosure with a heat-generating conducting body,Int.J.Comput.Methods2(3)(2005)431–450.

    [21]X.Xu,G.Sun,Z.Yu,Y.Hu,L.Fan,K.Cen,Numerical investigation of laminar natural convective heat transfer from a horizontal triangular cylinder to its concentric cylindrical enclosure,Int.J.Heat Mass Transf.52(2009)3176–3186.

    [22]K.M.Gangawane,Convective Flow and Heat Transfer Analysis by Using Thermal Lattice Boltzmann Method,Ph.D.Thesis,IIT Roorkee,India,2015.

    [23]H.F.Oztop,Z.Zhao,B.Yu,Fluid flow due to combined convection in lid-driven enclosure having a circular body,Int.J.Heat Fluid Flow30(2009)886–901.

    [24]M.A.H.Mamun,M.M.Rahman,M.M.Billah,R.Saidur,A numerical study on the effect of a heated hollow cylinder on mixed convection in a ventilated cavity,Int.Commun.Heat Mass Transfer37(2010)1326–1334.

    [25]S.H.Hussain,A.K.Hussein,Numerical investigation of natural convection phenomena in a uniformly heated circular cylinder immersed in square enclosure filled with air at different vertical locations,Int.Commun.Heat Mass Transfer37(2010)1115–1126.

    [26]M.M.Rahman,M.A.Alim,M.M.A.Sarker,Numerical study on the conjugate effect of joule heating and magnato-hydrodynamics mixed convection in an obstructed lid-driven square cavity,Int.Commun.Heat Mass Transfer37(2010)524–534.

    [27]V.A.F.Costa,A.M.Raimundo,Steady mixed convection in a differentially heated square enclosure with an active rotating circular cylinder,Int.J.Heat Mass Transf.53(2010)1208–1219.

    [28]S.H.Hussain,A.K.Hussein,Mixed convection heat transfer in a differentially heated square enclosure with a conductive rotating circular cylinder at different vertical locations,Int.Commun.Heat Mass Transfer38(2011)263–274.

    [29]Y.G.Park,H.S.Yoon,M.Y.Ha,Natural convection in square enclosure with hot and cold cylinders at different vertical locations,Int.J.Heat Mass Transf.55(2012)7911–7925.

    [30]K.Khanafer,S.M.Aithal,Laminar mixed convection flow and heat transfer characteristics in a lid driven cavity with a circular cylinder,Int.J.Heat Mass Transf.66(2013)200–209.

    [31]C.Choi,S.Jeong,M.Y.Ha,H.S.Yoon,Effect of a circular cylinder's location on natural convection in a rhombus enclosure,Int.J.Heat Mass Transf.77(2014)60–73.

    [32]D.Chatterjee,S.K.Gupta,B.Mondal,Mixed convective transport in a lid-driven cavity containing a nano fluid and a rotating circular cylinder at the center,Int.Commun.Heat Mass Transfer56(2014)71–78.

    [33]H.Xu,R.Xiao,F.Karimi,M.Yang,Y.Zhang,Numerical study of double diffusive mixed convection around a heated cylinder in an enclosure,Int.J.Therm.Sci.78(2014)169–181.

    [34]M.Sheikholeslami,M.Gorji-Bandpy,I.Pop,S.Soleimani,Numerical study of natural convection between a circular enclosure and a sinusoidal cylinder using control volume based finite element method,Int.J.Therm.Sci.72(2013)147–158.

    [35]H.Bararnia,S.Soleimani,D.D.Ganji,Lattice Boltzmann simulation of natural convection around a horizontal elliptic cylinder inside a square enclosure,Int.Commun.Heat Mass Transfer38(2011)1436–1442.

    [36]K.M.Gangawane,Computational analysis of mixed convection heat transfer characteristics in lid-driven cavity containing triangular block with constant heat flux:Effect of Prandtl and Grashof numbers,Int.J.Heat Mass Transf.(2016),http://dx.doi.org/10.1016/j.ijheatmasstransfer.2016.09.061.

    [37]Z.T.Yu,L.W.Fan,Y.C.Hu,K.F.Cen,Prandtl number dependence of laminar natural convection heat transfer in a horizontal cylindrical enclosure with an inner coaxial triangular cylinder,Int.J.Heat Fluid Flow53(7–8)(2010)1333–1340.

    [38]A.A.Mehrizi,M.Farhadi,S.Shayamehr,Natural convection flow of Cu–water nano fluid in horizontal cylindrical annuli with inner triangular cylinder using lattice Boltzmann method,Int.Commun.Heat Mass Transfer44(2013)147–156.

    [39]M.Kalteh,K.Javaherdeh,T.Azarbarzin,Numerical solution of nanofluid mixed convection heat transfer in a lid-driven square cavity with triangular heat source,Powder Technol.253(2014)780–788.

    [40]A.W.Islam,M.A.R.Sharif,E.S.Carlson,Mixed convection in a lid driven square cavity with an isothermally heated square blockage inside,Int.J.Heat Mass Transf.55(2012)5244–5255.

    [41]M.Lamsaadi,M.Na?mi,M.Hasnaoui,Natural convection heat transfer in shallow horizontal rectangular enclosures uniformly heated from the side and filled with non-Newtonian power law fluids,Energy Convers.Manag.47(15–16)(2006)2535–2551.

    [42]C.Y.Cheng,Natural convection heat transfer of non-Newtonian fluids in porous media from a vertical cone under mixed thermal boundary conditions,Int.J.Heat Mass Transf.36(7)(2009)693–697.

    [43]R.R.Huilgol,G.H.R.Kefayati,Natural convection problem in a Bingham fluid using the operator-splitting method,J.Non-Newtonian Fluid Mech.220(2015)22–32.

    [44]A.Bose,N.Nirmalkar,R.P.Chhabra,Effect of aiding-buoyancy on mixed-convection from a heated cylinder in Bingham plastic fluids,J.Non-Newtonian Fluid Mech.220(2015)3–21.

    [45]O.Turan,N.Chakraborty,R.J.Poole,Laminar natural convection of Bingham fluids in a square enclosure with differentially heated side walls,J.Non-Newtonian Fluid Mech.165(2010)901–913.

    [46]O.Turan,R.J.Poole,N.Chakraborty,Aspect ratio effects in laminar natural convection of Bingham fluids in rectangular enclosures with differentially heated side walls,J.Non-Newtonian Fluid Mech.166(2011)208–230.

    [47]O.Turan,A.Sachdeva,N.Chakraborty,R.J.Poole,Laminar natural convection of power-law fluids in a square enclosure with differentially heated side walls subjected to constant temperatures,J.Non-Newtonian Fluid Mech.166(2011)1049–1063.

    [48]N.B.Khelifa,Z.Alloui,H.Beji,P.Vasseur,Natural convection in a horizontal porous cavity filled with a non-Newtonian binary fluid of power-law type,J.Non-Newtonian Fluid Mech.169–170(2012)15–25.

    [49]C.Sasmal,R.P.Chhabra,Laminar natural convection from a heated square cylinder immersed in power-law liquids,J.Non-Newtonian Fluid Mech.166(2011)811–830.

    [50]A.Chandra,R.P.Chhabra,Laminar free convection from a horizontal semi-circular cylinder to power-law fluids,Int.J.Heat Mass Transf.55(2012)2934–2944.

    [51]M.H.Matin,I.Pop,S.Khanchezar,Natural convection of power-law fluid between two-square eccentric duct annuli,J.Non-Newtonian Fluid Mech.197(2013)11–23.

    [52]M.Nazari,L.Louhghalam,M.H.Kayhani,Lattice Boltzmann simulation of double diffusive natural convection in a square cavity with a hot square obstacle,Chin.J.Chem.Eng.23(2015)22–30.

    [53]F.Selimefendigil,H.F.Oztop,Natural convection and entropy generation of nanofluid filled cavity having different shaped obstacles under the influence of magnetic field and internal heat generation,J.Taiwan Inst.Chem.Eng.56(2015)42–56.

    [54]Q.Ren,C.L.Chan,Natural convection with an array of solid obstacles in an enclosure by lattice Boltzmann method on a CUDA computation platform,Int.J.Heat Mass Transf.93(2016)273–285.

    [55]A.Azimifar,S.Payan,Enhancement of heat transfer of con fined enclosures with free convection using blocks with PSO algorithm,Appl.Therm.Eng.101(2016)79–91.

    [56]B.P.Leonard,Order of accuracy of QUICK and related convection-diffusion schemes,Appl.Math.Model.19(2016)640–653.

    [57]R.P.Bharti,R.P.Chhabra,V.Eswaran,Steady forced convection heat transfer from a heated circular cylinder to power-law fluids,Int.J.Heat Mass Transf.50(2007)977–990.

    [58]G.De VahlDavis,Natural convection of air in a square cavity:a bench mark numerical solution,Int.J.Numer.Methods Fluids3(3)(1983)249–264.

    av电影中文网址| 亚洲精品在线观看二区| 日韩免费av在线播放| 好男人电影高清在线观看| 午夜免费成人在线视频| 97人妻精品一区二区三区麻豆 | 香蕉久久夜色| 香蕉久久夜色| 免费在线观看视频国产中文字幕亚洲| 日韩精品中文字幕看吧| 老司机深夜福利视频在线观看| 成人18禁在线播放| 美女扒开内裤让男人捅视频| 欧美日韩亚洲国产一区二区在线观看| 久久精品国产99精品国产亚洲性色 | 国产色视频综合| 男人操女人黄网站| 国产亚洲欧美98| 欧美日韩亚洲国产一区二区在线观看| 一卡2卡三卡四卡精品乱码亚洲| 法律面前人人平等表现在哪些方面| 欧美中文日本在线观看视频| 狂野欧美激情性xxxx| 国产亚洲欧美98| 夜夜夜夜夜久久久久| 国产成人欧美在线观看| 三级毛片av免费| 亚洲欧美激情综合另类| 国内精品久久久久久久电影| 香蕉国产在线看| 操美女的视频在线观看| 级片在线观看| 一级毛片精品| 亚洲欧美日韩另类电影网站| 日韩欧美国产一区二区入口| 日韩 欧美 亚洲 中文字幕| 精品人妻在线不人妻| 免费在线观看完整版高清| 一二三四社区在线视频社区8| 日韩三级视频一区二区三区| 国产亚洲欧美在线一区二区| 久久久久久大精品| 性色av乱码一区二区三区2| 性色av乱码一区二区三区2| 日韩欧美免费精品| 久久中文字幕一级| 一个人观看的视频www高清免费观看 | 国产高清videossex| 多毛熟女@视频| 亚洲国产中文字幕在线视频| 淫妇啪啪啪对白视频| 91老司机精品| 精品高清国产在线一区| 国产精品亚洲美女久久久| 欧美日本视频| 岛国在线观看网站| 国产午夜福利久久久久久| 久久精品91蜜桃| 国产xxxxx性猛交| 久久久精品欧美日韩精品| 久久久国产成人免费| 久久午夜综合久久蜜桃| 亚洲熟女毛片儿| 精品人妻在线不人妻| 亚洲性夜色夜夜综合| 1024视频免费在线观看| 国产精品久久久久久精品电影 | 757午夜福利合集在线观看| 咕卡用的链子| 看免费av毛片| 亚洲av成人一区二区三| 亚洲午夜精品一区,二区,三区| 一边摸一边抽搐一进一出视频| netflix在线观看网站| 亚洲中文日韩欧美视频| 成人欧美大片| 真人做人爱边吃奶动态| 国产精品一区二区免费欧美| 精品日产1卡2卡| 国产精品一区二区在线不卡| 美女扒开内裤让男人捅视频| 我的亚洲天堂| 看免费av毛片| 国产成人免费无遮挡视频| 欧美日韩瑟瑟在线播放| 18禁观看日本| 国产成人精品在线电影| 两人在一起打扑克的视频| 精品久久久精品久久久| 国产主播在线观看一区二区| 欧美乱妇无乱码| 50天的宝宝边吃奶边哭怎么回事| 自线自在国产av| or卡值多少钱| 精品日产1卡2卡| 国产乱人伦免费视频| 男女做爰动态图高潮gif福利片 | 亚洲欧美精品综合一区二区三区| 精品国产美女av久久久久小说| 亚洲精品久久国产高清桃花| 日日爽夜夜爽网站| 国产精品久久久av美女十八| 在线观看午夜福利视频| 成人欧美大片| 美女午夜性视频免费| 免费人成视频x8x8入口观看| 久久国产精品人妻蜜桃| 女人被躁到高潮嗷嗷叫费观| 黄频高清免费视频| 一卡2卡三卡四卡精品乱码亚洲| av免费在线观看网站| 亚洲视频免费观看视频| 69精品国产乱码久久久| 熟妇人妻久久中文字幕3abv| 18禁国产床啪视频网站| 好男人在线观看高清免费视频 | 大型黄色视频在线免费观看| 精品免费久久久久久久清纯| 亚洲情色 制服丝袜| 亚洲精品中文字幕在线视频| 国产精品爽爽va在线观看网站 | 天堂影院成人在线观看| 国产欧美日韩一区二区三区在线| 亚洲中文av在线| 久久久国产成人免费| x7x7x7水蜜桃| 国产单亲对白刺激| 一个人免费在线观看的高清视频| 久久久久久久午夜电影| 国产精品av久久久久免费| 夜夜爽天天搞| 视频在线观看一区二区三区| 亚洲人成电影观看| 亚洲成a人片在线一区二区| 国产激情久久老熟女| 色综合欧美亚洲国产小说| 日日夜夜操网爽| 午夜福利影视在线免费观看| 欧美日韩瑟瑟在线播放| 天天躁夜夜躁狠狠躁躁| 精品国产美女av久久久久小说| 国产亚洲av嫩草精品影院| 窝窝影院91人妻| av天堂在线播放| av视频在线观看入口| 在线播放国产精品三级| 国产精品影院久久| 又大又爽又粗| 91在线观看av| 国产xxxxx性猛交| 日韩高清综合在线| 亚洲国产日韩欧美精品在线观看 | 免费女性裸体啪啪无遮挡网站| 亚洲精品国产一区二区精华液| 国产亚洲av高清不卡| 天堂√8在线中文| 国产精品秋霞免费鲁丝片| 一本大道久久a久久精品| 一区福利在线观看| 国产欧美日韩综合在线一区二区| 国产乱人伦免费视频| 国产精品一区二区三区四区久久 | 成熟少妇高潮喷水视频| 纯流量卡能插随身wifi吗| 免费观看人在逋| 黄色丝袜av网址大全| ponron亚洲| 免费女性裸体啪啪无遮挡网站| 在线播放国产精品三级| 国产亚洲精品第一综合不卡| 亚洲成av人片免费观看| 亚洲自偷自拍图片 自拍| 国产精品久久久人人做人人爽| 69精品国产乱码久久久| 国语自产精品视频在线第100页| 国产精品一区二区三区四区久久 | 日韩欧美在线二视频| 999久久久精品免费观看国产| 久久久久久大精品| 九色国产91popny在线| 亚洲专区中文字幕在线| 国产精品久久视频播放| 好看av亚洲va欧美ⅴa在| 久热这里只有精品99| 久久久国产成人免费| 免费在线观看视频国产中文字幕亚洲| 在线天堂中文资源库| 亚洲成国产人片在线观看| 多毛熟女@视频| 亚洲成人久久性| 欧美日韩亚洲综合一区二区三区_| 日本在线视频免费播放| 国产野战对白在线观看| 夜夜爽天天搞| 搡老熟女国产l中国老女人| 久久精品国产亚洲av香蕉五月| 后天国语完整版免费观看| 国产熟女xx| 久久婷婷人人爽人人干人人爱 | 国产精品综合久久久久久久免费 | 久久青草综合色| 午夜影院日韩av| 又大又爽又粗| 免费一级毛片在线播放高清视频 | 给我免费播放毛片高清在线观看| 高潮久久久久久久久久久不卡| 少妇 在线观看| 中文字幕人成人乱码亚洲影| 久久久久久久精品吃奶| 黄色丝袜av网址大全| 岛国在线观看网站| 国产蜜桃级精品一区二区三区| 久久精品国产亚洲av高清一级| 免费在线观看亚洲国产| 香蕉久久夜色| 精品国产一区二区三区四区第35| 啪啪无遮挡十八禁网站| 成人18禁高潮啪啪吃奶动态图| 国产精品国产高清国产av| 丁香欧美五月| av超薄肉色丝袜交足视频| 成人三级做爰电影| av免费在线观看网站| 日韩精品免费视频一区二区三区| 亚洲精品中文字幕一二三四区| 久久久国产成人精品二区| 精品久久久精品久久久| 一区福利在线观看| 大香蕉久久成人网| 亚洲午夜理论影院| 一二三四社区在线视频社区8| 老司机福利观看| 悠悠久久av| 亚洲aⅴ乱码一区二区在线播放 | 日日夜夜操网爽| av欧美777| 热99re8久久精品国产| 黑人巨大精品欧美一区二区mp4| 91成年电影在线观看| 国产高清videossex| 欧美日本亚洲视频在线播放| 操出白浆在线播放| 不卡一级毛片| 日韩高清综合在线| 国语自产精品视频在线第100页| or卡值多少钱| 国产精品野战在线观看| 亚洲熟妇中文字幕五十中出| 国产片内射在线| 亚洲成人精品中文字幕电影| 国产麻豆69| 国产亚洲av高清不卡| 午夜亚洲福利在线播放| 色综合婷婷激情| 在线免费观看的www视频| 久久精品国产99精品国产亚洲性色 | 在线免费观看的www视频| 狂野欧美激情性xxxx| 欧美一区二区精品小视频在线| 人成视频在线观看免费观看| 不卡av一区二区三区| 国产亚洲欧美在线一区二区| 啦啦啦观看免费观看视频高清 | 久久 成人 亚洲| 在线永久观看黄色视频| 亚洲av电影不卡..在线观看| 亚洲国产精品999在线| 一级作爱视频免费观看| 中文亚洲av片在线观看爽| 一夜夜www| 麻豆一二三区av精品| 精品乱码久久久久久99久播| 校园春色视频在线观看| 精品久久久久久,| 久久香蕉激情| 麻豆久久精品国产亚洲av| 色精品久久人妻99蜜桃| 激情在线观看视频在线高清| 日本免费a在线| 日日摸夜夜添夜夜添小说| 18禁国产床啪视频网站| 亚洲avbb在线观看| 国产真人三级小视频在线观看| 一进一出抽搐动态| 国产欧美日韩一区二区三区在线| 久久久国产成人免费| 又紧又爽又黄一区二区| 中亚洲国语对白在线视频| 欧美乱色亚洲激情| av在线天堂中文字幕| 国产激情欧美一区二区| 亚洲熟妇熟女久久| 黄色a级毛片大全视频| 久久狼人影院| 国产97色在线日韩免费| 在线观看66精品国产| 99国产精品免费福利视频| 可以在线观看的亚洲视频| 精品国产国语对白av| 国产精品久久久av美女十八| 国产色视频综合| 波多野结衣巨乳人妻| 亚洲自拍偷在线| 夜夜爽天天搞| 91精品国产国语对白视频| 人人妻,人人澡人人爽秒播| 一夜夜www| 亚洲第一电影网av| 伊人久久大香线蕉亚洲五| 国产精品自产拍在线观看55亚洲| 中文字幕色久视频| 亚洲国产高清在线一区二区三 | 国产一区在线观看成人免费| 国产色视频综合| 真人做人爱边吃奶动态| 久久亚洲精品不卡| 欧美日韩福利视频一区二区| 岛国在线观看网站| 亚洲熟妇熟女久久| 亚洲av电影不卡..在线观看| 一区二区三区高清视频在线| 国产高清有码在线观看视频 | 黑人欧美特级aaaaaa片| 亚洲色图av天堂| 窝窝影院91人妻| 脱女人内裤的视频| 亚洲国产看品久久| 午夜精品久久久久久毛片777| 亚洲一区二区三区不卡视频| 亚洲视频免费观看视频| 操美女的视频在线观看| 国产精品综合久久久久久久免费 | 美女高潮到喷水免费观看| 久久精品国产99精品国产亚洲性色 | 成人永久免费在线观看视频| 色综合亚洲欧美另类图片| 久久性视频一级片| 丝袜美足系列| 老司机靠b影院| 国产亚洲欧美在线一区二区| 97超级碰碰碰精品色视频在线观看| 精品国内亚洲2022精品成人| 亚洲精品中文字幕在线视频| 国产91精品成人一区二区三区| 久99久视频精品免费| 91在线观看av| 男女做爰动态图高潮gif福利片 | 精品一区二区三区av网在线观看| 999久久久精品免费观看国产| 精品乱码久久久久久99久播| 十八禁人妻一区二区| 亚洲色图av天堂| 欧美成人性av电影在线观看| 757午夜福利合集在线观看| 午夜亚洲福利在线播放| 欧美日韩瑟瑟在线播放| 国产精品久久久久久人妻精品电影| 国产精品乱码一区二三区的特点 | 在线观看66精品国产| 夜夜爽天天搞| 黑人巨大精品欧美一区二区mp4| 亚洲精华国产精华精| 欧美精品亚洲一区二区| 日韩高清综合在线| 妹子高潮喷水视频| 久久久国产成人免费| 久久婷婷成人综合色麻豆| 如日韩欧美国产精品一区二区三区| 777久久人妻少妇嫩草av网站| 日韩av在线大香蕉| 国产乱人伦免费视频| 亚洲黑人精品在线| 精品少妇一区二区三区视频日本电影| 一区二区三区精品91| 夜夜看夜夜爽夜夜摸| 亚洲五月婷婷丁香| 国产成人欧美| bbb黄色大片| 97碰自拍视频| 国产精品亚洲av一区麻豆| 亚洲成人精品中文字幕电影| 老司机福利观看| av在线播放免费不卡| 精品午夜福利视频在线观看一区| av在线播放免费不卡| 亚洲男人的天堂狠狠| 日本 av在线| 多毛熟女@视频| 国产免费av片在线观看野外av| 成人18禁在线播放| 岛国视频午夜一区免费看| 亚洲在线自拍视频| 成人免费观看视频高清| 在线观看舔阴道视频| 国产精品av久久久久免费| 欧美乱色亚洲激情| 免费搜索国产男女视频| 亚洲中文字幕日韩| 中文字幕人妻熟女乱码| 亚洲伊人色综图| 女人精品久久久久毛片| 啦啦啦免费观看视频1| 校园春色视频在线观看| 国产私拍福利视频在线观看| 免费在线观看日本一区| 淫秽高清视频在线观看| 在线观看舔阴道视频| 日本免费a在线| 久久久久国产精品人妻aⅴ院| 两个人视频免费观看高清| 激情在线观看视频在线高清| 日韩大码丰满熟妇| 国产精品综合久久久久久久免费 | 动漫黄色视频在线观看| 日韩有码中文字幕| 人人妻人人澡欧美一区二区 | 在线免费观看的www视频| 亚洲欧美日韩高清在线视频| 久久久国产精品麻豆| 9191精品国产免费久久| 亚洲一码二码三码区别大吗| 国产亚洲av高清不卡| 香蕉国产在线看| 成人av一区二区三区在线看| 黄色片一级片一级黄色片| 精品人妻1区二区| 国产成人欧美在线观看| 婷婷六月久久综合丁香| 欧美一级a爱片免费观看看 | 88av欧美| 色综合婷婷激情| 欧美日韩乱码在线| 国产精品一区二区精品视频观看| 青草久久国产| www.自偷自拍.com| 亚洲欧美精品综合久久99| 亚洲片人在线观看| 淫妇啪啪啪对白视频| 亚洲精品粉嫩美女一区| 国产精品亚洲美女久久久| 欧美在线一区亚洲| 亚洲精品在线美女| 国产精品秋霞免费鲁丝片| 精品久久久久久,| 国产激情久久老熟女| 国产精品,欧美在线| 夜夜躁狠狠躁天天躁| 久久久久久大精品| 欧美激情久久久久久爽电影 | 丝袜在线中文字幕| 国产亚洲欧美98| 精品国产一区二区三区四区第35| 国产日韩一区二区三区精品不卡| 桃红色精品国产亚洲av| 少妇粗大呻吟视频| 亚洲一区高清亚洲精品| 国产精品久久久人人做人人爽| 国产人伦9x9x在线观看| 国产黄a三级三级三级人| 黑人操中国人逼视频| 男女之事视频高清在线观看| 一边摸一边做爽爽视频免费| 国产麻豆成人av免费视频| 色在线成人网| 亚洲欧洲精品一区二区精品久久久| 国产欧美日韩精品亚洲av| 精品国产国语对白av| 国产成人精品久久二区二区91| 欧美成人午夜精品| 无人区码免费观看不卡| 美女午夜性视频免费| 亚洲天堂国产精品一区在线| 国产三级黄色录像| 女性被躁到高潮视频| 男人的好看免费观看在线视频 | 亚洲电影在线观看av| 黄色片一级片一级黄色片| 97人妻精品一区二区三区麻豆 | 亚洲中文日韩欧美视频| 男女下面插进去视频免费观看| 国产精品亚洲一级av第二区| 在线观看免费午夜福利视频| 国产成人系列免费观看| 久久热在线av| 国产成人精品在线电影| 97人妻天天添夜夜摸| 亚洲色图 男人天堂 中文字幕| 97碰自拍视频| 亚洲avbb在线观看| 欧美日本中文国产一区发布| 日韩国内少妇激情av| 午夜福利在线观看吧| 每晚都被弄得嗷嗷叫到高潮| 老司机午夜福利在线观看视频| 宅男免费午夜| 精品高清国产在线一区| 夜夜爽天天搞| 日本免费一区二区三区高清不卡 | 国产91精品成人一区二区三区| 麻豆一二三区av精品| 大型av网站在线播放| 久久久久久免费高清国产稀缺| 非洲黑人性xxxx精品又粗又长| 男人舔女人下体高潮全视频| 人妻丰满熟妇av一区二区三区| 天天添夜夜摸| 黄色成人免费大全| 麻豆国产av国片精品| 日韩成人在线观看一区二区三区| 91麻豆av在线| 亚洲精品久久成人aⅴ小说| 国产精品久久电影中文字幕| 日韩精品免费视频一区二区三区| 成人亚洲精品av一区二区| 国产不卡一卡二| 国产成人啪精品午夜网站| 天堂√8在线中文| 国产亚洲精品综合一区在线观看 | 国产精品爽爽va在线观看网站 | 亚洲欧美激情在线| 日韩一卡2卡3卡4卡2021年| 久久婷婷人人爽人人干人人爱 | 成在线人永久免费视频| 一级a爱视频在线免费观看| 在线观看免费午夜福利视频| 99re在线观看精品视频| 夜夜躁狠狠躁天天躁| 午夜亚洲福利在线播放| 国产一级毛片七仙女欲春2 | 日韩国内少妇激情av| 黄片大片在线免费观看| 自线自在国产av| 日本三级黄在线观看| 变态另类丝袜制服| 欧美日韩瑟瑟在线播放| 中文字幕人妻熟女乱码| 国产xxxxx性猛交| 久久久久久久久免费视频了| 亚洲av成人av| 午夜日韩欧美国产| 免费看美女性在线毛片视频| 女性被躁到高潮视频| 国产精品久久久久久人妻精品电影| av免费在线观看网站| 这个男人来自地球电影免费观看| 麻豆成人av在线观看| 两性夫妻黄色片| 亚洲国产欧美日韩在线播放| 人人妻人人澡人人看| 波多野结衣高清无吗| 成人国语在线视频| 一a级毛片在线观看| 在线播放国产精品三级| 免费在线观看完整版高清| 久久青草综合色| 久久久久久国产a免费观看| 丝袜美足系列| 午夜亚洲福利在线播放| 97人妻天天添夜夜摸| 这个男人来自地球电影免费观看| 色哟哟哟哟哟哟| 在线观看免费午夜福利视频| 村上凉子中文字幕在线| 自拍欧美九色日韩亚洲蝌蚪91| 国产精品久久久人人做人人爽| 久久热在线av| 亚洲无线在线观看| 欧美 亚洲 国产 日韩一| 国产精品 国内视频| 妹子高潮喷水视频| 黄色片一级片一级黄色片| 欧美av亚洲av综合av国产av| 欧美日韩亚洲国产一区二区在线观看| 国产一卡二卡三卡精品| 怎么达到女性高潮| 成在线人永久免费视频| 又黄又爽又免费观看的视频| a在线观看视频网站| 黄色 视频免费看| e午夜精品久久久久久久| 久久精品国产亚洲av香蕉五月| 久久精品国产综合久久久| 一本久久中文字幕| 色在线成人网| 国产激情欧美一区二区| 一个人观看的视频www高清免费观看 | 午夜影院日韩av| 欧美日韩亚洲国产一区二区在线观看| 一区福利在线观看| 国产精品一区二区免费欧美| 亚洲中文av在线| 老司机深夜福利视频在线观看| 99国产精品一区二区蜜桃av| 老熟妇乱子伦视频在线观看| 又大又爽又粗| 黄色 视频免费看| 99久久国产精品久久久| 丁香欧美五月| 欧美国产日韩亚洲一区| 中文字幕人妻丝袜一区二区| 成人精品一区二区免费| 少妇裸体淫交视频免费看高清 | 久久中文字幕人妻熟女| 天天一区二区日本电影三级 | 亚洲五月色婷婷综合| 两性夫妻黄色片| 嫩草影视91久久| 熟女少妇亚洲综合色aaa.| 午夜福利一区二区在线看| 中文字幕久久专区| 99国产综合亚洲精品| 欧美精品啪啪一区二区三区| 亚洲成av人片免费观看| 97碰自拍视频| 免费久久久久久久精品成人欧美视频|