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

    Numerical analysis of bubble dynamics in the diffuser of a jet pump under variable ambient pressure*

    2017-06-07 08:22:46XinpingLong龍新平QingqingWang王晴晴LongzhouXiao肖龍洲JunqiangZhang章君強MaosenXu徐茂森WeifengWu吳偉烽BinJi季斌
    水動力學研究與進展 B輯 2017年3期
    關(guān)鍵詞:吳偉新平

    Xin-ping Long (龍新平), Qing-qing Wang (王晴晴), Long-zhou Xiao (肖龍洲), Jun-qiang Zhang (章君強), Mao-sen Xu (徐茂森), Wei-feng Wu (吳偉烽), Bin Ji (季斌)

    1. State Key Laboratory of Water Resources and Hydropower Engineering Science, School of Power and Mechanical Engineering, Wuhan University, Wuhan 430072, China, Email: xplong@whu.edu.cn

    2. Hubei Key Laboratory of Waterjet Theory and New Technology, Wuhan University, Wuhan 430072, China

    3. School of Energy and Power, Xi’an Jiaotong University, Xi’an 710049, China

    Numerical analysis of bubble dynamics in the diffuser of a jet pump under variable ambient pressure*

    Xin-ping Long (龍新平)1,2, Qing-qing Wang (王晴晴)1,2, Long-zhou Xiao (肖龍洲)1,2, Jun-qiang Zhang (章君強)1,2, Mao-sen Xu (徐茂森)1,2, Wei-feng Wu (吳偉烽)3, Bin Ji (季斌)1,2

    1. State Key Laboratory of Water Resources and Hydropower Engineering Science, School of Power and Mechanical Engineering, Wuhan University, Wuhan 430072, China, Email: xplong@whu.edu.cn

    2. Hubei Key Laboratory of Waterjet Theory and New Technology, Wuhan University, Wuhan 430072, China

    3. School of Energy and Power, Xi’an Jiaotong University, Xi’an 710049, China

    2017,29(3):510-519

    Recent studies have shown that the collapse of cavitation bubbles in a jet pump can generate an extremely high pressure with many potential applications. The dynamics of the bubble is governed by the Rayleigh-Plesset equation. With the bubble dynamics equation and the heat and mass transfer model solved with the Runge-Kutta fourth order adaptive step size method, the oscillations of the bubble in the diffuser of the jet pump are assessed under varied conditions. To obtain the pressure variation along the diffuser, the Bernoulli equation and the pressure measured in experiment are coupled. The results of simulation show that a transient motion of the bubbles can be obtained in the diffuser quantitatively, to obtain the pressure and temperature shock in the bubble. Moreover, increasing the outlet pressure coefficient would result in a more intense bubble collapsing process, which can be used in the subsequent studies of the cavitation applications. The predictions are compared with experiments with good agreement.

    Cavitation, jet pump, bubble collapse

    Introduction

    The cavitation occurs in industrial structures, such as the rotating machinery, the injectors, the jet pumps, and other hydraulic devices. Most of the time,it is accompanied with serious effects like the erosion, the noise and the decrease of the fluid-machinery efficiency. However, owing to the generation of hot spots, highly reactive free radicals and turbulence associated with liquid circulation currents, in industrial applications, biological engineering[1], water disinfection[2], and even food industry, like milk sterilization, the cavitation has attracted much of attention during past several decades, and our understanding of the cavita-tion has been greatly improved. Recent reviews were made by Gogate and Kabadi[3], Luo et al.[4]and Peng and Shimizu[5].

    To understand the fundamental mechanism of the cavitation phenomena[6], experiments were carried out for the growth and collapse of individual bubbles[7-12].

    For numerical modeling of gas bubbles oscillating in liquids, the first analysis of cavitation and bubble dynamics was made by Rayleigh. Plesset further considered the influence of the physical characteristics of fluid viscosity and surface tension and derived the Rayleigh-Plesset equation. For many real physical problems, the shape of a bubble is unlikely to remain to be spherical, especially during the collapse phase of its motion, or when the bubble is near the boundary of the fluid. In this case, the velocity potential of the bubble motion is usually determined by the boundary element method[13].

    In simulating the bubble collapse, the processing under the ambient pressure can be divided into two approaches, i.e., the constant and variable forms.

    The constant pressure model is considered as able to simplify the calculation in investigating the motionof a bubble. Rayleigh, as mentioned above, took the ambient pressure as a constant. Plesset and Chapman simulated two specific cases of initially spherical bubbles collapsing near a plane solid wall, a bubble initially in contact with the wall and a bubble initially a distance of half its radius away from the wall at the closest point, with the assumption of the constant ambient pressure. Other researchers such as Popinet and Zaleski[14], Dular and Coutier-Delgosha[15]and Li et al.[16]also assume that the ambient pressure is constant in the simulation domain.

    On the other hand, for bubbles under a variable ambient pressure condition, however, there is no direct, exact or accepted expression for the ambient pressure expression, p¥(t ), as a controversial issue. However, it is a crucial parameter for the mechanism of the bubble motion. Several simplified functions were suggested for the ambient pressure for convenience. A harmonic functional form is expressed as p¥= p0(1±es inw t ), where e is the amplitude of the oscillation,w is the frequency and t is the time, as suggested by Zhang and Li[17]. Furthermore, the expression for the pressure induced by the external sound field on the outside of the bubble wall,is reduced to. For small values of kR0, sin(k R0)? kR0and the pressure outside the bubble could be treated as uniform[18], whereAP is the amplitude of the acoustic wave, k is the wave number,w is the angular frequency of the acoustic wave. The equation, p¥(x,t) =p0+pasin (kx -wt + j)is employed by Wang and Manmi[19] with consideration of both the phase position and the wave velocity, where0p is the hydrostatic pressure, t is the time, and k,ap,w and j are the wave number, the pressure amplitude, the angular frequency, and the phase shift of the ultrasound, respectively. Since these pressure expressions are based neither on the measurement nor on the calculation for a real flow, the results obtained can not adequately explain the cavitation related phenomena accurately such as the thermal and pressure shock in the bubble. Qin[20]developed a numerical methodology for the calculation of the ambient pressure, i.e., the pressure distribution around the bubbles in a convergent-divergent nozzle by a CFD code, based on the RNG -ke model of the turbulent flow. In this paper, in the Ray leigh-Plesset equation, a position dependence is added, D x =v ×Dt , from the original time dependent form. Based on the above treatment of the ambient pressure, Qin et al.[21]employed a heat transfer model, including the effects of conduction plus radiation, to describe the thermodynamics of the non-condensable gas inside the bubble and the results of this heat transfer model match the previously published experimental data well. Kumar and Moholkar[22]used the approach of discrete calculations while modeling the pressure variation in the nozzle flow based on the Bernoulli equation between consecutive points. As shown in the above description, the calculation of the ambient pressure was improved by Qin and Kumar significantly. Nevertheless, the calculated ambient pressure can still not well reflect the real flow. For many cavitation applications, an accurate control of cavitation is highly desirable. For example, the milk sterilization based on the cavitation needs a high sterilization efficiency to reduce the protein inactivation. As shown in the above reviews, the ambient pressure in the simulation is still a crucial parameter, and there is no efficient method describing the pressure diffusion in a real flow.

    In the present paper, the gas bubble dynamics in the liquid under the hydraulic excitation in a jet pump is numerically studied to calculate the nonlinear oscillations of the bubbles. For a thermal analysis, the heat transfer model of Qin et al.[21,23]is adopted with consideration of the energy carried by the vapor. In addition, the effect of mass diffusion is considered based on the study of Kumar and Moholkar[22]. According to our experimental results, the cavitation is induced in the outlet of the nozzle, then develops in the throat shear layer, and collapses in the diffuser finally if the cavitation number is small enough. In other word, the upstream flow of the jet pump would influence the flow in the diffuser. However, in a common diffuser, the cavitation state is absolutely different, for example, the bubble might collapse at the boundary. To sum up, the paper, considers a single bubble in the diffuser of the jet pump. The description of the ambient pressure variation is based on the pressure measured in experiment and the Bernoulli equation, which can reflect the ambient pressure variation along the diffuser in a real flow accurately. On this premise, based on the Rayleigh-Plesset equation and combined with the mass-diffusion and thermal effects, we can evaluate the pressure and temperature variations in the bubble quantitatively, which can serve as a reference for the subsequent cavitation application studies in the jet pump.

    1. Mathematical formulation

    In this paper, based on the concept of the confined co-axial jet cavitation, we integrate the structure of the orifice and the venturi tube with the jet pump theory. Strong vortex cavitation will be induced by using a high-speed jet to entrain the low-speed secondary flow in a confined space. This sort of confined co-axial jet cavitation is different from the traditional form of cavitation, which can induce cavitation under a relatively high environmental pressure due to the vortex in the shear layer. Moreover, the containedbubbles or cavity clouds can be conveyed downstream immediately, collapse under a high pressure, and consequently induce strong pressure shocks and thermal shocks.

    Similar to the traditional central jet pump, the cavitation reactor also comprises of a nozzle, a suction chamber, a throat and a diffuser, as shown in Fig.1. The main structural parameters in Fig.1 are set as follows: Lt=96mm (the length of the throat), Dt= 16mm (the diameter of the throat), Dd=50mm (the diameter of the diffuser outlet), a=45o(the inclined angle of the suction chamber), b= 12o(the inclined angle of the diffuser) and the diameter of the nozzle outlet is 12 mm. With the above parameters, we can obtain the area ratio ()m, which is the ratio between the throat area and the nozzle outlet area and takes a value of 1.78. Besides the stainless steel suction duct and nozzle, the suction chamber, the throat and the diffuser are made of transparent Perspex material for observing the cavity cloud during the mixing process.

    Fig.1 Schematic diagram of jet pump

    Moreover, the red arrows in Fig.2 represent the mounting positions of the pressure transducers.

    Fig.2 The positions of pressure transducers

    In the experiment, the working pressure is controlled at the value of 398 kPa and various working conditions can be established by adjusting the outlet pressure. Based on the published analysis[7], the outlet pressure coefficient (Cpc) is employed to describe the dimensionless outlet pressure, which is defined as

    Here cp and sp are the outlet pressure and the static pressure of the secondary flow, respectively,lr is the density of the liquid andnV is the average velocity of the nozzle outlet.

    1.1 Rayleigh-Plesset equation for a spherical bubble Consider a spherical bubble of radius R, as a function of the time t, in the domain of the jet pump cavitation generator. Neglecting the gravity and with the assumption of the liquid incompressibility, the generalized Rayleigh-Plesset equation for the motion of the bubble wall is[21]

    where S and m represent the surface tension and the viscosity of the liquid,Bp is the pressure inside the bubble, p¥is the ambient pressure. According to Qin et al.[21,23], the Rayleigh-Plesset equation can be modified by adding the position dependence from its original time dependent form of Eq.(2). This can be achieved by using the relationship between the velocity (v) and the displacement (x), Dx =v ×Dt , and then the position, xi, can be expressed as a function of the time t, that is,

    Similarly, the ambient pressure, p¥, can be expressed as

    In order to estimateBp, it is assumed[24]that a bubble contains some quantity of contaminant gas whose quality and partial pressure are0M and0gp , respectively at a certain bubble reference size0R, and some water vapor. The water vapor pressure,vp, is constant at a constant temperature, and the gas is assumed to be incondensable. With the ideal gas state equation, the following expression for the pressure in a bubble can be obtained

    Substituting this into the Rayleigh-Plesset equation enables the numerical solution of the cavitation bubble radius with consideration of the thermodynamics.

    1.2 Energy equation

    According to the first law of thermodynamics, at any time interval dt, the net energy, UD , can be expressed as

    where dw is the work done on the bubble, dQcond, dQconvand dQradrepresent the heat losses by conduction, convection and radiation, respectively, hvd Mvis the energy in the bubble carried by the vapor, hvdM = hvrvdV =4 p hvrvR2dR , hvis the vapor enthalpy,vr is the vapor density.

    It is assumed that the bubble is spherical. At any instant of time,it, the pressure and the temperature inside the bubble areBp andBT, respectively. After a finite time step tD , the volume change of the bubble is VD . The work done on the bubble is

    Assuming that the velocity of the bubble relative to the surrounding liquid is small compared to the velocity of the bubble wall during the collapse, the convection effect can be ignored. According to Qin et al.[21]and Kumar and Moholkar[22], the heat losses by conduction and radiation are expressed as

    where k is the thermal conductivity of the water surrounding the bubble,thd is the distance between the two surfaces, e is the emissivity of the hot surface and kSBis the Stefan-Boltzmann constant, kSB= 5.670′10-8Wm-2K-4. When the surface of a bubble is water, the emissivity e is taken to be 0.95. A is the area of the bubble.

    The net energy will raise the temperature of the gas in the bubble by TD given by[21,22]

    where n is the number of moles andmixC is the specific heat of the mixed gas at constant volume. For a mixed gas,mixC can be expressed as

    where,viC is the specific heat of a certain gas at a constant volume andin is the corresponding number of moles.

    Substituting Eqs.(6), (7), (8) and (10) into Eq.(11), we have

    1.3 Mass-diffusion equation

    The mass-diffusion process across the bubbleliquid interface plays a significant role in the behavior of the gas bubbles because it may ultimately determine the presence or the absence of the bubbles in a liquid. The governing equation is

    where C is the gas concentration,diffD is the coefficient of diffusion of the gas in the liquid which is determined by Eq.(15). The problem is substantially simplified as

    where C¥is the dissolved mass concentration at a place a large distance away from the bubble,sC is the equilibrium (or saturation) concentration of the gas in the liquid andgr is the density of the gas within the bubble. The key to these processes is provided by Henry’s law, a connection between the partial pressure of the gas acting on a liquid surface,gp, andsC

    To sum up, the mass-diffusion equation in the bubble-liquid interface is

    1.4 Numerical solution method

    Equations (2), (13) and (17) are the complete formulation for the radial motion of the bubble with associated heat and mass transfer effects, which can besolved simultaneously using the Runge-Kutta (RK) fourth-order method, which is a kind of special high precision numerical algorithm for solving differential equation, and is widely used in engineering practice with the truncation error of O(hR4). Furthermore, although the RK method involves a larger calculation amount as compared with the methods of Euler and Heun, its precision is much higher than the other two methods in the same step length[25]. The algorithm is well described in Ref.[26].

    As a higher order ordinary differential equation cannot be solved by RK directly, Eq.(2) needs to be cast into two first-order equations first, i.e.

    Meanwhile, the energy equation and the massdiffusion equation are rearranged as

    The calculation equation set, which contains five variables, i.e., (t, s, R, T, M), consists of four equations, Eqs.(18)-(21). In the process of calculation, we denote the time step (hR) as hR=Dt . The selection of the time step D t influences the result of the calculation greatly. Especially when the slope of the bubble radius R becomes very steep in some places, too large time step D t will increase the calculation error or even lead to numerical overflow, however, too small tD can increase the computing burden. Therefore, a variable time step length is developed to optimize the iteration process.

    Based on the above considerations, a new adaptive algorithm is developed, in which the time step varies with the ratio of the bubble sizes to reduce the computing burden and at the same time ensure a high precision[27]. The details of this method are described inRef.[23].

    2. Modelling of pressure variation in diffuser

    To describe the ambient pressure along the diffuser accurately, Bernoulli equation is applied between the sole pressure tap (x/Dt=9.00) and any other point to obtain the ambient pressure variation. The method was described by Kumar and Moholkar[22], who divided the nozzle into 100 intervals with the assumption that the limiting case is at the throat, while the ambient pressure variation is discrete. In the present study, the continuous function of the ambient pressure can be obtained due to the fact that the flow can be obtained from the flowmeter. The velocity at the cross-section in the diffuser is calculated as

    where c is dimensionless, c= x/Dt. Particularly, we denote the position of the pressure tap (x/ Dt= 9.00) as c0, then the section radius at x/Dt=9.00 can be expressed as r(c0). r(c) is expressed as

    According to the Bernoulli equation with ignoring the loss of the liquid, the total pressure at an arbitrary point in the diffuser is calculated as

    Combining Eqs.(22), (23) and (24), ()pc is expressed as

    It is assumed that a bubble travels with the flow and shares the same velocity while it travels through the diffuser and the flow field is stable. The trace line of the bubble is coincident with that of the liquid, which is dc / dt= V(c). The equation of the trace line is expressed as followswhere C0is the integration constant, which depends on the location of the bubble at t =0. The collapse locations of the cavitation cloud are different under different outlet pressures, therefore the location of the bubble at t =0 is related with a specific condition.

    Taking Cpc=0.234 as an example, the relevant parameters are: p(c0) =104.8 kP a, Qc=3.938L/s and V(c0) =9.02 6m/s . Substituting the above values into Eq.(25), the relation between the pressure distribution in the diffuser and the location is obtained. Moreover, with c= 7.120 (x=113.93mm)after several iterations the pressure is shown to be approximately equal to the vapor pressure of the liquid. Although the observed collapse of the cavitation cloud is located at c=7.690(x=123.00mm)approximately, there is a distinction between the location of the cavitation cloud collapse and the beginning of the process, which is closer to the upstream. In this paper, the location of the cavitation cloud collapse is assumed to be at c=7.120(x=113.93mm). Therefore, C0can be calculated by

    According to the above method, we can obtain the position of the cavitation generation and0C under other conditions.

    Fig.3 Comparison of radius pulsation predicted by the numerical simulation with the experimental data

    3. Results and discussions

    3.1 Experimental validations

    In order to ensure the reliability of the simulation, the simulating result is compared with the experimental data. Figure 3 shows a comparison of the bubble radius predicted by means of the above methods with measurements by Popinet and Zaleski[28]. In the process of calculation, it is assumed that the bubble is under 100.0 kPa ambient pressure at c= 0.131 (x= 2.10 mm and the bubble contraction to the minimum radius for the first time is at the time =0t . For the first collapse (0 ms-0.4 ms), the numerical prediction fits the experimental data well. For the second collapse (0.4 ms-0.8 ms), the prediction gives a higher value. According to Zhang and Li[29], the acoustic damping mechanism is the fundamental issue for the gas bubbles oscillating in liquids. In any case, the effects of the fluid compressibility, in particular, the acoustic damping, should be taken into account. This implies that the real process of the bubble collapse is more rapid. The pulsation of the bubble cannot be depicted adequately, even with the consideration of the energy equation and the mass-transfer equation. Nevertheless, the numerical bubble collapse at the initial-time in this paper is quite accurate.

    3.2 Influences of variable pressures

    Using the previous analysis combined with our previous experimental data, the ambient pressure (())pc distributions along the diffuser can be obtained, as shown in Fig.4. In this figure, the ambient pressure in the diffuser increases along the x direction, while the slope of the ambient pressure decreases gradually. In addition, the higher the outlet pressure coefficient (Cpc), the higher the ambient pressure in the diffuser and the more violent the gradient of the ambient pressure. Moreover, the initial increase stage of the pressure shown in the figure is when the bubble is under a high ambient pressure. Obviously, the position of the bubble collapse and the initial increase stage of the pressure are related withpcC .

    Fig.4 The ambient pressure distributions under different conditions

    Under the above four conditions of ambient pressure distributions and the position of the bubble generation, the process of the bubble collapse is simulated based on the previous models. During the simulation, the parameters such as the radius variation, the pre-ssure and the temperature inside the bubble are also obtained with the assumptions that the initial radius of the bubble is 2 mm at 300 K. The variations of these parameters, may give some insight into the complicated process of the bubble collapse. The entire simulation is made with VB.net.

    3.2.1 The variation of bubble radius

    The simulated variations of the radii in the process of the bubble collapse under different outlet pressures are shown in Fig.5, where the bubble finishes the cycle of the collapse and starts the rebounding due to the non-equilibrium between the pressure in the bubble and the ambient pressure. However, the loss of energy and the mass diffusion of the soluble gas in the bubble result in the decaying of the bubble oscillation. The fluctuation of the bubble decays gradually until the bubble reaches the equilibrium position. By that time, the bubble is mainly filled with the soluble gas. In the experiment, it is observed that the bubble could not vanish absolutely even at the boundary of the bubble collapse but is dispersed into small bubbles, as consistent with the numerical simulation.

    Fig.5 The variation of bubble radius under different conditions (p1=398kPa , m =1.78)

    Obviously, the fluctuations of the bubble radii are similar for different outlet pressure coefficients (Cpc). However, the amplitude of the bubble pulsation and the damping rate see a great variation due to the different ambient pressures under various conditions. The effect ofpcC on the process of the first compression of the bubble in the diffuser can be seen by comparing the results of simulations for parameters in Table 1. Here D t1and D x1/Dtare the time of the first compression of the bubble and the corresponding variation of the position, RB1is the minimum size of the bubble, RB1/RB0is the ratio between the minimum and the initial sizes of the bubble. Dt1and Dx1/Dtare decreased to 0.462 ms and 0.476 from 0.813 ms and 0.557, respectively with the increase of pcC from 0.092 to 0.228. Apparently, a smallerpcC leads to a larger velocity of the bubble travelling to downstream. On the one hand, increasingpcC leads to a more intense increase of the pressure on the collapse boundary and a larger ambient pressure. On the other hand, the velocity of the bubble to downstream is larger due to a smaller cross section area on the collapse boundary and the ratio betweenBR and t is larger in view of Eq.(20). Therefore, increasing the outlet pressure coefficient, the vibration of the bubble in both space and time domains will be more intense and the pulsation amplitude of the bubble radius will increase significantly.

    Table 1 Comparison of parameters at the first compression of the bubble

    In addition, increasing Cpcleads to a smaller damping velocity of the bubble radius. For the pulse amplitude of the bubble radius to decrease to 0.25 mm at Cpc=0.092, it will take the bubble 5 pulse processes as shown in Fig.5(a). However, it takes the bubble 14 pulse processes at Cpc=0.288to achieve the same amplitude as shown in Fig.5(d). It demonstrates that a higherpcC will lead to a smaller oscillation damping. Since improvingpcC will increase the ambient pressure for the bubble and increase the pressure difference inside and outside of the bubble, which accounts for more energy for the bubble fluctuation. It demonstrates that increasingpcC leads to a smaller mass diffusion in the bubble motion. According to themass diffusion, the difference (C¥- CS) diminishes and one may have a stable equilibrium. Absolutely, the larger Cpcresults in a smaller decrease of (C¥- CS), therefore, it will take more oscillations to reach the equilibrium state.

    During the process of the bubble collapse in the simulation, one sees a violent oscillation and it decays gradually, while the real pressure of the cavitation flow is not homogeneous, as a result, the randomness due to the interactions between the bubbles will be enhanced. In addition, the bubble may generate centripetal explosions breaking small bubbles at the first compression of the bubble. It is of engineering significance investigating the pressure and temperature shock in the process of a bubble collapse. How to obtain an appropriate pressure and temperature shock is also an important part of the study for the cavitation applications. Therefore, the following two sections will address the pressure and temperature shock resulted from the process of the bubble collapse.

    3.2.2 Pressure shock

    In the above section, the bubble is simulated by using the conduction plus radiation model and the mass-diffusion equation. During the simulation, other parameters such as the pressure and the temperature inside the bubble in the bubble dynamics equation are also obtained.

    In this section, we will discuss the oscillation of the pressure in the bubble, as shown in Figure 6, under various outlet pressures corresponding to the radius variations shown in Fig.5. Obviously, when the bubble is under a high environmental pressure, the bubble radius fluctuates drastically accompanied with a strong pressure shock. When the radius of the bubble is reduced to its minimum size, the pressure in the bubble is increased to its maximum value. The pressure shock reaches up to 0.259 MPa, which is much higher than the ambient pressure, at the first compression of the bubble when Cpc=0.092, however, the maximal amplitude decreases to 0.221 MPa at the second compression as shown in Fig.6(a). It demonstrates that the pressure pulsation will decay constantly with the bubble traveling to the downstream, and eventually, the pressure in the bubble will slightly fluctuate around the environmental pressure. In Fig.6(a), when the bubble transfers to the place x/ Dt=10.5, the pressure in the bubble is approximately the environmental pressure, which is 58.1kPa, and the pulsation is weak.

    As the process of the pressure shock is synchronous with the pulsation of the bubble radius, namely, with the compression of the bubble, the pressure in the bubble increases violently and when the bubble expands, the pressure decreases dramatically. Therefore, whenpcC increases, the time interval between the pressure shocks is shrunk and the velocity of the amplitude decaying is smaller. Figure 6 demonstrates that the oscillating amplitude is of the order of 0.1 MPa when the bubble is under three pressure shocks at Cpc=0.092. And it would take the bubble seven pressure shocks to achieve the above results at Cpc= 0.157. Similar to the variation of the bubble radius, a larger Cpcresults in a smaller decrease of (C¥- CS), then it takes more pressure shocks to reach the equilibrium state.

    Fig.6 The variation of pressure under various conditions (1=p 398kPa, =1.78m )

    Moreover, the amplitude is increased apparently with the increase of Cpc. As shown in Fig.6, when Cpc=0.092, 0.157, 0.233 and 0.288, the maximum values of the pressures are 0.157 MPa, 0.932 MPa, 2.929 MPa and 5.312 MPa, respectively. However, unlike Figs.6(a) and 6(b), the maximums of the pressure amplitudes are not observed at the first pressure shock but at the second one in Figs.6(c) and 6(d). This suggests that not only the increase of the initial-time pressure but also the ambient pressure impacts the amplitude of the pressure shock. For Figs.6(c) and 6(d), the bubble is under a larger ambient pressure when travelling to the position of the second pressure shock, which results in a higher pressure amplitude (the maxi-mum one). Nevertheless, the amplitudes at the first pressure shock in Figs.6(c) and 6(d) are still larger than the previous two. In other words, increasingpcC leads to a more intense pressure shock.

    3.2.3 Temperature shock

    Figure 7 shows the variations of the temperature in the bubble under various conditions, which is similar to the variation of the pressure shown in Fig.6. The temperature shock occurs when the bubble collapses to its minimum size. As the bubble travels to downstream and the pulsation of the bubble weakens, the amplitude of the temperature shock decays. Eventually, the temperature fluctuates around the environment temperature.

    Fig.7 The variation of temperature under various conditions (p1=398kPa , m =1.78)

    We can deduce from Fig.7 that a larger Cpcleads to a higher temperature shock. As shown in Figure 7(a), the maximum temperature is 494.9 K at Cpc=0.092, however, the maximum temperatures are 768.1 K and 1117.1 K at Cpc=0.157 and 0.233, respectively. Particularly, when Cpcincreases to 0.288, the maximum temperature is 1368.1 K during the process of the bubble contraction. Similar to the pressure shock, when Cpcincreases, the time interval between the temperature shocks is shrunk and the velocity of the amplitude decaying is smaller. As shown in Fig.7, the oscillating amplitude is of the order of 400 K when the bubble is under three temperature shocks at Cpc= 0.092. It takes the bubble ten temperature shocks to achieve the above results at Cpc=0.288. Therefore, Cpccan increase not only the intensity of the pressure shock but also the temperature shock. The increase of the pressure difference inside and outside of the bubble drives the liquid do more work on the bubble during the bubble contraction. Moreover, more energies are transformed into the internal energy of the bubble with the increase of the temperature.

    In addition, the maximums of the temperature amplitudes are not observed at the first temperature shock but at the second one in Figs.7(c) and 7(d), similar to the pressure shock. This proves again that not only the increase of the initial-time pressure but also the ambient pressure impacts the bubble fluctuation.

    4. Conclusions

    The present work simulates the bubble collapse based on the Rayleigh-Plesset equation, which includes the thermal and mass diffusion. Moreover, the variation of the ambient pressure is calculated based on the Bernoulli equation combined with the experiment measurements. The results of the simulation show that a transient motion of the bubbles can be obtained in the diffuser, giving rise to the radial formation in the bubble to induce the pressure and temperature effects. The results can be a reference for subsequent studies of cavitation applications in the jet pump. In addition, due to the simple flow geometry of the jet pump, such a system is very easy to design and scale up for large-scale applications. By manipulating the outlet pressure of the flow, one can very easily control the intensity of the cavitation. In summary, the following conclusions are drawn.

    (1) When the bubble is under a high ambient pressure, the bubble radius fluctuates drastically accompanied with a strong pressure and temperature shock. Due to the dissipation of the energy and the diffusion of the contaminant gas, the oscillation of the bubble decays gradually until the bubble radius slightly fluctuating around a constant value. At this time, the main content of the bubble is the contaminant gas.

    (2) The amplitudes of the temperature and the pressure are decided by two factors: the velocity of the ambient pressure increase at the initial time and the value of the whole ambient pressure.

    (3) Moreover, increasing the outlet pressure coefficient leads to a more intense bubble collapsing process, where the amplitudes of the pressure and tempe-rature shock increase significantly and the time interval decreases between the adjacent shocks.

    [1] Balasundaram B., Pandit A. Selective release of invertase by hydrodynamic cavitation [J]. Biochemical Engineering Journal, 2001, 8(3): 251-256.

    [2] Gogate P. R. Cavitation: An auxiliary technique in wastewater treatment schemes [J]. Advances in Environmental Research, 2002, 6(3): 335-358.

    [3] Gogate P. R., Kabadi A. M. A review of applications of cavitation in biochemical engineering/biotechnology [J]. Biochemical Engineering Journal, 2009, 44(1): 60-72.

    [4] Luo X. W., Ji B., Tsujimoto Y. A review of cavitation in hydraulic machinery [J]. Journal of Hydrodynamics, 2016, 28(3): 335-358.

    [5] Peng G., Shimizu S. Progress in numerical simulation of cavitating water jets [J]. Journal of Hydrodynamics, 2013, 25(4): 502-509.

    [6] Ji B., Wang J., Luo X. et al. Numerical simulation of cavitation surge and vortical flows in a diffuser with swirling flow [J]. Journal of Mechanical Science and Technology, 2016, 30(6): 2507-2514.

    [7] Xiao L., Long X. Cavitating flow in annular jet pumps [J]. International Journal of Multiphase Flow, 2015, 71: 116-132.

    [8] Zhang A. M., Cui P., Cui J. et al. Experimental study on bubble dynamics subject to buoyancy [J]. Journal of Fluid Mechanics, 2015, 776: 137-160.

    [9] Zhang A. M., Wang S. P., Huang C. et al. Influences of initial and boundary conditions on underwater explosion bubble dynamics [J]. European Journal of Mechanics. B, Fluids, 2013, 42(2): 69-91.

    [10] Long X., Zhang J., Wang Q. et al. Experimental investigation on the performance of jet pump cavitation reactor at different area ratios [J]. Experimental Thermal and Fluid Science, 2016, 78: 309-321.

    [11] Xiao L., Long X., Zhang J. Shear cavitation in an annular jet pump under recirculation conditions [J]. Journal of Fluids Engineering, 2016, 138(6): 061303.

    [12] Peng X. X., Ji B., Cao Y. et al. Combined experimental observation and numerical simulation of the cloud cavitation with U-type flow structures on hydrofoils [J]. International Journal of Multiphase Flow, 2016, 79: 10-22.

    [13] Zhang A. M., Li S., Cui J. Study on splitting of a toroidal bubble near a rigid boundary [J]. Physics of Fluids, 2015, 27(6): 062102.

    [14] Popinet S., Zaleski S. Bubble collapse near a solid boundary: A numerical study of the influence of viscosity [J]. Journal of Fluid Mechanics, 2002, 464: 137-163.

    [15] Dular M., Coutier-Delgosha O. Thermodynamic effects during growth and collapse of a single cavitation bubble [J]. Journal of Fluid Mechanics, 2013, 736: 44-66.

    [16] Li B. B., Zhang H. C., Lu J. et al. Experimental investigation of the effect of ambient pressure on laser-induced bubble dynamics [J]. Optics and Laser Technology, 2011, 43(8): 1499-1503.

    [17] Zhang Y., Li S. Mass transfer during radial oscillations of gas bubbles in viscoelastic mediums under acoustic excitation [J]. International Journal of Heat and Mass Transfer, 2014, 69(2): 106-116.

    [18] Zhang Y., Du X. Influences of non-uniform pressure field outside bubbles on the propagation of acoustic waves in dilute bubbly liquids [J]. Ultrasonics Sonochemistry, 2015, 26: 119-127.

    [19] Wang Q. X., Manmi K. Three dimensional microbubble dynamics near a wall subject to high intensity ultrasound [J]. Physics of Fluids, 2014, 26(3): 132-143.

    [20] Qin Z. Investigation of the cavitation mechanism and erosion of submerged high pressure water jets [D]. Doctoral Thesis, Brisbane, Australia: The University of Queensland, 2004.

    [21] Qin Z., Bremhorst K., Alehossein H. Simulation of cavitation bubbles in a convergent-divergent nozzle water jet [J]. Journal of Fluid Mechanics, 2007, 573: 1-25.

    [22] Kumar K. S., Moholkar V. S. Conceptual design of a novel hydrodynamic cavitation reactor [J]. Chemical Engineering Science, 2007, 62(10): 2698-2711.

    [23] Alehossein H., Qin Z. Numerical analysis of Rayleigh-Plesset equation for cavitating water jets [J]. International Journal for Numerical Methods in Engineering, 2007, 72(7): 780-807.

    [24] Franc J. P., Michel J. M. Fundamentals of cavitation [J]. Fluid Mechanics and Its Applications, 2005, 76(11): 1-46. [25] Yang W. Y., Cao W., Chung T. S. Applied numerical methods using MATLAB [M]. New York, USA: Wiley, 2005.

    [26] Mathews J. H., Fink K. D. Numerical methods using Matlab [M]. 3rd Edition, Ney Jersey, USA: Pertince-Hall Inc, 2004.

    [27] Kiusalaas J. Numerical methods in engineering with Python [J]. Special Issueon Proceedings of Siggraph, 2010, 29(2): 4183-4209.

    [28] Popinet S., Zaleski S. Bubble collapse near a solid boundary: A numerical study of the influence of viscosity [J]. Journal of Fluid Mechanics, 2002, 464: 137-163.

    [29] Zhang Y. N., Li S. Improved formulas for thermal behavior of oscillating nanobubbles [J]. Journal of Hydrodynamics, 2016, 28(2): 325-328.

    10.1016/S1001-6058(16)60763-1

    October 22, 2015, Revised September 6, 2016)

    * Project supported by the National Natural Science Foundation of China (Grant Nos. 11472197, 51679169 and 51109174), the Fundamental Research Funds for the Central Universities (Grant No. 2042016gf0032).

    Biography:Xin-ping Long (1967-), Male, Ph. D., Professor

    Bin Ji, E-mail: jibin@whu.edu.cn

    猜你喜歡
    吳偉新平
    Transport properties of CrP
    幼兒園里歡樂多
    幼兒園(2021年18期)2021-12-06 02:45:42
    小螞蟻去游玩
    幼兒園(2021年16期)2021-12-06 01:06:48
    Scheme to measure the expectation value of a physical quantity in weak coupling regime?
    老腔唱新歌
    金秋(2021年22期)2021-03-10 07:59:16
    讓蘑菇
    幼兒園(2020年3期)2020-03-27 07:00:07
    劉新平 油畫作品
    你總是給我力量
    當代音樂(2018年11期)2018-06-30 14:43:16
    八字與命運
    特別文摘(2016年24期)2016-12-29 20:44:56
    再次相見
    国产免费一区二区三区四区乱码| av电影中文网址| 中文字幕精品免费在线观看视频| 美女视频免费永久观看网站| av视频免费观看在线观看| 波野结衣二区三区在线| av国产精品久久久久影院| 少妇人妻精品综合一区二区| 国产精品香港三级国产av潘金莲 | 国产av码专区亚洲av| 制服诱惑二区| 日韩成人av中文字幕在线观看| 夜夜骑夜夜射夜夜干| 韩国精品一区二区三区| 亚洲av电影在线进入| av免费观看日本| 久久久久人妻精品一区果冻| 国产亚洲av高清不卡| 精品视频人人做人人爽| 欧美精品一区二区大全| 午夜激情久久久久久久| a 毛片基地| 无遮挡黄片免费观看| 亚洲成色77777| 亚洲一卡2卡3卡4卡5卡精品中文| 国产片内射在线| 国产精品久久久久久精品电影小说| 黑人巨大精品欧美一区二区蜜桃| 久久鲁丝午夜福利片| 青草久久国产| 两个人看的免费小视频| 免费看av在线观看网站| 丁香六月天网| 日韩精品免费视频一区二区三区| 久久毛片免费看一区二区三区| 高清在线视频一区二区三区| 18禁动态无遮挡网站| 亚洲欧美日韩另类电影网站| 亚洲熟女毛片儿| 欧美乱码精品一区二区三区| 亚洲精品日本国产第一区| 久久热在线av| 国产日韩一区二区三区精品不卡| 午夜久久久在线观看| 欧美另类一区| 欧美精品av麻豆av| 成人漫画全彩无遮挡| 亚洲国产精品一区三区| 免费久久久久久久精品成人欧美视频| 久久久精品94久久精品| 亚洲国产精品一区三区| 久久人人爽av亚洲精品天堂| 人人妻人人添人人爽欧美一区卜| 亚洲五月色婷婷综合| 80岁老熟妇乱子伦牲交| 中文欧美无线码| 午夜久久久在线观看| 一区二区av电影网| 亚洲精品成人av观看孕妇| 狂野欧美激情性bbbbbb| 水蜜桃什么品种好| 亚洲精品第二区| 免费黄色在线免费观看| 热re99久久精品国产66热6| 国产日韩欧美亚洲二区| videos熟女内射| 天天躁夜夜躁狠狠躁躁| tube8黄色片| 亚洲欧美成人精品一区二区| 欧美日本中文国产一区发布| 一区二区av电影网| 国产老妇伦熟女老妇高清| 搡老岳熟女国产| 亚洲一码二码三码区别大吗| 久久影院123| 韩国高清视频一区二区三区| 女性生殖器流出的白浆| 一级片'在线观看视频| 黄色视频在线播放观看不卡| 亚洲成av片中文字幕在线观看| 人成视频在线观看免费观看| 久久热在线av| 三上悠亚av全集在线观看| 啦啦啦 在线观看视频| 欧美日韩福利视频一区二区| 我要看黄色一级片免费的| 亚洲国产毛片av蜜桃av| 男的添女的下面高潮视频| 狂野欧美激情性xxxx| 不卡av一区二区三区| 欧美97在线视频| 美女大奶头黄色视频| 欧美国产精品一级二级三级| 老司机靠b影院| 在线天堂中文资源库| 操出白浆在线播放| 欧美国产精品va在线观看不卡| 免费黄频网站在线观看国产| 欧美亚洲 丝袜 人妻 在线| 两个人看的免费小视频| 999精品在线视频| 女人久久www免费人成看片| 亚洲av成人精品一二三区| 热re99久久国产66热| 精品亚洲乱码少妇综合久久| 母亲3免费完整高清在线观看| 高清av免费在线| 成年人午夜在线观看视频| 精品一品国产午夜福利视频| 欧美精品一区二区免费开放| 我要看黄色一级片免费的| 日韩一区二区视频免费看| 亚洲精品日韩在线中文字幕| 久久精品久久精品一区二区三区| 亚洲精品一区蜜桃| 天堂中文最新版在线下载| 少妇的丰满在线观看| a 毛片基地| 亚洲第一av免费看| 午夜福利在线免费观看网站| 熟女少妇亚洲综合色aaa.| 久久国产亚洲av麻豆专区| 婷婷成人精品国产| 精品免费久久久久久久清纯 | 黄色怎么调成土黄色| 三上悠亚av全集在线观看| 欧美最新免费一区二区三区| 免费不卡黄色视频| 2021少妇久久久久久久久久久| 国产精品久久久av美女十八| 街头女战士在线观看网站| 卡戴珊不雅视频在线播放| 新久久久久国产一级毛片| 男女免费视频国产| 中国国产av一级| 国产欧美日韩一区二区三区在线| 成年人免费黄色播放视频| 日韩一卡2卡3卡4卡2021年| 99久国产av精品国产电影| 国产淫语在线视频| 哪个播放器可以免费观看大片| 午夜福利视频精品| 日韩一区二区视频免费看| 91国产中文字幕| 亚洲,一卡二卡三卡| 伊人久久国产一区二区| 日韩大码丰满熟妇| 国产成人免费无遮挡视频| 久久久国产欧美日韩av| 国产精品成人在线| 不卡av一区二区三区| 1024香蕉在线观看| 国产成人欧美在线观看 | 久热这里只有精品99| 国产一区二区三区av在线| 卡戴珊不雅视频在线播放| 国产日韩欧美在线精品| av免费观看日本| 最近手机中文字幕大全| 久热这里只有精品99| 日本色播在线视频| 91精品伊人久久大香线蕉| 免费少妇av软件| 日韩大码丰满熟妇| 国产精品久久久久成人av| 免费看不卡的av| 欧美日韩福利视频一区二区| 在现免费观看毛片| 黄色一级大片看看| 黑人猛操日本美女一级片| 美女视频免费永久观看网站| 在线天堂最新版资源| 欧美国产精品一级二级三级| 涩涩av久久男人的天堂| 两性夫妻黄色片| 一区二区三区乱码不卡18| 欧美变态另类bdsm刘玥| 成人国语在线视频| 飞空精品影院首页| 纵有疾风起免费观看全集完整版| 国产99久久九九免费精品| 欧美av亚洲av综合av国产av | 亚洲精品中文字幕在线视频| 成人毛片60女人毛片免费| 久久影院123| 97在线人人人人妻| 日本av免费视频播放| 咕卡用的链子| 99国产精品免费福利视频| 操美女的视频在线观看| 天堂俺去俺来也www色官网| 久久婷婷青草| 欧美xxⅹ黑人| 亚洲激情五月婷婷啪啪| 狂野欧美激情性xxxx| 欧美乱码精品一区二区三区| 精品国产超薄肉色丝袜足j| tube8黄色片| 黄片小视频在线播放| 中文字幕另类日韩欧美亚洲嫩草| a级毛片在线看网站| 国产极品粉嫩免费观看在线| 久久ye,这里只有精品| 中文字幕高清在线视频| 美女视频免费永久观看网站| 国产av一区二区精品久久| 中文字幕人妻丝袜制服| 亚洲成国产人片在线观看| 老司机亚洲免费影院| 高清视频免费观看一区二区| 2018国产大陆天天弄谢| 乱人伦中国视频| 一级,二级,三级黄色视频| 国产成人免费无遮挡视频| 午夜免费观看性视频| 国产极品粉嫩免费观看在线| 久热爱精品视频在线9| 美女高潮到喷水免费观看| 国产亚洲欧美精品永久| 国产乱人偷精品视频| 一级,二级,三级黄色视频| 大片免费播放器 马上看| 日韩欧美一区视频在线观看| 十分钟在线观看高清视频www| 亚洲五月色婷婷综合| 最近中文字幕2019免费版| 久久久精品94久久精品| 国产免费又黄又爽又色| 亚洲国产av影院在线观看| 少妇人妻精品综合一区二区| 香蕉丝袜av| 青青草视频在线视频观看| 国产一区有黄有色的免费视频| 最新的欧美精品一区二区| 女人爽到高潮嗷嗷叫在线视频| 午夜福利乱码中文字幕| av在线观看视频网站免费| 一级毛片我不卡| kizo精华| 国产男女内射视频| videos熟女内射| 如日韩欧美国产精品一区二区三区| 日韩免费高清中文字幕av| 亚洲男人天堂网一区| 国产男人的电影天堂91| 亚洲精品一区蜜桃| 黄片无遮挡物在线观看| 久久久久国产精品人妻一区二区| 最近中文字幕2019免费版| 9色porny在线观看| 久热这里只有精品99| e午夜精品久久久久久久| 秋霞在线观看毛片| 电影成人av| 精品国产国语对白av| 宅男免费午夜| 免费久久久久久久精品成人欧美视频| 十八禁高潮呻吟视频| 如日韩欧美国产精品一区二区三区| 精品少妇内射三级| 日本91视频免费播放| 欧美精品一区二区大全| 激情五月婷婷亚洲| 777久久人妻少妇嫩草av网站| 丁香六月天网| 国产av精品麻豆| 欧美国产精品一级二级三级| 欧美日韩av久久| 久久av网站| 精品国产乱码久久久久久小说| 两个人看的免费小视频| 久久久久久人人人人人| 色网站视频免费| 亚洲精品乱久久久久久| 啦啦啦 在线观看视频| 狂野欧美激情性xxxx| 在线 av 中文字幕| 久久狼人影院| 亚洲一区二区三区欧美精品| 中文精品一卡2卡3卡4更新| 欧美人与性动交α欧美软件| 国产亚洲精品第一综合不卡| 久久久国产精品麻豆| 99久久精品国产亚洲精品| 一级毛片我不卡| 亚洲成人免费av在线播放| 国产成人精品福利久久| 亚洲伊人久久精品综合| 咕卡用的链子| 男女午夜视频在线观看| 久久女婷五月综合色啪小说| 中文字幕人妻熟女乱码| 中文字幕av电影在线播放| 午夜激情av网站| 男人添女人高潮全过程视频| 一级毛片黄色毛片免费观看视频| 免费在线观看黄色视频的| 久久久久久久久免费视频了| 欧美日韩亚洲高清精品| 午夜福利在线免费观看网站| 中文精品一卡2卡3卡4更新| 久久久精品区二区三区| 亚洲精品视频女| av不卡在线播放| 九草在线视频观看| 日本爱情动作片www.在线观看| 亚洲国产看品久久| 自拍欧美九色日韩亚洲蝌蚪91| 午夜精品国产一区二区电影| av免费观看日本| 欧美成人午夜精品| 18在线观看网站| 水蜜桃什么品种好| 男女边摸边吃奶| 亚洲精品日本国产第一区| 国产av码专区亚洲av| 中文字幕高清在线视频| 国产精品一区二区在线观看99| 成人18禁高潮啪啪吃奶动态图| 巨乳人妻的诱惑在线观看| 菩萨蛮人人尽说江南好唐韦庄| 精品国产超薄肉色丝袜足j| 日韩av在线免费看完整版不卡| av电影中文网址| 久久女婷五月综合色啪小说| 一二三四中文在线观看免费高清| 午夜福利,免费看| 成人影院久久| 欧美黄色片欧美黄色片| 制服人妻中文乱码| 尾随美女入室| 亚洲情色 制服丝袜| 午夜免费男女啪啪视频观看| 黑人猛操日本美女一级片| 日韩精品免费视频一区二区三区| 精品午夜福利在线看| 女人高潮潮喷娇喘18禁视频| 色吧在线观看| 国产亚洲av片在线观看秒播厂| 日韩伦理黄色片| 黄色 视频免费看| 中文欧美无线码| 亚洲熟女毛片儿| 婷婷色麻豆天堂久久| 欧美精品亚洲一区二区| 国产精品 国内视频| 老司机影院成人| 女人久久www免费人成看片| 国产精品欧美亚洲77777| 欧美亚洲日本最大视频资源| 国产成人精品在线电影| 自拍欧美九色日韩亚洲蝌蚪91| 丝袜脚勾引网站| 精品一区二区三卡| 天天躁夜夜躁狠狠久久av| 欧美另类一区| 国产精品久久久久久久久免| 国产男女超爽视频在线观看| 日韩中文字幕视频在线看片| 亚洲精品中文字幕在线视频| 国产伦理片在线播放av一区| 亚洲欧美成人综合另类久久久| 综合色丁香网| 亚洲伊人色综图| 国产成人精品久久久久久| 青春草亚洲视频在线观看| 亚洲精品一区蜜桃| 色婷婷av一区二区三区视频| 国产一区有黄有色的免费视频| 99热国产这里只有精品6| 中文字幕亚洲精品专区| 亚洲欧美成人综合另类久久久| 欧美成人精品欧美一级黄| 两性夫妻黄色片| 国产无遮挡羞羞视频在线观看| 麻豆av在线久日| 国产精品一区二区精品视频观看| 欧美精品人与动牲交sv欧美| 韩国av在线不卡| 免费高清在线观看视频在线观看| 麻豆精品久久久久久蜜桃| 免费人妻精品一区二区三区视频| 国产亚洲av高清不卡| 国产极品天堂在线| 亚洲精品中文字幕在线视频| 又大又黄又爽视频免费| 日日啪夜夜爽| 最近最新中文字幕免费大全7| 九九爱精品视频在线观看| 精品国产乱码久久久久久小说| 国产一级毛片在线| 看非洲黑人一级黄片| 伊人久久大香线蕉亚洲五| 日韩不卡一区二区三区视频在线| 久久久久久人妻| 日韩不卡一区二区三区视频在线| 久久人人97超碰香蕉20202| 免费观看a级毛片全部| 国产精品久久久久久久久免| 国产一区二区在线观看av| 精品亚洲成国产av| 亚洲综合色网址| 亚洲精品久久成人aⅴ小说| 欧美中文综合在线视频| 国产成人午夜福利电影在线观看| 亚洲国产精品一区三区| 国产欧美亚洲国产| 国产1区2区3区精品| 国产高清不卡午夜福利| 老司机影院成人| 久久久亚洲精品成人影院| 亚洲五月色婷婷综合| 两性夫妻黄色片| 中文字幕亚洲精品专区| 桃花免费在线播放| 美女中出高潮动态图| 国产精品 欧美亚洲| 国产又色又爽无遮挡免| 亚洲色图综合在线观看| 久久国产精品大桥未久av| 免费在线观看视频国产中文字幕亚洲 | 老司机亚洲免费影院| 亚洲国产欧美一区二区综合| 超碰97精品在线观看| 欧美人与性动交α欧美精品济南到| 另类精品久久| 视频区图区小说| 国产乱人偷精品视频| 狂野欧美激情性xxxx| h视频一区二区三区| 成人18禁高潮啪啪吃奶动态图| 午夜福利,免费看| 亚洲精品av麻豆狂野| 最近中文字幕2019免费版| 日本欧美国产在线视频| 男女边摸边吃奶| 免费看av在线观看网站| 日韩一卡2卡3卡4卡2021年| 国产精品麻豆人妻色哟哟久久| 日韩中文字幕视频在线看片| 久久天堂一区二区三区四区| 免费观看性生交大片5| 午夜福利免费观看在线| 亚洲四区av| 亚洲精品视频女| 黄色 视频免费看| 性少妇av在线| 日韩电影二区| 人人妻人人添人人爽欧美一区卜| 女性生殖器流出的白浆| 免费日韩欧美在线观看| 超色免费av| 如何舔出高潮| 久久久久视频综合| 日本wwww免费看| 九九爱精品视频在线观看| 国产97色在线日韩免费| 秋霞伦理黄片| av国产精品久久久久影院| 最近中文字幕2019免费版| 曰老女人黄片| 黄片无遮挡物在线观看| 麻豆av在线久日| 亚洲精品在线美女| 国产精品秋霞免费鲁丝片| 巨乳人妻的诱惑在线观看| 国产人伦9x9x在线观看| 久久久国产精品麻豆| 在现免费观看毛片| 久久精品亚洲av国产电影网| 伊人久久国产一区二区| 免费av中文字幕在线| 精品人妻熟女毛片av久久网站| 黄色毛片三级朝国网站| 最近最新中文字幕大全免费视频 | av视频免费观看在线观看| 天美传媒精品一区二区| 国产av精品麻豆| 欧美精品亚洲一区二区| 欧美日韩福利视频一区二区| 国产又爽黄色视频| 国产日韩一区二区三区精品不卡| 黄色视频不卡| 老司机深夜福利视频在线观看 | 免费观看人在逋| 九色亚洲精品在线播放| 免费观看a级毛片全部| 男女午夜视频在线观看| 国产视频首页在线观看| 久久女婷五月综合色啪小说| 国产精品99久久99久久久不卡 | 丁香六月天网| 成人影院久久| 亚洲人成77777在线视频| 在线 av 中文字幕| 亚洲色图 男人天堂 中文字幕| 九九爱精品视频在线观看| 最新在线观看一区二区三区 | 亚洲欧美一区二区三区国产| 最近最新中文字幕大全免费视频 | 一级毛片我不卡| 黄色一级大片看看| 久久久久久久精品精品| 国产伦理片在线播放av一区| √禁漫天堂资源中文www| 极品少妇高潮喷水抽搐| av福利片在线| 久久久久久久大尺度免费视频| 最新的欧美精品一区二区| 国产精品久久久久久久久免| 精品久久蜜臀av无| 亚洲,欧美,日韩| 亚洲精品国产av成人精品| 亚洲精品久久久久久婷婷小说| 国产精品蜜桃在线观看| 日韩中文字幕视频在线看片| 欧美黑人欧美精品刺激| 亚洲综合精品二区| 国产黄频视频在线观看| 无遮挡黄片免费观看| 丰满迷人的少妇在线观看| 亚洲人成77777在线视频| 一二三四在线观看免费中文在| 欧美精品一区二区免费开放| 欧美精品高潮呻吟av久久| 欧美久久黑人一区二区| 在线天堂最新版资源| 国产精品一二三区在线看| 色婷婷av一区二区三区视频| 亚洲欧美成人精品一区二区| 精品久久久精品久久久| 91国产中文字幕| 国产亚洲av片在线观看秒播厂| svipshipincom国产片| www.av在线官网国产| 久久精品亚洲熟妇少妇任你| 国产黄频视频在线观看| 欧美成人精品欧美一级黄| 少妇被粗大猛烈的视频| 嫩草影视91久久| 免费观看a级毛片全部| 如何舔出高潮| 1024视频免费在线观看| 一级a爱视频在线免费观看| 午夜日本视频在线| 侵犯人妻中文字幕一二三四区| 亚洲国产精品一区二区三区在线| 精品久久久精品久久久| av.在线天堂| 中文字幕精品免费在线观看视频| 叶爱在线成人免费视频播放| 我要看黄色一级片免费的| 亚洲国产精品一区二区三区在线| 美国免费a级毛片| 免费观看人在逋| 中文乱码字字幕精品一区二区三区| 免费观看av网站的网址| 曰老女人黄片| 视频在线观看一区二区三区| 亚洲欧美日韩另类电影网站| 成人手机av| 三上悠亚av全集在线观看| 精品久久久久久电影网| 伊人久久大香线蕉亚洲五| 色婷婷av一区二区三区视频| 国产成人系列免费观看| 久久久久人妻精品一区果冻| 国产伦人伦偷精品视频| 女的被弄到高潮叫床怎么办| 男女边吃奶边做爰视频| 伦理电影大哥的女人| 亚洲欧美日韩另类电影网站| 亚洲一级一片aⅴ在线观看| 亚洲精品美女久久久久99蜜臀 | 日本91视频免费播放| 校园人妻丝袜中文字幕| 亚洲欧美一区二区三区黑人| 中文字幕亚洲精品专区| 男女国产视频网站| 你懂的网址亚洲精品在线观看| 国产亚洲av高清不卡| 一本色道久久久久久精品综合| 宅男免费午夜| 老司机亚洲免费影院| 视频区图区小说| 欧美激情 高清一区二区三区| 欧美乱码精品一区二区三区| 亚洲熟女精品中文字幕| 欧美黑人精品巨大| 一级毛片我不卡| 亚洲色图 男人天堂 中文字幕| 午夜影院在线不卡| 搡老乐熟女国产| 国产免费福利视频在线观看| 999久久久国产精品视频| 国产人伦9x9x在线观看| 欧美xxⅹ黑人| 三上悠亚av全集在线观看| 久久ye,这里只有精品| 精品免费久久久久久久清纯 | 亚洲成色77777| 黄网站色视频无遮挡免费观看| 亚洲一级一片aⅴ在线观看| 国产一区二区三区综合在线观看| 国产片特级美女逼逼视频| 国产精品秋霞免费鲁丝片| 人人妻人人爽人人添夜夜欢视频| 超碰成人久久| 日韩不卡一区二区三区视频在线| 久久久精品免费免费高清| 婷婷色综合www| 又大又黄又爽视频免费| 久久av网站| 午夜激情av网站| 一本—道久久a久久精品蜜桃钙片|