楊釗 ,潘曉明 ,余俊
(1. 同濟大學 地下建筑與工程系,上海,200092;2. 同濟大學 巖土及地下工程教育部重點實驗室,上海,200092;3. 中南大學 土木建筑學院,湖南 長沙,410075)
在盾構隧道中,在大數(shù)多情況下,二次襯砌不作為隧道結構的主要受力構件,其主要目的是修正施工“蛇行”、運行減糙、提高外管片耐久性和防水等維護作用,盾構隧道復合襯砌的內力計算在設計中一直沒有受到重視。近年來,我國擬采用盾構法修建多條大型輸水隧洞和城市大型下水管道,二次襯砌在輸水隧洞中將作為主要受力構件出現(xiàn)。因此,研究盾構隧洞復合襯砌計算模型已成為一項重要的課題?,F(xiàn)有的盾構隧洞復合襯砌計算模型多基于國際隧道協(xié)會發(fā)布的盾構隧道襯砌設計指南[1]。在設計指南中提出了 2種盾構隧洞復合襯砌計算模型:雙層框架模型和彈性方程模型。張厚美等[2]改進了第 1種模型中內、外層襯砌接觸界面模型,提出了關于接觸界面的3種模型:抗壓縮模型、局部抗彎模型和剪壓模型。以上模型均是基于荷載-結構法,沒有考慮土體開挖、回填注漿等施工過程對襯砌結構受力的影響。章青等[3]基于非連續(xù)介質變形理論提出了一種新的模型,模型通過放松離散單元交界面上的位移連續(xù)條件,來模擬螺栓連接和各種縫面的不連續(xù)變形特征。該模型仿真性強,計算效率高。本文作者從理論上分析雙層框架模型存在的不完善之處,提出一種基于地層?結構法的計算模型即實體疊合計算模型。該模型不僅能更精確地考慮土體與管片的相互作用,而且能考慮隧洞施工的各個階段對結構受力的附加影響。
考慮到外管片與內襯不是同期施工,且施工期和運營期的受力狀況不同,將襯砌的受力過程分成3種工況。
工況 1:外管片施工期,管片自重、施工期荷載及全部水土壓力由外管片環(huán)單獨承擔。
工況 2:運營初期,此時外部水土壓力已經(jīng)達到穩(wěn)定,且全部由外管片環(huán)承擔。內襯和外管片共同承擔內水壓力和內襯自重作用。
工況 3:運行穩(wěn)定期,內襯和外管片共同承擔工況2后繼增大的外部水土壓力。
由于工況3不為設計控制工況,雙層框架模型在控制工況下的受力模型如圖1所示。在模型中,外管片的最終受力、變位狀態(tài)由工況1與工況2的受力、變位狀態(tài)疊加而得。由圖1可知,雙層框架模型存在1個假定:工況1管片的受力與變位狀態(tài)不影響工況2管片和內襯的受力與變位,即外管片已發(fā)生的變形與受力對管片與內襯二者相互作用無影響。這一假定只有在管片接頭剛度、管片、內襯混凝土、管片周圍土體抗力均為線彈性以及管片周圍地基彈簧的分布范圍不變的情況下才能成立。而在實際工程中,管片接頭的抗彎剛度是非線性的,在不同工況下管片周圍地基彈簧分布范圍是變化的,管片、內襯混凝土、管片周圍土體抗力只能近似地看著是線彈性的。
在雙層框架模型中,管片與內襯接觸界面的力學行為可通過在接觸界面處設置壓縮彈簧、壓縮剪切彈簧或剛性鏈桿來實現(xiàn)。
考慮 2根長度為l的Euler-Bernouli疊合梁[4?5],彈性模量為E,截面彎曲慣性矩為I。假定2根梁的接觸界面之間只能傳遞接觸壓力,而不能傳遞彎矩與剪力。2根梁的接觸界面之間分別采用壓縮彈簧和goodman接觸單元來模擬接觸界面的接觸力學行為。
把梁用1個二節(jié)點的梁單元模擬。采用壓縮彈簧模擬界面力學效應,設壓縮彈簧的剛度系數(shù)為K,則壓縮彈簧單元對整個體系的附加剛度矩陣為:
圖1 雙層框架受力疊加模型Fig.1 Superposition model of double linings under loading
在界面之間設置只能傳遞壓力、不能傳遞剪力的goodman接觸單元,接觸單元上緣與下緣的垂直位移差可以表示為[6]:
式中:
單元內正應力與垂直位移差成正比,即
假定單元各結點上產(chǎn)生虛位移δ*e,那么,單元內虛位移差為:
在單位長度上,單元應力所做虛功為:
沿單元長度積分后,即得到單元應力所做的虛功,等于單元結點力所做的虛功,由此可得接觸單元的剛度矩陣(即為接觸單元對整個體系的附加剛度矩陣)為:
2種附加剛度矩陣只有在梁單元轉角為0°或梁單元長度足夠小,且梁單元兩節(jié)點的位移相等時,兩者的剛度矩陣利用K彈簧=0.5λtl換算,得到的計算結果才是完全等同的。
同理可以得到壓剪彈簧也只有在接觸界面劃分足夠多單元的情況下,才能較為精確地模擬2根梁之間的真實接觸狀態(tài)。
對于復合襯砌結構,第1階段由外管片承受外周水、土壓力,第2階段由內襯與外管片共同承擔內水壓力、內水與內襯自重。因此,復合襯砌結構可看作2個階段受力的疊合梁結構。現(xiàn)考慮一階段受力疊合梁,疊合梁由2根截面屬性與材料屬性均一樣的Euler-Bernouli梁組成,梁的彈性模量為2.0 GPa,截面尺寸為1 m×1 m,梁的長度為20 m,梁的一端固定,另一端承受集中力荷載10 kN。假設2根梁之間的接觸界面通過設置錨筋,表面鑿毛等方式使得2根梁在界面之間可以傳遞彎矩與剪力并保持界面之間的位移協(xié)調。若采用雙層框架模型計算,則應在2梁之間設置剛性鏈桿來模擬這種接觸界面性質。有限元模型如圖2所示。在實際計算過程中,可以將通過上述界面處理的2根梁看作1根整體梁進行計算。
當采用剛性鏈桿模擬界面力學性質,可以實現(xiàn) 2根梁在接觸界面間的位移協(xié)調。采用雙層框架模型計算所得梁的撓度與采用整梁計算所得的撓度相近,雙層框架模型計算所得最大撓度值為1.59 cm;整體梁模型計算所得最大撓度值為1.50 cm。但2種模型計算所得的梁端最大彎矩與采用整梁計算所得的最大彎矩值相差甚遠。采用雙層框架模型計算得到梁截面的最大彎矩為4.11×105N·m;將梁截面作為一個整體,計算所得梁截面最大彎矩為1.98×106N·m。這主要是因為當采用雙層框架模型計算時,上梁與下梁均具有中性軸,位于上梁與下梁的截面形心處,而實際上中性軸只有1條位于整體梁的截面形心處,這與實際情況不符。
在 ABAQUS中,雙層框架模型不精確模擬一階段受力疊合梁,更不可能精確模擬受力更加復雜的二階段受力疊合梁。因此,當采用 ABAQUS數(shù)值模擬軟件計算時,復合襯砌計算模型應采用實體單元模型。
隧道襯砌結構沿隧道縱向可以認為是一個無限長的結構體,在不考慮縱向變形的條件下,管片環(huán)的內力不沿縱向發(fā)生變化,因此,在計算過程中,可以將每一環(huán)管片作為一個平面應變狀況考慮。
圖2 有限元模型Fig.2 Model of FEM
在實體疊合模型中采用地層?結構法計算管片、內襯的受力與變位,考慮到土體材料的彈塑性,選用Drucker-Prager屈服準則和相關聯(lián)流動法則來模擬土體材料的非線性應力?應變特性。
工程實踐與理論研究表明[7?8]:當隧道上覆土層為砂性土時,采用水土分算,水壓力直接作用于管片上,土體的容重采用浮容重計算,土體的側壓力系數(shù)取為k0。
2.3.1 模型簡介
當管片采用梁單元模擬時,其縱向接頭的力學性能可采用旋轉彈簧、壓縮彈簧和剪切彈簧來分別模擬接頭的抗彎、抗壓和抗剪的力學狀態(tài)[9]。這種“梁?彈簧模型或殼彈簧模型”[10?13]已經(jīng)被廣泛應用于工程設計與計算中,并經(jīng)過實測驗證。當采用試驗給定的上述諸接頭力學參數(shù)后,計算所得的管片截面內力和變位值均能與實測值基本吻合。
當管片采用實體單元模擬時,其接頭的力學性能現(xiàn)多采用“接觸單元”來模擬[14?15]。根據(jù)接頭的實際受力與變形情況,管片接頭處的接觸單元必須綜合考慮:彈性密封墊與混凝土管片接觸、管片與管片的接觸、螺栓與螺栓孔接觸以及螺栓的錨固作用,這將使接頭處的建模問題顯得十分復雜;同時,在國內外眾多的接頭實驗資料中,很少提及管片接頭相關的接觸參數(shù)。這使得管片的接頭模型不僅在技術上不易實現(xiàn),而且在精度上也難以準確保證。為此,本文作者提出了一種新的實體簡化接頭模型以用于采用實體單元建模的管片分析。本文建議的接頭模型既能避免接頭部位上的復雜情況,又能借用過去從常規(guī)實驗所得或直接采用接頭剛度理論計算的接頭力學參數(shù)。實體簡化接頭模型的示意圖如圖3所示。
圖3 實體簡化接頭模型示意圖Fig.3 Simplified model for solid element
在模型中將接頭端肋與剛性體相連,并設置剛性體的參考點位于管片截面形心處。在兩剛性參考點之間設置剪切彈簧、壓縮彈簧和抗彎彈簧來模擬接頭的力學性質,這 3根彈簧的剛度取用“梁?彈簧模型或殼彈簧模型”中對應的彈簧剛度。
2.3.2 模型驗證
分別采用實體簡化接頭模型與“梁?接頭模型”計算如圖4所示模型。AB邊上受到400 kN/m的均布荷載,BC邊上受到300 kN/m的均布荷載。在AD和DC邊分別施加沿X和Y軸的對稱邊界條件。洞室外直徑為7.0 m,管片厚度為40 cm。洞室由4片管片組成,管片接頭分別位于 0°,90°,180°和 270°。0°角和90°角分別位于X軸正方向和Y軸正方向。
圖4 模型示意圖Fig.4 Simplified plot of model
分別采用實體單元與梁單元模擬管片,將兩者計算所得的管片內緣環(huán)向應力對比,如圖5所示。采用接頭模型計算所得的管片內緣環(huán)向應力相對誤差小于5%,實體簡化接頭模型能有效地應用于工程實際計算。
圖5 管片內緣環(huán)向應力對比Fig.5 Comparisons of circumferential stress between beam element and solid element
地層與管片之間采用無厚度的接觸單元模擬。接觸面的徑向力學行為采用“硬接觸”模擬,即接觸面之間可以傳遞無窮大的徑向壓應力,但不能傳遞徑向拉應力(在徑向拉力的作用下接觸面自動脫開)。接觸面的切向力學行為采用庫侖摩擦模型模擬[6]。庫侖摩擦模型用于判斷接觸面時是否發(fā)生相對錯動,同時也用于分析相對滑動對管片襯砌應力場的影響。
庫侖摩擦模型可以定義如下:
式中:τtrue為計算所得的真實剪應力;τlim為最大允許剪應力;μ為接觸界面綜合摩擦因數(shù);p為接觸面法向壓力。在接觸計算中,接觸面之間傳遞的最大剪應力不得超過接觸壓力乘以庫侖摩擦因數(shù)。當計算所得的剪應力大于接觸壓力乘以庫侖摩擦因數(shù)時,接觸面發(fā)生相對滑動,將接觸界面的狀態(tài)設為滑移進入下一步迭代計算。
復合襯砌接觸界面模型與內、外層襯砌間接觸界面的處理方式有關。根據(jù)接觸界面處理方式的不同,提出4種接觸界面模型。
2.5.1 無黏聚力模型
內襯施工前將管片內表面的螺栓手孔、注漿孔、起吊孔等凹槽用水泥充填抹平,然后澆注內襯混凝土。此時,內外層襯砌之間可以傳遞徑向壓力,且可以通過摩擦力的形式傳遞切向剪力。復合襯砌接觸界面模型采用與管土接觸相同的庫侖摩擦模型。
2.5.2 有黏聚力模型
澆筑內襯前抹平外管片內表面較大凹槽,其余部分通過鑿毛處理。此時,內、外層襯砌間可以傳遞壓力,且可以通過黏聚力與摩擦力的形式傳遞剪力。接觸面的徑向力學行為采用“硬接觸”模擬。接觸面的切向力學行為采用有黏聚力的庫侖摩擦模擬。
有摩擦的庫侖模型可以定義如下:
式中:τtrue為計算所得的真實剪應力;τlim為最大允許剪應力;μ為接觸界面綜合摩擦因數(shù);c為新老混凝土之間的黏聚力;p為接觸面法向壓力。
2.5.3 位移協(xié)調模型
在管片預制時設置錨筋孔,內襯澆筑時,在管片錨筋孔處設置錨筋,錨筋伸入到內襯混凝土內,并與內襯受力鋼筋點焊連接。此時,內、外層襯砌界面之間的位移可以認為是協(xié)調的,接觸界面模型采用位移協(xié)調模型。
2.5.4 局部位移協(xié)調模型
內襯混凝土或鋼筋伸入外管片手孔或凹槽中,這些部位的接觸界面的位移可以認為是協(xié)調的。其余部位應接觸界面的處理方式選用無黏聚力模型或有黏聚力模型。
為了考慮內襯施作之前外管片先已發(fā)生的變形與受力、對后續(xù)管片與內襯二者相互作用的附加影響,計算中應需計入外管片、內襯結構二者先后的不同受力作用。即在管片參與內襯結構共同受力時,要先計入已經(jīng)發(fā)生的管片內力與內凈變形的收斂值。
在數(shù)值模擬分析計算過程中,可以采用以下步驟來考慮并計入外管片和內襯的實際受力及其變形狀態(tài):(1) 建立有限元數(shù)值模型;(2) 平衡初始地應力;(3) 開挖土體單元,并施作外層管片;(4) 施加外水壓和管片自重;(5) 輸出按此時模型計算所得各單元的應力以及外層管片的變位值;(6) 輸入上面輸出的應力,用以平衡原先應力。此時的受力狀態(tài)為:土體已開挖,外層管片已經(jīng)施作,外周土水壓力作用在管片上。經(jīng)過應力平衡后,此時土體與外層管片的“計算變形”為0,應力狀態(tài)為第5步輸出的應力狀態(tài),即為二次襯砌施作之前的真實受力狀態(tài)。(7) 施作二次復合襯砌的內襯;(8) 施加內水壓和二襯自重;(9) 輸出此時的應力狀態(tài),即為真實的二襯復合襯砌受力狀態(tài)。輸出此時外層管片的變位值,并使之與第5步管片變位的輸出值相迭加,即為真實狀態(tài)下管片的變位值。
當物體處于小變形狀態(tài)時,其名義應變等于真實應變。因此,管片平衡應力后的應變(第9步輸出的應變)應即等于管片的真實應變減去管片平衡前的應變值(第5步輸出的應變),進而得出管片的真實變位值,其值可分解為平衡應力前的變位值(第 5步輸出的變位值)再加平衡后的變位值(第9步輸出的變位值)。因土體的塑性狀態(tài)是通過應力判斷,應力是通過增量應變計算,因此經(jīng)過第6步的平衡應力后,土體的應力狀態(tài)與塑性區(qū)與實際狀態(tài)相符。
青草沙原水過江隧洞工程位于上海長江隧道下游約80 m處,浦東側越江點在五號溝,長興島越江點在該島新開河附近,全長7.23 km。越江輸水管道設計擬采用 TBM 施工,采用有壓輸水,設計為圓形斷面,襯砌結構外直徑為6.8 m,管片厚為400 mm,內襯厚為 250 mm。隧道設計的控制截面上覆土為砂性土,覆土高度為 10.5 m,隧道中心線處內水壓為 0.309 MPa,隧道中心線處外水壓為0.300 MPa。內、外層襯砌間采用20 cm×20 cm分布的錨筋,故復合襯砌接觸界面模型采用位移協(xié)調模型計算。
有限元計算范圍為盾構隧道左右各取30 m,地表往下取45 m。對邊界,假設左右兩側存在水平約束,下部存在豎向約束,上部邊界自由。有限元計算模型如圖6所示??v向接頭的正、負抗彎、抗壓、抗剪剛度系數(shù)分別取為:kθ+=4.0×107N/(m·rad),kθ-=2.0×107N/(m·rad),Kn=5.0×1012N/m,KT=5.0×1011N/m。管片的壓縮模量E和泊松比μ取為:E=3.55×104MPa,μ=0.2;內襯混凝土的彈性參數(shù)取為:E=3.25×104MPa,μ=0.2。土層參數(shù)見表1。
圖6 有限元模型Fig.6 Model of finite element method
表1 土層參數(shù)Table 1 Parameter of soil layer
根據(jù)上述計算模型,采用 ABAQUS有限元軟件實現(xiàn)青草沙原水過江隧洞控制截面計算。采用實體單元計算時,可以通過預先定義截面得到管片與內襯的內力,但在后處理中無法圖形顯示內力。將內力與變位值輸出,編制后處理程序得到結構的內力圖與變位圖。
圖7所示為管片變位矢量。由圖7可知:內水壓施加以后,管片豎直方向壓縮,垂直方向伸長。管片的形狀由原先的圓形變成扁平的橢圓形。
圖8所示為管片截面內力。由圖8可見:彎矩以管片內側受拉為正。管片彎矩的峰值出現(xiàn)在管頂、管底和兩腰,其管頂、管底為正值,兩腰彎矩為負值。最大正彎矩位于管底處,其值為1.10×105N·m。最大負彎矩位于管腰處,其值為?0.87×105N·m。設軸力以受拉為正,受壓為負。管片全截面受壓。管頂、管底軸壓力小,而管腰軸壓力大。軸壓力最大值位于管腰處,最大為6.57×105N。
表2 管片控制截面變位Table 2 Displacement of segment in control section
圖7 管片變位矢量Fig.7 Displacement vector diagram of segment
圖8 管片截面內力Fig.8 Member force of segment section
內水壓施加以后,由于內襯質量與內水的質量,內襯主要發(fā)生向下的剛體位移,如圖9所示。
圖9 內襯變位矢量Fig.9 Displacement vector diagram of second linings
圖10所示為內襯截面內力。由圖10可見:彎矩以襯砌內側受拉為正。內襯彎矩的峰值出現(xiàn)在內襯頂部、內襯底部和兩腰,內襯頂部、底部為正,兩腰為負。最大正彎矩位于內部底部,為1.32×104N·m。最大負彎矩位于內襯腰處,其值為?0.97×104N·m。內襯全截面受拉。內襯頂部、底部軸拉力大,而腰部軸拉力小。軸拉力最大值位于內襯底部,最大值為?4.56×105N。
圖10 內襯截面內力Fig.10 Member force of seconding linings section
運營期內水壓作用,內、外襯砌聯(lián)合承擔內水壓力引起軸拉力。研究內、外襯砌對軸拉力的貢獻(見表3),對以后類似的工程有重大意義。當內水壓力作用時,內襯產(chǎn)生軸拉力,外襯軸力減小。因此,可以認為:內水壓力作用產(chǎn)生軸拉力=內襯軸拉力+外襯軸力減少量(其中,外襯軸力減小量=施工期外管片軸力?運營其管片軸力)。
內、外層襯砌軸力對比如表3所示。由表3可知:內水壓力引起的襯砌總拉力約40%由內襯承擔。
表3 內、外層襯砌軸力對比Table 3 Comparison of axial force on primary and seconding linings
(1) 雙層框架模型的接觸界面處理不能精確模擬內、外層襯砌間界面的力學行為,且不能考慮外管片先已發(fā)生的受力與變形對后續(xù)內、外襯砌之間相互作用的影響。
(2) 實體疊合模型內、外層襯砌接合面力學性態(tài)采用傳統(tǒng)的接觸單元模擬,理論體系清晰,能真實、準確地反映接觸界面的力學性態(tài)。
(3) 實體簡化接頭模型參數(shù)物理意義明確,參數(shù)可通過常規(guī)試驗或接頭剛度計算理論確定。該模型不僅簡化了當管片采用實體單元建模時接頭部位處理的復雜性,而且能有效地保證計算的精度。
(4) 實體疊合模型不僅可考慮隧洞施工的各個階段對結構受力的影響,而且能較好地反映外管片已發(fā)生的變形與受力、對管片與內襯二者相互作用的附加影響。實際工程應用結果表明:實體疊合模型對于分析盾構隧洞復合襯砌結構是可行且有效的。
[1] Takano Y H. Guidelines for the design of shield tunnel lining[J].Tunneling and Underground Space Technology, 2000, 15(3):303?331.
[2] 張厚美, 過遲, 呂國梁. 盾構壓力隧洞雙層襯砌的力學模型研究[J]. 水利學報, 2001(4): 28?34.ZHANG Hou-mei, GUO Chi, LU Guo-liang. Mechanical model for shield pressure tunnel with secondary linings[J]. Journal of Hydraulic Engineering, 2004(4): 28?34.
[3] 章青, 卓家壽. 盾構式輸水隧洞的計算模型及其工程應用[J].水利學報, 1999(2): 19?22.ZHANG Qing, ZHUO Jia-shou. A computational model of shield tunnel for water conveyance[J]. Journal of Hydraulic Engineering, 1999(2): 28?34.
[4] 馮紫良. 桿系結構的計算機分析[M]. 上海: 同濟大學出版社:1991: 50?79.FENG Zhi-liang. Trussing structure computer analysis[M].Shanghai: Tongji University Press, 1991: 50?79.
[5] 王勖成. 有限單元法[M]. 北京: 清華大學出版社, 2006:128?173.WANG Xu-cheng. The finite element method[M]. Beijing:Tsinghua University Press, 2006: 128?173.
[6] 朱伯芳. 有限單元法原理與應用[M]. 北京: 中國水利水電出版社, 2004: 72?156.ZHU Bo-feng. The finite element method theory and application[M]. Beijing: China Water Power Press, 2004:72?156.
[7] 鐘登華, 佟大威, 王帥, 等. 深埋 TBM 施工輸水隧洞結構的三維仿真分析[J]. 巖土力學, 2008, 29(3): 609?613.ZHONG Deng-hua, TONG Da-wei, WANG Shuai, et al. 3D simulation analysis of structure of deep water diversion tunnel with TBM constructing[J]. Rock and Soil Mechanics, 2008,29(3): 609?613.
[8] Msdhimo H, Ishimura T. Evaluation of the load on shield tunnel lining in gravel[J]. Tunneling and Underground Space Technology, 2003, 18(2): 233?241.
[9] ZHU He-hua, Hashimoto T. Back-analysis for shield tunnel using beam-joint model[R]. Shanghai: Tongji University.College of Civil Engineering, 1996.
[10] 朱合華, 丁文其, 李曉軍. 盾構隧道施工力學性態(tài)模擬及工程應用[J]. 土木工程學報, 2000, 33(3): 98?103.ZHU He-hua, DING Wen-qi, LI Xiao-jun. Construction simulation for the mechanical behavior of shield tunnel and its application[J]. China Civil Engineering Journal, 2000, 33(3):98?103.
[11] 朱偉, 黃正榮, 梁精華. 盾構襯砌管片的殼?彈簧設計模型研究[J]. 巖土工程學報, 2008, 28(8): 940?947.ZHU Wei, HUANG Zheng-rong, LIANG Jing-hua. Studies on shell-spring design model for segment of shield tunnel[J].Chinese Journal of Geotechnical Engineering, 2008, 28(8):940?947.
[12] 蘇宗賢, 何川. 盾構隧道管片襯砌內力分析的殼?彈簧?接觸模型及其應用[J]. 工程力學, 2007, 24(10): 131?136.SU Zong-xian, HE Chuan. Shell-spring-contact model for shield tunnel segmental lining analysis and its application[J].Engineering Mechanics, 2007, 24(10): 131?136.
[13] 朱合華, 劉庭金. 超淺埋盾構法隧道施工方案三維有限元分析[J]. 現(xiàn)代隧道技術, 2001, 38(6): 14?18.ZHU He-hua, LIU Ting-jin. 3D FEM analysis of shield driven tunnel with super shallow depth[J]. Modern Tunneling Technology, 2001, 38(6): 14?18.
[14] Blom C B M, Host E J, Jovanovic P S. Three-dimensional structural analyses of the shield-driven “Green Heart” tunnel of the high-speed line south[J]. Tunneling and Underground Space Technology, 1999, 14(2): 217?224.
[15] Klotz U, Vermeer P A. A 3D finite element simulation of a shield tunnel in weathered Singapore Bukit Timah Granite[J].Tunneling and Underground Space Technology, 2006, 21:272?283.