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

    Numerical Simulation of the Flow around Two-dimensional Partially Cavitating Hydrofoils

    2014-07-30 09:54:44FahriCelikYaseminArikanOzdenandSakirBal

    Fahri Celik, Yasemin Arikan Ozden and Sakir Bal

    1. Department of Naval Architecture and Marine Engineering, Y?ld?z Technical University, Istanbul 34349, Turkey

    2. Department of Naval Architecture and Marine Engineering, Istanbul Technical University, Istanbul 34469, Turkey

    1 Introduction1

    Cavitation appears to be an unavoidable phenomenon especially in water devices like pumps, marine propellers and hydrofoils where fast flow regimes exist. It can cause undesirable results like performance losses, structural failure,corrosion, noise and vibration. The prediction of performance losses caused by the sheet cavitation which is common particularly on hydrofoils is very important in their design stages. So the development of computational analyses methods is very important since the cavitation occurrence is extensively unavoidable for modern high speed marine vehicles with hydrofoils.

    In the past the two-dimensional (2D) cavitating hydrofoil flows were basically formulated under linear theories (Tulin,1953; Acosta, 1955; Geurst and Timman, 1956). The linear theory predicts an increase in the cavity size and volume with the increasing of the foil thickness in the same flow conditions. On the other hand the numerical nonlinear surface vorticity method developed by Uhlman (1987)where the cavity surface is obtained by an iterative manner until both the kinematic and dynamic boundary conditions are satisfied on the cavity surface, predicts that the cavity size should decrease with the increasing of the foil thickness.This defect of the linear theory has been corrected by the introduction of the nonlinear leading-edge correction(Kinnas, 1991).

    With the beginning of the use of boundary element methods (BEM), a fast approach for the cavitation analysis has been provided. After the first applications of the boundary element methods (Hess and Smith, 1966; Hess,1973; Hess, 1975), various studies using BEM have been conducted for the flow analysis of 2D non-cavitating hydrofoils (Katz and Plotkin, 2001; Chow, 1979). For the cavity prediction of partial and supercavitating 2D and 3D hydrofoils based on the potential theory through the use of boundary element methods, a nonlinear analysis method has been developed (Kinnas and Fine, 1993).

    Some important studies on the use of the boundary element methods for the flow analysis of 2D and 3D cavitating hydrofoils can be found (Fine and Kinnas, 1993;Kinnas and Fine, 1992; Kinnaset al., 1994; Kinnas, 1998;Kimet al., 1994; Kim and Lee, 1996; Dang and Kuiper,1998; Dang, 2001; Krishnaswamy, 2000). In these studies,non-linear analysis methods for viscous and inviscid flows around 2D and 3D cavitating hydrofoils and the prediction methods for unsteady sheet cavitation around 3D cavitating hydrofoils are presented. The effects of various cavity termination models on 2D and 3D cavitating hydrofoils have been also investigated. In addition the studies (Bal and Kinnas, 2003; Bal, 2007, 2008, 2011) where the iterative BEM methods are applied to the finite depth and wave tank effects on cavitating 2D and 3D hydrofoil analyses(including the surface piercing type of hydrofoils) are other important studies on this subject.

    Another approach for the cavitation simulation methods including the viscous effects is based on the numerical solution of the RANS and Euler equations. A method based on the solution of the Navier-Stokes equations (Deshpandeet al., 1993) and the Euler equations have been applied by Deshpandeet al. (1994). They have predicted that the presence of the viscous effects on a cavitating flow over 2D hydrofoils have little impact on the cavitating region in regard to the Euler analysis. A different study on the analysis of the partial cavity flow on 2D hydrofoils with the use of thek-εturbulence model has been carried out (Dupont and Avellan, 1991) and has been extended (Hirschiet al., 1997)for the cavity prediction on a twisted hydrofoil and a pump impeller.

    CFD methods also offer effective approaches to model the cavitating flows. The cavitating flow can be modeled as steady or unsteady phenomena. An unsteady flow approach has been studied (Berntsenet al., 2001) where the cavity length is calculated with a high level of accuracy. The cavitation erosion risk has been predicted on a 2D hydrofoil(Liet al., 2010) for steady and unsteady conditions and it is concluded that a RANS code, FLUENT, shows promising results with the prediction of unsteady cavitation phenomena.The effects of different turbulence models and wall treatments on the cavitating flow around hydrofoils have been compared by Huanget al. (2010). They have concluded that the RNG turbulence model with enhanced near-wall functioning offers more satisfactory results than those of the other models. To predict the impact of the early stages of the cavity development on the hydrofoils, a CFD model has been used in Hoekstra and Vaz (2009). They have concluded that the re-entrant jet model must be considered as an invalid approach for the modeling of steady partial cavities.

    In this study an iterative nonlinear method based on the potential theory was developed for the cavity prediction around partially cavitating 2D hydrofoils. The present analyses method is based on the method developed by the same authors in Ar?kanet al.(2012). For the 2D cavitation analysis, the NACA 16006, NACA 16012 and NACA 16015 sections were selected to predict the cavity shapes and pressure distributions, including lift and drag coefficients.The analyses were made for two different incoming flow angles of attack (4 and 6 degrees). The results obtained from the present method were compared with those of another potentially based BEM method (PCPAN) (Kinnas and Fine,1993a) and a commercial CFD code (FLUENT).

    2 Mathematical formulation for partially cavitating 2D hydrofoil analysis

    The flow over of the hydrofoil is assumed to be incompressible, inviscid and irrotational according to the potential theory. The submerged 2D hydrofoils are subject to a uniform inflow (U) (see Fig. 1).

    Fig. 1 2D hydrofoil

    The total potential for the flow around the hydrofoil can be described for the 2D case as (Bal and Kinnas, 2003):

    The total and the perturbation velocity potential must satisfy the Laplace equation:

    The following boundary conditions should be satisfied on the foil surface.

    1) Kinematic boundary condition:

    The flow must be tangential to the foil and the cavity surfaces,

    wherenis the unit vector normal to the foil surface and cavity surface.

    2) Dynamic boundary condition:

    The dynamic boundary condition requires that the pressure on the cavity surface is constant and equal to the vaporization pressure. This condition is satisfied by the deformation of the foil geometry so that on the cavity surface the pressure coefficient (-Cp) values are equal to the cavitation number (s),

    whereP∞is the total,Pvis the vaporization andPmis the static pressure.

    3) Kutta condition:

    The velocity in the trailing edge of the hydrofoil must be finite,

    4) Cavitation condition:

    This requires the cavity to close at the cavity trailing edge(Kinnas and Fine, 1993).

    3 Numerical implementation

    An iterative solution method is developed for the prediction of cavitation over 2D hydrofoils. The cavitating flow analyses are carried out in an iterative manner using an analysis method for non-cavitating hydrofoils. For the flow analysis around a 2D hydrofoil, a potential based boundary element method with a constant-strength source and doublet distributions (Dirichlet type of boundary condition) is used.The details of the method can be found in Katz and Plotkin,(2001). For a fixed cavity length and a constant cavitation number, the cavity shape is searched iteratively by deforming the foil surface in the cavity region with small intervals (Dh). The flow around the hydrofoil is analyzed for the new foil geometries (deformed by cavitation) in each of the iteration stages by the BEM method. When the dynamic boundary condition is provided over the cavity surface with an acceptable tolerance value, in the last iteration, the final cavity shape and the pressure distributions on the cavitating hydrofoil are obtained (Fig. 2). For a constant cavitation number, cavity shapes are searched for various fixed cavity lengths on the foil surface, and the appropriate cavity length and shape is selected according to the minimum error criteria where the sum of |-Cp+σ| along the cavity length is the minimum.

    Fig. 2 Cavity shapes at different iteration stages

    3.1 Cavity detachment

    An important parameter regarding the cavitation analysis is the location of the detachment point. While the cavitation begins in the leading edge in cases with a sharp leading edge,the cavity detaches from a point downstream of the laminar boundary layer separation point in cases with a smooth curved leading edge. The location of the cavity detachment point can be determined by using the smooth detachment condition (Brillouin-Villat condition). In this condition it is assumed that the cavity does not intersect with the foil surface at its leading edge and that the pressure on the wetted foil surface in front of the cavity is not lower than the vapor pressure (Krishnaswamy, 2000; Kinnas and Fine,1991). In the present study the detachment point is assumed to be located on the leading edge of the hydrofoil. The search algorithm similar to that of Kinnas and Fine (1991)for the cavity detachment point is under progress.

    3.2 Cavity termination

    The cavity closure region consists of a two phase turbulent zone where a very complex flow occurs. The difficulties in obtaining the flow characteristics in this region require the use of the cavity termination models.There are various cavity termination models to simulate the cavity closure region such as the end plate Riabouchinsky or termination wall model, the re-entrant jet model, Tulin’s spiral vortex model, the open wake model and the pressure recovery model. The most physically realistic cavity termination model is the re-entrant jet model due to its best representation of the streamlines around a cavitating body.Another physically realistic termination model is the termination wall (the Riabouchinsky) model which is similar to the re-entrant jet model. The wall termination model in which the cavity region is closed with a vertical wall is shown in Fig. 3. With both the termination wall model and the re-entrant jet models, a stagnation point which has been observed experimentally, occurs in the trailing edge of the cavity (Krishnaswamy, 2000). Another cavity termination model is the pressure recovery model, and its details are given in the following paragraph; the use of this model has been described by Kinnas and Fine (1993).

    In this study two different termination models are considered for the cavitation analysis of the 2D hydrofoil:First, the pressure recovery model and second the termination wall model.

    Fig. 3 Termination wall model

    The pressure recovery model requires that in the transition zone (T-L), on the cavity surface, the pressure is not equal to the vaporization pressure in order to close the cavity region (Fig. 4).

    Fig. 4 Pressure recovery model for partially cavitating the two-dimensional hydrofoil

    For the application of the pressure recovery termination model the equations below are employed (Kinnas and Fine,1993).

    here,sfis the arc-length of the foil measured beneath the cavity, measured from the cavity’s leading edge. The(0>lt;A>lt;1) andv(v>gt;0) are arbitrary constants. The cavity length is shown withland the abbreviationsDstands for the detachment point,Tfor the beginning of the transition zone andLfor the trailing edge of the cavity. In the present method the cavitation number (σT) is taken as:

    Through the comparison of the pressure coefficient (-Cp)with the new cavitation number (σT), the dynamic boundary condition on the foil surface (D-L) is provided.

    In the present method, for the 2D analysis, the termination wall model is applied automatically if thef(Sf)value in Eq. (7) in the termination model is selected as zero for the entire cavity region (D-L).

    4 CFD Analysis Method

    In the computational method, the governing equations are discretized using a first-order upwind interpolation scheme,and the discretized equations are solved using the SIMPLEC algorithm. The related equations are solved for the flow,vapor and turbulence. The convergence criterion for the solution is selected as 10-6. The mixture model is chosen as the multiphase model and the cavitation model is included in the calculations (Fluent 6.3 User’s Guide, 2006).

    For the mathematical model the incompressible,time-averaged, unsteady mean flow equations of continuity and momentum which are used in the mixture multiphase model are given below (Fluent 6.3 User’s Guide, 2006):

    whereθis the mass averaged velocity,ρmis the mixture density andαkis the volume fraction of phasek.

    herenis the number of the phases,Fis a body force,μmis the viscosity of the mixture andvdr.kis the drift velocity for secondary phasek.

    Thek-ε(RNG) turbulence model with enhanced wall treatment is selected. Turbulence kinetic energy,k(Eq. (11))and its dissipation,ε(Eq. (12)) for the mixture model are:

    In these equations,Gkrepresents the generation of the turbulence kinetic energy due to the mean velocity gradients,Gbis the generation of the turbulence kinetic energy due to buoyancy,YMis the contribution of the fluctuating dilatation in the compressible turbulence to the overall dissipation rate,αkandαεare the inverse effective Prandtl numbers forkandε,SkandSεare the source terms.

    The cavitation effects for the two-phase mixture flows are included by using the standard cavitation model in Fluent.The cavitation model implemented in Fluent is based on the“full cavitation model” (Singhalet al., 2001). The model can be used with all of the turbulence models and with the mixture multiphase model. For the calculations, the first phase is selected as liquid water and the second phase is selected as water vapor. The vapor transport equation with“f” as the vapor mass fraction is given as (Hoekstra and Vaz,2009):

    whereρis the mixture density,vvis the velocity vector of the vapor phase,γis the effective exchange coefficient,ReandRcare the vapor generation and condensation rate terms.The final phase rate expressions obtained after accounting for the effects of the turbulence-induced pressure fluctuations and non condensable gases are forp>lt;pv:

    The suffixes (”l” and “v”) denote the liquid and the vapor phases,σis the surface tension coefficient of the liquid,pvis the liquid vaporization pressure andCeandCcare empirical constants.

    5 Numerical results and discussions

    In this section some numerical results of the present method for the cavity prediction on 2D hydrofoils are presented. The calculations where the cavity shapes and pressure distributions are determined for the 2D hydrofoils are carried out for three different sections (NACA 16006,NACA 16012 and NACA 16015) at two angles of attack(α=4oandα=6o). For all the cases a constant cavity length(x/c=0.5) is considered, and the cavitation numbers (s) are obtained from another potential based BEM (PCPAN) with the pressure recovery model. Furthermore, the influences of the different cavity termination models on the cavity shapes and pressure distributions of the cavitating foil are investigated. The fluctuations on pressure distribution is avo?ided in a similar method (filtering technique) given in(Longuet-Higgins and Cookelet, 1976). The 2D analyses results are compared with those of PCPAN (Kinnas and Fine,1993a) and a CFD technique (FLUENT).

    Cavitation prediction for the NACA 16006, NACA 16012 and NACA 16015 hydrofoils

    With the 2D BEM analysis a total number of 129 panels are implemented on the entire hydrofoil (Fig. 5). 129 panels are distributed in equal numbers on the pressure side [DM],on the cavity surface [DL] and on the non-cavitating part of the suction side [LM]. The panels are placed with full cosine spacing. The panel numbers are determined according to the best converged results obtained from the numerical studies.For all of the applications of the present method, the deformation amount non-dimensionalized by the chord length of the section on the cavity surface in each iteration step is Dh/c=0.00001 and the convergence rate is reached in a maximum 1 000-iterations in each case investigated here.The pressure recovery termination and the wall termination models implemented in the analyses are explained in the previous sections. With the present method, when a termination model is not selected, the wall termination model is applied by default. For comparison with similar conditions, the constants in the pressure recovery termination model are selected as the same as in PCPAN(n= 1,l= 0.1,A=0.5 as suggested for NACA 16006 section)(Kinnas and Fine, 1993a).

    Fig. 5 Panel arrangement of the 2D hydrofoil (height magnified by 5)

    In Figs. 6-17 the comparisons of the cavity shapes and pressure distributions by the present method and PCPAN for NACA 16006 (σ=0.896 and 1.349;α=4oand 6o), NACA 16012 (σ=0.878 and 1.330;α=4oand 6o) and NACA 16015(σ=0.863 and 1.311;α=4oand 6o) foil sections (height magnified by 5) are demonstrated. Furthermore, the results obtained from the present method with two different cavity termination models are compared.

    Fig. 6 Pressure distributions and cavity shapes (pressure recovery termination model) (2D hydrofoil, height magnified by 5)

    Fig. 7 Pressure distributions and cavity shapes from present method for the pressure recovery and wall termination models

    Fig. 8 Pressure distributions and cavity shapes (pressure recovery termination model)

    Fig. 9 Pressure distributions and cavity shapes from present method for the pressure recovery and wall termination models

    Fig. 10 Pressure distributions and cavity shapes (pressure recovery termination model)

    Fig. 11 Pressure distributions and cavity shapes from present method for the pressure recovery and wall termination models

    Fig. 12 Pressure distributions and cavity shapes (pressure recovery termination model)

    Fig. 13 Pressure distributions and cavity shapes from present method for the pressure recovery and wall termination models

    Fig. 14 Pressure distributions and cavity shapes (pressure recovery termination model)

    Fig. 15 Pressure distributions and cavity shapes from present method for the pressure recovery and wall termination models

    Fig. 16 Pressure distributions and cavity shapes (pressure recovery termination model)

    Fig. 17 Pressure distributions and cavity shapes from present method for the pressure recovery and wall termination models

    In Figs. 6-17 the non-dimensional cavity lengths (x/c) are obtained approximately as 0.5 from the present method with the pressure recovery model, and as 0.475 with the wall termination model.

    The results obtained for these cavity lengths from the present method are in good agreement with the results of PCPAN by means of the cavity shapes and pressure distributions for all cases (Figs. 6, 10, 12, 14, 16) except for the case of NACA 16006 atα=6o(Fig. 8). It is considered that the differences seen in the transition zone in Fig. 8 are due to the section with a small thickness and a high flow angle of attack. It is observed from Figs. 7, 9, 11, 13, 15, and 17 where the cavity closure models of the present method are compared,that the pressure distributions are in accordance except for the cavity closure region. In the comparison of the present method with the wall termination model, the occurrence of a small slip on the peak values of the pressure coefficient (-Cp)is seen due to the characteristics of the termination model in the transition zone.

    In the CFD analysis a NACA 16006 foil is used for the sheet cavity prediction in the 2D calculations. The CFD solver FLUENT 6.3 is used and the grid is generated with GAMBIT 2.4.6. For the 2D case a structured mesh with 46,800 elements is created. The near foil surface region is meshed denser for a correct cavitation simulation. The boundary conditions and the grid over the foil are shown in Figs. 18 and 19, respectively. The following boundary conditions are considered for the calculations, the inlet and the side surfaces are set as the velocity inlet and the outlet is set as the pressure-outlet. A first order upwind scheme is selected for the discretization. The mixture fluids for the cavitation modelling are selected as liquid water and water vapour. The density of the liquid water is selected as 1 000 kg/m3and the density of the water vapour is selected as 0.025 58 kg/m3. The non-condensable mass fraction is taken as 1.6×10-5. As for the turbulence model, thek-ε(RNG)turbulence model is set. The cavitation analyses are made by adjusting the vaporisation pressure, the flow velocity and the gauge pressure in the outlet to provide the intended cavitation number. Calculations for the two dimensional hydrofoil are made forα=4oandα=6owithσ=0.896,σ=1.35,and for the Reynolds numberRe= 7.9×106.

    Fig. 18 Mesh view and boundary conditions

    Fig. 19 Mesh around the hydrofoil

    In Figs 20 and 22, the pressure distributions for the 2D hydrofoil (NACA 16006) are shown forα=4o,σ=0.896 andα=6o,σ=1.349, respectively. The results from the present method with the pressure recovery termination model are in comparison with the results of CFD.

    Fig. 20 Pressure distributions from the present method and CFD

    Fig. 21 Pressure distributions from the present method and CFD

    Fig. 22 Pressure distributions from the present method and CFD

    Fig. 23 Pressure distributions from the present method and CFD

    Fig. 25 Pressure distributions from the present method and CFD

    In Figs. 20-25, the accordance of the pressure distributions (-Cp) with the present method and CFD is observed apart from the region of the transition zone where the closure model is used in the present method. So the comparison of the results of CFD including the viscous effects with the BEM results is reasonable except for near the region of the cavity termination. Further results are supported by the results presented in Krishnaswamy (2000)and Huanget al. (2010).

    As it can be seen in Fig. 21, especially near the region of the cavity termination, the pressure distributions obtained from the present method differ from the results of CFD due to the fact that the thinner section and the higher angles of attack of the incoming flow give thicker cavity shapes. This can also be seen in the case of the cavity shapes given in

    Figs. 26-30. If the angle of attack increases, the cavity becomes thicker and the difference between the cavity lengths of the present method and those of CFD becomes larger (Figs. 26 and 27). If the thickness of the foil increases,the difference now between the cavity lengths of the present method and those of CFD becomes smaller (Figs. 26 and 30).

    The cavity shape from the present method with the pressure recovery termination model and CFD are demonstrated in Figs. 26 and 31.

    From Figs. 22 and 23 it is observed that the cavity length obtained by CFD is shorter than the length from the present method, but the cavity shape is in accordance except for near the zone of the cavity termination. Furthermore, from Fig. 22 (α= 4o) it can be seen that the cavity length obtained by CFD is approximately equal to the length obtained by the present method in regard to the pressure distribution (-Cp)on the hydrofoil. It is assumed that the shorter cavity lengths obtained by CFD (which is more realistic because of the including of the viscous effects) occur due to the effects of the break off cycle. These effects are presented in de Lange and de Bruin (1998) based on the results of the experiments.

    For the NACA 16006, NACA 16012 and NACA 16015 sections and for the angles of attackα=4oandα=6o, the lift and drag coefficients (CL,CD) from the present method,PCPAN with the pressure recovery cavity closure model and CFD are given in Table 1.

    Fig. 26 Cavity shapes, present method (solid line) and CFD(colored), NACA 16006, s= 0.89 and α= 4o

    Fig. 27 Cavity shapes, present method (solid line) and CFD(colored), NACA 16006, s= 1.349 and α= 6o

    Fig. 28 Cavity shapes, present method (solid line) and CFD(colored), NACA 16012, s= 0.878 and α= 4o

    Fig. 29 Cavity shapes, present method (solid line) and CFD(colored), NACA 16012, s= 1.330 and α= 6o

    Fig. 30 Cavity shapes, present method (solid line) and CFD(colored), NACA 16015, s= 0.864 and α= 4o

    Fig. 31 Cavity shapes, present method (solid line) and CFD(colored), NACA 16015, s= 1.311 and α= 6o

    The 2D hydrofoil’s lift coefficient (CL) and drag coefficient (CD) are as follows:

    Table 1 The lift and drag coefficients (CL, CD) from different methods

    It is shown in Table 1 that the lift coefficients (CL) from this method show good agreement with those of PCPAN and CFD while the differences in the drag coefficients are relatively larger. TheCLvalues obtained from the present method show an approximately 2% deviation for all sections and angles in comparison to those of PCPAN. While the lift coefficients obtained from CFD for the NACA 16006 section differ about 3% from the PCPAN, the error increases for the section thickness (NACA 16012 and NACA 16015).It is assumed that these results are caused by the viscous effects which are included in the CFD analysis and by the decrease of the cavity thickness while the section thickness increases.

    6 Conclusions

    An iterative method based on the potential theory has been presented for the prediction of cavities around partially cavitating 2D hydrofoils. For a specified cavitation number and cavity length, the surface part of the hydrofoil in the cavity length is deformed iteratively until the dynamic boundary condition has been satisfied with a fully wetted analysis within an acceptable accuracy. The appropriate cavity length and shape of the 2D hydrofoil have been determined according to the minimum error criterion among different cavity lengths.

    For the 2D analysis, the NACA 16006, 16012 and NACA 16015 sections have been selected. The analyses have been performed with incoming flow angles of attack of 4 and 6 degrees. The results obtained from the 2D analysis have been compared by means of cavity shapes and pressure distributions with the results obtained by PCPAN and by the CFD technique. The results have shown good consistency with the results of PCPAN and CFD.

    As for the cavity closure model, the pressure recovery and wall termination models are used. From the applications carried out in the study, it has been found that the cavity lengths with the wall termination models are less than those with the pressure recovery model for all cases as expected.The other termination models can also be included with the present method.

    Since the present method developed in this study is based on a fully wetted analysis, the dynamic boundary condition over the cavity surface is satisfied automatically with the kinematic boundary condition. So the present method gives more realistic cavity shapes. The present method gives a faster approach for the cavity prediction than especially the CFD or experimental methods in terms of time.

    In future studies the methods presented here for the analysis of 2D cavitating hydrofoils can be extended to cavitating 3D hydrofoils and marine propellers. Furthermore,these methods can also be applied to supercavitating cases and to cases where face cavitation occurs.

    Acosta AJ (1955). A note on partial cavitation of flat plate hydrofoils. Tech. Rep. E-19.9, California Institute of Technology, Hydrodynamics Lab. Pasadena, USA.

    Ar?kan Y, ?elik F, Do?rul A, Bal ? (2012). Prediction of cavitation on two- and three-dimensional hydrofoils by an iterative BEM,

    Proceedings of the 8th International Symposium on Cavitation,CAV2012, Singapore, 696-702.

    Bal S, Kinnas SA (2003). A numerical wave tank model for cavitating hydrofoils.Computational Mechanics, 32(4-6),259-268.

    Bal S (2007). A numerical method for the prediction of wave pattern of surface piercing cavitating hydrofoils.Proceedings of the Institution of Mechanical Engineers, Part C,Journal of Mechanical Engineering Sciences, 221(12), 1623-1633.

    Bal S (2008). Performance prediction of surface piercing bodies in numerical towing tank.International Journal of Offshore and Polar Engineering, 18(2), 106-111.

    Bal S (2011). The effect of finite depth on 2-D and 3-D cavitating hydrofoils.Journal of Marine Science and Technology, 16(2),129-142.

    Berntsen GS, Kjeldsen M, Arndt REA (2001). Numerical modeling of sheet and tip vortex cavitation with fluent 5.Proceedings of the 4th International Symposium on Cavitation, Paris, France,CAV2001: Session, B5.006.

    Chow C (1979).An introduction to computational fluid mechanics.John WileyandSons, Inc, Colorado.

    Dang J, Kuiper G (1998). Re-entrant jet modeling of partial cavity flow on three-dimensional hydrofoils.Proceedings of ASME FEDSM’98, Washington DC, USA.

    Dang J (2001). Numerical simulation of partial cavity flows. PhD Thesis, Technical University Delft, Delft, Netherlands.

    de Lange D, de Bruin G (1998). Sheet cavitation and cloud cavitation, re-entrant jet and three-dimensionality.Applied Scientific Research, 58, 91-114.

    Deshpande MD, Feng J, Merkle CL (1993). Navier–stokes analysis of 2-d cavity flows.Proceedings of Cavitation and Multi-Phase Flow Forum,ASME FED, 153-158.

    Deshpande MD, Feng J, Merkle CL (1994). Cavity flow predictions based on the euler equation.Transactions of the ASME 116,36-44.

    Dupont P, Avellan F (1991). Numerical computation of a leading edge cavity.Proceedings of Cavitation’91, ASME FED,Portland, Oregon, USA, 116, 47 - 54.

    Fine NE, Kinnas SA (1993). A boundary element method for the analysis of the flow around 3-D cavitating hydrofoils.Journal of Ship Research, 37(3), 213-224.

    Fluent 6.3 User’s Guide (2006).

    Geurst JA, Timman R (1956). Linearized theory of two-dimensional cavitational flow around a wing section.IX Intl Cong. Applied Mech, Brussels 1, 486.

    Hess J, Smith A (1966). Calculation of potential flow about arbitrary bodies.Progress in Aeronautical Sciences, 3,138.

    Hess J (1973). Higher order numerical solution of the integral equation for the two-dimensional Neumann problem.Computer Methods in Applied Mechanics and Engineering, 2, 1-15.

    Hess J (1975). The use of higher-order surface singularity distributions to obtain improved potential flow solutions for two-dimensional lifting airfoils.Computer Methods in Applied Mechanics and Engineering, 5, 11-35.

    Hirschi R, Dupont P, Avellan F (1997). Centrifugal pump performance drop due to leading edge cavitation: Numerical predictions compared with model tests.FEDSM’97, Vancouver,Canada, 3342.

    Hoekstra M, Vaz G (2009). The partial cavity on a 2-d foil revisited.

    Proceedings of the 7th International Symposium on Cavitation,CAV2009, Grenoble, France, No:43.

    Huang S, He M, Wang C, Chang X (2010). Simulation of cavitating flow around a 2-D hydrofoil.Journal of Marine Science and Application, 9(1), 63-68.

    Katz J, Plotkin A (2001).Low-speed aerodynamics, 2nd edition.Cambridge University Press, Cambridge, UK.

    Kinnas S (1991). Leading-edge corrections to the linear theory of partially cavitating hydrofoils.Journal of Ship Research, 35(1),15-27.

    Kinnas SA, Fine NE (1991). Analysis of the flow around supercavitating hydrofoils with midchord and face cavity detachment.Journal of Ship Research, 35(3), 198-209.

    Kinnas SA, Fine NE (1992). A nonlinear boundary element method for the analysis of unsteady propeller sheet cavitation.

    Proceedings of 19th ONR, Seoul, Korea.

    Kinnas SA, Fine NE (1993). A numerical nonlinear analysis of the flow around two- and three- dimensional partial cavitating hydrofoils.Journal of Fluid Mechanics, 254, 151-181.

    Kinnas SA, Fine NE (1993a). MIT-PCPAN and MIT-SCPAN(Partially cavitating and super cavitating 2-D panel methods)User’s Manual, Version 1.0.

    Kinnas SA, Mishima S, Brewer WH (1996). Non-linear analysis of viscous flow around cavitating hydrofoils.Proceedings of 20th ONR, Santa Barbara, USA, 446-465.

    Kinnas SA (1998). The prediction of unsteady sheet cavitation.Proceedings of the Third International Symposium on Cavitation, Grenoble, France, 16-36.

    Kim YG, Lee CS, Suh JC (1994). Surface panel method for prediction of flow around a 3-D steady or unsteady cavitating hydrofoil.Proceedings of the 2nd International Symposium on Cavitation, Tokyo, Japan, 113-120.

    Kim YG, Lee CS (1996). Prediction of unsteady performance of marine propellers with cavitation using surface-panel method.Proceedings of 21st Symposium on Naval Hydrodynamics,Trondheim, Norway, 913-929.

    Krishnaswamy P (2000). Flow modeling for partially cavitating hydrofoils. PhD Thesis, Technical University of Denmark,Lyngby, Denmark.

    Li Ziru, Pourquie M, Terwisga TJC (2010). A numerical study of steady and unseady cavitation on a 2D hydrofoil.9th International Conference on Hydrodynamics, Shanghai, China,728-735.

    Longuet-Higgins MS, Cokelet ED (1976). The deformation of steep surface waves on water (I)—a numerical method of computation.Proceedings of Royal Society of London,Series A, 350, 1-26.

    Singhal AK, Li HY, Athavale MM, Jiang Y (2001). Mathematical basis and validation of the full cavitation model.ASME FEDSM'01, New Orleans, Louisiana, USA.

    Tulin MP (1953). Steady two-dimensional cavity flows about slender bodies. DTMB Technical Report N. 834.

    Uhlman JJS (1987). The surface singularity method applied to partially cavitating hydrofoils.Journal of Ship Research, 31(2),107-124.

    Author biographies

    Fahri ?elik, graduated from the Department of Naval Architecture and Marine Engineering, ?stanbul Technical University, Turkey in 1994. He completed his M.Sc., PhD in 1997, 2005, respectively. He has been working as a Professor at Yildiz Technical University,Turkey.

    Yasemin Arikan ?zden, graduated from the Department of Naval Architecture and Marine Engineering, Yildiz Technical University, Turkey in 2008. Shecompleted her M.Sc in 2010 and has been studying for her PhD since 2010 at Yildiz Technical University, Turkey. She has also been working as a Research Assistant at the NAME of YTU.

    Sakir Bal was born in 1967. He is a full professor at the Department of Naval Architecture and Marine Engineering and superintendent of the Ata Nutku Ship Model Testing Laboratory,Istanbul Technical University. He was a visiting research scholar at the University of Texas, Austin, USA and the University of Newcastle upon Tyne (UK). He gave lectures at the University of Liege, Belgium. His current research interests include ship hydrodynamics, and marine propellers.

    国产在视频线在精品| 久久久久久久久中文| 欧美极品一区二区三区四区| 久久久久性生活片| 国产精品一区www在线观看 | 欧美丝袜亚洲另类 | 看十八女毛片水多多多| 内射极品少妇av片p| 可以在线观看毛片的网站| 中文字幕高清在线视频| 国产激情偷乱视频一区二区| 免费观看精品视频网站| 国产在线精品亚洲第一网站| 乱码一卡2卡4卡精品| 九九爱精品视频在线观看| 成人二区视频| 色视频www国产| 精品欧美国产一区二区三| 亚洲欧美清纯卡通| 日本一二三区视频观看| 久久国内精品自在自线图片| 嫁个100分男人电影在线观看| 一本一本综合久久| 日韩高清综合在线| 男人的好看免费观看在线视频| 我要搜黄色片| 国产精品女同一区二区软件 | 啦啦啦韩国在线观看视频| 国产又黄又爽又无遮挡在线| 中文字幕免费在线视频6| 欧美3d第一页| 欧美3d第一页| 亚洲avbb在线观看| 免费黄网站久久成人精品| 色尼玛亚洲综合影院| 欧美另类亚洲清纯唯美| 精品午夜福利视频在线观看一区| 免费看日本二区| 少妇人妻精品综合一区二区 | 国产精品不卡视频一区二区| 伦精品一区二区三区| 色5月婷婷丁香| 在线国产一区二区在线| av专区在线播放| 国产在线男女| 91在线观看av| 国产激情偷乱视频一区二区| 成人一区二区视频在线观看| 联通29元200g的流量卡| 尤物成人国产欧美一区二区三区| 又爽又黄a免费视频| 亚洲天堂国产精品一区在线| 欧美高清性xxxxhd video| 精品日产1卡2卡| 久久久久精品国产欧美久久久| 日本三级黄在线观看| av在线老鸭窝| 大又大粗又爽又黄少妇毛片口| 亚洲av成人精品一区久久| 久久久久久久久久久丰满 | 舔av片在线| 在线观看av片永久免费下载| 精品人妻1区二区| 日本五十路高清| 美女大奶头视频| 色在线成人网| 色噜噜av男人的天堂激情| 高清日韩中文字幕在线| 国产精品免费一区二区三区在线| 国产精品一区二区三区四区久久| 免费高清视频大片| 亚洲熟妇熟女久久| 亚洲精品国产成人久久av| 亚洲熟妇熟女久久| a级毛片a级免费在线| 国产精品爽爽va在线观看网站| 听说在线观看完整版免费高清| 亚洲无线在线观看| 国产精品女同一区二区软件 | 色5月婷婷丁香| 又黄又爽又免费观看的视频| 亚洲性久久影院| 欧美激情久久久久久爽电影| 日本成人三级电影网站| 日本-黄色视频高清免费观看| 国内少妇人妻偷人精品xxx网站| 久久精品人妻少妇| 日本五十路高清| 五月伊人婷婷丁香| 国产精品一及| 麻豆国产97在线/欧美| 亚洲国产欧美人成| 亚洲精品国产成人久久av| 一个人免费在线观看电影| av视频在线观看入口| 一区二区三区免费毛片| 最后的刺客免费高清国语| 最近在线观看免费完整版| 欧美色欧美亚洲另类二区| 亚洲久久久久久中文字幕| 国产色婷婷99| 日日摸夜夜添夜夜添av毛片 | 精品免费久久久久久久清纯| 久久午夜福利片| 国产淫片久久久久久久久| 成人鲁丝片一二三区免费| 天堂影院成人在线观看| 精品一区二区三区av网在线观看| 91在线观看av| 午夜福利高清视频| eeuss影院久久| 日本免费a在线| 热99re8久久精品国产| 国产又黄又爽又无遮挡在线| 日韩在线高清观看一区二区三区 | 91午夜精品亚洲一区二区三区 | 亚洲精品456在线播放app | 97超视频在线观看视频| 国产一区二区亚洲精品在线观看| 日本a在线网址| 日本一二三区视频观看| 日本撒尿小便嘘嘘汇集6| 内地一区二区视频在线| 日韩欧美在线乱码| 免费av不卡在线播放| 国产精品福利在线免费观看| 亚洲av美国av| 精品免费久久久久久久清纯| 在线免费十八禁| 欧美另类亚洲清纯唯美| 九色国产91popny在线| 久久久久久久久中文| 婷婷亚洲欧美| 人妻夜夜爽99麻豆av| 五月伊人婷婷丁香| 国产精品日韩av在线免费观看| 最近中文字幕高清免费大全6 | 国产高清三级在线| 成人永久免费在线观看视频| 成人特级黄色片久久久久久久| 午夜福利18| 欧美精品国产亚洲| 最近最新免费中文字幕在线| 69av精品久久久久久| 成人性生交大片免费视频hd| 久久精品国产亚洲av天美| 看黄色毛片网站| 国产视频一区二区在线看| 久久欧美精品欧美久久欧美| 精品一区二区免费观看| 国产精品自产拍在线观看55亚洲| 亚洲黑人精品在线| 18禁黄网站禁片午夜丰满| 国产 一区 欧美 日韩| 欧美绝顶高潮抽搐喷水| 国产视频一区二区在线看| 色在线成人网| 国产男人的电影天堂91| 久久亚洲真实| 在线看三级毛片| 久久国产精品人妻蜜桃| 老司机深夜福利视频在线观看| 国产精品一区二区性色av| 少妇被粗大猛烈的视频| 亚洲最大成人av| 人妻少妇偷人精品九色| 男女那种视频在线观看| 一个人看视频在线观看www免费| 搡老岳熟女国产| 99国产精品一区二区蜜桃av| 成熟少妇高潮喷水视频| 国产毛片a区久久久久| 欧美日韩精品成人综合77777| 久久久久久久精品吃奶| 色av中文字幕| 国产精品1区2区在线观看.| 99久久中文字幕三级久久日本| 99久久精品国产国产毛片| 精品久久久久久久久久免费视频| 亚洲一区二区三区色噜噜| 少妇的逼好多水| 精品无人区乱码1区二区| 国产又黄又爽又无遮挡在线| 亚洲av五月六月丁香网| 午夜福利18| 国产一区二区激情短视频| 欧美在线一区亚洲| 国产成人av教育| 亚洲在线自拍视频| 看免费成人av毛片| 精品人妻1区二区| 有码 亚洲区| 国产亚洲av嫩草精品影院| 女的被弄到高潮叫床怎么办 | 久久精品国产清高在天天线| 欧美一区二区亚洲| 久久久国产成人免费| 国产久久久一区二区三区| 搞女人的毛片| 国内少妇人妻偷人精品xxx网站| 色吧在线观看| 亚洲黑人精品在线| 亚洲精品456在线播放app | 亚洲18禁久久av| 一个人观看的视频www高清免费观看| 成年女人毛片免费观看观看9| 午夜免费男女啪啪视频观看 | 一进一出抽搐gif免费好疼| 久久久成人免费电影| 国产精品女同一区二区软件 | 久久久午夜欧美精品| 国产成人福利小说| 春色校园在线视频观看| 蜜桃亚洲精品一区二区三区| av在线亚洲专区| 中文资源天堂在线| 欧美日本视频| 啪啪无遮挡十八禁网站| 99久久久亚洲精品蜜臀av| 国产亚洲精品久久久久久毛片| 别揉我奶头~嗯~啊~动态视频| 国产淫片久久久久久久久| 一卡2卡三卡四卡精品乱码亚洲| 亚洲va日本ⅴa欧美va伊人久久| 国产 一区精品| 超碰av人人做人人爽久久| 成人二区视频| 少妇熟女aⅴ在线视频| 国产精品无大码| 亚洲av第一区精品v没综合| 欧美另类亚洲清纯唯美| 免费看光身美女| 91av网一区二区| 亚洲av第一区精品v没综合| 亚洲欧美精品综合久久99| 亚洲欧美日韩东京热| 久久精品国产自在天天线| videossex国产| 国产综合懂色| 亚洲一区二区三区色噜噜| 天堂动漫精品| 嫁个100分男人电影在线观看| 欧美一级a爱片免费观看看| 久久久久久九九精品二区国产| 国产国拍精品亚洲av在线观看| 麻豆成人午夜福利视频| 天堂av国产一区二区熟女人妻| av天堂在线播放| 嫩草影院新地址| 精品日产1卡2卡| 国产精品98久久久久久宅男小说| 黄色女人牲交| www.www免费av| 国产 一区 欧美 日韩| 欧美日韩黄片免| 亚洲18禁久久av| 真人一进一出gif抽搐免费| 国产色婷婷99| 日韩一本色道免费dvd| 露出奶头的视频| 乱系列少妇在线播放| 人妻丰满熟妇av一区二区三区| 精品一区二区三区人妻视频| 国产精品国产三级国产av玫瑰| 给我免费播放毛片高清在线观看| 成人美女网站在线观看视频| 悠悠久久av| 久久精品综合一区二区三区| 国产麻豆成人av免费视频| 国产伦精品一区二区三区四那| 国产乱人视频| 观看美女的网站| 成人二区视频| 日韩,欧美,国产一区二区三区 | 亚洲av日韩精品久久久久久密| 国内精品久久久久久久电影| 韩国av一区二区三区四区| 成人三级黄色视频| 欧美在线一区亚洲| 国产在视频线在精品| av.在线天堂| 国产国拍精品亚洲av在线观看| 精品无人区乱码1区二区| 国产伦在线观看视频一区| 国产三级在线视频| 天美传媒精品一区二区| 色av中文字幕| 又紧又爽又黄一区二区| 哪里可以看免费的av片| 好男人在线观看高清免费视频| 日韩欧美在线乱码| 日日啪夜夜撸| 99国产精品一区二区蜜桃av| 最近最新中文字幕大全电影3| 日本-黄色视频高清免费观看| 九色成人免费人妻av| 干丝袜人妻中文字幕| 久久久久国产精品人妻aⅴ院| 男女做爰动态图高潮gif福利片| 嫁个100分男人电影在线观看| 久久久久久久亚洲中文字幕| 国产极品精品免费视频能看的| 色综合站精品国产| 成年免费大片在线观看| 欧美激情久久久久久爽电影| 久久亚洲精品不卡| 性色avwww在线观看| 国内精品久久久久精免费| 91久久精品电影网| 欧美一区二区国产精品久久精品| 午夜激情欧美在线| av在线亚洲专区| 18禁黄网站禁片免费观看直播| 欧美在线一区亚洲| 男女那种视频在线观看| 中文字幕av成人在线电影| 俄罗斯特黄特色一大片| 亚洲av不卡在线观看| 国产高清不卡午夜福利| 亚洲国产欧洲综合997久久,| bbb黄色大片| 赤兔流量卡办理| 亚洲精品日韩av片在线观看| 亚洲黑人精品在线| 国产综合懂色| 免费在线观看影片大全网站| 欧美性感艳星| 精品日产1卡2卡| 动漫黄色视频在线观看| 琪琪午夜伦伦电影理论片6080| 国产亚洲欧美98| 亚洲专区国产一区二区| 变态另类丝袜制服| 婷婷丁香在线五月| av.在线天堂| 成人特级av手机在线观看| 桃色一区二区三区在线观看| 少妇高潮的动态图| 国产精品国产三级国产av玫瑰| 国产成人aa在线观看| 久久久久免费精品人妻一区二区| 成人三级黄色视频| 免费搜索国产男女视频| 午夜久久久久精精品| 日本-黄色视频高清免费观看| 欧美性猛交黑人性爽| 精品日产1卡2卡| 国产女主播在线喷水免费视频网站 | 舔av片在线| 国产久久久一区二区三区| 三级国产精品欧美在线观看| a级毛片免费高清观看在线播放| 欧美丝袜亚洲另类 | 天美传媒精品一区二区| 直男gayav资源| 成人毛片a级毛片在线播放| 黄色视频,在线免费观看| 国产久久久一区二区三区| 国产极品精品免费视频能看的| 亚洲狠狠婷婷综合久久图片| 内地一区二区视频在线| 两个人的视频大全免费| 国产黄色小视频在线观看| 日本在线视频免费播放| 午夜福利高清视频| 日日干狠狠操夜夜爽| 真实男女啪啪啪动态图| 成人国产一区最新在线观看| 国产高清有码在线观看视频| 国产高清三级在线| 亚洲av中文av极速乱 | 国产精品乱码一区二三区的特点| 99久久精品热视频| 精品人妻熟女av久视频| 久久久久久久久久成人| 国产高清激情床上av| 久99久视频精品免费| 欧美不卡视频在线免费观看| 国产av不卡久久| 内地一区二区视频在线| 色哟哟哟哟哟哟| 最近中文字幕高清免费大全6 | 亚洲人成网站在线播放欧美日韩| 欧美一区二区国产精品久久精品| 99热这里只有精品一区| av在线观看视频网站免费| 男人的好看免费观看在线视频| 欧美+日韩+精品| 色吧在线观看| 美女高潮的动态| 国产午夜福利久久久久久| 俄罗斯特黄特色一大片| 级片在线观看| 亚洲av中文av极速乱 | 看十八女毛片水多多多| 天堂动漫精品| 成年女人毛片免费观看观看9| 全区人妻精品视频| 色综合站精品国产| 九九爱精品视频在线观看| 国产爱豆传媒在线观看| 日本免费a在线| 久久久精品大字幕| 精华霜和精华液先用哪个| 在现免费观看毛片| 欧美精品啪啪一区二区三区| 欧美精品国产亚洲| 国产单亲对白刺激| 永久网站在线| 狠狠狠狠99中文字幕| 大又大粗又爽又黄少妇毛片口| 婷婷色综合大香蕉| 成人av在线播放网站| 婷婷亚洲欧美| 色综合站精品国产| av黄色大香蕉| 亚洲成人精品中文字幕电影| 国模一区二区三区四区视频| 伦精品一区二区三区| АⅤ资源中文在线天堂| 国产精品久久久久久久久免| 身体一侧抽搐| 99国产精品一区二区蜜桃av| 天堂av国产一区二区熟女人妻| 国产精品一及| 久久久久性生活片| 久久精品91蜜桃| 校园人妻丝袜中文字幕| 免费人成在线观看视频色| 国产不卡一卡二| 无人区码免费观看不卡| 五月玫瑰六月丁香| 欧美中文日本在线观看视频| 国产免费一级a男人的天堂| 欧美黑人巨大hd| 18+在线观看网站| 美女黄网站色视频| 国产成人aa在线观看| av福利片在线观看| 亚洲av一区综合| 99久久精品一区二区三区| 久久这里只有精品中国| 国产中年淑女户外野战色| 全区人妻精品视频| 欧美激情在线99| 亚洲午夜理论影院| 日本精品一区二区三区蜜桃| 老司机福利观看| 免费看a级黄色片| 热99re8久久精品国产| 12—13女人毛片做爰片一| 久久久色成人| 一区二区三区激情视频| 最近中文字幕高清免费大全6 | 日本黄大片高清| 男女视频在线观看网站免费| 色5月婷婷丁香| 成年女人永久免费观看视频| 级片在线观看| 亚洲成人中文字幕在线播放| 少妇人妻精品综合一区二区 | 琪琪午夜伦伦电影理论片6080| 乱码一卡2卡4卡精品| 一a级毛片在线观看| 久久99热6这里只有精品| 99久久久亚洲精品蜜臀av| 亚洲av.av天堂| 国产精品福利在线免费观看| 久久久色成人| 在线a可以看的网站| 淫秽高清视频在线观看| 中文字幕久久专区| 国产精品免费一区二区三区在线| 欧美极品一区二区三区四区| 亚洲自拍偷在线| 精品无人区乱码1区二区| 欧美一区二区精品小视频在线| 好男人在线观看高清免费视频| 动漫黄色视频在线观看| 国产亚洲精品久久久com| 亚洲精品一卡2卡三卡4卡5卡| 亚洲欧美日韩卡通动漫| 日日摸夜夜添夜夜添av毛片 | 日韩一区二区视频免费看| 五月玫瑰六月丁香| 黄色丝袜av网址大全| 少妇人妻一区二区三区视频| 在现免费观看毛片| 波野结衣二区三区在线| 日本免费a在线| 久久久久久久久中文| АⅤ资源中文在线天堂| 亚洲美女搞黄在线观看 | 99久久久亚洲精品蜜臀av| 久久天躁狠狠躁夜夜2o2o| 两性午夜刺激爽爽歪歪视频在线观看| 欧美日韩中文字幕国产精品一区二区三区| 又爽又黄a免费视频| 最新中文字幕久久久久| 久久99热6这里只有精品| 老熟妇仑乱视频hdxx| 午夜福利视频1000在线观看| 国产乱人伦免费视频| а√天堂www在线а√下载| 国产精品人妻久久久久久| 欧美精品国产亚洲| 国产色爽女视频免费观看| 丰满的人妻完整版| 一卡2卡三卡四卡精品乱码亚洲| 国产真实伦视频高清在线观看 | 成人毛片a级毛片在线播放| 亚洲va在线va天堂va国产| 亚洲国产精品合色在线| 最近视频中文字幕2019在线8| 中文字幕久久专区| 99久久中文字幕三级久久日本| 久久精品国产鲁丝片午夜精品 | 精品无人区乱码1区二区| 国产乱人视频| 国产人妻一区二区三区在| 亚洲第一区二区三区不卡| 少妇猛男粗大的猛烈进出视频 | 美女免费视频网站| 国产日本99.免费观看| 日韩欧美国产一区二区入口| 人妻少妇偷人精品九色| 国产高清视频在线播放一区| 真人一进一出gif抽搐免费| 三级国产精品欧美在线观看| 日本欧美国产在线视频| 欧美三级亚洲精品| 日日干狠狠操夜夜爽| 国产伦一二天堂av在线观看| av在线天堂中文字幕| 国产精品不卡视频一区二区| 欧美一区二区亚洲| 一个人看视频在线观看www免费| 91在线精品国自产拍蜜月| 黄色配什么色好看| 天美传媒精品一区二区| 日韩国内少妇激情av| 一区二区三区四区激情视频 | 俺也久久电影网| 国产视频一区二区在线看| 桃红色精品国产亚洲av| 夜夜看夜夜爽夜夜摸| 国产精品久久视频播放| 在线观看66精品国产| 嫩草影院入口| 免费av不卡在线播放| www.www免费av| 日本-黄色视频高清免费观看| 亚洲欧美清纯卡通| 日本黄大片高清| 中文字幕精品亚洲无线码一区| 人人妻,人人澡人人爽秒播| 两人在一起打扑克的视频| 欧美日本亚洲视频在线播放| 久久久国产成人精品二区| 无遮挡黄片免费观看| 成人亚洲精品av一区二区| 日本 av在线| 精品久久久久久久久久久久久| 精品乱码久久久久久99久播| 欧美一级a爱片免费观看看| 欧美黑人欧美精品刺激| 欧美潮喷喷水| 91久久精品国产一区二区三区| 美女高潮喷水抽搐中文字幕| 99热这里只有精品一区| 18禁黄网站禁片免费观看直播| 又爽又黄a免费视频| 天堂√8在线中文| 18+在线观看网站| 制服丝袜大香蕉在线| 特级一级黄色大片| 女人十人毛片免费观看3o分钟| 国产白丝娇喘喷水9色精品| 99久久精品国产国产毛片| 99久久久亚洲精品蜜臀av| 国产中年淑女户外野战色| 亚洲自拍偷在线| 亚洲美女搞黄在线观看 | 亚洲综合色惰| 国内精品久久久久精免费| 国产一区二区在线av高清观看| 国产在线男女| 超碰av人人做人人爽久久| 亚洲av五月六月丁香网| 久久久久九九精品影院| 麻豆久久精品国产亚洲av| 欧美xxxx性猛交bbbb| 精品人妻1区二区| 国产v大片淫在线免费观看| 99久久精品热视频| 亚洲av成人精品一区久久| 国产一区二区亚洲精品在线观看| 在线观看舔阴道视频| 午夜亚洲福利在线播放| 一个人看的www免费观看视频| 男插女下体视频免费在线播放| ponron亚洲| 成人鲁丝片一二三区免费| 麻豆精品久久久久久蜜桃| 少妇的逼水好多| 舔av片在线| 国产精品99久久久久久久久| 久久午夜福利片| 一区二区三区高清视频在线| 美女 人体艺术 gogo| 五月伊人婷婷丁香| 国产大屁股一区二区在线视频| 亚洲18禁久久av| 在线观看舔阴道视频| 天天躁日日操中文字幕|