• <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在线天堂中文字幕| 黑人欧美特级aaaaaa片| 亚洲国产日韩欧美精品在线观看 | 搡老岳熟女国产| 麻豆国产97在线/欧美| 丝袜人妻中文字幕| 精品久久久久久成人av| 香蕉久久夜色| 丁香欧美五月| 51午夜福利影视在线观看| 精品无人区乱码1区二区| 国产激情欧美一区二区| 人妻夜夜爽99麻豆av| 宅男免费午夜| 高清毛片免费观看视频网站| 亚洲男人的天堂狠狠| 不卡av一区二区三区| 日本五十路高清| 女生性感内裤真人,穿戴方法视频| 色精品久久人妻99蜜桃| 十八禁人妻一区二区| 欧美色视频一区免费| 桃色一区二区三区在线观看| 国产av麻豆久久久久久久| 成人特级黄色片久久久久久久| 久久精品夜夜夜夜夜久久蜜豆| 一夜夜www| 两人在一起打扑克的视频| 看片在线看免费视频| 在线观看免费视频日本深夜| 最新中文字幕久久久久 | 欧美最黄视频在线播放免费| 三级国产精品欧美在线观看 | 日韩欧美在线乱码| 亚洲美女黄片视频| 岛国在线免费视频观看| 国产久久久一区二区三区| 国产精品国产高清国产av| 亚洲国产精品sss在线观看| 黄色视频,在线免费观看| 国产成人一区二区三区免费视频网站| 国产亚洲精品一区二区www| 男人舔奶头视频| 久久久久久久午夜电影| 在线看三级毛片| 小蜜桃在线观看免费完整版高清| 18禁观看日本| 国产精品精品国产色婷婷| 在线免费观看不下载黄p国产 | 亚洲国产精品成人综合色| 又爽又黄无遮挡网站| 成人永久免费在线观看视频| 在线视频色国产色| 母亲3免费完整高清在线观看| 色av中文字幕| 网址你懂的国产日韩在线| 最新中文字幕久久久久 | 特级一级黄色大片| 国产主播在线观看一区二区| 老鸭窝网址在线观看| 法律面前人人平等表现在哪些方面| 成人国产综合亚洲| 色吧在线观看| 可以在线观看的亚洲视频| 国产av麻豆久久久久久久| 18禁观看日本| 日日干狠狠操夜夜爽| 成人三级做爰电影| 国产乱人伦免费视频| 黄色成人免费大全| 午夜免费观看网址| 1024手机看黄色片| 99久久99久久久精品蜜桃| 欧美又色又爽又黄视频| 国产高清有码在线观看视频| 午夜成年电影在线免费观看| 小说图片视频综合网站| 免费高清视频大片| xxxwww97欧美| 国产黄片美女视频| 长腿黑丝高跟| 最新在线观看一区二区三区| 免费看光身美女| 国产一区在线观看成人免费| 一个人观看的视频www高清免费观看 | 亚洲午夜精品一区,二区,三区| 国产又色又爽无遮挡免费看| 欧美xxxx黑人xx丫x性爽| 国产精品永久免费网站| 国产精品亚洲一级av第二区| 十八禁网站免费在线| 亚洲人成电影免费在线| 免费在线观看视频国产中文字幕亚洲| 国产高清激情床上av| 18美女黄网站色大片免费观看| av黄色大香蕉| 国产一级毛片七仙女欲春2| 亚洲精品美女久久久久99蜜臀| 国产欧美日韩一区二区精品| 欧美高清成人免费视频www| 无遮挡黄片免费观看| 日日干狠狠操夜夜爽| 亚洲精品乱码久久久v下载方式 | 天堂动漫精品| 日本撒尿小便嘘嘘汇集6| 亚洲一区二区三区不卡视频| 国产精品香港三级国产av潘金莲| 高清毛片免费观看视频网站| 窝窝影院91人妻| www.自偷自拍.com| 久久亚洲真实| 国产探花在线观看一区二区| 亚洲精品在线观看二区| 啪啪无遮挡十八禁网站| 国产精品电影一区二区三区| 亚洲真实伦在线观看| 男女之事视频高清在线观看| 亚洲欧洲精品一区二区精品久久久| 亚洲电影在线观看av| 长腿黑丝高跟| 精品午夜福利视频在线观看一区| 亚洲欧美精品综合一区二区三区| 中文字幕人妻丝袜一区二区| 午夜激情福利司机影院| 激情在线观看视频在线高清| 在线免费观看的www视频| 麻豆国产av国片精品| 波多野结衣巨乳人妻| 精品久久久久久,| 色播亚洲综合网| 国产av一区在线观看免费| 国产99白浆流出| 淫秽高清视频在线观看| 国产久久久一区二区三区| 精品一区二区三区视频在线观看免费| 听说在线观看完整版免费高清| 国产单亲对白刺激| 在线a可以看的网站| 国产成年人精品一区二区| 国产麻豆成人av免费视频| 亚洲专区字幕在线| 99热只有精品国产| 97超级碰碰碰精品色视频在线观看| www日本在线高清视频| 久久精品国产亚洲av香蕉五月| 久久久久九九精品影院| 成人欧美大片| 国产精品香港三级国产av潘金莲| 久久久色成人| 国产精品乱码一区二三区的特点| 美女黄网站色视频| 两人在一起打扑克的视频| 国语自产精品视频在线第100页| 欧美黑人欧美精品刺激| 又爽又黄无遮挡网站| 一级黄色大片毛片| 久久久久久国产a免费观看| 又爽又黄无遮挡网站| 日本成人三级电影网站| 露出奶头的视频| 午夜福利成人在线免费观看| 欧美日韩黄片免| 舔av片在线| 久久久久久人人人人人| 久久亚洲真实| 欧美日本视频| 岛国视频午夜一区免费看| 美女 人体艺术 gogo| 国产精品乱码一区二三区的特点| 国产精品一区二区三区四区久久| 欧美3d第一页| 国产一区二区三区视频了| 国产一区二区三区在线臀色熟女| 欧美国产日韩亚洲一区| 此物有八面人人有两片| 高清毛片免费观看视频网站| 亚洲无线观看免费| 久久天堂一区二区三区四区| 无人区码免费观看不卡| 九九在线视频观看精品| 免费观看人在逋| 制服人妻中文乱码| 国产三级中文精品| 午夜免费激情av| 99国产精品99久久久久| 91麻豆精品激情在线观看国产| 欧美日韩一级在线毛片| 国产精品1区2区在线观看.| 久久久水蜜桃国产精品网| 桃色一区二区三区在线观看| 一区二区三区激情视频| 观看美女的网站| 亚洲天堂国产精品一区在线| 日韩精品青青久久久久久| 欧美zozozo另类| 人妻丰满熟妇av一区二区三区| 欧美av亚洲av综合av国产av| 欧美色视频一区免费| 日韩国内少妇激情av| 99精品欧美一区二区三区四区| 国产日本99.免费观看| 日本成人三级电影网站| 亚洲七黄色美女视频| 天堂网av新在线| 丁香欧美五月| 很黄的视频免费| 亚洲精品一卡2卡三卡4卡5卡| 女生性感内裤真人,穿戴方法视频| 亚洲av电影不卡..在线观看| 日本免费一区二区三区高清不卡| 夜夜爽天天搞| 国产伦在线观看视频一区| 1024香蕉在线观看| 欧美色欧美亚洲另类二区| 小说图片视频综合网站| 欧美zozozo另类| 99久久精品一区二区三区| 男女下面进入的视频免费午夜| 别揉我奶头~嗯~啊~动态视频| 女同久久另类99精品国产91| 俄罗斯特黄特色一大片| 国产精品电影一区二区三区| 在线永久观看黄色视频| 性欧美人与动物交配| 动漫黄色视频在线观看| 51午夜福利影视在线观看| 男女那种视频在线观看| 日韩欧美精品v在线| 欧美乱妇无乱码| 亚洲欧美日韩东京热| 亚洲国产欧洲综合997久久,| 婷婷丁香在线五月| 免费看a级黄色片| 日本一二三区视频观看| 免费在线观看视频国产中文字幕亚洲| 一个人免费在线观看的高清视频| 真实男女啪啪啪动态图| 亚洲美女视频黄频| 成年女人毛片免费观看观看9| 亚洲中文字幕一区二区三区有码在线看 | 国产综合懂色| 久久精品aⅴ一区二区三区四区| 很黄的视频免费| 日本黄色视频三级网站网址| 成人午夜高清在线视频| 亚洲一区二区三区不卡视频| 久久香蕉国产精品| 黑人欧美特级aaaaaa片| 国产欧美日韩一区二区三| 最新在线观看一区二区三区| 国产爱豆传媒在线观看| 亚洲色图 男人天堂 中文字幕| www日本在线高清视频| 国产精品日韩av在线免费观看| 日韩欧美在线二视频| 亚洲真实伦在线观看| 后天国语完整版免费观看| 老司机在亚洲福利影院| 日日干狠狠操夜夜爽| 小蜜桃在线观看免费完整版高清| 深夜精品福利| 九色国产91popny在线| 欧美在线黄色| 女人被狂操c到高潮| 国产精品综合久久久久久久免费| 99riav亚洲国产免费| 少妇熟女aⅴ在线视频| 亚洲国产欧美人成| 欧美一级a爱片免费观看看| 午夜a级毛片| 国产日本99.免费观看| 亚洲一区二区三区色噜噜| 色播亚洲综合网| 欧美乱码精品一区二区三区| 国产成人福利小说| 欧美大码av| www日本黄色视频网| 国产成人一区二区三区免费视频网站| 精品一区二区三区视频在线观看免费| 啪啪无遮挡十八禁网站| 欧美日韩国产亚洲二区| 亚洲性夜色夜夜综合| 天堂网av新在线| 国产精品久久久人人做人人爽| 日本黄色片子视频| 久久久久亚洲av毛片大全| 一本一本综合久久| 变态另类成人亚洲欧美熟女| 99国产精品一区二区蜜桃av| 欧美乱色亚洲激情| 俄罗斯特黄特色一大片| 久久精品国产综合久久久| 亚洲熟妇中文字幕五十中出| av黄色大香蕉| 香蕉av资源在线| 亚洲精品456在线播放app | 国产精品永久免费网站| 草草在线视频免费看| 美女高潮的动态| 成人国产一区最新在线观看| 精品欧美国产一区二区三| 一个人看的www免费观看视频| 我要搜黄色片| 国产伦一二天堂av在线观看| av在线天堂中文字幕| 真人做人爱边吃奶动态| АⅤ资源中文在线天堂| 人妻丰满熟妇av一区二区三区| 久久天堂一区二区三区四区| 高清在线国产一区| 国产免费av片在线观看野外av| 精品久久蜜臀av无| 国产精品久久久av美女十八| 国产69精品久久久久777片 | 国产av一区在线观看免费| 午夜福利欧美成人| 99热这里只有精品一区 | 精品一区二区三区四区五区乱码| 久久久久久九九精品二区国产| 青草久久国产| 亚洲电影在线观看av| 欧美xxxx黑人xx丫x性爽| 成人国产综合亚洲| 一个人观看的视频www高清免费观看 | avwww免费| 中出人妻视频一区二区| 一级a爱片免费观看的视频| 老司机午夜十八禁免费视频| 日本五十路高清| 美女 人体艺术 gogo| 国内少妇人妻偷人精品xxx网站 | 国产成人系列免费观看| 精品久久蜜臀av无| 国产不卡一卡二| 国产成人欧美在线观看| 国产精品久久久av美女十八| 久久久久国产一级毛片高清牌| 午夜影院日韩av| 很黄的视频免费| 男人舔奶头视频| 国产成人aa在线观看| 久久久久久人人人人人| 亚洲人与动物交配视频| 日韩欧美国产一区二区入口| 午夜激情欧美在线| 午夜福利在线在线| 桃色一区二区三区在线观看| 日韩成人在线观看一区二区三区| 性欧美人与动物交配| 亚洲美女黄片视频| 国产v大片淫在线免费观看| 色在线成人网| 国产单亲对白刺激| 日韩 欧美 亚洲 中文字幕| 欧美黄色淫秽网站| 国产真人三级小视频在线观看| 女人被狂操c到高潮| 亚洲美女黄片视频| 欧美成人免费av一区二区三区| 97人妻精品一区二区三区麻豆| 桃色一区二区三区在线观看| 亚洲国产精品sss在线观看| 国产精品亚洲av一区麻豆| 久久久国产欧美日韩av| 国产成人欧美在线观看| 在线观看一区二区三区| 亚洲精品乱码久久久v下载方式 | 午夜精品在线福利| 日韩免费av在线播放| bbb黄色大片| 国产乱人视频| xxxwww97欧美| 美女免费视频网站| 校园春色视频在线观看| 国产精品av视频在线免费观看| 午夜成年电影在线免费观看| 宅男免费午夜| 中文字幕久久专区| 成年女人永久免费观看视频| 中文资源天堂在线| 黄频高清免费视频| 网址你懂的国产日韩在线| 最近最新免费中文字幕在线| 熟女电影av网| 久久久色成人| 男人舔奶头视频| 97超级碰碰碰精品色视频在线观看| 2021天堂中文幕一二区在线观| xxxwww97欧美| 中文字幕人成人乱码亚洲影| www.精华液| 亚洲美女视频黄频| 国产精品一区二区免费欧美| 天堂√8在线中文| 国产1区2区3区精品| a级毛片a级免费在线| 亚洲av成人一区二区三| 在线观看美女被高潮喷水网站 | 特大巨黑吊av在线直播| 色综合亚洲欧美另类图片| 亚洲专区国产一区二区| 国产亚洲av嫩草精品影院| 成年女人看的毛片在线观看| 欧美日韩一级在线毛片| 欧美日本亚洲视频在线播放| 精品日产1卡2卡| 国产不卡一卡二| 国产精品永久免费网站| 波多野结衣高清作品| 成人一区二区视频在线观看| or卡值多少钱| 亚洲最大成人中文| 母亲3免费完整高清在线观看| 日本成人三级电影网站| 久久久国产成人免费| 欧美+亚洲+日韩+国产| 色视频www国产| 搡老岳熟女国产| 法律面前人人平等表现在哪些方面| 久久天堂一区二区三区四区| 99精品在免费线老司机午夜| 天天添夜夜摸| 成人国产综合亚洲| 国产精华一区二区三区| 亚洲专区中文字幕在线| 欧美3d第一页| 久久99热这里只有精品18| 男女那种视频在线观看| 视频区欧美日本亚洲| 国产熟女xx| а√天堂www在线а√下载| 老熟妇乱子伦视频在线观看| 亚洲成av人片免费观看| 97碰自拍视频| 老司机福利观看| 久久天堂一区二区三区四区| 高潮久久久久久久久久久不卡| 色综合站精品国产| 国产精品99久久久久久久久| 最近最新免费中文字幕在线| 每晚都被弄得嗷嗷叫到高潮| 两人在一起打扑克的视频| 国产一区在线观看成人免费| 久久人妻av系列| 床上黄色一级片| 国产熟女xx| 99久久无色码亚洲精品果冻| 精品一区二区三区视频在线观看免费| 中文字幕精品亚洲无线码一区| 麻豆成人午夜福利视频| 伊人久久大香线蕉亚洲五| 欧美中文日本在线观看视频| 免费av毛片视频| 亚洲欧美日韩高清在线视频| 十八禁网站免费在线| 日韩欧美在线乱码| 91久久精品国产一区二区成人 | av天堂在线播放| 欧美一级毛片孕妇| 国产精品,欧美在线| 99久久精品一区二区三区| 日韩欧美免费精品| 变态另类成人亚洲欧美熟女| 日韩有码中文字幕| 97人妻精品一区二区三区麻豆| 伦理电影免费视频| 视频区欧美日本亚洲| 亚洲最大成人中文| 给我免费播放毛片高清在线观看| 又粗又爽又猛毛片免费看| 首页视频小说图片口味搜索| 日本 av在线| 亚洲在线自拍视频| 免费大片18禁| 精品不卡国产一区二区三区| 无人区码免费观看不卡| 国产精品香港三级国产av潘金莲| 日韩大尺度精品在线看网址| 国产精品精品国产色婷婷| 网址你懂的国产日韩在线| 1024手机看黄色片| 欧美av亚洲av综合av国产av| 久久热在线av| 国产精品av久久久久免费| 国产精品一区二区三区四区免费观看 | 无人区码免费观看不卡| 少妇丰满av| 日本一本二区三区精品| 精品久久久久久久毛片微露脸| 久久久久国产精品人妻aⅴ院| 免费观看的影片在线观看| 男人舔奶头视频| 欧美日韩亚洲国产一区二区在线观看| 免费高清视频大片| 天天一区二区日本电影三级| 国产亚洲精品久久久com| 国产av在哪里看| 久久伊人香网站| 麻豆国产97在线/欧美| 亚洲最大成人中文| 欧美色视频一区免费| 亚洲欧美日韩无卡精品| 中文字幕最新亚洲高清| a级毛片a级免费在线| 老司机在亚洲福利影院| 可以在线观看的亚洲视频| 日韩 欧美 亚洲 中文字幕| 夜夜躁狠狠躁天天躁| 99国产综合亚洲精品| 日本成人三级电影网站| x7x7x7水蜜桃| 国产欧美日韩一区二区三| 人人妻人人看人人澡| 免费搜索国产男女视频| 91字幕亚洲| 欧美极品一区二区三区四区| 亚洲国产欧洲综合997久久,| 亚洲va日本ⅴa欧美va伊人久久| 久久九九热精品免费| 午夜免费成人在线视频| av黄色大香蕉| 国产精品精品国产色婷婷| 在线观看免费视频日本深夜| 欧美色欧美亚洲另类二区| 精品国产美女av久久久久小说| 看片在线看免费视频| 男女视频在线观看网站免费| 午夜视频精品福利| 国产精品久久电影中文字幕| 色综合站精品国产| 三级男女做爰猛烈吃奶摸视频| 麻豆国产av国片精品| 五月伊人婷婷丁香| 精品国内亚洲2022精品成人| 日日干狠狠操夜夜爽| 免费在线观看影片大全网站| 欧美成人免费av一区二区三区| АⅤ资源中文在线天堂| 欧美激情久久久久久爽电影| 操出白浆在线播放| 久久99热这里只有精品18| 色综合婷婷激情| 久久久久国内视频| 人人妻人人澡欧美一区二区| 国产aⅴ精品一区二区三区波| or卡值多少钱| 全区人妻精品视频| 欧美3d第一页| 免费高清视频大片| 白带黄色成豆腐渣| 国产精品美女特级片免费视频播放器 | 午夜激情欧美在线| 久久天堂一区二区三区四区| 国产伦一二天堂av在线观看| 久久久色成人| 午夜a级毛片| 中文字幕av在线有码专区| 免费看十八禁软件| 亚洲美女黄片视频| 精品国产乱子伦一区二区三区| а√天堂www在线а√下载| 日韩欧美 国产精品| www.精华液| 国产一区二区三区在线臀色熟女| 亚洲成人免费电影在线观看| 国产精品亚洲av一区麻豆| 亚洲av成人av| 色老头精品视频在线观看| 99精品久久久久人妻精品| 草草在线视频免费看| 欧美日本视频| 久久午夜亚洲精品久久| 男人舔女人下体高潮全视频| 色综合婷婷激情| 亚洲五月婷婷丁香| 国产精品久久久av美女十八| 一级毛片精品| 久久久久国内视频| 在线视频色国产色| 国产伦精品一区二区三区四那| a级毛片在线看网站| 免费看a级黄色片| 午夜精品一区二区三区免费看| cao死你这个sao货| 人人妻,人人澡人人爽秒播| 男女那种视频在线观看| 岛国在线免费视频观看| 此物有八面人人有两片| xxx96com| 亚洲国产精品成人综合色| 人人妻,人人澡人人爽秒播| 亚洲人成网站高清观看| 99国产精品一区二区蜜桃av| 久久久久亚洲av毛片大全| 淫妇啪啪啪对白视频| 免费搜索国产男女视频| 亚洲真实伦在线观看| 啦啦啦观看免费观看视频高清| 99国产精品一区二区蜜桃av| 午夜福利在线在线| 高潮久久久久久久久久久不卡| 欧美黑人巨大hd| 村上凉子中文字幕在线| 人妻丰满熟妇av一区二区三区| 日本 欧美在线| 无遮挡黄片免费观看| 免费一级毛片在线播放高清视频| 在线观看午夜福利视频| 日韩欧美精品v在线| 免费一级毛片在线播放高清视频| 高清毛片免费观看视频网站| 国产一区在线观看成人免费| 9191精品国产免费久久| 久久九九热精品免费| 国产伦精品一区二区三区视频9 |