• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    基于真實細(xì)觀結(jié)構(gòu)的FFT數(shù)值方法對巖石材料非線性力學(xué)行為的研究

    2022-06-15 08:30:08李明耀彭磊左建平王智敏李紹金薛喜仁
    礦業(yè)科學(xué)學(xué)報 2022年4期
    關(guān)鍵詞:細(xì)觀方解石石英

    李明耀彭磊左建平王智敏李紹金薛喜仁

    中國礦業(yè)大學(xué)(北京)力學(xué)與建筑工程學(xué)院,北京 100083

    巖石是自然界常見的非均質(zhì)材料,其內(nèi)部包含許多細(xì)觀結(jié)構(gòu),如礦物夾雜、微裂隙和微孔洞等[1-2]。 這些細(xì)觀結(jié)構(gòu)特征改變了巖石的局部應(yīng)力分布,進(jìn)而影響其宏觀力學(xué)行為和破壞過程[3]。因此,在研究巖石的宏觀力學(xué)問題時,必須考慮其細(xì)觀結(jié)構(gòu)的影響。 由于巖石的細(xì)觀結(jié)構(gòu)具有很大的隨機(jī)性,而且不容易觀測,因此很難通過實驗直接建立巖石細(xì)觀結(jié)構(gòu)與宏觀力學(xué)性能的聯(lián)系。 通過數(shù)值模擬分析巖石細(xì)觀結(jié)構(gòu)與其宏觀力學(xué)性能的聯(lián)系已經(jīng)成為一種趨勢[4-5]。

    目前,數(shù)值方法被廣泛應(yīng)用于研究巖石宏觀性能和細(xì)觀結(jié)構(gòu)之間的聯(lián)系,常用的數(shù)值方法主要有基于連續(xù)體、不連續(xù)體的方法和混合方法等[6-7]。 有限元法(FEM)和有限差分法(FDM)都是基于連續(xù)體的數(shù)值方法,被廣泛應(yīng)用于模擬巖石的宏觀力學(xué)行為。 于慶磊等[8]采用有限元法對花崗巖進(jìn)行了單軸壓縮數(shù)值模擬實驗,研究了巖石細(xì)觀結(jié)構(gòu)對其變形和強(qiáng)度的影響。 朱澤奇等[9]采用有限差分法模擬了花崗巖的破壞過程,研究了巖石細(xì)觀結(jié)構(gòu)對其破壞過程的影響。 然而,有限元法和有限差分法在建立數(shù)值模型時需要生成大量的網(wǎng)格,巖石細(xì)觀結(jié)構(gòu)的幾何形狀是復(fù)雜多變的,這對網(wǎng)格生成技術(shù)造成了較大的困難。

    基于不連續(xù)體的數(shù)值方法,顆粒流分析法(PFC)最為常見。 Hu 等[10]采用了基于顆粒的離散元模型模擬了北山花崗巖的損傷過程,研究了北山花崗巖的宏觀力學(xué)和微裂紋行為。 韓振華等[11]采用基于顆粒流的離散元模型模擬了不同粒徑花崗巖的破壞過程,研究了礦物粒徑對其破壞形式和裂紋形態(tài)的影響。 PFC 方法的優(yōu)勢是可以有效還原巖石的細(xì)觀結(jié)構(gòu),并且可以通過統(tǒng)計化學(xué)鍵的斷裂數(shù)量來模擬巖石細(xì)觀結(jié)構(gòu)的破壞過程。 但是,PFC 方法沒有考慮顆粒本身的變形和失效過程,而且在模型參數(shù)標(biāo)定方面也存在一定的困難。

    近年來,基于快速傅里葉變換(Fast Fourier Transforms,FFT)的數(shù)值方法作為一種無網(wǎng)格方法,被逐漸應(yīng)用于巖石材料的數(shù)值模擬。 Li 等[12]采用基于FFT 的數(shù)值方法研究了脆性巖石的細(xì)觀結(jié)構(gòu)對其宏觀力學(xué)行為和破壞過程的影響。Escoda 等[13]利用基于FFT 的數(shù)值方法研究了混凝土細(xì)觀結(jié)構(gòu)對其局部應(yīng)力和宏觀力學(xué)性能的影響。FFT 方法作為一種無網(wǎng)格的數(shù)值方法,可以更真實地還原巖石的細(xì)觀結(jié)構(gòu),簡化了各種復(fù)雜細(xì)觀結(jié)構(gòu)的幾何形狀的網(wǎng)格化過程。 同時,作為一種可以直接進(jìn)行全場模擬的數(shù)值方法,有大量的開源數(shù)據(jù)庫的支持,而且對于求解具有復(fù)雜細(xì)觀結(jié)構(gòu)的問題仍具有較高的計算效率[14]。

    在過去的研究中,通常通過定義隨機(jī)參數(shù)將巖石的非均質(zhì)性納入到數(shù)值模型中[15]。 例如,將巖石的細(xì)觀結(jié)構(gòu)簡化為隨機(jī)分布的顆粒形狀來研究巖石非均質(zhì)性的影響。 Li 等[16]利用計算機(jī)重建的方法復(fù)制巖石初始的細(xì)觀結(jié)構(gòu),以此來考慮巖石的非均質(zhì)性,并探究了黏土巖的非均質(zhì)性對其破壞過程的影響。 這些方法雖然在一定程度上考慮了巖石非均質(zhì)性的影響,但是隨機(jī)參數(shù)的選取具有很大的主觀性,并且局限于統(tǒng)計分布。 隨著X 射線衍射(XRD)、掃描電子顯微鏡(SEM)和計算機(jī)斷層掃描(CT)等實驗技術(shù)的發(fā)展,為數(shù)字圖像處理技術(shù)(DIP)在數(shù)值模擬中的應(yīng)用提供了條件。 Yue 等[17]提出了利用DIP 技術(shù)將巖石細(xì)觀結(jié)構(gòu)納入到有限元模型的方法,并研究了巖石材料的非均質(zhì)性。 Wu等[18]提出了一種基于數(shù)字圖像處理的FDEM 模型,研究花崗巖細(xì)觀結(jié)構(gòu)對其宏觀力學(xué)行為和破壞過程的影響。 Tan 等[19]利用DIP 技術(shù)和Voronoi 鑲嵌過程將花崗巖的細(xì)觀結(jié)構(gòu)細(xì)分成隨機(jī)的子塊,通過控制子塊的尺寸和分布來考慮巖石的非均質(zhì)性,以此研究了花崗巖的非均質(zhì)性對其破壞過程的影響。

    上述研究考慮了巖石非均質(zhì)性的影響,并且提出了許多方法來研究巖石材料細(xì)觀結(jié)構(gòu)對其宏觀力學(xué)性能的影響。 但是,有限元法(FEM)和有限差分法(FDM)等[8-9]在建立模型時需要生成大量復(fù)雜的網(wǎng)格,顆粒流分析法(PFC)[20]無法考慮顆粒本身的變形和失效過程,而且對模型參數(shù)的定義存在一定的局限性。 采用隨機(jī)生成的方法[21-22]將巖石的細(xì)觀結(jié)構(gòu)等效為隨機(jī)分布不同形狀的顆粒,只能近似模擬,無法真實還原巖石的細(xì)觀結(jié)構(gòu),而且這種考慮巖石非均質(zhì)性的方法局限于統(tǒng)計分布的參數(shù)。

    鑒于此,為了更準(zhǔn)確地研究巖石材料的真實細(xì)觀結(jié)構(gòu)對其宏觀非線性力學(xué)行為的影響,本文提出了一種可以直接與數(shù)字圖像有機(jī)結(jié)合的FFT 數(shù)值方法,主要內(nèi)容包括:

    (1) 巖石細(xì)觀結(jié)構(gòu)的數(shù)字圖像表征。 對巖石表面圖像進(jìn)行數(shù)字圖像處理,獲取其細(xì)觀結(jié)構(gòu)表征圖像。

    (2) 建立真實細(xì)觀結(jié)構(gòu)的數(shù)值模型。 利用FFT 數(shù)值方法可直接導(dǎo)入像素信息的優(yōu)勢,得到還原巖石真實細(xì)觀結(jié)構(gòu)的數(shù)值模型。

    (3) 建立合理的細(xì)觀本構(gòu)模型。 根據(jù)巖石內(nèi)部不同礦物所表現(xiàn)出的局部力學(xué)性能,建立合理的局部本構(gòu)模型,如線彈性模型、彈塑性模型、損傷模型等。

    (4) 力學(xué)計算及分析應(yīng)用。 利用FFT 數(shù)值方法進(jìn)行力學(xué)計算和全場數(shù)值模擬,得到巖石細(xì)觀的應(yīng)力場的演化過程,并與實驗數(shù)據(jù)對比分析,得到細(xì)觀結(jié)構(gòu)與宏觀力學(xué)行為的內(nèi)在聯(lián)系。

    1 巖石細(xì)觀結(jié)構(gòu)的數(shù)字圖像表征

    數(shù)字圖像處理技術(shù)(Digital Image Processing,DIP)是利用計算機(jī)將圖像轉(zhuǎn)化為數(shù)字矩陣,然后利用圖像處理方法對數(shù)字矩陣進(jìn)行分析和處理的過程[23]。 相比于用隨機(jī)參數(shù)表征巖石非均質(zhì)性的方法[15,24],DIP 技術(shù)能夠更真實地模擬巖石細(xì)觀結(jié)構(gòu),是將巖石細(xì)觀結(jié)構(gòu)納入數(shù)值模型更合理的方法。 由于巖石類材料其內(nèi)部不同的礦物成分在數(shù)字圖像上表現(xiàn)為不同的顏色,因此DIP 技術(shù)在巖石表面圖像處理方面存在特定的優(yōu)勢,在獲取巖石內(nèi)部信息中得到了廣泛應(yīng)用[17]。

    1.1 巖樣的礦物成分與分布

    本文以法國ANDRA 地下研究實驗室的Callovo-Oxfordian(COx)黏土巖為研究對象[25],詳細(xì)闡述所提出的基于真實細(xì)觀結(jié)構(gòu)的FFT 數(shù)值方法在巖石力學(xué)領(lǐng)域的應(yīng)用。 圖1 為COx 黏土巖樣品在掃描電鏡下采集的灰度圖像,分辨率為250 μm。 COx 黏土巖中的礦物因化學(xué)成分不同而具有明顯不同的顏色和顆粒大小,主要由方解石(淺灰色大顆粒)、石英(深色大顆粒)和黏土(黑色散狀顆粒)3 種礦物組成。 此外,還有少量其他礦物和孔洞(因其含量很少,可以忽略),本文根據(jù)力學(xué)性質(zhì)將其歸于性質(zhì)較接近的區(qū)域。

    圖1 COx 黏土巖灰度圖像Fig.1 Gray image of the COx claystone

    圖1 所示COx 黏土巖的灰度圖像中,雖然各種礦物可以根據(jù)圖像中的顏色大致區(qū)分為3 類,但同一礦物成分的灰度在特定區(qū)域仍然有深有淺,而且有些區(qū)域不同礦物之間的灰度值較為接近,因此很難嚴(yán)格確定不同礦物的區(qū)域和邊界。 灰度圖像是由大量像素組成的二維矩陣,每個像素具有確定的位置信息和亮度信息,灰度值的取值范圍為[0,255],表示像素由黑色到白色的變化過程。 由圖1 可以看出,COx 黏土巖灰度圖像中,石英和黏土區(qū)的灰度值較為接近,而方解石區(qū)與其余區(qū)域差異明顯。 因此,為了準(zhǔn)確描述巖石內(nèi)部細(xì)觀結(jié)構(gòu),必須對灰度圖像進(jìn)行處理,從而嚴(yán)格區(qū)分出3 種不同礦物的區(qū)域和邊界。

    1.2 巖石細(xì)觀結(jié)構(gòu)的數(shù)字圖像表征過程

    基于上述黏土巖灰度圖像特點的分析,為了更好地區(qū)分黏土巖內(nèi)部的方解石、石英和黏土礦物,對其細(xì)觀結(jié)構(gòu)的數(shù)字圖像需分2 個步驟進(jìn)行處理:

    (1) 從整體圖像中表征方解石礦物;

    (2) 在其余區(qū)域表征石英礦物和黏土礦物。

    表征過程如圖2 所示。

    圖2 基于COx 黏土巖細(xì)觀結(jié)構(gòu)表征的流程Fig.2 Flow chart based on microstructural characterization of COx clay rock

    1.2.1 方解石礦物區(qū)域的表征

    如圖1 所示,在灰度圖像中,方解石區(qū)呈淺灰色,而石英和黏土區(qū)均呈深黑色。 方解石區(qū)與其他區(qū)域呈現(xiàn)出明顯不同的顏色,其灰度值差異較大,因此將方解石礦物從整體區(qū)域中區(qū)分出來是相對比較容易的。 本文采用閾值分割法對其進(jìn)行分割處理,將同一閾值分段里的像素歸為同一種礦物,用公式(1)描述[26]:

    式中,G(x)為表征圖像中像素點的值;f(x,y)為原始灰度圖像像素的灰度值;ζn為所選用的閾值。

    常用的閾值分割法主要有灰度直方圖分割法和基于掃描線確定閾值的分割法,灰度直方圖分割法比較適用于圖像灰度直方圖出現(xiàn)多個峰值,閾值在相鄰峰值之間的峰谷選擇[26]。 由于COx 黏土巖圖像的灰度直方圖峰值單一、沒有明顯的峰谷,所以本文采用基于掃描線上確定閾值的分割方法。 如圖3(a)所示,在COx 黏土巖的原始圖像上選擇一條白色掃描線,必須穿過方解石區(qū)和石英與黏土區(qū),并且兩個區(qū)域礦物的顆粒肉眼可以區(qū)分。 白色掃描線上灰度值的變化曲線如圖3(b)所示。 通過比較圖3(a)中掃描線穿過的區(qū)域變化與圖3(b)中曲線的變化情況可以發(fā)現(xiàn),石英和黏土區(qū)的灰度值在紅色線以下,而方解石區(qū)的灰度值在紅色線以上波動。 經(jīng)過多次實驗,最終確定灰度值99 為方解石區(qū)與石英黏土區(qū)的閾值。 因此,在灰度圖像矩陣中,凡是灰度值大于閾值99 的像素都?xì)w為方解石區(qū)。 由此,可以得到方解石的表征圖像,如圖3(c)所示的紅色區(qū)域。 通過對比灰度圖和表征圖,方解石區(qū)的表征圖像能很好地反映其真實細(xì)觀結(jié)構(gòu)。

    圖3 方解石礦物的表征過程Fig.3 Feldspar mineral characterization process

    由此可見,DIP 技術(shù)提供了一種更合理的方法來模擬巖石內(nèi)部復(fù)雜的初始細(xì)觀結(jié)構(gòu)。 經(jīng)過DIP技術(shù)處理后獲取的細(xì)觀結(jié)構(gòu)表征圖像,本質(zhì)上是由像素組成的二維矩陣。 每個像素都有確定的位置信息和顏色信息,用數(shù)字可以表示為(i,j,I)(其中i和j代表像素的位置,I代表像素的顏色屬性)。以COx 黏土巖為例,由于其細(xì)觀結(jié)構(gòu)表征圖像有3 種顏色屬性,本文采用數(shù)字來替換像素的顏色屬性,1 代表黏土,2 代表方解石,3 代表石英。 因此,黏土區(qū)、方解石區(qū)、石英區(qū)像素的數(shù)字形式分別為(i,j,1)、(i,j,2)、(i,j,3)。 這樣,就把200×200 像素的COx 黏土巖細(xì)觀結(jié)構(gòu)表征圖像轉(zhuǎn)化為由數(shù)字1、2、3 構(gòu)成的200×200 的二維數(shù)字矩陣。

    對方解石礦物表征后的圖像如圖3(c)所示,將局部區(qū)域放大可以看出其對應(yīng)的矩陣表,數(shù)字2代表方解石礦物區(qū),字母x 代表還未區(qū)分的黏土(數(shù)字1)和石英(數(shù)字3)。

    1.2.2 石英礦物區(qū)域的表征

    如圖1 所示,石英區(qū)和黏土區(qū)的灰度值較為接近,但是石英和方解石顆粒的尺寸大于黏土顆粒,它們形成了均勻的區(qū)域,該區(qū)域的局部灰度變化不大。 所以,本文采用在局部方差圖像上用分水嶺分割算法對石英區(qū)和黏土區(qū)進(jìn)行分割[25]。

    分水嶺分割算法把灰度圖像看成一幅地形浮雕,亮度強(qiáng)的區(qū)域灰度值較大,亮度暗的區(qū)域灰度值較小,通過尋找匯水盆地和分水嶺界線對圖像進(jìn)行分割。 分水嶺分割方法有時會產(chǎn)生過度分割的現(xiàn)象,為了減少這種情況,需要在分割前對圖像進(jìn)行高斯濾波處理。 首先,計算COx 黏土巖原始圖像的局部方差(在石英黏土區(qū)進(jìn)行,排除方解石區(qū)灰度值的影響),獲取原始圖像的局部方差圖像[圖4(a)]。 從圖4(a)中可以區(qū)分出與石英相關(guān)的均勻區(qū)域,從而可以將黏土與石英分割開。 其次,在濾波后的局部方差圖像上[圖4(b)]進(jìn)行分水嶺分割算法,得到局部方差的分水嶺分割結(jié)果圖[圖4(c)]。 由于局部方差圖像上顯示的石英區(qū)只是真實的一部分[圖4(d)],因此,將局部方差圖像上顯示的石英區(qū)的邊界與分水嶺分割算法結(jié)果圖相結(jié)合,重新識別石英區(qū)。 最后,保留分水嶺分割算法結(jié)果圖中與局部方差圖像中石英區(qū)相對應(yīng)的盆地,將其識別為石英區(qū)[圖4(e)]。 通過與原灰度圖像對比可以看出,經(jīng)過分水嶺算法分割圖像的石英區(qū)能準(zhǔn)確地表征其邊界。

    圖4 石英礦物的表征過程Fig.4 Process of quartz mineral characterization

    1.2.3 黏土礦物區(qū)域的表征

    經(jīng)過上述兩步,將閾值分割法確定的方解石區(qū)[圖3(c)]和分水嶺算法確定的石英區(qū)[圖4(e)]相疊加之后,剩余區(qū)域就是黏土礦物區(qū)。 最終,就可以得到COx 黏土巖整體細(xì)觀結(jié)構(gòu)的表征圖像,如圖5 所示。 與原灰度圖像(圖1)對比可見,通過DIP 技術(shù)得到的巖石細(xì)觀結(jié)構(gòu)表征圖像能較準(zhǔn)確地描述巖石的細(xì)觀結(jié)構(gòu)特征。

    圖5 COx 黏土巖的細(xì)觀結(jié)構(gòu)表征Fig.5 Microscopic characterization of COx argillaceous rocks

    2 數(shù)值方法

    2.1 基于快速傅里葉變換的數(shù)值方法

    基于快速傅里葉變換(FFT)的數(shù)值方法是文獻(xiàn)[27-28]提出的,用于計算復(fù)雜細(xì)觀結(jié)構(gòu)的非均質(zhì)材料力學(xué)行為的數(shù)值計算方法。 基于FFT 的數(shù)值計算方法可以直接利用圖像中的像素點作為材料點,與傳統(tǒng)的數(shù)值模擬方法相比無須再考慮復(fù)雜細(xì)觀結(jié)構(gòu)的邊界。 依據(jù)像素所屬的材料,賦予對應(yīng)材料點相應(yīng)的力學(xué)屬性,以此來將復(fù)雜的細(xì)觀結(jié)構(gòu)特征引入到FFT 模型中。 與FFT 方法相比,傳統(tǒng)的數(shù)值方法在表征巖石復(fù)雜的細(xì)觀結(jié)構(gòu)時,需要生成大量的網(wǎng)格;越復(fù)雜的細(xì)觀結(jié)構(gòu),網(wǎng)格的劃分越困難;大量的網(wǎng)格又影響了數(shù)值模擬的計算效率。 近年來,越來越多的學(xué)者將基于FFT 數(shù)值計算方法應(yīng)用到數(shù)值模擬領(lǐng)域中,FFT 方法得到了發(fā)展和進(jìn)步,并被用來研究孔洞和夾雜的問題[29-30]。

    以非均質(zhì)線彈性材料為例,現(xiàn)介紹基于FFT的數(shù)值計算方法。 材料單元在邊界?V上受到一個均勻的宏觀應(yīng)變?nèi)我庖稽c的局部應(yīng)變ε(x)可由Lippmann-Schwinger 方程求得

    式中,*為卷積符號;τ(x)為線彈性和均勻介質(zhì)中“非均質(zhì)部分”的極化應(yīng)力;σ(x)為局部應(yīng)力;C0為均勻介質(zhì)的彈性張量,與材料的局部位置x無關(guān)。

    格林算子Γ0在傅里葉空間中可以表示為

    式中,ξ為傅里葉空間中的離散頻率;λ0和μ0為材料介質(zhì)的拉梅系數(shù);δ為Kronecker's delta。

    由公式(2)求得代表單元內(nèi)任意一點的局部應(yīng)變ε(x)后,即可通過代表單元體內(nèi)的彈性張量C(x)求得對應(yīng)的局部應(yīng)力σ(x)。C(x)嚴(yán)格依賴于局部坐標(biāo)位置x。 材料的宏觀應(yīng)力可由整個單元局部應(yīng)力的體積平均值計算:

    針對非線性問題,本構(gòu)關(guān)系可以采用增量形式表示,即將總的加載過程劃分為有限的N步,從初始條件到n步的局部應(yīng)力σn,局部應(yīng)變以及內(nèi)部變量Vn都是已知的。 對于宏觀應(yīng)變增量在n+1 步時,局部非線性本構(gòu)關(guān)系可以用式(5)增量形式表示:

    式中,εp(x)為非線性局部應(yīng)變增量,如塑性應(yīng)變;Ctan為與加載歷史和加載方向有關(guān)的切線4 階張量。

    2.2 基于真實細(xì)觀結(jié)構(gòu)的FFT 數(shù)值方法的建立

    FFT 方法的主要優(yōu)勢在于避免了對非規(guī)則復(fù)雜細(xì)觀幾何形狀的網(wǎng)格劃分,取而代之的是一系列離散的材料點來模擬非均質(zhì)材料的細(xì)觀結(jié)構(gòu)特征。通過材料點的位置信息和材料屬性的對應(yīng)關(guān)系即可準(zhǔn)確獲得非均質(zhì)材料的細(xì)觀結(jié)構(gòu)特征和與其對應(yīng)的力學(xué)性質(zhì)。 因此,可以直接將經(jīng)過圖像處理的巖石真實細(xì)觀結(jié)構(gòu)表征的圖像像素信息(數(shù)字矩陣)作為FFT 方法的材料點,導(dǎo)入到FFT 程序,并根據(jù)數(shù)字所屬礦物的種類賦予相應(yīng)的細(xì)觀力學(xué)參數(shù),就可以進(jìn)行計算分析。

    如圖6 所示,首先獲取巖石的數(shù)字圖像[圖6(a)],然后提取圖像的像素信息得到每一個像素的位置信息和顏色信息[如圖6(b)],最后進(jìn)行圖像分割得到細(xì)觀結(jié)構(gòu)表征圖像[圖6(c)]。 圖6(c)中的像素值就是FFT 方法中的材料點,可以直接將其導(dǎo)入FFT 算法建立幾何模型,然后根據(jù)像素信息確定每個像素點的力學(xué)性質(zhì)。 因此,用于計算分析的FFT 模型與COx 黏土巖的細(xì)觀結(jié)構(gòu)表征圖像完全一致,真實完整地反映了巖石內(nèi)部的細(xì)觀結(jié)構(gòu)。 綜上所述,基于真實細(xì)觀結(jié)構(gòu)的FFT 數(shù)值計算方法的流程如圖7 所示。

    圖6 材料點生成過程[18]Fig.6 Material point generation process

    圖7 基于真實細(xì)觀結(jié)構(gòu)的FFT 數(shù)值計算方法流程Fig.7 Flow chart of the actual-mesostructure-based FFT numerical method

    基于真實細(xì)觀結(jié)構(gòu)的FFT 模型直接利用了像素原有的信息,不需要進(jìn)行任何處理,考慮了巖石的非均質(zhì)性。 同時,對于同一形狀的物體,像素越高,細(xì)觀結(jié)構(gòu)的邊界越接近真實的邊界(圖8)。 因此,基于真實細(xì)觀結(jié)構(gòu)的FFT 數(shù)值方法在模擬巖石細(xì)觀結(jié)構(gòu)上具有理想的效率和準(zhǔn)確性。

    圖8 像素對邊界的影響[31]Fig.8 The effect of the pixel on the boundary

    2.3 細(xì)觀本構(gòu)模型

    在分析COx 黏土巖礦物成分和細(xì)觀結(jié)構(gòu)的基礎(chǔ)上,依據(jù)黏土基質(zhì)、方解石和石英礦物的力學(xué)性能,選擇合理的局部本構(gòu)模型。 研究發(fā)現(xiàn),黏土礦物的拉伸和抗壓強(qiáng)度相對較低,在荷載作用下容易產(chǎn)生較大的變形。 因此,黏土礦物可以選用彈塑性本構(gòu)模型。 與黏土礦物相比,方解石和石英的礦物顆粒尺寸遠(yuǎn)大于黏土,相對堅硬,不易變形,可以采用線彈性模型[18]。 為了更準(zhǔn)確地描述本研究所使用的局部本構(gòu)模型,圖9 給出各種礦物的局部本構(gòu)曲線示意圖。

    圖9 細(xì)觀結(jié)構(gòu)和局部本構(gòu)模型的關(guān)系[12]Fig.9 Relationship between microstructure and local constitutive model

    石英和方解石采用線彈性模型,只需要確定彈性模量和泊松比即可。 而根據(jù)黏土礦物的力學(xué)性質(zhì),本文采用非關(guān)聯(lián)Drucker-Prager 準(zhǔn)則來描述其非線性力學(xué)行為

    式中,σp為平均應(yīng)力;σq為等效應(yīng)力;S為偏應(yīng)力張量;Cp為材料內(nèi)聚系數(shù);αp(γp)為摩擦系數(shù);γp為等效塑性剪切應(yīng)變,為αp的內(nèi)變量。

    本文采用了各向同性塑性硬化規(guī)律,用公式表示為

    式中,為摩擦系數(shù)的初始值;為摩擦系數(shù)的最終值;b為控制硬化速率的系數(shù)。

    為了更好地描述COx 黏土巖的塑性體積應(yīng)變,采用非關(guān)聯(lián)流動準(zhǔn)則,塑性勢能函數(shù)表示為

    式中,βp為塑性體積膨脹系數(shù)為塑性體積膨脹系數(shù)的初始值;為塑性體積膨脹系數(shù)的最終值;b1為控制塑性體積膨脹速率的系數(shù)。

    針對非線性問題,局部本構(gòu)關(guān)系以增量的形式可表示為

    式中,Cm為基質(zhì)的彈性剛度張量;Δεp為塑性應(yīng)變增量;Δλp為塑性乘子;dev 為偏量。

    Δεp和Δλp的值受當(dāng)前塑性狀態(tài)和加載/卸載準(zhǔn)則的影響,最終由塑性一致性條件確定。

    3 數(shù)值驗證與分析

    3.1 參數(shù)標(biāo)定與宏觀力學(xué)行為的預(yù)測

    將本文提出的基于真實細(xì)觀結(jié)構(gòu)的FFT 數(shù)值方法用于模擬COx 黏土巖的力學(xué)行為。 首先對模型參數(shù)進(jìn)行標(biāo)定。 方解石和石英彈性參數(shù)采用文獻(xiàn)[32]中測得的實驗值。 參照COx 黏土巖、方解石和石英的宏觀彈性參數(shù)的實驗值,采用Mori-Tanaka 方法反向求取了黏土基質(zhì)的彈性參數(shù)[16]。 FFT 模型彈性參數(shù)標(biāo)定的結(jié)果見表1。 文獻(xiàn)[32]對不同深度的COx 黏土巖進(jìn)行成分測量,不同深度所含礦物的含量是不同的。 為了提高數(shù)值模擬的準(zhǔn)確性,運(yùn)用數(shù)字圖像處理技術(shù)對每個深度巖石的表面圖像進(jìn)行了表征,將巖石的細(xì)觀結(jié)構(gòu)納入數(shù)值模型中。 統(tǒng)計了不同礦物所對應(yīng)的像素信息,最終得到FFT 模型中礦物含量與實驗測得含量的對比結(jié)果(表2)。

    表1 COx 黏土巖細(xì)觀介質(zhì)的彈性參數(shù)[32]Table 1 Elastic parameters of meso-medium of COx clay rock

    表2 不同深度各種礦物的含量[32]Table 2 The content of various minerals at different depths

    利用深度為482.2 m 的COx 黏土巖樣品,在圍壓為10 MPa 下的三軸壓縮實驗來標(biāo)定黏土基質(zhì)的塑性參數(shù),如圖10 所示。 將數(shù)值模擬曲線與實驗曲線擬合,取最優(yōu)的一組參數(shù),最終的非彈性參數(shù)結(jié)果如下-3.9,0.1,b=170,b1=300。

    圖10 參數(shù)標(biāo)定Fig.10 Parameter calibration

    文獻(xiàn)[32]對不同深度下的COx 黏土巖進(jìn)行了三軸壓縮實驗,本文對深度為451.4 m、469.0 m、469.1 m、482.2 m 和451.5 m 的樣品進(jìn)行了數(shù)值模擬實驗,并且預(yù)測了不同深度、不同圍壓下的宏觀力學(xué)行為。 首先,依照數(shù)字圖像處理的結(jié)果,將COx 黏土巖的細(xì)觀結(jié)構(gòu)納入數(shù)值模型中,并按照對應(yīng)的本構(gòu)關(guān)系和標(biāo)定的參數(shù)進(jìn)行計算,最終的對比結(jié)果如圖11 所示。

    圖11 數(shù)值模擬與三軸壓縮實驗數(shù)據(jù)的對比Fig.11 Comparison of numerical simulation and experimental data

    從數(shù)值模擬和實驗數(shù)據(jù)對比結(jié)果可以看出,本文提出的基于真實細(xì)觀結(jié)構(gòu)的FFT 數(shù)值計算方法能根據(jù)巖石不同細(xì)觀結(jié)構(gòu)預(yù)測其宏觀力學(xué)行為,不同深度、不同圍壓下COx 黏土巖的塑性階段軸線和側(cè)向曲線均擬合較好,峰值強(qiáng)度的預(yù)測準(zhǔn)確。FFT 模型真實地還原了巖石內(nèi)部的細(xì)觀結(jié)構(gòu),并且考慮了不同礦物的分布和相互作用。 但是模擬實驗中沒有考慮巖石的損傷階段,沒有考慮微裂紋、孔洞等缺陷對巖石宏觀力學(xué)行為的影響,這些問題將會在今后的研究工作中進(jìn)一步解決。

    3.2 細(xì)觀結(jié)構(gòu)特征與局部應(yīng)力分布

    為了進(jìn)一步揭示巖石細(xì)觀結(jié)構(gòu)特征與局部應(yīng)力的影響,取上述4 組數(shù)值模擬的細(xì)觀結(jié)構(gòu)圖及對應(yīng)的應(yīng)力云圖(圖12)。

    由圖12 可以看出,不同種類礦物的分布和形狀直接影響局部應(yīng)力場的分布。 應(yīng)力較大的區(qū)域(淺綠色),主要集中在不同礦物的邊界和細(xì)小的方解石區(qū)。 這是因為部分方解石礦物夾雜的形狀較為細(xì)小且不均勻,受其非均質(zhì)性的影響,這些區(qū)域容易產(chǎn)生局部應(yīng)力集中現(xiàn)象,使得該區(qū)域的應(yīng)力較大。 同時,在不同礦物的邊界也容易產(chǎn)生應(yīng)力集中現(xiàn)象。 石英區(qū)的應(yīng)力除了一些邊界會出現(xiàn)應(yīng)力集中現(xiàn)象,其他區(qū)域的應(yīng)力分布相對均勻(淺黃色),這是因為石英礦物夾雜是成塊存在的,并且這些塊狀區(qū)域內(nèi)不包含其他礦物,相對均質(zhì)。 黏土礦物的區(qū)域最大,應(yīng)力最小(深黃色區(qū)域),除去黏土內(nèi)部包含的其他礦物,黏土的應(yīng)力相對均勻分布于應(yīng)力云圖中。

    圖12 細(xì)觀結(jié)構(gòu)特征與局部應(yīng)力Fig.12 Microstructure and local stress

    綜上分析,非均質(zhì)巖石類材料的細(xì)觀結(jié)構(gòu)的分布和幾何形狀對其局部的應(yīng)力分布有重要的影響。 本研究提供的細(xì)觀結(jié)構(gòu)局部應(yīng)力的演化過程,可以合理地解釋細(xì)觀結(jié)構(gòu)特征對其應(yīng)力分布的影響規(guī)律。

    3.3 細(xì)觀結(jié)構(gòu)與宏觀力學(xué)行為的關(guān)系

    為了進(jìn)一步探索COx 黏土巖細(xì)觀結(jié)構(gòu)對其宏觀力學(xué)行為的影響,本文以深度482.2 m 的巖石樣品為例,利用基于真實細(xì)觀結(jié)構(gòu)的FFT 模型,給出了該樣品在10 MPa 三軸壓縮模擬試驗的應(yīng)力-應(yīng)變曲線與細(xì)觀應(yīng)力云圖之間的關(guān)系,如圖13 所示。

    圖13 細(xì)觀結(jié)構(gòu)應(yīng)力演化對宏觀力學(xué)行為的影響Fig.13 Influence of microstructure stress evolution on the macroscopic mechanical behavior

    從B 點應(yīng)力云圖可以看出,黏土基質(zhì)對應(yīng)的應(yīng)力值最小,方解石對應(yīng)的應(yīng)力值最大,石英的居中。 較大的應(yīng)力值主要分布在方解石和黏土基質(zhì)的分界線上,這是由于方解石分布不均勻,有較多的細(xì)小區(qū)域?qū)е戮植繎?yīng)力集中現(xiàn)象。 石英區(qū)相對均勻,應(yīng)力分布相對均勻。 A-C 階段云圖變化比較迅速,大面積的黏土基質(zhì)區(qū)由深黃色變?yōu)闇\黃色(淺黃色代表應(yīng)力較大的區(qū)域),每種礦物上升的幅度都很大,對應(yīng)的整體宏觀應(yīng)力-應(yīng)變曲線上升較迅速,應(yīng)力變化范圍較大。 當(dāng)應(yīng)力增加到圖中的C 點時,黏土基質(zhì)的應(yīng)力繼續(xù)增大,但是因為黏土基質(zhì)采用的是彈塑性本構(gòu)模型,此時大部分黏土礦物開始進(jìn)入塑性階段,應(yīng)力變化范圍較小,石英和方解石區(qū)的應(yīng)力繼續(xù)上升。 當(dāng)曲線到達(dá)D 點時,黏土基質(zhì)大部分區(qū)域的應(yīng)力已經(jīng)進(jìn)入屈服,只有少數(shù)區(qū)域的應(yīng)力還在增加,曲線應(yīng)力的上升主要由石英和方解石控制。 石英和方解石區(qū)域的應(yīng)力繼續(xù)增加,石英區(qū)域增加較均勻,方解石受分布不均勻的影響出現(xiàn)了一些應(yīng)力較大的區(qū)域,此階段宏觀曲線上升趨勢進(jìn)一步減緩。

    可見,本文提出的數(shù)值方法可以從細(xì)觀應(yīng)力云圖來可視化宏觀力學(xué)行為的變化過程,合理解釋細(xì)觀結(jié)構(gòu)特征與宏觀力學(xué)行為之間的內(nèi)在聯(lián)系。

    4 結(jié) 論

    利用FFT 數(shù)值方法可以直接獲取圖像像素信息的特點,與數(shù)值圖像處理技術(shù)有機(jī)結(jié)合,本文提出了一種考慮巖石內(nèi)部真實細(xì)觀結(jié)構(gòu)特征的FFT數(shù)值分析方法。 在此基礎(chǔ)上,以COx 黏土巖為例,研究了巖石內(nèi)部真實細(xì)觀結(jié)構(gòu)對其宏觀非線性力學(xué)行為的影響。 主要結(jié)論如下:

    (1) 針對巖石具有復(fù)雜細(xì)觀結(jié)構(gòu)特征的特點,提出了基于真實細(xì)觀結(jié)構(gòu)的FFT 數(shù)值方法。 該方法無須劃分網(wǎng)格,在研究巖石類材料時具有一定的優(yōu)勢。

    (2) 實現(xiàn)了非均質(zhì)性材料的力學(xué)分析,合理預(yù)測了不同賦存深度、不同圍壓下黏土巖的峰前非線性行為和峰值強(qiáng)度,以及細(xì)觀結(jié)構(gòu)特征對其宏觀非線性力學(xué)行為的影響。

    (3) 從細(xì)觀尺度上揭示了礦物的含量、形狀、大小和分布等特征對宏觀力學(xué)行為的影響機(jī)理。 數(shù)值計算結(jié)果表明,巖石材料內(nèi)部細(xì)觀結(jié)構(gòu)決定了材料的應(yīng)力分布,顯著影響材料的宏觀非線性行為。

    猜你喜歡
    細(xì)觀方解石石英
    F-在方解石表面的吸附及其對方解石表面性質(zhì)的影響
    硅酸鹽通報(2022年8期)2022-09-08 04:25:42
    石英及方石英的堿浸動力學(xué)研究
    氯化鈣和碳酸鈉對方解石浮選的影響及其機(jī)理研究
    基于細(xì)觀結(jié)構(gòu)的原狀黃土動彈性模量和阻尼比試驗研究
    地震研究(2021年1期)2021-04-13 01:05:24
    貴州重晶石與方解石常溫浮選分離試驗研究
    關(guān)于寬帶石英濾波器配套諧振器的選用
    電子制作(2018年14期)2018-08-21 01:38:10
    嚴(yán)重楔狀缺損石英纖維樁樹脂修復(fù)臨床分析
    石英云母片巖力學(xué)性質(zhì)各向異性的模擬方法探討
    基于四叉樹網(wǎng)格加密技術(shù)的混凝土細(xì)觀模型
    PBX炸藥的抗壓強(qiáng)度及抗拉強(qiáng)度細(xì)觀尺度的數(shù)值計算
    亚洲天堂国产精品一区在线| 久久久水蜜桃国产精品网| 亚洲国产精品成人综合色| 欧美日韩福利视频一区二区| 国产av一区在线观看免费| 黑人操中国人逼视频| 999精品在线视频| 国产精品一及| 18美女黄网站色大片免费观看| 可以免费在线观看a视频的电影网站| 亚洲av中文字字幕乱码综合| 91字幕亚洲| 午夜免费观看网址| 18禁裸乳无遮挡免费网站照片| 欧美日韩亚洲国产一区二区在线观看| 欧美性长视频在线观看| 女人高潮潮喷娇喘18禁视频| 琪琪午夜伦伦电影理论片6080| 一个人观看的视频www高清免费观看 | 亚洲va日本ⅴa欧美va伊人久久| 久久欧美精品欧美久久欧美| 亚洲国产欧美人成| 最近在线观看免费完整版| 国产av一区在线观看免费| 白带黄色成豆腐渣| 国内揄拍国产精品人妻在线| 国产高清有码在线观看视频 | 亚洲 欧美 日韩 在线 免费| 美女黄网站色视频| 久久性视频一级片| 日韩欧美国产在线观看| 精品国产亚洲在线| 长腿黑丝高跟| 精品熟女少妇八av免费久了| 天堂av国产一区二区熟女人妻 | 99riav亚洲国产免费| 色尼玛亚洲综合影院| 亚洲国产精品合色在线| 色在线成人网| 亚洲精品美女久久av网站| 一本一本综合久久| 国产在线精品亚洲第一网站| 三级国产精品欧美在线观看 | 制服诱惑二区| 欧美最黄视频在线播放免费| 日韩精品青青久久久久久| 亚洲国产欧美人成| 免费在线观看亚洲国产| 亚洲va日本ⅴa欧美va伊人久久| 亚洲自拍偷在线| 俄罗斯特黄特色一大片| 亚洲欧美精品综合久久99| 欧美乱色亚洲激情| 18禁黄网站禁片免费观看直播| 香蕉丝袜av| 三级毛片av免费| 麻豆一二三区av精品| 久久久久久九九精品二区国产 | 我的老师免费观看完整版| 少妇被粗大的猛进出69影院| 日韩大码丰满熟妇| 日韩大尺度精品在线看网址| 亚洲av成人不卡在线观看播放网| 舔av片在线| 在线看三级毛片| 亚洲av成人一区二区三| 国产精品久久久av美女十八| 啪啪无遮挡十八禁网站| 国产又黄又爽又无遮挡在线| 黄色 视频免费看| 老汉色av国产亚洲站长工具| 三级毛片av免费| 国产亚洲精品综合一区在线观看 | aaaaa片日本免费| 亚洲中文日韩欧美视频| 欧美又色又爽又黄视频| 俄罗斯特黄特色一大片| 国产成人啪精品午夜网站| 手机成人av网站| 国产主播在线观看一区二区| 1024视频免费在线观看| 日韩欧美国产在线观看| 午夜视频精品福利| 最近视频中文字幕2019在线8| 免费在线观看成人毛片| 国产成+人综合+亚洲专区| 日韩欧美精品v在线| 午夜亚洲福利在线播放| 亚洲第一电影网av| 老汉色∧v一级毛片| 亚洲熟妇熟女久久| 亚洲成a人片在线一区二区| 中文资源天堂在线| 18禁黄网站禁片免费观看直播| 亚洲av熟女| 法律面前人人平等表现在哪些方面| 国产乱人伦免费视频| 日韩精品免费视频一区二区三区| 给我免费播放毛片高清在线观看| 亚洲av中文字字幕乱码综合| 亚洲成人中文字幕在线播放| 日韩大码丰满熟妇| 中文字幕av在线有码专区| 久久欧美精品欧美久久欧美| 国产成人精品久久二区二区91| 精品久久久久久成人av| 久久精品91蜜桃| 一二三四社区在线视频社区8| 最近在线观看免费完整版| 一区二区三区国产精品乱码| 国产又色又爽无遮挡免费看| 制服诱惑二区| 男男h啪啪无遮挡| 亚洲真实伦在线观看| 一个人免费在线观看电影 | 亚洲va日本ⅴa欧美va伊人久久| 色噜噜av男人的天堂激情| 中亚洲国语对白在线视频| 日韩欧美国产在线观看| 男人舔奶头视频| 免费看美女性在线毛片视频| 国产片内射在线| 日本黄色视频三级网站网址| 欧美性猛交╳xxx乱大交人| 可以免费在线观看a视频的电影网站| 欧美日韩精品网址| 亚洲熟妇熟女久久| 成人三级做爰电影| 精品一区二区三区视频在线观看免费| 亚洲av电影在线进入| 国产97色在线日韩免费| 真人做人爱边吃奶动态| 国内精品一区二区在线观看| 成人国产一区最新在线观看| 天堂√8在线中文| 无人区码免费观看不卡| 国产精品国产高清国产av| 国产免费av片在线观看野外av| 两性夫妻黄色片| 日日摸夜夜添夜夜添小说| 国产三级在线视频| 这个男人来自地球电影免费观看| 国产成人av激情在线播放| 午夜福利高清视频| 最近最新中文字幕大全免费视频| 国产aⅴ精品一区二区三区波| 久久婷婷成人综合色麻豆| 90打野战视频偷拍视频| www.精华液| 成人特级黄色片久久久久久久| 免费在线观看完整版高清| 女警被强在线播放| 曰老女人黄片| 国产午夜精品久久久久久| 97超级碰碰碰精品色视频在线观看| 超碰成人久久| 丰满人妻一区二区三区视频av | 黄色视频不卡| 精品久久蜜臀av无| 每晚都被弄得嗷嗷叫到高潮| 亚洲一码二码三码区别大吗| 中文字幕久久专区| 五月伊人婷婷丁香| 99在线视频只有这里精品首页| 在线国产一区二区在线| 久久99热这里只有精品18| 精品欧美国产一区二区三| 亚洲av片天天在线观看| 精品国产美女av久久久久小说| 1024手机看黄色片| 日韩国内少妇激情av| 国产91精品成人一区二区三区| 女警被强在线播放| 三级男女做爰猛烈吃奶摸视频| 国产av不卡久久| 国产区一区二久久| 亚洲人成网站在线播放欧美日韩| 极品教师在线免费播放| 99热这里只有是精品50| 美女大奶头视频| 国产精华一区二区三区| 欧美 亚洲 国产 日韩一| 免费av毛片视频| 国产成人av教育| 久久久久国产精品人妻aⅴ院| 久久天堂一区二区三区四区| 国产单亲对白刺激| 一进一出抽搐gif免费好疼| 久久人妻福利社区极品人妻图片| 久久久精品大字幕| 亚洲国产看品久久| 亚洲人与动物交配视频| 在线十欧美十亚洲十日本专区| 久久中文看片网| 亚洲人与动物交配视频| 国产亚洲精品久久久久久毛片| 国产亚洲欧美在线一区二区| 精品久久久久久久久久久久久| 18禁观看日本| 久久天躁狠狠躁夜夜2o2o| 麻豆一二三区av精品| 久久精品夜夜夜夜夜久久蜜豆 | 一夜夜www| xxx96com| 日韩精品免费视频一区二区三区| 精品一区二区三区四区五区乱码| 在线观看www视频免费| 欧美日韩亚洲综合一区二区三区_| 高清在线国产一区| 97人妻精品一区二区三区麻豆| 俺也久久电影网| 亚洲国产日韩欧美精品在线观看 | 成人三级黄色视频| av福利片在线| 精品不卡国产一区二区三区| 51午夜福利影视在线观看| 亚洲精品国产一区二区精华液| 色哟哟哟哟哟哟| 18禁美女被吸乳视频| 日本一本二区三区精品| 中文字幕人成人乱码亚洲影| 狂野欧美白嫩少妇大欣赏| 午夜激情av网站| 国产精品日韩av在线免费观看| 熟女少妇亚洲综合色aaa.| 成人av一区二区三区在线看| 日韩大码丰满熟妇| 欧美不卡视频在线免费观看 | 岛国在线观看网站| 午夜福利免费观看在线| 国产aⅴ精品一区二区三区波| 国产野战对白在线观看| 亚洲成人久久爱视频| 免费无遮挡裸体视频| 久久香蕉国产精品| 亚洲在线自拍视频| 午夜a级毛片| bbb黄色大片| 久热爱精品视频在线9| 舔av片在线| 精品一区二区三区av网在线观看| 亚洲精品久久成人aⅴ小说| 精品一区二区三区四区五区乱码| 麻豆久久精品国产亚洲av| 欧美zozozo另类| 一级毛片高清免费大全| 波多野结衣高清作品| 色av中文字幕| 午夜福利18| 给我免费播放毛片高清在线观看| av福利片在线观看| 俄罗斯特黄特色一大片| 精品高清国产在线一区| www.999成人在线观看| 欧美成人免费av一区二区三区| 老司机靠b影院| 夜夜爽天天搞| 欧美日韩国产亚洲二区| 久久久国产精品麻豆| 成人18禁高潮啪啪吃奶动态图| 97人妻精品一区二区三区麻豆| 欧美日韩福利视频一区二区| 三级国产精品欧美在线观看 | 精品久久久久久久久久久久久| 亚洲美女黄片视频| 中文字幕精品亚洲无线码一区| 国产视频一区二区在线看| 俄罗斯特黄特色一大片| 国产成人精品久久二区二区免费| 国产欧美日韩精品亚洲av| 国产精品久久电影中文字幕| 欧美在线黄色| 免费看日本二区| 大型黄色视频在线免费观看| 制服诱惑二区| 日本黄色视频三级网站网址| 在线观看午夜福利视频| 国产精品98久久久久久宅男小说| 大型av网站在线播放| 亚洲av片天天在线观看| 亚洲国产欧美人成| 脱女人内裤的视频| 99精品欧美一区二区三区四区| 人妻久久中文字幕网| 欧美又色又爽又黄视频| 黄色片一级片一级黄色片| 国产精品久久久久久久电影 | 亚洲无线在线观看| 99热这里只有是精品50| 18禁国产床啪视频网站| 久久久久国产一级毛片高清牌| 99久久久亚洲精品蜜臀av| 精品国产乱码久久久久久男人| 高潮久久久久久久久久久不卡| 草草在线视频免费看| 99久久精品国产亚洲精品| av福利片在线观看| 国产午夜精品久久久久久| 国产aⅴ精品一区二区三区波| 日韩 欧美 亚洲 中文字幕| av片东京热男人的天堂| 亚洲真实伦在线观看| 欧美日韩一级在线毛片| 国产爱豆传媒在线观看 | 99国产极品粉嫩在线观看| 国产欧美日韩一区二区精品| 老汉色∧v一级毛片| 日韩高清综合在线| 日本成人三级电影网站| 日本 欧美在线| 亚洲欧美日韩东京热| av在线天堂中文字幕| 怎么达到女性高潮| 国产私拍福利视频在线观看| 日本三级黄在线观看| 亚洲国产精品999在线| 黄色片一级片一级黄色片| 三级国产精品欧美在线观看 | 国产成人av激情在线播放| 少妇被粗大的猛进出69影院| 国产av一区在线观看免费| 国产不卡一卡二| 性欧美人与动物交配| 男人舔奶头视频| 国产真实乱freesex| 亚洲精品在线美女| 禁无遮挡网站| 久久久久久免费高清国产稀缺| 欧美黑人巨大hd| 天天一区二区日本电影三级| 18禁美女被吸乳视频| 国产亚洲精品一区二区www| www.999成人在线观看| 最近最新中文字幕大全电影3| 午夜a级毛片| 特大巨黑吊av在线直播| 18禁黄网站禁片免费观看直播| 亚洲免费av在线视频| 波多野结衣高清无吗| 一级作爱视频免费观看| 亚洲全国av大片| 在线观看日韩欧美| 久久精品人妻少妇| 99久久国产精品久久久| 人人妻,人人澡人人爽秒播| 操出白浆在线播放| 伊人久久大香线蕉亚洲五| 久99久视频精品免费| 悠悠久久av| 制服诱惑二区| 欧美一区二区精品小视频在线| 麻豆国产av国片精品| 亚洲 欧美 日韩 在线 免费| 久久午夜亚洲精品久久| 每晚都被弄得嗷嗷叫到高潮| 亚洲人与动物交配视频| 精品人妻1区二区| 老司机午夜十八禁免费视频| 亚洲自拍偷在线| 久久久久久亚洲精品国产蜜桃av| 搡老岳熟女国产| 中出人妻视频一区二区| 日韩精品免费视频一区二区三区| 亚洲第一电影网av| 精品乱码久久久久久99久播| 淫妇啪啪啪对白视频| 亚洲成人国产一区在线观看| 精品一区二区三区av网在线观看| 床上黄色一级片| 日日爽夜夜爽网站| 欧美性猛交黑人性爽| 99国产精品99久久久久| 国产精品自产拍在线观看55亚洲| 一本大道久久a久久精品| a级毛片a级免费在线| 看黄色毛片网站| 看黄色毛片网站| 麻豆国产97在线/欧美 | 深夜精品福利| 国产精品电影一区二区三区| 露出奶头的视频| 在线国产一区二区在线| 在线观看美女被高潮喷水网站 | 黄片小视频在线播放| 午夜久久久久精精品| 久久天堂一区二区三区四区| 欧美一级a爱片免费观看看 | 色老头精品视频在线观看| 欧美在线一区亚洲| 老司机在亚洲福利影院| 99riav亚洲国产免费| 亚洲电影在线观看av| 国产精华一区二区三区| 午夜福利高清视频| 夜夜夜夜夜久久久久| 亚洲欧美日韩无卡精品| 国产69精品久久久久777片 | 99久久无色码亚洲精品果冻| 午夜两性在线视频| 久久九九热精品免费| 777久久人妻少妇嫩草av网站| 在线观看美女被高潮喷水网站 | 动漫黄色视频在线观看| 国产爱豆传媒在线观看 | 九色国产91popny在线| 男女午夜视频在线观看| 久久九九热精品免费| 国产精品亚洲美女久久久| 国产精品,欧美在线| 中文字幕精品亚洲无线码一区| 日本熟妇午夜| 免费电影在线观看免费观看| 亚洲第一欧美日韩一区二区三区| 天天躁夜夜躁狠狠躁躁| 免费看十八禁软件| 中文字幕精品亚洲无线码一区| 99在线视频只有这里精品首页| 国产av麻豆久久久久久久| 亚洲美女视频黄频| 欧美日韩一级在线毛片| 久热爱精品视频在线9| 亚洲av熟女| 国内精品久久久久久久电影| 好男人电影高清在线观看| 后天国语完整版免费观看| 全区人妻精品视频| 母亲3免费完整高清在线观看| 99精品在免费线老司机午夜| 亚洲男人天堂网一区| 久久午夜综合久久蜜桃| 中文字幕人成人乱码亚洲影| av片东京热男人的天堂| 最新在线观看一区二区三区| 日韩欧美 国产精品| 欧美成人一区二区免费高清观看 | 国产伦一二天堂av在线观看| 日韩三级视频一区二区三区| xxx96com| 在线永久观看黄色视频| www国产在线视频色| 人人妻,人人澡人人爽秒播| 欧美日韩一级在线毛片| 亚洲中文字幕一区二区三区有码在线看 | 亚洲激情在线av| 又爽又黄无遮挡网站| 国产精品久久久久久亚洲av鲁大| 日本 欧美在线| 欧美日韩瑟瑟在线播放| 舔av片在线| 丁香六月欧美| 免费在线观看完整版高清| 欧美成人免费av一区二区三区| 欧美激情久久久久久爽电影| 成人av一区二区三区在线看| 欧美在线一区亚洲| 久久久久免费精品人妻一区二区| 哪里可以看免费的av片| 国产片内射在线| avwww免费| 精品欧美一区二区三区在线| 久久久久亚洲av毛片大全| 99精品久久久久人妻精品| 熟女电影av网| 久久天躁狠狠躁夜夜2o2o| 又黄又爽又免费观看的视频| 免费观看人在逋| 91在线观看av| 久久久国产精品麻豆| 亚洲精品av麻豆狂野| www.999成人在线观看| 国产在线精品亚洲第一网站| 91麻豆精品激情在线观看国产| 欧美另类亚洲清纯唯美| 一级a爱片免费观看的视频| 女人被狂操c到高潮| 婷婷精品国产亚洲av在线| 亚洲成人久久爱视频| 久久婷婷成人综合色麻豆| 久久人人精品亚洲av| 欧美日韩中文字幕国产精品一区二区三区| 国产在线精品亚洲第一网站| 亚洲成人国产一区在线观看| 伦理电影免费视频| 亚洲,欧美精品.| 欧美一级a爱片免费观看看 | 婷婷丁香在线五月| 国产精品久久久久久亚洲av鲁大| 亚洲欧美日韩高清专用| 变态另类成人亚洲欧美熟女| 听说在线观看完整版免费高清| 日日摸夜夜添夜夜添小说| 变态另类丝袜制服| 国产在线精品亚洲第一网站| 免费在线观看完整版高清| 91大片在线观看| 亚洲激情在线av| 午夜久久久久精精品| 久久天堂一区二区三区四区| 亚洲第一电影网av| 久久精品国产亚洲av香蕉五月| 午夜福利成人在线免费观看| 国产真人三级小视频在线观看| 精品欧美一区二区三区在线| 美女午夜性视频免费| 老司机午夜福利在线观看视频| 国产在线观看jvid| 免费看日本二区| 亚洲国产精品999在线| 99久久综合精品五月天人人| 亚洲成人中文字幕在线播放| 成人三级做爰电影| 午夜福利欧美成人| 久久久久国产一级毛片高清牌| 午夜福利18| 日本 欧美在线| 好男人在线观看高清免费视频| 亚洲电影在线观看av| 两性夫妻黄色片| 成人av一区二区三区在线看| 制服诱惑二区| 午夜激情av网站| 亚洲aⅴ乱码一区二区在线播放 | 91麻豆精品激情在线观看国产| 特大巨黑吊av在线直播| 啦啦啦免费观看视频1| 可以免费在线观看a视频的电影网站| 97人妻精品一区二区三区麻豆| 免费看日本二区| 非洲黑人性xxxx精品又粗又长| 久久婷婷成人综合色麻豆| 欧美人与性动交α欧美精品济南到| 999精品在线视频| 日韩有码中文字幕| 97碰自拍视频| e午夜精品久久久久久久| 亚洲精品一区av在线观看| 欧美丝袜亚洲另类 | 亚洲中文字幕日韩| 亚洲五月婷婷丁香| 亚洲中文字幕一区二区三区有码在线看 | 日韩精品免费视频一区二区三区| 熟女少妇亚洲综合色aaa.| 黄色女人牲交| 欧美大码av| 一二三四社区在线视频社区8| 窝窝影院91人妻| 视频区欧美日本亚洲| netflix在线观看网站| 午夜亚洲福利在线播放| 婷婷亚洲欧美| 99riav亚洲国产免费| 色哟哟哟哟哟哟| 美女 人体艺术 gogo| 午夜福利在线在线| 成年女人毛片免费观看观看9| 一区福利在线观看| 一本精品99久久精品77| 十八禁网站免费在线| 一级毛片精品| 少妇裸体淫交视频免费看高清 | 宅男免费午夜| 久久 成人 亚洲| 久久天躁狠狠躁夜夜2o2o| 日本五十路高清| 最近最新中文字幕大全电影3| 一个人免费在线观看电影 | 国产三级黄色录像| 变态另类成人亚洲欧美熟女| 最新在线观看一区二区三区| 波多野结衣高清作品| 亚洲午夜精品一区,二区,三区| 99久久国产精品久久久| 1024香蕉在线观看| 国产伦在线观看视频一区| 国产av又大| 日韩 欧美 亚洲 中文字幕| 麻豆一二三区av精品| 在线观看午夜福利视频| 日韩成人在线观看一区二区三区| 免费在线观看日本一区| 亚洲精华国产精华精| 国产三级在线视频| 怎么达到女性高潮| 母亲3免费完整高清在线观看| 国产av在哪里看| 国产成人av激情在线播放| 欧美在线一区亚洲| 人人妻,人人澡人人爽秒播| 精品久久久久久久久久免费视频| 国产成人精品久久二区二区91| 又粗又爽又猛毛片免费看| 最新美女视频免费是黄的| 亚洲av第一区精品v没综合| 天天躁夜夜躁狠狠躁躁| 国产又黄又爽又无遮挡在线| 日日干狠狠操夜夜爽| 亚洲专区中文字幕在线| 丰满人妻熟妇乱又伦精品不卡| 精品国产乱码久久久久久男人| 两个人的视频大全免费| 丰满人妻熟妇乱又伦精品不卡| 日韩欧美免费精品| 国产爱豆传媒在线观看 | 一个人免费在线观看的高清视频| 国产精品久久久久久久电影 | 成年免费大片在线观看| 亚洲精品美女久久av网站| 在线观看66精品国产| 性欧美人与动物交配| 国产亚洲欧美在线一区二区| 一区二区三区激情视频| 熟女少妇亚洲综合色aaa.| 听说在线观看完整版免费高清| 老司机靠b影院| 日韩欧美国产在线观看|