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

基于模型試驗的船-冰碰撞載荷空間分布演變歷程研究

2021-03-17 05:55:12孫劍橋
振動與沖擊 2021年5期
關鍵詞:區域

孫劍橋,黃 焱

(1.天津大學 水利工程仿真與安全國家重點實驗室,天津 300350;2.天津大學 建筑工程學院,天津 300350;3.天津大學 港口與海洋工程天津市重點實驗室,天津 300350)

具有自破冰能力的高冰級極地船舶在極地海域的航道開辟、資源開發及科學考察等活動中扮演著重要的角色,同時也承擔著由船-冰碰撞所帶來的巨大的結構損傷風險。目前,國際船級社協會(IACS)頒發的極地船級規范(URI)是極地船舶結構強度設計遵循的主要標準。該規范以船首與巨型浮冰的碰擦(glancing impact)作為船體外板強度設計的控制載荷情形,并基于理論分析構建了船-冰碰撞載荷的直接計算方法[1-2]。

從便于設計的角度出發,IACS規范對船-冰碰撞載荷進行了簡化處理,將其以作用在一定范圍的均布壓力來表示。然而,已公開發表的實船測試工作表明,船-冰碰撞載荷具有突出的空間分布不均勻性與瞬時多變性。Vuorio等[3]基于波羅的海實船航行測試,指出了局部冰壓力的非均勻分布特征。Ritch等[4]通過破冰船與小型冰山碰撞的現場試驗,展示了單次碰撞事件中冰載荷局部空間分布的非均勻性及其隨時間的演變過程。對于高等級極地船舶的結構設計而言,船體結構塑性響應與損傷累積過程的合理評估,是與冰載荷空間分布動態演變歷程的準確把握緊密相聯的。盡管在已發表的實船測試中試圖對船-冰碰撞載荷的空間分布演變特征進行描述,但由于現場測量的限制,仍無法提供能夠用于評估結構損傷發展的動態載荷模型。由此,依靠冰水池內的物理模型試驗探究船-冰碰撞載荷的動態特征,就成為一種重要的研究手段。

國內對于船-冰碰撞的試驗研究尚處于起步階段。張健等[5]采用自由落體的方式進行了冰體碰撞舷側板架模型的試驗,獲得了冰體對板架結構的碰撞力、板架應力分布及板架變形等結果,并與數值仿真結果進行了對比。閆孟嬌等[6]采用水平沖擊試驗機開展了船體板-冰碰撞試驗,驗證了數值仿真模型的有效性,并研究了船體板-浮冰碰撞參數對結構動態響應與塑性變形的影響規律。在本文的先期研究中,黃焱等[7]通過多系列的冰水池物理模型試驗,探索了船-冰碰撞事件的冰水池試驗模擬方法,并發展了碰撞載荷動態特征的測試技術。在此基礎上,黃焱等[8]對船-冰碰撞載荷的時間歷程特征進行了研究,指出了以“單峰”型與“雙峰”型為主要特征的局部冰載荷時程曲線形態。本文所進行的船-冰碰撞載荷空間分布演變歷程研究,可視為上述研究的進一步延伸與補充。

1 模型試驗概述

1.1 試驗設施與相似準則

本文所進行的模型試驗是在天津大學冰力學與冰工程實驗室內完成的。實驗室低溫空間面積為320 m2,用于容納冰水池進行模型試驗。冰水池長40 m,寬6 m,深2 m,可制取的模型冰厚度為1.0~30.0 cm。試驗主拖車的最大水平驅動力為5 t,拖車車速可在0.001~1.0 m/s的范圍內無級調節。

傅汝德和柯西相似準則是冰水池模型試驗中所遵循的主要相似準則,適用于大多數冰與結構相互作用情形。在針對的船-冰碰撞情形中,船體以一定質量和速度撞擊浮冰,即慣性力的作用占主導地位,因此需遵循傅汝德準則;浮冰受到船體撞擊后發生擠壓和彎曲變形,即彈性力的作用占主導地位,因此需遵循柯西準則。上述二種相似準則已寫入國際拖曳水池會議(ITTC)的冰水池試驗技術規程中[9]。

在模型與原型同時滿足傅汝德數和柯西數相等的情況下,可得到模型條件下的冰強度、冰厚、冰彈性模量和壓力的比尺同為模型幾何縮尺比λ,時間和速度的模型縮尺比為λ1/2,質量和力的模型縮尺比為λ3。

1.2 試驗模型與模型冰

試驗中采用的船體模型是根據我國新一代極地科考船“雪龍2”號的船體線型,按照1∶40的幾何縮尺比制作的,船體的主要幾何參數如表1所示。

表1 新型極地科考船的船體主要幾何參數

試驗采用國際第二代低溫模型冰——尿素冰。在模型冰的制備過程中,采用“噴霧引晶”技術模擬天然海冰的初始結晶過程,冰晶格直徑控制在1 mm以下;然后,通過均勻的冷風風場使模型冰生長過程中的室溫維持在-25 °C;最終,當冰蓋生長至預定厚度時停止降溫,并通過回溫調整冰蓋的力學性質。如圖1(a)所示,室內生成的模型冰表層為細密的粒狀結晶層,下層為垂向分布的柱狀結晶層,在紋理結構上與北極地區海冰保持一致,如圖1(b)所示。通過對模型冰生長過程、結晶尺寸及紋理結構的控制,使其在冰的變形與破壞模式、冰載荷特征等關鍵性問題的模擬上與現實情況保持高度的相似性。

(a) 模型冰斷面照片(b) 北極海冰晶體結構圖1 模型冰與天然海冰的對比Fig.1 Comparison of the model ice and natural sea ice

本文涉及的原型船“雪龍2”號的結構強度滿足極地船級PC3級的要求。根據IACS規范,其對應的設計冰況為二年冰,具有較高的冰厚和冰強度。鑒于規范技術背景中將浮冰邊緣的彎曲破壞作為冰體最終的破壞形式,因此試驗中采用彎曲強度作為模型冰強度的指標?,F實情況下,海冰的擠壓強度σc與彎曲強度σf呈現一定的比例關系[10],試驗中制備的模型冰的強度同樣遵循了這一比例關系(σc/σf=3~4)。因此,試驗在對冰彎曲強度進行縮尺時,同時也實現了對冰擠壓強度的縮尺。根據IACS規范技術背景[11],PC3對應的設計冰厚為5.0 m,冰彎曲強度為1.2 MPa。根據傅汝德和柯西相似準則,模型試驗中的目標冰厚為12.5 cm,目標冰彎曲強度為30 kPa。

1.3 試驗測試儀器與碰撞載荷提取

試驗中采用觸覺式傳感器對船體表面的冰壓力進行直接測量。此傳感器由電阻式壓電傳感單元陣列組成,測試區覆蓋船模單側首柱至船肩的整個區域,共包含方形測試單元為1 024個,每個單元尺寸為14.5 mm,量程為1 000 kPa。

試驗中,船-冰碰撞載荷的提取過程可分為:

步驟1 針對每一幀測試數據,消除水壓力的干擾,突出局部冰壓力作用區域;

步驟2 針對每一幀測試數據,進行有效接觸區域邊界的平滑處理,及高壓力區的等值線勾勒(見圖2);

步驟3 針對每次碰撞事件,提取船-冰碰撞力Ft與有效平均壓力Pavg時程曲線;

步驟4 繪制冰載荷作用軌跡(見圖3)。

圖2 觸覺傳感器得到的冰載荷有效作用區域與高壓力區Fig.2 Effective load area and high-pressure zones by tactile sensor

圖3 觸覺傳感器測得的冰載荷沿船模表面的作用軌跡Fig.3 Ice loading trail along the hull surface by tactile sensor

以上船-冰碰撞載荷的提取過程在文獻[7]和[8]中均有詳細的介紹,本文便不再贅述。

1.4 試驗模擬方法與工況

為合理地模擬碰撞能量在船體與冰體間的轉換平衡,在前期研究中進行了一系列的探索性試驗,對船體和冰體在不同約束條件的組合下(船體固定拖曳或自由牽引,冰體自由漂浮或凍結于水池邊壁),碰撞過程的合理性和碰撞載荷的可應用性(即還原至原型后是否能達到規范要求的水平)分別進行了評估,最終確定采用船體固定牽引撞擊自由漂浮冰體的模式進行。在此模式下,船體與拖車進行剛性連接,航速由拖車提供,浮冰則處在自由漂浮的狀態。此模式與IACS規范假定的碰撞模式存在相似之處,即只約束了船體或冰體一方的自由度。

此外,為保證試驗得到的碰撞載荷水平與規范計算結果保持較高的吻合度,在前期研究中,針對不同的初始撞擊位置、浮冰質量和浮冰邊緣角度進行了試驗,最終確定在初始撞擊位置在首柱附近區域、浮冰邊緣角度為90°,且浮冰質量在達到1.5倍~2.5倍的船體排水量后,載荷水平基本達到規范的要求。

試驗中船模的拖曳速度為0.553 m/s,其原型航速為3.5 m/s,與規范技術背景文件中PC3等級下的標準設計航速相對應。由于碰撞載荷具有一定的隨機性,因此本文在上述同等條件下共進行了16組次的試驗,以保證試驗結果的合理性和有效性。

2 船-冰碰撞過程觀測

試驗過程中重點對浮冰的運動狀態的變化,及冰載荷作用位置及水平的變化進行了觀測。以下將結合試驗錄像與觸覺傳感器測試數據,對船-冰碰撞過程進行描述。

圖4~圖6展示的一次碰撞事件中不同時期的試驗場景及濾除水壓力干擾后的冰壓力分布。根據試驗錄像及冰壓力的演變,可以將碰撞事件分為以下三個階段:

(1) 初期——“接觸”階段:浮冰與船體首柱附近區域碰撞后,隨著船體的行進,接觸位置沿水線向后移動;

(2) 中期——“貫入”階段:在船體進一步的行進過程中,由于浮冰邊緣受到首部傾斜線型的引導而發生下壓彎曲變形,接觸面也隨之移動至水線下方的區域。此時,冰壓力也達到了最大水平;

(3) 末期——“分離”階段:隨著船體的繼續前進和浮冰的向外旋轉,浮冰重新上浮至水線處,并開始漂離,最終與船體脫離接觸。

圖4 碰撞初期的試驗場景及冰壓力分布Fig.4 Test scene and ice load at the initial stage of the impact

圖5 碰撞中期的試驗場景及冰壓力分布Fig.5 Test scene and ice load at the middle stage of the impact

圖6 碰撞末期的試驗場景及冰壓力分布Fig.6 Test scene and ice load at the final stage of the impact

上述碰撞過程中的船-冰接觸軌跡可通過圖7進行簡要示意。從圖7可以看出,在碰撞過程中,船體與浮冰保持了充分的接觸,接觸區域基本覆蓋了從首柱至船肩整個首部水線區域。

圖7 碰撞過程中的船-冰接觸軌跡示意Fig.7 Illustration of the ship-ice contact during the impact

3 船-冰碰撞載荷的整體空間移動

在圖4~圖6展示的船體與浮冰的碰擦過程中,由于船-冰之間的相對運動,及冰體的局部擠壓破壞與彎曲變形,導致冰載荷的大小、區域形狀及作用位置是處在不斷變化之中的。據此,在本文的研究中,船-冰碰撞載荷的空間分布演變特征被分為以下兩種:

(1) 載荷沿船體表面的空間移動:這一特征在本文研究中被定義為一種“整體”特征,其主要源于船-冰之間的相對運動及冰體的彎曲變形。

(2) 載荷作用區域的形狀及區域內冰壓力的空間分布:這一特征被定義為“局部”特征,其主要源于冰體的局部擠壓破壞。

在上述特征中,冰載荷的整體空間移動過程可通過圖3所示的作用軌跡進行揭示。通過本文進行的16組試驗得到的冰載荷作用軌跡擬合結果如圖8所示。由于船-冰碰撞過程中冰體下壓彎曲變形的存在,整體冰載荷沿船體外板的作用軌跡,并非簡單地從船首向船肩處沿水線掃略,而是呈現出一種近似拋物線的軌跡。冰載荷通常在冰體下壓彎曲變形最大的時刻,即拋物線形軌跡的最下方,達到其最大值。這種作用軌跡也從側面印證了已發表的船體冰區航行損傷統計中[12],在“冰帶”或冰加強區域下方仍出現了船體受損的情況,如圖9所示。

圖8 試驗得到的冰載荷作用軌跡拋物線擬合(原型尺度)Fig.8 Parabolic fits of the ice loading trails on the hull in full scale

圖9 基于波羅的海冰區航行數據的船體損傷位置統計[12]

4 碰撞載荷的局部空間分布

如前所述,在船-冰碰撞過程中,局部載荷作用區域的形狀及大小通常也處在不斷變化之中。這種特征可歸因于冰擠壓破壞過程中碎冰沫的剝落和擠出現象,如圖10所示。這一過程造成了冰載荷有效作用面積的持續變化,同時由于冰破壞的不均勻性,導致接觸面內若干高壓力區的存在。

圖10 冰與結構的擠壓破壞示意[13]Fig.10 Illustration of the crushing failure of ice

在本系列試驗中,冰載荷局部空間分布的演變過程同樣可通過觸覺傳感器測試數據進行呈現。圖11展示的是單次碰撞事件過程中8個時刻點下的冰壓力局部空間分布情況。在上述8個時刻點中,有5個處于平均冰壓力Pavg時程曲線的峰值,其余3個則處于谷值。根據圖11中Pavg的相對大小,結合前文所述的碰撞過程觀測,可將時程曲線所反映出的冰局部破壞過程與前述三個階段相對應。

圖11 單次碰撞事件中8個時刻點的冰壓力局部空間分布情況Fig.11 Local ice pressures at 8 timepoints during one impact

(1) 初始的接觸與局部擠壓(8.2~8.7 s):此階段中,船-冰接觸面內的碎冰沫剝落過程仍在發展,由此造成了接觸區域的逐漸縮小。同時,此階段的高壓力區分布也較為離散。

(2) 局部擠壓與彎曲變形(8.7~8.9 s):此階段為冰體下壓彎曲變形最大的階段。此時,船-冰接觸區域已縮減至最小,同時冰壓力分布也更為集中。

(3) 最終的分離(8.9~9.2 s):此階段中,冰體的彈性變形開始恢復,致使冰體向外旋轉、漂離。在此過程中,船-冰接觸區域重新向水線附近靠攏,接觸面積有所增大,同時高壓力區的分布也更為離散。

此外,從圖11還可以看出,接觸區域內高壓力區的個數與平均冰壓力Pavg的峰值或谷值存在聯系。在時程曲線的峰值點,通常對應一個高壓力區;而在谷值點,通常伴隨兩個高壓力區。顯然,在每一個“峰值-谷值-峰值”循環中均存在著動態、周期性的碎冰沫剝落與擠出過程。

從整個碰撞事件的角度,圖12展示了“單高壓力區”型與“雙高壓力區”型局部空間分布特征下平均冰壓力Pavg值的統計概率分布。從中可以看出,較大的Pavg值只存在于“單高壓力區”型局部空間分布特征下,而“雙高壓力區”型特征下的Pavg值一般較小。

圖12 “單高壓力區”(實線)與“雙高壓力區”(虛線)所對應平均冰壓力Pavg值的統計概率分布

5 移動冰載荷的構建方法

本文的試驗結果初步揭示了船-冰碰撞過程中冰載荷的“整體”空間移動與“局部”空間分布的演變過程。目前,國際上已發表的研究已初步揭示了移動載荷或非均勻載荷對結構的塑性響應有著顯著的影響[14-15]。因此,本文將基于模型試驗冰載荷測試結果,對上述特征進行構建,以期形成船-冰碰撞載荷空間分布演變歷程的數學描述。

5.1 整體空間移動特征構建

前文第三章通過冰載荷作用軌跡刻畫了船-冰碰撞過程中載荷沿船體表面的空間移動特性。在實際應用過程中,可將冰載荷的移動路徑概化為如圖13所示的軌跡。

圖13 概化冰載荷整體作用軌跡Fig.13 Generalized global ice loading trail

在圖13中,船首區域被劃分為4個等長的子區域。冰載荷軌跡的起點位于首柱區域的后邊界,最大冰載荷的水平位置出現在0.4B(B為船寬),垂向位置則位于水線以下0.5h(h為設計冰厚,本船為5 m)。載荷作用軌跡可延長至船肩,但一般終于船肩子區域中點的下方(見圖8)。

為考慮載荷的整體空間移動,一種簡捷易行的做法為:采用IACS規范的設計載荷板直接覆蓋冰載荷作用軌跡,如圖14所示。載荷板首先應保證對最大載荷作用點的覆蓋,且應互相連接。對于本船而言,規范設計載荷板的長度為3.6 m,高度為0.6 m;采用5個載荷板便可實現載荷軌跡的全覆蓋,其中最大冰載荷對應載荷板#4。

圖14 采用規范設計載荷板實現載荷空間移動的模擬Fig.14 Modeling of the load migration by design load patches

將每個載荷板投影到觸覺傳感器所對應的測試單元與區域,即可獲得該區域的載荷時間歷程。圖15展示的是單次碰撞試驗獲得的5個載荷板的時間歷程。從圖中可以看出,時程曲線多呈現“單峰”形態,且相鄰載荷板之間存在時間上的交叉與重疊。

圖15 單次碰撞試驗獲得的5個載荷板的時間歷程Fig.15 Time histories of the ice forces for the five load patches

據此,可構建載荷板時程曲線的簡化形式(見圖16),其中包含的參數有:i為載荷板編號,i=1,2,…,5;Fi為第i個載荷板的最大冰力;tb,i和te,i分別為載荷板i的冰力起始與終止時刻;TT,i為載荷板i冰力的總作用時間;TR,i為載荷板i冰力的加載時間;Δti為相鄰載荷板冰力的重疊時間。

圖16 載荷板時間歷程的簡化形式Fig.16 Simplified form of the load history on each patch

對于第1個載荷板而言,其冰力初始時刻可設為0,因此僅需獲得每個載荷板的F、TT、TR和Δt四個參數,即可完成對整個碰撞載荷時間歷程的構建。為便于應用,現將上述四個參數做無量綱化處理,即分別轉換為F/FURI、TT/Tspan、TR/TT和Δt/Tspan的形式,其中:FURI為IACS規范計算冰力,對于本船的原型值為15.92 MN;Tspan為載荷歷經單個載荷板的理論時間,通過下式進行計算

(1)

式中:w為載荷板長度;vship為設計航速;α為水線角。對于本船,w和vship的原型值分別為3.6 m和3.5 m/s。

上述無量綱參數的統計值在表2中列出。為消除樣本中個別極端數值的影響,本文在各無量綱參數的最終確定中選取了中位數為代表。據此,可實現5個載荷板上冰力時間歷程的最終構建,如圖17所示。其中,在最大冰力作用時刻4.3 s,載荷板#4將承受與規范設計冰力相當的冰力作用。

表2 無量綱參數的統計結果

5.2 局部空間分布特征構建

為考慮冰載荷局部空間分布的非均勻特性,一種較為簡便的做法即將載荷板進一步劃分為若干個子區域,考察各個子區域間冰壓力的相對大小。在本文研究中,擬將每個載荷板(長3.6 m、高0.6 m)沿長度方向劃分為8個子區域,而沿載荷板高度方向的冰壓力則假定為均勻分布。

前文第四章已初步展示了冰載荷時程曲線波動過程中“單-雙”高壓力區的周期性動態轉換過程,并指出了較大冰載荷只出現在“單高壓力區”型空間分布情形的基本規律。由于在結構強度設計中,設計人員更為關注極端冰載荷情形,因此本文將冰載荷的局部空間分布特征統一簡化為“單高壓力區”形態。這種形態可基于高斯函數進行表達

圖17 五個載荷板上的概化冰力時間歷程Fig.17 Time histories of the ice forces on five load patches

(2)

(3)

式中:x代表載荷板某一子區域的相對位置(見圖18);ls為載荷板子區域至載荷板前邊界的距離;w為載荷板的長度;y為代表載荷板子區域冰壓力相對大小的一無量綱參數,以P/Pmax表示;a,b,c為描述“單高壓力區”分布形態的高斯函數參數,其中,a與曲線峰值的幅值有關,b為峰值的中心位置,c控制峰度;d,e則共同決定函數的基線。

圖18 載荷板子區域的相對位置定義Fig.18 Definition of the local subpatch

將式(2)用于每個載荷板局部空間分布的擬合,擬合結果在表3中列出。圖19以瀑布圖的形式對5個載荷板局部空間分布形態的擬合曲線進行了展示。通過圖19可以看出,5個載荷板局部空間分布曲線的峰值中心與峰度是各不相同的。其中,承受最大冰力作用的載荷板#4所對應的冰壓力分布曲線更為陡峭,說明其冰壓力分布也更為集中。

通過圖19可獲得載荷板各子區域冰壓力的相對大小。而若要得到各子區域的冰壓力絕對值,則需引入圖17中載荷板的總冰力時程曲線,以保證在最大冰力作用時刻,施加在載荷板上的冰力與規范設計冰力在數值上是相等的。在載荷板#4劃分的8個子區域中,最大的子區域冰壓力值為19.7 MPa,約為規范設計平均冰壓力的2.7倍。

表3 冰壓力局部空間分布的高斯函數擬合結果

圖19 五個載荷板上冰壓力局部空間分布形態的擬合Fig.19 Fitted P/Pmax curves for the five load patches

6 結 論

本文針對船-冰碰撞載荷空間分布的演變歷程,通過一系列的冰水池模型試驗,分別從“整體”與“局部”的角度,對船-冰碰撞載荷沿船體表面的空間移動軌跡,及局部冰壓力空間分布形態隨時間的變化過程進行了研究。本文得到的結論如下:

(1) 由于碰撞過程中船-冰相對運動及冰體下壓彎曲變形的存在,冰載荷沿船體外板的整體作用軌跡,并非簡單地從船首向船肩處沿水線掃略,而是呈現出近似拋物線的軌跡。冰載荷通常在冰體下壓彎曲變形最大的時刻,即拋物線軌跡的最下方,達到其最大值。

(2) 冰體擠壓破壞過程中存在著動態、周期性的碎冰沫剝落與擠出過程,由此引發了冰載荷時程曲線中“峰-谷”循環的波動特征。與此相對應,冰壓力的局部空間分布也呈現出“單-雙”高壓力區的周期性動態演變過程。

(3) 極端冰載荷只存在于“單高壓力區”型的局部空間分布形態下。

在此基礎上,本文采用多個載荷板對船-冰接觸軌跡進行覆蓋的理念,提出一種移動冰載荷的構建方法。相鄰載荷板之間的載荷時程存在重疊以更好地模擬載荷的空間“移動”,而各載荷板上的非均布冰壓力則通過高斯函數進行擬合。本文所進行的上述探索性研究,可為規范冰載荷從“靜態均布”向“動態非均布”的擴展提供參考。

猜你喜歡
區域
分割區域
探尋區域創新的密碼
科學(2020年5期)2020-11-26 08:19:22
基于BM3D的復雜紋理區域圖像去噪
軟件(2020年3期)2020-04-20 01:45:18
小區域、大發展
商周刊(2018年15期)2018-07-27 01:41:20
論“戎”的活動區域
敦煌學輯刊(2018年1期)2018-07-09 05:46:42
區域發展篇
區域經濟
關于四色猜想
分區域
公司治理與技術創新:分區域比較
主站蜘蛛池模板: 国产一级α片| 亚洲中文在线看视频一区| 91亚洲精品国产自在现线| 成人一区专区在线观看| 亚洲大尺码专区影院| 久久性妇女精品免费| 国产成人精品亚洲日本对白优播| 91色在线观看| 伊人AV天堂| 成年av福利永久免费观看| 制服无码网站| 免费又爽又刺激高潮网址 | 国产资源免费观看| 国产人成在线观看| 国产欧美日韩一区二区视频在线| 亚洲一区二区无码视频| 久久精品人人做人人爽| 91黄色在线观看| 欧美日韩成人在线观看| 午夜视频免费一区二区在线看| 国产成人欧美| 欧美激情综合| 亚洲最新在线| 国产麻豆精品在线观看| 精品小视频在线观看| 日韩一二三区视频精品| 国产第一页亚洲| 国产成人一区在线播放| 国产一级毛片网站| 不卡无码网| 91国语视频| 思思热精品在线8| 国产美女精品一区二区| 国产永久免费视频m3u8| 国产综合欧美| 呦视频在线一区二区三区| 美女啪啪无遮挡| 中日无码在线观看| 国产精品爽爽va在线无码观看| 欧美日韩北条麻妃一区二区| 精品久久蜜桃| 亚洲av无码牛牛影视在线二区| 国产综合亚洲欧洲区精品无码| 久草视频中文| 亚洲最大看欧美片网站地址| 免费观看欧美性一级| 国产日韩精品一区在线不卡| 久久人搡人人玩人妻精品| 激情无码视频在线看| 91视频精品| 国产成人夜色91| 超薄丝袜足j国产在线视频| 丰满人妻久久中文字幕| 久久青草免费91线频观看不卡| 欧美五月婷婷| 青青草欧美| 国产一国产一有一级毛片视频| 久久精品一卡日本电影 | 欧美日韩国产高清一区二区三区| 精品無碼一區在線觀看 | 韩国福利一区| 亚洲一区二区三区香蕉| 亚洲精品天堂自在久久77| 国产精品自在在线午夜区app| 久久国语对白| 国产成a人片在线播放| 呦视频在线一区二区三区| 国产精品香蕉在线| 欧美笫一页| 免费无遮挡AV| 男女性午夜福利网站| 日韩欧美国产精品| 精品国产成人高清在线| a级毛片免费看| 国产成人精品一区二区免费看京| 超碰aⅴ人人做人人爽欧美| 免费观看亚洲人成网站| 亚洲毛片在线看| 三级毛片在线播放| 国产亚洲欧美在线人成aaaa| 国产丝袜第一页| 日韩精品一区二区三区免费|