第五屆海峽兩岸隧道與地下工程學術與技術研討會 2006 年 3 月 台北 A23 - 1
砂岩隧道彈塑性模式之分析與探討
蔡立盛1 翁孟嘉2 鄭富書3 林銘郎3 1中華顧問工程司工程師 2 中興工程顧問社地工中心研究員 3台灣大學土木工程學系教授摘要
台灣西部麓山帶主要係為中低強度之第三紀以上年輕地層。這些地層岩性主 要以沉積岩類之雜砂岩為主，根據室內實驗結果，顯示此系列之砂岩具非常重要 且值得注意的變形行為，且此行為為一般等向性彈性模式或完全彈塑性模式所無 法合適描述。本文除了建議可合適描述砂岩力學變形行為特徵之組成律模式外， 亦撰寫可應用於數值分析之程式。透過配合數值分析程式模擬隧道開挖時之地盤 行為，分析結果顯示本研究之模式可合適於定性地描述隧道開挖後擠壓變形特 徵。此外本模式分析所得變形量約為線彈性模式及Drucker-Prager模式之1.7~2.9 倍，相對於線彈性模式過於保守，本模式之分析結果可更接近於實測之變形量。 關鍵字：非線性彈性、塑性、砂岩。一、前言
台灣西部麓山帶岩層以第三紀的沉積岩為主，由於該岩層成岩作用時間短， 岩石膠結及壓密程度不佳；復以台灣既溫且濕的氣候環境，風化作用快速進行， 使得此類岩層常呈膠結不良、孔隙率高與變形性大等特性，形成大地工程自然之 不利因素。根據北二高之中和隧道災變調查，該隧道通過之岩層為中顆粒膠結疏 鬆的木山層砂岩，由於岩性軟弱，加上該處砂頁岩互層發達，使得岩體自持力減 低，造成隧道嚴重擠壓變形，產生高達36公分之頂拱沈陷量(鄭富書等人，1996)。 一般等向線彈性模式，其承受純剪應力時，體積並不會發生變化。然而，軟 弱岩石在受受剪應力作用下，產生膨脹的剪脹現象。軟弱岩石此種剪脹的特性， 對工程分析而言相當重要，若忽略掉此效應會低估岩石的變形性。 以隧道為例，運用線彈性模式分析時，隧道周圍元素需承受相當大之剪應 力，若無考慮材料之剪脹特性，勢必錯估隧道周圍之變形特性。故本研究首先建 議可適用於描述砂岩剪脹特徵的組成律模式後，更進一步將此模式撰寫為可搭配 數值分析軟體Abaqus使用之副程式，以便工程應用於實務分析上。蔡立盛等人 砂岩隧道彈塑性模式之分析與探討 A23 - 2
二、砂岩之變形特徵
2.1 台灣西部麓山帶砂岩變形特性 本研究機構林傑(1997)、唐孟瑜(1998)、蔡立盛(2000，2005)及翁孟嘉(2002) 曾針對台灣西部麓山帶第三紀與第四紀砂岩之力學行為進行系列力學試驗， 力學主要以三軸純剪應力路徑為主，並配合其它應力路徑。三軸純剪應力路 徑之實驗方法可參考翁孟嘉(2002)及上述其他作者之文獻內容。根據實驗之 結果，可歸納出台灣西部麓山帶砂岩之彈塑性行為特徵如下： （a）彈性變形方面：體積應力作用階段，隨著體積應力增加，岩材之彈性體 積模數隨之提高。在剪應力作用階段，彈性剪力模數幾與初始剪力模數 一致，而呈一定值，不同體積應力下之剪力模數隨著體積應力增加而增 加。純剪應力與彈性體積應變間存在明顯的偶合關係，且此關係亦與體 積應力有關。 （b）塑性變形方面：體積應力作用階段，低應力下即產生頗大塑性應變增量； 隨著應力增加，增量逐漸下降，呈現塑性應變閉鎖之現象。剪應力作用 階段，塑性剪應變一開始不顯著，趨近破壞時大量剪應變始產生。低剪 應力造成之塑性體積應變呈現壓縮狀態，隨著剪應力增加通過塑性剪脹 門檻後，體積開始大量膨脹直至破壞。 （c）根據受剪時之塑性流分佈，可推知砂岩之塑性勢能面初時呈橢圓帽狀， 趨近破壞時則幾與破壞包線重合；觀察與破壞包絡線交會之塑性流，可 發現其與破壞包絡線大致呈正交，此現象可大膽假定砂岩之塑性應變增 量符合諧和流動準則。三、砂岩彈塑性組成律模式
3.1 非線性彈性組成模式 3.1.1 高階格林彈性模式 砂岩之彈性體積變形與體積應力間存有非線性關係，且純剪應力與彈性 體積應變間存有偶合關係；然而在現有之等向性線彈性模式中其體積應 變僅與體積應力有關，而與剪應力無關。為反應上述現象，本文採用高 階格林彈性模式(Green Elastic Model；Hyperelastic Model)，建議一適 用於砂岩變形特徵之能量密度函度。格林彈性模式其應變張量式可表示 為：第五屆海峽兩岸隧道與地下工程學術與技術研討會 2006 年 3 月 台北 A23 - 3 ij ij ε σ ∂Ω = ∂ （1） 上式中Ω為一能量密度函數，可表示成三個應力不變量之型式，εij及σij 分別為應變與應力張量。翁孟嘉(2002)根據圍壓40MPa實驗觀察所得之 岩石彈性變形行為，選定合適台灣麓山帶砂岩之能量密度函數Ω為： 3/ 2 1 1 2 1 2 3 2 bI b I J b J Ω = + + （2） 上式中b1、b2及b3三項為材料常數，I1為第一應力不變量，J2 則為第二 偏差應力不變量。蔡立盛(2005)則進一步根據不同圍壓之實驗資料，將 能量密度函數建議修正為： 3/ 2 1 1 1 2 1 2 3 2 bI b I−J b J Ω = + + （3） 將式(3)代入式(1)，可得非線性彈性之應力–應變關係，如下所示： 1 2 1 2 3 11 2 1 2 2 1 3 ij 2 ( ) ( ) s ij ij ij b I b I J b I b ε δ σ − − ∂Ω = = + + + ∂ （4） 其中sij為偏差應力張量，而δij為Kronecker Delta張量。由式(4)可進一 步推求純剪應力路徑試驗，引致之體積應變及剪應變，如下所示： （a）體積應力作用階段 1 2 9 1 kk 2 1 1 ' I I = ε = b （5） 2 2 J' 0 γ= = （6） 上式中 ' 1 I 為第一應變不變量，代表體積應變。γ 為剪應變， ' 2 J 則為第二 偏差應變不變量。式(5)反應當體積應力增加，體積模數隨之遞增，體積 應變與應力呈一非線性關係。式(6)顯示體積應力作用並不導致剪應變產 生。 （b）剪應力作用階段 2 1 2 2 1 ' -3 (J /I ) I = b （7） 1 2 1 3 2 2(b I b) J γ = − + （8） 上式(7)顯示純剪應力作用下，彈性體積應變與剪應力及體積應力均存有 偶合關係，且為一非線性剪脹行為。式(8)則反應剪應力與彈性剪應變呈 一線性關係，隨著體積應力愈高，其彈性剪力模數愈高。模式充份反應 砂岩之彈性變形特性，此項特色已非一般線彈性模式可加以模擬。
蔡立盛等人 砂岩隧道彈塑性模式之分析與探討 A23 - 4 3.1.2 非線性彈性模式之增量型式 為了將上述之彈塑模式運用於有限元素分析中，本節介紹如何推導彈性 模式增量形式之應力應變關係。由塑性力學之定義，增量型式之應力應 變關係可表示為(Chen，1982)： ε• σ• ⎧ ⎫ ⎧ ⎫ ⎪ ⎪= ⎡ ⎤⎪ ⎪ ⎨ ⎬ ⎣ ⎦⎨ ⎬ ⎪ ⎪ ⎪ ⎪ ⎩ ⎭ Dt ⎩ ⎭ （9） 其中⎡ ⎤⎣ ⎦Dt 代表應力應變增量式間的切線柔度矩陣。以能量密度函數而 言，一般增量形式的應力-應變的關係為： 2 2 ij kl ijkl kl ij kl ij kl ijkl kl ij kl W H H σ ε ε ε ε ε σ σ σ σ ∂ = = ∂ ∂ ∂ Ω ′ = = ∂ ∂ （10） 其中σij為應力增量；εij為應變增量，Hijkl為應力-應變間的切線勁度矩 陣，H ′ijkl為應力應變間的切線柔度矩陣。首先由式(4)得到非線性彈性之 應力應變關係，令成如下式(11) 1 2 sij ij ij ε = φ δ + φ （11） 其中 0.5 2 1 1.5 b I1 1 b I2 1 φ = − − ； 1 2 b I2 1 b3 φ = − + 式(11)中的張量式對應力張量再偏微分一次，即可得到應力-應變增量的 關係，如下所示： ij 1 2 2 ij s s ij ij kl kl kl kl φ φ ε δ φ σ σ σ σ • ⎡∂ ∂ ∂ ⎤ • =⎢∂ + ∂ + ∂ ⎥ ⎣ ⎦ （12） 其中 sij 1 3 ik jl kl ij kl δ δ δ δ σ ∂ = − ∂ ； 1 1 2 11 2 3 4 kl kl kl kl kl b I b S A BS φ δ δ σ − ⎡ ⎤ ∂ = + = + ⎢ ⎥ ∂ ⎣ ⎦ ； 2 2 1 B= −b I− ； 1 3 2 11 2 1 2 3 2 4 A= b I− + b I− J ； 2 -2 12 kl kl kl b I S BS φ σ − ∂ = = ∂ 將式(12)展開得到增量形式的非線性彈性之應力-應變關係如下式所示： 2 2 1 3 ij A ij kl B ij klS ik jl ij klS B kl ijS kl ε• =⎢⎡ δ δ + δ +φ δ δ − φ δ + δ ⎤⎥σ• ⎣ ⎦ （13） 因此式(13)即可求得適用於描述砂岩彈性變形行為之切線柔度矩陣及切 線勁度矩陣，接著本研究即採用Fortran程式語言，撰寫增量型式應力應 變關係的程式。
第五屆海峽兩岸隧道與地下工程學術與技術研討會 2006 年 3 月 台北 A23 - 5 3.2 塑性組成律模式 3.2.1 降伏面 翁孟嘉(2002)採用之降伏面如式(14)所示，該降伏面係以Desai (1986, 1987a, 及1987b)所提出之
δ
0模式(或稱HiSS模式)為基礎所發展而得。 塑性模式之降伏面方程式如下：(
I1,J2)
=J2 −[
− (I1+T) + (I1+T)2]
=0 F αD n γD （14） 其中γD為一變數，由硬化準則所決定；γD為破壞包絡線斜率之平方； T為破壞包絡線與I1軸截距；而n則為控制降伏面形狀之常數。 3.2.2 塑性勢能面 此外翁孟嘉(2002)根據實驗結果，觀察砂岩受剪變形最終之塑性流幾與 破壞包絡線正交，因此假定塑性應變增量符合諧和流動準則(Associated Flow Rule)，亦即塑性勢能面G I( ,1 J2)=0與降伏面F I( ,1 J2)=0如式(15)。 ) , ( ) , (I1 J2 F I1 J2 G ≡ （15） 3.2.3 硬化準則 針對台灣麓山帶砂岩特性，其適用之硬化準則，如式(16)所示： 3 2 1 4 1 p ij D D d β ββ
ε
β α
α
⎛ ⎞ = + ⎜ ⎟ ⎝ ⎠ （16） 式(16)中右式之第一項，由β1、β2兩材料常數控制，其反應體積應力 作用下，塑性應變增量逐漸減小之情形；右式之第二項，為δ
0模式中所 建議之硬化準則，其反應剪應力作用下，塑性應變增量之變化情形，由 β3、β4兩常數控制。由於塑性模式已為增量型式，因此無需再推導其 增量型式的應力應變關係。 3.3 單元素驗證 為驗證程式之正確性，本研究將木山層砂岩純剪試驗實驗曲線資料與進行純 剪路徑的單元素數值模擬結果進行比較，一方面檢核此組成律模式是否可以 描述砂岩變形行為；另一方面驗證程式之正確性。 彈性行為方面，圖一至圖三為數值模擬結果與加解壓曲線試驗資料之比較。蔡立盛等人 砂岩隧道彈塑性模式之分析與探討 A23 - 6 圖一顯示程式可適當描述體積應力作用下之非線性行為。圖二顯示了不同圍 壓下，純剪應力與彈性體積應變關係之模擬結果甚佳，此結果亦為一般之彈 性模式所無法模擬。圖三則顯示程式亦可反應彈性剪力模數值會隨圍壓增加 之特性。 總變形之數值模擬結果如圖四至六所示。觀察各應力階段之結果顯示，無論 是體積應力或純剪應力作用階段，ABAQUS程式模擬曲線與實驗資料均相當 的一致。由此可知模式可適當地模擬砂岩變形行為，同時亦驗證了程式之正 確性。此外，隧道開挖過程中，砂岩並不會僅遵循純剪應力路徑，而是會各 種應力路徑。關於此項疑問，翁孟嘉(2002)曾針對不同應力路徑，將預測結 果與實驗資料進行比對，結果顯示預測值與實驗值相當的一致。這也說明了 此彈塑模式可適用於描述砂岩不同應力路徑之變形行為。 0 10 20 30 40 50 60 70 0 2000 4000 6000 8000 10000
Elastic volumetric strain (10-6)
Hydrostatic stress (MPa)
MS60-15 MS60-30 MS60-45 MS60-60 Regressional curve ABAQUS model , e v p ε R2=0.97 1 1119.9 , 3 I e v p ε = 0 10 20 30 40 50 60 -600 -500 -400 -300 -200 -100 0 100
Elastic volumetric strain (10-6)
S h ear st ress ( J2 ) 1/ 2 (MP a ) MS20 MS40 MS60 MS60-ABAQUS model MS40-ABAQUS model MS20-ABAQUS model 圖一 體積應力作用下之木山層(MS2)砂岩 之彈性體積應變模擬曲線 圖二 剪應力作用下木山層(MS2)砂岩之 彈性體積應變模擬曲線 0 10 20 30 40 50 60 0 1000 2000 3000 4000 5000
Elastic shear strain 2*(J'2)
1/2 (10-6) Shea r st ress ( J2 ) 1/2 (M Pa) MS20 MS40 MS60 MS60-ABAQUS model MS40-ABAQUS model MS20-ABAQUS model 0 5 10 15 20 25 30 35 40 45 0 2000 4000 6000 8000 10000 12000 14000 Volumetric strain (10-6) Volumetric stress (M Pa) Weng ( 2002) ABAQUS analysis MS2 sandstone Hydrostatic stress stage p=40 MPa
圖三 剪應力作用下木山層(MS2)砂岩之彈 性剪應變模擬曲線
圖四 不同體積應力作用下木山層砂岩 之彈性體積應變模擬曲線
第五屆海峽兩岸隧道與地下工程學術與技術研討會 2006 年 3 月 台北 A23 - 7 0 10 20 30 40 50 60 -5000 -4000 -3000 -2000 -1000 0 1000 Volumetric strain (10-6) Shear stress (J 2 ) 1/ 2 (MPa) Weng (2002) ABAQUS analysis MS2 sandstone TC test, p = 40 MPa 0 10 20 30 40 50 60 0 2000 4000 6000 8000 10000 12000 14000 Shear strain 2*(J'2) 1/2 (10-6) She ar str es s (J 2 ) 1/2 (MPa) Weng (2002) ABAQUS analysis MS2 sandstone TC test, p = 40 MPa 圖五 剪應力作用下木山層砂岩之彈性剪 應變模擬曲線 圖六 剪應力作用下木山層砂岩之彈性 體積應變模擬曲線
四、工程案例分析
4.1 隧道案例模型之建立 本研究採用台灣北二高甲隧道作為案例，以平面應變數值分析為基礎，探討 本 研 究 彈 塑 組 成 律 模 式 ( 簡 稱 NTU) 與 等 向 線 彈 性 模 式 ( 簡 稱 Elastic) 及 Drucker-Prager彈塑性模式(簡稱D-P)在隧道開挖變形行為上之差異性。隧道 案例之數值模型長度為200公尺，高度為160公尺，隧道覆土深度約為60公 尺，隧道斷面寬為16公尺，高為11.4公尺。 4.2 材料參數之選取 本研究分析採用之參數可分成直接由三軸試驗所得之參數(簡稱三軸參數)及 考量現場野外岩體為具有節理情況參數(簡稱岩體參數)。綜合整理各模式之 參數如表一所示。表一中之三軸參數係採木山層砂岩經由室內之單壓及純剪 三軸應力路徑求得。岩體參數值則是參考中興顧問公司分析木柵一號隧道與 中華顧問公司分析中和隧道所採行之岩體參數。 4.3 隧道案例分析結果討論 隧道案例分析之結果本文分為兩個主題進行討論之，包含(1)整體變形行為之 定性討論及(2)特定位置如頂拱、仰拱及側壁之變位量的定量討論。由此兩種 主題之比較結果來彰顯各模式之特色及優缺點。 定性分析整體變形行為方面，本研究觀察隧道開挖後的體積應變及剪應變之 分佈情形。其中剪應變部份，不同模式的剪應變集中區堿大致類似，彼此之 間之差異僅為剪應變值的大小。例如同為三軸參數時，本研究建議之NTU模 式求得之剪應變會大於D-P模式及Elastic模式。體積應變分佈情況方面，不 同模式則會有明顯的差異性。圖七為不同模式分析結果之總體積應變分佈情蔡立盛等人 砂岩隧道彈塑性模式之分析與探討 A23 - 8 況。由圖中可知，隧道頂拱及仰拱部份，各模式分析結果均為體積膨脹，亦 即擠壓的變形行為，然而在側壁部份，圖七(a)與(d)之側壁之膨脹量相對較 小，甚至有些區域之體積應變接近於零。圖七(b)為D-P模式配合岩體參數分 析結果，顯示隧道周圍均為擠壓行為，但是大多數集中在側壁。圖七(d)為NTU 模式分析結果顯示，隧道周圍為均勻的體積膨脹行為，且愈靠近隧道表面之 處，其體積膨脹量愈大，此現象與現地觀察到之隧道變形行為十分相似。由 以上討論可知，本研究建議之NTU模式應用在隧道分析上時，可以合適定性 地描述隧道開挖後周遭岩石的變形行為。由於各模式的彈性體積應變分佈情 形相當一致；因此不同模式總體積應變分佈情形之所會有如此大差異性，其 原因主要係歸咎於各模式之塑性應變值的計算結果不同所致。圖八為不同模 式分析結果之塑性體積應變的分佈圖。由圖中顯示圖八(a)與(c)之塑性體積應 變均為零，而(b)與(d)之圖則有明顯之塑性變形區域。 在定量分析隧道的變形量方面整理如表二所示。表中顯示岩體參數時，不管 是頂拱、仰拱或側壁，D-P模式之變位量均大於Elastic模式結果，本模式無 岩體之參數，因此僅由三軸參數比較不同模式時的變形量，本模式預估之變 形量比線彈性和彈塑性均較大，約為1.7~2.9倍，會更接近實測之隧道變形量。 表一 各模式之材料常數 模式 參數 三軸參數 岩體參數
E：彈性模數 7500 MPa 200 MPa 線彈性模式 ν：柏松比 0.21 0.21 k：破壞包絡線截距 11.6 Mpa 1.4 MPa φ：摩擦角 39 39 D-P 模式 ψ：剪脹角 39 39 b1：非線彈模式材料參數 143.7 µ⋅(MPa)−0.5 − b2：非線彈模式材料參數 1462.8 − b3：非線彈模式材料參數 30.2 µ⋅(MPa)−1 − γD：破壞包絡線斜率之平方 0.12 − T：破壞包絡線與 I1軸截距 40.7 MPa − n：降伏面形狀參數 3.47 − β1：硬化準則之參數 1.3E+5 − β2：硬化準則之參數 2.52 − β3：硬化準則之參數 0.64 − NTU 彈塑模式 β4：硬化準則之參數 7.34E-8 − 表二 不同模式分析隧道開挖變形量之結果 模式 頂拱擠壓 (cm) 仰拱擠壓(cm) 側壁擠壓(cm) 三軸參數 岩體參數 三軸參數 岩體參數 三軸參數 岩體參數 線彈性模式 0.168 6.28 0.335 12.55 0.162 6.07 D-P 彈塑模式 0.168 8.09 0.335 13.63 0.162 8.88 NTU 彈塑模式 0.41 − 0.59 − 0.48 −
第五屆海峽兩岸隧道與地下工程學術與技術研討會 2006 年 3 月 台北 A23 - 9 圖七 隧道開挖後，不同模式之總體積應變區分佈圖。其中Elastic代表採用線彈 性模式，D-P代表Drucker-Proager模式，NTU代表本研究之彈塑模式。 圖八 隧道開挖後之塑性體積應變區分佈圖。其中Elastic代表採用線彈性模式， D-P代表Drucker-Proager模式，NTU代表本研究之彈塑模式。
蔡立盛等人 砂岩隧道彈塑性模式之分析與探討 A23 - 10
五、結論
本研究建議之彈塑性模式不僅較一般線彈性模式更能合理描述砂岩材料的 彈性及塑性行為。利用單元素進行數值分析和試驗資料進行比對後，結果顯示兩 者之值相當一致，由此可驗證所撰寫的使用者自訂材料副程式可以有效地模擬砂 岩材料之彈性與塑性力學行為，反應大地材料特有之變形特徵。 隧道案例定性分析方面，本研究彈塑模式分析成果顯示隧道周圍為均勻體積 膨脹行為，且愈靠近隧道表面之處，其體積膨脹量愈大，此現象與現地觀察到之 隧道變形行為十分相似。在定量分析方面，岩體參數分析結果顯示，線彈性模式 及Drucker-Prager模式分析結果與實測值均有甚大之差距。以頂拱為例，線彈性 模式頂拱沉陷量約為6.28公分，Drucker-Prager模式約為8.09公分，此與現場實 測頂拱沉陷約36公分之差距仍有一段距離。而由三軸參數分析結果顯示，本研究 模式分析所得之頂拱變形量約為線彈性模式及Drucker-Prager模式之2.44倍，依 此倍率關係推估，本研究之模式如可得適當之岩體參數時，頂拱沉陷量估計約為 19.7公分。綜合本研究成果顯示，本文建議之彈塑模式無論由定性或定量方面， 模擬結果均應可以更近似於實際之隧道變形行為。參考文獻
1. 林傑 (1997)「軟弱砂岩剪脹及潛變行為研究」，碩士論文，台灣大學土木工程學研究所，台 北。 2. 唐孟瑜(1998)，「木山層砂岩變形行為研究」，碩士論文，台灣大學土木工程學研究所，台北。 3. 翁孟嘉(2002)，「麓山帶砂岩之力學特性及其與微組構關係研究」，博士論文，台灣大學土木 工程學研究所，台北。 4. 蔡立盛、鄭富書、林銘郎、黃燦輝、翁孟嘉 (2005)， 砂岩之非線性彈性變形及其組成律模式 研究，第十一屆大地工程研討會，台北縣萬里鄉。 5. 鄭富書、林銘郎、黃燦輝、洪如江 (1996)，「長隧道工程技術之研究—台灣北部隧道天然施 工障礙之研究」，國立台灣大學土木系研究報告，交通部。6. Chen, W.F. and Saleeb, A.F. (1982), Constitutive equations for engineering materials, Vol. I: Elasticity and modeling, John Wiley & Sons.
7. Desai, C.S., Somasundram, S. and Frantziskonis, G. (1986), “A hierarchical approach for constitutive modeling of geologic materials,” Int. J. for Numerical and Analytical Methods in Geomechanics, Vol. 10, pp. 225-257.
8. Desai, C.S. and Salami, M.R. (1987), “A constitutive model and associated testing for soft rock,” Int. J. Rock Mech. Min. Sci. & Geomech. Abstr,. Vol. 24, pp. 299-307.
9. Desai, C.S. and Varadarajan, A. (1987), “A constitutive model for quasi-static behavior of rock salt,” J. of Geophysical Research, Vol. 92, pp. 11445-11456.
