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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

式中：σ1為ζ1平面單位圓周上點的坐標(biāo)；t1為Z1坐標(biāo)系下孔洞1邊界上點的坐標(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)可以通過Laurent 級數(shù)表示。常數(shù)B、B′與C′反映了遠(yuǎn)場圍巖應(yīng)力的情況。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

以拱形孔洞為研究對象(圖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)原點為孔洞的形心。

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

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

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

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

由式(8)可得

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

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

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

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

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

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

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

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

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

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

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

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

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

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

在無側(cè)壓情況下，單拱形孔邊界以及雙拱形孔體系中孔洞1邊界環(huán)向應(yīng)力分布隨連接角θ的變化見圖8。圖8 中，位于孔洞外部的應(yīng)力為正，即為壓應(yīng)力，反之則為拉應(yīng)力；紅點和綠點分別表示最大壓力、拉應(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)可見：當(dāng)連接角θ為90°時，孔洞1邊界的切向應(yīng)力呈對稱布置；但在其他角度下，孔洞2的存在影響孔洞1周邊應(yīng)力的對稱性；當(dāng)連接角為0°～60°時，孔洞1環(huán)向最大壓應(yīng)力在孔洞左側(cè)拱底位置，而當(dāng)連接角為75°時，孔洞1 左側(cè)拱底環(huán)向壓力顯著降低，最大壓應(yīng)力出現(xiàn)在孔洞右側(cè)；隨著連接角的變化，孔周最大拉伸應(yīng)力位置會在孔洞頂、底中心位置附近變動，但無明顯變化規(guī)律；當(dāng)連接角為0°、45°、60°、75°時，最大拉應(yīng)力在孔底中心位置附近，而在其他角度下，最大拉應(yīng)力在孔頂中心位置。此外，從圖8還可見：當(dāng)2個孔洞沿斜線布置時，靠近相鄰孔洞一側(cè)與遠(yuǎn)離相鄰孔洞一側(cè)的邊界應(yīng)力呈現(xiàn)明顯的差異；在靠近孔洞2 的側(cè)邊(如左側(cè)邊墻和孔底)，環(huán)向應(yīng)力隨連接角變化明顯，而在遠(yuǎn)離孔洞2 的側(cè)邊(右側(cè)邊墻和孔頂)區(qū)域，應(yīng)力變化幅度相對較小，說明相鄰2個孔洞之間的相互影響主要集中在相鄰一側(cè)的區(qū)域。

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

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

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

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

離散元數(shù)值仿真可通過微觀尺度研究孔洞周邊的微裂紋的萌生和演化，從而反映出孔周應(yīng)力場的分布情況。為研究連接角對雙拱孔洞體系承載特性與破裂模式的影響，本文采用離散元軟件PFC2D開展含雙拱形孔洞體系巖體數(shù)值模擬?？锥葱螤睢⒊叽缂捌洳季志c理論求解中的模型一致，模型長×寬為100 m×100 m，側(cè)壓系數(shù)為0(如圖10所示)，豎直方向以恒定速率加載直至模型完全被破壞。采用PFC2D模擬時，首先以某砂巖為參照對象，基于該砂巖標(biāo)準(zhǔn)單軸壓縮試驗結(jié)果對數(shù)值模型微觀參數(shù)進(jìn)行標(biāo)定，各參數(shù)取值如表1所示。圖11 所示為砂巖室內(nèi)試驗單軸壓縮試驗與數(shù)值模擬結(jié)果。從圖11可見：數(shù)值模型除了無壓密階段外，其強度和彈性模量都與室內(nèi)試驗結(jié)果一致；同時，在破壞模式上，兩者都屬于斜剪切破壞，這說明可利用表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ī)律見圖12。從圖12可見：隨著θ增加，試樣的峰值應(yīng)力在15°時增加，隨后降低；當(dāng)θ超過30°時，試樣的峰值應(yīng)力持續(xù)增加。雙拱形孔洞體系峰值應(yīng)力與應(yīng)力集中系數(shù)的關(guān)系見圖13。從圖13 可見：隨著壓應(yīng)力集中系數(shù)增加，峰值應(yīng)力呈單調(diào)降低的趨勢；紅色虛線為峰值應(yīng)力與壓應(yīng)力集中系數(shù)線性擬合曲線，較高的確定系數(shù)(0.907)表明可用線性函數(shù)表征兩者的聯(lián)系；而拉應(yīng)力集中系數(shù)，與峰值應(yīng)力無明顯關(guān)系，這說明含雙孔洞巖體的破壞主要受壓應(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模擬中含拱形孔洞圍巖破壞模式見圖14。圖14 中，紅色點代表裂紋，藍(lán)色圓圈表示巖石起裂位置。含單孔洞圍巖裂紋萌生于最大拉應(yīng)力集中的底板中心位置(圖14(a))，主裂紋從孔洞頂?shù)撞恐行难刂Q直方向延伸。對于雙孔洞圍巖，破裂主要發(fā)生于2個孔洞連接處，起裂位置、破壞程度均與連接角有關(guān)；當(dāng)θ為0°或90°時，2個孔洞間的裂紋方向與加載方向平行，其貫通模式屬于拉伸破壞。然而，當(dāng)2個孔洞傾斜布置時，其破壞屬于剪切破壞或拉/剪混合破壞。值得注意的是，當(dāng)θ為30°時，初始破壞出現(xiàn)在理論計算中孔洞1 邊界最大壓應(yīng)力所在位置附近，但該區(qū)域的破壞并未隨著初始破壞的出現(xiàn)而繼續(xù)發(fā)展，而是在2個孔洞的頂板或底板的拉伸裂紋發(fā)展到一定程度后才出現(xiàn)剪切貫通破壞。而對于其他角度，起裂位置均發(fā)生在孔洞頂?shù)撞?即拉應(yīng)力集中處)，但與中心位置略有偏移，這是第2孔洞的存在對應(yīng)力分布造成的影響(見圖8)。

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

3 結(jié)論

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

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

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