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

    Numerical investigation of entropy generation and heat transfer of pulsating flow in a horizontal channel with an open cavity*

    2017-09-15 13:55:51FatmaZamzariZouhaierMehrezAfifElCafsiAliBelghithPatrickLeQu
    水動力學研究與進展 B輯 2017年4期

    Fatma Zamzari, Zouhaier Mehrez, Afif El Cafsi, Ali Belghith, Patrick Le Quéré

    1.Faculty of Sciences of Tunis, Laboratory of Energy, Heat and Mass Transfer (LETTM), Department of Physics, El Manar University, El Manar 2092, Tunisia, E-mail: fatmazamzari@live.fr

    2.Computer Science Laboratory for Mechanics and Engineering Sciences (LIMSI), CNRS Bat 508, BP 133, 91403 Orsay Cedex, France

    Numerical investigation of entropy generation and heat transfer of pulsating flow in a horizontal channel with an open cavity*

    Fatma Zamzari1, Zouhaier Mehrez1, Afif El Cafsi1, Ali Belghith1, Patrick Le Quéré2

    1.Faculty of Sciences of Tunis, Laboratory of Energy, Heat and Mass Transfer (LETTM), Department of Physics, El Manar University, El Manar 2092, Tunisia, E-mail: fatmazamzari@live.fr

    2.Computer Science Laboratory for Mechanics and Engineering Sciences (LIMSI), CNRS Bat 508, BP 133, 91403 Orsay Cedex, France

    In this study, the entropy generation and the heat transfer of pulsating air flow in a horizontal channel with an open cavity heated from below with uniform temperature distribution are numerically investigated. A numerical method based on finite volume method is used to discretize the governing equations. At the inlet of the channel, pulsating velocity is imposed for a range of Strouhal numbersStpfrom 0 to 1 and amplitudeApfrom 0 to 0.5. The effects of the governing parameters, such as frequency and amplitude of the pulsation, Richardson number, Ri, and aspect ratio of the cavity,L/H, on the flow field, temperature distribution, average Nusselt number and average entropy generation, are numerically analyzed. The results indicate that the heat transfer and entropy generation are strongly affected by the frequency and amplitude of the pulsation and this depends on the Richardson number and aspect ratio of the cavity. The pulsation is more effective with the aspect ratio of the cavityin terms of heat transfer enhancement and entropy generation minimization.

    Pulsating flow, entropy generation, mixed convection heat transfer, open cavity

    Introduction

    Heat transfer in open cavities has received considerable attention, due to its presence in various engineering applications such as cooling of electronic devices, electronic chips and solar receivers. Thus, the knowledge of this flow type is crucial from an economic viewpoint or in terms of environmental acceptability. For this reason, a large number of papers have been recently published on this specific subject. In this context, Manca et al.[1]studied numerically the mixed convection in an open cavity heated from different sides with constant flux. They found that the opposing forced flow configuration has the highest thermal performance in terms of heat transfer rate. In Ref.[2], the same authors carried out an experimental investigationon mixed convection in an open cavity heated from the left wall. They showed that the maximum dimensional temperature rise values decrease as the Reynolds and Richardson numbers decrease. Leong et al.[3]analyzed the mixed convection in an open cavity in a horizontal channel. They found that the heat transfer rate was reduced and the flow becomes unstable in the mixed convection regime. Aminossadati and Ghasemi[4]presented a numerical study of mixed convection in a horizontal channel with a discrete heat source in an open cavity. They showed that the best cooling performance is given when the heat source is located on the right wall. Rahman et al.[5]studied the conjugated effect of joule heating and magneto-hydrodynamic on double-diffusive mixed convection in a horizontal channel with an open cavity. They showed that the flow behavior, heat and mass transfers are affected by applying external magnetic field. Mehrez et al.[6]performed a numerical investigation for the effect of non-isothermal heater on nanofluid flow in an open cavity. They showed thatthe flow behavior, heat transfer and entropy generation are strongly affected by the heater position. Abdelmassih et al.[7]conducted a three-dimensional numerical simulation for mixed convection flow inside a cubic open cavity. They showed that the flow becomes unstable and complicated for a higher Reynolds and Richardson numbers. Stiriba et al.[8]investigated heat transfer and fluid flow characteristics of laminar flow past a three-dimensional open cavity heated from below. They showed that the flow becomes steady and the heat diffusion is predominant at low Reynolds and Richardson numbers.

    To improve heat transfer, several techniques have been introduced which can be classified into two methods: active method (which requires external power source) and passive method (requires no external power source). Recently, the pulsating flow (active technique of control) received a great deal of research interest because it is frequently encountered in large application areas: biological systems (arterial flows in the human body), engineering systems (pulsed heat exchanger, cooling system of nuclear reactor, electronic cooling, pulse combustor for both civil and military uses, etc.). A numerical study of mixed convection of pulsating flow and heat transfer characteristics over a backward-facing step under laminar regime was numerically investigated by Khanafer et al.[9]. They showed that the average Nusselt number increases by increasing the pulsation frequency. Khanafer et al.[10]performed a numerical simulation of unsteady mixed convection in a driven cavity using an externally exited sliding lid. They showed that the average Nusselt number increases by increasing the lid frequency. Velazquez et al.[11]investigated a two-dimensional laminar pulsating flow in a heated rectangular cavity with different aspect ratios. They showed that the pulsation enhances heat transfer in the cavity. Sourtiji et al.[12]have numerically investigated the oscillating fluid flow and heat transfer characteristics in a square cavity with inlet and outlet ports. They showed that the heat transfer is enhanced for all the investigated Strouhal numbers in comparison to its steady state case. Selimefendigil and Oztop[13]have numerically investigated the forced convection in pulsating nanofluid flow over a backward facing step with a stationary circular cylin- der. They showed that the heat transfer is enhanced by increasing the frequency of the oscillation. Unal at al.[14]performed a numerical study of heat transfer characteristics and nanofluid flow in a wavy minichannel under pulsating inlet flow conditions. They showed that the heat transfer performance is enhanced at high nanoparticle volume fraction and at low pulsating frequency. In the case of open cavity geometry, Selimefendigil and Yurddas[15]performed a numerical study of pulsating flow in a channel with a cavity heated from below and a left vertical side. They showed that the heat transfer enhancement is affected by the pulsation and this depends on the Reynolds and Richardson numbers.

    The second law of thermodynamics is considerably important because it presents the basis of assessing the irreversibility (entropy generation). The minimization of entropy generation has been considered as an effective tool to improve the performance of heat transfer process and consequently enhances the efficiency of energy utilization. Entropy generation has recently been the topic of great interest in fields such as heat exchangers, turbo machinery, electronic cooling, porous media and combustion. The entropy generation in the convective flows is caused by the variations of irreversibilities due to the heat transfer (due to the conduction) and friction (due to the viscosity). Many studies have been published on entropy generation. Abu-Nada[16]investigated the entropy generation of laminar flow over a backward facing step. He showed that the entropy generation varies by variation of Reynolds number, Brinkman number and dimensionless temperature parameter. In the same context, Abu-Nada[17]showed that the entropy generation depends on the expansion ratio. Khorasanizadeh et al.[18]investigated the mixed convection and entropy generation of Cu-water nanofluid and pure water in a lid-driven square cavity. They showed that it is possible to obtain a configuration giving maximum heat transfer against minimum entropy production. Nayak et al.[19]have numerically studied the mixed convection of copper-water nanofluid inside a differentially heated skew enclosure. They showed that the average rate of heat transfer and entropy generation is sensitive to the skew angle and the nanoparticle volume fraction. Recently, a critical review of contributions to the theory and application of entropy generation analysis to different types of engineering systems is performed by Sciacovelli et al.[20].

    In the case of open cavity geometry, Mehrez et al.[21]numerically studied the heat transfer and entropy generation of nanofluid flow in an open cavity for various types of nanoparticles (Cu, CuO, Al2O3, TiO2). They showed that the heat transfer and the entropy generation rates are affected by varying the volume fraction and type of nanoparticles, Reynolds and Richardson numbers and aspect ratio of the cavity. For the same configuration, Mehrez et al.[22]analyzed the MHD effects on the flow structure, heat transfer and entropy generation of Cu-water nanofluid flow. They showed that the entropy generation and the heat transfer are strongly affected by the presence of external magnetic field. Recentely in Ref.[23], Mehrez et al. studied the effect of inclination angle of the channel on nanofluid flow. They showed that the heat transfer phenomenon and irreversibility mechanism depend mainly on the inclination angle. Zamzari etal.[24]studied the entropy generation and the heat transfer of air-flow in a horizontal channel with an open cavity, they showed that, with an appropriate choice of pertinent parameters of the flow and aspect ratio, it can obtain minimum entropy generation with high heat transfer.

    To our knowledge, the effect of pulsation on the flow behavior, heat transfer and entropy generation in an open cavity is not yet studied; this motivated us to study this problem. Therefore, the main objective of the current investigation is to understand how heat transfer and entropy generation are influenced by the presence of pulsation at the inlet. The results are graphically presented by streamlines, isotherms, spatial and spatial-temporal averaged Nusselt numbers, entropy generation and Bejan numbers for different Richardson numbers, aspect ratios of the cavity and various frequencies, phases and amplitudes of the pulsation.

    1. Physical problem

    A schematic description of the physical problem is shown in Fig.1. It consists of a two-dimensional horizontal channel including an open cavity.LandHare the length and height of the cavity respectively. The channel lengthlis set to 7H. The bottom wall of the cavity is heated by a uniform temperature distribution while the other walls are insulated. The Richardson number is varied in the range of 0.25≤Ri≤1. Three aspect ratios of the cavity,L/H=1, 1.5 and 2 are considered in this study. At the inlet the cold fluid is maintained at a uniform temperatureTc. The bottom wall of the cavity is heated with a constant temperatureTH. A sinusoidal velocity is imposed in the entrance of the cavity:wherefpandApare the frequency and the amplitude of the pulsation respectively.

    Fig.1 Geometric configuration

    2. Governing equations

    A two-dimensional, incompressible and laminar flow is considered. The fluid properties are assumed to be constant. The dimensional form equations, continuity (1), momentum (2), (3) and energy (4), describing the flow under Boussinesq approximation are as follows:

    These equations are used in previous published works[3,4,7].

    The following dimensionless parameters are used to obtain the dimensionless form of the above equations:

    The governing equations, in dimensionless form, are expressed as follows:

    The relevant physical parameters are:

    Richardson number

    whereGris the Grashof number,

    Reynolds number

    This expression ofRe, based on the cavity height and velocity at the inlet, is already used by previous works: Manca et al.[1], Abu-Nada et al.[16,17]and Mehrez et al.[21].

    Prandtl number

    in all the work has a fixed value of 0.71 correspond to the air.

    The boundary condition in dimensionless form can be expressed as:

    At the channel inlet: A velocity with a sinusoidal time dependent part is imposedU=[1+Asin(2πStp·ta)], whereAis the non-dimensional amplitude,

    is the Strouhal number based on the pulsation frequency

    On the bottom wall of the cavity, temperature is constantθ=1 and the remaining walls are insulated. (10)

    On the solid surfaces, a No-slip conditions are prescribed:U=V=0 (11)

    At the exit, convective boundary conditions are applied

    Heat exchanges for fluid flow within the cavity between the hot wall and the fluid are characterized by local Nusselt number given by the relation

    The spatial average Nusselt number,Nus, is obtained by integrating the above local Nusselt number

    The spatial and temporal average Nusselt number,Num, is obtained by integratingNusfor one period of the pulsation,τ, and is given by

    The normalized Nusselt numberis defined as

    whereis the non pulsating case.gives the rate of heat transfer variation by applying pulsation.

    Entropygeneration rate for a two-dimensional flow in Cartesian coordinates is given by Bejan[25]

    While dividingSgenby

    and by using the same dimensionless parameters given above, the dimensionless entropy generation is given as follows:

    where the irreversibility factorφis given as

    φcan be expressed as a function of non-dimensional numbers as follows

    whereis the Eckert number,is the Prandtl number,is the dimensionless temperature difference and

    The first term on the right-hand side of Eq.(19) corresponds to the heat transfer irreversibility (Sgen,c) whereas the second term corresponds to the irreversibility due to viscous dissipation (Sgen,v).

    The spatial average entropy generation over the entire flow domain is calculated by using the following equation

    where

    andθis the volume of computational domain.

    The spatial-temporal average entropy generation is obtained by integratingNssduring one period of pulsation

    The normalized entropy generationis defined as

    whereSt=0is the non pulsating case.gives the rate of entropy generation variation.

    The Bejan number,Be, defined as the ratio between the entropy generation due to heat transfer by the total entropy generation, is expressed as

    It is known that whenBe?0.5 the irreversibility due to heat transfer dominates. WhenBe?0.5 the irreversibility due to the viscous effects dominates the processes and ifBe=0.5 the entropy generation due to the viscous effects and the heat transfer effects are equal.

    3. Numerical procedure

    3.1 Numerical method

    The unsteady Navier-Stokes equations are discretized in space by the finite volume method. The time integration is performed using a time-splitting algorithm, also known as a prediction projection algorithm, which allows one to decouple pressure from velocity[26,27]. Assuming all quantities known at timenΔT, the solution at time (n+1)ΔTis obtained as follows:

    An intermediate velocity fieldV*is first computed using a second-order time scheme. This time stepping combines a second order backward Euler scheme for the diffusion terms, with an explicit second-order Adams-Bashforth extrapolation for the non-linear terms, taking into account known pressure field. This step reads

    In the second step, this intermediate velocity field is projected on to the subspace of divergence free vector field using the Helmholtz decomposition theorem. This step reads

    where

    It is accomplished by taking the divergence of equation giving rise to a Poisson’s type equation for the incremental pressure

    This equation is solved with a multigrid algorithm, in which the presence of internal blockings is automatically taken into account.

    Once the pressure field is obtained, the newquantities atn+1 are given by

    A second order upwind finite difference method by means of a QUICK scheme, proposed by Leonard[28]is used for the convective terms and a centered scheme for the diffusive ones.

    The computations are carried out by using a FORTRAN house code developed in collaboration between LIMSI (France) and LETTM (Tunisia). This code has been successfully used for different types of problems[26,27,29,30].

    3.2 Grid testing

    The effect of grid size is examined by running different combinations ofRi,L/H,ApandStpon different uniform grids. (Figs.2, 3) represent a sample of results giving the effect of grid refinement on the time and space averaged Nusselt and entropy generation, numbers forRe=200,Ri=1,L/H=1,Ap=0.1 andStp=0.4 (Fig.2) and local Nusselt numberNuforRe=200,Ri=1,L/H=1 andAp=0 (Fig.3). The grid of 132×124 is chosen for the present investigation because it makes a good compromise between the accuracy of the solution and the cost of computation time. It is noted that in the present wok, the aspect ratio of the cavityL/His varied by varying the length of the cavityLwhile fixing the heightH. For higher aspect ratios of the cavity, the number of grid points in theXdirection is increased accordingly.

    3.3 Solution verification

    It is known that the main source of numerical errors is the discretization due to limited number of grid points. To evaluate the numerical errors and associated uncertainties, the solution verification is performed following the method of safety factor proposed by Xing et al.[31-33]with three different grids: fine (132×124), medium (110×103) and coarse (92×86) forRe=200,Ri=1,L/H=1,Ap=0.1 andStp=0.4 using two parameters: time and space averaged Nusselt numberNumand time and space averaged entropy generation numberNsm(Table 1).The three chosen grids are generated using a constant grid refinement ratior=1.2 in the two spatial directions (XandY). If 3, 2, and 1 represent the coarse, medium, and fine grids with spacing Δx3, Δx2and Δx1,ris given as follows

    Fig.2 Variation of time- spatial averaged Nusselt and entropy generation numbers with the grid size

    Fig.3 Local Nusselt number profiles for various grid sizes

    Table 1 Solution verification for time and space averaged Nusselt and entropy generation numbers

    If the solutions for the fine, medium and coarse grids areS1,S2andS3, respectively, solution changesεfor medium-fine and coarse-medium solutions and the convergence ratioRare defined by:

    The monotonic convergence is achieved when 0<R<1. Then the three grid solutions can be used to compute the estimated order of accuracypRE, error

    When solutions are in the asymptotic range, i.e.,however, in many circumstances, especially for industrial applications, solutions are far from the asymptotic range such thatpREis greater or smaller thanpth[34]. In this work all the solutions are calculated in the asymptotic range. The ratio ofpREtopthis used as the distance metric

    The grid uncertainty is estimated by:

    As shown in Table 1, monotonic convergence is achieved with a low grid uncertainty of 0.38%S1forNumand 1.28%S1forNsm. This shows that the current fine grid resolution is sufficient for the present problem. It is used for all simulations.

    3.4 Code validation

    In order to testing the performance of the present numerical code, the steady-state solution represented by maximum temperatureθmaxand dimensionless minimum stream functionΨmin, which is obtained in a two-dimensional channel with a cavity heated from different sides with uniform heat flux atRe=100,Pr=0.71, is compared with the results of Manca et al.[1](see Table 2). Also, Fig.4 compares the temperature profile at midsection of the enclosure for assisting forced flow case given by our code and that of Manca et al.[1], =100ReandRi=0.1. In Table 2 and Fig.4 a good agreement between the present results and those of Ref.[1] is observed.

    Table 2 A comparison between the present study and Manca et al.[1]

    Fig.4 Comparison of temperature profile at midsection between the present study and Manca et al.[1]

    4. Results and discussion

    In the present work, investigations were performed within the range of Strouhal number:, amplitude of pulsation:Richardson number 0.25≤Ri≤1 and aspect ratio of the cavity 1≤L/H≤2 and for Reynolds numberRe=200, the Prandtl number is kept constant at 0.71 to guarantee constant fluid physical properties for mo-derate and small values of temperature difference. The irreversibility factorφis set equal to 0.01. To analyze the flow behavior and the effect of the mixed convection, the streamlines and the isotherms during a periodic-steady cycle for three values of aspect ratios:L/H=1,L/H=1.5 andL/H=2 are illustrated in Fig.5 forAp=0.1,Stp=0.3,Ri=0.5.

    Fig.5 On the top: Streamlines, on the bottom: Isotherms at different pulsation phase and variousL/H,Ap=0.1,Stp=0.3 and

    It is noted that, for all the consideredL/H, the presence of a large recirculation cell which fills the entire cavity. The intensity of this eddy varies during a pulsation period due to the temporal variation of the velocity on the entrance. It is also noted that the intensification of the flow field by increasingL/Hduring the whole period except tota=τ/2 when the flow forL/H=1.5 is slightly more intense thanL/H=2.

    Isotherm plots show the formation of thermal flow forL/H=1.5 is slightly more intense than layer near the heated surface (thermal stratification), where the temperature gradually varied from the heated wall to the cold part indicating the domination of heat transfer conduction mechanism. It is clearly seen that the thickness of the thermal layer near the heated surface and the clustering of the isotherms are changed during a period. It is noted also that, by increasingL/Hthe isotherms become closer near the hot wall showing the increase of temperature gradients.

    Fig.6 Variations of spatial-averaged Nusselt number versus time for different amplitudes of pulsation atRi=0.5

    Temporal variations of spatial-averaged Nusselt numberNus(Eq.(14)) for different pulsation amplitudesAp, different Strouhal numbersStpand various aspect ratiosL/Hare depicted in Fig.6. As can be seen in this Figure a sinusoidal shape of theNuswith time can be observed. The oscillation ofNusaround a mean value is obtained with a time period equal to that of the pulsation. The maximum value ofNusincreases with increasing the pulsation amplitude. The shape of the oscillation is preserved for low frequency, but by increasingStpthere appears distortion of the pure sinusoids.

    In order to evaluate the effects of the pulsation flow on the heat transfer enhancement, Fig.7 displays the time-spatial averaged Nusselt numberNum(Eq.(15)) versus Strouhal number, for different Richardson numbers (Ri=0.25,Ri=0.5 andRi= 1), for various aspect ratios of the cavity atAp=0.1. It is clearly seen that the all curves have the same increase-decrease profile with a peak at a certainStp.

    This phenomenon clearly illustrates the existence of an optimum frequency for the maximum heat transfer enhancement ratio. Starting fromStp=0 (steady flow)Numincreases by increasing Strouhal number until an optimalStpheat transfer rate is reached, and then decreases afterward. This enhancement of heat transfer rate is more effective by increasingRi. Thus, the heat transfer into the cavity is enhanced by increasing the buoyancy effects. It isnoted also that the location of the maximum heat transfer shifts towards the right (low value ofStp) by increasing the aspect ratio of the cavity. The maximum enhancement of heat transfer is obtained atStp=0.4 forL/H=1, atStp=0.4 forL/H=1.5 and forStpbetween 0.2 and 0.3 forL/H=2. This figure shows also that the flow with the aspect ratio of the cavityL/H=1.5 gives higher values ofNumcompared with the other cases.

    Fig.7 Time- spatial averaged Nusselt number against pulsation frequency for various Richardson numbers forAp=0.1

    The time-spatial averaged Nusselt numberNumagainst pulsation amplitudeApis plotted in Fig.8 for various Richardson numbers, and various aspect ratios (Stp=0.3). This figure shows the effect of increasing Richardson number on improvement of heat transfer for allApandL/Hvalues, and this is due to the intensification of buoyancy effect.

    Fig.8 Time-spatial averaged Nusselt number against pulsation amplitude for various Richardson numbers forStp=0.3

    Furthermore it is observed that, for all Richardson numbers, the time-spatial averaged Nusselt number increases sharply with the velocity amplitudes and this can be explained by the increase of the temperature gradient due to the increase ofAp. This shows the importance of the pulsation to improve heat transfer. It is noted also that the highest heat trans-fer is obtained withL/H=1.5.

    Fig.9 Variations of spatial-averaged entropy generation versus time for different amplitudes of pulsation atRi=0.5

    Temporal variations of spatial-averaged entropy generationNss(Eq.(21)) for different velocity amplitudes, different Strouhal numbers and various aspect ratios are depicted in Fig.9.Nssoscillates periodically around a mean value with a time period equal to that of pulsation. The sinusoidal aspect ofNssappears by increasing the Strouhal number. In all cases, the peak of the entropy generation increases generally with increasing the amplitude of the pulsation periodically around a mean value with a time period equal to that of pulsation. The sinusoidal aspect ofNssis more visible for low value ofStpand a distortion appears by increasing the Strouhal numfor various Richardson numbers atis displayed in Fig.10. Obviously, the figure shows a resonant type behavior similar to Fig.11.Nsmincreases by increasing pulsation frequency, attains a maximum and then decreases. This diagram not only verifies the effect ofStpon average entropy generation but also shows that by increasing Richardson numberNsmincreases due to augmentation of the ber. In all cases, the peak of the entropy generation increases generally with increasing the amplitude of the pulsation.

    To testing the influence of the pulsation on entropy generation, the time-spatial averaged entropy generationNsm(Eq.(23)) versus Strouhal numberbuoyancy effect. This increases the velocity and temperature gradients and therefore leads to an increase in viscous and conduction contribution to the entropy generation.

    Fig.10 Time- spatial averaged entropy generation against pulsation frequency for various Richardson number for

    These effects are more pronounced atL/H= 1.5 forRi=0.25 andRi=0.5 which gives a maximum entropy generation, but forRi=1 entropy generation becomes minimal. It is noted also that the maximum ofNsmis shifted toward the right by increasing the aspect ratio of the cavity. Thus this maximum is obtained atStp=0.4 forL/H=1,Stp=0.4 forL/H=1.5 andStpbetween 0.2 and 0.3 forL/H=2.

    Fig.11 Time- spatial averaged entropy generation against pulsation amplitude for various Richardson numbers for

    Figure 12 illustrates the variations of Bejan numberBe(Eq.(25)) versus Strouhal number for various aspect ratios atAp=0.1. It is seen that, by comparing the curves of different cases, the Bejan number increases by increasing aspect ratio of the cavity and this can be explained by the intensification of heat transfer irreversibility relative to the friction one. Results show also that the increase ofStpreducesBecompared with non-pulsed case. This is due to the intensification of viscous irreversibility compared with that of conduction.

    The Bejan number against pulsating amplitude for various aspect ratios, various Richardson numbersatStp=0.3 is plotted in Fig.13. It is observed, for all the cases, that the Bejan number decreases with the increment of the pulsation amplitudeAp, except atL/H=1.5 forRi=0.25 andRi=0.5 where Bejan number increases slightly with increasingApand then decreases. It can be observed from this figure that for low value ofAp,Beincreases by increasingL/H, but from a certain value of pulsation amplitude (Ap=0.25forRi=0.25,Ap=0.35 forRi=0.5 andAp=0.5 forRi=1), the values ofBeforL/H=1.5 are greater than those forL/H=2.

    Fig.12 Bejan number against Strouhal number for various aspect ratios atAp=0.1

    Fig.13 Bejan number against pulsation amplitude for various aspect ratios atStp=0.3

    It is noted that the values of Bejan number are high (close to 0.9). This proves that entropy generation contributed by conduction,Nsm,c, is dominant relative to that contributed by viscous irreversibility,Nsm,yin the entropy generation.

    Fig.14against Strouhal number for various Richardson numbers at

    To complete our investigation it is very important to see the effectiveness of the pulsation. So we choose to study the variations of the quotient between the normalized Nusselt number(Eq.(16)) and the normalized entropy generation(Eq.(24)) versus Strouhal number for various Richardson number atAp=1 (Fig.14). As shown in this figure, the quotientdecreases as a function ofStpattains a minimum and then increases and this for a value of aspect ratioL/H=1 andL/H=2. However forL/H=1.5, the reverse trend was observed,increases as a function ofStpattains a maximum forStp=0.4 and then decreases. In this case, the maximum of this quotient becomes more apparent by increasingRi. By comparing the quotient with the unit, it can be concluded that for configurations withL/H=1 andwhich means that the rate of entropy generation is more important than the rate of heat transfer enhancement. But for the case ofL/H=1.5, this ratio is greater than 1. This shows that the rate of heat transfer enhancement is greater than the rate of entropy generation. It should be noted that the mentioned variations ofare generally small, since the maximum variation is about 4 %.

    Finally, it can be concluded that, for all cases studied the pulsation is relatively more effective forL/H=1.5,Ri=1 and for a pulsation frequencySt=0.4 with a pulsation amplitudeAp=0.1.

    5. Conclusions

    In this paper, heat transfer and entropy generation of pulsating flow in an open cavity heated from the bottom wall are numerically studied. The effects of the pulsation parameters such as phase, frequency and amplitude as well as pertinent parameters such as Richardson number and the aspect ratio of the cavity are studied. The main conclusions obtained are summarized as follows:

    (1) The flow behavior and the temperature distribution are strongly affected by the pulsation.

    (2) The space-averaged Nusselt and entropy generation numbers oscillate periodically versus the time with a period equal to that of pulsation.

    (3) The pulsation increases heat transfer and entropy generation. The maximum heat transfer rate is observed for an optimum Strouhal number which depends on Richardson number and aspect ratio of the cavity.

    (4) The heat transfer and the entropy generation rates during a pulsation cycle increase considerably by increasing the pulsation amplitude.

    (5) The contribution rate of each irreversibility type in the entropy generation varies by varying the frequency and the amplitude of the pulsation. It depends on Richardson number and aspect ratio of the cavity.

    (6) The configuration with aspect ratio of the cavityL/H=1.5 is more effective because it gives, by applying a pulsation at the inlet, an enhancement rate of heat transfer greater than the increase rate of the entropy generation.

    [1] Manca O., Nardini S., Khanafer K. Effect of heatedwall position on mixed convection in a channel with an open cavity [J].Numerical Heat Transfer, 2003, 43(3): 259-282.

    [2] Manca O., Nardini S., VafaI K. Experimental investigation of mixed convection in a channel with an open cavity [J].Experimental Heat Transfer, 2006, 19(1): 53-68.

    [3] Leong J. C., Brown N. M., Lai F. C. Mixed convection from an open cavity in a horizontal channel [J].Inter-national Communications Heat Mass Transfer, 2005, 32(5): 583-592.

    [4] Aminossadati S. M., Ghasemi B. A numerical study of mixed convection in a horizontal channel with a discrete heat source in an open cavity [J].European Journal of Mechanics B/Fluids, 2009, 28(4): 590-598.

    [5] Rahman M. M., Saidur R., Rahimn A. Conjugated effect of joule heating and magneto-hydrodynamic on doublediffusive mixed convection in a horizontal channel with an open cavity [J].International Journal of Heat and Mass Transfer, 2011, 54(15): 3201-3213.

    [6] Mehrez Z., Cafsi A. E., Belghith A. et al. Effect of heated wall position on heat transfer and entropy generation of Cu water nanofluid flow in an open cavity [J].Canadian Journal of Physics, 2015, 93(12): 1615-1629.

    [7] Abdelmassih G., Vernet A., Pallares J. Numerical simulation of incompressible laminar flow in a three-dimensional channel with a cubical open cavity with a bottom wall heated [J].Journal of Physics: Conference series, 2012, 395(1): 012099

    [8] Stiriba Y., Ferré J. A., Grau F. X. Heat transfer and fluid flow characteristics of laminar flow past an opencavity with heating from below [J].International Communications in Heat and Mass Transfer, 2013, 43(4): 8-15.

    [9] Khanafer K., Al-Azmi B., Al-Shammar A. et al. Mixed convection analysis of laminar pulsating flow and heat transfer over a backward-facing step [J].International Journal of Heat and Mass Transfer, 2008, 51(25): 5785-5793

    [10] Khanafer K. M., Al-Amiri A. M., Pop I. Numerical simulation of unsteady convection in a driven cavity using an externally exited sliding lid [J].European Journal of Mechanics B/Fluids, 2007, 26(5): 669-687.

    [11] Velazquez A., Arias J., Montanes J. Pulsating flow and convective heat transfer in a cavity with inlet and outlet sections [J]International Journal of Heat and Mass Transfer, 2009, 52(3-4): 647-654.

    [12] Sourtiji E., Hosseinizadeh S. F., Gorji-Bandpy M. et al. Heat transfer enhancement of mixed convection in a square cavity with inlet and outlet ports due to oscillation of incoming flow [J].International Communications in Heat and Mass Transfer, 2011, 38(6): 806-814.

    [13] Selimefendigil F., Oztop H. Identification of forced convection in pulsating flow at a backward facing step with a stationary cylinder subjected to nanofluid [J].International Communications in Heat and Mass Transfer, 2013, 45(7): 111-121.

    [14] Unal A., Selma A., Dogan D. Heat transfer enhancement with laminar pulsating nanofluid flow in a wavy channel [J].International Communications in Heat and Mass Transfer, 2014, 59: 17-23

    [15] Selimefendigil F., Yurddas A. Numerical analysis of mixed convection heat transfer in pulsating flow for a horizontal channel with a cavity heated from vertical side and below [J].Heat Transfer Research, 2012, 43(6): 509-525.

    [16] Abu-Nada E. Numerical prediction of entropy generation in separated flows [J].Entropy, 2005, 7(4): 234-252.

    [17] Abu-Nada E. Entropy generation due to heat and fluid flow in backward facing step flow with various expansion ratios [J].International Journal of Exergy, 2006, 3(4): 419-435.

    [18] Khorasanizadeh H., Nikfar M., Amani J. Entropy generation of Cu-water nanofluid mixed convection in a cavity [J].European Journal of Mechanics, 2013, 37(28): 143-152.

    [19] Nayak R. K., Bhattacharyya S., Pop I. Numerical study on mixed convection and entropy generation of Cu-water nanofluid in a differentially heated skewed enclosure [J].International Journal of Heat and Mass Transfer, 2015, 85(6): 620-634.

    [20] Sciacovelli A., Verda V., Sciubba E. Entropy generation analysis as a design tool-A review [J].Renewable and Sustainable Energy Reviews, 2015, 43: 1167-1181.

    [21] Mehrez Z., Bouterra M., Cafsi A. E. et al. Heat transfer and entropy generation analysis of nanofluids flow in an open cavity [J].Computers and Fluids, 2013, 88(4): 363-373

    [22] Mehrez Z., Cafsi A. E., Belghith A. et al. MHD effects on heat transfer and entropy generation of nanofluidflow in an open cavity [J].Magnetism and Magnetic Materials, 2015, 374: 214-224.

    [23] Mehrez Z., Cafsi A. E., BELGHITH A. et al. The entropy generation analysis in the mixed convective assisting flow of Cu-water nanofluid in an inclined open cavity [J].Advanced Powder Technology, 2015, 26(5): 1442-1451.

    [24] Zamzari F., Mehrez Z., Cafsi A. E. et al. Entropy generation and mixed convection in a horizontal channel with an open cavity [J].International Journal of Exergy, 2015, 17(2): 219-239.

    [25] Bejan A. A study of entropy generation in fundamental convective heat transfer [J].Journal of Heat Transfer, 1979, 101: 718-725.

    [26] Mehrez Z., Bouterra M., Cafsi A. E. et al. Mass transfer control of a backward-facing step flow by local forcingeffect of Reynolds number [J].Thermal Science, 2011, 15(2): 367-378.

    [27] Nasri Z., Laatar A. H., Balti J. Natural convection enhancement in an asymmetrically heated channelchimney system [J].International Journal of Thermal Sciences, 2015, 90: 122-134.

    [28] LEONARD B. P. Simple high accuracy resolution program for convection modeling of discontinuities [J].International Journal for Numerical Methods in Fluids, 1988, 8(10): 1291-1318.

    [29] Mehrez Z., Bouterra M., Cafsi A. E. et al. The influence of the periodic disturbance on the local heat transfer in separated and reattached flow [J].Heat Mass Transfer, 2009, 46(1): 107-112.

    [30] Taieb S., Laatar A. H., Balti J. Natural convection in an asymmetrically heated vertical channel with an adiabatic auxiliary plate [J].International Journal of Thermal Sciences, 2013, 74(6): 24-36.

    [31] Xing T., Stern F. Factors of safety for Richardson extrapolation [J].Journal of Fluids Engineering, 2010, 132(6): 061403.

    [32] Xing T., Stern F. Closure to “Discussion of “Factors of Safety for Richardson Extrapolation” 2011, J. Fluids Eng., 133, 11550 [J].Journal of Fluids Engineering, 2011, 133(11): 115502.

    [33] Xing T. Direct numerical simulation of Open Von Kármán Swirling Flow [J].Journal of Hydrodynamics, 2014, 26(2): 165-177.

    [34] Xing T., Carrica P., Stern F. Computational towing tank procedures for single run curves of resistance and propulsion [J].Journal of Fluids Engineering, 2008, 130(10): 101102.

    (Received May 2, 2016, Revised May 23, 2017)

    *Biography:Fatma Zamzari (1986-), Female, Ph. D., Researcher

    Zouhaier Mehrez,

    E-mail: zouhaier.mehrez@yahoo.fr

    国产一区二区三区视频了| 黄色视频,在线免费观看| 日韩av在线大香蕉| 亚洲男人的天堂狠狠| 在线观看舔阴道视频| 亚洲中文字幕日韩| 久久精品国产亚洲av香蕉五月| 日韩欧美免费精品| 国产精品伦人一区二区| 搡老熟女国产l中国老女人| 老熟妇仑乱视频hdxx| 香蕉av资源在线| 日本a在线网址| 欧美色欧美亚洲另类二区| 看免费av毛片| 午夜福利在线观看免费完整高清在 | 男人舔奶头视频| 亚洲美女视频黄频| 日本五十路高清| 亚洲熟妇熟女久久| 久久久久九九精品影院| 天堂网av新在线| 日韩中文字幕欧美一区二区| 丰满的人妻完整版| 十八禁人妻一区二区| 久久伊人香网站| 亚洲国产精品合色在线| 国产高清激情床上av| 欧美zozozo另类| 极品教师在线视频| 免费看光身美女| 日本与韩国留学比较| 国产欧美日韩精品亚洲av| 高清毛片免费观看视频网站| 亚洲欧美日韩高清在线视频| 精品人妻偷拍中文字幕| 我要看日韩黄色一级片| 成年免费大片在线观看| 又黄又爽又免费观看的视频| 久久99热这里只有精品18| 女人十人毛片免费观看3o分钟| 国产淫片久久久久久久久 | 少妇被粗大猛烈的视频| 午夜福利在线观看吧| 国产久久久一区二区三区| 伊人久久精品亚洲午夜| 91字幕亚洲| 国产精品久久久久久亚洲av鲁大| 精品久久久久久久久久久久久| 级片在线观看| www.www免费av| xxxwww97欧美| 久久草成人影院| 免费在线观看成人毛片| 国产单亲对白刺激| 国产精品电影一区二区三区| 少妇熟女aⅴ在线视频| 1000部很黄的大片| 国产精品久久久久久亚洲av鲁大| 少妇人妻一区二区三区视频| 美女高潮喷水抽搐中文字幕| 一夜夜www| 亚洲欧美清纯卡通| 亚洲一区高清亚洲精品| 看免费av毛片| 波多野结衣巨乳人妻| 日本成人三级电影网站| 国产精品综合久久久久久久免费| 一级a爱片免费观看的视频| 最近视频中文字幕2019在线8| 欧美成人一区二区免费高清观看| 精品熟女少妇八av免费久了| 精品99又大又爽又粗少妇毛片 | 老司机午夜福利在线观看视频| 日韩欧美在线乱码| 日韩中文字幕欧美一区二区| 简卡轻食公司| 国产日本99.免费观看| 波多野结衣巨乳人妻| 亚洲人成网站在线播放欧美日韩| 国产不卡一卡二| 极品教师在线免费播放| 男人和女人高潮做爰伦理| 精品日产1卡2卡| 久久久久免费精品人妻一区二区| 欧美一区二区精品小视频在线| 亚洲成人中文字幕在线播放| 丰满乱子伦码专区| 青草久久国产| 中文在线观看免费www的网站| 中文字幕av在线有码专区| 九九在线视频观看精品| 久久久久久久久大av| 亚洲最大成人av| a在线观看视频网站| 欧美黄色淫秽网站| 国产视频内射| 午夜日韩欧美国产| 夜夜看夜夜爽夜夜摸| 精品人妻视频免费看| 少妇的逼好多水| av视频在线观看入口| 国产中年淑女户外野战色| 尤物成人国产欧美一区二区三区| 精品福利观看| 国产av在哪里看| 91麻豆精品激情在线观看国产| 男女之事视频高清在线观看| xxxwww97欧美| 欧美高清成人免费视频www| 亚洲人成网站在线播放欧美日韩| 久久欧美精品欧美久久欧美| 亚洲综合色惰| 日韩欧美免费精品| 桃红色精品国产亚洲av| 免费在线观看成人毛片| 精品一区二区三区av网在线观看| 日本熟妇午夜| 国产高清视频在线播放一区| 久久这里只有精品中国| 深爱激情五月婷婷| 1024手机看黄色片| 日本黄色片子视频| 脱女人内裤的视频| 国产精品1区2区在线观看.| 亚洲精品在线观看二区| 亚洲av免费在线观看| 久久久久久久精品吃奶| 最新在线观看一区二区三区| av专区在线播放| 日本五十路高清| 久久6这里有精品| 在线播放国产精品三级| 一进一出抽搐gif免费好疼| 亚洲av电影在线进入| 免费观看人在逋| 蜜桃久久精品国产亚洲av| 在线a可以看的网站| 长腿黑丝高跟| 欧美丝袜亚洲另类 | 大型黄色视频在线免费观看| 久久久久久久精品吃奶| 真人做人爱边吃奶动态| 亚洲片人在线观看| 69av精品久久久久久| 18禁黄网站禁片午夜丰满| 免费在线观看亚洲国产| 99久久成人亚洲精品观看| 小说图片视频综合网站| 婷婷色综合大香蕉| 深爱激情五月婷婷| 男人的好看免费观看在线视频| 成人亚洲精品av一区二区| 真人做人爱边吃奶动态| 永久网站在线| 最新中文字幕久久久久| 日韩欧美免费精品| 午夜激情福利司机影院| 欧美日本视频| 欧美+亚洲+日韩+国产| 麻豆国产97在线/欧美| 亚洲欧美清纯卡通| 国产av一区在线观看免费| 亚洲精品一区av在线观看| 少妇的逼好多水| 国产精品嫩草影院av在线观看 | 午夜精品在线福利| 久久久久久久久久成人| 五月玫瑰六月丁香| 一卡2卡三卡四卡精品乱码亚洲| 动漫黄色视频在线观看| 国产伦人伦偷精品视频| 日日夜夜操网爽| 精品久久久久久,| 久久天躁狠狠躁夜夜2o2o| 看免费av毛片| 亚洲片人在线观看| 久久久久性生活片| 精品国产亚洲在线| 亚洲欧美日韩卡通动漫| 国产欧美日韩精品一区二区| 久久亚洲精品不卡| 久久久久性生活片| 亚洲精品一卡2卡三卡4卡5卡| 亚洲av成人av| 成熟少妇高潮喷水视频| 亚洲av成人不卡在线观看播放网| 99久久无色码亚洲精品果冻| 国产伦精品一区二区三区四那| 91av网一区二区| 18禁裸乳无遮挡免费网站照片| 日日摸夜夜添夜夜添小说| 亚洲精品成人久久久久久| 国产精华一区二区三区| 成人国产综合亚洲| 久久久精品欧美日韩精品| 国产精品永久免费网站| av在线观看视频网站免费| 性插视频无遮挡在线免费观看| 看免费av毛片| 一进一出好大好爽视频| 久久精品国产99精品国产亚洲性色| 国产视频一区二区在线看| 成人特级av手机在线观看| 久久精品国产亚洲av天美| 中文字幕av成人在线电影| 亚洲黑人精品在线| 国产伦在线观看视频一区| 9191精品国产免费久久| 久久99热这里只有精品18| 国产成人aa在线观看| 国产白丝娇喘喷水9色精品| 床上黄色一级片| 色综合欧美亚洲国产小说| 麻豆成人av在线观看| 午夜福利免费观看在线| 精品一区二区三区av网在线观看| 国产主播在线观看一区二区| 亚洲欧美日韩高清专用| 国产色婷婷99| 精品人妻1区二区| 男人的好看免费观看在线视频| 亚洲不卡免费看| 国产老妇女一区| 亚洲av.av天堂| 午夜精品在线福利| 色综合婷婷激情| 亚洲乱码一区二区免费版| 一级毛片久久久久久久久女| 高清在线国产一区| 日本黄大片高清| 国产一区二区三区在线臀色熟女| 校园春色视频在线观看| 免费av不卡在线播放| 色吧在线观看| 国产欧美日韩一区二区精品| 国产伦精品一区二区三区四那| 亚洲乱码一区二区免费版| 久久午夜福利片| 99热精品在线国产| 日韩欧美国产一区二区入口| 观看美女的网站| 老司机午夜十八禁免费视频| 内射极品少妇av片p| 少妇人妻一区二区三区视频| av视频在线观看入口| 99久久成人亚洲精品观看| 亚洲一区二区三区色噜噜| 深爱激情五月婷婷| 直男gayav资源| 欧美三级亚洲精品| 国产一区二区激情短视频| 久久99热6这里只有精品| 欧美不卡视频在线免费观看| 91麻豆精品激情在线观看国产| 99精品久久久久人妻精品| 国产精品不卡视频一区二区 | 日韩亚洲欧美综合| 看黄色毛片网站| 丁香欧美五月| 男人和女人高潮做爰伦理| 亚洲av中文字字幕乱码综合| 亚洲精品在线观看二区| 欧美zozozo另类| 在线播放国产精品三级| 亚洲 国产 在线| 久久6这里有精品| 国产精品美女特级片免费视频播放器| 天堂影院成人在线观看| 永久网站在线| 欧美性猛交╳xxx乱大交人| 久久伊人香网站| 九九在线视频观看精品| 一二三四社区在线视频社区8| 欧美黄色淫秽网站| 99久久九九国产精品国产免费| 99热只有精品国产| 欧美丝袜亚洲另类 | 999久久久精品免费观看国产| 五月玫瑰六月丁香| 精品一区二区三区视频在线| 中文字幕av在线有码专区| 午夜免费成人在线视频| 搡老妇女老女人老熟妇| 日韩欧美在线二视频| 亚洲 欧美 日韩 在线 免费| 有码 亚洲区| 欧美黄色淫秽网站| 亚洲经典国产精华液单 | 99热这里只有是精品50| 午夜福利在线在线| 欧美成人性av电影在线观看| 天堂影院成人在线观看| 国产爱豆传媒在线观看| 欧美日韩乱码在线| 国产av一区在线观看免费| 一区二区三区高清视频在线| 欧美性感艳星| 国产精品永久免费网站| 免费在线观看影片大全网站| a级毛片a级免费在线| 国产69精品久久久久777片| 999久久久精品免费观看国产| 亚洲精华国产精华精| 少妇的逼水好多| 久久精品国产自在天天线| 亚洲欧美日韩卡通动漫| 亚洲美女黄片视频| 久久6这里有精品| 亚洲欧美精品综合久久99| 又紧又爽又黄一区二区| 亚洲18禁久久av| 女人被狂操c到高潮| 亚洲成人免费电影在线观看| 观看美女的网站| 国产色婷婷99| 最新中文字幕久久久久| 麻豆成人午夜福利视频| 91九色精品人成在线观看| 宅男免费午夜| 91麻豆av在线| 男人舔奶头视频| 在现免费观看毛片| 一级黄色大片毛片| 亚洲精品粉嫩美女一区| 极品教师在线视频| 免费在线观看影片大全网站| 嫩草影院新地址| 蜜桃久久精品国产亚洲av| 免费在线观看影片大全网站| 少妇裸体淫交视频免费看高清| 国产精品久久电影中文字幕| 国产免费一级a男人的天堂| 一个人免费在线观看电影| 国产 一区 欧美 日韩| 成人毛片a级毛片在线播放| 制服丝袜大香蕉在线| 精品久久久久久久久久免费视频| 最好的美女福利视频网| 男人的好看免费观看在线视频| 中文资源天堂在线| www.色视频.com| 免费在线观看亚洲国产| 亚洲专区国产一区二区| a级毛片免费高清观看在线播放| 内地一区二区视频在线| 国产乱人视频| 久久午夜福利片| 国产熟女xx| 亚洲av日韩精品久久久久久密| 国产欧美日韩精品一区二区| a级毛片免费高清观看在线播放| 国产精品电影一区二区三区| 亚洲无线在线观看| 91狼人影院| 极品教师在线免费播放| 两性午夜刺激爽爽歪歪视频在线观看| 国产黄色小视频在线观看| 好看av亚洲va欧美ⅴa在| 乱人视频在线观看| 国产免费男女视频| 两个人视频免费观看高清| 亚洲欧美日韩东京热| 观看美女的网站| 国产精品99久久久久久久久| 两个人视频免费观看高清| 欧美精品啪啪一区二区三区| 天天躁日日操中文字幕| 一级黄色大片毛片| 欧美成狂野欧美在线观看| 亚洲中文字幕日韩| 麻豆国产av国片精品| 99久久久亚洲精品蜜臀av| 午夜福利视频1000在线观看| 成人国产综合亚洲| 欧美日韩中文字幕国产精品一区二区三区| 在线播放国产精品三级| 久久精品夜夜夜夜夜久久蜜豆| 91麻豆av在线| 九九久久精品国产亚洲av麻豆| 欧美高清性xxxxhd video| 国产主播在线观看一区二区| 国产大屁股一区二区在线视频| 精品乱码久久久久久99久播| 国产精品电影一区二区三区| 我要看日韩黄色一级片| 亚洲欧美日韩东京热| 欧美zozozo另类| 国产白丝娇喘喷水9色精品| 精品久久久久久久末码| 亚洲国产高清在线一区二区三| 久久九九热精品免费| 国产免费一级a男人的天堂| 欧美日韩福利视频一区二区| 人人妻人人澡欧美一区二区| 国产精品伦人一区二区| 久久这里只有精品中国| 99热这里只有是精品50| 国产精品日韩av在线免费观看| 午夜免费男女啪啪视频观看 | 日韩欧美国产一区二区入口| 国产成人欧美在线观看| 久久亚洲精品不卡| 日本一本二区三区精品| 久久久久久国产a免费观看| 欧美高清成人免费视频www| 国内毛片毛片毛片毛片毛片| 国内精品一区二区在线观看| 久久久久国内视频| 成人av一区二区三区在线看| 免费在线观看成人毛片| 婷婷色综合大香蕉| 亚洲天堂国产精品一区在线| 天堂√8在线中文| 精品久久久久久久人妻蜜臀av| 国产精品自产拍在线观看55亚洲| 久久久久国产精品人妻aⅴ院| 日本 欧美在线| 老司机深夜福利视频在线观看| 最近中文字幕高清免费大全6 | 1000部很黄的大片| 久久久久性生活片| 日韩中字成人| 国产在视频线在精品| 精品午夜福利在线看| 日日夜夜操网爽| 亚洲狠狠婷婷综合久久图片| 欧美激情久久久久久爽电影| 高清毛片免费观看视频网站| 亚洲专区中文字幕在线| 国内毛片毛片毛片毛片毛片| 91在线观看av| 免费搜索国产男女视频| 午夜福利在线观看免费完整高清在 | 久久久久久久久久成人| 日本 欧美在线| 亚洲第一区二区三区不卡| 欧美日韩福利视频一区二区| 亚洲欧美日韩无卡精品| 亚洲国产精品久久男人天堂| 亚洲av成人av| 亚洲一区二区三区不卡视频| 色精品久久人妻99蜜桃| 色哟哟哟哟哟哟| 国产欧美日韩一区二区精品| 国产亚洲精品久久久久久毛片| 狠狠狠狠99中文字幕| 内射极品少妇av片p| 天天躁日日操中文字幕| 12—13女人毛片做爰片一| aaaaa片日本免费| 精品久久久久久久久亚洲 | 国产老妇女一区| av女优亚洲男人天堂| av中文乱码字幕在线| 成熟少妇高潮喷水视频| 久久亚洲精品不卡| 免费观看精品视频网站| 我要看日韩黄色一级片| 一个人看的www免费观看视频| 免费观看的影片在线观看| 国内少妇人妻偷人精品xxx网站| 免费在线观看亚洲国产| 久久亚洲精品不卡| 日日夜夜操网爽| 制服丝袜大香蕉在线| 亚洲专区国产一区二区| 可以在线观看的亚洲视频| 亚洲欧美日韩无卡精品| 国产综合懂色| 亚洲狠狠婷婷综合久久图片| av国产免费在线观看| 国产精品一区二区性色av| 日本 欧美在线| 国内精品一区二区在线观看| 观看美女的网站| 成人永久免费在线观看视频| a级毛片免费高清观看在线播放| 丁香六月欧美| 天堂√8在线中文| 99热精品在线国产| 波多野结衣高清作品| 夜夜躁狠狠躁天天躁| 精品午夜福利在线看| 亚洲人成网站高清观看| 91麻豆精品激情在线观看国产| 亚洲第一电影网av| 亚洲综合色惰| 国产又黄又爽又无遮挡在线| 九九久久精品国产亚洲av麻豆| 乱码一卡2卡4卡精品| 国产三级黄色录像| 婷婷亚洲欧美| 国产真实乱freesex| 91av网一区二区| 十八禁国产超污无遮挡网站| 俄罗斯特黄特色一大片| 九九在线视频观看精品| 久久久久性生活片| 欧美午夜高清在线| 9191精品国产免费久久| 一区二区三区免费毛片| 波多野结衣高清无吗| 亚洲avbb在线观看| 免费高清视频大片| 日韩欧美在线乱码| 亚洲av一区综合| 老司机深夜福利视频在线观看| 一区二区三区四区激情视频 | 国产黄色小视频在线观看| 欧美国产日韩亚洲一区| 欧美日本亚洲视频在线播放| 国产69精品久久久久777片| 亚洲电影在线观看av| 亚洲 国产 在线| 久久午夜亚洲精品久久| 精品一区二区三区人妻视频| 美女 人体艺术 gogo| 一级a爱片免费观看的视频| 人妻久久中文字幕网| 精品一区二区免费观看| 日本在线视频免费播放| 久久婷婷人人爽人人干人人爱| 亚洲不卡免费看| 九九久久精品国产亚洲av麻豆| 高清在线国产一区| 黄片小视频在线播放| 老师上课跳d突然被开到最大视频 久久午夜综合久久蜜桃 | 美女黄网站色视频| 国产色爽女视频免费观看| 亚洲激情在线av| av黄色大香蕉| 宅男免费午夜| 99热这里只有是精品50| АⅤ资源中文在线天堂| 欧美潮喷喷水| 90打野战视频偷拍视频| 日日摸夜夜添夜夜添av毛片 | 啪啪无遮挡十八禁网站| 亚洲色图av天堂| 日韩高清综合在线| 色哟哟哟哟哟哟| 国产美女午夜福利| 亚洲真实伦在线观看| 亚洲精品456在线播放app | 国产成年人精品一区二区| 啦啦啦韩国在线观看视频| 精品久久国产蜜桃| 国产精品不卡视频一区二区 | 日韩国内少妇激情av| 免费搜索国产男女视频| 国产亚洲精品久久久久久毛片| 麻豆国产av国片精品| 神马国产精品三级电影在线观看| 欧美中文日本在线观看视频| 国内揄拍国产精品人妻在线| 欧美激情在线99| 男女那种视频在线观看| 精品乱码久久久久久99久播| 日本一本二区三区精品| 久久人人爽人人爽人人片va | 亚洲精品一卡2卡三卡4卡5卡| 一边摸一边抽搐一进一小说| 久久香蕉精品热| 男插女下体视频免费在线播放| 我的老师免费观看完整版| 99热6这里只有精品| 97人妻精品一区二区三区麻豆| 校园春色视频在线观看| 国产91精品成人一区二区三区| 久久久精品大字幕| 欧美性猛交黑人性爽| 欧美激情国产日韩精品一区| 夜夜爽天天搞| 欧美三级亚洲精品| 51国产日韩欧美| 村上凉子中文字幕在线| 成熟少妇高潮喷水视频| 国产欧美日韩一区二区三| 免费观看精品视频网站| 舔av片在线| 激情在线观看视频在线高清| 18禁裸乳无遮挡免费网站照片| 亚洲第一区二区三区不卡| 夜夜夜夜夜久久久久| 亚洲av成人av| 亚洲精品一卡2卡三卡4卡5卡| 亚洲成a人片在线一区二区| 国产成年人精品一区二区| 婷婷亚洲欧美| 欧美色欧美亚洲另类二区| 久久中文看片网| 天天一区二区日本电影三级| 亚洲人与动物交配视频| 国产欧美日韩一区二区三| 两性午夜刺激爽爽歪歪视频在线观看| 能在线免费观看的黄片| 91麻豆精品激情在线观看国产| 伊人久久精品亚洲午夜| 国产成人av教育| 亚洲av美国av| 亚洲色图av天堂| 成年免费大片在线观看| 亚洲内射少妇av| 免费av不卡在线播放| 国产成人欧美在线观看| 国产精品日韩av在线免费观看| 亚洲中文字幕日韩| 欧美日韩亚洲国产一区二区在线观看| or卡值多少钱| 国内精品美女久久久久久| 日本在线视频免费播放| 日本免费a在线| 嫩草影视91久久| 亚洲在线观看片|