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

基于ANSYS和JC法的滑坡抗滑樁結構可靠度分析

2016-02-25 08:29:26羅麗娟夏香波
災害學 2016年1期

羅麗娟,熊 帆,陳 悅,夏香波,王 瑞

(1. 長安大學 建筑工程學院,陜西 西安 710061; 2. 陜西師范大學 旅游與環境學院,陜西 西安 710062;3. 長安大學 地下結構與工程研究所,陜西 西安 710061)

?

基于ANSYS和JC法的滑坡抗滑樁結構可靠度分析

羅麗娟1,2,3,熊帆1,陳悅1,夏香波1,王瑞1

(1. 長安大學 建筑工程學院,陜西 西安 710061; 2. 陜西師范大學 旅游與環境學院,陜西 西安 710062;3. 長安大學 地下結構與工程研究所,陜西 西安 710061)

摘要:在深入分析可靠度指標的幾何意義和JC法的基本原理的基礎上,構造了滑坡抗滑樁結構的功能函數。分析了陜西省吳起縣大路溝巨型黃土滑坡的工程地質條件,根據滑坡的工程地質勘察資料和滑坡治理工程的設計資料,建立了基于ANSYS的抗滑樁樁土體系有限元模型,經過自重作用下樁土結構數值的模擬,驗證了該有限元模型的有效性。在此基礎上,選取滑體土重度γ、滑面土內聚力C和內摩擦角φ、滑面以下土的壓縮模量Espan四個隨機變量,建立了基于ANSYS軟件和JC法的抗滑樁結構可靠度計算模型,并研制了相應算法,算例結果表明:基于ANSYS和JC法的抗滑樁結構的失效概率pspan=0.49%,可靠度指標β=2.58,失效概率略大于蒙特卡洛(MCS)抽樣算法結果和CCD響應面法(RSM)結果,但是基于ANSYS和JC法的效率最高,從而驗證了該可靠度模型及算法的有效性。

關鍵詞:可靠度;JC法;抗滑樁;黃土滑坡;ANSYS

在我國,滑坡是威脅人民生命財產安全的主要地質災害之一。據中國地質環境監測院統計數據,2014年全國共發生地質災害10 907起,共造成349人死亡、218人受傷、51人失蹤,直接經濟損失54.1億元,其中滑坡8 128起,占地質災害總數的74.5%。據統計,近20年來,我國每年因為地質災害造成的經濟損失高達200多億元[1];滑坡治理一直是我國地質災害防治的重要任務,抗滑樁作為一種常見的滑坡治理方案得到廣泛應用[2]。

在完善輔助工程措施的基礎上,滑坡防治工程向輕型化、復合型和機械化施工方向發展,基于可靠度的滑坡防治是這種發展方向的重要實踐。國內外在研究滑坡(邊坡)基于矩方法的穩定性可靠度和支護方案優化設計方面做了大量工作,并取得了許多重要成果。在邊坡(滑坡)穩定性的可靠度計算方面取得的主要成果如下所示。楊令強等[2]利用反演分析和考慮殘差的最小二乘擬合的方法確定斷層粘聚力,反分析確定出抗滑樁的抗力統計參數,利用離散化降維解法求得邊坡穩定的可靠度指標。還有許多科技工作者[3-5]用一次二階矩方法評價邊坡穩定性可靠度。吳明坤等[6]基于強度折減,建立位移迭代有限元模型模擬樁-土間的相互作用,基于強度折減的位移迭代有限元法對邊坡抗滑樁加固工程進行可靠性分析。

近年來,許多學者將結構可靠度的矩方法應用于地下結構和巖土工程構筑物的可靠度分析,主要有將一次二階矩與拉普拉斯漸進法相結合在驗算點處作二次展開[7]、一次二階矩驗算點法及其改進[8-9]、二次二階矩法[10]、剩余推力法與JC 法的相結合[11]等等。范雷、唐輝明等[12]引入二元函數插值逼近的一次二階矩方法對邊坡穩定性可靠度進行了分析。有的學者在巖土工程參數敏感性方面進行了研究,取得了許多成果[13-15],還有學者對樁基[16]、擋土墻結構[17]、基坑支護結構[18]的可靠度計算方法進行了研究。

滑坡的剩余下滑力是作用在抗滑樁結構上的荷載,剩余下滑力的大小主要由滑體自重、滑面的抗剪強度及滑面形態決定,荷載的大小及分布決定了結構的尺寸和強度要求。從可靠度的角度來看,滑面土的抗剪強度指標及其統計分布對大型滑坡的治理成本和抗滑結構的可靠性起決定作用。怎樣在滑坡治理工程的造價經濟性與結構可靠性之間尋找最優解是工程界和廣大學者共同關心的問題。

基于以上背景,作者對某黃土滑坡滑帶土的抗剪強度指標、滑體土的重度等開展了大量的室內實驗和統計分析工作,得到了滑體和滑面土的有關隨機變量的統計分布特征,在此基礎上,基于ANSYS和JC法對抗滑樁結構進行了可靠度分析。

1基于JC法的抗滑樁可靠度理論

1.1 JC法的基本原理

定義抗滑樁結構極限狀態功能函數Z為:

Z=g(X1,X2,…,Xn)。

(1)

(2)

則可定義結構可靠指標β為:

(3)

若將基本隨機變量Xi(i=1,2,…,n)進行標準正態化,即:

(4)

則原極限狀態方程Z=g(X1,X2,…,Xn)=0轉化為Z′=f(Y1,Y2,…,Yn)=0,極限狀態超曲面Z=0上任意一點P在n+1維正態空間(Z,X1,X2,…,Xn)中的坐標向量為(z,x1,x2,…,xn)在n+1維標準正態空間(X′,Y1,Y2,…,Yn)中的坐標向量為(z′,y1,y2,…,yn)它們是一一對應的,即:

(5)

(6)

P′*到坐標原點O′的距離為:

(7)

因此,可靠度指標β的幾何意義就是n+1維標準正態空間(Z′,Y1,Y2,…,Yn)的坐標原點O′到極限狀態超曲面Z′=0的最短距離,即

β=min(d)。

(8)

由此可得到結構失效概率Pf=Φ-1(β),見圖1。

圖1 失效邊界與中心點的關系

在抗滑樁可靠度分析中,通常有多個基本隨機變量,當隨機變量為非正態變量時,采用R-F變換方法可以將非正態變量當量正態化[19]。

在一次二階矩驗算點法中,當基本隨機變量Xi、Xj之間相互獨立時,則Z與Z′的均值、標準差分別表示為:

(9)

(10)

(11)

(12)

當基本隨機變量Xi、Xj之間不相互獨立時,要考慮相關系數的影響[19],可靠度指標β仍然用下式計算:

(13)

1.2 抗滑樁結構的功能函數

實際工程中,鋼筋混凝土抗滑樁失效模式多數由于抗滑樁結構受拉側拉應變超過允許值,本文算例中混凝土強度等級為C30,熱軋鋼筋為HRB400,分別根據混凝土和鋼筋的抗拉強度設計值、彈性模量計算得到混凝土極限拉應變εcu=3.3×10-3、鋼筋極限拉應變εy=1.8×10-3。綜合考慮結構抗力的安全性,取抗滑樁結構允許拉應變為[ε]=1.0×10-3,當樁身結構受拉側最大拉應變εy,max=1.0×10-3時,認為樁身結構達到極限狀態,因此,定義滑樁結構功能函數為:

Z=max(ε)-1.0×10-3。

(14)

式中:max(ε)為抗滑樁結構最大正應變εy, max。

2抗滑樁樁土體系有限元計算

2.1 工程背景

現以陜西省吳起縣大路溝滑坡治理工程抗滑樁結構為工程背景。大路溝滑坡屬老滑坡,受坡腳路塹工程開挖(滑坡體中部、上部為多級旱地,中下部被人工開挖為11級臺階)的影響形成體積約856.8×104m3的大型滑坡。地形陡峻、南低北高,溝谷切割深度大,呈深窄的“V”型斷面,屬于黃土梁峁谷坡地貌;滑坡總體縱向坡度約30°。滑坡體海拔介于1 505~1 630 m之間,相對高差125 m;滑坡體長約350 m,寬約400 m,平面約15.2×104m2。大路溝滑坡全貌見圖2。

該滑坡治理方案的主要擋土結構為抗滑樁,樁截面分為1.5×2.0 m、2 m×3 m、2.5 m×3.5 m等幾種形式,樁長分別有20 m、24 m、26 m等三種,縱斷面工程結構布置見圖3。

圖2 大路溝滑坡全貌圖(鏡向65°)

圖3  大路溝滑坡縱斷面工程結構布置圖

大路溝滑坡各地層(巖性)由新到老依次描述如下:

(1)第四系(Q)——馬蘭黃土;離石黃土;淺灰色、黃色沖積粉土、粉質粘土,局部含有卵礫石;

(2)新生界新近系保德組(N2b)——主要為深紅或紫紅、棕紅色粘土巖(或砂質粘土巖),含鈣質結核,但層理不明顯;

(3)中生界白堊系華池組(K1h)——上部為強~中等風化棕紅色塊狀細粒長石砂巖;中部可見棕褐色泥巖、咖啡色頁巖、砂質泥巖(或泥巖)互層;下部為棕紅~暗紅色厚層塊狀細粒長石砂巖,夾有咖啡色粉砂質泥巖。各土層物理指標的統計平均值如表1所示。

2.2 基于ANSYS的有限元模型

現以樁長為20m,截面尺寸為1.5m×2.0m的抗滑樁為例,采用大型通用有限元軟件ANSYS建立滑坡抗滑樁樁土體系有限元模型[20]。以樁身結構、滑體及滑面為主要研究對象,不考慮鋼筋混凝土材料的非線性,假定鋼筋混凝土樁為均質線彈性體,采用鋼筋混凝土均質體等效模型建模;不考慮地下水滲流、地震作用的影響,為方便后續可靠度的計算,暫考慮二維問題;抗滑樁截面為矩形,數值模型中樁截面高度和寬度根據樁間距換算得到,樁長取工程樁實際長度;滑帶土體、滑床、滑體采用D-P彈塑性模型;考慮樁側、樁端與土的接觸,樁與土變形協調,交界面無滑移;滑體與滑床之間設置為接觸單元。

因大路溝滑坡表面形態較復雜,為便于有限元模型的建立,對滑坡縱斷面表面形態(見圖3)作適當簡化。根據勘察、設計資料,該滑坡的幾何模型及結構尺寸等參數見表2;該有限元模型中各土層物理指標見表1;樁體、滑體及滑床材料計算參數見表3;樁側與樁側土體、樁端與樁端土體、滑體與滑床界面處接觸單元和參數設置見表4。據此所構建的有限元模型共劃分2 713個節點、63個接觸單元、2 661個單元,有限元模型單元網格見圖4。

表1 土的物理指標統計平均值

表2 滑坡幾何模型及結構尺寸參數

表3 樁體、滑體及滑床材料計算參數

注:表中“-”表示沒有該項數據。

表4 接觸單元和參數設置表

表5 隨機變量統計特征值

圖4 有限元模型單元網格圖

2.3 計算分析

為了計算土體自重作用下滑坡體的變形,設置如下邊界條件:左右兩側限制水平方向位移,取Ux=0;底部邊界限制豎向位移,取Vy=0;頂部邊界設為自由邊界。有限元模型計算得到的Y方向位移云圖見圖5。

圖5 Y方向位移Vy云圖

由圖5可知:自重作用下,滑體和滑床沿滑面的相對滑移非常突出,且剪出口和滑坡后壁相對位置明顯,顯示整體滑動明顯。

樁體Y方向正應力σy云圖見圖6。

由圖6可知:樁背側受拉明顯,最大拉應力達18.6MPa,樁前側受壓明顯,最大壓應力為19.8MPa,抗滑樁符合受彎變形特征,但受彎受壓區范圍有限,顯示樁長安全儲備偏大。圖5和圖6顯示了該抗滑樁樁土有限元模型的有效性。

圖6 樁體Y方向正應力σy云圖

3基于JC法的抗滑樁可靠度計算

3.1 抗滑樁樁土結構體系的隨機變量特征

在抗滑樁樁土結構體系中,樁身結構所受荷載為滑坡體剩余下滑力,樁身結構抗力的發揮由樁身嵌固段圍巖反力提供。因滑坡巖土體的物理力學性質指標的獲得具有不確定性,導致抗滑樁結構設計成果具有一定的失效風險。滑坡的下滑力主要是滑體自重γ提供的,抗滑力主要由滑帶土內摩擦角φ和內聚力C控制(不考慮滑坡形態帶來的影響),抗滑樁抗力的發揮與滑面以下土層的壓縮模量Es有關,因此,本文主要考慮此4個隨機變量對抗滑樁結構可靠度的影響。采用SPSS軟件對此四個隨機變量作統計分析,經擬合、K-S假設檢驗后得到滑面及滑體各隨機變量的統計分布及特征值見表5。

3.2 基于ANSYS和JC法的抗滑樁可靠度分析

基于ANSYS和JC法求解抗滑樁可靠度主要思路為:可靠度計算采用基于JC法的FORTRAN源程序(Relia-JC.FOR),抗滑樁有限元模型是基于ANSYS建立的,FORTRAN源程序(Relia-JC.FOR)所需要調用的功能函數值是通過調用ANSYS命令流文件PR2D.txt實現返回樁身最大正應變εy,max。

其相應的過程:源程序(Relia-JC.FOR)將隨機變量取值寫入PRINPUT.TXT,FORTRAN源程序(Relia-JC.FOR)執行至調用ANSYS命令流文件,通過執行*VREAD,HMCF(1),′PRINPUT′,′TXT′,實現對幾何模型和有限元模型參數(隨機變量)賦值,ANSYS計算結束后,通過*VWRITE,MAXSTRY將樁身最大正應變(可靠度計算中功能函數組成部分)MAXSTRY(εy,max)寫入至PROUT.TXT文件中,供FORTRAN源程序(Relia-JC.FOR)調用傳遞抗滑樁結構功能函數值。

基于ANSYS和JC法計算可靠度流程見圖7。

圖7 基于JC法和ANSYS的抗滑樁可靠度計算流程圖

FORTRAN語言源程序(Relia-JC.FOR)對ANSYS命令流文件PR2D.txt主要的調用語句如下:

LOGICAL(4)result

result=SYSTEMQQ('C:ProgramFilesAnsysIncv90ANSYSinintelANSYS90-b-pane3fl-iD:programPILERELIAPR2D.txt-oD:programPILERELIAPROUT.TXT')

基于ANSYS軟件計算抗滑樁結構功能函數的命令流文件PR2D.txt的主要命令流如下:

!(1)定義相關變量參數

/CLEAR

*DIM,HMCF,ARRAY,4!定義四個隨機變量組成的數組

!定義HUATI,MOLIANG,COHENSION,FRICTION變量

*CREATE,TEMP1 !創建臨時宏文件TEMP1

*VREAD,HMCF(1),'PRINPUT','TXT' !隨機變量讀入

(F7.2,/,E7.2,/,F8.2,/,F5.2)!格式同PRINPUT.TXT格式

*END!結束宏文件

/INPUT,TEMP1 !四個隨機變量賦值

/PREP7

!(2)定義單元類型、實常數、材料屬性

ET,2,PLANE182 !周圍土體實體單元

MP,EX,1,3.4122E10 !抗滑樁結構彈摸、泊松比、密度MP,PRXY,1,0.18

MP,DENS,1,2500

MP,EX,2,1.5E7!滑床及樁前土壓縮模量、泊松比、密度

MP,PRXY,2,0.30

MP,DENS,2,1800

TB,DP,2!DP模型

TBDATA,1,41.0,13.7,10.0

MP,EX,3,1.5E7 !滑體土體壓縮模量、泊松比、密度

MP,PRXY,3,0.32

MP,DENS,3,HUATI!滑體密度(隨機變量)

TB,DP,3!DP模型

TBDATA,2,28.3,26.3,15.0

MP,EX,4,MOLIANG!滑帶土壓縮模量(隨機變量)

MP,PRXY,4,0.35

MP,DENS,4,HUATI!滑帶土密度(隨機變量)

TB,DP,4!滑帶土體內聚力和內摩擦角(隨機變量)

TBDATA,3,COHENSION,FRICTION,5.0 !

!(3)創建邊坡及抗滑樁模型關鍵點(略)

!(10)后處理

/POST1

ESEL,S,MAT,,1!提取樁身結構

PLNSOL,S,Y,0,1!應力云圖

SET,FIRST

NSORT,EPTO,Y!將樁身節點正應變 排序

*GET,MAXSTRY,SORT,,MAX!提取最大正應變

*CREATE,TEMP2 !定義臨時宏文件TEMP2

*CFOPEN,'PROUT','TXT',' '!功能函數值存儲文件

*VWRITE,MAXSTRY!寫入 至PROUT.TXT文件

(E11.3) !數據格式E11.3

*CFCLOS!關閉PROUT.TXT文件

*END!宏結束標志

/INPUT,TEMP2 !

FINISH

3.3 可靠度計算結果分析

基于ANSYS和JC法計算結構的失效概率Pf=0.49%,可靠度指標β=2.58,與蒙特卡洛法和響應面法計算結果對比如表6所示。

表6 幾種可靠度計算方法結果對比

由表6可知,相同條件下計算抗滑樁結構失效概率從小到大的方法依次為:蒙特卡洛法(MCS)、響應面(RSM)中心復合設計法、基于ANSYS和JC法,但三種方法算得的失效概率相差不大。但是,計算效率從低到高依次為:蒙特卡洛法(MCS)、響應面(RSM)中心復合設計法(試驗點數量N=25)、基于ANSYS和JC法(無需抽樣,僅通過數值方法實現功能函數對隨機變量求偏導數)。

4結論

本文以大路溝滑坡抗滑樁工程為算例,基于ANSYS建立了抗滑樁樁土體系有限元模型和抗滑樁可靠度模型, 基于ANSYS和JC法進行抗滑樁可靠度計算,可得出以下幾點認識。

(1)ANSYS建模時,采用鋼筋混凝土均質體等效模型建立抗滑樁結構有限元模型是可行的,建立的基于ANSYS和JC法的抗滑樁可靠度模型是有效的。

(2)基于各方法的抗滑樁可靠度計算結果表明:蒙特卡洛法(MCS)計算得到的失效概率最小,響應面(RSM)中心復合設計法次之,基于ANSYS和JC法計算得到的失效概率最大,但三種方法算得的失效概率相差不大。蒙特卡洛法(MCS)計算效率最低,響應面(RSM)中心復合設計法次之,基于ANSYS和JC法效率最高。

參考文獻:

[1]羅麗娟,趙法鎖. 滑坡防治工程措施研究現狀與應用綜述[J].自然災害學報,2009,18(4): 158-164.

[2]楊令強,馬靜,張社榮. 抗滑樁加固邊坡的穩定可靠度分析[J].巖土工程學報,2009,31(8):1299-1302.

[3]王建華,易朋瑩,鄧繼輝,等. 滑坡穩定性可靠度分析[J]地下空間與工程學報,2007, 3(4):668-672.

[4]唐亞明. 基于可靠度的黃土斜坡穩定性分析[J].地質通報,2008,27(8):1217-1222.

[5]張智洪,鄭紅娟. 土坡穩定性可靠度分析[J].重慶交通大學學報:自然科學版,2007,26(Supp.1):65-69.

[6]吳明坤,王建國,譚曉慧. 基于強度折減的邊坡抗滑樁可靠性分析[J].工業建筑,2011,41(5):99-103.

[7]姚澤良,李寶平,周雪峰. 結構可靠度分析的一次二階矩方法與二次二階矩方法[J].西北水力發電,2005,21(3):20-23.

[8]李景奎,張義民. 正態分布連續體結構可靠性拓撲優化設計[J].東北大學學報:自然科學版,2011,32(9):1304- 1307.

[9]吳雪莉,姜進,張樂民,等. 改進的一次二階矩法在隧道可靠度中的應用[J].工業建筑,2008,38(Supp.1):694-697.

[10]宋云連,焦同戰,趙海芳. 部分二次二階矩法在公路邊坡可靠度分析中的應用[J].公路交通科技,2007,24(9):1-5.

[11]羅麗娟,趙法鎖,胡江洋,等. 基于剩余推力法的黃土高邊坡穩定性可靠度分析[J].長安大學學報:自然科學版,2008,28(4):27-31.

[12]范雷,唐輝明,胡斌,等. 二元函數逼近在斜坡可靠度分析中的應用研究[J].巖土力學,2008,29(3):624-628.

[13]周艷. 重力式擋墻可靠度及其敏感度分析[J].國外建材科技,2008,29(3):116-118.

[14]馬德云,左勇志,霍達,等. 結構可靠性對隨機參數的敏感性研究[J].工程抗震與加固改造,2008,30(2):33-38.

[15]胡曉軍,吳延枝. 抗滑樁錨固深度的可靠度與參數敏感性[J].河海大學學報:自然科學版,2011,39(1):49-53.

[16]丁繼輝,袁滿,王建國. 多樁型組合樁復合地基承載力的可靠度分析[J].工程力學,2008,25(Supp.2):168-172.

[17]杜永峰,余鈺,李慧. 重力式擋土墻穩定性的結構體系可靠度分析[J].巖土工程學報,2008,30(3):349-353.

[18]尹盛斌,丁紅巖. 基坑支護結構的漸近積分可靠度分析方法及應用[J].巖土力學,2011,32(12):3728-3732.

[19]胡志平,羅麗娟. 管片襯砌結構可靠度分析的優化方法[J],巖石力學與工程學報,2005,24(22):4145-4150.

[20]王秀麗,董文燕,金兆鑫. 品字型抗滑樁對滑坡的加固效果研究[J]. 災害學,2015,30(2):8-10.

劉曉梅,劉國忠,趙金彪. 廣西暴雨過程雨團時空分布特征研究[J].災害學, 2016,31(1):39-44.[ Liu Xiaomei, Liu Guozhong and Zhao Jinbiao. Temporal and Spatial Distributions of Rain Cluster about the Rainstorm Process in Guangxi[J].Journal of Catastrophology, 2016,31(1):39-44.]

羅麗娟,熊帆,陳悅,等. 基于ANSYS和JC法的滑坡抗滑樁結構可靠度分析[J].災害學, 2016,31(1):33-38.[ Luo Lijuan, Xiong Fan,Chen Yue,et al.Reliability Analysis of Anti-slide Pile Based on JC Method and ANSYS[J].Journal of Catastrophology, 2016,31(1):33-38.]

Reliability Analysis of Anti-slide Pile Based on JC Method and ANSYS

Luo Lijuan1, 2, 3, Xiong Fan1, Chen Yue1, Xia Xiangbo1, and Wang Rui1

(1.SchoolofCivilEngineering,ChanganUniversity,Xi’an710061,China;

2.TourismandEnvironmentCollegeofShaanxiNormalUniversity,Xi’an710062,China;

3.InstituteofUndergroundStructureandEngineering,ChanganUniversity,Xi’an710061,China)

Abstract:Based on the in-depth analysis of the geometric meaning of reliability index and the basic theory of JC method, the performance function of anti-slide pile soil system is proposed. On the basis of analysis of engineering geological conditions of Dalugou loess landslide in Wuqi County Shaanxi Province, according to the landslide geological survey data and design information of landslide remediation engineering, the finite element model-based ANSYS software platform of anti-slide pile soil system is established. The numeric simulation results of pile and soil system under own weight condition proves the validity of the finite element model. Moreover, the density γ of landslide, cohesion c and internal friction angle φ of slide surface, and compression modulus Es, four random variables of sliding bed, the reliability calculation model of anti-slide pile based on ANSYS software and JC method is established, and the corresponding algorithm is developed. The sample calculation shows that, the probability of failure of anti-slide pile structure Pf=0.49% and β=2.58, which is little larger than that of MCS Sampling method, and that of CCD Responding Surface method. The example verified the validity of the reliability model and corresponding algorithm.

Key words:reliability; JC method; anti-slide pile; loess landslide; ANSYS

作者簡介:羅麗娟(1973-),女,山西平遙人,博士(后),副教授,主要從事邊坡(滑坡)工程與樁基工程等方面的教學和科研研究. E-mail:luojuan@chd.edu.cn

基金項目:國家自然科學 (41202190);長安大學中央高校基本科研業務費專項資金資助(2013G3282014;2014G2260009)

收稿日期:2015-07-01修回日期:2015-08-21

中圖分類號:P642.22;X43

文獻標志碼:A

文章編號:1000-811X(2016)01-0033-06

doi:10.3969/j.issn.1000-811X.2016.01.008

主站蜘蛛池模板: 亚洲无线国产观看| 制服丝袜亚洲| 国产喷水视频| 亚洲男人的天堂在线| 亚洲成aⅴ人片在线影院八| 一级毛片中文字幕| 国产成人高清精品免费| 综合社区亚洲熟妇p| 97se亚洲综合不卡| 国产高清无码第一十页在线观看| 国产精品视屏| 国产后式a一视频| 久久久久人妻一区精品色奶水| 本亚洲精品网站| 特级精品毛片免费观看| 白浆免费视频国产精品视频| 國產尤物AV尤物在線觀看| 精品無碼一區在線觀看 | 中文字幕亚洲综久久2021| 国产精品自在在线午夜区app| 日韩专区欧美| 99热这里只有精品免费| 国产视频a| 精品久久久久久成人AV| 国产丝袜一区二区三区视频免下载| 午夜性爽视频男人的天堂| 成年人国产视频| 欧美日本视频在线观看| 欧美一级专区免费大片| 国产欧美网站| 99久久国产综合精品2023| 欧美午夜在线视频| 国产精品午夜电影| 日韩高清无码免费| 亚洲黄色高清| 亚洲日韩图片专区第1页| 午夜啪啪福利| yjizz国产在线视频网| 亚洲天堂色色人体| 欧美特黄一免在线观看| 在线欧美一区| 中文字幕在线播放不卡| 狠狠综合久久| 日韩无码视频播放| 91在线激情在线观看| 在线看片免费人成视久网下载| 这里只有精品在线| 日韩东京热无码人妻| 青青草原偷拍视频| 亚洲国产精品成人久久综合影院| 91区国产福利在线观看午夜| 美女免费黄网站| 无码中字出轨中文人妻中文中| 国产成人精品视频一区二区电影| 国产亚洲美日韩AV中文字幕无码成人 | 国产又爽又黄无遮挡免费观看| 黄色福利在线| 无遮挡国产高潮视频免费观看| 精品国产自在在线在线观看| 成人亚洲天堂| 国产精品爽爽va在线无码观看 | 热热久久狠狠偷偷色男同| 国产精品乱偷免费视频| 无码精油按摩潮喷在线播放| 成人亚洲国产| 国产精品yjizz视频网一二区| 国产成人艳妇AA视频在线| 精品久久香蕉国产线看观看gif| AV不卡在线永久免费观看| 精品自窥自偷在线看| 成人免费午夜视频| 中文一级毛片| 91色国产在线| 国产又大又粗又猛又爽的视频| 国产精品成人免费视频99| 久久精品国产免费观看频道| 91免费观看视频| 午夜福利免费视频| 国产黑人在线| 国产美女精品人人做人人爽| Jizz国产色系免费| 中文字幕 欧美日韩|