劉善琪，李永兵，田會全，劉旭耀，朱伯靖，石耀霖

中國科學(xué)院計算地球動力學(xué)重點實驗室，中國科學(xué)院大學(xué)，北京 100049

1 引 言

巖石熱導(dǎo)率是巖石主要的物性參數(shù)之一，對地質(zhì)基礎(chǔ)研究（巖石圈熱結(jié)構(gòu)、地球內(nèi)部熱傳導(dǎo)過程、地球深部熱狀態(tài)、深部巖石的變質(zhì)作用和熱演化歷史等）和應(yīng)用研究（礦山采掘、石油開發(fā)和地熱能利用）均具有重要的意義.在巖石圈內(nèi)，熱傳導(dǎo)是巖石間熱傳遞的最主要方式.巖石圈內(nèi)的各種地球物理過程（如巖石的流變狀態(tài)、熱結(jié)構(gòu)、大地熱流等）都與巖石的熱導(dǎo)率緊密相關(guān).

目前取得巖石熱導(dǎo)率的方法主要有兩種：實驗室測量和原位測量.實驗室測量一般是在實驗室對經(jīng)過風干（少部分為烘干）的巖芯樣品進行測定.這種方法對于低孔隙率的巖石，不會引起大的誤差；對高孔隙率的碎屑沉積巖類，則會產(chǎn)生顯著的影響.大多數(shù)巖石都具有孔隙或裂隙，即使致密的結(jié)晶巖類，也有一定程度的孔隙或裂隙.泥巖、頁巖、砂巖等巖石的孔隙率可高達20%～30%，對于這些含孔隙的巖石，干試樣的熱導(dǎo)率與含水試樣的熱導(dǎo)率存在不同程度的差異［1］.忽略巖石孔隙中的水對熱導(dǎo)率的影響是實驗結(jié)果與實際情況不符的主要原因.因此，最理想的熱導(dǎo)率測量方法應(yīng)該在原位進行，現(xiàn)已發(fā)展了一系列的原位測量方法，但是原位測量不僅方法復(fù)雜，而且一些巖石無法進行原位測量.同時，上述兩種方法都必須耗費大量時間.

實際上，孔隙巖石是一種多孔介質(zhì)［2］，并且往往含有水.因此，從巖石的性質(zhì)出發(fā)，建立孔隙巖石的數(shù)學(xué)模型或經(jīng)驗公式來預(yù)測巖石的熱導(dǎo)率一直是令人感興趣的問題［3］.前人在這方面已開展大量研究，主要是基于孔隙巖石的特性建立簡化多孔介質(zhì)導(dǎo)熱模型，即：（1）統(tǒng)計模型［4］；（2）局部結(jié)構(gòu)模型［5］；（3）半經(jīng)驗?zāi)Ｐ停?］；（4）容積平均模型［7］；（5）分形模型［8－9］.其中，前四種經(jīng)典模型無法準確描述多孔介質(zhì)內(nèi)孔隙分布情況及結(jié)構(gòu)特征，只能近似地在大尺度范圍內(nèi)描述多孔介質(zhì)中的熱量傳導(dǎo)過程，無法揭示真實的熱量傳遞規(guī)律和溫度場的分布情況，模型的計算結(jié)果與實際測量值有較大偏差［10］；而分形模型對分析傳遞過程本身尚未取得有意義的進展，也沒能定量化建立導(dǎo)熱和結(jié)構(gòu)之間的關(guān)系.

針對含濕多孔介質(zhì)有效熱導(dǎo)率的數(shù)學(xué)模型研究主要有：Bogaty［11］和 Hollies［12］對含有不同纖維、孔隙率及飽和度的介質(zhì)進行了實驗，發(fā)現(xiàn)織物的有效熱導(dǎo)率主要受厚度、孔隙率和飽和度的影響，纖維的排列方式對它的影響不大；Singh等［13］在Pande等［14－15］模型的基礎(chǔ)上，考慮了有效連續(xù)介質(zhì)近似值及其所有可能的相互作用，計算了含濕土壤在不同溫度下的有效熱導(dǎo)率，其結(jié)果與實驗值吻合較好；Zhang等［16－18］提出了用多組分體系的隨機混合模型預(yù)測含濕多孔介質(zhì)有效熱導(dǎo)率的方法，該模型對于孔隙率小于0.6的含濕多孔介質(zhì)準確度很高，而對于孔隙率大于0.6的含濕多孔介質(zhì)有一定程度的偏離；Wang等［19］提出了基于三維微觀方法來預(yù)測有效熱導(dǎo)率的隨機多孔介質(zhì)模型；Tiak等［20］根據(jù)傅里葉定律和能量守恒定律提出含濕多孔介質(zhì)微元體熱導(dǎo)率模型；Aurangzeb等［21］采用美國材料實驗協(xié)會的ASTM標準測量密度和孔隙率，在保證實驗精確度的前提下提出了指數(shù)衰減模型.這些模型雖然考慮了多孔介質(zhì)中液相的存在，便于計算確定熱導(dǎo)率.但是，上述方法對孔隙分布情況的處理與自然界的真實情況還有一些差距，得出的經(jīng)驗公式僅在一定程度上可用.

前人的研究表明，應(yīng)用數(shù)值分析方法獲得準確的含濕孔隙巖石有效熱導(dǎo)率的前提是建立符合實際情況的巖石數(shù)值模型.本文在景惠敏對雙礦物熱導(dǎo)率數(shù)值模擬的基礎(chǔ)上［22］，將其發(fā)展至多組分體系，采用隨機網(wǎng)格劃分和材料指定的方法，建立接近實際情況的巖石模型，用有限元方法進行含濕孔隙巖石有效熱導(dǎo)率的三維數(shù)值模擬，并將數(shù)值模擬所得的大量數(shù)據(jù)結(jié)果進行統(tǒng)計分析，最終獲得巖石有效熱導(dǎo)率.含濕孔隙巖石采用圓柱體模型，上下表面施加不同的溫度，側(cè)面絕熱.通過有限元法計算出總熱流q，結(jié)合上下表面的溫度梯度，計算巖石的有效熱導(dǎo)率.

2 方法原理

2.1 基本原理

含濕孔隙巖石可以看作由礦物、空氣和水構(gòu)成的三組分復(fù)合介質(zhì)，分別用下標m、a、w標記.本文用孔隙率f和飽和度s來描述其組成.孔隙率f為水和空氣所占的巖樣體積的體積分數(shù)，飽和度s為孔隙中水所占的體積分數(shù).這樣，礦物、水和空氣的體積分數(shù)分別為：1－f、f×s、f×（1－s）.

在無熱源條件下，對于三維穩(wěn)態(tài)熱傳導(dǎo)問題［23］，其偏微分控制方程為

邊界條件是

kx、ky、kz是材料分別沿物體三個主方向（x、y、z方向）的熱導(dǎo)率，單位為 W·m－1·K－1，φ＝是在Γ1邊界上的給定溫度.

對于微分方程（1）和邊界條件（2），可推導(dǎo)出變分原理中的泛函表達式

利用δΠ（φ）＝0可以建立有限元的求解方程.

將求解域Ω離散為有限個單元體，單元內(nèi)各點的溫度φ近似地由單元的節(jié)點溫度φi插值得到

其中ne是每個單元的節(jié)點數(shù)；Ni（x，y，z）是C0型插值函數(shù)，具有下述性質(zhì)：

將（4）式代入有限元離散后的泛函，從δΠ（φ）＝0，可以得到穩(wěn)態(tài)熱傳導(dǎo)問題的有限元求解方程：

式中K是熱傳導(dǎo)矩陣，它是對稱矩陣，在引入給定溫度條件以后，K是對稱正定的；φ＝ ［φ1φ2…φn］T是節(jié)點溫度列陣.矩陣K元素表示如下：

（8）式中右側(cè)是各單元對熱傳導(dǎo)矩陣的貢獻，可將其改寫成單元集成的形式：

以上就是三維穩(wěn)態(tài)熱傳導(dǎo)問題的有限元方程.三維穩(wěn)態(tài)熱傳導(dǎo)的傅里葉定律為

其中，q表示熱流密度.有限元模型的等效熱導(dǎo)率為

式中，Keff為模型的等效熱導(dǎo)率，qi為節(jié)點的熱流密度，ΔT為模型上下表面的溫度差，h為模型垂直方向的高度，N為節(jié)點總數(shù).

2.2 模型及網(wǎng)格剖分方法

本數(shù)值分析模型基于四面體網(wǎng)格，并假設(shè)由若干網(wǎng)格組成一個礦物顆粒，在給定孔隙率及飽和度的前提下，固液氣三相的空間分布完全隨機確定.巖石的有效熱導(dǎo)率與孔隙的分布有密切關(guān)系［24］.圖1為巖石中只有一種流體相時，巖石中孔隙分布的極端情況［25］.巖石可能包含比例不同的多種礦物，但在初步模擬中我們暫時先假定它們的熱導(dǎo)率相差不大，近似視為一種礦物，稱為礦物的熱導(dǎo)率.如圖1a所示，當熱流與孔隙平行時，巖石的有效熱導(dǎo)率最大，可表示為

其中Keff為巖石的有效熱導(dǎo)率，Kf為流體相的熱導(dǎo)率，Km為巖石骨架的熱導(dǎo)率，f為孔隙率.當熱流與孔隙垂直時（圖1b），巖石的有效熱導(dǎo)率最小，可表示為

圖1 巖石中孔隙分布的極端情況示意［25］Fig.1 Sketch of the extreme situation of pore distribution［25］

考慮到固體礦物成分、孔隙率及飽和度完全相同的巖石，由于內(nèi)部孔隙分布和材料分布的不同，其導(dǎo)熱性質(zhì)會有很大的不同，因此為獲得最符合客觀實際的含濕孔隙巖石的熱導(dǎo)率，模型采用隨機剖分單元、隨機指定材料屬性的方法來模擬巖石內(nèi)部各種孔隙分布情況下的巖石導(dǎo)熱性.基于此，先生成四面體網(wǎng)格，然后由若干網(wǎng)格組成一個顆粒，并按一定混合比例，隨機對各個顆粒指定材料性質(zhì)，模型內(nèi)部的孔隙分布與每次試驗所產(chǎn)生的一個偽隨機數(shù)相關(guān)，在固定孔隙率及飽和度的條件下，每次試驗都可以產(chǎn)生不同的孔隙分布，而且孔隙的形狀、大小、方向不規(guī)則，更接近于真實情況，使數(shù)值模擬結(jié)果更加客觀真實.

本文的網(wǎng)格劃分方法可任意控制組成一個顆粒的網(wǎng)格數(shù)，并且可以按照研究需要為不同的顆粒賦予相應(yīng)的材料屬性.首先對模型進行四面體網(wǎng)格劃分.假設(shè)生成的網(wǎng)格總數(shù)為maxelem，需要生成的顆粒數(shù)為M，則組成一個顆粒的網(wǎng)格數(shù)為n＝int（maxelem／M）.

其中，x為由計算機產(chǎn)生（0～1）之間的隨機數(shù)，i為單元號，1、2、3為材料編號，依次代表礦物、水、空氣.

從第一個節(jié)點開始按節(jié)點號進行單元拓撲，搜索以此節(jié)點為頂點的單元，按單元號從小到大的順序選定n個單元，并據(jù)（15）式賦材料屬性mate（i）.如果在某次循環(huán)中，進行單元拓撲時，搜索到的單元已被賦屬性，則跳過繼續(xù)搜索其它單元；如果以此節(jié)點為頂點的單元都已被賦屬性，則進行下一個節(jié)點的單元拓撲，如此循環(huán)M次.如果循環(huán)M次后，還有單元未賦屬性，則再循環(huán)一次.

圖2 孔隙率很小時部分模型中水、空氣的分布Fig.2 Distribution of water and air in some models when porosity is very small

為驗證文中采用方法的正確性及合理性，先建立一個網(wǎng)格數(shù)量較少的模型（3960個網(wǎng)格），6個網(wǎng)格組成一個顆粒.為清楚顯示顆粒的組成及分布情況，此模型的孔隙率設(shè)定為0.01、飽和度為0.5，因此代表水、空氣的顆粒數(shù)都為3個.限于篇幅，挑選三次隨機生成的模型來進行說明，如圖2所示，綠色代表的是礦物，藍色代表的是水，粉紅色代表的是空氣.圖2a是某次隨機生成的整體模型，圖2b是與之對應(yīng)的水顆粒、空氣顆粒的分布情況，圖2c、圖2d是其它兩次隨機模型中水顆粒、空氣顆粒的分布.圖3為數(shù)值模擬中建立的含濕孔隙巖石模型，上側(cè)為巖石整體模型，下側(cè)為礦物、水、空氣的分布.

3 結(jié)果及討論

圖3 含濕孔隙巖石模型Fig.3 Wet porous rock model and the parts

如上所述，孔隙的分布會影響試驗數(shù)值的分布，因此對于相同的一組輸入?yún)?shù)，多次試驗的結(jié)果會有分散性.所以，單次的計算結(jié)果不具有代表性，要用Monte Carlo法生成多個模型進行計算.在每一種固定的孔隙率及飽和度下，隨機生成200個模型進行計算，每次試驗得到的Keff值均繪出.圖4（a—d）分別是網(wǎng)格為4萬、8萬、12萬、15萬時的Keff值分布情況，圖中的f軸表示孔隙率，數(shù)值依次取0.05、0.1、0.15、0.2、0.25、0.3、0.4；s軸表示飽和度，數(shù)值由0遞增至1；Keff軸表示巖石有效熱導(dǎo)率.Km、Kw、Ka的取值分別為3.3、0.6、0.0274W·m－1·K－1.

由圖4可見：（1）當孔隙率和飽和度固定時，Keff的分布有起伏（即Keff不同），說明孔隙的分布、結(jié)構(gòu)等對有效熱導(dǎo)率有影響；（2）巖石的孔隙率對有效熱導(dǎo)率Keff有顯著影響，Keff隨著孔隙率的增大而減??；（3）網(wǎng)格的數(shù)量對計算結(jié)果有影響，網(wǎng)格較少時，計算的Keff值起伏較大、分布較分散，網(wǎng)格較多時，計算的Keff值起伏較小，分布相對集中.圖5更好地說明了網(wǎng)格數(shù)目對Keff集中程度的影響.

圖5只列出了孔隙率為0.3，飽和度分別取0、0.5、1時的Keff值分布圖.由于采用了隨機網(wǎng)格模型，試驗所得數(shù)值較為符合高斯分布，高斯分布的均值對應(yīng)不同比例下的Keff的平均值，圖中顯示固定孔隙率及飽和度下的數(shù)值分布與高斯曲線擬合較好.對比圖5（a—d）四幅圖可見：數(shù)值試驗中在固定孔隙率及飽和度的情況下，網(wǎng)格數(shù)越多可以獲得相對越集中的數(shù)值，峰值越明顯.

圖4 （a）4萬網(wǎng)格時 Keff值分布；（b）8萬網(wǎng)格時 Keff值分布；（c）12萬網(wǎng)格時 Keff值分布；（d）15萬網(wǎng)格時 Keff值分布Fig.4 （a）Distribution of Keff，40000grids；（b）Distribution of Keff，80000grids；（c）Distribution of Keff，120000grids；（d）Distribution of Keff，150000grids

圖5 孔隙率為0.3時，不同飽和度及網(wǎng)格數(shù)目下的Keff值分布Fig.5 Gaussian distribution of Keffwhen porosity is 0.3with different degree of saturation and the number of grids

本文采取固定孔隙率和飽和度的方法來分析不同情況下網(wǎng)格數(shù)目對誤差的影響.圖6a為只含固體礦物時，數(shù)值計算結(jié)果的相對誤差隨網(wǎng)格數(shù)的變化圖，圖6（b—d）是孔隙率為0.3、飽和度分別為0、0.5、1時有效熱導(dǎo)率的標準偏差隨網(wǎng)格數(shù)目的變化圖.從圖6a可以看出相對誤差隨網(wǎng)格數(shù)的增多而減小，當網(wǎng)格數(shù)目大于10萬時，相對誤差雖然還有減小的趨勢，但變化已相對平緩；網(wǎng)格為15萬時，數(shù)值計算出的Keff的相對誤差為0.007%，這同時也驗證了程序的可靠性及正確性.圖6（b—d）中標準偏差的分布顯示了類似的趨勢：標準偏差隨網(wǎng)格數(shù)目的增多而減小，即網(wǎng)格數(shù)目越多數(shù)值計算出的Keff的分散度越小、試驗穩(wěn)定性越好；當網(wǎng)格增多到一定程度時，標準偏差變化趨于平緩.如圖所示，網(wǎng)格數(shù)目對誤差影響的臨界值約為1×105.以上分析說明，為保證誤差穩(wěn)定在一定范圍內(nèi)，網(wǎng)格數(shù)不能少于1×105，而網(wǎng)格繼續(xù)增多雖然可以進一步減小誤差，但效果已不明顯.因此，在數(shù)值計算中既為了節(jié)省計算時間又保證精度，可以把網(wǎng)格劃為10萬左右.

圖7為網(wǎng)格數(shù)為12萬時，孔隙率為0.3、飽和度分別為0、0.5、1時，顆粒數(shù)為500、1000、2000、5000、10000、20000、30000時的標準偏差分布圖.圖7（a—c）顯示了相同的趨勢：標準偏差隨顆粒數(shù)目的增多而減小；當顆粒數(shù)目大于1×104時，標準偏差隨顆粒數(shù)目的變化已趨于平緩.以上分析說明，巖石試件體積與礦物顆粒平均體積之比對標準偏差影響的臨界值約為1×104.據(jù)此可解釋對小尺寸試樣實際測量時的尺度效應(yīng)，試件要足夠大、或者說包含了足夠多的礦物顆粒和孔隙時，測量結(jié)果才有代表性.

用熱導(dǎo)率儀測量試件的熱導(dǎo)率時，由于試件尺寸有限，因此不同試件測量結(jié)果會有差異.一般由于時間、經(jīng)費的限制，實際測量數(shù)目有限，而數(shù)值計算方法，如同圖7所顯示的，不但能夠得到熱導(dǎo)率的平均值，而且能對有限尺寸（包含有限礦物顆粒和孔隙）的試件測量的誤差分布有所了解.因此計算方法提供了更多的信息，也為我們提供了消除試件大小尺度效應(yīng)的途徑.

圖6 誤差隨網(wǎng)格數(shù)目變化的分布圖Fig.6 Error of Keffversus number of grids

圖7 標準偏差隨顆粒數(shù)目變化的分布圖Fig.7 Standard deviation of Keffversus number of grains

用熱導(dǎo)率儀測量試件的熱導(dǎo)率時，對試件的尺寸有嚴格要求.實驗試件最大一般為7.85×104mm3（美國EKO公司生產(chǎn)的 HC－110（190）熱導(dǎo)率儀對試件的尺寸要求為直徑51～63.5mm，厚度為0～25mm，最大約78.5×104mm3，精度為5%）.據(jù)圖7可知，測量含濕孔隙巖石的有效熱導(dǎo)率時，為了保證標準偏差穩(wěn)定在5%以內(nèi)，巖石試件的體積與礦物顆粒的體積之比應(yīng)大于2000.因此，只有當?shù)V物粒徑小于4mm時，熱導(dǎo)率儀的測量誤差才可達5%以下.而對于礦物粒徑大于4mm的巖石來說，試件存在尺度效應(yīng)，單個試件的測量結(jié)果則可能出現(xiàn)較大的偏差，要對大量試件測量和求平均值才能減小誤差.而我們的計算方法則不受試件大小或礦物顆粒大小的限制.

圖8 數(shù)值計算結(jié)果和實驗數(shù)據(jù)［26］的對比Fig.8 Comparison of the experimental results［26］and the numerical calculation

圖8給出了本數(shù)值計算結(jié)果和實驗數(shù)據(jù)［26］的對比，各組分的熱導(dǎo)率Km、Kw、Ka按文獻分別取1.856、0.574、0.025W·m－1·K－1.由圖可知，本文數(shù)值計算的結(jié)果和實驗數(shù)據(jù)基本一致，各測點的平均偏差為2.4%，最大偏差為7.37%.這說明本文所述的孔隙巖石模型是合理的，計算方法可行.

圖9給出了孔隙率、飽和度對有效熱導(dǎo)率的影響.整體而言，隨著孔隙率的增加，有效熱導(dǎo)率逐漸減??；當孔隙率較小時（約小于0.15），飽和度對有效熱導(dǎo)率的影響較小，孔隙率對有效熱導(dǎo)率的影響特別大；當孔隙率較大時，飽和度對有效熱導(dǎo)率的影響趨于明顯；孔隙率不變時，有效熱導(dǎo)率隨飽和度的增大而增大.由圖9可知，在已知模型孔隙率和飽和度的情況下，可直接在圖中確定出巖石的有效熱導(dǎo)率.模型中礦物、水、空氣的熱導(dǎo)率分別取3.3、0.6、0.0274W·m－1·K－1.

圖9 孔隙率、飽和度對Keff的影響Fig.9 Effects of the porosity and the degree of saturation on Keff

這樣，為快速方便確定巖石的有效熱導(dǎo)率，可以模擬計算得到含濕孔隙巖石的孔隙率及飽和度與有效熱導(dǎo)率的關(guān)系圖；固體礦物、水、空氣所占比例與有效熱導(dǎo)率關(guān)系圖，然后根據(jù)巖石的孔隙率及飽和度、或者是固體礦物、水、空氣所占比例，從圖中確定出巖石有效熱導(dǎo)率的范圍.

此有限元方法可以較好地模擬含濕孔隙巖石內(nèi)部的熱流情況.圖10（a—d）分別給出四個不同位置（h／5、2h／5、3h／5、4h／5）截面的x、y、z 三個方向的熱流分布情況，可以看出熱流的分布顯示出了明顯的不均勻性，這說明熱流的分布同樣受孔隙分布情況的影響.根據(jù)能量平衡，流過每個截面的總熱流是一個定值，但由于孔隙的存在，導(dǎo)致每個截面的熱流分布不均勻.

4 結(jié) 論

本文應(yīng)用隨機剖分單元、隨機指定材料屬性的有限元前處理方法，對含濕孔隙巖石的有效熱導(dǎo)率進行了數(shù)值模擬計算，并與實驗數(shù)據(jù)進行了對比.結(jié)果表明：（1）用有限元方法和Monte Carlo方法進行含濕孔隙巖石的熱力學(xué)性質(zhì)數(shù)值模擬是高效可行的；（2）本數(shù)值方法可以應(yīng)用于各類多孔介質(zhì)的有效熱導(dǎo)率預(yù)測上，也可以推廣到更為復(fù)雜的非均勻多孔介質(zhì)的研究上，從而進一步認識多孔介質(zhì)中的傳熱過程，為工程計算提供指導(dǎo)；（3）數(shù)值計算結(jié)果表明，孔隙巖石的傳熱不僅取決于礦物、水、空氣本身的熱物理特性，而且在很大程度上還與孔隙分布情況和孔隙結(jié)構(gòu)有關(guān)，在局部熱性質(zhì)是不均勻的，但整個巖體宏觀可以用一個熱導(dǎo)率描述；（4）試樣尺度會對巖體宏觀熱導(dǎo)率測量或網(wǎng)格數(shù)目對巖體宏觀熱導(dǎo)率計算造成誤差，這種誤差隨著實測時試件的增大或計算時網(wǎng)格數(shù)目的增加而減??；（5）在數(shù)值計算中，為了節(jié)省計算時間、保證計算精度，可以把網(wǎng)格劃為10萬左右；（6）通過我們的計算可知，用熱導(dǎo)率儀測量巖石試件熱導(dǎo)率時，只有當?shù)V物粒徑小于4 mm時，才能保證誤差在5%以內(nèi)，當?shù)V物粒徑大于4mm時，要對大量試件測量和求平均值才能減小誤差.本文僅簡單假定孔隙巖石骨架各種礦物的熱導(dǎo)率相差不大，可以近似視為一個常數(shù)，不過沿用本文的方法，不難實現(xiàn)對組分已知的多種礦物復(fù)合孔隙巖石的模擬.

本方法結(jié)合CT＿LBM技術(shù)［27］可以推廣到巖石的力學(xué)、電學(xué)等其它物理性質(zhì)的研究，也可以推廣到人工復(fù)合材料的研究、設(shè)計，成為確定復(fù)合材料性質(zhì)的重要手段.本文的方法在常溫常壓下得到了證實，進一步，如果我們知道各種礦物高壓高溫下的物性，以及巖石中各種礦物所占的比例，那么也不難計算多礦物巖石在高溫高壓下的性質(zhì).高溫高壓下實驗測試成本更高，以礦物物性測量為基礎(chǔ)，計算組分已知的巖石的物性，本文的方法應(yīng)該有廣闊的應(yīng)用天地.

（References）

［1］楊淑貞，張文仁，沈顯杰.孔隙巖石熱導(dǎo)率的飽水試驗研究.巖石學(xué)報，1986，2（4）：83－91.Yang S Z，Zhang W R，Shen X J.Experimental research on the thermal conductivity of water－saturated porous rocks.Acta Petrologica Sinica （in Chinese），1986，2（4）：83－91.

［2］Bear J著.李競生，陳崇希譯.多孔介質(zhì)流體動力學(xué).北京：中國建筑工業(yè)出版社，1983.Bear J.Translated by Li J Z，Chen C X.Dynamics of Fluids in Porous Media（in Chinese）.Beijing：China Architecture and Building Press，1983.

［3］李大心.國外巖石熱導(dǎo)率的研究動向.地質(zhì)科技情報，1985，4（2）：12－20.Li D X.Research trends of thermal conductivity of rocks abroad.Geological Science and Technology Information （in Chinese），1985，4（2）：12－20.

［4］Miller M N.Bounds for effective electrical，thermal，and magnetic properties of heterogeneous materials.Mathematical Physics，1969，12（10）：1988－2004.

［5］Carbonell R G，Whitaker S.Heat and mass transfer in porous media.／／Proceedings of the NATO Advanced Study Institute on Mechanics of Fluids in Porous Media.Newark，USA，1982；121－198.

［6］Hadley G R.Thermal conductivity of packed metal powders.Int.J.Heat and Mass Transfer，1986，29（6）：909－920.

［7］Kaviany M.Principle of Heat Transfer in Porous Media.New York：Springer－Verlag，1995.

［8］Katz A J，Thompson A H.Fractal sandstone pores：implications for conductivity and pore formation.Physical Review Letters，1985，54（12）：1325－1328.

［9］Thompson A H，Katz A J，Krohn C E.The microgeometry and transport properties of sedimentary rock.Advances in Physics，1987，36（5）：625－694.

［10］李小川，施明恒，張東輝.非均勻多孔介質(zhì)中導(dǎo)熱過程.東南大學(xué)學(xué)報（自然科學(xué)版），2005，35（5）：761－766.Li X C，Shi M H，Zhang D H.Heat conduction process in non－uniform porous media.Journal of Southeast University（Natural Science Edition）（in Chinese），2005，35（5）：761－766.

［11］Bogaty H，Hollies N R S，Milton H M.Some thermal properties of fabrics：Part Ⅰ.The effect of arrangement.Text.Res.J.，1957，27（6）：445－449.

［12］Hollies N R S，Bogarty H.Some thermal properties of fabrics PartⅡ.The effect of water content.Text.Res.J.，1965，35（2）：187－190.

［13］Singh A K，Singh R，Chaudhary D R.Heat conduction through moist soils at different temperatures.Pramana，1988，31（6）：523－528.

［14］Pande R N，Kumar V，Chaudhary D R.Thermal conduction in a homogeneous two－phase system.Pramana，1984，22（1）：63－70.

［15］Pande R N，Chaudhary D R.Thermal conduction through loose and granular two－phase materials at normal pressure.Pramana，1984，23（5）：599－605.

［16］Zhang H F，Ge X S，Ye H.Effectiveness of the heat conduction reinforcement of particle filled composites.Modelling Simul.Mater.Sci.Eng.，2005，13（3）：401－412.

［17］Zhang H F，Ge X S，Ye H.Randomly mixed model for predicting the effective thermal conductivity of moist porous media.J.Phys.D：Appl.Phys.，2006，39（1）：220－226.

［18］Zhang H F，Ge X S，Ye H，et al.Heat conduction and heat storage characteristics of soils.Applied Thermal Engineering，2007，27（2－3）：369－373.

［19］Wang M R，Wang J L，Pan N，et al.Three－dimensional effect on the effective thermal conductivity of porous media.Phys.D：Appl.Phys.，2007，40（1）：260－265.

［20］Tiak D，Delkumburewatte G B.The influence of moisture content on the thermal conductivity of a knitted structure.Measurement Science and Technology，2007，18（5）：1304－1314.

［21］Aurangzeb，Maqsood A.Modeling of the effective thermal conductivity of consolidated porous media with different Saturants：A test case of Gabbro rocks.International Journal of Thermophysics，2007，28（4）：1371－1386.

［22］景惠敏.數(shù)值模擬在地震海嘯、港口風浪等地學(xué)問題上的應(yīng)用［博士論文］.北京：中國科學(xué)院研究生院地球科學(xué)學(xué)院，2011.Jing H M.Numerical simulation applied in evaluation of tsunami and wave hazard in a harbor ［Ph.D.thesis］（in Chinese）.Beijing：College of Earth Science，Graduate University of Chinese Academy of Sciences，2011.

［23］王勖成.有限單元法.北京：清華大學(xué)出版社，2003.Wang X C.Finite Element Method （in Chinese）.Beijing：Tsinghua University Press，2003.

［24］Bouguerra A.Prediction of effective thermal conductivity of moist wood concrete.J.Phys.D：Appl.Phys.，1999，32（12）：1407－1414.

［25］Vermia L S，Shrotriya A K，Singh R，et al.Prediction and measurement of effective thermal conductivity of three phase systems.J.Phys.D：Appl.Phys.，1991，24（9）：1515－1526.

［26］Somerton W H. Thermal Properties and Temperature Related Behavior of Rock／Fluid Systems. New York：Elsevier，1992.

［27］Zhu B J，Cheng H H，Qiao Y C，et al.Porosity and permeability evolution and evaluation in anisotropic porosity multiscale－multiphase－multicomponent structure.Chinese Science Bulletin，2012，57（4）：320－327.