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

電場-滲流場-應力場耦合的電滲固結數值分析

2012-01-08 07:12:24王柳江劉斯宏汪俊波
巖土力學 2012年6期
關鍵詞:理論

王柳江,劉斯宏,汪俊波,朱 豪

(1. 河海大學 水利水電學院,南京 210098;2. 西安熱工研究院有限公司,西安 710032;3. 南水北調中線干線建設管理局,北京 100038)

1 引 言

電滲是一種采用電能進行地基加固的技術[1]。Casagrande 于1939 年首次將該技術應于巖土工程中[2-3],現已應用到斜坡、堤壩、軟土地基加固等巖土工程中[4-6]。近幾年,國內主要將其應用到圍海造地、袋裝吹填淤泥土加固、基坑降水、高速公路地基處理等工程。

目前,對于電滲法,尤其當電滲法與其他工法聯合使用時,大多數還是依靠工程技術人員的經驗,尚缺乏合適的設計計算理論。為克服電滲法理論上的空白,國內外眾多學者對電滲固結理論展開了研究。其中,Esrig[7]最早對電滲固結理論進行了研究,然而他僅考慮了電滲的作用;Wan 等[8]通過進一步完善使其適用于電滲和堆載共同作用的情況,但只適用于一維情況。Lewis[9]在Esrig 理論的基礎上 考慮了電滲過程中的電流變化。近年來,國內學者也對電滲固結理論進行了深入研究。莊艷峰,王 釗[10-11]研究了電滲過程中的電學問題,考慮了電極界面電阻的變化情況,從能量與電荷守恒角度分別推導了電滲的能級梯度理論和電荷累積理論,使其能夠較好地解釋電滲過程中電勢的分布、電流以及土體電導率的變化規律,然而其適用性還有待于進一步驗證。蘇金強等[12]完善了Esrig 電滲固結理論,使之適用于二維情況,能夠用于電極排形布置的情況,但對復雜電極布置的情況卻不適用。李瑛等[13]通過電場等效的方法建立了等應變條件下考慮堆載-電滲聯合作用的軸對稱電滲固結理論,該理論雖然已經具備了較強的工程實際意義,但還是局限在飽和土的理論框架內。由于電滲是一個涉及水分、電子、離子在土體中遷移的復雜問題,在電勢梯度的作用下產生負孔隙水壓力,同時滲流過程中水流與土骨架變形是相互作用的,屬于非飽和土范疇的滲流-變形耦合問題。因此,電滲固結理論的推導應建立在非飽和土理論的基礎上,并采用有限元方法對電滲以及電滲與其他工法聯合使用的情況進行計算。

本文在綜合已有研究成果的基礎上,基于非飽和土多孔介質力學理論,假設等溫條件,忽略孔隙水體的化學作用以及各組分發生相變的影響,采用電荷守恒原理、質量守恒原理、歐姆定律以及達西定律,建立能夠反映電場、滲流場以及應力場相互作用的電滲固結方程;然后,將該方程通過Galerkin法離散得到有限元計算列式,編制計算程序,對電滲的室內模型試驗進行數值模擬,并與試驗結果進行了比較。

2 電滲固結理論的建立

2.1 基本假定

為推導能夠反映電滲過程中電場-滲流場-應力場耦合的控制方程,需做如下一些假設:

(1)假設土體內各點等溫,不考慮土體內化學、相變等因素的影響;

(2)忽略土體內細顆粒發生電泳而產生的電流;

(3)土體內電荷傳遞滿足歐姆定律;

(4)土體內水體流動滿足達西定律;

(5)由電勢梯度和水頭梯度產生的流量可以疊加;

(6)由土體中離子移動而產生的電場相對外加電場忽略不計;

(7)固體顆粒不可壓縮、水可壓縮且體積壓縮模量不變,孔隙氣體與大氣相通。

2.2 基本控制方程的建立

根據歐姆定律可得

根據電荷守恒原理,同時假設電流處于穩定狀態,則得到如下方程:

式中:eR 為電流匯源項。

將式(1)代入式(2)得到電勢分布的控制方程為

根據達西定律,水頭梯度引起的流速為

式中:hK 為土體的滲透系數(m/s);H 為水頭。

由電勢梯度引起的流速為

式中:eoK 為土體的電滲透系數(m2?V-1?s-1)。

根據假設(5),電滲過程中流過土體單元的流速為水頭梯度和電勢梯度引起的流速之和,則土體中孔隙水的滲流流速可表示為

同時根據假設(7),得到土體單元的質量守恒方程(即滲流連續性方程)為

式中:C 為非飽和土的容水度( C =?θw/?uw),θw為土體的體積含水率,θw= nSr;uw為孔隙水壓力;Sr為飽和度;n 為孔隙率;u 為位移矢量; Kw為水的體積壓縮模量(MPa)。

將式(6)代入式(7),得到滲流場的控制方程為

土體單元的平衡微分方程為

式中:ρ、sρ 、wρ 分別為土體、土體中固體顆粒、土體中水的密度。根據假設(7),土體中的孔隙氣壓等于大氣壓,則此時土體的有效應力可表示為

式中:σi′j為有效應力張量;δij為Kroneker 符號;χ為與飽和度 Sr有關的非線性常數,當土體處于飽和狀態時, χ = 1;非飽和狀態時, χ = χ( Sr),為簡化該關系式也可近似等于rS ,即 rSχ≈ 。

式(3)、(8)、(9)分別為電場、滲流場以及應力場耦合作用下的電滲固結理論的控制方程。為求解方程,將其初始條件及邊界條件表示如下:

Dirichlet 邊界條件:

Neumann 邊界條件:

2.3 土體的非線性本構模型

本文在非飽和土多孔介質理論的基礎上研究電滲作用下土體中電場、滲流場、應力場的相互作用,計算采用的參數除eoK 、hK 和eσ 之外,還包括土水特性以及力學本構模型參數。其中,土體的電導率與含水率、飽和度以及孔隙率的關系密切。根據文獻[14]的試驗結果,土體的電導率采用推廣的阿爾奇公式換算得到:

式中:K、a、b 為與土性有關的參數;wσ 為孔隙水的電導率。

非飽和土的飽和度隨吸力的變化而變化,兩者之間的關系式可以通過擬合土體的特征曲線而得到,本文采用van-Genuchten 模型[15]:

式中:s 為吸力( ua- uw),根據假設(7),土體中的孔隙氣壓等于大氣壓,則吸力等于負孔隙水壓力,其中規定吸力最大不超過104kPa;p0為進氣值;n、m 為曲線擬合參數,m = 1 - 1/n,Se為有效飽和度:

式中:rwS 為殘余水飽和度;swS 為土體的最大飽和度。

土體的滲透系數采用van-Genuchten-Mualem模型[16]計算,主要考慮了非飽和土體的滲透系數與飽和度的關系:

式中:rk 為相對滲透系數,將其與土體飽和狀態時的滲透系數 swK 相乘即為土體的滲透系數,其余參數意義與式(16)、(17)中的相同。

電滲通常用于超軟黏土的排水加固,則土體的變形預測采用能夠較好地反映黏土彈塑性變形特性的修正劍橋模型,其屈服函數見式(19),該模型具有嚴格的理論依據,適用于正常固結和欠固結黏土的應力變形計算。

式中:p 為平均主應力;q 為廣義剪應力;λ、κ 分別為土體的壓縮和回彈系數;e0為初始孔隙比;M為臨界狀態線的斜率;為塑性體積應變,也是該模型的硬化參數。

3 控制方程的離散

Galerkin 加權余量法可以方便地運用于固體力學、流體力學和物理學問題,并且已經發展出了相對標準的推導過程,現使用加權余量法對控制方程在空間域上進行離散,離散后得到的有限元計算矩陣列式為

其中

將矩陣表達式(20)以簡化形式表示為

X 為主要變量,變量上方的點代表時間導數。采用向后歐拉法求解該方程組:

式中:1nX+、1nX+˙ 、1nF+為 1nt t+= 時刻對應的值,采用Newton-Raphson 迭代法對其進行求解。首先,令每一時步計算時變量的初值為

則:

由式(22)可知,當右端項和左端項之差滿足設定的收斂精度時迭代結束:

由于該問題的高度非線性,將式(24)代入式(22)后,得到如下形式:

上標i 為迭代步數,下一迭代步采用的變量更新為

將更新后的變量代入式(25),若誤差滿足收斂要求,則結束該時步的迭代;若不滿足,則代入式(27)繼續進行迭代計算。由于式(20)左端的系數矩陣不對稱,則采用非對稱方程組求解器求解該方程組。

4 算例驗證及分析

根據以上介紹的有限元計算列式和方程計算方法,自主開發了有限元計算程序。為驗證該程序的正確性和適用性,本文進行了電滲聯合真空排水的室內模型試驗,并對其進行了數值模擬。

4.1 電滲聯合真空排水模型試驗

模型試驗箱尺寸和布置情況如圖1 所示。試驗在一個80 cm×80 cm×60 cm 的模型箱中進行,采用的土樣取自大連大窯灣的海相淤泥質軟土,其含水率高達92.3%,土樣填筑高度為50 cm;在模型箱的4 個對角處布置陽極,陰極和塑料排水板布置在模型箱中間位置;陰、陽極電極材料為直徑1 cm 的鋼筋,塑料排水板的截面尺寸為100 mm×4 mm;為防止土體中水分蒸發或入滲,在土體表面鋪設一層土工膜進行隔水;試驗中設置了位移計和振弦式孔壓計(量測范圍:-100~100 kPa)分別監測圖1中s1,s2,s3 處的沉降和k1 處的孔壓。試驗開始后,在初始72 h 內先進行真空排水,爾后接通直流電源進行電滲降水,電滲電壓為10 V,采用穩壓方式通電,試驗總共持續30 d,試驗過程中每隔4 h監測一下相關數據。由于試驗中排水板附近土體干裂經常發生漏氣現象,真空度在試驗中是變化的,由真空表測得其變化曲線如圖2 所示;由于該土為海相淤泥土,具有含鹽量高、電導率大的特點,且通過Miller soil box 進行試驗得到土體電導率在電滲過程中確實會減小,然而在含水率大于82%時,土體電導率的減小并不十分明顯;其次,電極材料采用的是鋼筋,在通電條件下,陽極在該土中極易腐蝕,使其表面形成了一層保護膜,導致陽極的導電率迅速降低,增加了陽極和土樣之間的接觸電阻,進而使有效作用在兩極間土體上的電壓在試驗過程中不斷減小的,因此,可以認為接觸電阻是電滲計算中必須考慮的因素。本文采用萬用表測定有效電壓,將其一探針與陰極連接,另一探針與50 cm 長的銅線連接后作為加長的探針,插入陽極邊上的土體中。試驗中,探針與陽極之間的距離設為4 mm,該值相對試驗土樣尺寸來說較小,因此,探針與陽極之間的電勢差可近似認為由接觸電阻引起。圖3為試驗過程中作用在兩極間土體上有效電壓的實測值。整個試驗相關的結果示于圖5~9 中。

4.2 計算模型及條件

圖1 電滲聯合真空排水模型試驗示意圖 Fig.1 Schematic diagrams of model test of electroosmosis combined with vacuum drainage

圖2 試驗過程中真空度隨時間變化 Fig.2 Variation of vacuum degree along test process

圖3 試驗過程中有效電壓隨時間變化 Fig.3 Variation of effective voltage along time of model test

圖4 模型試驗有限元網格 (單位:m) Fig.4 Finite element meshes of model test (unit: m)

圖5 電流的實測值與計算值比較 Fig.5 Comparison between calculated and measured current

圖6 兩極間土體電阻的實測值 Fig.6 The evolution of soil resistance between two poles

圖7 表面沉降的實測值與計算值比較 Fig.7 Comparison between calculated and measured surface settlements

圖8 點k1 處相對孔隙水壓力實測值與計算值比較 Fig.8 Comparison between calculated and measured relative pore water pressure at node k1

圖9 試驗總排水量的實測值與計算值比較 Fig.9 Comparison between calculated and measured overall drainage

圖4 為模型試驗的有限元網格。為簡化網格,通過綜合陰極導電特性和排水板的滲透特性,將其作為一種材料進行模擬。計算模型的邊界條件為:模型箱底部及4 側為已知位移邊界,其中底部為固定邊界,4 側只約束其法向位移;由于土體表面鋪有土工膜,除排水板頂部單元之外,模型的底部,頂部以及4 側均為不排水邊界,而排水板頂部設為水頭邊界,且該邊界上的孔壓約束值等于圖2 所示的實測真空度;電極頂部單元為電勢邊界,其電勢值輸入如圖3 所示。初始條件包括初始應力和初始水頭,初始應力由土體自重決定;試驗前土體處于完全飽和狀態,初始水頭為表面高度,等于0.5 m。計算所采用的參數包括土體、陽極以及陰極和排水板復合單元的電導率、電滲透系數、滲透系數、土水特征參數以及修正劍橋模型中的參數。其中土體 的電導率,土水特征關系以及力學本構關系如式(15)~(19)所示,其參數通過試驗得到如下:σw=1.0 S/m, Keo=2.5×10-9m2/s?V,K=3.6,a= 2.051,b=5.531, p0=5 kPa,n=1.84, Srw=0.6,Ssw=1.0, Kw=1.95×106kPa, Ksw=2.5×10-9m/s,λ=4.179,κ=0.052,M=0.65,ν=0.3,e0=2.96。陽極以及陰極的電導率為5.0×106S/m,其電滲透系數相對土體可以忽略,取為1.0×10-20m2/ s?V;陽極和陰極的滲透系數也取為常數,由于陰極與排水板在模擬時作為一種材料處理,則將排水板的滲透系數等效到陰極后為2.0×10-4m/s,陽極的滲透系數遠小于陰極,假設等于淤泥土的飽和滲透系數為2.5×10-9m/s;由于陰極和陽極的尺寸相對較小,對土體的變形影響可以忽略,因此,假設陰極和陽極采用的力學本構模型參數與土體一致。

4.3 計算結果分析

圖5 為電滲過程中電流隨時間變化的計算值與實測值的對比,兩者十分接近。由圖可知,電流隨電滲時間的增長而減小,且其變化趨勢與圖3 中有效電壓的變化趨勢一致。這是由于電滲過程中陽極腐蝕以及陽極附近土體含水率急劇降低導致接觸電阻增大,而由于土體初始含水率較高,且試驗過程中平均含水率變化幅度較小,測得土體總電阻變化也較小,如圖6 所示,由此說明了在該試驗中有效電壓的減小主要受接觸電阻的影響。根據歐姆定律中電壓和電流的關系可知,電流隨電壓的減小而減小;在每次間歇期之后通電,電流較通電前會有顯著增大,這與間歇期土體中含水率、離子重分布以及累積電荷消散有關[17]。

圖7 為試驗過程中土體表面s1-s3 測點(參見圖1)沉降的計算值與實測值比較。由于s2 處位移計失效,所以缺少了s2 處的實測值。由圖可見,試驗結束后,陰極表面沉降的計算值和實測值分別為5.07 cm 和4.98 cm,而陽極處表面沉降的計算值和實測值為5.29 cm 和5.24 cm,計算值與實測值基本吻合。在試驗初始72 h 內,由于未開始電滲,陰極處沉降明顯大于陽極;待土體中開始通電后,由于在電滲作用下陽極附近土體中的水分往陰極流動,引起陽極處的沉降速度加快,且試驗結束后其最終沉降量大于陰極。

圖8為k1處減去初始孔壓值后的相對孔隙水壓力的實測值與計算值的對比。在試驗開始的72 h,由于曼德爾效應,k1 處孔壓略有增大;電滲開始后,k1 處孔壓逐漸減小且轉變為負孔壓,而在電滲間歇期,孔壓略有增大。該現象說明了電滲能夠使陽極附近產生負孔壓,從而起到加固效果。

圖9 為試驗過程中排水總量的計算結果與實測結果的對比,兩者也相當接近。由于試驗的排水量與真空度、有效電壓、土樣和塑料排水板的滲透系數有關,是電滲法或電滲與其他工法聯合的綜合效果,圖9 的計算結果表明了本文建立的電滲固結理論方程以及所采用的計算條件和參數的合理性。

圖10 為計算得到的試驗結束時刻對角線A-A剖面兩極間土體中孔隙水壓力、沉降、最大和最小有效主應力分布。可見,在電滲和真空排水聯合作用下,土體中的最大負孔壓分布在兩極附近,而最小負孔壓分布在兩極中間略偏向陰極處。由于水流與土體介質的耦合關系,孔壓減小,土體的有效應力增大,則沉降也增大。因此,孔壓分布決定了土體中沉降和應力的分布,最大沉降發生在最大負孔壓處,在陽極和陰極表面,表面沉降呈兩極大、中間小的分布形式;土體中的應力分布同樣呈兩極大、中間小形式,說明了電滲和真空抽水均能使土體中產生負孔壓,從而增大土體的有效應力,達到土體的加固效果。

圖10 試驗結束時刻A-A 剖面孔壓、位移、 應力分布等值線 Fig.10 Isoline distribution of pore pressure, displacement, and stress in section A-A at the end of test

圖11 試驗300 h 后電勢沿兩極間分布情況 Fig.11 Potential distribution between anode and cathode after 300 hours test

圖11 為試驗結束時刻兩極間土體中電勢分布的計算結果。可見,兩極間的電勢并非呈直線分布,而是近似由3 條不同斜率的直線構成。靠近電極附 近的土體中電勢梯度較大,而兩極中間土體中的電勢梯度較小。該分布規律與界面電阻有關,若電極與土體的接觸面積增大,則界面電阻減小,電極附近土體中的電勢梯度也將減小。該現象與文獻[18]中的試驗結果一致。

5 結 語

電滲固結問題是一個典型的非飽和土固結問題,它的解決必須充分考慮電場、滲流場以及應力場之間的相互作用關系。本文在非飽和土多孔介質力學理論的基礎上,利用電荷守恒原理、質量守恒原理、應力平衡方程、達西定律以及歐姆定律,推導出基于三場耦合的電滲固結理論方程,提出了相應的數值解法,且自主開發了有限元程序,并對電滲聯合真空排水的室內模型試驗進行了數值模擬。算例表明,該模型的數值預測值與試驗實測值無論在量值還是趨勢上都吻合地較好,說明了本文所建立的電滲固結數值模型及其計算方法是正確、有效的,它不僅適用于工程中單獨的電滲固結問題,還可以解決電滲與其他工法聯合使用下的土體固結問題,從而為工程中涉及電滲固結的數值模擬提供了理論支持。由于電滲是一個涉及多個物理場相互作用的復雜問題,本文的工作尚處于初級階段,還有待于深入研究。

[1] 汪聞韶. 汪聞韶院士土工問題論文選集[M]. 北京: 中國水利水電出版社, 1999.

[2] CASAGRANDE L. Electric-osmosis on soil[J]. Geotechnique, 1949, (1): 159-177.

[3] CASAGRANDE L. Electro-osmotic stabilization of soils[J]. Journal of the Boston Society of Civil Engineers, 1952, (39): 51-83.

[4] CASAGRANDE L. Stabilization of soils by means of electroosmosis-state-of-the-art[J]. Journal of the Boston Society of Civil Engineers, 1983, 69(2): 255-302.

[5] MILLIGAN V. First application of electro-osmosis to improve friction pile capacity-three decades later[J]. Proceedings of the Institution of Civil Engineers, Geotechnical Engineering, 1995, 113:112-117.

[6] PERRY W. Electro-osmosis dewaters large foundation excavation[J]. Construction Methods and Equipment, 1963, 45(9): 116-119.

[7] ESRIG M I. Pore pressures, consolidation and electrokinetics[J]. Journal of the Soil Mechanics and Foundations Division, ASCE, 1968, 94: 899-921.

[8] WAN T Y, MITCHELL J K. Electro-osmotic consolidation of soil[J]. Journal of the Geotechnical Engineering Division, 1976, 102(GTS): 473-491.

[9] LEWIS W R, HUMPHESON C. Numerical analysis of electroosmotic flow in soils[J]. Journal of the SMFD, ASCE, 1973, 95 (SM4): 603-616.

[10] 莊艷峰, 王釗, 林清. 電滲的能級梯度理論[J]. 哈爾濱工業大學學報, 2005, 37(2): 283-286. ZHUANG Yan-feng, WANG Zhao, LIN Qing. Energy level gradient theory for electro-osmotic consolidation[J]. Journal of Harbin Institute of Technology, 2005, 37(2): 283-286.

[11] 莊艷峰, 王釗. 電滲的電荷累積理論[J]. 巖土力學, 2005, 26(4): 629-632. ZHUANG Yan-feng, WANG Zhao. Electric charge accumulation theory for electric- osmotic consolidation[J]. Rock and Soil Mechanics, 2005, 26(4): 629-632.

[12] 蘇金強, 王釗. 電滲的二維固結理論[J]. 巖土力學, 2004, 25(1): 125-131. SU Jin-qiang, WANG Zhao. Theory of two-dimensional electric-osmotic consolidation of soils[J]. Rock and Soil Mechanics, 2004, 25(1): 125-131.

[13] 李瑛, 龔曉南, 盧萌盟, 等. 堆載-電滲聯合作用下的耦合固結理論[J]. 巖土工程學報, 2010, 32(1): 77-81. LI Ying, GONG Xiao-nan, LU Meng-meng, et al. Coupling consolidation theory under combined action of load and electro-osmosis[J]. Chinese Journal of Geotechnical Engineering, 2010, 32(1): 77-81.

[14] 劉國華, 王振宇, 黃建平. 土的電阻率特性及其工程應用研究[J]. 巖土工程學報, 2004, 26(1): 83-87. LIU Guo-hua, WANG Zhen-yu, HUANG Jian-ping. Research on electrical resistivity feature of soil and it’ s application[J]. Chinese Journal of Geotechnical Engineering, 2004, 26(1): 83-87.

[15] VAN GENUCHTEN T H M. A closed form equation predicting the hydraulic conductivity of unsaturated soils[J]. Soil Science Society of America Journal, 1980, 44:892-898.

[16] MUALEM Y. A new model for predicting the hydraulic conductivity of unsaturated porous media[J]. Water Resources Research, 1976, 12(3): 513-522.

[17] 汪俊波, 劉斯宏, 徐偉, 等. 電滲法處理大連大窖灣超軟土室內試驗[J]. 水運工程, 2010, 437(1): 15-19.

[18] 莊艷峰, 王釗. 電滲固結中的界面電阻問題[J]. 巖土力學, 2004, 25(1): 117-120. ZHUANG Yan-feng, WANG Zhao. Study on interface electric resistance of electroosmotic consolidation[J]. Rock and Soil Mechanics, 2004, 25(1): 117-120.

猜你喜歡
理論
堅持理論創新
當代陜西(2022年5期)2022-04-19 12:10:18
神秘的混沌理論
理論創新 引領百年
相關于撓理論的Baer模
多項式理論在矩陣求逆中的應用
基于Popov超穩定理論的PMSM轉速辨識
大電機技術(2017年3期)2017-06-05 09:36:02
十八大以來黨關于反腐倡廉的理論創新
“3T”理論與“3S”理論的比較研究
理論宣講如何答疑解惑
學習月刊(2015年21期)2015-07-11 01:51:44
婦女解放——從理論到實踐
主站蜘蛛池模板: 美女视频黄又黄又免费高清| 国产第一色| 国产性爱网站| 五月天久久婷婷| 欧美亚洲国产日韩电影在线| 制服丝袜 91视频| 九九精品在线观看| 青青草国产免费国产| 国产精品漂亮美女在线观看| 中文字幕欧美日韩| 国产精品3p视频| 1769国产精品免费视频| 国产精品网拍在线| 久久久久青草大香线综合精品| 直接黄91麻豆网站| 久久a毛片| 中文字幕色在线| 69视频国产| 久久综合亚洲鲁鲁九月天| 亚洲aⅴ天堂| 亚洲天堂精品在线观看| 国产一区免费在线观看| 久久无码高潮喷水| 久久久久亚洲精品成人网| 一本大道东京热无码av| 亚洲娇小与黑人巨大交| 国内老司机精品视频在线播出| 婷婷综合缴情亚洲五月伊| 91精品人妻一区二区| 亚洲精品第一页不卡| 国产精品美人久久久久久AV| 久久www视频| 99久久人妻精品免费二区| 亚洲欧美日韩另类在线一| 亚洲区视频在线观看| 99re经典视频在线| 91亚洲精品国产自在现线| 久久久久国产精品熟女影院| 香蕉视频在线观看www| 亚洲天堂久久| 国产特级毛片| 18禁影院亚洲专区| 国产成人无码综合亚洲日韩不卡| 99热这里只有精品国产99| 精品夜恋影院亚洲欧洲| 国产婬乱a一级毛片多女| 午夜啪啪福利| 97色婷婷成人综合在线观看| 九月婷婷亚洲综合在线| 极品国产一区二区三区| 性色一区| 婷婷色中文网| 午夜三级在线| 波多野结衣亚洲一区| 欧美午夜网| 久久综合亚洲鲁鲁九月天| 国产欧美视频综合二区| 亚洲精品动漫| 天天综合网亚洲网站| 久久狠狠色噜噜狠狠狠狠97视色| 午夜国产在线观看| 亚洲人成人无码www| 91毛片网| 思思热在线视频精品| 国产丝袜91| 国产午夜一级毛片| 国产麻豆另类AV| 天天操精品| 亚洲中文在线看视频一区| 日韩高清成人| 亚洲国产成人精品青青草原| 亚洲欧美精品在线| 秘书高跟黑色丝袜国产91在线| 成人免费视频一区二区三区 | 国产v精品成人免费视频71pao| 91小视频版在线观看www| 97亚洲色综久久精品| 日韩AV手机在线观看蜜芽| 亚洲国产日韩在线成人蜜芽| 国产对白刺激真实精品91| 久久中文无码精品| 日本欧美午夜|