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

        ?

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

        2022-08-17 11:24:10儲(chǔ)紀(jì)龍
        船舶力學(xué) 2022年8期
        關(guān)鍵詞:穩(wěn)性螺旋槳波浪

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

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

        0 引 言

        國際海事組織(IMO)于2020年12月發(fā)布了第二代完整穩(wěn)性衡準(zhǔn)臨時(shí)指南[1],目前,正在制定第二代完整穩(wěn)性衡準(zhǔn)直接評(píng)估方法解釋性文件。二代完整穩(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)性評(píng)估的高難度,目前急需一種準(zhǔn)確有效的騎浪/橫甩數(shù)值預(yù)報(bào)方法。所謂騎浪,就是船舶在隨浪或尾斜浪中高速航行時(shí),被波浪捕獲并以波速前進(jìn)的現(xiàn)象。通常船舶在波浪的下坡段發(fā)生騎浪,處于騎浪狀態(tài)下的船舶,多數(shù)會(huì)因航向不穩(wěn)定性而不可控制地轉(zhuǎn)向,發(fā)生橫甩,橫甩產(chǎn)生的離心力嚴(yán)重時(shí)可以導(dǎo)致船舶傾覆。

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

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

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

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

        1.1 船舶運(yùn)動(dòng)坐標(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)處初始波谷的縱向相對(duì)位置。

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

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

        船舶發(fā)生騎浪/橫甩時(shí)船舶遭遇頻率遠(yuǎn)小于升沉和縱搖的固有頻率,升沉和縱搖運(yùn)動(dòng)可通過靜平衡法求解穩(wěn)定平衡點(diǎn)來近似[11]。本文基于MMG 模型構(gòu)建縱蕩-橫蕩-首搖-橫搖四自由度數(shù)學(xué)模型進(jìn)行騎浪運(yùn)動(dòng)數(shù)值預(yù)報(bào),該數(shù)學(xué)模型滿足IMO 提案中對(duì)騎浪/橫甩直接穩(wěn)性評(píng)估的相關(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 提案中對(duì)騎浪/橫甩直接穩(wěn)性評(píng)估的相關(guān)要求指出,騎浪/橫甩數(shù)值模擬方法應(yīng)該正確地模擬由于船體漩渦脫落而產(chǎn)生的水動(dòng)力,也就是除了靜水操縱力,還應(yīng)包括因波浪粒子速度和船舶速度共存產(chǎn)生的水動(dòng)升力。因此,本文在雙槳雙舵的螺旋槳推力模型和舵力模型里考慮波浪粒子速度的影響,同時(shí)還進(jìn)一步考慮了螺旋槳和舵實(shí)時(shí)出入水的影響。

        雙螺旋槳推進(jìn)時(shí),左右兩側(cè)螺旋槳受波浪粒子速度和實(shí)時(shí)出入水情況的影響,產(chǎn)生的推力不同,同時(shí)還會(huì)產(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í)時(shí)浸水面積。

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

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

        1.6 波浪力

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

        1.7 非線性橫搖阻尼

        橫搖阻尼是預(yù)報(bào)橫搖運(yùn)動(dòng)的關(guān)鍵,尤其大幅橫搖運(yùn)動(dòng),船舶橫甩時(shí)通常伴有大幅橫搖,甚至傾覆,本文研究中采用了線性和立方項(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è)垂直或外飄的相似船型,因此該船在惡劣海況中高速航行時(shí)容易發(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òng)力導(dǎo)數(shù)對(duì)內(nèi)傾船騎浪運(yùn)動(dòng)的影響。試驗(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)動(dòng)區(qū)域和騎浪運(yùn)動(dòng)區(qū)域劃分如圖3和圖4中的左圖所示,計(jì)算獲得的騎浪邊界與試驗(yàn)結(jié)果對(duì)比如圖3和圖4中的右圖所示。

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

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

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

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

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

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

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

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

        3.2 縱蕩繞射影響

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

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

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

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

        分析非線性水動(dòng)力導(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對(duì)船舶運(yùn)動(dòng)的影響,分別計(jì)算考慮非線性水動(dòng)力導(dǎo)數(shù)和不考慮非線性水動(dòng)力導(dǎo)數(shù)兩種情況下內(nèi)傾船運(yùn)動(dòng)模式,計(jì)算結(jié)果如圖9 所示,獲得的騎浪運(yùn)動(dòng)和周期運(yùn)動(dòng)的邊界重合,說明非線性水動(dòng)力導(dǎo)數(shù)對(duì)內(nèi)傾船騎浪運(yùn)動(dòng)預(yù)報(bào)影響較小,在缺乏試驗(yàn)數(shù)據(jù)時(shí)可忽略。

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

        4 結(jié) 論

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

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

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

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

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

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

        猜你喜歡
        穩(wěn)性螺旋槳波浪
        波浪谷和波浪巖
        船舶穩(wěn)性控制系統(tǒng)研究
        基于CFD的螺旋槳拉力確定方法
        波浪谷隨想
        去看神奇波浪谷
        3800DWT加油船螺旋槳諧鳴分析及消除方法
        廣東造船(2015年6期)2015-02-27 10:52:46
        絞吸式挖泥船的穩(wěn)性計(jì)算
        廣東造船(2015年6期)2015-02-27 10:52:45
        螺旋槳轂帽鰭節(jié)能性能的數(shù)值模擬
        波浪中并靠兩船相對(duì)運(yùn)動(dòng)的短時(shí)預(yù)報(bào)
        中國航海(2014年1期)2014-05-09 07:54:24
        我的三大絕招
        高清国产日韩欧美| 日韩一区二区三区无码影院| 狠狠色狠狠色综合| 91麻豆国产香蕉久久精品| 久久精品国产av大片| 中文字幕女同人妖熟女| 放荡的少妇2欧美版| 女人做爰高潮呻吟17分钟| 国产精品久久这里只有精品| 国产av精选一区二区| 夹得好湿真拔不出来了动态图| 无套内谢的新婚少妇国语播放| 精品国产品欧美日产在线| 在线观看视频国产一区二区三区 | 国产香蕉一区二区三区| 青青草成人免费在线视频| 欧美日韩精品久久久免费观看| 亚洲影院天堂中文av色| 中文字幕日本人妻一区| 97超碰国产成人在线| 亚洲色欲色欲大片www无码| 亚洲黄色一级毛片| 国产自拍精品视频免费观看| 精品高朝久久久久9999| 国产成年女人特黄特色毛片免| 青青草视频网站免费观看| 水蜜桃网站视频在线观看| 日本高清视频wwww色| 熟妇五十路六十路息与子| 搡老女人老妇女老熟妇69| 日本人妻免费一区二区三区| 少妇无码av无码专区| 亚洲精品123区在线观看| 国产一区二区三区尤物| 精品久久久久久久久午夜福利| 老熟女多次高潮露脸视频| 中文少妇一区二区三区| 久久亚洲精品中文字幕| 国产自偷亚洲精品页65页| 亚洲午夜无码久久久久软件| 国产一区二区三区亚洲avv|