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

    Numerical simulation and analysis of confined turbulent buoyant jet with variable source*

    2015-12-01 02:12:23ELAMINMohamedALGHAMDIAbdulmajeedSALAMAAmgadSUNShuyu1MathematicsDepartmentFacultyofScienceAswanUniversityAswan8158EgyptComputationalTransportPhenomenaLaboratoryCTPLDivisionofPhysicalSciencesandEngineeringPSEKingAbdullahUni

    EL-AMIN Mohamed F., AL-GHAMDI Abdulmajeed, SALAMA Amgad, SUN Shuyu1. Mathematics Department, Faculty of Science, Aswan University, Aswan 8158, Egypt. Computational Transport Phenomena Laboratory (CTPL), Division of Physical Sciences and Engineering(PSE), King Abdullah University of Science and Technology (KAUST), Thuwal 955-6900, Jeddah, Kingdom of Saudi Arabia, E-mail: mohamed.elamin@kaust.edu.sa. College of Engineering and Islamic Architecture, Umm Al-Qura University, Makkah, Kingdom of Saudi Arabia

    Numerical simulation and analysis of confined turbulent buoyant jet with variable source*

    EL-AMIN Mohamed F.1,2, AL-GHAMDI Abdulmajeed3, SALAMA Amgad2, SUN Shuyu2
    1. Mathematics Department, Faculty of Science, Aswan University, Aswan 81528, Egypt
    2. Computational Transport Phenomena Laboratory (CTPL), Division of Physical Sciences and Engineering(PSE), King Abdullah University of Science and Technology (KAUST), Thuwal 23955-6900, Jeddah, Kingdom of Saudi Arabia, E-mail: mohamed.elamin@kaust.edu.sa
    3. College of Engineering and Islamic Architecture, Umm Al-Qura University, Makkah, Kingdom of Saudi Arabia

    (Received August 8, 2014, Revised December 16, 2014)

    In this work, experimental and numerical investigations are undertaken for confined buoyant turbulent jet with varying inlet temperatures. Results of the experimental work and numerical simulations for the problem under consideration are presented. Four cases of different variable inlet temperatures and different flow rates are considered. The realizable k-εturbulence model is used to model the turbulent flow. Comparisons show good agreements between simulated and measured results. The average deviation of the simulated temperature by realizable k-εturbulent model and the measured temperature is within 2%. The results indicate that temperatures along the vertical axis vary, generally, in nonlinear fashion as opposed to the approximately linear variation that was observed for the constant inlet temperature that was done in a previous work. Furthermore, thermal stratification exits,particularly closer to the entrance region. Further away from the entrance region the variation in temperatures becomes relatively smaller. The stratification is observed since the start of the experiment and continues during the whole course. Numerical experiments for constant, monotone increasing and monotone decreasing of inlet temperature are done to show its effect on the buoyancy force in terms of Richardson number.

    realizable k-εmodel, turbulent jet, stratification, heat stores, CFD

    Introduction0F

    Although stratification of fluids due to the existence of temperature gradients is not desirable in many processes that require homogenization, it is, in other processes (e.g., heat storage tanks) desirable because of the low mixing mechanisms involved which help to maintain the required temperature distribution. The problem, however, in thermally stratified heat storage tanks is its sensitivity to external disturbance. That is the state of thermal stratification could be destroyed once sufficient disturbance is introduced. In particular,the conditions at the inlet are considered as one example of such disturbance sources. Therefore, it is important to study the effect of inlet conditions on stratification behavior of such systems. Since, in heat storage tanks and in many other applications, fluids enter the tank in the form of buoyant jet, great deal of works on jet flow have been published either experimentally or numerically. However, the problem of buoyant heated jet with variable source temperature as can be found in many industrial and environmental applications, has received relatively little attention. Jet flow can be divided mainly into three types, pure jet, pure plume and forced plume. In pure jets, fluids enter the domain with high momentum fluxes that essentially cause high intensity of turbulent mixing. In pure plume, on the other hand, buoyancy fluxes cause local acceleration leading to turbulent mixing. In the general case of a forced plume a combination of initial momentum and buoyancy fluxes are responsible for turbulent mixing. Several techniques were developed to study jets in confined spaces as will be explained later. Lately, with the increase of computers efficiencies andcapacities, computational fluid dynamics (CFD) became one of the essential tools to exploring the fluids' behavior of such fundamental importance. In turbulent flows the Reynolds average Navier-Stokes (RANS)technique are usually adopted in order to make the system amenable to solution. The problem in using RANS approach however, is that till now, there is no uniform set of equations to model all kinds of turbulent flows and heat transfer scenarios. Therefore, it is important to choose the model which suites the case under investigation and even to calibrate its coefficients in order to fit experimental results. El-Amin et al.[1]investigated the 2-D upward, axisymmetric turbulent confined jet and developed several models to describe flow patterns using realizable k-εturbulence model. Furthermore, CFD analysis of the flow structure of a horizontal water injected jet into a rectangular tank has been done by El-Amin et al.[2]. Their findings were later used by Panthalookaran et al.[3]to calibrate both realizable and RNG k-εturbulence models so that they may be used for simulating stratified hot water storage tanks.

    Many experimental works were conducted to highlight the interesting patterns and the governing parameters pertinent to this kind of flows including the work of Agrawal and Prasad[4]and O'Hern et al.[5]and many other authors have investigated several phenomena pertinent to buoyant jets. For example, the problem of bifurcation in a buoyant horizontal laminar jet was studied by Arakeri et al.[6]. With respect to the kind of fluids used in buoyant jet studies, several researchers have considered different fluids for either experimental or numerical investigations. For example O'Hern et al.[5]performed experimental work on a turbulent buoyant helium plume. El-Amin and Kanayama[7,8]studied buoyant jet resulting from hydrogen leakage. They developed the similarity formulation and obtained solutions of the centerline quantities such as velocity and concentration. Furthermore,El-Amin[9]introduced a numerical investigation of a vertical axisymmetric non-Boussinesq buoyant jet resulting from hydrogen leakage in air as an example of injecting a low-density gas into high-density ambient. On the other hand, the mechanics of buoyant jet flows issuing from a general 3-D geometry into an unbounded ambient environment with uniform density or stable density stratification and under stagnant or steady sheared current conditions was investigated by Jirka[10]. An integral model is formulated for the conservation of mass, momentum, buoyancy and scalar quantities in the turbulent jet flow. Furthermore Jirka[11]extended also this work to account for the plane buoyant jet dynamics resulting from the interaction of multiple buoyant jet effluxes spaced along a diffuser line. In the previous work by El-Amin et al.[1],analyses of the components of 2-D axisymmetric vertical unheated/heated turbulent confined jet using realizable k-εturbulence model were conducted. Moreover experimental work was elaborated for temperature measurements of such system to provide verification of the models used. In that work, several models were considered to describe axial velocity, centerline velocity, radial velocity, dynamic pressure, mass flux,momentum flux and buoyancy flux for both unheated(non-buoyant) and heated (buoyant) jet. In that work inlet temperatures were considered fixed. However, in many applications inlet temperatures are not, generally fixed. Many recent related papers can be found in Refs.[12-16].

    Fig.1 Schematic diagram of the experimental setup

    Fig.2 Schematic representation of the calculation domain

    Table 1 Summary of the experimental data with variable inlet temperature

    In this work, analysis of vertical hot water jet entering a cylindrical tank filled with cold water with variable inlet temperature is conducted. The inlet temperature of the buoyant jet is allowed to change within a small range and is considered as a function of time. Numerical investigations under the above mentioned conditions are performed in order to obtain fields of pressure, velocity, temperature and turbulence. 2-D axisymmetric simplification is assumed to reduce the grid size in the solution domain and the realizable k-εmodel is used to model turbulent flow.

    1. Measurements

    Schematic diagram of the experimental setup is shown in Fig.1. The cylindrical tank, made of 0.005 m thick galvanized iron sheets, has a diameter of 0.36 m and a height of 0.605 m is shown in Fig.2(a). The inlet pipe is located at center of the bottom of the tank with an inner diameter of 0.02 m. the inlet pipe is inserted in the tank up to a height of 0.06 m above the base of the tank. The outlet pipe, located at center of the top of the tank, has an inner diameter of 0.02 m with a depth in the tank of 0.055 m from the top plate. The uncertainty in the diameter of the inlet and outlet pipes is ±0.001m. This geometry suggests that the tank, the inlet, and the outlet may be modeled as axisymmetric around the vertical axis. Thermal effects are measured by thermocouples of K-type that were calibrated against a standard PT-25 resistance thermometer with an average calibration error of ±0.25 K. The temperature is recorded in Kelvin using a data acquisition system connected with a personal computer. The data acquisition system has an error about ±1K. Moreover, the zero-order error in terms of resolution error for the thermocouples is ±0.005 K. Flow rate is measured using a magnetic-type flowmeter. The operation principle of a magnetic flowmeter is based on the Faraday's law of electromagnetic induction. When water flows through the magnetic field, a small voltage is induced which is proportional to the velocity of the flow and is accurately measured by two stainless steel electrodes mounted opposing each other inside the metering pipe. The resolution error of the voltages read by the flowmeter is negligible because it reads until 6 decimal digits. The instrumental error of the flowmeter is ±3.5%at various values of temperature.

    Fig.3 Variable inlet (a) temperature, (b) Reynolds number and(c) turbulence intensity as a function of time for Cases 1-4

    The above estimated errors are included in the measured data. Thermal effects are measured by inserting a vertical rod with 9 stainless steel sheathed K-type thermocouples. The nine sensors are distributed at different heights as, 0.06 m, 0.12 m, 0.18 m,0.24 m, 0.30 m, 0.36 m, 0.42 m, 0.48 m and 0.54 m measured from the bottom. The distance between the symmetry axis and the thermocouples' rod is 0.09 m,i.e., in the middle between the symmetry axis and the tank wall. The inlet temperature was measured using another thermocouple which is located at the inlet pipe. It is important to indicate that all precautions have been taken to make sure the geometrical symmetry is achieved in the sense of the alignment of the inlet and outlet tubes, the smoothness of tube materials, the inlet geometry, etc..

    For the variable inlet temperature, the used parameters are listed in Table 1. The duration for measurements for each case was approximately 30 min. The inlet temperature and the initial temperature are given in column 2 and 3 of the Table 1, respectively, and the flow rate is given in column 4. The inlet velocity,Reynolds number and turbulence intensity at the inlet nozzle are calculated from the given data.

    The best fitting for the given curves of the inlet variable temperature, Fig.3(a), can be represented by a 5th order polynomial as a function of time, as given by Eqs.(A1)-(A4) in the Appendix. The corresponding Reynolds numbers, Fig.3(b), and inlet turbulence intensity, Fig.3(c), are calculated as functions of inlet temperature which in turn have been represented by 5th order polynomials of time, which are given in Eqs.(A5)-(A8) and (A9)-(A12), respectively, with the aid of Eqs.(13)-(14). These functions may be represented by other polynomials with less order but the deviation from the measured data will increase.

    The following empirical relation (Fluent User's Guide, Fluent Inc.) is used to describe the turbulent intensity at the inlet nozzle as a function of the Reynolds number,I0=0.16(Re)-0.125. The Reynolds number, upon considering the inlet diameter as a length scale, is given by,Re=u0d/n.

    The measured temperature profiles for variable inlet temperature are plotted as a function of time, at different sensors positions (Cases 1-4) in Figs.4(a)-4(d). It is apparent that the temperatures along the vertical axis vary in nonlinear fashion with time, especially, at larger heights z≥0.18 m(plume region). Also, the figures indicate that the temperature increases as time increases. Thermal stratification is observed looking at the difference in temperature from top (higher)to bottom (lower). The degree of stratification, however, seems to be more pronounced in the lower half of the tank than in the top half. The stratification is verified from the beginning of the experiment and continues during the whole time.

    Fig.4 Profiles of measured temperature as a function of time, at different sensors positions, for variable inlet temperature,Cases 1-4

    2. Mathematical formulation

    A comparison was done by El-Amin et al.[1]totest different turbulence models when simulating confined buoyant jet, and they reported that the best model is the realizable k-εmodel. Therefore, in this work we use this model to simulate the problem under consideration. The realizable k-εmodel developed by Shih et al.[17]involves a new eddy-viscosity formula and a new model equation for dissipationεbased on the dynamic equation of the mean-square vorticity fluctuations. The RANS are given in Eqs.(2)-(3), and the energy equation is represented by Eq.(4). The governing equations of mass, momentum and turbulence are given below.

    Continuity equation:

    Radial momentum equation:

    Energy equation:

    Turbulent kinetic energy(k )equation:

    Turbulence dissipation rate(ε)equation:

    In the above equations,uandvare the mean axial and radial velocity components, respectively. The other quantities are time,t , densityρ, acceleration due to gravity,g , pressure,p , dynamic viscosity,μ, kinematics viscosity,n, thermal conductivity,λ, turbulent Prandtl number fork,σk, turbulent Prandtl number forε,σk,T is the temperature,T0is the reference operating temperature.

    The turbulence dissipation rate is denoted byε, whilek is the turbulent kinetic energy of the turbulent fluctuations per unit mass. The turbulent viscosity μtis defined as

    where the coefficient Cμis a new variable defined in the realizable k-εmodel and given by the relation:

    where

    The eddy viscosity formulation is based on the realizability constraints, the positivity of the normal stress and Schwarz inequality for turbulent shear stresses.

    Furthermore, in Eq.(20),C1defined by the form

    and

    is the generation of turbulent kinetic energy due to the mean velocity gradient.

    The velocity componentu parallels to the gravitational vector andv is the component of the velocity perpendicular to the gravitational vector. In this way, C3ε=1for buoyant shear layers in which the main flow direction is aligned with the direction of gravity(the case under study). For buoyant shear layers those are perpendicular to the gravitational vector,C3ε=0.

    The relation gives the generation of turbulence due to buoyancy

    where Prtis the turbulent Prandtl number for energy andβis the thermal expansion coefficient. The model constants of the k-εmodel are established to ensure that the model performs well for certain canonical flows such as pipe flow, jet flow, and boundary layer flow, in this paper C1ε=1.44,C2=1.9,σk= 1.0,σε=1.2and Prt=0.85will be used.

    Shih et al.[17]have compared their model (realizable k-εturbulence model) with experimental data as well as with the standard k-εmodel for a round jet flow and other flows. The comparison shows a good matching between their model and the experimental data than the standard model. The realizable k-ε model implies that the model satisfies specific constraints on the Reynolds' stresses that make the model more consistent with the physics of turbulent flows and hence more accurate than the other turbulent models. This model contains a new transport equation for the turbulent dissipation rate. Also, a critical coefficient of the model,Cμ, is expressed as a function of mean flow and turbulence properties, rather than assumed to be constant as in the standard model. This allows the model to satisfy certain mathematical constraints on the normal stresses consistent with the physics of turbulence (realizability). Additionally, the realizable k-εmodel uses different sources and sinks terms in the transport eddy dissipation. The modified prediction ofεalong with the modified calculation of μt, makes this model superior to the other k-εmodels. For the jet flow this model does better in predicting the spreading rate especially, for near region z<0.35(see El-Amin et al.[1]).

    Boundary conditions need to be specified on all surfaces of the computational domain. Boundaries presented in this study include inflow (Inlet), outflow(outlet), solid wall and axis of symmetry as shown in Fig.2(b). with ?in,?outand ?walldenote the boundary of the inlet, outlet and wall, respectively. In addition to the non-realistic boundary on the axis of symmetry ?axis. The velocity-inlet boundary conditions imposed at the nozzle are

    here Tinis defined in Eqs.(A1)-(A4) for the cases of variable inlet temperature.

    Due to the stagnant conditions of water inside the tank before the beginning of the influx, all velocity components were initially set to zero. Heat transfer through the walls of the tank is not taken into consideration (adiabatic walls). Density of water, specific heat,thermal conductivity and viscosity are formulated as a piecewise-linear profile of temperature. The turbulence intensity and hydraulic diameter characterizes the turbulence at the inlet boundary. The following equation of empirical correlation for pipe flows is used to describe the turbulence intensity at the inlet boundary as a function of the Reynolds number,

    Alternatively, one can use the followingkandεon the inlet boundary as (see Kadem et al.[18])

    whered is the nozzle diameter.

    The boundary condition on the axis of symmetry is represented by free-slip condition, that is, a non-realistic wall with no-friction when velocity and other components near the wall are not retarded. Unlike the no-slip boundary condition for which flow has zero velocity in the wall, free-slip flow is tangent to the surface. On the axis of symmetry, the radial velocity componentv , and the gradient of the other dependent variables (u,ε,k,T ) were equal to zero. So, one may write

    Solid wall boundary conditions are imposed along the solid walls, the no-slip boundary condition for velocities, zero value for turbulent kinetic energy, and zero gradients for temperature and energy dissipation rate were used.

    where n1is the outward normal to the wall.

    Finally, the outlet boundary at which water discharged freely, can be formulated as

    where n2is the outward normal to the outlet boundary.

    Also, static pressure can be defined at a known given point in the domain and Fluent will extrapolate all other conditions from the interior of the domain. Initial conditions are described as follows

    In fact, very small values are given as initial conditions forkandεinstead of zero to speed up convergence of the solution.

    3. Numerical investigations

    Figure 2(b) shows the computational domain with dimensions of: radius is 0.18 m and height is 0.605 m. The radius of the inlet and outlet pipes is 0.01 m, while the inlet height inside the tank is 0.06 m and the outlet depth in the tank is 0.055 m. The meshes are built up of Quadratic Submap cells. The number of grid elements used for all calculations is 7 984. Fluent[19]is used to model the flow in the tank by solving the continuity, momentum, turbulence and energy equations.

    Fig.5 Grid independence test

    Table 2 Time step comparison for temperature, axial velocity, turbulent kinetic energy and turbulent dissipation rate of kinetic energy, for Case 4

    In order to prove grid independence, numerical experiment for Case 4, is repeated on the systematically refined grids in quadrilateral cells of sizes 7 984(grid-1), 9 240 (grid-2), 12 880 (grid-3) and 23 560(grid-4) respectively. The minimum distances between the nodes in the respective grids are 0.00125 m,0.002 m, 0.0014285 m and 0.001 m and the maximum distances between the nodes are 0.0055 m,0.004583 m, 0.00366 m and 0.00275 m respectively in the order of refinement. Figures 5(a), 5(b) show the results of the grid refinement studies for the axial velocity and temperature, respectively. The maximumdeviation caused by grid is about 3% for the velocity,and 0.14 % for the temperature that could be negligible.

    Fig.6 Comparison between measured and simulated temperature as function of tank height for variable inlet temperature (Cases 1-4)

    In order to achieve convergence, under-relaxation is applied to pressure, velocities, energy, turbulent viscosity, turbulence kinetic energy and turbulent dissipation rate calculations. Body force weighted discretization is used for pressure and the velocity-pressure coupling is treated using the SIMPLE algorithm. A second-order upwind scheme is used in the equations of momentum, energy, turbulence kinetic energy and turbulence dissipation rate. Segregated implicit solver with the implicit second-order scheme is used.

    In order to use a suitable time step, we performed a comparison for one case with different time steps as 0.01 s, 0.1 s, 0.5 s and 1 s, shown in Table 2. This comparison includes temperature, axial velocity, turbulent kinetic energy and turbulent dissipation rate of kinetic energy. From this table one can note that the differences are negligible. Then, to reduce the time of calculation we have to use the time step of 1 s.

    Both measured and simulated temperatures as a function of the tank height for various times and variable inlet temperatures (Cases 1-4) are plotted in Figs.6(a)-6(d), respectively. Good agreement between the experimental and numerical data is observed. The maximum error observed is 0.35 K, however, for Cases 1, 2, 3 and 4 the maximum error is 0.2 K, 0.2 K,0.35 K and 0.35 K, respectively, occurs after 30 min of charging process.

    Fig.7 Normalized axial velocity as a function of r/ czat different positions ofz of Case 4 at t=15 min

    3.1 Axial and radial velocities

    The mean positive axial velocity,u(excluding the reflected velocity with the negative values) is normalized with reference to the centerline velocity ucagainst r/ cz(r normalized by the jet width b=cz, cis the lateral spread rate of the jet) with different heights, for the Case 4 at 15 min, is plotted in Fig.7. It can be seen from this figure that axial velocity profile shows self-similar behavior. Therefore, axial velocity may be represented by a Gaussian distribution with centerline velocity,uc, height,z , and width,b, as parameters. The Gaussian function takes the form

    This empirical model is plotted in Fig.8 along with the comparison with the simulated axial velocity. In this study, the parameter of lateral spread rate of the jet is as reported by Fischer et al.[20]. One can note a relatively large error at small velocities at the both ends of the bell-shaped curve. c=0.11which lies in the range of the standard values

    Fig.8 Comparison between simulated and empirical Gaussian model of axial velocity as a function ofr for different positions ofz in the Case 4 at t=15 min

    Using the definition of axial velocity, Eq.(24),centerline axial velocity (velocity on the axis of symmetry) can be given as

    such that uc=u( r→0).

    Alternatively, the centerline velocity may be written in the form

    as to be comparable with the common formula of the centerline velocity given in literature. It is notable that Au=Bud,Buspecifies the decay rate of the time averaged centerline velocity.

    Dimensional arguments together with experimental observations suggest that the mean flow variables,which are known as similarity solutions, are conforming to Eqs.(25) (Fisher et al.[20]). The continuity equation, Eq.(1), for the time-averaged velocities can be solved by substituting the axial velocity form into Eq.(1) to obtain the cross-stream radial velocity in the form

    where

    3.2 Dynamic pressure

    The dynamic pressure behaves similar to the axial velocity but of course it is a scalar quantity so we do not see negative beaks of the curve. The dynamic pressure can be defined according to the equation

    or

    where Pcis the centerline dynamic pressure, and h=c /.

    Fig.9 Comparison between the simulated dynamic pressure and its Gaussian fitting as a function ofr for different positions ofz in the Case 4 at t=15 min

    Figure 9 illustrates a comparison between the simulated dynamic pressure and its Gaussian fitting using Eq.(28) as a function ofr for different positio-ns ofz of the unheated jet at t=5 min(Case 4). This figure shows a good matching for this Gaussian distribution of the dynamic pressure.

    Fig.10 Temperature as a function ofz on the axis of symmetry(r =0)with varying times, Case 1

    Fig.11 Turbulence intensity as a function ofz on the axis of symmetry(r =0)with varying times, Case 1

    Fig.12 Velocity magnitude as a function ofr with varying values of z at t=10 min, Case 3

    3.3 Selected simulated results

    In Fig.10 temperature profiles are plotted against zat different times. One may notice relatively higher temperatures closer to the inlet up to, approximately, z≈0.12 m, and then it decreases aszincreases. As the time proceeds, the temperature closer to the inlet decreases as shown in the figure while it increases further away. The turbulence intensity as a function of the coordinatez along the axis of symmetry for various times is shown in Fig.11. The turbulence intensity is defined as the ratio of the root-mean-square of the turbulent velocity fluctuations and the mean velocity. Apparently closer to the inlet the velocity fluctuations increases due to the impingement of the jet in the relatively quiescent fluid in the tank. However,away from the inlet the intensity of turbulence decreases because of the decrease in velocity as the jet spreads laterally as manifested in Fig.12. For z≤0.14 m the turbulence intensity is the same all the time, while for z>0.14 mm the turbulence intensity decreases with time. It is interesting to note that inside the outlet pipe the turbulence intensity increases as manifested by the sharp increase in turbulence intensity towards the outlet pipe as a result of the influence of pipe wall.

    Fig.13 Temperature as a function ofr with varying values of zat t =10 min, Case 3

    The velocity magnitude as a function of radial axis distance,r , at different positions ofzat t= 10 min is plotted in Fig.12. The velocity magnitude in bottom part of the tank is larger closer to the axis of symmetry while it has smaller values far from it (i.e., asr increases). Aszincreases velocity magnitude decreases closer to axis of symmetryz while it increases asr increases. This behavior may be explained by the fact that the jet leaves the inlet with a higher velocity and disperses laterally as it moves far from the source. Figure 13 shows temperature as a function ofr for various values ofzat t=10 min. At the bottom of the tank (i.e., smallz ), the temperature is higher closer to the symmetry axis and it sharply decreases far from it (i.e., as r increases). Therefore, asz increases and the jet disperses more laterally, the temperature closer to the axis of symmetry decreases while increasing as rincreases.

    Figure 14 displays the contours of static pressure,velocity magnitude, temperature, and turbulence intensity at 1 min, 10 min, 20 min and 30 min. The contours of the stream function are shown in Fig.15. These figures provide a complete understanding of the flow structures under the current conditions especially in the two regions around inlet and outlet.

    Fig.14 Contours of, static pressure (1st row), velocity magnitude (2nd row), temperature (3rd row) and turbulence intensity (4th row) at 1 min (1st column), 5 min (2nd column), 15 min (3rd column) and 25 min (4th column),for Case 4

    Fig.15 Contours of stream function at 30 min for Case 2

    For the bottom part of the tank the flow profiles behave as most of jet flows acts in accordance to the Gaussian distribution with one peak, while the flow seems to be quiet with small velocity around the inlet and it may vanishes under the level of inlet. In the top part of the tank and after impinging the ceiling wall of the tank, three peaks can be seen. One represents the positive main flow inside the jet body and the other two are small close to the wall in the opposite direction (with negative velocity) and it increases as z increases. When the jet impinges the ceiling wall the flow opposites and causes motion from top to bottom close to the walls to form vortex-like structures.

    3.4 Jet Richardson number

    Richardson number is defined as the ratio of buoyancy to inertia forces. But for more convenience we will define the Richardson number according to the local centerline velocity. Richardson number is calculated using buoyancy-related terms (density difference)and the velocity at the same point. In jet flows,Richardson number takes the form[17]

    Richardson number is plotted in Fig.16 against the heightz , at different times for Case 3. From this figure it can be seen that Richardson number is reduced in the region closer to the nozzle inlet, and then it increases with the height. In the bottom part the inertia effect dominates the buoyancy effect (jet-like zone),therefore Richardson number decreases. In the top part of the tank, on the other hand, the buoyancy effect dominates the inertia (plume-like zone) as manifested by the increase of Richardson number. Also, in this zone Richardson number increases with time because temperature increases with time and enhances the buoyancy while it decreases closer to the inlet as the inlet temperature is set to decrease.

    Fig.16 Richardson number as a function of the heightz , at different times, with constant inlet temperature

    In order to examine the effect of varying the inlet temperature on Richardson number we perform three numerical experiments, one of them with constantinlet temperature, and two with monotone increasing and monotone decreasing inlet temperature, respectively. The inlet temperatures for these numerical experiments are defined as:

    Tin=308.9539 K, for the constant inlet temperature,

    T=301.5519+0.01t2, for the monotone in

    increasing inlet temperature,

    T=308.9539-0.01t2, for the monotone de

    increasing inlet temperature, where,1min≤t≤30 min, for monotone increasing inlet temperature,301.552≤Tin≤309.962and for monotone decreasing inlet temperature,308.954≤ Tin≤300.554.

    Fig.17 Richardson number as a function of the heightz , at different times, with monotone increasing inlet temperature

    Fig.18 Richardson number as a function of the heightz , at different times, with monotone decreasing inlet temperature

    Figure 16 shows Richardson number for the case of constant inlet temperature. From this figure it can be seen that Richardson number's behavior is similar for all times closer to the inlet (i.e., in the jet-like region). In the plume-like region Richardson number increases with time because of the increase in temperature. Figures 17 and 18 illustrate Richardson number for the monotone increasing and monotone decreasing inlet temperature, respectively. One can notice that Richardson number in the plume-like region in the case of monotone increasing inlet temperature increases with time as a manifestation of the increased buoyancy, Fig.17. On the other hand, for the monotone decreasing inlet temperature, Fig.18, Richardson number decreases in the Jet-like region as a manifestation of the decreased temperature.

    4. Conclusion

    This work is devoted to investigate the problem of non-uniform inlet temperature of buoyant jet. An analysis for vertical hot water jets entering a cylindrical tank filled with cold water under the condition of variable inlet temperature is introduced. The variable inlet temperature is considered as a function of time of charging process. Experimental measurements are performed for the different cases in sequential time steps for both constant and variable inlet temperature. Numerical investigations under the above mentioned conditions are performed. The realizable k-εturbulence model is used to simulate turbulent flow for this problem. Comparisons between the measured and simulated temperature show good agreements. Selected empirical Gaussian model with standard parameters are used to represent the simulated results. Selected simulated quantities such as velocity magnitude, temperature and turbulence intensity are investigated. The results indicate that temperature varies, approximately,linearly with time for the constant inlet temperature cases, while, it seems to be, approximately, polynomial or logarithms functions of time for the variable inlet temperature, especially, for plume region. Also,thermal stratification exits, however thinner the thermal layers in top part of the tank are than that in the bottom part. The stratification is verified from the beginning of experiment and continues during the whole time.

    Acknowledgements

    The second author would like to thank the Institute of Scientific Research and Revival of Islamic Heritage, Umm Al-Qura University, for the support of the project (Grant No. 43308012).

    [1] EL-AMIN M. F., SUN S. and HEIDEMANN W. et al. Analysis of a turbulent buoyant confined jet modeled using realizable k-εmodel[J]. Heat Mass Transfer, 2010, 46(8): 943-960.

    [2] EL-AMIN M. F., SUN Shuyu and SALAMA Amgad. Simulation of buoyancy-induced turbulent flow from a hot horizontal jet[J]. Journal of Hydrodynamics, 2014,26(1): 104-113.

    [3] PANTHALOOKARAN V., EL-AMIN M. F. and HEI-DEMANN W. et al. Calibrated models for simulation of stratified hot water heat stores[J]. International Journal of Energy Research, 2008, 32(7): 661-676.

    [4] AGRAWAL A., PRASAD A. K. Integral solution for the mean flow profiles of turbulent jets, plumes, and wakes[J]. Journal of Fluids Engineering, 2003, 125(5):813-822.

    [5] O'HERN T. J., WECKMAN E. J. and GERHART A. L. et al. Experimental study of a turbulent buoyant helium plume[J]. Journal of Fluid Mechanics, 2005, 544:143-171.

    [6] ARAKERI J. H., DAS D. and SRINIVASAN J. Bifurcation in a buoyant horizontal laminar jet[J]. Journal of Fluid Mechanics, 2000, 412: 61-73.

    [7] EL-AMIN M. F., KANAYAMA H. Integral solutions for selected turbulent quantities of small-scale hydrogen leakage: A non-buoyant jet or momentum-dominated buoyant jet regime[J]. International Journal of Hydrogen Energy, 2009, 34(3): 1607-1612.

    [8] EL-AMIN M. F., KANAYAMA H. Similarity consideration of the buoyant jet resulting from hydrogen leakage[J]. International Journal of Hydrogen Energy,2009, 34(14): 5803-5809.

    [9] EL-AMIN M. F. Non-Boussinesq turbulent buoyant jet resulting from hydrogen leakage in air[J]. International Journal of Hydrogen Energy, 2009, 34(18): 7873-7882.

    [10] JIRKA G. H. Integral model for turbulent buoyant jets in unbounded stratified flows. Part 1: Single round jet[J]. Environmental Fluid Mechanics, 2004, 4(1): 1-56.

    [11] JIRKA G. H. Integral model for turbulent buoyant jets in unbounded stratified flows, Part 2: Plane jet dynamics resulting from multiport diffuser jets[J]. Environmental Fluid Mechanics, 2006, 6(1): 43-100.

    [12] CHEN Jian-gang, ZHANG Jian-min and XU Wei-lin et al. Particle image velocimetry measurements of vortex structures in stilling basin of multi-horizontal submerged jets[J]. Journal of Hydrodynamics, 2013, 25(4):556-563.

    [13] BARATIAN-GHORGHI Z., KAYE N. B. and KHAN A. A. et al. The merging of two unequal axisymmetric parallel turbulent jets[J]. Journal of Hydrodynamics,2012, 24(2): 257-262.

    [14] DU Yu-kun, WANG Rui-he and NI Hong-jian et al. Dynamical analysis of high-pressure supercritical carbon dioxide jet in well drilling[J]. Journal of Hydrodynamics, 2013, 25(4): 528-534.

    [15] WANG X. L., LEE J. H. and LU T. J. et al. A comparative study of single-/two-jet crossflow heat transfer on a circular cylinder[J]. International Journal of Heat and Mass Transfer, 2014, 78(7): 588-598.

    [16] SHIM Y. M., RICHARDS P. J. and SHARMA R. N. Turbulent structures in the flow field of plane jet impinging on a circular cylinder[J]. Experimental Thermal and Fluid Science, 2014, 57(3): 27-39.

    [17] SHIH T. H., LIOU W. W. and SHABBIR A. et al. A new k-εeddy-viscosity model for high Reynolds number turbulent flows-model development and validation[J]. Computers and Fluids, 1995, 24(3): 227-238.

    [18] KADEM K., MATAOUI A. and SALEM A. et al. Numerical simulation of heat transfer in an axisymmetric turbulent jet impinging on a flat plate[J]. Advanced Modeling and Optimization, 2007, 9(2): 207-217.

    [19] FLUENT INC. Fluent user’s guide[M]. Lebanon, NH,USA: Fluent Inc., 1998.

    [20] FISCHER H. B., LIST E. J. and KOH R. C. Y. et al.

    Mixing in inland and coastal waters[M]. San Diego,USA: Academic Press, 1979.

    Appendixes

    The inlet variable temperature may be given as a function of time for each case as follows:

    These polynomials are plotted in Fig.3(a). The ranges of variation of the inlet temperature are:

    The corresponding Reynolds numbers are:

    These polynomials are shown in Fig.3(b). The ranges of variation of the Reynolds numbers are:

    Also, the nozzle inlet turbulence intensity can be represented in a polynomial form of time:

    These profiles are illustrated in Fig.3(c). The ranges of variation of the inlet turbulence intensities are:

    * Biography: EL-AMIN Mohamed F. (1971-), Male, Ph. D.,Professor

    又黄又粗又硬又大视频| av在线天堂中文字幕| 一本综合久久免费| 久久精品国产亚洲av高清一级| 免费在线观看影片大全网站| 亚洲国产中文字幕在线视频| 男女之事视频高清在线观看| 一进一出抽搐gif免费好疼| 欧美中文综合在线视频| 夜夜爽天天搞| av欧美777| 久久国产精品影院| 亚洲avbb在线观看| 欧美老熟妇乱子伦牲交| 在线十欧美十亚洲十日本专区| 十分钟在线观看高清视频www| 88av欧美| 精品国产亚洲在线| 国内精品久久久久久久电影| 88av欧美| 999精品在线视频| 亚洲情色 制服丝袜| 99riav亚洲国产免费| 久久亚洲精品不卡| 欧美不卡视频在线免费观看 | 午夜久久久久精精品| 两个人视频免费观看高清| 美女午夜性视频免费| 757午夜福利合集在线观看| 最新在线观看一区二区三区| 国产一区二区三区综合在线观看| 久久久久久免费高清国产稀缺| 中文字幕最新亚洲高清| 欧美色视频一区免费| 日韩欧美国产一区二区入口| 99精品在免费线老司机午夜| 久久久久国内视频| 久久精品国产亚洲av香蕉五月| 国产精品一区二区精品视频观看| 久久伊人香网站| 成人手机av| 午夜福利欧美成人| 国产精品爽爽va在线观看网站 | 十八禁网站免费在线| 狠狠狠狠99中文字幕| 亚洲专区字幕在线| 国产成人一区二区三区免费视频网站| 久久久精品国产亚洲av高清涩受| 狂野欧美激情性xxxx| 午夜福利欧美成人| 久久久国产成人精品二区| 18禁观看日本| 国产成人啪精品午夜网站| 国产私拍福利视频在线观看| 欧美一级a爱片免费观看看 | 国产色视频综合| 国产真人三级小视频在线观看| 日日夜夜操网爽| 我的亚洲天堂| 青草久久国产| 国产午夜精品久久久久久| 动漫黄色视频在线观看| 久久精品国产99精品国产亚洲性色 | 国产亚洲av嫩草精品影院| 久久久久久久久中文| 欧美成人性av电影在线观看| 天天一区二区日本电影三级 | 亚洲精品美女久久av网站| 性欧美人与动物交配| 国产av在哪里看| 18禁黄网站禁片午夜丰满| 黄色女人牲交| 日韩欧美国产一区二区入口| 久久午夜综合久久蜜桃| 亚洲国产精品sss在线观看| 亚洲欧美激情综合另类| 自线自在国产av| 久久这里只有精品19| 久久精品人人爽人人爽视色| 波多野结衣高清无吗| 国产激情欧美一区二区| 免费在线观看亚洲国产| 欧美在线一区亚洲| 日韩 欧美 亚洲 中文字幕| 在线观看www视频免费| svipshipincom国产片| 国产熟女午夜一区二区三区| 国内毛片毛片毛片毛片毛片| 午夜免费观看网址| 天天躁狠狠躁夜夜躁狠狠躁| 侵犯人妻中文字幕一二三四区| 精品国产乱子伦一区二区三区| 热99re8久久精品国产| 亚洲av熟女| 久久香蕉国产精品| 午夜福利欧美成人| 一区二区三区高清视频在线| 女人爽到高潮嗷嗷叫在线视频| 国产人伦9x9x在线观看| 精品欧美国产一区二区三| av网站免费在线观看视频| 88av欧美| 免费av毛片视频| 一区二区三区国产精品乱码| 日本a在线网址| 动漫黄色视频在线观看| 国产麻豆成人av免费视频| 国产亚洲欧美在线一区二区| 国产又爽黄色视频| 日韩国内少妇激情av| 老司机福利观看| 欧美人与性动交α欧美精品济南到| 99国产综合亚洲精品| 国产欧美日韩一区二区三| 免费无遮挡裸体视频| 国产成人欧美在线观看| 两个人视频免费观看高清| 亚洲国产精品sss在线观看| 岛国在线观看网站| 99久久综合精品五月天人人| 中文亚洲av片在线观看爽| 国产伦人伦偷精品视频| 变态另类成人亚洲欧美熟女 | 久久国产精品影院| 两性夫妻黄色片| 亚洲精品美女久久av网站| 咕卡用的链子| 成年女人毛片免费观看观看9| 欧美激情极品国产一区二区三区| 精品国产超薄肉色丝袜足j| 丝袜美足系列| 国产激情久久老熟女| 国产精品久久电影中文字幕| 精品欧美国产一区二区三| 国产私拍福利视频在线观看| 日本免费一区二区三区高清不卡 | 在线观看66精品国产| 国产aⅴ精品一区二区三区波| 人人妻人人澡人人看| 最近最新中文字幕大全免费视频| 国产成年人精品一区二区| 国产av在哪里看| 9热在线视频观看99| 亚洲色图综合在线观看| 国产极品粉嫩免费观看在线| 中文字幕另类日韩欧美亚洲嫩草| 久久精品影院6| 丰满的人妻完整版| 巨乳人妻的诱惑在线观看| 国产高清视频在线播放一区| 美国免费a级毛片| 国产精品 欧美亚洲| bbb黄色大片| 亚洲七黄色美女视频| 国产91精品成人一区二区三区| 97碰自拍视频| 精品不卡国产一区二区三区| 好看av亚洲va欧美ⅴa在| 久久久久国产精品人妻aⅴ院| 国产精品一区二区在线不卡| 日本撒尿小便嘘嘘汇集6| 久久久久九九精品影院| 久久中文字幕一级| 啦啦啦免费观看视频1| 亚洲五月色婷婷综合| 国产精品野战在线观看| 美女扒开内裤让男人捅视频| 日本vs欧美在线观看视频| 亚洲七黄色美女视频| 久久精品国产亚洲av香蕉五月| 国产男靠女视频免费网站| 精品不卡国产一区二区三区| 看片在线看免费视频| 高清在线国产一区| 亚洲第一欧美日韩一区二区三区| 久久性视频一级片| 黄片小视频在线播放| 99国产极品粉嫩在线观看| 国产亚洲av高清不卡| 久久影院123| www.www免费av| 亚洲少妇的诱惑av| 免费观看精品视频网站| 欧美日韩福利视频一区二区| 国产免费av片在线观看野外av| 亚洲国产欧美网| 人人妻人人澡人人看| 在线天堂中文资源库| 身体一侧抽搐| 黄色毛片三级朝国网站| 欧美成狂野欧美在线观看| 国产私拍福利视频在线观看| 亚洲精品在线美女| 欧美另类亚洲清纯唯美| 久久午夜亚洲精品久久| 成年人黄色毛片网站| 香蕉国产在线看| 久久久国产成人免费| 禁无遮挡网站| 最近最新免费中文字幕在线| 男女下面插进去视频免费观看| 淫秽高清视频在线观看| 99国产综合亚洲精品| 精品熟女少妇八av免费久了| 午夜免费鲁丝| 黄网站色视频无遮挡免费观看| 天天一区二区日本电影三级 | 久久欧美精品欧美久久欧美| 亚洲电影在线观看av| 国产视频一区二区在线看| 欧美亚洲日本最大视频资源| 亚洲精品一卡2卡三卡4卡5卡| 免费无遮挡裸体视频| 欧美激情久久久久久爽电影 | 久久久久久大精品| 国产一区二区三区综合在线观看| 久久精品国产清高在天天线| 搡老岳熟女国产| 午夜a级毛片| 日韩欧美免费精品| videosex国产| 国产精品av久久久久免费| 免费观看人在逋| 欧美成人性av电影在线观看| 精品免费久久久久久久清纯| 69av精品久久久久久| 黄色视频,在线免费观看| 国产精品自产拍在线观看55亚洲| 亚洲性夜色夜夜综合| 免费看美女性在线毛片视频| 久久中文字幕一级| 高清黄色对白视频在线免费看| 亚洲国产欧美日韩在线播放| 亚洲天堂国产精品一区在线| 免费看十八禁软件| 欧美在线一区亚洲| 老鸭窝网址在线观看| 婷婷六月久久综合丁香| 激情在线观看视频在线高清| 最近最新免费中文字幕在线| 中文字幕久久专区| 国产精品久久久久久精品电影 | 日本 欧美在线| 亚洲视频免费观看视频| 日韩欧美一区二区三区在线观看| 久久欧美精品欧美久久欧美| 一进一出抽搐gif免费好疼| 日韩精品中文字幕看吧| 日日爽夜夜爽网站| 一级毛片女人18水好多| 日韩欧美国产在线观看| 女人精品久久久久毛片| 欧美黄色淫秽网站| 久久久国产成人精品二区| 男人的好看免费观看在线视频 | 精品久久蜜臀av无| 欧美人与性动交α欧美精品济南到| 久99久视频精品免费| 亚洲五月天丁香| 国产一区二区三区在线臀色熟女| 1024香蕉在线观看| 久久精品aⅴ一区二区三区四区| 久久狼人影院| 黄色毛片三级朝国网站| svipshipincom国产片| 亚洲avbb在线观看| 亚洲av熟女| 村上凉子中文字幕在线| 777久久人妻少妇嫩草av网站| 9色porny在线观看| 麻豆国产av国片精品| 午夜精品在线福利| 成熟少妇高潮喷水视频| 在线av久久热| 欧美人与性动交α欧美精品济南到| 日日干狠狠操夜夜爽| 亚洲人成77777在线视频| 91成年电影在线观看| 嫩草影院精品99| 淫秽高清视频在线观看| cao死你这个sao货| av网站免费在线观看视频| 1024视频免费在线观看| 大型av网站在线播放| 亚洲国产精品成人综合色| 亚洲精华国产精华精| 男人的好看免费观看在线视频 | 久久精品国产综合久久久| 色综合婷婷激情| 久久精品国产综合久久久| 超碰成人久久| 韩国av一区二区三区四区| 日韩欧美在线二视频| 免费久久久久久久精品成人欧美视频| 亚洲欧美日韩高清在线视频| 中文字幕av电影在线播放| 亚洲一区二区三区不卡视频| 亚洲av五月六月丁香网| 丝袜在线中文字幕| 97人妻精品一区二区三区麻豆 | 亚洲专区字幕在线| 不卡av一区二区三区| 不卡av一区二区三区| 久久久久久免费高清国产稀缺| 18禁观看日本| 久久狼人影院| 啪啪无遮挡十八禁网站| 国产色视频综合| 欧美丝袜亚洲另类 | 好看av亚洲va欧美ⅴa在| 曰老女人黄片| 激情视频va一区二区三区| 亚洲精品在线美女| 无遮挡黄片免费观看| 国产av一区二区精品久久| 日韩视频一区二区在线观看| 久久精品91无色码中文字幕| 亚洲成人国产一区在线观看| 国产精品爽爽va在线观看网站 | 少妇粗大呻吟视频| 成人国产一区最新在线观看| 亚洲自拍偷在线| 国产精品,欧美在线| 久久久国产精品麻豆| 久久精品成人免费网站| 可以在线观看的亚洲视频| 一进一出抽搐动态| 国产精品久久电影中文字幕| 50天的宝宝边吃奶边哭怎么回事| 国产免费男女视频| 每晚都被弄得嗷嗷叫到高潮| 国产又色又爽无遮挡免费看| 美女高潮到喷水免费观看| 久久中文字幕人妻熟女| 中国美女看黄片| 久久这里只有精品19| 国产激情久久老熟女| 亚洲中文日韩欧美视频| 老熟妇乱子伦视频在线观看| 动漫黄色视频在线观看| 国产欧美日韩一区二区三区在线| 一a级毛片在线观看| 黄色视频不卡| av超薄肉色丝袜交足视频| 一进一出抽搐动态| 9色porny在线观看| 老司机在亚洲福利影院| 久久婷婷人人爽人人干人人爱 | 多毛熟女@视频| 久久午夜综合久久蜜桃| 亚洲av片天天在线观看| 日韩精品中文字幕看吧| 亚洲一区中文字幕在线| 亚洲中文字幕日韩| 久久国产精品男人的天堂亚洲| 久久香蕉国产精品| 大陆偷拍与自拍| 欧美日韩亚洲综合一区二区三区_| 国产精品日韩av在线免费观看 | 午夜福利视频1000在线观看 | 国产av一区在线观看免费| 精品欧美一区二区三区在线| 国产成人av教育| 日韩有码中文字幕| 在线观看日韩欧美| 18禁黄网站禁片午夜丰满| 成年人黄色毛片网站| 制服诱惑二区| 免费搜索国产男女视频| 91精品三级在线观看| 成人三级黄色视频| 午夜精品国产一区二区电影| 大型黄色视频在线免费观看| 久久欧美精品欧美久久欧美| 欧美一区二区精品小视频在线| 亚洲中文字幕一区二区三区有码在线看 | 少妇粗大呻吟视频| 欧美黑人欧美精品刺激| 国产成人影院久久av| 精品久久久久久成人av| 国产成人av激情在线播放| 不卡一级毛片| 欧美黄色片欧美黄色片| 在线观看免费日韩欧美大片| 成人18禁在线播放| 久久精品亚洲精品国产色婷小说| 久久国产精品影院| avwww免费| 啪啪无遮挡十八禁网站| 欧美中文综合在线视频| 老熟妇乱子伦视频在线观看| 国产三级在线视频| 精品久久久久久久毛片微露脸| 青草久久国产| 精品一区二区三区av网在线观看| 欧美精品亚洲一区二区| 男人舔女人下体高潮全视频| 91av网站免费观看| 精品久久蜜臀av无| 曰老女人黄片| 他把我摸到了高潮在线观看| 亚洲七黄色美女视频| 国产精品亚洲av一区麻豆| 国产激情欧美一区二区| 精品一区二区三区视频在线观看免费| 亚洲av美国av| 首页视频小说图片口味搜索| 国产精品1区2区在线观看.| 日韩 欧美 亚洲 中文字幕| 夜夜爽天天搞| 9色porny在线观看| 一区二区三区激情视频| 99热只有精品国产| 欧美成狂野欧美在线观看| 窝窝影院91人妻| 波多野结衣高清无吗| 亚洲av日韩精品久久久久久密| 一级作爱视频免费观看| xxx96com| 国产精品1区2区在线观看.| 久热这里只有精品99| 亚洲国产精品成人综合色| 美女扒开内裤让男人捅视频| 欧美不卡视频在线免费观看 | 久久中文字幕人妻熟女| 欧美成人性av电影在线观看| 黄片大片在线免费观看| 久久香蕉激情| 欧美色欧美亚洲另类二区 | 曰老女人黄片| 身体一侧抽搐| 日本免费一区二区三区高清不卡 | 十八禁网站免费在线| 欧美日韩中文字幕国产精品一区二区三区 | 日日摸夜夜添夜夜添小说| 欧美日本视频| 九色亚洲精品在线播放| 午夜影院日韩av| 亚洲成av片中文字幕在线观看| 免费观看人在逋| 国产av又大| 国产精品永久免费网站| 国产成人免费无遮挡视频| 97碰自拍视频| 一区福利在线观看| 91国产中文字幕| 免费一级毛片在线播放高清视频 | 国产精品 国内视频| 黄色a级毛片大全视频| 久久久久亚洲av毛片大全| 天天添夜夜摸| √禁漫天堂资源中文www| 丝袜在线中文字幕| 99国产综合亚洲精品| 久久伊人香网站| 一区二区三区高清视频在线| 精品免费久久久久久久清纯| 精品久久久久久成人av| 身体一侧抽搐| 欧美日韩乱码在线| 亚洲国产精品久久男人天堂| 日韩免费av在线播放| 满18在线观看网站| 正在播放国产对白刺激| 黄网站色视频无遮挡免费观看| 亚洲熟妇熟女久久| 午夜精品久久久久久毛片777| 午夜成年电影在线免费观看| 国产成人精品久久二区二区免费| 无限看片的www在线观看| 一卡2卡三卡四卡精品乱码亚洲| 色在线成人网| 激情在线观看视频在线高清| 日韩欧美三级三区| 久9热在线精品视频| 欧美乱色亚洲激情| 国产精品,欧美在线| 热99re8久久精品国产| 久久精品亚洲熟妇少妇任你| 色综合亚洲欧美另类图片| 国产熟女午夜一区二区三区| 色播在线永久视频| 久久亚洲精品不卡| 日韩有码中文字幕| 女性被躁到高潮视频| 免费av毛片视频| 99在线视频只有这里精品首页| 成人永久免费在线观看视频| 色综合亚洲欧美另类图片| av福利片在线| 午夜福利影视在线免费观看| 两个人免费观看高清视频| av福利片在线| 色老头精品视频在线观看| 宅男免费午夜| 操美女的视频在线观看| 麻豆成人av在线观看| 日韩三级视频一区二区三区| 91麻豆av在线| 亚洲 欧美 日韩 在线 免费| 国产精品爽爽va在线观看网站 | 亚洲国产精品999在线| 亚洲人成电影观看| 热re99久久国产66热| 日韩欧美一区视频在线观看| 男人操女人黄网站| 两性午夜刺激爽爽歪歪视频在线观看 | 午夜福利一区二区在线看| 99在线视频只有这里精品首页| 啦啦啦韩国在线观看视频| 亚洲国产毛片av蜜桃av| 国产色视频综合| 曰老女人黄片| videosex国产| 久久久久久久久免费视频了| 亚洲国产精品999在线| 一级,二级,三级黄色视频| 国产成人av激情在线播放| 波多野结衣高清无吗| 欧美在线一区亚洲| 成人精品一区二区免费| 日韩欧美一区视频在线观看| 精品国产一区二区久久| 久久香蕉精品热| 中文字幕人妻丝袜一区二区| 老司机午夜福利在线观看视频| 日韩中文字幕欧美一区二区| 欧美成人性av电影在线观看| 丁香六月欧美| 宅男免费午夜| 亚洲欧美精品综合久久99| 老汉色∧v一级毛片| 中国美女看黄片| 熟妇人妻久久中文字幕3abv| 国产成年人精品一区二区| 欧美另类亚洲清纯唯美| 亚洲自拍偷在线| 欧美乱妇无乱码| 又紧又爽又黄一区二区| 日韩三级视频一区二区三区| 亚洲avbb在线观看| 亚洲美女黄片视频| 给我免费播放毛片高清在线观看| 99国产精品一区二区三区| 90打野战视频偷拍视频| 久久热在线av| 精品人妻在线不人妻| 夜夜夜夜夜久久久久| 久久久久久人人人人人| 亚洲人成电影观看| 一边摸一边抽搐一进一出视频| 免费看美女性在线毛片视频| 精品国内亚洲2022精品成人| 深夜精品福利| 亚洲精品国产色婷婷电影| 午夜a级毛片| www国产在线视频色| 人人妻人人澡欧美一区二区 | 无限看片的www在线观看| АⅤ资源中文在线天堂| 我的亚洲天堂| 午夜福利18| 欧美色欧美亚洲另类二区 | 青草久久国产| 麻豆成人av在线观看| 亚洲七黄色美女视频| 在线观看舔阴道视频| 成人亚洲精品一区在线观看| 国产日韩一区二区三区精品不卡| 在线观看免费视频网站a站| 精品久久蜜臀av无| 美女大奶头视频| 成人国语在线视频| 亚洲精品在线观看二区| tocl精华| 日本黄色视频三级网站网址| 中亚洲国语对白在线视频| 女性被躁到高潮视频| 免费女性裸体啪啪无遮挡网站| 欧美久久黑人一区二区| 宅男免费午夜| 九色国产91popny在线| 久久中文看片网| av天堂久久9| 亚洲精品久久成人aⅴ小说| 校园春色视频在线观看| 男女下面进入的视频免费午夜 | 国产私拍福利视频在线观看| 欧洲精品卡2卡3卡4卡5卡区| 国产一区在线观看成人免费| 久久久国产成人精品二区| 久久婷婷成人综合色麻豆| 亚洲一卡2卡3卡4卡5卡精品中文| 欧美日韩黄片免| 在线观看日韩欧美| 妹子高潮喷水视频| 久久性视频一级片| 999久久久精品免费观看国产| 亚洲五月婷婷丁香| 亚洲电影在线观看av| 男人操女人黄网站| 色哟哟哟哟哟哟| 男人操女人黄网站| 午夜视频精品福利| 久99久视频精品免费| 成人国产一区最新在线观看| 国产欧美日韩综合在线一区二区| 麻豆久久精品国产亚洲av| 久久久久久久精品吃奶| 成人欧美大片| av在线播放免费不卡| 国产精品美女特级片免费视频播放器 | 在线观看66精品国产| 欧美日韩中文字幕国产精品一区二区三区 | 国产成年人精品一区二区|