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

樁靴連續貫入過程的動態模擬方法研究①

2015-06-09 12:36:02郭紹曾彭碧瑤
地震工程學報 2015年2期
關鍵詞:變形影響方法

周 龍, 劉 潤,2, 郭紹曾, 彭碧瑤

(1.天津大學水利安全與仿真國家重點試驗室,天津 300072; 2.港口巖土工程技術交通行業重點實驗室,天津 300072)

樁靴連續貫入過程的動態模擬方法研究①

周 龍1, 劉 潤1,2, 郭紹曾1, 彭碧瑤1

(1.天津大學水利安全與仿真國家重點試驗室,天津 300072; 2.港口巖土工程技術交通行業重點實驗室,天津 300072)

采用CEL方法(耦合的歐拉-拉格朗日分析法)對比研究樁靴下沉速度不同對計算結果的影響,分析不同樁靴下沉速度下附近樁基礎的響應,研究質量放大方法的適用性。結果表明:無樁時樁靴貫入速度不同對土體阻力影響較小,土體的破壞形態和剪應力水平略有不同;樁靴貫入速度對樁身水平位移影響較小,樁靴貫入速度越大,樁身最大應力越小;質量放大系數的增加對樁身最大位移的影響較小,但對樁身最大應力有較大影響,因此建議謹慎使用質量放大方法。

樁靴; 樁基礎; 貫入速度; 質量放大

0 引言

隨著海洋油氣資源的不斷開發,海洋鉆井平臺的應用日益增多。海洋鉆井平臺是海上油氣資源勘察開發的主要設備,其主要作用是為鉆井設備提供一個能夠正常工作的平臺。在所有形式的海洋鉆井平臺中,自升式鉆井平臺因作業時自身穩定性好且浮運方便而獲得了較多應用。自升式鉆井平臺的樁靴貫入深度大,貫入過程中對土體的破壞和擾動劇烈,可能發生穿刺破壞現象導致工程事故發生,同時對臨近結構基礎也可能產生不良影響。因此,準確模擬樁靴連續貫入的整個動態過程可為分析自升式鉆井平臺的安全性提供參考。

國內外學者已針對樁靴貫入過程進行了一系列試驗研究。Hossain等[1-4]利用室內實驗研究了樁靴貫入正常固結黏土和軟黏土時的承載特性和土體破壞模式,并將實驗數據與工程實測數據進行了對比,取得較好一致性。Xie等[5]開展了離心模型試驗對樁靴貫入軟黏土時的土體破壞過程進行模擬,并運用PIV圖像識別技術精確反映了樁靴臨近土體的運動機制,重點分析了樁靴附近的孔壓變化,得到橫向及縱向的土體位移情況。張浦陽等[6]利用流固耦合模型對均質及非均質黏土中樁靴貫入深度與承載力關系進行研究,得到了樁靴上部孔穴高度未完全形成時不同土質對樁靴承載力的影響情況。

由于樁靴貫入過程的復雜性,目前的試驗研究多采用室內小比尺模型試驗,只能得出定性的結論。有限元模擬方法由于本身是基于連續體力學的,無法模擬出土體的破壞擠壓過程,更難達到模擬整個動態過程的目的。因此,眾多學者采用耦合的歐拉-拉格朗日分析法(CEL)對樁靴連續貫入過程進行了模擬。2010年Qiu等[7]采用亞塑性本構關系定義沙土,結合CEL法對樁靴貫入導致的沙土大變形過程進行數值模擬;2013年Zheng等[8-9]利用CEL法分析了三層黏土中樁靴貫入和靜力觸探的探入過程;2014年Haydar等[10]利用CEL法模擬了樁靴貫入過程中土體及臨近導管架樁基的變形響應規律,并與傳統拉格朗日分析法結果進行比對,驗證了CEL法的優越性。

以上研究多只針對自升式平臺樁靴自身的插樁過程,而樁靴貫入位置通常距導管架平臺較近,貫入引起的土體大幅破壞和擾動會對臨近樁基的穩定性產生較大影響。1990年Siciliano等[11]開展離心機實驗,針對軟黏土中樁靴貫入對臨近樁基橫向位移影響展開研究;2009—2010年吳永韌等[12-14]利用模型實驗、abaqus非線性有限元法及ALE數值模擬法得到臨近主平臺基礎承載力變化規律;2013年Tho等[15]采用歐拉有限元模擬對比室內實驗方法得到樁靴貫入與鄰近樁基間作用模式。

目前CEL方法在樁靴連續貫入問題的研究中方興未艾,但在樁靴下沉速度的設定和質量放大方法的適用性等方面尚有很多疑問,需要進一步的研究驗證。本文首先對比研究無樁時樁靴下沉速度不同對計算結果的影響,然后對比分析不同樁靴下沉速度下附近樁基礎的響應,最后研究質量放大方法的適用性。

1 耦合的歐拉-拉格朗日算法(CEL)

1.1 物質變形的描述方法

根據坐標系選取的不同,可以將物質運動的描述分為Lagrange描述(又稱物質描述)和Euler描述(又稱空間描述)兩種方法。

物體在t=0時刻所占據的空間區域稱為初始構形,記為V0。為了度量物體的運動,需要選取一個特定的構形作為基準,稱為參考構型。一般取初始構形作為參考構形,如圖1所示。

圖1 初始構形和現時構形Fig.1 Initial configuration and current configuration

在參考構形中質點的矢徑X為

(1)

式中:ei為直角坐標系的基矢量;Xi為參考構形中質點的矢徑X的分量。

質點的矢徑X不隨時間t變化。Xi稱為物質坐標或拉格朗日坐標,它可以作為該質點的標記。

在現時構形中質點的矢徑x為

(2)

式中:xi為矢徑x的分量,其余符號同式(1)。

坐標x給出了質點在空間中的位置,稱為空間坐標系或歐拉坐標。

在Lagrange描述的方法中,取物質坐標Xi和時間t作為獨立坐標,即借助于運動著的質點來考察物體的運動和變形。在這種描述中質點的位移為

(3)

式中:ui為質點位移的分量;其余符號同式(1)和式(2)。

在Euler描述的方法中空間點在不同時刻被物質點取空間坐標xi和時間t作為獨立坐標,描述同一空間點在不同時刻被物質點占據的情況。在這種描述中,質點的位移為

(4)

式中:符號同式(4)。

1.2 ABAQUS軟件中的CEL算法

傳統的拉格朗日分析中節點與材料相互綁定,材料隨單元的變形而變形。該方法依賴于網格的變形,通過求解單元節點位移得到力學結果,但對于樁靴下沉這類大變形問題會存在極大的單元變形奇異,造成收斂困難,因而無法模擬。純粹的歐拉分析方法主要應用于流體力學分析中,網格節點被固定在空間中,材料在不變形的網格內自由流動。一旦歐拉單元流出歐拉網格之外,它就從分析中徹底流失了。純粹的歐拉分析方法通過求解速度場得到整個流場的流態,但無法求解結構準確的應力-應變響應,亦不能完全滿足模擬樁靴下沉對附近樁基礎影響的要求。

CEL技術是ABAQUS軟件中計算流固耦合的關鍵技術,吸取了歐拉網格和拉格朗日網格的優點,采用歐氏網格中網格固定而材料可在網格中自由流動的方式建立模型,有效地解決了有關大變形、材料破壞和流體材料等諸多問題;同時建立拉式網格與歐氏網格的接觸算法,利用拉式網格得到結構準確的應力-應變響應。

CEL技術應用強大簡易的Abaqus/Explicit通用接觸算法,可以較好地解決流體和固體結構的接觸問題以及穿刺問題。樁靴的下沉問題就是穿刺問題的一個特例:樁靴在外力的作用下從外部逐漸侵入土中,將土體擠密推走,占據原有土體的位置。土看成可以自由流動的流體,采用歐拉網格進行模擬,解決土體的大變形問題;樁靴可以采用剛體模型,樁基礎采用三維拉氏單元進行模擬,建立土與樁靴和樁之間的接觸關系,利用Explicit動力學求解器得到樁靴下沉過程中對附近樁基礎的影響。

1.3 質量放大方法

在ABAQUS有限元軟件中,CEL方法是基于中心差分的動態顯式求解算法計算的。利用中心差分算法求解具體問題時,時間步長Δt必須限定在由該問題求解方程性質所決定的某個穩定極限值Δtcr之內,否則算法不收斂。中心差分法解的穩定條件為:

(5)

(6)

式中:Δtcr為穩定極限值,又稱為臨界步長;Le為穩定單元長度,一般為最小單元尺寸;Cd為材料波速,與材料的彈性模量E和密度ρ有關。

人為地將材料的密度ρ增加f2倍,則波在該材料內的傳播速度會降低f倍,從而將穩定時間增量提高f倍。這樣進行同樣的分析所需要的增量步就會減少,達到提高計算速度的目的。但是質量不能進行無限的放大,否則計算結果將不符合實際情況。

2 模型建立

某自升式鉆井平臺樁靴直徑18m,最大預壓載10 000t,貫入深度為10m。不考慮樁靴本身的變形與內力,將之視為剛體,使用離散剛體單元(R3D4)進行模擬。土體視為可以自由流動的液體,使用歐拉網格單元(EC3D8R)進行模擬,土體模型采用Mohr-Coulomb彈塑性模型,具體參數如表1所示。

表1 土體參數表

Qiu[7]和Zheng等[8-9]利用CEL法分析樁靴貫入問題時將貫入速度設定在米每秒的量級,這與樁靴的實際作業速度相去甚遠。根據工程資料,樁靴的實際貫入速度大概在厘米每秒甚至毫米每秒的量級。為研究人為提高樁靴貫入速度對計算結果的影響,建立在樁靴貫入深度相同時貫入速度分別為1 m/s和0.01 m/s的兩個無樁模型。有限元模型圖如圖2所示。

在使用CEL方法模擬時設定的樁靴貫入速度不同也有可能影響附近樁基礎的響應結果。為分析這種影響,建立兩個樁靴貫入速度分別為1 m/s和0.25 m/s的有樁模型。海洋工程中使用的樁基礎一般為大直徑鋼管樁基礎,由于圓形樁基礎的壁厚較模型的整體幾何尺度小幾個數量級,且必須考慮樁體本身的內力和變形,不能將之簡單地簡化為剛體。根據式(5)和式(6)將大大降低求解所需的臨界步長,增大計算時間。為提高計算速度,有學者提出可以對樁體進行質量放大。本文為驗證質量放大方法的有效性,在保證其他條件都不變的前提下,建立質量放大系數分別為1、10、100和1 000的計算模型。為提高計算效率,將問題進行簡化——圓形薄壁鋼管樁基礎簡化為邊長為0.5 m的實心方樁基礎,樁與樁靴的平面位置關系圖如圖3所示,樁的具體參數如表2所示。

圖2 無樁時有限元模型圖Fig.2 FE model of the spudcan with no piles

圖3 方樁基礎與樁靴的平面位置示意圖Fig.3 Plan of the spuncan and pile foundation

表2 方樁參數表

Table2 Parameters of the square pile

截面長度/m長度/m入泥深度/m彈性模量/GPa泊松比0.5115962100.3

3 模擬結果

3.1 無樁時不同樁靴貫入速度結果對比

圖4為樁靴貫入速度分別為1 m/s和0.01 m/s時土體阻力與貫入深度的關系曲線。由圖可知,雖然樁靴貫入速度相差了100倍,但兩種速度下的土體阻力與貫入深度曲線基本一致。在貫入深度小于30 m時,速度慢的土體阻力較大,當貫入深度大于30 m時,速度快的土體阻力較大。在相同深度處兩者的土體阻力最大相差20%,說明樁靴貫入速度不會影響土體阻力的大小。這是因為土體阻力主要由樁靴下部土體的力學屬性決定,不會因樁靴貫入速度的不同發生變化。

圖5為樁靴貫入速度分別為1 m/s和0.01 m/s時土體的等效塑性云圖。等效塑性云圖表征樁靴貫入過程中土體的擾動和破壞范圍。由圖可知貫入速度為0.01 m/s的土體擾動破壞范圍較貫入速度1 m/s略大,且插樁速度較慢時土體回流現象明顯,速度較大時土體基本沒有回流。這說明樁靴貫入速度對土體的擾動和破壞范圍影響不大,但貫入速度越大對土體的水平向擠壓效應越明顯,土體向四周排開而不會回流。

圖4 不同貫入速度下土體阻力與貫入深度的關系曲線Fig.4 The soil resistance vs.penetration depth at different penetration speeds

圖5 不同貫入速度下土體等效塑性云圖Fig.5 The equivalent plastic contours of soil at different penetration speeds

圖6為樁靴貫入速度分別為1 m/s和0.01 m/s時土體的剪應力云圖。由圖可知,貫入速度不同時土體的剪應力分布形態基本一致,最大剪應力分布在樁靴下部,沿圓弧形向四周遞減。貫入速度較快時土體的最大剪應力也較大,說明樁靴貫入速度越快,土體的剪應力水平越高。

圖6 不同貫入速度下土體Tresca剪應力云圖 (單位:Pa)Fig.6 The Tresca shear stress contours of soil at different penetration speeds (unit:Pa)

雖然在不同的樁靴貫入速度情況下土體的破壞形態和剪應力水平有所區別,但實際工程中更加關注土體阻力沿深度的變化曲線,為提高計算效率,可以人為設定較高的樁靴貫入速度。

3.2 有樁時不同樁靴貫入速度結果對比

圖7為樁靴貫入速度分別為1 m/s和0.25 m/s時樁身的水平位移云圖。由圖可知,當貫入速度為1 m/s時樁身的最大水平位移為1.11 m,當貫入速度為0.25 m/s時為0.99 m,兩者相差不到12%。可近似認為樁靴貫入速度對計算樁的水平變形影響較小。

圖7 不同貫入速度下樁身的水平位移云圖 (單位:m)Fig.7 The horizontal displacement contours of the pile with different penetration speeds (unit:m)

圖8 不同貫入速度下樁身的Mises應力云圖 (單位:Pa)Fig.8 The Mises stress contours of the pile with different penetration speeds (unit:Pa)

圖8為樁靴貫入速度分別為1 m/s和0.25 m/s時樁身的Mises應力云圖。由圖可知,樁身的最大Mises應力當貫入速度為1 m/s時為128 MPa,當貫入速度為0.25 m/s時為163 MPa,表明樁靴貫入速度越大,計算得到的樁身最大Mises應力越少。這是由于樁靴的貫入速度越大,在貫入相同深度處所需的時間越小,土的動能不能有效地轉化為樁的應變能,使得樁身應力較小。

3.3 樁質量放大的影響

圖9為其他條件相同,質量放大系數分別為1倍、10倍、100倍和1 000倍時樁身的水平位移云圖。由圖可知,當質量放大系數分別為1倍、10倍、100倍和1 000倍時樁身的最大水平位移分別為9.9 m、9.7 m、1.19 m和1.18 m。說明隨著樁的質量放大系數增加,樁身最大水平位移略有增大,最大增幅在12%左右。可近似認為樁身最大水平位移對質量放大系數不敏感。

圖9 不同質量放大系數下樁身的水平位移云圖 (單位:m)Fig.9 The horizontal displacement contours of pile shaft with different mass scaling factors (unit:m)

圖10為其他條件都相同質量放大系數分別為1倍、10倍、100倍和1 000倍時樁身的Mises應力云圖。

圖10 不同質量放大系數下樁身的Mises應力云圖 (單位:Pa)Fig.10 The Mises stress contours of the pile shaft different mass scaling factors (unit:Pa)

由圖可知當質量放大系數分別為1倍、10倍、100倍和1 000倍時樁身的最大Mises應力分別為163 MPa、157 MPa、49 MPa和42 MPa。說明質量放大系數對樁身最大Mises應力有較大影響,且兩者關系為非線性的。筆者還對質量放大系數為30倍、50倍和80倍的模型進行模擬(限于篇幅不一一列舉),結果表明隨著質量放大系數的增加,樁身最大Mises應力有下降的趨勢,但無規律可循,因此很難找到不影響計算結果的合適的放大系數。

綜上所述,質量放大方法雖然能顯著提高計算速度,對樁身變形計算影響也較小,但在樁身應力的計算上缺乏準確性,并且很難評價質量放大系數的大小對計算結果的影響。因此在無可靠對照資料的情況下建議慎用質量放大方法。

4 結論與建議

CEL動態模擬方法在樁靴連續貫入問題的研究中被廣泛應用,但現有計算中設定的樁靴貫入速度與實際工程有較大差異。是否可以使用質量放大方法提高計算速度還有待進一步的研究。本文首先對比研究在無樁時不同樁靴貫入速度對計算結果的影響,然后對比分析不同樁靴貫入速度對附近樁基礎響應的影響,最后研究使用質量放大方法對計算結果的影響。結果表明:

(1) 無樁時不同樁靴貫入速度對土體阻力影響較小,土體的破壞形態和剪應力水平略有不同,因此為提高計算效率可設定較高的樁靴貫入速度。

(2) 樁靴貫入速度對樁身水平位移影響較小,貫入速度為1 m/s時樁身的最大水平位移較貫入速度為0.25 m/s時大12%;樁靴貫入速度越大,在貫入相同深度處所需的時間越短,土的動能轉化為樁的應變能也越小,從而使得樁身最大應力越小。

(3) 樁身最大位移對質量放大系數的增加不敏感。隨著質量放大系數的增加,樁身最大Mises應力有下降的趨勢,但并無規律可循,并且很難評價質量放大系數的大小對計算結果的影響,因此在無可靠對照資料的情況下建議慎用質量放大方法。

References)

[1] Hossain M S,Randolph M F,White D T.Punch-through of Spudcan Foundations in Two-layer Clay[C]//Frontiers in Offshore Geotechnics,Taylor and Francis Balkema,The Netherlands,2005,535-540.

[2] Hossain M S,Hu Y,Randolph M F.Limiting Cavity Depth for Spudcan Foundations Penetrating Clay[J].Geotechnique,2005,55(9):679-690.

[3] Hossain M S,Randolph M F.Experimental Investigation of Punch-through Potential for Spudcan Foundations [C]// Proceedings of the 7th International Conference on Physical Modeling in Geotechnics,Zurich,Switzerland.2010:1025-1031.

[4] Hossain M S,Zheng J,Menzies D,et al.Spudcan Penetration Analysis for Case Histories in Clay[J].Journal of Geotechnical and Geoenvironmental Engineering,2014,140(7):04014034.

[5] Xie Y,Leung C F,Chow Y K.Study of Soil Responses Adjacent to Jack-up Spudcan Foundation[C]// Proceedings of the 29thInternational Conference on Ocean,Offshore and Arctic Engineering,Shanghai,China,2010:1-7.

[6] 張浦陽,于曉洋,丁紅巖.海上自升式鉆井平臺插樁階段樁靴承載力計算[J].石油勘探與開發,2011,38(5):613-619.ZHANG Pu-yang,YU Xiao-yang,DING Hong-yan.Spudcan Bearing Capacity Calculation of the Offshore Jack-up Drilling Platform During the Preloading Process[J].Petroleum Exploration and Development,2011,38(5):613-619.(in Chinese)

[7] Qiu G,Henke S,Grabl J.3D Analysis of the Installation Process of Spudcan Foundations[C]//Proceedings of the Frontiers in Offshore Geothchnics Ⅱ.2010:685-690.

[8] Zheng J B,Hossain M S,Wang D.3D Large Deformation FE Analysis of Spudcan and Cone Penetration on Three-layer Clay[C]//Proceedings of the Twenty-third International Offshore and Polar Engineering,Anchorage,Alaska,USA.2013.

[9] Zheng J B,Hossain M S,Wang D.3D Large Deformation FE Analysis of Spudcan Foundations on Layered Clays Using CEL Approach[J].Springer Series in Geomechanics and Geoengineering,2013:803-810.

[10] Arslan H,Wang P C.Advanced Continuum Modeling of Pile Response to Jack-up Spudcan Penetration[C]//Proceedings of the Twenty-fourth International Ocean and Polar Engineering Conference.2014:506-513.

[11] Siciliano R,Hanmilton J,Murff J,et al .Effect of Jackup Spud Cans on Piles[C]//Offshore Technology Conference.1990:OTC 6467.

[12] 吳永韌,魯曉兵,王淑云,等.Spudcan基礎貫入對固定平臺基礎影響[J].中國海洋平臺, 2008, 23(1):35-38.WU Yong-ren,LU Xiao-bing,WANG Shu-yun,et al.Effects Of Spuncan On Pile Foundations Of Fixed Platform[J].China Offshore Platform, 2008, 23(1):35-38(in Chinese)

[13] 吳永韌,王淑云,魯曉兵.樁靴壓入對固定平臺基礎擾動的數值模擬[C]//第三屆全國巖土與工程學術大會論文集.2009:46-49.WU Yong-ren,WANG Shu-yun,LU Xiao-bing.The Effects to the Fixed Platform Foundation by Spuncan Penetration[C]//Proceedings of the 3rdNational Conference on Soil Mechanics and Geotechnical Engineering of China,2009:46-49.(in Chinese)

[14] Wu Y R,Lu X,Zhang X.Effects of the Spudcan Penetration on the Adjacent Foundations[J].Open Ocean Engineering Journal,2010:38-44.

[15] Tho K K, Leung C F,Chow Y K,Swaddiwudhipong S.Eulerian Finite Element Simulation of Spudcan-pile Interaction[J].Canadian Geotechnical Journal,2013,50(6):595-608.

A Dynamic Simulation Method for Continuous Spudcan Penetration

ZHOU Long1, LIU Run1,2, GUO Shao-zeng1, PENG Bi-yao1

(1.StateKeyLaboratoryofHydraulicEngineeringSimulationandSafety,TianjinUniversity,Tianjin300072,China;2.KeyLaboratoryofPortGeotechnicalEngineering,MinistryofTransport,Tianjin300072,China)

The CEL method,which is widely used in the research on spudcan penetration,is a new simulation method.However,there are many questions which have not been answered,especially about the effect of penetration velocity and the mass scaling method.This paper firstly simulates two models with no piles which have different penetration speeds.Then,two models with square piles are also simulated with different penetration speeds.Additionally,the mass scaling method is verified.The simulation results show that the penetration speed has little effect on the soil resistance and more effect on the soil failure mode and shear stress.When the model has a pile,the horizontal displacement is relatively unaffected by the penetration speed,and the Mises stress of the pile increases with the penetration speed.The mass scaling factor also has little effect on the horizontal displacement of the pile.The Mises stress has an irregular relation to the mass scaling factor.Therefore,caution should be carefully taken when the mass scaling method is used.

spudcan; pile foundation; penetration speed; mass scaling

2014-08-20

國家重點基礎研究發展計劃(973計劃)(2014CB046800);國家自然科學基金優秀青年基金(51322904);天津市自然科學基金面上項目(12JCYBJC4700)

周 龍(1988-),男,安徽淮北人,碩士研究生,主要從事海洋樁基礎方面的研究.E-mail:zhoulonglmn@126.com

劉 潤(1974-),女,天津人,教授、博士生導師,主要從事海洋土力學和土與海洋結構物相互作用方面的研究.E-mail:liurun@tju.edu.cn

TU43

A

1000-0844(2015)02-0460-007

10.3969/j.issn.1000-0844.2015.02.0460

猜你喜歡
變形影響方法
是什么影響了滑動摩擦力的大小
哪些顧慮影響擔當?
當代陜西(2021年2期)2021-03-29 07:41:24
談詩的變形
中華詩詞(2020年1期)2020-09-21 09:24:52
“我”的變形計
例談拼圖與整式變形
會變形的餅
擴鏈劑聯用對PETG擴鏈反應與流變性能的影響
中國塑料(2016年3期)2016-06-15 20:30:00
用對方法才能瘦
Coco薇(2016年2期)2016-03-22 02:42:52
四大方法 教你不再“坐以待病”!
Coco薇(2015年1期)2015-08-13 02:47:34
捕魚
主站蜘蛛池模板: 曰韩人妻一区二区三区| 久久久久久国产精品mv| av大片在线无码免费| 亚洲AV成人一区二区三区AV| AV不卡国产在线观看| 激情网址在线观看| 欧日韩在线不卡视频| 国产清纯在线一区二区WWW| 中文字幕久久波多野结衣| 在线日本国产成人免费的| 成人国产精品一级毛片天堂| 欧美日韩在线成人| 亚洲h视频在线| Aⅴ无码专区在线观看| 华人在线亚洲欧美精品| 一区二区理伦视频| 国产日本一线在线观看免费| 日本色综合网| 狠狠色香婷婷久久亚洲精品| 91无码人妻精品一区| 国产美女精品人人做人人爽| 国产麻豆91网在线看| 欧美日韩国产高清一区二区三区| 四虎国产在线观看| 熟妇丰满人妻av无码区| 91青青视频| 精品国产免费观看一区| 亚洲男人在线天堂| 亚洲日本中文字幕乱码中文| www.99精品视频在线播放| 欧美曰批视频免费播放免费| 国产成人久久777777| 国产91蝌蚪窝| 国产成人h在线观看网站站| 久久国产精品电影| 国产真实二区一区在线亚洲| 成人小视频网| 五月天婷婷网亚洲综合在线| 香蕉国产精品视频| 中文字幕在线不卡视频| 白浆视频在线观看| 亚洲人成日本在线观看| 欧美激情二区三区| 浮力影院国产第一页| 一区二区影院| 成人午夜视频在线| 亚洲中文字幕在线观看| 1769国产精品免费视频| 999国内精品久久免费视频| 在线免费a视频| 亚洲嫩模喷白浆| 漂亮人妻被中出中文字幕久久 | 在线无码九区| 在线观看亚洲成人| 中文字幕天无码久久精品视频免费 | 丰满人妻一区二区三区视频| 波多野结衣国产精品| 亚洲美女一区二区三区| 国产天天射| 成人午夜久久| 天堂岛国av无码免费无禁网站| 中文字幕永久在线观看| 老司机aⅴ在线精品导航| 久久黄色影院| 99久久精品国产自免费| 永久成人无码激情视频免费| 综合色88| 亚洲精品无码不卡在线播放| 野花国产精品入口| 精品国产一区二区三区在线观看 | av一区二区无码在线| 国产欧美精品专区一区二区| 黄色不卡视频| 99在线观看免费视频| 一级毛片无毒不卡直接观看| 亚洲美女高潮久久久久久久| 欧美h在线观看| 一级毛片无毒不卡直接观看 | 亚洲综合18p| 亚洲日韩AV无码精品| 国产在线拍偷自揄拍精品| 亚洲综合18p|