亚洲免费av电影一区二区三区,日韩爱爱视频,51精品视频一区二区三区,91视频爱爱,日韩欧美在线播放视频,中文字幕少妇AV,亚洲电影中文字幕,久久久久亚洲av成人网址,久久综合视频网站,国产在线不卡免费播放

        ?

        超固結(jié)土模型平面應(yīng)變分叉理論解與數(shù)值模擬

        2011-01-31 06:10:24孫德安甄文戰(zhàn)
        關(guān)鍵詞:本構(gòu)主應(yīng)力傾角

        孫德安, 段 博, 甄文戰(zhàn)

        (上海大學(xué)土木工程系,上海200072)

        土體應(yīng)變局部化最常見的現(xiàn)象就是巖土工程中宏觀剪切帶的出現(xiàn).實(shí)際工程中,基坑坍塌、山體滑坡、地基失穩(wěn)破壞等具有共同的破壞機(jī)理,即土體在特定的外荷載作用下產(chǎn)生局部變形,伴隨著應(yīng)變局部化的發(fā)展逐漸形成剪切帶,最后導(dǎo)致失穩(wěn)破壞.因此,研究土體應(yīng)變局部化的形成和發(fā)展對(duì)評(píng)價(jià)和預(yù)防各類工程事故具有重要意義.

        巖土材料應(yīng)變局部化觸發(fā)剪切帶的現(xiàn)象早已在室內(nèi)試驗(yàn)中觀察到.例如,文獻(xiàn)[1-3]通過真三軸試驗(yàn),分析了應(yīng)力狀態(tài)對(duì)砂土變形局部化的影響,研究了砂在不同應(yīng)力路徑下的分叉特性;Khalid等[4]對(duì)黏性土局部化分叉進(jìn)行了平面應(yīng)變及軸對(duì)稱條件下的試驗(yàn)研究.另外,從理論上對(duì)材料局部化分叉加以分析也越來越受到重視.Rudnicki等[5]和Vardoulakis等[6]針對(duì)不同性質(zhì)巖土材料提出了應(yīng)變局部化準(zhǔn)則;Anand等[7-8]應(yīng)用分叉理論研究了平面應(yīng)變條件下彈塑性材料的初始剪切帶行為;錢建固等[9]理論推導(dǎo)了平面應(yīng)變條件下Mohr-Coulomb彈塑性模型局部化分叉的理論解.隨著計(jì)算機(jī)計(jì)算能力的迅速發(fā)展,數(shù)值模擬為巖土材料應(yīng)變局部化研究提供了一個(gè)有利途徑,其中Huang等[10]采用亞塑性本構(gòu)模型研究了砂的剪切帶問題;徐連民等[11-12]采用高精度歐拉向后應(yīng)力積分算法,模擬分析了平面應(yīng)變時(shí)正常固結(jié)土和超固結(jié)土剪切帶的形成過程.

        本工作主要針對(duì)基于伏斯列夫(Hvorslev)面超固結(jié)黏土三維彈塑性本構(gòu)模型[13-14],推導(dǎo)平面應(yīng)變條件下局部化分叉理論解,并據(jù)此理論解給出常平均主應(yīng)力和常最小主應(yīng)力平面應(yīng)變條件下黏土局部化分叉的理論預(yù)測(cè);然后,采用與文獻(xiàn)[15]相似的方法,編寫材料子程序,把基于伏斯列夫面的超固結(jié)黏土三維本構(gòu)模型嵌入到大型非線性有限元軟件ABAQUS中,并在平面應(yīng)變路徑下對(duì)局部化分叉問題進(jìn)行數(shù)值分析、計(jì)算,給出分叉點(diǎn)的數(shù)值解;最后,對(duì)比分析分叉點(diǎn)的數(shù)值解和理論解,驗(yàn)證理論推導(dǎo)的可靠性.

        1 超固結(jié)黏土本構(gòu)模型

        本構(gòu)模型的彈塑性矩陣在很大程度上決定了應(yīng)變局部化的形成與發(fā)展.由于復(fù)雜的生成條件和地質(zhì)歷史作用,超固結(jié)土體受到的前期固結(jié)壓力大于當(dāng)前的固結(jié)壓力時(shí),其力學(xué)特性不同于正常固結(jié)土,主要表現(xiàn)為低孔隙率、強(qiáng)度高、壓縮性小、具有剪脹和應(yīng)變軟化等變形特性.因此,合理選擇既能反映應(yīng)變局部化的變形特點(diǎn),又能反映土體在數(shù)值模擬中的漸進(jìn)性破壞的超固結(jié)土本構(gòu)模型,顯得尤為重要.

        本工作采用姚仰平等[13-14]提出的基于伏斯列夫面的超固結(jié)黏土三維彈塑性本構(gòu)模型,該模型可以反映土體硬化、軟化、剪縮及剪脹等變形特性.通過變換應(yīng)力的方法[16]將SMP(spatially mobilized plane)準(zhǔn)則與該模型相結(jié)合,實(shí)現(xiàn)該模型的三維化,從而可以合理地反映土體在三維應(yīng)力下的變形和強(qiáng)度特性,并且三維化后模型的塑性流動(dòng)方向與屈服面不正交,為非相關(guān)聯(lián)本構(gòu)模型.模型的屈服函數(shù)表達(dá)式為

        式中,

        H為硬化參數(shù),λ為壓縮指數(shù),κ為回彈指數(shù),e0為初始孔隙比,為塑性體積應(yīng)變,M為臨界狀態(tài)時(shí)的應(yīng)力比,Mf為潛在強(qiáng)度,為變換應(yīng)力空間中對(duì)應(yīng)的應(yīng)力比,及分別為變換應(yīng)力空間中對(duì)應(yīng)的初始平均應(yīng)力、平均應(yīng)力及廣義剪應(yīng)力,為變換應(yīng)力空間中對(duì)應(yīng)的前期固結(jié)壓力,R為超固結(jié)參數(shù).

        2 平面應(yīng)變局部化分叉解析

        應(yīng)變局部化現(xiàn)象是材料失穩(wěn)破壞的先兆,它的發(fā)生點(diǎn)對(duì)應(yīng)于材料的局部化分叉點(diǎn),而分叉點(diǎn)的預(yù)測(cè)又強(qiáng)烈依賴于本構(gòu)模型的彈塑性矩陣.基于伏斯列夫面的超固結(jié)黏土三維本構(gòu)模型能夠合理描述黏土材料應(yīng)力應(yīng)變關(guān)系、硬化、軟化、剪縮及剪脹等變形特性,并且能夠合理預(yù)測(cè)應(yīng)變局部化的發(fā)生.以下就該模型采用Rudnicki等[5]提出的聲學(xué)張量法對(duì)平面應(yīng)變條件下的局部化分叉進(jìn)行解析(見圖1).

        圖1 剪切帶示意圖Fig.1 Chart of shear bands

        聲學(xué)張量法認(rèn)為:均質(zhì)各向同性的材料,在應(yīng)變局部化平面P上的速度場(chǎng)保持連續(xù),而速度梯度場(chǎng)產(chǎn)生跳躍,且速度及速度梯度方向在P面內(nèi)保持均勻.因此,有

        可知

        因速度梯度場(chǎng)產(chǎn)生跳躍,則有

        式中,εi,j和ε·i,j分別為應(yīng)變張量和應(yīng)變率張量,ui為位移張量.

        剪切帶邊界上還滿足靜力平衡邊界條件,即

        把式(9)代入式(10),得

        若應(yīng)變局部化發(fā)生,則分叉條件表示為

        對(duì)于平面應(yīng)變問題,有

        則分叉的表達(dá)式為

        式中,

        若材料發(fā)生平面應(yīng)變局部化分叉,則式(14)有實(shí)數(shù)解.由根的判別式,可知

        剪切帶法線傾角為

        本工作采用的基于伏斯列夫面的超固結(jié)黏土三維本構(gòu)模型的彈塑性矩陣為

        式中,

        其中p為平均主應(yīng)力,δij為Kronecker張量,ν為泊松比,I1,I2及I3為第一、第二和第三應(yīng)力不變量,skl為偏應(yīng)力張量.

        用式(17)中的彈塑性矩陣Dijkl替換式(12)中的,得到的式(15)和(16)即為基于伏斯列夫面的超固結(jié)黏土三維本構(gòu)模型開始觸發(fā)平面應(yīng)變分叉的必要條件和相應(yīng)剪切帶法線傾角表達(dá)式.

        基于本工作中推導(dǎo)的關(guān)于基于伏斯列夫面的超固結(jié)黏土三維本構(gòu)模型在平面應(yīng)變條件下的局部化分叉理論解,下面給出藤森黏土[13-14]局部化分叉的理論預(yù)測(cè).藤森黏土的材料特性參數(shù)如表1所示.

        表1 藤森黏土的模型參數(shù)Table 1 Model parameters for Fujinomori clay

        根據(jù)表1的參數(shù)值,可以預(yù)測(cè)常最小主應(yīng)力條件下不同超固結(jié)比(over-consolidation ratio,OCR)時(shí)的應(yīng)力應(yīng)變關(guān)系和相應(yīng)的應(yīng)力比峰值時(shí)的應(yīng)變.圖2給出了最小主應(yīng)力為200 kPa時(shí)平面應(yīng)變條件下不同超固結(jié)比藤森黏土的理論解分叉點(diǎn)和峰值點(diǎn),其中縱坐標(biāo)對(duì)應(yīng)最大主應(yīng)變.可以看出,不同超固結(jié)比土體的局部化分叉均發(fā)生在應(yīng)力應(yīng)變塑性硬化階段,即各分叉點(diǎn)皆在對(duì)應(yīng)的峰值點(diǎn)之前;并且超固結(jié)比越大,分叉點(diǎn)越接近峰值點(diǎn).

        圖3從理論上分析了最小主應(yīng)力為200 kPa時(shí),平面應(yīng)變條件下不同超固結(jié)比對(duì)剪切帶傾角α的影響,α=90°-θ,其中θ為剪切帶法線傾角(見圖1).由圖3可知,當(dāng)OCR為2時(shí),剪切帶傾角α=47.70°,當(dāng)OCR為12時(shí),θ=51.51°;隨著OCR值增大,局部化分叉觸發(fā)的剪切帶傾角呈現(xiàn)緩慢增大的趨勢(shì).

        圖2 平面應(yīng)變下不同超固結(jié)比的藤森黏土的分叉理論解Fig.2 Analytical solutions for bifurcation of Fujinomori clay with different OCRs under plane strain condition

        圖3 不同超固結(jié)比的藤森黏土的剪切帶傾角Fig.3 Inclination angles of shear band on Fujinomori clay with different OCRs

        圖4為理論預(yù)測(cè)的平均主應(yīng)力為200 kPa時(shí),平面應(yīng)變條件下剪切帶傾角α隨最大主應(yīng)變?chǔ)?變化的過程.當(dāng)ε1=0~0.9%時(shí),局部化不明顯,沒能形成剪切帶;在ε1=0.9%~2.5%時(shí),開始形成明顯的應(yīng)變局部化平面,且其傾角從產(chǎn)生就開始隨著最大主應(yīng)變的增加而顯著減小;局部化繼續(xù)發(fā)展,應(yīng)變局部化平面傾角趨于穩(wěn)定,直至ε1=2.9%時(shí),發(fā)生分叉,應(yīng)變局部化平面轉(zhuǎn)變?yōu)槊黠@的剪切帶.可以看出,分叉發(fā)生后,剪切帶傾角α減小很緩慢,趨于穩(wěn)定.

        3 平面應(yīng)變分叉數(shù)值預(yù)測(cè)

        3.1 數(shù)值解析條件

        為驗(yàn)證本工作推導(dǎo)的局部化分叉理論解的可靠性,采用與文獻(xiàn)[15]相似的方法編寫材料子程序,把基于伏斯列夫面的超固結(jié)黏土三維本構(gòu)模型嵌入到大型非線性有限元軟件ABAQUS中,并對(duì)超固結(jié)黏土在平面應(yīng)變條件下的局部化分叉問題進(jìn)行數(shù)值預(yù)測(cè)分析.在計(jì)算過程中,因整體剛度矩陣出現(xiàn)奇異,即負(fù)特征值,導(dǎo)致計(jì)算停止.負(fù)特征值表示在應(yīng)力應(yīng)變關(guān)系曲線的硬化階段出現(xiàn)負(fù)斜率,與原應(yīng)力應(yīng)變關(guān)系不同,控制偏微分方程由橢圓型變?yōu)殡p曲線型,即分叉.從數(shù)學(xué)上分析,負(fù)特征值意味著剛度矩陣不正定,方程組存在多組解,從而可求得負(fù)特征值出現(xiàn)時(shí)對(duì)應(yīng)的應(yīng)力所在位置,即為數(shù)值計(jì)算的分叉點(diǎn).

        圖4 最大主應(yīng)變與剪切帶傾角的關(guān)系Fig.4 Relationship between maximum principal strain and inclination angle of shear band

        材料的失穩(wěn)分叉分為彌散連續(xù)分叉和不連續(xù)失穩(wěn)分叉,前者指失穩(wěn)不是以出現(xiàn)明顯的不連續(xù)帶為標(biāo)志,而是由大量彌散分布的微破壞構(gòu)成導(dǎo)致的.研究認(rèn)為,這種分叉對(duì)應(yīng)于材料控制方程橢圓性的喪失,即切線剛度矩陣出現(xiàn)奇異,det(Dep)=0.本工作中數(shù)值分析時(shí)涉及的分叉就屬于彌散連續(xù)分叉.不連續(xù)分叉指宏觀上發(fā)生分叉,即出現(xiàn)不連續(xù)剪切帶,在剪切帶上應(yīng)變場(chǎng)產(chǎn)生跳躍,即det(nDepn)=0.本工作在理論上給出了不連續(xù)分叉的解.不連續(xù)分叉在彌散連續(xù)分叉的基礎(chǔ)上考慮了剪切帶方向的影響,二者結(jié)果相差不會(huì)很大.如果本構(gòu)模型為關(guān)聯(lián),則彌散連續(xù)分叉解對(duì)應(yīng)的硬化模量大于或等于不連續(xù)分叉解對(duì)應(yīng)的硬化模量,則不連續(xù)分叉先于彌散連續(xù)分叉出現(xiàn);而對(duì)非關(guān)聯(lián)模型,彌散連續(xù)分叉解對(duì)應(yīng)的硬化模量可能大于、小于或等于不連續(xù)分叉解對(duì)應(yīng)的硬化模量,則彌散連續(xù)分叉有可能先于或落后于不連續(xù)分叉的出現(xiàn).

        因?yàn)閼?yīng)變局部化是材料隨著均勻加載由均勻變形場(chǎng)過渡到非均勻變形場(chǎng)的一個(gè)現(xiàn)象,所以初始均勻各向同性材料就成為本工作數(shù)值研究的對(duì)象.這里,分析對(duì)象為邊長(zhǎng)10 cm的均勻各向同性立方體,采用8節(jié)點(diǎn)六面體線性完全積分單元類型,將試樣劃分為2 744個(gè)單元,有限元網(wǎng)格劃分如圖5所示.坐標(biāo)軸x方向設(shè)為最大主應(yīng)力方向,在x方向上,設(shè)定試樣上下兩端面光滑,給上端面施加豎向強(qiáng)制位移;剪切過程中固定住y方向上兩端面在y方向上的位移,同時(shí)固定下端面在x方向;z方向?yàn)閼?yīng)力邊界.試樣假定為藤森黏土,本構(gòu)模型參數(shù)如表1所示.

        圖5 有限元網(wǎng)格Fig.5 Meshes for finite element analysis

        3.2 結(jié)果分析

        圖6為數(shù)值模擬平面應(yīng)變條件下的2條不同的應(yīng)力路徑,分別為常最小主應(yīng)力(σ3=200 kPa)和常平均主應(yīng)力(p=200 kPa).

        圖6 應(yīng)力路徑Fig.6 Stress path

        圖7為2條不同應(yīng)力路徑的應(yīng)力應(yīng)變關(guān)系,σ1/σ3為最大主應(yīng)力與最小主應(yīng)力之比.圖中對(duì)比了理論解分叉點(diǎn)與數(shù)值解分叉點(diǎn).當(dāng)σ3=200 kPa時(shí),理論解分叉點(diǎn)對(duì)應(yīng)最大主應(yīng)變?yōu)?.4%,峰值點(diǎn)對(duì)應(yīng)的最大主應(yīng)變?yōu)?.1%,數(shù)值解分叉點(diǎn)對(duì)應(yīng)的最大主應(yīng)變?yōu)?.9%;當(dāng)p=200 kPa時(shí),峰值點(diǎn)對(duì)應(yīng)的最大主應(yīng)變?yōu)?.6%,理論解分叉點(diǎn)對(duì)應(yīng)的最大主應(yīng)變?yōu)?.9%,數(shù)值解分叉點(diǎn)對(duì)應(yīng)的最大主應(yīng)變?yōu)?.7%.上述結(jié)果表明,基于伏斯列夫面的超固結(jié)黏土三維彈塑性本構(gòu)模型,在常最小主應(yīng)力和常平均主應(yīng)力的剪切應(yīng)力路徑下,應(yīng)變局部化分叉均發(fā)生在土體的應(yīng)變硬化階段.

        圖7 主應(yīng)變和應(yīng)力比的關(guān)系Fig.7 Relationship between principal strain and stress ratio

        圖8給出了在σ3=200 kPa時(shí),平面應(yīng)變條件下不同超固結(jié)比的藤森黏土應(yīng)變局部化分叉理論解與數(shù)值解的比較.可以發(fā)現(xiàn),各超固結(jié)比下土體的數(shù)值解分叉點(diǎn)均發(fā)生在峰值前,進(jìn)一步證實(shí)了平面應(yīng)變條件下土體應(yīng)變局部化發(fā)生在應(yīng)力應(yīng)變塑性硬化階段的理論解結(jié)果(見圖2);并且隨著超固結(jié)比的增大,土體強(qiáng)度越高,數(shù)值解分叉點(diǎn)越接近峰值點(diǎn),這與理論解分叉點(diǎn)隨固結(jié)比的增長(zhǎng)而接近峰值點(diǎn)的特點(diǎn)吻合得很好.不同應(yīng)力路徑與超固結(jié)比條件下分叉的理論解與數(shù)值解的應(yīng)力比誤差如表2所示,不難看出,二者應(yīng)力比相差較小,這說明了平面應(yīng)變局部化分叉理論解的合理性.

        圖8 不同超固結(jié)比的藤森黏土的應(yīng)力比應(yīng)變關(guān)系Fig.8 Stress ratio-strain relationships of Fujinomori clay with different OCRs

        表2 局部化分叉理論解和數(shù)值解應(yīng)力比的比較Table 2 Comparison of stress ratio between analytical and numerical solutions for localization bifurcation

        4 結(jié)束語

        本工作基于聲學(xué)張量法系統(tǒng)地對(duì)基于伏斯列夫面的超固結(jié)黏土三維彈塑性本構(gòu)模型在平面應(yīng)變條件下局部化分叉進(jìn)行了解析,并運(yùn)用理論解對(duì)常最小主應(yīng)力平面應(yīng)變條件下不同超固結(jié)比土體的分叉點(diǎn)和剪切帶傾角進(jìn)行了對(duì)比分析.理論分析表明,超固結(jié)比越大,分叉點(diǎn)發(fā)生得越早,對(duì)應(yīng)的最大主應(yīng)變?cè)叫?,剪切帶傾角越大;同時(shí)理論分析了常平均主應(yīng)力平面應(yīng)變條件下局部化分叉形成的剪切帶的傾角隨最大主應(yīng)變變化的關(guān)系,理論證明了剪切帶形成后,剪切帶傾角隨最大主應(yīng)變?cè)龃笞兓⑿?,趨于穩(wěn)定;最后,本工作借助非線性有限元軟件ABAQUS分別在常最小主應(yīng)力和常平均主應(yīng)力平面應(yīng)變條件下對(duì)局部化分叉進(jìn)行了數(shù)值模擬分析,得出了分叉點(diǎn)數(shù)值解與理論解基本一致的結(jié)論,從而驗(yàn)證了平面應(yīng)變條件下局部化分叉理論解的可靠性.

        [1] CHUJ,LOS C R,LEEI K.Strain softening and shear band formation ofsand in multi-axialtest[J].Geotechnique,1996,46(1):63-82.

        [2] WANGQ,LADEP V.Shear banding in true triaxial tests and its effect on failure in sand[J].Journal of Engineering Mechanics,2001,127(8):754-761.

        [3] SUND A,HUANGW X,YAOY P.An experimental study of failure and softening in sand under threedimensional stress condition[J].Granular Matter,2008,10(3):187-195.

        [4] KHALIDA,ALSHIBLII S A.Strain localization in clay:plane strain versus triaxial loading conditions[J].Geotechnical and Geological Engineering,2007,25:45-55.

        [5] RUDNICKIJ W,RICEJ R.Conditions for the localization of the deformation in pressure sensitive dilatant materials[J].Journal of the Mechanics and Physics of Solids,1975,23:371-394.

        [6] VARDOULAKISI,SULEMJ.Bifurcation analysis in geomechanics[M].New York:Blackie Academic&Professional,1995.

        [7] ANANDL,SPITZIGW A.Initiation of localized shear bands in plane-strain[J].Journal of the Mechanics and Physics of Solids,1980,28(2):113-128.

        [8] ANANDL,SPITZIGW A.Shear-band orientations in plane-strain[J].Acta Metallurgica,1982,30(2):553-561.

        [9] 錢建固,黃茂松.土體應(yīng)變局部化現(xiàn)象的理論解析[J].巖土力學(xué),2005,26(3):432-436.

        [10] HUANGW X,SUND A,SLOANS W.Analysis of the failure mode and softening behaviour of sands in true triaxial tests[J].International Journal of Solids and Structures,2007,44:1423-1437.

        [11] 徐連民,王興然.用有限變形理論研究黏性土試樣中變形的局部化問題[J].巖土工程學(xué)報(bào),2004,26(2):125-129.

        [12] 徐連民,朱合華,中井照夫,等.超固結(jié)粘土的剪切帶數(shù)值模擬[J].巖土力學(xué),2006,27(1):61-66.

        [13] 姚仰平,侯偉,周安楠.基于Hvorslev面的超固結(jié)黏土本構(gòu)模型[J].中國(guó)科學(xué):E輯(技術(shù)科學(xué)),2007,37 (11):1417-1429.

        [14] YAOY P,LIZ Q,HOUW,et al.Overconsolidated clay model based on revised Hvorslev envelope[C]∥Proceeding of the 3rd Sino-Japan Geotechnical Symposium. Beijing: China Communication Press,2007:698-705.

        [15] 孫德安,甄文戰(zhàn).不同應(yīng)力路徑下剪切帶的數(shù)值模擬[J].巖土力學(xué),2010,31(7):2253-2258.

        [16] MATSUOKAH,YAOY P,SUND A.The Cam-clay model revised by the SMP criterion[J].Soils and Foundations,1999,39(1):81-95.

        猜你喜歡
        本構(gòu)主應(yīng)力傾角
        地球軸傾角的改斜歸正
        激光傾角儀在CT引導(dǎo)下經(jīng)皮肺穿刺活檢中的應(yīng)用
        車輪外傾角和前束角匹配研究
        北京汽車(2021年3期)2021-07-17 08:25:42
        離心SC柱混凝土本構(gòu)模型比較研究
        復(fù)合斷層對(duì)地應(yīng)力的影響研究
        鋸齒形結(jié)構(gòu)面剪切流變及非線性本構(gòu)模型分析
        一種新型超固結(jié)土三維本構(gòu)模型
        深部沿空巷道圍巖主應(yīng)力差演化規(guī)律與控制
        考慮中主應(yīng)力后對(duì)隧道圍巖穩(wěn)定性的影響
        探討大傾角綜采面“三機(jī)”防倒防滑技術(shù)控制研究
        河南科技(2014年23期)2014-02-27 14:18:49
        国产18禁黄网站免费观看| 在线毛片一区二区不卡视频| 久久国产高潮流白浆免费观看| 粉嫩人妻91精品视色在线看| 国产激情艳情在线看视频| 久久久久久久久888| 国产三级精品美女三级| 日本不卡不二三区在线看| 日本一区二区三区小视频| 韩国日本一区二区在线| 日本少妇春药特殊按摩3| 躁躁躁日日躁| japanese色国产在线看视频| 日本一区二区高清精品| 少妇性俱乐部纵欲狂欢电影| 精品人体无码一区二区三区| 国产日韩午夜视频在线观看| 91精品久久久中文字幕| 亚洲中文字幕无码一久久区| 亚洲自偷自拍另类图片小说| 一区二区三区国产偷拍| 日本国产一区在线观看| 久久久久人妻一区二区三区| 婷婷成人基地| 亚洲综合久久一本久道| h视频在线播放观看视频| 亚洲色爱免费观看视频| 国产成人亚洲综合无码精品| 国产熟妇一区二区三区网站| 国产特级毛片aaaaaa高潮流水| 国产精品va在线观看无码| 麻豆国产成人AV网| 美女免费视频观看网址| 亚洲h在线播放在线观看h| 国产福利小视频在线观看| 亚洲国产91高清在线| 少妇高潮无套内谢麻豆传| 亚洲精品高清你懂的| 国产精品午夜福利天堂| 亚洲中文字幕无码av永久| 精品国产制服丝袜高跟|