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

        ?

        水面比降對(duì)船舶航行阻力影響的數(shù)值研究

        2023-01-03 08:12:24舒麟棹錢志鵬
        關(guān)鍵詞:船模船體航行

        冀 楠,楊 光,舒麟棹,錢志鵬,楊 春

        (重慶交通大學(xué) 航運(yùn)與船舶工程學(xué)院,重慶 400074)

        0 引 言

        船舶航行過程中經(jīng)常會(huì)遇到有水面比降的情況,例如:明渠河道、受洪水漲落或潮汐影響的河道(長(zhǎng)江南京至鎮(zhèn)江河段[1])或者河底高程沿流向逐漸降低的河道水面均會(huì)產(chǎn)生比降;海洋中長(zhǎng)波波面也可以認(rèn)為是有比降水面;山區(qū)河流因河道高程和寬度多變, 水面比降現(xiàn)象更常見。船舶航行于帶有比降的水面時(shí)會(huì)產(chǎn)生不同于一般水面的水動(dòng)力特性,導(dǎo)致船舶航行姿態(tài)改變,進(jìn)而影響航速甚至操縱性能的惡化,容易造成事故。因此,為了提高船舶航行在有比降水面時(shí)的安全性,有必要對(duì)水面比降對(duì)船舶的航行阻力與浮態(tài)影響進(jìn)行研究。

        對(duì)水面比降影響船舶航行問題的研究主要有基于經(jīng)驗(yàn)公式的理論分析方法和模型試驗(yàn)法,研究范圍基本集中在內(nèi)河航道的上灘問題上?;诮?jīng)驗(yàn)公式的理論分析方面,曹民雄等[2]利用急流灘模型的船模上灘試驗(yàn)資料,對(duì)比分析了幾種內(nèi)河船舶上灘阻力計(jì)算方法;周曉杰等[3]根據(jù)航道規(guī)范和理論計(jì)算,對(duì)流速-比降組合值對(duì)應(yīng)的船舶上灘指標(biāo)進(jìn)行計(jì)算分析,得到了研究河段對(duì)應(yīng)的規(guī)范設(shè)計(jì)尺度;董曉韜[4]采用船舶推力-航行阻力平衡的方法,計(jì)算多種功率、載重條件下船舶自航上灘的臨界流速-比降組合,從而獲取多種工況下的急灘通航水力指標(biāo);張鵬等[5]結(jié)合船舶在洪水急流灘自航上灘的實(shí)測(cè)資料,確定了三峽庫(kù)區(qū)洪水急流灘河段船舶自航上灘的水力指標(biāo);李慶寧等[6]研究了三峽-葛洲壩航段航行的3種船型在泄洪期安全上水過灘的通航能力;童思陳等[7]以山區(qū)河流瀾滄江為代表,提出了適合于瀾滄江代表船舶的水流阻力計(jì)算式,并且結(jié)合坡降阻力計(jì)算, 建立了瀾滄江急流灘的通航水力指標(biāo)。模型試驗(yàn)方面,曹民雄等[8]進(jìn)行了多種水力條件下的船模上灘試驗(yàn),發(fā)現(xiàn)水面比降、水流速度及船舶水深吃水比是船舶上灘的主要影響因素;劉曉菲等[9]開展了烏江小幺灘河段航道整治物理模型試驗(yàn),通過調(diào)整河段流速和比降,改善水流條件,以滿足船舶安全航行要求。

        截至目前,國(guó)內(nèi)外學(xué)者對(duì)涉及船舶航行于比降水面的問題主要集中在評(píng)估水力特性和通航指標(biāo)上,對(duì)航行阻力的計(jì)算采用的是經(jīng)驗(yàn)公式或模型試驗(yàn)。筆者采用計(jì)算流體力學(xué)(CFD)方法,基于商業(yè)軟件STAR-CCM+,以某型集裝箱船為研究對(duì)象,通過改變河道斷面縱向收縮角和入口流速獲得不同的水面比降工況,分析了水面比降對(duì)船舶航行阻力與浮態(tài)的影響,為船舶航行于帶有比降水面時(shí)的安全決策提供一定的參考。

        1 船舶航行于比降水面阻力分析

        船舶航行于比降水面時(shí),由于水面坡度的影響,船舶的浮態(tài)及阻力都會(huì)區(qū)別于一般水面。如圖1,根據(jù)文獻(xiàn)[5]~文獻(xiàn)[7]的分析,船舶處于比降水面臨界狀態(tài)時(shí),會(huì)產(chǎn)生縱傾,船舶重力沿船體縱向的分力被定義為比降阻力。此時(shí)船舶航行的總阻力和比降阻力的經(jīng)驗(yàn)計(jì)算公式為:

        R=Rv+Rj

        (1)

        式中:Rv為水流阻力,N;Rj為比降阻力,N。

        比降阻力Rj的計(jì)算公式為:

        Rj=βWJ

        (2)

        式中:W為船舶排水量,N;J為船長(zhǎng)范圍內(nèi)平均水面比降;β為考慮船舶上灘時(shí)水面比降局部增大的修正系數(shù),一般取β=1.1~1.2。

        圖1 船舶航行于比降水面Fig. 1 The ship sailing on the water surface gradient

        船舶航行時(shí),由于水面存在一定比降,船首尾的水位差導(dǎo)致船舶傾斜。船傾斜后,浮心位置會(huì)發(fā)生偏移,致使浮心和重心不在同一鉛垂線上,于是作用方向相反的重力和浮力會(huì)促使船舶縱傾角增大,增大后的縱傾角設(shè)為θ′,θ′理論上要大于比降對(duì)應(yīng)的水面坡度角θ。因此,實(shí)際上船舶重力在船體縱向的分力應(yīng)該為Wsinθ′,這才是比降阻力。由此可知,采用式(2)計(jì)算得到的比降阻力會(huì)小于實(shí)際值,從而導(dǎo)致利用經(jīng)驗(yàn)公式計(jì)算得到的船舶航行于比降水面時(shí)的總阻力會(huì)偏小。

        2 數(shù)值方法

        2.1 控制方程

        在空間固定坐標(biāo)系下,數(shù)值模擬采用的是不可壓縮流動(dòng)RANS方程,控制方程為連續(xù)性方程和雷諾平均方程,如式(3)和式(4):

        (3)

        (4)

        2.2 湍流模型

        考慮計(jì)算時(shí)間和計(jì)算精度,選用Realizablek-ε湍流模型,湍流動(dòng)能k及湍流耗散率ε方程如式(5)和式(6):

        (5)

        (6)

        式中:k為湍流動(dòng)能;ε為湍流耗散率;μt為湍流黏性系數(shù);C1=1.44;C2=1.9;σk=1.0;σε=1.2。其中μt為:

        (7)

        式中:Cμ=0.09。

        3 計(jì)算模型與驗(yàn)證

        3.1 計(jì)算模型

        以國(guó)際研討會(huì)標(biāo)準(zhǔn)集裝箱船KCS為研究對(duì)象,如圖2。船模主要參數(shù)為:垂線間長(zhǎng)Lpp=7.278 6 m,型寬B=1.019 m,型深D=0.341 8 m,設(shè)計(jì)吃水T=0.34 m,濕表面積S=9.438 m2,方形系數(shù)Cb=0.65,模型縮尺比λ=31.67。

        KCS雖為海船船型,但其船型的長(zhǎng)、寬、吃水與長(zhǎng)江過閘船型集-18[10]的相應(yīng)參數(shù)基本呈2∶1的尺度比。雖然KCS船型的方形系數(shù)0.65比集-18的方向系數(shù)0.53大,但根據(jù)船舶原理[11]中關(guān)于方形系數(shù)對(duì)船舶阻力影響的敘述可知,在低于1.93 m/s的航速范圍內(nèi),方形系數(shù)差異造成的船型阻力變化僅為10%左右。

        圖2 KCS船幾何模型Fig. 2 KCS ship geometric model

        3.2 船舶浮態(tài)與阻力之間關(guān)系的驗(yàn)證

        船舶在比降水面航行時(shí)會(huì)由于水面坡度而產(chǎn)生不同程度的尾傾,進(jìn)而造成船舶阻力的變化。為了驗(yàn)證文中數(shù)值模型的合理性,將船舶在不同初始縱傾角下的阻力數(shù)值計(jì)算結(jié)果與文獻(xiàn)[12]進(jìn)行了對(duì)比。

        表1 網(wǎng)格基本參數(shù)信息Table 1 Grid basic parameter information

        表2 船??v傾阻力數(shù)值計(jì)算與試驗(yàn)對(duì)比Table 2 Comparison between numerical calculation and experiment of trim resistance of ship model

        由表2可以發(fā)現(xiàn),三套網(wǎng)格的計(jì)算結(jié)果與試驗(yàn)的擬合度較好,最大誤差僅為3.21%。網(wǎng)格1在縱傾角小于0.5°后,與試驗(yàn)結(jié)果存在一定偏離,特別是在縱傾角為0.25°和-0.5°時(shí)的誤差較大;網(wǎng)格3在縱傾角0.25°時(shí)誤差最大,大于0.25°時(shí)的計(jì)算結(jié)果較網(wǎng)格2偏差也較大。在綜合考慮計(jì)算精度與時(shí)間成本后,選擇網(wǎng)格2的屬性設(shè)置參數(shù)作為數(shù)值計(jì)算的網(wǎng)格劃分參數(shù)。

        3.3 計(jì)算工況設(shè)置

        計(jì)算水深設(shè)定為8 m,水深吃水比為h/T=23.5。通過改變河道斷面縱向收縮角和入口流速V來獲得不同的水面比降。河道斷面縱向收縮角是邊坡與x方向的夾角α,如圖3。筆者設(shè)置夾角α分別為1°、2°、3°、5°,入口流速V分別為0.5、1.25、2 m/s,以α-V來表示對(duì)應(yīng)工況。數(shù)值計(jì)算一共設(shè)置3組工況,第一組工況:1-0.5、2-0.5、3-0.5、5-0.5;第二組工況:1-1.25、2-1.25、3-1.25、5-1.25;第三組工況:1-2、2-2、3-2、5-2。

        船體坐標(biāo)系如圖4,原點(diǎn)O位于船模重心,x軸指向船首,y軸指向左側(cè)岸壁,z軸為垂直水面向上。船模首傾時(shí)縱傾角為正值,尾傾時(shí)縱傾角為負(fù)值。

        圖3 計(jì)算域平面Fig. 3 Plan of calculation domain

        圖4 船體坐標(biāo)系Fig. 4 Ship hull coordinate system

        3.4 邊界條件及網(wǎng)格劃分

        假定船舶沿航道中線航行,且航道關(guān)于航道中線對(duì)稱,因此,航道和船舶均可采用半模進(jìn)行設(shè)置和計(jì)算。邊界條件的設(shè)置如圖5,計(jì)算域的進(jìn)口、頂部及底部均設(shè)置為速度進(jìn)口,出口設(shè)置為壓力出口,側(cè)邊及船體表面設(shè)置為壁面,計(jì)算域中縱剖面設(shè)置為對(duì)稱平面。計(jì)算域中,進(jìn)口距船首2Lpp,出口距船尾3Lpp,出口的寬度固定為0.7Lpp,邊坡長(zhǎng)度固定為3.3Lpp。設(shè)定船模距岸壁的距離超過一倍船長(zhǎng),水深吃水比h/T=23.5,認(rèn)為船舶不受岸壁和淺水效應(yīng)的影響,只在帶有比降水面的敞水水域航行。數(shù)值計(jì)算放開了船體升沉和縱傾兩個(gè)自由度;求解時(shí)間離散采用一階;擴(kuò)散項(xiàng)采用二階中心差分格式離散;對(duì)流項(xiàng)采用二階迎風(fēng)格式,以提高計(jì)算精度;壓力-速度耦合方程采用SIMPLEC算法進(jìn)行求解。

        網(wǎng)格采用切割體網(wǎng)格,對(duì)自由液面及船體周圍進(jìn)行體加密,如圖6。船體首尾型線曲度較大,需進(jìn)行體加密,如圖7。經(jīng)過多次試算,最終確定無(wú)量綱化的壁面距離y+值的范圍在30~60之間,說明船體邊界層的網(wǎng)格設(shè)置是合理的。

        圖5 計(jì)算域邊界條件Fig. 5 Boundary conditions of calculation domain

        圖6 計(jì)算域整體網(wǎng)格Fig. 6 Global grid of calculation domain

        圖7 船首尾網(wǎng)格Fig. 7 Bow and stern grid

        4 結(jié)果分析

        4.1 阻力分析

        表3匯總了3組工況計(jì)算的阻力值,其中水流阻力作為對(duì)比項(xiàng),是在提取到對(duì)應(yīng)比降工況下船體附近平均流速數(shù)值后,另外計(jì)算的船舶在敞水中的阻力值。比降阻力的數(shù)值模擬值由總阻力與水流阻力的差值得到,比降阻力的經(jīng)驗(yàn)公式值由式(2)求得,并且也給出了經(jīng)驗(yàn)公式值與數(shù)值模擬值的差值百分比。從表3可知:每組工況中,隨著比降的增加,船舶的總阻力逐漸增大,意味著航行難度增加;當(dāng)入口流速一定,每組工況的水流阻力相差并不大,水流阻力主要由船速?zèng)Q定,或者說是固定船模數(shù)值計(jì)算中的入口水流流速?zèng)Q定。在每組工況中,船舶周圍的流速與入口流速相比變化也不大,這說明式(1)將船舶在有比降水面上航行時(shí)的總阻力分解為水流阻力和比降阻力兩個(gè)部分的思路是合理的。

        表3給出了3組工況的比降阻力數(shù)值模擬值和經(jīng)驗(yàn)公式值隨水面比降變化的對(duì)比。由表3可以發(fā)現(xiàn),采用經(jīng)驗(yàn)公式計(jì)算的比降阻力明顯偏小,這與第1節(jié)理論分析得出的結(jié)論一致。而且,在每組工況中,隨著比降的增大,相對(duì)誤差又在縮小,說明式(2)在大比降時(shí)計(jì)算準(zhǔn)確度高,而并不適用于小比降時(shí)的情況。究其原因,應(yīng)該是因?yàn)楸冉递^小時(shí),船舶重力浮力力矩造成的船舶縱傾角調(diào)整相對(duì)影響較大,而在大比降時(shí),由于比降數(shù)值大,因此船舶重力浮力力矩造成的船舶縱傾角調(diào)整相對(duì)較小,比降阻力的數(shù)值趨向于直接用比降坡度角計(jì)算的結(jié)果。

        圖8是3組工況的摩擦阻力Rf和壓阻力Rp隨不同比降變化的曲線。

        圖8 不同比降下的摩擦阻力和壓阻力曲線Fig. 8 The frictional resistance and pressure resistance curves under different gradient

        從圖8中可以看出,阻力的主要成分是壓阻力Rp(粘壓和興波阻力之和),且壓阻力Rp會(huì)隨著水面比降的增大而明顯增大;摩擦阻力Rf隨水面比降的變化較小,在高流速工況下,甚至有一定的減小。摩擦阻力Rf與船體附近流速以及船模的濕表面積Sm有關(guān)。表4給出了3組工況船模在不同比降下的濕表面積Sm。由表4可以發(fā)現(xiàn):每組工況內(nèi),濕表面積Sm變化很?。还r之間,由于入口流速變大,船模濕表面積Sm有一定的增加,但增加數(shù)值仍然很小。由此可知,船舶在有比降水面航行時(shí),船舶摩擦阻力Rf的變化主要由船體附近流速?zèng)Q定,而與濕表面積Sm的變化關(guān)系不大。

        表4 不同比降下KCS船模濕表面積Table 4 The wet surface area of KCS ship model under different gradient

        圖9展示了第一組工況船體底部動(dòng)壓分布情況,從圖9中可知:比降為0.1‰ 時(shí),低壓區(qū)覆蓋船底首部到尾部的全部區(qū)域;比降為0.2‰ 時(shí),船首部壓力增大,船底低壓區(qū)縮減到僅覆蓋尾部區(qū)域;比降為0.4‰ 時(shí),船底低壓區(qū)進(jìn)一步縮減到僅覆蓋尾肩部;比降為0.6‰ 時(shí),雖然流速幾乎沒變,但比降的增大仍然使得船底的壓力數(shù)值整體增大。由此可見,比降的增大,不但會(huì)增大船體首尾的壓力差,也會(huì)增大船體所受水壓力的數(shù)值。

        圖9 船體底部壓力隨比降變化分布Fig. 9 Distribution of hull bottom pressure changing with gradient

        4.2 自由液面波形和船舶浮態(tài)分析

        圖10為第三組工況的自由液面波形等值線。圖11為3組工況距船中縱剖面1倍的船寬B處波高H沿船長(zhǎng)X方向的分布。由圖10、圖11可知:每組工況下,隨著比降的增大,波面的平均傾角會(huì)增大。3組工況間的對(duì)比規(guī)律為:隨著船速的提高,船體興波對(duì)水面的擾動(dòng)作用逐漸顯著,船首壅水導(dǎo)致首波峰的幅值逐漸增大,船尾后部的波浪起伏也逐漸明顯。從能量的角度看,波浪起伏越大,需要消耗越多的能量,這一部分能量最終以阻力的形式體現(xiàn),這也是船舶航行于比降水面時(shí)阻力增大的原因之一。

        圖10 自由液面波形等值線Fig. 10 Free surface waveform contour

        圖12是3組工況下船體縱傾角隨水面比降變化的曲線,并且給出了水面坡度角作為對(duì)比,因?yàn)槲矁A時(shí)船體縱傾角為負(fù)值,為了便于比較,圖12中的水面坡度角取負(fù)值。由圖12可發(fā)現(xiàn),船舶在上灘航行時(shí)船身會(huì)產(chǎn)生尾傾現(xiàn)象。隨著水面比降的增大,船身尾傾更嚴(yán)重,并且尾傾角度大于水面坡度角,這也進(jìn)一步證明采用式(2)計(jì)算出的比降阻力值會(huì)偏小。

        圖13是不同流速下,船尾自由面隨水面比降變化分布,數(shù)值仿真中采用VOF方法捕捉船體自由液面。從圖13中可以明顯看到,隨著水面比降增大,自由液面高度在逐漸抬升,說明船體尾傾角度在增大。

        圖11 不同比降下的波高分布Fig. 11 Wave elevation distribution under different gradient

        圖12 縱傾隨比降的變化曲線Fig. 12 Variation curves of trim changing with gradient

        圖13 船尾自由面隨比降變化分布Fig. 13 Distribution of free surface of stern changing with gradient

        4.3 流場(chǎng)特性分析

        槳盤處伴流場(chǎng)是計(jì)算螺旋槳推力及主機(jī)功率的重要參數(shù),因此需研究水面比降對(duì)槳盤處伴流場(chǎng)的影響[14]。單槳船的槳盤面關(guān)于中縱剖面左右對(duì)稱,為方便對(duì)比,采用圖14進(jìn)行分析,圖14中Wx為槳盤處標(biāo)稱伴流分?jǐn)?shù)。由圖14可知:圖14(a)中,槳盤面中心處的高伴流區(qū)隨著比降的增大呈現(xiàn)擴(kuò)大的趨勢(shì),槳盤面下部的伴流隨著比降的增大有向外擴(kuò)展的趨勢(shì),即隨著比降的增大,槳盤處的平均伴流增大;圖14(b)中,槳盤面下部的伴流隨著比降增大仍呈現(xiàn)向外擴(kuò)展的趨勢(shì),但上部水流卻隨著比降的增大呈收縮趨勢(shì),說明高伴流區(qū)在向下移動(dòng),伴流的變化必將對(duì)螺旋槳的性能造成影響。

        圖14 第三組工況下標(biāo)稱伴流分?jǐn)?shù)Wx隨比降變化的分布Fig. 14 Distribution of nominal wake fraction changing with gradient in the third group of working conditions

        圖15~圖17分別展示了第三組工況船底、船首、船尾渦量細(xì)節(jié)分布,渦量顯示根據(jù)Q準(zhǔn)則得到。從圖15中可以看到船首舭渦、側(cè)渦、肩渦和船尾舭渦、船尾渦等紊亂的渦流。對(duì)比4個(gè)工況下的船底渦量分布可以發(fā)現(xiàn),隨著水面比降的增大,船底渦量強(qiáng)度會(huì)逐漸減弱。水面比降為1.7‰ 時(shí),船底渦量強(qiáng)度最高,船首、尾舭渦附著在船底;水面比降增大至7.4‰ 時(shí),舭渦渦量減少,強(qiáng)度也變?nèi)?。從圖16可以發(fā)現(xiàn),水面比降為1.7‰ 時(shí),船首側(cè)渦和肩渦清晰可見,并且呈現(xiàn)融合狀態(tài);比降增大至3.3‰ 時(shí),側(cè)渦、肩渦分離,且側(cè)渦強(qiáng)度已開始減弱;比降增大至4.6‰ 時(shí),肩渦強(qiáng)度開始減弱;當(dāng)水面比降為7.4‰ 時(shí),船首側(cè)渦分離成兩小段,兩側(cè)翼狀結(jié)構(gòu)消失。從圖17可以發(fā)現(xiàn),水面比降增加時(shí),船尾渦的強(qiáng)度也會(huì)減弱,水面比降增大到4.6‰ 時(shí),船尾渦的尾翼狀結(jié)構(gòu)消失。

        圖15 船底渦系分布Fig. 15 Diseribution of vontex system at bottem of hull

        圖16 船首渦量Fig. 16 Bow vorticity

        圖17 船尾渦量Fig. 17 Stern vorticity

        5 結(jié) 論

        筆者以實(shí)驗(yàn)船模為研究對(duì)象,對(duì)船舶航行于有比降水面進(jìn)行了數(shù)值模擬,分析了船舶的阻力、興波、浮態(tài)和黏性流場(chǎng)的變化。得出了以下主要結(jié)論:

        1)計(jì)算流體力學(xué)(CFD)方法可以用來模擬船舶遭遇有比降水面航行時(shí)的阻力變化和航行特性,評(píng)估船舶航行能力。

        2)船舶航行在帶有比降的水面時(shí),水面比降越大,船舶所受到的總阻力也越大,增加了航行難度;水面比降對(duì)摩擦阻力的影響小,對(duì)壓阻力的影響較大。

        3)船舶在有比降水面航行會(huì)產(chǎn)生一定量的尾傾。隨著水面比降的增大,船身尾傾會(huì)更嚴(yán)重,并且尾傾角度會(huì)大于比降對(duì)應(yīng)水面坡度角,這將導(dǎo)致比降阻力的真實(shí)結(jié)果比經(jīng)驗(yàn)公式計(jì)算出的數(shù)值大。

        4)船舶航行在帶有比降的水面時(shí)的流場(chǎng)特性表現(xiàn)為:船尾槳盤處的伴流分?jǐn)?shù)會(huì)隨著水面比降的增大而增大,影響螺旋槳的運(yùn)轉(zhuǎn)效率,增加航行難度;水面比降的變化會(huì)影響船體渦量強(qiáng)度和渦系的發(fā)展,渦量強(qiáng)度隨著水面比降的增大會(huì)減弱。

        猜你喜歡
        船模船體航行
        基于模糊PID的船模航向控制研究
        船體行駛過程中的壓力監(jiān)測(cè)方法
        到慧骃國(guó)的航行
        小舟在河上航行
        航行
        青年歌聲(2017年6期)2017-03-13 00:57:56
        船模靜水橫搖試驗(yàn)的不確定度分析
        焊接殘余應(yīng)力對(duì)船體結(jié)構(gòu)疲勞強(qiáng)度的影響分析
        焊接(2015年9期)2015-07-18 11:03:51
        赴美軍“仁慈”號(hào)醫(yī)院船駐船體會(huì)
        西洋船模王——童鑑良
        航海(2014年6期)2014-12-12 10:36:03
        水下爆炸氣泡作用下船體總縱強(qiáng)度估算方法
        狠狠摸狠狠澡| 日韩精品中文字幕人妻系列| 日本激情一区二区三区| 亚洲AV秘 无码一区二区三| 国产日韩AV无码免费一区二区| 黑丝美女喷水在线观看| 国产午夜精品av一区二区三| 日韩av一区二区三区激情在线| 色777狠狠狠综合| 欧美aaaaaa级午夜福利视频| 国产精品无码不卡一区二区三区| 国产男女猛烈无遮挡免费视频| 亚洲色无码中文字幕| 亚洲婷婷久久播66性av| 国产精品毛片无遮挡| 精品国产乱码久久久久久口爆网站| 国产h视频在线观看网站免费| 亚洲av影片一区二区三区| 日本久久精品视频免费| 人禽杂交18禁网站免费| 日本牲交大片免费观看| 女人与牲口性恔配视频免费| 国产亚洲美女精品久久久2020| 中文字幕有码在线亚洲| 亚洲国产成人av二区| 国产精品无码久久久久| 久久亚洲AV成人一二三区| 亚洲天堂无码AV一二三四区| 日本在线一区二区在线| 天天色天天操天天日天天射| 亚洲av综合av成人小说| 国产成人亚洲精品无码mp4| 色综合久久加勒比高清88| 亚洲熟女少妇精品久久| 999国产精品999久久久久久| 午夜福利院电影| 男人的天堂在线无码视频| 日本av第一区第二区| 亚洲色精品三区二区一区| 久久中文字幕人妻熟av女蜜柚m| 一本无码人妻在中文字幕|