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

        ?

        基于蒙特卡洛隨機(jī)采樣的土石混合體數(shù)值模型構(gòu)建方法

        2020-06-29 04:09:34
        湖南交通科技 2020年2期
        關(guān)鍵詞:混合體土石蒙特卡洛

        (深圳市綜合交通設(shè)計(jì)研究院,廣東 深圳 518003)

        0 引言

        土石混合體(SRM)是一種非常復(fù)雜的不連續(xù)介質(zhì)材料,它是由具有一定尺寸的強(qiáng)度較高的巖塊、強(qiáng)度較低的土體和孔隙等組成的多相體系,是一種遍布全國(guó)乃至世界許多地方且對(duì)工程建設(shè)有重要影響的地質(zhì)材料[1-6]。在工程實(shí)際中,構(gòu)成土石混合體的各種組分在外荷載作用下的力學(xué)性質(zhì)有很大的差異,組分間又存在著極其復(fù)雜的相互作用。這種巖土材料的力學(xué)性能(如應(yīng)力傳遞、破壞模式、裂紋擴(kuò)展、承載能力等)與均質(zhì)巖土體有著較大的差別,且塊石的空間分布、級(jí)配和含石量等關(guān)鍵因素對(duì)SRM的宏觀力學(xué)性質(zhì)有重要影響[7-8]。因此,有必要對(duì)土石混合體的力學(xué)特性進(jìn)行更為準(zhǔn)確的分析。最初,學(xué)者們多采用試驗(yàn)的方法獲取土石混合體的力學(xué)性質(zhì)[9-11],但是在工程尺度上進(jìn)行土石混合體的原位試驗(yàn)和實(shí)驗(yàn)室試驗(yàn)較為困難且耗時(shí)[7],數(shù)值分析為解決這一問(wèn)題提供了有效的方法[12]。構(gòu)建與工程實(shí)際符合程度較高的數(shù)值模型是目前開展土石混合體數(shù)值計(jì)算的關(guān)鍵問(wèn)題[13]。

        塊石顆粒的生成是構(gòu)建土石混合體模型的第一步,以往學(xué)者多用一些簡(jiǎn)單的形體來(lái)模擬塊石。TSESARSKY等[14]基于球形巖石塊,研究巖石的彈性模量和各向同性。GRAZIANI A等[15]基于規(guī)則形狀的巖石塊(如圓形、三角形、矩形等)對(duì)二維SRM隨機(jī)模型進(jìn)行了一系列的數(shù)值模擬來(lái)研究SRM的力學(xué)行為。然而上述方法過(guò)度簡(jiǎn)化了顆粒的真實(shí)形態(tài),對(duì)具有更為復(fù)雜外形的顆粒材料適用性不足?;诖?,XU W J等[16-17]和DU C等[18]應(yīng)用任意凸多邊形來(lái)模擬塊石,這相對(duì)于簡(jiǎn)單形體的塊石有很大的改進(jìn)。羅偉等[19]則基于AutoCAD 二次開發(fā)開展了土石混合體隨機(jī)結(jié)構(gòu)模型生成與直剪強(qiáng)度數(shù)值試驗(yàn)研究。喻江武等[6]運(yùn)用傅里葉變換方法生成特定形狀的塊石并構(gòu)建土石混合體模型。趙鑫曜等[20]通過(guò)MATLAB軟件對(duì)大量的塊石樣本的數(shù)字圖像進(jìn)行邊界提取、平滑、規(guī)則化等處理,得到塊石形狀坐標(biāo)和塊石形狀參數(shù),構(gòu)建了塊石庫(kù),從而建立土石混合體模型。這與實(shí)際情況比較貼合,但構(gòu)建塊石庫(kù)需要耗費(fèi)大量的精力。WANG Z M等[8]提出了一種基于蒙特卡洛隨機(jī)采樣原理生成任意形狀隨機(jī)骨料結(jié)構(gòu)的程序,開展混凝土的細(xì)觀研究。本文在WANG Z M等[8]的基礎(chǔ)上加以改進(jìn),分別基于圓形、矩形、三角形采用蒙特卡洛隨機(jī)采樣來(lái)進(jìn)行任意形狀二維塊石顆粒的構(gòu)建。

        在構(gòu)建土石混合體模型的過(guò)程中,塊石間的位置關(guān)系判定是極其重要的一個(gè)步驟。2個(gè)塊石在二維平面上的位置關(guān)系有相離、接觸、相交、內(nèi)含4種情況。根據(jù)自然界的實(shí)際情況,在進(jìn)行土石混合體模型的構(gòu)建時(shí),需要保證塊石之間彼此不能相交和內(nèi)含[13,16,21]。簡(jiǎn)單形體間(圓形、三角形、矩形等)的位置關(guān)系判定較為簡(jiǎn)單。XU W J等[16]采用隨機(jī)凸多邊形來(lái)模擬塊石并基于Area criterion和Sum of the included angle criterion方法進(jìn)行隨機(jī)凸多邊形塊石的相交判定;MENG Q X等[13]基于蒙特卡洛隨機(jī)采樣原理生成任意形狀的塊石顆粒,但并未對(duì)塊石顆粒間的位置關(guān)系判定方法進(jìn)行說(shuō)明,也沒(méi)有其他學(xué)者對(duì)任意形狀塊石位置關(guān)系的判定算法進(jìn)行研究。

        本文基于蒙特卡洛隨機(jī)采樣構(gòu)建的塊石,提出了一種任意形狀塊石的位置關(guān)系判定算法,突破了以往在進(jìn)行塊石位置關(guān)系判定時(shí)對(duì)塊石“凸性”的限制;介紹了一種構(gòu)建懸浮結(jié)構(gòu)土石混合體模型的方法,并對(duì)其進(jìn)行了直剪試驗(yàn)分析。

        1 任意形狀塊石顆粒的隨機(jī)生成

        本文基于蒙特卡洛隨機(jī)采樣方法進(jìn)行二維塊石的構(gòu)建,基本方法是:首先規(guī)定一個(gè)簡(jiǎn)單的初始形狀(如矩形);然后在初始形狀的輪廓上進(jìn)行點(diǎn)的隨機(jī)采樣,并通過(guò)改變采樣點(diǎn)與初始形狀形心之間的距離來(lái)得到一系列不規(guī)則形狀的塊石。具體實(shí)現(xiàn)過(guò)程如下:

        1)生成一組隨機(jī)角度增量:

        {Δθ1,Δθ2,Δθ3,…,Δθn-1}

        其中:

        (i=1,2,…,n-1)

        (1)

        (2)

        式中:ζ、η均是介于0和1之間的系數(shù)。當(dāng)η很小時(shí),表示角度增量很小,隨機(jī)角度更均勻。

        在隨機(jī)角度增量生成的過(guò)程中,需保證:

        (3)

        3) 求得從原點(diǎn)到隨機(jī)采樣點(diǎn)的距離為{d1,d2,d3,…,dn},對(duì)于每一條極徑,在di-Δd到di+Δd范圍內(nèi)生成一組新的隨機(jī)長(zhǎng)度序列{D1,D2,D3,…,Dn},則新的極徑可以描述為:

        Di=di+(2ζ-1)Δd(i=1,2,…,n)

        (4)

        式中:ζ是0到1之間的隨機(jī)數(shù);di是簡(jiǎn)單形體極徑的長(zhǎng)度;Di是新的極徑長(zhǎng)度;Δd是極徑的變化幅度。

        新的隨機(jī)長(zhǎng)度序列{D1,D2,D3,…,Dn}對(duì)應(yīng)一組新的采樣點(diǎn){P1,P2,P3,…,Pn}。

        4)在極坐標(biāo)系中使用{D1,D2,D3,…,Dn}和{Δθ1,Δθ2,Δθ3,…,Δθn}來(lái)生成塊石輪廓,并使用式(5)將極坐標(biāo)系轉(zhuǎn)換為笛卡爾坐標(biāo)系。

        (5)

        蒙特卡洛隨機(jī)采樣生成塊石的過(guò)程如圖1所示,其中:

        (6)

        2 塊石位置關(guān)系判定

        2個(gè)塊石的位置關(guān)系有相離、接觸、相交、內(nèi)含4種。本節(jié)通過(guò)構(gòu)造包圍盒來(lái)進(jìn)行塊石位置關(guān)系的初步判定,若包圍盒相離,則2個(gè)塊石相離;若包圍盒相交、接觸或內(nèi)含,則需通過(guò)判斷構(gòu)成2個(gè)塊石輪廓的線段之間的位置關(guān)系來(lái)進(jìn)一步判定2個(gè)塊石的位置關(guān)系。

        2.1 基于包圍盒的塊石位置關(guān)系初步判定

        包圍盒是一個(gè)能將復(fù)雜形態(tài)的塊石包圍起來(lái)的盒子。為了提高算法的效率,迅速排除相離及內(nèi)含的塊石,通過(guò)構(gòu)造包圍盒初步判定塊石的位置關(guān)系是一種較好的方法。

        2.1.1包圍盒的構(gòu)建

        較為常用的包圍盒有AABB(Axis Aligned Bounding Box)包圍盒[22]、OBB(Oriented Bounding Box)包圍盒[23]和球包圍盒[24-25]。AABB包圍盒是能夠包圍原物體的最小軸向包圍盒;OBB包圍盒則在AABB包圍盒的基礎(chǔ)上考慮了物體的朝向;球包圍盒則是能夠包圍原物體的最小圓形。圖2直觀地顯示了三者的區(qū)別。

        圖2 3種包圍方式對(duì)比

        就構(gòu)造方式來(lái)看,AABB包圍盒最簡(jiǎn)單,球包圍盒次之,OBB包圍盒最復(fù)雜;就包圍盒與塊石的符合程度來(lái)看,OBB包圍盒最為緊密,球包圍盒次之,AABB包圍盒最為疏松。考慮到算法的效率,選用AABB包圍盒進(jìn)行塊石初步位置關(guān)系的判定。

        AABB包圍盒的構(gòu)建方法如圖3所示,可以得到包圍盒1的4個(gè)頂點(diǎn)坐標(biāo):A1、B1、C1、D1和兩邊長(zhǎng):d1x=x1max-x1min,d1y=y1max-y1min。同理可以得到包圍盒2的4個(gè)頂點(diǎn)坐標(biāo):A2、B2、C2、D2和兩邊長(zhǎng):d2x=x2max-x2min,d2y=y2max-y2min。

        圖3 AABB包圍盒構(gòu)建

        2.1.2包圍盒位置關(guān)系初步判定

        得到2個(gè)塊石的包圍盒之后,通過(guò)構(gòu)造能將2個(gè)塊石同時(shí)包圍的公共包圍盒(如圖4所示),來(lái)判定2個(gè)包圍盒相離、接觸、相交、內(nèi)含4種位置關(guān)系。

        圖4 AABB包圍盒位置關(guān)系判斷

        具體方法如下:首先提取出公共包圍盒的2個(gè)邊長(zhǎng):dx=xmax-xmin,dy=ymax-ymin,其中:xmin=min(x1min,x2min);ymin=min(y1min,y2min);xmax=max(x1max,x2max);ymax=max(y1max,y2max)。然后采用以下條件進(jìn)行2個(gè)包圍盒位置關(guān)系判定。

        1) 當(dāng)d1x+d2x>dx且d1y+d2y>dy時(shí),則這2個(gè)包圍盒相交。

        2) 當(dāng)滿足3個(gè)條件之一時(shí),2個(gè)包圍盒接觸:

        ①d1x+d2x=dx且d1y+d2y=dy;

        ②d1x+d2x=dx且d1y+d2y>dy;

        ③d1x+d2x>dx且d1y+d2y=dy。

        3) 當(dāng)d1x+d2x

        4) 當(dāng)滿足下列情況之一時(shí),2個(gè)包圍盒內(nèi)含。

        ①x1min>x2min,x1maxy2min,y1max

        ②x2min>x1min,x2maxy1min,y2max

        2.2 塊石位置關(guān)系精確判定

        塊石位置關(guān)系有相離、相交、接觸、內(nèi)含4種情況,具體判定方法如下:

        1) 若2個(gè)包圍盒相離,則2個(gè)塊石必定相離。

        2) 若2個(gè)包圍盒接觸且接觸面上2個(gè)塊石也接觸,則2個(gè)塊石相接觸,如圖5a所示;否則,2個(gè)塊石相離,如圖5b所示。

        3) 若2個(gè)包圍盒相交,先求出2個(gè)包圍盒的相交包圍盒(定義為2個(gè)包圍盒的相交部分,圖6中虛線部分),然后遍歷相交包圍盒內(nèi)構(gòu)成2個(gè)塊石輪廓的線段進(jìn)行相交判斷[25]。若相交包圍盒構(gòu)成2個(gè)塊石的線段均不相交,則這2個(gè)塊石相離(見圖6a);若相交包圍盒構(gòu)成2個(gè)塊石的線段存在相交的情況,且某塊石上存在一頂點(diǎn)在另一塊石的內(nèi)部,則這2個(gè)塊石相交(見圖6b);否則,這2個(gè)塊石接觸(見圖6c)。

        圖5 包圍盒接觸塊石位置關(guān)系判斷示意圖

        4) 若2個(gè)包圍盒內(nèi)含,同樣先構(gòu)造出2個(gè)包圍盒的相交包圍盒,即為較小的包圍盒;若包圍盒相交則用“3)”中的方法可以判定出塊石的相交(圖7b)和接觸(圖7c);若2個(gè)包圍盒不相交則判斷塊石2(較小的塊石)是否有一點(diǎn)在塊石1(較大的塊石)內(nèi),若是,則2個(gè)塊石內(nèi)含(圖7d),否則2個(gè)塊石相離(圖7a)。

        圖6 包圍盒相交時(shí)塊石相離、相交或接觸示意圖

        圖7 包圍盒內(nèi)含時(shí)塊石位置關(guān)系判斷

        3 土石混合體模型構(gòu)建

        本文基于MATLAB隨機(jī)生成方法進(jìn)行土石混合體模型的構(gòu)建,將蒙特卡洛隨機(jī)采樣生成的塊石按照特定的級(jí)配和含石量隨機(jī)投放進(jìn)指定的區(qū)域內(nèi)[26]。在投放的過(guò)程中,采用第2節(jié)的方法判定此塊石與區(qū)域內(nèi)已有塊石的位置關(guān)系,當(dāng)新投的塊石與已有塊石相離時(shí),判定該塊石有效,否則舍棄此塊石,重新投放,直至含石量達(dá)到要求或投放次數(shù)達(dá)到閾值[13]。

        3.1 塊石的隨機(jī)投放

        本文按照特定的級(jí)配曲線和含石量在一定區(qū)域內(nèi)進(jìn)行塊石的投放,未被塊石覆蓋的區(qū)域作為土體。在投放時(shí),首先隨機(jī)生成一個(gè)粒徑為d(d在限定的范圍內(nèi))的塊石,然后從塊石庫(kù)中隨機(jī)取出1個(gè)塊石M(xi,yi), 在前面塊石庫(kù)建立的塊石中心o的坐標(biāo)為:

        (7)

        其中xi、yi為邊界上各點(diǎn)的坐標(biāo),將各塊石以中心o隨機(jī)縮放,縮放系數(shù)為η,然后將塊石的中心點(diǎn)坐標(biāo)o(xo,yo)隨機(jī)插在區(qū)域內(nèi)的一點(diǎn)s(xs,ys)上,則最終塊石在土石混合體區(qū)域內(nèi)的坐標(biāo)為:

        (8)

        用第2節(jié)的方法進(jìn)行塊石的相交、非相交判定,若新插入的塊石與已有的塊石不相交,則存儲(chǔ)此塊石,否則,舍棄此塊石,重新投放。

        3.2 投放終點(diǎn)的確定

        在投放過(guò)程中,首先設(shè)定一個(gè)含石率P(定義為塊石面積與土石混合體面積的比值)。則對(duì)于每一粒徑范圍內(nèi)的塊石Mi,當(dāng)Mi達(dá)到此粒徑范圍內(nèi)塊石的含石率Pi時(shí),Mi投放結(jié)束,開始投放下一粒徑的塊石Mi+1。

        由于每次塊石的投放點(diǎn)都是隨機(jī)的,當(dāng)邊坡內(nèi)的塊石數(shù)量達(dá)到一定值時(shí),可能會(huì)出現(xiàn)每次新投放的塊石總會(huì)與邊坡內(nèi)已有的塊石相交的情況。因此,針對(duì)這種情況,本文設(shè)定了一個(gè)閾值,當(dāng)投放次數(shù)達(dá)到此閾值時(shí),即使沒(méi)有達(dá)到規(guī)定的含石率,依然終止投放。

        本文建立懸浮結(jié)構(gòu)土石混合體邊坡的詳細(xì)流程如圖8所示。

        圖8 土石混合體模型構(gòu)建流程圖

        3.3 有限差分?jǐn)?shù)值模型的生成

        本文使用MIDAS GTS來(lái)對(duì)懸浮結(jié)構(gòu)的土石混合體進(jìn)行網(wǎng)格劃分。首先將MATLAB生成的土石混合體模型導(dǎo)入CAD并保存為dxf文件;然后將dxf文件導(dǎo)入MIDAS進(jìn)行有限元網(wǎng)格劃分;網(wǎng)格劃分好后,將模型在MIDAS中沿垂直模型平面的方向拉伸1 mm,形成1個(gè)假三維模型。

        4 土石混合體直剪試驗(yàn)分析

        本節(jié)共構(gòu)造3個(gè)土石混合體模型(見圖9),3個(gè)模型級(jí)配不同,但尺寸與含石量均相同(含石量為40%)。模型級(jí)配如圖10a所示,級(jí)配3土石混合體模型中大粒徑塊石顆粒的含量最多,級(jí)配2次之,級(jí)配1最少。模型尺寸、邊界條件及加載方式見圖10b,計(jì)算用土石材料參數(shù)見表1。

        采用FLAC3D對(duì)土石混合體與均質(zhì)土體進(jìn)行直剪試驗(yàn),其計(jì)算結(jié)果如圖11所示。計(jì)算結(jié)果表明:均質(zhì)土體的剪切帶為一條規(guī)則的直線;土石混合體的剪切帶則更加寬厚和曲折,且土石混合體在剪切過(guò)程中存在著明顯的繞石效應(yīng),即塑性區(qū)沿著塊石之間的土體部分進(jìn)行擴(kuò)展。

        a)模型1

        b)模型

        c)模型3

        圖10 級(jí)配曲線與模型試驗(yàn)條件

        a)均質(zhì)模型b)土石混合體模型1

        c)土石混合體模型

        2d)土石混合體模型3

        實(shí)驗(yàn)所得土石混合體和均質(zhì)土體的剪切應(yīng)力與剪切位移關(guān)系曲線如圖12所示。

        圖12 土石混合體與均質(zhì)土體剪切應(yīng)力-應(yīng)變曲線

        由圖12可得:在剪切的初始階段,均質(zhì)土體的剪切應(yīng)力隨剪切應(yīng)變線性增長(zhǎng),此時(shí)土體處于彈性變形階段;在彈性變形階段之后、到達(dá)峰值強(qiáng)度之前,曲線增長(zhǎng)速率不斷減小,曲線段逐漸變得平緩,即進(jìn)入初始屈服階段;在初始屈服階段,土體的剪切應(yīng)力繼續(xù)增加,達(dá)到峰值后,土體發(fā)生破壞,其剪切應(yīng)力開始降低,并逐漸達(dá)到相應(yīng)的殘余應(yīng)力。這與XU W J等[1]和GONG J等[4]的研究結(jié)果一致。相比均質(zhì)土體,土石混合體剪切過(guò)程中還會(huì)存在若干個(gè)局部剪切階段,即應(yīng)力-應(yīng)變曲線呈現(xiàn)多次上升、下降,這是因?yàn)橥潦旌象w中存在著塊石與塊石、塊石與土體之間的摩擦、咬合、滑動(dòng)等情況。由圖12可知,局部剪切階段可能會(huì)發(fā)生在峰值應(yīng)力出現(xiàn)前,也可能發(fā)生在峰值應(yīng)力出現(xiàn)之后,并且不同土石混合體模型局部剪切階段出現(xiàn)的時(shí)間和次數(shù)均不相同,這可能與塊石的空間分布有關(guān)。

        由于塊石的存在,土石混合體內(nèi)部的相互作用較為復(fù)雜,因此,土石混合體的剪切過(guò)程較均質(zhì)土體更為緩慢,到達(dá)殘余應(yīng)力時(shí)土石混合體的水平位移更大。在剪切的過(guò)程中,土石混合體的剪切應(yīng)力一直大于均質(zhì)土體,且峰值應(yīng)力遠(yuǎn)大于均質(zhì)土體,這說(shuō)明塊石的存在使得土體的抗剪強(qiáng)度增大,塊石對(duì)土體強(qiáng)度起到了一種增強(qiáng)作用。另外,通過(guò)對(duì)不同級(jí)配的土石混合體模型進(jìn)行對(duì)比分析,發(fā)現(xiàn)在含石量相同的情況下,級(jí)配3的土石混合體模型的峰值應(yīng)力和殘余應(yīng)力最大,級(jí)配2次之,級(jí)配1最小。這說(shuō)明在含石量相同的情況下,大粒徑塊石顆粒的含量越多,土石混合體的抗剪強(qiáng)度越高,即大粒徑塊石顆粒的存在有助于提高土石混合體的抗剪強(qiáng)度。

        5 結(jié)論

        本文采用蒙特卡洛隨機(jī)采樣進(jìn)行塊石隨機(jī)生成,并基于AABB包圍盒提出了一種任意形狀塊石位置關(guān)系判定算法,可以有效處理以往研究中在塊石位置關(guān)系判定方面對(duì)塊石凸性的限制?;谠撍惴ㄟM(jìn)行土石混合體模型構(gòu)建并對(duì)3種級(jí)配不同、含石量相同的模型進(jìn)行直剪數(shù)值試驗(yàn)分析,結(jié)果表明:

        1) 均質(zhì)土體剪切帶為一條規(guī)則直線,土石混合體剪切帶由于塊石的存在而變得更加寬厚和曲折。土石混合體剪切帶上存在明顯繞石效應(yīng),即塑性區(qū)沿塊石之間的土體部分進(jìn)行擴(kuò)展。

        2) 相比均質(zhì)土體,土石混合體剪切過(guò)程剪應(yīng)力-應(yīng)變曲線除了具有近似線性的彈性階段、達(dá)到峰值強(qiáng)度前的初始屈服階段、剪切破壞面逐漸形成的剪切破壞階段、剪切變形持續(xù)增加的殘余變形階段4個(gè)階段之外,還會(huì)存在若干個(gè)局部剪切階段,即應(yīng)力-應(yīng)變曲線呈現(xiàn)多次上升、下降。

        3) 土石混合體抗剪強(qiáng)度較均質(zhì)土體有較大提升,表明塊石存在對(duì)土體強(qiáng)度起到了一種增強(qiáng)作用。并且在含石量相同的情況下,大粒徑塊石顆粒的含量越多,土石混合體的抗剪強(qiáng)度越高,即大粒徑顆粒的增強(qiáng)作用更明顯。此外,土石混合體剪切過(guò)程較均質(zhì)土體更為緩慢,到達(dá)殘余應(yīng)力時(shí)的水平位移更大。

        猜你喜歡
        混合體土石蒙特卡洛
        強(qiáng)風(fēng)化地層上土石圍堰的截流結(jié)構(gòu)設(shè)計(jì)及施工
        征服蒙特卡洛賽道
        市政道路高填方區(qū)域土石混填技術(shù)研究
        壓實(shí)度對(duì)筑堤用土石混合料力學(xué)特性影響研究
        治淮(2018年1期)2018-02-09 09:15:01
        利用控制變量方法縮減蒙特卡洛方差
        不同水環(huán)境下土石混填地基水平推剪試驗(yàn)研究
        克制的宣泄——鋼琴家普列特涅夫:浪漫主義、現(xiàn)代審美和個(gè)人風(fēng)格的混合體
        獨(dú)領(lǐng)風(fēng)騷的夏日
        蒙特卡洛模擬法計(jì)算電動(dòng)汽車充電負(fù)荷
        基于蒙特卡洛的非線性約束條件下的優(yōu)化算法研究
        午夜免费观看一区二区三区| 国产免费专区| 国产女人体一区二区三区| 丝袜美腿亚洲综合一区| 插鸡网站在线播放免费观看| 97精品国产手机| 亚洲AV无码精品色午夜超碰| 日本一区二区在线播放| 一区二区三区国产精品| 粉嫩人妻91精品视色在线看| 国产av激情舒服刺激| 亚洲av午夜国产精品无码中文字| 亚洲色欲久久久久综合网| 国产资源在线视频| 久久99免费精品国产| (无码视频)在线观看| 亚洲欧美日韩在线一区| 丝袜美女污污免费观看的网站| 男女视频网站免费精品播放| 国产精品亚洲av无人区一区香蕉 | 亚洲一区二区在线视频播放| 男人一插就想射的原因| 欧美国产激情二区三区| 欧美国产日本高清不卡| 杨幂Av一区二区三区| 日本午夜a级理论片在线播放| 亚洲av丰满熟妇在线播放| 欧美黑人群一交| 揄拍成人国产精品视频| 亚洲国产精品自产拍久久蜜AV | 久久精品国产夜色| 中文字幕一区二区三在线| 水野优香中文字幕av网站| 国产成人午夜精华液| 亚洲欧美日韩国产综合专区| 日韩中文字幕不卡在线| 国产人妻丰满熟妇嗷嗷叫| 欧美亚洲高清日韩成人| 亚洲蜜臀av一区二区三区漫画| 亚洲欧洲成人精品香蕉网| 精品无码国产污污污免费网站|