小波理論於管結構損壞識別之應用

國立交通大學

土木工程研究所

碩士論文

小波理論於管結構損壞識別之應用

The Application of Wavelet Transform in Damage Identifications

of Pipes

研究生：謝濬鴻

指導教授：洪士林 博士

小波分析應用非破壞性檢測於管損壞識別

The Application of Wavelet Transform in Damage Identifications

of Pipes

研 究 生：謝濬鴻 Student: Chun-Hung Hsieh 指導教授：洪士林 博士 Advisor: Dr. Shih-Lin Hung

國立交通大學 土木工程研究所

碩士論文

A Thesis

Submitted to Department of Civil Engineering College of Engineering

National Chiao Tung University In Partial Fulfillment of the Requirements

for the Degree of Master of Science in Civil Engineering September 2005 HsinChu,Taiwan,Republic of China 中華民國九十四年九月

小波理論於管結構損壞識別之應用

學生：謝濬鴻 指導教授：洪士林

國立交通大學土木工程學系

摘要

近年來小波轉換應用於各個研究領域之中，由於其對訊號具有多重解析的特 性，以及提供以往傅力葉轉換只於頻率域上的解析外，訊號隨著時間變化關係的 解析能力，表現出其對訊號處理方面上的優越性。 本研究提出一個應用於管結構的非破壞性損壞識別方法，先建立識別參數資 料庫，將衝擊檢測獲得的訊號以小波分析後與資料庫進行比對，以診斷出損壞發 生位置與損壞的程度。研究過程使用數值方法模擬管勁度衰減損壞下的加速度反 應，利用小波分析的多重解析技術將加速度訊號分解成不同尺度頻率區間的子訊 號，利用特定尺度的子訊號能量變化判斷損壞發生位置及損壞程度。使用四個案 例對管損壞識別系統進行測試，顯示出除了靠近邊界時，對於損壞發生位置的識 別以及損壞程度上的識別都有相當不錯的表現。

The Application of Wavelet Transform in Damage

Identifications of Pipes

Student: Chun-Hung Hsieh Advisor: Dr. Shih-Lin Hung

Department of Civil Engineering

College Engineering

National Chaio Tung University

Abstract

This thesis proposes an application of NDT method for damage identification of pipes. A numerical model is employed to process simulate accelerative responses of pipes which are subjected to the stiffness degradation damage in circular.First,a finite-element software is employed to establish the database of identification parameters.The wavelet transform is then adopted to process the signal from impact tests performed in FEM model.Finally, the damage location and the degree of damage is diagnosed by comparing analyzed results with parameters in the database.The accelerative responses are decomposed into several different frequency regions by the multi-resolution technique of wavelet analsys.It uses the energy variation of sub-signals in specific scale to diagnose damage.Testing the identify system by four cases,it reveal that the proposed model has excellent performances in damage identification except that damage is close to two supported ends.

誌謝

經過了兩年的研究生生涯，到這裡能完成這本論文首先感謝我的指導老師洪 士林教授，提點指導我在研究生活中所遇到的挫折與困難，以及黃炯憲老師、鄭 復平老師、林昌佑老師在口試中的寶貴意見，以補足本論文的不足之處。 感謝洪門大師兄詹君志學長，在我修改與完成論文這段期間給予很多的建議 與幫忙，接著感謝多日在研究室中帶來許多歡樂的同仁，哈寶、小龜、洪安以及 學弟心農、忠錦與宏宇，讓我在生活中時常圍繞著歡樂的氣氛。 最後感謝家人的支持，尤其是父母與祖母，有了他們我才能在求學生涯中不 用顧慮生活中的瑣事，可以專心來完成自己的學業，並感謝許多在生活中鼓勵與 關心我的朋友，謝謝你們。

目錄

摘要... I ABSTRACT... II 誌謝... III 目錄... IV 表目錄 ...VII 圖目錄 ... IX 第一章 緒論 ...1 1.1 前言...1 1.2 研究動機與目的...2 1.3 文獻回顧...3 1.4 研究方法與步驟...6 1.5 論文章節與架構...7 第二章 小波理論 ...8 2.1 起源簡介...8 2.2 傅立葉轉換...9 2.2.1 短時傅立葉轉換...10 2.3 小波轉換...10

2.3.2 連續小波轉換...11

2.3.3 離散小波轉換...12

2.4 多尺度分析（MULTI-RESOLUTION ANALYSIS，MRA）...13

2.4.1 近似空間與細部空間...14 2.4.2 Mallat 演算法 ...15 第三章 有限元素法動力分析 ...17 3.1 簡介...17 3.2 有限元素與結構動力...18 3.2.1 中心差積分法...21 3.2.2 Newmark 積分法...22 3.2.3 六面體實體元素...23 3.3 有限元素分析軟體－ANSYS...25 3.3.1 軟體基本架構與分析流程...25 第四章 系統設計與案例分析 ...27 4.1 有限元素分析方法建立數值模型...27 4.1.1 實體數值模型...27 4.1.2 模型元素...28 4.1.3 衝擊試驗動力分析...29 4.2 訊號分析...30

4.2.1 訊號前處理...30 4.2.2 小波分析...30 4.3 分析案例與歸納結果...31 4.3.1 損壞發生位置識別...32 4.3.2 損壞程度識別...34 4.3.3 歸納整理...36 4.4 架構管損壞識別系統...37 4.4.1 案例分析...38 第五章 結論與建議 ...42 5.1 結論與建議...42 5.2 未來展望...43 參考文獻 ...44 附表...47 附圖...65 附錄A...127

表目錄

表3.1 ξ_j、η_j、ζ_j於六面體元素各結點的對應值...47 表4.1 數值模型基本幾何與性質...47 表4.2 量測點之編號與位置...48 表4.3 無損壞各測點之E 值表 ...48 ₅ 表4.4 單環損壞50%時損壞位置與各測點之E 值表...49 ₅ 表4.5 單環損壞40%時損壞位置與各測點之E 值表...50 ₅ 表4.6 單環損壞30%時損壞位置與各測點之E 值表...51 ₅ 表4.7 單環損壞20%時損壞位置與各測點之E 值表...52 ₅ 表4.8 單環損壞10%時損壞位置與各測點之E 值表...53 ₅ 表4.9 雙環損壞20%時損壞位置與各測點之E 值表...54 ₅ 表4.10 雙環損壞30%時損壞位置與各測點之E 值表...55 ₅ 表4.11 雙環損壞40%時損壞位置與各測點之E 值表...56 ₅ 表4.12 損壞於0.5 時不同尺度下損壞程度與Ce值關係表...57 表4.13 損壞於0.4 時不同尺度下損壞程度與Ce值關係表...58

表4.15 本研究MEYER WAVELET離散尺度與頻率對照表（當Δt=4×10-5）....59 表4.16 損壞發生於0.5 時前 35 個模態的頻率...60 表4.17 檢測案例一...61 表4.18 檢測案例二...62 表4.19 檢測案例三...63 表4.20 檢測案例四...64

圖目錄

圖2.1 短時傅立葉轉換在時間與尺度域平面上的解析度...65 圖2.2 小波轉換在時間與尺度域平面上的解析度...65 圖2.3 多尺度一維小波分解結構圖...66 圖2.4 一維小波分解演算流程圖［22］...67 圖2.5 一維小波重構演算流程圖［22］...68 圖3.1 8 結點六面體元素［21］...68 圖4.1 動力有限元素分析流程圖...69 圖4.2 管模型邊界條件與受力點...70 圖4.3 經由對稱面後簡化之管數值模型...70 圖4.4 簡化模型網格示意圖（側面）...71 圖4.5 激振力歷時圖...71 圖4.6 管模型量測點位置...72 圖4.7 訊號處理與損壞識別流程圖...73 圖4.8 量測點加速度歷時舉例圖...74 圖4.9 雙環損壞第一組對稱示意圖...74 圖4.10 雙環損壞第二組對稱示意圖...75 圖4.11 雙環損壞第三組對稱示意圖...75

圖4.13 損壞環位置0.4 損壞 50%之細部重建訊號 ...76 圖4.14 損壞環位置0.5 損壞 50%之細部重建訊號 ...77 圖4.15 損壞環位置0.6 損壞 50%之細部重建訊號 ...77 圖4.16 損壞環位置0.8 損壞 50%之細部重建訊號 ...78 圖4.17 損壞環位置0.2 損壞 40%之細部重建訊號 ...78 圖4.18 損壞環位置0.4 損壞 40%之細部重建訊號 ...79 圖4.19 損壞環位置0.5 損壞 40%之細部重建訊號 ...79 圖4.20 損壞環位置0.6 損壞 40%之細部重建訊號 ...80 圖4.21 損壞環位置0.8 損壞 40%之細部重建訊號 ...80 圖4.22 損壞環位置0.2 損壞 30%之細部重建訊號 ...81 圖4.23 損壞環位置0.4 損壞 30%之細部重建訊號 ...81 圖4.24 損壞環位置0.5 損壞 30%之細部重建訊號 ...82 圖4.25 損壞環位置0.6 損壞 30%之細部重建訊號 ...82 圖4.26 損壞環位置0.8 損壞 30%之細部重建訊號 ...83 圖4.27 損壞環位置0.2 損壞 20%之細部重建訊號 ...83 圖4.28 損壞環位置0.4 損壞 20%之細部重建訊號 ...84 圖4.29 損壞環位置0.5 損壞 20%之細部重建訊號 ...84 圖4.30 損壞環位置0.6 損壞 20%之細部重建訊號 ...85 圖4.31 損壞環位置0.8 損壞 20%之細部重建訊號 ...85

圖4.32 損壞環位置0.2 損壞 10%之細部重建訊號 ...86 圖4.33 損壞環位置0.4 損壞 10%之細部重建訊號 ...86 圖4.34 損壞環位置0.5 損壞 10%之細部重建訊號 ...87 圖4.35 損壞環位置0.6 損壞 10%之細部重建訊號 ...87 圖4.36 損壞環位置0.8 損壞 10%之細部重建訊號 ...88 圖4.37 單損壞環於不同位置損壞50%時各量測點訊號E 值...89 ₅ 圖4.38 單損壞環於不同位置損壞40%時各量測點訊號E 值...90 ₅ 圖4.39 單損壞環於不同位置損壞30%時各量測點訊號E 值...91 ₅ 圖4.40 單損壞環於不同位置損壞20%時各量測點訊號E 值...92 ₅ 圖4.41 單損壞環於不同位置損壞10%時各量測點訊號E 值...93 ₅ 圖4.42 雙損壞環於0.2 與 0.5 位置損壞 20%時之細部重建訊號 ...94 圖4.43 雙損壞環於0.5 與 0.8 位置損壞 20%時之細部重建訊號 ...94 圖4.44 雙損壞環於0.25 與 0.65 位置損壞 20%時之細部重建訊號 ...95 圖4.45 雙損壞環於0.35 與 0.75 位置損壞 20%時之細部重建訊號 ...95 圖4.46 雙損壞環於0.2 與 0.8 位置損壞 20%時之細部重建訊號 ...96 圖4.47 雙損壞環於0.2 與 0.5 位置損壞 30%時之細部重建訊號 ...96 圖4.48 雙損壞環於0.5 與 0.8 位置損壞 30%時之細部重建訊號 ...97

圖4.50 雙損壞環於0.35 與 0.75 位置損壞 30%時之細部重建訊號 ...98 圖4.51 雙損壞環於0.2 與 0.8 位置損壞 30%時之細部重建訊號 ...98 圖4.52 雙損壞環於0.2 與 0.5 位置損壞 40%時之細部重建訊號 ...99 圖4.53 雙損壞環於0.5 與 0.8 位置損壞 40%時之細部重建訊號 ...99 圖4.54 雙損壞環於0.25 與 0.65 位置損壞 40%時之細部重建訊號 ...100 圖4.55 雙損壞環於0.35 與 0.75 位置損壞 40%時之細部重建訊號 ...100 圖4.56 雙損壞環於0.2 與 0.8 位置損壞 40%時之細部重建訊號 ...101 圖4.57 雙損壞環於0.2 與 0.5 位置損壞 20%時各量測點訊號E 值...102 ₅ 圖4.58 雙損壞環於0.25 與 0.65 位置損壞 20%時各量測點訊號E 值...102 ₅ 圖4.59 雙損壞環於0.2 與 0.8 位置損壞 20%時各量測點訊號E 值...102 ₅ 圖4.60 雙損壞環於0.2 與 0.5 位置損壞 30%時各量測點訊號E 值...103 ₅ 圖4.61 雙損壞環於0.25 與 0.65 位置損壞 30%時各量測點訊號E 值...103 ₅ 圖4.62 雙損壞環於0.2 與 0.8 位置損壞 30%時各量測點訊號E 值...103 ₅ 圖4.63 雙損壞環於0.2 與 0.5 位置損壞 40%時各量測點訊號E 值...104 ₅ 圖4.64 雙損壞環於0.25 與 0.65 位置損壞 40%時各量測點訊號E 值...104 ₅ 圖4.65 雙損壞環於0.2 與 0.8 位置損壞 40%時各量測點訊號E 值...104 ₅ 圖4.66 損壞程度檢測損壞位置與施力點、量測點示意圖...105

圖4.67 損壞位置於0.5 時不同損壞程度之尺度_{2 小波轉換係數...106} 5 圖4.68 損壞位置於0.4 時不同損壞程度之尺度_{2 小波轉換係數...107} 5 圖4.69 損壞位置於0.2 時不同損壞程度之尺度_{2 小波轉換係數...108} 5 圖4.70 損壞環位置0.2 損壞 50%各測點訊號頻譜圖 ...109 圖4.71 損壞環位置0.4 損壞 50%各測點訊號頻譜圖 ...109 圖4.72 損壞環位置0.5 損壞 50%各測點訊號頻譜圖 ...110 圖4.73 損壞環位置0.6 損壞 50%各測點訊號頻譜圖 ...110 圖4.74 損壞環位置0.8 損壞 50%各測點訊號頻譜圖 ... 111 圖4.75 損壞環位置0.2 損壞 40%各測點訊號頻譜圖 ... 111 圖4.76 損壞環位置0.4 損壞 40%各測點訊號頻譜圖 ...112 圖4.77 損壞環位置0.5 損壞 40%各測點訊號頻譜圖 ...112 圖4.78 損壞環位置0.6 損壞 40%各測點訊號頻譜圖 ...113 圖4.79 損壞環位置0.8 損壞 40%各測點訊號頻譜圖 ...113 圖4.80 損壞環位置0.2 損壞 30%各測點訊號頻譜圖 ...114 圖4.81 損壞環位置0.4 損壞 30%各測點訊號頻譜圖 ...114 圖4.82 損壞環位置0.5 損壞 30%各測點訊號頻譜圖 ...115 圖4.83 損壞環位置0.6 損壞 30%各測點訊號頻譜圖 ...115 圖4.84 損壞環位置0.8 損壞 30%各測點訊號頻譜圖 ...116 圖4.85 損壞環位置0.2 損壞 20%各測點訊號頻譜圖 ...116

圖4.86 損壞環位置0.4 損壞 20%各測點訊號頻譜圖 ...117 圖4.87 損壞環位置0.5 損壞 20%各測點訊號頻譜圖 ...117 圖4.88 損壞環位置0.6 損壞 20%各測點訊號頻譜圖 ...118 圖4.89 損壞環位置0.8 損壞 20%各測點訊號頻譜圖 ...118 圖4.90 損壞環位置0.2 損壞 10%各測點訊號頻譜圖 ...119 圖4.91 損壞環位置0.4 損壞 10%各測點訊號頻譜圖 ...119 圖4.92 損壞環位置0.5 損壞 10%各測點訊號頻譜圖 ...120 圖4.93 損壞環位置0.6 損壞 10%各測點訊號頻譜圖 ...120 圖4.94 損壞環位置0.8 損壞 10%各測點訊號頻譜圖 ...121 圖4.95 無損壞管各測點訊號頻譜圖...121 圖4.96 損壞程度檢測於位置0.2 不同程度時量測點訊號頻譜圖...122 圖4.97 損壞程度檢測於位置0.4 不同程度時量測點訊號頻譜圖...122 圖4.98 損壞程度檢測於位置0.5 不同程度時量測點訊號頻譜圖...123 圖4.99 不同損壞位置於各種損壞程度之Ce 值...123 ₅ 圖4.100 管的損壞識別建構與應用流程圖...124 圖4.101 檢測案例一各量測點訊號E 值 ...125 ₅ 圖4.102 檢測案例二各量測點訊號E 值 ...125 ₅ 圖4.103 檢測案例三各量測點訊號E 值 ...126 ₅

第一章 緒論

1.1 前言

現時存在大量過去建構出的結構物，在經過歲月的洗禮後，可能已經存在許 多外來或內在的損壞。大多數的狀況其損壞程度很難由外觀判斷，為了防止突然 的破壞與防範未然，需要有方便與可靠的檢驗方法以診斷出其潛伏的損壞。 非破壞性檢測（Nondestructive Testing，NDT）為一門運用廣泛之應用科學， 其主要之概念在於不破壞受測結構物的情況下，探測隱藏於內部或潛伏於表面的 缺陷。藉由給予結構物足以反映出其物理現象的外力媒介（力、電磁、光…等）， 先量測不同情況下的物理狀態變化，接著收集所需之資料並予以研判，最後決定 結構之完整性與使用之可靠性。非破壞性檢測的目的在於系統結構的破壞識別， 破壞識別的判定可分成幾個層級（Rytter 1993、Doebling et al. 1998），第一階段 為判斷是否有損壞的發生，第二階段為判斷損壞的位置，第三階段為判斷損壞的 程度，第四階段為預測結構物能繼續使用的可能性。 非破壞性檢測方法研究裡，往往先建立數學模型來探討分析的步驟與方法， 並藉由修改數值模型做一個先行的模擬分析，以減少在真實實驗中的錯誤與浪 費，並提供了分析者較細微的角度去觀察，並可得到許多在真實實驗中無法輕易 量得的數據或現象，例如整個結構的模態，以及波動的傳遞等。 新數值分析方法的出現使得曾經研究分析過的題目得以有新的研究角度與 不同的結果討論，在不斷進步的理論發展其各種不同的應用方式，如小波理論的 出現讓人們以不同的方式觀察訊號的特徵，對於訊號、影像處理上有著優良的表 現，在非破壞檢測的領域中，藉由小波理論提供了訊號處理的新工具以發展出更 方便與準確的檢測方法。

1.2 研究動機與目的

運輸管主要是負責流體的運送以節省運輸資源的浪費，對於其長期暴露在環 境當中以及遭受運送物質的侵蝕與壓力，可能造成某部份的細小損壞並逐步演變 成較嚴重的損壞，甚至造成整個運輸管線的停擺，為了避免運輸管線發生突發性 的破壞狀況，需要一個非破壞性的檢驗方法去檢測管的損壞狀況，來確保其使用 上的可靠性。 非破壞性檢測的檢測方式有許多樣的方法來達成，如利用受測物的靜力行為 上的形狀曲線、動力行為的速度、加速度歷時反應與阻抗行為如壓電反應等，結 構物的動力反應行為，包含著結構物本身的性質，諸如幾何形狀、材料性質等， 在力學的行為上，局部勁度衰減的損壞管結構物在受靜力行為下的行為變形不是 相當的明顯，則對於局部損壞所造成變形曲線不連續部份之偵測效果較為不彰。 動力反應加速度訊號是損壞識別重要且可行的媒介之ㄧ，藉由分析加速度訊號在 不同損壞位置與損壞程度的差異性作為損壞識別的判斷依據，而動力行為下所產 生的加速度訊號可藉由不同的數值分析工具以找出其中的差異性，例如將訊號轉 至頻率域上已觀察訊號頻率特性的傅力葉轉換，以及近年來的小波理論分析，其 在訊號處理上的強大功能，如訊號壓縮與濾除雜訊、偵測訊號於頻率或二次以上 不連續處，和它可以同時觀察訊號在頻率域與時間域上特性的能力。 本研究是探討管結構體，在其發生損壞時在動力行為上的差異性，利用數值 模型模擬局部環狀的勁度衰減損壞，並藉由有限元素動力分析模擬此受損管結構 物在受一衝擊激振力後的動力加速度反應，接著利用小波分析的多重解析技術來 解析此加速度訊號，以其對動力訊號進行去除雜訊以及分解訊號的能力，將訊號 分解成不同的頻率區間以檢視在不同損壞情況之下訊號的組成成分有何差異 性，藉以觀察損壞發生於不同位置與程度時，所量測訊號分析後結果之間的關

較佳之檢測方法與小波分析方法，找出較簡易的損壞發生識別參數以及損壞程度 識別參數，並建立管結構損壞的識別方法以建構整個管損壞識別的程序與一套可 行的檢測方式。

1.3 文獻回顧

對於文獻的收集與參考主要以兩個方向去整理，首先就是對於非破壞性檢測 的方法，從中學習非破壞性檢測中數值模型的模擬方法與所分析量測值之間的可 行性與優缺點，再來就是訊號的處理方法，對於現今訊號處理的方法中學習有哪 些較為合適於模型訊號的分析方法。而近年來小波轉換（Wavelet Transform， WT）的研究發展，提供了對訊號處理的另一種工具，由於其對於時間域上的敏 感性，在工程的應用上時常被拿來偵測訊號事件的發生，如機械震動，齒輪的局 部承載缺失，梁、板的裂縫位置判定等，以及利用訊號組成成分的分解來去除雜 訊、訊號壓縮與純頻率分解等，所以對於訊號處理與分析文獻則傾向於小波理論 之應用。 楊黎熙（1994）使用數值方法模擬振動訊號，以梁結構的自然頻率與加速度 值極值並配合類神經網路學習以達到梁結構裂縫的位置與深度進行檢測［1］。 Z. Sun 與 C.C Chang（2002）使用小波包分析（Wavelet Packets method）以滿足 小波分析於高頻多解析能力的不足之處，將結構物之振動行為以小波分解之各個 階層的能量並輔以類神經網路學習來衡量結構物的破壞［2］。H. Melhem and H. Kim 比較快速傅立葉轉換與連續小波轉換於道路結構體損壞的偵測，其研究指出 使用連續小波轉換後對於偵測道路結構上的梁與板裂縫皆有很好的分析能力［3 ］。Y.J. Yan，H. N. Hao 與 L.H. Yam（2004）以自然頻率與使用小波理論於動力 反應的結果在複合板結構的裂縫檢測，其結果顯示出小波轉換的能量譜值，比使 用自然頻率在小破壞的偵測上更具敏感性［4］。Quan Wang 與 Xiaomin Deng （1998）使用空間小波（Spatial Wavelets）分析來偵測梁與板裂縫的位置，即以

小波轉換分析結構物的模態形狀，藉由小波分析偵測訊號不連續性的有著敏感反 應特性，以顯露出空間中形狀的不連續處，此研究中指出，對於較低的訊號解析 度即測點的取樣頻率而言也有不錯的表現，此方法也不需藉由與健康結構物的資 料進行比對，純粹從損壞結構物所得之模態形狀即可進行分析［5］。以針對訊號 處理與轉換的文獻裡S.K. Tnag（2000）使用短時傅立葉轉換與諧和小波分析振 動訊號的衰減行為，研究指出對於單頻訊號時對於萃取衰減的指數常數兩者皆有 不錯的表現，而多頻率訊號時短時傅力業轉換則有較佳的表現［6］。Z. Su 與 L. Ye （2004）發展出一套的方法，其為將感應網路所得之訊號經連續小波轉換後，經 由訊號過濾與壓縮、特徵擷取與再用類神經網路來做片段識別，以達到損壞的評 估［7］。Biswajit Basu 利用小波轉換並計算其係數能量來識別多自由度系統當部 分系統勁度減弱時之動力行為，並且指出其也可應用於對損壞時間發生點的偵測 ［8］。J.N. Watson 與 P.S. Addison（2002）使用小波轉換於非破壞性檢測裡的超 音波回聲訊號，先分析出時間－頻率域平面上係數最大值的線，保留其係數並與 以分析［9］。C. Junsheng，Y. Dejie 與 Y. Yu（2005）使用小波轉換振動訊號後， 計算不同頻率波段的能量於時間上的密度分佈來判定滾輪受載時的缺陷［10］。 Ｂ. Liu(2003)的研究顯示出使用適當的諧和小波比小波包方法有更加的頻率解析 能力［11］。S. Loutridis(2004)使用小波的局部能量密度以偵測齒輪系統的缺陷。 ［12］Y. Zou，L. Tong and G. P. Steven(2004)使用有限元素模型分析夾層板得出 板的自然頻率與動力反應，其研究指出分解後小波訊號的能量譜對於小裂縫有不 錯的測能力［13］。蘇威智(2003)使用小波包分析對離散化的結構系統進行系統 識別，進而求出系統的自然頻率、阻尼比等資訊。［14］

模型方法（Model-Based Method）一般是以有限元素分析來實行結構分析， 藉由修改模型的材料性質或外在幾何條件來模擬結構破壞，並可由各種分析得到 結構物不同的特徵反應，如模態分析（Modal Analysis Method）、時間反應、頻

性質（質量、阻尼、勁度）的方程式所得，而結構物的破壞正是這些物理性質的 改變，一般而言損壞通常會導致質量與勁度的降低，並局部的增加阻尼比，就以 上三種性質而言，結構損壞時阻尼的改變最為敏感，其次是勁度改變再來才是質 量的改變。而由於不同的偵測技術，模態分析方法主要能分為下列幾種，如模態 形狀改變法（Modal Shape Change Method）、模態形狀曲線法（Modal Shape Curve Method）、應用敏感性更新法（Sensitivity-Based Update Method）、最佳矩陣更新 法（Optimal Matrix Update Method）、頻率反應方程式法（Frequency Response Function method）、綜合模態參數法（Combined Modal Parameters Method）等。 以上方法由於是採用模態中較低頻的部份，故顯現出的是結構物整體行為，也因 如此在量測點的選擇上也比較有彈性。然而對於模態方法也有些缺點，首先是其 只能偵測出固定形式的破壞，以及由於必須要有一個健康完整的結構物資訊做比 對，所以通常在複雜結構物中需要保留下龐大的資料，再來它對於局部的小損壞 較難偵測得出來，而且對於局部損壞的位置也很難加以判定。［15］ 而完整的研究方法中訊號處理也是相當重要的一環，如何消除訊號中的雜訊 與突顯夾雜在訊號裡的特徵，是訊號處理的主要工作。由於損壞檢測的對象多是 基本的結構要件，所以在土木結構研究上也多以梁、板等基本結構物件來做探 討，進而至複合板、深梁，橋梁等較複雜的結構物做研究。而在模型模擬後所需 分析的量值，也因量測的物理量與方法而有所不同，所擷取物理量如位移、速度、 加速度、電壓等改變，或是結構模態形狀，阻抗反應等。在訊號處理中，多是將 資料轉至不同的值域來觀察與比較訊號的性質和特徵。對於訊號在頻率域 （Frequency Domain）上做分析，此類的方法是基於結構物勁度的減少將會導致 自然頻率的降低，藉此只以結構物的頻率反應而診斷損壞，由於局部性的損壞將 導致在某固定頻率上的改變，但是對於小破壞而言，在局部的區域上偵測高頻會 有較好的辨識度，但如此一來卻又對於判定損壞的位置顯得不切實際，也有研究 指出回聲的頻率在辨識上有著較好的表現。對於訊號的時間域（Time Domain）

而言，基本上所有的方法都是來自時間的歷時資料做分析，藉由輸入的外力與震 動反應來衡量破壞程度與破壞位置，在時間域內對於事件發生之時間點較具敏感 性。而阻抗域（Impedance Domain）上的分析是基於結構物的每一個部份，都將 提供結構物阻抗，所以當結構物的部份性質改變時，其阻抗也跟著改變，在力學 上阻抗可定義成 v F ，F 為結構物某點所受之外力，而v則為此點的速度，而電學 上 則 是 電 壓 的 比 率 ， 此 類 方 法 很 適 合 應 用 於 即 時 監 測 （Real-Time Health Monitoring）。［15］

1.4 研究方法與步驟

1. 資料蒐集 蒐集相關文獻資料如論文、期刊等，尋找有關於非破壞檢測的數值模擬方 法，以及瞭解小波理論的優點特性與其應用於損壞識別的研究狀況與成果， 從中學習與思考適用的模擬與分析方法來建構整個損壞識別系統。學習研究 時所需之工具，如分析軟體工具的合適性與使用方法。 2. 建立分析資料 建立管結構的數值模型，模擬不同損壞位置與損懷程度的管結構案例，並逐 步修正不同的檢測方式以得到小波分析時所需較合適的動力訊號。 3. 尋找管結構的損壞識別參數 利用小波分析於不同尺度下分解或轉換有限動力分析所得之訊號，觀察並找 出某一尺度對於損壞發生與程度較具有關聯性，並以識別參數定量的呈現此 種現象。 4. 建構損壞識別系統 整合數值模型案例所分析之結果，架構整個損壞識別流程與方法，並以其分

5. 研究成果討論與論文撰寫 彙整本研究之成果結論與相關理論撰寫成論文。

1.5 論文章節與架構

本論文的研究架構分為五個章節，第一章為緒論，其內容說明了研究背景、 研究動機與目的、文獻回顧以及研究方法與步驟。第二章為小波理論的理論基 礎，介紹小波理論的發展、特性與適用性，以及理論的實質內容如小波轉換的種 類與基本公式、多尺度分析的介紹。第三章為動力有限元素分析的理論，介紹有 限元素法的發展與數值運算上的公式，最後介紹有限程式分析的步驟。第四章為 系統設計與案例分析，在此以有限元素模型建立不同損壞模式的管結構，利用小 波分析其動力訊號與資料整理後，找出識別損壞參數與識別損壞的方法，接著建 立整個管的損壞識別系統，第五章為結論與建議，本章將就本研究所做出之方法 做一個結論並分析其不足待改善與精進之處。

第二章 小波理論

2.1 起源簡介

1910 年 Alfred Haar 提出了最早的小波單範正交基底，但當時並沒有出現「小 波」這名詞。Y. Meyer 認為，小波分析思想萌芽於 1930 年至 1980 年，但真正最 初研究小波是法國地質物理學家Jean Morlet。20 世紀 60 年代，由於工業發展需 要，尋找地下石油成為法國一項重要目的，如何從反射信號中提取有用的石油信 息是Morlet 的主要工作，1981 年 Morlet 仔細研究 Gabor 變換的方法，對 Fourier 變換與加窗 Fourier 變換的異同、特點及函數構造做了創造性研究，首次提出了 「小波分析」概念，建立以他名字命名的 Morlet 小波。小波分析的方法主要是 由Y. Meyer 與其同僚所發展，主要的演算法可追朔到 1988 由 Stephane Mallat 所 提出，從此以後小波的研究工作就推展到國際上。諸如Ingrid Daubechies、Ronald Coifman、Chui 和 Victor Wickerhauser 等人的研究工作聯合奠定了小波分析的基 礎。［16］ 傳統的信號分析是建立在傅立葉轉換（Fourier Transform）的基礎之上，對 於許多訊號而言，傅立葉分析是一種很有用的工具，因為訊號的頻率包含了許多 重要的訊息，但由於其轉換是將訊號轉至頻率域，是一種全域的轉換，在時間局 部域的資訊卻被遺漏，然而對於訊號中短暫出現的特徵，傅立葉分析並不是相當 適合去偵測它們的存在，而這種特徵剛好是非平穩信號最根本和最關鍵的性質。 為了改善以上的缺點，Dennis Gabor（1946）在傅立葉轉換上採用了窗口（window） 的觀念，稱之為短時程傅立葉轉換（Short-Time Fourier Transform，STFT），同時 支援了訊號在時間與頻率上的看法，並提供了一些訊號事件發生的時間與在哪一 個頻率的資訊，但其資訊的精度受限於窗口的大小，而且一但決定了窗口的大 小，此窗口將固定的在所有頻率都將被使用。為了在訊號分析上有著更加彈性的

低頻資訊時使用長時間的區間，與高頻資訊時使用短時間的區間，得到更高的精 度。近年來，已有學者提出使用小波當作系統識別的工具，而以小波架構的解析 工具能同時具有良好的時間域以及頻率域解析能力，非常適合探測瞬態反常現 象。小波分析是傅立葉分析的重要發展，它不僅保留了傅氏理論的優點，又能克 服其不足之處。小波分析尚擁有許多優良的特性使其不僅應用在單一領域中，其 涵蓋的應用範圍，諸如：數學、物理、訊號處理、影像處理、地震勘探、海洋工 程、機械故障診斷與及時監控…等領域。［16］

2.2 傅立葉轉換

訊號的轉換主要是在使用不同的基底下轉換，傅立葉轉換是以三角函數 cos(x) 、sin(x) 或指數函數_{e 為基底，將一訊號由時間域轉換到頻率域下，故為}ix 了探討訊號之頻率特性，傅立葉轉換是相當有用的分析工具，其轉換式表示如下 式

∫

− = f x e dx f f iωx π ω ( ) 2 1 ) ( （2.2.1） 而其逆轉換為 dx e f x f f i x

∫

= _ω ω π ( ) 2 1 ) ( （2.2.2） 對於處理一離散時間點之訊號，傅立葉轉換與逆轉換式則修改如（2.2.1）式 與（2.2.2）式，將訊號轉至離散之頻率域，稱之為離散傅立葉轉換（Discrete Fourier transform，DFT），若訊號向量長度為N，則

∑

− = − = 1 0 2 ) ( ) ( N n nk N i f _{k f n e} f π ，k =0,1,...,N −1 （2.2.3）

∑

− = = 1 0 2 ) ( 1 ) ( N k nk N i f _{k e} f N n f π ，n=0,1,...,N −1 （2.2.4）

2.2.1 短時傅立葉轉換

傅立葉轉換後從頻率域上可知有暫態訊號產生，但卻不知在何時發生， Dennis Gabor(1946)將視窗(windowing)的觀念引進傅立葉轉換中，把信號劃分成 許多的時間間隔，用傅立葉轉換分析每一個時間區間，以確定該時間區間內存在 的頻率，改善了傅立葉分析無法偵測時間區段內頻率特性的不足。取傅立葉轉換 前先乘上一個窗口函數，則至少可知暫態發生在哪一個窗口範圍內，設此窗口函 數 為 g(x) ， 則 此 種 傅 立 葉 轉 換 稱 為 短 時 傅 立 葉 轉 換(Short-Time Fourier Transform，STFT)，可寫成 dx x e g x f f f_(ω,τ)=

∫

_{( )} ∗ −iωx_{( −τ) （2.2.5）} ) (x−τ g 為窗口移動函數，τ 則為平移參數，當窗口函數g(x)為高斯函數時，STFT 則稱之為Gabor 轉換。STFT 有固定持續時間的視口函數，以及一個固定的頻率 解析度，在時頻域平面上的解析度可由【圖2.1】所示。但由於 STFT 的視窗大 小為固定，有一些需要更佳解析度時，並無法進行調整，對於較細微的部分訊息 被忽略，造成訊息模糊不清的問題。而小波轉換理論則是更進一步的將的利用函 式的縮放改變窗口大小，具備了STFT 的優點更提升了其解析度。

2.3 小波轉換

小波轉換則是提供一個不一樣的基底，使的訊號在某些基底下更容易的展現 出其特質。其基底是由一個原型(prototype)函數的縮張(dilation)與平移(translation) 所形成的。這些基底函數具有短持續時間且高頻率，以及長持續時間且低頻率的 特性，因此相當適合表達高頻的突發暫態訊號或是長時間變化緩慢的訊號。正是 這種特性，使小波轉換具有對信號的自適應性［16］。小波轉換在時間與尺度域 平面上的解析度如【圖2.2】所示。

2.3.1 小波函數

對於一個有限函蓋函數ψ(x)，符合以下三點： 其全域積分值為零，即 （1）

∫

ψ(x)dx=0 且滿足 ψ(x) =1，也就是 （2）

_∫

ψ(x)2dx=1 （3）

{

ψ(2jx-k)| j,k∈Ζ

}

_形成L2_{(R)的一組基底} 由於（1）式，此函數ψ x 的圖形必定上下起伏，且由（2）式可知( ) ψ(x)下 所圍之面積亦不會很大，故名之小波或凌波。例如Harr 函數就是這樣一個函數， 且為眾多小波函數族類中最基本的一個，亦稱之為一階小波函數，而所謂的小波 函數族類包括了，正交、雙正交、半正交，區間中的小波函數等等，如Daubechies wavelets、Coiflets、Symlets、Biorthogonal wavelets，這些族類內的函數隨著階數 而變化。 在做小波轉換時，首先必須挑選一個小波函數ψ x 一般稱之為母小波( ) （mother waverlet），而母小波經由尺度的縮張與平移後之函數 ( - ) a b x ψ 一般稱之 為子小波，其中a稱為尺度參數，b稱為平移參數。

2.3.2 連續小波轉換

假設訊號為平方可積(即能量有限的信號空間表示為 _f(_x)_∈L2(_R)_{)，則所謂} 的連續小波轉換（Continuous Wavelet Transform，CWT），即將一個單變數函數

) (x

dx a b x ψ x f a b a f w( , )= 1

∫

( ) ( − ) _{a,b∈ℜ，a≠0 （2.3.1）} 而其逆轉換為 dadb a x-b ψ b a f a C x f ψ ) ( ) , ( 1 1 ) ( R w 2

∫∫

₊ = （2.3.2） 其中C 是一個隨著_ψ ψ 而變的常數 ω ω ω d ψ C_ψ =

∫

2 ) ( （2.3.3） ) (ω ψ 為小波函數的傅立葉轉換 dx x ψ ψ _{( )e}-iwx 2 1 ) ( =

∫

π ω （2.3.4）

2.3.3 離散小波轉換

離散小波轉換（Discrete Wavelet Transform，DWT）其轉換式類似於連續小 波轉換一樣，不同處在於尺度參數a與平移參數b為離散值，如此可以避免在連 續小波轉換時產生的多餘值（Redundancy），其轉換式如下 dx a b x ψ x f a k j f w( , )= 1

∫

( ) _j,k( − ) _a=cj_,b=kcj_;(j,k)∈Z2 _（2.3.5） 通常c值取2，則尺度a參數為2 的冪次方，即 20、21、22…，如此一來在做轉換 時能更加的有效率。 其中，ψ 必須為單範正交，其必需滿足 _j,k ⎩ ⎨ ⎧ = = =

∫

_j,k t _j'_,k' t dt if_elsej j' k k' 0 , 1 ) ( ) ( ψ ψ （2.3.6） 這表示在同一尺度 j ，k下，它都是單範正交的。

以上不論是連續小波轉換或者是離散小波轉換，其原始訊號在時間域上均為 連續的，對於一般所量測的訊號或者是數位訊號而言，其在時間域上是不連續， 本研究所處理的訊號就是屬於這種離散訊號，遇到這種離散訊號時，則將式 （2.3.1）與式（2.3.5）的積分式換為總和 ) ( ) ( 1 ) , ( w ₌

∑

− a b x ψ x f a b a f （2.3.7）

2.4 多尺度分析（Multi-resolution Analysis，MRA）

1988 年 S. Mallat 在發展正交小波基底時提出了多尺度分析(Multi-resolution Analysis，MRA)的概念，以空間的概念上抽象地說明了小波多分辨率特性，將 此之前的所有正交小波基底的構造法統一起來，並發展出快速小波轉換的演算法 （Fast Wavelet Transform Algorithm），亦即所謂的 Mallat 演算法。

在提到多尺度分析之前，先介紹多層解析空間（Multi-resolution space）的概 念，若一個擁有無限解析度的空間V ，其所包含的任意V 子空間與_j V 子空間，_j+1 1 j+ V 的維度比V 高且_j V 包含於_j V_j+1（即V_j ⊂ V_j+1），當 j 無限大時即為空間V ， 也就是說整個空間V 是一層套著一層的多層空間，如此一來當 j 越來越大時，V_j 中的函數 f_j(x)就越能逼近V 的函數 f(x)，則空間族

{ }

V 為多層解析空間［17］。 _j 多尺度分析的應用，主要是為了表現出訊號在某一尺度的有著明顯的特徵， 多尺度分析是將訊號 f(x)層層分解為不同尺度的分量，由不同尺度的分量進而 觀察原始訊號的特性藉此進行分析。其多尺度一維小波分解結構的組織形式如 【圖2.3】所示，在圖中 cA 為經分解後的近似係數，cD 為分解後的細節係數， 而下標數字則表示分解的尺度階數。

2.4.1 近似空間與細部空間

近似空間V 與細部空間_m W 之間關係如下： _m j j j j j V W V W V ₊₁ = ⊕ , ⊥ （2.4.1） 其中符號⊕稱為兩空間V 與_m W 的正交直和(Orthogonal directsummation)。 _m 1 + j V ：較高解析度的近似訊號函數空間 j V ：較低解析度的近似訊號函數空間 j W ：V 與_j+1 V 相對應近似訊號函數之細節訊號函數空間 _j 在訊號處理中，一個較高解析度和一個較低解析度的訊號，其差別在於較高 解析度的訊號比較低解析度的訊號含有更多的資訊，即所謂的細節部份。這些細 節部份則是落在W 內。將式（2.4.1）展開可寫成 _m 0 0 j j j W W W V V = ₋₁⊕ ₋₂ ⊕_L⊕ ⊕ （2.4.2） 當j→∞時，所有W 的正交直和即為_{j L}2(_R)_。 j V 是層層相套的空間，全是由 一個尺度函數（scaling function）φ(x)的縮張平移所造成；W 是兩兩正交的空間，_j 它是V 中，層與層之間的細節部分，全是由一個小波函數_j ψ(x)的縮脹平移所造 成，也就是說φ(x)與ψ_jk(x)為_L2(_R)_{的一組基底。［17］} 1} -2 k ,0 0 ) ( { ) ( j jk x j x ∪ ψ ≤ <∞ ≤ < φ （2.4.4） 若將（2.4.2）是推廣至j< ，則 0

L L⊕ ⊕ ⊕ ⊕ ⊕ = _{m− m− 0 -1} m W W W W V _{1 2} （2.4.5） ) ( 2 _R L 的一組基底（2.4.4）式即可改寫成 Z k Z, j | ) x ( jk ∈ ∈ ψ （2.4.6） 其中 ) -(2 2 ) (x j jx k jk ψ ψ = （2.4.7）

2.4.2 Mallat 演算法

對於一離散的訊號，Mallat 演算法應用在快速離散小波分解（Decomposing） 與重構(Reconstruction)中，是類似於訊號處理中將訊號分成兩個頻率波段一樣， 可將訊號一階一階進行分解成近似係數與細節係數，即所謂的高頻與低頻部份， 以及藉由這些近似係數與細節係數還原成近似訊號以至原始訊號，在分解演算時 首先將原始訊號分解成近似係數和細節係數兩部分，再往下一階分解時則是提取 近似部份做分解，一樣得到近似係數和細節係數兩部分，以此類推，其演算流程 如【圖2.4】。分解過程是將訊號經由濾波器進行分解成兩相正交的低頻部分和高 頻部份，訊號經過低通濾波器(low-pass filters)分解成低頻部份，經過高通濾波器 (high-pass filters) 分解成高頻部分，低通濾波器由尺度係數所組合，高通濾波器 則是由小波係數所組合。每一次分解後的數據量減半，因此分析後得到的低頻成 分和高頻成分的時域解析度比分解前訊號減低一半。重構演算則是分解演算的逆 過程，經由每一層重構之後，訊號的數據量增加一倍，因此可提高訊號時域解析 度，重構演算之流程如【圖2.5】。 【圖2.4】與【圖 2.5】中，符號↓ 是做降取樣(downsampling)，符號 22 ↑ 是 做升取樣(upsampling)；降取樣是將分解後的訊號用每兩個點取一點的方式，重 新排序，其觀念同跳取，升取樣是在前面做次取樣後每兩個點之間插入一點零， 再重新排序，其觀念同零內插，舉例說明如下：

假設有一訊號序列為［0 1 2 3 4 5 6 7 8］ 跳取：

若跳取因子為2，則新的序列為［0 2 4 6 8］。 零內插(Zero interpolation)：

第三章 有限元素法動力分析

3.1 簡介

有限元素法（finite element method）是用來解物理與工程微分方程的一種解 析方法，源於1950 年代早期的太空工業，首先由 Turner、Clough、Martin 與 Topp 先後提出發表，而後其他研究人員如Szabo 與 Lee(1969)，Zienkiewice(1971)等先 後提出有關結構力學、熱傳遞與流體力學的元素方程式，亦能用加權剩餘程序 （weighted residual procedure）如 Galerkin method 或最小平方法來導出，擴大了 有限元素法的應用範圍，也因此使得有限元素法能應用到任何微分方程式。［18 ］ 此方法的基本觀念是將連體視為次區域，作為結構解析之模擬，任何連續 量，如：溫度、壓力、位移等，均可用一不連續函數的形式做為近似的表示。此 型式為有限區域的集合分段連續函數所組成，使用連續量的值，以定義分段連續 函數在其有限數次域。［19］ 將連續的求解域離散為一組有限元素的組合體，這樣的組合體能解析的模擬 或逼近求解區域。有限元素作為一種數值分析方法的另一重要步驟，是利用在每 一個元素內假設的近似函數來表示全求解區域上待求的未知場函數。元素內的近 似函數通常由未知場函數在各個元素結點上的數值以及形狀函數表達。如此一 來，一個問題的有限元素分析中，未知場函數的結點值就成為新的未知量，從而 使一個連續的無限自由度問題變成離散的有限自由度問題。一經求解出這些未知 量，就可以利用形狀函數確定元素組合體上的場函數。顯然，在理論上，隨著元 素數目的增加，亦即元素尺寸的縮小，若元素是滿足收斂要求的，近似解最後將 收斂於精確解。［20］

有限元素法在解題的過程需要儲存極多的數值資料及進行大量的計算，隨著 個人電腦的運算能力提昇，使用有限元素分析進行研究的代價越來越低，並促使 有限元素法的應用領域越來越大，不論是力學、電磁學、聲學等方面研究，有限 元素法都扮演著極重要角色。

3.2 有限元素與結構動力

由於在結構動力分析時，結構在隨時間變化的載重之下，位移、速度、加速 度，以至於應力、應變都是時間的函數。由於引入了時間變數，故所處理的是四 維的問題

(

x,y,z,t

)

，在有限元素法當中一般採用部分離散的方法，即只對空間進 行離散，則元素內的位移u，v，w的形狀表示式 ) ( ) , , ( ) , , , ( _i 1 i i x y z u t N t z y x u

∑

n = = ) ( ) , , ( ) , , , ( _i 1 i i t v z y x N t z y x v

∑

n = = （3.2.1） ) ( ) , , ( ) , , , ( _i 1 i i x y z w t N t z y x w

∑

n = = 亦可寫成矩陣形式 e NX u = （3.2.2） 其中 ⎪ ⎭ ⎪ ⎬ ⎫ ⎪ ⎩ ⎪ ⎨ ⎧ = ) , , , ( ) , , , ( ) , , , ( u t z y x w t z y x v t z y x u

[

N N Nn

]

N = _{1 2 L} ， N_i = N_iI_3×3 i =1,2,3_L,n

⎪ ⎪ ⎪ ⎪ ⎭ ⎪⎪ ⎪ ⎪ ⎬ ⎫ ⎪ ⎪ ⎪ ⎪ ⎩ ⎪⎪ ⎪ ⎪ ⎨ ⎧ = n 2 1 e X X X X M ， ⎪ ⎪ ⎭ ⎪ ⎪ ⎬ ⎫ ⎪ ⎪ ⎩ ⎪ ⎪ ⎨ ⎧ = ) ( ) ( ) ( X i i i i t w t v t u i=1,2,3_Ln 則 元 素 的 應 變 與 應 力 可 藉 由 所 得 之 結 點 位 移 u ， 與 應 變 － 位 移 矩 陣 （strain-displacement matrix）B 及彈性矩陣（elasticity matrix）D 相乘而得

e BX = ε （3.2.3） e DBX D = = ε σ （3.2.4） 為了建立結構動力學方程式，可利用Lagrange 方程式 0 X R X L -X L ₌ ∂ ∂ + ∂ ∂ ∂ ∂ & & dt d (3.2.5) 式（3.2.5）中，L 是 Lagrange 函數，R 為耗散函數 P -T L= ∏ （3.2.6） T 為系統之動能，∏_P是系統的位能

∫

= Ve T e _{u u} 2 1 T ρ_{& &}dV （3.2.7） dS dV dV - u f - u T D 2 1 S T Ve T Ve T e p =

∫

∫

∫

∏ σ ε ε （3.2.8） （3.2.7）式中ρ代表元素的密度，（3.2.8）式中第一項為元素的應變能，二、 三項代表外力的位能，f 是作用於元素的體積力，T 則是作用於元素的面積力。 至於元素的耗散係數R，如果元素的阻尼力為各質點速度的比例時，可表示為

∫

= Ve T e _{u u} 2 1 R μ_{& &}dV （3.2.9） 其中μ代表阻尼係數，將（3.2.2）式引入（3.2.7）、（3.2.8）和（3.2.9）中，則得 到

e e T e e _{(X ) M X} 2 1 T = & & e T e e e T e e p (X ) K X -(X ) Q 2 1 = ∏ （3.2.10） e e T e e _{(X ) C X} 2 1 R = & & 其中為元素的質量矩陣：M N NdV Ve T e

∫

= ρ 勁度矩陣：K B DBdV Ve T e ₌

∫

阻尼矩陣：C N NdV Ve T e

∫

= μ 元素荷載向量：Q N fdV N TdS S T Ve T e ₌

∫

₊

∫

σ 集合各個元素的矩陣和向量，以形成整個系統的運動方程式 X M X 2 1 T T T e e _{= & &} =

∑

Q X -KX X 2 1 T T e e p p = ∏ = ∏

∑

（3.2.11） X C X 2 1 R R T e e _{= & &} =

∑

M、K、C、R 分別是系統的質量矩陣、勁度矩陣、阻尼矩陣、元素荷載向量， 將（3.2.11）式代入 Lagrange 方程式（3.2.5）即可得到系統的運動方程式 ) Q( ) KX( ) ( X C ) ( X M&& t + & t + t = t （3.2.12） 此為一個二階常微分方程組，在有限元素動力分析中，一般採用直接積分法 或者是振型疊加法，直接積分法是對運動方程式直接進行積分，而振型疊加法則 是先求解一無阻尼系統之自由震動，即（3.2.12）式中忽略則尼 C，與等式右側 之荷載力為零時所得之式

0 ) KX( ) ( X M&& t + t = （3.2.13） 這是屬於一個矩陣特徵植問題，利用所得的特徵向量及系統之振型對方程式 （3.2.12）進行變換，最後對各個自由度的運動方程式進行積分並與以疊加，以 得到整個系統問題之解。 直接積分法為對運動方程式積分前，不進行對方程式形式的變換，而直接進 行逐步數值積分，其基本概念為，在求解時間域內0< t<T，任何時間點t 均滿 足運動方程式，或者是在一相隔時間點Δt的各個離散時間點上，需滿足運動方 程 式 。 假 定 將 求 解 時 間 域 0< t<T 離 散 成 n 個 時 間 區 間 ， 則 Δt T/= n ， t t t_k+1 = _k +Δ ，若時間點t 之位移、速度與加速度均已知，且₀ t₁、t₂、……、t_k時 間點之解都已求得，則計算的目的在於求得系統於時間點t_k+1之解，並以此過程 建構整個離散求解時間域內的解。

3.2.1 中心差積分法

在直接積分法中，中心差積分法是一種顯示的積分方法，其方法為使用位移 來表示加速度與速度，假設時間間格Δt皆相同，於時間點t_k的加速度X&& 與速度_k k X& 可表示為 ) X 2X -(X 1 X _{2 -1} + ₊₁ Δ = _{k k k} k t && （3.2.14） ) X (-X 2 1 X _-1 + ₊₁ Δ = _{k k} k t & （3.2.15） 由之前的運動方程式（3.2.12）可知，時間點t_k之運動方程式滿足 k k k k CX KX Q X M&& + & + = （3.2.16） 為求時間t_k+1的位移X_k+1可經由將（3.2.14）式與（3.2.15）式帶入（3.2.16）並整 理後可得

1 -2 2 1 2 ₂ C)X 1 -M 1 ( -M)X 2 -(K -Q C)X 2 1 M 1 ( _{k t k k} t t t t t + Δ = Δ Δ Δ Δ + （3.2.17） 由於積分起步時須需要一較特別的已知值X_-1，可由下式求得 0 2 0 0 1 - X 2 X -X X = Δt& + Δt && （3.2.18） 由（3.2.17）式可知，若質量矩陣 M 與阻尼矩陣 C 為對角矩陣，或阻尼矩陣 C 為零矩陣時（無阻尼系統），求解運動方程式並不需要進行逆矩陣的轉換，且 求解下一時間點之位移X_k+1是基於時間t_k的運動方程式，此種積分方法一般稱之 為顯式（explicit）積分法，此積分法的缺點是其為有條件穩定（conditionally stable） 的算法，也就是時間步長（time step）必須小於一臨界值Δ ，解才會趨向穩定。 t_cr max cr 2 ω = Δ ≤ Δt t （3.2.19） 其中，ω_max為系統之最大自然頻率。若是為阻尼系統，且阻尼比（Damping ratio） 為ξ，則（3.2.19）修改為 ) -1 ( 2 2 max cr =_ω +ξ ξ Δ ≤ Δt t （3.2.20）

3.2.2 Newmark 積分法

1959 年 Newmark 推導一系列數值積分公式，稱之為 Newmark Method，相 較於顯式（explicit）積分法，Newmark 積分法求解時間t_k+1的位移X_k+1是基於時 間t_k+1的運動方程式，這種積分方法稱之為隱式（implicit）積分法。Newmark 積 分法是線性加速度法的延伸，其採用下列假設 t k k k k+ =X +[(1- )X + X + ]Δ

2 1 1 - )X X ] 2 1 [( X X

X_k+ = _k + &_kΔt + β &&_k +β &&_k+ Δt （3.2.22） 由上述二式可得到加速度a&& 之求解式 _k+1 k k k k k t t 2 -1)X 1 ( -X 1 -) X -(X 1

X&& _{1 2 1} & &&

β β βΔ Δ = ₊ + （3.2.23） 其中α 、β為 Newmark 積分法之參數，經由不同的參數組合可得出不同之加速 度法，例如當 6 1 = β 、 2 1 = α 時，此積分法即所謂的線性加速度法。為了求解X_k+1， 時間點t_k+1之運動方程式 1 1 1 1 CX KX Q X M&&_k+ + &_k+ + _k+ = _k+ （3.2.24） 將（3.2.21）式導入（3.2.19）式再代入（3.2.22）式後可得 ]M X 1 X 1 X 1 [ Q C)X M 1 (K _{2 k 1 k 1 2 k k k} t t t t β β β & β && α βΔ + Δ = + Δ + Δ + + _{+ +} ]C X 1) -2 ( X 1) -( X [ _{k k} t _k t + & + Δ && Δ + β α β α β α （3.2.25） Zienkiewicz 指出，Newmark 積分法若要成為非條件穩定（unconditionally stable），必須滿足下列條件 1. ) 2 1 ( 4 1 _α β ≥ + 2. 2 1 ≥ α （3.2.26） 3. 0 2 1 > + +β α

3.2.3 六面體實體元素

實體元素為三維的元素，其包含了四面體元素、三角錐元素、三角形菱柱元 素與六面體元素，其中六面體元素如【圖3.1】，線性元素為一個8 結點六面體元 素，若自然座標系統為ξ、η、ζ ，則其 j 結點的形狀函數可表示為

) )(1 )(1 (1 8 1 j j j j N = +ξξ +ηη +ζζ （3.3.1） 其中ξ_j、η_j、ζ_j於各結點的值如【表3.1】所示，則形狀函數矩陣 N 為 3×24 的 矩陣可表示為 ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ = 8 1 8 2 1 2 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 N N N N N N N N L L L （3.3.2） 應變－位移矩陣B 則為 6×24 矩陣，表示為 N 0 0 0 0 0 0 0 0 0 B × ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ = x z y z x y z y x （3.3.3） 經由偏微分的規則，形狀函數矩陣對自然座標的偏導數可以表示成 ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ ∂ ∂∂ ∂∂ ∂ × = ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ ∂ ∂∂ ∂∂ ∂ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ = ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ ∂ ∂∂ ∂∂ ∂ z Ny Nx N z Ny Nx N z y x z y x z y x N N N i i i i i i i i i J ζ ζ ζ η η η ξ ξ ξ ζ η ξ （3.3.4） 式中，J 為 Jacobian matrix，為了求得形狀函數矩陣對場座標的偏導數，可將 J 之反矩陣

⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ ∂ ∂∂ ∂∂ ∂ × = ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ ∂ ∂∂ ∂∂ ∂ ζ η ξ i i i i i i N N N z Ny Nx N 1 -J （3.3.5）

3.3 有限元素分析軟體－ANSYS

ANSYS 是一個有限元素分析軟體，其應用可於結構力學、流體力學、熱傳 導、電磁場之有限元素分析，並可結合不同分析間的互相影響，做一整體工程上 的 分 析 。LS-DYNA 原 本 是 顯 式 有 限 元 素 分 析 程 式 （ explicit finite element program），ANSYS LS-DYNA 則是整合了 ANSYS 程式裡原本有的前處理器與後 處理器功能，用來建構分析模型以及呈現分析後的資訊。在ANSYS 的內建元素 裡，提供了許多種類的元素可做分析，如梁元素、連結元素、實體元素、流體元 素、殼元素與管元素等，不同元素又可提供不同種類的分析，如力學分析、電磁 分析、熱傳導分析與流體分析等。

3.3.1 軟體基本架構與分析流程

ANSYS 基本的分析流程分為下述幾個步驟 1. 環境設定：設定環境參數，如檔案名稱，環境單位制。 2. 前處理器：在前處理器裡建立出分析結構物的實體模型，選定材料性質 參數，定義元素後，建立出有限元素模型。 3. 解題處理器：在此定義分析方式與解題方法，設立邊界條件與載重，設 定輸出入控制並求解。 4. 後處理器：整合設定收集所需之分析資料，可輸出至檔案或輸出文字圖 形資訊，並可驗證求解之結果。

若是以背景模式執行時，分析過程中可以加入些變數與流程控制，如此對於 較大複雜的分析過程提供了方便的選擇。

第四章 系統設計與案例分析

先利用有限元素動力分析程式建立並分析管結構損壞發生在不同損壞位置 與損壞程度時的加速度反應，在藉由小波分析處理所得之訊號歸納出識別參數， 最後利用此識別參數來架構整個管結構損壞識別系統。

4.1 有限元素分析方法建立數值模型

本研究有限元素結構動力分析是使用有限元素分析程式 ANSYS，模擬管結 構受一衝擊載重後自由振動之動力行為。數值模型的分析過程分為以下幾個步 驟，首先必須建立分析結構物的實體數學模型，再來決定元素與其基本性質，接 著提供外力與邊界條件，最後求解，其流程圖如【圖4.1】。

4.1.1 實體數值模型

利用ANSYS 的前置處理器可以很方便地建立實體之數學模型，本研究的實 體數值模型為一長2 公尺，管外半徑為 8 公分，管壁厚 0.8 公分，邊界條件為兩 端各取兩點做全束制之鋼管，如此在邊界的行為上提供了管壁振動的可能，其基 本的幾何與邊界條件外觀如【圖4.2】。 為了減低少計算上的成本，合理的簡化模型是可行的一種策略，理論上當管 結構物在中點受一載重時，管的縱方向與施力方向所形成的平面將會是管結構物 的在動力分析時的對稱面，利用這個對稱面可以將數學模型簡化至整體管模型的 一半，如此一來將大幅的減少分析所耗費的時間，在實際操作上在經過與未利用 對稱特性的整體數學模型做比較，假設網格大小均為0.5 公分時其節點數、元素 量則相差了約四倍，經過對稱面簡化後的實體數學模型如【圖4.3】。

對於管的損壞模式，假設是呈一環狀的損壞，以局部的降低損壞環勁度並改 變其位置以模擬之，損壞的程度即損壞環元素勁度的衰減百分比，衰減越多表示 損壞越嚴重，損壞的位置為在管的兩支承端間變動。

4.1.2 模型元素

元素挑選對有限元素分析是否正確有著關鍵性的決定，在管結構之模擬可用 殼元素（shell element）或者實體元素（solid element），兩者使用的結點數與形 狀函數都不同，由於殼元素為二維元素，所以殼元素的實體數學模型為平面模 型，而使用實體元素為三維元素，所以實體元素的實體數學模型則為三維的立體 模型，一般若是只求管結構較低的自然頻率，利用殼元素模擬之即足夠，且在相 同的元素邊長的時間計算量上亦較實體元素來的少，但對於局部破壞的識別上， 在以往的參考文獻中指出，振動訊號低頻的部份改變並不是那麼的明顯，反而是 在較高頻率的細節上有所反應，為了表現出激振力對結構物產生之高頻部份細 節，在網格大小的控制也必須要小到足以表現出健康與損壞之管結構物在高頻的 差異性，由於本研究網格化後之元素邊長將比管壁厚度小，如此一來若使用殼元 素將會有不合理的狀況，所以挑選較為合適的實體元素。實體元素中四面體元素 與六面體元素相比較之下，以同樣的元素邊長，四面體元素所需之數量比六面體 元素多，但因為分析模擬的關鍵在於激振力所造成的彈力波為沿著管長方向傳 遞，所以相較於六面體元素規則網格化而言，四面體元素對於管結構高頻細節部 份並沒有很大的幫助，但卻會造成較多的耗費，所以挑選六面體實體元素為本研 究所使用的有限元素。 元素的物理性質即管結構的材料性質，諸如楊氏係數、包松比、密度，由於 採用近似於鋼的材料，材料性質具有等向性，楊氏係數在三維空間中互相垂直的 方向皆相同，模擬管健康部分的元素楊氏係數取 2.1×1010 公斤/平方公尺，包松

比取 0.3，密度為 7800 公斤/立方公尺。對於損壞部分的元素，將楊氏係數做修 折減以模擬損壞部份勁度的改變。 網格大小的挑選，主要是依據理論上網格切的越細小則有限元素分析將會越 準確，但實際上受限於電腦的效能與精度是不可能將元素切的無限小，但還是依 據以上的觀念選擇一個較適合的網格大小，由於使用的是中間插積分法，由 （3.2.20）式可知，當元素切的越細小時，臨界時間間距也將越小以滿足計算結 果的穩定，當然同時結點數與元素數也會增加，如此雖然可以得到較精確的解（於 更高頻率的部份精確），但會對整個運算時間將大幅的增加，對於分析的主要目 的並沒有幫助太多，且訊號的取樣頻率並沒有以臨界頻率作取樣，如此在較高頻 率的振動同樣的將於取樣時被忽略掉，本研究中網格大小挑選為 0.5×0.5×0.4 公 分的六面體元素，即於管壁厚度方向為0.4 公分，其它方向為 0.5 公分，如此管 壁厚度有2 個元素，而在管長方向劃分 400 個元素，管的一半周長有 49 個元素 組成，其簡化模型經由網格化後之圖形如【圖4.4】。

4.1.3 衝擊試驗動力分析

模擬管結構所受的激振外力為一個衝擊載重，總受載時間為 0.01 秒內，在 整個分析歷程中其時間與載重關係圖為【圖4.5】，給予適當的系統阻尼後，利用 顯式結構動力法進行直接積分，整個分析時程由0 至 0.1 秒。動力分析後的結果， 加速度取樣點取7 個測點來取樣，位置為由管的兩支承端等間距分散在管上方， 其編號與位置如【表4.2、圖 4.6】。如此取點的動機是，當損壞發生在量測點上 或是發生兩者的中間時，該量測點所得之動力訊號在某些頻率區間內的能量將會 改變，藉由小波分析後判別損壞發生的位置。

4.2 訊號分析

經由動力分析後所得到的加速度歷時資料後，便可將其做數值上的分析，以 求出或比較出訊號中的差異性。使用離散小波轉換，以其多層解析的特性，將訊 號層層分解成許多不同的子訊號，將不同損害情況下各測點子訊號能量予以比 較，了解在哪一階子訊號對於損壞位置之間有著明顯的相關性，訊號分析流程圖 如【圖4.7】，本研究中訊號處理的部份是使用 MATLAB 所架構的程式語言進行 運算。

4.2.1 訊號前處理

為了減少管結構在承受衝擊時的瞬態反應干擾，以增加分析後辨識能力的準 確度。且使用離散小波轉換時，訊號取樣點數取2 的冪次方較佳，所以捨棄前面 一小段的取樣值，原本的訊號取樣點數為2500 點，將訊號取最為接近的_{2 取樣}11 數，從結束點往前擷取至第2048 個訊號點停止，如此一來可得較大部份未受載 重作用下之訊號，減低瞬態反應對訊號之干擾，取樣時間點之加速度歷時其中一 取樣訊號表示如【圖4.8】。

4.2.2 小波分析

使用小波轉換時，首先必須選擇一個母小波，對於訊號的特性選擇適當的母 小波，將造成在訊號轉換時有更佳的解析能力，但如何挑選最佳的母小波仍是值 得探討的研究，本研究中採用Meyer wavelet 做分析，由於其尺度函數φ(x)與小 波函數ψ(x)皆是定義在頻率域裡，如此經由離散小波轉換後各階層的係數將落 於某一段頻率區間內，且Meyer wavelet 也可以使用快速離散小波演算法進行快 速轉換，由於使用快速轉換後的各個子訊號可完整的重建原始訊號，如此可方便 於計算各個子訊號的能量百分比，關於Meyer wavelet 其餘性質可參考【附錄 A】。

本研究利用小波分析方法分為兩種，第一種是使用快速離散小波轉換將訊號 分解並重建至不同頻率區間的子訊號，第二種為使用離散尺度小波轉換並觀察轉 換後的小波轉換係數。利用小波對於頻率解析的特性，以偵測損壞發生時損壞位 置或附近所得之加速度訊號頻率的能量分布改變，並將此現象定量的表現出反應 變化來達到損壞辨識的能力。

4.3 分析案例與歸納結果

本研究分析案例分為以下形式： 1. 單環損壞的模式，將損壞環的寬度固定為 10 公分寬，對損壞環的勁度進行 折減，分別折減勁度為10%、20%、30%、40%與 50%，而損壞環的位置從 受力端點沿著管長軸方向移動，損壞發生的位置為（0.2）、（0.4）、（0.5）、（0.6） 與（0.8）這五種狀況總共 25 個分析案例。 2. 雙環損壞的模式裡，同樣的將每個損壞環的寬度固定為 10 公分寬，模擬與 單一環損壞情況相同的邊界條件與激振力之下所造成的加速度反應，並取同 樣位置的量測點訊號，損壞程度為兩個損壞環的損壞程度相同，模擬案例的 損壞程度為勁度衰減 20%、30%與 40%。案例模擬雙環損壞位置發生的型 式，兩個損壞環位置以其分布在管上的對稱型式分成三組討論。第一組為損 壞位置發生於同一側，於（0.2）、（0.5）位置與（0.5）、（0.8）位置【圖 4.9】， 第二組為損壞位置發生在管兩側但不對稱，於（0.25）、（0.65）位置與（0.35）、 （0.75）位置【圖 4.10】，第三組則為損壞發生管的兩側但是對稱在（0.2） 與（0.8）位置【圖 4.11】。 案例分析中的對於損壞的位置均採無因次化，即座標位置除以管長，以方便 閱讀與容易了解內容，而在量測點後將會有其無因次位置之座標，如C1（0.125）。