温度场和渗流场耦合对恐龙化石风化破坏的影响

刘凤臣，于学峰，贾超，杜圣贤，宋香锁，杨斌，张尚坤

(1.山东省地质科学研究院，国土资源部金矿成矿过程与资源利用重点实验室，山东省金属矿产成矿地质过程与资源利用重点实验室，山东 济南 250013；2.山东大学土建与水利学院，山东 济南 250013)

0 引言

古生物化石保护得到人们越来越多的关注，恐龙化石在古生物化石中最为重要[1]，然而它同样经历着大自然的风化破坏，研究发现，温度和水是其中的两大主要因素[2-4]。

岩体温度场与渗流场耦合分析是岩体水力学科中重要的研究内容，在大量的人类工程及开采活动中，如地下水开采、石油开采、地热开发、地质灾害治理、大型水利工程等领域得到广泛的应用[5-18]。恐龙化石及围岩受日照影响，在其表面和内部形成温度场，同时化石及围岩内部本身存在一定孔隙水，形成渗流场。该文将岩体水力学中温度场与渗流场的耦合关系运用到恐龙化石及围岩中，研究恐龙化石及围岩中温度场与渗流场如何相互影响、相互作用对化石及围岩造成风化破坏。

1 温度场与渗流场耦合基本方程

该文渗流场采用多孔介质布里克曼方程[19-24]，与温度场相结合。

(1)多孔介质模型传热的数学模型为：

(1)

(ρCp)eq=θpρpCρ·p+(1-θp)ρCρ

(2)

keq=θp+(1-θp)k

(3)

式中：ρ—液体密度(kg/m3)；Cp—液体等压热容(J/K)；k—热传导系数(W/(m2K))；u—速度场；Y—比热比，等压热容与等容热容之比；(ρCp)eq—等效体积等压热容；θp—材料的体积分数；kp—固体材料的热传导系数；ρp—固体材料密度(kg/m3)；Cρ·p—固体材料的比热容(J/(kg·K))。

(2)多孔介质流动的数学模型是由连续性方程和动量平衡方程组合而成，它们一起构成了Brinkman方程：

(4)

(5)

式中：μ—流体的动态粘滞系数(Pa·s)；u—流体的速度；ρ—流体的密度(kg/m3)；p—压力(Pa)；εp—孔隙率；k—渗透系数(m/d)；Qbr—质量源；F—体积力。

(3)温度场与渗流场的耦合：在温度场的计算中，速度采用布里克曼方程计算得出的速度；而在渗流场计算中，体积力计算公式为：

F=ρgβ(T-Tc)

(6)

式中：T—温度场计算得出的温度；Tc—初始温度。

2 温度场与渗流场耦合数值模型建立与条件设定

根据恐龙化石的简化模型和围岩的简化模型，建立热流耦合数值模型如图1所示[2]。其边界条件为：内部初始温度采用诸城地区多年平均气温12℃，左右两边界对称，前后两边界对称，上表面裸露于空气中。

图1 热流耦合数值模型 (据杜圣贤等，2015)

随着时间变化，外界温度变化过程如图2所示。其中包括温度上升和下降两个阶段。0～30000s为温度上升阶段，30000～60000s为温度下降阶段，之后温度又上升。

图2 时间-温度变化曲线

热流耦合模型参数如表1所示。

3 温度场与渗流场耦合数值模型结果分析

3.1 随外界温度变化化石及围岩表面渗流场变化

随着外界温度升高，分别模拟了0.5h(1800s)，1h(3600s)，2h(7200s)，3h(10800s)，6h(21600s)，8h(28800s)的变化情况(图3)。图中等值线为渗流场流速梯度的变化情况，矢量图为速度矢量的变化情况。可以看出，在外界温度变化的前8个小时(0～28800s)，随着外界温度不断升高，表面流场速度梯度迅速变化，速度矢量也不断变化。在0.5个小时(1800s)时，化石表面流场速度梯度最大，流速也是化石两端最大，而此时围岩表面流速梯度普遍变化较小。这是因为当外界温度变化时，由于化石裸露面积较大，化石内部孔隙水运动剧烈导致的；第1个小时(3600s)时，明显可以看出，在恐龙化石和围岩接触处，流场流速梯度和流速均明显大于其他部分，这是由于化石和围岩的各力学性质和水力学性质不同，导致它们之间存在相互作用，故在外界条件变化时，化石和围岩接触处受影响最大，此处风化速度会比其他地方快。

随着外界温度温度开始降低时，分别模拟了13h(46800s)，15h(54000s)，17h(61200s)，24h(86400s)的变化情况(图4)。可以看出，外界温度开始降低后，最先受影响的也是化石裸露区，然后扩展到围岩部分。15(54000s)～17h(61200s)，通过流速梯度和流速大小变化图形可以看出，外界温度降低的情况下，化石两端和内部速度梯度和速度大小非常不均匀，这种情况更会加速化石的风化。根据总体变化趋势可以发现，围岩部分内部空隙水随外界温度变化而变化持续的时间更长，故当外界条件，特别是水和温度持续变化时，围岩的风化速率可能会大于化石。总之，随外界温度的不断变化，化石和围岩的表面流场剧烈变化，孔隙水的运动会破坏化石和围岩的内部结构，导致其强度降低；水的作用会带动内部化学离子的运动，加速化学风化的进行。

表1 热流耦合模型参数

图3 温度上升时恐龙化石及围岩表面渗流场变化图

图4 温度降低时恐龙化石及围岩表面渗流场变化图

3.2 随外界温度变化化石及围岩内部渗流场变化

随着外界温度升高时，分别模拟了1h(3600s)，3h(10800s)的变化情况(图5)。可以看出，外界温度变化1h(3600s)时，内部流场在恐龙化石和围岩接触处变化最大，在围岩内部也有几处变化较大区域，但比例较小；随着外界温度的不断升高，到第3个小时(10800s)时，内部流场在化石和围岩中流速均较大，即内部空隙水运动剧烈。

随着外界温度开始降低直至再次升高时，分别模拟了12h(43200s)，18h(64800s)的变化情况(图6)。可以看出，随着外界温度的降低，第12个小时(43200s)时，恐龙化石及围岩内部流速都有所降低，但部分区域仍有较大的流速梯度；同时，在由温度升高到温度降低的过程中，其内部孔隙水的运动方向发生了变化。之后，随外界温度又回升，内外温差较大，孔隙水运动更加剧烈，到第18个小时(64800s)时，整个化石内部空隙水流动剧烈，围岩内孔隙水运动也开始深入内部。

3.3 随外界温度变化化石整体内部渗流场变化

随外界温度变化化石内部渗流场变化，分别模拟了9h(32400s)，24h(86400s)的变化情况(图7)。可以看出，随着外界温度变化，化石内部流速梯度较大，在化石裸露区两侧，空隙水向外界流出；而在化石与围岩接触区，孔隙水有穿过围岩区向上表面运动的趋势，而这种运动会加速化石与围岩接触区化石和围岩的风化；在化石中间部分也出现了较大的流速区，表明化石内部随着外界温度的变化，内部空隙水流动会出现集中区，循环往复，这部分化石的强度会降低，出现微小的裂缝，随后裂缝部分水的运动会更加集中，风化速度加快。

图5 温度升高时恐龙化石及围岩内部渗流场变化图

图6 温度降低至升高时恐龙化石及围岩内部渗流场变化图

图7 化石内部渗流场分布图

图8 化石及围岩热流耦合内部温度及流速变化图

3.4 化石及围岩温度场和渗流场耦合内部温度及流速变化

图8为恐龙化石及围岩热流耦合的二维切片，分别模拟了1h(3600s)，20h(72000s)的变化情况。可以看出，初始阶段温度场和渗流场相互作用，使得化石和围岩表面产生孔隙水运动；随着外界温度的持续变化，化石和围岩内部孔隙水运动范围加大，但可以看出，在化石和围岩接触处孔隙水运动速度远大于其他部分；随着外界温度的升降变化，化石和围岩内部空隙水运动的方向会发生变化；由于孔隙水的运动，影响了温度分布，20h(72000s)时的温度云图非常明显，这会更加加剧化石和围岩内部的温度差，使得内部产生温度应力，加速风化进程；同时温度的分布也影响了内部空隙水的运动速度。

4 结论

(1)从恐龙化石及围岩的内部结构出发，以恐龙化石及围岩作为岩体，介绍了恐龙化石风化过程中温度场与渗流场耦合作用的基本理论，并且建立了化石及围岩热流耦合作用的多孔介质模型，模拟了温度场与渗流场在恐龙化石及围岩内部及表面的相互作用，描述了恐龙化石在温度和水的相互作用下的破坏机制。

(2)实验结果显示，化石及围岩受日照影响产生的温度会影响其内部孔隙水的变化，影响渗流场的变化；同时，化石及围岩的孔隙水也会参与到化石及围岩内部的热量交换中，影响化石及围岩的温度场的分布。两者相互作用，共同对恐龙化石的破坏产生影响。

(3)通过对恐龙化石及围岩进行热流耦合作用的模拟研究，对恐龙化石的保护提供了一定的科学依据。今后可借助岩石力学、流体力学、传热学等相关学科，开展化石及围岩温度场—渗流场—应力场三场之间耦合作用研究，建立三场耦合作用的数学模型，为化石保护建立完整的物理场防护体系。