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

飽和粘土中球孔擴張問題彈塑性解析

2014-06-15 17:18:12李鏡培唐劍華張亞國鐘光玉同濟大學土木工程學院地下建筑與工程系20092上海巖土及地下工程教育部重點實驗室同濟大學200092上海上海南匯建工建設集團有限公司20199上海
哈爾濱工業大學學報 2014年12期
關鍵詞:模型

李鏡培,唐劍華,張亞國,鐘光玉(1.同濟大學土木工程學院地下建筑與工程系,20092上海;2.巖土及地下工程教育部重點實驗室(同濟大學),200092上海;.上海南匯建工建設(集團)有限公司,20199上海)

飽和粘土中球孔擴張問題彈塑性解析

李鏡培1,2,唐劍華1,2,張亞國1,2,鐘光玉3
(1.同濟大學土木工程學院地下建筑與工程系,20092上海;2.巖土及地下工程教育部重點實驗室(同濟大學),200092上海;3.上海南匯建工建設(集團)有限公司,201399上海)

為了研究靜力觸探試驗及沉樁擴孔等工程問題,基于修正劍橋模型,推導了不排水條件下球孔擴張問題的半解析解.將擴張球孔周圍土體分為臨界狀態區、塑性區以及彈性區三個區域.彈性區內,利用彈性理論得到應力和孔隙水壓力的解答;臨界狀態區及塑性區內,利用相關聯的流動法則、拉格朗日分析法建立了關于應力的一階非線性常微分方程組,以彈塑性界面處的應力分量作為初值,求解微分方程組可得到應力和孔隙水壓力的解答.研究結果表明:各向同性超固結比對擴孔壓力、土體應力、超孔隙水壓力以及塑性區范圍均具有顯著影響,且擴孔過程中土體剪切模量并非常量,其隨擴孔半徑、各向同性超固結比的變化而變化;同時通過與已有解答進行比較,對本文方法的可靠性進行了驗證.

球孔擴張;剪切模量;修正劍橋模型;各向同性超固結比

球孔擴張理論在旁壓試驗、靜力觸探、沉樁及壓力注漿、樁基承載力等巖土工程問題中有著廣泛應用[1].然而,由于采用的本構關系不同,所得到的結果也各不相同.Vesic[2]考慮塑性區土體的壓縮性,給出了土體服從M-C屈服準則狀態下球孔擴張的近似解;Yu等[3]假定土體為線性理想彈塑性體,采用非關聯流動法則的M-C屈服準則,在考慮土體剪脹性的條件下求解了大應變情況下球形孔的擴張;Banerjee等[4]把理想的剛塑性模型應用于正常固結粘性土的擴孔問題;然而理想的彈塑性模型與剛塑性模型均不能考慮土體應力歷史的影響,而劍橋模型則可以克服這些缺陷. Collins等[5]采用臨界狀態模型推導了不排水情況下球形孔擴張的大應變解析解,分析了超固結比對擴孔壓力以及對超孔壓的影響;Cao等[6]采用修正劍橋模型對不排水狀態下球孔擴張問題進行解析,但其在求解過程中假設塑性區內的偏應力為極限狀態的偏應力,剪切模量為常量,使得結果存在一定誤差.胡偉等[7]基于劍橋模型推導了球孔不排水擴張的解析解,但其對偏應力做了線性插值的近似處理,與土體塑性區內偏應力非線性變化的特征不相吻合.Chen等[8]在不對偏應力和平均有效應力進行任何簡化假設的情況下,利用修正劍橋模型給出了圓柱孔擴張問題的精確半解析解;然而,研究沉樁擴孔,靜力觸探等問題時,往往認為樁端或探頭周圍土體呈球形擴張[9],故而基于平面應變假設的圓柱形擴孔理論并不能直接使用.

以往的球孔擴張理論假設在擴孔過程中,剪切模量為常量,基于此,本文在以往研究基礎上,利用修正的劍橋模型理論,對球孔擴張問題進行解析.得到了飽和粘土中球形孔擴張后,彈性區、臨界狀態區、塑性區內的應力及超孔隙水壓力分布特征,給出了球孔周圍土體剪切模量的變化規律,分析了各向同性超固結比對球孔周圍應力以及超孔壓的影響.并將球孔擴張理論應用到靜力觸探實驗中,本文的研究結果可為靜力觸探等試驗提供一定的理論基礎.

1 力學模型及基本假設

采用如圖1所示球孔擴張模型,球孔初始孔徑為a0,初始孔隙水壓力為u0,初始平均總應力為p0,初始平均有效應力為p0′,經過擴張后,最終孔徑為a.假設壓應力和壓應變為正,任意一點平均總應力為p,有效平均總應力為p′,偏應力為q,孔隙水壓力為u.rf為塑性區與臨界狀態區交界面的半徑,rp為彈塑性區交界處的半徑,其初始半徑為rp0,對于塑性區內任意一點rx,其初始位置為rx0.對于a≤r≤rf區域內的土體,此區域內的土體的應力狀態均處于臨界狀態CSL線上,土體達到此狀態后,應力狀態不再變化,該區域為臨界狀態區;當rf≤r≤rp,土體處于塑性狀態,服從修正的劍橋模型;當r≥rp,土體處于彈性狀態,服從虎克定律.假定土體飽和、均質、各向同性、不可壓縮,將球孔擴張看成不排水過程.

根據彈塑性理論可得

球孔擴張前土體處于各向同性狀態,可得到初始徑向應力σr0,初始切向應力σθ0:

σ0為初始應力,球孔擴張過程中,土體內部的任意一點都滿足下列平衡方程:

若用有效應力表示,則方程(3)可化為

式中:σr、σθ分別為徑向、切向總應力,σ′r、σ′θ分別為徑向、切向有效應力,u為孔隙水壓力.

圖1 球孔擴張模型

2 彈性區應力分析

彈性區內,采用小應變理論,假設壓應變為正,因此徑向應變增量與切向應變增量可由式(5)表示:

d ur為徑向位移增量,由虎克定律可得彈性區的應力應變關系:

式中:ν為泊松比,E為彈性模量,σr′、σθ′分別為徑向和切向有效應力.對于修正的劍橋模型,泊松比ν為常數.

G為剪切模量,可由比體積υ、平均有效應力p′表示[10]:

式中κ為υ-ln p′平面上卸載-再加載曲線的斜率.由于球孔擴張過程可看成不排水擴張,故彈性區的體積變化為零,因此有:

由式(9)可得,彈性區內的有效應力p′和比體積υ保持不變,因此在彈性區內,剪切模量G為常數.

由平衡方程、應力應變關系可得彈性區內的應力和位移[11]:

式中:σp為彈塑性邊界處的總徑向應力,G0為初始剪切模量.由式(10)可得彈性區的平均應力保持不變,又因為彈性區的平均有效應力保持不變,因此彈性區內的孔壓也保持不變.

3 塑性區應力分析

3.1 有效應力

塑性區內,土體單元服從修正劍橋模型,其屈服函數[10-12]為

式中:M為p-q′平面中CSL線的斜率,pc′為各向同性狀態條件下的屈服應力.

利用彈塑性邊界處應力連續性條件可得

由式(11)、(12)可知,在彈塑性邊界處,偏應力為[13]

式中R為各向同性超固結比,R=pc′/p0′[14].

聯立式(1)、式(10)可得彈塑性邊界處的有效應力:

根據塑性理論的相關聯流動法則可得塑性應變增量:

式中:ψ為塑性流動因子,η=q/p′,定義為應力比.

塑性區土體的應變服從大變形理論,由塑性理論可知

球孔周圍土體任意一點的體應變為零,即dεv= 0,因此有

根據chen等[8]的方法,結合式(16)、(18)可得到塑性區內關于應力的微分方程:

式(19)是運用拉格朗日分析法建立的微分方程組,該方程組適用于塑性區內任意一點rx,若要通過式(19)求得rx處的應力,只需求得rx處土體單元由彈性狀態變為塑性狀態的位置rxp,以及在rxp處的應力初始值.下面論述中將給出應力初始值.

rxp處的應力初始值與徑向距離r無關[8],可由式(14)、(15)確定,即

球孔擴張過程可看作不排水過程,球孔周圍的土體體積應變為零,結合式(10)可得rxp表達式

以彈塑性界面處一點為研究對象,則rx、rxp等于rp,由式(22)可得

因此,塑性區內球孔擴張問題可歸結為求解一系列具有初值條件的非線性常微分方程,其中式(19)為控制方程,式(21)、(23)為初值條件;該非線性微分方程組可通過數值方法求解.

3.2 超孔隙水壓力

擴孔過程中,在彈性區內,由于孔壓保持不變,不產生超孔壓;在塑性區內,超孔壓可通過對式(4)積分得到,積分區間為[rxrp],則

由于無法得到超孔隙水壓力的解析解,故可利用數值積分求解式(24).

4 應力應變關系

與偏應力相對應,偏應變可由下式確定[15]:

在彈性區內,偏應力q可由式(10)確定,偏應變可由式(26)確定,因而聯立式(10)、式(26)可確定彈性區內的應力應變關系;在塑性區內,通過求出微分方程(19)的數值解后,可得塑性區內土體的偏應力,結合式(30)可得塑性區土體的應力應變關系.

5 算例分析

圖2為擴孔半徑a/a0=2,各向同性超固結比R=1.001、2、3、10時,偏應力q的變化規律,隨著各向同性超固結比的增大R,臨界狀態區的半徑逐漸減小.當R<2,隨著徑向距離的減小,偏應力q逐漸增大,當土體達到臨界狀態時,偏應力q保持不變,土體在屈服之后表現硬化的性質;當R= 2,土體一旦屈服,偏應力q就保持不變,土體在屈服之后表現理想彈塑性的性質;當R>2,隨著徑向距離的減小,偏應力q先增大后減小,當土體達到臨界狀態時,偏應力q保持不變,土體在屈服之后表現軟化的性質.在塑性區內,各向同性超固結比對偏應力q的影響較大,隨著各向同性超固結比的增大,偏應力q增大,在彈性區內各向同性超固結比對偏應力q的影響較小,可以忽略不計.與文獻[6]的計算方法進行對比,可以看出兩種方法得到的偏應力吻合較好.

圖2 偏應力q/cu的徑向分布規律

當球孔擴展到某一孔徑時,土體中不同位置處的有效應力大小也不同,由式(8)可知剪切模量G會隨之發生變化.圖3為a/a0=2時,孔周土體剪切模量G的分布規律,可以看出,臨界狀態區與彈性區內土體的剪切模量幾乎不變,而塑性區內的剪切模量G變化較大.各向同性超固結比R對剪切模量G有較大影響,當R<2時,隨著徑向距離的增大,剪切模量先保持不變,然后逐漸增大,最后保持不變;當R=2時,隨著徑向距離的增大,剪切模量G保持不變;當R>2時,隨著土體徑向距離的增大,剪切模量G先保持不變,然后逐漸減小,最后保持不變;從圖3還可得到,在臨界狀態區、塑性區內,各向同性超固結比越大,剪切模量G越大.當R<2時,隨著擴孔半徑的增大,剪切模量G先保持不變,然后逐漸減小,直至穩定;當R=2時,隨著擴孔半徑的增大,剪切模量G保持不變;當R>2時,隨著擴孔半徑的增大,剪切模量G先保持不變,然后逐漸增大,最后保持不變.由此可見,擴孔過程中土體的剪切模量G并非常量,而是隨著球孔的擴孔半徑、土體的各向同性超固結比變化而變化.

圖3 剪切模量G/G0與徑向距離r/a的關系

圖4 中,當1≤a/a0≤2時,隨著擴孔半徑的增加,擴孔壓力急劇增大;當a/a0≥2時,隨著擴孔半徑的增加,擴孔壓力幾乎不變,表明擴孔壓力達到極限值.各向同性超固結比對擴孔壓力的影響也比較大,隨著各向同性超固結比的增大,擴孔壓力也急劇增加.

圖4 擴孔壓力σa/cu與擴孔半徑a/a0的關系

圖5 為擴孔半徑a/a0=2,各向同性超固結比R=1.001、3、10時,超孔壓Δu徑向分布規律.隨著各向同性超固結比的增大,孔壁處的超孔壓逐漸增大,而彈塑性邊界處的超孔壓減小,當各向同性超固結比增大到一定程度時,在彈塑性邊界處將出現負孔壓,當R=10時,彈塑性邊界處超孔壓為負.在彈性區內,土體的超孔壓為零.文獻[6]的方法計算得到的超孔壓與本文結果有一定差異.當R=1.001時,兩者差異較大.這是因為擴孔過程中的剪切模量為初始剪切模量,塑性區的偏應力為臨界狀態偏應力,而本文認為擴孔過程中剪切模量與有效應力成正比,且對塑性區內的偏應力沒有簡化,故本文得到的超孔壓更準確.同時也說明文獻[6]的假設對超孔壓的影響較大.

根據圖6,當1≤a/a0≤2時,隨著擴孔半徑的增加,孔壁處的超孔壓Δu急劇增加;當a/a0≥2時,隨著擴孔半徑的增加,孔壁處的超孔壓保持不變.各向同性超固結比對擴孔壓力的影響也比較大,隨著各向同性超固結比的增大,孔壁處的超孔壓急劇增加.當各向同性超固結比較大時,孔壁出現負孔壓.這主要是因為當各向同性超固結比較大時,土體表現出剪脹的特性.

圖5 超孔壓Δu/cu與徑向距離r/a的關系

圖6 超孔壓Δu/cu與擴孔半徑a/a0的關系

各向同性超固結比R對塑性區的半徑有一定的影響.圖7表明:當1≤R≤3時,隨著各向同性超固結比的增大,塑形區的半徑急劇減小;當R>3時,塑性區的半徑基本保持不變.

圖7 塑性區半徑rp/a與R的關系

圖8 為當擴孔半徑a/a0=2時,徑向應力σr與切向應力σθ隨徑向距離的變化規律.孔壁附近土體的徑向應力與切向應力基本保持不變,這表明球孔周圍土體已達到臨界應力狀態.塑性區內的應力急劇增大或者減小.在彈性區內,隨著徑向距離的增大,應力趨于穩定.此外,各向同性超固結比對徑向應力、切向應力也有顯著的影響,隨著各向同性超固結比的增大,臨界狀態區內的應力急劇增大,然而在彈性區內,徑向應力與切向應力幾乎不受各向同性超固結比的影響.

圖8 徑向應力σr/cu、切向應力σθ/cu的徑向分布規律

圖9 為球孔周圍土體中應力應變關系.當偏應變εq較小時,此時偏應力q與偏應變εq呈線性關系,此時土體發生彈性變形;隨著球孔的擴張,偏應變εq逐漸變大,偏應力q與偏應變εq呈非線性關系,此時土體已發生塑性變形.此外,各向同性超固結比對土體的應力應變關系具有顯著的影響,當各向同性超固結比R=1.001時,土體在屈服之后表現出硬化的性質;當R=2時,土體在屈服之后表現理想彈塑性的性質;當R=10時,土體在屈服之后表現軟化的性質.

圖9 偏應力q與偏應變εq的關系

6 應用分析

球孔擴張理論可用于樁基承載力、旁壓實驗、靜力觸探實驗.本文以靜力觸探實驗為例,分析靜力觸探實驗錐頭的極限阻力.文中不考慮錐頭的粗糙程度,因此根據錐頭的靜力平衡方程可得錐頭的極限阻力為

式中:qt為錐頭的極限阻力,σu為球孔的極限擴孔壓力,基于本文理論分析,只需令a/a0→∞,即可得到極限擴孔壓力σu.

為了驗證該理論模型在實際應用中的可行性,以下將以具體的靜探實驗為例.本文選取文獻[16]的數據:M=1.37、v=0.3、R=1.34~3.00、p0′=23.2~104 kPa,實測圓錐阻力qc=204~763 kPa.錐頭貫入過程中,由于超孔壓的影響,會使實測超孔壓偏小,因此應考慮孔壓影響,修正后的圓錐阻力q[17]為

式中:qc為實測圓錐阻力,qt為修正后的錐頭極限阻力,u為錐頭孔壓,α為凈面積比值,根據文獻[17]的研究,α可取為0.84.

從圖10可看出,利用本文理論方法得到的數據與文獻[16]的試驗數據雖然有一定的誤差,但總體趨勢一致,在一定程度上,可利用該方法預測試驗數據,因此本文理論模型有一定的實用價值.

圖10 錐頭極限阻力深度分布規律

7 結 論

采用修正劍橋模型,在不對偏應力做任何假設的條件下,得到了不排水條件下球孔周圍土體應力和孔隙水壓力的半解析解.并通過與已有解答的對比分析說明了本文研究方法的正確性以及結果的可靠性.研究結果表明:

1)各向同性超固結比R對球孔周圍土體剪切模量的影響顯著,R<2時,隨徑向距離的增大,剪切模量先保持不變,后逐漸增大直至穩定;R=2時,剪切模量為定值;當R>2時,剪切模量先保持不變,后逐漸減小直至穩定;R對臨界狀態區、塑性區內的應力及超孔隙水壓力影響較大,R越大,切向和徑向應力越大;彈性區內,R對孔隙水壓力、徑向和切向應力幾乎無影響.

2)各向同性超固結比對臨界狀態區、塑性區內的應力、超孔壓影響較大,各向同性超固結比越大,切向應力、徑向應力越大;在彈性區內,各向同性超固結比對孔壓、徑向應力、切向應力幾乎無影響.此外,各向同性超固結比越大,塑性區的半徑、臨界狀態區的半徑越小.

3)將球孔擴張理論應用到靜力觸探原位實驗中,對球孔擴張理論的應用性進行了說明,可以在一定程度上促進靜力觸探原位實驗在實際工程中的應用.

[1]宋勇軍,胡偉,王德勝,等.基于修正劍橋模型的擠密樁擠土效應分析[J].巖土力學,2011,32(3):811-814.

[2]VESICA S.Expansion of cavity in infinite soilmass[J]. Journal of the Soil Mechanics and Foundation Division,ASCE,1972,98(SM3):265-290.

[3]YU H S,HOULSBY G T.Finite cavity expansion in dilatant soils:loading analysis[J].Geotechnique,1991,41(2):173-183.

[4]BANERJEE P K,DAVIES T G,FATHALLAH R C. Behavior of axially loaded driven piles in saturated clay from model studies[M].England:Applied Science Publishers Ltd,1982.

[5]COLLINS IF,YU H S.Undrained cavity expansions in criticalstate soils[J].Int JNumer Anal Meth Geomech,1996,20(7):489-516.

[6]CAO L F,THE C I,CHANG M F.Undrained cavity expansion in modified Cam clay I:theoretical analysis[J].Geotechnique,2001,51(4):323-334.

[7]胡偉,黃義,劉增榮,等.飽和粘土中擠土樁球形孔擴張的彈塑性分析[J].工程力學,2008,25(8):180-187.

[8]CHEN S L,ABOUSLEIMAN Y N.Exact undrained elasto-plastic solution for cylindrical cavity expansion in modified Cam Clay soil[J].Geotechnique,2012,62(5):447-456.

[9]CHANG M F,TEH C I,CAO L F.Undrained cavity expansion in modified Cam clay II:application to the interpretation of the piezocone test[J].Geotechnique,2001,51(4):335-350.

[10]WOOD D M.Soil behaviour and critical state soil mechanics[M].Cambridge:Cambridge University Press,1990.

[11]YU H S.Cavity expansion methods in geomechanics[M].Netherlands:Springer,2000.

[12]李廣信.高等土力學[M].北京:清華大學出版社,2004.

[13]WROTH C P.The interpretation of in situ tests[R]. Twenty Fourth Rankine Lecture Geotechnique,1984,34(4):449-489.

[14]CHANG M F,TEH C I,CAO L F.Critical state strength parameters of saturated clays from the modified Cam clay model[R].Canadian Geotechnical Journal,1999,36(5):876-890.

[15]SILVESTRI V,ABOUSAMRA G.Application of the exact constitutive relationship of modified Cam clay to the undrained expansion of a spherical cavity[J]. International Journal for Numerical and Analytical Methods in Geomechanics,2011,35(1):53-66.

[16]HIGHTDW,BOND A J,LEGGE JD.Characterization of the Bothkennaar clay:an overview[J].Géotechnique,1992,42(2):303-347.

[17]張誠厚,施健.孔壓靜力觸探試驗的應用[J].巖土工程學報,1997,19(1):50-57.

(編輯 趙麗瑩)

Elastic-plastic solution of sphere cavity expansion in saturated clay

LIJingpei1,2,TANG Jianhua1,2,ZHANG Yaguo1,2,ZHONG Guangyu3
(1.Department of Geotechnical Engineering,Tongji University,200092 Shanghai,China;2.Key Laboratory of Geotechnical and Underground Engineering(Tongji University),Ministry of Education,201399 Shanghai,China;3.Shanghai Nanhui Construction Group CO.LTD.,201399 Shanghai,China)

An exact semi-analytical solution in the undrained cavity expansion can be obtained on the basis of the MCCmodel to research the cone penetration test and the pile driving.The field around the cavity can be divided into three zones:critical zone,plastic deformation zone and elastic deformation zone.In the elastic zone,an analytical solution for the distributions of stress and excess pore pressure is deduced according to the elastic theory.In the critical and plastic zone,a set of first-order nonlinear ordinary differential equations concerning stress can be obtained according to the associated flow rule and the lagrangian analysismethod.The stressss and pore pore pressure can be solved as an initial value problem starting at the elastic-plastic boundary.The results show that the isotropic over consolidation ratio has a significant influence on the stresses and the excess pore pressure.The shearmodulus vary significantly with the cavity radius and the isotropic over consolidation in the course of cavity expansion.

sphere cavity expansion;shearmodulus;modified cambridgemodel;isotropic over consolidation ratio

TU473

A

0367-6234(2014)12-0071-07

2014-05-05.

國家自然科學基金(41272288);浦東新區科技發展基金創新資金項目(PKJ2013-C08).

李鏡培(1960—),男,教授,博士生導師.

李鏡培,lijp2773@tongji.edu.cn.

猜你喜歡
模型
一半模型
一種去中心化的域名服務本地化模型
適用于BDS-3 PPP的隨機模型
提煉模型 突破難點
函數模型及應用
p150Glued在帕金森病模型中的表達及分布
函數模型及應用
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 亚洲无码视频图片| 亚洲三级色| 乱人伦视频中文字幕在线| 亚洲成人网在线播放| 成年人国产网站| 91精品亚洲| 欧美成人一区午夜福利在线| 亚洲最大看欧美片网站地址| 日韩精品无码一级毛片免费| 成人91在线| 亚洲无码高清视频在线观看| 日韩精品无码不卡无码| 91亚洲免费视频| 91国内外精品自在线播放| 毛片a级毛片免费观看免下载| 2021国产精品自产拍在线| 午夜精品久久久久久久无码软件| 人妻无码中文字幕第一区| 国禁国产you女视频网站| 日韩无码视频网站| 免费在线视频a| 国产成人亚洲精品色欲AV| 免费一级成人毛片| 国产永久在线视频| 日韩区欧美区| 毛片网站免费在线观看| 亚洲视频四区| 久久国产亚洲欧美日韩精品| 国模沟沟一区二区三区| 色婷婷天天综合在线| 成人无码区免费视频网站蜜臀| 亚洲中文精品人人永久免费| 成人日韩视频| 精品午夜国产福利观看| 超清人妻系列无码专区| 青青久久91| 国产成人精品男人的天堂下载| 久久黄色免费电影| 呦视频在线一区二区三区| 福利在线免费视频| 国产区福利小视频在线观看尤物| 熟妇人妻无乱码中文字幕真矢织江| 国产一级小视频| 日韩天堂网| 亚洲一区免费看| 曰韩人妻一区二区三区| 免费又爽又刺激高潮网址| 呦女亚洲一区精品| 青青操视频免费观看| www.亚洲色图.com| 亚洲中久无码永久在线观看软件| 日本爱爱精品一区二区| 青青草一区| 久久香蕉国产线| 日韩大片免费观看视频播放| 一级香蕉人体视频| 免费看美女毛片| 免费国产福利| 国产乱子伦一区二区=| 亚洲欧美成人在线视频| 日韩精品毛片人妻AV不卡| 久久特级毛片| 波多野结衣中文字幕一区二区| 亚洲成aⅴ人片在线影院八| 欧美综合区自拍亚洲综合绿色| 57pao国产成视频免费播放| 日韩精品无码免费专网站| 夜夜操天天摸| 无码免费的亚洲视频| 欧美天堂久久| 色悠久久综合| 丝袜无码一区二区三区| 国产精品美人久久久久久AV| 九九这里只有精品视频| 日本一区二区不卡视频| 日韩精品一区二区深田咏美| 丁香亚洲综合五月天婷婷| 一级毛片在线播放免费| 18黑白丝水手服自慰喷水网站| 国产黄在线免费观看| 狠狠色成人综合首页| 亚洲第一成年网|