999精品在线视频,手机成人午夜在线视频,久久不卡国产精品无码,中日无码在线观看,成人av手机在线观看,日韩精品亚洲一区中文字幕,亚洲av无码人妻,四虎国产在线观看 ?

非常規(guī)內(nèi)傾船型騎浪數(shù)值預(yù)報方法研究

2022-08-17 11:24:10儲紀(jì)龍
船舶力學(xué) 2022年8期
關(guān)鍵詞:船舶

儲紀(jì)龍,魯 江,顧 民

(中國船舶科學(xué)研究中心 水動力學(xué)重點(diǎn)實(shí)驗(yàn)室,江蘇 無錫 214082)

0 引 言

國際海事組織(IMO)于2020年12月發(fā)布了第二代完整穩(wěn)性衡準(zhǔn)臨時指南[1],目前,正在制定第二代完整穩(wěn)性衡準(zhǔn)直接評估方法解釋性文件。二代完整穩(wěn)性衡準(zhǔn)改變以往只依靠經(jīng)驗(yàn)公式制定衡準(zhǔn)的方法,涵蓋了五種波浪穩(wěn)性失效模式:參數(shù)橫搖、純穩(wěn)性喪失、癱船穩(wěn)性、過度加速度和騎浪/橫甩,其中騎浪/橫甩涉及波浪中船舶操縱性、耐波性、快速性和波浪穩(wěn)性等多學(xué)科的交叉耦合,是一種強(qiáng)非線性的復(fù)雜穩(wěn)性失效模式。由于騎浪/橫甩現(xiàn)象的復(fù)雜性和直接穩(wěn)性評估的高難度,目前急需一種準(zhǔn)確有效的騎浪/橫甩數(shù)值預(yù)報方法。所謂騎浪,就是船舶在隨浪或尾斜浪中高速航行時,被波浪捕獲并以波速前進(jìn)的現(xiàn)象。通常船舶在波浪的下坡段發(fā)生騎浪,處于騎浪狀態(tài)下的船舶,多數(shù)會因航向不穩(wěn)定性而不可控制地轉(zhuǎn)向,發(fā)生橫甩,橫甩產(chǎn)生的離心力嚴(yán)重時可以導(dǎo)致船舶傾覆。

第二代完整穩(wěn)性衡準(zhǔn)中每種失效模式衡準(zhǔn)都由三個層次的評估方法組成,分別為第一層薄弱性衡準(zhǔn)、第二層薄弱性衡準(zhǔn)和直接穩(wěn)性評估,三層評估方法的計(jì)算復(fù)雜性依次遞增,評估的準(zhǔn)確性也依次提高。最新提案[1]中給出了騎浪/橫甩直接穩(wěn)性評估的相關(guān)要求:對于騎浪/橫甩失效模式,船舶運(yùn)動數(shù)值模擬至少包括縱蕩、橫蕩、橫搖和首搖四自由度運(yùn)動,不考慮的自由度運(yùn)動應(yīng)該滿足靜平衡假定;船舶運(yùn)動數(shù)值模擬應(yīng)該考慮由于船體漩渦脫落而產(chǎn)生的水動力,除了靜水操縱力,還應(yīng)包括因波浪粒子速度和船舶速度共存而產(chǎn)生的水動升力;采用合適的舵控制方法;初始航速應(yīng)該設(shè)置得足夠小以避免假騎浪(artificial surf-riding)的發(fā)生。

針對騎浪/橫甩數(shù)值預(yù)報方法的研究,國外學(xué)者Umeda 等[2]采用考慮線性波浪力的縱蕩-橫蕩-橫搖-首搖四自由度MMG 模型定性預(yù)報了ITTC A2 漁船的騎浪/橫甩現(xiàn)象;Hashimoto 等[3]拓展已有的縱蕩-橫蕩-橫搖-首搖四自由度MMG 模型預(yù)報了內(nèi)傾船的騎浪/橫甩現(xiàn)象;Matsubara 等[4]基于縱蕩-橫蕩-橫搖-首搖四自由度MMG模型,結(jié)合臨界波法和直接計(jì)數(shù)法計(jì)算了短峰不規(guī)則波中船舶由橫甩導(dǎo)致大角度橫傾的概率;Angelou 等[5]基于六自由度數(shù)學(xué)模型數(shù)值模擬風(fēng)帆游艇的騎浪/橫甩現(xiàn)象,獲得了風(fēng)帆游艇的動態(tài)穩(wěn)定性圖譜和穩(wěn)性失效的臨界波陡;Carrica 等[6-7]應(yīng)用CFD 方法再現(xiàn)了船舶在規(guī)則波和不規(guī)則波中橫甩的發(fā)生過程。國內(nèi)對于騎浪/橫甩數(shù)值模擬方法的研究相對較少,于立偉等[8-9]采用操縱性與耐波性耦合的六自由度弱非線性統(tǒng)一模型,對一艘漁船隨浪條件下的騎浪/橫甩運(yùn)動進(jìn)行數(shù)值模擬,探究漁船的騎浪/橫甩界限,并改進(jìn)舵出水模型,分析舵出水對橫甩運(yùn)動的影響;張寶吉等[10]以四自由度操縱性方程為基礎(chǔ),數(shù)值模擬了漁船在隨浪、尾斜浪下騎浪/橫甩現(xiàn)象,分析了不同波浪條件下漁船的運(yùn)動規(guī)律。

國內(nèi)對于騎浪/橫甩的研究主要集中于漁船的騎浪/橫甩現(xiàn)象數(shù)值預(yù)報,而對于內(nèi)傾船這樣具有雙槳雙舵、細(xì)長體、穿浪船艏、水線以上舷側(cè)內(nèi)傾、大型上層建筑等特性的典型非常規(guī)高速船的騎浪/橫甩數(shù)值模擬幾乎沒有。因此,本文針對非常規(guī)內(nèi)傾船型,建立考慮波浪中雙槳雙舵實(shí)時出入水變化的縱蕩-橫蕩-橫搖-首搖四自由度運(yùn)動耦合的數(shù)學(xué)模型,采用FORTRAN 語言編寫計(jì)算程序,進(jìn)行運(yùn)動仿真,預(yù)報非常規(guī)內(nèi)傾船的騎浪現(xiàn)象,可為騎浪/橫甩直接穩(wěn)性評估奠定基礎(chǔ)。

1 數(shù)學(xué)模型

1.1 船舶運(yùn)動坐標(biāo)系

本文采用兩種坐標(biāo)系,如圖1 所示:空間固定坐標(biāo)系O-ξηζ,原點(diǎn)O位于水平面,起始于波谷,ζ軸向下為正,用來描述波浪;船體坐標(biāo)系G-xyz,以船舶重心G為原點(diǎn),x軸在中線面內(nèi),平行于基面,指向船首為正,z軸向下為正。圖中χ為航向角,φ為橫搖角,δ為舵角,ξG為船舶重心G在空間固定坐標(biāo)系中與原點(diǎn)處初始波谷的縱向相對位置。

圖1 坐標(biāo)系Fig.1 Coordinate systems

1.2 數(shù)學(xué)模型構(gòu)建

船舶發(fā)生騎浪/橫甩時船舶遭遇頻率遠(yuǎn)小于升沉和縱搖的固有頻率,升沉和縱搖運(yùn)動可通過靜平衡法求解穩(wěn)定平衡點(diǎn)來近似[11]。本文基于MMG 模型構(gòu)建縱蕩-橫蕩-首搖-橫搖四自由度數(shù)學(xué)模型進(jìn)行騎浪運(yùn)動數(shù)值預(yù)報,該數(shù)學(xué)模型滿足IMO 提案中對騎浪/橫甩直接穩(wěn)性評估的相關(guān)要求[1]。數(shù)學(xué)模型如下所示:

1.3 船體力

船體力XH,YH,NH和KH的表達(dá)式為

1.4 靜水阻力

靜水阻力R(u)的表達(dá)式為

式中,ρ為水密度,SF為船舶濕表面積,CT為船舶總阻力系數(shù),L為船舶垂線間長,g為重力加速度。

1.5 螺旋槳推力和舵力

IMO 提案中對騎浪/橫甩直接穩(wěn)性評估的相關(guān)要求指出,騎浪/橫甩數(shù)值模擬方法應(yīng)該正確地模擬由于船體漩渦脫落而產(chǎn)生的水動力,也就是除了靜水操縱力,還應(yīng)包括因波浪粒子速度和船舶速度共存產(chǎn)生的水動升力。因此,本文在雙槳雙舵的螺旋槳推力模型和舵力模型里考慮波浪粒子速度的影響,同時還進(jìn)一步考慮了螺旋槳和舵實(shí)時出入水的影響。

雙螺旋槳推進(jìn)時,左右兩側(cè)螺旋槳受波浪粒子速度和實(shí)時出入水情況的影響,產(chǎn)生的推力不同,同時還會產(chǎn)生首搖力矩,這里左舷推力用TP表示,右舷推力用TS表示,首搖力矩用NP表示。雙螺旋槳推力XP和首搖力矩NP的表達(dá)式如下:

式中:nP為螺旋槳轉(zhuǎn)速;DP為螺旋槳直徑;tP為螺旋槳推力減額分?jǐn)?shù);KT為螺旋槳推力系數(shù);JPP、JPS分別為左右側(cè)螺旋槳的進(jìn)速系數(shù);wp為螺旋槳伴流分?jǐn)?shù);yPP、yPS分別為左右側(cè)螺旋槳橫向位置;uWPP、uWPS分別為左右側(cè)螺旋槳處波浪粒子速度;β0P、β0S分別為左右側(cè)螺旋槳有效推力系數(shù),這里設(shè)定為β0P=APP/AP0,β0S=APS/AP0,AP0為螺旋槳面積,APP,APS分別為左右側(cè)螺旋槳實(shí)時浸水面積。

與螺旋槳一樣,左右兩側(cè)舵受波浪粒子速度和實(shí)時出入水情況的影響,產(chǎn)生的舵力不同,同時還會產(chǎn)生首搖力矩,舵力模型在參考文獻(xiàn)[12]提供的舵力模型基礎(chǔ)上修改考慮波浪粒子速度和實(shí)時出入水的影響,各方向舵力的表達(dá)式如下:

式中:tR為舵阻力減額系數(shù);xR為舵的縱向位置;zR為舵的垂向位置;aH為操舵誘導(dǎo)船體橫向力的修正因子;xHR為操舵誘導(dǎo)船體橫向力作用點(diǎn)縱向位置;zHR為操舵誘導(dǎo)船體橫向力作用點(diǎn)垂向位置;yRP、yRS分別為左右舵的橫向位置;ARP、ARS分別為左右舵實(shí)時浸水面積;URP、URS分別為左右舵來流速度;αRP、αRS分別為左右舵攻角;fαP、fαS分別為左右舵力系數(shù);ΛP、ΛS分別為左右舵實(shí)時展弦比;ε為槳舵伴流分?jǐn)?shù)比;η為螺旋槳直徑與舵展比;κ為槳舵相互作用系數(shù);uWPP、uWPS分別為左右舵位置處的波浪粒子速度;AR為舵面積;Λ為舵展弦比;hRP、hRS分別為左右舵相對于波面的瞬時垂向高度。

1.6 波浪力

波浪力主要包括Froude-Krylov 力和繞射力,其中縱蕩波浪力的繞射效應(yīng)通過IMO 提案SDC 2/INF.10附件35[13]中日本提出的經(jīng)驗(yàn)修正因子μx進(jìn)行修正。各方向波浪力的計(jì)算公式如下所示:

1.7 非線性橫搖阻尼

橫搖阻尼是預(yù)報橫搖運(yùn)動的關(guān)鍵,尤其大幅橫搖運(yùn)動,船舶橫甩時通常伴有大幅橫搖,甚至傾覆,本文研究中采用了線性和立方項(xiàng)阻尼系數(shù),公式如下:

式中,α和γ分別為線性項(xiàng)和立方項(xiàng)阻尼系數(shù)。

2 模型介紹

本文選取ONR 內(nèi)傾船開展騎浪數(shù)值模擬研究。ONR 內(nèi)傾船是典型的非常規(guī)高速船,具有細(xì)長體、穿浪船艏、水線以上舷側(cè)內(nèi)傾,以及雙槳雙舵等特性。由于舷側(cè)內(nèi)傾,ONR內(nèi)傾船的橫搖恢復(fù)力矩小于舷側(cè)垂直或外飄的相似船型,因此該船在惡劣海況中高速航行時容易發(fā)生穩(wěn)性失效。ONR 內(nèi)傾船的主要參數(shù)和橫剖面圖分別如表1和圖2所示。

表1 船型主要參數(shù)Tab.1 Main parameters of the ship

圖2 船體橫剖面圖Fig.2 Body plan of the ship

3 結(jié)果與分析

基于以上建立的騎浪/橫甩數(shù)學(xué)模型,分別分析舵和螺旋槳出水、縱蕩繞射效應(yīng)和非線性水動力導(dǎo)數(shù)對內(nèi)傾船騎浪運(yùn)動的影響。試驗(yàn)結(jié)果來自參考文獻(xiàn)[3]。

3.1 舵和螺旋槳出水

計(jì)算工況如下:波長與船長比λ/LPP為1.25;波陡H/LPP為0.05;航向?yàn)?°、15°、22.5°、30°和37.5°;Fn從0.25 到0.5,間隔為0.01。分別計(jì)算不考慮舵和螺旋槳出水和考慮舵和螺旋槳出水兩種情況,數(shù)值計(jì)算的內(nèi)傾船周期運(yùn)動區(qū)域和騎浪運(yùn)動區(qū)域劃分如圖3和圖4中的左圖所示,計(jì)算獲得的騎浪邊界與試驗(yàn)結(jié)果對比如圖3和圖4中的右圖所示。

圖3 船舶運(yùn)動模式數(shù)值計(jì)算結(jié)果與試驗(yàn)結(jié)果[3]對比(不考慮舵和螺旋槳出水)Fig.3 Comparison of boundaries of ship motion modes between numerical results and experimental results without rudders and propellers emersion

圖4 船舶運(yùn)動模式數(shù)值計(jì)算結(jié)果與試驗(yàn)結(jié)果[3]對比(考慮舵和螺旋槳出水)Fig.4 Comparison of boundaries of ship motion modes between numerical results and experimental results with rudders and propellers emersion

從圖3 和圖4 中可以看出,不考慮舵和螺旋槳出水時,固定航向下內(nèi)傾船隨著Fn的增大,船舶運(yùn)動從周期運(yùn)動過渡到騎浪運(yùn)動再到周期運(yùn)動,而考慮舵和螺旋槳出水后計(jì)算的騎浪區(qū)域擴(kuò)大;同時考慮舵和螺旋槳出水時,隨著航向的增大,內(nèi)傾船發(fā)生騎浪的臨界弗勞德數(shù)Fncr增大。通過與試驗(yàn)結(jié)果對比可以看出,除了試驗(yàn)中橫甩情況,考慮舵和螺旋槳出水計(jì)算的周期運(yùn)動和騎浪運(yùn)動區(qū)域與試驗(yàn)結(jié)果更為一致。而試驗(yàn)中內(nèi)傾船發(fā)生橫甩的情況,數(shù)值計(jì)算中只發(fā)生了騎浪,通常騎浪被看作是發(fā)生橫甩的先兆,而橫甩是騎浪平衡的危險分岔,數(shù)值計(jì)算中內(nèi)傾船發(fā)生騎浪后達(dá)到穩(wěn)定,舵和螺旋槳出水等非線性因素產(chǎn)生的擾動沒有打破內(nèi)傾船的騎浪平衡,還需要進(jìn)一步研究。同時橫甩運(yùn)動與船舶初始狀態(tài)息息相關(guān),數(shù)值計(jì)算中未設(shè)置初始擾動,而試驗(yàn)中無法排除初始擾動。

選取的航向?yàn)?°時,不考慮舵和螺旋槳出水時縱蕩速度u的時間歷程曲線如圖5所示。當(dāng)內(nèi)傾船的Fn為0.25 時,內(nèi)傾船做穩(wěn)定的周期運(yùn)動,穩(wěn)定狀態(tài)的縱蕩速度u呈規(guī)則的正余弦。隨著航速增大,內(nèi)傾船發(fā)生騎浪,如圖5中Fn為0.35的情況,內(nèi)傾船的航速增大到波速,穩(wěn)定狀態(tài)的縱蕩速度u不變。當(dāng)航速繼續(xù)增大時,內(nèi)傾船很快加速到大于波速,做穩(wěn)定的周期運(yùn)動,如圖中Fn為0.45 的情況所示,穩(wěn)定狀態(tài)的縱蕩速度u呈規(guī)則的正余弦。

圖5 縱蕩速度u時間歷程(不考慮舵和螺旋槳出水)Fig.5 Time histories of u without rudders and propellers emersion

圖6 為航向?yàn)?°時考慮舵和螺旋槳出水后內(nèi)傾船的縱蕩速度u的時間歷程曲線。考慮舵和螺旋槳出水后,內(nèi)傾船發(fā)生騎浪的Fn范圍增大。以Fn為0.45 的工況為例,考慮舵和螺旋槳出水后,螺旋槳推力效率降低,舵效降低,在縱蕩波浪力的作用下,內(nèi)傾船被加速到波速發(fā)生騎浪。圖7 為該工況下內(nèi)傾船左右舵實(shí)時浸水面積ARP、ARS和螺旋槳有效推力系數(shù)β0P、β0S的時間歷程曲線,從圖中可以看出內(nèi)傾船發(fā)生騎浪后,左右舵的浸水面積ARP、ARS由原來的24 m2分別降到7.8 m2和8.3 m2,左右螺旋槳有效推力系數(shù)β0P、β0S由1.0分別降到0.66和0.69,舵和螺旋槳效率均降低很多。

圖6 縱蕩速度u時間歷程(考慮舵和螺旋槳出水)Fig.6 Time histories of u with rudders and propellers emersion

圖7 騎浪運(yùn)動舵和螺旋槳出水動態(tài)時間歷程Fig.7 Time histories of rudders and propellers emersion of surf-riding motion

3.2 縱蕩繞射影響

分析波浪繞射效應(yīng)對船舶運(yùn)動模式的影響,由于切片法無法計(jì)算縱蕩方向的波浪繞射力,這里縱蕩繞射效應(yīng)采用經(jīng)驗(yàn)公式(32)對縱蕩Froude-Krylov力進(jìn)行修正。分別計(jì)算考慮縱蕩方向的繞射效應(yīng)和不考慮縱蕩方向的繞射效應(yīng)兩種情況,獲得的內(nèi)傾船周期運(yùn)動和騎浪運(yùn)動區(qū)域的邊界如圖8所示。

圖8 考慮縱蕩繞射和不考慮縱蕩繞射計(jì)算的船舶運(yùn)動模式邊界對比Fig.8 Boundaries of ship motion modes with and without surge diffraction

從圖中可以看出,考慮縱蕩繞射效應(yīng)計(jì)算結(jié)果與試驗(yàn)值吻合得更好。因?yàn)榭v蕩方向波浪力只考慮Froude-Krylov 力時,波浪力估算值偏大,使得內(nèi)傾船提前達(dá)到騎浪平衡狀態(tài)。而且航向角為5°時,不考慮縱蕩繞射效應(yīng)所計(jì)算的騎浪邊界對應(yīng)的Fn小于0.3,與騎浪橫甩第一層薄弱性衡準(zhǔn)Fn≤0.3 不易發(fā)生騎浪橫甩的結(jié)論不一致。因此,縱蕩方向的繞射效應(yīng)對船舶騎浪橫甩運(yùn)動影響較大,直接數(shù)值模擬時不可忽略。

3.3 非線性水動力導(dǎo)數(shù)

分析非線性水動力導(dǎo)數(shù)X'vv、X'vr、X'rr、Y 'vvv、Y 'vvr、Y 'vrr、Y 'rrr、N'vvv、N'vvr、N'vrr、N'rrr對船舶運(yùn)動的影響,分別計(jì)算考慮非線性水動力導(dǎo)數(shù)和不考慮非線性水動力導(dǎo)數(shù)兩種情況下內(nèi)傾船運(yùn)動模式,計(jì)算結(jié)果如圖9 所示,獲得的騎浪運(yùn)動和周期運(yùn)動的邊界重合,說明非線性水動力導(dǎo)數(shù)對內(nèi)傾船騎浪運(yùn)動預(yù)報影響較小,在缺乏試驗(yàn)數(shù)據(jù)時可忽略。

圖9 考慮非線性水動力導(dǎo)數(shù)和不考慮非線性水動力導(dǎo)數(shù)計(jì)算的船舶運(yùn)動模式邊界對比Fig.9 Boundaries of ship motion modes with and without nonlinear hydrodynamic derivatives

4 結(jié) 論

本文以雙槳雙舵的ONR 內(nèi)傾船為例,開展了非常規(guī)船型的騎浪數(shù)值預(yù)報方法研究,得出了以下結(jié)論:

(1)本文構(gòu)建的四自由度縱蕩-橫蕩-首搖-橫搖運(yùn)動數(shù)學(xué)模型,可以預(yù)報隨浪和尾斜浪中雙槳雙舵非常規(guī)內(nèi)傾船的騎浪運(yùn)動,用于騎浪直接穩(wěn)性評估。

(2)舵和螺旋槳出入水是影響非常規(guī)內(nèi)傾船騎浪運(yùn)動的一個關(guān)鍵因素,進(jìn)行數(shù)值預(yù)報時不可忽略。

(3)縱蕩方向的繞射效應(yīng)對非常規(guī)內(nèi)傾船騎浪運(yùn)動影響較大,進(jìn)行數(shù)值預(yù)報時不能忽略。

(4)非線性水動力導(dǎo)數(shù)對非常規(guī)內(nèi)傾船騎浪運(yùn)動預(yù)報影響較小,在缺乏試驗(yàn)數(shù)據(jù)時可忽略。

該力學(xué)模型需要進(jìn)一步考慮橫甩的潛在影響因素,如大幅橫搖對雙槳雙舵出水的影響等,為雙槳雙舵船型的橫甩直接評估提供可靠力學(xué)模型。

猜你喜歡
船舶
船舶避碰路徑模糊控制系統(tǒng)
計(jì)算流體力學(xué)在船舶操縱運(yùn)動仿真中的應(yīng)用
CM節(jié)點(diǎn)控制在船舶上的應(yīng)用
基于改進(jìn)譜分析法的船舶疲勞強(qiáng)度直接計(jì)算
《船舶》2022 年度征訂啟事
船舶(2021年4期)2021-09-07 17:32:22
船舶!請加速
BOG壓縮機(jī)在小型LNG船舶上的應(yīng)用
船舶 揚(yáng)帆奮起
軍工文化(2017年12期)2017-07-17 06:08:06
船舶壓載水管理系統(tǒng)
中國船檢(2017年3期)2017-05-18 11:33:09
小型船舶艉軸架設(shè)計(jì)
船海工程(2015年4期)2016-01-05 15:53:30
主站蜘蛛池模板: 欧美日韩激情在线| 91在线丝袜| 大香网伊人久久综合网2020| 自拍欧美亚洲| 日日拍夜夜嗷嗷叫国产| 日韩一区精品视频一区二区| 精品国产一区二区三区在线观看| 精品无码国产自产野外拍在线| 久久黄色小视频| 久久人搡人人玩人妻精品一| 亚洲综合精品第一页| 国产福利影院在线观看| av一区二区三区高清久久| 日韩二区三区| 国产黑丝一区| 91免费国产在线观看尤物| 噜噜噜久久| 日韩毛片免费| 不卡午夜视频| 亚州AV秘 一区二区三区| 中文字幕天无码久久精品视频免费| 欧美一区二区自偷自拍视频| 在线日本国产成人免费的| 一区二区三区精品视频在线观看| 99精品欧美一区| jizz在线免费播放| 丁香六月综合网| 国产真实乱子伦视频播放| 欧美第一页在线| 国产精品浪潮Av| 国产成人高清精品免费| 国内精品久久人妻无码大片高| 亚洲VA中文字幕| 亚洲熟妇AV日韩熟妇在线| 青青草原偷拍视频| 九九热免费在线视频| 成人韩免费网站| 日韩欧美一区在线观看| 在线五月婷婷| 色精品视频| 国产精品午夜福利麻豆| 91精品国产福利| 欧美视频在线播放观看免费福利资源 | 人妻无码中文字幕一区二区三区| 国产美女91视频| 国产日韩欧美一区二区三区在线| 亚洲婷婷丁香| 久草国产在线观看| 嫩草影院在线观看精品视频| 国产日韩精品一区在线不卡| 亚洲一级毛片在线观播放| 一级成人a毛片免费播放| 国产喷水视频| 国产日本欧美在线观看| 国产视频欧美| 亚洲成人精品在线| 国产激爽大片在线播放| av在线无码浏览| 日韩成人高清无码| 亚洲bt欧美bt精品| 亚洲不卡无码av中文字幕| 日韩美一区二区| 欧美翘臀一区二区三区| 波多野结衣在线se| 一本色道久久88亚洲综合| 中文字幕在线播放不卡| 久久窝窝国产精品午夜看片| 国产成人无码Av在线播放无广告| 特级精品毛片免费观看| 在线观看91香蕉国产免费| 亚洲精品无码成人片在线观看| 国产精品污污在线观看网站| 在线播放国产99re| 99热这里只有精品免费| 凹凸国产熟女精品视频| 亚洲天堂精品视频| 国产成人亚洲毛片| 日韩黄色大片免费看| 影音先锋丝袜制服| 国产麻豆精品久久一二三| 一区二区三区在线不卡免费| 亚洲精品成人片在线观看|