水工圓形隧洞圍巖襯砌摩擦滑動接觸的新解法1)

尹崇林 呂愛鐘

(華北電力大學水電與巖土工程研究所，北京 102206)

引言

隧洞是廣泛應用于水電、交通、礦山及軍事工程領域的常用地下結構.近來，隨著我國經濟的持續發展及高新技術的不斷應用，我國的隧洞工程得到了前所未有的迅速發展，支護-圍巖相互作用，以及如何保證隧洞圍巖的穩定性是我們要研究的基本問題[1-2].

利用解析法求解隧洞工程問題是一種基本的分析方法.而由Muskhelishvili[3]提出的平面彈性復變函數方法是求解隧洞圍巖和襯砌中應力以及位移的常用解析法，通過這種方法容易獲得無支護隧洞的應力和位移解[4-5].而在實際工程中，常在隧洞內設置襯砌支護以保證隧洞的安全，因此研究支護-圍巖的相互作用更是十分重要，何川等[6]以及吳順川等[7]分別利用D-P準則和莫爾庫倫屈服準則，獲得了圍巖和襯砌相互作用的彈塑性解析解.

Hasanpour 等[8]、Ramoni 等[9]以及Son[10]都對支護和巖土材料之間的相互作用問題進行了研究.在求解有支護的深埋隧洞問題的解時，可以將其簡化為平面應變無限域問題，于學馥等[11]、Wang 等[12]以及Lu 等[13-14]已經獲得了各向同性彈性巖體中圓形隧洞襯砌和圍巖相互作用的解析解.Lu等[15-16]考慮了更復雜的孔形，以直墻半圓拱形隧洞為例獲得了非圓形隧洞的解析解.Bobet[17]和王少杰等[18]考慮了更復雜的巖體材料，分別視巖體為橫觀各向同性和正交各向異性的彈性體，獲得了非圓形隧洞的應力和位移的解析解.

隧洞在施加襯砌支護之后，圍巖和襯砌接觸并相互作用，在以往的解析研究中，大都考慮了完全接觸和光滑接觸兩種理想的接觸方式.文獻[12-13,15,17-18]都假定圍巖與襯砌之間的接觸為完全接觸，即認為接觸界面上法向的徑向應力和徑向位移及環向的剪切應力和切向位移都是連續的.這種接觸假定圍巖與襯砌之間非常粗糙，接觸面可以承受很大的摩擦力，不允許圍巖與襯砌之間產生相對滑動；而文獻[11,14,16]在求解過程中將圍巖與襯砌之間的接觸視為光滑接觸，即認為接觸界面上法向的徑向應力和徑向位移仍然是連續的，即圍巖和襯砌在法線方向不能相互脫開，但假定接觸面充分光滑，在接觸面的切線方向不能承受摩擦力，即認為環向方向的剪切應力為零，這種接觸允許圍巖與襯砌之間產生相對滑動.Atkinson 和Eftaxiopoulos[19]在求解套管井和膠結井水力壓裂問題的解析解時，在水泥和巖石的接觸面上考慮了這種光滑接觸條件.高永濤等[20]考慮雙層厚壁圓筒之間的光滑接觸并獲得了非均布荷載作用下內外壁的應力解析解.儲昭飛等[21]在求解非靜水應力場中圓形隧道黏彈性解析解時，其中也考慮了襯砌和圍巖之間的光滑接觸.

完全接觸和光滑接觸是人為假定的兩種理想情形，實際上兩種材料的接觸很可能并非為這兩種理想情形，例如，在建立連接結構接觸界面的非線性力學模型時，王東等[22]就將名義的光滑平面視作凹凸不平的粗糙面，考慮了微凸體的黏滑摩擦行為.圍巖和襯砌接觸時，它們之間并非完全光滑，也并非可以承受任意大的摩擦力，當圍巖與襯砌之間的剪應力大于所能承受的最大靜摩擦力后，將發生切向滑動.庫侖摩擦模型作為剛體的摩擦模型，因其簡單和實用而被廣泛應用于工程分析中.該模型認為切向摩擦力的數值不能超過法向壓力和摩擦系數的乘積.Meschke 等[23]建立計算模型以模擬樁土相互作用時隧道的機械化開挖，其中就利用了庫侖摩擦模型來模擬樁-土的相互作用.Cavalieri 等[24-25]利用增廣拉格朗日法求解涉及摩擦的接觸問題時，也利用庫侖摩擦定律建立了計算模型.蘇宗賢等[26]和楊釗等[27]利用庫侖摩擦模型分別模擬了管片和圍巖以及復合襯砌中內外層襯砌之間的接觸.呂愛鐘等[28]在庫侖摩擦模型的基礎上利用復變函數方法和最優化方法獲得了圓形隧洞圍巖和襯砌摩擦滑動接觸的解析解，并與有限元軟件ANSYS 所得結果進行了對比，驗證了其方法合理性.

本文同樣將通過庫侖摩擦模型模擬接觸面上的切向力學行為，但不同于文獻[28]，本文在優化過程中減少了設計變量的個數，極大地簡化了優化模型.在彈性接觸分析中，庫侖摩擦模型可表述如下

式中，τrθ接觸面上的剪應力，σr為接觸面法向壓力，fr為接觸界面綜合摩擦系數.該表達實質上是靜摩擦力的屈服條件或接觸面的滑移條件，將切向力與徑向壓力聯系起來.這種接觸條件最符合實際情況，即為摩擦滑動接觸.本文基于平面應變假定，將圍巖和襯砌視為各向同性線彈性體，考慮支護滯后效應，利用復變函數方法推導出深埋圓形隧洞在原始地應力和襯砌內部靜水壓力共同作用下圍巖與襯砌在摩擦滑動接觸情況下的解析解，為隧洞圍巖和襯砌接觸的支護設計計算提供理論基礎.

1 基本原理及方法

如圖1所示，為在原始地應力和靜水壓力共同作用下圓形襯砌支護隧洞.在開挖之前，圍巖無窮遠處水平和豎直方向的原始地應力分別為：，λ 為側壓力系數，定義拉應力為負，壓應力為正.隧洞襯砌的內外半徑分別為R0和R1，定義徑向和x方向的夾角為θ，以逆時針為正.

圖1 原始地應力和靜水壓力共同作用下深埋圓形襯砌隧洞Fig.1 Lined circular tunnel under in situ stresses and uniform hydrostatic pressure

1.1 摩擦滑動接觸的準則及數學模型

當圍巖和襯砌間的接觸為摩擦滑動接觸時，由于圍巖和襯砌之間的摩擦作用，接觸面可以傳遞剪力，但允許圍巖與襯砌之間的切向相對滑動，本文將通過最優化方法來處理這樣的問題.當圍巖與襯砌之間的剪應力大于所能承受的最大靜摩擦力后，圍巖和襯砌間接觸面間將產生相對滑動，接觸面兩側，即圍巖一側和襯砌一側的切向位移將產生間斷，即在r=R1的接觸面上，會存在uθ1uθ2的區域.本文作者提出在保證圍巖與襯砌之間剪應力小于或者等于接觸面所能提供的摩擦力前提下，將接觸面間產生最小滑移量的狀態作為襯砌的真實工作狀態.這樣的準則可以用最優化的數學模型表示為

其中，F(X)=|uθ1?uθ2|，bj(X)0 為庫侖摩擦模型表示的不等式組.

1.2 僅開挖作用時引起的圍巖位移

如圖2所示在隧洞開挖之后，假設沒有支護，圍巖會完成相應的位移.此時圍巖內任一點的徑向位移ur和切向位移uθ可以表示為[13]

圖2 開挖引起的圍巖位移Fig.2 Surrounding rock displacement caused by excavation

式中，R1為無支護隧洞的半徑，z是一個復數且z=R1eiθ，κ1=3-4μ1，G1=E1/[2(1+μ1)]，μ1，E1和G1分別是巖石的泊松比、彈性模量和剪切模量.式(3)即隧洞開挖后無支護時圍巖的全部位移.

1.3 求解圍巖和襯砌解析函數的基本方程

對于圓形隧洞，其應力分量和位移分量極易求得[3]

其中，?0(z)和ψ0(z)是關于圍巖或者襯砌的兩個解析函數.若能求出滿足邊界條件的圍巖和襯砌的兩個解析函數，則應力分量和位移分量均可以利用上述結果得出.

隧洞開挖后進行支護，圍巖和襯砌相互作用.圍巖對應的兩個解析函數為[13]

同時，襯砌對應的兩個解析函數可以用洛朗級數表示為

式(9)和式(10)中的系數ck，dk，ek，fk，gk和hk都是待求的實常數.

通過襯砌內邊界L1的應力邊界條件和圍巖襯砌接觸面L2的接觸條件可以建立求解解析函數系數的基本方程.

用σr2，σθ2，τrθ2分別表示襯砌的徑向應力、環向應力和剪應力.在襯砌內邊界L1上，作用有大小為p0的靜水壓力，則其應力邊界條件可以表示為(式中z=R0eiθ)

用σr1，σθ1，τrθ1分別表示圍巖的徑向應力，環向應力和剪應力.在圍巖襯砌接觸面L2上，圍巖和襯砌的徑向應力和剪應力連續，其應力連續條件可以表示為

本文考慮支護滯后的隧洞開挖過程，η 為位移釋放系數，如果η=0，即隧洞在開挖后立即進行支護，圍巖在支護之前沒有發生位移.明顯這不符合工程實際，即使在開挖后立即安裝支護，也不可避免地在圍巖中發生了一部分變形.假設隧洞位移在完成了最大位移ur+iuθ的η(0η1)倍后，再進行支護，此時圍巖完成的位移為η(ur+iuθ).設置襯砌以后，圍巖和襯砌相互作用，支護限制了部分圍巖位移的產生，此時圍巖的位移表示為ur1+iuθ1.而在圍巖作用下，襯砌的位移可以表示為ur2+iuθ2.

在圍巖襯砌接觸面L2上，根據法向位移連續條件可得，其中z=R1eiθ

式中κ2=3-4μ2，G2=E2/[2(1+μ2)].μ2，E2和G2分別是襯砌的泊松比、彈性模量和剪切模量.

將相應的z值代入式(9)～式(12)，式(15)～式(17)，利用冪級數解法，比較θ的同冪次項系數，可得到關于系數ck，dk+2，ek，fk+2，gk+2和hk的方程，其中k1.聯立所有方程，整理后發現當k2時，ck，dk+2，ek，fk+2，gk+2和hk都為0，只有c1，d1，d3，e1，f1，f3，g1，g3和h1為非零實數，可得方程(18)～式(25)

這樣利用邊界條件就得到了式(18)～式(25)共8 個線性方程.8 個方程，9 個未知量c1，d1，d3，e1，f1，f3，g1，g3，h1，只通過這8 個方程，是不可能求解9 個未知量的，如何求解這樣的問題是本文的關鍵所在.

本文將用圍巖一側和襯砌一側的切向位移間斷值uθ1?uθ2的大小來衡量接觸面間的滑移量，在式(2)中用目標函數F(X)=|uθ1?uθ2| 表示.由式(8)可得(式中z=R1eiθ)

由式(28)和式(29)可得

由式(30)可以看出，對任意θ，當|{}|內的值達到最小值時，所產生的位移間斷值都最小，只是對于不同θ，位移間斷值的最小值不同而已.所以取目標函數F為

使式(31)達到最小的一組c1，d1，d3，e1，f1，f3，g1，g3和h1，可以保證接觸面間產生最小的滑移量.在接觸面，必須滿足庫侖摩擦模型，式中τrθ和σr分別為接觸面上的剪應力和徑向正應力，由于在接觸面上應力的連續性，因此τrθ，σr既可取τrθ1，σr1，也可取τrθ2，σr2.本文取前者，則有

對任意的θ，式(32)都應該成立.由于問題的對稱性，只考慮θ∈[0,π]的情形即可.將[0,π]分為m等分(本文取180)，則式(32)可化歸為m+1 個約束條件

2 摩擦滑動接觸的解法

2.1 最優化問題的求解

文獻[28]使用混合罰函數方法來進行優化計算[29-30]，其數學模型為

按照文獻[28]的求解過程，其將式(18)～式(25)組成的8 個方程作為優化模型中的等式約束函數aj(X)，而將式(33)組成的不等式組作為不等式約束條件，然后將等式約束函數aj(X)和不等式約束函數bj(X)分別構成懲罰項加到目標函數F(X)構成一個新的無約束的稱為罰函數的目標函數p(x,r)，其中r稱為罰因子.這樣，就把原目標函數F(X)的約束極小值問題轉為求罰函數p(x,r)的無約束極小值問題.其罰函數的具體形式為

式中下標集合I1，I2定義為

其中X(0)是給出的原問題的計算初始點.在該模型中，X=[c1,d1,d3,e1,f1,f3,g1,g3,h1]，即所有的變量都需要參與到優化的過程中，計算初始點X(0)=我們知道，設計變量的增加會優化的數學模型更加復雜，從而增加優化的難度，有時甚至得不到滿足約束的最優解.如果可以減少設計變量的個數，這將大大地簡化優化的數學模型，從而使求解變得容易.

本文同樣使用混合罰函數法來進行優化計算，不同的是本文將設計變量的個數從9 個變為1 個，極大地簡化了優化模型.具體來說，本文可以將c1，d1，d3，e1，f1，f3，g1，g3和h1其中的任意一個作為設計變量，以h1為例，賦予h1初值，根據式(18)～式(25)組成方程組即可以求解出c1，d1，d3，e1，f1，f3，g1，g3這8 個變量，再將c1，d1，d3，e1，f1，f3，g1，g3和h1代入如下的優化模型中，即可獲得優化結果

式中，bj(X)0 表示的是不等式約束條件，它由式(33)構成，而F(X)即為式(31).在該優化模型中，其罰函數的具體形式簡化為

可以看出，較之文獻[28]中的罰函數的具體形式，本文的罰函數形式沒有二次損失項中的并且計算初始點X(0)僅有，這極大地簡化了優化的數學模型.模型中只有h1為設計變量，改變h1的值，使式(36)達到最小值的即為該優化問題的解，此時其他8個變量可以由解方程組得出且分別表示為為了比較兩種模型的精確程度，我們以摩擦滑動接觸為例，利用新舊兩種模型分別計算出不同摩擦系數fr對應的目標函數的值，其對比如表1所示.

表1 改變摩擦系數時兩種優化模型所得目標函數值比較Table 1 Objective function values obtained by two optimization models when changing friction coefficient fr

從表1 可以看出，比較不同摩擦系數fr計算出的目標函數值，本文提出的新優化模型所得結果更小，這是因為本文通過解方程的方法可以使式(34)中的aj(X)=0 精確滿足，較之文獻[28]中的優化模型，該方法更加精確.

同時可以通過訪問目標函數的次數來對比兩種方法的迭代速度(見表2)，可以看出新優化模型較之舊模型有很大地提升.

表2 兩種優化模型訪問目標函數的次數Table 2 Times to access F(X)for two optimization models

利用新優化模型求解出的這樣一組X=[c1,d1,d3,e1,f1,f3,g1,g3,h1]即為本文問題的解，從而完成了圍巖支護后，圍巖和襯砌中復勢函數的求解.從以上的分析和求解過程可以看出，這種摩擦滑動接觸的解法具有一般性，它包含了完全接觸和光滑接觸兩種極限情形.

當fr=0 時，由式(32)表示的不等式約束條件將變為

若令τrθ1=0，則可以獲得與式(38)完全相同的方程，即式(38)表示的就是光滑接觸的情形.

當fr的值大于某個值時，式(32)表示的不等式約束條件總會滿足，式(36)達到最小值(其值為零)時的解為

則式(39)所表示的實際就是uθ1=uθ2，這由式(26)、式(27)可以清晰地看出.由后面的算例可知，滿足完全接觸條件所需要的fr值大小與側壓力系數λ 密切相關，λ 越接近于1，所需要的fr值越小，理論上當λ=1 時，fr=0 就可以滿足完全接觸條件.

2.2 襯砌和圍巖中的應力

給定R0，R1，p，λ，p0，G1，μ1，G2，μ2，η，可以求出待定的c1，d1，d3，e1，f1，f3，g1，g3，h1.

襯砌中的應力可表示為[13]

圍巖中的應力應由開挖前，開挖后，支護后三部分的應力進行迭加.迭加后的應力仍用σr1，σθ1，τrθ1符號表示，可求得[13]

3 分析和討論

3.1 圍巖和襯砌接觸面上的接觸方式

取計算參數為：p=10.0 MPa，p0=2.0 MPa，R1=3.0 m，R0=2.7 m，μ1=0.25，μ2=0.20，E2/E1=10.0，η=0.20.對于不同的側壓力系數λ以及摩擦系數fr，獲得的解析函數的系數也不同.由于問題的對稱性，只需要分析θ∈[0?,90?]內的結果.

由前面的推導我們知道，本文作者提出的摩擦滑動接觸可以計算光滑接觸和完全接觸兩種極限的接觸工況，并且這與摩擦系數fr的取值相關.取側壓力系數λ=0.5，分析摩擦系數fr取不同值時接觸面上接觸應力σr，τrθ和切向位移間斷值|uθ1?uθ2|的變化規律，并將其與完全接觸[13]和光滑接觸[14]的結果進行對比.

從圖3 可以看出，接觸面上一旦發生滑動，除了對稱點0?和90?，其余各點都會發生相對滑動，且在45?處有最大值.當摩擦系數fr為零時，接觸上發生最大的相對滑動，當fr0.64 時，接觸面上的切向位移間斷值|uθ1?uθ2|基本為零，即uθ1≈uθ2，并且隨著fr的增大，|uθ1?uθ2|逐漸減小.

圖3 接觸面上切向位移間斷值|uθ1?uθ2|的分布Fig.3 Distributions of|uθ1?uθ2|on the interface

如圖4(a)所示，接觸面上徑向應力σr在45?處為一定值，而當fr=0 即光滑接觸時，σr的絕對值分別在0?和90?處取得極小值和極大值.當fr較大時，σr的絕對值在[0?,90?]范圍內單調遞減.圖4(b)表明隨著fr的增大，接觸面上的剪應力τrθ的絕對值也隨之增大，且都在45?處取得最大值.當fr0.64時，圖4 中的σr，τrθ不再變化，且與完全接觸的結果一致，結合圖3，fr0.64 時，uθ1≈uθ2，此時得到的解即為完全接觸的解.

圖4 接觸面上接觸應力隨摩擦系數fr 變化的分布規律Fig.4 Distributions of the contact stresses on the interface for different fr

我們已經確定利用摩擦滑動接觸解法可以求解完全接觸和光滑接觸兩種極限接觸問題，而且通過推導已知當fr的值大于某個值時，總能滿足完全接觸的條件，稱這個值為閾值確定閾值對判斷接觸面上的接觸方式很有幫助.文獻[28]通過反復試算的方式討論了不同彈性模量和側壓力系數所確定的閾值，但是這種方法不夠精確，本文嘗試利用一種精確的方法來確定閾值

表3 改變位移釋放系數η 時不同側壓力系數λ所確定的閾值Table 3 Values of for various η and λ

表3 改變位移釋放系數η 時不同側壓力系數λ所確定的閾值Table 3 Values of for various η and λ

表4 改變襯砌厚度R1?R0時不同側壓力系數λ所確定的閾值Table 4 Values of for various R1?R0and λ

表4 改變襯砌厚度R1?R0時不同側壓力系數λ所確定的閾值Table 4 Values of for various R1?R0and λ

表5 改變水壓力p0時不同側壓力系數λ 所確定的閾值Table 5 Values of for various p0and λ

表5 改變水壓力p0時不同側壓力系數λ 所確定的閾值Table 5 Values of for various p0and λ

同時，為了同文獻[28]中所得結果進行比較，取R0=2.5 m，側壓力系數分別為0.2，0.5，0.8，1.0，在不同彈模比值下所得閾值如表6 所示.

表6 改變彈模比值E2/E1時不同側壓力系數λ所確定的閾值Table 6 Values of for various E2/E1and λ

表6 改變彈模比值E2/E1時不同側壓力系數λ所確定的閾值Table 6 Values of for various E2/E1and λ

3.2 襯砌和圍巖各邊界上切向應力的變化規律

在隧洞的支護設計中，各邊界上的切向應力變化規律十分重要，所以在這一節，我們取不同的摩擦系數fr，分別分析它們的變化對襯砌內外邊界的切向應力以及圍巖開挖邊界上切向應力(σθ1)的影響.3.1 節已經分析過當λ=0.5，滿足完全接觸條件的閾值=0.64，因此，在本節的討論中，摩擦系數fr的最大取值為0.64，其余參數不變化.

由圖5 可以看出，當λ=0.5，不管摩擦系數如何改變，所有的切向應力都為壓應力.由圖5(a)和圖5(b)可以看出，隨著fr的增大，襯砌內邊界的切向應力最大壓應力增大，最小壓應力減小，在[0?,90?]內，壓應力隨著角度的增大而減小，其變化范圍增大.在襯砌外邊界上，隨著fr的增大，切向應力變化范圍減小.只是當fr較小且趨于0時，從0?到90?，壓應力隨著角度的增大而增大，分別在0?和90?取得最小和最大壓應力；當fr較大時，結論恰好相反.由圖5(c)可得，保持λ=0.5 不變，在[0?,90?]內，隨著角度的增大，σθ1為壓應力單調減小.隨著fr的增大，σθ1在90?出現的最小壓應力增大，在0?出現的最大壓應力減小.圖5 表明摩擦滑動接觸條件下所得襯砌內外邊界上的切向應力以及圍巖開挖邊界上的切向應力σθ1都介于光滑接觸和完全接觸的結果之間.

圖5 各邊界上切向應力隨摩擦系數fr 變化的分布規律Fig.5 Distributions of the tangential stresses on the 3 boders for different fr

4 結論

本文提出了更加符合實際情況的摩擦滑動接觸條件來描述實際工程中圍巖和襯砌之間的接觸問題.以庫侖摩擦模型模擬圍巖襯砌之間的接觸，當接觸面上發生相對滑動之后，將接觸面上產生最小滑移量的狀態視為襯砌的真實工作狀態，并以此為基礎，利用最優化方法建立了摩擦滑動接觸情形下新的解法.在利用混合罰函數法求解解析函數的系數時，減少了設計變量的個數，去掉了罰函數二次損失項中的極大地簡化了最優化解法的數學模型.

同時提出了能夠精確獲得滿足完全接觸條件的摩擦系數閾值的公式，并且得到了不同位移釋放系數，不同襯砌厚度，不同水壓力大小以及不同彈模比值情況下改變側壓力系數時所獲得的閾值為判斷圍巖和襯砌接觸方式提供了理論基礎.摩擦滑動接觸條件下所得的解都介于完全接觸和光滑接觸之間，并且能同時可以求解光滑接觸和完全接觸兩種極限情況，具有一般性.特別地，在襯砌內邊界上，摩擦滑動接觸所得最大切向應力介于兩種極限接觸之間，且光滑接觸所得最大切向應力最小，這說明在實際工程中，盡量減小襯砌和圍巖之間的摩擦，可以增加襯砌的承載能力.