亚洲免费av电影一区二区三区,日韩爱爱视频,51精品视频一区二区三区,91视频爱爱,日韩欧美在线播放视频,中文字幕少妇AV,亚洲电影中文字幕,久久久久亚洲av成人网址,久久综合视频网站,国产在线不卡免费播放

        ?

        Prediction of Ship-Ship Interaction Forces in Shallow Water Using a High-order Panel Method

        2016-05-16 02:41:50XUHufuZOUZojinLIUXioyn
        船舶力學(xué) 2016年12期
        關(guān)鍵詞:會(huì)遇上海交通大學(xué)元法

        XU Hu-fu,ZOU Zo-jin,LIU Xio-yn

        (a.School of Naval Architecture,Ocean and Civil Engineering;b.State Key Laboratory of Ocean Engineering,Shanghai Jiao Tong University,Shanghai 200240,China)

        Prediction of Ship-Ship Interaction Forces in Shallow Water Using a High-order Panel Method

        XU Hua-fua,b,ZOU Zao-jiana,b,LIU Xiao-yana

        (a.School of Naval Architecture,Ocean and Civil Engineering;b.State Key Laboratory of Ocean Engineering,Shanghai Jiao Tong University,Shanghai 200240,China)

        A three-dimensional high order panel method based on Non-Uniform Rational B-Spline (NURBS)is developed for predicting the ship-ship hydrodynamic interaction forces during meeting and overtaking in shallow water.A NURBS surface is used to precisely represent the body geometry. Velocity potential on the body surface is described by B-spline after the source density distribution on the boundary surface is determined.A collocation approach is applied to the boundary integral equation discretization.Under the assumption of low passing speed,the effect of free surface elevation is neglected in the numerical calculation,and the infinite image method is used to deal with the finite water depth effect.The time stepping method is used to solve the velocity potential at each time step.The present results of hydrodynamic interaction forces are compared with those obtained by slender body theory and CFD-based RANS method to show the validity of the proposed numerical method.Calculations are conducted for different water depths,lateral distances between ships and ships’speeds,and the detailed results are presented to demonstrate the effects of these factors.The results show that interaction forces are much larger and much more susceptible to water depth,but less susceptible to lateral distance in meeting case than those in overtaking case.

        meeting;overtaking;ship-ship hydrodynamic interaction;NURBS; high-order panel method

        0 Introduction

        Ship-to-ship hydrodynamic interactions involving significant hydrodynamic forces and moments will occur when two ships are moving in close proximity.Particular cases of special interest are the unsteady-state situation of two ships moving in tandem in re-fueling or lighting operation,and situations of two ships meeting or overtaking in restricted water.In each of these cases,hydrodynamic interaction forces and moments may lead to accidents such as collision.

        Many previous studies are based on experimental method since it typically produces reliable and realistic results.Vantorre et al(2002)[1]carried out an extensive experimental study onship-to-ship interaction in shallow water involving various ship types.

        Theoretical study of ship-ship hydrodynamic interaction was traditionally based on the slender-body theory.Tuck and Newman(1974)[2]extended the slender-body theory to deal with the hydrodynamic interactions between two ships moving with different speed in deep water or with same speed in shallow water.Yeung(1978)[3]undertook theoretical studies on ship-ship hydrodynamic interaction by using the slender-body theory and the method of matched asymptotic expansions.A similar approach was used by Yeung and Tan(1980)[4]to study the hydrodynamic interaction of ships with some fixed obstacles.This study was extended by Hsiung and Gui(1988)[5]to a larger variety of obstacles.In all the studies mentioned above,the free surface effects were ignored under the assumption of low ship speed.

        During the last decades,there has also been a growth in numerical study of ship-ship hydrodynamic interaction.Korsmeyer et al(1993)[6]studied the body-to-body interaction in a rectangular canal by using a rectangular Green function.Parallel motion was studied for two spheroids,two Mariner ships and two Panamax bulk-carriers,and comparisons with experimental data showed a fairly good agreement.Yasukawa(2003)[7]derived motion equations for two ships maneuvering in close proximity within potential theory,and simulated the motion of two ferries when one overtakes the other and passes away.A 3D panel method was applied to calculate the added mass and the interaction forces,and the collision caused by the hydrodynamic interaction between ships in shallow water was demonstrated.Xiang and Faltinsen (2010)[8],Zhou et al(2012)[9]carried out the calculation of hydrodynamic interaction forces in restricted waters using potential theory,both encountering and overtaking between target ship and own ship were considered.Sutulo et al(2012)[10]made validation of potential-flow estimation of interaction forces between ships,and affecting factors were analyzed in detail.In all these numerical studies,the free surface effects were ignored.Moreover,all these numerical researches use constant panel method,which is difficult to satisfy the C1-continuity on the body surface.Thus,a high-order panel method is needed.

        During the last decades,high-order panel method based on B-spline or Non-Uniform Rational B-spline(NURBS)has been developed rapidly.B-spline and NURBS can provide a more precise description of the body geometry of the ship surface for the purpose of hydrodynamic calculation;moreover,they can be used to represent the source density or velocity potential distribution for potential flow problem,so that the continuity of higher order derivatives of velocity potential on the ship surface can be ensured precisely.Zhang et al(2003)[11]developed a numerical computation method based on B-splines for the hydrodynamic interaction forces between two ships,where the free surface effects were ignored.

        In the present paper,a three-dimensional high-order panel method based on NURBS is developed to calculate the ship-to-ship interaction force both in meeting and overtaking conditions.A NURBS surface is used to precisely represent the body geometry.The velocity potential distribution on the body surface is described by B-spline after the source density distribution on the body surface is determined.The collocation method is adopted and the collo-cation points are distributed on the vertices of NURBS surface,and the singularity sources are distributed on the Gaussian points of each panel on the NURBS surface.Under the assumption of low ship speed,the effect of free surface elevation is ignored,and the infinite image method is used to deal with the finite water depth effect.Calculations are carried out for different water depths,lateral distances between ships and ships’speeds to analyze the influences of these factors.

        1 Mathematical formulation

        We consider two ships moving along parallel courses in proximity of each other,one ship (Ship 1)navigates with a constant speed U1,and another ship(Ship 2)navigates with a constant speed U2,as shown in Fig.1.The coordinate systems o-xyz fixed in the space and oi-xiyizifixed to Ship i(i=1,2)are adopted,where the xi-axis is pointing from ship stern to the bow, yi-axis to port side and zi-axis vertically upward;the oi-xiyiplane coincides with the undisturbed free surface(z=0).The water depth is assumed to be constant and expressed by z=-h. The longitudinal and transverse distances between the ships are denoted by ST and Sp,respectively.

        Fig.1 Sketch of ship-to-ship interaction

        It is assumed that(1)The fluid is inviscid and the flow is irrotational;(2)The ship speeds are very small,so that the effects of free-surface elevation can be ignored.Furthermore,the shedding of vortices due to viscous effect and flow separation is neglected.The neglect of free-surface effects and shed vortices removes all the memory effects and allows the solution to be stepped through time as a series of independent hydrodynamic calculations(Korsmeyer,et al, 1993)[6].

        The perturbation velocity potential representing ship flow is defined as φ( t,x,y,z),it satisfies the Laplace equation in the fluid domain and the following boundary conditions:

        The velocity potential φ is expressed by source distribution on the hull surfaces as follows:

        where σ(i)denotes the strength of source distributed on the hull surface of Ship i.P( x,y,z)is a field point andis a singular point on the surface S which consists of the wetted surfaces S1and S2.The Green function takes the following form:

        With application of this Green function formula of infinite mirror-image,the boundary conditions(3)and(4)on the undisturbed free surface and the sea bottom are satisfied automatically;while by satisfying the boundary conditions(1)and(2)on the hull surfaces,it follows the following Fredholm integral equation of second kind:

        After Eq.(7)is solved,the disturbance velocity potential and disturbance velocity can be expressed through the density σ by Eq.(5)and Eq.(8):

        The pressure distribution on Ship i is determined from the unsteady Bernoulli equationin the following form:

        The hydrodynamic force and moment can be also represented in the standard component form.The lateral force Yi=Fi2and yaw moment Ni=Mi3on Ship i can be expressed as:

        2 Solution procedure of the panel method based on NURBS

        In present study,not only the body geometry,but also the source density is described by NURBS.The expressions of an arbitrary point P and the distributed source density σ on the body surface at a given moment take the forms:

        where u,v are the parameters representing two directions of the body surface and u,v∈[0,1]; Dijis the control point of the body surface,Sijis the control point of the distributed source at a given moment;wijis the weight in relation to the control points;m,n are the numbers of the control points in the u,v directions,respectively;Ni,k(u),Nj,l(v)are the B-spline basis functions;k,l are the degrees of B-spline basis functions.In the present study,the degree of B-spline basis functions is chosen as 3,which can provide a C2-continuous representation of the body surface.

        Substituting Eq.(13)into Eq.(7),it follows

        where E=Pu·Pu,F=Pu·Pv,G=Pv·Pv.

        The discrete form of the second kind of Fredholm integral equation is satisfied on the collocation points P( u0,v0),and the source is distributed on the Gauss points of each panel on the body surfaces.

        From Eq.(14),the source density control points Sijcan be obtained;then by substituting Sijinto Eq.(13),the source density on the body surfaces is obtained,and from Eq.(5),the velocity potential of arbitrary point in the fluid field can be determined.Finally,the pressure and hydrodynamic interaction forces can be calculated from Eqs.(9),(10)and(11).

        3 Numerical results and discussions

        Taking two Wigley hulls of same size as an example,calculations are conducted.The Wigley hull is a mathematical ship hull,whose expression is as follows:

        where L,B and D are the ship length,breadth and draught,respectively.The main parameters of the Wigley hull are given in Tab.1.

        ξ=2.0ST/(L1+L2),the normalized longitudinal distance between ships,is used to replace t in the time stepping method,where Liis the length of Ship i.Considering the reality that the interaction hydrodynamic forces vanish to zero when the longitudinal distance between ships is larger than two times of ship length,the calculation starts at an initial position of ξ=-2 and ends at the position of ξ=2.In order to satisfy the assumption of low ship speed,the speeds of the ship models in the meeting case are U1=U2=0.358 9 m/s,corresponding to Fn=0.066 2;while in the overtaking case,the speeds of the ship models are U1=0.358 9 m/s and U2=0.538 4 m/s,the corresponding Froude number is Fn= 0.066 2 and 0.10,respectively.

        The numerical results of hydrodynamic lateral force and yaw moment are normalized as:

        Tab.1 Main parameters of the Wigley hull

        According to the convergence study carried out by Xu and Zou[12]under the condition of water depth to draught ratio h/D=1.15 and the lateral distance between ships Sp=1.5B,5×16 panels on each of the ship hull and time step of 0.261 2 s are appropriate.Also,for the infinite image Green function truncated terms,a=±6 is sufficient.

        3.1 Comparison of the results

        The present numerical results of the hydrodynamic force and moment on Ship 1 are compared with those of the slender-body theory(Tuck and Newman[2],Yeung[3],Varyani and Krishnankutty[13]),as shown in Fig.2 and Fig.3.In these figures,the numerical results obtained by the authors using the CFD-based RANS method are also compared.

        Fig.2 Comparison of lateral force and yaw moment on Ship 1 in meeting condition, h/D=1.3,Sp=1.25 m(4.17B),U1=0.358 9 m/s,U2=0.358 9 m/s

        Fig.3 Comparison of lateral force and yaw moment on Ship 1 in overtaking condition, h/D=2.0,Sp=2.5 B,U1=0.358 9 m/s,U2=0.538 4 m/s

        It is obvious that the interaction forces and moments emerge at the distance between two ships within 2 times of ship length and vanish to zero as the distance is becoming larger and larger both in meeting and overtaking conditions,and the interaction forces and moments are much larger in meeting case than those in overtaking case.For both conditions,the lateral forces are characterized by initial repulsion,followed by attraction and finally repulsion again; and the first repulsion peak occurs when the longitudinal distance between the bows of the two ships is zero for meeting condition or the longitudinal distance between the bow of the overtaking ship(Ship 2)and the stern of the overtaken ship(Ship 1)is zero for overtaking con-dition.The peak values of attraction occur when the two ships are abeam.Just like the first repulsion peak,the second repulsion peak occurs when the longitudinal distance between the stern of Ship 1 and the stern of Ship 2 is zero in meeting case or the longitudinal distance between the stern of the overtaking ship(Ship 2)and the bow of the overtaken ship(Ship 1)is zero in overtaking case.The change tendencies of the yaw moment are rather different from those of lateral forces.There are four phases can be distinguished for meeting condition,which are characterized by consecutive actions of bow repulsion,bow attraction,bow repulsion and bow attraction,while the yaw moment shows only two peaks,that is,consecutive bow repulsion and bow attraction for overtaking condition.In meeting case,the peak value of yaw moment occurs not only when the longitudinal distance between two ships equal to one ship length,the same position when lateral forces extremes occur,it also occurs when the distance between two ships equals to half ship length in overtaking case.

        From these figures it can be seen that the present numerical method can give a qualitatively correct prediction of the change tendency of ship-ship hydrodynamic interaction forces, and the agreements between the numerical results are better in the overtaking case than in the meeting case,which means that the viscous effects do not play an important role in the overtaking case.

        3.2 Influence of water depth

        To investigate the influence of water depth on the hydrodynamic interaction forces,calculations are conducted for different water depth to draught ratio h/D=1.15,1.30,1.50,1.80 and 2.00,while the lateral distance between ships Sp=1.5B keeps unchanged.Fig.4 and Fig.5 show the results in meeting and overtaking conditions,respectively.

        It can be seen that the change tendency is same for different water depth,and the magnitudes of the hydrodynamic interaction force and moment decrease with the increase of water depth in both meeting and overtaking conditions.The peak values of attraction lateral forces are about 2 times of the peak values of repulsion forces both in meeting and overtaking conditions at various water depths.Although all the 3 phases and change tendency of the lateral forces in both meeting and overtaking cases are the same,the peak values of correspondingly phases of overtaking case are about 13.5%smaller than those of meeting case at h/D=1.15,while with the increase of water depth,the difference is not so obvious,which means that the lateral force of meeting is more susceptible to water depth than that of overtaking.On the contrary,the yaw moment of meeting is much smaller than that of overtaking,only nearly one half.

        Fig.4 Hydrodynamic force and moment coefficients of Ship 1 at different water depth during meeting,Sp=1.5B,U1=0.358 9 m/s,U2=0.358 9 m/s

        Fig.5 Hydrodynamic force and moment coefficients of Ship 1 at different water depth during overtaking,Sp=1.5B,U1=0.358 9 m/s,U2=0.538 4 m/s

        3.3 Influence of lateral distance between ships

        To investigate the influence of the lateral distance between ships on the hydrodynamic interaction forces,calculations are conducted for different lateral distance between ships Sp= 1.5B,2.0B,2.5B and 3.0B,where the water depth to draught ratio h/D=1.15 keeps unchanged. Fig.6 and Fig.7 show the numerical results.

        Fig.6 Hydrodynamic force and moment coefficients of Ship 1 at different lateral distance between ships during meeting,h/D=1.15,U1=0.358 9 m/s,U2=0.358 9 m/s

        Fig.7 Hydrodynamic force and moment coefficients of Ship 1 at different lateral distance between ships during overtaking,h/D=1.15,U1=0.358 9 m/s,U2=0.538 4 m/s

        It can be seen that the change tendency is same for different lateral distance between ships,and the magnitudes of the hydrodynamic interaction force and moment decrease with theincrease of lateral distance between ships both in meeting and overtaking conditions.The similar tendency can be found,compared to the results of water depth effects.Particularly,the peak values of both lateral force and yaw moment are changing rapidly when Sp=1.5B,and the peak values of them are particularly large.It is obvious that the interaction forces are more susceptible to lateral distance in overtaking case than that in meeting case.

        3.4 Influence of ship speeds during meeting

        To investigate the influence of ship speeds,calculations are carried out for different ratio of ship speeds U1/U2=2.0,1.5,1.0,1/1.5 and 1/2,under h/D=1.15 and Sp=1.5B.The results are shown in Fig.8.

        Fig.8 Hydrodynamic force and moment coefficients of Ship 1 at different ratio of ship speeds,h/D=1.15,Sp=1.5B

        The change tendency of the hydrodynamic force and moment is the same as the water depth effect as described in Fig.4.The smaller the ratio of ship speeds U1/U2,the larger the force and moment on Ship 1 are.The same tendency can be found in overtaking condition.

        3.5 Interaction forces on both ships during overtaking

        In Fig.9,the hydrodynamic interaction force and moment acting on both ships during overtaking are compared.

        Fig.9 Hydrodynamic force and moment coefficients of both ships during overtaking, h/D=1.15,Sp=1.5B,U1=0.358 9 m/s,U2=0.538 4 m/s

        As it can be seen in Fig.9,the change tendency of the lateral force and yaw moment acting on both ships is same,but the peak values of lateral force and yaw moment coefficients of the slower ship(overtaken ship)are about 2 times those of the faster ship(overtaking ship) when the speed ratio U2/U1=1.5.

        4 Concluding remarks

        The ship-ship interaction forces both in meeting and overtaking conditions have been numerically studied by using a high-order panel method.The boundary integral equation is solved at each time step.Under the assumption of low ship speed,the effect of free surface elevation is ignored.Comparisons with slender body theory and CFD method show that the present method is effective.

        The influences of finite water depth and lateral distance between ships are analyzed in detail.It is concluded that the change tendencies of the interaction forces are qualitatively same at different water depth and lateral distance between ships,but the peak values of interaction forces decrease with the increase of these parameters.Moreover,it can be found that lateral force and yaw moment are more susceptible to water depth in meeting case than those in overtaking case.On the contrary,the lateral force and yaw moment are more susceptible to lateral distance during overtaking than those during meeting.

        The interaction forces and moments acting on the overtaken ship are much larger than those acting on the overtaking ship.Typically,under the condition of h/D=1.15,Sp=1.5B,and the speeds of overtaken and overtaking ships U1=0.358 9 m/s and U2=0.538 4 m/s,the interaction forces and moments acting on the overtaken ship are almost 2 times of those acting on the overtaking ship.

        The present high-order panel method based on NURBS is able to predict qualitatively the tendency of the interaction forces.The further study will focus on extending the method to predict the hydrodynamic interaction between real ships with linear boundary conditions on the free surface being taken into consideration.

        [1]Vantorre M,Verzhbitskaya E,Laforce E.Model test based formulations of ship-ship interaction forces[J].Ship Technology Research,2002,49:124-139.

        [2]Tuck E O,Newman J N.Hydrodynamic interactions between ships[C]//10th Symposium on Naval Hydrodynamics.Cambridge,Mass.,USA,1974:35-70.

        [3]Yeung R W.On the interactions of slender ships in shallow water[J].Journal of Fluid Mechanics,1978,85(Part 1):143-159.

        [4]Yeung R W,Tan W T.Hydrodynamic interactions of ships with fixed obstacles[J].Journal of Ship Research,1980,24(1): 50-59.

        [5]Hsiung C C,Gui Q Y.Computing interaction forces and moments on a ship in restricted waterways[J].International Shipbuilding Progress,1988,35(403):219-254.

        [6]Korsmeyer F T,Lee C H,Newman J N.Computation of ship interaction forces in restricted waters[J].Journal of Ship Research,1993,37(4):298-306.

        [7]Yasukawa H.Simulation of ship collision caused by hydrodynamic interaction between ships[C]//MARSIM’03,International Conference on Marine Simulation and Ship Maneuverability.Kanazawa,Japan,2003.

        [8]Xiang X,Faltinsen O M.Maneuvering of two interacting ships in calm water[C].11th International Symposium on Practical Design of Ships and Other Floating Structures.Rio de Janeiro,Brazil,2010.

        [9]Zhou X Q,Sutulo S,Guedes Soares C.Computation of ship hydrodynamic interaction forces in restricted waters using potential theory[J].J Marine Sci.Appl.,2012(11):265-275.

        [10]Sutulo S,Guedes Soares C,Otzen J F.Validation of potential-flow estimation of interaction forces acting upon ship hulls in parallel motion[J].Journal of Ship Research,2012,56(3):129-145.

        [11]Zhang X T,Teng B,Liu Z Y,et al.Computation of hydrodynamic interaction forces between ships based on B-splines[J]. Journal of Hydrodynamics,Ser.B,2003(2):83-88.

        [12]Xu H F,Zou Z J.Prediction of hydrodynamic forces on a moored ship induced by a passing ship in shallow water by using a high-order panel method[J]Journal of Shanghai Jiaotong University,2016,21(2):129-135.

        [13]Varyani K S,Krishnankutty P.Modification of ship hydrodynamic interaction forces and moment by underwater ship geometry[J].Ocean Engineering,2006,33(8-9):1090-1104.

        基于高階面元法的淺水域中船—船水動(dòng)力相互作用數(shù)值預(yù)報(bào)

        徐華福a,b,鄒早建a,b,劉曉艷a

        (上海交通大學(xué)a.船舶海洋與建筑工程學(xué)院;b.海洋工程國家重點(diǎn)實(shí)驗(yàn)室,上海200240)

        文章開發(fā)了一種基于非均勻有理B樣條(NURBS)的三維高階面元法對(duì)淺水域中兩船會(huì)遇和追越時(shí)的船-船水動(dòng)力相互作用進(jìn)行預(yù)報(bào)。該方法采用NURBS精確表達(dá)船體曲面,采用配置點(diǎn)法滿足離散的邊界積分方程,在邊界面上的源強(qiáng)分布確定之后采用B樣條表達(dá)物面速度勢(shì)分布。基于低速假設(shè),忽略了自由面的興波效應(yīng),采用無窮鏡像法處理剛性自由面和淺水效應(yīng),采用時(shí)間步進(jìn)法在時(shí)域內(nèi)求解速度勢(shì)。將文中計(jì)算結(jié)果與基于RANS方程求解的CFD方法和細(xì)長體理論方法計(jì)算結(jié)果比較,驗(yàn)證了該方法的有效性。在此基礎(chǔ)上,分別在不同水深、船間橫向間距和船速下進(jìn)行了系列的計(jì)算,分析了這些因素對(duì)船-船水動(dòng)力相互作用影響。結(jié)果表明,在相同的工況下,船-船會(huì)遇時(shí)相互作用力比追越時(shí)更大,對(duì)水深變化更敏感,而對(duì)橫向間距不如追越時(shí)敏感。

        會(huì)遇;追越;船—船水動(dòng)力相互作用;NURBS;高階面元法

        U661.1

        A

        徐華福(1984-),男,上海交通大學(xué)博士研究生;鄒早建(1956-),男,上海交通大學(xué)教授,博士生導(dǎo)師;劉曉艷(1987-),女,上海交通大學(xué)博士研究生。

        U661.1 < class="emphasis_bold">Document code:A

        A

        10.3969/j.issn.1007-7294.2016.12.004

        1007-7294(2016)12-1535-12

        Received date:2016-06-29

        Foundation item:Supported by National Natural Science Foundation of China(51179019,51309152)

        Biography:XU Hua-fu(1984-),male,Ph.D.candidate of Shanghai Jiao Tong University,E-mail:huafuxu@sjtu.edu.cn; ZOU Zao-jian(1956-),male,professor/tutor,corresponding author,E-mail:zjzou@sjtu.edu.cn.

        猜你喜歡
        會(huì)遇上海交通大學(xué)元法
        上海交通大學(xué)
        基于AIS數(shù)據(jù)的船舶會(huì)遇挖掘與分析
        基于AIS數(shù)據(jù)的交匯水域船舶會(huì)遇態(tài)勢(shì)辨識(shí)
        中國航海(2021年1期)2021-03-10 13:31:32
        換元法在解題中的運(yùn)用
        基于離散元法的礦石對(duì)溜槽沖擊力的模擬研究
        上海交通大學(xué)參加機(jī)器人比賽
        一種基于DBSCAN的船舶會(huì)遇實(shí)時(shí)識(shí)別方法
        換元法在解題中的應(yīng)用
        “微元法”在含電容器電路中的應(yīng)用
        基于AIS的特定船舶會(huì)遇實(shí)況分布
        中國航海(2014年3期)2014-11-28 11:17:08
        久久久久99精品成人片直播| 中文字幕av人妻一区二区| 国产一区二区三区中出| 色欲av永久无码精品无码蜜桃| 日躁夜躁狠狠躁2001| 国产精品高清视亚洲乱码有限公司 | 国产精品 无码专区| 国产一级二级三级在线观看视频| 国产啪亚洲国产精品无码| 日产精品久久久久久久| 日韩午夜在线视频观看| 狠色人妻丝袜中文字幕| 日韩aⅴ人妻无码一区二区| 亚洲乱码视频在线观看| 初尝人妻少妇中文字幕在线| 国内自拍色第一页第二页| 亚洲日本一区二区一本一道| 欧美精品videosse精子| 免费观看又污又黄的网站| 国产成人福利在线视频不卡| 全国一区二区三区女厕偷拍| 夜夜爽日日澡人人添| √天堂中文官网8在线| 亚洲人成网站在线播放小说| 少妇高潮久久蜜柚av| 精品免费久久久久久久| 国产污污视频| 精品蜜桃在线观看一区二区三区 | 天美麻花果冻视频大全英文版| 激情五月婷婷久久综合| 日韩av毛片在线观看| 国产亚洲精品久久久ai换| 亚洲VA中文字幕欧美VA丝袜| 国产福利一区二区三区在线观看| 成人精品一区二区三区电影| 无码人妻精品一区二区三区66| 久久精品国产亚洲av桥本有菜 | 国产成人无码综合亚洲日韩| 5级做人爱c视版免费视频| 日韩av免费在线不卡一区| 亚洲综合色区一区二区三区|