劉鳳臣,于學峰,賈超,杜圣賢,宋香鎖,楊斌,張尚坤
(1.山東省地質(zhì)科學研究院,國土資源部金礦成礦過程與資源利用重點實驗室,山東省金屬礦產(chǎn)成礦地質(zhì)過程與資源利用重點實驗室,山東 濟南 250013;2.山東大學土建與水利學院,山東 濟南 250013)
古生物化石保護得到人們越來越多的關(guān)注,恐龍化石在古生物化石中最為重要[1],然而它同樣經(jīng)歷著大自然的風化破壞,研究發(fā)現(xiàn),溫度和水是其中的兩大主要因素[2-4]。
巖體溫度場與滲流場耦合分析是巖體水力學科中重要的研究內(nèi)容,在大量的人類工程及開采活動中,如地下水開采、石油開采、地熱開發(fā)、地質(zhì)災(zāi)害治理、大型水利工程等領(lǐng)域得到廣泛的應(yīng)用[5-18]??铸埢皣鷰r受日照影響,在其表面和內(nèi)部形成溫度場,同時化石及圍巖內(nèi)部本身存在一定孔隙水,形成滲流場。該文將巖體水力學中溫度場與滲流場的耦合關(guān)系運用到恐龍化石及圍巖中,研究恐龍化石及圍巖中溫度場與滲流場如何相互影響、相互作用對化石及圍巖造成風化破壞。
該文滲流場采用多孔介質(zhì)布里克曼方程[19-24],與溫度場相結(jié)合。
(1)多孔介質(zhì)模型傳熱的數(shù)學模型為:
(1)
(ρCp)eq=θpρpCρ·p+(1-θp)ρCρ
(2)
keq=θp+(1-θp)k
(3)
式中:ρ—液體密度(kg/m3);Cp—液體等壓熱容(J/K);k—熱傳導(dǎo)系數(shù)(W/(m2K));u—速度場;Y—比熱比,等壓熱容與等容熱容之比;(ρCp)eq—等效體積等壓熱容;θp—材料的體積分數(shù);kp—固體材料的熱傳導(dǎo)系數(shù);ρp—固體材料密度(kg/m3);Cρ·p—固體材料的比熱容(J/(kg·K))。
(2)多孔介質(zhì)流動的數(shù)學模型是由連續(xù)性方程和動量平衡方程組合而成,它們一起構(gòu)成了Brinkman方程:
(4)
(5)
式中:μ—流體的動態(tài)粘滯系數(shù)(Pa·s);u—流體的速度;ρ—流體的密度(kg/m3);p—壓力(Pa);εp—孔隙率;k—滲透系數(shù)(m/d);Qbr—質(zhì)量源;F—體積力。
(3)溫度場與滲流場的耦合:在溫度場的計算中,速度采用布里克曼方程計算得出的速度;而在滲流場計算中,體積力計算公式為:
F=ρgβ(T-Tc)
(6)
式中:T—溫度場計算得出的溫度;Tc—初始溫度。
根據(jù)恐龍化石的簡化模型和圍巖的簡化模型,建立熱流耦合數(shù)值模型如圖1所示[2]。其邊界條件為:內(nèi)部初始溫度采用諸城地區(qū)多年平均氣溫12℃,左右兩邊界對稱,前后兩邊界對稱,上表面裸露于空氣中。
圖1 熱流耦合數(shù)值模型 (據(jù)杜圣賢等,2015)
隨著時間變化,外界溫度變化過程如圖2所示。其中包括溫度上升和下降兩個階段。0~30000s為溫度上升階段,30000~60000s為溫度下降階段,之后溫度又上升。
圖2 時間-溫度變化曲線
熱流耦合模型參數(shù)如表1所示。
隨著外界溫度升高,分別模擬了0.5h(1800s),1h(3600s),2h(7200s),3h(10800s),6h(21600s),8h(28800s)的變化情況(圖3)。圖中等值線為滲流場流速梯度的變化情況,矢量圖為速度矢量的變化情況??梢钥闯?,在外界溫度變化的前8個小時(0~28800s),隨著外界溫度不斷升高,表面流場速度梯度迅速變化,速度矢量也不斷變化。在0.5個小時(1800s)時,化石表面流場速度梯度最大,流速也是化石兩端最大,而此時圍巖表面流速梯度普遍變化較小。這是因為當外界溫度變化時,由于化石裸露面積較大,化石內(nèi)部孔隙水運動劇烈導(dǎo)致的;第1個小時(3600s)時,明顯可以看出,在恐龍化石和圍巖接觸處,流場流速梯度和流速均明顯大于其他部分,這是由于化石和圍巖的各力學性質(zhì)和水力學性質(zhì)不同,導(dǎo)致它們之間存在相互作用,故在外界條件變化時,化石和圍巖接觸處受影響最大,此處風化速度會比其他地方快。
隨著外界溫度溫度開始降低時,分別模擬了13h(46800s),15h(54000s),17h(61200s),24h(86400s)的變化情況(圖4)??梢钥闯?,外界溫度開始降低后,最先受影響的也是化石裸露區(qū),然后擴展到圍巖部分。15(54000s)~17h(61200s),通過流速梯度和流速大小變化圖形可以看出,外界溫度降低的情況下,化石兩端和內(nèi)部速度梯度和速度大小非常不均勻,這種情況更會加速化石的風化。根據(jù)總體變化趨勢可以發(fā)現(xiàn),圍巖部分內(nèi)部空隙水隨外界溫度變化而變化持續(xù)的時間更長,故當外界條件,特別是水和溫度持續(xù)變化時,圍巖的風化速率可能會大于化石??傊?,隨外界溫度的不斷變化,化石和圍巖的表面流場劇烈變化,孔隙水的運動會破壞化石和圍巖的內(nèi)部結(jié)構(gòu),導(dǎo)致其強度降低;水的作用會帶動內(nèi)部化學離子的運動,加速化學風化的進行。
表1 熱流耦合模型參數(shù)
圖3 溫度上升時恐龍化石及圍巖表面滲流場變化圖
圖4 溫度降低時恐龍化石及圍巖表面滲流場變化圖
隨著外界溫度升高時,分別模擬了1h(3600s),3h(10800s)的變化情況(圖5)??梢钥闯?,外界溫度變化1h(3600s)時,內(nèi)部流場在恐龍化石和圍巖接觸處變化最大,在圍巖內(nèi)部也有幾處變化較大區(qū)域,但比例較?。浑S著外界溫度的不斷升高,到第3個小時(10800s)時,內(nèi)部流場在化石和圍巖中流速均較大,即內(nèi)部空隙水運動劇烈。
隨著外界溫度開始降低直至再次升高時,分別模擬了12h(43200s),18h(64800s)的變化情況(圖6)??梢钥闯觯S著外界溫度的降低,第12個小時(43200s)時,恐龍化石及圍巖內(nèi)部流速都有所降低,但部分區(qū)域仍有較大的流速梯度;同時,在由溫度升高到溫度降低的過程中,其內(nèi)部孔隙水的運動方向發(fā)生了變化。之后,隨外界溫度又回升,內(nèi)外溫差較大,孔隙水運動更加劇烈,到第18個小時(64800s)時,整個化石內(nèi)部空隙水流動劇烈,圍巖內(nèi)孔隙水運動也開始深入內(nèi)部。
隨外界溫度變化化石內(nèi)部滲流場變化,分別模擬了9h(32400s),24h(86400s)的變化情況(圖7)。可以看出,隨著外界溫度變化,化石內(nèi)部流速梯度較大,在化石裸露區(qū)兩側(cè),空隙水向外界流出;而在化石與圍巖接觸區(qū),孔隙水有穿過圍巖區(qū)向上表面運動的趨勢,而這種運動會加速化石與圍巖接觸區(qū)化石和圍巖的風化;在化石中間部分也出現(xiàn)了較大的流速區(qū),表明化石內(nèi)部隨著外界溫度的變化,內(nèi)部空隙水流動會出現(xiàn)集中區(qū),循環(huán)往復(fù),這部分化石的強度會降低,出現(xiàn)微小的裂縫,隨后裂縫部分水的運動會更加集中,風化速度加快。
圖5 溫度升高時恐龍化石及圍巖內(nèi)部滲流場變化圖
圖6 溫度降低至升高時恐龍化石及圍巖內(nèi)部滲流場變化圖
圖7 化石內(nèi)部滲流場分布圖
圖8 化石及圍巖熱流耦合內(nèi)部溫度及流速變化圖
圖8為恐龍化石及圍巖熱流耦合的二維切片,分別模擬了1h(3600s),20h(72000s)的變化情況??梢钥闯觯跏茧A段溫度場和滲流場相互作用,使得化石和圍巖表面產(chǎn)生孔隙水運動;隨著外界溫度的持續(xù)變化,化石和圍巖內(nèi)部孔隙水運動范圍加大,但可以看出,在化石和圍巖接觸處孔隙水運動速度遠大于其他部分;隨著外界溫度的升降變化,化石和圍巖內(nèi)部空隙水運動的方向會發(fā)生變化;由于孔隙水的運動,影響了溫度分布,20h(72000s)時的溫度云圖非常明顯,這會更加加劇化石和圍巖內(nèi)部的溫度差,使得內(nèi)部產(chǎn)生溫度應(yīng)力,加速風化進程;同時溫度的分布也影響了內(nèi)部空隙水的運動速度。
(1)從恐龍化石及圍巖的內(nèi)部結(jié)構(gòu)出發(fā),以恐龍化石及圍巖作為巖體,介紹了恐龍化石風化過程中溫度場與滲流場耦合作用的基本理論,并且建立了化石及圍巖熱流耦合作用的多孔介質(zhì)模型,模擬了溫度場與滲流場在恐龍化石及圍巖內(nèi)部及表面的相互作用,描述了恐龍化石在溫度和水的相互作用下的破壞機制。
(2)實驗結(jié)果顯示,化石及圍巖受日照影響產(chǎn)生的溫度會影響其內(nèi)部孔隙水的變化,影響滲流場的變化;同時,化石及圍巖的孔隙水也會參與到化石及圍巖內(nèi)部的熱量交換中,影響化石及圍巖的溫度場的分布。兩者相互作用,共同對恐龍化石的破壞產(chǎn)生影響。
(3)通過對恐龍化石及圍巖進行熱流耦合作用的模擬研究,對恐龍化石的保護提供了一定的科學依據(jù)。今后可借助巖石力學、流體力學、傳熱學等相關(guān)學科,開展化石及圍巖溫度場—滲流場—應(yīng)力場三場之間耦合作用研究,建立三場耦合作用的數(shù)學模型,為化石保護建立完整的物理場防護體系。