基于改進(jìn)迭代求解的雙孔洞圍巖應(yīng)力場(chǎng)解析研究

蔡鑫，陳立業(yè)，周子龍，譚力海，陽(yáng)俊

(1. 中南大學(xué) 資源與安全工程學(xué)院，湖南 長(zhǎng)沙，410083；2. 南華大學(xué) 資源環(huán)境與安全工程學(xué)院，湖南 衡陽(yáng)，421000；3. 湖南有色金屬職業(yè)技術(shù)學(xué)院 資源環(huán)境系，湖南 株洲，412000)

地下工程中常包含不同形態(tài)的孔洞，如硐室、隧(巷)道、溶洞以及采空區(qū)等。由于孔洞周圍的應(yīng)力集中引起巖體變形與裂紋萌生擴(kuò)展，極易誘發(fā)巖體整體失穩(wěn)[1-3]。在地下工程中，局部區(qū)域內(nèi)往往存在多個(gè)孔洞，這些地下孔洞幾何形態(tài)與空間分布特征各異，當(dāng)孔洞間距較小時(shí)，彼此會(huì)相互影響，致使其周圍應(yīng)力場(chǎng)分布特征十分復(fù)雜[4-5]。為揭示外荷載作用下多孔洞圍巖系統(tǒng)的損傷破裂特征及機(jī)制，國(guó)內(nèi)外學(xué)者進(jìn)行了大量研究，如YANG 等[6-9]開(kāi)展了外荷載作用下圓形、橢圓、矩形雙孔洞的室內(nèi)試驗(yàn)與數(shù)值仿真研究，獲得了不同形狀孔洞的方位角及相對(duì)位置對(duì)巖石力學(xué)性質(zhì)及裂紋演化的影響規(guī)律。然而，以上研究局限于孔洞對(duì)巖石試樣力學(xué)性質(zhì)與破壞模式的影響，對(duì)相鄰孔洞相互作用機(jī)理研究較少，更缺乏對(duì)孔洞體系周邊圍巖應(yīng)力分布的定量分析。實(shí)際上，相較于室內(nèi)實(shí)驗(yàn)和數(shù)值模擬，應(yīng)力解析方法可靠性高，能對(duì)地下孔洞相關(guān)的工程問(wèn)題進(jìn)行快捷計(jì)算，并對(duì)相關(guān)設(shè)計(jì)參數(shù)進(jìn)行分析，可為地下工程設(shè)計(jì)與施工提供重要參考，具有重要的理論意義與應(yīng)用價(jià)值。在含孔洞圍巖應(yīng)力分析研究中，對(duì)非圓類復(fù)雜幾何形狀，常用復(fù)變函數(shù)方法求解以獲得高精度的解析解[10-14]。目前，對(duì)單個(gè)孔洞問(wèn)題的復(fù)變函數(shù)解析理論較完善[14-17]，但對(duì)于2個(gè)或多個(gè)孔洞的情況，由于鄰近孔洞之間的相互影響，此時(shí)圍巖的應(yīng)力分布與單個(gè)孔洞情況存在顯著區(qū)別，其應(yīng)力函數(shù)也變得更加復(fù)雜。為此，JEFFERY[18]采用雙極坐標(biāo)法研究了含雙圓孔無(wú)限域的應(yīng)力場(chǎng)，并給出了雙坐標(biāo)下雙調(diào)和應(yīng)力函數(shù)的一般形式。在此基礎(chǔ)上，LING[19]對(duì)含雙圓孔無(wú)限域的應(yīng)力場(chǎng)進(jìn)行了研究。KOOI等[20]對(duì)深埋雙平行巷道的應(yīng)力場(chǎng)進(jìn)行了計(jì)算。MINDLIN[21]對(duì)孔洞應(yīng)力場(chǎng)進(jìn)行了初始化處理，并推導(dǎo)了考慮重力場(chǎng)的應(yīng)力函數(shù)。朱大勇等[22-23]基于應(yīng)力疊加原理，將含多個(gè)孔洞的無(wú)限域中應(yīng)力場(chǎng)當(dāng)作各個(gè)單獨(dú)孔洞在洞口內(nèi)邊界上作用虛擬面力產(chǎn)生的應(yīng)力場(chǎng)的疊加。此外，Schwarz交替法也是解決含多孔洞的多連通域問(wèn)題的一個(gè)重要方法。早期受限于當(dāng)時(shí)的計(jì)算條件，只能完成1 次迭代，其求解精度較低[24]。近年來(lái)，隨著計(jì)算機(jī)計(jì)算能力提升，張路青等[25]利用Schwarz 交替法對(duì)雙圓形孔洞體系區(qū)域應(yīng)力場(chǎng)求解，通過(guò)2次迭代，顯著提高了求解精度。

基于迭代求解雙孔洞體系應(yīng)力場(chǎng)的核心問(wèn)題是通過(guò)求解每步迭代過(guò)程中的應(yīng)力解析函數(shù)得到在孔洞邊界產(chǎn)生的多余面力，而這涉及復(fù)變應(yīng)力解析函數(shù)或孔洞邊界坐標(biāo)點(diǎn)在不同坐標(biāo)系的轉(zhuǎn)換問(wèn)題。對(duì)于圓形孔洞，其映射函數(shù)和逆映射函數(shù)表達(dá)式簡(jiǎn)單，可直接對(duì)應(yīng)力解析函數(shù)進(jìn)行上述轉(zhuǎn)換。而對(duì)非圓孔洞，其形狀需要通過(guò)復(fù)雜的映射函數(shù)進(jìn)行描述，對(duì)應(yīng)的應(yīng)力解析函數(shù)也十分復(fù)雜。對(duì)于雙孔洞體系這種多連通域問(wèn)題，Schwarz交替法局限于圓形孔洞體系問(wèn)題的求解，而對(duì)于復(fù)雜形狀多孔洞問(wèn)題缺少有效的求解方法。為此，張路青等[26]利用孔洞邊界坐標(biāo)點(diǎn)的轉(zhuǎn)換代替復(fù)應(yīng)力解析函數(shù)的轉(zhuǎn)換，根據(jù)不同孔洞的映射函數(shù)、逆映射函數(shù)以及孔洞相對(duì)位置關(guān)系對(duì)孔洞邊界多余面力進(jìn)行求解。Lü 等[27]針對(duì)雙直墻半圓拱形孔洞體系，通過(guò)傅里葉級(jí)數(shù)形式將多余面力的函數(shù)展開(kāi)，并進(jìn)行20 余次迭代，從而獲得高精度的應(yīng)力解析解。ZHANG等[28]借鑒此方法獲得了無(wú)限域中含2個(gè)以上橢圓孔洞體系的應(yīng)力場(chǎng)解析解。卞曉琳等[29]進(jìn)一步結(jié)合復(fù)變函數(shù)解析拓延法，對(duì)2個(gè)相鄰深埋隧道的應(yīng)力場(chǎng)和位移場(chǎng)進(jìn)行求解。以上研究極大推動(dòng)了復(fù)雜形狀雙孔洞體系應(yīng)力解析理論的發(fā)展，但由于非圓孔洞的逆映射函數(shù)十分復(fù)雜，特別是當(dāng)孔洞形狀與圓形差異較大時(shí)，運(yùn)用現(xiàn)有方法求解應(yīng)力函數(shù)時(shí)求解精度較低。

本文在Schwarz交替迭代法的基礎(chǔ)上，提出一種新的理論解析方法。通過(guò)引入最優(yōu)化技術(shù)，直接完成不同坐標(biāo)系平面中孔洞邊界坐標(biāo)點(diǎn)的變換，無(wú)需在求解過(guò)程中引入逆映射函數(shù)，可避免逆映射函數(shù)求解造成的誤差，在理論上可顯著提升求解精度。基于此方法，以拱形孔洞為例，構(gòu)建無(wú)限平面域含單拱形孔洞和雙拱形孔洞體系的數(shù)學(xué)模型，獲得拱形孔洞周邊應(yīng)力場(chǎng)分布規(guī)律，并利用數(shù)值模擬方法對(duì)解析結(jié)果進(jìn)行驗(yàn)證與拓展。

1 基于交替迭代法求解雙空洞體系應(yīng)力場(chǎng)

1.1 交替迭代法求解雙孔洞體系應(yīng)力解基本原理

含任意形狀雙孔洞模型平面示意圖如圖1 所示。考慮彈性半空間內(nèi)存在2個(gè)任意形狀孔洞(孔洞1與孔洞2)且分布在任意位置，以孔洞1、孔洞2的形心為坐標(biāo)系原點(diǎn)分別建立直角坐標(biāo)系Z1(x1O1y1)和Z2(x2O2y2)，對(duì)應(yīng)的映射函數(shù)分別為z1=ω1(ζ1)與z2=ω2(ζ2)(其中，ζ1與ζ2分別為含孔洞1 平面與含孔洞2平面在含單位圓映射平面上點(diǎn)的坐標(biāo)，計(jì)算時(shí)使用復(fù)數(shù)形式ρeiθ，ρ為實(shí)部，i 為虛數(shù)，下同)。令遠(yuǎn)場(chǎng)水平應(yīng)力、垂直應(yīng)力與剪切應(yīng)力分別為σx∞、σy∞與τx∞y，根據(jù)交替迭代法原理，孔洞2 應(yīng)力求解過(guò)程如下。

圖1 含任意形狀雙孔洞模型平面示意圖Fig.1 Schematic of a plate model containing two holes with arbitrary shapes

1) 假設(shè)當(dāng)孔洞1 開(kāi)挖完成時(shí)，孔洞2 尚未開(kāi)挖，在遠(yuǎn)場(chǎng)應(yīng)力作用下，只存在孔洞1的應(yīng)力。

2) 計(jì)算坐標(biāo)系Z1下孔洞1開(kāi)挖后虛擬的孔洞2邊界的多余應(yīng)力。

3) 下一步開(kāi)挖孔洞2時(shí)，必須給孔洞2邊界施加外力以平衡孔洞1 在孔洞2 邊界產(chǎn)生的多余面力。假設(shè)此時(shí)孔洞1不存在，將坐標(biāo)系轉(zhuǎn)為Z2，則可根據(jù)孔洞2邊界應(yīng)力平衡條件計(jì)算相應(yīng)的應(yīng)力函數(shù)，而此時(shí)孔洞2 又會(huì)在孔洞1 邊界產(chǎn)生多余面力，需要再轉(zhuǎn)換到Z1坐標(biāo)系，計(jì)算新的邊界平衡條件下關(guān)于孔洞1的應(yīng)力函數(shù)。

隨著迭代次數(shù)增加，某一孔洞在相鄰孔洞邊界產(chǎn)生的多余面力會(huì)越來(lái)越小，求解結(jié)果逐漸逼近真實(shí)結(jié)果。一般地，當(dāng)求解精度滿足實(shí)際要求時(shí)，便可停止迭代計(jì)算，完成應(yīng)力求解。

1.2 逆映射函數(shù)未知情況下多余面力的求解

針對(duì)非圓形孔洞映射函數(shù)與逆映函數(shù)求解復(fù)雜、誤差大等問(wèn)題[30]，本研究提出在不考慮逆映射函數(shù)的情況下，采用交替迭代求解復(fù)雜形狀雙孔洞體系外部應(yīng)力解析解的改進(jìn)方法。

假設(shè)只存在孔洞1時(shí)，令解析函數(shù)為φ10(ζ1)與ψ10(ζ1)，則孔洞應(yīng)力邊界條件為

式中：σ1為ζ1平面單位圓周上點(diǎn)的坐標(biāo)；t1為Z1坐標(biāo)系下孔洞1邊界上點(diǎn)的坐標(biāo)，t1=ω1(σ1)/′(σ1)，其中′(σ1)為ω1(σ1)的共軛導(dǎo)數(shù)；(σ1)為φ10(σ1)的共軛導(dǎo)數(shù)；(σ1)為ψ10(σ1)的共軛。

φ10(ζ1)與ψ10(ζ1)可表示為

式中：φ01(ζ1)與ψ01(ζ1)可以通過(guò)Laurent 級(jí)數(shù)表示。常數(shù)B、B′與C′反映了遠(yuǎn)場(chǎng)圍巖應(yīng)力的情況。

求出解析函數(shù)φ10(ζ1)與ψ10(ζ1)的解后，需消除孔洞2 邊界的多余面力，其核心問(wèn)題是利用φ10(ζ1)與ψ10(ζ1)得到多余面力分布函數(shù)。在此過(guò)程中，需對(duì)ζ2平面的單位圓上點(diǎn)σ2進(jìn)行坐標(biāo)系變換，以得到孔洞2邊界多余面力f112(σ2)：

式中：γ1為點(diǎn)σ2經(jīng)過(guò)坐標(biāo)系變換到ζ1平面中的對(duì)應(yīng)點(diǎn)的坐標(biāo)。

早期研究表明，對(duì)類圓形孔洞，由其逆映射函數(shù)得到的近似單位圓與實(shí)際單位圓基本吻合。然而，對(duì)于矩形、梯形、拱形等含棱角形狀，其逆映射函數(shù)求解精度較低[29]。為此，本文引入低維變量最優(yōu)化技術(shù)，通過(guò)不同平面與不同坐標(biāo)系中采樣點(diǎn)坐標(biāo)轉(zhuǎn)換，替代逆映射函數(shù)對(duì)σ2點(diǎn)對(duì)應(yīng)的γ1點(diǎn)以及σ1點(diǎn)對(duì)應(yīng)的γ2點(diǎn)進(jìn)行坐標(biāo)求解，最終得到對(duì)應(yīng)的多余面力表達(dá)式，避免了逆映射函數(shù)的繁瑣求解及與理論解之間的誤差。下面以計(jì)算σ2點(diǎn)對(duì)應(yīng)γ1點(diǎn)的坐標(biāo)為例，對(duì)該方法步驟進(jìn)行簡(jiǎn)要說(shuō)明。

1) 確定孔洞2 對(duì)應(yīng)單位圓上σ2點(diǎn)采樣個(gè)數(shù)n，則采樣間距(Δσ2)為Δσ2=2π/n。第j(j=1，2，3，…，n)個(gè)采樣點(diǎn)的坐標(biāo)為σj2=eiΔσ2×j。

2) 根據(jù)映射函數(shù)z2=ω2(σ2)確定σ2采樣點(diǎn)在Z2坐標(biāo)系中孔洞2邊界上點(diǎn)t2的坐標(biāo)。

3) 計(jì)算點(diǎn)t2在Z1坐標(biāo)系中的對(duì)應(yīng)點(diǎn)坐標(biāo)η1=t2+c。

4) 利用最優(yōu)化技術(shù)求解點(diǎn)η1在ζ1平面對(duì)應(yīng)的映射點(diǎn)γ1，如圖2所示。令初始域(圖2中著色扇形區(qū)域)角度范圍為[αa，αb]，半徑范圍為[ra，rb]。將初始域角度以0.5(αa+αb)為增量，半徑以0.5(ra+rb)為增量，分為4 個(gè)子域，通過(guò)比較每個(gè)區(qū)域的端點(diǎn)值，確定最優(yōu)值所在的子域(圖中紅色區(qū)域)。之后，再將子區(qū)域按照前面的方法分為4 個(gè)子區(qū)域，如此反復(fù)迭代，直至區(qū)域中心值達(dá)到設(shè)定的精度為止。

圖2 利用最優(yōu)化技術(shù)求解映射點(diǎn)坐標(biāo)Fig. 2 Determination of mapping point coordinates by using optimization method

5) 依次確定采樣點(diǎn)在ζ1平面中的坐標(biāo)值，計(jì)算各點(diǎn)對(duì)應(yīng)的多余面力。

同理，可根據(jù)以上步驟計(jì)算σ1點(diǎn)經(jīng)過(guò)坐標(biāo)系變換到ζ2平面中的對(duì)應(yīng)點(diǎn)γ2，僅需改變相應(yīng)的映射函數(shù)表達(dá)式即可。

當(dāng)確定了孔洞采樣點(diǎn)在映射單位圓上的多余面力之后，可利用傅里葉復(fù)數(shù)級(jí)數(shù)形式表示多余面力的分布。以f112(σ2)為例，令其傅里葉復(fù)數(shù)級(jí)數(shù)表達(dá)式為

將采樣點(diǎn)坐標(biāo)值代入式(5)，則可以建立高維方程組，求解可得f112(σ2)的傅里葉級(jí)數(shù)表達(dá)式。

1.3 復(fù)應(yīng)力函數(shù)迭代求解

在確定f112(σ2)表達(dá)式之后，在第1次迭代計(jì)算中，孔洞2對(duì)應(yīng)單位圓的應(yīng)力邊界條件可表示為

由式(6)求解得到的孔洞2 的應(yīng)力解析函數(shù)φ12(ζ2)、ψ21(ζ2)與孔洞1 的應(yīng)力解析函數(shù)φ10(ζ1)、ψ10(ζ1)疊加即可得到第1次迭代計(jì)算得到的求解結(jié)果，至此，任意形狀雙孔洞的第1次迭代完成。

由于迭代求解是從孔洞1開(kāi)始的，故在第1次迭代后，孔洞2邊界已滿足零面力的邊界條件，但孔洞1邊界會(huì)產(chǎn)生新的多余面力f211(σ1)，故需要繼續(xù)在孔洞1的邊界施加新的反面力進(jìn)行平衡，其計(jì)算過(guò)程與首次的迭代過(guò)程一致，在此后的每次迭代過(guò)程中，只需根據(jù)其上一步的迭代求解結(jié)果對(duì)雙孔洞應(yīng)力解析函數(shù)進(jìn)行更新，直至兩洞周邊附加面力在要求范圍之內(nèi)。

2 算例分析與數(shù)值模擬驗(yàn)證

以巖體工程中常見(jiàn)的拱形隧道為例，利用提出的改進(jìn)迭代求解方法計(jì)算雙拱形孔洞圍巖應(yīng)力場(chǎng)，并采用有限元模擬對(duì)理論解析結(jié)果進(jìn)行對(duì)比驗(yàn)證。隨后，利用離散元軟件PFC2D 對(duì)含雙拱形孔洞圍巖破壞模式進(jìn)行分析。

2.1 拱形孔洞映射函數(shù)求解

以拱形孔洞為研究對(duì)象(圖3)，其具體尺寸為：R1=6.30 m，R2=8.80 m，R3=1.60 m，R4=15.60 m，L1=12.25 m，L2=2.50 m，L3=4.64 m，L4=0.96 m，坐標(biāo)原點(diǎn)為孔洞的形心。

圖3 拱形孔洞示意圖Fig. 3 Schematic of an arched hole

Z平面上孔洞外區(qū)域變化到ζ平面上單位圓外區(qū)域的映射函數(shù)表達(dá)式為[30]

式中：R為實(shí)數(shù)，與Z平面上孔洞尺寸有關(guān)。當(dāng)孔洞關(guān)于x軸對(duì)稱時(shí)，Ck為實(shí)數(shù)，c為孔洞邊界的積分路徑。一般地，孔洞形狀越復(fù)雜，就需要越大的k來(lái)提高映射形狀的精度，即確定R和Ck后便可得到孔洞的映射函數(shù)。

令孔洞周邊任一點(diǎn)Aj極坐標(biāo)為(rj，αj)，其在ζ平面上單位圓上對(duì)應(yīng)點(diǎn)Bj極坐標(biāo)為(1，θj)，則由式(7)可得

由式(8)可得

根據(jù)歐拉公式，將式(9)中實(shí)部與虛部分開(kāi)得：

對(duì)平面上單位圓上均勻分布的m個(gè)點(diǎn)進(jìn)行采樣，當(dāng)孔洞自身關(guān)于x軸對(duì)稱時(shí)，可只對(duì)范圍為[0，π]的單位圓進(jìn)行采樣。令孔洞中心為坐標(biāo)原點(diǎn)，Bj極坐標(biāo)為(1，θj)，其通過(guò)式(11)確定的Z平面對(duì)應(yīng)點(diǎn)Aj坐標(biāo)為(rj，αj)，實(shí)際對(duì)應(yīng)點(diǎn)A′j坐標(biāo)為(r′j，αj)。此時(shí)，采樣點(diǎn)與目標(biāo)點(diǎn)的平均離差平方和f為

本文采用Box復(fù)合型最優(yōu)化算法求解映射函數(shù)的系數(shù)。該算法的核心在于在n維空間可行域內(nèi)構(gòu)建具有N(N≥n+1)個(gè)頂點(diǎn)的多面體(又稱復(fù)合形)，然后逐一比較復(fù)合形各頂點(diǎn)目標(biāo)函數(shù)值。在每一步計(jì)算中，獲取并刪除復(fù)合形中使目標(biāo)函數(shù)值最大的頂點(diǎn)(壞點(diǎn))，同時(shí)生成使目標(biāo)函數(shù)值更小的新點(diǎn)。如此反復(fù)迭代，使得復(fù)合形不斷收縮，其頂點(diǎn)不斷移向最優(yōu)點(diǎn)，直至滿足中止迭代條件，本文中設(shè)定的精度為5×10-4[31]。復(fù)合形法求解過(guò)程中孔洞映射形狀的演化見(jiàn)圖4。在迭代求解過(guò)程中，隨著平均離差平方和f不斷降低，根據(jù)復(fù)合形中當(dāng)前最優(yōu)點(diǎn)Xb得到的映射孔洞形狀逐漸趨近于目標(biāo)拱形孔洞形狀。

圖4 復(fù)合形法求解過(guò)程中孔洞映射形狀的演化Fig. 4 Evolution of the mapping shape during the solving process using complex method

最終得到的拱形映射函數(shù)表達(dá)式為

根據(jù)式(13)得到的拱形孔洞映射形狀與目標(biāo)形狀進(jìn)行對(duì)比，見(jiàn)圖5。從圖5 可以看到：拱形孔洞的映射形狀與實(shí)際形狀高度吻合，證明了復(fù)合型最優(yōu)化算法的適用性。

圖5 孔洞映射形狀與目標(biāo)形狀比較Fig. 5 Comparison between the mapping shape and original shape of the hole

2.2 雙拱形孔洞體系應(yīng)力場(chǎng)解析結(jié)果與數(shù)值驗(yàn)證

雙拱形孔洞體系布局示意圖如圖6所示。雙拱形孔洞體系以孔洞1形心為中心，兩孔洞形心距為9 m，兩孔洞形心連線與水平線夾角(下稱連接角)為θ。為研究?jī)煽锥催B接角的影響，取θ的范圍為0°～90°，以15°間隔遞增。令σx∞=p，σy∞=τx∞y=0，利用本文提出的改進(jìn)迭代求解方法計(jì)算雙拱形孔洞周邊應(yīng)力分布。

圖6 雙拱形孔洞體系布局示意圖Fig. 6 Schematic of the two arched holes system

為驗(yàn)證應(yīng)力場(chǎng)解析解的準(zhǔn)確性，采用有限元分析軟件ANSYS建立二維平面應(yīng)變模型，獲得側(cè)壓系數(shù)為0時(shí)的雙拱形孔洞應(yīng)力分布規(guī)律，并將其與理論解析結(jié)果進(jìn)行比較。模型中，圍巖材料選用線彈性本構(gòu)模型，單元類型為PLANE183 平面單元，巖層彈性模量取5 GPa，泊松比取0.3，垂直方向遠(yuǎn)場(chǎng)應(yīng)力σx∞取1 MPa。為減少邊界效應(yīng)的影響，模型尺寸需大于孔洞尺寸的10 倍，故平面數(shù)值模型的尺寸設(shè)為300 m×300 m。

當(dāng)θ為0°時(shí)，孔洞1周邊環(huán)向應(yīng)力的分布規(guī)律(以壓應(yīng)力為正)見(jiàn)圖7。從圖7 可見(jiàn)：整體上，理論解析與數(shù)值模擬結(jié)果變化趨勢(shì)近乎一致，這說(shuō)明本文提出的雙孔洞應(yīng)力解析解方法是準(zhǔn)確、可靠的。值得注意的是，在局部位置，特別是左右邊墻底端拐角處解析結(jié)果與模擬結(jié)果存在很小的差異，其原因是：一方面，在數(shù)值模擬中，孔洞并不是處于無(wú)限域之中，邊界對(duì)孔洞周邊應(yīng)力場(chǎng)可能會(huì)產(chǎn)生一定影響；另一方面，受計(jì)算精度影響，映射形狀與目標(biāo)拱形孔洞形狀也并不完全吻合，但這種差別很小，可以忽略不計(jì)。

圖7 孔洞1周邊環(huán)向應(yīng)力分布理論解析與數(shù)值模擬結(jié)果比較Fig.7 Comparison between analytical and numerical results of stress around the hole 1

2.3 連接角對(duì)孔周應(yīng)力分布的影響規(guī)律

在無(wú)側(cè)壓情況下，單拱形孔邊界以及雙拱形孔體系中孔洞1邊界環(huán)向應(yīng)力分布隨連接角θ的變化見(jiàn)圖8。圖8 中，位于孔洞外部的應(yīng)力為正，即為壓應(yīng)力，反之則為拉應(yīng)力；紅點(diǎn)和綠點(diǎn)分別表示最大壓力、拉應(yīng)力所在位置。

圖8 單孔洞與不同布局雙孔洞體系中孔洞1邊界環(huán)向應(yīng)力解析解Fig.8 Hoop stress results around the boundary of a single hole and hole 1 in the two-hole system with different layouts

由圖8(h)可見(jiàn)：當(dāng)連接角θ為90°時(shí)，孔洞1邊界的切向應(yīng)力呈對(duì)稱布置；但在其他角度下，孔洞2的存在影響孔洞1周邊應(yīng)力的對(duì)稱性；當(dāng)連接角為0°～60°時(shí)，孔洞1環(huán)向最大壓應(yīng)力在孔洞左側(cè)拱底位置，而當(dāng)連接角為75°時(shí)，孔洞1 左側(cè)拱底環(huán)向壓力顯著降低，最大壓應(yīng)力出現(xiàn)在孔洞右側(cè)；隨著連接角的變化，孔周最大拉伸應(yīng)力位置會(huì)在孔洞頂、底中心位置附近變動(dòng)，但無(wú)明顯變化規(guī)律；當(dāng)連接角為0°、45°、60°、75°時(shí)，最大拉應(yīng)力在孔底中心位置附近，而在其他角度下，最大拉應(yīng)力在孔頂中心位置。此外，從圖8還可見(jiàn)：當(dāng)2個(gè)孔洞沿斜線布置時(shí)，靠近相鄰孔洞一側(cè)與遠(yuǎn)離相鄰孔洞一側(cè)的邊界應(yīng)力呈現(xiàn)明顯的差異；在靠近孔洞2 的側(cè)邊(如左側(cè)邊墻和孔底)，環(huán)向應(yīng)力隨連接角變化明顯，而在遠(yuǎn)離孔洞2 的側(cè)邊(右側(cè)邊墻和孔頂)區(qū)域，應(yīng)力變化幅度相對(duì)較小，說(shuō)明相鄰2個(gè)孔洞之間的相互影響主要集中在相鄰一側(cè)的區(qū)域。

雙拱形洞體系中孔洞1最大解析拉、壓應(yīng)力集中系數(shù)隨著θ的變化規(guī)律見(jiàn)圖9。由圖9可見(jiàn)：

1) 隨著θ增加，雙拱形洞體系孔1周邊最大拉伸應(yīng)力集中系數(shù)先增大后減小；當(dāng)θ為45°時(shí)，孔洞1的頂板區(qū)域拉應(yīng)力達(dá)到最大。

2) 最大壓應(yīng)力隨著θ增加表現(xiàn)出先降低、后增加再快速降低的變化趨勢(shì)；當(dāng)θ為30°時(shí)，孔洞1邊界切向壓應(yīng)力達(dá)到最大，此后逐漸降低，直到θ為90°時(shí)，切向壓應(yīng)力集中系數(shù)最小，僅為2 個(gè)孔水平布置時(shí)(θ=0)的64.2%。

2.4 連接角對(duì)雙拱孔洞體系承載特性與破壞模式的影響

離散元數(shù)值仿真可通過(guò)微觀尺度研究孔洞周邊的微裂紋的萌生和演化，從而反映出孔周應(yīng)力場(chǎng)的分布情況。為研究連接角對(duì)雙拱孔洞體系承載特性與破裂模式的影響，本文采用離散元軟件PFC2D開(kāi)展含雙拱形孔洞體系巖體數(shù)值模擬。孔洞形狀、尺寸及其布局均與理論求解中的模型一致，模型長(zhǎng)×寬為100 m×100 m，側(cè)壓系數(shù)為0(如圖10所示)，豎直方向以恒定速率加載直至模型完全被破壞。采用PFC2D模擬時(shí)，首先以某砂巖為參照對(duì)象，基于該砂巖標(biāo)準(zhǔn)單軸壓縮試驗(yàn)結(jié)果對(duì)數(shù)值模型微觀參數(shù)進(jìn)行標(biāo)定，各參數(shù)取值如表1所示。圖11 所示為砂巖室內(nèi)試驗(yàn)單軸壓縮試驗(yàn)與數(shù)值模擬結(jié)果。從圖11可見(jiàn)：數(shù)值模型除了無(wú)壓密階段外，其強(qiáng)度和彈性模量都與室內(nèi)試驗(yàn)結(jié)果一致；同時(shí)，在破壞模式上，兩者都屬于斜剪切破壞，這說(shuō)明可利用表1中參數(shù)很好地模擬砂巖的力學(xué)行為。

表1 PFC2D離散元數(shù)值模型微觀參數(shù)Table 1 Micro-parameters for the PFC2D numerical model

圖10 離散元數(shù)值模型Fig.10 Discrete element numerical models

圖11 數(shù)值模型標(biāo)定結(jié)果Fig.11 Calibration results for numerical model

數(shù)值模擬中含雙拱形孔洞試樣峰值應(yīng)力隨連接角θ的變化規(guī)律見(jiàn)圖12。從圖12可見(jiàn)：隨著θ增加，試樣的峰值應(yīng)力在15°時(shí)增加，隨后降低；當(dāng)θ超過(guò)30°時(shí)，試樣的峰值應(yīng)力持續(xù)增加。雙拱形孔洞體系峰值應(yīng)力與應(yīng)力集中系數(shù)的關(guān)系見(jiàn)圖13。從圖13 可見(jiàn)：隨著壓應(yīng)力集中系數(shù)增加，峰值應(yīng)力呈單調(diào)降低的趨勢(shì)；紅色虛線為峰值應(yīng)力與壓應(yīng)力集中系數(shù)線性擬合曲線，較高的確定系數(shù)(0.907)表明可用線性函數(shù)表征兩者的聯(lián)系；而拉應(yīng)力集中系數(shù)，與峰值應(yīng)力無(wú)明顯關(guān)系，這說(shuō)明含雙孔洞巖體的破壞主要受壓應(yīng)力集中系數(shù)控制。

圖12 雙孔洞模型峰值應(yīng)力隨θ的變化Fig.12 Variation in the peak stress of numerical model containing two holes against θ

圖13 雙拱形孔洞體系峰值應(yīng)力與應(yīng)力集中系數(shù)的關(guān)系Fig.13 Relationship between the peak stress and stress concentration factor in the twin arched hole system

PFC模擬中含拱形孔洞圍巖破壞模式見(jiàn)圖14。圖14 中，紅色點(diǎn)代表裂紋，藍(lán)色圓圈表示巖石起裂位置。含單孔洞圍巖裂紋萌生于最大拉應(yīng)力集中的底板中心位置(圖14(a))，主裂紋從孔洞頂?shù)撞恐行难刂Q直方向延伸。對(duì)于雙孔洞圍巖，破裂主要發(fā)生于2個(gè)孔洞連接處，起裂位置、破壞程度均與連接角有關(guān)；當(dāng)θ為0°或90°時(shí)，2個(gè)孔洞間的裂紋方向與加載方向平行，其貫通模式屬于拉伸破壞。然而，當(dāng)2個(gè)孔洞傾斜布置時(shí)，其破壞屬于剪切破壞或拉/剪混合破壞。值得注意的是，當(dāng)θ為30°時(shí)，初始破壞出現(xiàn)在理論計(jì)算中孔洞1 邊界最大壓應(yīng)力所在位置附近，但該區(qū)域的破壞并未隨著初始破壞的出現(xiàn)而繼續(xù)發(fā)展，而是在2個(gè)孔洞的頂板或底板的拉伸裂紋發(fā)展到一定程度后才出現(xiàn)剪切貫通破壞。而對(duì)于其他角度，起裂位置均發(fā)生在孔洞頂?shù)撞?即拉應(yīng)力集中處)，但與中心位置略有偏移，這是第2孔洞的存在對(duì)應(yīng)力分布造成的影響(見(jiàn)圖8)。

圖14 孔洞周圍區(qū)域破壞模式Fig.14 Failure patterns around the hole

3 結(jié)論

1) 利用最優(yōu)化技術(shù)進(jìn)行不同平面和坐標(biāo)系下孔洞邊界坐標(biāo)點(diǎn)的變換，直接求解孔洞邊界的多余面力分布函數(shù)，實(shí)現(xiàn)了地下2個(gè)任意形狀孔洞在任意相對(duì)位置下的應(yīng)力求解。該方法無(wú)需進(jìn)行逆映射函數(shù)的繁瑣求解，避免了逆映射函數(shù)求解造成的誤差，適用于任意雙孔洞形狀和布置情況，降低了計(jì)算難度并提高了求解精度。

2) 鄰近孔洞的存在會(huì)改變現(xiàn)有孔洞邊界切向應(yīng)力分布，加劇應(yīng)力集中程度，其主要影響區(qū)域?yàn)榭拷噜徔锥吹囊粋?cè)。對(duì)于拱形孔洞，在單軸壓縮條件下，當(dāng)連接角θ到達(dá)30°時(shí)，孔洞邊界切向應(yīng)力的壓應(yīng)力集中程度達(dá)到最大；而θ為45°時(shí)，其拉應(yīng)力集中程度達(dá)到最大。

3) 理論計(jì)算與數(shù)值模擬結(jié)果表明，在單軸壓縮條件下，孔洞體系的初始萌生裂紋為孔洞頂?shù)装謇瓚?yīng)力集中區(qū)域形成的拉伸裂紋，但最終破壞表現(xiàn)為雙孔洞間的貫通破壞。不同連接角的雙孔洞體系峰值應(yīng)力與壓應(yīng)力集中系數(shù)呈線性相關(guān)關(guān)系，說(shuō)明雙孔洞巖石的破壞取決于孔洞邊界的壓應(yīng)力集中程度。