張尚坤,于學峰，賈超,杜圣賢,宋香鎖，陳誠，劉鳳臣，陳軍，陳文芳

(1.山東省地質(zhì)科學研究院，國土資源部金礦成礦過程與資源利用重點實驗室，山東省金屬礦產(chǎn)成礦地質(zhì)過程與資源利用重點實驗室，山東 濟南 250013;2.山東大學土建與水利學院，山東 濟南 250013)

1 概述

熱應力是引起巖石破壞的一個重要因素，由熱應力引起的巖石破裂稱為熱破裂。根據(jù)熱傳導理論，對于同樣的物體，在相同加熱條件下，如果初始條件和邊界條件不同，即加熱方式不同，由于熱傳導在時間和空間上的遲滯效應，在物體表面和內(nèi)部的溫度分布將存在較大的差異，這將導致物體中應力的集中和分布的差異[1-3]。在研究陶瓷材料的熱性能過程中研究人員發(fā)現(xiàn)熱沖擊在物體表面產(chǎn)生的熱應力要大于均勻加熱產(chǎn)生的熱應力。Moon V等人[4-6]在研究蒸汽壓力對沉積巖石熱破裂的影響實驗中，得到了與此相同的結論?？铸埢療崞茐闹饕窃谔柟庹盏臈l件下，由于受熱不均勻及日照長久往復循環(huán)作用，使得化石表面可能產(chǎn)生裂縫，從而加速了恐龍化石風化損害。

引起熱應力的基本條件是在約束下有溫度的變化。產(chǎn)生熱應力的約束條件大致可以歸納為：外部變形的約束，相互變形的約束，內(nèi)部各區(qū)域之間變形的約束。就巖石而言，各礦物之間的物理力學性質(zhì)有較大的差異，在溫度變化時，由于內(nèi)部各區(qū)域熱變形的不協(xié)調(diào)而產(chǎn)生熱應力；當有地溫梯度存在時，即使沒有外部約束，但由于相互變形的約束，也將產(chǎn)生較大的熱應力。由這些因素引起的熱應力很容易導致巖石中弱介質(zhì)的破壞，從而改變巖體力學性質(zhì)和力學行為[7-8]。恐龍化石挖掘過程中，大量的化石暴露地表，很多化石的一部分出露地表，一部分仍埋藏于地表下[9-18]。由于空氣和地表下巖石介質(zhì)不同，受熱條件不同，比熱容和熱傳導特性各異，造成化石暴露在空氣中的部分與埋藏在地表下的部分產(chǎn)生溫度差，受熱膨脹不均勻而產(chǎn)生破壞應力。

恐龍骨骼化石是被埋葬的恐龍骨骼經(jīng)過自然地質(zhì)作用，經(jīng)鈣化、碳化、硅化、礦化形成的巖石，而保留了原骨骼的形狀和結構。礦物成分包含多種礦物顆粒，在受熱作用下，由于這些礦物顆粒在熱力學上表現(xiàn)出的差異(熱膨脹各向異性、熱膨脹不均勻性等)使得化石內(nèi)部結構發(fā)生相應的變化，當這些變化超過化石自身的內(nèi)部結合力限度時，就會在化石內(nèi)部產(chǎn)生微裂縫，從而引起化石物性參數(shù)(如孔隙度、滲透率等)的變化，加速恐龍化石的風化。

該試驗主要針對日照產(chǎn)生在恐龍化石件各處不同溫度差異分布情況下，恐龍化石試件的物理力學受力變形特性進行研究。通過有限差分軟件FLAC3D進行數(shù)值模擬[19-23]，揭示日光輻射造成的恐龍化石溫度差異分布對化石風化的影響。

2 研究的理論基礎

通常熱力耦合方程可以描述為：

(1)

σij,ii+Fbi=ρiii

(2)

εij=(ui,j+uj,i)/2

(3)

σij=λεmmδij+2Gεij-βΔΤδij

(4)

β=(3λ+2G)α

(5)

式中：σij和εij為應力和應變項；Fbi為體力；β為熱應力系數(shù)；ΔΤ為溫度改變量，即ΔΤ=T-T0；δij為Kronecker函數(shù)；Q，λ，G，α，k，ρ和c分別為熱生成率、拉梅常數(shù)、剪切模量、熱膨脹系數(shù)、熱傳導系數(shù)、密度和比熱。

嚴格地說，這些方程是相互關聯(lián)的，應該一起求解。然而，對大多數(shù)的工程問題來說，應變對溫度分布的影響比較小，且慣性變量項也非常小，即式(1)中的和式(2)中的可以忽略不計。這樣一來，溫度場的計算就純屬熱傳導范圍，應力場的計算中包含了熱應力項。所以在計算時就可以先計算出溫度場，然后再計算由溫度和外荷載共同引起的應力和變形等。

基于上述分析，該數(shù)值模型中有關熱力耦合的求解就可以分兩步進行，即先根據(jù)溫度條件計算出模型中的溫度場，然后根據(jù)力學條件求解應力場并進行相應的破壞處理。即每個單元滿足如下的熱傳導方程和應力場方程。

熱傳導方程為：

(6)

Τ(p,t)|p∈S1=φ(p,t) (邊界條件為S1)

(7)

(8)

(9)

Τ|t=t0=Τ(p,t0) (在t0時刻)

(10)

應力-應變場方程為：

σij,j+Fbi=0 (平衡方程)

(11)

εij=(ui,j+uj,i)/2 (幾何方程)

(12)

σij=λεkkδi,j+2Gεij-βΔΤδij(物理方程)

(13)

(14)

(15)

β=(3λ+2G)α(邊界條件為S1)

(16)

巖石材料在外載作用下，當應力達到其所受力強度極限后，其應力會隨著變形的增加而降低,這種現(xiàn)象通常被稱之為“應變軟化”。由于應變軟化的影響，結構變形穩(wěn)定發(fā)展到一定程度之后，會突然失去變形的穩(wěn)定性，導致結構的動態(tài)破壞失效。另一方面，破裂的巖石仍具有一定的殘余承載力。針對巖石應變軟化破壞特性，F(xiàn)LAC3D提供了一種具有應變軟化/硬化特性的本構模型—Strain-Hardening/ Softening Mohr-Coulomb Model[20,24]，基于摩爾庫倫的破壞準則，破壞屈服后摩擦角、粘聚力、剪脹角、抗拉強度等能以一定的規(guī)律遞減，最終達到一個穩(wěn)定的強度值。該構模型如圖1所示；屈服后粘聚力、摩擦角變化如圖2所示。

圖1 應變軟本構模型應力應變曲線

圖2 屈服后粘聚力(a)、摩擦角(b)的變化曲線示意圖

3 模擬試驗條件的設定

模擬試驗的主要目的是研究不同溫度場下，恐龍化石試件單軸壓縮變形特性[25]。為了方便研究，采用Strain-Hardening/Softening Mohr-Coulomb Model本構模型，將恐龍化石試件設定為直徑50mm，長100mm的圓柱體試塊，模型共計6144個單元，6897個節(jié)點(圖3)。

圖3 模擬計算模型

由于地表介質(zhì)不同，吸收的陽光輻射也不同，導致日照作用下地表溫度與氣溫存在較大差異，實測顯示：一般情況下，巖石的地表溫度要大大高于氣溫，在日照直射的區(qū)域，當空氣溫度為40℃時，巖石地表最高溫度可以達到50℃～60℃。而埋藏在地下1m范圍的物體，其最低溫度可以只有十幾攝氏度。因此，開挖出露的恐龍化石在日照輻射作用下，化石暴露在空氣中的部分與埋藏在土層中的部分存在較大溫差，化石內(nèi)部存在溫度梯度場，在溫度不均勻的熱膨脹作用下，會使得試塊中產(chǎn)生易導致化石破壞的溫度拉應力,巖石的熱物理性質(zhì)如表1所示。根據(jù)恐龍化石賦存情況，設計采用7種不同的工況進行對比試驗研究，具體情況如表2所示。

表1 巖石熱物理性質(zhì)

表2 恐龍化石地下與地表部分溫度比較

4 試驗模擬結果與分析

(1)溫度應力分布場

不同溫度梯度場作用下，圓柱體試件的溫度分布情況和其產(chǎn)生的溫度應力分布情況如圖4—圖10所示，對比溫度應力可以發(fā)現(xiàn)：隨著溫差逐漸增大，圓柱體試塊的上半部分的四周表面出現(xiàn)拉應力，在溫度差異值大于30℃后變得比較明顯。而在試件內(nèi)部卻為壓應力，應力分布不均勻容易導致塊體的表面的開裂。

(a)—溫度場分布；(b)—溫度拉應力分布圖4 工況1(溫差0℃)溫度場與溫度拉應力分布云圖

(a)—溫度場分布；(b)—溫度拉應力分布圖5 工況2(溫差10℃)溫度場與溫度拉應力分布云圖

(a)—溫度場分布；(b)—溫度拉應力分布圖6 工況3(溫差20℃)溫度場與溫度拉應力分布云圖

(a)—溫度場分布；(b)—溫度拉應力分布圖7 工況4(溫差30℃)溫度場與溫度拉應力分布云圖

(a)—溫度場分布；(b)—溫度拉應力分布圖8 工況5(溫差40℃)溫度場與溫度拉應力分布云圖

(a)—溫度場分布；(b)—溫度拉應力分布圖9 工況6(溫差50℃)溫度場與溫度拉應力分布云圖

(a)—溫度場分布；(b)—溫度拉應力分布圖10 工況7(溫差60℃)溫度場與溫度拉應力分布云圖

(2)不同工況單軸壓縮破壞結果

對不同溫度應力分布的試件進行單軸壓縮模擬試驗，試驗結果如圖11—圖17所示，試驗結果顯示隨溫度梯度增長塑性破壞區(qū)域變化不大，但是塑性流動區(qū)域，即破壞后變形滑動的區(qū)域，逐漸從中部擴展到試件的兩端表面。這個過程在溫度差大于30℃以后，表現(xiàn)得較為明顯。

(3)應力應變曲線變化分析

在對既有溫度應力的恐龍化石單軸加載試驗過程中，在模型上設置了4個應力監(jiān)測關鍵點(圖18)，用于觀察在試驗加載過程中，關鍵點的應力隨模型的軸向應變的變化關系，以及加載應力隨軸向應變的變化關系(圖19—圖25)。

(a)—塑性破壞區(qū)；(b)—剪切流動變形區(qū)圖11 工況1(溫差0℃)試件單軸壓縮剪切破壞單元分布

(a)—塑性破壞區(qū)；(b)—剪切流動變形區(qū)圖12 工況2(溫差10℃)試件單軸壓縮剪切破壞單元分布

(a)—塑性破壞區(qū)；(b)—剪切流動變形區(qū)圖13 工況3(溫差20℃)試件單軸壓縮剪切破壞單元分布

(a)—塑性破壞區(qū)；(b)—剪切流動變形區(qū)圖14 工況4(溫差30℃)試件單軸壓縮剪切破壞單元分布

(a)—塑性破壞區(qū)；(b)—剪切流動變形區(qū)圖15 工況5(溫差40℃)試件單軸壓縮剪切破壞單元分布

(a)—塑性破壞區(qū)；(b)—剪切流動變形區(qū)圖16 工況6(溫差50℃)試件單軸壓縮剪切破壞單元分布

(a)—塑性破壞區(qū)；(b)—剪切流動變形區(qū)圖17 工況7(溫差60℃)試件單軸壓縮剪切破壞單元分布

圖20 工況2監(jiān)測點應力應變曲線

圖21 工況3監(jiān)測點應力應變曲線

圖22 工況4監(jiān)測點應力應變曲線

從以上實驗結果可以看出，溫度應力場的存在對化石應力峰值抗壓強度的影響不大，但是對化石應力峰后強度有較為顯著的影響。對比溫度差異較小的工況1、工況2與工況3,4,5,6,7的應力應變曲線可知，溫度應力場較小時，監(jiān)測點應力峰后應力應變曲線較為平滑一致，而溫度應力場較大工況下，監(jiān)測點應力峰后的應力應變曲線會有較大的起伏波動，據(jù)此分析，波動是由于單元體產(chǎn)生較大塑性流動破壞造成的，說明塑性流動區(qū)隨溫度梯度增大而擴展。

圖23 工況5監(jiān)測點應力應變曲線

圖24 工況6監(jiān)測點應力應變曲線

圖25 工況7監(jiān)測點應力應變曲線

5 結論

(1)熱應力是引起化石體破裂風化的一個重要原因，化石作為一種特殊巖石其礦物成分、結構構造等存在很大的差別。在受熱作用下，化石內(nèi)部發(fā)生熱開裂現(xiàn)象，導致化石內(nèi)部結構發(fā)生變化。

(2)溫度梯度差異對化石的影響，主要表現(xiàn)在溫度的熱膨脹作用下，化石的表面出現(xiàn)拉應力區(qū)，由于化石本身是脆性材料，若溫度不斷循環(huán)變化，在化石表面就會形成長久往復的拉壓循環(huán)應力，極易造成化石表面產(chǎn)生微小裂縫，并使裂縫不斷擴張，形成空氣進出通道，加速化石風化破壞。

(3)溫度分布不均勻，往往還會加重一定初始地應力場中化石的受力狀況。埋藏在地層中的化石是在一定的地應力場中生成的，當部分被發(fā)掘暴露后，其本身有一定的應力釋放，內(nèi)部的應力分布是不均勻的。因此，在受外界溫度影響下，溫度熱應力的作用，很可能造成應力分布情況的惡化，使得局部應力過大，造成化石開裂。

(4)在化石發(fā)現(xiàn)挖掘中，暴露后的化石風化破壞明顯加快，通過實驗證明埋藏在地下的化石部分與暴露在空氣中的化石部分的溫度差異，是造成化石內(nèi)部熱應力分布不均勻的重要原因。另外，在受外部熱源(如太陽輻射)作用下，產(chǎn)生在化石表面與內(nèi)部溫度梯度分布，也是造成化石應力分布不均勻的原因，因此化石保存要盡量放在恒溫環(huán)境下，避免化石溫度往復變化。

參考文獻：

[1] 張淵,張賢,趙陽升.砂巖的熱破裂過程[J].地球物理學報,2005,48(3):656-659.

[2] 秦本東,何軍,諶倫建.石灰?guī)r和砂巖高溫力學特性的試驗研究[J].地質(zhì)力學學報,2009,15(3):253-261.

[3] 楊俊杰,黃思靜,張文正，等.表生和埋藏成巖作用的溫壓條件下不同組成碳酸鹽巖溶解過程的實驗模擬[J].沉積學報,1995,13(4):49-54.

[4] Moon V,Jajawardane J.Geomechanical and geochemical changes during early stages of weathering of Karamu Basalt,New Zealand[J].Engineering Geology,2004,74:57-72.

[5] White A F,Bullen T D.Differential rates of feldspar weathering in granitic regoliths[J]．Geochimica et Cosmochimica，2001,65(6):847-869.

[6] Matsukara Y,Hashizume K,Oguchi C T．Effect of microstructure and weathering on the strength anisotropy of porous rhyolite[J].Engineering Geology,2002,62:39-47.

[7] 肖洪天,周維垣.脆性巖石變形與破壞的細觀力學模型研究[J].巖石力學與工程學報，2001,20(2):151-155.

[8] 劉元高,周維垣,趙吉東，等.裂隙巖體損傷局部化破壞分岔模型及其應用[J].力學學報,2003,35(4):411-418.

[9] 杜圣賢，于學峰，宋香鎖，等.山東省古生物化石保護規(guī)劃研究[M].北京：地質(zhì)出版社，2017：1-84.

[10] 王申娜，陳樹清，張艷霞，等．山東“恐龍澗”創(chuàng)造三個“世界之最”[J]．中國地理雜志，2010(6):88-93.

[11] 陳樹清．山東諸城恐龍化石發(fā)掘新成果[J]．化石,2010(1):9-11．

[12] 趙喜進，李敦景，韓崗，等．山東的巨大諸城龍[J]．地球?qū)W報，2007,2(4):111-122．

[13] 季燕南．巨型山東龍的系統(tǒng)分類、生活習性與生態(tài)環(huán)境研究[J]．地學前緣，2010,17(1):378-385.

[14] 韓振玉,王經(jīng),王秀靜,等.萊陽恐龍谷地層古生物現(xiàn)狀與保護建議[J].山東國土資源,2015,31(10):36-39．

[15] 夏文杰,李秀華．恐龍埋藏環(huán)境及巖相古地理特征[M]．成都:四川科學技術出版社，1988．80-82．

[16] 柳永清,曠紅偉，姬書安，等．魯東諸城地區(qū)晚白堊世恐龍集群埋葬地沉積相與埋葬學初步研究[J]．地質(zhì)論評，2010.56(4):4457-4468.

[17] 楊群，王怡林,李朝真，等．原子光譜和X射線衍射研究祿豐恐龍化石與圍巖特征[J]．光譜學與光譜分析，2005，25(2)：299-301.

[18] 杜圣賢,張尚坤,于學峰,等.恐龍化石風化效應的TM耦合分析研究[J].山東國土資源,2015,31(10):65-70.

[19] 劉波,韓彥輝.FLAC原理實例與應用指南[M].北京:人民交通出版社，2004:79-90.

[20] 陳育民,徐鼎平.FLAC/FLAC3D基礎與工程實例[M].北京：中國水利水電出版社,2009：1-20.

[21] Fairhurst C,Cook N G W.The phenomenon of rock splitting parallel to the direction of maximum compression in the neighborhood of a surface crack[A].In:Proc.First Congress International Society for Rocks Mechanics[C].Lisbon:Lisbon Press,1966:15-65.

[22] Kemeny J,Cook N G W.Effective moduli no-linear deformation and strength of a cracked elastic solid[J].Int. J .Rock Mechanics. Min. Sci. & Geomech. Abstr.,1986,23(2):107-118.

[23] Kemeny J,Cook N G W.Micromechanics of deformation in rock[A].In:Toughening Mechanism in Quasi-Brittle Materials[C].Netherlands:Kluwer Academic Publishers,1991:62-78.

[24] 孫衛(wèi)軍,周維垣.裂隙巖體彈塑性—損傷本構模型[J].巖石力學與工程學報，1990,9(2):108-119.

[25] 傅晏,劉新榮,張永興，等．水巖相互作用對砂巖單軸強度的影響研究[J]．水文地質(zhì)工程地質(zhì)，2009(6)：54-58．