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

    Numerical Modeling and Analysis of Gas Entrainment for the Ventilated Cavity in Vertical Pipe*

    2014-07-24 15:40:13XIANGMinJIANGZhenyu江振宇ZHANGWeihua張為華andTUJiyuan屠基元
    關(guān)鍵詞:基元

    XIANG Min (向 敏), JIANG Zhenyu (江振宇), ZHANG Weihua (張為華)and TU Jiyuan (屠基元)

    1Collage of Aerospace Science and Engineering, National University of Defense Technology, Changsha 410073, China

    2Institute of Nuclear and New Energy Technology, Tsinghua University, Beijing 100084, China

    3School of Aerospace, Mechanical and Manufacturing Engineering, RMIT University, Victoria 3083, Australia

    Numerical Modeling and Analysis of Gas Entrainment for the Ventilated Cavity in Vertical Pipe*

    XIANG Min (向 敏)1,2,**, JIANG Zhenyu (江振宇)1, ZHANG Weihua (張為華)1and TU Jiyuan (屠基元)2,3

    1Collage of Aerospace Science and Engineering, National University of Defense Technology, Changsha 410073, China

    2Institute of Nuclear and New Energy Technology, Tsinghua University, Beijing 100084, China

    3School of Aerospace, Mechanical and Manufacturing Engineering, RMIT University, Victoria 3083, Australia

    A semi-empirical gas entrainment model was proposed for the ventilated cavity in vertical pipe, based on which, a complete numerical scheme was established by coupling with the Eulerian-Eulerian two-fluid model to predict the multiscale flow field created by ventilated cavity. Model predictions were validated against experimental measurements on void fraction and bubble size distributions. Simulations were carried out to explore the effect of ventilation rate and inlet turbulence intensity on the macroscale cavity shape and the bubbly flow downstream of the ventilated cavity. As the ventilation rate increasing, a reverse trend was observed for the void fraction and bubble size distributions. It is concluded that the average void fraction in the pipe flow region is determined by the volumetric ratio between liquid and gas. However, the bubble size evolution is dominated by the breakage effect induced by turbulence in the vortex region. Furthermore, simulations were conducted to analyze geometric scale effect based upon Froude similitude. The results imply that the velocity distributions were properly scaled. Slight scale effect was seen for the void fraction caused by faster dispersion of bubbles in the larger size model. The comparatively greater bubble size was predicted in the smaller model, implying significant scale effects in terms of turbulence and surface tension effect. It reveals that empirical correlations valid in wide range are required for the extrapolation from small-size laboratory models.

    ventilated cavity, gas entrainment, bubbly flow, simulation

    1 INTRODUCTION

    The flow phenomenon of a ventilated cavity is frequently encountered in practical situations such as U-tube fermenters and airlift bioreactors [1]. The ventilated cavity refers generally to a large Taylor bubble attached to a gas-injecting nozzle in a downward flow pipe as depicted in Fig. 1. Gas entrainment, by which gas is introduced into the liquid in the form of gas bubbles due to the action of turbulence and surface tension, occurs at the bottom of the ventilated cavity. The bubbly flow region containing numerous and diverse bubbles is generated downstream of the ventilated cavity which dominates the mass and momentum transfer characteristics between gas and liquid. Three dynamical multiphase phenomena are coupled: (1) macroscale cavity formation, (2) gas entrapment and breakup at the cavity bottom, and (3) dispersion of bubbles and its effect on the liquid phase turbulence, posing a tremendous challenge in comprehending this kind of multiscale two-phase flow.

    Numerous studies have been undertaken to better understand the basic phenomenas involved in the gas entrainment process. Some research [1-3] was focused on establishing relationships for the cavity length and bubble entrainment rate and correlating them with gas and liquid flow rates, and concluded that the gas entrainment rate was affected by the turbulence intensity, jet velocity, liquid surface tension and viscosity etc. Other paper [4-6] concentrated on the flow patterns in the bubbly flow region generated by gas entrainment. In particular, Su [4] observed three different regions including the vortex region, transitional region and pipe flow region behind the cavity as shown in Fig. 1, validating the complexity and diversity of the generated bubbly flow. Most of the previous experimental research was carried out in the vertical pipe with diametersranging between 50 mm and 100 mm, providing primary knowledge on the gas entrainment process. However, considering the large number of relevant parameters affecting the gas entrainment process including fluid properties, pipe geometry, inflow conditions, and entrained air bubble characteristics, it is hard to perform parameter sensitivity study for this phenomenon.

    Figure 1 Schematic of the flow field below ventilated cavity [5]

    Furthermore, due to the multiscale characteristics imbedded in this problem, an extreme difficulty is induced in carrying out similarity analysis of this two-phase flow. The Froude number similitude [7, 8] is needed for the macroscale cavity while a Weber similitude is required for the entrapment of air bubbles which is dominated by surface tension effects. Furthermore, the momentum exchanges between cavity and streaming film flow which is dominated by viscous effects implies the need for a Reynolds similitude. For geometrically-similar models, it is impossible to satisfy simultaneously Froude, Reynolds and Weber similarities. Therefore, the gas entrainment process for a ventilated cavity may be affected by significant scale effects.

    To date, little research has been conducted at the parameter sensitivity study on the gas entrainment for ventilated cavity due to the limit knowledge on the microscale nature of the slug flow. In our previous study [9], a numerical model, based on the Eulerian-Eulerian two-fluid model coupled with the population balance model, has been successfully established to predict the turbulent two-phase flow structure as well as the bubble size distribution downstream of a ventilated cavity. In this study, the previous work is extended by developing a gas entrainment model based on the analysis of the gas entrainment mechanisms. A complete numerical scheme which is capable of predicting the gas entrainment process for a wide range of flow situations is then established by integrating the gas entrainment model into the two-fluid model. The proposed models are validated against the experimental work carried out by Su [4] and Delfos et al [5]. Particular emphasis is then directed towards investigating the effect of ventilation condition, inlet turbulence intensity and geometric scale acting on the bubbly flow downstream of the ventilated cavity. The numerical results will provide us deeper information on the flow field generated by gas entrainment, helping us to condense valuable insights to enhance this two-phase system performance.

    2 SEMI-EMPIRICAL GAS ENTRAINMENT MODEL

    According to Chanson [10], gas entrainment at a ventilated cavity tail is governed by the free surface gas entrainment caused by plunging jet. It is well established that when the water jet velocity Ujat the cavity base is less than a certain threshold velocity Ue, no gas entrainment takes place. When it exceeds the onset velocity, the gas entrainment is clarified into two regimes. (1) When the jet velocity is small (lower than 4-6 m·s?1), gas may be captured in the induction trumpet caused by the surface roughness [11-13] as shown in Fig. 2. Subsequently, the entrained gas pockets are broken up into small bubbles which are then dragged into the liquid bulk. (2) When the water film jet velocity increases further, the gas entrainment is accomplished dominantly by the boundary layer gas entrainment mechanism [14, 15]. Presently, we restrict ourselves to the first regime in this paper.

    Figure 2 Schematic of gas entrainment at the cavity base

    The volumetric rate of gas entrainment per unit peripheral length at the cavity base is considered to be proportional to the jet velocity Ujand the surface roughness hw[16]:

    The surface roughness is proved to be proportional to the kinetic of the liquid jet [11, 12, 17]:

    According to the experimental results [16], an intermittence factor Iwhas to be introduced to consider the probability of the surface fluctuation. Based on the assumption that the increase of intermittency on the surface of the falling film shows a similar behavior as the growth of turbulent spots in a boundary layer [18], Iwfor ventilated cavity in vertical pipe is calculated as the following formula:

    where Lcrefers to the cavity length, Lonis the cavity length for the onset of gas entrainment which is calculated according to the onset velocity for gas entrainment. This formula was firstly validated in a vertical pipe with 100 mm diameter and the cavity length vary from 1D to 15D. In our research, this formula is used in small diameter pipe ranging from 50 mm to 200 mm, and the cavity length in the pipe varies from 1D to 20D. The final volumetric gas entrainment rateis expressed as

    where Dcrepresents the cavity base diameter, and C1is an adjustable constant according to the experimental condition [12]. In this paper, a value of 2.5×10?4was adopted for C1which is in accordance with the plunging jet. The water jet velocity at the cavity tail is obtained according to the mass conservation:

    where Qlis the inlet liquid volumetric rate, and K is a proportionality constant with a value of 1.2 adopted [19].

    It is imperative that the cavity boundary is appropriately defined in order to calculate the cavity diameter and jet velocity at the cavity base. The cavity profile proposed by Dumitrescu [20] is used in this study:

    where z refers to the axial distance from the cavity top where the cavity diameter is zero, r represents the radial coordinate and R is the pipe radius. In the case of a long confined cavity, if the gravity acting on the film is balanced by the shear force, the film thickness will reach a constant value which can be calculated from the following correlation given by Henstock and Hanratty (1979):

    where Ulis the liquid velocity at the inlet,lν is the kinematic viscosity and D represents the pipe diameter.

    Figure 3 Comparison between predicted gas entrainment rate with experimental data

    Due to the limited information on the size of the inception bubbles, the entrained bubble size is taken to be the maximum stable bubble diameter rmaxin the mixing zone behind the cavity. The bubble diameter is calculated using the following formula which was validated through comparison with experiments carried out in a vertical pipe with 50 mm diameter and the inlet liquid velocity range from 0.5 to 1 m·s?1[21]:

    where Wecis the critical Weber number for breakage with a value of 4.7 adopted. σ refers to the surface tension coefficient and Rcis the cavity radius.

    The experimental data on vertical pipes with different size by Su [4] and Delfos et al. [5] were employed to assess the model predictions. Considering the similar low turbulence level existing in the experiments, a constant onset velocity of 1 m·s?1for gas entrainment was adopted according to empirical formula for the vertical circular plunging jet [22]. Given different liquid rates, the variation of gas entrainment rate along with cavity length was obtained as shown in Fig. 3. It was observed that higher inlet liquid rate induced shorter cavity length for the onset of gas entrainment. When the gas entrainment was inspired, the gas entrainment rate climbed up slowly at the initial stage due to the low probability of the surface fluctuation. Subsequently, the gas entrainment rate increased linearly along with the cavity length until reaching the maximum value gradually. The gas entrainment model is validated by good agreement with the experimental results.

    3 TWO-FLUID MODEL

    3.1 Mass and momentum conservation equations

    The two-fluid model based on the Eulerian-Eulerianframework [23] solves the ensemble-averaged of mass, momentum and energy equations for each phase, whereby the liquid is considered as the continuum phase and the gas as disperse phase. Interactions between phases are effected via interfacial transfer terms for heat, mass and momentum exchange. Since there is no interfacial mass or heat transfer between the phases in the present study, the energy equation is not needed.

    In the absence of interfacial mass transfer, the continuity equation of the two-phases can be written as where α, ρ and u are the void fraction, density and velocity vector of each phase. Subscripts i=l or g denote the liquid and gas phase respectively. The momentum equation can be expressed as:

    where ρlis adopted as the reference density ρrefto calculate the buoyancy force. Firepresents the interfacial forces for the interfacial momentum transfer and g is the gravity acceleration vector. It is noted that the interfacial forces appearing in the momentum equation strongly govern the distribution of the liquid and gas phases within the flow volume. Details of the different forces acting within the two-phase flow can be found in our previous work [9].

    3.2 Population balance model

    The population balance method is adopted to predict the size distribution of the poly-dispersed bubbles. The population balance equation is solved with application of the MUSIG model [9] which employs multiple discrete bubble size groups to represent the population balance of bubbles. Assuming each bubble class travel at the same mean algebraic velocity, individual number density of bubble class k can be expressed as [24]

    where nkis the average bubble number density of the kth group, the source terms BC, BB, DCand DBare the birth rates due to coalescence and break-up, and the death rate due coalescence and break-up of bubbles respectively. Details of the break-up and coalescence rate can be found in our previous work [9].

    3.3 Turbulence model

    Turbulence model is required to close the terms of effective viscosity in the Navier-Stokes equations. In this paper, separate turbulence model is adopted for gas and liquid. The effective viscosity of the liquid phase is considered as being composed of the molecular viscosity μl, the turbulent viscosity μt,land the bubble induced turbulent viscosity μt,b:

    The two-equation k-ε model is applied for the liquid turbulent viscosity, as it offers a good compromise between numerical effort and computational accuracy, giving:

    where k represents the turbulent kinetic energy of liquid and ε describes its dissipation rate, Cμis a model constant, k and ε are calculated from their transport equations [25].

    The characteristics of the gas phase turbulence are deduced from the ones of the liquid phase as follows:

    The extra bubble induced turbulent viscosity in Eq. (11) is evaluated according to the model by Sato et al. [26]:

    where Dsrepresents the Sauter mean bubble diameter and is calculated according to the diameter and number density of every bubble class:

    4 EXPERIMENTAL AND NUMERICAL DETAILS

    The experiments [4] carried out in a 3 m long vertical pipe with an internal diameter of 0.058 m by Su [4] are used for further investigation. A gas injection pipe was used to produce a metered flow of gas into the down flowing liquid stream to form a large attached cavity. The superficial liquid velocity Ulvaries from 0.27 to 0.9 m·s?1and the superficial gas velocity Ugfrom 0.01 to 0.1 m·s?1to form cavity with different length. Three cases as depicted in Table 1 are selected as the basis of simulation, where Qland Qgrefer to the liquid flow rate and gas ventilation rate respectively. Fr is the Froude number defined by inlet liquid velocity Uland pipe diameter D. Re is the Reynolds number defined by the liquid inlet velocity and pipe diameter.

    The generic computational fluid dynamics (CFD) code ANSYS CFX12 was utilized to investigate the liquid flow around the ventilated cavity and the two-phase flow behavior downstream of the cavity. Numerical calculations were performed on a 45° radialsector of the pipe with symmetry boundary conditions imposed at both vertical sides of the pipe. Taking the liquid inlet as the reference datum, the outlet boundary was placed at a distance of 1.7 m downstream. The ventilated cavity generated at the injection pipe was placed 100 mm downstream of the pipe inlet. The cavity length and boundary profile which were calculated through the air entrainment model and assumed to be fixed during the simulation. Therefore, the gas field inside the cavity was excluded from the computational domain. A velocity inlet is set up at the pipe inlet, while a mass flow inlet boundary condition is specified at the circumference of the cavity base for the entrained bubbles. A low turbulence intensity of 1% based on velocity magnitude was specified at the liquid inlet. The bubble inlet thickness is determined based on the initial bubble diameter which is estimated using Eq. (7).

    Table 1 Experimental cases description

    Figure 4 Visualization and mesh distribution of computational model

    Figure 5 Water velocity vector distribution for Case 2

    Figure 6 Dimensionless water velocity distribution downstream cavity base for Case 2 (z: axial distance from the cavity base)

    The mesh structure and boundary conditions were presented in Fig. 4. In order to carry out mesh sensitivity study, three levels of computational mesh consisting of around 30000 nodes (coarse), 60000 nodes (medium) and 120000 nodes (fine) were adopted. Fig. 5 gives the water vector distribution throughout the computational domain. The predicted water velocity distribution downstream of the cavity base for Case 2 using different mesh size were compared in Fig. 6 where the maximum differences between the mediumand fine mesh were less than 3%, indicating the numerical prediction was mesh independent. Therefore, the medium mesh with the maximum grid size taken as 4 mm and the minimum grid size taken as 0.6 mm was adopted for the following research. The transport equations of the two-fluid and population balance models were discretized by the finite volume approach. The convection terms were approximated by the high resolution scheme which uses a slope limiter to limit the gradient around shocks or discontinuities, while the diffusion terms were approximated by the second-order central difference scheme. These different schemes are detailed in the technical documentation for CFX 11 [27]. A fixed physical time scale of 0.001 s was adopted to achieve steady state solutions.

    5 RESULTS

    5.1 Influence of ventilation rate on the multiscale flow field

    Figure 7 shows the cavity shape and the void fraction distributions for different cases. Resulted from a higher ventilation rate, the stable cavity length was increased by 50% for Case 2 compared with Case 1. The bubbly flow downstream the cavity was therefore affected. The size for the vortex regions characterized by high void fraction distributions and the nature of the swirling streamlines just below the base of the ventilated cavity was also presented in the figure. In comparison with Case 1, the vortex region length was expanded by 16% for Case 2, which was mainly caused by the higher jet velocity at the cavity tail. The detail distributions of the bubbly flow parameters including the void fraction and the Sauter mean bubble diameter at different axial positions are given in Fig. 8. A good agreement was observed between the predicted parameters and the experimental results. However, disagreement was found for the void fraction distribution in the near wall region, indicating the underestimate of the free-bubble region thickness. This maybe resulted from the limitation of the radial force model. The flow turbulence and volumetric rate ratio between gas and liquid were considered as the main parameter to determine the void fraction in the fully developed pipe flow region. However, smaller bubbles were generated in Case 2, showing a reversed variation trend with the void fraction. This was resulted from the stronger turbulence in the vortex region which enhanced the breakage of bubbles, demonstrating the dominant effect of the turbulence in the vortex region on the bubble size evolution.

    5.2 Effect of inlet turbulence intensity

    The onset velocity for gas entrainment is observed to reduce with the increasing of inlet turbulence intensity, approaching around 0.5-1.0 m·s?1for rough cavity interface [27-29], therefore the gas entrainment rate and the bubbly flow downstream the cavity is influenced as well. In application of the empirical formula calculating the onset velocity for the vertical circular plunging jet [22], the variation of gas entrainment rate along with the cavity length under different inlet turbulence intensity represented by Itwas obtained as shown in Fig. 9. The liquid volumetric rate was maintained as 1.17×10?3m3·s?1. Resulted from a shorter cavity length for the onset of gas entrainment, the initial part of the variation curve was shifted left when the turbulence intensity was enhanced. However, as the increasing of the cavity length, the effect of turbulence intensity gradually diminished,leading to a coincidence of the two curves. As shown in the figure, with a stable cavity length of 400 mm, the gas entrainment rate was raised by 50% when the turbulence intensity was increased from 1% to 3%.

    Figure 7 Void fraction distributions and vortex region size under different ventilation rate

    Figure 8 Bubbly flow parameters distributions in the fully developed pipe flow region under different ventilation rate

    Figure 9 Variation of gas entrainment rate along with the cavity length under different inlet turbulence intensity

    Figure 10 presents detail distributions of the bubbly flow parameters for further investigation on the inlet turbulence intensity effect. The average void fraction at different sections was all augmented with higher turbulence intensity. The local maximum void fraction at the exit section of the vortex region was increased by 50% which was in accordance with the raising magnitude of the gas entrainment rate, indicating the aforementioned dominant effect of the volumetric ratio between liquid and gas on the void fraction. Nevertheless, the bubble size was only slightly increased at the exit of the vortex region due to high turbulence. In the pipe flow region, the high void fraction effect gradually dominated, causing the average bubble diameter increased by 6%. When the cavity was expanded beyond 600 mm, the gas entrainment rate and the resulted bubbly flow were unaffected by the inlet turbulence intensity.

    5.3 Geometric scale effect analysis

    Simulations were conducted with a scaled up model with the geometric scaling ratio rl=2. The simulated liquid flow conditions (Table 2) were chosen based upon the Froude similitude in comparison with the experimental condition of Cases 2 and 3. The ventilation rate required to maintain the same geometric scaling ratio for the cavity length was calculated according to the air entrainment model. For identical fluids in all simulated cases, the Froude similitude implies that the Weber number differs between simulations that surface tension-dominated processes may not be properly scaled. Comparative analysis of dimensionless water velocity, void fraction and bubble size were given in Fig. 11. The distributions of velocity were properly scaled with a Froude similitude for the investigated flow conditions, which was illustrated by the good agreement in terms of dimensionless distributions of water velocity. Caused by faster dispersion of bubbles, a slightly lower void fraction was observed in the larger size model. The entrained bubbles, as well as the bubbles downstream the ventilated cavity were comparatively larger in the smaller size model,implying significant scale effects in terms of turbulence and surface tension effect. In dimensional terms, the size distributions of the smallest bubbles in the pipe flow region were about the same for both models.

    6 CONCLUSIONS

    Figure 10 Bubbly flow parameters distributions in the fully developed pipe flow region

    Table 2 Computational condition for scale effect analysis

    Figure 11 Comparison of dimensionless distributions of water velocity, void fraction and bubble size

    A numerical scheme for the ventilated cavity was established by coupling the Eulerian-Eulerian two-fluid model with a semi-empirical gas entrainment model. Simulations were carried out to explore the effect of ventilation rate, inlet turbulence intensity and geometric scale effect on the bubbly flow downstream of the ventilated cavity. As the ventilation rate was doubled,the void fraction in the pipe flow region was raised up to the same degree. However, the bubble size was reduced by 20%, indicating that the bubble size evolution process is dominated by the breakage effect induced by turbulence in the vortex region. When the inlet turbulence intensity was enhanced, the gas entrainment rate and the resulted bubbles were only affected at the initial stage of gas entrainment where the intermittence factor remained below 1. In conclusion, the proposed model provides a useful methodology to predict the multiscale multiphase flow field created by ventilated cavity, which will lead a valuable insight in designing and controlling of the two phase systems with ventilated cavity existing.

    NOMEMCALATURE

    BBbirth rate due to breakage, m?3·s?1

    BCbirth rate due to coalescence, m?3·s?1

    C1adjustable model constant

    Cμk-ε turbulence model constant

    Cμbbubble induced turbulent viscosity constant

    D diameter, m

    DBdeath rate due to breakage, m?3·s?1

    DCdeath rate due to coalescence, m?3·s?1

    DsSauter mean bubble diameter, m

    Fitotal interfacial force, kg·m·s?2

    g gravitational acceleration, m·s?2

    hWsurface roughness, m

    K proportionality constant

    k turbulent kinetic energy, m2·s?2

    Lccavity length, m

    Loncavity length for the onset of gas entrainment, m

    nkaverage number density of the kth class, m?3

    P pressure, kg·m?1·s?2

    Q volumetric rate, m3·s?1

    q gas volumetric rate per unit peripheral width of the cavity base, m2·s?1

    Rccavity radius, m

    Re Reynolds number

    r distance from the pipe centerline, m

    rlgeometric scaling ratio

    rmaxmaximum stable bubble diameter, m

    U superficial velocity, m·s?1

    Ueonset velocity for gas entrainment, m·s?1

    Ujwater jet velocity at the cavity base, m·s?1

    u velocity, m·s?1

    Weccritical Weber number

    z axial distance, m

    α void fraction

    δ liquid film thickness, m

    ε turbulent kinetic energy dissipation rate, m2·s?3

    μ molecular viscosity, kg·m?1·s?1

    μt,bbubble induced turbulent viscosity, kg·m?1·s?1

    μeeffective viscosity, kg·m?1·s?1

    ν kinematic viscosity, m2·s?1

    ρ density, m3·kg?1

    σ surface tension coefficient, kg·s?2

    Subscripts

    c cavity

    g gas

    i index of gas/liquid phase

    l liquid

    t turbulence induced

    REFERENCES

    1 Bacon, R.P., Scott, D.M., Thorpe, R.B., “Large bubbles attached to spargers in downwards two-phase flow”, Int. J. Multiphase Flow, 21, 949-959 (1995).

    2 Riiser, K., Fabre, J., Suzanne, C., “Gas entrainment at the rear of a Taylor bubble”, In: Proceedings of the European Two-phase Flow Group Meeting, Stockholm (1992).

    3 Lee, Y.H., Scott, D.M., Thorpe, R.B., “The scale-up of large bubbles attached to spargers in downward two-phase flow”, Chem. Eng. Sci, 52, 3797-3809 (1997).

    4 Su, C., “Gas entrainment at the bottom of a Taylor bubble in vertical gas-liquid slug flows”, Ph.D. Thesis, University of Houston, USA (1995).

    5 Delfos, R., Wisse, C.J., Oliemans, R.V.A., “Measurement of air-entrainment from a stationary Taylor bubble in a vertical tube”, Int. J. Multiphase Flow, 27, 1769-1787 (2001).

    6 Sotiriadis, A.A., Thorpe, R.B., “Liquid re-circulation in turbulent vertical pipe flow behind a cylindrical bluffbody and a ventilated cavity attached to a sparger”, Chem. Eng. Sci., 60, 981-994 (2005).

    7 Chanson, H., Toombes, L., “Hydraulics of stepped chutes: the transition flow”, J. Hyd. Res. IAHR, 42, 43-54 (2004).

    8 Henderson, F.M., Open Channel Flow, MacMillan Co., New York (1966).

    9 Xiang, M., Cheung, S.C.P., Yeoh, G.H., Zhang, W.H., Tu, J.Y., “On the numerical study of bubbly flow created by ventilated cavity in vertical pipe”, Int. J. Multiphase Flow, 37, 756-768 (2011).

    10 Chanson, H., Air Bubble Entrainment in Free-Surface Turbulent Shear Flows. Academic Press, London (1996).

    11 Chanson, H., Brattberg, T., “Air entrainment by two-dimensional plunging jets: The impingement region and the very-near flow field”, In: Proceedings of the ASME Fluids Engineering Conference, Washington DC (1998).

    12 Ma, J., Oberai, A.A., “A quantitative sub-grid air entrainment model for bubbly flows—Plunging jets”, Comput. Fluids, 39, 77-86 (2010).

    13 Chanson, H., Manasseh, R., “Air entrainment processes in a circular plunging jet: Void-fraction and acoustic measurements”, J. Fluids Eng., 125, 910-921 (2003).

    14 Fernandes, R.C., Semiat, R., “Hydrodynamic model for gas-liquid slug flow in vertical tubes”, AlChE J., 29 (6), 981-989 (1983).

    15 Fu, X., “Interfacial area measurement and transport modeling in air-water two-phase flow”, Ph.D. Thesis, Purdue University, USA (2001).

    16 Kockx, J.P., Nieuwstadt, F.T.M., “Gas entrainment by a liquid film falling around a stationary Taylor bubble in a vertical tube”, Int. J. Multiphase Flow, 31, 1-24 (2005).

    17 Sene, K.J., “Air entrainment by plunging jets”, Chem. Eng. Sci., 43, 2615-2623 (1988).

    18 Johnson, M.W., Fashifar, A., “Statistical properties of turbulent bursts in transitional boundary layers”, Int. J. Heat Fluid Flow, 15, 283-290 (1994).

    19 Evans, G.M., Machniewski, P.M., Bin, A.K., “Bubble size distribution and void fraction in the wake region below a ventilated gas cavity in downward pipe flow”, Chem. Eng. Res. Des., 82 (9), 1095-1104 (2004).

    20 Dumitrescu, D.T., “Stromung an einer Luftblase im Senkrechten Rohr”, Z. angew., Math. Mech., 23, 139-149 (1943).

    21 Thorpe, R.B., Evans, G.M., Zhang, K., “Liquid recirculation and bubble breakup beneath ventilated gas cavities in downward pipe flow”, Chem. Eng. Sci., 56, 6399-6409 (2001).

    22 Ervine, D., Mckeogh, E., Elsawy, E., “Effect of turbulence intensity on the rate of air entrainment by plunging water jets”, ICE Proc., 69, 425-450 (1980).

    23 Shah, A., Chughtai, I.R., Inayat, M.H., “Numerical simulation of direct-contact condensation from a supersonic steam jet in subcooled water”, Chin. J. Chem. Eng., 18 (4), 577-587 (2010).

    24 Kumar, S., Ramkrishna, D., “On the solution of population balance equations by discretization—I. A fixed pivot technique”, Chem. Eng. Sci., 51, 1311-1332 (1996).

    25 Launder, B., Salding, D., Lectures in Mathematical Models for Turbulence, Academic Press, New York (1972).

    26 Sato, Y., Sadatomi, M., Sekoguchi, K., “Momentum and heat transfer in two-phase bubbly flow”, Int. J. Multiphase Flow, 7, 167-178 (1981).

    27 Cummings, P.D., “Aeration due to breaking waves”, Ph.D. Thesis, University of Queensland, Australia (1996).

    28 Mckeogh, E.J., Ervine, D.A., “Air entrainment rate and diffusion pattern of plunging liquid jets”, Chem. Eng. Sci., 36, 1161-1172 (1981).

    29 Wood, I.R., IAHR Hydraulic Structures Design Manual, Balkema Press, Rotterdam, 1-149 (1991).

    FLUID DYNAMICS AND TRANSPORT PHENOMENA

    Chinese Journal of Chemical Engineering, 22(3) 252—260 (2014)

    10.1016/S1004-9541(14)60033-1

    2012-10-25, accepted 2013-04-23.

    *Supported by the Research Project Foundation of National University of Defense Technology (JC12-01-04) and the National Science Foundation for Post-doctoral Scientists of China (2012M520268).

    **To whom correspondence should be addressed. E-mail: xiangmin333@hotmail.com

    猜你喜歡
    基元
    基于化學(xué)學(xué)科理解視角,建構(gòu)判斷“基元反應(yīng)”的思維模型
    關(guān)注基元反應(yīng)的考查
    面向游戲場景生成的細(xì)分插槽WFC算法研究
    圖像聲納換能器指向性優(yōu)化研究
    基于多重示范的智能車輛運(yùn)動(dòng)基元表征與序列生成
    人體細(xì)胞內(nèi)存在全新DNA結(jié)構(gòu)
    基元樹建筑物圖像偽造組件檢測(cè)算法
    一種彩色圖像的主基元圖模型和提取算法
    基于圓弧基元的工件實(shí)時(shí)定位與匹配方法
    基于二維壓差矢量傳感器的定向誤差分析
    久久这里有精品视频免费| 晚上一个人看的免费电影| 亚洲国产精品国产精品| 大香蕉久久网| 最新的欧美精品一区二区| 成人国语在线视频| 国产在视频线精品| 伊人久久国产一区二区| 老司机影院毛片| 观看av在线不卡| 精品一区在线观看国产| 侵犯人妻中文字幕一二三四区| 国产亚洲av片在线观看秒播厂| 国产白丝娇喘喷水9色精品| 国产精品无大码| 一级片'在线观看视频| av电影中文网址| 天天躁夜夜躁狠狠躁躁| 国产伦理片在线播放av一区| 国产精品无大码| 青春草亚洲视频在线观看| 水蜜桃什么品种好| 免费看不卡的av| 国产成人精品久久久久久| 精品久久久久久电影网| 一区二区三区乱码不卡18| 中文字幕亚洲精品专区| 国产精品.久久久| freevideosex欧美| 亚洲情色 制服丝袜| 最近2019中文字幕mv第一页| 免费少妇av软件| 香蕉丝袜av| 精品一区二区免费观看| 免费在线观看视频国产中文字幕亚洲 | 90打野战视频偷拍视频| 国产精品秋霞免费鲁丝片| 国产1区2区3区精品| 欧美成人精品欧美一级黄| 亚洲男人天堂网一区| av网站免费在线观看视频| 美女大奶头黄色视频| 性色avwww在线观看| 色视频在线一区二区三区| 一本久久精品| 999精品在线视频| 亚洲伊人色综图| 男人添女人高潮全过程视频| 午夜福利,免费看| 国产欧美亚洲国产| 精品酒店卫生间| 熟女电影av网| 欧美日韩成人在线一区二区| 免费不卡的大黄色大毛片视频在线观看| 国产一区二区 视频在线| 性色avwww在线观看| 99久久精品国产国产毛片| 久久精品久久精品一区二区三区| 精品一区在线观看国产| 97在线视频观看| 国产精品亚洲av一区麻豆 | 激情五月婷婷亚洲| 欧美成人午夜精品| av有码第一页| 色吧在线观看| 中文字幕精品免费在线观看视频| 少妇被粗大的猛进出69影院| 91国产中文字幕| 精品第一国产精品| 女的被弄到高潮叫床怎么办| 99国产综合亚洲精品| 久久久久久伊人网av| 天天躁日日躁夜夜躁夜夜| 极品少妇高潮喷水抽搐| 永久免费av网站大全| 亚洲国产av影院在线观看| 亚洲内射少妇av| 国产精品免费大片| 久久精品久久精品一区二区三区| 久久精品久久久久久久性| 少妇被粗大的猛进出69影院| 在线天堂最新版资源| 亚洲精品视频女| 国产一区二区在线观看av| 日本av手机在线免费观看| 人妻系列 视频| 久久av网站| 国产精品偷伦视频观看了| 日韩成人av中文字幕在线观看| 亚洲av综合色区一区| 99国产精品免费福利视频| 成人黄色视频免费在线看| 啦啦啦中文免费视频观看日本| 欧美 日韩 精品 国产| www.av在线官网国产| 日韩伦理黄色片| 国产精品久久久久久av不卡| 久久久国产欧美日韩av| 制服丝袜香蕉在线| 在线观看国产h片| 亚洲精品中文字幕在线视频| 欧美最新免费一区二区三区| 精品午夜福利在线看| 2021少妇久久久久久久久久久| 久久久久久久亚洲中文字幕| 国产精品国产三级专区第一集| 色视频在线一区二区三区| 中文字幕制服av| 91精品三级在线观看| 99国产综合亚洲精品| 女人久久www免费人成看片| 国产成人精品福利久久| 女人久久www免费人成看片| 国产成人午夜福利电影在线观看| 国产精品国产三级专区第一集| 观看美女的网站| 欧美日韩亚洲高清精品| 丝袜美腿诱惑在线| 国产精品三级大全| 超碰成人久久| 亚洲精品第二区| 韩国av在线不卡| 日韩 亚洲 欧美在线| 在线免费观看不下载黄p国产| 少妇被粗大的猛进出69影院| 制服诱惑二区| 欧美 亚洲 国产 日韩一| 久久精品aⅴ一区二区三区四区 | 色婷婷久久久亚洲欧美| 国产成人午夜福利电影在线观看| 亚洲av成人精品一二三区| 一级毛片我不卡| 国产精品麻豆人妻色哟哟久久| 日本91视频免费播放| 亚洲欧洲日产国产| 99久久精品国产国产毛片| 亚洲四区av| 性色av一级| 久久久a久久爽久久v久久| 精品酒店卫生间| 午夜日本视频在线| 日本vs欧美在线观看视频| 欧美日韩视频高清一区二区三区二| 色婷婷久久久亚洲欧美| 久久精品国产亚洲av天美| 另类亚洲欧美激情| 精品国产乱码久久久久久男人| 91国产中文字幕| 国产精品欧美亚洲77777| 亚洲欧美清纯卡通| videos熟女内射| 交换朋友夫妻互换小说| 卡戴珊不雅视频在线播放| av在线老鸭窝| 青草久久国产| 春色校园在线视频观看| 人妻一区二区av| 一区二区av电影网| 9热在线视频观看99| 免费在线观看完整版高清| 日韩不卡一区二区三区视频在线| 国产精品99久久99久久久不卡 | 国产免费福利视频在线观看| 精品少妇久久久久久888优播| 亚洲图色成人| 1024视频免费在线观看| 欧美日韩国产mv在线观看视频| 国产色婷婷99| 26uuu在线亚洲综合色| 日本欧美国产在线视频| 岛国毛片在线播放| 欧美国产精品一级二级三级| 久久精品国产亚洲av天美| 欧美精品一区二区大全| 精品少妇内射三级| av视频免费观看在线观看| 哪个播放器可以免费观看大片| 丰满少妇做爰视频| 国产毛片在线视频| 桃花免费在线播放| 少妇 在线观看| 国产成人精品福利久久| √禁漫天堂资源中文www| 国产欧美日韩一区二区三区在线| 亚洲精品日韩在线中文字幕| 一本大道久久a久久精品| 热99久久久久精品小说推荐| 亚洲伊人久久精品综合| 久久久久久免费高清国产稀缺| 国产精品一区二区在线观看99| 国产成人精品婷婷| 老熟女久久久| 91aial.com中文字幕在线观看| 最近手机中文字幕大全| 捣出白浆h1v1| 伦理电影免费视频| 天天躁狠狠躁夜夜躁狠狠躁| 少妇被粗大猛烈的视频| 水蜜桃什么品种好| 亚洲经典国产精华液单| 欧美精品av麻豆av| 亚洲三级黄色毛片| 精品久久久久久电影网| 亚洲精品成人av观看孕妇| 成人免费观看视频高清| 成人18禁高潮啪啪吃奶动态图| 亚洲少妇的诱惑av| 91午夜精品亚洲一区二区三区| 国精品久久久久久国模美| 久久狼人影院| 精品一区二区三区四区五区乱码 | 国产片内射在线| 老熟女久久久| 汤姆久久久久久久影院中文字幕| 好男人视频免费观看在线| 午夜免费男女啪啪视频观看| 成人国语在线视频| www.自偷自拍.com| 91国产中文字幕| 亚洲四区av| 日本黄色日本黄色录像| 搡老乐熟女国产| 日日啪夜夜爽| 亚洲美女黄色视频免费看| 午夜日韩欧美国产| 精品少妇黑人巨大在线播放| 丰满乱子伦码专区| 国产成人精品一,二区| 香蕉国产在线看| 国产精品久久久久成人av| 99re6热这里在线精品视频| 99热网站在线观看| 晚上一个人看的免费电影| 美女大奶头黄色视频| 王馨瑶露胸无遮挡在线观看| 国产探花极品一区二区| 男女高潮啪啪啪动态图| 精品国产一区二区三区四区第35| 一区二区av电影网| 欧美日韩亚洲高清精品| 中文欧美无线码| av不卡在线播放| 亚洲国产日韩一区二区| 欧美另类一区| 国产色婷婷99| 成人国语在线视频| 乱人伦中国视频| 免费高清在线观看视频在线观看| 亚洲国产色片| 大陆偷拍与自拍| 久久久久久久久免费视频了| 精品一区二区三卡| 亚洲一级一片aⅴ在线观看| 久热久热在线精品观看| 边亲边吃奶的免费视频| 2021少妇久久久久久久久久久| 日韩,欧美,国产一区二区三区| 午夜久久久在线观看| 国产深夜福利视频在线观看| 国产黄色免费在线视频| 亚洲国产最新在线播放| 精品酒店卫生间| 中国国产av一级| 精品国产一区二区三区四区第35| 各种免费的搞黄视频| 美女视频免费永久观看网站| 一区二区三区精品91| av不卡在线播放| videosex国产| 国产免费视频播放在线视频| 久久久久精品人妻al黑| 在线观看美女被高潮喷水网站| 成年动漫av网址| 日韩制服骚丝袜av| 久久久久久人妻| 国产黄频视频在线观看| 最近的中文字幕免费完整| 你懂的网址亚洲精品在线观看| 看免费av毛片| 免费在线观看完整版高清| 国产男女内射视频| 国产成人精品一,二区| 国产精品秋霞免费鲁丝片| 超碰97精品在线观看| 边亲边吃奶的免费视频| 中文乱码字字幕精品一区二区三区| 老司机影院毛片| xxxhd国产人妻xxx| 国产探花极品一区二区| 亚洲欧洲精品一区二区精品久久久 | 激情视频va一区二区三区| 成人毛片a级毛片在线播放| 嫩草影院入口| 人人妻人人澡人人爽人人夜夜| 久久久久久久久免费视频了| 国产毛片在线视频| 国产日韩一区二区三区精品不卡| 欧美变态另类bdsm刘玥| 国产男女内射视频| 在线观看免费视频网站a站| 一区二区av电影网| 丝瓜视频免费看黄片| 国产精品成人在线| 妹子高潮喷水视频| 亚洲av中文av极速乱| 99久国产av精品国产电影| 国产男女超爽视频在线观看| 午夜免费鲁丝| 国产探花极品一区二区| 日本爱情动作片www.在线观看| 亚洲视频免费观看视频| 丰满乱子伦码专区| 久久毛片免费看一区二区三区| 又粗又硬又长又爽又黄的视频| av线在线观看网站| 极品人妻少妇av视频| 中文字幕亚洲精品专区| 精品国产乱码久久久久久小说| 久久精品久久久久久久性| 1024视频免费在线观看| 精品少妇一区二区三区视频日本电影 | 女人久久www免费人成看片| 国产综合精华液| 另类亚洲欧美激情| 国产精品欧美亚洲77777| 精品一区二区三卡| 日本猛色少妇xxxxx猛交久久| 久久久a久久爽久久v久久| 国产麻豆69| 天天躁夜夜躁狠狠躁躁| 搡女人真爽免费视频火全软件| 女的被弄到高潮叫床怎么办| 免费在线观看黄色视频的| 人妻系列 视频| 1024视频免费在线观看| 麻豆乱淫一区二区| 2022亚洲国产成人精品| 丰满饥渴人妻一区二区三| a级毛片在线看网站| 免费日韩欧美在线观看| 亚洲一码二码三码区别大吗| 国产精品熟女久久久久浪| 精品国产一区二区三区久久久樱花| 热re99久久国产66热| 美女高潮到喷水免费观看| 免费观看性生交大片5| 99国产精品免费福利视频| 亚洲情色 制服丝袜| 国产精品不卡视频一区二区| 久久久久久久久免费视频了| 欧美人与善性xxx| 久久ye,这里只有精品| 男人舔女人的私密视频| 母亲3免费完整高清在线观看 | 大香蕉久久网| 日韩制服骚丝袜av| 观看av在线不卡| 免费在线观看黄色视频的| 一区二区三区乱码不卡18| 国产成人午夜福利电影在线观看| 一级毛片 在线播放| 亚洲av综合色区一区| 高清黄色对白视频在线免费看| 亚洲,一卡二卡三卡| 久久久精品免费免费高清| 久久久精品国产亚洲av高清涩受| 国产日韩欧美视频二区| 日本免费在线观看一区| 两个人看的免费小视频| 久久久久久久亚洲中文字幕| 国产精品偷伦视频观看了| 欧美 亚洲 国产 日韩一| www.av在线官网国产| 国产精品欧美亚洲77777| 国产成人aa在线观看| 精品国产一区二区三区久久久樱花| 中文字幕最新亚洲高清| 日日啪夜夜爽| 国产福利在线免费观看视频| 中文字幕av电影在线播放| 一区二区三区四区激情视频| 欧美日韩精品成人综合77777| av一本久久久久| 国产精品99久久99久久久不卡 | 美女国产视频在线观看| 久久久久国产精品人妻一区二区| 亚洲成人av在线免费| 9色porny在线观看| 国产精品麻豆人妻色哟哟久久| 青草久久国产| 男人爽女人下面视频在线观看| 免费大片黄手机在线观看| 多毛熟女@视频| 精品亚洲乱码少妇综合久久| 午夜福利网站1000一区二区三区| 中文字幕av电影在线播放| 久久久欧美国产精品| 一级毛片 在线播放| 国产精品久久久久久精品古装| 午夜福利在线观看免费完整高清在| 久久久久久久精品精品| 免费少妇av软件| 香蕉国产在线看| 欧美日韩精品网址| 国产av码专区亚洲av| 国产1区2区3区精品| 成人亚洲精品一区在线观看| 国产精品久久久久成人av| 国产成人精品无人区| 成年人午夜在线观看视频| 老鸭窝网址在线观看| 在线免费观看不下载黄p国产| 免费人妻精品一区二区三区视频| 久久精品久久久久久久性| 免费观看a级毛片全部| 成年av动漫网址| 一级毛片 在线播放| 日本欧美视频一区| 久久ye,这里只有精品| 伦精品一区二区三区| 啦啦啦在线免费观看视频4| 亚洲成色77777| 一本大道久久a久久精品| 亚洲一码二码三码区别大吗| 一边摸一边做爽爽视频免费| 久久精品国产自在天天线| 男人舔女人的私密视频| 爱豆传媒免费全集在线观看| 欧美变态另类bdsm刘玥| 久久精品久久精品一区二区三区| 欧美日韩综合久久久久久| 国产精品熟女久久久久浪| 国产av码专区亚洲av| 丝袜美足系列| 老汉色av国产亚洲站长工具| 亚洲国产看品久久| 丝袜人妻中文字幕| 天天躁日日躁夜夜躁夜夜| 热99久久久久精品小说推荐| 一级片'在线观看视频| av线在线观看网站| 这个男人来自地球电影免费观看 | 一级片'在线观看视频| 少妇被粗大的猛进出69影院| 91在线精品国自产拍蜜月| 久久久国产精品麻豆| 一级毛片我不卡| 永久免费av网站大全| 人妻系列 视频| av福利片在线| 一级爰片在线观看| 伦理电影免费视频| 伦理电影大哥的女人| 一级毛片我不卡| 中文天堂在线官网| 美女视频免费永久观看网站| 国产97色在线日韩免费| 午夜福利乱码中文字幕| 啦啦啦中文免费视频观看日本| 黄色配什么色好看| 精品久久久精品久久久| 亚洲激情五月婷婷啪啪| 久热久热在线精品观看| 国产精品亚洲av一区麻豆 | 1024视频免费在线观看| 久久久久久人人人人人| 久久久久国产一级毛片高清牌| 午夜福利乱码中文字幕| 免费观看性生交大片5| 人体艺术视频欧美日本| 免费黄色在线免费观看| 日韩不卡一区二区三区视频在线| 男人添女人高潮全过程视频| 国产成人av激情在线播放| 久久久久久久久久久免费av| 亚洲国产色片| av免费观看日本| 亚洲少妇的诱惑av| av一本久久久久| 国产在线免费精品| 超碰成人久久| 两性夫妻黄色片| 女人高潮潮喷娇喘18禁视频| 日本猛色少妇xxxxx猛交久久| 中文字幕人妻丝袜一区二区 | videos熟女内射| 欧美少妇被猛烈插入视频| 亚洲久久久国产精品| 精品久久久精品久久久| 18在线观看网站| 亚洲美女搞黄在线观看| 国产亚洲欧美精品永久| 欧美黄色片欧美黄色片| 天天躁日日躁夜夜躁夜夜| 精品人妻偷拍中文字幕| 亚洲国产成人一精品久久久| 午夜激情久久久久久久| 老汉色av国产亚洲站长工具| 99久国产av精品国产电影| 国产精品偷伦视频观看了| 高清在线视频一区二区三区| 精品一区在线观看国产| 人人妻人人澡人人爽人人夜夜| 国产精品久久久久久精品电影小说| 制服诱惑二区| 国产日韩欧美在线精品| 午夜免费鲁丝| 大话2 男鬼变身卡| 亚洲国产欧美日韩在线播放| 亚洲成国产人片在线观看| 黄片小视频在线播放| 嫩草影院入口| 亚洲精品视频女| 九九爱精品视频在线观看| 午夜激情久久久久久久| 黄色配什么色好看| 人妻少妇偷人精品九色| 久久久久精品久久久久真实原创| 日韩一区二区视频免费看| 精品少妇久久久久久888优播| 免费女性裸体啪啪无遮挡网站| 天堂俺去俺来也www色官网| 精品一区二区三区四区五区乱码 | 欧美成人午夜精品| 久久久久久伊人网av| 国产不卡av网站在线观看| 91精品国产国语对白视频| 永久免费av网站大全| 欧美日韩亚洲高清精品| 免费久久久久久久精品成人欧美视频| 岛国毛片在线播放| 午夜激情av网站| 人人妻人人澡人人看| 男女无遮挡免费网站观看| 国精品久久久久久国模美| 人成视频在线观看免费观看| 纯流量卡能插随身wifi吗| 国产男人的电影天堂91| 在线天堂中文资源库| 不卡av一区二区三区| 一区在线观看完整版| 青春草视频在线免费观看| 欧美精品av麻豆av| 精品国产露脸久久av麻豆| 天堂俺去俺来也www色官网| 亚洲一级一片aⅴ在线观看| 伦精品一区二区三区| 欧美国产精品va在线观看不卡| 国产熟女午夜一区二区三区| 国产xxxxx性猛交| 亚洲精品久久午夜乱码| 天天操日日干夜夜撸| 在线亚洲精品国产二区图片欧美| √禁漫天堂资源中文www| 色94色欧美一区二区| 国产精品嫩草影院av在线观看| 午夜av观看不卡| 国产xxxxx性猛交| 高清不卡的av网站| 中文天堂在线官网| 亚洲av在线观看美女高潮| 在线观看国产h片| 中文字幕人妻丝袜制服| 超碰97精品在线观看| 国产成人精品婷婷| videosex国产| 亚洲三区欧美一区| xxx大片免费视频| 黄色一级大片看看| 水蜜桃什么品种好| 99热网站在线观看| 亚洲在久久综合| 曰老女人黄片| 人人妻人人澡人人看| 蜜桃在线观看..| 国语对白做爰xxxⅹ性视频网站| 欧美精品一区二区免费开放| av在线老鸭窝| 啦啦啦在线免费观看视频4| 男女啪啪激烈高潮av片| 高清在线视频一区二区三区| 我的亚洲天堂| 国产成人精品无人区| 丝瓜视频免费看黄片| 精品第一国产精品| 国产成人精品无人区| av线在线观看网站| 久久久久久久久免费视频了| 成人毛片60女人毛片免费| 丝瓜视频免费看黄片| 黑人欧美特级aaaaaa片| 国产女主播在线喷水免费视频网站| 国产午夜精品一二区理论片| 美女福利国产在线| 在线观看免费日韩欧美大片| 国产精品久久久av美女十八| 国产1区2区3区精品| 亚洲人成电影观看| 美女午夜性视频免费| 91精品国产国语对白视频| 欧美日韩成人在线一区二区| 精品国产一区二区三区四区第35| 丝袜美足系列| 两性夫妻黄色片| 日韩一区二区视频免费看| av电影中文网址| 看非洲黑人一级黄片| 宅男免费午夜| 不卡av一区二区三区| 国产精品秋霞免费鲁丝片| 欧美+日韩+精品| 男人爽女人下面视频在线观看| 日本vs欧美在线观看视频| 久久久久精品人妻al黑| 亚洲精品第二区|