梵凈山紅云金頂形成演化流固耦合數(shù)值模擬分析

梵凈山紅云金頂形成演化流固

耦合數(shù)值模擬分析［摘要］根據(jù)紅云金頂?shù)牡刭|(zhì)背景，采用三維有限差分軟件對(duì)其形成演化進(jìn)行兩組對(duì)比分析：一組只受應(yīng)力場(chǎng)作用，另一組受滲流場(chǎng)和應(yīng)力場(chǎng)耦合作用，通過開挖命令來模擬不同構(gòu)造時(shí)期對(duì)研究區(qū)形成演化的影響。結(jié)果表明，考慮巖土體內(nèi)部流固耦合作用能更好實(shí)現(xiàn)對(duì)紅云金頂形成演化較真實(shí)的還原，可為類似地貌的形成演化研究提供參考。［關(guān)鍵詞］形成演化；有限差分；流固耦合；梵凈山0引言梵凈山，武陵山脈主峰,得名于“梵天凈土”，位于貴州省銅仁市的江口、印江、松桃交界處，最高山峰海拔2572m,地形最大高差2100m。梵凈山總面積約960km2，是第42屆世界遺產(chǎn)大會(huì)認(rèn)定的世界自然遺產(chǎn)，這里不僅有著大自然鬼斧神工造就的地質(zhì)奇觀，還有被譽(yù)為“基因庫(kù)'、“人類的寶貴遺產(chǎn)”的豐富生物資源，同時(shí)這里是西南地區(qū)著名的千年佛教名山。梵凈山紅云金頂是武陵山脈的最高峰，因早晨朝陽(yáng)照射常見紅云瑞氣環(huán)繞而聞名。隨著貴州旅游業(yè)的迅速發(fā)展,越來越多的目光將聚焦于此。開展對(duì)梵凈山紅云金頂形成演化的研究，不僅對(duì)分析貴州高原地質(zhì)環(huán)境變遷具有重要的科學(xué)意義,而且對(duì)景區(qū)發(fā)展規(guī)劃產(chǎn)生深遠(yuǎn)影響。流固耦合是指巖土體內(nèi)部的滲流場(chǎng)和應(yīng)力場(chǎng)之間相互影響、相互作用的過程國(guó)內(nèi)有許多學(xué)者對(duì)梵凈山地質(zhì)開展了大量研究，主要是從構(gòu)造地質(zhì)學(xué)、環(huán)境地質(zhì)學(xué)、水文地質(zhì)學(xué)等傳統(tǒng)方向入手，但結(jié)合數(shù)值模擬進(jìn)行分析的卻很少，考慮流固耦合影響的更是匱乏。因此，本文根據(jù)其地質(zhì)歷史背景，通過有限差分軟件進(jìn)行紅云金頂形成演化的流固耦合分析,力求能最大限度地還原形成演化中滲流場(chǎng)和應(yīng)力場(chǎng)相互作用帶來的影響。研究結(jié)論有助于深化對(duì)梵凈山地質(zhì)構(gòu)造發(fā)育規(guī)律的認(rèn)識(shí)，有效保護(hù)好這來之不易的世界自然遺產(chǎn)。1地質(zhì)背景1.1地形地貌梵凈山是云貴高原向湘西丘陵過渡斜坡上的第一高峰（相對(duì)高度達(dá)2000m），它不僅是烏江與沅江分水嶺，還是橫亙于貴州、重慶、湖南、湖北4省（區(qū)）的武陵山脈的最高主峰。紅云金頂位于梵凈山頂部山脊，由一塊高約104m的錐柱狀山體組成，整個(gè)山體高聳直立，上半段地形坡度達(dá)85°-90°,許多區(qū)段形成“凸嘴'、&凹巖腔”，下部?jī)蓚?cè)為脈狀山體的陡坡地帶，地形坡度達(dá)30°~55°。1.2地層巖性研究區(qū)出露地層有第四系松散堆積層（Q4）與元古界板溪群烏葉組第一段（Pt3w1）、甲路組（Pt3j）及梵凈山群（Pt2fj）0第四系松散堆積層主要為殘坡積碎石土、黏土夾碎石土等，廣泛分布于下方斜坡和溝谷中，厚度0.2~3.5m。元古界板溪群烏葉組第一段（Pt3wl）主要為淺灰及深灰色厚層變余細(xì)砂巖及變余硅質(zhì)粉砂巖，夾淺灰色粉砂質(zhì)板巖；淺灰、灰、灰綠色薄層、板狀絹云母板巖及粉砂質(zhì)板巖,中上部偶夾變余絹云母粉砂巖⑵。完整巖體較為堅(jiān)硬，錘擊聲脆。表層巖體多風(fēng)化嚴(yán)重，常具劈理構(gòu)造，將巖體按一定方向分割成平行密集的薄片狀、薄板狀。巖體中節(jié)理裂隙發(fā)育。元古界板溪群甲路組（Pt3j）巖性為灰、紫灰色薄層含鈣質(zhì)絹云母板巖夾千枚狀構(gòu)造，該層厚度較薄，總體厚度約10余米，主要分布于烏葉組形成的陡崖之下的基座位置,巖體節(jié)理裂隙發(fā)育。梵凈山群（Pt2fj）巖性為淺灰微帶綠色厚層狀變余砂巖、塊狀變質(zhì)細(xì)碧-角斑巖及變質(zhì)基性火山巖-蝕變巖等，間夾粉砂質(zhì)板巖,為梵凈山景區(qū)的主要地層。1.3地質(zhì)構(gòu)造區(qū)內(nèi)總體大地構(gòu)造處于羌塘-揚(yáng)子-華南板塊揚(yáng)子陸塊江南復(fù)合造山帶黔南坳陷區(qū)銅仁復(fù)式褶皺變形區(qū)⑵，區(qū)內(nèi)構(gòu)造復(fù)雜，主要發(fā)育北北東向構(gòu)造及南北向扭曲性構(gòu)造形跡。歷史上經(jīng)歷了五陵、雪峰、燕山等多期構(gòu)造運(yùn)動(dòng)，褶皺、斷裂發(fā)育，區(qū)內(nèi)主要構(gòu)造有梵凈山穹狀背斜、郭栗壓性斷裂等。1.4地質(zhì)演化研究區(qū)處于揚(yáng)子地層區(qū)，發(fā)育一套新元古代變質(zhì)巖系，基底梵凈山群為一套巖漿侵入碎屑巖發(fā)育的粉砂質(zhì)板巖與變余砂巖系，與上部淺變質(zhì)板溪群呈角度不整合。該區(qū)經(jīng)歷了以武陵構(gòu)造旋回為特征的造山運(yùn)動(dòng)、以雪峰-加里東構(gòu)造旋回為特征的內(nèi)陸裂谷拉張運(yùn)動(dòng)、以海西-印支-燕山構(gòu)造旋回為特征的地層褶曲斷裂、以喜馬拉雅及新構(gòu)造旋回期為特征的整體抬升。該區(qū)古地理演化受大地構(gòu)造演化控制。武陵構(gòu)造旋回期、雪峰-加里東構(gòu)造旋回期、海西-印支-燕山構(gòu)造旋回期均為由淺及深再變淺的海侵-海退過程，構(gòu)造旋回期末均因造山運(yùn)動(dòng)而成為陸地⑶。武陵構(gòu)造旋回期（新元古代中期），從起始到結(jié)束呈現(xiàn)的巖漿活動(dòng)性質(zhì),早期表現(xiàn)為夭折陸間裂谷（發(fā)育初始洋殼，紅海型）的島弧-弧后拉斑-鈣堿性玄武質(zhì)巖漿噴溢和侵位性質(zhì)，晚-末期表現(xiàn)為后碰撞高溫高分異A2型酸性巖漿侵位及后碰撞-板內(nèi)S型酸性巖漿侵位。雪峰-加里東構(gòu)造旋回期（新元古代晚期-早古生代），在武陵運(yùn)動(dòng)之后，黔東地區(qū)發(fā)育造山期后的磨拉石盆地相快速堆積的底礫巖,總體反映該時(shí)期屬濱岸-淺海陸棚環(huán)境，該地區(qū)在沉積環(huán)境上已經(jīng)擺脫了武陵造山運(yùn)動(dòng)的影響而基本夷平。隨后揚(yáng)子地塊與華夏地塊在此基礎(chǔ)上再次發(fā)生裂解，其間南華裂谷盆地形成，由此形成板溪群、下江群和丹洲群反映出的臺(tái)地、斜坡、盆地的沉積格局。處于該構(gòu)造旋回期造山帶邊部的黔東地區(qū)表現(xiàn)為輕微變質(zhì)。海西-印支-燕山構(gòu)造旋回期（晚古生代-早白堊世），貴州已進(jìn)入陸內(nèi)裂陷演化發(fā)展階段，沉積格局出現(xiàn)重大變化，黔東地區(qū)主要受江南復(fù)合造山帶欽杭帶晚古生代欽防海槽發(fā)展演化的遠(yuǎn)程影響。使該地區(qū)從泥盆-早白堊世,經(jīng)歷了從裂谷盆地至前陸盆地的演化歷程。貴州處于濱太平洋與特提斯構(gòu)造域的聯(lián)合作用。發(fā)生于早白堊世末的燕山運(yùn)動(dòng)產(chǎn)生強(qiáng)烈而廣泛的褶皺和斷裂，使若干先期斷裂復(fù)活。在強(qiáng)大的區(qū)域擠壓應(yīng)力作用下，受遠(yuǎn)程影響效應(yīng)的黔東地區(qū)的老變質(zhì)巖系遭區(qū)域低溫動(dòng)力變質(zhì)疊加影響,產(chǎn)生與燕山期區(qū)域褶皺軸向一致的透入性劈理及定向組構(gòu)，沿脆性斷裂帶的動(dòng)力碎裂變質(zhì)作用遍及全省。喜馬拉雅及新構(gòu)造旋回期（晚白堊世-第四紀(jì)），貴州處于板內(nèi)隆升活動(dòng)階段,具有陸內(nèi)活動(dòng)性質(zhì)。新構(gòu)造活動(dòng)主要表現(xiàn)為區(qū)域性隆升背景下的斷塊活動(dòng)，具有明顯的掀斜性、間歇性隆升和差異性隆升等特征，而且現(xiàn)代仍處在隆升趨勢(shì)之中。對(duì)貴州現(xiàn)今發(fā)育的河谷階地、第四系殘坡積物、溫泉、地震及地貌和水系格式的分布具有重要影響。2數(shù)值模擬流固耦合的基本理論2.1基本理論對(duì)于流固耦合的基本理論［4-5］，前人已經(jīng)作了大量研究說明，本文就不再贅述。使用三維有限差分軟件進(jìn)行流固耦合計(jì)算時(shí)，流體在孔隙介質(zhì)中的流動(dòng)服從各向同性Darcy定律，孔隙介質(zhì)可被看成是可變形體，同時(shí)滿

足Biot方程，主要包括運(yùn)動(dòng)方程、平衡方程、本構(gòu)方程和幾何方程，方程中通過建立孔隙壓力p與流體滲流速度6及巖體應(yīng)力"、應(yīng)變+之間的關(guān)系實(shí)現(xiàn)流體與固體的耦合計(jì)算。在流固耦合計(jì)算中，有4種類型的滲流計(jì)算邊界條件:①給定孔隙水壓力；②給定邊界外法線方向流速分量；③不透水邊界;④透水邊界。系統(tǒng)默認(rèn)為不透水邊界⑹。1.1運(yùn)動(dòng)方程液體的流動(dòng)規(guī)律通過Darcy定律來描述。對(duì)于恒定均勻流，該定律可表達(dá)如下：6,=-,；!,（s）[p-pf-jgjh （1）式中:6；為單位流量向量;p為孔隙水壓力；，初為滲透系數(shù)張量;（（s）為與飽和度s相關(guān)的滲透系數(shù);pf為流體密度;g.（；=1，3）為重力加速度的3個(gè)分量,m/s2。2.1.2平衡方程「刃1）流體質(zhì)量平衡方程。公式如下：%）6，+67=彼(2)式中:6頃為流量，m@/s;g.為體積流體源強(qiáng)度，(2)）為孔隙介質(zhì)體內(nèi)單位體積的液體體積改變量。2） Cauchy運(yùn)動(dòng)方程。公式如下：dv-皿+"7="式中:g為重力加速度;%為速度對(duì)時(shí)間的導(dǎo)數(shù);"方為應(yīng)力張量。本構(gòu)方程體積應(yīng)變的改變引起流體孔隙壓力的變化;反過來,孔隙壓力的變化也會(huì)導(dǎo)致體積應(yīng)變的發(fā)生囪。孔隙介質(zhì)本構(gòu)方程的增量形式為：+a-p.=）（,-+礦） （4）式中:為應(yīng)力增量;）為給定函數(shù);+初為總應(yīng)變。相容方程應(yīng)變率和速度梯度之間的關(guān)系為：=（I,+1，）/22.2建模步驟建模步驟簡(jiǎn)化見圖1,主線表示應(yīng)力場(chǎng)計(jì)算的過程，左右兩邊表示考慮滲流場(chǎng)后的新增步驟。圖1建模流程圖3研究區(qū)形成演化過程模擬1三維有限差分軟件建模紅云金頂出露的地層主要是板溪群，在武陵構(gòu)造旋回期（新元古代中期）形成的是紅云金鼎基底的梵凈山群巖層，板溪群主要是雪峰-加里東構(gòu)造旋回期（新元古代晚期-早古生代）沉積變質(zhì)而來。在經(jīng)歷了海西-印支-燕山構(gòu)造旋回期（晚古生代-早白堊世）和喜馬拉雅及新構(gòu)造旋回期（晚白堊世-第四紀(jì)）的一系列復(fù)雜的斷裂、隆升、冰蝕等內(nèi)外力地質(zhì)作用下才形成了現(xiàn)在紅云金頂?shù)钠嬗^。因此，根據(jù)研究區(qū)的地質(zhì)背景，建立三維模型進(jìn)行數(shù)值分析。在建模時(shí)分為4層，分別代表研究區(qū)發(fā)生的4次構(gòu)造旋回期，并且根據(jù)不同時(shí)期對(duì)研究區(qū)形成演化的影響，將每層的厚度作出相應(yīng)的調(diào)整。在建立數(shù)值模型時(shí)考慮模型邊界條件、卸荷等地質(zhì)現(xiàn)象，模型的范圍必須要足夠

大，以覆蓋開挖可能影響的區(qū)域，同時(shí)又要兼顧計(jì)算效率。研究區(qū)的實(shí)際尺寸為27mx25mx104m,參考已有的研究成果，為盡可能消除模型范圍過小而導(dǎo)致邊界效應(yīng)的影響，四周邊界距金頂水平距離不小于5倍金頂直徑,模型的尺寸設(shè)計(jì)為444mx400mx342m。為減少不必要的網(wǎng)格，網(wǎng)格劃分原則為金頂附近密集，遠(yuǎn)處稀疏，模型共有節(jié)點(diǎn)126511個(gè)，單元337863個(gè)。本模擬進(jìn)行3次開挖來反映不同構(gòu)造旋回期對(duì)研究區(qū)的影響，劃分4個(gè)階段（初始階段#萌芽階段#發(fā)展階段#完成階段）進(jìn)行具體的形成演化描述,對(duì)比只受應(yīng)力場(chǎng)作用和應(yīng)力場(chǎng)與滲表1巖體名稱體積模量/MPa剪切模量巖體名稱體積模量/MPa剪切模量/MPa抗拉強(qiáng)度/kPa黏聚力/MPa內(nèi)摩擦角/（"）密度kg-(m3)_：滲透系數(shù)1/m2,(Pa-sec)_1孔隙率混合巖1103821000.532.520003.06X10110.35混合巖23031823001.242.526703.06X10-100.25巖體參數(shù)表3.3.1初始階段經(jīng)過武陵構(gòu)造旋回期（新元古代中期）的巖漿活動(dòng)后,在雪峰-加里東構(gòu)造旋回期（新元古代晚期-早古生代），相對(duì)深水滯留缺氧的環(huán)境中進(jìn)行沉積，水平層理比較發(fā)育，該地區(qū)在沉積環(huán)境上已經(jīng)擺脫了武陵造山運(yùn)動(dòng)的影響而基本夷平。所以，初始模型以此為地質(zhì)背景，分為在應(yīng)力場(chǎng)作用下的初始地應(yīng)力場(chǎng)平衡（圖2）和在力場(chǎng)與滲流場(chǎng)耦合作用下的初始孔壓場(chǎng)平衡（圖3）03.3.2萌芽階段（第一次開挖）雪峰-加里東構(gòu)造旋回期（新元古代晚期-早古生代），基本夷平之后揚(yáng)子地塊與華夏地塊在此基礎(chǔ)上再次發(fā)生裂解,其間南華裂谷盆地形成，由此形成板溪群反映出的臺(tái)地沉積格局，經(jīng)歷了從裂谷盆地？（初始）洋盆-前陸盆地-磨拉石盆地的演化歷程。所以，原始地面模型以此為地質(zhì)背景,這一時(shí)期以沉積為主，發(fā)育著局部構(gòu)造變形跡象，在內(nèi)力和外力的作用下地表不斷剝蝕。第一次開挖后，在只考慮應(yīng)力場(chǎng)的作用下，由圖4、圖5可知,原始山體形成，兩側(cè)出現(xiàn)卸荷變形區(qū)，中間出現(xiàn)應(yīng)力集中區(qū)，山體未見明顯塑性區(qū)。在考慮應(yīng)力場(chǎng)與滲流場(chǎng)耦合作用下，由圖6可知，山體表面出現(xiàn)大量塑性屈服和拉裂破壞單元，這是由于開挖后孔隙水壓力的下降流場(chǎng)耦合作用下的開挖模擬結(jié)果,分析流固耦合條件下的變形破壞特征是否更符合現(xiàn)實(shí)。3.2參數(shù)設(shè)置巖體采用摩爾-庫(kù)侖彈塑性力學(xué)模型，滲流場(chǎng)模擬采用各向同性滲流模型，底部和兩側(cè)邊界為不透水邊界，頂部邊界為自由邊界。為了能體現(xiàn)海西-印支-燕山構(gòu)造旋回期（晚古生代-早白堊世）的變質(zhì)作用帶來的巖體力學(xué)參數(shù)的變化，設(shè)置混合巖1和混合巖2進(jìn)行區(qū)別。具體參數(shù)通過巖石物理力學(xué)試驗(yàn)的成果，結(jié)合一些規(guī)范、手冊(cè)以及以往工程經(jīng)驗(yàn)綜合取值。具體參數(shù)見表1。圖2模型初始應(yīng)力分布云圖ContourOfGpPorePressure7.5522E+06I7.5000E+06I7.0000E+066.5000E+066.0000E+065.5000E+065.0000E+064.5000E+064.0000E+063.5000E+063.0000E+062.5000E+062.0000E+061.5000E+06U1.0000E+06I5.0000E+05|0.0000E+00圖3模型初始孔壓場(chǎng)分布云圖3.3.3 發(fā)展階段（第二次開挖）海西-印支-燕山構(gòu)造旋回期（晚古生代-早白堊世），以發(fā)育侏羅山式褶皺、日爾曼式褶皺、逆沖推覆斷層、淺層滑脫構(gòu)造、平行走滑斷層為代表，構(gòu)造運(yùn)動(dòng)方向由東向西,變形以脆性變形為主要特點(diǎn)，具淺部層次構(gòu)造變形特征。區(qū)域褶皺伴隨區(qū)域低溫動(dòng)力變質(zhì)-脆性斷裂伴隨動(dòng)力碎裂變質(zhì)。所以，原始地面模型以此為地質(zhì)背景,多次構(gòu)造運(yùn)動(dòng)，既有激烈的褶皺運(yùn)動(dòng)，又有緩和的升降運(yùn)動(dòng)。以海相碎屑地層和火山巖經(jīng)區(qū)域變質(zhì)作用形成的地層為主。第二次開挖后，在只考慮應(yīng)力場(chǎng)的作用下，由圖7、圖8可知，原始山體繼續(xù)發(fā)展，不斷的卸荷作用，在兩側(cè)出現(xiàn)卸荷變形區(qū)，中間出現(xiàn)應(yīng)力集中區(qū)，山體有零星的塑性區(qū)。在考慮應(yīng)力場(chǎng)與滲流場(chǎng)耦合作用下，由圖9可知，山體周圍局部出現(xiàn)較密集塑性屈服和拉裂破壞單元。ContourOfX-DisplacementContourOfX-Displacement圖4第一次開挖X方向位移圖5第一次開挖塑性區(qū)圖6第一次開挖流固耦合塑性區(qū)3.3.4完成階段（第三次開挖）喜馬拉雅及新構(gòu)造旋回期（晚白堊世-第四紀(jì)），由海變陸并隆升成為高原，該構(gòu)造階段貴州受青藏高原隆升影響，垂直運(yùn)動(dòng)特征明顯，典型構(gòu)造樣式為隆升背景下的地壘-地塹式構(gòu)造組ContourOfX-DisplacementContourOfX-Displacement5.5685E-02I5.0000E-024.0000E-023.0000E-022.0000E-021.0000E-020.0000E+00-1.0000E-02-2.0000E-02-3.0000E-02I-4.0000E-02|-4.7787E-02圖7第二次開挖X方向位移圖8第二次開挖塑性區(qū)圖9第二次開挖流固耦合塑性區(qū)圖10