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

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

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

1 研究背景

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

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

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

2 已有圍巖縱剖面曲線(xiàn)

2.1 彈性力學(xué)解

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

(1)

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

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

對(duì)于X*≤0情況

(2)

對(duì)于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 基于實(shí)測(cè)值的擬合

Hoek對(duì)Chern等[13]實(shí)測(cè)結(jié)果進(jìn)行擬合后，得到如下LDP曲線(xiàn)最佳擬合公式[14]：

(4)

同時(shí)，Hoek在其RocSupport軟件中又給出如下圍巖縱剖面曲線(xiàn)：

(5)

2.3 數(shù)值模擬獲取

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

對(duì)于X*≤0情況

(6)

對(duì)于X*≥0情況

(7)

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

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

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

圖1 隧洞三維計(jì)算模型

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

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

4 LDP曲線(xiàn)影響及處理方法

4.1 開(kāi)挖步長(zhǎng)

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

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

4.2 開(kāi)挖模擬方式

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

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

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

4.3 其他因素

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

圖2 不同開(kāi)挖步尺寸下LDP曲線(xiàn)

圖3 兩種命令計(jì)算的不同開(kāi)挖方式下LDP曲線(xiàn) 圖4 兩種命令計(jì)算的不同開(kāi)挖方式下LDP曲線(xiàn)

圖5 巖體自重對(duì)LDP曲線(xiàn)的影響

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

5 結(jié) 論

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

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

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

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

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