LDP曲線數(shù)值求解中波浪形發(fā)展問題的解決

李昌博， 姚磊華， 劉立鵬， 孫興松

(1.中國地質(zhì)大學(xué)(北京) 工程技術(shù)學(xué)院， 北京 100083； 2. 中國水利水電科學(xué)研究院， 北京 100048)

1 研究背景

隨著新奧法在地下工程施工中的廣泛采用，基于彈塑性分析的收斂約束方法在圍巖支護設(shè)計中逐漸受到重視[1]。獲得合理可用的圍巖縱向變形曲線(LDP曲線)是采用收斂約束方法設(shè)計地下洞室圍巖支護參數(shù)的基礎(chǔ)，目前業(yè)界主要是通過數(shù)值計算獲得該曲線。

國內(nèi)外很多學(xué)者均對LDP曲線展開過研究，Panet[2]給出了彈性條件下LDP曲線的擬合公式。Vlachopoulos等[3]采用二維軸對稱有限元模型，考慮圍巖塑性區(qū)半徑對于掌子面支撐效應(yīng)的影響，模擬分析得到考慮圍巖塑性區(qū)與洞室半徑間關(guān)系的LDP曲線公式。吳順川等[4]基于Hoek-Brown屈服準則，采用FLAC3D模擬軟件，研究了不同應(yīng)力條件和巖體質(zhì)量水平下隧道開挖的LDP曲線，并得到擬合公式。張傳慶等[5]以Hoek經(jīng)驗公式為基礎(chǔ)，基于現(xiàn)場圍巖監(jiān)測收斂位移，研究了隧道開挖過程中前期損失位移和LDP曲線的求取方法。孫鈞等[6]和蘇國韶[7]對應(yīng)力釋放進行了早期研究，隨后楊友彬等[8]基于應(yīng)力釋放率可近似等效為位移釋放率的原理，提出通過體積損失率求解應(yīng)力釋放率的方法。郭瑞等[9]研究了隧道開挖過程中應(yīng)力釋放與位移釋放的相關(guān)關(guān)系。杜雷功[10]對水工隧洞的施工關(guān)鍵技術(shù)進行了詳細闡述。趙大洲等[11]考慮隧洞施工方式和圍巖支護形式等因素，采用數(shù)值模擬分析了圍巖與襯砌的相互作用關(guān)系。

數(shù)值計算中，特別是地下洞室開挖圍巖穩(wěn)定性分析中，開挖步長與剖分精度間組合關(guān)系對計算結(jié)果存在較為明顯的影響，目前對采用數(shù)值計算方法獲得LDP曲線過程中開挖步長與剖分精度對于曲線的影響規(guī)律及如何消除影響未見相關(guān)研究，同時對于巖體開挖模擬方式以及巖體自重等因素對于LDP曲線最終形態(tài)等亦研究較少?；诖?，本文在收集整理現(xiàn)有圍巖縱剖面曲線的基礎(chǔ)上，采用數(shù)值模擬軟件對地下洞室開挖過程進行了模擬分析，詳盡研究開挖步長與剖分精度間不同組合關(guān)系對于圍巖縱剖面曲線形態(tài)的影響，并論證了巖體開挖模擬方式以及自重等對于LDP曲線形態(tài)以及圍巖最終收斂位移等的影響，給出了如何消除LDP曲線存在波浪形發(fā)展問題的方法，研究結(jié)果可為設(shè)計人員采用收斂-約束方法進行地下洞室支護設(shè)計中LDP曲線的求解提供一定的借鑒和參考。

2 已有圍巖縱剖面曲線

2.1 彈性力學(xué)解

Panet基于彈性力學(xué)理論分析給出了LDP曲線彈性解[2]：

(1)

式中：X*=X/R0，X為洞壁點距離掌子面的距離,m；R0為洞室半徑，m；uR為洞壁點的位移,cm；umax為洞壁最大位移,cm。

Unlu等[12]認為圍巖縱剖面曲線并不單純只與掌子面距離相關(guān)，同時與巖體泊松比ν亦相關(guān)，并提出不同經(jīng)驗公式：

對于X*≤0情況

(2)

對于X*≥0情況

(3)

其中:

Aa=-0.22ν-0.19;Ba=0.73ν+0.81

Ab=-0.22ν+0.81;Bb=0.39ν+0.65

2.2 基于實測值的擬合

Hoek對Chern等[13]實測結(jié)果進行擬合后，得到如下LDP曲線最佳擬合公式[14]：

(4)

同時，Hoek在其RocSupport軟件中又給出如下圍巖縱剖面曲線：

(5)

2.3 數(shù)值模擬獲取

Vlachopoulos等[3]研究發(fā)現(xiàn)圍巖塑性區(qū)對于圍巖縱剖面曲線具有較大影響，并提出不同塑性區(qū)RP與洞徑RT比值(R*=RP/RT)時LDP曲線：

對于X*≤0情況

(6)

對于X*≥0情況

(7)

由上述公式可知，基于不同的假設(shè)所獲得的圍巖縱剖面變形公式并不一致，在具體工程中難以直接確定采用哪個公式更為合適，另外由于每個工程所賦存的環(huán)境不同，如巖體強度、地應(yīng)力量級及方向等，使得理論推導(dǎo)公式或基于某個工程測量值的擬合公式均難以直接運用，通常仍需采用數(shù)值模擬計算的方式獲得具體工程的LDP曲線以用于后續(xù)支護設(shè)計中。

3 數(shù)值模型及邊界條件

此處以西南某隧洞工程為研究對象，建立數(shù)值模擬模型，并給定相同邊界條件，以研究開挖步長、開挖方法等的影響。利用FLAC3D的自帶建模命令建立隧洞三維計算模型，如圖1所示。

圖1 隧洞三維計算模型

模型沿隧洞軸線方向(y向)長100 m、寬(x向)100 m、高(z向)100 m，隧洞直徑D為10 m，沿隧洞軸線方向網(wǎng)格長度為1 m(0.1D)，共剖分281 600個單元，含289 365個節(jié)點。隧洞處于20 MPa的靜水應(yīng)力場中，x、y向施加相應(yīng)x、y向位移邊界限制，模型底面施加z向位移限制。巖體為理想彈塑性材料，屈服準則采用Mohr-Coloumb準則，物理力學(xué)參數(shù)如表1所示。

表1 巖體物理力學(xué)參數(shù)

4 LDP曲線影響及處理方法

4.1 開挖步長

鉆爆法開挖施工中，一般根據(jù)地下洞室斷面面積大小，采用臺階法分布開挖，而在數(shù)值模擬中通?？紤]盡量與施工方法較為接近，所以在求解LDP曲線時一般也采用分段開挖的方式進行巖體開挖，此處設(shè)置開挖長度分布為0.1D、0.2D、1D及一次性開挖完畢(8D)4種情況模擬地下洞室開挖(均采用zone cmodel assign null命令)。不同開挖步長下拱頂圍巖收斂變形以及歸一化處理后圍巖縱剖面變形曲線(LDP曲線)如圖2所示。

由圖2可知，巖體開挖步長對于LDP曲線形態(tài)影響很大，如采用0.2D、1D開挖長度時，曲線存在明顯波浪形態(tài)，這勢必使得采用收斂約束方法進行圍巖支護設(shè)計時無法采用該曲線來選擇合理的支護距離和支護體強度；開挖體長度為0.1D和8D時，即分別對應(yīng)1個剖分網(wǎng)格和一次性開挖完所有待開挖巖體時，LDP曲線形態(tài)較為光滑，無明顯波浪形態(tài)存在。故而，在采用數(shù)值模擬方法求解地下洞室?guī)r體開挖時圍巖LDP曲線，應(yīng)注意開挖步長對于圍巖縱向變形曲線形態(tài)的影響，可采用開挖長度為單個剖分單元體長度或一次性開挖完的方式，此時所獲得的歸一化圍巖LDP曲線較為理想，但需要注意的是采用1個單元長度或一次性開挖方式時，仍需要考慮該方式對于圍巖最終收斂變形量值的影響。

4.2 開挖模擬方式

在現(xiàn)有數(shù)值模擬分析軟件中，對于巖體開挖的模擬一般是通過將需要挖去的單元體直接刪除或“殺死”等方式進行，如FLAC3D軟件中采用delete或zone cmodel assign null命令消除待開挖巖體，再將模型計算至平衡狀態(tài)，以此模擬隧洞開挖過程。FLAC3D軟件中一直存在大體積巖體開挖慣性效應(yīng)導(dǎo)致的塑性區(qū)體積過大問題，最新版本的FLAC3D軟件中為防止在靜力求解過程中出現(xiàn)非真實慣性效應(yīng)導(dǎo)致的洞壁圍巖非真實性的破壞，采用zone relax命令給待開挖單元體乘以0～1之間的某一系數(shù)，對稱性地逐漸降低開挖單元剛度、應(yīng)力以及密度，當(dāng)系數(shù)降為0時單元本構(gòu)模型變?yōu)閚ull，同時移除施加于單元體上所有力或位移邊界，以此模擬開挖過程。圖2顯示采用無應(yīng)力釋放方式開挖時，0.1D和8D開挖步長對于LDP曲線形態(tài)基本無影響，但最終位移收斂量值存在明顯不同，此處采用zone cmodel assign null和zone relax兩種不同的命令來模擬0.1D和8D兩種不同開挖長度下圍巖LDP曲線形態(tài)，具體如圖3所示。

由圖3可知，LDP曲線在距離掌子面5倍洞徑外逐漸趨于穩(wěn)定，F(xiàn)LAC3D軟件常規(guī)的zone cmodel assign null命令或者新版軟件中zone relax命令所計算得到的圍巖縱剖面變形曲線基本相同。采用不同模擬方式下洞壁圍巖收斂變形量值如圖4所示。

由圖4可知，同一開挖步長下采用zone relax和zone cmodel assign null兩種不同命令計算時，圍巖收斂位移表現(xiàn)出不同的變化形式。一次性挖完時(8D)，采用zone relax命令計算時所得到的圍巖收斂位移在掌子面后與采用zone cmodel assign null命令計算值具有較大差異，并隨著距離增加而逐漸加大。這將給采用收斂-約束方法進行支護設(shè)計時帶來很大的困惑，即相同支護系統(tǒng)、相同安裝時機卻得到不同的支護安全系數(shù)。而采用0.1D長度開挖時(每次開挖一個單元長度)，兩種模擬方式得到的圍巖收斂位移在任何時刻均較為一致，且所得到圍巖最終收斂位移與彈塑性力學(xué)芬納解更為接近，所以推薦采用0.1D開挖方式獲得更為合理的曲線形態(tài)和圍巖最終收斂位移，以便于支護設(shè)計選擇。

4.3 其他因素

LDP曲線提出時假設(shè)巖體賦存于靜水應(yīng)力場中，但是多數(shù)學(xué)者在利用數(shù)值模擬方法求解圍巖變形時，通常會設(shè)定巖體密度(重度)，用以考慮在自重應(yīng)力作用下巖體的變形，此處對相同應(yīng)力場情況下，仿真模擬了有、無巖體自重兩種不同的情況，以分析自重對LDP曲線和圍巖最終位移的影響，結(jié)果如圖5所示。

圖2 不同開挖步尺寸下LDP曲線

圖3 兩種命令計算的不同開挖方式下LDP曲線 圖4 兩種命令計算的不同開挖方式下LDP曲線

圖5 巖體自重對LDP曲線的影響

由圖5可知，在數(shù)值模擬中考慮自重應(yīng)力作用時圍巖變形略大于不設(shè)置巖體密度(自重)的靜水壓力場情況，較高地應(yīng)力作用圍巖屈服范圍更大時這種差異會變大，但之間量值相差總體較小，LDP曲線形態(tài)及發(fā)展型式基本相近。所以，采用數(shù)值計算的方法獲取LDP曲線時，是否設(shè)置自重應(yīng)力場對于LDP曲線的影響較小，基本可忽略不計。

5 結(jié) 論

本文對采用數(shù)值模擬技術(shù)獲取圍巖縱剖面變形曲線時如何避免曲線波浪形變化以及開挖模擬方法、設(shè)置巖體自重等對于LDP曲線的影響進行了研究，主要得出以下結(jié)論：

(1)采用數(shù)值模擬方法求解圍巖開挖縱剖面變形曲線時，巖體開挖步長對曲線形態(tài)影響較大，曲線易于存在明顯波浪形態(tài)。開挖長度為單個剖分單元體長度或一次性開挖完的方式可有效避免該問題，但仍需考慮不同方式對于圍巖最終收斂變形的影響。

(2)開挖步數(shù)越長，采用直接刪除單元或應(yīng)力釋放法圍巖最終收斂變形計算值差距越大，開挖長度與剖分單元長度相同時，兩種模擬方式得到的收斂位移時刻一致，且彈塑性力學(xué)解更為接近，推薦采用該方法獲取LDP曲線，以便于支護設(shè)計選擇。

(3)采用數(shù)值計算獲取LDP曲線時，是否設(shè)置自重應(yīng)力場對于曲線影響較小，基本可忽略不計。

本文對獲得合理可用LDP曲線進行了研究，可為隧洞支護設(shè)計提供一定的參考，但由于隧洞工程多穿越不同巖性及質(zhì)量地層，如何考慮巖體參數(shù)隨機分布對于LDP曲線的影響需要進一步開展研究。