趙延林 ，曹平，萬文 ，王衛(wèi)軍 ，張盛國 ，鄒聲華

(1. 湖南科技大學(xué) 煤礦安全開采技術(shù)湖南省重點實驗室，湖南 湘潭，411201；2. 湖南科技大學(xué) 能源與安全工程學(xué)院，湖南 湘潭，411201；3. 中南大學(xué) 資源與安全工程學(xué)院，湖南 長沙，410083)

突水事故造成的直接經(jīng)濟損失一直排在各類煤礦災(zāi)害前列。南方煤礦突水事故的水源主要有巖溶水、老窯積水和地表水，而且三者之間連通關(guān)系復(fù)雜，往往其中一類突水導(dǎo)致其他突水事故并發(fā)[1-3]。在各類水源中, 地下巖溶突水壓力大, 突發(fā)性特別強, 造成的破壞損失最嚴重。礦井開挖過程中，當擊穿大型充水巖溶體(如暗河、填充溶洞及巖溶破碎帶等)時，會引發(fā)礦井涌水突水地質(zhì)災(zāi)害[4-5]。如2003 年湖南漣源市七一煤礦石壩井的“4·16”特大突水災(zāi)害(17 人遇難)、2010 年湖南張家界市鑫茂煤礦“9.28”溶洞透水事故、2004 年貴州思南縣天池煤礦“12.12”特大透水事故(21人死亡，15 人失蹤)都是晚二疊系茅口灰?guī)r內(nèi)高承壓溶洞水體突破巖層引起的突水事故。對于隱伏溶洞突水機理的研究主要集中在西南巖溶隧道突水機理的研究上。劉招偉等[6-8]分別對圓梁山隧道、滬蓉西高風(fēng)險巖溶隧道、季家坡隧道的溶洞突水致災(zāi)機制和圍巖、襯砌破壞機理進行了研究，總結(jié)了巖溶突水的不同地質(zhì)模式、基本力學(xué)機制、突水流態(tài)演化模型和襯砌壓裂致災(zāi)機制。目前，巖溶突水的力學(xué)機理研究主要針對巖溶隧道，而對于巖溶礦井在采掘過程中由于溶洞承壓水擊穿巖柱而產(chǎn)生工作面巖溶突水的力學(xué)機理和突水過程分析研究較少。在巖溶礦井中，一般通過布置探放水鉆孔來探放巖溶水，通過防突巖柱的留設(shè)來防治巖溶突水。而《煤礦防治水規(guī)定》[9]沒有具體給出針對隱伏承壓溶洞防突巖柱留設(shè)的安全厚度計算方法。進行隱伏承壓溶洞防突巖柱的穩(wěn)定性及安全厚度研究，對于防治礦井巖溶突水具有重要意義[10-11]。本文作者從巖石力學(xué)和滲流力學(xué)出發(fā)采用流固耦合分析方法探討巷道前伏承壓溶洞防突巖柱在采動應(yīng)力和高水力梯度作用下的流固耦合效應(yīng)，研究高水壓下防突巖柱的失穩(wěn)演化，同時引入強度折減法理論提出防突巖柱安全系數(shù)的概念，分析巷道前伏承壓溶洞防突巖柱的安全厚度，為合理設(shè)計和預(yù)留防突巖柱提供理論依據(jù)。

1 承壓溶洞突水的流固耦合機理

當巷道前伏承壓溶洞時，采動應(yīng)力誘發(fā)圍巖應(yīng)力場與滲流場的急劇變化，在溶洞與工作面之間的圍巖部位將會顯現(xiàn)應(yīng)力場和滲流場的強耦合效應(yīng)。這種耦合作用表現(xiàn)在：高水力梯度(高動水壓力)作用于圍巖導(dǎo)致應(yīng)力場的改變；在高靜水壓力和應(yīng)力場共同作用下圍巖裂隙擴展，其滲透特性大大加強。從巖體流體力學(xué)角度出發(fā)，承壓溶洞突水可認為是采動應(yīng)力和高動、靜水壓力作用下，溶洞與工作面之間巖柱的突變失穩(wěn)并伴隨巖溶水大規(guī)模涌出。防突巖柱在采動應(yīng)力與滲流的雙重作用下強度不斷降低，當系統(tǒng)狀態(tài)瀕臨臨界狀態(tài)時，微小的擾動誘發(fā)防突層整體破斷突水。

1.1 非線性滲流分析

防突巖柱滲流分析微分控制方程為[12]：

式中：p 為滲透壓；k 為滲透系數(shù)；M 為比奧模量；ε為體積應(yīng)變；α 為有效應(yīng)力的比奧系數(shù)。

滲透壓和應(yīng)力對巖體滲透性的影響如下[13]

式中：k0為巖體滲透系數(shù)初值； k(Θ , p)為耦合分析下巖體的滲透系數(shù)；Θ = σ1+ σ2+ σ3，為體積應(yīng)力；σ1，σ2， σ3分別為第一、第二、第三有效主應(yīng)力；ξ 為滲透突跳系數(shù)。

在耦合分析中，對于處于彈性狀態(tài)的單元，認為巖體滲透系數(shù)是有效體積應(yīng)力的負指數(shù)函數(shù)，取ξ=1.0，β=0.5。

對于塑性屈服單元，滲透系數(shù)大大增加，滲透突跳系數(shù)顯著提高，取ξ=5 000，β =1.0。

式(2)反映了巖體應(yīng)力對滲流的耦合作用，尤其反映了塑性單元對滲流的強耦合效應(yīng)。

1.2 彈塑性力學(xué)分析

應(yīng)力場與滲流場耦合作用下巖體的彈塑性力學(xué)本構(gòu)關(guān)系[13-14]：

式中：(α p),j為滲透水力梯度作為等效體積力作用于巖體骨架，反映了滲流場對應(yīng)力場的耦合作用。

屈服準則采用Mohr-Coulomb 準則。該準則的剪切破壞判據(jù)為：

拉伸破壞判據(jù)為：

式中： σt為巖體抗拉強度；φ 為巖體內(nèi)摩擦角；c 為黏聚力。

1.3 流固耦合程序設(shè)計

巷道承壓溶洞突水的流固耦合分析，采用間接耦合手段，首先計算出ti時刻巖柱的彈塑性應(yīng)力場，由式(2)計算出ti時刻巖柱滲透系數(shù)，將該滲透系數(shù)傳遞給滲流場計算模塊，實現(xiàn)應(yīng)力場對滲流的耦合；再將滲流場計算得到的滲透體積力施加給彈塑性應(yīng)力計算單元，完成滲流場對應(yīng)力場的耦合。如此對時間循環(huán)，直至達到預(yù)定的滲流計算時刻。為研究防突巖柱的流固耦合效應(yīng)，在FLAC3D2.1 的平臺下，通過FISH 語言成功開發(fā)了巷道前伏承壓溶洞防突巖柱流固耦合分析程序。該程序分為彈塑性應(yīng)力計算、滲流計算和耦合分析3 個模塊。

基于FLAC3D開發(fā)的流固耦合間接分析方法，其優(yōu)點體現(xiàn)在各耦合參數(shù)(如滲透系數(shù)、滲透體積力等)隨計算時步時空演化和各耦合子系統(tǒng)(滲流子系統(tǒng)和彈塑性力學(xué)子系統(tǒng))的耦合信息隨計算時步準確傳遞。

1.4 流固耦合-強度折減法聯(lián)動分析

采用強度折減法思想研究防突巖柱的穩(wěn)定性。

強度折減法的基本原理是將巖體的抗剪參數(shù)c,tanφ 同時除以某個折減系數(shù)Fi, 得到一組新的參數(shù)c′，φ′, 將新的材料參數(shù)作為計算參數(shù)進行試算。

將承壓溶洞突水的流固耦合理論和強度折減法結(jié)合起來，形成流固耦合-強度折減法聯(lián)動分析方法。研究在防突巖柱強度參數(shù)不斷折減的情況下，巖柱失穩(wěn)突出和滲流突變的演化過程。

流固耦合-強度折減法聯(lián)動分析中，首先進行彈塑性應(yīng)力-滲流的流固耦合分析，判斷計算是否收斂，若收斂則對防突巖柱的強度參數(shù)進行折減，并記錄折減系數(shù)，用折減后強度參數(shù)作為防突巖柱的強度，進行流固耦合分析，再判斷計算是否收斂，如此重復(fù)，直至計算不收斂，此時的折減系數(shù)即為防突巖柱的安全系數(shù)。流固耦合-強度折減法聯(lián)動分析流程圖如圖1所示。

圖1 流固耦合-強度折減法聯(lián)動分析流程圖Fig.1 Flow diagram of Linkage analysis between fluid-solid coupling and strength reduction method

2 突水案例

以湖南省婁底漣源市七一煤礦石壩井-160 m 水平老水倉擴容掘進工作面發(fā)生巖溶突水事故為案例，進行防突巖柱的穩(wěn)定性分析。該礦井水文地質(zhì)條件復(fù)雜，煤層下伏茅口灰?guī)r層，在茅口灰?guī)r中溶洞往往成串珠狀存在，溶洞呈半充填狀態(tài)，充填物含泥砂。23采區(qū)-160 m 水平老水倉擴容掘進工作面位于茅口灰?guī)r之中，2003 年4 月15 日在掘進過程中發(fā)生了承壓溶洞突水的重大事故，16 人遇難，28 h 共計涌水約12 500 m3，涌泥5 440 m3。清碴后，在工作面前方3～4 m 處露出1 個10 m×6 m×10 m 左右的溶洞。

2.1 計算模型

圖2 所示為七一煤礦石壩井-160 m水平老水倉擴容掘進防突巖柱穩(wěn)定性分析的計算模型，模型尺寸x×y×z=50 m×15 m×30 m，掘進巷道位于模型中部，半圓拱形巷道，模型右邊界設(shè)有承壓溶洞，為簡化假定溶洞為近橢球形，老水倉擴容巷道埋深450 m。

力學(xué)邊界及網(wǎng)格劃分見圖2，其中，y=0 平面為固定約束，σ1=12.5 MPa，σ3=10 MPa。滲流邊界條件：茅口灰?guī)r內(nèi)溶洞水壓p 分別取2，3 和4 MPa，工作面水壓為0.1 MPa。

圖2 中，1～5 為模擬分步掘進。在數(shù)值計算中，假定每步進尺5 m。當工作面距溶洞只有10 m 時，為研究不同厚度防突巖柱的穩(wěn)定性，假定每步進尺1 m，以防突巖柱厚度分別為10，9，8，7，6，5 和4 m 為例研究不同厚度防突巖柱的穩(wěn)定性。

對石壩井-160 m 水平的茅口灰?guī)r層進行取樣，并進行巖石力學(xué)和滲透性實驗，得到該類巖石的強度參數(shù)和滲流力學(xué)參數(shù)見表1。

采用流固耦合-強度折減法聯(lián)動分析法對承壓溶洞突水過程展開研究。以溶洞內(nèi)壓p=4 MPa，防突巖柱h=8 m 為例來研究不同強度下防突巖柱失穩(wěn)演化的突水突變過程。

圖2 計算模型Fig.2 Numerical calculation model

2.2 聯(lián)動分析結(jié)果

隨著折減系數(shù)的增大，防突巖柱強度參數(shù)削減，在采動應(yīng)力和高水壓力作用下，防突巖柱的塑性區(qū)擴展演化，導(dǎo)致巖柱的滲透系數(shù)顯著增加。圖3 所示為折減系數(shù)Fi=1.10 和1.54 時防突巖柱的滲透系數(shù)分布圖。Fi=1.10 時，巖柱強度參數(shù)為c=1.36 MPa，tan φ=0.52，巖柱滲透系數(shù)主要受體積應(yīng)力和孔隙壓的影響，巖柱滲透系數(shù)自工作面向溶洞形成了工作面附近的滲透增強區(qū)、巖柱中部的未擾動滲透區(qū)及溶洞附近的滲透增強區(qū)。距工作面1.5 m 左右的滲透增強區(qū)，滲透系數(shù)為3.25×10-10～15.21×10-10m/s；巖柱中部的未擾動滲透區(qū)，其滲透系數(shù)為0.70×10-10～1.79×10-10m/s；溶洞附近的滲透增強區(qū)，滲透系數(shù)為2.97×10-10～5.68×10-10m/s。Fi=1.54 時，巖柱強度參數(shù)為c=0.97 MPa，tan φ=0.39，此強度下巖柱內(nèi)85%左右的單元處于塑性狀態(tài)，防突巖柱滲透系數(shù)成數(shù)量級增加，為1.47×10-6～2.80×10-6m/s。在此強度下的巖柱不具備防突水的能力。

圖3 防突巖柱的滲透系數(shù)分布Fig.3 Distribution of hydaulic conductivity of rock pillar for preventing water inrush

表1 巖體力學(xué)計算參數(shù)Table 1 Mechanical parameters of rock mass

圖4 所示為折減系數(shù) Fi=1.10 時，計算模型的滲流場分布。巷道隱伏承壓溶洞，掘進擾動了圍巖滲流場，巖溶水在高水力梯度作用下向臨空面排出，滲流場不斷調(diào)整以消除承壓溶洞高勢能不穩(wěn)定狀態(tài)。在巖柱內(nèi)形成了自溶洞向工作面逐漸降低的滲透壓，聚集于巖柱的平均滲透壓梯度高達0.5 MPa/m，高滲透壓梯度產(chǎn)生指向工作面的水平體積力是導(dǎo)致防突巖柱突出失穩(wěn)的動力。

圖4 折減系數(shù)Fi =1.10 時滲透場分布Fig.4 Distribution of seepage field when Fi =1.10

圖5 不同折減系數(shù)下的工作面的滲流量Fig.5 Discharge of working face under different reduction factors

在流固耦合分析中，通過記錄工作面的流量，得到不同折減系數(shù)下工作面的滲流量見圖5。由圖5 可知：伴隨著防突巖柱強度的減弱，巖柱工作面滲流量也由小流量演化成大規(guī)模的涌水量。Fi=1.10 時，工作面流量為6.158×10-4m3/s，此時該強度下的巖柱穩(wěn)定，巖柱內(nèi)塑性區(qū)范圍不大，工作涌水流量少；Fi=1.54～1.58 時，巖柱內(nèi)塑性單元占85%～87%，與塑性單元強耦合的巖柱滲透系數(shù)突跳增長，導(dǎo)致工作涌水量急劇增加，工作面涌水量高達0.753～ 0.811 m3/s，當Fi=1.59 時，防突巖柱全部處于塑性狀態(tài)，計算不收斂，工程上則表現(xiàn)為承壓溶洞突水。

在采動應(yīng)力和滲透體積力共同作用下巖柱向工作面產(chǎn)生位移，隨著折減系數(shù)的增加，其位移量也隨之增大。圖6 所示為當Fi=1.54 時，防突巖柱水平位移等值線圖，耦合計算得出：當Fi=1.54 時，防突巖柱向工作面位移高達14.58 cm。而當Fi=1.59 時，力學(xué)計算不收斂，此時防突巖柱向工作面突出失穩(wěn)。

圖6 折減系數(shù)Fi =1.54 時防突巖柱中剖面上水平位移云圖Fig.6 Horizontal displacement diagrams of rock pillar when Fi =1.54

2.3 防突巖柱安全系數(shù)

采用流固耦合-強度折減法聯(lián)動分析，改變?nèi)芏磧?nèi)壓p 和巖柱厚度h，得到水壓p=2，3 和4 MPa 時，防突巖柱計算厚度h 分別為3，4，5，6，7，8，9 和10 m 時，防突巖柱的安全系數(shù)和溶洞水壓、防突巖柱厚度的關(guān)系曲線(圖7)。

由圖7 可知：在一定溶洞水壓下，隨著防突巖柱厚度的增大，巖柱安全系數(shù)也逐漸增大，采用指數(shù)函數(shù)對其進行擬合：

式中：Fs為防突巖柱安全系數(shù)；h 為防突巖柱計算厚度；a，b，d 為待定系數(shù)。

由式(9)可擬合得到不同溶洞水壓下防突巖柱的安全系數(shù)與巖柱計算厚度之間關(guān)系式，見表2。將Fs=1.5 作為評價防突巖柱穩(wěn)定性的判斷標準，可得出如下判斷標準[15]：

(1) Fs≥1.5 時防突巖柱穩(wěn)定；

(2) Fs＜1.5 時防突巖柱不穩(wěn)定。

圖7 防突巖柱安全系數(shù)和巖柱計算厚度的關(guān)系曲線Fig.7 Variation relation curves between safe factor of rock pillar for preventing water inrush and rock pillar thickness

表2 不同溶洞水壓下防突巖柱安全系數(shù)與巖柱計算厚度的關(guān)系Table 2 Relationship between safety factor of rock pillar for preventing water inrush and rock pillar thickness under different karst cave water pressures

由表2 可見：隨溶洞水壓力增大，安全防突巖柱計算厚度要相應(yīng)增加，如溶洞水壓為2 MPa 和4 MPa時，安全防突巖柱計算厚度分別為5.8 m 和8.0 m，后者的安全防突巖柱厚度要比前者多2.2 m。這是由于溶洞水壓的增加，聚集在防突巖柱內(nèi)的滲透體積力增大，巖柱受到向工作面方向的側(cè)向推力增加，因此安全防突巖柱要相應(yīng)增加留設(shè)厚度。

在七一煤礦石壩井-160 m 水平老水倉擴容掘進工作面發(fā)生的溶洞高承壓水體突破巖柱涌出的重大突水事故中發(fā)現(xiàn)：10 m×6 m×10 m 的溶洞距工作面3～4 m，即工作面與溶洞之間只有3～4 m 的防突巖柱，且?guī)r溶水與地表水溝通良好，溶洞內(nèi)水壓為3.5～4.0 MPa。防突巖柱工程留設(shè)厚度至少應(yīng)該達到炮眼長度(約2.2 m)和安全防突巖柱計算厚度(7～8 m)之和，為9.2～10.2 m。3～4 m 厚的防突巖柱完全不具備有防突功能。

3 結(jié)論

(1) 巷道隱伏承壓溶洞突水的流固耦合模型揭示突水機理在于采動應(yīng)力和高滲流體積力共同作用于防突巖柱，導(dǎo)致巖柱塑性區(qū)擴展而突出失穩(wěn)，伴隨巖柱的力學(xué)失穩(wěn)，防突巖柱發(fā)生滲流突變，工作面發(fā)生大規(guī)模突水。

(2) 采用承壓溶洞突水的流固耦合-強度折減法聯(lián)動分析方法，研究高水力梯度作用下防突巖柱的安全儲備，提出防突巖柱安全系數(shù)的概念，為安全防突巖柱留設(shè)厚度提供了一種新的定量方法。

(3) 將具有安全系數(shù)為1.5 的防突巖柱厚度作為防突巖柱的計算安全厚度。隨溶洞水壓的增加，防突巖柱的計算安全厚度相應(yīng)增大。

(4) 七一煤礦石壩井-160 m 水平老水倉擴容掘進工作面發(fā)生的溶洞突水是由于防突巖柱留設(shè)過小，不具有安全儲備，導(dǎo)致溶洞高承壓水體突破巖柱涌出。

[1] 李敬敬, 常聚才, 陳陽. 厚松散含水層下特厚煤層綜放開采相似模擬[J]. 礦業(yè)工程研究, 2013, 28(4): 47-50.LI Jingjing, CHANG Jucai, CHEN Yang. Similar simulation research on caving mining of extra thick coal seam under thick loose aquifer[J]. Mineral Engineering Research, 2013, 28(4):47-50.

[2] 李文敏. 承壓水上采煤的流固耦合數(shù)值模擬[J]. 礦業(yè)工程研究, 2013, 28(2): 37-41.LI Wenmin. The fluid-solid coupling numerical simulation of confined water mining[J]. Mineral Engineering Research, 2013,28(2): 37-41.

[3] 李建明, 朱斌, 武強. 趙各莊礦大傾角煤層綜放開采突水條件的力學(xué)機制分析[J]. 煤炭學(xué)報, 2007, 32(5): 453-457.LI Jianming, ZHU Bin, WU Qiang. Mechanical mechanism analysis of water inrushcondition on mechanized top-coal caving in steep seam in Zhaogezhuangcoal mine[J]. Journal of China Coal Society, 2007, 32(5): 453-457.

[4] 楊天鴻, 唐春安, 譚志宏, 等. 巖體破壞突水模型研究現(xiàn)狀及突水預(yù)測預(yù)報研究發(fā)展趨勢[J]. 巖石力學(xué)與工程學(xué)報, 2007,26(2): 268-277.YANG Tianhong, TANG Chunan, TAN Zhihong, et al. State of the art of inrush models in rock mass failure and developing trend for prediction and forecast of groundwater inrush[J].Chinese Journal of Rock Mechanics and Engineering, 2007,26(2): 268-277.

[5] Zhang J C. Investigations of water inrushes from aquifers under coal seam[J]. International Journal of Rock Mechanics and Mining Sciences, 2005, 42(2): 350-360.

[6] 劉招偉, 何滿潮, 王樹仁. 圓梁山隧道巖溶突水機理及防治對策研究[J]. 巖土力學(xué), 2006, 27(2): 228-232.LIU Zhaowei, HE Manchao, WANG Shuren. Study on karst water burst mechanism and prevention countermeasures in Yuanliangshan tunnel[J]. Rock and Soil Mechanics, 2006, 27(2):228-232.

[7] 李利平, 李術(shù)才, 張慶松. 巖溶地區(qū)隧道裂隙水突出力學(xué)機制研究[J]. 巖土力學(xué), 2010, 31(2): 523-528.LI Liping, LI Shucai, ZHANG Qingsong. Study of mechanism of water inrush induced by hydraulic fracturing in karst tunnels[J]. Rock and Soil Mechanics, 2010, 31(2): 523-528.

[8] 許振浩, 李術(shù)才, 李利平, 等. 一種典型的巖溶隧道襯砌壓裂突水災(zāi)害成因與防治[J]. 巖石力學(xué)與工程學(xué)報, 2011, 30(7):1396-1404.XU Zhenhao, LI Shucai, LI Liping, et al. Cause, disaster prevention and controlling of a typical kind of water inrush and lining fracturing in karst tunnels[J]. Chinese Journal of Rock Mechanics and Engineering, 2011, 30(7): 1396-1404.

[9] 國家安全生產(chǎn)監(jiān)督管理總局, 國家煤礦安全監(jiān)檢局. 煤礦防治水規(guī)定[M]. 北京: 煤炭工業(yè)出版社, 2009.Sate Administration of Work Safety, Sate Administration of Coal Safety. Rule of mine prevent ion and cure water disaster[M]. Beijing: China Coal Industry Publishing House,2009.

[10] 高保彬, 高佳佳, 袁東升. 基于UDEC 的大采高覆巖破裂的模擬與分析[J]. 湖南科技大學(xué)學(xué)報(自然科學(xué)版), 2013, 28(2):83-87.GAO Banbin, GAO Jiajia, YUAN Dongshen. Simulation and analysis of overlying strata fracture in lame mining height top-coal caving face of coal seams based on UDEC[J]. Journal of Hunan University of Science & Technology (Natural Science Edition), 2013, 28(2): 83-87.

[11] 劉超, 張義平, 齊飛, 等. 地面直流電法在煤礦防治水中的應(yīng)用[J]. 礦業(yè)工程研究, 2012, 27(3): 64-67.LIU Chao, ZHANG Yiping, QI Fei, et al. The application of ground DC method in mine water prevention and control[J].Mineral Engineering Research, 2012, 27(1): 54-58.

[12] Itasca Consulting Group Inc. FLAC3D users manual[R].Minneapolis: Itasca Consulting Group Inc., 2004.

[13] 王衛(wèi)軍, 趙延林, 李青鋒, 等. 礦井巖溶突水災(zāi)變機理[J]. 煤炭學(xué)報, 2010, 35(3): 443-448.WANG Weijun, ZHAO Yanlin, LI Qingfeng, et al. Disaster mechanism of karst water bursting in mine[J]. Journal of China Coal Society, 2010, 35(3): 443-448.

[14] 趙延林, 彭青陽, 萬文, 等. 高水壓下巖體裂紋擴展的滲流-斷裂耦合機制與數(shù)值實現(xiàn)[J]. 巖土力學(xué), 2014, 35(2): 556-564.ZHAO Yanlin, PENG Qingyang, WAN Wen, et al. Seepagefracture coupling mechanism of rock masses cracking propagation under high hydraulic pressure and numerical verification[J]. Rock and Soil Mechanics, 2014, 35(2): 556-564.

[15] 趙延林, 吳啟紅, 王衛(wèi)軍, 等. 基于突變理論的采空區(qū)重疊頂板穩(wěn)定性強度折減法及應(yīng)用[J]. 巖石力學(xué)與工程學(xué)報, 2010,29(7): 1412-1513.ZHAO Yanlin, WU Qihong, WANG Wenjun, et al. Strength reduction method to study stability of goaf overlapping roof based on catastrophe theory[J]. Chinese Journal of Geotechnical Engineering, 2010, 29(7): 1412-1513.