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

    EFFECT OF A PROPELLER AND GAS DIFFUSION ON BUBBLE NUCLEI DISTRIBUTION IN A LIQUID*

    2012-08-22 08:32:14HSIAOChaoTsungCHAHINEGeorges
    水動力學研究與進展 B輯 2012年6期

    HSIAO Chao-Tsung, CHAHINE Georges L.

    DYNAFLOW, INC., Jessup, Maryland 20794, USA, E-mail: ctsung@dynaflow-inc.com

    (Received December 26, 2011, Revised April 28, 2012)

    EFFECT OF A PROPELLER AND GAS DIFFUSION ON BUBBLE NUCLEI DISTRIBUTION IN A LIQUID*

    HSIAO Chao-Tsung, CHAHINE Georges L.

    DYNAFLOW, INC., Jessup, Maryland 20794, USA, E-mail: ctsung@dynaflow-inc.com

    (Received December 26, 2011, Revised April 28, 2012)

    A multi-bubble dynamics code accounting for gas diffusion in the liquid and through the bubble wall was developed and used to study the modification of a bubble nuclei population dynamics by a propeller. The propeller flow field was obtained using a Reynolds-Averaged Navier-Stokes (RANS) solver and bubble nuclei populations were propagated in this field. The numerical procedure enabled establishment of the possibility of production behind the propeller of relatively large visible bubbles starting from typical ocean nuclei size distributions. The resulting larger bubbles are seen to cluster in the blade wakes and tip vortices. Parametric investigations of the initial nuclei size distribution, the dissolved gas concentration, and the cavitation number were conducted to identify their effects on bubble entrainment and the resultant void fractions and bubble distribution modifications downstream from the propeller. Imposed synthetic turbulence-like fluctuations unto the average RANS flow field were also used to study the effect averaging in the RANS procedure has on the results.

    gas diffusion, bubble entrainment, propeller

    Introduction

    The generation and entrainment of bubbles by ship motion raises many two phase flow fundamental questions and is of interest as bubble generation could in hostile environments affect the safety of the ship. These bubbles are transported by the flow field to the stern area, captured in the ship wake, trapped in the large vortical structures, and then act as tracers of the ship wake. In littoral warfare, such bubbly wakes provide an excellent opportunity to homing devices enabling them to easily find their target because of the large acoustic cross sections and high acoustic response of bubbles to acoustic waves.

    It is believed that bubble entrainment has two main sources. One source is from bubbles generated by air entrainment at the surface due to free surface and hull interactions resulting in breaking and spilling waves at the bow and stern and bubble entrainment in the boundary layer along the hull. The other potential source is bubble “production” by the propellers. While there have been some research efforts to understandand observe bubble generation and entrainment in breaking waves[1-4], much less has been done to investigate the second means of bubble generation/entrainment due to the propeller[5], which is investigated here.

    One hypothesized scenario for generation and entrainment of bubbles by the propeller is rooted in the fact that natural waters always contain suspended microscopic bubble nuclei in addition to dissolved gas. These nuclei, when subjected to variations in the local liquid pressures, will respond dynamically by changing volume, oscillating, and eventually growing from sub-visual to become visible due to cumulative gas transfer into the bubbles. Several techniques have been used to measure the distribution and sizes of these nuclei both in the ocean and in laboratories. These include Coulter counter, holography, light scattering methods, cavitation susceptibility meters and acoustic methods such as the ABS ACOUSTIC BUBBLE SPECTROMETER?[6]. Franklin[7]had summarized some earlier nuclei size distribution measurements done by different researchers in different waters.

    In order to investigate this type of bubble size enhancement, we have considered a numerical approach combining propeller flow filed viscous solution coupled with multi-bubble dynamics incorporation mode-ling of gas diffusion to account for the transfer of gases across the bubble walls. Such an approach has been summarized recently in Chahine[8]. 3DYNAFSDSM?is a Lagrangian multi-bubble tracking and dynamics code[9]using a bubble motion equation with a modified Rayleigh-Plesset bubble dynamics equation incorporating bubble Surface Averaged Pressures (SAP) to account for flow field non-uniformities, and a bubble slip velocity pressure term[10]. This numerical model allows consideration of the spatiotemporal variations of dissolved gas concentration around the bubble and enables full tracking of a realistic bubble nuclei population.

    To investigate bubble “production” by a propeller, we have utilized this model to propagate a nuclei distribution in the flow field of a marine propeller (NSWCCD Propeller 5168). This flow field was simulated earlier by Hsiao and Pauley[11]using a Reynolds-Averaged Navier-Stokes (RANS) solver and validated with experimental measurements.

    1. Numerical approach

    1.1 Bubble dynamics model

    The dynamics of a spherical bubble has been extensively studied following the original works of Rayleigh[12]and Plesset[13]which included for an incompressible liquid, the effects of inertia, bubble content compressibility, and ambient pressure variations. Since then, a very large number of studies have been conducted to include a host of other physical phenomena. In the present study we have considered the following form of the equation for the bubble radius, R(t ), which accounts for liquid and gas compressibility, liquid viscosity, surface tension, and non-uniform pressure fields, and is based on Keller-Herring equation[14]. where c is the sound speed in the liquid,ρ its density, μ its viscosity, and pvits vapor pressure. pgis the pressure of the gas in the bubble, γ is the surface tension, and ubis the bubble travel velocity. In Eq.(1), which we dub the SAPbubble dynamics equation[9], we have accountedfor a slip velocity between the bubble and the host liquid, and for a nonuniform pressure field along the bubble surface.pencand uencare the liquid pressure and velocity respectively.pencand uencare defined as the averages of the liquid pressures and velocities over the bubble surface. The use ofpencresults in a major improvement over the classical spherical bubble model which uses the pressure at the bubble center in its absence. The gas pressure, pg, is obtained, as described in the next section, from the solution of the gas diffusion problem and the assumption that the gas is an ideal gas.

    The bubble trajectory is obtained using the following motion equation

    where ρbis the air density, CDis the drag coefficient given by an empirical equation from Haberman and Morton[15]and ? is the deformation tensor. The 1st right hand side term is a drag force. The 2nd and 3rd terms account for the added mass. The 4th term accounts for the presence of a pressure gradient, while the 5th term accounts for gravity and 6th term is a lift force.

    1.2 Gas diffusion model

    Water can contain gas not only in the form of nuclei with a particular distribution, but also as a gas dissolved in the liquid with a concentration C Over time, dissolved gas will diffuse from high concentration regions into low concentration regions following a gas transport equation for the time and space dependent dissolved gas concentration, C(x,t ), given by

    wheregD is the molar diffusivity of the gaseous component in the liquid (in practice the turbulent diffusivity in high turbulence areas).

    For an oscillating bubble in a liquid containing dissolved gas, Eq.(3) needs to be solved while satisfying the following initial and far field boundary conditions:

    where ris the radial distance from the bubble center and C∞is the dissolved gas concentration far away from the bubble surface.

    At the bubble surface Henry’s law, which relates the concentration of gas in a liquid to the partial pressureof the gas above the liquid, applies. Hence, we have at r=R

    whereH is the Henry constant.

    As the bubble oscillates, non-condensable gas is transported across the bubble wall. The rate of transport of gas moles into the bubble, n˙g, can be related to the gradient of gas concentration in the liquid at the bubble wall

    where Sis thesurface of the bubble, ?C/?n is the normal derivative to the bubble surface of the liquid gas concentration. The integration on the handside of Eq.(6) is over the bubble surface. Time integration of Eq.(6) determines at every instant the total number of moles of gas, ng, in the bubble.

    The two components of the bubble content: vapor and gas, are both assumed to b e ideal gases which follo w an ideal gas law

    where Vbis the volume o f bubble, pvis the vapor pressure,nvis number of moles of vapor within the bubble, Ruis the universal gas constant and Tbis the absolutetemperature of the gas and vapor mixture. Due to therelatively short vaporization time compared to bubble dynamics and gas diffusion characteristic times, the vapor is considered to instantaneously flow in and out of the bubble, and pvis taken equal to the equilibrium vapor pressure of the liquid at the bubble wall temperature. One consequence of this assumption is that the amount ofng, and nvin the bubble are directly proportional to the ratio of their respective partial pressures.

    Use of the idea gas equation of state introduces another unknown Tb. This problem can be closed by considering the thermodynamics of the contents of the bubble using the first law energy balance of the ideal gas mixture within the bubble. We consider the bubble wall to constitute a deformable and permeable control surface and write the first law energy balance for the control volume bounded by that surface

    where dUis the change in internal energy, dW is the work done on the control volume, hiis the specific enthalpy, cVis the specific heat at constant volume, cPis the specific heat at constant pressure and Tiis the liquid temperature.

    Combining Eq.(8) and (9) one can obtain

    where the superscripts Tband Tlfor cVand cPindicate that the specificheats are evaluated at corresponding bubble or liquid temperature. Using the fact that cP-cV=Ruone can rearrange Eq.(10) to become:

    Integration of Eq.(11) provides the instantaneous gas pressure to be used in Eq.(1) and the boundary condition, Eq.(5).

    A general, detailed solution of Eq.(3) would involve a time consuming numerical procedure, such as a space and time dependent finite difference scheme.To avoid this, an d since we want to consider very large number of nuclei/bubbles, we adopt the thin boundary layer approximation introduced by Plesset and Zwick[16]to obtain a solution to Eq.(3) for the case of an isolated spherical bubble. In this approach large gradients in gas concentration are concentrated in a bubble wall boundary layer which is small compared to the bubble radius. An analytical solution then exists and relates the gas concentration at the bubble wall to the concentration at “infinity”. This expression requires integration over the whole history of the bubble dynamics in order to enable computation of the amount of gas inside the bubble, and has the following form

    Although solving Eq.(12) instead of Eq.(3)reduces numerical complexity, some loss of solution generality is expected due to the thin layer approximation and the neglect of convection term.

    Equation (11) requiresgn andto be integrated in time. To do so we rewrite Eq.(12) by first assigning the time interval integration part to be

    The interval time integration[0,t] can be split into two intervals to separate out the current time step as

    where Δt is the time step size. The second terms in Eq.(14) is evaluated analytically to avoid the (tτ)-1/2sing ularity att=τ

    where the tildes denote average values over the interval [t-Δt,t]. If we approximatethese average vales by≈(t)and≈[R(t)+R(t-Δt )]/2, then Eq.(15) becomes

    After substituting Eq.(17) and (13) into Eq.(15), we can rearrange Eq.(15) to become

    To determine each bubble motion andvolume variation, the set of four differential Eq.(1), (2), (12) and (18) are solved using a Runge-Kutta fourth-order scheme to integrate through time.

    Fig.1 Comparison of bubble radius versus time and resulting bubble growth from rectified diffusion between the current numerical model and Crum and Mao’s analytical model[17]

    The current numerical model is compared to the analytical model derived by Crum and Mao[17]for a rectification diffusion bubble driven by an acoustic pressure field as shown in Fig.1. It is seen that the comparison between our current numerical model and Crum’s analytical model shows a good agreement on the bubble growth rate except at the very beginning. Also, due to the neglect of all but the first-order term in the Crum analytical model, the bubble maxima andminima of the analytical model are seen to be a few percent smaller than in the present numerical model.

    Fig.2 Axial velocity Vxfrom experiment and computationat x/Rp=0.2386

    Fig.3 Radial velocity Vrfrom experiment and computation at x/Rp=0.2386

    2. Results and discussion

    The above method was applied to study nuclei motion and size variations in a selected propeller flow field. Nuclei in natural water with a prescribed air void fraction were releasedupstream of a rotating propell?er and their dynamics computed using 3DYNAFSVIS. The propeller model was the David Taylor Propeller 5 168 propeller which is a five-blade propeller with a 15.856 inch (0.4 m) diameter. The flow field around this propeller was simulated earlier for three different advance coefficients by Hsiao and Pauley[11]using a RANS solver and were validated against experimental measurements conducted by Chesnakas and Jessup[18]. Figure 2 and Fig.3 illustrate the good comparison between simulation and measurement for the axial and radial components of velocity field. Additional more detailed comparisons can be found in Hsiao and Pauley[11]. The propeller flow field condition used in the study corresponded to an advance coefficient J=U∞/nD =1.1, where the liquid velocity U∞=10.7 m/s, the propeller number of revolutions per minutes, n=1 450 rpm, and D=0.4 m the propeller diameter.

    Fig.4 Release location of the nuclei upstream from the rotating propeller for the “synthetic” case study presented here to illustrated effect of propeller on nuclei space and size distribution

    3. Gas diffusion effects

    To study t he gas diffusion effects on the bubbledynamics and the size distribution downstream from the propeller, we compare the results for a “synthetic”bubble nuclei distributionof including or not gas diffusion effects. Nuclei with a uniform radius of 50 μm were released from a preset grid located upstream from the propeller at x=–0.08 m as illustrated in Fig.4. The total number of nuclei released was 4 500 and the nuclei were tracked from x=–0.08 m to x= 0.3 m. For the computations with gas diffusion, the water was assumed to be saturated with gas (i.e., the dissolved gas concentration was 100%, i.e., C= 0.66 mol/m3).

    Fig.5 Bubble size and location distribution at theexit of the computational domain for simulations without gas diffusion at the exit from the computational domain,x= 0.22 m

    Figure 5 and Fig.6 show the resulting bubble size andspatial distributions at x=0.22 m. for thetwo conditions where gas diffusion was not taken into account (Fig.5) and where itwas (Fig.6). These simulations were conducted at a cavitation number, σ= 1.75, where σ is defined as

    where p∞and pvare the liquid pressure and velocity upstreamfrom the propeller.

    From the comparison of the bubble size distributions it is seen that the propeller modifies the size range of the bubbles from the uniform 50 μm nuclei sizes to a range of sizes: 50μm<R<70μm when therewas no gas diffusion. When gas diffusion is taken into account a much larger range of bubble sizes is generated: 50μm<R<190μm . In both cases larger sized bubbles are seen to coll ect mainly in two regions: the tip vortices areas and the wakes of the blades, which emanate from its blade suction side.

    Fig.6 Bubble size and location distribution at theexit of the computational domain for simulations withgas diffusion at the exit from the computational domain, x=0.22 m

    By observing the behavior of individual bubbles, we are able to determine where the larger bubble originated from. Figure 7 shows three typical scenarios of bubble paths through the propeller flow field.

    (1) In the first (Fig.7(a)) the bubble has a helicoidal motion and gets entrapped in a tip vortex. It grows to its maximum size while approaching the vortex core and feeling the lowest pressure region. The bubble then shrink in volume while oscillating as the pressure along the vortex center recovers. However, the bubble does not return to its initial size (50 μm), instead it retains a much larger size (80 μm) due to a net increase of the amount of gas in the bubble accumulated during the volume fluctuations.

    (2) In the second scenario, the bubble travels over the blade surface on the suction side (Fig.7(b)). Here the bubble grows to its maximum size as it encounters the minimum pressure on the blade surface. It then collapses as the pressure recovers, but here again retains a much larger size (110 μm) than the original after it passes the blade area.

    (3) The third scenario is the most common, and concerns bubbles passing through the blade areas away from the tip vortices and blade suction side regions. For these bubbles the size remains very close to the original size during and after crossing the prope-ller f low field, as shown in Fig.7(c).

    Fig.7 Bubble trajectories along the propeller and bubble radius along the trajectories

    Fig.8 Illustration of the fictitious volume used to feed nuclei into the inlet (or release area) of the computational domain

    4. Modeling of a real nuclei field

    In order to simulate realistic water conditions a known size distribution of nuclei was selected and used toseed the water feeding the propeller. The nuclei were considered to be distributed randomly in a fictitious supply volume feeding the computational domain inlet plane (release area). The fictitious volume size was determined as the product of the release area by the sought physical duration of the simulation and the inlet velocity, U∞, as illustrated in Fig.8. The nuclei size density distribution function, n(R), is defined as the number of nuclei per unit volume having radii in the range[R,R+dR]

    where N(R) is the number of nuclei of radiu s R in aunit volume. This function can be expressed as a discrete distribution of M selected nuclei sizes. Thus, the total void fractio n, α, in the liquid can be obtained from

    where Niis the discrete number of nuclei of radius Riused in the computations.

    Fig.9(a)Thenuclei size density distribution selected to resemblethefield measurement byMedwin[7]

    Fig.9(b) The resultant number of nuclei released

    For the results shown below we have selected the following characteristic values, which correspond to typical field measurement by Medwin[7]size range of 10 μm to 200 μm and void fraction, α=3.23×10–5. Figure 9(a) shows the selected nuclei size numberdensi ty and Fig.9(b) shows the numbers of discrete bubble sizes and bandwidths selected.

    5. Time-Averaged void fraction distribution

    In order to minimize computations CPU time, a time of release, Δt, during which actual bubble compu tations are conducted was selected and correspond ed to a meaningful bubble population distribution in the fictitious bubble supply volume. The resulting bubble behavior was then assumed to be repeated for mΔt and the effects of m on the results were analyzed.

    Fig.10Effect of the bubble release duration on the total number of bubbles within the computational domain. Computations were conducted for Δt=3.7×10–3s and results were then repeated to cover much longer release times, mΔt,m=1-25

    Fig.11 Snapshot of the bubbles in the propeller flow field of the bubbles in the computational domain ata given time during the computation. Bubble sizes are to-scale in volume II, they are enhanced by a factor of5 in the volumes I and III

    Figure 10 shows the time variation of the total number of bubbles within the computational domain (-0.08m<0.3m) for Δt=3.7×10–3s and for different number of repeats, m. The results clearly indicate that a steady sta teis achieved after t>0.03 ms whena minimum of 10 repeats is effectuated. For a smaller value of the release duration the solution is only transient. Increasing the number of repeats, m, over 10 extends the quasi-steady state period in the computational domain.

    Figure 11 illustrates the sizes and locations of the bubbles on the propeller blades and in the tip vortices. A snapshot of both front and side views are shown for σ=1.75. It is noted that in the side view the bubbles in the volume II are shown to scale, while they are amplified 5 times (to make them more “visible” in print) in the volumes I and III. Figure 11 illustrates the marked increase in “visible” bubbles downstream of thepropeller, traveling bubble cavitation, tip vortex cavitation and the initiation of sheet cavitation.

    Fig.12 The distribution of the time-averaged void fraction,, along the propeller axialdirectionfor avalueof void fraction at the inlet, α=3.23×10–5and α=1.75

    In order to illustrate the results in a simp lified manner, some average quantities are introduced. We can define a section-averaged time dependent void fraction, αΔx(x,t), as follows

    where d is the shaft diameter,D, the propeller diameter,and N(x,t)is the number of bubble within the interval [x,x+Δx]at thegiven time, t. The time-averaged void fraction ofα(x,t)att1, over a time period dt can also be defined as

    Figure 12 shows the distributionofαalong x. It is seen thatthe void fraction increases significantly as a result ofcavitation on the blades (first hump in the figure) and in the tip vortices (second hump). The void fraction decreases asthe pressure encountered by the bubble increases. However, the downstream void fraction is in this case still becomes 1.4 times that of upstream.

    Fig.13 Comparison of the initial nuclei size distributionentering the computational domain and that resulting at the domain exit at x=0.2 m

    The effect of propeller cavitation on the number of bubbles having a large size is much more important than on. Figure 13 shows the bubble size distribution at x=0.2 m. It is seen that a large number of bubbles of 400 μm radius are now present in the field and would be seen, while all bubbles were originally smaller than 200 μm upstream when they entered the propeller flow field.

    Fig.14 Comparison of the time-averaged void fraction distribution,, along the axial direction between two different durations of the nuclei release

    6. Effect of the duration of nuclei release

    In order to investigate whether the duration of bubble release selected in the previous section, Δt= 3.7×10–3s, is long enough to be statistically meaningful, comparison with a computation with a 10fold longer release time, Δt=3.7×10–2s is considered here. Comparison ofthe void fractionsαisshown in Fig.14. It is seen that differences between the two release durations are minimal with less than 2.5% difference in the time-averaged void fraction downstream. This justifies the use of the shorter duration of 3.7×10–3s to conduct parametric studies.

    7. Parametric studies

    The effect of three important physical parameters: the cavitation number, the initial nuclei size distribution, and the initial dissolved gas concentration, on the bubble nuclei entrainment and dynamics are addressed in this section.

    Fig.15Comparison of the void fraction variations along x for three cavitation numbers, σ=1.75, 1.65 and 1.5

    7.1Cavitation number

    The time-averaged void fraction and downstream nucleisize distribution obtained at three different cav itation numbers, σ=1.5, 1.65 and 1.75, are compared in Fig.15 and Fig.16. As can be seen in the figures, in all three cases the voidfaction increases significantly inthe propeller liquid volume II (Fig.11(b)), th en plateaus. As the cavitation number decreases, the void fractions within the cavitation areas increase and the relative importance of the tip vortices cavitation volume diminishesrelative to the propeller blade cavitation. With smaller values ofσ, the bubbles retain a much larger size downstream of the propeller due to enhanced net gain of gas mass diffusing into the bubbles. The downstream void fractionαatx=0.2 m is 50 times the upstream value for σ=1.5, 2.35 times for σ=1.65 and 1.4 times for σ=1.75. The largest downstream bubble radii observed are 800 μm for σ=1.5, 600 μm for σ=1.65 and 400 μm for σ=1.75, as compared to 200 μm at release.

    Fig.17 Comparison of the void fraction variation along the propeller streamwise direction for two initial nuclei size distributions

    Fig.18 Bubble size distribution at x=0.2 m for two initial nuclei size distributions

    7.2Nuclei distribution

    To investigate the effect of the initial nuclei size distribution on the results a bubble size distribution ranging from 10 μm-200 μm is compared in Fig.17 to thatof a uniform distribution of 100 μm bubbles.Both distributions are selected to have the same initial void fraction. Figure 17 showsα(x) in both cases. Higher bubble void fractions are seen near the propeller blade region for the initial uniform nuclei size case. However, this does not result in any significant difference in the void fraction observed downstream. Figure 18 compares the resulting downstream nuclei size distributions at x=0.2 m. It is seen that the largest bubble radius is 400 μm for the non-uniform nuclei distribution, while it is 300 μm for the uniform case. This is related to the existence of larger bubbles in the initial non-uniform bubble nuclei distribution.

    7.3 Dissolved gas concentration

    The effect of the dissolved gas concentration on the bubble dynamics is illustrated in Fig.19 by c3onsidering disso3lved gas concentration of 0.66 mol/m and 1.32 mol/m. As seen in Fig.19, this results in further increase in the void fraction downstream. However, although the dissolved gas concentration was doubled, the void fraction increased by only 11%.

    Fig.19Comparison of the void fraction variation along the propeller streamwise direction for two initial gas concentrations

    8. Eff ect of artificial turbulence-like fluctuations

    The computations shown above used the flow field obtained by a RANS computation.As a result, the pressures and velocities used were quantities averaged by the RANS numerical procedure. On a real propeller, fluctuating turbulent quantities exists and modulate the RANS solution. In order to investigate the influence of turbulence on the results shown above, artificial turbulence-like fluctuations, (u',p'), were randomly added to the velocity and pressure field, (,), obtained from RANS

    To reproduce a realistic turbulent flow, the fluctuation velocity field, u', was composed of both spatial and temporal fluctuations. The pressure fluctuation field were then approximated by

    The following expressions give the fluctuating velocity field for theaxial component only (for simplicity of the presentation):

    To make the artificial fluctuations physics-based, their intensity Au(Ω)was related to the RANS determined vorticity distrib ution, Ω, andmade maximum at

    where A is the imposed amplitude of the oscillations and Ωmaxis the maximum vorticity in a given y-z plane.

    Fig.20 Contour plots of the imposed velocity fluctuations amplitude Au(Ω) in the y-z plane located at x= 0.05 m

    Figure 20 illustrates the contours of Au(Ω)in the y-z plane downstream of the propeller at x= 0.05 m.Two flow passages only of the five-bladed propeller are shown. As expected, the highest amplitude oscillations are in the tip vortices and along the wake regions of the propeller blades. In Eq.(25), f(x), f(y)andf(z) are spatial fluctuation functions, madeofΝx, Νyand Νzcomponents with wave lengths and phases, λxand φ(λx), λyand φ(λy), λz, and φ(λz) respectively. The phase shiftsφ are randomly distributed between [0,1]. Ax, Ayand Azare factors used to ensure that the fluctuation functions have unit values at the highest oscillations amplitude.

    Fig.21 Contour plots of the axial velocity fluctuations obtained with A=0.05 and composed of 400 different wave lengths ranging from 0.001 m to 1 m in boththey and z direction

    By multiplying f(x), f(y) and f(z )with Au(Ω), one can obtain the spatial velocity field fluctuations. The axial velocity fluctuation co ntours shown in Fig.21 were obtained with A=0.05 and 400 wave lengths ranging from 1 m to 0.001 m in both the y and z directions.

    Fig.22 Example fluctuations function with ωpeak=8 kHz, and composed of frequencies between 1 kHz and 15kHz

    Equation (25) also includes an imposed temporal fluctua tion, f(t), which is composed of Ntdiffe-re nt frequencies,ω, with phase shift φ(ω)randomly distributed between [0,1]. At is afactor to ensure that the fluctuation function has a unit value at the highest amplitude of oscillations. The exponential term included in the temporal function allows one to generate bias fluctuations with the highest energy centered in the specified peak frequency, ωpeak. Figure 22 illustratesanexample signal composed of frequencies between 1 kHz to 15 kHz with ωpeak=8 kHz shown in both the time and frequency domains.

    Fig.23 Encountered pressure (black) and bubble radius (grey) for a 50 μm nuclei which grew over a propeller blade suction surface, σ=1.75

    With this imposed artificial turbulence, weinvestigate the effects of the fluctuations amplitudeand frequency on the bubble dynamics and entrainment. Figure 23 shows in black the bubble encountered pressure (i.e., pressure felt by the bubble during its travel) and in grey the bubble radius versus time for a 50 μm bubble in absence of imposed fluctuations. This bubble grew to more than 0.004 m near a blade surface and then stabilized to ab out 200 μm downstream ofthe propeller. The cavitation number was σ= 1.75.

    Figure 24 shows similar curves as in Fig.23 but when artificial turbulent flow field like fluctuations were imposed. The results for three conditions are shown in the figure with different peak frequencies: ωpeak=2 kHz, 6 kHz and 10 kHz and A=0.4. It is clear from the figures that the bubble has the strongest oscillations and the largest downstream radius when the peak frequency is equal to 6 kHz. This can beexplained by considering the bubble natural frequency, ωR, given by

    where k is the polytrophic gas constant (=k1.4 was used in the current study). Not surprisingly, for the bubble shown in Fig.23, the natural oscillation frequencyis about 6 kHz since it has a radius of 200 μm and encounters a pressure of around 30 000 Pa. As expected, oscillations with a driving pressure field resonating with the bubble natural frequency enhance the gas diffusion. This resonating oscillations are expected to increase further the size of the resulting bubble as we have shown in Fig.1.

    Fig.24 Encountered pressure (black) and bubble radius (grey) for a 50 μm nuclei subjected to imposed artificial flow field fluctuations with three different peak frequencies, ωpeak=2 kHz, 6 kHz and 10 kHz, A=0.4

    Also, the effect of the amplitude of the fluctuations on the results is shown in Fig.25 for threedifferent values of the amplitude,A, and with ωpeak= 6kHz. As expected, the bubble oscillations increase when the amplitude of the fluctuations increases.

    Fig .25 Encountered pressure (black) and bubble radius (grey) for a 50 μm nuclei subjected to imposed artificial flow field fluctuations with three different fluctuation amplitudes: A=0.1, 0.2 and 0.4, ωpeak=6 kHz

    9.Conclusions

    Application of multi-bubble dynamics, accounting for gas diffusion effects to the study of bubble entrainment and dynamics in a propeller flow enables studyof the modification of natural waters nuclei distribution by a propeller. Comparison of the results between inclusion or neglect of gas diffusion shows that gas diffusion plays an important role on the resulting bubble sizes downstream from the propeller. Bubble sizes become larger than the original upstream sizes, downstream of the propeller due to a net influx of dissolved gas into the bubble. Bubble explosive stream void fraction. Different upstream nuclei distributions with the same void fraction have similar effegrowth and collapse, are also seen to be an essential“catalyst” to enable significant diffusion as bubbles which do not experience intense dynamics do not change size significantly.

    The model enabled simulations of bubble entrainment by a rotating propeller with a realistic ocean nuclei size distribution. Large visible bubbles were seen to collect in the blades vortices and wake regions. As the cavitation number decreases the downstream void fraction increases and the nuclei size distribution shifts towards larger sizes. Increase in dissolved gas concentration was also found to increase the downcts on the downstream void fraction but have a more significant effect on the downstream nuclei size distribution. Changes in the bubble size are indeed the most important effect observed in the study, with the propeller flow field significantly modifying the size distribution by both broadening the size range and increasing significantly the number of larger bubbles.

    Imposing artificial turbulence on the flow field results in enhanced bubble activity, which increases with the fluctuations amplitude. The magnitude of the bubble size oscillations reaches a maximum as the peak frequency of the imposed turbulent fluctuations matches the resonant natural bubble oscillation frequency.

    Acknowledgements

    This work was conducted at Dynaflow, Inc. (www.dynaflow-inc.com) and was supported by the Office of Naval Research (Grant No. N00014-05-C-0170) monitored by Dr. Patrick L. Purtell. We would also like to thank Dr. Kim Ki-Han at ONR for his support of this program.

    [1] RAICHLEN F. A., BRENNEN C. E. Measurement of air entrainment by bow waves[J]. Journal of Fluids Engineering, 2001, 123(1): 57-63.

    [2] DUNCAN J. H. Spilling breakers[J]. Annual Review of Fluid Mechanics, 2001, 33: 519-547.

    [3] MUSCARI R., MASCIO A. D. Numerical modeling of breaking waves generatedby a ship’s hull[J]. Journal of Marine Science and Technology, 2004, 9(4): 158- 170.

    [4] MORAGAF. J., CARRICA P. M. and DREW D. A. et al. A sub-grid air entrainment model for breaking bow waves and naval surface ships[J]. Computes and Fluids, 2008, 37(3): 281-298.

    [5] PEREIRA F. Transport of gas nuclei in a propeller[C]. Second International Symposium on Marine Propu- lsors SMP’11. Hamburg, Germany, 2011.

    [6] WU Xiong-jun, CHAHINE Georges L. Development of an acoustic instrument for bubble size distribution measurement[J]. Journal of Hydrodynamics, 2010, 22(5

    Suppl.): 330-336.

    [7] FRANKLIN R. E. A note on the radius distribution function for microbubbles of gas in water[C]. ASME Cavitation and Multiphase Flow Forum. Los Angeles, California, USA,1992, 77-85.

    [8] CHAHINE Georges L. Numerical simulation of bubble flow interactions[J]. Journal of Hydrodynamics, 2009, 21(3): 316-332.

    [9] HSIAO C.-T., CHAHINE G. L. Scaling of tip vortex cavitation inception noise with a statistic bubble dynamics model accounting for nuclei size distribution[J]. Journal of Fluid Engineering, 2005, 127(1): 55-64.

    [10] HSIAO C.-T., CHAHINE G. L. and LIU H.-L. Scaling effects on prediction of cavitation inception in a line vortex flow[J]. Journal of Fluid Engineering, 2003, 125(1): 53-60.

    [11] HSIAO C.-T., PAULEY L. L. Numerical computation of the tip vortex flow generated by a marine propeller[J]. Journal of Fluids Engineering, 1999, 121(3): 638-645. [12] RAYLEIGH L. On the pressure developed in a liquid during collapse of a spherical cavity[J]. Philosophical Magazine Series 6, 1917, 34(200): 94-98.

    [13] PLESSET M. S., DAVIES R. Cavitation in real liquids[M]. Elsevier Publishing Company, 1964, 1-17.

    [14] VOKURKA K. Comparison of Rayleigh’s, Herring’s, and Gilmore’s models of gas bubbles[J]. Acustica, 1986, 59(3): 214-219.

    [15] HABERMAN W. L., MORTON R. K. An experimental investigation of the drag and shape of air bubbles rising in various liquids[R]. Report 802, DTMB, 1953.

    [16] PLESSET M. S., ZWICK S. A. A non-steady heat diffusion problem with spherical symmetry[J]. Journal of Applied Physics, 1952, 23(1): 95-98.

    [17] CRUM L. A., MAO Y. Acoustically enhanced bubble growth at low frequencies and its implications for human diver and marine mammal safety[J]. Journal of the Acoustical Society of America, 1996, 99(5): 2898- 2907.

    [18] CHESNAKAS C., JESSUP S. Experimental characterization of propeller tip flow[C]. 22th Symposium on Naval Hydrodynamics. Washington D.C., 1998, 156-170.

    10.1016/S1001-6058(11)60308-9

    * Biography: HSIAO Chao-Tsung (1964-), Male, Ph. D., Principal Research Scientist

    天天添夜夜摸| 久久热在线av| 少妇熟女aⅴ在线视频| 国产亚洲欧美98| 亚洲欧美激情综合另类| а√天堂www在线а√下载| 免费高清视频大片| 亚洲男人天堂网一区| 人人妻人人澡欧美一区二区 | 12—13女人毛片做爰片一| 亚洲色图av天堂| 午夜日韩欧美国产| 欧美 亚洲 国产 日韩一| 大又大粗又爽又黄少妇毛片口| 亚洲av成人精品一区久久| 欧美色欧美亚洲另类二区| 亚洲av不卡在线观看| 联通29元200g的流量卡| 亚洲色图av天堂| 91午夜精品亚洲一区二区三区 | 99riav亚洲国产免费| 久99久视频精品免费| 蜜桃亚洲精品一区二区三区| 久久久久久久久久成人| 久久精品人妻少妇| 色综合亚洲欧美另类图片| 变态另类丝袜制服| 一夜夜www| 最近在线观看免费完整版| 如何舔出高潮| 中文字幕高清在线视频| 欧美激情国产日韩精品一区| 在线天堂最新版资源| 精品人妻视频免费看| 亚洲精品亚洲一区二区| 超碰av人人做人人爽久久| 欧美激情国产日韩精品一区| 国产极品精品免费视频能看的| 精品久久久久久,| 久久精品国产清高在天天线| 亚洲av一区综合| 精品久久久久久久末码| 特大巨黑吊av在线直播| 国产午夜精品论理片| 免费观看在线日韩| 99热网站在线观看| 99九九线精品视频在线观看视频| 久久九九热精品免费| 欧美另类亚洲清纯唯美| 国产精品人妻久久久久久| 国产视频一区二区在线看| 国产三级在线视频| 成人特级av手机在线观看| 午夜福利在线观看吧| 麻豆精品久久久久久蜜桃| 九色国产91popny在线| 搞女人的毛片| 免费看a级黄色片| 欧美黑人欧美精品刺激| 午夜老司机福利剧场| 午夜爱爱视频在线播放| 国产探花在线观看一区二区| 波野结衣二区三区在线| 亚洲国产精品sss在线观看| АⅤ资源中文在线天堂| 亚洲欧美清纯卡通| 国产av一区在线观看免费| 中文资源天堂在线| 午夜精品一区二区三区免费看| 精品国内亚洲2022精品成人| 悠悠久久av| 香蕉av资源在线| 日本爱情动作片www.在线观看 | 九九热线精品视视频播放| 国产亚洲91精品色在线| 午夜免费激情av| 蜜桃亚洲精品一区二区三区| 久久精品国产亚洲av涩爱 | 一级黄片播放器| 级片在线观看| 国产毛片a区久久久久| 亚洲真实伦在线观看| 国产亚洲91精品色在线| 一进一出抽搐动态| 联通29元200g的流量卡| 大型黄色视频在线免费观看| 男女边吃奶边做爰视频| 精品人妻熟女av久视频| 国内精品久久久久精免费| h日本视频在线播放| 国产三级在线视频| 超碰av人人做人人爽久久| 国产精品一区www在线观看 | 91av网一区二区| 国产精品一区二区三区四区久久| 午夜福利在线观看吧| 欧美色视频一区免费| 国产单亲对白刺激| 国产成人一区二区在线| 波多野结衣高清无吗| 国产精品一区二区三区四区久久| 99riav亚洲国产免费| 欧美另类亚洲清纯唯美| 免费无遮挡裸体视频| 精品久久久噜噜| 人人妻人人澡欧美一区二区| av在线老鸭窝| 亚洲av电影不卡..在线观看| 免费观看精品视频网站| 亚洲欧美日韩高清在线视频| 中文字幕免费在线视频6| 最近最新中文字幕大全电影3| 欧美极品一区二区三区四区| av在线老鸭窝| 国产伦在线观看视频一区| 一个人看视频在线观看www免费| 人妻制服诱惑在线中文字幕| 黄色配什么色好看| 亚洲精品亚洲一区二区| 99久久无色码亚洲精品果冻| 成人鲁丝片一二三区免费| 成人高潮视频无遮挡免费网站| 日韩中字成人| 毛片女人毛片| 国产成人一区二区在线| 免费av毛片视频| 日本一二三区视频观看| 色综合婷婷激情| 国产三级在线视频| 午夜免费激情av| 99热精品在线国产| 成人三级黄色视频| 一本久久中文字幕| 91狼人影院| av国产免费在线观看| 精品久久久久久久末码| 精品久久久久久久人妻蜜臀av| 午夜福利在线在线| 不卡视频在线观看欧美| 在线观看一区二区三区| 欧美zozozo另类| 精品午夜福利视频在线观看一区| 成年女人看的毛片在线观看| 日本精品一区二区三区蜜桃| 国产69精品久久久久777片| 黄色一级大片看看| 哪里可以看免费的av片| 欧美成人免费av一区二区三区| 波野结衣二区三区在线| a级毛片a级免费在线| 免费看av在线观看网站| 国产精品久久久久久久电影| 露出奶头的视频| 日本撒尿小便嘘嘘汇集6| 国产精品99久久久久久久久| 日韩欧美三级三区| 51国产日韩欧美| 中国美白少妇内射xxxbb| 亚洲国产精品sss在线观看| 国产精品久久久久久av不卡| 小蜜桃在线观看免费完整版高清| 热99在线观看视频| 国产精华一区二区三区| 亚洲无线观看免费| 亚洲一级一片aⅴ在线观看| 精品一区二区三区视频在线| 亚洲,欧美,日韩| 日韩欧美一区二区三区在线观看| 欧美色视频一区免费| 亚洲国产日韩欧美精品在线观看| 久久久成人免费电影| 久久亚洲真实| 真实男女啪啪啪动态图| 亚洲精品日韩av片在线观看| 国模一区二区三区四区视频| 91在线观看av| 国产精品一及| 成人鲁丝片一二三区免费| 国产高清不卡午夜福利| 亚洲性夜色夜夜综合| 又爽又黄无遮挡网站| 亚洲色图av天堂| АⅤ资源中文在线天堂| 九九爱精品视频在线观看| 1000部很黄的大片| 成人特级av手机在线观看| 人妻久久中文字幕网| 熟女电影av网| 最近中文字幕高清免费大全6 | 日本熟妇午夜| 国产男人的电影天堂91| 国产v大片淫在线免费观看| 九九在线视频观看精品| 免费观看人在逋| 色尼玛亚洲综合影院| 日韩欧美 国产精品| 午夜亚洲福利在线播放| 可以在线观看毛片的网站| 搡老岳熟女国产| 男插女下体视频免费在线播放| 真人做人爱边吃奶动态| 他把我摸到了高潮在线观看| 一级av片app| 久久香蕉精品热| 综合色av麻豆| 精品久久久久久久久久免费视频| 最新中文字幕久久久久| xxxwww97欧美| 永久网站在线| 国产精品女同一区二区软件 | 男女做爰动态图高潮gif福利片| 国模一区二区三区四区视频| 免费一级毛片在线播放高清视频| 免费av毛片视频| 亚洲人成伊人成综合网2020| 亚洲成人免费电影在线观看| 亚洲av第一区精品v没综合| 老女人水多毛片| 精品人妻偷拍中文字幕| 春色校园在线视频观看| 老熟妇乱子伦视频在线观看| 最近视频中文字幕2019在线8| 久久久久久久久久久丰满 | 美女cb高潮喷水在线观看| 中国美白少妇内射xxxbb| 国产私拍福利视频在线观看| 91麻豆av在线| 亚洲精品在线观看二区| 给我免费播放毛片高清在线观看| 婷婷亚洲欧美| 日韩欧美免费精品| 国产在线男女| 亚洲av免费在线观看| 国产日本99.免费观看| 女同久久另类99精品国产91| 国产伦人伦偷精品视频| 成人精品一区二区免费| 丰满人妻一区二区三区视频av| 18禁黄网站禁片午夜丰满| 亚洲av中文字字幕乱码综合| 波多野结衣高清无吗| 一区二区三区四区激情视频 | 久久精品影院6| 欧美性猛交黑人性爽| 69人妻影院| 色播亚洲综合网| 午夜老司机福利剧场| 色在线成人网| 精品人妻偷拍中文字幕| 人妻夜夜爽99麻豆av| 深夜a级毛片| 99久久精品热视频| 久久久久久国产a免费观看| 毛片一级片免费看久久久久 | 香蕉av资源在线| 深夜精品福利| 99国产精品一区二区蜜桃av| 成熟少妇高潮喷水视频| 国产精品亚洲美女久久久| 男人舔奶头视频| 麻豆国产97在线/欧美| 99九九线精品视频在线观看视频| 少妇裸体淫交视频免费看高清| 久久久精品欧美日韩精品| 男女视频在线观看网站免费| 99久国产av精品| 我的老师免费观看完整版| 亚洲av熟女| 中文字幕高清在线视频| 男人舔女人下体高潮全视频| 欧美丝袜亚洲另类 | 久99久视频精品免费| 波野结衣二区三区在线| 岛国在线免费视频观看| 久久久久久久久久成人| 亚洲最大成人av| 国产高清三级在线| 热99在线观看视频| 精品日产1卡2卡| 亚洲精品亚洲一区二区| 99精品久久久久人妻精品| 又紧又爽又黄一区二区| 91久久精品国产一区二区成人| 中文字幕高清在线视频| 亚洲av美国av| 高清在线国产一区| 午夜a级毛片| 99热这里只有是精品在线观看| 国产精品久久久久久久电影| 欧美xxxx黑人xx丫x性爽| 亚洲无线观看免费| a级毛片a级免费在线| 国产成人a区在线观看| 久久精品国产自在天天线| 免费在线观看成人毛片| 亚洲专区中文字幕在线| 好男人在线观看高清免费视频| 伦理电影大哥的女人| 乱人视频在线观看| 成人精品一区二区免费| 成年女人永久免费观看视频| 啦啦啦观看免费观看视频高清| 国产伦精品一区二区三区视频9| 最新在线观看一区二区三区| 波多野结衣高清无吗| 男人和女人高潮做爰伦理| 国产黄片美女视频| 一区二区三区四区激情视频 | 精品久久久久久久久久免费视频| 香蕉av资源在线| 欧美日韩亚洲国产一区二区在线观看| 亚洲自拍偷在线| 人妻少妇偷人精品九色| 我的老师免费观看完整版| 中亚洲国语对白在线视频| 国产亚洲精品av在线| 床上黄色一级片| 日本在线视频免费播放| 日本 av在线| 男人狂女人下面高潮的视频| 成人一区二区视频在线观看| 天天一区二区日本电影三级| 美女大奶头视频| АⅤ资源中文在线天堂| 最近中文字幕高清免费大全6 | 两人在一起打扑克的视频| 午夜久久久久精精品| 日韩在线高清观看一区二区三区 | 日本与韩国留学比较| 欧美一区二区国产精品久久精品| 日韩 亚洲 欧美在线| 国产大屁股一区二区在线视频| 国产精品乱码一区二三区的特点| 成年女人永久免费观看视频| av专区在线播放| 国产成年人精品一区二区| 日本在线视频免费播放| 国产免费av片在线观看野外av| 久久草成人影院| 欧美日本亚洲视频在线播放| 国产69精品久久久久777片| 男人舔奶头视频| 国产精品美女特级片免费视频播放器| 12—13女人毛片做爰片一| 毛片女人毛片| 国产免费男女视频| 3wmmmm亚洲av在线观看| 久久久久久久久久黄片| 亚州av有码| 国产视频内射| 熟女电影av网| 偷拍熟女少妇极品色| 在线播放国产精品三级| 日本免费一区二区三区高清不卡| 观看美女的网站| 亚洲第一电影网av| 国产精品一区二区性色av| 一级毛片久久久久久久久女| 色尼玛亚洲综合影院| 两个人的视频大全免费| 精品免费久久久久久久清纯| 亚洲在线自拍视频| 欧美又色又爽又黄视频| 日本a在线网址| 亚洲图色成人| 国产精品久久久久久久电影| 国产一区二区亚洲精品在线观看| 亚洲成人久久爱视频| 熟妇人妻久久中文字幕3abv| 嫩草影院新地址| 久久精品国产亚洲av涩爱 | 亚洲欧美日韩无卡精品| 国产v大片淫在线免费观看| 少妇人妻一区二区三区视频| 婷婷丁香在线五月| 草草在线视频免费看| 亚洲人成网站在线播| 亚洲精品成人久久久久久| 免费看日本二区| 亚洲美女黄片视频| 婷婷丁香在线五月| 亚洲美女视频黄频| 精品人妻一区二区三区麻豆 | 精品久久国产蜜桃| 国产亚洲精品综合一区在线观看| 国内毛片毛片毛片毛片毛片| 欧美日韩瑟瑟在线播放| 一个人观看的视频www高清免费观看| 少妇熟女aⅴ在线视频| 亚洲七黄色美女视频| 黄色配什么色好看| 国模一区二区三区四区视频| 午夜福利在线观看免费完整高清在 | 国产高清视频在线播放一区| av女优亚洲男人天堂| 老熟妇仑乱视频hdxx| 国产精品一区www在线观看 | 99久久成人亚洲精品观看| 国产伦精品一区二区三区视频9| 精品久久久久久久久久久久久| 色吧在线观看| 99精品在免费线老司机午夜| 国产成人福利小说| 老司机深夜福利视频在线观看| a级毛片免费高清观看在线播放| 午夜日韩欧美国产| 天美传媒精品一区二区| 国产乱人伦免费视频| 国模一区二区三区四区视频| 免费在线观看影片大全网站| aaaaa片日本免费| 村上凉子中文字幕在线| 性欧美人与动物交配| 91av网一区二区| 免费黄网站久久成人精品| 色综合亚洲欧美另类图片| 国产精品一区二区免费欧美| 欧美一区二区国产精品久久精品| 欧美人与善性xxx| 欧美zozozo另类| 亚洲av五月六月丁香网| 校园人妻丝袜中文字幕| 国产精品1区2区在线观看.| 久久久久九九精品影院| 久久热精品热| 久久草成人影院| 人人妻人人澡欧美一区二区| 亚洲欧美精品综合久久99| avwww免费| 亚洲天堂国产精品一区在线| 色哟哟·www| av在线亚洲专区| 亚洲av熟女| 一个人看的www免费观看视频| 99国产精品一区二区蜜桃av| 亚洲国产高清在线一区二区三| 在线免费十八禁| 国产精品久久久久久亚洲av鲁大| 91麻豆av在线| 丰满人妻一区二区三区视频av| 日日夜夜操网爽| 国产精品自产拍在线观看55亚洲| 久久久久国产精品人妻aⅴ院| 亚洲精品影视一区二区三区av| 国产精品亚洲一级av第二区| 欧美成人一区二区免费高清观看| 免费在线观看影片大全网站| 啦啦啦啦在线视频资源| 婷婷六月久久综合丁香| 国产在线精品亚洲第一网站| 色综合站精品国产| 亚洲性久久影院| 长腿黑丝高跟| 亚洲国产日韩欧美精品在线观看| 亚洲av熟女| 免费av不卡在线播放| 国产精品嫩草影院av在线观看 | 男人和女人高潮做爰伦理| 成人国产一区最新在线观看| 一级av片app| 久久久成人免费电影| 国内精品久久久久精免费| 嫩草影院精品99| 久久久久国产精品人妻aⅴ院| 99riav亚洲国产免费| 欧美一区二区精品小视频在线| 欧美xxxx性猛交bbbb| 亚洲成人中文字幕在线播放| 国产高清三级在线| 我要看日韩黄色一级片| 国产精品1区2区在线观看.| 在现免费观看毛片| 国产三级在线视频| av黄色大香蕉| 欧美一区二区国产精品久久精品| 搡老熟女国产l中国老女人| 欧美成人性av电影在线观看| 精品一区二区三区视频在线观看免费| 亚洲人成伊人成综合网2020| 成人三级黄色视频| 日本一本二区三区精品| 99九九线精品视频在线观看视频| av天堂中文字幕网| 欧美xxxx性猛交bbbb| 久久精品国产自在天天线| 日本一本二区三区精品| 俄罗斯特黄特色一大片| 日本成人三级电影网站| 99热6这里只有精品| 一区福利在线观看| 免费在线观看日本一区| 观看美女的网站| 国产亚洲精品久久久久久毛片| 精品人妻熟女av久视频| 男人舔女人下体高潮全视频| 精品国内亚洲2022精品成人| 最近中文字幕高清免费大全6 | 夜夜爽天天搞| ponron亚洲| 老司机福利观看| 亚州av有码| 在线观看舔阴道视频| 国产精品久久久久久av不卡| 国内少妇人妻偷人精品xxx网站| 国产精品久久久久久久久免| 精品午夜福利在线看| 国产伦一二天堂av在线观看| 久久久久九九精品影院| 免费看a级黄色片| 国产av一区在线观看免费| 床上黄色一级片| 成人高潮视频无遮挡免费网站| 久久午夜亚洲精品久久| 91av网一区二区| 亚洲人成网站高清观看| 一个人免费在线观看电影| 亚洲人成网站在线播放欧美日韩| 精品无人区乱码1区二区| 男人舔女人下体高潮全视频| 美女被艹到高潮喷水动态| 欧美zozozo另类| 99久久久亚洲精品蜜臀av| 一级av片app| 国产成年人精品一区二区| 成人永久免费在线观看视频| 桃红色精品国产亚洲av| 国产精品福利在线免费观看| 国产aⅴ精品一区二区三区波| 乱人视频在线观看| 在线天堂最新版资源| 桃色一区二区三区在线观看| 亚洲一级一片aⅴ在线观看| 亚洲无线观看免费| 春色校园在线视频观看| 国产 一区精品| 欧美一级a爱片免费观看看| 一个人免费在线观看电影| 999久久久精品免费观看国产| 男人狂女人下面高潮的视频| 十八禁网站免费在线| 成人高潮视频无遮挡免费网站| 精品一区二区三区视频在线| 在线观看一区二区三区| 国产在线男女| 一边摸一边抽搐一进一小说| 2021天堂中文幕一二区在线观| 老女人水多毛片| 日韩欧美免费精品| 身体一侧抽搐| 极品教师在线视频| 变态另类丝袜制服| 亚洲av二区三区四区| 国产美女午夜福利| 色综合亚洲欧美另类图片| 亚洲国产日韩欧美精品在线观看| 春色校园在线视频观看| 搡老岳熟女国产| 窝窝影院91人妻| 亚洲精品影视一区二区三区av| 日日啪夜夜撸| 国产aⅴ精品一区二区三区波| 亚洲av一区综合| 国产精品一区二区性色av| 亚洲av电影不卡..在线观看| 精品久久久久久久久av| 日韩欧美精品免费久久| 免费看a级黄色片| av在线老鸭窝| 91在线精品国自产拍蜜月| 一区二区三区高清视频在线| 蜜桃久久精品国产亚洲av| 综合色av麻豆| 一区二区三区激情视频| 久久草成人影院| 色哟哟·www| 变态另类成人亚洲欧美熟女| 日本黄色片子视频| 天堂动漫精品| 又黄又爽又免费观看的视频| 午夜激情欧美在线| 在线免费观看的www视频| 熟妇人妻久久中文字幕3abv| 99精品久久久久人妻精品| 久久婷婷人人爽人人干人人爱| 在线国产一区二区在线| 亚洲精品久久国产高清桃花| 久久人人精品亚洲av| 联通29元200g的流量卡| 色精品久久人妻99蜜桃| av在线观看视频网站免费| 精品久久久噜噜| 女人十人毛片免费观看3o分钟| 嫩草影视91久久| 日本一本二区三区精品| 丝袜美腿在线中文| 日本五十路高清| 亚洲人成网站高清观看| 亚洲国产精品成人综合色| 日本五十路高清| 国产色婷婷99| 欧美黑人巨大hd| 人人妻人人澡欧美一区二区| 免费在线观看日本一区| 欧美xxxx黑人xx丫x性爽| 亚洲欧美日韩无卡精品| 免费一级毛片在线播放高清视频| 国产成人影院久久av| 欧美日韩精品成人综合77777| 色综合婷婷激情| 亚洲av免费高清在线观看| 亚洲美女视频黄频| 欧美人与善性xxx| 日韩 亚洲 欧美在线| 久久久久久久久久黄片| 成年版毛片免费区|