• <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ù)值計算
    精品卡一卡二卡四卡免费| 亚洲,一卡二卡三卡| √禁漫天堂资源中文www| 少妇的逼好多水| 国产精品欧美亚洲77777| 搡老乐熟女国产| 成人影院久久| 综合色丁香网| 看免费av毛片| 国产精品女同一区二区软件| 久久青草综合色| 大片免费播放器 马上看| 91精品国产国语对白视频| 一区二区三区四区激情视频| 日本av手机在线免费观看| 最近中文字幕高清免费大全6| 99久久综合免费| 大香蕉97超碰在线| 亚洲人与动物交配视频| 在线免费观看不下载黄p国产| 韩国高清视频一区二区三区| 蜜桃国产av成人99| 中文字幕亚洲精品专区| 女人被躁到高潮嗷嗷叫费观| 女人被躁到高潮嗷嗷叫费观| 亚洲伊人色综图| 久久久a久久爽久久v久久| 精品一区二区三卡| 国产精品不卡视频一区二区| 成年人免费黄色播放视频| 99精国产麻豆久久婷婷| av.在线天堂| 一级,二级,三级黄色视频| 观看av在线不卡| 菩萨蛮人人尽说江南好唐韦庄| 五月玫瑰六月丁香| 欧美国产精品一级二级三级| 日韩成人av中文字幕在线观看| 高清不卡的av网站| 人妻人人澡人人爽人人| www日本在线高清视频| 最近2019中文字幕mv第一页| 亚洲五月色婷婷综合| 国产精品久久久久久久久免| tube8黄色片| av视频免费观看在线观看| 人人妻人人爽人人添夜夜欢视频| 免费大片黄手机在线观看| 人体艺术视频欧美日本| 精品国产国语对白av| av不卡在线播放| 亚洲精品第二区| 香蕉精品网在线| 久久精品国产鲁丝片午夜精品| h视频一区二区三区| 高清视频免费观看一区二区| 亚洲经典国产精华液单| 国产免费现黄频在线看| 男女下面插进去视频免费观看 | 日韩精品有码人妻一区| 爱豆传媒免费全集在线观看| 少妇的丰满在线观看| 制服人妻中文乱码| 啦啦啦视频在线资源免费观看| 国产精品久久久久久精品电影小说| 性色av一级| 国产亚洲精品第一综合不卡 | 丰满迷人的少妇在线观看| 好男人视频免费观看在线| 男人爽女人下面视频在线观看| 国产精品国产三级国产专区5o| 亚洲精品一区蜜桃| 久久婷婷青草| 国产精品.久久久| 成人毛片a级毛片在线播放| av免费观看日本| 日本vs欧美在线观看视频| 日韩av免费高清视频| 久久99热这里只频精品6学生| 国产亚洲精品久久久com| 久久精品久久久久久久性| 国产精品蜜桃在线观看| 午夜福利视频在线观看免费| 欧美变态另类bdsm刘玥| 成人亚洲精品一区在线观看| 亚洲av欧美aⅴ国产| 婷婷色综合大香蕉| 日日撸夜夜添| 午夜久久久在线观看| 日日爽夜夜爽网站| 五月天丁香电影| 一边摸一边做爽爽视频免费| 亚洲精品一区蜜桃| 国产av一区二区精品久久| 久久精品人人爽人人爽视色| 男的添女的下面高潮视频| 一区二区av电影网| 菩萨蛮人人尽说江南好唐韦庄| 激情五月婷婷亚洲| 美女福利国产在线| 最近的中文字幕免费完整| 中文字幕免费在线视频6| 两个人免费观看高清视频| 夫妻性生交免费视频一级片| 久久免费观看电影| 日韩大片免费观看网站| 成人漫画全彩无遮挡| 在现免费观看毛片| 婷婷色av中文字幕| 午夜福利视频在线观看免费| 色5月婷婷丁香| 国产白丝娇喘喷水9色精品| 中文天堂在线官网| 一区二区av电影网| 狂野欧美激情性xxxx在线观看| 咕卡用的链子| 韩国av在线不卡| 亚洲美女搞黄在线观看| 美女中出高潮动态图| www.av在线官网国产| 国产永久视频网站| 精品久久久久久电影网| 国产日韩一区二区三区精品不卡| 亚洲av免费高清在线观看| 韩国精品一区二区三区 | 伦理电影免费视频| 亚洲国产最新在线播放| 久久99热6这里只有精品| 亚洲伊人色综图| 亚洲在久久综合| 伊人亚洲综合成人网| 精品一区二区免费观看| 国产日韩一区二区三区精品不卡| 亚洲高清免费不卡视频| 丝袜脚勾引网站| 日日啪夜夜爽| 日韩成人伦理影院| 寂寞人妻少妇视频99o| 亚洲美女视频黄频| 久久热在线av| 日韩 亚洲 欧美在线| 极品少妇高潮喷水抽搐| 成人午夜精彩视频在线观看| 亚洲人成77777在线视频| 免费人妻精品一区二区三区视频| 成人国产av品久久久| 久热这里只有精品99| 狠狠婷婷综合久久久久久88av| 最新中文字幕久久久久| 国产精品一区二区在线不卡| 亚洲精品日本国产第一区| 久久精品国产鲁丝片午夜精品| 亚洲,欧美,日韩| 男女下面插进去视频免费观看 | 亚洲国产看品久久| 少妇高潮的动态图| 少妇 在线观看| 久久精品国产a三级三级三级| 99re6热这里在线精品视频| 视频在线观看一区二区三区| 亚洲精品色激情综合| 久久久久久久精品精品| 在线免费观看不下载黄p国产| 午夜影院在线不卡| 久久久久久久久久久免费av| 欧美人与性动交α欧美精品济南到 | 日本色播在线视频| 久久ye,这里只有精品| 久久久久网色| 亚洲国产毛片av蜜桃av| 99国产精品免费福利视频| 一级片免费观看大全| 日本猛色少妇xxxxx猛交久久| 一二三四中文在线观看免费高清| 亚洲五月色婷婷综合| 丁香六月天网| 久久国产精品大桥未久av| 亚洲熟女精品中文字幕| 久久ye,这里只有精品| 97精品久久久久久久久久精品| 一区二区日韩欧美中文字幕 | 久久综合国产亚洲精品| 青春草国产在线视频| 日韩精品免费视频一区二区三区 | 国产亚洲最大av| 美国免费a级毛片| 亚洲伊人色综图| 90打野战视频偷拍视频| 如日韩欧美国产精品一区二区三区| 亚洲精品国产色婷婷电影| 一本大道久久a久久精品| 伊人久久国产一区二区| 少妇高潮的动态图| 制服诱惑二区| 久久人人爽av亚洲精品天堂| 久久久久网色| 免费少妇av软件| videosex国产| 欧美bdsm另类| 91精品国产国语对白视频| 晚上一个人看的免费电影| 国产在线视频一区二区| 美女主播在线视频| 国产成人av激情在线播放| 少妇人妻精品综合一区二区| 久久精品熟女亚洲av麻豆精品| 满18在线观看网站| 精品久久国产蜜桃| 中文乱码字字幕精品一区二区三区| av在线app专区| 免费观看av网站的网址| 亚洲av.av天堂| 久久精品国产综合久久久 | 日本wwww免费看| av卡一久久| 日韩制服丝袜自拍偷拍| 国产探花极品一区二区| 亚洲激情五月婷婷啪啪| 亚洲欧洲日产国产| 欧美国产精品va在线观看不卡| 亚洲欧美色中文字幕在线| 狠狠婷婷综合久久久久久88av| 视频区图区小说| 嫩草影院入口| 精品99又大又爽又粗少妇毛片| 亚洲中文av在线| 日本午夜av视频| 人人妻人人添人人爽欧美一区卜| 丝袜在线中文字幕| 99re6热这里在线精品视频| 中文欧美无线码| 在线观看国产h片| 51国产日韩欧美| 赤兔流量卡办理| av在线观看视频网站免费| 国产综合精华液| 美女视频免费永久观看网站| 日本黄色日本黄色录像| 国产一区二区在线观看日韩| 91精品国产国语对白视频| 高清毛片免费看| 国产av码专区亚洲av| 草草在线视频免费看| 美女福利国产在线| 亚洲国产精品国产精品| videosex国产| 久久久久久伊人网av| 99久久精品国产国产毛片| 午夜精品国产一区二区电影| 看免费成人av毛片| 国产极品粉嫩免费观看在线| 9191精品国产免费久久| 国产在线一区二区三区精| 国产一区二区激情短视频 | freevideosex欧美| 激情五月婷婷亚洲| 久久久精品区二区三区| 日本黄色日本黄色录像| 欧美bdsm另类| 精品国产国语对白av| 免费黄频网站在线观看国产| 永久网站在线| 久久久精品94久久精品| 人妻系列 视频| 亚洲国产精品999| 欧美人与性动交α欧美精品济南到 | 免费观看av网站的网址| 日韩av免费高清视频| 少妇熟女欧美另类| 热99久久久久精品小说推荐| 日产精品乱码卡一卡2卡三| 亚洲精品久久午夜乱码| 日韩制服丝袜自拍偷拍| 欧美日韩综合久久久久久| 国产精品国产三级国产av玫瑰| 夜夜骑夜夜射夜夜干| 成人综合一区亚洲| 欧美变态另类bdsm刘玥| 欧美精品人与动牲交sv欧美| 精品第一国产精品| 国产精品蜜桃在线观看| 免费高清在线观看日韩| 成人国语在线视频| av不卡在线播放| 亚洲欧美成人精品一区二区| 欧美精品一区二区免费开放| 男男h啪啪无遮挡| 男女午夜视频在线观看 | 国产精品熟女久久久久浪| 午夜福利乱码中文字幕| 欧美精品一区二区免费开放| 黄色怎么调成土黄色| 中文字幕免费在线视频6| 亚洲欧美精品自产自拍| www.熟女人妻精品国产 | 少妇被粗大猛烈的视频| 久久午夜福利片| 十八禁高潮呻吟视频| 亚洲精品aⅴ在线观看| 丰满迷人的少妇在线观看| 欧美成人精品欧美一级黄| 99视频精品全部免费 在线| 国产男人的电影天堂91| 一边摸一边做爽爽视频免费| 十八禁网站网址无遮挡| 极品人妻少妇av视频| 大香蕉97超碰在线| 午夜视频国产福利| 美女主播在线视频| 午夜老司机福利剧场| 赤兔流量卡办理| 少妇被粗大猛烈的视频| 国产亚洲精品久久久com| 国产精品人妻久久久久久| 国产一区二区在线观看日韩| 欧美日韩国产mv在线观看视频| 亚洲丝袜综合中文字幕| 亚洲美女黄色视频免费看| 亚洲国产最新在线播放| 午夜免费观看性视频| 亚洲熟女精品中文字幕| 美女内射精品一级片tv| 久久婷婷青草| 欧美精品一区二区免费开放| 91精品国产国语对白视频| videos熟女内射| 亚洲精品日韩在线中文字幕| 免费观看无遮挡的男女| av在线老鸭窝| 2022亚洲国产成人精品| xxxhd国产人妻xxx| tube8黄色片| 欧美日韩成人在线一区二区| 午夜91福利影院| 久久精品熟女亚洲av麻豆精品| 亚洲精品久久久久久婷婷小说| 国精品久久久久久国模美| 在线观看人妻少妇| 日韩av不卡免费在线播放| 看免费av毛片| 成年动漫av网址| 少妇被粗大猛烈的视频| 五月天丁香电影| 日韩成人伦理影院| 最近最新中文字幕免费大全7| 少妇熟女欧美另类| 久久久久久久久久久免费av| 国产激情久久老熟女| 极品少妇高潮喷水抽搐| 日产精品乱码卡一卡2卡三| 欧美人与性动交α欧美软件 | 欧美少妇被猛烈插入视频| 久久久久久人人人人人| 国产精品 国内视频| 久久青草综合色| av在线老鸭窝| 五月伊人婷婷丁香| 九色亚洲精品在线播放| 免费高清在线观看日韩| 黄色怎么调成土黄色| 精品人妻在线不人妻| 日韩制服丝袜自拍偷拍| 黑人巨大精品欧美一区二区蜜桃 | 亚洲国产精品成人久久小说| 日韩熟女老妇一区二区性免费视频| 两个人看的免费小视频| 免费人成在线观看视频色| 精品一区二区三区四区五区乱码 | 精品亚洲成国产av| 欧美日韩成人在线一区二区| 亚洲综合精品二区| 国产又色又爽无遮挡免| 最近最新中文字幕免费大全7| 中文字幕亚洲精品专区| 2022亚洲国产成人精品| 中文乱码字字幕精品一区二区三区| 99国产精品免费福利视频| 国产成人精品婷婷| 免费黄网站久久成人精品| 午夜av观看不卡| 亚洲av免费高清在线观看| 免费av不卡在线播放| 国产精品女同一区二区软件| 日韩电影二区| 日韩精品有码人妻一区| 国产黄色视频一区二区在线观看| 久久人妻熟女aⅴ| 美女福利国产在线| 国产xxxxx性猛交| 欧美日韩精品成人综合77777| 99热6这里只有精品| 男女啪啪激烈高潮av片| 亚洲四区av| 日韩,欧美,国产一区二区三区| 搡老乐熟女国产| 精品一品国产午夜福利视频| 亚洲国产av影院在线观看| 另类精品久久| 久久久亚洲精品成人影院| 大香蕉97超碰在线| 久久久久久久久久人人人人人人| 黄色毛片三级朝国网站| 欧美变态另类bdsm刘玥| 五月玫瑰六月丁香| 亚洲av在线观看美女高潮| 亚洲一区二区三区欧美精品| 男男h啪啪无遮挡| 高清视频免费观看一区二区| 啦啦啦在线观看免费高清www| 乱人伦中国视频| freevideosex欧美| 一区二区三区乱码不卡18| 精品视频人人做人人爽| 一区在线观看完整版| 少妇熟女欧美另类| 国产探花极品一区二区| 国产精品欧美亚洲77777| 在线亚洲精品国产二区图片欧美| 男女下面插进去视频免费观看 | 午夜免费观看性视频| 美女大奶头黄色视频| av国产久精品久网站免费入址| 一本一本久久a久久精品综合妖精 国产伦在线观看视频一区 | 婷婷色综合大香蕉| 一二三四在线观看免费中文在 | 国产精品人妻久久久影院| 久久久久久人妻| 69精品国产乱码久久久| 国产亚洲一区二区精品| 在线观看免费日韩欧美大片| 美国免费a级毛片| 色网站视频免费| 久久 成人 亚洲| 欧美老熟妇乱子伦牲交| 日日撸夜夜添| 18禁动态无遮挡网站| 亚洲精品久久久久久婷婷小说| 欧美精品亚洲一区二区| 亚洲欧美一区二区三区黑人 | 国产黄色视频一区二区在线观看| 在线看a的网站| 热99国产精品久久久久久7| 全区人妻精品视频| 美女福利国产在线| 亚洲中文av在线| 国产成人av激情在线播放| 午夜福利乱码中文字幕| 久久久国产一区二区| 国产在线一区二区三区精| 国产一区二区在线观看av| 女人被躁到高潮嗷嗷叫费观| 久久久亚洲精品成人影院| 国产精品一区www在线观看| 欧美最新免费一区二区三区| 国产精品久久久久久av不卡| 久久久久国产网址| av天堂久久9| 97超碰精品成人国产| 在线观看三级黄色| 亚洲伊人色综图| 国产成人91sexporn| 久久精品国产a三级三级三级| 久久久久精品久久久久真实原创| 热re99久久精品国产66热6| 91国产中文字幕| 免费看av在线观看网站| 亚洲精品,欧美精品| 国产精品99久久99久久久不卡 | 涩涩av久久男人的天堂| 最黄视频免费看| 麻豆乱淫一区二区| 寂寞人妻少妇视频99o| 欧美精品av麻豆av| 亚洲成人一二三区av| 香蕉丝袜av| 久久97久久精品| 国产日韩欧美亚洲二区| 久久久亚洲精品成人影院| av国产久精品久网站免费入址| 五月开心婷婷网| 亚洲综合色惰| 精品久久久精品久久久| 99国产精品免费福利视频| 黑人高潮一二区| 内地一区二区视频在线| 亚洲欧洲国产日韩| 自线自在国产av| 欧美bdsm另类| 999精品在线视频| 精品人妻偷拍中文字幕| 丰满少妇做爰视频| 高清视频免费观看一区二区| 一本—道久久a久久精品蜜桃钙片| 超碰97精品在线观看| av在线观看视频网站免费| 日韩一区二区三区影片| 国产免费福利视频在线观看| 亚洲精品国产av成人精品| 高清黄色对白视频在线免费看| 免费日韩欧美在线观看| 久久久国产精品麻豆| 亚洲高清免费不卡视频| 性高湖久久久久久久久免费观看| 嫩草影院入口| 免费女性裸体啪啪无遮挡网站| 国产精品久久久久成人av| 日本-黄色视频高清免费观看| 国产成人a∨麻豆精品| 久久女婷五月综合色啪小说| 国产免费一区二区三区四区乱码| 色94色欧美一区二区| 免费观看无遮挡的男女| 免费大片黄手机在线观看| 天堂8中文在线网| 九草在线视频观看| 一二三四中文在线观看免费高清| 国产精品不卡视频一区二区| 国产亚洲午夜精品一区二区久久| 久久久久久久亚洲中文字幕| 国产亚洲欧美精品永久| 国产午夜精品一二区理论片| 大码成人一级视频| 亚洲丝袜综合中文字幕| 日韩大片免费观看网站| 国产精品久久久av美女十八| 自拍欧美九色日韩亚洲蝌蚪91| 巨乳人妻的诱惑在线观看| 2018国产大陆天天弄谢| 80岁老熟妇乱子伦牲交| 黑人欧美特级aaaaaa片| 国产成人一区二区在线| 啦啦啦视频在线资源免费观看| 九草在线视频观看| 9191精品国产免费久久| 人人妻人人澡人人爽人人夜夜| 99九九在线精品视频| 欧美日韩综合久久久久久| 色婷婷久久久亚洲欧美| 亚洲av国产av综合av卡| 99热全是精品| 一区二区三区精品91| 久久久久久久大尺度免费视频| 视频区图区小说| 国产爽快片一区二区三区| 亚洲av.av天堂| 涩涩av久久男人的天堂| 久久久久精品久久久久真实原创| 午夜福利视频在线观看免费| 一级a做视频免费观看| h视频一区二区三区| 两个人看的免费小视频| 久久久国产精品麻豆| 最近手机中文字幕大全| 永久免费av网站大全| 久久久久精品性色| av在线观看视频网站免费| 内地一区二区视频在线| 亚洲成国产人片在线观看| 青春草视频在线免费观看| 国产亚洲最大av| 亚洲精品一二三| 国产成人免费无遮挡视频| 一级毛片黄色毛片免费观看视频| 国产精品一区www在线观看| 国产精品 国内视频| 建设人人有责人人尽责人人享有的| 国产精品国产三级专区第一集| 青春草国产在线视频| 五月伊人婷婷丁香| 免费黄网站久久成人精品| 中文精品一卡2卡3卡4更新| 免费看不卡的av| 黄片播放在线免费| 免费av不卡在线播放| 亚洲国产最新在线播放| 亚洲av日韩在线播放| 色婷婷久久久亚洲欧美| 男女免费视频国产| 大香蕉97超碰在线| 日韩一本色道免费dvd| 在线观看www视频免费| 亚洲熟女精品中文字幕| 久久精品久久久久久久性| 亚洲美女视频黄频| 极品少妇高潮喷水抽搐| 99热网站在线观看| 新久久久久国产一级毛片| 激情五月婷婷亚洲| 成年人免费黄色播放视频| 毛片一级片免费看久久久久| 国产成人精品福利久久| 欧美 亚洲 国产 日韩一| 亚洲一区二区三区欧美精品| 飞空精品影院首页| 久久久久视频综合| 热99久久久久精品小说推荐| 国产精品久久久久久av不卡| 亚洲欧美精品自产自拍| 麻豆精品久久久久久蜜桃| 日韩精品有码人妻一区| 国产淫语在线视频| 欧美日韩视频高清一区二区三区二| 国产一区二区激情短视频 | 一区二区三区四区激情视频| 亚洲人成77777在线视频| 成年人午夜在线观看视频| 亚洲av中文av极速乱| 国产69精品久久久久777片| 一区二区av电影网| 国产日韩欧美在线精品| 黄片播放在线免费| 伦理电影大哥的女人| 亚洲欧美清纯卡通| 伊人亚洲综合成人网| 日韩一区二区三区影片|