梵净山红云金顶形成演化流固耦合数值模拟分析

王健阳

(贵州大学 喀斯特环境与地质灾害防治重点实验室，贵阳 550025)

0 引 言

梵净山，武陵山脉主峰，得名于“梵天净土”，位于贵州省铜仁市的江口、印江、松桃交界处，最高山峰海拔2 572 m，地形最大高差2 100 m。梵净山总面积约960 km2，是第42届世界遗产大会认定的世界自然遗产，这里不仅有着大自然鬼斧神工造就的地质奇观，还有被誉为“基因库”、“人类的宝贵遗产”的丰富生物资源，同时这里是西南地区著名的千年佛教名山。梵净山红云金顶是武陵山脉的最高峰，因早晨朝阳照射常见红云瑞气环绕而闻名。随着贵州旅游业的迅速发展，越来越多的目光将聚焦于此。开展对梵净山红云金顶形成演化的研究，不仅对分析贵州高原地质环境变迁具有重要的科学意义，而且对景区发展规划产生深远影响。

流固耦合是指岩土体内部的渗流场和应力场之间相互影响、相互作用的过程[1]。国内有许多学者对梵净山地质开展了大量研究，主要是从构造地质学、环境地质学、水文地质学等传统方向入手，但结合数值模拟进行分析的却很少，考虑流固耦合影响的更是匮乏。因此，本文根据其地质历史背景，通过有限差分软件进行红云金顶形成演化的流固耦合分析，力求能最大限度地还原形成演化中渗流场和应力场相互作用带来的影响。研究结论有助于深化对梵净山地质构造发育规律的认识，有效保护好这来之不易的世界自然遗产。

1 地质背景

1.1 地形地貌

梵净山是云贵高原向湘西丘陵过渡斜坡上的第一高峰(相对高度达2 000 m)，它不仅是乌江与沅江分水岭，还是横亘于贵州、重庆、湖南、湖北4省(区)的武陵山脉的最高主峰。红云金顶位于梵净山顶部山脊，由一块高约104 m的锥柱状山体组成，整个山体高耸直立，上半段地形坡度达85°～90°，许多区段形成“凸嘴”、“凹岩腔”，下部两侧为脉状山体的陡坡地带，地形坡度达30°～55°。

1.2 地层岩性

研究区出露地层有第四系松散堆积层(Q4)与元古界板溪群乌叶组第一段(Pt3w1)、甲路组(Pt3j)及梵净山群(Pt2fj)。

第四系松散堆积层主要为残坡积碎石土、黏土夹碎石土等，广泛分布于下方斜坡和沟谷中，厚度0.2～3.5 m。元古界板溪群乌叶组第一段(Pt3w1)主要为浅灰及深灰色厚层变余细砂岩及变余硅质粉砂岩，夹浅灰色粉砂质板岩；浅灰、灰、灰绿色薄层、板状绢云母板岩及粉砂质板岩，中上部偶夹变余绢云母粉砂岩[2]。完整岩体较为坚硬，锤击声脆。表层岩体多风化严重，常具劈理构造，将岩体按一定方向分割成平行密集的薄片状、薄板状。岩体中节理裂隙发育。元古界板溪群甲路组(Pt3j)岩性为灰、紫灰色薄层含钙质绢云母板岩夹千枚状构造，该层厚度较薄，总体厚度约10余米，主要分布于乌叶组形成的陡崖之下的基座位置，岩体节理裂隙发育。梵净山群(Pt2fj)岩性为浅灰微带绿色厚层状变余砂岩、块状变质细碧-角斑岩及变质基性火山岩-蚀变岩等，间夹粉砂质板岩，为梵净山景区的主要地层。

1.3 地质构造

区内总体大地构造处于羌塘-扬子-华南板块扬子陆块江南复合造山带黔南坳陷区铜仁复式褶皱变形区[2]，区内构造复杂，主要发育北北东向构造及南北向扭曲性构造形迹。历史上经历了五陵、雪峰、燕山等多期构造运动，褶皱、断裂发育，区内主要构造有梵净山穹状背斜、郭栗压性断裂等。

1.4 地质演化

研究区处于扬子地层区，发育一套新元古代变质岩系，基底梵净山群为一套岩浆侵入碎屑岩发育的粉砂质板岩与变余砂岩系，与上部浅变质板溪群呈角度不整合。该区经历了以武陵构造旋回为特征的造山运动、以雪峰-加里东构造旋回为特征的内陆裂谷拉张运动、以海西-印支-燕山构造旋回为特征的地层褶曲断裂、以喜马拉雅及新构造旋回期为特征的整体抬升。该区古地理演化受大地构造演化控制。武陵构造旋回期、雪峰-加里东构造旋回期、海西-印支-燕山构造旋回期均为由浅及深再变浅的海侵-海退过程，构造旋回期末均因造山运动而成为陆地[3]。

武陵构造旋回期(新元古代中期)，从起始到结束呈现的岩浆活动性质，早期表现为夭折陆间裂谷(发育初始洋壳，红海型)的岛弧-弧后拉斑-钙碱性玄武质岩浆喷溢和侵位性质，晚-末期表现为后碰撞高温高分异A2型酸性岩浆侵位及后碰撞-板内S型酸性岩浆侵位。

雪峰-加里东构造旋回期(新元古代晚期-早古生代)，在武陵运动之后，黔东地区发育造山期后的磨拉石盆地相快速堆积的底砾岩，总体反映该时期属滨岸-浅海陆棚环境，该地区在沉积环境上已经摆脱了武陵造山运动的影响而基本夷平。随后扬子地块与华夏地块在此基础上再次发生裂解，其间南华裂谷盆地形成，由此形成板溪群、下江群和丹洲群反映出的台地、斜坡、盆地的沉积格局。处于该构造旋回期造山带边部的黔东地区表现为轻微变质。

海西-印支-燕山构造旋回期(晚古生代-早白垩世)，贵州已进入陆内裂陷演化发展阶段，沉积格局出现重大变化，黔东地区主要受江南复合造山带钦杭带晚古生代钦防海槽发展演化的远程影响。使该地区从泥盆-早白垩世，经历了从裂谷盆地至前陆盆地的演化历程。贵州处于滨太平洋与特提斯构造域的联合作用。发生于早白垩世末的燕山运动产生强烈而广泛的褶皱和断裂，使若干先期断裂复活。在强大的区域挤压应力作用下，受远程影响效应的黔东地区的老变质岩系遭区域低温动力变质叠加影响，产生与燕山期区域褶皱轴向一致的透入性劈理及定向组构，沿脆性断裂带的动力碎裂变质作用遍及全省。

喜马拉雅及新构造旋回期(晚白垩世-第四纪)，贵州处于板内隆升活动阶段，具有陆内活动性质。新构造活动主要表现为区域性隆升背景下的断块活动，具有明显的掀斜性、间歇性隆升和差异性隆升等特征，而且现代仍处在隆升趋势之中。对贵州现今发育的河谷阶地、第四系残坡积物、温泉、地震及地貌和水系格式的分布具有重要影响。

2 数值模拟流固耦合的基本理论

2.1 基本理论

对于流固耦合的基本理论[4-5]，前人已经作了大量研究说明，本文就不再赘述。

使用三维有限差分软件进行流固耦合计算时，流体在孔隙介质中的流动服从各向同性Darcy定律,孔隙介质可被看成是可变形体,同时满足Biot方程,主要包括运动方程、平衡方程、本构方程和几何方程,方程中通过建立孔隙压力p与流体渗流速度q及岩体应力σ、应变ε之间的关系实现流体与固体的耦合计算。在流固耦合计算中，有4种类型的渗流计算边界条件:①给定孔隙水压力；②给定边界外法线方向流速分量；③不透水边界；④透水边界。系统默认为不透水边界[6]。

2.1.1 运动方程

液体的流动规律通过Darcy定律来描述。对于恒定均匀流，该定律可表达如下：

(1)

2.1.2 平衡方程[7]

1) 流体质量平衡方程。公式如下：

(2)

式中：qi,j为流量，m3/s；qu为体积流体源强度，L/s；ζ为孔隙介质体内单位体积的液体体积改变量。

2) Cauchy运动方程。公式如下：

(3)

2.1.3 本构方程

体积应变的改变引起流体孔隙压力的变化；反过来，孔隙压力的变化也会导致体积应变的发生[8]。孔隙介质本构方程的增量形式为：

(4)

2.1.4 相容方程

应变率和速度梯度之间的关系为：

εij=(vi,j+vj,i)/2

(5)

2.2 建模步骤

建模步骤简化见图1，主线表示应力场计算的过程，左右两边表示考虑渗流场后的新增步骤。

图1 建模流程图

3 研究区形成演化过程模拟

3.1 三维有限差分软件建模

红云金顶出露的地层主要是板溪群，在武陵构造旋回期(新元古代中期)形成的是红云金鼎基底的梵净山群岩层，板溪群主要是雪峰-加里东构造旋回期(新元古代晚期-早古生代)沉积变质而来。在经历了海西-印支-燕山构造旋回期(晚古生代-早白垩世)和喜马拉雅及新构造旋回期(晚白垩世-第四纪)的一系列复杂的断裂、隆升、冰蚀等内外力地质作用下才形成了现在红云金顶的奇观。

因此，根据研究区的地质背景，建立三维模型进行数值分析。在建模时分为4层，分别代表研究区发生的4次构造旋回期，并且根据不同时期对研究区形成演化的影响，将每层的厚度作出相应的调整。在建立数值模型时考虑模型边界条件、卸荷等地质现象，模型的范围必须要足够大，以覆盖开挖可能影响的区域，同时又要兼顾计算效率。研究区的实际尺寸为27 m×25 m×104 m，参考已有的研究成果，为尽可能消除模型范围过小而导致边界效应的影响，四周边界距金顶水平距离不小于5倍金顶直径，模型的尺寸设计为444 m×400 m×342 m。为减少不必要的网格，网格划分原则为金顶附近密集，远处稀疏，模型共有节点126 511个，单元337 863个。

本模拟进行3次开挖来反映不同构造旋回期对研究区的影响，划分4个阶段(初始阶段→萌芽阶段→发展阶段→完成阶段)进行具体的形成演化描述，对比只受应力场作用和应力场与渗流场耦合作用下的开挖模拟结果，分析流固耦合条件下的变形破坏特征是否更符合现实。

3.2 参数设置

岩体采用摩尔-库仑弹塑性力学模型，渗流场模拟采用各向同性渗流模型，底部和两侧边界为不透水边界，顶部边界为自由边界。为了能体现海西-印支-燕山构造旋回期(晚古生代-早白垩世)的变质作用带来的岩体力学参数的变化，设置混合岩1和混合岩2进行区别。具体参数通过岩石物理力学试验的成果，结合一些规范、手册以及以往工程经验综合取值。具体参数见表1。

表1 岩体参数表

3.3 形成演化分析

3.3.1 初始阶段

经过武陵构造旋回期(新元古代中期)的岩浆活动后，在雪峰-加里东构造旋回期(新元古代晚期-早古生代)，相对深水滞留缺氧的环境中进行沉积，水平层理比较发育，该地区在沉积环境上已经摆脱了武陵造山运动的影响而基本夷平。所以，初始模型以此为地质背景，分为在应力场作用下的初始地应力场平衡(图2)和在力场与渗流场耦合作用下的初始孔压场平衡(图3)。

3.3.2 萌芽阶段(第一次开挖)

雪峰-加里东构造旋回期(新元古代晚期-早古生代)，基本夷平之后扬子地块与华夏地块在此基础上再次发生裂解，其间南华裂谷盆地形成，由此形成板溪群反映出的台地沉积格局，经历了从裂谷盆地-(初始)洋盆-前陆盆地-磨拉石盆地的演化历程。所以，原始地面模型以此为地质背景，这一时期以沉积为主，发育着局部构造变形迹象，在内力和外力的作用下地表不断剥蚀。第一次开挖后，在只考虑应力场的作用下，由图4、图5可知，原始山体形成，两侧出现卸荷变形区，中间出现应力集中区，山体未见明显塑性区。在考虑应力场与渗流场耦合作用下，由图6可知，山体表面出现大量塑性屈服和拉裂破坏单元，这是由于开挖后孔隙水压力的下降导致。

图2 模型初始应力分布云图

图3 模型初始孔压场分布云图

3.3.3 发展阶段(第二次开挖)

海西-印支-燕山构造旋回期(晚古生代-早白垩世)，以发育侏罗山式褶皱、日尔曼式褶皱、逆冲推覆断层、浅层滑脱构造、平行走滑断层为代表，构造运动方向由东向西，变形以脆性变形为主要特点，具浅部层次构造变形特征。区域褶皱伴随区域低温动力变质-脆性断裂伴随动力碎裂变质。所以，原始地面模型以此为地质背景，多次构造运动，既有激烈的褶皱运动，又有缓和的升降运动。以海相碎屑地层和火山岩经区域变质作用形成的地层为主。第二次开挖后，在只考虑应力场的作用下，由图7、图8可知，原始山体继续发展，不断的卸荷作用，在两侧出现卸荷变形区，中间出现应力集中区，山体有零星的塑性区。在考虑应力场与渗流场耦合作用下，由图9可知，山体周围局部出现较密集塑性屈服和拉裂破坏单元。

图4 第一次开挖X方向位移

图5 第一次开挖塑性区

图6 第一次开挖流固耦合塑性区

图7 第二次开挖X方向位移

图8 第二次开挖塑性区

图9 第二次开挖流固耦合塑性区

图10 第三次开挖X方向位移

3.3.4 完成阶段(第三次开挖)

喜马拉雅及新构造旋回期(晚白垩世-第四纪)，由海变陆并隆升成为高原，该构造阶段贵州受青藏高原隆升影响，垂直运动特征明显，典型构造样式为隆升背景下的地垒-地堑式构造组合，以脆性变形为主要特点，具浅表层次构造变形特征。所以，原始地面模型以此为地质背景，第三次开挖后，由图10、图11可知，在只考虑应力场的作用下，原始山体发育基本完成，经过多期的卸荷作用，风化作用的影响下，山体顶端因应力集中现象出现塑性区，其方向沿着构造面的走向随着时间不断发展，这也是红云金顶主裂缝的形成原因。在考虑应力场与渗流场耦合作用下，由图12可知，山体上部周围出现较密集塑性屈服和拉裂破坏单元。这也是红云金顶四周节理裂隙发育，不断剥落的原因。

图11 第三次开挖塑性区

图12 第三次开挖流固耦合塑性区

综上所述，通过分析可知考虑流固耦合作用更能体现这几次构造旋回期不断的造山运动所带来的影响，其结果也更符合现状。此外，梵净山地区第四纪历史上曾发生过多次冰川流行，从山上到山下保存着一套完整的冰蚀地貌、冰碛台地及冰碛泥砾等形迹[9]，对当今红云金顶地貌的形成也起到控制性作用。

4 结 论

1) 梵净山的大地质背景，经历了多期复杂的洋陆转换和板内活动，是一个处于充水环境下，不断沉积、变质、固结成岩，隆起成山的过程。用三维有限差分软件进行数值模拟分析，可以较好地反映红云金顶形成演化的每一个构造旋回期。

2) 采用流固耦合分析，考虑形成演化中应力场和渗流场之间的相互影响，对比只受应力场作用，模拟可知流固耦合的结果更符合当今地貌的发展特征。

3) 红云金顶目前处于整体稳定，但由于山体节理裂隙比较发育，存在发生崩塌的风险。作为世界自然遗产，应该进行针对性的保护措施治理。