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

聚焦波作用下平面網衣結構的水動力特性研究

2022-10-18 10:13:02俞嘉臻張顯濤
海洋工程 2022年5期

俞嘉臻,張顯濤,李 欣

(1.上海交通大學 海洋工程國家重點實驗室,上海 200240; 2.上海交通大學三亞崖州灣深海科技研究院,海南 三亞 572000)

隨著人類對海鮮需求急劇擴大,網箱養殖在近幾年迅猛發展,然而由于近海的地域限制以及魚類養殖對環境的污染,將網箱養殖轉移到深遠海,尋求更加優越的養殖環境,成為網箱養殖發展的必然趨勢[1]。但是深遠海有著更為惡劣、開放的環境,特別是可能出現的極端波浪,對網箱造成巨大的載荷,給養殖裝備的結構安全帶來全新的挑戰。網衣作為養殖裝備的關鍵結構之一,其水動力特性直接影響了養殖裝備的動力響應。因此,研究極端波浪與網衣結構的相互作用對于養殖裝備的設計與性能評估有重要的意義。

過去十年間,大多數研究聚焦于流作用下的網衣水動力特性,通過試驗和數值模擬來確定流作用下網衣受力與各參數之間的關系。Fridman和Danilov[2]研究發現網衣阻力系數與雷諾數和密實度有關。Milne[3]比較了有結和無結網衣的水動力特性,發現網衣的阻力系數受結節形式的影響。Kristiansen和Faltinsen[4]提出了一種Screen載荷模型來計算定常流下網箱的水動力。Lader等[5]研究了生物附著對網衣阻力的影響,研究發現生物附著效應導致網衣阻力顯著增加。Tang等[6]提出了新的網衣阻力系數計算公式,綜合考慮了網衣密實度、雷諾數、網的材料以及網衣的結節形式等因素的影響。

最近,國內外的研究人員開始研究規則波下網衣結構的水動力特性,并提出了適用于網衣結構的多孔介質模型。Song等[7]研究了柔性網衣在規則波下的受力狀態,發現網衣受到的水平波浪力取決于KC數和網衣參數。Lader等[8]研究了規則波和網衣的相互作用,發現網衣的水平波浪力大約是垂向力的10倍。Patursson等[9]和Zhao等[10-11]應用多孔介質模型分別模擬了在定常流作用下不同來流攻角、多網衣結構以及網箱布置下的流場變化。Bi等[12]和Zhao等[13]進一步將多孔介質模型用于模擬波浪與網衣結構的相互作用。Chen和Christensen[14]基于已有的多孔介質模型,通過與Morison模型計算的網衣受力等效分析,獲得多孔介質阻力系數的直接估計方法,通過將數值模擬結果與試驗結果對比,證明該模型具有一定準確性。Dong等[15]對規則波中不同類型的網衣進行了一系列試驗,并提出了一種基于Morison公式的經驗模型用來計算規則波下的網衣受力。

上述研究集中在研究定常流以及規則波下網衣結構的水動力特性,并建立了多孔介質模型和經驗模型來求得網衣的受力。然而對于極端波浪作用下網衣結構水動力特性的研究目前還十分有限。基于waves2Foam建立數值波浪水池,極端波浪采用基于NewWave理論的聚焦波模擬,網衣結構采用多孔介質模型模擬,通過與Morison模型計算的網衣受力等效分析,獲得多孔介質模擬網衣結構的阻力系數的直接估計方法。開展聚焦波與網衣結構相互作用的數值模擬,與網衣水槽試驗結果作對比,驗證此數值模擬的準確性。在此基礎上,考慮不同的波浪參數和網衣密實度,系統研究聚焦波作用下網衣結構的水動力特性。

1 數學模型

1.1 聚焦波理論

Tromans等[16]提出了NewWave理論,用來模擬給定波峰值的聚焦波,該方法通過預定的波峰高度,分配每個頻率的波浪譜能量,最后通過式(1)求得各組成波的波幅。

(1)

式中:ai為聚焦波各組分波浪波幅,Ac為聚焦波最大波幅,S(fi)為組成波對應的波浪譜(本文采用Jonswap波浪譜),δf為能量譜的頻率間隔,N為組成波數目。則基于NewWave理論的聚焦波時歷計算公式為:

(2)

式中:η為波面升高,ki、ωi分別為第i個組成波的波數和圓頻率,Xt、Tt分別為聚焦波的聚焦位置和聚焦時刻。如圖1所示,聚焦波是一種瞬時極端波浪,屬于非平穩過程,現采用局部波浪參數對其進行描述,Ac為靜水面到最大波峰的垂向距離;聚焦波譜峰周期Tp為波浪組分中能量最大的波浪對應的周期;譜峰波長λp通過譜峰周期計算得到;波陡s定義為聚焦波波峰Ac與譜峰波長λp的比值。

圖1 聚焦波參數定義

1.2 多孔介質模型

1.2.1 控制方程

針對多孔介質模型,Jensen等[17]得到了體積平均、雷諾平均的流體運動控制方程:

(3)

(4)

(5)

式中:C是非線性多孔介質阻力系數矩陣。

網衣受到的瞬時力是數值模擬中最重要的輸出部分,式(5)描述的是多孔介質對流體的作用力,而流體對多孔介質的作用力為多孔介質所受阻力即所需要輸出的部分,兩者互為反作用力。多孔介質所受瞬時力Qi最終表達式為:

(6)

式中:PV是多孔介質控制區域。

另外,Van Gent[18]將多孔介質的慣性力系數Cm定義為:

(7)

式中:γp是無量綱經驗系數,一般取值為0.34。

1.2.2 多孔介質阻力系數直接估計方法

獲得網衣水動力系數是網衣水動力研究中關鍵的問題之一。一般而言,將網衣水動力系數分為阻力系數(順流方向)和升力系數(垂直于來流方向),也可以分為法向力(垂直于網衣平面)和切向力(位于網衣平面內)。當來流方向垂直于網衣平面時,網衣阻力與其法向力保持一致。文中主要關心網衣的阻力系數。由于數值模擬中采用多孔介質模擬網衣結構,因此,多孔介質的阻力即為所關心的網衣結構的阻力。

圖2為網衣結構示意,網衣由若干網目組成,網目根據形狀可劃分為菱形、矩形和六邊形結構,根據結節可以分為有結節和無結節兩種情況。每個網目由目腳和結節構成,對于菱形和矩形網目,一般認為單個網目由4根目腳和4個結節組成。文中用到的網衣均為無結節矩形網衣。在網衣水動力計算中涉及到一個重要的參數即密實度Sn,代表的意義為網線投影面積與網衣整體投影面積的比值。

圖2 網衣結構

(8)

式中:dw為網線直徑,lw為網衣目腳長度。

通過與Morison載荷模型等效來得到多孔介質非線性阻力系數(即式(5)中的Cij)。假設網衣處于定常流中,且為無結節方形網衣。以網衣平面法線方向為x方向建立坐標系,網衣固定在y-z平面內,因此無需對系數進行坐標變換。水流可分解為x-y平面和z軸的分量,圖3即為速度矢量的分解,θ、β表示來流速度與網衣平面之間的夾角,網衣的傾斜角通過這兩個參數體現。

圖3 水流速度矢量分解和單個網目網線簡化示意

Tauti[19]在分析網衣水動力時提出了假設:網衣在水中運動時,作用在每個網目的目腳和結節上的水動力彼此獨立,互不相關,因此認為作用在網衣上的總力是每個網線所受力的疊加, 單個網目有水平、豎直兩類網線(如網線1和網線2),網目在x方向的受力Fx(即沿網衣平面法線方向上的受力)是網線1和網線2在x方向受力的疊加,即Fx=Fx,1+Fx,2。

(9)

(10)

(11)

式中:Fx,Fy,Fz分別為網衣在x,y,z方向上的受力,M是平行于網線2的網線總數,N是平行于網線1的網線總數,A1,A2是網線1和網線2的投影面積,Cd,twine是網線的阻力系數。基于式(6)作用在多孔介質上的力可表示為:

(12)

(13)

(14)

式中:Qx,Qy,Qz分別為網衣在x,y,z方向上所受流體的瞬時力(在下文中,Qx即為網衣阻力Fd,Qz即為網衣升力Fl),V為多孔介質區域體積,U∞為無窮遠處來流速度。由于將網衣看作多孔介質,所以作用在多孔介質上的力等于每個網線受力的總和,即F=Q,計算得到C1,C2,C3:

(15)

(16)

(17)

此外,通過正確地降低Morison模型中系數的夾角依賴性,仍然可以得到合理的結果。因此文中引入a、b兩個系數來描述網線與網線之間的相互作用效應,并將角度從公式中去掉,以減少多孔介質的角度依賴性,其中的影響已經包含在a、b兩個系數中。阻力系數最終表達式為:

(18)

(19)

(20)

式中:S1,S2分別為網衣內豎直網線和水平網線的投影面積。從Rudi等[20]針對不同網衣沖角、密實度的系列試驗中選取試驗數據對系數a、b進行擬合,通過線性插值的方法得到一定密實度下的a、b表達式:

(21)

(22)

式中:Sn為網衣密實度。

2 試驗和數值模型

2.1 物理試驗模型

試驗在波浪水槽中進行,如圖4(a)所示,該波浪水槽總長18 m,其中過渡區3 m,工作區12 m,消波區3 m;水槽寬1 m,水深可調節,試驗中水深為0.86 m。造波機搖板精度誤差在0.01°以下,搖擺最大幅值可達±25°,最小造波周期可達0.3 s,波高最大可達20 cm,完全滿足試驗的精度要求。

圖4 物理試驗模型

圖4(b)為網衣試驗裝置,由網衣、夾具、細鐵桿、頂部和底部固定裝置及三分力傳感器組成。傳感器固定于夾具上,上下夾具分別和頂部鋁型材、底部壓鐵固定連接,細鐵桿固定于兩個夾具上,網衣通過上下兩根細鐵桿與夾具連接并展開,通過調節頂部鋁型材上直角鋼的高度來控制網衣豎直方向上的張緊程度。其中三分力傳感器量程為50 N,精度0.01 N,滿足試驗的精度要求。

試驗利用電阻式浪高儀測量波浪升高,利用三分力傳感器來測量網衣在波浪作用下受到的阻力和升力。試驗中共用到兩個三分力傳感器,量程50 N,精度0.05%,安裝在網衣裝置的頂部和底部。當網衣受到波浪力時,傳感器在z方向和x方向受力,數據采集系統采集傳感器數據得到實時的網衣受力。此處說明在試驗安裝過程中,網衣的阻力方向為傳感器的z方向,網衣的升力方向為傳感器的x方向。在試驗前需對傳感器進行標定,分別得到三分力傳感器在3個方向上的標定系數。

2.2 試驗工況

模型試驗測試了網衣在不同聚焦波作用下的水動力結果。具體的網衣參數及波浪工況如表1和表2所示。

表1 網衣參數

表2 波浪工況

在試驗開始之前,首先進行了空白對照試驗,即把網衣拆除。由于整個裝置只有底部傳感器單元浸沒于水中,所以對底部的三分力傳感器的受力進行監測。對在若干個波浪下的數據進行采集,發現三分力傳感器在無網衣的情況下基本不受力,其幅值與無外力作用時的幅值接近,所以認為除網衣結構之外的部件對傳感器的受力影響可以忽略不計。

每個工況重復3次以排除試驗的偶然性并驗證試驗結果的可靠性。如圖5所示,對比了Exp_Wave2工況下3次試驗的結果。從圖5中可以看出,3次試驗的結果吻合度非常好,說明在波浪作用下該試驗裝置測得的數據重復性好,可靠性高。在下文與數值模擬結果對比時,將3次重復試驗結果的峰值做對比,取中間值的那次試驗作為對比數據。

圖5 Exp_Wave2工況重復性試驗結果對比

2.3 數值模型及準確性驗證

在數值模擬計算中,基于式(3)~(22)開發適用于計算網衣結構水動力的多孔介質模型求解器,在原有的造波工具waves2Foam中嵌入多孔介質區域。相關的參數輸入包括聚焦波的波浪參數:各波浪組分的頻率、波幅、波數、相位等;以及網衣的物理參數Sn、Cm、V、Cd,twine,系數a,b以及網衣中水平、豎直網線的總投影面積S1、S2。在數值計算過程中,多孔介質區域外的有限體積單元采用waves2Foam默認求解器求解波浪運動。多孔介質區域內的有限體積單元帶有porosity標識,當讀取到該單元時會切換到多孔介質求解器進行求解,并將每個時間步長的計算結果進行保存。為得到想要的網衣受力結果,對多孔介質區域內的有限體積單元計算得到的阻力和升力進行積分并實時輸出到結果文檔,得到實時的網衣受力結果,具體計算過程如圖6所示,多孔介質區域內的有限體積元e1,e2……en,每個有限體積元處的速度場都已知,根據式(6)得到每個有限體積元的瞬時力,最終將多孔介質區域的所有有限體積元進行積分得到整個多孔介質區域即網衣的瞬時力。為得到波浪在網衣結構下的演化特征,在數值水槽若干處設置“浪高儀”,浪高的計算原理是基于局部水體積分數α(在流體中α=1,在空氣中α=0)積分求得該監測點的水深。下文中的結果都基于上述的輸出文檔數據進行分析。

圖6 數值模擬結果輸出

所建數值波浪水池參照實際水槽尺寸設計,基于waves2Foam來開展聚焦波與網衣結構相互作用的數值模擬。如圖7所示,水槽總長14 m,高1.1 m,水深0.86 m,消波區域長2 m,多孔介質區域厚50 mm,多孔介質區域中心距離速度入口7 m,聚焦波聚焦位置位于多孔介質區域中心后方50 mm處,即距離速度入口7.05 m(避免多孔介質區域內部對浪高儀監測的影響)。

圖7 數值波浪水槽

試驗中網衣的參數及波浪工況同實際水槽試驗,如表1和表2所示,其中聚焦波聚焦位置距造波板7 m。

網格密度按最大波高和波長進行劃分,對波浪自由面的網格進行加密,對消波區等區域的網格進行加疏。文中對3種網格類型進行分析,每種網格的具體劃分精度見表3。用這3類網格密度對工況Exp_Wave1進行數值模擬,驗證其收斂性,數值模擬結果如圖8所示。

表3 網格劃分精度 Tab.3 Mesh scheme 單位:cm

圖8 不同網格的模擬對比

由圖8可以得到,隨著網格精度的提高,Mesh1和Mesh2的結果有一定差距,但是隨著網格繼續細化,Mesh3精度并沒有顯著提高,所以Mesh2滿足數值模擬的精度要求,在下文中就以Mesh2的網格精度對聚焦波—網衣結構進行數值模擬。

在數值模擬中,需要與試驗一致的參數有密實度Sn(與孔隙率關系為n=1-Sn),以及在各波浪下網線對應的阻力系數Cd,twine,根據試驗的已有數據,數值模擬多孔介質區域所需的參數見表4和表5。

表4 多孔介質阻力系數

表5 多孔介質參數

圖9至圖12為4個波浪工況下數值模擬和試驗結果對比,從中可以得到,數值模擬和試驗結果所得數據趨勢相同且結果相似,但是可以發現,數值模擬得到的網衣阻力總是大于試驗阻力結果,而升力總是小于試驗升力結果。這是因為在實際試驗中,網衣張緊程度對于網衣受力有影響。在試驗中,當聚焦波到達網衣時,網衣由于波浪力發生變形,與波浪傳播方向形成一定沖角,從而一定程度上減小了阻力,增大了升力。其次,由于網衣材料的可變性,在網衣張緊狀態下其網衣密實度會隨網線直徑減小而減小。但是在數值模擬過程中,網衣一直處于理想的張緊狀態,在波浪下不會發生變形,所以密實度也始終保持一致。

圖9 Exp_Wave1波浪工況下數值模擬與試驗結果對比

圖10 Exp_Wave2波浪工況下數值模擬與試驗結果對比

圖11 Exp_Wave3波浪工況下數值模擬與試驗結果對比

圖12 Exp_Wave4波浪工況下數值模擬與試驗結果對比

綜上,數值模擬結果中網衣阻力偏大且升力偏小是合理的。對比4個波浪條件下的結果,可以驗證該數值模型的有效性,能夠較準確地模擬網衣在聚焦波下的受力情況。

3 結果與討論

3.1 網衣密實度對平面網衣水動力特性的影響

本節通過該數值方法研究分析極端波浪下網衣密實度對平面網衣結構水動力特性的影響。以網衣密實度、網線直徑、目腳長度為變量設計了12種網衣,具體參數如表6。

表6 網衣參數

其中Net1至Net6通過改變目腳長度來控制密實度;Net7至Net12通過改變網線直徑來改變密實度。橫向分析兩組的水動力特性來比較兩種改變密實度的方式對網衣水動力的影響。

在數值模擬中,選取上節的Exp_Wave2聚焦波為模擬工況,表7為具體波浪參數。

表7 聚焦波參數

李玉成等[21]在研究波浪條件下漁網的水動力特性時,分析得到網衣結構在波浪下Morison方程的阻力系數受KC數影響不大,所以在計算網線的阻力系數時,忽略KC數的影響,由此得到多孔介質模型所需的參數,具體如表8所示。

表8 多孔介質模型參數

基于以上獲得的波浪參數以及多孔介質參數,在wave2Foam波浪水槽中進行平面網衣的數值模擬,獲得各類網衣的水動力特性結果,現對結果數據進行分析。圖13(a)和13(b)展示了不同網衣密實度下,經過平面網衣結構后的聚焦波波浪時歷,由圖可以得到網衣結構的存在影響聚焦波的時空演化,且波浪波幅峰值隨著網衣密實度的增大而減少,不過總體上衰減幅度不明顯,如圖13(a),當密實度Sn從0.05增至0.30,波幅峰值僅從5.04 cm降至4.34 cm,減少了14%。

圖13 Net1~Net12的水動力特性

圖13(c)和13(d)及圖13(e)和13(f)證明網衣密實度對網衣的受力有極為顯著的影響,總體趨勢為網衣受力隨著密實度的增大而明顯升高。縱向比較圖13(c)和圖13(e),發現阻力和升力時歷上峰值出現時刻不一致,阻力峰值出現時刻與聚焦波聚焦時刻較為吻合,且阻力曲線左右對稱,升力曲線大致為中心對稱。這可以根據Morison公式中升力的計算原理進行解釋,由于網衣升力大小取決于波浪中水粒子豎直方向上的平均速度,當水粒子水平速度達到最大時,網衣阻力時歷的斜率絕對值最小,即此時豎直速度為0;當水粒子水平速度為0時,網衣阻力時歷的斜率絕對值最大,此時豎直速度達到極值。比較圖13(c)和圖13(e)可以證明,升力每個波峰和波谷出現時刻都與阻力為0的時刻相對應,總體上,網衣阻力曲線大致是以聚焦時刻為軸的軸對稱圖形,升力曲線大致是以聚焦時刻為中心點的中心對稱圖形,而且當密實度增大時,聚焦時刻前后升力的極值開始不相等,聚焦時刻后峰值絕對值大于聚焦時刻前的峰值絕對值,說明水粒子豎直速度總是在聚焦波到達之后達到最大值。觀察圖13(e)和圖13(f),發現當網衣密實度增加時,升力時歷在受力周期里出現了多個極值,曲線開始不光順,這是因為聚焦波為強非線性波,當密實度增大,升力整體提高時,高諧波成分對升力的影響逐漸明顯。

為了具體分析阻力、升力和密實度之間的關系,就網衣密實度—受力兩者進行分析,如圖14所示。由圖14可以得到,阻力、升力和密實度之間大致成線性關系,如圖14(a)所示當密實度從0.05增加至0.30,其阻力峰值從0.68 N增值5.14 N,升力峰值從0.15 N增至0.68 N,阻力升力比也從4.5增至7.5。所以阻力較升力對密實度的變化更加敏感。

圖14 密實度—網衣受力曲線

3.2 聚焦波波幅對平面網衣水動力特性的影響

對平面網衣—聚焦波結構中的聚焦波進行研究,以聚焦波的波幅為控制變量來分析不同聚焦波下平面網衣的水動力特性,探究聚焦波波幅與網衣受力之間的數值關系。網衣選用上節的Net3,根據波幅大小設計6種聚焦波,表9為具體參數。

表9 波浪及多孔介質參數

同理,將各參數代入到數值模型中得到結果,現對結果數據進行處理分析。圖15(a)展示了不同聚焦波的波浪時歷,波幅從0.83 cm增加至6.61 cm。由圖15(b)和15(c)可得,聚焦波波幅的大小對網衣受力有極其顯著的影響,總體趨勢為聚焦波越大波幅越大,網衣阻力和升力越大。

圖15 不同聚焦波作用下網衣的水動力特性分析

為探究網衣受力與波幅之間的具體關系,現對網衣受力峰值和譜峰波幅進行分析,具體如圖16所示。由圖16(a)可以發現,阻力、升力和最大波幅之間呈非線性關系,且受力增長速率隨最大波幅的增加而升高;圖16(b)可以得出,阻力與升力之比隨著最大波幅的增加而增大,由此推出,最大波幅對阻力的影響較升力更加明顯。

圖16 波幅峰值—網衣受力曲線和波幅峰值—阻力與升力比曲線

4 結 語

通過網衣試驗與數值模擬結果進行對比,驗證了數值模型計算的準確性和可行性。其中極端波浪數值水槽基于waves2Foam建立,采用多孔介質模型模擬網衣結構,其多孔介質阻力系數通過與Morison模型計算的網衣受力等效分析的直接估計方法獲得。對不同網衣密實度及不同波浪參數下網衣結構的升阻力特性以及網衣結構對波浪場的擾動規律進行分析。得到以下結論:

1)在相同聚焦波下,網衣密實度參數對網衣受力有極為顯著的影響,網衣所受阻力和升力隨著密實度的增大而明顯升高,與密實度大致呈線性關系,且阻力較升力對密實度的變化更加敏感。阻力和升力峰值出現時刻不同,阻力峰值出現在升力斜率最大時刻,升力同理。

2)聚焦波波幅的大小對網衣受力有極其顯著的影響,總體趨勢為聚焦波最大波幅越大,網衣阻力和升力越大,阻力、升力和最大波幅之間呈非線性關系,且受力增長速率隨最大波幅的增加而升高。

3)網衣結構對聚焦波時空演化特征影響相對較小,改變了聚焦波波形,且隨著網衣密實度的增大,波浪波幅峰值逐漸減少,不過總體上衰減幅度不明顯。

主站蜘蛛池模板: 99视频在线观看免费| 国产精女同一区二区三区久| 凹凸国产分类在线观看| 77777亚洲午夜久久多人| 91年精品国产福利线观看久久| 无遮挡一级毛片呦女视频| 欧美国产综合视频| 91青青视频| 欧美成人a∨视频免费观看| 国产电话自拍伊人| 丰满人妻中出白浆| 人妻无码AⅤ中文字| 国产亚洲精品自在久久不卡| 日韩无码视频专区| 97国产精品视频自在拍| 99精品国产自在现线观看| 无码区日韩专区免费系列| 欧美激情视频一区| 在线免费观看AV| 在线中文字幕网| 国产一级二级在线观看| yjizz视频最新网站在线| 在线观看亚洲人成网站| 日韩毛片在线播放| 日韩免费无码人妻系列| 亚洲女同欧美在线| 欧美高清三区| 国产丰满大乳无码免费播放| 国产毛片片精品天天看视频| 日韩在线永久免费播放| 亚洲欧美成人| 欧美日韩国产一级| 国产真实乱子伦视频播放| 国产成人免费观看在线视频| 日本一区二区三区精品AⅤ| 国产无码制服丝袜| 国产精品美女免费视频大全| 日韩精品资源| 激情综合激情| 亚洲第一黄色网址| 激情五月婷婷综合网| 国产精品55夜色66夜色| 午夜不卡视频| 亚洲精品第一页不卡| 精品無碼一區在線觀看 | 99热这里只有精品在线观看| 再看日本中文字幕在线观看| 日韩精品一区二区深田咏美| 国产成人一区免费观看| 日韩不卡高清视频| 日本www色视频| 高潮毛片免费观看| 黄色福利在线| 国产熟睡乱子伦视频网站| 亚洲视频免费在线| 国产人碰人摸人爱免费视频| 亚洲综合网在线观看| 亚洲人成日本在线观看| 亚洲国模精品一区| 一级毛片高清| 九色在线观看视频| 国产99免费视频| 国产精品九九视频| 三上悠亚在线精品二区| 91精品国产麻豆国产自产在线| 日韩欧美国产区| 国产一区二区三区免费| 亚洲午夜天堂| 欧美日韩一区二区三区在线视频| 久久国产热| 54pao国产成人免费视频| 97se亚洲综合在线| 久久精品波多野结衣| 国产精品制服| 免费观看成人久久网免费观看| 伦伦影院精品一区| 国产成人免费观看在线视频| 又爽又大又黄a级毛片在线视频| 狠狠色噜噜狠狠狠狠奇米777| 伊人成人在线视频| 久久综合激情网| 国产精选小视频在线观看|