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

FLAC3D在滑坡穩(wěn)定性分析評(píng)價(jià)中的應(yīng)用

2016-02-28 02:37:22趙忠海
城市地質(zhì) 2016年4期
關(guān)鍵詞:分析

趙忠海

(北京市地質(zhì)礦產(chǎn)勘查開發(fā)局,北京 100195)

FLAC3D在滑坡穩(wěn)定性分析評(píng)價(jià)中的應(yīng)用

趙忠海

(北京市地質(zhì)礦產(chǎn)勘查開發(fā)局,北京 100195)

本文以高治村滑坡為例,論述了FLAC3D在滑坡穩(wěn)定性分析評(píng)價(jià)中的應(yīng)用。作者首先介紹了FLAC3D的基本原理和主要流程,然后采用摩爾-庫倫模型,對(duì)高治村滑坡在天然狀態(tài)和暴雨?duì)顟B(tài)下的變形和位移情況進(jìn)行了數(shù)值模擬,模擬結(jié)果直觀地顯示了該坡體在上述兩種狀態(tài)下的塑性變形及水平方向位移的分布特征,清楚地顯示了坡體內(nèi)潛在的滑動(dòng)面。文章根據(jù)數(shù)值模擬結(jié)果,對(duì)該滑坡的穩(wěn)定性進(jìn)行了分析和評(píng)價(jià),對(duì)該滑坡的變形和破壞形式進(jìn)行了研究和探討。通過對(duì)該坡體的塑性變形與位移情況分析認(rèn)為,該滑坡在天然狀態(tài)下基本穩(wěn)定;在連續(xù)高強(qiáng)度降雨?duì)顟B(tài)下欠穩(wěn)定,局部失穩(wěn)滑塌的可能性很大,這與現(xiàn)場(chǎng)調(diào)查分析的結(jié)果一致。這個(gè)結(jié)果表明,利用FLAC3D進(jìn)行滑坡穩(wěn)定性分析和評(píng)價(jià)簡單、可行,且更加直觀、方便,具有傳統(tǒng)極限平衡法所無法比擬的優(yōu)勢(shì)。

滑坡;穩(wěn)定性分析評(píng)價(jià);FLAC3D數(shù)值模擬

0 前言

滑坡是山區(qū)和丘陵區(qū)比較常見的一種地質(zhì)災(zāi)害,是斜坡上部的巖土體,受自然及人為等因素的影響,在重力作用下,沿一定的軟弱面或滑動(dòng)帶,整體地或分散地順坡向下滑動(dòng)的自然現(xiàn)象。滑坡的危害性極大,是僅次于地震和洪水的一種嚴(yán)重的地質(zhì)災(zāi)害。滑坡穩(wěn)定性分析是滑坡隱患治理工作的基礎(chǔ),它是在勘查工作的基礎(chǔ)上,運(yùn)用數(shù)學(xué)和力學(xué)等分析手段,對(duì)滑動(dòng)面以上的滑坡體的穩(wěn)定狀態(tài)進(jìn)行評(píng)價(jià)。目前有關(guān)滑坡穩(wěn)定性的分析方法大致可分為兩類,即極限平衡法和數(shù)值分析法。其中,極限平衡法主要根據(jù)邊坡破壞的邊界條件,應(yīng)用力學(xué)分析的方法,對(duì)可能發(fā)生的滑動(dòng)面,在各種荷載作用下進(jìn)行理論計(jì)算和力學(xué)分析。該方法的理論模型和計(jì)算過程比較簡單,計(jì)算結(jié)果的物理意義非常明確,但由于采用的是剛體極限平衡理論,難以動(dòng)態(tài)反映滑坡堆積體在外界誘發(fā)因素作用下應(yīng)力場(chǎng)的改變以及由其引起的變形情況;加之在模型建立時(shí)假設(shè)的前提條件往往與實(shí)際情況有出入,導(dǎo)致其計(jì)算結(jié)果與實(shí)際情況存在一定差異(宣世進(jìn)等,2009)。而數(shù)值分析法主要是通過數(shù)值模型來模擬邊坡的應(yīng)力分布和變形情況,研究巖土體的應(yīng)力和應(yīng)變的變化過程,求解邊坡的安全系數(shù),并據(jù)此判斷該邊坡是否穩(wěn)定。同極限平衡法相比,數(shù)值分析法克服了邊坡形狀不規(guī)則和材料不均勻的限制,在求解過程中,既考慮了單元體的變形協(xié)調(diào),又考慮了巖土體的破壞準(zhǔn)則,故其計(jì)算結(jié)果更加精確、合理,結(jié)果分析更加方便、直觀(連金芳,2009)。目前常用的數(shù)值分析方法主要有離散元法、有限元法和FLAC法。其中FLAC法也稱連續(xù)介質(zhì)快速拉格朗日差分法(Fast Lagrangian Analysis of Continua),是由美國明尼蘇達(dá)大學(xué)(University of Minnesota)和美國Itasca咨詢集團(tuán)公司(Itasca Consulting Group Inc.)于1986年在拉格朗日算法基礎(chǔ)上開發(fā)的一種顯式有限差分程序,有二維(FLAC2D)和三維(FLAC3D)兩個(gè)版本,該程序有效綜合了離散元法和有限元法的優(yōu)點(diǎn),通過運(yùn)用離散模型方法、動(dòng)態(tài)松弛方法以及有限差分方法,將連續(xù)介質(zhì)的動(dòng)態(tài)演化過程轉(zhuǎn)化為離散節(jié)點(diǎn)的運(yùn)動(dòng)方程和離散單元的本構(gòu)方程進(jìn)行計(jì)算,計(jì)算過程簡單、快速,且能較好地模擬不同介質(zhì)在達(dá)到其強(qiáng)度極限或屈服極限時(shí),所發(fā)生的變形破壞或塑性流動(dòng)等力學(xué)行為,故廣泛應(yīng)用于邊坡穩(wěn)定性分析、支護(hù)工程設(shè)計(jì)、地下洞室開挖、路基填筑以及采礦區(qū)地面沉陷的預(yù)測(cè)等多個(gè)工程地質(zhì)領(lǐng)域。本文以高治村滑坡為例,論述了FLAC3D在滑坡穩(wěn)定性分析評(píng)價(jià)中的應(yīng)用(李朝暉等,2011)。

1 計(jì)算理論

1.1 方法原理

FLAC3D軟件是一種基于快速拉格朗日算法的三維顯式有限差分法程序,是用顯式差分法來求解運(yùn)動(dòng)方程和動(dòng)力方程。該程序?qū)⒂?jì)算區(qū)域內(nèi)的介質(zhì)離散成若干個(gè)計(jì)算單元,每個(gè)單元在給定的邊界條件下都遵循事先指定的線性或非線性本構(gòu)關(guān)系,各單元之間用節(jié)點(diǎn)相互連接,每個(gè)節(jié)點(diǎn)都要受到來自其周圍區(qū)域其它節(jié)點(diǎn)的合力影響,當(dāng)某個(gè)節(jié)點(diǎn)受到其周圍其它節(jié)點(diǎn)的合力后,該節(jié)點(diǎn)就會(huì)失穩(wěn)而發(fā)生運(yùn)動(dòng),其運(yùn)動(dòng)方程可以寫成時(shí)間步長△t的有限差分的形式,從而可以在一個(gè)時(shí)步長中求得速度和位移的增量。對(duì)于任一區(qū)域,均可根據(jù)其周圍節(jié)點(diǎn)的運(yùn)動(dòng)速度而求得它的應(yīng)變率,然后根據(jù)材料的本構(gòu)關(guān)系求得應(yīng)力的增量,再由應(yīng)力的增量求出t和t+△t時(shí)刻各節(jié)點(diǎn)的不平衡力和各節(jié)點(diǎn)在t+△t時(shí)的加速度和積分加速度,進(jìn)而求出各節(jié)點(diǎn)新的位移值和新的應(yīng)力值,如此循環(huán)迭代和計(jì)算,直至不平衡力足夠小或者各節(jié)點(diǎn)的位移趨于平衡(聶廣煥等,2012;叢凱,2012)。

FLAC3D軟件內(nèi)有一條自動(dòng)查找安全系數(shù)的內(nèi)嵌命令-solve fos,該命令就是利用強(qiáng)度折減法,通過不斷內(nèi)插逼近的辦法,來求解邊坡的安全系數(shù)。安全系數(shù)可以定義為使巖土體剛好達(dá)到臨界破壞狀態(tài)時(shí),對(duì)巖土體剪切強(qiáng)度進(jìn)行折減的程度(鄭文棠,2010)。應(yīng)用到有限差分法中可以表述為∶

保持外荷載不變,逐步減小抗剪強(qiáng)度指標(biāo),即將抗剪強(qiáng)度參數(shù)c和tan φ同時(shí)除以折減系數(shù)F,從而得到一組新的強(qiáng)度參數(shù)值c'和tanφ' ,然后將折減后的c'、φ'作為材料新的參數(shù)代入,進(jìn)行有限差分計(jì)算分析,不斷循環(huán)迭代計(jì)算,直至巖土體達(dá)到臨界破壞狀態(tài),此時(shí)采用的折減系數(shù)F即為該巖土體的安全系數(shù)(歐湘萍等,2009;張晶等,2013)。

待求得滑坡的安全系數(shù)后,再依據(jù)《滑坡防治工程勘察規(guī)范》(DZ/T0218-2006)中的評(píng)價(jià)標(biāo)準(zhǔn),判斷該滑坡是否穩(wěn)定。

1.2 計(jì)算流程

利用FLAC3D軟件進(jìn)行滑坡穩(wěn)定性分析的主要流程為∶ ①根據(jù)坡體的主位移剖面,建立其幾何計(jì)算模型;②輸入強(qiáng)度參數(shù),運(yùn)行FLAC3D程序;③根據(jù)事先確定的坡體失穩(wěn)判據(jù),判斷坡體是否達(dá)到臨界破壞狀態(tài),如果達(dá)到臨界破壞狀態(tài)就結(jié)束程序運(yùn)行,輸出計(jì)算結(jié)果,否則就繼續(xù)將折減后的強(qiáng)度參數(shù)作為新的參數(shù)重新輸入,繼續(xù)運(yùn)行程序,直到坡體達(dá)到臨界破壞狀態(tài)為止;④結(jié)合其它分析方法或?qū)崪y(cè)結(jié)果對(duì)運(yùn)行結(jié)果進(jìn)行評(píng)價(jià)。具體流程見圖1。

圖1 利用FLAC3D進(jìn)行滑坡穩(wěn)定性分析工作流程圖Fig.1 Work fowchart of analyzing landslide stability by FLAC3D

2 工程算例

2.1 算例描述

本文選取四川德陽高治村滑坡作為實(shí)例。該滑坡為一小型推移式的土質(zhì)滑坡,滑坡體長約130m,寬約180m,其前緣高程約570m,后緣高程約590m,總體坡向220°~ 235°,坡度11°~16°。滑坡體主要由崩坡積和殘坡積的粉質(zhì)粘土夾碎石組成,厚度約0.4~4.5m,平均厚度約2.0m。前緣(剪出口)厚度約0.4~1.3m,中部厚度約2.0~4.5m,后緣厚度約2.4~3.2m。據(jù)鉆孔及探槽揭露,該滑坡的滑床為白堊系下統(tǒng)古店組(K1g)紅褐色泥巖和砂巖,巖層單層厚度0.15~1.3m,巖層,傾向?yàn)镾185°~ 230°W,傾角為5°~ 10°,傾向于坡向一致。該滑坡的滑動(dòng)帶(面)受基巖頂界面控制,至溝邊臨空,沿殘坡積層與基巖接觸的薄弱帶滑動(dòng)剪出,滑動(dòng)帶(面)深度約0.6~4.2m,傾角約12°~ 15°,滑帶土主要為棕褐色—黃白色的粘性土。根據(jù)實(shí)地調(diào)查,該滑坡是于一次強(qiáng)降雨后開始出現(xiàn)下滑和地表變形、破壞現(xiàn)象,變形和破壞主要發(fā)生在坡體的前緣和中部地帶,并導(dǎo)致了房屋、圍墻以及引水渠等地表建筑物的變形和破壞(四川省廣漢地質(zhì)工程勘察院,2009;北京市地質(zhì)研究所;2009)。

2.2 模型建立

由于該滑坡規(guī)模不大,故本次計(jì)算選取一條沿主位移方向的D—D′剖面(圖2)進(jìn)行簡化處理抽象出一個(gè)二維模型,以巖性界面及滑動(dòng)面為分區(qū)邊界分層概化。

(1)建立計(jì)算模型時(shí),以滑坡水平運(yùn)動(dòng)方向?yàn)閤軸正方向,豎直向上方向?yàn)閥軸正方向。

(2)計(jì)算采用的是比較常見的彈塑性模型,屈服準(zhǔn)則為摩爾-庫倫準(zhǔn)則。模型采用的巖土體力學(xué)參數(shù)主要有內(nèi)聚力、內(nèi)摩擦角、容重、泊松比、彈性模量以及剪切模量。

(3)初始條件中,未考慮構(gòu)造應(yīng)力,僅考慮自重產(chǎn)生的初始地應(yīng)力場(chǎng)。

(4)根據(jù)水文、鉆探資料,本區(qū)地下水位在滑動(dòng)面以下,故計(jì)算時(shí),巖土體強(qiáng)度在天然狀態(tài)下取天然強(qiáng)度,持續(xù)暴雨?duì)顟B(tài)下取飽和強(qiáng)度。

圖2 高治村滑坡D—D'剖面示意圖Fig.2 The D—D'section sketch of GaoZhi Village landslide

(5)由于本次計(jì)算主要是重點(diǎn)分析滑坡在重力作用下的失穩(wěn)變形情況,故模型中未施加水平構(gòu)造應(yīng)力,對(duì)模型左、右側(cè)緣邊界及前、后邊界采用單向水平約束;對(duì)模型底部邊界采用固定約束;只令模型的頂部(即坡面)為自由邊界。

(6)本次計(jì)算將D—D′剖面剖分為4350個(gè)單元,共計(jì)為9042個(gè)節(jié)點(diǎn)。

2.3 參數(shù)選取

本次計(jì)算所選取的巖土體力學(xué)參數(shù)由前期勘察和室內(nèi)試驗(yàn)并參考類似條件下、相同巖性的巖土體的相關(guān)數(shù)據(jù)來確定,具體取值情況見表1。由于影響滑坡穩(wěn)定性系數(shù)的指標(biāo)主要為滑帶土的抗剪強(qiáng)度,故本次計(jì)算主要是通過調(diào)整滑帶土的強(qiáng)度參數(shù)來進(jìn)行滑坡穩(wěn)定性的數(shù)值模擬和分析評(píng)價(jià)。

2.4 結(jié)果分析

本次數(shù)值模擬計(jì)算結(jié)果如圖3—圖5所示。其中圖3和圖4分別為FLAC3D軟件模擬的高治村滑坡在天然狀態(tài)及暴雨?duì)顟B(tài)下D—D′剖面的塑性區(qū)云圖,圖5為高治村滑坡在暴雨?duì)顟B(tài)下D—D′剖面的X方向位移等值線圖。由圖3可以看出,在天然狀態(tài)下,從滑坡的后緣至前緣,塑性區(qū)(圖中紅色部分)只是在坡體的淺表處零星分布,這表明該滑坡在天然狀態(tài)下的塑性變形較弱,坡體比較穩(wěn)定;而由圖4可以看出,在暴雨?duì)顟B(tài)下,從滑坡的后緣至前緣,塑性區(qū)(圖中紅色部分)在坡體的第四系土層(滑坡體)中大量分布,沿土層與基巖的分界面(滑動(dòng)面)微弱貫通,坡體的前緣和中部地帶的塑性變形尤為強(qiáng)烈,這表明該滑坡在暴雨?duì)顩r下處于欠穩(wěn)定狀態(tài),雖然坡體發(fā)生整體失穩(wěn)滑動(dòng)的可能性較小,但在坡體前緣及其中部陡坎處發(fā)生局部失穩(wěn)破壞的可能性極大。

圖3 天然狀態(tài)下D—D′剖面的塑性區(qū)分布圖Fig.3 The plastic zone distribution of D—D section under the natural state

圖4 暴雨然狀態(tài)下D—D′剖面的塑性區(qū)分布圖Fig.4 The plastic zone distribution of D—D′section under the heavy rain state

圖5 暴雨然狀態(tài)下D—D′剖面X方向位移云圖Fig.5 The X-displacement contours of D - D 'section under the heavy rain state

表1 D—D'剖面巖土體特性參數(shù)Tab.1 The characteristic parameters of rock-soil mass of D - D 'section

從圖5中可以看出,在暴雨?duì)顟B(tài)下,坡體的前緣及其中部陡坎處發(fā)生了位移,位移以水平方向?yàn)橹鳎怪狈较虻奈灰撇淮蟆_@表明在暴雨條件下,該坡體的前緣和中部地帶發(fā)生了水平方向的滑動(dòng)變形,變形部位較淺,最大滑動(dòng)位移發(fā)生在坡體中部的陡坎處,這與實(shí)地調(diào)查中所見到的地表變形現(xiàn)象是一致的。

根據(jù)對(duì)D—D′剖面采用強(qiáng)度折減法進(jìn)行計(jì)算所得到的結(jié)果,該坡體在暴雨?duì)顟B(tài)下的穩(wěn)定系數(shù)為 ,依據(jù)《滑坡防治工程勘察規(guī)范》(DZ/T0218-2006)中的評(píng)價(jià)標(biāo)準(zhǔn)判斷,該坡體在暴雨條件下欠穩(wěn)定,易發(fā)生局部的滑動(dòng)變形。這與該坡體實(shí)際發(fā)生的滑動(dòng)變形情況完全吻合(楊波等,2013)。

通過對(duì)高治村滑坡穩(wěn)定性的數(shù)值模擬分析可知,該滑坡在天然狀態(tài)下是穩(wěn)定的,在暴雨?duì)顟B(tài)下是欠穩(wěn)定的,暴雨的作用導(dǎo)致該滑坡巖土體的物理力學(xué)指標(biāo)降低,進(jìn)一步弱化了滑帶土的強(qiáng)度指標(biāo),從而誘發(fā)了坡體局部地帶(前緣及中部陡坎處)的滑移失穩(wěn)(李巖等,2012)。

3 結(jié)語

本文以德陽高治村滑坡為例,論述了FLAC3D在滑坡穩(wěn)定性分析評(píng)價(jià)中的應(yīng)用。實(shí)踐表明,通過FLAC3D對(duì)滑坡穩(wěn)定性進(jìn)行數(shù)值模擬和分析評(píng)價(jià)是合理、可行的。該方法不僅計(jì)算簡單,可迅速求得滑坡穩(wěn)定性系數(shù),而且直觀性強(qiáng),可方便模擬任何地貌、任意形狀滑動(dòng)面滑坡的失穩(wěn)破壞過程,且可采用不同的本構(gòu)關(guān)系,考慮巖土體的非線性,有效克服了傳統(tǒng)極限平衡分析和有限元分析方法的不足。隨著數(shù)值計(jì)算方法和計(jì)算機(jī)技術(shù)的飛速發(fā)展,該方法在工程地質(zhì),特別是巖土工程的力學(xué)分析方面的應(yīng)用更加廣泛,在滑坡穩(wěn)定性的分析和評(píng)價(jià)方面,也必將發(fā)揮更大的作用。

[1]宣世進(jìn),沈萬里. 滑坡穩(wěn)定性評(píng)價(jià)方法研究[J]. 科技資訊,2009,7(1):101.

[2]連金芳. 滑坡穩(wěn)定性評(píng)價(jià)方法評(píng)述[J]. 長江大學(xué)學(xué)報(bào)(自然科學(xué)版),2009,6(4):305~307.

[3]李朝暉,陳建平,馬艷玲. 基于強(qiáng)度折減法的邊坡穩(wěn)定性分析[J]. 采礦技術(shù),2011,11(4):64~65.

[4]聶廣煥,張志增,楊子勝,等. 基于FLAC的南平市紅桃山邊坡穩(wěn)定性數(shù)值分析[J]. 中原工學(xué)院學(xué)報(bào),2012,23(5):51~56.

[5]叢凱. 基于FLAC的牽引式滑坡穩(wěn)定性分析[J]. 甘肅地質(zhì),2012,21(1):69~72.

[6]鄭文棠. 基于FLAC3D的強(qiáng)度折減法和點(diǎn)安全系數(shù)法對(duì)比[J]. 水利與建筑工程學(xué)報(bào),2010,8(4):54~57.

[7]歐湘萍,白楷,朱云升,等. 基于FLAC3D的強(qiáng)度折減法邊坡穩(wěn)定性分析[J]. 武漢理工大學(xué)學(xué)報(bào),2009,31(9):59~61.

[8]張晶,陳劍平,何隆祥,等. 基于FLAC3D的新巖滑坡穩(wěn)定性分析[J]. 中國水運(yùn),2013,13(4):263~265.

[9]四川省廣漢地質(zhì)工程勘察院. 汶川地震災(zāi)區(qū)震后德陽市旌陽區(qū)地質(zhì)災(zāi)害詳細(xì)調(diào)查與區(qū)劃報(bào)告[R]. 2009:1~300.

[10]北京市地質(zhì)研究所. 四川省德陽市旌陽區(qū)和新鎮(zhèn)高治村堰塘庫壩滑坡應(yīng)急治理工程勘查報(bào)告[R]. 2009:1~25.

[11]楊波,鄭杰. 基于FLAC3D的層狀巖質(zhì)邊坡失穩(wěn)數(shù)值模擬[J]. 露天采礦技術(shù),2013,29(1):3~5.

[12]李巖,武磊,李海平. 基于FLAC的北碚滑坡穩(wěn)定性及防治措施分析[J]. 公路與汽運(yùn),2012,28(3):140~143.

Application of FLAC3D in Landslide Stability Analysis and Evaluation

ZHAO Zhonghai
(Beijing Geology Prospecting and Developing Bureau, Beijing 100195)

Taking the Gaozhi Village landslide as an example, this paper discussed the application of FLAC3D in landslide stability analysis and evaluation. The author frst introduced the basic principle and main process of FLAC3D, and then simulated the Gaozhi Village landslide stability under natural state and heavy rain state by Mohr-Coulomb of FLAC3D. The simulation result intuitively showed the distribution characteristics of plastic deformation and horizontal displacement of the slope in the above two conditions, and clearly showed the potential sliding surface inside the slope. According to the result of numerical simulation, the article analyzed and evaluated the landslide stability, studied and discussed the deformation and failure forms of the landslide. By analyzing the plastic deformation and displacement of slope, the author considered that the landslide is basically stable under the natural state and is unstable and likely to slide and lose its stability in the local area under the continuous high intensity rainfall state. This is consistent with the investigation and analysis. This result shows it is feasible and more simple and more intuitive and convenient to analyze and evaluate the stability of landslide by FLAC3D. Compared with the traditional limit equilibrium method, FLAC3D has the incomparable advantage.

Landslide; Stability analysis and evaluation; FLAC3D simulation.

P642.22

A

1007-1903(2016)04-0094-05

10.3969/j.issn.1007-1903.2016.04.018

趙忠海(1970- ),男,教高,從事區(qū)域地質(zhì)調(diào)查、礦產(chǎn)資源勘查、工程物探勘察以及地質(zhì)災(zāi)害評(píng)估、勘查、防治工作。

猜你喜歡
分析
禽大腸桿菌病的分析、診斷和防治
隱蔽失效適航要求符合性驗(yàn)證分析
電力系統(tǒng)不平衡分析
電子制作(2018年18期)2018-11-14 01:48:24
電力系統(tǒng)及其自動(dòng)化發(fā)展趨勢(shì)分析
經(jīng)濟(jì)危機(jī)下的均衡與非均衡分析
對(duì)計(jì)劃生育必要性以及其貫徹實(shí)施的分析
GB/T 7714-2015 與GB/T 7714-2005對(duì)比分析
出版與印刷(2016年3期)2016-02-02 01:20:11
中西醫(yī)結(jié)合治療抑郁癥100例分析
偽造有價(jià)證券罪立法比較分析
在線教育與MOOC的比較分析
主站蜘蛛池模板: 国产欧美日韩专区发布| 欧美在线导航| 久操线在视频在线观看| a级毛片免费看| 成人午夜天| 六月婷婷激情综合| 在线国产91| 国产网站免费看| 国产精品九九视频| 亚洲综合18p| 亚洲人成影院午夜网站| 国产成人艳妇AA视频在线| 日本高清视频在线www色| 亚洲欧美国产五月天综合| 萌白酱国产一区二区| 国产精品制服| 区国产精品搜索视频| 免费人成网站在线观看欧美| 久久99热66这里只有精品一| 欧美另类精品一区二区三区| 在线毛片免费| 亚洲免费福利视频| 五月婷婷导航| 婷婷亚洲综合五月天在线| 午夜激情福利视频| 精品第一国产综合精品Aⅴ| 色有码无码视频| 人妻中文久热无码丝袜| 久久人人爽人人爽人人片aV东京热| 91在线丝袜| 亚洲无码91视频| 91在线高清视频| 欧美日本在线播放| 综合色天天| 亚洲系列中文字幕一区二区| 97青青青国产在线播放| 国产黄色爱视频| 成人字幕网视频在线观看| 亚洲欧美日韩天堂| 国产一级精品毛片基地| 乱人伦视频中文字幕在线| 国产小视频a在线观看| 久久久波多野结衣av一区二区| 91po国产在线精品免费观看| 国产欧美日韩另类| 欧美天堂在线| 一级一毛片a级毛片| AV片亚洲国产男人的天堂| 国产高清国内精品福利| 国产 日韩 欧美 第二页| 亚洲黄色片免费看| 久久久精品无码一区二区三区| 欧美啪啪一区| 亚洲视频在线青青| 中文字幕第1页在线播| a免费毛片在线播放| 狠狠做深爱婷婷综合一区| 国产精品无码AⅤ在线观看播放| 91久久精品日日躁夜夜躁欧美| 久久这里只精品国产99热8| 国产精品99久久久久久董美香| 中文字幕66页| 青青操视频在线| 国产96在线 | 亚洲成av人无码综合在线观看| 多人乱p欧美在线观看| 欧美yw精品日本国产精品| 免费国产小视频在线观看| 99久久国产精品无码| 欧美v在线| 精品一区国产精品| 色综合日本| 欧美日韩理论| 香蕉视频在线精品| 国产精品视频观看裸模| 欧美www在线观看| 国产成人综合久久| 欧美a在线| 日韩午夜伦| 亚洲高清国产拍精品26u| 亚洲五月激情网| 欧美中文一区|