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

    Large eddy simulation of turbulent attached cavitating flow with special emphasis on large scale structures of the hydrofoil wake and turbulence-cavitation interactions*

    2017-03-09 09:09:20BinJi季斌YunLong龍云XinpingLong龍新平ZhongdongQian錢忠東JiajianZhou周加建
    水動力學研究與進展 B輯 2017年1期
    關鍵詞:龍云新平

    Bin Ji (季斌), Yun Long (龍云), Xin-ping Long (龍新平), Zhong-dong Qian (錢忠東),Jia-jian Zhou (周加建)

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

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

    3.Science and Technology on Water Jet Propulsion Laboratory, Shanghai 200011, China

    (Received June 9, 2016, Revised December 6, 2016)

    Large eddy simulation of turbulent attached cavitating flow with special emphasis on large scale structures of the hydrofoil wake and turbulence-cavitation interactions*

    Bin Ji (季斌)1,2,3, Yun Long (龍云)1,2, Xin-ping Long (龍新平)1,2, Zhong-dong Qian (錢忠東)1,Jia-jian Zhou (周加建)3

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

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

    3.Science and Technology on Water Jet Propulsion Laboratory, Shanghai 200011, China

    (Received June 9, 2016, Revised December 6, 2016)

    In this paper, the turbulent attached cavitating flow around a Clark-Y hydrofoil is investigated by the large eddy simulation (LES) method coupled with a homogeneous cavitation model. The predicted lift coefficient and the cavity volume show a distinctly quasi-periodic process with cavitation shedding and the results agree fairly well with the available experimental data. The present simulation accurately captures the main features of the unsteady cavitation transient behavior including the attached cavity growth, the sheet/cloud cavitation transition and the cloud cavitation collapse. The vortex shedding structure from a hydrofoil cavitating wake is identified by theQ-criterion, which implies that the large scale structures might slide and roll down along the suction side of the hydrofoil while being further developed at the downstream. Further analysis demonstrates that the turbulence level of the flow is clearly related to the cavitation and the turbulence velocity fluctuation is much influenced by the cavity shedding.

    Cavitation, large eddy simulation (LES), vortex structure, turbulence-cavitation interactions

    Introduction

    The cavitation involves a very complex twophase flow. It can have negative impacts on the hydraulic machinery such as through vibrations, erosions, noises and performance break-down[1]. Therefore, it is necessary to investigate its physical mechanism for controlling cavitation in engineering applications.

    The fundamental physics of the cavitation phenomenon was widely studied. Much attention was paid on the evolution of both sheet and cloud cavitations on the hydrofoil. The typical quasi-periodic phenomenon was a focus, including the sheet cavitation inception,growth and shedding[2-4]. Li et al.[5]reproduced the quasi-periodic phenomenon by using the large eddy simulation (LES) method and the VOF model. The VOF technique is widely used for the simulation of the free surface flow[6]. The cavity interface can be predicted accurately[7]. Gnanaskandan and Mahesh[8]studied the sheet to cloud cavitation transition over a wedge. It is found that the frequency of the shedding process obtained by the LES is in good agreement with that of the experiments. These studies show that the recent methods can accurately capture the detailed features of the cavitation. Studies about the sheet to cloud cavitation show that the re-entrant jet is a very important source factor for the cavity shedding. Le et al.[9]studied partial cavities. They observed that the cavitation periodic shedding was closely associated with the re-entrant jet. Pham et al.[10]confirmed the conclusions of Le et al.[9]by analyzing the frequencies of the re-entrant jet surges and the cloud shedding. Recent studies confirm that the re-entrant jet is the main cause for the development and the shedding of the attachedcavity[2]. The studies of the cylindrical objects[11], the head-form body[12]and others[13-15]demonstrate that the re-entrant jets are responsible for the cavity shedding.

    The unsteady cavity shedding often involves a complex cavity vortex structure. Therefore, the cavitation-vortex interactions are important research issues. Gopalan and Katz[16]demonstrated by means of the PIV and the high speed photography that the collapse of the vapor cavity in the downstream region was mainly responsible for the generation of the vortex. The LES method is a useful tool to reproduce the unsteady cavitating flow, with great advantages in accurately capturing the complex vortex structures. It was firstly proposed by Smagorinsky[17], and now it is widely applied for the simulations[18,19]. Roohi et al.[7]combined the 2-D LES with the VOF and the cavitation model and showed that it was useful to simulate the shape of the cloud cavity and its dynamics. Ji et al.[13]and Luo et al.[20]studied the 3-D cavitation and vortex structures based on the LES around a twisted hydrofoil. Ji et al.[21]indicated that the mass transfer along the cavity surface induces an increase in the vortex dilatation. They found that the baroclinic torque term is important along the liquid-vapor interface but has a negligible effect inside the attached cavity region. Dreyer et al.[22]studied the tip leakage vortex on a NACA0009 hydrofoil using the stereo-PIV. The results revealed that it is prone to generate the cavitation when the vortex length is the maximum. All these studies have enhanced the understanding of the cavitation-vortex interactions, but the large scale vortex structures remain a difficult issue. More recently, Peng et al.[23]performed a series of experimental observations around hydrofoils in the cavitation tunnel. They observed that the U-type fl ow structures were common in cloud cavities. The side reentrant jets were identified as the main cause of the sheet cavity shedding and the formation of the U-type vortex structure was controlled in the cloud cavitation, which usedQ-criterion[24]to visualize the vortex structure evolution in numerical results. This helps us gain a deeper understanding of the cavitation shedding dynamics and cavitation-vortex interactions.

    Besides these fundamental studies of the dynamics of the sheet and cloud cavitation and the cavitation-vortex interactions, various kinds of cavitations were also studied, among others, the cavitating fl ows around a semi-circular leading-edge fl at plate[25], the axisymmetric projectile[26], the 3-D hemispherical head-form body[12], and behind a 3-D disk[27]by numerical and experimental methods. Wang et al.[28]studied the effect of free surface on the unsteady behavior of the cloud cavity. They found that the cavity evolution was remarkably different when the axisymmetric projectile was near the free surface. Park and Rhee[29]performed a comparative study of incompressible and compressible flow solvers for the cavitating flow. They concluded that the compressibility effects can enhance the reproduction of the cavitation. Various features of the cavitation at different cavitators are thus revealed, and the understanding and the ability to control and utilize the cavitation are enhanced.

    The pressure and velocity fluctuations were also attracted some attention. Chen et al.[30]found that the pressure fl uctuations are closely related to the transient cavitation behaviors, and they are excited by the acceleration due to the changes in the cavity volume. Ji et al.[31]illustrated the relationship between the cavitation evolution and the pressure fluctuations by analyzing the cavity volume fluctuation. They came to a similar conclusion that the cavity volumetric acceleration is the main source of the cavitation excited pressure. Huang et al.[32]introduced a filter-based corrected model. They found that the transient cavitation behaviors significantly enhanced the turbulence velocity fluctuations to induce a thicker turbulence boundary layer. More researches show that the streamwise velocity fluctuations generally are dominant inside the cavity[33], and the streamwise and spanwise fluctuations are equally important at the cavity closure and in the wake[8]. The turbulence-cavitation interactions continue to attract many research interests, but their mechanisms remain not well understood[34,35].

    Although a great progress has been made on the cavitating flow and the cavity vortex structures by both numerical and experimental methods, the large scale structures in the hydrofoil wake are still a worthwhile topic. Furthermore, the turbulence-cavitation interactions remain elusive due to the unsteadiness and the complex flow structures. This paper studies the cavitating flow around a Clark-Y hydrofoil by investigating the turbulent attached cavitating flow with a special emphasis on the large scale structures in the hydrofoil wake and the turbulence-cavitation interactions with the LES approach.

    1. Cavitation model and numerical methods

    The LES method is employed in this paper to solve the unsteady Navier-Stokes equations coupled with a mass transfer cavitation model. The main features of the solver are as follows.

    1.1Physical cavitation model

    Schnerr and Sauer[36]developed the cavitation model adopted in this paper to describe the mass transfer through the interface of the vapor and the liquid. The cavitation model is described by the follow equation

    whereαvis the vapor volume fraction. The source termsm˙+andm˙?in Eq.(1) are deduced from the bubble dynamics equation and they denote the evaporation and condensation rates, which are defined as follows

    whereRbis the bubble radius which is computed by

    whereNbis the bubble number density, which is set as 1013, as recommended by Schnerr and Sauer[36].

    1.2Governing equations and the LES approach

    In the basic governing equations in the mixture model adopted in this study, the same velocity and pressure for the vapor and the liquid are assumed, and the mass and momentum conservation equations are as follows

    whereuiis the velocity component in theidirection andpis the mixture pressure which is assumed taking the same value for the vapor and the liquid, andμis the laminar viscosity andρis the mixture density, and they are defined as follows

    The LES equations are obtained by applying a Favre-filtering operation to Eqs.(5) and (6). The equations are as follow:

    where the over-bars denote the filtered quantities andτijmeans the subgrid-scale stresses which are the unknown terms and defined as

    whereis the rate of the strain tensor in the resolved scale and is defined as

    whereμtis the subgrid-scale turbulence viscosity closed by applying the wall-adopting local eddy-viscosity (WALE) model[37]in this paper. The WALE model can return a zero turbulence viscosity for the laminar shear flows to allow the correct treatment of the laminar zones in the domain. In the WALE model,μtis modeled as

    wherekis the von Karman’s constant,dis the distance to the closest wall,Vis the volume of the computational cell andCSis the default WALE constant with the value of 0.5[37].

    1.3Simulation setup

    A Clark-Y hydrofoil with a chord length ofC=0.07 m is studied in this paper. The hydrofoil is placed horizontally in the computational domain with a fixed angle of attack ofα=8oas shown in Fig.1. Note that it is difficult to use the experimental spanwise size in view of the huge computational cost. Therefore, the spanwise length is chosen as0.3Cin Fig.1(a). The adopted length and height in the computational section is10Cand2.7C, respectively, as the height is usually chosen at least as twice of the thickness of the hydrofoil for capturing the streamwise vortex accurately[38]. The boundary inlet is arranged an imposed velocity and the outlet is with a fixed static pressure and a no-slip boundary on the surface of the hydrofoil. The free slip wall is imposed on the upper and lower walls. The side wall boundary is arranged as the periodicity interface. An O-H type grid as shown in Fig.2 is generated with a sufficient refinement near the foil surface for the LES method. The calculation mesh is a 3-D mesh and the 2-D mesh in Fig.2 is just its 2-D side view. The grid sizes near the hydrofoil of Fig.2 in the stream-wise, wall-normal and span-wise directions are Δx=2× 10?4m,Δy=2× 10?5mand Δz=2.1× 10?4m, respectively.

    Fig.1 Computational domain and boundary conditions

    The time-dependent governing equations are discretized in both space and time domains. The second order central difference scheme is adopted for the diffusion term, with a high order resolution scheme[39]for the convective term in the governing equations. The second-order implicit formulation is adopted for the transient term. The simulation setup is similar to that described in our previous paper[21]. The simulation is carried out initially under the steady non-cavitation condition and then, the unsteady cavitating flow is modeled with the steady non-cavitation results as the initial conditions. The time step isTref/200and theYplus is about 1. The residual is 10?3for the continuity equation, the momentum equations and the mass transport equation. The maximum internal iteration steps for the unsteady cavitating flow in each time step is 50. The calculation is conducted in a computer of 64G memory and 8-core parallel, and with about 10 000 time steps, and it takes about 1 week. The convergence process is smooth and the results are as follows.

    Fig.2 Structural mesh around the hydrofoil surface

    2. Results and discussions

    2.1Periodic behaviors of the hydrodynamic performance and the cavity volume fluctuationsThe lift coefficient is defined as

    Figure 3 shows the predicted lift coefficient against the time. The lift coefficient clearly shows periodic behaviors against the time. It varies dramatically as the pressure distribution changes sharply, which is reasonable in view of the sheet cavity shedding and the cloud cavity, as illustrated in the followed discussion. Figure 4 presents the numerical total cavity volume variation over three typical cycles. The cavity volume obtained by simulation is defined as

    whereNis the total number of the control volumes in the computational domain,αiis the vapor volumefraction in each control volume andViis the volume of each cell.Scavin Fig.4 is the hydrofoil profile area. The reference time isTref=40 msbased on Huang et al.[32]. The experimental cavity volume is obtained from Huang et al.[32]. It is indicated that the fluctuation of the cavity volume occurs during the cavity shedding. Wang et al.[40]confirmed that the large fluctuation of the total cavity volume occurs when the cavitation exists. It is observed that the cavity volume variation is a distinctly quasi-periodic process with the cavity shedding. The numerical largest cavity volume is smaller than that in the experiment when the LES method is used, but the predicted total cavity volume is in good agreement with the experimental data in general. In Huang et al.[32], the numerical largest cavity volume is a little bit higher than the experimental one obtained with the fi lter-based density corrected model (FBDCM), but the variation of the cavity volume is indeed captured. Overall, the lift coefficient and the cavity volume show periodic behaviors. The large fluctuation of the lift coefficient and the cavity volume occurs during the cavitation shedding. Next we will have detailed discussions about the time evolution of the cavity, the unsteady cavity and the vortex structures and their interactions.

    Fig.3 Predicted lift coefficients against the time.σ=0.8,α=8o,U=10 m/s∞

    Fig.4 (Color online) Cavity volume variation during cavitation shedding (Exp.[32]vs. Cal.).σ=0.8,α=8o,U∞= 10 m/s

    2.2Time evolution of the unsteady cavitation and discussions of large scale structures in the hydrofoil wake

    To show the calculated time evolution of the unsteady cavitation, eight numerical snapshots are shown in Fig.5. The experimental results are from the doctoral thesis of Huang[41]. Figure 5 shows the comparison between the experimental and predicted time evolutions of the cavitation shedding in one typical cycle. The numerical snapshots show that the experimental and predicted cavitation shedding patterns within one typical cycle are in good agreement, except for a few differences. The numerical results obtained with the LES method have captured the transient cavitation features visualized in the experiment, including the inception, the growth and the shedding of the sheet cavity. The detailed analyses are as follows. Att=1/8T, a thin attached cavity appears near the foil leading edge (L.E.) in Fig.5(b1), and the cloud cavity near the trailing edge (T.E.) resulted by the breakup of the previous sheet cavity is still seen apparently. Then, the numerical cavitation pattern shows that the length of the attached cavity increases as shown in Figs.5(b2) and 5(b3). The length of the attached cavity grows to its maximum att=4/8Tand the cavity surface near the T.E. becomes a bubble from now on. In Fig.5(b5), a re-entrant jet appears apparently while the adverse pressure gradient is strong enough to conquer the momentum near the foil surface. The re-entrant jet will be visualized by the distribution of the axial velocity along the foil later. Then, the re-entrant jet develops towards the L.E. as shown in Figs.5(b6) and 5(b7). With the development of the re-entrant jet, the attached cavity surface becomes an unstable turbulent vapor cloud and begins to shed from the foil surface. The unstable turbulent vapor cloud finally leads to a highly unstable cloudy cavitation in Figs.5(b5)-5(b8), resulting in a highly turbulent flow which will be discussed later on in Section 2.3. Meanwhile, the cavity volume decreases to the minimum as shown in Fig.4. Overall, the numerical results predicted with the LES method are in a reasonable agreement with the experimental pictures to show that the LES method can accurately capture the features of the unsteady cavitation transient behavior.

    The vortex structures in one typical cycle visualized based on theQ-criterion andα=10%isosurface are shown in Fig.6. It is useful to apply theQ-criterion for presenting the vorticity due to the vortex rather than the shear. TheQ-criterion is introduced by Hunt et al.[24], and it is defined as

    Fig.5 (Color online) Time evolution of cavitation shedding in one typical cycle.σ=0.8,α=8o,U∞=10 m/s

    Fig.6 (Color online) Vortex structures in one typical cycle forσ=0.8,α=8o

    where?is the vorticity tensor andSis the rate of the strain tensor. The value of theQ-criterion is set as 20 000 in Fig.6 to ensure that the vorticity is mainly resulted by the vortex. As is shown in Fig.6, the vortex shape based on theQ-criterion is much more complicated than the vapor volume fraction. The cloud cavitation collapses downstream and then becomes a large number of small bubbles, but the vortex is closely related to the cavitation and it can still slide and develop into the downstream. The vortex structurewhich is visually based on theQ-criterion can be seen in Figs.6(b1)-6(b8).

    Fig.7 (Color online) The predicted vapor volume fraction, vortex stretching term, vortex dilatation and baroclinic torque term in one typical cycle.σ=0.8,α=8o,U∞=10 m/s

    Att=1/8Tand 2/8Tnear the T.E. in Fig.6, the shed vapor cloud has a clearly horseshoe-shaped vortex structure indicated by the arrows. In Figs.6(b1)-6(b7), the horseshoe-shaped structure rolls down alongthe suction face of the hydrofoil and slides towards the hydrofoil wake as the red dash line shows, and the evolution of the same large-scale structure can be seen. The large-scale structure intensity becomes weaker as it progresses further downstream. Note that by the isosurface of theQ-criterion, the horseshoe-shaped structure is generated in Fig.6(b1) after the cavity shedding due to the re-entrant jet, and the main flow transports the horseshoe-shaped structure towards the downstream. Gopalan and Katz[16]show that the horseshoe-shaped vortex is formed from the initial rounded cavity structure which shrinks due to the increase of the local pressure. The head of the horseshoe-shaped vortex structure travels further downstream than the two legs and therefore, the horseshoeshaped vortex structure is at an oblique angle with the horizontal plane. The initial horseshoe-shaped vortex structure is thicker and complicated above the face near the trailing edge, but with the time passing, the vortex structure becomes more and more weaker. Finally, the vortex structure collapses initially at the head due to the local high pressure. The sheet cavity shedding and the cloud cavity have a substantial influence on the vortex structure. The large-scale vortex structure gives rise to a very complex and turbulent flow, which will be shown in Section 2.3.

    The transient evolution of the cavity shedding and the horse-shaped vortex were also observed and illustrated in other investigations of different hydrofoils. Ji et al.[13]studied the twisted hydrofoil and found that the cavitation is also evolved due to the effects of both the re-entrant jet and a pair of sideentrant jets. It is noted that the vapor cloud will generate a more clear and dramatical horse-shaped vortex at the twist and large attack angle in the middle area. Another study of the NACA66 hydrofoil by Ji et al.[21]is also carried out with the LES method. Both the numerical and experimental investigations demonstrate that the highly turbulent vapor cloud with a 3-D structure is developed from a 2-D smooth pocket of vapor.

    As is shown in Fig.6, the highly unsteady flow in the cloud cavitation can result in a very complicated vortex structures. In order to investigate the cavitationvortex interaction, the vorticity transport equation in a variable density flow is employed as follows

    In this equation, the vorticity consists of three terms. The vortex stretching term (ω??)Vdenotes the stretching and the tilting of a vortex due to the velocity gradients, the vortex dilatation term(ω??)Vcomes from the volumetric expansion/contraction, and the baroclinic torque (?ρm×?p)/ρm2represents the generation of vorticity due to the misaligned pressure and the density gradients. The last term on the right side indicates the viscous diffusion of the vorticity, and it can be neglected in a high Reynolds number flow[34].

    Figure 7 presents the predicted vapor volume fraction and the contours of the terms of the vortex stretching, the vortex dilatation and the baroclinic torque at one typical cycle as in Fig.4. Att=1/8T, the three terms show very different levels of the vorticity transportation. The vortex stretching term is dominant near the L.E., as a thin attached cavity appears near the L.E. With the attached cavity development and the cavity shedding occurring as shown by other pictures, the vortex stretching term is always an important part due to the velocity gradients. In Figs.7(b) and 7(c), the vortex becomes thicker with the increase of the length of the attached cavity, and the interface of the vapor and the liquid has a larger magnitude than in the other area of the vortex dilatation term. When the re-entrant jet develops towards the L.E., a cavity shedding occurs. At the same time, in Figs.7(e)-7(h), the vortex dilatetion term becomes much more significant than that in Figs.7(b) and 7(c), due to the fact that the local density changes dramatically. It is equal to zero in the non-cavitation region without the volumetric expansion/contraction. Although it is not so remarkable as compared with the vortex stretching and dilatation terms, the baroclinic torque term is very important that it might be the main mechanism for the production of the cavitating vortex as illustrated by Ji et al.[21]. These pictures show the strong correlation between the cavitation and the vortex, and more details and the mechanism might be found in a future research.

    2.3The cavitating flow structure and the turbulencecavitation interactions

    The normalized ensemble averaged velocity and Reynolds stress may illustrate the features of the turbulence cavitating flow structure more directly. The measured and predicted normalized ensemble averaged amplitudes of the axial velocity are compared at the selected monitoring locations along the foil, as shown in Fig.8. It can be seen that the numerical predictions are in good agreement with the experimental results, although the differences increase near the trailing edge (x/C=0.8)due to the sheet cavity shedding and the complicated vortex. The vertical location indicates that the velocity increases to the normal value. It can be seen that the vertical location moves upward with the selected monitoring locations approaching the trailing edge, which shows that the thickness of the boundary layer significantly increases due to the unsteady cavitation. Note that the normalized ensemble averaged velocity is negative near the foil face in Figs.8(a) and 8(b), and the height that the velocityis negative increases clearly in Figs.8(c) and 8(d), namely, there is an apparent re-entrant jet due to the adverse pressure gradient.

    Fig.8 Comparisons of the measured and predicted normalized ensemble averaged amplitudes of the axial velocity at the selected monitoring locations along the foil.σ=0.8,α=8o,=10 m/s

    The running average of the Reynolds stress is shown in Fig.9, where the Points A and B denote the locationsx/C=0.4andx/C=0.8, respectively, as shown in Fig.10(a). It can be seen that the convergence is achieved after 6 000 samples. These data also suggest that the magnitudes of the Reynolds stress remain little fluctuanted after the convergence, but this does not affect the turbulent flow features to be determined in Fig.10. Figure 10 shows the distributions of the normal Reynolds stressesu′u′andv′v′, and the shear Reynolds stressu′v′, whereu′andv′are the horizontal and vertical components of the turbulent velocity fluctuations.

    Fig.9 (Color online) Running average of Reynolds stress at random sample points at selected locations along the foil.σ=0.8,α=8o,U∞=10 m/s

    As seen in Fig.10, the magnitudes of the normal and shear Reynolds stresses are significantly higher around the suction side than in other area. This indicates that the turbulence velocity fluctuations mainly locate in the region where the cavity exists. The turbulence in the region away from the hydrofoil is small. It is interesting to note that the horizontal and vertical Reynolds stresses show completely different distributions along the foil face and in the wake. In Fig.10(a), the horizontal Reynolds stress′is dominant above the foil face, as is demonstrated in Fig.9 that the magnitude of Reynolds stress′at Point A is clearly larger than that at Point B. It can be seen that the turbulence velocity fluctuations are restricted to a narrow height range in the vertical direction. This shows that the detached cavity is only formed for a short time and is not transformed into a large-scale cavity structure as shown in Fig.5. Therefore, the turbulent fluctuation is restricted to a small area. Conversely, the vertical Reynolds stress shows its predominance in the wake over other area in Fig.10(c). This indicates that the shedding of the sheet cavity and the horse-shaped vortex structure produce a high turbulence fluctuation in the wake. Note that in the cavityclosure and the downstream regions, there is a remarkable increase of the turbulence as is seen in Figs.10(a)-10(c). Again, this strongly illustrates that the cloud cavity and the large-scale vortex structure in the wake are the cause of the large turbulent fluctuation. The similar conclusions were reached by Laberteaux and Ceccio[42]. Therefore, the turbulence level of the flow is clearly related to the cavitation. The turbulence velocity fluctuation is greatly influenced by the cavity shedding and the vortex structures.

    Fig.10 (Color online) Distributions of Reynolds stresses.σ= 0.8,α=8o,=10 m/s

    3. Conclusions

    The LES coupled with a homogeneous cavitation model is carried out for the unsteady turbulent attached cavitating flow around a Clark-Y hydrofoil. The quasi-periodic behaviors of the lift coefficient and the cavity volume, the evolution of the unsteady cavitation, the cavity vortex structure, the cavitating flow structure and the turbulence-cavitation interactions are examined. The main conclusions are as follows:

    (1) The predicted lift coefficient and the cavity volume show a distinct quasi-periodic process with the cavitation shedding and they vary dramatically as the pressure distribution changes sharply by the sheet cavity shedding. The LES results of the cavity volume are generally in good agreement with the experiment measured values.

    (2) The time evolution of the unsteady cavitation predicted with the LES method is in good agreement with the experimental results, which shows that the LES method can accurately capture the features of the unsteady cavitation transient behavior. The vortex shape based on theQ-criterion is much more complicated than the vapor volume fraction. The horseshoe-shaped structure is generated after the cavity shedding due to the re-entrant jet, and the horseshoeshaped structure can roll down along the suction face of the hydrofoil and slide towards the hydrofoil wake. The sheet cavity shedding and the cloud cavity have a substantial influence on the vortex structure. The largescale vortex structure gives rise to a very complex and turbulent flow. The vortex stretching term becomes much more predominant when the highly unstable cavitation becomes more severe. The baroclinic torque term is very important and it might be the main mechanism for the production of the cavitating vortex.

    (3) The predicted normalized ensemble averaged amplitudes of the axial velocity are in good agreement with the measured ones. The ensemble averaged velocity shows that the attached cavitation and the cloud cavitation shedding are largely responsible for the increase of the thickness of the boundary layer. The distribution of the normalized ensemble averaged velocity shows the existence of an apparent re-entrant jet. There is a remarkable increase of the turbulence due to the cloud cavitation and the cavity vortex structure in the cavity closure and wake regions. The turbulence level of the flow is clearly related to the cavitation, and the turbulence velocity fluctuation is greatly influenced by the cavity shedding and the cloud cavitation.

    The Verification is an important part in the CFD. However, there is no clear guide about the Verification in the cavitation research. Our previous paper[21]indicates that a sufficient accuracy can be obtained with a reasonable mesh resolution in a fully wetted model with the LES method. The grid effects on the cavitation remain unknown, and an inconsistent conclusion might be reached if we compare the results obtained with three different grids with the process as the ITTC standard for verification and validation in CFD (V&V in CFD)[43,44]about ship hydrodynamics. What we should do is to investigate the relationship between the subgrid-scale turbulence model and the grids. Therefore, the grid independency for the LES method is a difficult and challenging issue for us[45]. We are trying to conduct more investigations on the verification research in the cavitation study with the LES method in our future study.

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

    [2] Huang B., Zhao Y., Wang G. Large eddy simulation of turbulent vortex-cavitation interactions in transient sheet/ cloud cavitating flow [J].Computers and Fluids, 2014, 92(3): 113-124.

    [3] Wang G., Ostoja-Starzewski M. Large eddy simulation of a sheet/cloud cavitation on a NACA0015 hydrofoil [J].Applied Mathematical Modelling, 2007, 31(3): 417-447.

    [4] Dai S., Younis B. A., Sun L. Large-eddy simulations of cavitation in a square surface cavity [J].Applied Mathematical Modelling, 2014, 38(23): 5665-5683.

    [5] Li L., Li B., Hu Z. et al. Large eddy simulation of unsteady shedding behavior in cavitating flows with time-average validation [J].Ocean Engineering, 2016, 125: 1-11.

    [6] Passandideh-Fard M., Roohi E. Coalescence collision of two droplets: bubble entrapment and the effects of important parameters [C].14th Annual (International) Mechanical Engineering Conference. Isfahan, Iran, 2006.

    [7] Roohi E., Zahiri A. P., Passandideh-Fard M. Numerical simulation of cavitation around a two-dimensional hydrofoil using VOF method and LES turbulence model [J].Applied Mathematical Modelling, 2013, 37(9): 6469-6488.

    [8] Gnanaskandan A., Mahesh K. Large eddy simulation of the transition from sheet to cloud cavitation over a wedge [J].International Journal of Multiphase Flow, 2016, 83: 86-102.

    [9] Le Q., Franc J. P., Michel J. M. Partial cavities: Global behavior and mean pressure distribution [J].Journal of Fluids Engineering, 1993, 115(2): 243-248.

    [10] Pham T. M., Larrarte F., Fruman D. H. Investigation of unsteady sheet cavitation and cloud cavitation mechanisms [J].Journal of Ffluids Engineering, 1999, 121(2): 289-296.

    [11] Zhang L. X., Khoo B. C. Computations of partial and super cavitating flows using implicit pressure-based algorithm (IPA) [J].Computers and Fluids, 2013, 73: 1-9.

    [12] Pendar M. R., Roohi E. Investigation of cavitation around 3D hemispherical head-form body and conical cavitators using different turbulence and cavitation models [J].Ocean Engineering, 2016, 112: 287-306.

    [13] Ji B., Luo X. W., Peng X. X. et al. Three-dimensional large eddy simulation and vorticity analysis of unsteady cavitating flow around a twisted hydrofoil [J].Journal of Hydrodynamics, 2013, 25(4): 510-519.

    [14] Park S., Rhee S. H. Numerical analysis of the three-dimensional cloud cavitating flow around a twisted hydrofoil [J].Fluid Dynamics Research, 2012, 45(1): 015502.

    [15] Sedlar M., Ji B., Kratky T. et al. Numerical and experimental investigation of three-dimensional cavitating flow around the straight NACA2412 hydrofoil [J].Ocean Engineering, 2016, 123: 357-382.

    [16] Gopalan S., Katz J. Flow structure and modeling issues in the closure region of attached cavitation [J].Physics of Fluids, 2000, 12(4): 895-911.

    [17] Smagorinsky J. General circulation experiments with the primitive equations: I. The basic experiment [J].Monthly Weather Review, 1963, 91(3): 99-164.

    [18] Moin P. Advances in large eddy simulation methodology for complex flows [J].International Journal of Heat and Fluid Flow, 2002, 23(5): 710-720.

    [19] Piomelli U. Large-eddy simulation: Achievements and challenges [J].Progress in Aerospace Sciences, 1999, 35(4): 335-362.

    [20] Luo X. W., Ji B., Peng X. X. et al. Numerical simulation of cavity shedding from a three-dimensional twisted hydrofoil and induced pressure fluctuation by large-eddy simulation [J].Journal of Fluids Engineering, 2012, 134(4): 379-389.

    [21] Ji B., Luo X. W., Arndt R. E. A. et al. Large eddy simulation and theoretical investigations of the transient cavitating vortical flow structure around a NACA66 hydrofoil [J].International Journal of Multiphase Flow, 2015, 68: 121-134.

    [22] Dreyer M., Decaix J., Münch-Alligné C. et al. Mind the gap: A new insight into the tip leakage vortex using stereo-PIV [J].Experiments in Fluids, 2014, 55(11): 1-13.

    [23] 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.

    [24] Hunt J. C. R., Wray A. A., Moin P. Eddies, streams, convergence zones in turbulent flows [C].Proceedings of the Summer Program 1988 in its Studying Turbulence Using Numerical Simulation Databases. California, USA, 1988, 2: 193-208.

    [25] Kravtsova A. Y., Markovich D. M., Pervunin K. S. et al. High-speed visualization and PIV measurements of cavitating flows around a semi-circular leading-edge flat plate and NACA0015 hydrofoil [J].International Journal of Multiphase Flow, 2014, 60: 119-134.

    [26] Yu X., Huang C., Du T. et al. Study of characteristics of cloud cavity around axisymmetric projectile by large eddy simulation [J].Journal of Fluids Engineering, 2014, 136(5): 051303.

    [27] Roohi E., Pendar M. R., Rahimi A. Simulation of threedimensional cavitation behind a disk using various turbulence and mass transfer models [J].Applied Mathematical Modelling, 2016, 40(1): 542-564.

    [28] Wang Y., Wu X., Huang C. et al. Unsteady characteristics of cloud cavitating flow near the free surface around an axisymmetric projectile [J].International Journal of Multiphase Flow, 2016, 85: 48-56.

    [29] Park S., Rhee S. H. Comparative study of incompressible and isothermal compressible flow solvers for cavitating flow dynamics [J].Journal of Mechanical Science and Technology, 2015, 29(8): 3287-3296.

    [30] Chen G., Wang G., Hu C. et al. Combined experimental and computational investigation of cavitation evolution and excited pressure fluctuation in a convergent–divergent channel [J].International Journal of Multiphase Flow, 2015, 72: 133-140.

    [31] Ji B., Luo X., Wu Y. et al. Numerical investigation of three-dimensional cavitation evolution and excited pressure fluctuations around a twisted hydrofoil [J].Journal of Mechanical Science and Technology, 2014, 28(7): 2659-2668.

    [32] Huang B., Young Y. L., Wang G. et al. Combined experimental and computational investigation of unsteady structure of sheet/cloud cavitation [J].Journal of Fluids Engineering, 2013, 135(7): 071301.

    [33] Iyer C. O., Ceccio S. L. The influence of developed cavitation on the flow of a turbulent shear layer [J].Physics of Fluids, 2002, 14(10): 3414-3431.

    [34] Dittakavi N., Chunekar A., Frankel S. Large eddy simulation of turbulent-cavitation interactions in a Venturi nozzle [J].Journal of Fluids Engineering, 2010, 132(12): 121301.

    [35] Decaix J., Goncalvès E. Investigation of three-dimensional effects on a cavitating Venturi flow [J].International Journal of Heat and Fluid Flow, 2013, 44: 576-595.

    [36] Schnerr G. H., Sauer J. Physical and numerical modelingof unsteady cavitation dynamics[C].Fourth International Conference on Multiphase Flow. New Orleans, USA, 2001.

    [37] Nicoud F., Ducros F. Subgrid-scale stress modelling based on the square of the velocity gradient tensor [J].Flow Turbulence and Combustion, 1999, 62(3):183-200.

    [38] Sagaut P., Lee Y. T. Large eddy simulation for incompressible Flows: An introduction. scientific computation series [J].Applied Mechanics Reviews, 2002, 55(6): 1745-1746.

    [39] Barth T. J., Jespersen D. C. The design and application of upwind schemes on unstructured meshes [C].27th Aerospace Sciences Meeting, Reno, NV, USA, 1989.

    [40] Wang G., Senocak I., Shyy W. et al. Dynamics of attached turbulent cavitating flow [J].Progress in Aerospace Sciences, 2001, 37(6): 551-581.

    [41] Huang B. Physical and numerical investigation of unsteady cavitating fl ows [D]. Doctoral Thesis, Beijing, China: Beijing Institute of Technology, 2012(in Chinese).

    [42] Laberteaux K. R., Ceccio S. L. Partial cavity flows. Part 1. Cavities forming on models without spanwise variation [J].Journal of Fluid Mechanics, 2001, 431: 1-41.

    [43] ITTC QM Procedure. Uncertainty analysis in CFD verification and validation methodology and procedures [R]. 7.5-03-01-01, 2002.

    [44] ITTC QM Procedure. Uncertainty analysis in CFD, examples for resistance and flow [R]. 7.5-03-02-01, 2002.

    [45] Xing T. A general framework for verification and validation of large eddy simulations [J].Journal of Hydrodynamics, 2015, 27(2): 163-175.

    * Project supported by the National Natural Science Foundation of China (Grant Nos. 51576143, 11472197).

    Biography:Bin Ji (1982-), Male, Ph. D., Associate Professor

    Xin-ping Long, E-mail:xplong@whu.edu.cn

    猜你喜歡
    龍云新平
    幼兒園里歡樂多
    幼兒園(2021年18期)2021-12-06 02:45:42
    小螞蟻去游玩
    幼兒園(2021年16期)2021-12-06 01:06:48
    老腔唱新歌
    金秋(2021年22期)2021-03-10 07:59:16
    出滇抗戰(zhàn)時期龍云對滇軍的治理研究
    創(chuàng)造(2020年6期)2020-11-20 05:58:42
    讓蘑菇
    幼兒園(2020年3期)2020-03-27 07:00:07
    劉新平 油畫作品
    新疆藝術(2019年3期)2019-07-17 11:26:44
    An integral calculation approach for numerical simulation of cavitating flow around a marine propeller behind the ship hull *
    URANS simulations of the tip-leakage cavitating flow with verification and validation procedures *
    你總是給我力量
    當代音樂(2018年11期)2018-06-30 14:43:16
    Some notes on numerical simulation and error analyses of the attached turbulent cavitating flow by LES *
    汤姆久久久久久久影院中文字幕| 久久99蜜桃精品久久| 黄片无遮挡物在线观看| av免费观看日本| www.色视频.com| 欧美精品国产亚洲| 女人十人毛片免费观看3o分钟| 色视频在线一区二区三区| 欧美成人a在线观看| 精品一区二区三区视频在线| a级毛片免费高清观看在线播放| av在线播放精品| 最近中文字幕高清免费大全6| 男人和女人高潮做爰伦理| 啦啦啦啦在线视频资源| 特大巨黑吊av在线直播| 成人午夜精彩视频在线观看| 黑丝袜美女国产一区| 亚洲精品自拍成人| 一区二区av电影网| 大码成人一级视频| 国国产精品蜜臀av免费| 一级毛片我不卡| 国内精品宾馆在线| 色视频www国产| 涩涩av久久男人的天堂| 高清视频免费观看一区二区| 永久网站在线| 国产精品成人在线| 亚洲av男天堂| 午夜激情久久久久久久| 亚洲欧美日韩卡通动漫| 一级a做视频免费观看| 久久久成人免费电影| 久久精品久久久久久噜噜老黄| 3wmmmm亚洲av在线观看| 亚洲欧洲国产日韩| 有码 亚洲区| 欧美性感艳星| 亚洲精品日本国产第一区| av国产精品久久久久影院| 中文字幕精品免费在线观看视频 | 国产 精品1| 国产亚洲精品久久久com| 亚洲综合精品二区| 乱系列少妇在线播放| 亚洲精品成人av观看孕妇| 成人二区视频| 女的被弄到高潮叫床怎么办| 精品人妻一区二区三区麻豆| 亚洲av中文字字幕乱码综合| 久久女婷五月综合色啪小说| 国产黄色免费在线视频| 黑人高潮一二区| 国产av国产精品国产| 男人和女人高潮做爰伦理| 美女高潮的动态| 亚洲va在线va天堂va国产| 青青草视频在线视频观看| 亚洲欧美日韩卡通动漫| 熟女人妻精品中文字幕| 亚洲精品国产av蜜桃| 午夜免费观看性视频| 国产成人一区二区在线| 精品国产乱码久久久久久小说| 久久精品久久精品一区二区三区| 97在线人人人人妻| 1000部很黄的大片| 2018国产大陆天天弄谢| 韩国av在线不卡| 久久亚洲国产成人精品v| 汤姆久久久久久久影院中文字幕| 91精品国产国语对白视频| 高清毛片免费看| .国产精品久久| 久久久成人免费电影| av国产精品久久久久影院| 丰满迷人的少妇在线观看| 日本一二三区视频观看| av免费在线看不卡| 欧美精品国产亚洲| 少妇的逼好多水| 久久97久久精品| 精品人妻熟女av久视频| 国产91av在线免费观看| 久久精品熟女亚洲av麻豆精品| 精品人妻偷拍中文字幕| 人人妻人人添人人爽欧美一区卜 | 欧美激情极品国产一区二区三区 | 久久国产精品男人的天堂亚洲 | 如何舔出高潮| 老司机影院毛片| 国产精品一及| 少妇 在线观看| 久久久精品免费免费高清| 亚洲精品视频女| 欧美一区二区亚洲| 尤物成人国产欧美一区二区三区| 午夜福利网站1000一区二区三区| 午夜免费鲁丝| 极品少妇高潮喷水抽搐| 国内少妇人妻偷人精品xxx网站| 国产亚洲最大av| 99热这里只有是精品在线观看| 亚洲性久久影院| 亚洲精品日韩av片在线观看| 一本久久精品| 大又大粗又爽又黄少妇毛片口| 国内少妇人妻偷人精品xxx网站| 舔av片在线| 亚洲成人一二三区av| www.av在线官网国产| 亚洲综合色惰| 一本一本综合久久| 国产精品99久久久久久久久| 各种免费的搞黄视频| 中文乱码字字幕精品一区二区三区| 成人午夜精彩视频在线观看| 视频区图区小说| 性色avwww在线观看| 国产色婷婷99| av免费观看日本| 国模一区二区三区四区视频| 蜜桃在线观看..| 国产 一区精品| 欧美成人a在线观看| 国产成人精品一,二区| 国内少妇人妻偷人精品xxx网站| 一级二级三级毛片免费看| 亚洲成人一二三区av| 国产成人freesex在线| 国产精品一区二区性色av| 亚洲欧洲日产国产| a 毛片基地| 欧美日韩亚洲高清精品| 国产精品久久久久成人av| 国产精品99久久久久久久久| 中国美白少妇内射xxxbb| 麻豆成人av视频| 国产男女超爽视频在线观看| 国产一区二区在线观看日韩| 高清av免费在线| 美女xxoo啪啪120秒动态图| 一本久久精品| 极品教师在线视频| 国产精品人妻久久久久久| 国产成人精品福利久久| 国产在线一区二区三区精| 美女中出高潮动态图| 网址你懂的国产日韩在线| 国产精品一区二区在线观看99| 寂寞人妻少妇视频99o| 日本av免费视频播放| 亚洲国产欧美在线一区| 日日摸夜夜添夜夜爱| 最近的中文字幕免费完整| 国产成人精品一,二区| 成年美女黄网站色视频大全免费 | 伦理电影大哥的女人| 永久网站在线| 成人国产av品久久久| 亚洲欧美精品自产自拍| 日产精品乱码卡一卡2卡三| 国产成人精品婷婷| 亚洲国产最新在线播放| 欧美成人精品欧美一级黄| 日韩一区二区视频免费看| 你懂的网址亚洲精品在线观看| av在线观看视频网站免费| 大码成人一级视频| 简卡轻食公司| 18+在线观看网站| 人体艺术视频欧美日本| 精品国产乱码久久久久久小说| 久久久国产一区二区| 成人18禁高潮啪啪吃奶动态图 | 国产一区有黄有色的免费视频| 不卡视频在线观看欧美| 日本色播在线视频| 男人爽女人下面视频在线观看| 亚洲国产精品999| 亚洲综合色惰| 香蕉精品网在线| 永久网站在线| 天天躁日日操中文字幕| 日韩成人伦理影院| 国产精品一及| 18禁动态无遮挡网站| 午夜福利在线观看免费完整高清在| 久久99热这里只有精品18| 国产精品偷伦视频观看了| 国产精品蜜桃在线观看| 国产深夜福利视频在线观看| 少妇的逼水好多| 高清黄色对白视频在线免费看 | 深夜a级毛片| 人人妻人人添人人爽欧美一区卜 | 看十八女毛片水多多多| 激情五月婷婷亚洲| 国产高清三级在线| 又大又黄又爽视频免费| 久久精品久久久久久久性| 99视频精品全部免费 在线| 国产视频首页在线观看| 日韩视频在线欧美| 国产伦精品一区二区三区视频9| 国产人妻一区二区三区在| 日本黄色日本黄色录像| 亚洲精品国产成人久久av| 国产精品久久久久久久久免| 蜜桃在线观看..| 国产成人freesex在线| 国产精品麻豆人妻色哟哟久久| 久久人妻熟女aⅴ| 精品久久国产蜜桃| 18禁裸乳无遮挡免费网站照片| 亚洲国产毛片av蜜桃av| 寂寞人妻少妇视频99o| 爱豆传媒免费全集在线观看| av卡一久久| 欧美高清性xxxxhd video| 三级国产精品欧美在线观看| 色网站视频免费| 国产亚洲欧美精品永久| 久久精品夜色国产| 久久鲁丝午夜福利片| 日韩伦理黄色片| 美女内射精品一级片tv| 国产色爽女视频免费观看| 99精国产麻豆久久婷婷| 中文精品一卡2卡3卡4更新| 麻豆精品久久久久久蜜桃| 国产精品一区二区三区四区免费观看| www.色视频.com| 国产成人aa在线观看| 女人久久www免费人成看片| 亚洲av男天堂| 久久久久精品性色| 午夜视频国产福利| 精华霜和精华液先用哪个| 啦啦啦中文免费视频观看日本| 高清在线视频一区二区三区| 亚洲av欧美aⅴ国产| av国产精品久久久久影院| 又黄又爽又刺激的免费视频.| 寂寞人妻少妇视频99o| 久久综合国产亚洲精品| 国产v大片淫在线免费观看| 波野结衣二区三区在线| 五月玫瑰六月丁香| 最近的中文字幕免费完整| 不卡视频在线观看欧美| av播播在线观看一区| 日本av手机在线免费观看| 女的被弄到高潮叫床怎么办| www.av在线官网国产| 日韩三级伦理在线观看| 欧美丝袜亚洲另类| 日韩国内少妇激情av| 国产真实伦视频高清在线观看| 精品久久久久久久末码| tube8黄色片| 午夜福利网站1000一区二区三区| 99九九线精品视频在线观看视频| 精品99又大又爽又粗少妇毛片| .国产精品久久| 日本猛色少妇xxxxx猛交久久| 欧美3d第一页| 天天躁夜夜躁狠狠久久av| 亚洲精品乱久久久久久| 免费在线观看成人毛片| 五月开心婷婷网| 色吧在线观看| 在线免费观看不下载黄p国产| 国产精品一区二区三区四区免费观看| 亚州av有码| 久久99热6这里只有精品| 久久久久久久久久人人人人人人| 亚洲精品国产成人久久av| 国产欧美日韩一区二区三区在线 | 国产精品免费大片| 欧美3d第一页| 亚洲经典国产精华液单| www.色视频.com| 能在线免费看毛片的网站| 国产欧美日韩一区二区三区在线 | 欧美xxxx性猛交bbbb| 人人妻人人看人人澡| 久久99蜜桃精品久久| 能在线免费看毛片的网站| 久久久精品免费免费高清| 交换朋友夫妻互换小说| 噜噜噜噜噜久久久久久91| 麻豆乱淫一区二区| 成人亚洲精品一区在线观看 | 日韩 亚洲 欧美在线| 在线播放无遮挡| 国产精品一及| 青春草国产在线视频| h日本视频在线播放| 免费久久久久久久精品成人欧美视频 | 亚洲精品亚洲一区二区| 欧美性感艳星| 香蕉精品网在线| 人人妻人人澡人人爽人人夜夜| av国产免费在线观看| 美女内射精品一级片tv| 国产免费又黄又爽又色| 亚洲欧美精品专区久久| 亚洲高清免费不卡视频| 久久久a久久爽久久v久久| 日韩一本色道免费dvd| 国产欧美另类精品又又久久亚洲欧美| 男人狂女人下面高潮的视频| 国产免费一级a男人的天堂| 久久ye,这里只有精品| 成人二区视频| 久久人妻熟女aⅴ| 久久ye,这里只有精品| 欧美高清性xxxxhd video| 亚洲精品日韩在线中文字幕| 少妇猛男粗大的猛烈进出视频| 免费少妇av软件| 免费高清在线观看视频在线观看| 王馨瑶露胸无遮挡在线观看| 午夜免费男女啪啪视频观看| 中文欧美无线码| 性色avwww在线观看| 少妇的逼水好多| 观看av在线不卡| 国产精品久久久久成人av| 久热这里只有精品99| 婷婷色综合大香蕉| 亚洲四区av| 亚洲国产精品999| 日日撸夜夜添| 国产成人aa在线观看| 成人亚洲欧美一区二区av| 国产精品三级大全| 欧美高清性xxxxhd video| 久久久亚洲精品成人影院| 久久97久久精品| 国产精品麻豆人妻色哟哟久久| 午夜日本视频在线| 91久久精品国产一区二区成人| 又黄又爽又刺激的免费视频.| 免费观看a级毛片全部| 久久国内精品自在自线图片| 人妻 亚洲 视频| 最黄视频免费看| 少妇熟女欧美另类| 精品国产一区二区三区久久久樱花 | 国产伦精品一区二区三区视频9| 亚洲国产色片| 亚洲无线观看免费| 黄片wwwwww| 97热精品久久久久久| 另类亚洲欧美激情| 婷婷色综合大香蕉| 亚洲av免费高清在线观看| 精品一品国产午夜福利视频| 国产 一区精品| 亚洲国产最新在线播放| 亚洲三级黄色毛片| 国产精品成人在线| 午夜老司机福利剧场| av在线app专区| 精品熟女少妇av免费看| 伦理电影免费视频| 国产色爽女视频免费观看| 91午夜精品亚洲一区二区三区| 久久久精品94久久精品| 亚洲精品成人av观看孕妇| 国产一区二区三区综合在线观看 | av国产免费在线观看| 国产免费一区二区三区四区乱码| 一个人免费看片子| 久久精品国产自在天天线| 久久韩国三级中文字幕| 国产高清三级在线| 欧美日韩精品成人综合77777| 国产在视频线精品| 亚洲av中文字字幕乱码综合| 久久久亚洲精品成人影院| 黄色怎么调成土黄色| 亚洲人成网站高清观看| 国产永久视频网站| 午夜老司机福利剧场| 观看免费一级毛片| 天堂中文最新版在线下载| 亚洲精品国产成人久久av| 如何舔出高潮| 高清黄色对白视频在线免费看 | 超碰av人人做人人爽久久| 极品教师在线视频| 久久女婷五月综合色啪小说| 亚洲熟女精品中文字幕| 高清av免费在线| 青春草国产在线视频| 最近中文字幕2019免费版| 亚洲国产欧美在线一区| 久久久久精品久久久久真实原创| 国产精品一区www在线观看| 亚洲性久久影院| 成人国产麻豆网| av不卡在线播放| 蜜桃在线观看..| 日韩制服骚丝袜av| 狂野欧美激情性bbbbbb| 亚洲国产日韩一区二区| 美女国产视频在线观看| 下体分泌物呈黄色| 在线精品无人区一区二区三 | 最后的刺客免费高清国语| 日日啪夜夜爽| .国产精品久久| 多毛熟女@视频| 香蕉精品网在线| 久久久亚洲精品成人影院| 一级毛片我不卡| 伊人久久精品亚洲午夜| 夫妻午夜视频| 国产日韩欧美亚洲二区| 国产精品偷伦视频观看了| 五月天丁香电影| 91精品一卡2卡3卡4卡| 在线观看三级黄色| 国产爽快片一区二区三区| 亚洲一区二区三区欧美精品| 亚洲成人手机| 日韩视频在线欧美| 久久久久久久久大av| 国产亚洲一区二区精品| 亚洲欧美精品自产自拍| 一区二区三区免费毛片| 成人午夜精彩视频在线观看| 国产精品成人在线| av国产久精品久网站免费入址| 免费人成在线观看视频色| 国产一区二区在线观看日韩| 最近2019中文字幕mv第一页| 欧美日韩一区二区视频在线观看视频在线| 亚洲激情五月婷婷啪啪| 综合色丁香网| 又黄又爽又刺激的免费视频.| 免费高清在线观看视频在线观看| 国产亚洲精品久久久com| 国产精品一区二区三区四区免费观看| 麻豆国产97在线/欧美| 涩涩av久久男人的天堂| 日韩伦理黄色片| 国产精品蜜桃在线观看| 最黄视频免费看| 国产69精品久久久久777片| 身体一侧抽搐| 久久午夜福利片| 久久精品熟女亚洲av麻豆精品| 又大又黄又爽视频免费| 亚洲精品国产av成人精品| 国产男女超爽视频在线观看| 欧美性感艳星| 欧美变态另类bdsm刘玥| 菩萨蛮人人尽说江南好唐韦庄| 狠狠精品人妻久久久久久综合| 精品人妻偷拍中文字幕| 一二三四中文在线观看免费高清| 国产91av在线免费观看| 哪个播放器可以免费观看大片| 99热全是精品| 亚洲精品亚洲一区二区| 国产 精品1| 精品国产一区二区三区久久久樱花 | 国产一区二区在线观看日韩| 欧美老熟妇乱子伦牲交| 国产精品爽爽va在线观看网站| 看免费成人av毛片| 一本久久精品| 2022亚洲国产成人精品| 亚洲不卡免费看| 又粗又硬又长又爽又黄的视频| 日日摸夜夜添夜夜添av毛片| 亚洲成人手机| 天天躁夜夜躁狠狠久久av| 欧美3d第一页| 99久久中文字幕三级久久日本| 80岁老熟妇乱子伦牲交| 一二三四中文在线观看免费高清| 在线看a的网站| 久久久久性生活片| 亚洲伊人久久精品综合| 日韩欧美一区视频在线观看 | 亚洲国产成人一精品久久久| 亚洲国产欧美在线一区| 国产精品99久久久久久久久| 亚洲欧美清纯卡通| 亚洲av国产av综合av卡| 精品99又大又爽又粗少妇毛片| 在线观看三级黄色| 一个人看视频在线观看www免费| 精品久久久精品久久久| 亚洲精品自拍成人| 亚洲图色成人| 男女国产视频网站| 亚洲成色77777| 黄色视频在线播放观看不卡| 婷婷色综合大香蕉| 黄色欧美视频在线观看| 午夜免费鲁丝| 亚洲欧美精品自产自拍| 国产亚洲一区二区精品| 少妇人妻久久综合中文| 国产精品人妻久久久久久| 亚洲成人手机| 久久精品人妻少妇| 老熟女久久久| 国产亚洲最大av| 黄色日韩在线| 一级黄片播放器| 最黄视频免费看| 男人添女人高潮全过程视频| 国产黄片视频在线免费观看| 性高湖久久久久久久久免费观看| 久久久久久久精品精品| av在线老鸭窝| 欧美激情国产日韩精品一区| 久久人人爽av亚洲精品天堂 | 在线天堂最新版资源| 欧美成人a在线观看| 五月开心婷婷网| 在线观看免费高清a一片| 国产在线免费精品| 观看美女的网站| 国产探花极品一区二区| av在线播放精品| 色综合色国产| 好男人视频免费观看在线| 国产av精品麻豆| 日本一二三区视频观看| 成人毛片60女人毛片免费| av国产久精品久网站免费入址| 国产欧美亚洲国产| 六月丁香七月| 91狼人影院| 国国产精品蜜臀av免费| kizo精华| 久久影院123| 日韩成人伦理影院| 精品一区二区免费观看| 热re99久久精品国产66热6| 亚洲欧美日韩无卡精品| 日本欧美国产在线视频| 欧美激情极品国产一区二区三区 | 青春草视频在线免费观看| 国产黄片视频在线免费观看| 日韩欧美 国产精品| 久久久久久人妻| 国产男女内射视频| 大香蕉久久网| 日韩 亚洲 欧美在线| 菩萨蛮人人尽说江南好唐韦庄| 日韩在线高清观看一区二区三区| 亚洲欧美清纯卡通| 男女啪啪激烈高潮av片| 韩国高清视频一区二区三区| 精华霜和精华液先用哪个| 中文字幕精品免费在线观看视频 | 亚洲欧美中文字幕日韩二区| 日韩大片免费观看网站| 少妇精品久久久久久久| 久久人人爽人人爽人人片va| 一级爰片在线观看| 一级a做视频免费观看| 亚洲av中文字字幕乱码综合| 国产综合精华液| 性高湖久久久久久久久免费观看| 2021少妇久久久久久久久久久| 亚洲精品亚洲一区二区| 啦啦啦啦在线视频资源| 国产精品无大码| 国产日韩欧美亚洲二区| 亚洲欧洲日产国产| 在线精品无人区一区二区三 | 国产在线视频一区二区| 国语对白做爰xxxⅹ性视频网站| 亚洲欧美精品自产自拍| 一级黄片播放器| 国产真实伦视频高清在线观看| 亚洲国产成人一精品久久久| 哪个播放器可以免费观看大片| 国产精品麻豆人妻色哟哟久久| 成人18禁高潮啪啪吃奶动态图 | av在线播放精品| 永久网站在线| 秋霞在线观看毛片| 精品国产乱码久久久久久小说| 女的被弄到高潮叫床怎么办| 自拍欧美九色日韩亚洲蝌蚪91 | 国产成人freesex在线| 免费观看的影片在线观看| 亚洲av欧美aⅴ国产| 99热6这里只有精品| av国产免费在线观看| 97在线人人人人妻| 九色成人免费人妻av| 亚洲性久久影院| 18禁在线无遮挡免费观看视频| av在线蜜桃| 天堂俺去俺来也www色官网| 在线观看一区二区三区| 国产一区有黄有色的免费视频| 激情五月婷婷亚洲| 内射极品少妇av片p| 99久久精品国产国产毛片| 精品人妻视频免费看| 日本欧美视频一区|