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

基于模型試驗及耦合數值分析的XCC樁擠土效應研究

2022-11-21 02:57:10李康樂邱恩喜孫希望聶文峰
土木工程與管理學報 2022年5期
關鍵詞:模型

李康樂, 邱恩喜,2, 孫希望, 張 蕊, 聶文峰

(1.西南石油大學 土木工程與測繪學院,四川 成都 610500;2.成都理工大學 地質災害防治與地質環境保護國家重點實驗室,四川 成都 610059;3.中鐵二院昆明勘察設計研究院有限責任公司,云南 昆明 650200)

在我國工程建設中,擠土樁因其可以改善地基土性能,起到補強作用,提高地基承載力,故被廣泛地運用于地基加固中[1],而靜壓施工法,以其施工速度快、噪音小、擾動小等特點已經逐漸成為一種重要的施工手段。但靜壓沉樁時,由于樁體入土將排開相應體積的土體,從而對樁身土體產生側向擠壓,樁端部土體產生側向下的壓縮剪切等行為,使樁身周圍土體的應力狀態發生變化,將不可避免地對相鄰建筑物及環境產生影響[2]。因此,研究靜壓樁沉樁機理特性具有重要的現實意義。

國內外學者針對沉樁擠土效應開展了廣泛的研究,目前擠土效應的研究理論主要有圓孔擴張理論及應變路徑法,研究手段主要有數值分析、模型試驗及現場試驗。羅占友[3]就沉樁擠土效應室內模型試驗研究進展進行了論述,指出現有的模型實驗中存在土體內部位移變化獲得較難,土樣制備難以保持原有特性及群樁相互作用機理的研究缺乏等問題。針對上述存在的缺點,不少學者采用調配的透明土結合測量系統開展室內模型實驗研究。Lehane和Gill[4]利用透明土樣中放置示蹤點的技術來研究沉樁擠土效應;White[5~9]采用半模試驗結合自主研發的PIV技術捕捉對稱面位置處的整個土體位移場來研究沉樁擠土效應位移場分布規律;Ni和Hird等[10,11]則結合透明土和PIV的技術來研究圓形截面樁、螺旋樁以及管樁的沉樁擠土效應問題。上述研究雖然實現了對土體內部位移可視化,但其調配的透明土樣與真實土體存在一定差別。為此,本文選擇制備重塑土樣,并采用在玻璃箱上打孔,借助全站儀觀察插于孔中的大頭針頭部位移變化,以實現對沉樁擠土半模型試驗的土體內部位移的可視化。

沉樁過程中樁土相互作用,因其涉及到幾何非線性、接觸非線性、材料非線性、大變形、本構模型等一系列的問題,僅憑借理論公式的推導和試驗研究難以滿足越來越精確的工程技術要求。有限單元法的出現,使有限元分析成為了解決這些非線性復雜問題的一個有效途徑[12]。葉成銀[13]、陳振華[14]等依托空間網格搜索法,建立FLAC-PFC耦合沉樁數值模型,從宏細觀的角度分析了沉樁機理,但計算過程繁瑣。譚鑫等[15]則通過wall(墻)建立FLAC-PFC之間的信息交流通道,實現FLAC-PFC的直接耦合,模擬了碎石樁單樁受荷室內模型試驗的變形及破壞全過程。

近年來,自劉漢龍教授科研團隊研發了現澆X形混凝土樁(XCC樁)后,因其與傳統的圓樁、方樁相比,增大了樁的摩擦面積,大幅提高了承載能力,故極大地推動了XCC樁在工程中的應用與研究。為研究XCC樁荷載傳遞機理和沉降特性,劉芝平[16]、王智強[17]、張敏霞[18,19]等先后利用大型土工模型試驗槽開展了XCC樁和圓樁的足尺模型試驗,對XCC樁承載力、端阻力和側摩阻力進行分析得到了比較類似的結果。曹兆虎等[20]通過透明土材料及PIV測量技術開展了對XCC樁極限承載力的可視化研究,得到了不同的XCC樁樁端破壞模式。

上述研究主要針對整體的承載性能,開展了包括理論分析、足尺模型試驗、現場試驗和數值模擬研究等,對XCC樁的沉樁擠土效應研究較少。周航等[21,22]通過模型試驗提出了可用于XCC異形樁及矩形樁的基于圓孔擴張理論的修正解。Kong等[23]研究了XCC樁在軟土地基貫入過程中樁周土體位移和應力等行為,并提出了一種適用于XCC樁的改進圓孔擴張理論,同時研究得出XCC樁的貫入效果與圓形樁類似,相同截面面積的(等效)圓形樁半徑可以代替XCC樁的半徑。劉漢龍等[24]通過現場試驗分析了X形樁的擠土效應,得出沿尖角方向的擠土壓力大于凹弧方向的擠土壓力。

上述研究表明,異形樁有著區別于普通圓形樁的擠土效應,因此借助有限差分-離散元耦合數值模擬方法,開展對異形樁與普通圓形樁的比較研究,通過類比普通圓形樁的沉樁擠土規律,得到XCC樁沉樁擠土位移場的分布規律,對工程應用具有重要的研究意義。

1 室內沉樁模型試驗

1.1 室內模型試驗方案

(1)模型箱尺寸。模型箱尺寸為60 cm×60 cm×60 cm,兩側各設4個排水孔,在沉樁側設位移觀測孔。

(2)土樣制備。試驗土樣為重塑土樣,土樣基本物理性質見表1,圖1為土樣的級配曲線。

表1 土樣基本物理性質

圖1 試驗土樣顆粒級配

模型箱底部鋪設60 mm厚細砂隔層,并于細砂及排水孔處鋪設透水土工布。將土樣按略大于液限配置成飽和土,分層搗實,土體頂面鋪設透水土工布后加壓固結,直至孔隙水壓力趨于穩定,60 d后可認為固結良好。固結完成后,關閉排水孔,加入適量水,使其達到飽和,靜置15 d后進行沉樁試驗。

(3)模型樁制備。根據相似理論,模型樁幾何比例n=25,樁體采用樹脂abs材料制作。試驗采用半樁模型,見圖2(圖中,rd為XCC樁截面等效半徑;t為圓形樁截面半徑,t=rd;D為圓形樁截面直徑;r0為XCC樁各角度下的截面半徑)。由圖2可知,在0~45°范圍內,隨著θ角度的變化,XCC樁的θ角對應下的截面半徑r0隨之變化,樁體幾何參數詳見表2。

圖2 模型樁實物/mm

表2 樁體參數

(4)試驗測量內容:土體內部位移和土表隆起量。

1)土體內部位移。試驗過程中借助全站儀觀察模型箱側細針的頭部位移,確定土體內部的徑向和豎向位移。

2)土表隆起量。試驗過程中通過土體表面設置的百分表觀測土表隆起量。

1.2 試驗結果分析

(1)土體內部徑向位移分析。由圖2知,r0在0~45°范圍內隨角度而變,為了研究XCC樁的沉樁擠土效應,分別進行了沿45°角(θ=45°)的凹弧段和沿0°角(θ=0°)的十字段進行半樁沉樁模擬,并與與其截面等效半徑rd相等的普通圓形樁對比研究。圖3,4為半樁沉樁完成后位移圖(圖中,z為沉樁深度;r為距樁軸線的距離,即擴孔后土體的徑向位置;θ=0°(45°)表示XCC樁的0°角或45°角;rd=1.2 cm表示與XCC樁截面等效半徑rd相等的普通圓形樁)。

圖3 徑向位移分布

圖3為XCC樁十字段、凹弧段及與XCC樁相同截面面積下的圓形樁沉樁后徑向位移分布圖。由圖3a可見,沉樁后,隨著距樁軸線距離的增大,樁周土體徑向位移逐漸減小且呈對數衰減,試驗結果與文獻[25]研究規律基本一致;XCC異形樁的徑向擠土位移約為圓形樁的1.0~1.4倍,且在r<1.25D范圍內,XCC樁的截面半徑越大,其對應角度下的徑向擠土位移越大;而在r≥1.25D范圍內,XCC樁的徑向擠土位移與其各角度下的截面半徑成反比。

由圖3b可見,隨著深度的增加,距樁軸線同一距離的徑向位移逐漸增大,接近樁端時快速減小。主要是由于土體較軟,樁土界面摩阻力較大,土體被拖拽朝下,致使隨著沉樁深度的增加,樁端下部土體更加密實,進而徑向擠土位移呈逐漸增大后迅速減小的規律。

(2)土體內部豎向位移分析。由圖4可知,沉樁后,豎向位移約為徑向位移的2~5倍。隨著距樁軸線距離的增大,樁周土體豎向位移逐漸減小;隨著深度的增加,距樁軸線同一距離的豎向位移呈先增大后逐漸減小的分布規律,其豎向位移最大值出現在距土體表面約5D處。

圖4 豎向位移分布

(3)土體隆起量。圖5為土體表面豎向位移圖,自沉樁開始至沉樁結束,樁周土體以凹陷為主,在近樁區1.25D范圍內的土體凹陷量為1.5~3.5 mm。這主要是由于室內模型試驗采用飽和軟黏土,在沉樁擠土過程中軟黏土隆起量較小,室內模型試驗結果表現為以凹陷為主,而實際工程中地表土體通常為一層硬塑狀黏土,在擠土樁施工過程中土體隆起量較大,與室內試驗結果存在一定差異。

圖5 土體表面位移量

2 FLAC-PFC耦合數值分析

2.1 模型建立

采用FLAC-PFC耦合數值分析方法,根據室內沉樁模型試驗結果進行計算參數反演。依據反演參數研究XCC樁樁周土體位移場,提出XCC樁擠土位移公式。

(1)三軸參數標定

由于離散的顆粒集合體是基于顆粒間接觸模型,而連續介質是基于土體本構模型,所以首先要進行離散單元的宏細觀參數標定,以使得離散單元與連續介質在宏觀表現上趨于一致。分別建立有限差分和離散集合體的三軸模型試驗進行反演,不斷調整顆粒的細觀參數(顆粒半徑Rm;顆粒孔隙率e;顆粒接觸模量E′和顆粒間摩擦系數ν)以匹配土體的宏觀參數變形模量Es和泊松比μ,使離散元數值試樣的宏觀力學表現逼近實際巖土體的真實行為,從而建立起能夠反應實際巖土體材料強度與變形特征的離散顆粒細觀力學模型[26]。

對離散顆粒和連續介質進行三軸壓縮實驗,來確定離散顆粒的細觀參數。模擬試件的尺寸為直徑6.18 cm,高度12.5 cm。根據試驗數據反演后確定的數值模型計算參數見表3,4。

表3 連續介質(FLAC)土體參數

表4 離散介質(PFC)土體參數

通過PFC和FLAC模擬三軸壓縮下的應力-應變曲線,如圖6所示。標定的細觀參數能夠有效地反映物體實際的宏觀特征。

圖6 三軸參數標定

(2)沉樁模型建立

采用連續-離散耦合方法對沉樁擠土效應進行研究。變形較大區域選擇球體顆粒,變形較小區域選擇連續介質。模型如圖7所示。模型樁長度為L=0.4 m;模型寬度為1.5L=0.6 m;高度為2L=0.8 m;離散區域半徑為0.03 m,高度為2L=0.8 m。

圖7 沉樁模型示意

模型四周約束水平位移,底部采用固定支座,模型頂面為自由面。按0.1 mm/s的位移速率模擬沉樁。

2.2 數值模擬結果分析

(1)各角度下XCC樁與圓形樁徑向位移比較分析

圖8為沉樁結束后徑向位移對比圖,可見數值分析結果與室內模型試驗數據變化規律基本一致,兩者數值差值在10%以內。

圖8a為距土體表面20 cm處樁周土體徑向位移,隨著距樁軸線距離的增大,樁周土體徑向位移逐漸減小且服從冪指函數uXr=rX(θ)(rX(θ)/(kr))n(θ)[21](式中:uXr為XCC樁在距離樁心距離為r時的徑向位移;n(θ)為給出方程中關于θ的變量;rX(θ)為XCC樁截面邊界的半徑(它是關于極角的一個函數);k為通過試驗及數值模擬結果得出的修正系數)衰減規律,與試驗結果一致;在同一截面面積下,XCC樁的徑向擠土位移也約為圓形樁的1.0~1.4倍。

由圖8b可知,在r<1.25D范圍內,XCC樁的截面半徑越大,其對應角度下的徑向擠土位移越大,且徑向擠土位移場呈十字形不均勻擴張;而在1.25D≤r<2.5D范圍內,XCC樁的徑向擠土位移與其各角度下的截面半徑成反比,且徑向位移場呈現出方形不均勻擴張;在r≥2.5D范圍內,XCC樁各角度下的徑向擠土位移逐漸減小趨于一致,基本呈圓形均勻擴張。

圖8 沉樁完成后徑向位移對比圖

(2)各角度下XCC樁與圓形樁豎向位移比較分析

圖9為沉樁完成后豎向位移對比。由圖9可知,沉樁結束后,其豎向位移變化規律與室內模型試驗變化規律基本一致,兩者數值差值在10%以內。

從圖9a,9b可知,隨著距樁軸線距離的增大,樁周土體豎向位移逐漸減小;隨著深度的增加,距樁軸線同一距離的豎向位移呈先增大后逐漸減小的分布規律,其豎向位移最大值出現在距土體表面約5D處;在相同截面面積下,XCC樁的豎向位移為圓形樁的1.0~1.4倍,與室內模型試驗結果基本一致。

圖9 沉樁完成后豎向位移對比

圖9c為XCC樁豎向位移等值線圖,其分布規律與圖8b徑向位移等值線圖分布規律基本一致,相同位置處土體豎向位移約為徑向位移的2~5倍。

2.3 近樁區球體顆粒運動規律

由圖10可知,沉樁完成后,近樁區土顆粒在豎向上呈明顯的波帶狀運動規律,在樁土界面處土體顆粒運動尤為明顯;在樁端底部,土體顆粒位移較小,顆粒運動規律與室內試驗一致。

圖10 近樁區顆粒運動規律

在圓形樁沉樁結束后,樁周土體顆粒以徑向運動為主,存在一定的切向運動;而XCC樁沉樁完成后,其樁周土體顆粒切向運動較圓形樁顯著,主要是土顆粒受十字段擠壓,致使其向兩側凹弧段運動,表現為1.25D≤r<2.5D范圍內凹弧段土顆粒徑向位移大于十字段土顆粒徑向位移,與室內試驗結果一致。

3 理論對比分析

3.1 傳統圓孔擴張理論

采用柱孔擴張理論分析,假設土體是各向同性的理想彈塑體,具有各項同性的初始應力,只有徑向位移,無切向位移,體力不計。以彈性力學中極坐標下的物理方程和邊界條件,推導在Mohr-Coulomb材料屈服條件下的位移場解答。

平衡方程:

(1)

式中:σr為徑向擠土應力;σθ為環向擠土應力;r為擴孔后土體的徑向位置。

幾何方程:

(2)

物理方程(胡克定律):

(3)

式中:εr為徑向正應變;εθ為環向正應變;E為楊氏彈性模量;μ為泊松比。

塑性區服從Mohr-Coulomb屈服準則:

σr-σθ=(σr+σθ)sinφ+2ccosφ

(4)

可得彈性區位移場ur:

(5)

可得塑性區邊界位移:

(6)

其中最大塑性區半徑:

(7)

式中:c為土體粘聚力;φ為土體內摩擦角;Ru為圓孔最終半徑;Rp為塑性區最大半徑;R0為圓孔初始半徑;a為圓孔擴張過程中孔半徑;up為影響區邊界的徑向位移;p為相應的孔內最終壓力值;Δ=0.015[28]。

采用圓孔擴張理論計算,假設土體從小孔半徑為R0開始擴張到半徑為2R0的過程模擬半徑為R的樁的沉樁擠土過程。R與R0的關系可以從沉樁過程中排開土體體積相等的關系推出,即:

(8)

在進行XCC樁擠土效應計算時,可將其轉化為同一等效直徑下的圓形樁按照現有圓孔擴張理論計算,但是其擠土位移偏于保守,需要對其進行修正。

3.2 修正擴張理論

周航等[21]基于室內透明土沉樁模型試驗得出的沉樁擠土位移場分布規律給出了關于XCC樁徑向擠土位移的修正圓孔擴張理論解答。但其制備的透明土樣與真實土體存在一定差異,其試驗規律與真實土體中沉樁擠土位移場分布規律存在一定出入。周航給出的XCC樁徑向擠土理論解答僅是考慮了在彈性情況下XCC樁各角度下截面半徑與徑向位移之間的關系,其計算結果較大,與實際工程中的沉樁擠土效應存在一定的差異。

由圖8可知,XCC樁徑向擠土位移場在四個象限內位移情況相同,且沿45°,-45°對稱,因此,重點研究了在0°~45°內位移場變化情況,研究發現:

(1)r<1.25D范圍內,XCC樁徑向擠土位移與其各角度下的截面半徑成正比,且位移場呈由十字形擴張逐漸轉為圓形擴張;

(2)在1.25D≤r<2.5D范圍內,XCC樁徑向擠土位移與其各角度下的截面半徑成反比,且由圓形擴張逐漸轉為方形擴張(r=1.6D),又由方擴張逐漸轉為圓形均勻擴張;

(3)在r≥2.5D范圍內,XCC樁各角度下的徑向擠土位移基本保持一致,呈圓形均勻擴張。

上述位移場表明,在XCC樁擴張過程中,其徑向擠土位移場存在2個過渡區(過渡區半徑Rt=1.25D, 2.5D),即存在2個圓形擴張過程,而非周航[21]的一個過渡區(Rt=1.5倍外包圓直徑)。鑒于XCC樁在飽和軟土中的擠土擴張過程,文章基于試驗和數值分析數據的統計分析,同時引入圓孔擴張理論,給出了更加符合真實土體的XCC樁徑向擠土位移經驗公式。

通過研究發現在一定范圍(r≤1.6D)時,各角度下的徑向位移與XCC樁截面半徑服從如下的變化規律,現基于XCC樁擠土機理推導[21,22]并結合數理統計給出XCC樁在各角度下的徑向位移:

(9)

式中:k為通過試驗及數值模擬結果得出的修正系數,即將式(9)計算結果與室內試驗數據、數值模擬試驗數據同時擬合后獲得擬合精度最佳時確定的修正系數k值,經擬合后,取k=3時擬合效果最佳。

下面將給出該函數的推導:

(1)r0求解:

區域1及區域3,分別如式(10),(11)所示:

(10)

(11)

圖11 XCC樁截面邊界幾何關系[21]

區域2:

r0=a[k2(cosθ+sinθ)-

(12)

(13)

(14)

(15)

(2)對n(θ)的求解

(16)

基于圓孔擴張理論的體積守恒原則,可用相同截面面積的圓形樁代替XCC樁截面面積Axcc,故可寫成如下表達式:

(17)

式中:rd為XCC樁截面等效半徑;Rt為過渡區半徑,Rt=1.25D;urt為過渡區半徑的徑向位移;Δ為塑性區平均體積應變,取0.06。

(18)

再將上式帶入式(9),兩邊化簡可得:

(19)

再對兩邊取對數化簡可得:

(20)

便可得:

(21)

這樣便求出了在r≤1.6D范圍內的關于n(θ)的XCC樁在各個角度下的擠土位移情況。

通過對圖8的研究發現,對于r>1.6D范圍內的XCC樁在各個角度下的徑向擠土位移場已經由方形擴張逐漸轉為圓形擴張,為了方便計算,本文認為在r>1.6D范圍內的XCC樁各角度下的擠土位移均服從一同種衰變曲線,但由于數值的不斷減小,導致各角度下的徑向位移基本一致,最終成圓形擴張(略去Rt=2.5D的過渡區)。鑒于此,給出了如下的表達式:

(22)

通過既有研究表明,距樁軸線20D時徑向擠土位移約為1%rd,如姚孟洋[28]通過理論計算、數值模擬表明距樁軸線20D時,徑向擠土位移約為1%rd,因此對于常數n的確定,本文以距離樁軸線20D時的徑向擠土位移為1%rd,故可確定式(22)中n=1.15。

由圖12可見,經修正擴孔理論計算的XCC樁各角度下的徑向擠土位移與模型計算結果和試驗數據基本相同,與兩者相差在5%以內,表明擬合推導的XCC樁擠土效應理論公式可行,可用于工程實踐。

圖12 距土體表面20 cm處徑向位移對比

4 結 論

本文根據室內沉樁模型試驗結果進行計算參數反演,并依據反演參數據建立了FLAC-PFC耦合數值模形,研究XCC樁樁周土體位移場,提出XCC樁擠土位移公式,得到了如下結果結論:

(1)文章提出的uXr=rX(θ)[rX(θ)/(kr)]n(θ),r≤1.6D;uXr=uXr=1.6D(1.6D/r)n,r>1.6D修正擴孔理論與沉樁數值模擬結果及試驗結果吻合度較高,三者數據相差在5%以內,表明提出的修正擴孔理論可行。

(2)XCC樁沉樁數值模擬結果表明 XCC 樁的沉樁擠土位移場可以分為三塊區域: 1)距樁軸線r<1.25D范圍內的近樁段樁周土體呈明顯不均勻的十字形擴張區;2)1.25D≤r<2.5D范圍內的方形過渡區,在這一區域,樁周土體呈現出均勻與不均勻交替擴張出現;3)r≥2.5D范圍內的圓形擴張區,在這一區域,樁周土體保持均勻的圓形擴張。

(3)在同一截面面積下,XCC樁的徑向擠土位移約為圓形樁的1.0~1.4倍,且XCC樁的擠土位移場與圓形樁的擠土位移場在沿徑向分布上存在差異。

(4)在飽和軟土中,由于樁土界面摩擦力較大,土體較軟,樁周土體以凹陷為主,在近樁區1.25D范圍內的土體凹陷量為1.5~3.5 mm。

(5)為了完善XCC樁的擠土效應研究,后續將進一步研究XCC樁的群樁擠土效應和XCC樁沉樁過程中超孔隙水壓力產生及消散規律。

猜你喜歡
模型
一半模型
一種去中心化的域名服務本地化模型
適用于BDS-3 PPP的隨機模型
提煉模型 突破難點
函數模型及應用
p150Glued在帕金森病模型中的表達及分布
函數模型及應用
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 成人午夜免费观看| 白浆免费视频国产精品视频| 亚洲男人在线| 美女毛片在线| 亚洲综合中文字幕国产精品欧美| 国产成人毛片| 精品久久久无码专区中文字幕| 久久网欧美| 亚洲无限乱码| 国产一级毛片高清完整视频版| 日韩AV手机在线观看蜜芽| 成人一级免费视频| 国产美女91呻吟求| 亚州AV秘 一区二区三区| 综合久久五月天| 99精品国产电影| 免费欧美一级| 最新国产麻豆aⅴ精品无| 国产综合色在线视频播放线视| 午夜日b视频| 欧美日韩第三页| 日韩午夜伦| 亚洲an第二区国产精品| 国产女同自拍视频| 国产精品女主播| 91外围女在线观看| 国产亚洲视频中文字幕视频| 中文字幕在线观看日本| 久久99热这里只有精品免费看| 国产91小视频| 亚洲国产中文在线二区三区免| 成年av福利永久免费观看| 狠狠干欧美| 欧美黄色网站在线看| 成人亚洲国产| 亚洲91精品视频| 2020国产精品视频| 极品国产一区二区三区| 国产毛片基地| 日本精品视频一区二区| 67194亚洲无码| 18禁黄无遮挡网站| 成人看片欧美一区二区| 久久人体视频| 538国产在线| 一级毛片在线免费视频| 狠狠亚洲五月天| aaa国产一级毛片| 日本精品中文字幕在线不卡| 99久久国产综合精品2020| 久久香蕉国产线看观看式| 免费无码网站| 国产永久免费视频m3u8| 国产v欧美v日韩v综合精品| 在线日韩一区二区| 99热这里只有精品免费| 欧美国产视频| 91精品国产自产在线老师啪l| 国产男人天堂| 国产精品分类视频分类一区| 欧美日本在线一区二区三区| 不卡视频国产| 免费一看一级毛片| 亚洲欧美在线综合一区二区三区| 日韩第九页| 国产福利拍拍拍| 欧美一区二区啪啪| 国产高清国内精品福利| 国产一区二区影院| 一级毛片在线免费视频| 欧美日韩综合网| 亚洲aⅴ天堂| 亚洲av无码片一区二区三区| 国产成人91精品| a毛片免费看| 天天做天天爱天天爽综合区| 91成人在线观看视频| 欧美a在线看| 美女无遮挡拍拍拍免费视频| 欧美三級片黃色三級片黃色1| 2021国产乱人伦在线播放| 欧美a网站|