鞠 楊，任張瑜，鄭江韜，毛靈濤，王 凱，周宏偉，謝和平

(1.中國(guó)礦業(yè)大學(xué)(北京) 煤炭資源與安全開采國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 100083；2.中國(guó)礦業(yè)大學(xué)(北京) 應(yīng)急管理與安全工程學(xué)院，北京 100083；3.中國(guó)礦業(yè)大學(xué)(北京) 能源與礦業(yè)學(xué)院，北京 100083；4.深圳大學(xué) 深地科學(xué)與綠色能源研究院，廣東 深圳 518060)

煤炭對(duì)我國(guó)經(jīng)濟(jì)社會(huì)發(fā)展與能源安全作出了巨大貢獻(xiàn)，在我國(guó)能源體系中起著壓艙石和兜底保障作用。我國(guó)能源資源稟賦特征是“富煤、貧油、少氣”，已探明化石能源儲(chǔ)量中煤炭占比達(dá)94%。習(xí)近平總書記在近期的中央經(jīng)濟(jì)工作會(huì)議上指出“要立足以煤為主的基本國(guó)情，抓好煤炭清潔高效利用，增加新能源消納能力，推動(dòng)煤炭和新能源優(yōu)化組合。要狠抓綠色低碳技術(shù)攻關(guān)”。2020年，習(xí)近平總書記在第75屆聯(lián)大一般性辯論上提出“二氧化碳排放力爭(zhēng)于2030年前達(dá)到峰值，努力爭(zhēng)取2060年前實(shí)現(xiàn)碳中和”。在“碳達(dá)峰、碳中和”國(guó)家重大戰(zhàn)略決策背景下，深入研究如何充分發(fā)揮煤炭的保障能源與戰(zhàn)略資源的作用、實(shí)現(xiàn)煤炭綠色低碳開發(fā)與潔凈利用迫在眉睫。

地球深部蘊(yùn)藏著豐富的煤炭資源，巖石是煤炭資源開發(fā)與工程建設(shè)的主要對(duì)象，狹義上，巖石指實(shí)驗(yàn)室尺度的巖塊，廣義上，巖石指工程尺度的巖體。煤炭開發(fā)與工程建設(shè)活動(dòng)打破了初始地應(yīng)力場(chǎng)的平衡狀態(tài)，引起巖石應(yīng)力場(chǎng)重分布，誘發(fā)沖擊地壓、瓦斯突出、礦震、生態(tài)環(huán)境破壞等重大災(zāi)害。然而，人們對(duì)煤炭賦存狀態(tài)、巖石構(gòu)造演化以及生產(chǎn)開發(fā)活動(dòng)的影響認(rèn)知還不夠，傳統(tǒng)開采理論與技術(shù)難以精準(zhǔn)探測(cè)和描述巖石復(fù)雜結(jié)構(gòu)及開采引發(fā)的演化過程，也缺乏對(duì)開采巖石力學(xué)行為與致災(zāi)機(jī)理的科學(xué)解析和準(zhǔn)確表征，因而難以對(duì)煤炭生產(chǎn)與工程建設(shè)活動(dòng)可能引發(fā)的工程災(zāi)害進(jìn)行超前預(yù)警和科學(xué)防控。解決這些基礎(chǔ)科學(xué)與技術(shù)難題已成為創(chuàng)新發(fā)展煤炭資源綠色智能安全開發(fā)理論與技術(shù)的重大基礎(chǔ)性前沿課題。

煤炭綠色智能安全開采基礎(chǔ)理論之一是巖石力學(xué)。然而，傳統(tǒng)的巖石力學(xué)理論在定量分析和描述生產(chǎn)開發(fā)活動(dòng)引發(fā)的巖石結(jié)構(gòu)與力學(xué)行為變化及災(zāi)變機(jī)理時(shí)存在諸多困難，例如，① 天然巖石賦含大量跨尺度、幾何非規(guī)則的孔隙、裂隙、節(jié)理或?qū)永淼确沁B續(xù)結(jié)構(gòu)，傳統(tǒng)的巖石力學(xué)方法，在連續(xù)均勻介質(zhì)假定的基礎(chǔ)上，在代表性體元(Representation Volume Element,RVE)尺度上，研究巖石的應(yīng)力、應(yīng)變、滲流與災(zāi)變等行為，建立巖石的宏觀本構(gòu)模型、強(qiáng)度準(zhǔn)則和失效判據(jù)等，難以準(zhǔn)確描述和定量解析巖石內(nèi)部非連續(xù)結(jié)構(gòu)特征及其對(duì)巖石力學(xué)行為與孕災(zāi)-致災(zāi)機(jī)理的影響。同時(shí)，由于人們對(duì)RVE尺度內(nèi)不連續(xù)結(jié)構(gòu)引發(fā)的應(yīng)力場(chǎng)、變形場(chǎng)和滲流場(chǎng)的變化認(rèn)識(shí)不足，導(dǎo)致確定合理的RVE尺度非常困難，基于連續(xù)均勻假設(shè)所建立的本構(gòu)模型、強(qiáng)度準(zhǔn)則和失效判據(jù)可能失去代表性，造成采用傳統(tǒng)的理論模型與數(shù)值方法解決工程實(shí)際問題時(shí)存在較大偏差。這些困難制約了構(gòu)建考慮內(nèi)部非連續(xù)結(jié)構(gòu)效應(yīng)的巖石力學(xué)與災(zāi)變分析體系；② 開采擾動(dòng)引發(fā)的巖石內(nèi)部非連續(xù)結(jié)構(gòu)演化、物質(zhì)-力學(xué)相互作用與多物理場(chǎng)耦合致災(zāi)是一個(gè)“看不見、摸不著”的“黑箱”過程，現(xiàn)有的巖石力學(xué)理論和方法難以直觀顯示和準(zhǔn)確描述，缺乏定量解析這種“黑箱”內(nèi)部應(yīng)力場(chǎng)、變形場(chǎng)和滲流場(chǎng)等多物理場(chǎng)的方法、理論和模型；③ 由于難以準(zhǔn)確獲取、直觀顯示和定量解析生產(chǎn)與建設(shè)活動(dòng)引發(fā)的巖石非連續(xù)結(jié)構(gòu)演化、多物理場(chǎng)耦合致災(zāi)機(jī)理及災(zāi)變力學(xué)特征，現(xiàn)有的礦山災(zāi)害防控多采用“隨采隨治”策略，難以超前研判不同開采方式和工藝可能引發(fā)的工程災(zāi)害的時(shí)間、位置和量級(jí)，難以實(shí)現(xiàn)開采災(zāi)害的超前預(yù)警和科學(xué)防控。巖石力學(xué)的這些基礎(chǔ)難題已成為制約煤炭綠色智能安全開采理論與災(zāi)害防控技術(shù)發(fā)展的瓶頸，迫切需要建立和發(fā)展顛覆性理論與變革性技術(shù)。

當(dāng)前新興技術(shù)的發(fā)展為直觀描述和定量表征巖石內(nèi)部復(fù)雜非連續(xù)結(jié)構(gòu)、應(yīng)力場(chǎng)、變形場(chǎng)、滲流場(chǎng)與災(zāi)變機(jī)理創(chuàng)造了條件。針對(duì)巖石內(nèi)部非連續(xù)結(jié)構(gòu)、物理場(chǎng)演化及災(zāi)變行為難以透明顯示和定量解析的難題，筆者采用高分辨μCT成像、3D數(shù)字重構(gòu)、3D打印、光彈性法、數(shù)字圖像相關(guān)法和三維流場(chǎng)示蹤等技術(shù)和方法，提出了巖石內(nèi)部非連續(xù)結(jié)構(gòu)、應(yīng)力場(chǎng)、應(yīng)變場(chǎng)、滲流場(chǎng)演化與多場(chǎng)耦合致災(zāi)機(jī)理的透明解析與透明推演方法，基本思路是：構(gòu)建含復(fù)雜非連續(xù)結(jié)構(gòu)RVE尺度的巖石數(shù)字模型和透明物理模型，研制具有與天然巖石相同或相近物理力學(xué)性能、具有應(yīng)力敏感性的透明3D打印材料，研究確定合理尺寸的非連續(xù)巖石RVE單元；采用透明可視化方法，在RVE尺度上，透明顯示和定量表征巖石非連續(xù)結(jié)構(gòu)特征及其對(duì)巖石內(nèi)部應(yīng)力場(chǎng)、應(yīng)變場(chǎng)、滲流場(chǎng)的影響與耦合致災(zāi)機(jī)理；基于巖石多物理場(chǎng)的透明定量解析，建立關(guān)聯(lián)內(nèi)部非連續(xù)結(jié)構(gòu)及物理力學(xué)效應(yīng)的巖石本構(gòu)模型、強(qiáng)度準(zhǔn)則及災(zāi)變判據(jù)；構(gòu)建RVE尺度的巖石透明模型與天然巖石非線性物理量之間的轉(zhuǎn)換關(guān)系及相似律模型，形成基于RVE尺度多物理場(chǎng)透明解析的非連續(xù)巖石力學(xué)與災(zāi)變分析理論的基礎(chǔ)框架。在此基礎(chǔ)上，將實(shí)驗(yàn)室透明物理模型研究與工程現(xiàn)場(chǎng)探測(cè)相結(jié)合，透明顯示并定量分析工程尺度下巖石結(jié)構(gòu)演化、能源物質(zhì)運(yùn)移、力學(xué)耦合作用與致災(zāi)過程。

近年來，國(guó)內(nèi)外學(xué)者圍繞巖石復(fù)雜結(jié)構(gòu)的三維數(shù)字重構(gòu)方法、3D打印材料研制、透明物理模型制備、尺度效應(yīng)、三維模型應(yīng)力場(chǎng)、變形場(chǎng)和滲流場(chǎng)演化以及災(zāi)變機(jī)理的透明定量表征開展了大量研究。筆者簡(jiǎn)要總結(jié)了國(guó)內(nèi)外學(xué)者及筆者在上述方面的研究進(jìn)展，分析了未來發(fā)展與面臨的挑戰(zhàn)。本文意在為發(fā)展非連續(xù)巖石力學(xué)與災(zāi)變分析理論、實(shí)現(xiàn)煤炭綠色智能安全開采及災(zāi)害超前預(yù)警防控提供理論研究與技術(shù)應(yīng)用參考。

1 巖石非連續(xù)結(jié)構(gòu)的數(shù)字模型

1.1 非連續(xù)結(jié)構(gòu)探測(cè)方法

巖石物理力學(xué)性質(zhì)與災(zāi)變行為很大程度上取決于其內(nèi)部非連續(xù)結(jié)構(gòu)。近年來，工業(yè)X射線CT掃描(X-ray CT)、磁共振成像(MRI)和聚焦離子束電鏡掃描(FIB-SEM)等先進(jìn)成像技術(shù)的發(fā)展，為識(shí)別和準(zhǔn)確表征巖石內(nèi)部跨尺度非連續(xù)結(jié)構(gòu)提供了有力手段。筆者基于X射線CT掃描數(shù)據(jù)構(gòu)建了巖石非連續(xù)結(jié)構(gòu)的三維數(shù)字模型及其RVE單元模型。該方法可直觀顯示和定量刻畫巖石內(nèi)部微細(xì)觀尺度的非連續(xù)結(jié)構(gòu)，已廣泛應(yīng)用于巖石力學(xué)研究。與CT成像技術(shù)相比，F(xiàn)IB-SEM是一種破壞性方法，但具有更高的分辨率，可達(dá)1～10 nm，為研究巖石內(nèi)部納米尺度非連續(xù)結(jié)構(gòu)提供了有效方法，但觀測(cè)區(qū)域小和破壞性限制了該技術(shù)在某些情況下的應(yīng)用。MRI在空間和時(shí)間分辨率方面相對(duì)靈活，具有快速識(shí)別內(nèi)部結(jié)構(gòu)變化的優(yōu)勢(shì)，但分辨率相對(duì)較低。將上述方法相結(jié)合可實(shí)現(xiàn)巖石內(nèi)部跨尺度非連續(xù)結(jié)構(gòu)的直觀定量觀測(cè)，有利于在RVE尺度范圍內(nèi)定量刻畫巖石復(fù)雜的內(nèi)部結(jié)構(gòu)并構(gòu)建反映內(nèi)部跨尺度結(jié)構(gòu)的巖石RVE模型。

1.2 非連續(xù)結(jié)構(gòu)的重構(gòu)方法與模型

實(shí)驗(yàn)室及現(xiàn)場(chǎng)原位探測(cè)提供了巖石內(nèi)部非連續(xù)結(jié)構(gòu)的詳細(xì)信息，但考慮到地層巖心鉆探的難度、高精度成像的高成本以及巖石非均質(zhì)非連續(xù)性帶來的物理力學(xué)量對(duì)比分析的困難，發(fā)展一種便捷、經(jīng)濟(jì)和可重復(fù)的巖石非連續(xù)結(jié)構(gòu)的直觀定量描述手段十分重要。近些年，利用有限地質(zhì)探測(cè)數(shù)據(jù)和計(jì)算機(jī)重構(gòu)算法來重構(gòu)巖石三維非連續(xù)結(jié)構(gòu)的研究取得了飛速發(fā)展。隨機(jī)生成、多點(diǎn)統(tǒng)計(jì)信息、模擬退火、基于過程的方法和機(jī)器學(xué)習(xí)等方法被廣泛用于重構(gòu)和刻畫巖石內(nèi)部復(fù)雜非連續(xù)結(jié)構(gòu)。

隨機(jī)生成方法可以簡(jiǎn)單、直接地構(gòu)建出巖石內(nèi)部非連續(xù)結(jié)構(gòu)。例如，基于統(tǒng)計(jì)信息，WANG和PAN提出了一種多參數(shù)隨機(jī)生成-生長(zhǎng)的方法來生成不同類型的多孔結(jié)構(gòu)；筆者改進(jìn)了隨機(jī)球及裂隙填充法，通過控制球及裂隙的大小和孔隙度參數(shù)來生成頁巖孔隙結(jié)構(gòu)；為了研究粗糙裂縫的流體輸送特性，筆者應(yīng)用分形理論直接生成了粗糙裂縫。這些隨機(jī)生成方法通常具有較高的重構(gòu)效率，準(zhǔn)確性取決于重構(gòu)過程中所使用的關(guān)鍵信息。

多點(diǎn)統(tǒng)計(jì)信息法使用模板掃描訓(xùn)練圖像來提取局部多點(diǎn)特征，并將這些特征模式運(yùn)用到新的重構(gòu)圖像中。重構(gòu)過程中考慮了多點(diǎn)統(tǒng)計(jì)信息，可以較好地再現(xiàn)多孔介質(zhì)的幾何形狀和連通性。

模擬退火法是一種被廣泛使用的重構(gòu)方法，它通過不斷調(diào)換代表孔隙相與固體相的數(shù)據(jù)點(diǎn)對(duì)的位置來重構(gòu)模型，這種方法可以將任意類型的目標(biāo)函數(shù)納入重構(gòu)。但由于需要進(jìn)行大量的位置調(diào)換，并生成大量的中間結(jié)構(gòu)，重構(gòu)速度較慢。近年來，為了提高巖石數(shù)字模型的重構(gòu)精度和重構(gòu)效率，各國(guó)學(xué)者在模擬退火算法的基礎(chǔ)上，采用高階相關(guān)函數(shù)、方向相關(guān)函數(shù)、優(yōu)化的模擬退火算法、不同權(quán)值相關(guān)函數(shù)、多尺度重構(gòu)和并行計(jì)算等方法來構(gòu)建巖石非連續(xù)結(jié)構(gòu)三維RVE數(shù)字模型。筆者基于該原理，提出了改進(jìn)的并行模擬退火算法，對(duì)巖石孔隙結(jié)構(gòu)進(jìn)行了高效準(zhǔn)確的三維重構(gòu)(圖1)。

基于沉積巖成巖過程來重建沉積巖孔隙結(jié)構(gòu)的方法由?REN和BAKKE率先提出，該方法考慮了壓實(shí)、顆粒沉積和成巖過程，具有廣泛應(yīng)用前景。

圖1 巖石孔隙結(jié)構(gòu)三維重構(gòu)模型Fig.1 3D reconstruction model of rock pores

機(jī)器學(xué)習(xí)方法通過大量真實(shí)的巖石掃描圖像來建立和訓(xùn)練神經(jīng)網(wǎng)絡(luò)，利用訓(xùn)練好的網(wǎng)絡(luò)來生成巖石三維復(fù)雜結(jié)構(gòu)。這種方法需要大量的數(shù)據(jù)樣本，但神經(jīng)網(wǎng)絡(luò)訓(xùn)練完成后可以在較短的時(shí)間內(nèi)生成具有較好相似度的真實(shí)結(jié)構(gòu)(圖1(b))。受限于神經(jīng)網(wǎng)絡(luò)的能力，目前大多數(shù)重建工作主要局限于64和128的RVE尺度規(guī)模。特別是，由于神經(jīng)網(wǎng)絡(luò)原理的難解釋性和泛化能力差等問題，該方法在重建結(jié)果的準(zhǔn)確性和普適性等方面面臨諸多挑戰(zhàn)，有待進(jìn)一步改進(jìn)。

需要指出的是，盡管研究表明數(shù)字化方法可以經(jīng)濟(jì)有效地獲取和透明顯示巖石內(nèi)部非連續(xù)結(jié)構(gòu)，但目前研究主要集中在利用有限微觀結(jié)構(gòu)信息進(jìn)行巖石RVE模型重構(gòu)。使用低階統(tǒng)計(jì)信息提高了RVE模型重構(gòu)效率，但犧牲了模型精度。隨著計(jì)算機(jī)運(yùn)算能力進(jìn)一步提高，有望采用更多的高階統(tǒng)計(jì)信息來提高RVE模型的重構(gòu)精度，這將是巖石內(nèi)部非連續(xù)結(jié)構(gòu)三維數(shù)字化重構(gòu)與透明表征研究的發(fā)展方向。

2 3D打印模型材料研制

3D打印技術(shù)已廣泛應(yīng)用于巖石力學(xué)研究，該技術(shù)為重復(fù)制備內(nèi)部包含復(fù)雜非連續(xù)結(jié)構(gòu)的巖石物理模型并定量分析非連續(xù)結(jié)構(gòu)對(duì)巖石災(zāi)變力學(xué)行為的影響提供了可能性和新途徑。然而，除內(nèi)部非連續(xù)結(jié)構(gòu)因素外，巖石物性也是影響巖石物理力學(xué)性質(zhì)與災(zāi)變行為的重要因素，研制能夠反映天然巖石物性的3D打印材料是利用3D打印模型開展非連續(xù)巖石災(zāi)變力學(xué)行為面臨的重要挑戰(zhàn)。

針對(duì)此問題，國(guó)內(nèi)外學(xué)者采用不同3D打印技術(shù)，例如，光固化成型技術(shù)、熔積成型法、選擇性燒結(jié)成型法等，研究了硅砂、石膏、光敏和熱固性樹脂等常用打印材料的物理力學(xué)性質(zhì)，探究了這些材料用于制備天然巖石模型的適用性。例如，PERRAS等對(duì)比分析了天然砂巖和3D打印模型的抗壓強(qiáng)度、抗拉強(qiáng)度、剛度、裂紋生成、裂紋損傷閾值和應(yīng)變行為等，發(fā)現(xiàn)Sand-Furan 3D打印試樣的力學(xué)性能與天然砂巖最接近(圖2)；ZHU和ZHOU研究了陶瓷、石膏、PMMA、SR20丙烯酸共聚物和光敏樹脂5種可打印材料在單軸壓縮下的力學(xué)性能和破壞模式，指出采用立體光固化成型法(Stereo Lithography Appearance，SLA)的光敏樹脂材料最適合模擬脆性和堅(jiān)硬巖石，他們通過靜態(tài)和動(dòng)態(tài)巴西劈裂試驗(yàn)發(fā)現(xiàn)光敏樹脂模型表現(xiàn)出與天然巖石相同的破裂行為。GELL等研究了用于制備巖石試件的水泥、樹脂和含砂混料的力學(xué)性質(zhì)，指出3D打印試件可以實(shí)現(xiàn)巖樣內(nèi)部非連續(xù)結(jié)構(gòu)的精細(xì)化復(fù)制，比傳統(tǒng)澆注制樣更具優(yōu)勢(shì)。

除了如何使3D打印模型具有與天然巖石一致的非連續(xù)結(jié)構(gòu)和物性外，如何能夠準(zhǔn)確反映和量化非連續(xù)巖石內(nèi)部應(yīng)力場(chǎng)也是3D打印材料研究面臨的重要挑戰(zhàn)。為了獲得高應(yīng)力敏感性、可用于透明顯示和定量表征巖石內(nèi)部應(yīng)力的3D打印材料，筆者利用美國(guó)Stratasys公司生產(chǎn)的Object Connex 500和J750 3D打印機(jī)，制備了用于打印材料物理力學(xué)性能測(cè)試的標(biāo)準(zhǔn)試件，詳細(xì)研究了3D打印材料VeroClear的組成成分、基本物理力學(xué)性能、打印成型方式及模型后處理的影響、應(yīng)力雙折射效應(yīng)(圖3(a),(b))、光熱曲線(圖3(c))、應(yīng)力凍結(jié)性質(zhì)(圖3(d))以及高低溫下的塑性及脆性性質(zhì)等。其中，Object Connex 500和J750打印機(jī)的最小成型厚度分別為16和14 μm。實(shí)驗(yàn)結(jié)果表明：VeroClear具有良好的應(yīng)力雙折射效應(yīng)和應(yīng)力凍結(jié)性能，這種打印材料在一定條件下具有與煤、砂巖相似的物理力學(xué)性質(zhì)與破壞模式，而且可用于透明顯示和定量表征受載破壞過程中巖石非連續(xù)結(jié)構(gòu)內(nèi)部應(yīng)力場(chǎng)的演化，為透明解析巖石RVE模型內(nèi)部三維應(yīng)力場(chǎng)提供了條件和基礎(chǔ)。

圖2 天然砂巖樣品及3D打印砂巖樣品[35]Fig.2 Natural sandstone samples and 3D printed sandstone samples[35]

圖3 3D打印材料的物理力學(xué)性能Fig.3 Physical and/or mechanical properties of 3D printing materials

需要指出的是，盡管國(guó)內(nèi)外學(xué)者做了大量工作，但從材料物性及物理力學(xué)性能角度來看，3D打印模型材料與天然巖石尚有不小差距，還需進(jìn)一步研究和改性現(xiàn)有的打印材料。根據(jù)目前的研究結(jié)果，改善3D打印材料性能可有以下途徑：① 在分子結(jié)構(gòu)層面設(shè)計(jì)打印材料；② 調(diào)整模型打印方式與參數(shù)；③ 增添不同成分的微觀結(jié)構(gòu)來模擬天然巖石并改進(jìn)打印模型的物理力學(xué)性質(zhì)；④ 優(yōu)化打印模型后處理等，希望應(yīng)用這些技術(shù)使打印模型能夠代替或更好地模擬天然巖石。

3 巖石非連續(xù)結(jié)構(gòu)的透明物理模型

利用現(xiàn)場(chǎng)和實(shí)驗(yàn)室探測(cè)技術(shù)及數(shù)字重構(gòu)方法可以建立直觀定量描述巖石RVE單元內(nèi)跨尺度非連續(xù)結(jié)構(gòu)的數(shù)字透明模型。然而，由于人們對(duì)RVE尺度內(nèi)不連續(xù)結(jié)構(gòu)引發(fā)的應(yīng)力場(chǎng)、變形場(chǎng)和滲流場(chǎng)的變化及其內(nèi)稟物理機(jī)制認(rèn)識(shí)不足，在數(shù)字模型基礎(chǔ)上，將基于連續(xù)性假設(shè)所建立的本構(gòu)模型、強(qiáng)度準(zhǔn)則和失效判據(jù)用于巖石災(zāi)變數(shù)值模擬分析，預(yù)測(cè)結(jié)果往往與工程實(shí)際存在較大偏差。因此，開展物理模型實(shí)驗(yàn)揭示巖石災(zāi)變“黑箱”過程與物理本質(zhì)、建立準(zhǔn)確的體現(xiàn)內(nèi)部非連續(xù)結(jié)構(gòu)效應(yīng)的本構(gòu)模型、強(qiáng)度準(zhǔn)則與失效判據(jù)是解決工程實(shí)際問題的必由之路。然而，如何重復(fù)制備含復(fù)雜內(nèi)部結(jié)構(gòu)的物理模型一直是巖石物理模型實(shí)驗(yàn)研究的一個(gè)難題。

近年來，快速發(fā)展的3D打印技術(shù)為解決上述難題提供了途徑。人們基于數(shù)字模型與3D打印模型材料，采用不同的3D打印技術(shù)制備了與天然巖石物性相似且內(nèi)部結(jié)構(gòu)一致的物理模型。例如，為了研究打印制備的孔隙模型的力學(xué)性質(zhì)和滲流特征，ISHUTOV等利用CT成像和3D打印技術(shù)制備了天然砂巖孔隙結(jié)構(gòu)的物理模型(圖4)。JIANG等打印制備了與天然巖石內(nèi)部節(jié)理分布一致的物理模型，指出3D打印技術(shù)能夠制備含有復(fù)雜孔隙和裂隙結(jié)構(gòu)的巖石模型，可用于巖石力學(xué)性質(zhì)研究。然而，由于采用不透明3D打印材料，上述方法制備的巖石物理模型無法直觀顯示和定量表征巖石內(nèi)部物理場(chǎng)變化與災(zāi)變過程。為了透明顯示和定量解析巖石內(nèi)部非連續(xù)結(jié)構(gòu)及物理場(chǎng)演化與災(zāi)變過程，筆者采用自主研發(fā)的原位三軸加載CT成像系統(tǒng)(美國(guó)、英國(guó)與中國(guó)發(fā)明專利)與分形數(shù)字重構(gòu)方法，應(yīng)用研制的透明3D打印材料，制備了具有與天然巖石一致非連續(xù)裂隙、孔隙結(jié)構(gòu)及相近物理力學(xué)性質(zhì)的巖石透明物理模型(圖5)。

圖4 3D打印孔隙巖石模型[47]Fig.4 3D-printed models of porous rocks[47]

為了模擬賦含非均質(zhì)非均布顆粒與礦物充填的巖石，作者基于多材料打印技術(shù)，采用不同力學(xué)性質(zhì)的材料復(fù)制天然巖石中的不同成分，例如，筆者采用兩種不同打印材料模擬骨料和基體，并按照骨料和孔隙的實(shí)際分布打印制備了混凝土透明模型(圖6)。

值得關(guān)注的是，天然巖石含有大量跨尺度的孔隙、裂縫、顆粒和節(jié)理等微結(jié)構(gòu)，在CT掃描識(shí)別精度范圍內(nèi)，部分微結(jié)構(gòu)在CT圖像中可能僅占有1個(gè)體像素，可采用體像素打印技術(shù)制備這類模型。體像素打印是一種在體像素尺度上使用多種材料進(jìn)行3D打印的技術(shù)，可以對(duì)每個(gè)體像素賦予不同的材料屬性。對(duì)于內(nèi)部組成極其復(fù)雜、跨尺度且難以精確構(gòu)建數(shù)字模型的非連續(xù)結(jié)構(gòu)，可以直接在CT圖像的基礎(chǔ)上，采用體像素打印方式制備透明物理模型，該方法將大幅度提高復(fù)雜物理模型的制備精度。

圖5 巖石透明物理模型[14,42]Fig.5 Transparent models of rocks[14,42]

圖6 混凝土透明模型[52-53]Fig.6 Transparent model of concrete [52-53]

此外，為了模擬天然巖石孔隙或裂隙中的流體介質(zhì)(油、水、氣等)，筆者采用固液兩相打印技術(shù)，實(shí)現(xiàn)了包含流體介質(zhì)的三維RVE模型一體化打印制備。圖7為打印制備的內(nèi)嵌流體的方板圓孔模型和孔隙結(jié)構(gòu)模型，液體頂部留有微小氣泡。結(jié)合體像素3D打印技術(shù)，通過設(shè)計(jì)固液材料的體像素空間分布，可以實(shí)現(xiàn)任意位置處孔隙、裂隙包含不同液體的三維模型制備。

圖7 內(nèi)嵌液體的方板圓孔和孔隙結(jié)構(gòu)的透明3D打印模型Fig.7 3D-printed models of a square plate embedded a circle pore and a porous rock, filled with liquids

4 RVE模型尺度

三維數(shù)字重構(gòu)和3D打印技術(shù)為實(shí)現(xiàn)巖石RVE模型內(nèi)部非連續(xù)結(jié)構(gòu)的透明解析、發(fā)展考慮內(nèi)部非連續(xù)結(jié)構(gòu)特征的巖石力學(xué)理論提供了條件和基礎(chǔ)。利用RVE尺度模型研究巖石的應(yīng)力、應(yīng)變、滲流與災(zāi)變等行為，建立反映內(nèi)部非連續(xù)結(jié)構(gòu)特征的巖石本構(gòu)模型、強(qiáng)度準(zhǔn)則和災(zāi)變判據(jù)，首先需要厘清巖石內(nèi)部非連續(xù)結(jié)構(gòu)特征及其對(duì)巖石力學(xué)行為與災(zāi)變行為的影響，并確定有代表性的RVE模型尺度。BELL和SIEGMUND開展了不同厚度3D打印模型的三點(diǎn)彎曲斷裂實(shí)驗(yàn)，分析了1～5 mm內(nèi)模型強(qiáng)度和斷裂韌性的尺寸效應(yīng)；AKIYAMA等通過四點(diǎn)彎實(shí)驗(yàn)研究了模型尺寸對(duì)基于FDM 技術(shù)的打印模型力學(xué)性能的影響。筆者測(cè)試了采用不同光敏打印材料制備的不同尺寸的含復(fù)雜孔隙和裂隙結(jié)構(gòu)的3D打印模型的抗壓強(qiáng)度、抗拉強(qiáng)度、彈性模量、泊松比、屈服強(qiáng)度、斷裂韌性等，分析了巖石打印模型的尺寸效應(yīng)(圖8),圖8(a)數(shù)據(jù)是根據(jù)至少3個(gè)重復(fù)試件實(shí)測(cè)結(jié)果計(jì)算獲得。實(shí)驗(yàn)結(jié)果表明： 3D打印RVE模型體現(xiàn)出與天然巖石一致的尺寸效應(yīng)，這為定量分析巖石內(nèi)部非連續(xù)結(jié)構(gòu)效應(yīng)及其對(duì)巖石災(zāi)變力學(xué)行為的影響，確定含不同類型非連續(xù)結(jié)構(gòu)的RVE單元尺寸提供了基礎(chǔ)和依據(jù)。

5 巖石內(nèi)部應(yīng)力場(chǎng)的透明解析

5.1 二維模型應(yīng)力場(chǎng)透明解析

為了揭示復(fù)雜非連續(xù)結(jié)構(gòu)對(duì)巖石災(zāi)變行為的影響機(jī)制，筆者采用3D打印技術(shù)制備了天然巖石非連續(xù)結(jié)構(gòu)的二維透明模型，采用光彈性實(shí)驗(yàn)方法，通過改進(jìn)相移和解包裹算法，定量解析了荷載作用下非連續(xù)結(jié)構(gòu)全場(chǎng)應(yīng)力分布特征與演化規(guī)律，為認(rèn)知和解析巖石災(zāi)變行為的內(nèi)在物理機(jī)制創(chuàng)造了條件。例如，筆者基于巖石裂隙結(jié)構(gòu)的CT圖像，制備了天然裂隙的二維透明RVE模型，發(fā)展改進(jìn)了相移和解包裹算法，獲得了外部荷載作用下裂隙結(jié)構(gòu)全場(chǎng)剪應(yīng)力和主應(yīng)力差(圖9)。此外，筆者還模擬巖石粗糙斷層結(jié)構(gòu)，打印制備了粗糙斷層的二維透明RVE模型，設(shè)計(jì)了斷層剪切變形的光彈性實(shí)驗(yàn)，開發(fā)了非連續(xù)結(jié)構(gòu)光彈性應(yīng)力場(chǎng)的定量提取算法，定量解析了斷層滑移時(shí)斷層區(qū)域全場(chǎng)主應(yīng)力差和剪應(yīng)力分布與演化規(guī)律(圖10)，建立了粗糙面起伏角與斷層周邊應(yīng)力場(chǎng)的定量關(guān)系。采用類似方法，作者打印制備了隧道襯砌與圍巖結(jié)構(gòu)的二維平面模型，提出了圍巖襯砌結(jié)構(gòu)全場(chǎng)應(yīng)力的實(shí)驗(yàn)提取與定量解析方法，獲得了外載作用下圍巖襯砌結(jié)構(gòu)全場(chǎng)主應(yīng)力差與剪應(yīng)力分布及演化規(guī)律(圖11)，為定量表征襯砌結(jié)構(gòu)及管片界面應(yīng)力分布、驗(yàn)證數(shù)值模擬結(jié)果提供了方法和途徑。類似的，MISSERONI等通過可視化內(nèi)含剛性顆粒模型周邊光彈性條紋的分布，分析了剛性顆粒模型周邊的應(yīng)力集中和應(yīng)力場(chǎng)分布特征；DANIELS和HAYMAN通過顆粒光彈模型，解析了不同節(jié)理結(jié)構(gòu)特征對(duì)巖石內(nèi)部應(yīng)力場(chǎng)的影響；BIGONI和NOSELLI通過建立墻體結(jié)構(gòu)的光彈模型，分析了承載梁壓縮引起的墻體內(nèi)部壓應(yīng)力分布的特征；GHISLAIN 等建立了地質(zhì)斷層的光彈模型，利用光彈性條紋和主應(yīng)力跡線，分析了斷層周圍應(yīng)力場(chǎng)(圖12)；SOLIVA 等和LUNINA等利用地質(zhì)斷層的光彈性模型，定量解析了斷層周邊應(yīng)力場(chǎng)的分布與演化。

圖8 巖石3D打印模型尺寸效應(yīng)的實(shí)驗(yàn)結(jié)果Fig.8 Experimental investigation on the size effects of 3D-printed rock models

圖9 巖石裂隙結(jié)構(gòu)2D透明模型及不同軸向載荷作用下裂隙區(qū)域應(yīng)力場(chǎng)演化的解析結(jié)果[62]Fig.9 2D transparent model of rock fractures and the evolution of stress fields under different uniaxial compressive loads[62]

圖10 巖石粗糙斷層2D透明模型及剪切位移變化時(shí)斷層周邊應(yīng)力場(chǎng)分布演化的解析結(jié)果[64]Fig.10 2D transparent model of a rough fault and the distribution and evolution of stress fields under different shear deformations[64]

圖11 襯砌結(jié)構(gòu)應(yīng)力場(chǎng)的透明解析結(jié)果[65]Fig.11 Transparentized solutions to the stress fields of lining structures[65]

5.2 三維模型應(yīng)力場(chǎng)透明解析

巖石非連續(xù)結(jié)構(gòu)二維模型應(yīng)力透明解析為認(rèn)知和量化非連續(xù)結(jié)構(gòu)對(duì)巖石災(zāi)變行為的影響、構(gòu)建非連續(xù)巖石力學(xué)理論提供了途徑和支撐。然而，天然巖石包含復(fù)雜的三維結(jié)構(gòu)，透明顯示和定量解析外載作用下巖石內(nèi)部三維應(yīng)力場(chǎng)分布及演化對(duì)于闡釋巖石結(jié)構(gòu)整體災(zāi)變機(jī)理和預(yù)測(cè)工程災(zāi)害發(fā)生的位置與量級(jí)具有重要意義。

圖12 斷層周邊應(yīng)力場(chǎng)定量解析結(jié)果[70]Fig.12 Quantitative solutions to the stress fields near faults

光彈性應(yīng)力凍結(jié)是一種獲取和分析三維連續(xù)均質(zhì)模型內(nèi)部應(yīng)力場(chǎng)的有效方法，該方法是在溫箱中對(duì)光彈性模型施加固定荷載，經(jīng)歷溫度升降循環(huán)，光彈性材料在溫升過程中分別出現(xiàn)液態(tài)相和彈性相，降溫時(shí)液態(tài)相固化，模型內(nèi)部應(yīng)力被“凍結(jié)”下來。對(duì)“凍結(jié)”后模型進(jìn)行不同位置的切片，通過觀測(cè)切片的光彈性條紋來解析三維模型內(nèi)部全場(chǎng)應(yīng)力分布。然而，天然巖石結(jié)構(gòu)復(fù)雜，傳統(tǒng)的光彈性應(yīng)力凍結(jié)法難以制備包含眾多跨尺度孔隙、裂隙、顆粒等不規(guī)則結(jié)構(gòu)的三維模型，加之缺乏有效的針對(duì)非連續(xù)結(jié)構(gòu)應(yīng)力場(chǎng)的提取算法，實(shí)現(xiàn)非連續(xù)巖石內(nèi)部三維應(yīng)力場(chǎng)定量解析一直是個(gè)難題。

針對(duì)上述問題，筆者采用高精度CT成像、3D數(shù)字重構(gòu)和3D打印技術(shù)，在二維模型應(yīng)力解析方法的基礎(chǔ)上，通過研制高應(yīng)力敏感3D打印材料、改進(jìn)光彈性應(yīng)力凍結(jié)法和發(fā)展針對(duì)復(fù)雜結(jié)構(gòu)的應(yīng)力場(chǎng)提取算法等措施，提出了非連續(xù)巖石內(nèi)部三維應(yīng)力場(chǎng)的透明定量解析方法，為實(shí)現(xiàn)巖石RVE模型內(nèi)部復(fù)雜三維應(yīng)力場(chǎng)的透明解析、構(gòu)建基于RVE內(nèi)部非連續(xù)結(jié)構(gòu)效應(yīng)的巖石力學(xué)理論提供了實(shí)驗(yàn)基礎(chǔ)與理論支撐。例如，圖13展示了裂隙煤巖三維透明模型內(nèi)部不同位置處的應(yīng)力條紋；圖14給出了孔隙巖石三維透明模型內(nèi)部不同位置處的應(yīng)力條紋與定量解析結(jié)果；圖15為非均質(zhì)砂礫巖三維透明模型內(nèi)部不同位置處的應(yīng)力分布解析結(jié)果，該結(jié)果直觀顯示和量化了相鄰顆粒間的高應(yīng)力分布特征，為闡釋非均質(zhì)砂礫巖變形破壞的內(nèi)部非均質(zhì)非均布顆粒效應(yīng)及應(yīng)力場(chǎng)演化提供了實(shí)驗(yàn)證據(jù)及依據(jù)；圖16直觀顯示和定量表征了軸向壓縮荷載作用下不同荷載階段裂隙煤巖內(nèi)部三維應(yīng)力場(chǎng)的分布與演化特征。

圖13 裂隙煤巖三維透明模型內(nèi)部不同位置處的應(yīng)力分布特征[14]Fig.13 Stress distribution at various positions inside a 3D transparent model of fractured coals[14]

為了揭示開采過程誘發(fā)巷道圍巖災(zāi)變機(jī)理和超前預(yù)測(cè)采動(dòng)圍巖破壞的位置和量級(jí)，筆者基于礦區(qū)現(xiàn)場(chǎng)三維地質(zhì)勘探信息，利用3D打印技術(shù)，制備了回采巷道和回采工作面的三維縮尺透明物理模型(圖17)，模型中的煤層及覆巖采用與其物理力學(xué)性質(zhì)相近的3D打印模型材料制備。采用上述應(yīng)力解析方法，透明顯示和定量表征了不同掘進(jìn)位置處回采巷道圍巖的應(yīng)力分布大小與演化規(guī)律，為超前預(yù)警不同開采方式可能誘發(fā)的巷道與工作面圍巖災(zāi)變提供了研究與定量分析手段。

圖14 孔隙巖石模型內(nèi)部應(yīng)力場(chǎng)的透明解析[42]Fig.14 Transparentized solutions to the stress fields inside porous rock models [42]

圖15 非均質(zhì)砂礫巖模型內(nèi)部應(yīng)力場(chǎng)的透明解析[53]Fig.15 Transparentized solutions to the stress fields inside heterogeneous glutenite models[53]

5.3 動(dòng)態(tài)應(yīng)力場(chǎng)的透明解析

準(zhǔn)確獲取和量化巖石全場(chǎng)應(yīng)力連續(xù)動(dòng)態(tài)演化與非連續(xù)結(jié)構(gòu)效應(yīng)對(duì)于揭示礦山巖爆、瓦斯突出、地震等動(dòng)力災(zāi)害的內(nèi)稟物理機(jī)理至關(guān)重要。國(guó)內(nèi)外眾多學(xué)者發(fā)展了動(dòng)態(tài)光彈性實(shí)驗(yàn)方法來獲取和分析巖石應(yīng)力場(chǎng)的動(dòng)態(tài)演化特征，例如，相移法、RGB法、載荷步法、載波法和焦散線法等。但由于巖石非連續(xù)結(jié)構(gòu)模型制備困難以及動(dòng)態(tài)應(yīng)力場(chǎng)的連續(xù)高速變化特征，傳統(tǒng)的巖石動(dòng)態(tài)應(yīng)力場(chǎng)表征主要集中于簡(jiǎn)單結(jié)構(gòu)模型以及應(yīng)力場(chǎng)演化的定性或半定量分析，難以直觀顯示和定量解析復(fù)雜巖石結(jié)構(gòu)的動(dòng)態(tài)應(yīng)力場(chǎng)演化特征。

圖16 軸向壓縮荷載作用下不同階段裂隙煤巖內(nèi)部應(yīng)力場(chǎng)的分布與演化特征Fig.16 Distribution and evolution of the stress fields inside a fractured coal under different compressive loads

圖17 回采巷道圍巖應(yīng)力場(chǎng)的定量解析與透明推演[79]Fig.17 Quantitative solution and transparentized interpretation of the stress field in mining roadways[79]

針對(duì)上述問題，筆者利用3D打印技術(shù)制備了內(nèi)嵌梭形裂紋的透明圓盤模型，通過對(duì)徑壓縮光彈性實(shí)驗(yàn)以及改進(jìn)動(dòng)態(tài)光彈條紋提取算法，獲取并定量解析了裂紋尖端應(yīng)力場(chǎng)特征及其與裂紋擴(kuò)展速度之間的對(duì)應(yīng)關(guān)系；對(duì)于復(fù)雜非連續(xù)結(jié)構(gòu)內(nèi)部動(dòng)態(tài)應(yīng)力場(chǎng)定量解析難題，作者提出了一種定量解析連續(xù)載荷作用下復(fù)雜非連續(xù)結(jié)構(gòu)模型內(nèi)部全場(chǎng)應(yīng)力演化的方法，獲得了連續(xù)加載下復(fù)雜孔隙模型內(nèi)部的主應(yīng)力差分布與演化規(guī)律(圖18)，該方法適用于定量解析沖擊或爆炸載荷作用下復(fù)雜巖石結(jié)構(gòu)模型內(nèi)部應(yīng)力場(chǎng)，但前提是能夠清晰地獲取快速連續(xù)變化的動(dòng)態(tài)光彈條紋。

由于存在大量不規(guī)則天然裂隙，巖石內(nèi)部裂隙三維擴(kuò)展行為及起裂判據(jù)一直是巖石斷裂與災(zāi)變理論研究的熱點(diǎn)和難點(diǎn)問題。巖石三維裂紋擴(kuò)展過程中的應(yīng)力場(chǎng)特征對(duì)于理解和定量分析裂紋起裂、擴(kuò)展、轉(zhuǎn)向和匯聚行為、建立裂紋非穩(wěn)定擴(kuò)展失效判據(jù)至關(guān)重要。然而，傳統(tǒng)實(shí)驗(yàn)方法在制備三維裂紋模型以及獲取裂紋三維擴(kuò)展過程中的應(yīng)力變化等方面存在困難，三維裂紋連續(xù)擴(kuò)展過程中的應(yīng)力場(chǎng)透明顯示與定量描述尚未得到解決。針對(duì)此問題，筆者采用3D打印技術(shù)和所研制的應(yīng)力敏感3D打印材料，制備了含內(nèi)嵌非貫穿型裂紋和貫穿型裂紋的透明模型，通過模型受壓的光彈性應(yīng)力凍結(jié)實(shí)驗(yàn)，利用改進(jìn)的光彈條紋分析技術(shù)和應(yīng)力場(chǎng)提取算法，獲得并定量解析了裂紋三維擴(kuò)展過程中模型應(yīng)力全場(chǎng)特別是裂紋尖端區(qū)域主應(yīng)力差和剪應(yīng)力的分布及演化特征(圖19)。

圖18 載荷連續(xù)作用下不同加載階段復(fù)雜孔隙結(jié)構(gòu)模型內(nèi)部主應(yīng)力差全場(chǎng)分布與演化規(guī)律[54]Fig.18 Distribution and evolution of the principal stress difference inside a porous rock model under different continuous compressive loads[54]

圖19 三維內(nèi)嵌裂隙和貫穿裂隙擴(kuò)展過程中內(nèi)部不同位置處主應(yīng)力差的透明解析[101-102]Fig.19 Transparentized solutions to the principal stress differences at different positions inside the 3D models containing embedded and penetrated fractures[101-102]

6 巖石內(nèi)部應(yīng)變場(chǎng)的透明解析

直觀顯示和定量表征巖石內(nèi)部變形破壞程度和集中分布區(qū)域?qū)τ诶斫鈳r石整體變形破壞機(jī)理具有重要意義。然而，目前常用的數(shù)字圖像相關(guān)法(Digital Image Correlation,DIC)和表面三維DIC技術(shù)只能測(cè)量巖石表面和離面的變形，難以用于量測(cè)巖石內(nèi)部變形，迫切需要發(fā)展內(nèi)部三維變形場(chǎng)的直觀定量表征的方法。借助三維成像技術(shù)，如微焦點(diǎn)CT、共聚焦顯微鏡、磁共振成像等技術(shù)，可以獲得巖石樣品的三維圖像，基于此，BAY等提出了數(shù)字體圖像相關(guān)法(Digital Volume Correlation,DVC)來量測(cè)受壓骨骼內(nèi)部應(yīng)變場(chǎng)。近年來，DVC方法得到迅速發(fā)展和應(yīng)用。但如何提高內(nèi)部微小應(yīng)變的測(cè)量與計(jì)算精度一直是DVC法面臨的主要挑戰(zhàn)。根據(jù)配準(zhǔn)算法的不同，將DVC方法分為基于子集的局部DVC (Local DVC， L-DVC)和基于有限元的全局DVC (Global DVC,G-DVC)兩種主要方法。相比之下，L-DVC計(jì)算效率高，在計(jì)算參數(shù)相同的情況下，G-DVC的精度略高于L-DVC，尤其是當(dāng)子集尺寸較小時(shí)。針對(duì)巖石內(nèi)部非連續(xù)結(jié)構(gòu)特征，筆者提出了多尺度子塊體搜索和亞體素平移方法MS-DVSP(Multi-scale and Shifting Digital Volumetric Speckle Photography)，顯著提高了復(fù)雜結(jié)構(gòu)內(nèi)部三維變形的測(cè)量精度，透明顯示和定量解析了荷載作用下煤巖(圖20)、紅砂巖和混凝土等內(nèi)部三維變形的分布與演化規(guī)律，為構(gòu)建基于內(nèi)部應(yīng)力場(chǎng)和應(yīng)變場(chǎng)透明解析的非連續(xù)巖石本構(gòu)關(guān)系與災(zāi)變準(zhǔn)則提供了途徑和基礎(chǔ)。

值得關(guān)注的是，DVC法雖然能夠應(yīng)用于量測(cè)巖石內(nèi)部變形場(chǎng)，但天然巖石結(jié)構(gòu)復(fù)雜，每一個(gè)用于實(shí)驗(yàn)量測(cè)的巖石樣品結(jié)構(gòu)差異性很大，難以保證被測(cè)巖石RVE結(jié)構(gòu)的一致性以及所獲得的應(yīng)變場(chǎng)具有足夠代表性。同時(shí)， 由于DVC法要求被測(cè)巖石RVE內(nèi)部必須有滿足計(jì)算要求的高質(zhì)量散斑，而很多類別巖石內(nèi)部并不存在這樣的散斑，這使得DVC法難以直接應(yīng)用于測(cè)量此類巖石的內(nèi)部變形場(chǎng)。為解決此難題，筆者采用多材料3D打印技術(shù)制備了內(nèi)置散斑載體與裂隙結(jié)構(gòu)的三維透明模型(圖21)，發(fā)展了三維散斑數(shù)字圖像相關(guān)法，獲得并定量分析了含裂隙模型的內(nèi)部變形場(chǎng)，為透明定量解析三維非連續(xù)結(jié)構(gòu)內(nèi)部變形分布演化特征與變形破壞機(jī)理提供了新方法和新途徑。

圖20 不同CO2吸附時(shí)間下煤巖內(nèi)部的體積應(yīng)變[116]Fig.20 Volumetric strain inside a coal at different CO2 sorption time [116]

圖21 內(nèi)嵌裂紋模型位移場(chǎng)的透明解析Fig.21 Transparentized solution to the displacement field of a 3D fractured model

7 巖石內(nèi)部滲流場(chǎng)的透明解析

準(zhǔn)確認(rèn)識(shí)和定量表征巖石內(nèi)流體輸運(yùn)行為對(duì)于研究和解決煤與瓦斯共采、油氣資源開發(fā)、土壤污染物運(yùn)移、CO地質(zhì)封存及核廢料處置等一系列工程問題具有重要意義。然而，目前常用的RVE尺度巖心驅(qū)替實(shí)驗(yàn)重點(diǎn)關(guān)注流體宏觀輸運(yùn)性質(zhì)，如滲透率、相滲曲線和毛細(xì)壓力曲線等特征，缺乏對(duì)這些宏觀輸運(yùn)性質(zhì)演變起控制作用的孔尺度流體運(yùn)移機(jī)理的深入分析和準(zhǔn)確描述，導(dǎo)致宏觀輸運(yùn)性質(zhì)演化規(guī)律難以預(yù)測(cè)。同時(shí)，由于巖心驅(qū)替實(shí)驗(yàn)中巖石內(nèi)部流體流動(dòng)過程不可見，且孔隙結(jié)構(gòu)存在顯著非均質(zhì)性，全場(chǎng)流動(dòng)模式極為復(fù)雜，僅憑宏觀認(rèn)識(shí)難以準(zhǔn)確反映和描述真實(shí)流動(dòng)規(guī)律。常用的微觀模型實(shí)驗(yàn)在模型制備及施加溫度壓力條件等方面受限，且模型一般為小尺度二維模型，導(dǎo)致微觀模型實(shí)驗(yàn)結(jié)果難以直接應(yīng)用于實(shí)際工程。此外，常用的流體動(dòng)力學(xué)計(jì)算方法(CFD)在處理復(fù)雜幾何結(jié)構(gòu)及多相滲流行為受到較大限制，且受限于計(jì)算資源，模擬計(jì)算難以準(zhǔn)確反映真實(shí)地層非均質(zhì)性孔隙結(jié)構(gòu)內(nèi)的多相輸運(yùn)行為。針對(duì)常規(guī)巖心實(shí)驗(yàn)、微觀模型實(shí)驗(yàn)及傳統(tǒng)數(shù)值分析存在的問題，筆者借助3D打印技術(shù)、自主設(shè)計(jì)研發(fā)的原位三軸加載CT成像系統(tǒng)(圖22)以及改進(jìn)的并行格子玻爾茲曼方法(LBM)，透明顯示和定量解析了復(fù)雜巖石復(fù)雜結(jié)構(gòu)內(nèi)流體輸運(yùn)行為及非均質(zhì)孔隙結(jié)構(gòu)中的非混相界面不穩(wěn)定性，探討了復(fù)雜孔隙系統(tǒng)內(nèi)非混相水油驅(qū)替動(dòng)力學(xué)特征、優(yōu)勢(shì)路徑及其對(duì)孔隙結(jié)構(gòu)幾何參數(shù)、流體黏度等物理參數(shù)的敏感性，通過直觀展示流場(chǎng)、識(shí)別和分析孔隙結(jié)構(gòu)內(nèi)流體流動(dòng)路徑，闡釋了巖石復(fù)雜結(jié)構(gòu)內(nèi)流體輸運(yùn)控制機(jī)理。這些研究為揭示巖石復(fù)雜結(jié)構(gòu)內(nèi)流體滲流機(jī)制及定量分析流體宏觀輸運(yùn)特性、解決實(shí)際工程問題提供了基礎(chǔ)和手段。

圖22 自主設(shè)計(jì)研發(fā)的原位三軸加載CT成像系統(tǒng)Fig.22 Self-developed in-situ triaxial loading and CT imaging system

針對(duì)長(zhǎng)期注水開發(fā)老油田含水率大幅上升、儲(chǔ)層水竄導(dǎo)致采收率下降、注水無效循環(huán)以及廢水處理成本升高等工程難題，筆者采用高精度微米CT成像方法，獲得了真實(shí)儲(chǔ)層不同水驅(qū)條帶代表性三維孔隙結(jié)構(gòu)，依據(jù)天然巖石真實(shí)結(jié)構(gòu)信息，采用3D打印技術(shù)制備了巖石孔隙結(jié)構(gòu)的透明模型，利用透明模型的多相驅(qū)替實(shí)驗(yàn)，直觀展示和定量解析了油、水、堵劑、聚合物、表面活性劑驅(qū)替過程(圖23)，揭示了條帶結(jié)構(gòu)滲透率級(jí)差及注入黏度差異是高滲透條帶形成水竄通道主要原因，堵劑運(yùn)移路徑及原位成膠是保證后續(xù)注入轉(zhuǎn)向進(jìn)入驅(qū)替高含油區(qū)域、進(jìn)而提高剩余油采收率的重要因素。

油氣注液驅(qū)替開采會(huì)造成儲(chǔ)層地應(yīng)力重分布，引發(fā)巖石結(jié)構(gòu)變形，進(jìn)而改變巖石內(nèi)流體的驅(qū)替特性。然而，儲(chǔ)層巖石深埋地下，驅(qū)替開采引發(fā)的巖石孔隙結(jié)構(gòu)三維形變及驅(qū)替行為“看不見、摸不著”，難以直觀顯示和定量描述，同時(shí)也缺乏三維孔隙變形過程中孔尺度非混相驅(qū)替動(dòng)力學(xué)行為與機(jī)理的科學(xué)認(rèn)知與理論解析模型，傳統(tǒng)實(shí)驗(yàn)方法無法觀測(cè)孔隙變形時(shí)的三維非混相驅(qū)替行為演化，難以揭示驅(qū)替過程的主控因素。為解決該問題，筆者在天然巖石三維孔隙結(jié)構(gòu)CT圖像的基礎(chǔ)上，采用3D打印技術(shù)，制備了與天然巖石相同孔隙結(jié)構(gòu)的三維物理模型(圖24，其中，為歸一化的圍壓；為驅(qū)替效率)，采用自主研發(fā)的原位三軸加載CT成像系統(tǒng)，開展了不同圍壓條件下三維孔隙結(jié)構(gòu)變形與變形孔隙結(jié)構(gòu)內(nèi)三維油水驅(qū)替演化的原位CT掃描實(shí)驗(yàn)，透明顯示并定量分析了不同圍壓條件下油水驅(qū)替過程中驅(qū)替前緣運(yùn)移行為及指進(jìn)現(xiàn)象、優(yōu)勢(shì)路徑以及殘余油空間分布的動(dòng)態(tài)演化(圖23)，建立了有效應(yīng)力、孔隙變形特征以及油水驅(qū)替效率之間的定量關(guān)系，為科學(xué)認(rèn)知和定量描述油氣資源開采引起的儲(chǔ)層巖石孔隙結(jié)構(gòu)形變與水油驅(qū)替機(jī)理提供了新手段、新方法，取得了新發(fā)現(xiàn)。

直接觀測(cè)與定量表征煤層開采過程中煤巖內(nèi)部三維裂隙網(wǎng)絡(luò)演化以及裂隙網(wǎng)絡(luò)內(nèi)部流體流動(dòng)是預(yù)測(cè)導(dǎo)水裂隙帶演化和防治礦井水污染的重要基礎(chǔ)。然而，采動(dòng)應(yīng)力作用下煤巖內(nèi)部三維裂隙網(wǎng)絡(luò)發(fā)育演化以及流體流動(dòng)難以直接觀測(cè)，現(xiàn)有理論模型和數(shù)值模擬方法難以計(jì)算和預(yù)測(cè)這一復(fù)雜過程。近年來，國(guó)內(nèi)外學(xué)者嘗試使用CT成像技術(shù)來研究巖石內(nèi)部三維裂隙演化與流體流動(dòng)規(guī)律，但受限于實(shí)驗(yàn)設(shè)備和測(cè)試方法，目前研究大多采用對(duì)變形破壞后的巖石樣品進(jìn)行CT掃描來觀測(cè)巖石內(nèi)部的三維裂隙網(wǎng)絡(luò)及流體分布。然而，這種方法難以觀測(cè)樣品在承載條件下的內(nèi)部結(jié)構(gòu)變形與流動(dòng)行為，相較于原位應(yīng)力狀態(tài)下三維裂隙結(jié)構(gòu)和流體流動(dòng)行為會(huì)產(chǎn)生較大誤差。針對(duì)上述挑戰(zhàn)，筆者利用所研發(fā)的原位三軸加載CT成像系統(tǒng)，開展了不同開采應(yīng)力路徑下煤巖三維裂隙網(wǎng)絡(luò)演化與水流特征的原位CT觀測(cè)實(shí)驗(yàn)，直觀展現(xiàn)并定量表征了不同應(yīng)力階段煤巖內(nèi)部三維裂隙網(wǎng)絡(luò)的拓?fù)浣Y(jié)構(gòu)、連通性和裂隙生長(zhǎng)閉合，并分析了其對(duì)裂隙煤巖滲透率的影響(圖25，為連通性系數(shù))，該研究提供了一種直接觀測(cè)和定量表征采動(dòng)應(yīng)力下煤巖三維裂隙網(wǎng)絡(luò)演化和滲透率改變的實(shí)驗(yàn)方法，為預(yù)測(cè)煤礦開采過程中導(dǎo)水裂隙帶的分布與演化以及礦井水污染防治提供了基礎(chǔ)和參考。

圖23 高中低滲非均質(zhì)孔隙結(jié)構(gòu)內(nèi)調(diào)剖堵水透明可視化實(shí)驗(yàn)[126]Fig.23 Visualization of water flooding,polymer gel injection,antidilution polymer flooding,and surfactant agent flooding processes in a 3D printed heterogeneous porous model[126]

值得關(guān)注的是，煤炭資源開發(fā)過程誘發(fā)的巖石災(zāi)變常常是固體變形與流體流動(dòng)耦合作用的結(jié)果，如何透明顯示和定量解析巖石災(zāi)變中的流固耦合效應(yīng)一直是巖石力學(xué)研究的前沿課題。基于作者建立的固體應(yīng)力場(chǎng)透明解析方法以及流體雙折射法，筆者制備了內(nèi)含流體介質(zhì)的巖石復(fù)雜結(jié)構(gòu)的3D打印模型(圖6(b))，提出災(zāi)變過程中巖石固體與流體應(yīng)力場(chǎng)耦合作用的透明解析方法。該方法同時(shí)捕捉外載作用下巖石模型固體和流體內(nèi)部的應(yīng)力光彈條紋(圖26)，利用作者改進(jìn)的光彈條紋解析方法，提取固體和流體介質(zhì)的應(yīng)力場(chǎng)，實(shí)現(xiàn)巖石變形過程中內(nèi)部固體與流體應(yīng)力場(chǎng)耦合作用及致災(zāi)機(jī)理的透明解析。其中流體雙折射法(Flow-Induced Birefringence，F(xiàn)IB)是指利用具有雙折射效應(yīng)的流體來直觀顯示流場(chǎng)應(yīng)力分布的方法，例如，PHILIPPOFF，PRADOS和PEEBLES等學(xué)者研究發(fā)現(xiàn)了流體的應(yīng)力雙折射效應(yīng)；MULLER和VERGNES研究了流體應(yīng)力雙折射條紋與流體內(nèi)部應(yīng)力場(chǎng)的對(duì)應(yīng)關(guān)系；SENGUPTA等利用流體應(yīng)力雙折射條紋研究了圓柱繞流條件下流體內(nèi)部的應(yīng)力場(chǎng)演化(圖27)；AUHL等借助流體雙折射條紋分析了不同應(yīng)變率條件下流體流動(dòng)過程中的應(yīng)力場(chǎng)，并與數(shù)值計(jì)算結(jié)果進(jìn)行了對(duì)比(圖28,其中，為應(yīng)變率；為時(shí)間)。上述研究為透明解析巖石內(nèi)部流體介質(zhì)與固體結(jié)構(gòu)內(nèi)部應(yīng)力場(chǎng)耦合效應(yīng)提供了途徑和方法。

圖24 不同圍壓條件下水驅(qū)路徑與殘余油空間分布[131]Fig.24 Water flow paths and 3D oil trapping under various confining pressures[131]

圖25 不同開采應(yīng)力路徑下煤巖內(nèi)部三維裂隙演化特征[135]Fig.25 Evolution of 3D fractures inside a coal under various mining-induced stresses[135]

圖26 豎向壓縮時(shí)內(nèi)含流體介質(zhì)的三維孔隙結(jié)構(gòu)固體與流體應(yīng)力條紋Fig.26 Stess patterns of the solid and liquid in a 3D porous model under vertical compressive loads

圖27 流體繞過圓柱障礙物時(shí)的流動(dòng)情況[139]Fig.27 Fluid birefringence fringe during fluid flow through one cylinder[139]

圖28 不同應(yīng)變率下流體內(nèi)部應(yīng)力實(shí)測(cè)值與數(shù)值結(jié)果的對(duì)比[140]Fig.28 Comparison of fluid birefringence fringes and numerical solutions under different strain rates[140]

8 結(jié)論與展望

(1)三維數(shù)字重構(gòu)方法可以經(jīng)濟(jì)有效地獲取和透明顯示巖石內(nèi)部非連續(xù)結(jié)構(gòu)，但目前研究主要集中在利用有限微觀結(jié)構(gòu)信息進(jìn)行巖石RVE模型重構(gòu)，在重構(gòu)模型尺度、重構(gòu)效率和模型精度方面仍需要進(jìn)一步提升。

(2)3D打印技術(shù)為制備巖石復(fù)雜非連續(xù)結(jié)構(gòu)物理模型并量化非連續(xù)結(jié)構(gòu)對(duì)巖石災(zāi)變力學(xué)行為的影響提供了新途徑。但研制反映天然巖石物性并能夠透明量化巖石多物理場(chǎng)行為的3D打印材料仍面臨諸多挑戰(zhàn)，尚需從多方面進(jìn)一步研究和改性現(xiàn)有的打印材料。

(3)基于具有自主知識(shí)產(chǎn)權(quán)的原位三軸加載CT成像系統(tǒng)、三維數(shù)字重構(gòu)方法及研制的透明3D打印材料，制備了具有與天然巖石一致非連續(xù)裂隙、孔隙結(jié)構(gòu)及相近物理力學(xué)性質(zhì)的巖石透明物理模型，為通過RVE物理模型揭示并量化巖石災(zāi)變的內(nèi)部非連續(xù)結(jié)構(gòu)及多物理場(chǎng)耦合效應(yīng)創(chuàng)造了條件。

(4)初步形成了巖石非連續(xù)結(jié)構(gòu)應(yīng)力場(chǎng)、應(yīng)變場(chǎng)、滲流場(chǎng)等多物理場(chǎng)定量透明解析理論與方法，為破解開采誘發(fā)的巖石災(zāi)變“黑箱”機(jī)理、建立基于RVE尺度多物理場(chǎng)透明解析的非連續(xù)巖石力學(xué)與災(zāi)變分析理論、實(shí)現(xiàn)開采災(zāi)害的透明推演、超前預(yù)警和科學(xué)防控提供了途徑和研究基礎(chǔ)。

需要指出的是，巖石災(zāi)變內(nèi)部非連續(xù)結(jié)構(gòu)及多物理場(chǎng)效應(yīng)的透明解析與透明推演研究尚屬剛剛起步，當(dāng)前研究面臨諸多挑戰(zhàn)。例如，除需進(jìn)一步研發(fā)反映天然巖石物性并能夠透明量化巖石多物理場(chǎng)行為的3D打印材料外，在透明解析巖石內(nèi)部塑性變形與應(yīng)力場(chǎng)特征、建立內(nèi)部塑性變形解析量與宏觀塑性力學(xué)響應(yīng)之間的定量關(guān)系、透明解析多物理場(chǎng)之間的耦合效應(yīng)、分析和量化多相(固、液、氣)多物理場(chǎng)之間的相互作用仍存在困難，迫切需要通過多學(xué)科交叉融合開展更深入的研究。

[1] 2050世界與中國(guó)能源展望(2020版)[R]. 中國(guó)石油經(jīng)濟(jì)技術(shù)研究院,2020.

[2] 中國(guó)政府網(wǎng). 中央經(jīng)濟(jì)工作會(huì)議[EB/OL]. http://www. gov. cn/xinwen/2021-12/10/content_5659796. html,2021-12-10.

[3] 新華社. 習(xí)近平在第七十五屆聯(lián)合國(guó)大會(huì)一般性辯論上的講話[EB/OL]. http://www. xinhuanet. com/politics/leaders/2020-09/22/c_1126527652. htm，2020-09-22.

[4] 潘一山,代連朋,李國(guó)臻,等. 煤礦沖擊地壓與冒頂復(fù)合災(zāi)害研究[J]. 煤炭學(xué)報(bào),2021,46(1):112-122.

PAN Yishan,DAI Lianpeng,LI Guozhen,et al. Study on compound disaster of rock burst and roof falling in coal mines[J]. Journal of China Coal Society,2021,46(1):112-122.

[5] 袁亮. 深部采動(dòng)響應(yīng)與災(zāi)害防控研究進(jìn)展[J]. 煤炭學(xué)報(bào),2021,46(3):716-725.

YUAN Liang. Research progress of mining response and disaster prevention and control in deep coal mines[J]. Journal of China Coal Society,2021,46(3):716-725.

[6] 齊慶新,潘一山,李海濤,等. 煤礦深部開采煤巖動(dòng)力災(zāi)害防控理論基礎(chǔ)與關(guān)鍵技術(shù)[J]. 煤炭學(xué)報(bào),2020,45(5):1567-1584.

QI Qingxin,PAN Yishan,LI Haitao,et al. Theoretical basis and key technology of prevention and control of coal-rock dynamic disasters in deep coal mining[J]. Journal of China Coal Society,2020,45(5):1567-1584.

[7] 李鐵,蔡美峰,蔡明. 采礦誘發(fā)地震分類的探討[J]. 巖石力學(xué)與工程學(xué)報(bào),2006,25(S2):3679-3686.

LI Tie,CAI Meifeng,CAI Ming. A discussion on classification of mining-induced seismicity[J]. Chinese Journal of Mechanics and Engineering,2006,25(S2):3679-3686.

[8] DEKOCK Tim,BOONE Marijn A,DE Schryver Thomas,et al. A pore-scale study of fracture dynamics in rock using X-ray micro-CT under ambient freeze-thaw cycling[J]. 2015,49(5):2867-2874.

[9] JIANG Zeyun,VAN Dijke Marinusizaak Jan,GEIGER Sebastian,et al. Pore network extraction for fractured porous media[J]. Advances in Water Resources,2017,107:280-289.

[10] JU Yang,REN Zhangyu,ZHENG Jiangtao,et al. Quantitative visualization methods for continuous evolution of three-dimensional discontinuous structures and stress field in subsurface rock mass induced by excavation and construction-An overview[J]. Engineering Geology,2020,265:105443.

[11] MING C,HORII H. A constitutive model of highly jointed rock masses[J]. Mechanics of Materials,1992,13(3):217-246.

[12] LI Shihai,ZHOU Dong. Progressive failure constitutive model of fracture plane in geomaterial based on strain strength distribution[J]. International Journal of Solids and Structures,2013,50(3):570-577.

[13] 謝和平. 深部巖體力學(xué)與開采理論研究進(jìn)展[J]. 煤炭學(xué)報(bào),2019,44(5):1283-1305.

XIE Heping. Research review of the state key research development program of China:Deep rock mechanics and mining theory[J]. Journal of China Coal Society,2019,44(5):1283-1305.

[14] JU Yang,XIE Heping,ZHENG Zemin,et al. Visualization of the complex structure and stress field inside rock by means of 3D printing technology[J]. Chinese Science Bulletin,2014,59(36):5354-5365.

[15] KETCHAMRICHARD A,CARLSON William D. Acquisition,optimization and interpretation of X-ray computed tomographic imagery:Applications to the geosciences[J]. Computers and Geosciences,2001,27(4):381-400.

[16] BLUNTMARTIN J,BIJELJIC Branko,DONG Hu,et al. Pore-scale imaging and modelling[J]. Advances in Water Resources,2013,51:197-216.

[17] WILDENSCHILD Dorthe,SHEPPARD Adrian P. X-ray imaging and analysis techniques for quantifying pore-scale structure and processes in subsurface porous medium systems[J]. Advances in Water Resources,2013,51(1):217S246.

[18] CNUDDE V,BOONE M N. High-resolution X-ray computed tomography in geosciences:A review of the current technology and applications[J]. Earth-Science Reviews,2013,123(4):1-17.

[19] WERTHCHARLES J,ZHANG Changyong,BRUSSEAU Mark L,et al. A review of non-invasive imaging methods and applications in contaminant hydrogeology research[J]. Journal of Contaminant Hydrology,2010,113(1-4):1-24.

[20] XIONG Q R,BAYCHEV T G,JIVKOV A P. Review of pore network modelling of porous media:Experimental characteris-ations,network constructions and applications to reactive transport[J]. Journal of Contaminant Hydrology,2016,192:101-117.

[21] NIASHAHRAM Farhadi,DASANI Devang,TSOTSIS Theodore T,et al. An integrated approach for the characterization of shales and other unconventional resource materials[J]. Industrial and Engineering Chemistry Research,2016,55(12):3718-3728.

[22] WANG Moran,PAN Ning. Modeling and prediction of the effective thermal conductivity of random open-cell porous foams[J]. International Journal of Heat and Mass Transfer,2008,51(5-6):1325-1331.

[23] ZHENG Jiangtao,WANG Ziyan,GONG Wenbo,et al. Characterization of nanopore morphology of shale and its effects on gas permeability[J]. Journal of Natural Gas Science and Engineering,2017,47:83-90.

[24] JU Yang,ZHANG Qinggang,YANG Yongming,et al. An experimental investigation on the mechanism of fluid flow through single rough fracture of rock[J]. Science China-Technological Sciences,2013,56(8):2070-2080.

[25] JU Yang,ZHANG Qingang,ZHENG Jiangtao,et al. Fractal model and lattice Boltzmann method for characterization of non-Darcy flow in rough fractures[J]. Scientific Reports,2017,7:1-9.

[26] OKABE H,BLUNT M J. Pore space reconstruction using multiple-point statistics[J]. Journal of Petroleum Science and Engineering,2005,46(1):121-137.

[27] JU Yang,ZHENG Jiangtao,EPSTEIN Marcelo,et al. 3D numerical reconstruction of well-connected porous structure of rock using fractal algorithms[J]. Computer Methods in Applied Mechanics and Engineering,2014,279(9):212-226.

[28] JIAO Y,STILLINGER F,TORQUATO S. A superior descriptor of random textures and its predictive capacity[J]. Proceedings of the National Academy of Sciences,2009,106(42):17634-17639.

[29] JU Yang,HUANG Yaohui,GONG Wenbo,et al. 3-D reconstruction method for complex pore structures of rocks using a small number of 2-D X-ray computed tomography images[J]. IEEE Transactions on Geoscience and Remote Sensing,2018,57(4):1873-1182.

[30] GERKEKIRILL M,KARSANINA Marina V. Improving stochastic reconstructions by weighting correlation functions in an objective function[J]. Europhysics Letters,2015,111(5):1-5.

[31] KARSANINA Marina,GERKE Kirill. Hierarchical optimization:fast

and robust multiscale stochastic reconstructions with rescaled correlation functions[J]. Physical Review Letters,2018,121(26):1-6.

[32] JU Yang,HUANG Yaohui,ZHENG Jiangtao,et al. Multi-thread parallel algorithm for reconstructing 3D large-scale porous structures[J]. Computers and Geosciences,2017,101:10-20.

[33] ?REN P E,BAKKE S. Process based reconstruction of sandstones and prediction of transport properties[J]. Transport in Porous Media,2002,46(2):311-343.

[34] MOSSER Lukas,DUBRULE Olivier,BLUNT MARTIN J. Reconstruction of three-dimensional porous media using generative adversarial neural networks[J]. Physical Review E,2017,96(4):1-17.

[35] PERRASMATTHEW A,VOGLER Daniel. Compressive and tensile behavior of 3D-printed and natural sandstones[J]. Transport in Porous Media,2018,129(2):559-581.

[36] ZHOU Tao,ZHU Jianbo. Identification of a suitable 3D printing mate-rial for mimicking brittle and hard rocks and its brittleness enhancements[J]. Rock Mechanics and Rock Engineering,2018,51(3):765-777.

[37] ZHOU Tao,ZHU Jianbo. An experimental investigation of tensile fracturing behavior of natural and artificial rocks in static and dynamic brazilian disc tests[J]. Procedia Engineering,2017,191:992-998.

[38] GELLE Leanor M,WALLEY Stephen M,BRAITHWAITE christo-

pher H. Review of the validity of the use of artificial specimens for characterizing the mechanical properties of rocks[J]. Rock Mechanics and Rock Engineering,2019,52(9):2949-2961.

[39] WANG Li,JU Yang,XIE Heping,et al. The mechanical and photoelastic properties of 3D printable stress-visualized materials[J]. Scientific Reports,2017,7(1):10918.

[40] JU Yang,GUO Wang,REN Zhangyu,et al. Experimental study on

mechanical and optical properties of printable photopolymer used for visualising hidden structures and stresses in rocks[J]. Optical Materials,2021,111:110691.

[41] JU Yang,REN Zhangyu,HU Xuan,et al. Critical phase-transition temperature for freezing stress in thermo-sensitive photopolymers used for visualizing stress fields in solids[J]. Optics and Lasers in Engineering,2021,139:106458.

[42] JU Yang,REN Zhangyu,LI Xiaolan,et al. Quantification of hidden whole-field stress inside porous geomaterials via three-dimensional printing and photoelastic testing methods[J]. Journal of Geophysical Research Solid Earth,2019,124(6):5408-5426.

[43] JING L. A review of techniques,advances and outstanding issues in numerical modelling for rock mechanics and rock engineering[J]. International Journal of Rock Mechanics and Mining Sciences,2003,40(3):283-353.

[44] KATO Naoyuki. Interaction of slip on asperities:Numerical simulation of seismic cycles on a two-dimensional planar fault with nonuniform frictional property[J]. Journal of Geophysical Research:Solid Earth,2004,109(B12):1-17.

[45] ISHUTOV Sergey,HASIUK Franciszek J,HARDING Chris,et al. 3D printing sandstone porosity models[J]. Interpretation,2015,3(3):SX49-SX61.

[46] ISHUTOV Sergey,HASIUK Franciszek J,FULLMER Shawn M,et al. Resurrection of a reservoir sandstone from tomographic data using three-dimensional printing[J]. AAPG Bulletin,2017,101(9):1425-1443.

[47] ISHUTOV S,HASIUK F J,JOBE D,et al. Using resin-based 3D printing to build geometrically accurate proxies of porous sedimentary rocks[J]. Ground Water,2018,56(3):482-490.

[48] JIANG Quan,FENG Xiating,GONG Yanhua,et al. Reverse modelling of natural rock joints using 3D scanning and 3D printing[J]. Computers and Geotechnics,2016,73:210-220.

[49] JU Yang,WANG Jianqiang,PENG Ruidong,et al. Industrial CT scanning test system[P]. US Patent：15539835，2018-02-13.

[50] JU Yang,WANG Jianqiang,PENG Ruidong,et al. Industrial CT scanning test system and fluid pressure loading apparatus[P]. US Patent：15539835,2018-01-30.

[51] 鞠楊,王建強(qiáng),彭瑞東,等. 一種工業(yè)CT掃描試驗(yàn)方法及試驗(yàn)系統(tǒng)[P]. 中國(guó)專利：CN106706676B,2019-07-26.

[52] JU Yang,WANG Li,XIE Heping,et al. Visualization and transparentization of the structure and stress field of aggregated geomaterials through 3D printing and photoelastic techniques[J]. Rock Mechanics and Rock Engineering,2017,50(6):1383-1407.

[53] JU Yang,WANG Li,XIE Heping,et al. Visualization of the three-

dimensional structure and stress field of aggregated concrete materials through 3D printing and frozen-stress techniques[J]. Construction and Building Materials,2017,143:121-137.

[54] JU Yang,REN Zhangyu,MAO Lingtao,et al. Quantitative visualisation of the continuous whole-field stress evolution in complex pore structures using photoelastic testing and 3D printing methods[J]. Optics Express,2018,26(5):6182-6201.

[55] HILLER Jonathan,LIPSON Hod. Design and analysis of digital materials for physical 3D voxel printing[J]. Rapid Prototyping Journal,2009,15(2):137-149.

[56] SOSSOU Germain,DEMOLY Frédéric,BELKEBIR Hadrien,et al. Design for 4D printing:A voxel-based modeling and simulation of smart materials[J]. Materials & Design,2019,175:107798.

[57] BELL Darren,SIEGMUND Thomas. 3D-printed polymers exhibit a strength size effect[J]. Additive Manufacturing,2018,21:658-665.

[58] AKIYAMA Hiroyuki,UCHIDA Makoto,KANEKO Yoshihisa. Evaluation of effect of sample size and layer direction on mechanical property of specimen manufactured by FDM-type 3D printer[J]. Key Engineering Materials,2019,794:324-332.

[59] THEOCARISPERICLES S,GDOUTOS Emmanuel E. Matrix theory of photoelasticity[M]. Berlin:Springer,2013.

[60] FROCHTMAX Mark. Photoelasticity:The selected scientific papers of MM Frocht[M]. Pergamon：1969.

[61] JU Yang,REN Zhangyu,WANG Li,et al. Photoelastic method to quantitatively visualise the evolution of whole-field stress in 3D printed models subject to continuous loading processes[J]. Optics and Lasers in Engineering,2018,100:248-258.

[62] JU Yang,WANG Yating,REN Zhangyu,et al. Optical method to quantify the evolution of whole-field stress in fractured coal subjected to uniaxial compressive loads[J]. Optics and Lasers in Engineering,2020,128:106013.

[63] YIN Peijie,ZHAO Can,MA Jianjun,et al. Experimental study of non-linear fluid flow through rough fracture based on fractal theory and 3D printing technique[J]. International Journal of Rock Mechanics and Mining Sciences,2020,129:104293.

[64] JU Yang,WAN Changbing,REN Zhangyu,et al. Quantification of

continuous evolution of full-field stress associated with shear deformation of faults using three-dimensional printing and phase-shifting methods[J]. International Journal of Rock Mechanics Mining Sciences,2020,126:104187.

[65] JU Yang,WAN Changbing,FU Guoming,et al. Optical method for identification and quantification of full-field stress distributions and evolution in assembled lining structures based on additively printed models and phase-shifting methods[J]. Tunnelling and Underground Space Technology,2022,119:104228.

[66] MISSERONI D,DAL CORSO F,SHAHZAD S,et al. Stress concentration near stiff inclusions:Validation of rigid inclusion model and boundary layers by means of photoelasticity[J]. Engineering Fracture Mechanics,2014,121-122:87-97.

[67] DANIELSKAREN E,HAYMAN Nicholas W. Force chains in seismogenic faults visualized with photoelastic granular shear experiments[J]. Journal of Geophysical Research:Solid Earth,2008,113(B11):B11411.

[68] BIGONI Davide,NOSELLI Giovanni. Localized stress percolation

through dry masonry walls. Part I - Experiments[J]. European Journal of Mechanics-A/Solids,2010,29(3):291-298.

[69] BIGONI Davide,NOSELLI Giovanni. Localized stress percolation thr-

ough dry masonry walls. Part II-Modelling[J]. European Journal of Mechanics-A/Solids,2010,29(3):299-307.

[70] GHISLAIN de Joussineau,PETIT JeanPierre,GAUTHIERBERTRAND D M. Photoelastic and numerical investigation of stress distributions around fault models under biaxial compressive loading conditions[J]. Tectonophysics,2003,363(1-2):19-43.

[71] SOLIVA Roger,MAERTEN Frantz,PETIT JeanPierre,et al. Field evidences for the role of static friction on fracture orientation in extensional relays along strike-slip faults:Comparison with photoelasticity and 3-D numerical modeling[J]. Journal of Structural Geology,2010,32(11):1721-1731.

[72] LUNINA O V,GLADKOV A S,AFONKIN A M,et al. Deformation style in the damage zone of the Mondy fault:GPR evidence(Tunka basin,southern East Siberia)[J]. Russian Geology and Geophysics,2016,57(9):1269-1282.

[73] CERNOSEK Jan. Three-dimensional photoelasticity by stress freezing

[J]. Experimental Mechanics,1980,20(12):417-426.

[74] KUSKE Albrecht,ROBERTSON George Stewart. Photoelastic stress analysis[M]. New Jersey :Wiley-Interscience,1974.

[75] FROCHTMAX Mark,GUERNSEY JR R. A special investigation to develop a general method for three-dimensional photoelastic stress analysis[J]. National Advisory Committee for Aeronautics,1952.

[76] GUERNSEY R. A special investigation to develop a general method for three-dimensional photoelastic stress analysis[M]. Oxford:Pergamon,1969,233-269.

[77] SWAIN Digendranath,PHILIP Jeby,PILLAI Annamala,et al. A revisit to the frozen stress phenomena in photoelasticity[J]. Experimental Mechanics,2016,56:903-917.

[78] HEYWOOD R B. Photoelasticity for Designers [M]. Oxford :Pergamon,1969.

[79] 謝和平,高峰,鞠楊,等. 深地煤炭資源流態(tài)化開采理論與技術(shù)構(gòu)想[J]. 煤炭學(xué)報(bào),2017,42(3):547-556.

XIE Heping,GAO Feng,JU Yang,et al. Theoretical and technological conception of the fluidization mining for deep coal resources[J]. Journal of China Coal Society,2017,42(3):547-556.

[80] LEE S,RAVICHANDRAN G. An investigation of cracking in brittle solids under dynamic compression using photoelasticity[J]. Optics and Lasers in Engineering,2003,40(4):341-352.

[81] JU Yang,XIE Heping,ZHAO Xi,et al. Visualization method for stress-field evolution during rapid crack propagation using 3D printing and photoelastic testing techniques[J]. Scientific Reports,2018,8(1):1-10.

[82] XIA Kaiwen,ROSAKISARES J,HIROO Kanamori. Laboratory ear-

thquakes:The sub-Rayleigh-to-supershear rupture transition[J]. Science,2004,303(5665):1859-1861.

[83] XIA Kaiwen,ROSAKISARES J. Laboratory earthquakes along faults

with a low velocity zone:Directionality and pulse-like ruptures[J]. Extreme Mechanics Letters,2021,46:101321.

[84] GABUCHIAN Vahe,ROSAKIS Ares,BHAT Harsha S,et al. Experimental evidence that thrust earthquake ruptures might open faults[J]. Nature,2017,545:336-340.

[85] GABUCHIAN Vahe,ROSAKIS Ares J,LAPUSTA Nadia,et al. Experimental investigation of strong ground motion due to thrust-fault earthquakes[J]. Journal of Geophysical Research Solid Earth,2014,119(2):133-143.

[86] ASUNDI Anand,TONG Liu,BOAY Chai Gin. Dynamic phase-shifting photoelasticity[J]. Applied Optics,2001,40(22):3654-3658.

[87] ASUNDI Anand,SAJAN M R,LIU Tong. Dynamic photoelasticity using TDI imaging[J]. Optics & Lasers in Engineering,2002,38(1-2):3-16.

[88] AJOVALASIT A,PETRUCCI G,SCAFIDI M. Review of RGB photoelasticity[J]. Optics & Lasers in Engineering,2015,68:58-73.

[90] AJOVALASIT A,PETRUCCI G. Developments in RGB Photoelasticity[J]. Applied Mechanics & Materials,2005,3-4:205-210.

[91] RAMAKRISHNAN Vivek,RAMESH K. Scanning schemes in white light Photoelasticity-Part I:Critical assessment of existing schemes[J]. Optics and Lasers in Engineering,2017,92:129-140.

[92] RAMAKRISHNAN Vivek,RAMESH K. Scanning schemes in white light photoelasticity-Part II:Novel fringe resolution guided scanning scheme[J]. Optics and Lasers in Engineering,2017,92:141-149.

[93] NG T W. Photoelastic stress analysis using an object step-loading method[J]. Experimental Mechanics,1997,37(2):137-141.

[94] ZUCCARELLO B,TRIPOLI G. Photoelastic stress pattern analysis using Fourier transform with carrier fringes:Influence of quarter-wave plate error[J]. Optics & Lasers in Engineering,2002,37(4):401-416.

[95] RAMESH K,KASIMAYAN T,NEETHI Simon B. Digital photoelasticity - A comprehensive review[J]. The Journal of Strain Analysis for Engineering Design,2011,46:245-266.

[96] MANOGG Peter. Investigation of the rupture of a plexiglas plate by means of an optical method involving high-speed filming of the shadows originating around holes drilled in the plate[J]. International Journal of Fracture Mechanics,1966,2(4):604-613.

[97] QIU Peng,YUE Zhongwen,YANG Renshu. Mode I stress intensity factors measurements in PMMA by caustics method:A comparison between low and high loading rate conditions[J]. Polymer Testing,2019,76:273-285.

[98] GAO Guiyun,LI Zheng,NEGAHBAN Mehrdad. Dynamic fracture analysis of polycarbonate by the optical method of caustics[J]. Procedia Materials Science,2014,3:165-176.

[99] PAPADOPOULOS G A. New formula of experimental stress intensity factor evaluation by caustics[J]. International Journal of Fracture,2011,171(1):79.

[100] REN Zhangyu,XIE Huimin,JU Yang. Quantification of photoela-stic fringe orders using polarized light camera and continuous loading[J]. Optics and Lasers in Engineering,2020,134:106263.

[101] JU Yang,LIU Peng,REN Zhangyu,et al. Characterization of stress field evolution during 3D internal fracture propagation using additively printed models and frozen stress techniques[J]. Theoretical and Applied Fracture Mechanics,2021,111:102870.

[102] LIU Peng,JU Yang,FU Guoming,et al. Visualization of full-field stress evolution during 3D penetrated crack propagation through 3D printing and frozen stress techniques[J]. Engineering Fracture Mechanics,2020,236:107222.

[103] BAYBRIAN K,Smith TAIT S,Fyhrie DAVID P,et al. Digital volume correlation:Three-dimensional strain mapping using X-ray tomography[J]. Experimental Mechanics，1999,39(3):217-226.

[104] SMITHTAIT S,BAY BRIAN K,RASHID Mark M. Digital volume correlation including rotational degrees of freedom during minimization[J]. Experimental Mechanics,2002,42(3):272-278.

[105] LIU Li,MORGAN Elise. Accuracy and precision of digital volume correlation in quantifying displacements and strains in trabecular bone[J]. Journal of Biomechanics,2007,40(15):3516-3520.

[106] OBATA Yoshihiro,BALE HRISHIKESH A,BARNARD HAROLD S,et al. Quantitative and qualitative bone imaging:A review of synchrotron radiation microtomography analysis in bone research[J]. 2020,110:103887.

[107] WANG Long,ZHANG Wei,LI Haibo,et al. 3D in-situ characterizations of damage evolution in c/sic composite under monotonic tensile loading by using X-Ray computed tomography[J]. Applied Composite Materials,2020,27(1):1-12.

[108] RETHORE J,LIMODIN N,BUFFIERE J Y,et al. Digital volume correlation analyses of synchrotron tomographic images[J]. Journal of Strain Analysis for Engineering Design,2011,46(7):683-695.

[109] BULJAC A,JAILIN C,MENDOZA A,et al. Digital volume correlation:Review of progress and challenges[J]. Experimental Mechanics,2018,58(5):1-48.

[110] XU Feng. Quantitative characterization of deformation and damage process by digital volume correlation:A review[J]. Theoretical and Applied Mechanics Letters,2018,8(2):83-96.

[111] XU Feng. Quantitative characterization of deformation and damage process by digital volume correlation:A review[J]. Theoretical and Applied Mechanics Letters,2018,8(2):83-96.

[112] BULJAC A,JAILIN C,MENDOZA A,et al. Digital volume correlation:Review of progress and challenges[J]. Experimental Mechanics,2018,58(5):661-708.

[113] PAN Bing,WANG Bo. Some recent advances in digital volume correlation[J]. Optics and Lasers in Engineering,2020,135:106189.

[114] MAO Lingtao,LIU Haizhou,ZHU Ying,et al. 3D strain mapping of opaque materials using an improved digital volumetric speckle photography technique with X-ray microtomography[J]. Applied Sciences,2019,9(7):1418.

[115] MAO Lingtao,ZHU Ying,WANG Yuwei,et al. An improved digital volumetric speckle photography technique with X-ray microtomography and its applications to investigating strain localization in red sandstone[J]. Applied Sciences,2020,53(3):1457-1466.

[116] MAO Lingtao,HAO Nai,AN Liqian,et al. 3D mapping of carbon dioxide-induced strain in coal using digital volumetric speckle photography technique and X-ray computer tomography[J]. International Journal of Coal Geology,2015,147-148:115-125.

[117] MAO Lingtao,CHIANG Fupen. 3D strain mapping in rocks using digital volumetric speckle photography technique[J]. Acta Mechanica,2016,227(11):3069-3085.

[118] MAO Lingtao,YUAN Zexun,YANG Min,et al. 3D strain evolution in concrete using in situ X-ray computed tomography testing and digital volumetric speckle photography[J]. Measurement,2019,133:456-467.

[119] LENOIR N,BORNERT M,DESRUES J,et al. Volumetric digital image correlation applied to X-ray microtomography images from triaxial compression tests on argillaceous rock[J]. Strain,2007,43(3):193-205.

[120] GANDOMKAR Asghar,RAHIMPOUR Mohammad Reza. Investigation of low-salinity waterflooding in secondary and tertiary enhanced oil recovery in limestone reservoirs[J]. Energy and Fuels,2015,29:7781-7792.

[121] HerringANNA L,ANDERSSON Linnéa,NEWELL D L,et al. Pore-

scale observations of supercritical COdrainage in Bentheimer sandstone by synchrotron X-ray imaging[J]. International Journal of Greenhouse Gas Control,2014,25:93-101.

[122] SKAUGE Arne,ZAMANI Nematollah,GAUSDAL Jacobsen J?rgen,et al. Polymer flow in porous media:Relevance to enhanced oil recovery[J]. Colloids and Interfaces,2018,2:27.

[123] JU Yang,GONG Wenbo,ZHENG Jiangtao. Characterization of immiscible phase displacement in heterogeneous pore structures:Parallel multicomponent lattice Boltzmann simulation and experimental validation using three-dimensional printing technology[J]. International Journal of Multiphase Flow,2019,114:50-65.

[124] JU Yang,GONG Wenbo,CHANG Wei,et al. Effects of pore characteristics on water-oil two-phase displacement in non-homogeneous pore structures:A pore-scale lattice Boltzmann model considering various fluid density ratios[J]. International Journal of Engineering Science,2020,154:103343.

[125] LIU Jinhui,JU Yang,ZHANG Yingqi,et al. Preferential paths of air-water two-phase flow in porous structures with special consideration of channel thickness effects[J]. Scientific Reports,2019,9(1):16204.

[126] ZHENG Jiangtao,WANG Zenglin,JU Yang,et al. Visualization of water channeling and displacement diversion by polymer gel treatment in 3D printed heterogeneous porous media[J]. Journal of Petroleum Science and Engineering,2021,198:108238.

[127] CNUDDE V,BOONE M N. High-resolution X-ray computed tomography in geosciences:A review of the current technology and applications[J]. Earth-Science Reviews,2013,123:1-17.

[128] CONNOLLY Paul R J,VOGT Sarah J,IGLAUER Stefan,et al. Capillary trapping quantification in sandstones using NMR relaxometry[J]. Water Resources Research,2017,53:7917-7932.

[129] ISHAKOGLU Ali,BAYTAS A Filiz.The influence of contact angle on

capillary pressure-saturation relations in a porous medium including various liquids[J]. International Journal of Engineering Science,2005,43:744-755.

[130] ISHIBASHI Takuya,MCGUIRE Thomas P,WATANABE Noriaki,et al. Permeability evolution in carbonate fractures:Competing roles of confining stress and fluid pH[J]. Water Resources Research,2013,49:2828-2842.

[131] JU Yang,XI Chaodong,ZHENG Jiangtao,et al. Study on three-dimensional immiscible water-oil two-phase displacement and trapping in deformed pore structures subjected to varying geostress via in situ computed tomography scanning and additively printed models[J]. International Journal of Engineering Science,2022,171:103615.

[132] ALHAMMADIAMER M,GAO Ying,AKAI Takashi,et al. Pore-scale X-ray imaging with measurement of relative perm-eability,capillary pressure and oil recovery in a mixed-wet micro-porous carbonate reservoir rock[J]. Fuel,2020,268:117018.

[133] GAO Ying,QASEMINEJAD Raeini Ali,BLUNT Martin J,et al. Pore occupancy,relative permeability and flow intermittency measurements using X-ray micro-tomography in a complex carbonate[J]. Advances in Water Resources,2019,129:56-69.

[134] HERRING Anna L,ANDERSSON Linnéa,SCHLüTER Steffen,et al. Efficiently engineering pore-scale processes:The role of force dominance and topology during nonwetting phase trapping in porous media[J]. Advances in Water Resources,2015,79:91-102.

[135] JU Yang,XI Chaodong,WANG Shaojie,et al. 3-D fracture evolution and water migration in fractured coal under variable stresses induced by fluidized mining:In situ triaxial loading and CT imaging analysis[J]. Energy Reports,2021,7:3060-3073.

[136] PHILIPPOFF Wladimir. Flow-Birefringence and Stress[J].Nature,

1956,178:811-812.

[137] PRADOS J W,PEEBLES F N. Two-dimensional laminar-flow anal-

ysis,utilizing a doubly refracting liquid[J]. 1959,5(2):225-234.

[138] MULLER R,VERGNES B. Validity of the stress optical law and application of birefringence to polymer complex flows[M]. Amsterdam:Elsevier,1996,257-284.

[139] SENGUPTA Anupam,PIEPER Christoph,ENDERLEIN J?rg,et al. Flow of a nematogen past a cylindrical micro-pillar[J]. Soft Matter,2013,9(6):1937-1946.

[140] AUHL D,HOYLE D M,HASSELL D,et al. Cross-slot extensional rheometry and the steady-state extensional response of long chain branched polymer melts[J]. Journal of Rheology,2011,55(4):875-900.