高鐵通過竹科之波傳研究

(1)

國立交通大學

土木工程研究所

碩士論文

高鐵通過竹科之波傳研究

Study of Wave Propagation Problem at Hsin-chu Science

Park due to Passing of High Speed Train

研究生：曹治本

指導教授：劉俊秀

(2)

高鐵通過竹科之波傳研究

研究生：曹治本指導教授：劉俊秀

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

論文摘要

本文為利用有限元素法，來分析高速鐵路通過新竹科學園區所產生的波傳振動問題。首先以有限元素法軟體建立理想化的數值模型，將分析模型數值解與理論解析解互相比較，來證明以有限元素法軟體分析波傳問題的精確性，接著討論如何選用適當的參數來模擬實際案例，對於外力作用所產生的振動量加以分析。研究結果顯示，在選用適當參數的前提下，以有限元素法分析複雜系統半無限域的波傳問題，是一種很好的方式。

(3)

Study of Wave Propagation Problem at Hsin-chu

Science Park due to Passing of High Speed Train

Thesis Advisor :Gin-Show Liou Graduate Student : Chih-Pen Tsao

ABSTRACT

The main goal of this research lies on analyzing wave propagation problem at Hsin-chu Science Park. Firstly we set up an ideal finite element modeland analyzed it to obtain same special numerical solutions. Then we compared the numerical solutions with the theoretical solution in order to examine the accuracy of the finite element method. After that, we discussed how to pick proper parameter such as element size to model the real situations. And then, wave propagation in half space medium caused by external forces are analyzed. The results show that with proper selection of parameters, analyzing the wave propagation in half space medium using the finite element method is appropriate.

(4)

誌

謝

在研究所的生活中，感謝恩師劉俊秀教授在課業上的指導與觀念上的啟發，使我能完成論文研究，在此表達我最高的敬意與謝忱；同時，在本文研究期間，承蒙交通大學黃炯憲教授與宜蘭科技大學李洋傑教授惠予寶貴意見和指正，使本論文更加完善，在此表達誠摯的謝意；另外，感謝瑜隆學長、博智學長、柏濤同學在研究學習經驗上的分享，友人國修、宗成、明憲、丞凱、在課業上的切磋砥礪與生活上的照顧，學弟紀戌、冠杰在研究期間的協助，在此一併致謝。最後，由衷感激長年來一直在精神上支持我和鼓勵我的父母親，遠你們永遠平安順遂。並以本文獻給所有關心我的朋友與家人。

(5)

目

_錄

中文摘要 --- Ⅰ 英文摘要 --- Ⅱ 誌謝 --- Ⅲ 目錄 --- Ⅳ 表目錄 --- Ⅶ 圖目錄 --- Ⅷ 第一章緒論 --- 1 1-1 研究動機--- 1 1-2 研究方法--- 2 1-3 內容大綱---3 第二章理論簡介 --- 4 2-1 彈性波傳理論--- 4 2-1-1 彈性波傳簡介---4 2-1-2 彈性波傳控制方程---5 2-1-3 實體波---9 2-1-4 表面波--- 10 2-2 有限元素法理論 --- 17

(6)

2-2-1 有限元素法分析原理 --- 17 2-2-2 靜力學的有限元素法分析--- 18 2-2-3 結構動力學的有限元素法分析--- 22 第三章半無限域波傳分析 --- 28 3-1 前言--- 28 3-2 高速鐵路高架橋振動反應分析--- 29 3-2-1 載重模型--- 29 3-2-2 結構模型--- 30 3-2-3 動態反力歷時分析結果--- 31 3-3 1/3倍八分貝頻帶圖 --- 33 3-4 地質參數判斷--- 35 3-4-1 震測地質參數--- 35 3-4-2 重錘實驗之等值地質參數--- 36 3-5 理想化分析--- 38 3-5-1 理想化分析模式--- 38 3-5-2 理想化模式分析結果--- 39 3-6 有限元素法分析--- 40 3-6-1 有限元素法分析模式--- 40 3-6-2 有限元素法分析結果--- 42

(7)

第四章結論與建議 --- 44 4-1 結論--- 44 4-2 建議--- 45 參考文獻

(8)

表

目

錄

表3-1 高架橋墩柱之斷面參數與力學性質--- 50 表3-2 高架橋大樑之斷面參數與力學性質--- 50 表3-3 1/3倍頻八分貝頻帶之上下限及中心頻率(Hz)--- 51 表3-4 土壤參數(damping ratioξ=0.01) --- 51 表3-5 重錘試驗回歸_α假設參數條件--- 52 表3-6 重錘試驗之等值土壤參數(ξ=0.05)--- 52 表3-7 高鐵列車到達各橋墩的時間--- 52 表3-8 時間延遲作用外力計算--- 52

(9)

圖

目

錄

圖2-1 壓力波波傳遞方向與介質振動方向的關係圖--- 53 圖2-2 剪力波波傳遞方向與介質振動方向的關係圖--- 53 圖2-3 雷利波波傳遞方向與介質振動方向的關係圖--- 53 圖2-4 元素表面應力分佈圖--- 54 圖2-5 半無限域表面平面波傳遞於X-Z平面 --- 54 圖2-6 各彈性波與柏松比所對應之比例關係--- 55 圖2-7 半無限域表面三種波傳攜帶能量百分比--- 55 圖3-1 THSR新幹線列車示意圖 --- 56 圖3-2 高速鐵路高架橋振動分析之結構模型--- 56 圖3-3 76K+465墩柱基礎之動態反力歷時--- 57 圖3-4 76K+495墩柱基礎之動態反力歷時--- 58 圖3-5 76K+535墩柱基礎之動態反力歷時--- 59 圖3-6 76K+595墩柱基礎之動態反力歷時--- 60 圖3-7 76K+635墩柱基礎之動態反力歷時--- 61 圖3-8 76K+665墩柱基礎之動態反力歷時--- 62 圖3-9 76K+695墩柱基礎之動態反力歷時--- 63

(10)

圖3-10 76K+725墩柱基礎之動態反力歷時--- 64 圖3-11 76K+755墩柱基礎之動態反力歷時--- 65 圖3-12 振動衰減之關係--- 66 圖3-13 3-D分析模型 --- 67 圖3-14 頻率域中X方向之外力(Rx) --- 68 圖3-15 頻率域中Z方向(垂直)之外力(Rz)--- 68 圖3-16 頻率域中Y方向之彎矩外力(My)--- 68 圖3-17 土壤Cs=816m/sec 距離300m Y方向(水平)之dB值 --- 69 圖3-18 土壤Cs=816m/sec 距離300m Z方向(垂直)之dB值 --- 69 圖3-19 土壤Cs=816m/sec 距離300m X方向(水平)之dB值 --- 69 圖3-20 兩層土壤(上層Cs=408m/sec,下層Cs=816m/sec) 距離300m Y方向(水平) 之dB值 --- 70 圖3-21 兩層土壤(上層Cs=408m/sec,下層Cs=816m/sec) 距離300m Z方向(垂直) 之dB值--- 70 圖3-22 兩層土壤(上層Cs=408m/sec,下層Cs=816m/sec) 距離300m X方向(水平) 之dB值 --- 70 圖3-23 土壤Cs=408m/sec,ξ=0.05,距離300m Y方向(水平)之dB值 --- 71

(11)

圖3-24 土壤Cs=408m/sec,ξ=0.05,距離300m Z方向(垂直)之dB值 --- --- 71 圖3-25 土壤Cs=408m/sec,ξ=0.05,距離300m X方向(水平)之dB值--- --- 71 圖3-26 有限元素法分析模型 --- 72 圖3-27 分析模型橋墩相對位置圖 --- 72 圖3-28 1-3Hz有限元素法分析模型 --- 73 圖3-29 3.5~6.0Hz有限元素法分析模型 --- 73 圖3-30 6-15Hz有限元素法分析模型 --- 74 圖3-31 二十個節點的實體元素C3D20 --- 74 圖3-32 十二個節點的無限元素CIN3D12R --- 75 圖3-33 Z方向(垂直)ABAQUS結果與解析解結果之比較(300M)---- 75 圖3-34 Y方向(水平)ABAQUS結果與解析解結果之比較(300M) --- 76 圖3-35 考慮時間延遲Pier1-Pier9外力共同作用下晶圓廠位置之三個方向之振動dB值 --- 76

(12)

第一章緒論

1-1 研究動機

隨著科技的進步，台灣的半導體產業蓬勃迅速的發展。由於半導體產品的精密度較高，製造過程中容易因為微小的振動而降低產品的良率，因此對於製程中環境振動量大小的分析與評估，成為一項相當重要的課題。高速鐵路是鐵路現代化的重要指標，目前國內正發展興建由台北至高雄的高速鐵路，其行經的路線鄰近新竹科學園區，而由高速鐵路行車所引致之振動必然會對精密產業之良率造成影響，所以有必要對這股巨大能量所產生的波傳振動問題加以研究。波傳所產生的振動問題，常用現地實驗、理論數學求解、數值模擬計算等不同的方式加以分析，就各種分析方法的比較方面；現地實驗的結果較為符合實際情況，但所費不貲且結果不適用於其他區域。理論數學求解具有理論基礎，可是必須藉由許多的假設來簡化求解的難度，如此過多的簡化步驟常使推算的波傳行為偏離真實現象。數值模擬計算經濟簡便且適用性良好，雖然其中也存在著數值解的精確性及處理速度受到計算能力限制的問題，但是隨著近年來計算機性能的

(13)

大幅提昇，已逐漸克服了此項缺點，因此也使得數值模擬成為較常被採用的分析方法。

目前應用最普遍的數值模擬方法為有限差分法(finite difference method) 及有限元素法(finite element method)。有限差分法為較早被採用的數值方法，主要是將控制方程式依差分原理離散成一組聯立方程式，再進行求解，可是對於具不規則邊界條件的問題，有限差分法就會變的較複雜且無法系統化。有限元素法的分析則是將整個問題直接以有限元素建立數值模型模擬加以求解，由於有限元素法具有高度的操作彈性，特別適用於複雜的幾何形狀與土壤性質等實際情況，因此本文選擇以有限元素法來分析半無限域之波傳問題，並探討如何選用適當的有限元素法參數，才能使分析結果能更加的完整及正確。

1-2 研究方法

本文主題在於討論高鐵通過竹科之波傳問題，並以有限元素法 (finite element method)程式來建立實際案例模型加以分析。

首先利用有限元素軟體ABAQUS，建立理想化的三維數值模型，並將分析出來之數值解，與相同條件以土壤阻抗矩陣計算出之理論解析解比較，來證實以有限元素法計算波傳問題的精確性及可行性。

(14)

接著以「高速鐵路行經新竹科學園區引致電子廠廠區的振動」為數值模擬分析案例；同樣利用有限元素軟體 ABAQUS，先依廠區附近的高程資料建立幾何模型，由地質鑽探報告判斷土壤的材料參數，由以上資料建立與實際情形相仿之數值模型後，再將高鐵列車行經時頻率域外力之大小輸入，經由有限元素法軟體的運算而對半無限域之波傳行為產生之振動量加以分析。

1-3 內容大綱

本論文的內容共分為四章。第一章為緒論，說明研究動機及研究方法。第二章為理論簡介，簡述半無限域受力彈性波數學理論及有限元素法基本理論。第三章為半無限域波傳分析，首先介紹動態反力分析方式，分貝轉換理論，以及地質參數的判斷方法。接著介紹理想化分析模式，並將理論計算所得之解析解與有限元素法數值模型計算所得之數值解相互比較，驗證有限元素法分析之精確性，最後經由以上的結論與現地高程資料來建立實際案例的數值模型加以分析。第四章為結論與建議，將本研究之分析結果，歸納出結論，並對未來相關研究提出建議。

(15)

第二章理論簡介

2-1 彈性波傳理論

2-1-1

彈性波傳簡介

波傳可視為一種能量的傳遞，而能量向外傳遞必需仰賴介質，當傳遞的介質為彈性時，此時即稱為彈性波，彈性波又因其傳遞路徑的不同，區分為實體波（body wave）和表面波（surface wave）兩種形式。

實體波在物體內部和表面傳播，又依波傳遞方向與介質振動方向的關係而分為壓力波（Primary wave）及剪力波（Secondary wave）；壓力波通過物體時傳播介質的振動方向和波傳方向一致或平行，如圖 2-1。剪力波通過物體時，傳播介質的振動方向和波傳方向垂直，如圖2-2 所示。表面波運動範圍主要在自由表面附近，又稱為雷利波 (Rayleigh wave) ，其傳播介質的振動方式是在垂直面上的橢圓運動，如圖 2-3。由於表面波的傳遞速度以及傳遞時能量隨距離的衰減均較實體波慢，因此常被視為造成地表振動的主要來源之一，於相關的研究中相當受到重視。

(16)

2-1-2 彈性波傳控制方程

彈性波傳遞控制方程式的推導已有百年以上的歷史，Richard【1】等曾於著作中，利用材料的組成律、應變與位移的幾何關係，加上力平衡原理，來推得彈性波傳遞控制方程式，以下就其推導過程加以說明

(1) 材料的組成律：首先由廣義虎克定律（General Hooke’ s law）得知，介質元素受力產生變形，應力與應變的關係可由 [σ]=[D][ε]來決定，而當介質為均質、等向、線彈性的材料時，材料的勁度矩陣[D] 可用兩個彈性常數（Lame’s constants）λ、μ 表示。 ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ = ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ xz yz xy zz yy xx xz yz xy zz yy xx D D D D D D D D D D D D D D D D D D D D D D D D D D D D D D D D D D D D γ γ γ ε ε ε σ σ σ σ σ σ 66 65 64 63 62 61 56 55 54 53 52 51 46 45 44 43 42 41 36 35 34 33 32 31 26 25 24 23 22 21 16 15 14 13 12 11 (2-1)

[ ]

⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ + + + = μ μ μ μ λ λ λ λ μ λ λ λ λ μ λ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 0 0 0 2 0 0 0 2 D (2-2) 其中

(17)

G ) 2(1 E ) 2 1 )( 1 ( = + = − + = υ μ υ υ υ λ E λ、μ= Lame’s constants E= Young’s Modulus G= Shear modulus υ= Poisson ratio σ=作用於介質元素上的應力 ε=介質元素的軸向應變 γ=介質元素的剪應變 (2)應變與位移的幾何方程關係式：在彈性及均質等向的材料中，位移與應變的幾何方程關係式定義為 ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ ∂ ∂ + ∂ ∂ ∂ ∂ + ∂ ∂ ∂ ∂ + ∂ ∂ ∂ ∂∂ ∂ ∂ ∂ = ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ z u x w y w z v x v y u z wy v x u γ γ γ ε ε ε xz yz xy zz yy xx (2-3) 其中u、v、w為介質元素 x、y、z 三方向的變位將式(2-2)、式(2-3)代入式(2-1)，得

(18)

⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ ∂ ∂ + ∂ ∂ ∂ ∂ + ∂ ∂ ∂ ∂ + ∂ ∂ ∂ ∂∂ ∂ ∂ ∂ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ + + + = ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ z u x w y w z v x v y u z w y v x u μ μ μ λ λ λ λ μ λ λ λ λ μ λ σ σ σ σ σ σ xz yz xy zz yy xx 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 0 0 0 2 0 0 0 2 μ (2-4) 式(2-4)展開後可得波傳介質元素應力與應變的關係式 ) )( ( ) )( ( ) )( 2 ( z w y v x u xx _∂ ∂ + ∂ ∂ + ∂ ∂ + = λ μ λ λ σ (2-4a) ) )( ( ) )( 2 ( ) )( ( z w y v x u yy _∂ ∂ + ∂ ∂ + + ∂ ∂ = λ λ μ λ σ (2-4b) ) )( 2 ( ) )( ( ) )( ( z w y v x u zz _∂ ∂ + + ∂ ∂ + ∂ ∂ = λ λ λ μ σ (2-4c) ) )( ( x v y u xy _∂ ∂ + ∂ ∂ = μ σ (2-4d) ) )( ( y w z v yz _∂ ∂ + ∂ ∂ = μ σ (2-4e) ) )( ( x w z u xz _∂ ∂ + ∂ ∂ = μ σ (2-4f) (3)力平衡方程式：根據牛頓第二運動定律，一個邊長為 dx、dy、 dz，密度為 ρ 之微小元素，若忽略徹體力(body force)之作用，如圖 2-4，可求得直角座標系三個方向的力平衡控制方程式

(19)

2 2 ) ( ) ( ) ( 0 t u dxdydz dxdy dxdy dz z dxdz dxdz dy y dydz dydz dx x F zx zx zx yx yx yx xx xx xx X ∂ ∂ = − ∂ ∂ + + − ∂ ∂ + + − ∂ ∂ + =

∑

ρ σ σ σ σ σ σ σ σ σ 整理可得 z y x t u xx yx _zx ∂ ∂ + ∂ ∂ + ∂ ∂ = ∂ ∂ σ σ σ ρ ₂ 2 (2-5) 同理，由

∑

Fy=0可得 z y x t v xy yy zy ∂ ∂ + ∂ ∂ + ∂ ∂ = ∂ ∂ σ σ σ ρ ₂ 2 (2-6) 由

∑

Fz=0可得 z y x t w xz yz _zz ∂ ∂ + ∂ ∂ + ∂ ∂ = ∂ ∂ σ σ σ ρ ₂ 2 (2-7) 將波傳介質應力與應變的關係式 (2-4 a-f)，分別代入力平衡控制方程式(2-5)、(2-6)、(2-7) ，可得 z x w y x v z u y u x u t u ∂ ∂ ∂ + + ∂ ∂ ∂ + + ∂ ∂ + ∂ ∂ + ∂ ∂ + = ∂ ∂ 2 2 2 2 2 2 2 2 2 2 ) ( ) ( ) ( ) ( ) 2 (λ μ μ μ λ μ λ μ ρ (2-8a) z y w y x u z v y v x v t v ∂ ∂ ∂ + + ∂ ∂ ∂ + + ∂ ∂ + ∂ ∂ + + ∂ ∂ = ∂ ∂ 2 2 2 2 2 2 2 2 2 2 ) ( ) ( ) ( ) 2 ( ) (μ λ μ μ λ μ λ μ ρ (2-8b) z y v z x u z w y w x w t w ∂ ∂ ∂ + + ∂ ∂ ∂ + + ∂ ∂ + + ∂ ∂ + ∂ ∂ = ∂ ∂ 2 2 2 2 2 2 2 2 2 2 ) ( ) ( ) 2 ( ) ( ) (μ μ λ μ λ μ λ μ ρ (2-8c) 其中若定義體積應變 z w y v x u zz yy xx _∂ ∂ + ∂ ∂ + ∂ ∂ = + + =ε ε ε ε 、拉普拉斯運算子

(20)

2 2 2 2 2 2 2 y z x ∂ ∂ + ∂ ∂ + ∂ ∂ = ∇ ，式(2-8a、b、c)可改寫為 u G x G t u 2 2 2 ) ( + ∇ ∂ ∂ + = ∂ ∂ ε λ ρ (2-9a) v G y G t v 2 2 2 ) ( + ∇ ∂ ∂ + = ∂ ∂ ε λ ρ (2-9b) w G z G t w 2 2 2 ) ( + ∇ ∂ ∂ + = ∂ ∂ ε λ ρ (2-9c) 式(2-9a)、(2-9b)及式(2-9c)即為常見之三維彈性波傳運動控制方程式。

2-1-3 實體波

三維彈性波傳運動控制方程式於無限域中，若考慮不同的運動形式，可求得兩個解。一個解描述純體積變化所產生的傳遞，即為壓力波，而另一個解描述純旋轉波的傳遞，為剪力波，兩者又因其主要於介質內部傳遞，而統稱為實體波。壓力波傳遞的方式為介質粒子疏密相間造成擠壓，藉由軸向彈力恢復的方式傳遞能量，若將式(2-9a)、(2-9b)及式(2-9c)分別對 x、y、z 微分後相加可得 ε λ ε ρ 2 2 2 ) 2 ( + ∇ = ∂ ∂ G t 或 ε 2 2ε 2 2 ∇ = ∂ ∂ p v t (2-10)

(21)

其中 ρ λ G v_p = +2 =壓力波波速 (2-11) 式(2-10)可視為壓力波波傳介質的膨脹量ε以速度ν_p來傳遞。剪力波的能量藉由介質粒子間切向剪應力來進行傳遞，若以式 (2-9b) 對 z 微分及式(2-9c) 對 y 微分，兩式相減後可得 x x G t ω ω ρ 2 2 2 ∇ = ∂ ∂ 或 x _v_s _x t ω ω 2 ₂ 2 2 ∇ = ∂ ∂ (2-12) 其中 ρ G v_s = =剪力波波速 (2-13) x ω 亦可用ω_y、ω_z 來表示， ω_x、ω_y 、ω_z 為繞 x 軸、y 軸、z 軸之旋轉角，與變位的關係為 z v y w x _∂ ∂ − ∂ ∂ = ω 2 x w z u y _∂ ∂ − ∂ ∂ = ω 2 y u x v z _∂ ∂ − ∂ ∂ = ω 2 式(2-12)顯示剪力波波傳介質的旋轉量ω 是以速度ν_s來傳遞。

2-1-4 表面波

波傳運動控制方程式於半無限域中，因為受到自由表面邊界條件的改變，除了壓力波及剪力波兩個解外，尚可求得第三種波傳形式的解，即為雷利波。雷利波為Lord Rayleigh【2】首先提出，其後 Lamb 【3】曾詳細描述其性質，又因其傳遞範圍受限於半無限域自由表面

(22)

附近，而稱為表面波。雷利波在求解時，考慮一作用於彈性體表面的平面波動，如圖 2-5。於此系統中，平面波朝 x 方向傳遞，座標 x、z 方向的位移為 u、 w，座標 y 方向的位移為 0，若定義壓力場勢能函數Φ、剪力場勢能函數Ψ來描述 x、z 方向的位移，其中 z x u ∂ ∂ + ∂ ∂ = Φ Ψ x z w ∂ ∂ − ∂ ∂ = Φ Ψ 則相關式可用勢能函數Φ、Ψ表示成體積應變 Φ ) Ψ Φ ( ) Ψ Φ ( ₌_∇2 ∂ ∂ − ∂ ∂ ∂ ∂ + ∂ ∂ + ∂ ∂ ∂ ∂ = ∂ ∂ + ∂ ∂ = + = x z z z x x z w x u zz xx ε ε ε 旋轉角 Ψ ) Ψ Φ ( ) Ψ Φ ( 2 ₌_∇2 ∂ ∂ − ∂ ∂ ∂ ∂ − ∂ ∂ + ∂ ∂ ∂ ∂ = ∂ ∂ − ∂ ∂ = x z x z x z x w z u y ω 二維運動控制方程式則可寫成 ) ( ) Φ ( ) 2 ( ) Ψ ( ) Φ ( 2 2 2 2 2 2 ψ λ ρ ρ ∇ ∂ ∂ + ∇ ∂ ∂ + = ∂ ∂ ∂ ∂ + ∂ ∂ ∂ ∂ z G x G t z t x (2-14a) ) ( ) Φ ( ) 2 ( ) Ψ ( ) Φ ( 2 2 2 2 2 2 ψ λ ρ ρ ∇ ∂ ∂ − ∇ ∂ ∂ + = ∂ ∂ ∂ ∂ − ∂ ∂ ∂ ∂ x G z G t x t z (2-14b) 二維彈性波傳運動控制方程式亦可求得壓力波與剪力波兩個解。控制方程式解出的第一個解為壓力波。將式(2-14a) 對 x 做微分，式

(23)

(2-14b)對 z 做微分，再進行相加可得 Φ ∇ + = ∂ Φ ∂ 2 2 2 ) 2 ( G t λ ρ 或 = ∇ Φ ∂ Φ ∂ 2 2 2 2 p t ν (2-15) 其中 ρ λ G v_p = +2 =壓力波波速 (2-16) 控制方程式解出的第二個解為剪力波。將式(2-14a) 對 z 微分及式(2-14b) 對 x 微分，兩式相減後可得 Ψ ∇ = ∂ Ψ ∂ 2 2 2 ) (G t ρ 或 = ∇ Ψ ∂ Ψ ∂ 2 2 2 2 s t ν (2-17) 其中 ρ G v_s = =剪力波波速 (2-18) 若再考慮半無限域自由表面邊界條件時，控制方程式可以解出第三個解，即為雷利波。假設雷利波為簡諧波，波移動的方向為正x 方向， F(z)、G(z)為描述雷利波擴張、旋轉隨深度變化之振幅函數，N 為波數，L 為波長，ω為頻率，勢能函數Φ、Ψ可表示成 )] ( F(z)exp[ Φ= i ωt−Nx (2-19a) )] ( G(z)exp[ Ψ = i ωt−Nx (2-19b) 其中 L N = 2π

(24)

將式(2-19a、b) 代入式(2-15)、(2-17)控制方程，可得 (z) F F(z) F(z) 2 2 2 ′′ + − = − N p υ ω (z) G G(z) G(z) 2 2 2 ′′ + − = − N s υ ω 經由整理可得 0 F(z) ) ( (z) F_′′ ₋ 2 ₋ 2₂ ₌ P N υ ω _(2-20a) 0 G(z) ) ( (z) G ₂ 2 2 ₋ ₌ − ′′ S N υ ω _(2-20b) 若令 ( ₂) 2 2 2 P N q υ ω − = 、 ( ₂) 2 2 2 S N s υ ω − = ，則式(2-20a)、(2-20b)可寫為 F(z)、 G(z)之二階常微分方程式 0 F(z) (z) F_′′ ₋_q2 ₌ _(2-21a) 0 G(z) (z) G_′′ ₋_s2 ₌ _(2-21b) 二階常微分方程式解的形式為 ) exp( B ) exp( A F(z)= ₁ −qz + ₁ qz (2-22a) ) exp( B ) exp( A G(z)= ₂ −sz + ₂ sz (2-22b) 式中等號右側第二項顯示，雷利波振幅函數隨著深度Z 的增加會變成無限大，與真實現象不符，因此B₁、B₂需為 0，最後振幅函數的解為 ) exp( A F(z)= ₁ −qz (2-23a)

(25)

) exp( A G(z)= ₂ −sz (2-23b) 將式(2-23a)、(2-23b)代回式(2-19a)、(2-19b)，於是勢能函數Φ、Ψ最後可表示成 )] ( exp[ A Φ= ₁ −qz+i ωt−Nx (2-24a) )] ( exp[ A Ψ = ₂ −sz+i ωt−Nx (2-24b) 接著由自由表面的邊界條件可知，在半無限域的自由表面Z=0 處應力 0 ) )( G 2 ( ) )( ( = ∂ ∂ + + ∂ ∂ = z w x u zz λ λ σ (2-25a) 0 ) )( G ( = ∂ ∂ + ∂ ∂ = x w z u xz σ (2-25b) 由位移 z x u ∂ ∂ + ∂ ∂ = Φ Ψ 、 x z w ∂ ∂ − ∂ ∂ = Φ Ψ的定義及勢能函數式(2-24a)、(2-24b) 之表示式，可將邊界條件式(2-25a)、(2-25b) 改寫為 0 G A 2 ] G) 2 [( A 2 2 ₂ 1 0 = + − − = = q N i Ns z zz λ λ σ 0 ) ( A A 2 2 2 2 1 0 = + + = = i Nq s N z xz σ 整理後可得 0 1 G 2 2G) ( A A 2 2 2 1 + − ₋ ₌ Ns i N q λ λ _(2-26a) 0 1 ) ( 2 A A 2 2 2 1 ₊ ₌ + N s qiN _(2-26b) 以上兩式相加可得 2 2 2 2 ₂ G 2 2G) ( N s qiN Ns i N q + − = − + λ λ

(26)

整理得 ] G) 2 )[( ( 4_qGsN2 ₌ _s2 ₊_N2 _λ₊ _q2 ₋_λ_N2 _(2-27) 將_q2_、_s2_{代入式(2-27)} 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 4 2 ₍ ₎₍ ₎ _[( ₂_G)( ₎ _] _[ ₍ _)] G 16 S P S P N N N N N N N ν ω λ ν ω λ ν ω ν ω ₋ ₌ ₊ ₋ ₋ ₊ ₋ − 並將等號兩邊同除_{G N}2 8_得 2 2 2 2 2 2 2 2 2 2 2 2 2 2 ) 2 ( ] ) )( G G 2 ( 2 [ ) 1 )( 1 ( 16 N N N N _S _P _S P ν ω ν ω λ ν ω ν ω ₋ ₌ ₋ + ₋ − (2-28) 由前述雷利波頻率為ω、波數為N 、波長為L、 L N = 2π 、 R N ν ω = =雷利波波速等之定義，式(2-28)可改寫成 2 2 2 2 2 2 2 2 2 2 ) 2 ( ] ) )( G G 2 ( 2 [ ) 1 )( 1 ( 16 S R P R S R P R ν ν ν ν λ ν ν ν ν ₋ ₌ ₋ + ₋ − (2-29) 若令 2 2 2 2 K P R _α ν ν = ， 2 2 2 K S R ₌ ν ν ， ₂ ₂ 2 1 G G 2 G G 2 α λ ρ ρ λ ν ν = + = + = S P _{則式(2-29)可寫為} 2 2 2 2 2 2 2 2 2 ₎₍₁ ₎ ₍₂ 1 ₎ ₍₂ ₎ 1 ( 16 − K −K = − α K −K α α (2-30) 展開可得 0 ) 1 ( 16 ) 16 24 ( 8 4 2 2 2 6 ₋ _K ₊ ₋ _{α K} ₊ _α ₋ ₌ K (2-31) 將 ) 2 1 )( 1 ( υ υ υ λ − + = E 、 ) 2(1 E G υ + = 代入 G G 2 1 2 + = λ α ，得 υ υ α 2 2 2 1 2 − − = 若介質柏松比υ已知，可求得_α2_{後代回式(2-31)，最後可解得}_K_之正實

(27)

數解，即可求出雷利波波速ν_R =K ν_S。比較上述提到的三種波波速，壓力波波速最快，其次為剪力波，雷利波則最慢；一般而言，三種波的波速為有條件之比例關係，此關係視土壤之柏松比（Poisson ratio）而定，如圖 2-6 所示。就傳遞能量而言，Miller 和 Pursey【4】對一均勻分佈之圓形基礎，在一均質、等向、彈性半無限空間地表上，以垂直方向震動所得之能量分配結果顯示，三種彈性波的能量分配為雷利波 67%、剪力波 26%、壓力波 7%，如圖 2-7。由於雷利波的能量佔總能量的67%，而且其傳遞主要集中於介質表面，是以在考慮地表振動時，雷利波將是一項很重要的影響因素。

(28)

2-2 有限元素法理論

2-2-1

有限元素法分析原理

有限元素法是一種新興分析結構問題的數值計算方法，這個名稱是由美國加州大學教授克拉夫Clough,R.W.【5】所首次引用。有限元素法原理是將原本複雜、龐大的彈性連續體結構，分割成各個簡單的計算元素或單元(Element)，這個人為地劃分元素的過程稱為彈性連續體的離散化，若是把離散化後各個元素組合成計算模型，並假設元素與元素之間僅在有限的節點(Node Point)上連接，來近似的替代原來的彈性連續體，如此就可以將原來分析彈性連續體上無限個質點變形和應力的問題，轉化成分析離散化計算模型中有限個節點變形和應力的問題。理論上來說，當元素劃分的越小、元素和節點的數目越多，則離散化模型越接近真實彈性連續體結構，計算分析的結果也與實際越接近，對於極限狀況，當將元素劃分無限小、元素和節點的數目無限多時，則離散化模型就變成與彈性連續體完全相同，計算分析的結果就與實際完全一樣了。但對於真實問題而言，元素和節點的數目總是有限的，因而此種方法稱之為有限元素法。

(29)

2-2-2

靜力學的有限元素法分析

以有限元素法求解結構問題，相關的理論及著作都已被編寫成參考書籍，Richard H【6】、王新榮【7】等；有限元素法分析問題常見的步驟為先將結構體離散化，建立元素的特性矩陣，接著將元素特性矩陣加以組集成整體結構矩陣進而求解。以下先就靜力學有限元素法分析的原理加以介紹。一、元素勁度矩陣的推導 1. 設定元素位移函數若假設元素的局部座標系為oxyz，可以依據元素的性質，將元素內任意點x,y,z三方向的位移分量u,v,w，寫成以多項式形式表示的位移函數u(x,y,z)、v(x,y,z)、w(x,y,z)，或寫成矩陣方程的形式為

{ }

[ ( , , )]{ } ) , , ( ) , , ( ) , , ( α z y x S z y x w z y x v z y x u d = ⎪ ⎭ ⎪ ⎬ ⎫ ⎪ ⎩ ⎪ ⎨ ⎧ = (2-32) } {d 是元素體內任意點u,v,w的位移分量 } {α 為待定係數 2. 以節點位移來表示元素內任意點的位移函數將元素上各個節點的座標值代入式(2-32)中，寫成矩陣的形式則得

(30)

{ }

α ] [ } {q = C 求解出待定係數得

{ }

q C_] 1 [ } {α = − (2-33)

{ }

T 2 2 2 1 1 1 ] [u v w u v w " q = 為元素的節點位移將式(2-33)代回式(2-32)得

{ }

q N

{ }

q C S d} [ ][ ] [ ] { ₌ −1 ₌ _(2-34) 其中[N]稱為形狀函數矩陣式(2-34)即是用節點位移

{ }

q 來表示元素體內任意點位移函數{d}的關係式。 3. 根據幾何方程式由位移函數求應變{ε}和應變矩陣[B] 將式(2-34)以節點位移表示之u、v、w，代入幾何方程關係式(2-3) 中，計算元素的應變

{ }

ε ，把結果寫成矩陣形式，則有

{ }

ε =[B(x,y,z)]

{ }

q (2-35) 其中[B]稱為元素的應變矩陣 4. 依據物理方程，透過應變{ε}求應力{σ} 將式(2-35)代入應力應變的關係式(2-1)中，計算元素的應力

{ }

σ 得

{ }

σ =[D]{ε}=[D][B]

{ }

q (2-36) 5. 以虛功方程式求出元素勁度矩陣就靜力學而言，建立元素勁度矩陣常用虛功方程代替平衡微分方

(31)

程，而三維問題的虛功方程為

{ }

q

{ }

dxdydz V ε σ T * T * _F ₌

∫∫∫

_(2-37) 式中

{ }

F =作用於元素各節點的外力

{ }

_q* _{=元素各節點處的虛位移}

{ }

σ =外力作用下元素任意點上的應力

{ }

_ε* _{=元素任意點上與虛位移相對應的虛應變} dx、dy、dz為微小元素體的邊長，它們的乘積表示微小元素體的體積dV，由於

{ }

_ε* _與

{ }

_q* _{相協調，它們之間也可以寫出與式(2-35)類似的} 關係式

{ }

_ε* ₌_[_B_]

{ }

_q*

{ } { }

_ε* T ₌ _q* T_[B_]T_(2-38) 將式 (2-38)代入式(2-37)，虛功方程可改寫成

{ }

=

_∫∫∫

{ }

=

{ }

_∫∫∫

{ }

V V q B dV q B dV q* T _F * T_[ _]T σ * T _[ _]T σ 由於虛位移為任意，所以可以由等號兩端消去

{ }

_q* T_得

{ }

=

_∫∫∫

{ }

V B σ dV T ] [ F (2-39) 最後將式(2-36)代入式(2-39)中，可得元素勁度矩陣的一般公式

{ }

B D B dV

{ }

q K

{ }

q V [ ] [ ][ ] [ ] F ₌ T _⋅ ₌

∫∫∫

(2-40)

∫∫∫

= B D B dV K] [ ] [ ][ ] [ T

(32)

[ ]

K 稱為元素勁度矩陣二、整體結構分析 1. 建立整體結構的靜力平衡方程式將局部座標系中各個元素的勁度矩陣

[ ]

K ，經由座標轉換並加以組集，可得結構的總勁度矩陣K，把各個節點外力對應組集，可得結構的總外力行矩陣F，最後可以求出整體結構的平衡方程式 Kq F = (2-41) 2. 進行邊界條件處理將結構節點位移的支承條件引入式(2-41)，使之成為正定的線性方程組。 3. 解方程組，求節點位移針對結構的總勁度矩陣K的特點，選用合適的線性方程組求解方法，對邊界條件處理後的方程組F =Kq求解，可求得結構各節點的位移q。

4.

根據節點位移求應力

若要計算應力，則在計算出節點位移q後，從中挑選出欲計算元素對應的節點位移，組成

{ }

q ，應用

{ }

ε =[D]

{ }

q 和

{ }

σ =[D]

{ }

ε =[D][B]

{ }

q ，即可求出相應節點的應力。

(33)

2-2-3

結構動力學的有限元素法分析

對於複雜結構的動力學計算分析，有限元素法也是一種很有效的工具。動力學問題的有限元素分析方法也是要把彈性連續體結構離散為有限個元素群的組合體。首先進行每個元素的特性分析，其中包含了元素勁度矩陣、元素質量矩陣、元素阻尼矩陣的計算，接著再把各個元素的特性矩陣組集起來，組成整個結構的總勁度矩陣、總質量矩陣、總阻尼矩陣，從而形成結構的動力學方程式之後再進行求解。以下就動力學的有限元素法分析的原理加以介紹。一、元素特性矩陣的推導由虛功原理的物理意義可知，外力在虛位移上所做的虛功等於彈性體內各點的應力在虛應變上所做的虛功的總和。在動力學問題中，元素外力項除了施加於節點的動態外力，還包括了元素的慣性力及阻尼力。若元素節點上的激振力為

{ }

F ，則

{ }

F 所做的虛功為

{ }

_q* T

{ }

F _(2-42a) 元素的慣性力與加速度成正比，方向相反。設元素材料的密度為ρ，則單位體積慣性力可寫成

{ }

f _ρ =−ρ

{ }

d } {d 為元素體內任意點加速度函數

(34)

將上式乘以元素任意點上的虛位移

{ }

_d* _{，並對元素體積域積分，可得} 慣性力所做的虛功為

{ }

d

{ }

d dV V T * _ρ

∫∫∫

− (2-42b) 假設阻尼力與速度成正比，方向相反，即假設為線性黏性阻尼，其阻尼係數為υ，則單位體積阻尼力為

{ }

f _υ =−υ

{ }

d } {d 為元素體內任意點速度函數同理，阻尼力所做的虛功為

{ }

d

{ }

d dV V T * _υ

∫∫∫

− (2-42c) 由式(2-42a)、(2-42b)、(2-42c)，可將動力學問題的虛功方程寫為

{ }

q

{ }

d

{ }

d dV

{ }

d

{ }

d dV

{ }

dV V V V ρ υ ε σ T * T * T * T * _F

∫∫∫

− = − (2-43) 因為有限元素法中施加的外力必須先轉化成作用於節點上的等值力，因此須經由關係式_{_d*_}₌_[_N_]

{ }

_q* _、_{_d_}₌_[_N_]

{ }

_q _、_{_d_}₌_[_N_]

{ }

_q _等，將式(2-43)改寫成節點位移表示式

{ }

q B D B

{ }

q dV

{ }

q N N

{ }

q dV V V [ ] [ ][ ] [ ] [ ] T T * T T *

∫∫∫

_ρ

∫∫∫

+

{ }

_q* T [ _N]T [ _N]

{ }

_q _dV

{ }

_q* T

{ }

F V = +

∫∫∫

υ (2-44) 其中虛位移

{ }

_q* _和

{ }

_q _、

{ }

_q _、

{ }

_q _{均可看成常量，可以提到積分號外，} 又

{ }

_q* _{為任意值上均成立，可以把它從等式兩端消去，這樣式(2-44)}

(35)

可寫為

{ }

q N N dV

{ }

q dV B D B V V ⋅ +

∫∫∫

⋅

∫∫∫

[ ]T[ ][ ] [ ]T_ρ[ ]

{ } { }

F ] [ ] [ T _⋅ ₌ +

_∫∫∫

N N dV q V υ (2-45) 式(2-45)的物理意義為動力學問題中元素節點上的彈性力、慣性力、和阻尼力與外加激振力相平衡。其中若令元素質量矩陣

[ ]

M N N dV V[ ] [ ] T

∫∫∫

= ρ (2-46) 元素阻尼矩陣

[ ]

=

_∫∫∫

V N N dV C [ ]T_υ[ ] _(2-47) 元素勁度矩陣

∫∫∫

= V B D B dV K] [ ] [ ][ ] [ T _(2-48) 則式(2-45)可改寫成

[ ]

M

{ }

q +

[ ]

C

{ }

q +

[ ]

K

{ } { }

q = F (2-49) 式(2-49)稱為元素的動力學方程式。二、整體結構分析 1. 元素矩陣的座標轉換在建立元素的動力學方程式之後，可經由座標轉換來得到整體彈性結構方程式。已知元素的座標轉換矩陣

[ ]

T ，可寫出元素節點虛位移

(36)

{ }

_{q 、}* _節點位移

{ }

_q _{、節點加速度}

{ }

_q _{在統一座標系及局部座標系的轉換} 關係式

{ }

q =

[ ]

T

{ }

q (2-50)

{ }

q =

[ ]

T

{ }

q (2-51)

{ }

_q* ₌

[ ]

_T _q

{ }

* _(2-52) 同理，元素節點力的座標變換關係為

{ }

F =

[ ]

T

{ }

F (2-53)

{ }

F =

[ ]

T T

{ }

F (2-54) 上式中上方加橫線的量是相對於統一座標系的，而不加橫線的量是相對於局部座標系。首先介紹元素勁度矩陣的座標轉換。將式(2-40)代入式(2-54)可得

{ }

F =

[ ] [ ]

T T K

{ }

q 再將式(2-50) 代入得

{ }

F =

[ ] [ ][ ]

T T K T

{ }

q (2-55) 式(2-55)左端

{ }

F 為統一座標系下元素的節點力，右端

{ }

q 為統一座標系下元素的節點位移，因此統一座標系下元素的勁度矩陣

[ ]

K 可知為

[ ]

K =

[ ] [ ][ ]

T T K T (2-56) 也就是說，只要先求得局部座標系下的元素勁度矩陣

[ ]

K ，以及元素的座標轉換矩陣

[ ]

T ，之後再進行

[ ] [ ][ ]

T T K T _{的乘積計算，就可求得統一}

[ ]

(37)

接著介紹元素質量矩陣的座標轉換。於兩種座標系中，元素慣性力所做的虛功應相等，即

{ } [ ]

q* T

(

₋ M

{ }

q

)

₌

{ }

q* T

(

₋

[ ]

M

{ }

q

)

_(2-57) 將式(2-51)、 (2-52)代入等號右項，可得

{ } [ ]

q* T

(

₋ M

{ }

q

)

₌

{ }

q* T

[ ] [ ][ ]

T T

(

₋ M T

{ }

q

)

_(2-58) 因為虛位移

{ }

_q* T_{為任意，所以可以由等號兩端消去得}

[ ]

M

{ }

q =−

[ ] [ ][ ]

T T M T

{ }

q − 比較上式兩端，就可以得出統一座標系下的質量矩陣

[ ]

M =

[ ] [ ][ ]

T T M T (2-59) 若將式(2-46)代入得

[ ]

M

[ ]

T N N

[ ]

T dV N N dV V V [ ] [ ] [ ] [ ] T T T

∫∫∫

= = ρ ρ (2-60) 式中

[ ]

N 稱為統一座標系下的形狀函數矩陣，若已知

[ ]

N ，則統一座標系下的質量矩陣

[ ]

M 可由式(2-60)求得。最後由相同原理，可得兩種座標系中元素阻尼矩陣的變換公式

[ ]

C =

[ ] [ ][ ]

T T C T (2-61) 2. 整體結構動力學方程式的建立與靜力學問題相同，於動力學中建立元素特性矩陣之後，也要把元素群綜合成為結構的計算模型，再建立起該計算模型的動力學方程式，並以求解計算模型的動力學方程式來近似替代求解原來的彈性連

(38)

續體結構的動力學方程。先將各元素座標變換後的

[ ]

K 、

[ ]

M 、

[ ]

C 進行組集，則得總勁度矩陣 K、總質量矩陣 M 和總阻尼矩陣 C，再將各元素的

{ }

F 進行組集，則得總外力行矩陣 F，因此對於整個結構而言，就有 F Kq q C q M+ + = (2-62) 式(2-62)即為整個彈性結構的動力學方程式，也就是用有限元素法來求解動力學問題的基本方程式。與靜力學分析相同，式(2-62)尚未考慮結構的支承條件，因此也要先進行支承條件處理後再加以求解。

(39)

第三章半無限域波傳分析

3-1 前言

為滿足台灣地區快速且大量運輸的需求，興建高速鐵路的計畫因應而生。高速鐵路為求高速行駛下的安全與舒適，全程大部份都是以高架橋方式設計，由於鐵路列車相當重，而且又是高速行駛，對於沿線附近必然會產生額外的振動影響。高速鐵路雖然沒有直接從新竹科學園區穿過，但是離最近的電子公司廠房只有 285 公尺，而晶圓廠之生產線為 0.11μm 製程，對於振動產生的影響很敏感，因此必須在高速鐵路尚未完工運作之前，事先進行評估，避免高鐵通車後對產品良率產生嚴重的影響，以致降低產品的競爭力及獲利能力。園區內的電子公司以前雖曾委託台安工程技術顧問股份有限公司於高鐵沿線基樁預定處進行重錘落錘實驗【8】，量測估算土層波傳消散特性，但此資料無法直接用於判斷高鐵列車通過竹科段所引致振動，因此本文以台灣高鐵行經新竹科學園區廠區附近路段為三維分析案例，首先利用套裝軟體實際模擬列車移動載重於高架橋基礎所產生的動態反力，劉俊秀黃炯憲李洋傑【9】，接著結合以上分析的

(40)

成果，利用理論解析與電腦數值模擬等方式，來分析高鐵列車通過竹科段所引致之振動，並探討該振動傳至園區內，對晶圓廠產生的影響。

3-2 高速鐵路高架橋振動反應分析

台灣高鐵行經新竹科學園區廠區附近，最近距離僅僅約為285 公尺，該路段主要是由高架橋與部分隧道及路堤路段所組成。就高架橋路段而言，高鐵行車移動載重將激發高架橋結構行車方向（縱向）的擺動振動模態，高架橋結構的振動能量透過高架橋基礎傳遞至附近地盤，再經由地盤土層的應力波傳遞至園區地盤，最後引起園區廠區之振動反應。而隧道及路堤路段則是因為連接著廣大半無限土層，沒有明顯的振動模態可供激發，所產生的地盤振動反應相較於高架橋路段之地盤振動反應之比例相當小，因此於振動分析時對隧道及路堤路段的貢獻均忽略不計。

3-2-1 載重模型

高鐵THSR 新幹線列車由 12 個車廂所組成，列車的 48 個輪軸載重大小及排列方式如圖 3-1 所示，劉俊秀黃炯憲李洋傑【9】，載

(41)

重移動之速度則以台灣高鐵設計速度300km/h 計。高鐵列車車廂對振動的要求非常嚴格，最後所求出輪軸下之可動支承反力歷時曲線相當接近於常數值水平線，對高鐵橋樑而言，這表示了 48 個移動載重相當接近擬靜態載重，即載重大小與方向不隨時間改變，但作用點隨時間移動改變，若分析的重點在於橋樑下方基礎地盤的動態行為時，則將 48 個移動載重考慮為擬靜態移動載重，已足以作出準確的分析結果，因此以下的分析中，高鐵列車所引起的載重一律考慮為擬靜態移動載重。

3-2-2 結構模型

台灣高鐵行經新竹科學園區廠區附近路段，樁號為 76K+435 至 76K+785 之高架橋共有 10 跨，其中包含 7 個簡支樑型式與 3 個連續樑型式，高架橋面接近於水平，而這些高架橋的9 個墩柱高則依地形高低不同而改變。高速鐵路高架橋振動分析之結構模型如圖 3-2，劉俊秀黃炯憲李洋傑【9】，其主要組成部分說明如下： (1) 高架橋 9 個墩柱之基礎模擬為固定支承，即基礎的勁度假設為無限大。在新竹科學園區附近地盤相當堅硬之環境下，若分析標的在於高架橋基礎之動態反力，這應屬一個合理的假設。

(42)

(2) 高架橋 9 個墩柱分別以樑元素模擬，各段樑元素之性質依其實際斷面計算，如表 3-1 所示。 (3) 各墩柱頂之帽樑則依其質量以二至三根之勁度無限大樑元素模擬。 (4) 高架橋之 10 個橋跨大樑則以每根 1m 長之樑元素模擬，樑元素個數視大樑長度而定，各段樑元素之性質則依其實際斷面計算，如表 3-2 所示。由於大樑中包含一些機電設施、道版與欄杆等非結構元件，這些部分以大樑之附加質量₂_.₇₉_×₁₀3_{kg/m 處} 理。 (5) 軌道則只考慮其軸向勁度，四根鋼軌合起來以一個軸向桿件（E=2.1×1010 kg/m2，A=4×7.7×10-3 m2）處理，扣件對於軌道與道版之連接也只考慮其軸向勁度，以每單位長度之軸向勁度 18.367×106 kg/m2納入模型中。

3-2-3 動態反力歷時分析結果

由載重模型及結構模型配合套裝軟體，即可進行頻率域動態分析，劉俊秀黃炯憲李洋傑【9】，其分析的頻率範圍為 0~15Hz。在求出每一個高架橋基礎動態反力六個分量的頻率域反應後，經由FFT

(43)

轉換到時間域，就可以得到高架橋基礎之動態反力六個分量的歷時。台灣高鐵行經新竹科學園區廠區附近路段9 個墩柱之動態反力歷時，如圖 3-3~3-11 所示。其中各反力分量的座標系以 X 向為行車方向，而Z 向為垂直方向，各圖時間軸已平移使原點對準高鐵列車頭通過各該墩柱中心時刻。就各反力分量之間的比較而言，由圖3-3~3-11 可以得到以下之觀察結論： (1)各反力分量中R_x、R_z與M_y三個分量相對於其他三個分量大相當多。 (2) R_x與M_y兩個分量為前述高架橋結構的縱向擺動振動模態被高鐵行車之移動載重激發後的結果，其主要頻率大約在 3.3Hz，此為高鐵行車移動載重之第一個主要頻率。 (3) Rz分量則為高鐵行車的移動載重透過高架橋結構系統傳遞至基礎的垂直向反力，其時間變化主要是移動載重擬靜力影響線的結果，高架橋結構之振動的影響對於此分量的貢獻相對較小。整個列車通過各該墩柱的時段可以由Rz分量歷時之高原區段清楚識別。就各墩柱基礎反力之間比較而言，由圖3-3~3-11 可以得到以下之觀察結論：

(44)

(1)三個連續樑間的兩根墩柱之基礎反力比其他墩柱之基礎反力大相當多，其主要理由為三個連續樑的跨距較其他橋跨跨距（皆為30M）要大。 (2)就整個列車通過後的時段，即對應各圖中Rz分量歷時之高原區段過後之時段而言，76K+595 墩柱基礎之反力R_x與M_y兩個分量皆可看到明顯接近單頻自由振動的反應，這表示 76K+595 墩柱的縱向擺動振動模態被高鐵行車之移動載重激發至接近共振的結果。

3-3 1/3 倍八分貝頻帶圖

由於晶圓廠的設計規範，有關微動部分皆以1/3 倍八分貝頻帶圖 ( one third octave band)來表示其容忍值之標準，因此相關的分析結果皆須先轉換成 1/3 倍八分貝頻帶圖與規範進行比較。1/3 倍八分貝頻帶圖轉換相關式首先介紹如下

歷時x(t)之能譜密度函數( power spectral density)

2 n 1 k k d xx d ) T , f ( X T n 2 ) f ( Sˆ

∑

= = (3-1) 其中，Sˆ_xx(f)為微動量歷時x(t)之能譜密度函數估算值 ) T , f ( X_k 為歷時x(t)之 Fourier 轉換

(45)

T 為歷時 x(t)之延時，在本研究，歷時 x(t)只取主振動段 nd為分段數。

然後，再依均方根(root-mean-square)與能譜密度函數之間之關係，估算1/3 倍八分貝頻帶( one third octave band) 所對應之速度均方根： 2 / 1 ) ( ) ( ( ) ˆ ) ( = _⎢⎣⎡

∫

f i _⎥⎦⎤ i f xx rms high low df f S i V (3-2) 其中，flow(i)與 fhigh(i)為各 1/3 倍八分貝頻帶之上下限(參考表 3-3)。

由式(3-2)所得之速度均方根，將透過下式轉換成分貝（dB）值 o rms V i V ( ) log 20 (3-3) 其中V0取為1×10−6in/sec=2.54×10−8m/sec 由於本章分析中計算所得的振動量值，皆已經過1/3 倍八分貝頻帶圖轉換相關式計算改為以分貝（dB）值來表示其成果，因此以下就其實際轉換的程序加以說明

1. 將下限頻率 flow(i)、中心頻率為 fmid(i)、上限頻率 fhigh(i)之頻帶，

上、下限頻率範圍內劃分成若干段。 2. 以理論或有限元素法方式，求出外力作用下於各分段點頻率所產生之速度振動量值，並由其求各分段內速度振動量平均值 ) T , f ( X_k 。 3. 由X_k(f,T)及式(3-1)，可得能譜密度函數Sˆ_xx(f)，並由各分段對應的

(46)

劃分頻率帶寬，求出各分段能譜密度函數面積。

4. 由式(3-2)，求下限頻率 flow(i)至上限頻率 fhigh(i)內所有分段能譜密

度函數面積和之均方根，即為外力作用下，中心頻率為 fmid(i)時所對應之速度均方根V_rms值。 5. 最後將速度均方根V_rms經由式(3-3)的轉換，即可求得外力作用下，於中心頻率為fmid(i)頻帶所產生振動量分貝 dB 值。

3-4 地質參數判斷

本研究之地質參數主要參考竹科晶圓廠建廠之地質探勘及震測資料與高鐵公司沿線之地質探勘資料來分析判別。至於地形之變化，以有限元素法方式的分析中，將依據現有地形之實測資料來建立模型，以符合廠區附近與高鐵沿線之實際地形。

3-4-1 震測地質參數

由中基土壤技術顧問有限公司【10】及亞新工程顧問股份有限公司【11】之相關資料可了解，園區之地盤主要是由頁岩上面覆蓋礫石與砂質土壤所形成，其覆蓋層厚度約為1m~10m，地質條件良好，壓

(47)

力波速可達2000m/sec 以上，因此於分析時地盤可假設為二層土壤；第一層厚為5m，其壓力波速為 1000m/sec。第二層為彈性半空間，其壓力波速為 2000m/sec。又依據宏偉工程顧問有限公司茂矽建廠地質探勘與震測資料【12】得知，其土壤柏松比υ=0.4，及土壤單位重 ρ=2.1ton/m3，所以第一層土壤之剪力波速為 408m/sec，第二層土壤剪力波速為816m/sec。因為分析對象為高鐵列車通過高架橋所產生之振動，其在土壤所產生之變形為小變形，所以考慮土層之動態參數時，皆以小應變之動態參數考慮，又由於土壤之應變小，故土壤之阻尼比只取 1%。由地質探勘資料所得，分析時使用之土壤動態參數表示於表3-4。

3-4-2 重錘實驗之等值地質參數

台安工程技術顧問公司曾受委託執行重錘落錘實驗以求得高鐵沿線至晶圓廠房之地盤振動之衰減，此實驗所假設之衰減公式為振動全由表面波(Surface Waves)所引起，並假設地盤為均質土壤，因此其振幅衰減關係可表示為 ) r r ( 2 1 1 2 e 2 1 r r v v = −α − (3-4) 其中

(48)

v1、v2分別為點1、點 2 之振幅 r1、r2分別為點1、點 2 至振源之距離 2 r r1 可視為幾何衰減 ) r r (2 1 e−α − 可視為土壤阻尼衰減。台安工程顧問試驗報告並以4Hz、8Hz、16Hz、32Hz 為例，建議式(3-4) 之衰減常數α=0.01。由振幅衰減關係式及試驗報告建議之結果，可用理想化模式回歸 α 值，決定出與重錘實驗等值的地質參數，其計算方式如下 (1)以理想化模式計算施加單位力時，距施力點不同距離各點理論振幅解析解。解析解之分析方法請參考 Gin-Show Liou 【13-16】，分析時所假設之各種參數條件表示於表 3-5。 (2) 由各點解析解及已知與原點距離，以式(3-4)求出各點理論衰減常數α_i，並由α_i計算所有點位理論衰減常數平均值α_avg。 (3) 將不同假設參數計算所得之理論衰減常數平均值α_avg，與試驗建議之衰減常數 α 做比較，可回歸出分析時等值的地質參數。上述的回歸程序中，若先將理論振幅解析解 v1、v2 由式(3-3)轉換成dB 值後，代入式(3-4)可得

(49)

) ( 10 ln 20 log 20 2 1 2 1 1 2 r r r r dB dB − − − = α (3-5) 其中 1 dB 、dB2分別為點 1、點 2 之分貝值 r1、r2分別為點1、點 2 至振源之距離由 (3-5)式以 dB 值回歸 α 值之結果，建議土壤之剪力波速為 408m/sec 及阻尼比ξ=0.05，分析所得之土壤參數表示於表 3-6。回歸所得以分貝dB 值表示之地盤振動衰減關係如圖 3-12 所示，其中橫軸為距施力點之距離，圖 3-12(a)中 Y 方向為與高鐵列車行駛方向垂直之方向的振動，圖3-12(b)中 Z 方向為垂直方向之振動。Y、 Z 或 α 括號中之數值為阻尼比(ξ) 。

3-5 理想化分析

3-5-1 理想化分析模式

首先以理想化之3-D 模式分析高鐵列車通過時之微振動。其中理想化之3-D 模式如圖 3-13 所示，基本假設為(1)地盤為水平層狀土壤； (2)土壤力學性質為線彈性；(3)基礎為表面基礎，而理想化模式解析

(50)

解之分析方法請參考Gin-Show Liou【13-16】。由表3-4 及表 3-6 之土壤地質參數可得二種地盤，這二種地盤的第一種為一5m 厚土層覆蓋於彈性半空間之地盤，第二種為彈性半空間地盤。又由於橋墩椿基礎皆座落於堅硬地盤上，因此亦考慮表 3-4 中無第一層較軟土壤的情況下之彈性半空間地盤。因此在理想化模式分析中，共考慮了三種地盤參數。由於振動之求解方法為在頻率域( Frequency Domain)中求解，圖 3-3~3-11 中較重要之外力 Rx、Rz、My三方向做Fourier Transform，可得各頻率外力之大小如圖3-14~3-16。圖中 Pier1~9 所代表之橋墩位置為76k+465、76k+495、……76k+755 等 9 個橋墩。

3-5-2 理想化模式分析結果

圖 3-17~3-19 為圖 3-14~3-16 之外力作用下，距離橋墩 300m 處三個方向(二個水平方向，一個垂直方向)振動之 dB 值，此三圖的地盤土壤參數為剪力波速Cs=816m/sec，及阻尼比 ξ=0.01。圖 3-20~3-22 為圖 3-14~3-16 外力作用下，假設二層土壤地盤，距離橋墩 300m 處三個方向振動之 dB 值。此二層地盤之參數分別為上層土壤厚 5m， Cs =408m/sec，下層彈性半空間之 Cs=816m/sec，兩層土壤之阻尼比

(51)

皆為 ξ=0.01，如表 3-4 所示。圖 3-23~3-25 為圖 3-14~3-16 外力作用下，假設地盤土壤參數Cs=408m/sec 及阻尼比 ξ=0.05，距離橋墩 300m 位置三個方向振動之dB 值，如表 3-6 所示。其中必須注意圖 3-17~3-25 之結果，均假設單一橋墩外力作用下之振動dB 值，然而於實際分析時必須模擬9 個橋墩之外力在高鐵列車通過時間延遲之效應下，作用於地盤所產生之振動dB 值，而此效應將與地形起伏變化之效應同時考慮。在單橋墩外力作用的情況下，距離高鐵最近廠房之地盤最大振動量由圖3-20 可知，Y 方向(與高架橋垂直之水平方向)之最大振動值為 37dB，其發生在中心頻率 6.3Hz 處。由圖 3-21 得知，Z 方向(垂直方向)之最大振動值為 42dB，亦發生在中心頻率 6.3Hz 處。由圖 3-22 得知，X 方向(與高架橋平行之水平方向)之最大振動值為 48dB，其發生在中心頻率10Hz 處。

3-6 有限元素法分析

3-6-1 有限元素法分析模式

由於地盤為半無限空間區域，用有限元素法分析，首先必須先決

(52)

定分析區域( Domain)。有限元素分析模式之分析區域如圖 3-26，橋墩之相對位置如圖3-27 所示。根據改變分析區域大小之有限元素法數值實驗結果【17】，得知用 1/2 模式所需之區域大小；頻率 3Hz 以下時為 880mX880mX445m 如圖3-28，頻率介於 3.5Hz~6.0Hz 之間使用 700mX500mX185m 如圖 3-29，頻率大於 6.0Hz 以上時使用 500mX400mX115m 如圖 3-30，並加上邊界使用無限元素，其數值結果即可收斂。利用有限元素法解波傳問題時，元素種類及尺寸之大小亦會影響數值的精度。有限元素法三維分析所使用的元素種類，近域的分析採用三維、二十個節點的實體元素C3D20 如圖 3-31 所示，在遠域的土壤採用三維、十二個節點的無限元素CIN3D12R 如圖 3-32 所示，來模擬吸能邊界。就元素尺寸大小的選擇方面，基本上如果波長愈短，元素尺寸應愈小，而根據數值實驗，使用以下之元素尺寸即可達到所要求之精度。頻率 3Hz 以下時，使用元素尺寸為 40mX40mX40m，頻率介於 3.5Hz~6.0Hz 之間使用 20mX20mX20m 之元素，頻率大於 6Hz 以上時，使用 10mX10mX10m 之元素。有限元素分析使用之程式為商用 ABAQUS 程式【18】。為証明前述之模式區域大小及元素尺寸大小可行，圖3-33 及 3-34 為比較距離高鐵 300m 處 ABAQUS 之結果與解析解之結果，解析解之分析方

(53)

法請參考Gin-Show Liou【13-16】。圖 3-33 及 3-34 中 Ab 之曲線為 ABAQUS 之結果，I 之曲線為解析解之結果。圖 3-33 及 3-34 之外力為 Pier4(76k+595)之外力，且假設無地型變化及土壤剪力波速Cs=816m/sec。由圖3-33 及 3-34 可知，前述所選之區域大小及元素尺寸合理。

3-6-2 有限元素法分析結果

由理想化模式之結果可知，地盤為二層土壤時之振動值較大，因此在有限元素法之分析模式中，地盤以二層土壤模擬，上層厚為5m，剪力波速Cs=408m/sec，下層為彈性半空間，Cs=816m/sec。兩層土壤之阻尼比(ξ)皆為 0.01。理想化模式分析結果皆為單一 Pier 之外力作用下且假設無地形變化，距離高鐵高架橋 300m 處之振動 dB 值，但實際情形為高鐵列車通過時 Pier1 到 Pier9 將依次產生振動外力，而且新竹科學園區為高低起伏的丘陵地，地形產生的影響亦不可輕忽，因此在實際案例之分析中，將依據現有地形之實測資料及 Pier 位置來建立模型，以符合園區與高鐵沿線之實際情況，並考慮每一Pier 之外力依高鐵列車通過 (時速 300km/h)之時間延遲作用在地盤上，求取晶圓廠位置處(距高鐵

(54)

高架橋約300m)之地盤振動。

時間延遲作用下外力的計算方式如表 3-7 及表 3-8 所示；其中表 3-7 為假設高鐵列車時速為 300km/h，車行方向由 Pier1 到 Pier9 時，列車到達各Pier 的時間 t，將各單一 Pier 的外力與列車到達各 Pier 的時間 t，經由如表 3-8 的計算，就可以求得單一 Pier 在時間延遲作用下的外力。表3-8 中ω =2πf 為外力的頻率，r、i、R、I 分別為頻率域外力的實部值與虛部值。在同時考慮9 個橋墩外力依列車行進之時間延遲作用之情況下，晶圓廠之地盤最大振動量由圖3-35 可知，在中心頻率 10Hz 處，X、 Y、Z 三個方向之最大振動值皆約為 40dB，在中心頻率 6.3Hz 處，三個方向之最大振動值皆約為35dB，在中心頻率 3.15Hz 處，在 X 方向之最大振動值約為36db。

(55)

第四章結論與建議

本文利用有限元素法為工具，建立各種不同的數值模型，來模擬分析震波於介質中的傳遞行為，以下將對本文分析所得到的結論與建議加以說明。

4-1 結論

1. 以 ABAQUS 分析波傳問題，所有的模擬與計算均是在電腦上完成，由於受限於系統記憶體及CPU 處理能力的不足，因此在決定模型分析範圍及網格分割大小時必須審慎加以考慮，本文以羅博智【17】無因次化數值實驗的結果，來決定分析案例模型模擬的關鍵參數，經過證實可以節省運算資源而得到令人滿意的結果。 2. 實際案例高架橋基礎的動態反力，因為受到高鐵行車移動載重的影響，所以在Rx、Rz 與 My（X 為行車方向、Z 為垂直方向、Y 為與高架橋垂直之水平方向）三個方向上的反力分量會較其他分量要大，此應屬合理現象，而當高架橋的跨距較大時，墩柱之基礎反力也會明顯的增大。 3. 廠區附近環境振動量的大小，通常都是由儀器現場量測所得，但

(56)

是若遇到必須事先評估振動量的案例，可以利用理論或有限元素法方式，先求出外力作用下所產生的位移或速度值，再經由與本文相同的程序，將其轉換成分貝（dB）值後與規範進行比較。 4. 由振動量的分析所得到的結論，在單橋墩外力作用下、土壤介質為兩層土壤，於晶圓廠地盤的X 方向會產生最大振動量值 48dB，若考慮9 個橋墩的外力依列車行進之時間延遲作用時，晶圓廠處的地盤在X、Y、Z 三個方向，均會產生約為 40dB 的最大振動值。

4-2 建議

1. ABAQUS 數值模型的建立相當繁複，分析時往往需要配合修改大量的輸入條件求解，為了節省分析時間與減少錯誤發生的機率，建議以篆寫程式的方式來產生ABAQUS 批次輸入檔，並配合ABAQUS/CAE 分析模組的前處理產生網格及偵錯的功能來輔助分析，而為了提昇計算能力及數值結果的準確性，則建議利用國家高速電腦中心的設備來幫助研究的進行。 2. 本研究三維分析模型的土壤參數以及地形變化，主要是參考實際量測的資料來分析判別，有鑑於現地實驗往往所費不貲

(57)

而且結果取得不易，因此建議分析時若無相關資料可供佐證，可將土壤介質視為遵守彈性原則之微小變形，阻尼比以小應變之動態參數考慮取1%來分析，而當分析區域廣大或地形起伏不明顯時，建議模型地表高程以平面來代替高低起伏之模擬。 3. 本研究分析振動量分貝（dB）值的大小，是於頻帶寬內取點，計算外力所產生的速度值平均而得，若是取的點位過少，則容易忽略外力的峰值而產生誤差，因此建議於外力極值附近增加計算的點位，使分析的成果更能反映出真實的現象。 4. 建議日後高速鐵路通車後，可至園區附近進行振動量的現地量測，並將實際量測的結果與本文互相比較，藉以研判以數值方法來分析波傳問題的精確性。

高鐵通過竹科之波傳研究

國 立 交 通 大 學

土 木 工 程 研 究 所

碩 士 論 文