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

凍結黏土非整數階修正S-M損傷蠕變模型研究

2025-02-19 00:00:00姚兆明賴龍輝湯海東陳軍浩
人民長江 2025年1期
關鍵詞:模型

摘要:

凍土蠕變特性是判斷凍土結構是否穩定的決定性因素。以山東張集地區黏土為研究對象,對凍結黏土進行單軸抗壓試驗與蠕變試驗,得到了在不同凍結溫度和應力加載等級下的凍結黏土蠕變規律。通過在Singh-Mitchell(S-M)模型中引入雙曲線函數和非整數階微積分,建立了能反映凍結溫度和應力加載等級因素影響的凍結黏土非整數階修正S-M蠕變模型。進一步引入損傷變量,并結合試驗數據建立了損傷變量與溫度關系公式,提出了非整數階修正S-M損傷蠕變模型,該模型能較好地描述凍結黏土非穩定蠕變階段。研究結果表明:隨著凍結溫度的降低,凍結黏土強度明顯提高,蠕變變形減小;應力加載等級變化對凍結黏土蠕變變形影響顯著。與經典S-M模型計算值和試驗值的對比表明,所建模型能夠準確反映溫度效應下不同應力加載等級的凍結黏土蠕變規律,并且擬合優度更高。該模型參數同時具有物理和數學意義,且模型參數較少,適用于實際工程應用。

關" 鍵" 詞:

凍土; 蠕變; 本構模型; Singh-Mitchell模型; 非整數階導數

中圖法分類號: TU445

文獻標志碼: A

DOI:10.16232/j.cnki.1001-4179.2025.01.024

收稿日期:2024-01-29;接受日期:2024-05-04

基金項目:

福建省自然科學基金面上項目(2022J01925);礦山地下工程教育部工程研究中心開放研究項目(JYBGCZX2021104)

作者簡介:

姚兆明,男,教授,博士,主要從事土體本構理論和巖土數值分析研究。E-mail:zhmyaoaust@126.com

Editorial Office of Yangtze River. This is an open access article under the CC BY-NC-ND 4.0 license.

文章編號:1001-4179(2025) 01-0180-07

引用本文:

姚兆明,賴龍輝,湯海東,等.凍結黏土非整數階修正S-M損傷蠕變模型研究

[J].人民長江,2025,56(1):180-186.

0" 引 言

人工凍結法是在軟土、流砂層中進行工程建設常用的一種施工方法。隨著城市地下交通大量建設,地下工程的施工深度不斷增加,大大增加了凍結法施工的風險。為了保障地下施工安全,研究凍結土體的蠕變特性十分關鍵。許多學者對凍土蠕變進行了深入的研究,建立了大量凍土蠕變模型,主要包括經驗模型、元件模型等[1-3]。

傳統經驗模型有對數型、雙曲線型、指數型及多項式型等。Yin等[4]使用對數函數建立模型,描述土體在一維應力條件下的彈塑性和蠕變行為;胡惠華等[5]指出雙曲線函數更適合描述一維固結蠕變試驗中的應力-應變曲線,以此代替Singh-Mitchell 模型中的指數函數,并建立了能夠更好反映應力水平影響的蠕變模型;余湘娟等[6]提出采用雙曲線形式對試驗數據進行擬合,建立了一維次固結的經驗模型,該模型考慮了壓力、超固結狀態對次壓縮的影響,能夠較好地計算次固結沉降量;黃海峰等[7]通過使用雙曲線函數描述應力-應變關系,并結合經驗模型理論,建立了一種新的經驗蠕變模型。經驗模型雖然缺少完善的物理概念,但具有擬合精度高、參數少、易于確定的優點,能夠方便地解決工程實際問題。

元件模型是由彈性體、塑性體、黏性體元件通過串并聯方式組合而成。咸冬冬等[8]通過大量分級加載和循環加卸載三軸蠕變試驗,提出了基于元件模型建立的黏彈-塑性變形的損傷蠕變模型。研究表明,元件模型的擬合精度與元件模型個數呈正比,但過多的元件個數會導致模型復雜,計算難度高,不利于工程實際使用。

以上模型主要用于常溫下巖土蠕變方面的研究,而對于凍土蠕變方面的模型研究較為少。凍結黏土的蠕變特性處于理想流體和理想固體之間,與普通巖土蠕變特性具有較大差異,通常的蠕變模型無法較好地描述溫度效應下凍結黏土的蠕變特性。研究發現通過引入非整數階導數理論,可以解決整數階模型不能較好處理蠕變時間記憶性的問題[9-11],對蠕變非線性性質描述更具優勢。陳軍浩等[12]基于廣義開爾文模型引入非整數階理論,提出能夠描述凍土蠕變特性的模型;Hou等[13]將彈黏塑性理論和非整數階導數相結合,構建出一種新的凍土蠕變本構模型;姚兆明等[14]在傳統經驗模型的基礎上,通過優化模型算法,得到了參數具有物理意義的凍土蠕變模型;楊歲橋等[15]進行了不同溫度、荷載及含冰量條件下的凍土蠕變試驗,并建立能反映各影響因素的凍土蠕變模型。目前大多數凍土蠕變模型未考慮在高應力下凍土損傷對蠕變的影響,也無法描述凍結黏土加速蠕變階段的蠕變特性。

本文基于Singh-Mitchell(S-M)模型并引入雙曲線函數和非整數階微積分,建立了凍結黏土非整數階S-M蠕變模型。在所建模型的基礎上,引入損傷因子,建立在高應力水平下能夠反映溫度效應的非整數階修正S-M損傷蠕變模型。通過將所建模型計算值與不同凍結溫度下的試驗值、S-M模型計算值進行對比,驗證了所建模型的準確性與適用性。

1" 凍結黏土單軸抗壓、蠕變試驗

1.1" 試驗方案

試驗黏土土樣取自山東張集礦,取樣深度為228.21~253.27 m,顏色為棕色,結構致密,可塑性差。按照GB/T 50123—2019《土工試驗方法標準》測得原狀土的物理性質指標見表1,試驗方案見表2~3。

由現場取土后,密封包裝送往實驗室,按照要求將原狀土加工成50 mm(直徑)×100 mm(高度)的圓柱體試樣,如圖1(a)所示。試驗前將試樣用薄膜密封包裹,防止水分流失。在-10,-15,-20 ℃共3種溫度下養護24 h,每個溫度養護3個試樣。養護后的試樣如圖1(b)所示。

1.2" 試驗結果及分析

單軸抗壓強度試驗在安徽理工大學WDT-100型凍土試驗機進行。試驗采用應變控制加載方式,應變加載速率為1%/min。破壞的土樣如圖2所示。

從圖2可以看出,在荷載作用下,凍結黏土內部出現滑裂面。試樣破壞主要有兩種形態,一種是由幾條上下貫通大裂縫產生的破壞,另一種是由多條細碎裂縫產生的破壞,兩者都表現出脆性破壞的特點。

從表2可見,在試驗凍結溫度范圍內,單軸抗壓強度平均值隨著溫度的降低而增大。凍土內冰晶主要起到顆粒之間的黏結作用。隨著溫度的降低,凍土中未凍水的含量降低,冰晶含量增加,提高了土體顆粒之間的黏結力,進而使凍土單軸抗壓強度明顯提高。

蠕變試驗采用分級連續加載,分級加載參照煤炭行業標準MT/T 593—2011《人工凍土物理力學性能試驗》,以 σ=0.3σs、0.5σs、0.7σs 進行連續加載。

圖3為不同凍結溫度各級加載下的蠕變曲線。

從圖3可以看出,在低應力階段(0.3σs、0.5σs)凍結黏土蠕變過程主要由衰減蠕變和穩定蠕變兩個階段組成。

衰減蠕變主要表現在蠕變的初期階段,在應力的作用下,凍結黏土的應變在短時間內迅速增大。相同凍結溫度下,應力加載等級的增加,會明顯延長衰減蠕變階段的時間,而凍結溫度對其影響較小。

穩定蠕變階段主要表現在蠕變值增長緩慢,逐漸趨于穩定。相同凍結溫度下,應變值受應力加載等級影響較大,應力加載等級的增加使應變值增大。

當處于高應力加載水平(0.7σs)時,蠕變初期會經歷短暫的衰減蠕變。此過程中凍結黏土的內部會產生結構性損傷并緩慢地累積,進而導致承載力降低,進入加速蠕變階段,直到土樣完全破壞。

2" 非整數階修正S-M蠕變模型

2.1" 修正S-M 經驗蠕變模型

S-M模型是廣泛應用的模型之一[16],能夠很好地反映常溫下土體的蠕變特性。S-M 蠕變模型方程為

ε=A1eb1σ(t/t1)λ(1)

式中:σ=η·σs,為蠕變加載應力值;η為應力加載等級;t1為單位時間;A1,b1,λ為蠕變試驗材料常數。

由式(1)可知,當σ=0時,ε≠0。為此,對式(1)進行修正:

ε=A1(eb1σ-1)(t/t1)λ(2)

式(2)克服了應力為零時ε不為零的缺限。

凍結黏土蠕變特征不僅取決于加載的應力大小,還受到加載時間的影響,具有較強的時間依賴性,但式(2)無法較好地體現蠕變的時間依賴性,而雙曲線函數更能體現凍結黏土蠕變的時間依賴性。圖4是不同凍結溫度下應變-應力等時曲線。

由圖4可見,在相同溫度下,選取的時間點越大,應力應變曲線的彎曲程度也愈加明顯,呈現出一簇相似的雙曲線。凍結黏土的蠕變特性也受凍結溫度的影響,凍結溫度越低,等時曲線的整體曲度明顯變大,雙曲線性蠕變特征更加明顯。

同時,由公式(2)可知,當t→∞時,ε→∞,這不符合在低應力蠕變加載時ε趨于穩定的性質,而雙曲線函數能描述當t→∞時,ε趨于常數的性質。鑒于此,對公式(2)進行如下修正:

ε=A(ebσ-1)tt+B(3)

式中:A,B,b為待定參數。

2.2" 非整數階修正S-M蠕變模型建立

非整數階理論能夠很好地反映凍結黏土蠕變處于理想流體和理想固體之間的特性。Riemann-Liouville非整數階導數常常被引用于本構模型中,用來構建非整數階本構蠕變模型。它的優點是對材料本身的歷史具有記憶性,能夠很好地對非線性特性進行描述[17]。凍結黏土蠕變特征適合使用非整數階導數來表達[18-20]。

當取t=v0τ時,它的非整數階微分為[21]

dβtdτβ=v0Γ(2-β)τ1-β(4)

式中:t為時間;v0為加載階段恒加載應變速率;τ為分數階積分后的自變量;βgt;0,且n-1lt;β≤n(n為正整數);Γ為Gamma函數。

式(4)中:τ1-β=(t/v0)1-β,即有:

dβtdτβ=vβ0Γ(2-β)t1-β(5)

將式(3)進行非整數階求導,并將式(5)代入式(3)中,可得到:

ε=A(ebσ-1)v0Γ(2-β)t1-βv0Γ(2-β)t1-β+B(6)

式中:A為與溫度有關的參數;B,b,β為待定參數。

2.3" 模型參數β的確定

令式(6)中

Y=v0Γ(2-β)εt1-βX=v0Γ(2-β)t1-β(7)

式中:參數v0=1。

將式(7)代入式(6)可得:

Y=1A(ebσ-1)X+BA(ebσ-1)(8)

殷德順等[22]認為當材料為理想流體時,β=1;當材料為理想固體時,β=0。對于凍土來說,β的取值范圍應該在0~1之間。

圖5是通過蠕變試驗數據來計算不同β值,得到的X-Y關系曲線。從圖5中可以看出,隨著β值的變小,X-Y關系曲線弧度減小,更接近于直線,線性關系增強。當β=0.1時,線性關系最明顯,線性擬合誤差最小,因此對于參數β=0.1最合適。

2.4" 參數A和B的確定

凍結溫度和應力加載等級是影響凍結黏土蠕變的主要因素。為考慮溫度和應力加載等級的影響,通過試驗數據的分析,建立關于凍結溫度和應力加載等級的模型參數公式。

將不同溫度下的 0.3σs,0.5σs應力加載等級的試驗值代入式(8),可得β=0.1時,不同溫度下的X-Y關系曲線,如圖6所示。

從圖6可以得到X-Y關系曲線的線性擬合方程,根據其擬合方程斜率和截距可以求得不同溫度下,0.3σs,0.5σs應力加載等級的參數A,B,b,見表4。

從表4可以發現,參數A和b同時受到溫度T影響,并且呈線性的關系。通過參數A、b與溫度T的擬合可以得到二者的函數關系為

A=0.0139T+1.002(9)

b=0.0176T+0.987(10)

由表4可知,溫度T和應力加載等級η對參數B的影響較小,且各條件下的參數B大小相近,所以對于參數B可取平均值B=0.61。

將各參數及其函數代入式(6),由此可以得出與溫度T和應力加載等級η有關的非整數階修正S-M蠕變模型為

ε=(0.0139T+1.002)·(ebσ-1)·1.04t0.91.04t0.9+0.61(11)

式中:b=0.0176T+0.987。

2.5" 所建蠕變模型與S-M蠕變模型對比

依據上節相同的思路,通過試驗數據求得傳統S-M蠕變模型為

ε=(-0.002T+0.17)e(0.04T+1.58)σt0.15(12)

根據式(11)建立的非整數階修正S-M蠕變模型與式(12)建立的S-M模型分別計算不同凍結溫度的黏土蠕變值。將兩模型計算出的結果與試驗值進行對比,并對模型進行擬合優度分析,如圖7所示。

從圖7可以看出,低應力加載水平(0.3σs、0.5σs)下,相同溫度的試驗曲線與計算曲線整體走勢相似且整體吻合度較高。

表5為所建模型與S-M模型擬合優度的對比。由表5可以看出,所建模型擬合優度整體更高,表明非整數階修正S-M蠕變模型可以很好地反映溫度效應下黏土的蠕變規律,同時驗證了該模型對凍結黏土的適用性。

3" 非整數階修正S-M損傷蠕變模型及驗證

3.1" 非整數階修正S-M損傷蠕變模型建立

上文非整數階修正S-M蠕變模型未考慮凍結黏土內部的損傷,無法準確描述在高應力水平(0.7σs)下凍結黏土蠕變特性。鑒于此,對上述模型引入損傷變量進一步進行修正。

式(13)是常用的損傷變量:

D=1-∫α0P(x)dx=exp(t/α)n(13)

式中:D為損傷變量;P(x) 為微元體概率密度函數;α為損傷變量參數;t為時間;

n為材料常數,取n=1。損傷變量D取值范圍為0~1。當巖土內部結構處于沒有任何裂縫的狀態時,D=0;當巖土處于完全破壞狀態時,D=1。

在式(6)中引入損傷變量D,則有:

ε=A(ebσ-1)v0Γ(2-β)t1-βv0Γ(2-β)t1-β+B·etα(14)

通過對不同溫度下0.7σs加載等級的蠕變數據進行擬合,得到-10,-15,-20 ℃下參數α分別為153.62,118.31,87.53。

可以看出參數α與溫度線性相關,溫度越低,參數α越小。由此可以發現在高應力水平下,溫度越低,損傷變量1/α也越大。對參數α進行擬合,建立與溫度T有關的公式:

α=6.61T+218.96(15)

高應力水平下考慮損傷效應的非整數階修正S-M損傷蠕變模型為

ε=(0.0139T+1.002)·(ebσ-1)·1.04t0.91.04t0.9+0.61etα(16)

將不同溫度下的試驗參數代入模型可以得到模型計算曲線,并將其與試驗曲線對比,如圖8所示。

由圖8可見,在高應力加載水平(0.7σs)下,不同溫度的試驗曲線與計算曲線均吻合較為理想,能夠較好地反映溫度效應下凍結黏土的損傷對蠕變曲線的影響。

3.2" 非整數階修正S-M損傷蠕變模型驗證

為進一步驗證溫度效應下所建模型的準確性和有效性,使用不同凍結溫度(-12℃,-18℃)的黏土蠕變試驗數據進行驗證分析,得到計算值與試驗值對比曲線,如圖9所示。

由圖9可以看出,所建模型較S-M模型總體上的擬合優度更好,并且能夠反映在高應力下凍結黏土的損傷蠕變特性。模型對于不同凍結溫度下的凍結黏土蠕變曲線擬合效果較好,能夠很好地反映不同溫度和應力加載等級的凍結黏土蠕變規律,驗證了所建模型的準確性和有效性。

4" 結 論

以山東張集礦的深部黏土為研究對象,開展不同凍結溫度下單軸抗壓試驗和蠕變試驗。通過分析不同凍結溫度的黏土單軸蠕變曲線,運用非整數階微積分函數和雙曲線函數,建立了能反映不同凍結溫度和應力加載等級因素影響的凍結黏土非整數階修正S-M損傷蠕變模型。主要結論如下:

(1) 雙曲線函數能夠較好地描述凍結黏土的非線性特征。非整數階微積分函數能夠反映凍結黏土材料處于理想流體與理想固體之間的特性,二者相結合建立了一種形式新穎的凍結黏土蠕變模型。

(2) 通過對不同凍結溫度下蠕變試驗數據進行驗證分析,模型計算值和試驗值一致性高,所建模型平均擬合優度為0.95,表明模型擬合精度較高。與經典S-M模型進行對比驗證,經典S-M模型平均擬合優度為0.84,表明所建模型能很好地反映凍結黏土的非穩定蠕變過程。

(3) 所建模型參數不僅具有物理概念,而且具有一定數學意義,模型參數較少,計算方便,便于工程應用。

參考文獻:

[1]" 魏堯,楊更社,葉萬軍.凍結溫度對凍融黃土力學特性的影響規律研究[J].長江科學院院報,2018,35(8):61-66.

[2]" 羅飛,張元澤,朱占元,等.一種青藏高原凍結砂土蠕變本構模型[J].哈爾濱工業大學學報,2020,52(2):26-32.

[3]" 董曉宏,張愛軍,連江波,等.非飽和凍融黃土固結蠕變特性研究[J].人民長江,2010,41(3):88-91.

[4]" YIN J H,GRAHAM J.Equivalent times and one-dimensional elastic viscoplastic modelling of time-dependent stress-strain behaviour of clays[J].Canadian Geotechnical Journal,1994,31(1):42-52.

[5]" 胡惠華,賀建清,聶士誠.洞庭湖砂紋淤泥質土一維固結蠕變模型研究[J].巖土力學,2022,43(5):1269-1276.

[6]" 余湘娟,殷宗澤,高磊.軟土的一維次固結雙曲線流變模型研究[J].巖土力學,2015,36(2):320-324.

[7] "黃海峰,巨能攀,周新,等.紅層滑坡滑帶土經驗型蠕變模型研究[J].人民長江,2017,48(5):91-95.

[8]" 咸冬冬,楊圣奇,柏正林,等.閃長玢巖蠕變試驗及黏彈-塑性損傷模型研究[J].人民長江,2023,54(12):225-232,240.

[9]" 李德建,劉校麟,韓超.基于等效黏彈性的變階分數階巖石損傷蠕變模型[J].巖土力學,2020,41(12):3831-3839.

[10]曹建軍,胡斌,王澤祺,等.基于分數階積分的軟弱夾層蠕變損傷模型研究[J].巖土力學,2024,45(2):454-464,476.

[11]吳斐,劉建鋒,武志德,等.鹽巖的分數階非線性蠕變本構模型[J].巖土力學,2014,35(增2):162-167.

[12]陳軍浩,姚兆明,徐穎,等.人工凍土蠕變特性粒子群分數階導數模型[J].煤炭學報,2013,38(10):1763-1768.

[13]HOU F,LI Q,LIU E,et al.A fractional creep constitutive model for frozen soil in consideration of the strengthening and weakening effects[J].Advances in Materials Science and Engineering,2016:5740292.

[14]姚兆明,李南,郭夢圓.人工凍結黏土經驗模型參數確定及蠕變規律驗證[J].金屬礦山,2022(11):58-63.

[15]楊歲橋,王寧寧,張虎.高溫凍土的蠕變特性試驗及蠕變模型研究[J].冰川凍土,2020,42(3):834-842.

[16]SINGH A,MITCHELL J K.General stress-strain-time function for clay[J].Journal of the Clay Mechanics and Foundation Division,1968,94(SM1):21-46.

[17]程愛平,付子祥,鄧代強,等.高低應力水平下膠結充填體蠕變特征及分數階本構模型[J].中國礦業大學學報,2023,52(2):276-285,299.

[18]沈朝輝,吳禮舟.考慮溫度效應的花崗巖分數階蠕變模型研究[J].人民長江,2019,50(8):167-171.

[19]姚兆明,蹇膨遠,郭夢圓.凍結黏土分數階損傷蠕變模型參數確定及驗證[J].科學技術與工程,2023,23(31):13507-13514.

[20]慕煥東,鄧亞虹,赫佳,等.基于分數階導數理論的Q_3黃土流變本構模型研究[J].工程地質學報,2022,30(4):1077-1086.

[21]姚兆明,蹇膨遠,郭夢圓.凍結黏土分數階損傷蠕變模型參數確定及驗證[J].科學技術與工程,2023,23(31):13507-13514.

[22]殷德順,和成亮,陳文.巖土應變硬化指數理論及其分數階微積分理論基礎[J].巖土工程學報,2010,32(5):762-766.

(編輯:鄭 毅)

Analysis of Singh-Mitchell damage creep model for frozen clay with non-integer order modification

YAO Zhaoming1,2,LAI Longhui1,2,TANG Haidong1,2,CHEN Junhao3

(1.College of Civil Engineering,Anhui University of Science and Technology,Huainan 232001,China;

2.Engineering Research Center of Underground Mine Construction of Ministry of Education,Huainan 232001,China;

3.School of Civil Engineering,Fujian University of Technology,Fuzhou 350118,China)

Abstract:

The creep characteristics of frozen soil are the determinant factors for judging the stability of frozen soil structures.In this study,the clay in a mining area of Shandong Province was taken as the research object.Uniaxial compression tests and creep tests were conducted on the frozen clay,through the tests we obtained the creep law of frozen clay at different freezing temperatures and stress load levels.The Singh-Mitchell (S-M) model was modified with the introduction of a hyperbolic function and non-integer order calculus,creating a non-integer order modified S-M creep model that can reflect the influence of freezing temperature and stress load levels.Further a damage variable was introduced,we established a formula reflecting the relation between damage variable and the temperature,so we proposed a non-integer order modified S-M creep-damage model that can describe the non-stable creep stage of frozen clay more accurately.Research results showed that as the freezing temperature fell,the strength of the frozen clay significantly increased,and the creep deformation reduced.The change of stress load levels also dramatically affected the creep deformation of the frozen clay.Compared with the calculated values of the classical S-M model and the experimental values,the established model could accurately reflect the creep law of frozen clay under different stress load levels in the context of temperature effects,and demonstrated a higher degree of fit.The parameters of this model have both physical and mathematical meanings,and the relative less number of parameters makes it more suitable for practical engineering applications.

Key words:

frozen soil; creep; constitutive model; Singh-Mitchell model; non-integer order derivative

猜你喜歡
模型
一半模型
一種去中心化的域名服務本地化模型
適用于BDS-3 PPP的隨機模型
提煉模型 突破難點
函數模型及應用
p150Glued在帕金森病模型中的表達及分布
函數模型及應用
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 久久国产乱子| 国产高潮流白浆视频| 91精品人妻一区二区| 国产精品黑色丝袜的老师| 国产97视频在线| 在线精品亚洲国产| 四虎在线观看视频高清无码| 广东一级毛片| 青青草国产在线视频| 国产精品亚洲精品爽爽 | 最新国语自产精品视频在| 免费国产小视频在线观看| 国产理论一区| 国产精品人人做人人爽人人添| 国内毛片视频| 国产在线91在线电影| 日韩精品免费在线视频| 日本高清免费一本在线观看| 欧美成人二区| 亚洲福利片无码最新在线播放| 熟妇无码人妻| 在线视频精品一区| 成人午夜精品一级毛片| 米奇精品一区二区三区| 国产在线视频导航| 欧美亚洲另类在线观看| 伊人成人在线视频| 青青草原偷拍视频| 国产后式a一视频| 精品中文字幕一区在线| 无码日韩精品91超碰| 日本一区高清| 亚洲AV无码乱码在线观看代蜜桃 | 女人毛片a级大学毛片免费| 欧美国产另类| 欧美黄网在线| 久久国产精品麻豆系列| 伊人天堂网| 鲁鲁鲁爽爽爽在线视频观看| 曰韩人妻一区二区三区| 免费无码一区二区| 亚洲黄网在线| 国产成人精品高清不卡在线| 一级做a爰片久久毛片毛片| 美美女高清毛片视频免费观看| 日本免费一区视频| 午夜a级毛片| 午夜影院a级片| 中文成人在线视频| 久久精品人人做人人爽97| a亚洲视频| 99精品免费欧美成人小视频 | 在线人成精品免费视频| 国产白丝av| 国产综合色在线视频播放线视| 狠狠色婷婷丁香综合久久韩国| 夜夜操天天摸| 免费人成网站在线观看欧美| 亚洲精品中文字幕无乱码| 色亚洲成人| 久久精品人妻中文视频| 国产精品久久久久久影院| 久久综合亚洲色一区二区三区| 在线毛片网站| 国产精品男人的天堂| 无码专区第一页| 国产成人久久777777| 国产女同自拍视频| 欧美成人aⅴ| 久久国产精品波多野结衣| 免费A级毛片无码免费视频| 国产在线一二三区| P尤物久久99国产综合精品| 国产全黄a一级毛片| 国产乱肥老妇精品视频| 国产精品亚洲αv天堂无码| 91在线播放免费不卡无毒| 中文无码精品A∨在线观看不卡| 亚洲高清在线天堂精品| 欧美精品高清| 中文精品久久久久国产网址| 一本大道香蕉中文日本不卡高清二区 |