高效能大地工程有限元素分析之研究(III)

行政院國家科學委員會專題研究計畫 成果報告

高效能大地工程有限元素分析之研究(3/3) 研究成果報告(完整版)

計 畫 類 別 ： 個別型

計 畫 編 號 ： NSC 96-2221-E-011-135-

執 行 期 間 ： 96 年 08 月 01 日至 97 年 07 月 31 日 執 行 單 位 ： 國立臺灣科技大學營建工程系

計 畫 主 持 人 ： 謝佑明

計畫參與人員： 碩士班研究生-兼任助理人員：孫煥期 碩士班研究生-兼任助理人員：宋俊宜 碩士班研究生-兼任助理人員：林弘倫 博士班研究生-兼任助理人員：潘茂森

處 理 方 式 ： 本計畫涉及專利或其他智慧財產權，1 年後可公開查詢

中 華 民 國 97 年 10 月 31 日

目錄

目錄... I  表目錄... II  圖目錄... III 

一、前言及目的... 1 

二、文獻探討... 4 

2.1 有限元素法與網格調整... 4 

2.2 內插方法之研究... 5 

2.3 分析後誤差評估（a posteriori error estimation） ... 7 

2.4 METIS... 9 

三、網格演化法... 9 

3.1 網格演化之研究方法... 9 

3.2 網格演化之系統分析流程... 11 

3.2.1 網格輪廓之擷取... 13 

3.2.2 應力場之解映射... 17 

3.3 程式架構... 18 

3.4 案例分析... 19 

3.4.1 二維度深開挖分析... 19 

3.4.2 三維度深開挖分析... 25 

3.4.2 三維度隧道開挖分析... 28 

3.5 效能探討... 29 

3.5.1 二維度深開挖分析之效能探討... 29 

3.5.2 三維度深開挖及隧道開挖分析之效能探討... 30 

四、網格誤差評估方法之研究... 33 

4.1 誤差評估研究方法... 33 

4.2 誤差評估程式架構... 34 

4.3 誤差評估成果... 35 

五、METIS 網格切割... 38 

六、結論... 40 

參考文獻... 41 

計畫成果自評... 47 

表目錄 

表 1 土壤及混凝土參數... 24 

表 2 執行平台構架及選用參數表... 24 

表 3 二維傳統有限元素分析之計算時間... 30 

表 4 二維網格演化有限元素分析之各階段計算之時間... 30 

表 5 三維深開挖之傳統有限元素分析之計算時間... 31 

表 6 三維深開挖之網格演化有限元素分析之各階段計算之時間... 31 

表 7 三維隧道開挖之傳統有限元素分析之計算時間... 31 

表 8 三維隧道開挖之網格演化有限元素分析之各階段計算之時間.. 32 

表 9 網格切割後元素與節點數量表... 39 

圖目錄 

圖 1 網格密度與應力值或變位值... 3 

圖 2 各元素較佳應力的位置... 9 

圖 3 網格演化法之分析流程概念... 10 

圖 4 一二向度深開挖施工順序... 11 

圖 5 網格演化法之整體分析流程... 13 

圖 6 節點移動之區域... 14 

圖 7 單一方向調整之流程... 15 

圖 8 節點單一方向調整方式... 15 

圖 9 以內插方法調整方式... 16 

圖 10 以內插方法做調整之流程... 17 

圖 11 不同材料參數設定(A)~(C)之區域 ... 18 

圖 12 網格模型調整類別核心程式架構... 18 

圖 13 解映射類別核心程式架構... 19 

圖 15 二向度深開挖案例模型... 21 

圖 16 傳統分析方法第九階段位移增量及破壞面... 22 

圖 17 網格演化分析方法第九階段位移增量及破壞面... 22 

圖 18 擋土牆位移及地表之沉陷變化量... 23 

圖 19 三維空間下之深開挖模型... 26 

圖 20 第 20 階段開挖區域之變位增量... 26 

圖 21 第 20 階段開挖之擾動與總變位量... 27 

圖 22 第 20 階段開挖每面擋土牆之變位增量... 27 

圖 23 第 20 階段開挖每面擋土牆之擾動量... 27 

圖 24 三維下之隧道開挖模型... 28 

圖 25 第 30 階段開挖之變位增量... 29 

圖 26 第 30 階段開挖之擾動量及總變位量... 29 

圖 27 誤差關聯建立流程圖... 34 

圖 28 誤差評估程式架構... 35 

圖 29 平面應變問題... 36 

圖 30 不同網格密度下之元素相對誤差... 37 

圖 31 開挖邊上的最大橫向位移的收斂... 37 

圖 32 地中最大位移的收斂... 37 

圖 33 以 SPR 計算誤差評估值 ... 38 

圖 34 以 SEEELEM 計算誤差評估值... 38 

圖 35 以 SPR 預測最大地中變位的 實際誤差值 ... 38 

圖 36 以 SEEELEM 預測最大地中變位的 實際誤差值... 38 

一、前言及目的 

人口密集的都會區裡，建築物林立，新建工程與鄰近建築物或設施之距離往 往十分有限，而在施工過程中難免對於鄰房及鄰近設施有所影響，此影響必須在 工程開始前進行評估，並透過成功的工程設計以確保在施工的過程中以及完工後 不致影響鄰近設施的服務性與安全性。

傳統大地工程分析方法如極限平衡法(limit equilibrium)著重於穩定性分析，

並透過安全係數的方式以克服土壤性質之不確定性以及考慮對鄰近設施的影響。

此設計方式並非直接地考慮深開挖工程對於鄰近設施的影響，而安全係數的選擇 主要仰賴過去工程經驗而無法以依客觀而合乎物理原則的方法決定，因此可能造 成過度設計，付出不必要的成本；亦有設計不足的可能，使鄰近設施受損，而必 須付出額外的成本。

有限元素法(finite element method)為目前廣泛應用於各領域之研究分析的數 值方法，大地工程常需面對複雜的幾何形狀、材料性質的變異性、材料的非線性 行為等等，有限元素法均可合理地考量這些特性，並可由分析的結果得到完整的 變位場(displacement field)與應力場( stress field )等相關資訊；但目前無論學術界 或是業界所廣泛使用的有限元素分析多為使用平面應變假設之二向度分析，但由 於實際工程問題多為三向度的問題，二向度之有限元素分析必須對實際工程問題 做出工程判斷或是利用經驗以決定適當的二向度幾何形狀、支撐勁度或施工順序 等等。而這些判斷或經驗也降低了有限元素分析方法的嚴謹性與可靠性。因此最 完整、所需假設最少的分析方式乃採用三向度有限元素分析。然而大地工程的三 向度有限元素模型具有大量的自由度，因此需要冗長的計算時間。

根據經驗，深開挖分析上側向邊界使用開挖深度的 5~10 倍距離為一適當的 距離，因此若一深開挖工程為 30 米深，則有限元素模型必須在側向上跨越 300~600 米的距離；若所考量的工程為隧道開挖，則有限元素模型之隧道向長度 可能為數百米至千米，因此大地工程之有限元素模型需涵括廣大的範圍。而工程

分析上所需考慮的工程行為與工程設施相對微小，如前例，一 30 米之深開挖可 能假設 5 至 10 階段進行數值模型的開挖，則每一階段僅為 3 至 6 米，而所用來 支撐的連續壁系統可能僅有 0.3 米至 1 米的厚度，因此，相對於有限元素分析網 格所涵括的範圍，其所需描述的工程行為與設施極為微小。因此，大地工程之有 限元素模型實際上需跨越多個長度尺度，而造成三向度有限元素分析，由於所需 求解的自由度數量極為龐大，即使進行線彈性分析亦需要冗長的計算時間，而大 地工程問題之非線性更使得計算時間成倍數成長。因此三向度有限元素分析無論 國內外都極少被應用於大地工程研究或案例分析上，而文獻中少數所能看到的三 向度有限元素分析常使用及為粗糙的網格，而降低了其分析結果之準確性與可靠 性。

因此為了解決大量自由度的問題，本研究將根據大地工程有限元素分析網格 幾何形狀變動之特性進行改善，提出「網格演化法」[1][2]之分析程序，以期能 降低大地工程三維度有限元素分析所需使用的自由度量，進而改善龐大計算量所 需耗費的時間。

本研究並延伸先前開發之軟體系統，配合電腦叢集系統使用平行運算之技巧 以減少三向度有限元素分析所須計算時間，此平行計算的方式乃採用無重疊之領 域分解法，但在此軟體系統中之領域分解步驟，是採用人工進行之二向度平面分 解並配合 extrusion 已完成三向度的有限元素模型之分解，此平面人工分解的方 式雖然對於切割出來的子領域(subdomain)形狀與介面性質有較佳的掌握，但相當 費工時。本研究使用 METIS 之公開程式庫，自動化領域分割演算法，並與先前 研究之平行化有限元素分析系統進行整合，使整個平行化有限元素分析系統自動 化。

本研究除了在效能改善上提出解決方案外，另外在網格品質方面也進行探討 [3]。有限元素法為一種用來評估和分析實際的問題的數值方法，其分析的誤差 會隨著網格密度的增加而減小，即收斂行為。傳統上，研究者利用此ㄧ特性來找

分析結果中的值（如應力值或變位值）也隨之改變，當網格密度達到一定程度時，

分析結果中的應力值或變位值會維持定值或改變量極小，此即達到分析的最精確 解（收斂解）。

圖 1 網格密度與應力值或變位值

上述方法最後可完全消除網格密度所造成的誤差，但為了得到分析的最精確 解，研究者可能會花費許多時間與計算成本來尋找最佳的網格密度，而且工程上 並非所有問題都需找到有限元素法的最精確解，大部分的情況為網格密度與分析 結果符合工程師的要求即可。若是要找出分析問題中各個位置的網格最佳密度

（如開挖問題中的擋土壁需較高的網格密度），則需花費更多時間。

為了使網格演化法之網格變動找出最合適的分析結果，本研究將針對深開挖 與隧道的問題根據有限元素法的分析結果計算誤差值，並建立誤差值與實際誤差 值之間的關係，可使工程師在短時間內根據可容忍的實際誤差值建立適當的網格 密度，並使網格演化法之不同階段網格品質穩定。

h

σ

) (u

二、文獻探討 

2.1 有限元素法與網格調整

有限元素法 (Finite Element Method, FEM) 為目前廣泛應用於各領域之研究 分析的數值方法，對於大地工程常需面對複雜的幾何形狀、材料性質的變異性、

材料的非線性行為等等，有限元素模型均可合理地考量這些特性，並可由分析的 結果得到完整的變位場 (displacement field) 與應力場 (stress field) 等相關資訊，

以做為工程設計之依據與參考。此實為解析方法 (analytical solution) 或極限平衡 法等其他分析方法所無法達成的，解析分法目前僅能考慮簡單的幾何形狀與材料 行為 (常是線彈性行為)；而極限平衡法則僅考慮應力狀況而忽略了土壤與擋土 牆結構變型所造成之影響。因此，有限元素法已被廣泛應用於深開挖問題的探討 [4][5][6][7] 與隧道分析[8][9][10][11][12][13][14] 等大地工程問題上。

網格調整法為有限元素法中之一特殊技巧。有限元素法雖可以應用於分析工程上 的許多問題，但對於網格中具有局部大變形 (large deformation) 上之處理，有限 元素法有其先天上的限制存在：1) 其會因網格中部份的有限元素扭曲過大，而 造成在計算元素積分時之不精確，更甚而造成計算結果無法收歛；2) 應力集中 的發生，因此應力集中處需要高精度的網格以避免應力描述之不良；3) 當分析 中考量到材料之塑性行為時，會有局部化問題 (localization problems)，造成有限 元素解相依於有限元素網格之切割而有截然不同的答案。而在有限元素法架構下 處理此類問題的方法目前可主要分為兩類：1) ALE (Arbitrary Lagrangian-Eulerian) [15][16] 法，與 2) 網格調整法。

網格調整法在處理大變形問題分析上是一個相當直觀的方法。此方法的概念 在於在有限元素法分析進行當中，「監視」有限元素網格之品質，當發現有限元 素網格中有元素發生過度之扭曲時，即中斷此次的有限元素分析，並根據目前之 分析結果將網格進行調整以產生一新的、元素形狀良好的網格，以「恢復」有限 元素網格之品質，以期能維持有限元素解之正確性。當新的網格產生了之後，則

必需將之前中斷的有限元素分析結果以「解映射」(solution mapping) 的過程映 射至新產生的網格上，而此解映射實質上為一內插或外插的過程，內插計算方式 將於 2.4 節討論。

早期所提出之網格調整法主要以人工的主觀判斷網格調整之時機。而後有學 者提出所謂適應性網格調整法 (adaptive remeshing) 的概念。此觀念在於訂定出 一 量 化 之 網 格 調 整 指 標 ， 例 如 網 格 之 應 力 誤 差 評 估 子 (error estimator) [17][18][19]，爾後由電腦程式根據所使用之量化標準自動調整網格。

網 格 調 整 法 常 被 使 用 於 機 械 工 程 上 之 金 屬 成 形 (metal forming) 分 析 [20][21] 及衝擊模擬 [22]，在大地工程上鮮少被應用。打擊樁 (driven pile) 貫入 地層屬於大變形問題，在樁身附近會造成土層中樁身附近網格高度扭曲，然並未 能於文獻中找到此分析方法之應用。目前唯有基腳承載力問題 [23][24][25] 有應 用此技巧進行研究分析。其它土木上應用網格調整之技巧者僅有動態載重於樑上 局部化問題 [26]。

而為了滿足使用者對於網格調整之需求，商業化之有限元素分析軟體也常會 提 供 網 格 調 整 之 分 析 能 力 。 以 ABAQUS 為 例 ， ABAQUS 提 供 了 *MAP SOLUTION 之功能進行解映射的動作，並於 ABAQUS/CAE 之前處理程式提供 了幾何形狀擷取與再劃分之功能。於 ABAQUS 範例問題中，即有以此方法，

對一鋼模壓縮的問題來進行模擬 (cylindrical billet upsetting)，此模型之標準試驗 根據 Lippmenn (1979)[27] 所定義，而其利用網格調整之分析結果與 Taylor (1981) [28] 的 結 果 近 似 。 而 最 新 之 ABAQUS 6.6[29] 版 更 加 入 了 adaptive remeshing 之功能，以方便大變形問題之模擬。

2.2 內插方法之研究

在 2.3 節中，為了解決有限元素之大變形問題，使用調整網格的方式，讓網 格重新產生，並將上一階段積分點上的解 (常為應力場)，傳遞到下一階段去，

其中傳遞的方式就是使用解映射 (solution mapping)的方法，而解映射的方法就

是內插 (interpolation) 或外插 (extrapolation)。

Hu 與 Randolph (1998)[24] 研究以網格調整法分析基腳之極限承載力，並 評估了以反轉距離法 (IDA, Inverse Distance Algorithm)、任意線性內插法 (ALI, Arbitrary Linear Interpolation)、以及單一元素內插法 (UEM, Unique Element Method) 之解的比較。他們的研究認為，單一元素內插法所得之結果較佳，網格 調整對於解之影響較規則而建議使用。

內 插 與 外 插 為 常 見 之 數 值 資 料 處 理 方 法 ， 因 此 一 般 的 數 值 分 析 書 籍 [30][31][33] 皆會介紹基本的數值內插方法，如拉格朗治 (Lagrange) 差值法、牛 頓 (Newton) 差值法、葉特肯( Aitken)差值法、與曲線近似法等。其它文獻中亦 有提及資料交換的方法[20][21] 等。

曲線近似法 (curve-fitting) 乃先假設一數學函數 (或稱基底函數，常使用多 項式、指數函數形式、與楔形曲線 spline) 以近似已知之資料點，而待此數學函 式決定後，即可利用來進行內插。在曲線近似中，基底函數中之待定係數常藉由 已知點的資料以最小平方法來決定之，故本文中將以曲線近似求取內插之方法稱 為最小平方近似法 (Least Square Method, LS)。此方式是以單一個數學函數來近 似與描述整體之資料場，故內插結果完全取決於所假設之數學函數是否適用於整 體地描述該資料。

Lancaster 與 Salkauskas (1981)[33] 提出變動式最小平方近似法 (Moving Least Square Method, MLS)，其與最小平方近似法不同處在於以多個相同形式但 係數不同之基底函數乘上一權重函數後疊加以近似已知之資料點，並可利用來進 行內插。MLS 中之權重函式的選擇也許多種，如指數型態加權函式、高階多項 式等，Krysl 與 Belytschko (2001)[34] 之文獻，提供了許多的權重函式。

變動式最小平方近似理論的內插方法，由於其內插出來之資料場 (data field) 可輕易透過其基底函數與權重函數進行一次或多次的微分，且其內插出來的資料 場具有高度連續性，因此近年來常被應用於各類型之無網格 (mesh free)分析方法，

節點內插法(Point Interpolation Method, PIM)，為 Liu(1999)[39]提出了進一步 改良的內插計算方法，其中依基底函式型態的不同又分為 Polynomial PIM 和 Radial PIM。

前述以曲線近似為基礎之內插方法常以近似曲線與已知資料點之誤差量最 小化以決定待定係數。但亦有研究者提出以場之連續性為優先之內插方法，如 Lekien 與 Marsden (2005)[40] 針對海洋洋流資料之內插，以整個內插資料場在 三維度空間中都可以保持是 C¹連續作為條件，而內插出海岸洋流的狀態。

2.3 分析後誤差評估（a posteriori error estimation）

數值方法的分析前誤差評估（a priori error estimation）長期以來一直是研究 者以數值方法分析問題必要工作之ㄧ，這些評估在網格參數作適當變化的時候能 提供不同分析軟體收斂性及穩定性的資訊，同時也能提供概略的誤差行為資訊。

傳統上進行數值模擬時如果預測到誤差的存在並不會去量化其值，而是靠研究者 的經驗去控制參數以及數值模擬結果的品質，但依賴研究者的經驗以及分析前對 誤差的預測仍有許多不確定因素存在，所以在分析後針對結果作誤差評估（a posteriori error estimation）變成一種改進分析品質較可靠的方式。在有限元素法 的分析中，研究者較為關心的通常是分析結果的導數（derivatives）而非分析結 果本身，以彈性問題為例，應力與應變是研究者較為關心的對象，而不是變位。

另外，分析結果的導數在元素之間通常為非連續的值（彈性問題中的應力即是如 此），因此研究者所得到的結果為一個不連續的應力場。針對這個問題許多有限 元素分析軟體在將最後結果呈現給使用者之前會先將不連續的參數作平滑

（smooth）的動作，使得使用者能得到連續的分析結果。在某些條件下，平滑過 後的結果會比直接分析的結果有更高的準確性也更接近真解，評估兩者之間的差 即 Recovery-based 誤差評估方法的基本概念，而如何取得比分析結果更佳的解則 是這類方法的重點。

為求得最佳應力的位置，因此有 Superconvergence 概念。元素中應力值較佳 的位置最早是由 J. Barlow (1976)[41]所提出，這些位置的應力值比有限元素法分 析結果有更高的準確性，根據 O.C. Zienkiewicz 和 J.Z. Zhu (1992)[42]兩位學者在 SPR 的建議，各元素較佳應力值的位置如圖 2 所示。

Recovery 方法有：1) Superconvergent Patch Recovery（SPR）、 2)內插方法配 合 superconvergence 的概念、 3) Simple Error Estimator for Element（SEEELEM）

及 4) Simple Error Estimator for Node（SEENode）。其中 SPR 為使用最小二乘法

（least square method）平滑區域內不連續參數的 recovery 方法；使用內插方合配 合 superconvergence 的概念，如變動式最小平方近似法 MLS 是一種使用最小二 乘法求解系統方程式以找出近似曲線或曲面方程式來求得未知特定點的資料值，

配合 superconvergence 的概念也能求出有限元素法中節點上更準確的應力值，與 SPR 不同的地方在於 MLS 不需要網格的資訊來決定影響未知點(求解點)的範圍，

MLS 並加上權重函數來決定對未知點影響的大小；而 SEEELEM 此方法為學生 的指導教授所建議的簡單的節點應力值 recover 方法。將元素內的積分點外插到 節點上取得節點上的應力值；最後，SEENode 此方法與 SEEELEM 方式近似，

為學生的指導教授所建議的簡單誤差評估方法，與先前方法不同的地方在於直接 將誤差的值計算在節點上。

經過 recovery 的計算之後更準確的應力值是在節點上，在經由 ZZ 誤差評估 (ZZ error estimato)[42]的計算，以元素為單位將應力的差值對元素作積分求得誤 差評估值，所以需利用形狀函數將元素內節點上的應力值內插到積分點上，積分 點就會同時有兩個不同的應力值，如此才可在元素上作 ZZ 誤差評估子的計算。

4 node element 8 node element,

reduce integration 8 node element

3 node element 6 node element

圖 2 各元素較佳應力的位置

2.4 METIS

METIS 為一套高效率、免費而公開之程式庫，可用來切割網格，1995 年，

Karypis and Kumar 將領域演算法予以自動化。詳細 METIS 說明可參考 METIS 手冊 4.0 版[43][44]。

三、網格演化法 

3.1 網格演化之研究方法

網格演化法此一概念乃從有限元素法中之一特殊技巧—網格調整法所演變 出來。網格調整法在處理大變形問題分析上是一個相當直觀的方法。在先前的研 究中，網格演化法的分析流程已被提出，參考圖 3之分析流程圖，此分析方法主要 的概念乃根據不同階段開挖的幾何形狀，分別產生不同的有限元素網格進行分析，並將 每一階段之分析結果，以內插的方式將計算結果(應力、應變、位移等…)分別解映射 (solution mapping)到下一階段所產生新網格元素的節點或積分點上，此種新分析模式同 等於模擬工程在時間上的連續性。如圖 1 所示：1) 初始 i=1；2) 建立第 1 階段之模型 與網格；3) 進行有限元素分析而得到第 i 階段分析結果；4) 從第 i 階段分析結果擷取其 外觀輪廓而得到第 i + 1 階段模型之初始幾何形狀並建立網格；5) 將第 i 階段的計算結

果解映射到第 i + 1 階段之網格上；6) 重覆步驟 3) - 5)直至所有步驟分析完為止。

圖 3 網格演化法之分析流程概念

圖 4 (a)為一以平面應變假設進行之深開挖分析模型示意圖及其網格切割方 式，圖中之數字表示開挖順序。目前實務上傳統的大地工程有限元素分析，將所 有不同階段的幾何形狀盡納入單一有限元素網格中，如圖 4 (a)所示，並在分析 不同階段時將對應於被挖除材料的有限元素移除，此開挖分成五個階段完成，以 傳統有限元素分析時，網格的產生須包括所有階段的幾何形狀，且必須沿著所有 階段的幾何形狀邊緣，如沿著圖 4 (a)中所有的線條產生網格分割(discretization)。

然而在進行計算時，事實上真正需要細微切割的部份為靠近開挖的網格部份，距 離開挖較遠處所需之網格分割可較為粗糙而不會影響到所得數值解的準確性。因 此本研究所提出的網格演化法乃改變此部份的分析流程而產生此新的方式。以圖 4 (a)之開挖施工順序為例，網格演化法所需考慮的開挖步驟之幾何形狀及其網格 切割步驟如圖 4 (b)~(f)所示，分別表示此開挖的五個階段，而在進行有限元素分 析時，根據步驟(b)階段的幾何形狀產生此階段的有限元素網格，之後根據傳統 的有限元素分析模式，進行數值分析，最後取得第一階段分析結果。若第一階段 結果順利求得，則此時將根據第二階段幾何形狀，先產生第二階段有限元素網格，

利用解映射的方式，將第一階段所得結果內插到第二階段的元素節點或積分點上，

如此依序求解有限元素分析、解映射，直到開挖完成或有限元素分析無法平衡為

網格演化法的計算程序上與傳統之網格調整法並無不同，然而本研究以「網 格演化法」命名此分析流程的目的在於避免與網格調整法混淆。傳統上，網格調 整法的目的在於調整網格以減少網格之過度扭曲，並增進解之精確程度；且其調 整網格之依據常為有限元素分割誤差之估計。然而，本研究所提出之網格調整法，

其調整依據為各開挖階段之幾何形狀而非網格中的誤差估計，而減少網格之過度 扭曲為網格演化之「副作用」(side effect)而非主要目的。因此，網格演化法並非 一完全新的分析流程，但是為網格調整法的一個全新應用。

圖 4 一二向度深開挖施工順序 (a)為單一網格所需通過之幾何形狀邊界；

(b)~(f)為使用多重網格所需考慮之不同階段開挖之幾何形狀邊界

3.2 網格演化之系統分析流程

本研究採用一泛用型之有限元素分析軟體

ABAQUS

詳細說明下：

由 ABAQUS 進行初始階段之有限元素分析，在此階段建立初始分析模型，

並設定相關參數，之後進行初始階段分析；分析後之結果透過 ABAQUS 產生輸 出檔(*.fil)，作為下一階段分析之初始應力場；之後由 ABAQUS 建立下一階段之 近似分析網格模型；有限元素模型和網格建立後，由 ABAQUS 產生輸出檔(*.inp)，

包 含 網 格 模 型 和 各 項 參 數 的 資 料 ； 之 後 執 行 自 行 開 發 之 模 型 調 整 函 式 庫 Modify.exe，將本階段的分析模型根據上一階段分析後之模型，進行變動；送出 變動後之網格模型予 ABAQUS 進行有限元素分析，並且將前一階段分析後之 整體應力狀態解映射到此階段各元素積分點上， ABAQUS 會主動呼叫 sigini 副程式用來對積分點上的應力值進行初始化；由 sigini 副程式呼叫自行開發之內 插函式庫 MyLib.lib；從*.fil 讀取前一階段輸出檔之應力場資訊；傳回此應力場 資訊，由 MyLib.lib 進行解映射計算；之後傳回計算後的數值送至 Sigini.lib；最 後傳回此階段有限元素模型上某積分點映射後之應力值；當所有積分點都拿取到 解映射後之應力值，則進行此階段分析，之後判斷是否有下一階段分析；如果要 進行下一階段分析，則產生此階段之輸出檔；並且由 ABAQUS 產生下一階段有 限元素分析之網格模型；如果沒有下一階段分析，則產生有限元素分析後之相關 輸出檔，並結束整體分析流程。

圖 5 網格演化法之整體分析流程

3.2.1 網格輪廓之擷取

網格演化法之分析流程中，在擷取分析模型輪廓可採用 ABAQUS 中的 Python 函式擷取，但 ABAQUS 只提供二維的輪廓擷取，而三維問題僅能得到 orphan mesh part，並無法得到模型輪廓，因此無法藉以產生下一階段的網格。在 此，本研究提出網格模型輪廓擷取之方式，主要先產生一近似的網格模型，爾後 根據前一階段分析結果與各設定對應區域內各節點變位大小，分別移動近似模型 上對應區域內之節點位置，使該近似模型與前一階段分析結果之外型一致。

以二維度四階段深開挖工程為例，每階段網格模型輪廓依所設定區域的節點，

根據前一階段同ㄧ區域內之節點變化作移動，而所選取區域如圖 6 中(a)~(d)所示，

在進行有限元素分析時，根據圖 6(a)之幾何形狀產生有限元素網格，再依傳統 有限元素分析模式，進行數值分析，最後取得 A 區域開挖之分析結果。在進行 下一階段分析時，先建立下一階段幾何形狀和其有限元素網格，並設定節點移動 之區域如圖 6(b)，之後拿取前一階段分析結果之變位數據，對下一階段之有限 元素模型做輪廓變化，之後各階段分析在網格模型輪廓拿取部份，重覆以上之步 驟。

(a) A 區域開挖之模型 (b) B 區域開挖之模型

5

1 2

3 4

C

(c) C 區域開挖之模型 (d) D 區域開挖之模型

圖 6 節點移動之區域

如此分析流程，不但保留住原網格演化法之優點，且每階段有限元素模型及 其網格為重新建立，再根據上一階段分析結果之變位數據作輪廓改變，故可避免 因拿取上一階段分析後之輪廓所產生的問題：1)下一階段開挖區域因模型輪廓已 改變，故在區域劃分 (如土壤分層) 上只能估計其位置，使整體模型建構上之精 確度較低。2)因每階段之分析，會造成模型輪廓的變動量逐漸增加，在有限元素 網格之產生時，容易造成 ABAQUS CAE 認定部份網格品質不好。

因前後兩階段設定區域上之節點數量不一致，無法將節點直接移動到前一階 段的節點位置上，故本研究針對節點移動此提出兩種方法：

(1) 將需移動之節點根據前一階段節點位置，做單一方向之調整

此方法之調整方式，為將每區域內之節點根據前一階段對應區域之節點主要 走向或主要分佈方向，延其正交方向做移動，並將每區域內端點上節點直接移動 至前一階段區域端點位置上，圖 7 為詳細之流程。以圖 8 二維空間下為例，空 心圓為前一階段節點位置，藍色圓為現階段節點移動前之位置，根據前一階段空 心圓位置之主要走向 (X 方向)，沿其正交之 Y 方向做移動，紅色圓為現階段節 點移動後之位置，而端點上的藍色點則移動到前一階段端點位置上。

圖 7 單一方向調整之流程

圖 8 節點單一方向調整方式

(2) 將前一階段之節點利用內插方法投影到相對應之幾何面上

此方法之調整方式，以三維空間做說明。1)先根據前一階段之節點利用最小 二乘法找出ㄧ趨勢平面如圖 9(a)； 2)將前一階段之節點投影至平面上如圖 9(b)，

得到屬於此趨勢面上之座標值； 3)從所有投影點中找出內積最小之兩個向量，

作為一組新的參考座標軸，並將所有節點對此一組向量作投影，使所有節點原本 為三維之座標轉換為二維之座標； 4)將現階段需移動之節點以步驟 2 的方式對 同一趨勢面做投影如圖 9(c)； 5)同步驟 3 對依據前一階段節點得到內積最小的 兩組向量作投影，使原三維座標值轉換為二維座標表示； 6)進行內插計算，將 前一階段節點轉換之二維座標作為已知點座標，而原三維座標值作為已知點的三 個數值，現階段節點轉換之二維座標作為內插點座標，去計算其內插點的三個內 插值，即為現階段節點調整後之三維座標如圖 9(d)，詳細之流程說明如圖 10。

a)趨勢面之表示 (b)前一階段節點對趨勢面做投影

(c)現階段節點對趨勢面做投影 (d)內插計算後之現階段節點位置 圖 9 以內插方法調整方式

圖 10 以內插方法做調整之流程

3.2.2 應力場之解映射

解映射主要在於分析過程中前後兩階段應力場使用內插方法做傳遞，不過在 網格模型中每個區域有著不同的材料參數設定下，如圖 11 所示，有三組不同的 材料設定(A)~(C)，根據前一階段分析後之整體應力場值，使用內插方法計算下 一階段各點位之應力值時，由於內插方法其本上假設連續狀況，而不同材料每個 區域邊界的點位可能會因鄰近不同參數設定區域之應力場影響其內插計算之應 力結果，故對此提出一解決方法。

當進行應力解映射，利用內插計算下一階段點位應力值時，將根據所要計算 的點位其所在之區域，拿取此區域前一階段之應力場值作為內插計算之參考數值。

以圖 11 為例，某一計算點位在區域(3)中，為(A)材料參數設定，故拿取前一階 段區域(3)且同為(A)材料參數設定之應力場值，計算其應力值。

圖 11 不同材料參數設定(A)~(C)之區域

3.3 程式架構

網格演化法之程式架構採用物件導向之程式語言 C⁺⁺來實作開發，自行撰寫 程式架構主要分三個主要部份，Modify.exe、Sigini.o 與 MyLib.lib。其中 Modify.exe 主要將新建立的有限元素模型，根據上一階段分析後之模型形狀作調整，將有限 元素網格模型根據前一階段模型外框之點位變化資訊，利用內插方法做計算，調 整下一階段模型輪廓上之點位資訊，已達到模型輪廓上之調整。整體核心程式架 構如圖 12 所示；而 Sigini.o 主要為提供 ABAQUS 呼叫自行撰寫的內插程式 MyLib.lib 的中繼介面；MyLib.lib 主要處理前後兩個相鄰階段應力場的傳遞，

ABAQUS 透過 Sigini.o 中的函式呼叫 MyLib.lib，MyLib.lib 讀取上一階段分析後 整體應力場值，根據 ABAQUS 所傳送過來的點位相關資訊，利用所選定的內插 方法，計算出此點位其相關應力值。MyLib.lib 中一樣包含讀檔類別以讀取有限 元素分析之輸出檔、內插類別、權重計算類別以及點搜尋類別等。其程式核心架 構如圖 13 解映射類別核心程式架構

所示。

圖 12 網格模型調整類別核心程式架構

圖 13 解映射類別核心程式架構

3.4 案例分析

3.4.1 二維度深開挖分析

如圖 14(a)所示，此二向度深開挖模型分十階段開挖，每階段開挖 3 米，各 階段開挖幾何形狀如圖 14(b)所示，擋土牆貫入深度 30 米，厚度 1.1 米，考慮之 開挖所影響之範圍為中央斷面左右各 100 米，岩盤位置假設深度為 50 米。土壤 及混凝土參數如表 1 所示，土壤組成律模式假設為完全彈塑性模式，並以莫爾

－庫倫降伏準則進行不排水分析。

根據上述開挖之模型，分別以傳統有限元素分析方法與網格演化有限元素分 析方法，進行分析，並比較圖 14(a)中擋土牆(A-A 斷面、B-B 斷面)及地表(C-C 斷面、D-D 斷面)上之變形，最後再比較兩種方法的計算時間。圖 15、圖 16 分 別為新舊方法第九階段開挖模型之位移增量分佈，由其位移之分佈可看出該開挖 之可能破壞機制與破壞面。此兩種分析方法於皆於第九階段分析步驟上產生大量 之變形而無法達到力平衡狀態，因此判斷破壞之發生，且由圖中可看出明顯之破 壞機制。另因第九階段所產生之位移量太大，若將與前八階段做比較將無法看出 前八階段之差異，故以下將使用前八階段擋土牆與地表的變形比較之。此分析所 執行的平台及選用參數如表 2 所示。

圖 17(a)、圖 17(b)、圖 17(c)、圖 17(d)分別為分析模型中 A-A、B-B、C-C、

D-D 斷面上擋土牆位移及地表之沉陷變化圖，從圖中可發現在開挖ㄧ到六階段時，

新舊方法之變形差異，幾乎都相似，但在開挖第七、八階段時，牆的變形及地表

之沉陷差異量有較明顯的的差異存在，其差異量約為 0.2 公分左右，造成差異較 大的原因可能是在此二階段應力場從上ㄧ階段解映射後，應力重新分配所造成位 移場的擾動。

(a) 一二向度深開挖分析模型

(b) 一二向度深開挖分析模型各階段開挖之幾何形狀

圖 14 二向度深開挖案例模型

200 m 100 m

50 m

1.1 m

25 m

Wall Wall

Soil

Excavate Soil

C-C D-D

50 m

1.1 m

Wall Wall

6

1 2

5

8 7

9 10

3 4

22 m 6 m

22 m

Excavate

Soil

圖 15 傳統分析方法第九階段位移增量及破壞面

圖 16 網格演化分析方法第九階段位移增量及破壞面

(a) A-A 斷面上各階段開挖擋土牆變形分佈 (b) B-B 斷面上各階段開挖擋土牆變形分佈

(c) C-C 斷面上各階段開挖地表變形分佈 (d) D-D 斷面上各階段開挖地表變形分佈 圖 17 擋土牆位移及地表之沉陷變化量

表 1 土壤及混凝土參數

Soil  Concrete 

Unit weight, γ (t/m

) 

Young’s modulus,    E (t/m2) 

5000  2.51E+6 

Poission’s ratio, υ 

cohesion, c (t/m2)

K

表 2 執行平台構架及選用參數表

CPU Intel Pentium D 3.0GHz (DualCore, 1MB L2 Cache)

Memory 1 GB*2 ( Dual Channel)

OS Linux 64bit ( Gentoo )

Compiler Intel C

Compiler 9.0 ( icpc )

Math Library Intel MKL 8.0.2

內插函式選用 MLS + 背景網格，3 階 (二維 m=10，三維 m=20)

元素型態 CPE8 (八個節點九個積分點之平面應變元素)

3.4.2 三維度深開挖分析

如圖 18 所示，此為三向度深開挖模型，此深開挖模型共分五層開挖區域，

每層之深度分別為 3 米，每層劃分為四區塊，且為模擬現地情況，每區塊之周圍 呈現一斜坡共分 20 個階段進行開挖，土壤及混凝土參數如表 1 所示，土壤組成 律模式假設為完全彈塑性模式，並以莫爾－庫倫降伏準則進行不排水分析。

在塑性材料設定下使用傳統有限元素方法和網格演化法，此兩種方法進行分 析，之後以第 20 階段開挖分析之變位增量、位移擾動量和總變位量上差異做探 討。圖 19 為兩者進行第 20 階段開挖分析之變位增量分佈；圖 20(a)中可知其擾 動量約為 0.078 公分至 0.9 公分之間；從圖 20(b)可得知最大總變位量約為 9.1 公 分，可得知最大擾動量和最大總變位量之相對誤差為 9.8%，從相對誤差值來看，

此階段之擾動量可能為兩者產生差異之原因，不過從圖中可以看到最大變位量幾 乎都在開挖底部處，而傳統分析方法和網格演化分析法上在土體挖除上有所差異，

故也可能造成兩者變位差異之原因。在圖 21 中為第 20 階段各個擋土牆之變位 增量分佈；在圖 22 中可看到各個擋土牆位移擾動量分佈，其範圍約為 0.024 公 分至 0.29 公分之間，而擋土牆之最大總變位量為 3.54 公分，其最大擾動量和最 大總變位量之相對誤差為 8.19%，故在此階段分析，其擾動量使兩者之分析結果 有所差異。

(a) 側視圖(XZ 平面)

(b)上視圖(XY 平面) (c)開挖區域之分割

圖 18 三維空間下之深開挖模型

(a)傳統分析方法 (b)網格演化方法 圖 19 第 20 階段開挖區域之變位增量

(a)擾動量 (b)總變位量 圖 20 第 20 階段開挖之擾動與總變位量

(a)傳統分析(Wall-1) (b)傳統分析(Wall-2) (c)傳統分析(Wall-3) (d)傳統分析(Wall-4)

(e)網格演化(Wall-1) (f)網格演化(Wall-2) (g)網格演化(Wall-3) (h)網格演化(Wall-4) 圖 21 第 20 階段開挖每面擋土牆之變位增量

(a) Wall-1 (b) Wall-2 (c) Wall-3 (d) Wall-4 圖 22 第 20 階段開挖每面擋土牆之擾動量

3.4.2 三維度隧道開挖分析

模型如圖 23 所示，此隧道開挖模型共分 30 個階段進行開挖，開挖區為一 直徑為 之圓形，並每階段開挖深度為 ，其土壤參數為γ = 16 kN/m3，E = 100 MPa，υ = 0.49，c =67 KPa 進行不排水分析，並假設土壤滿足 von Mises 準則。

在塑性材料參數設定下之隧道開挖模型以傳統有限元素方法和網格演化有 限元素分析方法進行分析，之後將第 30 階段開挖之變位增量、擾動量和總變位 量之分析結果做探討。在第 30 階段開挖之變位增量如圖 24 所示，由圖中可知 主要變位之發生在隧道開挖處，兩者之分佈情況有所差異，從圖 25(a)得知此階 段分析之擾動量約為 0.2 公分至 2.4 公分之間，而此階段之最大總變位量如圖 25(b)所示，約為 27 公分，而最大擾動量和最大總變位量之相對誤差約為 8.9%，

從相對誤差值可知此階段之擾動量為兩者產生差異之原因所在。

(a)上視圖(XY 平面)

60 m 6 m

3 m

(b)前視圖(XZ 平面) 圖 23 三維下之隧道開挖模型

(a)傳統分析方法 (b)網格演化方法 圖 24 第 30 階段開挖之變位增量

(a)擾動量 (b)總變位量 圖 25 第 30 階段開挖之擾動量及總變位量

3.5 效能探討

3.5.1 二維度深開挖分析之效能探討

表 3 為傳統有限元素分析 3.4.1 節模型的計算時間，開挖ㄧ到九階段之總計 算時間為 132 秒，表 4 為網格演化法各階段之計算時間，從開始到第九階段分 析結束，計算時間約為 78 秒，比使用傳統分析之方式快約 54 秒或是 40.9%。由 此可知，使用網格演化之方式應用於幾何形狀較為複雜之二維開挖模型的情況下，

會有較快的計算時間。以此案例而言，使用網格演化法的開挖模型，其各階段所 需之有限元素個數與節點個數均約為使用單一網格分析方法之ㄧ半。

表 3 二維傳統有限元素分析之計算時間

項 目 元素個數 節點個數 ABAQUS 計算時間

開挖ㄧ到九階段 4695 13828 132 sec.

表 4 二維網格演化有限元素分析之各階段計算之時間

項 目 元素個數 節點個數 ABAQUS 計算時間

放置擋土牆 2172 6381 5 sec.

開挖第一階段 1896 5847 8 sec.

開挖第二階段 1925 5950 7 sec.

開挖第三階段 1991 6148 7 sec.

開挖第四階段 1978 6109 7 sec.

開挖第五階段 2009 6208 7 sec.

開挖第六階段 2002 6191 7 sec.

開挖第七階段 1958 6063 8 sec.

開挖第八階段 1882 5843 8 sec.

開挖第九階段 1828 5685 17 sec.

總時間 78 sec.

3.5.2 三維度深開挖及隧道開挖分析之效能探討

表 5、表 6、表 7 和表 8 為傳統分析方法和網格演化法分別進行三維塑性 深開挖和隧道開挖之有限元素分析時間，在深開挖方面，從網格演化法在每階段 分析中，元素個數平均約為傳統有限元素分析之 9.1%，節點個數平均約為傳統 有限元素分析之 9.2%，在總時間花費上只需要傳統分析方法的 11.75%；而在隧 道開挖方面，元素個數平均約為傳統有限元素分析之 18.5%，節點個數平均約為 傳統有限元素分析之 19.7%，在總時間花費上只需要傳統分析方法的 35.28%。

從此數據上可明顯看出網格演化分析快速之優點。

表 5 三維深開挖之傳統有限元素分析之計算時間

項 目 元素個數 節點個數 ABAQUS 計算時間

開挖 1 到 20 階段 87626 119390 29663 sec.

表 6 三維深開挖之網格演化有限元素分析之各階段計算之時間

項 目 元素個數 節點個數 ABAQUS 計算時間

放置擋土牆 8278 10705 25

開挖第 1 階段 8379 12247 170

開挖第 2 階段 8118 11880 166

開挖第 3 階段 8174 11959 163

開挖第 4 階段 7610 11174 151

開挖第 5 階段 7974 11750 184

開挖第 6 階段 7964 11712 177

開挖第 7 階段 7923 11669 161

開挖第 8 階段 7991 11753 163

開挖第 9 階段 8721 12830 192

開挖第 10 階段 8003 11820 183

開挖第 11 階段 7709 11431 172

開挖第 12 階段 7258 10771 171

開挖第 13 階段 7098 10588 163

開挖第 14 階段 7958 11814 184

開挖第 15 階段 7601 11324 161

開挖第 16 階段 7580 11278 172

開挖第 17 階段 7493 11243 186

開挖第 18 階段 7541 11299 184

開挖第 19 階段 7379 11081 180

開挖第 20 階段 7323 10989 176

總時間 3484 sec.

表 7 三維隧道開挖之傳統有限元素分析之計算時間

項 目 元素個數 節點個數 ABAQUS 計算時間

開挖 1 到 30 階段 53797 75957 20387 sec.

表 8 三維隧道開挖之網格演化有限元素分析之各階段計算之時間

項 目 元素個數 節點個數 ABAQUS 計算時間

開挖第 1 階段 8775 12995 51

開挖第 2 階段 7887 11601 254

開挖第 3 階段 6581 9659 213

開挖第 4 階段 8138 11798 238 開挖第 5 階段 8250 12113 287

開挖第 6 階段 7349 10744 204

開挖第 7 階段 7247 10579 178

開挖第 8 階段 8262 12121 214

開挖第 9 階段 9007 13119 211 開挖第 10 階段 9388 13617 238 開挖第 11 階段 8276 12093 188

開挖第 12 階段 8865 13010 221

開挖第 13 階段 8502 12542 218

開挖第 14 階段 8356 12373 201

開挖第 15 階段 8237 13432 218 開挖第 16 階段 11016 15952 270 開挖第 17 階段 11126 16161 156

開挖第 18 階段 10324 14924 261

開挖第 19 階段 10444 15293 220

開挖第 20 階段 10868 15874 244

開挖第 21 階段 10820 15808 237 開挖第 22 階段 10518 15417 231 開挖第 23 階段 11130 16174 255

開挖第 24 階段 12905 18729 326

開挖第 25 階段 13553 19634 352

開挖第 26 階段 12075 17486 267

開挖第 27 階段 11049 16117 242 開挖第 28 階段 13577 19893 329 開挖第 29 階段 13524 19739 345

開挖第 30 階段 13870 20563 323

總時間 7192 sec.

四、網格誤差評估方法之研究 

4.1 誤差評估研究方法

本研究將針對深開挖與隧道的問題建立誤差評估值與實際誤差值的關係，其 建立的主要流程如圖 26 所示，1)決定典型分析問題：由於本研究是針對深開挖 與隧道的問題評估誤差，首先需蒐集國內深開挖與隧道的工程相關案例以了解國 內在深開挖與隧道工程的施工規模與狀況，再決定一個典型的施工規模與尺寸用 來代表分析的問題，並利用文獻回顧中的每個誤差評估方法對其評估誤差；2) 決定分析問題的收斂解：為建立每個網格密度分析結果的實際誤差值（分析結果 與真解的差值），我們必須決定有限元素法在此分析問題所能得到之最佳解（即 收斂解）以當作真解，其決定方法為觀察分析問題的特徵值在每一個網格密度的 變化情形，以開挖問題為例，特徵值可包括：最大擋土牆變位、最大擋土牆彎矩、

最大地表角變量、最大地表沉陷量、最大地中變位等，當網格達到一定密度時，

上述之特徵值將維持定量或改變量極小，此網格密度的分析結果即可作為此問題 之真解；3)誤差評估方法與特徵值：圖 26 中 3、4、5、6 步驟為一個將誤差評估 值與實際誤差值做關聯的重複性動作，在本研究中的每一個誤差評估方法所產生 的值將會與深開挖問題中的每個特徵值做關聯，如此在下一個步驟中我們可以判 斷哪一個誤差評估方法對哪一個特徵值能得到較好的誤差預測結果；4)判斷合適 的誤差評估方法：建立誤差評估值與實際誤差值的關聯後，一個特徵值會有來自 不同誤差評估方法的數個關聯結果。而理想的關聯結果資料點的分布集中使得趨 勢線可以很容易的被決定，且誤差評估值與實際誤差值的關係為接近線性並通過 原點。線性關係的優點在於每一個誤差評估值所對應到的實際誤差值較為準確，

代表此誤差評估方法在此分析問題與此位置可提供準確的誤差預測。5)建立參數 與誤差評估之關連；6)改變誤差評估方法與特徵植；7)找出合適的誤差評估方 法。

圖 26 誤差關聯建立流程圖

4.2 誤差評估程式架構

本研究之程式使用物件導向之程式語言 C⁺⁺自行開發撰寫。誤差程式之架構 主要可分為三個部分，1)有限元素資料庫程式、2)誤差評估程式、3)視覺化類別。

第一部分為用來讀取有限元素法分析結果的資訊，並將此資訊做儲存與管理；誤 差評估類別則透過公開函式成員存取第一部分儲存於記憶體之計算結果，並計算 誤差；視覺化類別則將誤差評估之結果做視覺化之處理（如 3D 呈現，動畫等）。

整體程式架構如圖 27 所示。

1.決定分析問題(開 挖、隧道)

2.找出分析問題 的收斂解

4.決定關聯特 徵值

5.建立參數與誤差 評估值之關聯

7.判斷合適的 誤差評估方法 3.決定誤差

評估方法

6.改變誤差評估 方法與特徵值

F.E.data

map<unsigned int, string*>

map<unsigned int, element*>

map<unsigned int, node*>

abaqusFormat elementLibrary

map<unsigned int, list<unsigned int>*>

set<unsigned int>

fileReader element

node

node element

C3D20R CPE6

CPE8 CPE8R

CPE4

element

C3D20 C3D10

SPR

SEEELEM myMLS

SEENode

SEENode2

ZZ

has

is

圖 27 誤差評估程式架構

4.3 誤差評估成果

以 ABAQUS 分析一平面應變問題，並在最後一步驟開挖右上四分之一的模 型大小，如圖 28 所示。在這裡選定開挖邊上的最大橫向位移與地中最大位移作 為觀察 ABAQUS 分析結果是否收斂的特徵值，如圖 30 與圖 31 所示，X 軸為元 素的大小。分析過程中所有網格的元素數量為：100、400、900、1600、3600、

6400、10000、14400、19600、25600，以元素數量為 25600 的網格的分析結果作 為此分析問題的真解，可求得每一個網格分析結果特徵值的實際誤差值（與真解 的差值），這裡將做最大誤差評估值與特徵值的關聯。

以誤差評估方法 SPR 與 SEEELEM 計算此分析問題所有網格的最大誤差評 估值，結果如圖 32 與圖 33 所示。最後將每個網格分析結果所得到的最大誤差

評估值與最大地中變位的實際誤差值做關聯，如圖 34 與圖 35 所示，為誤差評 估值與實際誤差值的關係，X 軸為實際誤差值，Y 軸為誤差評估值。可提供研究 者或工程人員對不同網格密度的實際誤差值作查詢。

圖 28 平面應變問題

(a) 元素數量:100 (b) 元素數量:400

(c) 元素數量:900 (d) 元素數量:1600

(e) 元素數量:3600 (f) 元素數量:6400

(g) 元素數量:10000 (h) 元素數量:14400

(i) 元素數量:19600 (j) 元素數量:25600 圖 29 不同網格密度下之元素相對誤差

圖 30 開挖邊上的最大橫向位移的收斂 圖 31 地中最大位移的收斂

0.207 0.208 0.209 0.21 0.211 0.212 0.213

0.01 0.1

1

log h

side displacement

0.3874 0.3875 0.3876 0.3877 0.3878 0.3879 0.388

0.01 0.1

1

log h

max earth displacement

圖 32 以 SPR 計算誤差評估值 圖 33 以 SEEELEM 計算誤差評估值

圖 34 以 SPR 預測最大地中變位的 實際誤差值

圖 35 以 SEEELEM 預測最大地中變位的 實際誤差值

五、METIS 網格切割 

本研究以一三維度隧道開挖網格為列，使用 METIS 分別分割四個及八個部 份，如圖 36 (c) (d)所示，原始網格之切割如圖 36 (a) (b)所示，切割完後之節點 及元素數量如表 9，每部的元素及節點數量約為近似，在用於平行計算上，可達 負載平衡之效能。

0.2 0.21 0.22 0.23 0.24 0.25 0.26 0.27

0.01 0.1

1

log h

error estimation(SPR)

0.12 0.13 0.14 0.15 0.16 0.17 0.18 0.19 0.2 0.21

0.01 0.1

1

log h

error estimation(SEEELEM)

0.21 0.22 0.23 0.24 0.25 0.26

0 0.0001 0.0002

0.0003 0.0004

0.0005

actual error

error estimation(SPR)

0.13 0.14 0.15 0.16 0.17 0.18 0.19 0.2 0.21

0 0.0001 0.0002

0.0003 0.0004

0.0005

actual error

error estimation(SEEELEM)

(a) 原始網格 (b) 原始網格

(c) 切割成 4 份之網格 (d) 切割成 8 份之網格 圖 36 一隧道開挖網格網格與 METIS 分割成 4 及 8 份之結果

表 9 網格切割後元素與節點數量表

4 Parts 8 Parts

Partitions Number of Elements Number of Nodes Number of Elements Number of Nodes

1 1940 2488 932 1294

2 2062 2652 1031 1408

3 1959 2510 975 1324

4 2049 2649 1031 1406

5 - - 948 1368

6 - - 1031 1357

7 - - 1031 1424

8 - - 1031 1421

六、結論 

• 本研究根據網格重置 (remeshing) 的概念發展出適用於大地工程問題之有 限元素分析程序，發展了與有限元素分析軟體 ABAQUS 結合之技術，以實 現本研究所提出之分析程序，此新的分析程序可大幅降低具複雜幾何形狀之 有限網格之建立。

• 本研究回顧了多種內插方法，對選取其中適用於不規則取樣、內插點的六種 方法：IDA、LS、MLS、PIM、RPIM-MQ 及 RPIM-EXP 進行詳細的誤差與 效能評估[1][2]，根據本研究之評估結果，RPIM-EXP 和 RPIM-MQ 在應力 計算上，其誤差和計算效能上皆有不錯的表現，建議使用於應力場傳遞之內 插計算。同時並使用物件導向技巧開發了一泛用的類別，可用以實行不同的 內插方式，並提供一致 (consistent) 的介面供其它程式使用各種不同的內插 方法。

• 對於網格演化中拿取前一階段分析後之輪廓，本研究提出兩種對幾何形狀調 整方式：單一方向調整和以內插方法做調整，並以所提出的內插方法做結合 驗證。根據評估結果，以內插方法做調整在平面對曲面之情況下效果較佳，

而在其他情況下以單一方向調整效果較佳；而無論那一種調整方式，皆需配 合一內插方法，此內插方法以 MLS 內插方法效果最好，故建議使用 MLS 為主要幾何形狀調整方法。

• 在二維深開挖案例下，本研究根據一假設案例比較舊有方法與本研究提出之 新程序比較其求得之應力場與位移場。結果發現整體趨勢完全一致，而誤差 方面由所選取的斷面上得知位移量相對差異小於 2.5%，而應力相對差異則 小於 10%。由於有限元素法計算所求得之應力位於積分點上，並以 weak formulation 以滿足力平衡條件，因此不同的網格分割將得到不同的應力值，

也因此計算結果之應力值有較大的差異是可以預期的。而與原始分析方法之 間的位移差異小於 2.5%應屬於合理的差異。

• 在三維隧道開挖案例下，網格演化方法在建構網格模型時，在現階段進行開 挖區塊，網格切割上會比較密，而其他部分則較鬆散，故在網格鬆散部分可 能產生較大之擾動量，且和傳統分析方法上也有所差異是可以預期的。

• 就執行效能上而言，網格演化法之計算效能在二維及三維分析上均優於傳統 之有限元素分析方法，且在三維開挖案例下，其效能表現更是遠優於傳統有 限元素分析(三維深開挖提升 88%、三維隧道開挖提升 64.%)，其主要原因是 網格演化法能大幅降低分析所需之自由度量，加速整體分析效率。

• 本文所提之誤差評估假設案例中使用一特徵值建立誤差評估值與實際誤差 值的關連，但在實際案例中需考慮數個特徵值，以大地工程之開挖案例來說，

需考慮最大擋土牆變位、最大擋土牆彎矩、最大地表沉陷、最大地表角變量 及最大地中變位等特徵值，因此需針對這些特徵值建立其關連性，並在此關 聯中提出最具影響的結果。

• 本研究成功的使用 METIS 之公開程式庫自動化領域之切割，並節省人力的 花費在切割網格上，且 METIS 分割的每塊區域之元素為連續，分割之後元 素量與節點量幾乎近似，使用叢集電腦於平行化分析時，可使達到較好平衡 負載節省計算時間。

參考文獻 

[1] 潘茂森(2006)，”網格演化有限元素分析於大地工程之應用” ，國立台灣科技 大學營建工程所碩士論文

[2] 孫煥期(2008)，” 三維網格演化有限元素分析於大地工程之應用” ，國立台 灣科技大學營建工程所碩士論文

[3] 宋俊宜(2007)，”有限元素網格對開挖分析結果影響之估計方法初探”，國立 台灣科技大學營建工程所碩士論文

[4] Georgiadis, M., and Anagnostopoulos, C. (1998). "Effect of berms on sheet-pile