行政院國家科學委員會專題研究計畫 成果報告
總計畫暨子計畫:地聲與電磁波監測方法之應用（III）
計畫類別： 整合型計畫 計畫編號： NSC92-2625-Z-002-024- 執行期間： 92 年 08 月 01 日至 93 年 07 月 31 日 執行單位： 國立臺灣大學水工試驗所 計畫主持人： 劉格非 報告類型： 完整報告 處理方式： 本計畫可公開查詢中 華 民 國 93 年 12 月 17 日
行政院國家科學委員會補助專題研究計畫成果報告
土石流監測與預報系統之研究（三）
Debris Flow Monitoring and Forcasting System (III)
計畫類別：整合型計畫
計畫編號：NSC 92－2625－Z －002－ 024 －
執行期間： 92 年 8 月 1 日至 93 年 7 月 31 日
計畫主持人：劉格非
計畫參與人員： 李欣輯、黃登賢、陳盈蓁、張婉真
成果報告類型(依經費核定清單規定繳交)：完整報告
處理方式： 立即可公開查詢
執行單位：國立台灣大學土木工程系
中 華 民 國 93 年 10 月 31 日
目錄
第一章 前言 3 第二章 土石流監測與預報系統 10 第三章 地聲探測器 17 第四章 與電磁波探測器 51 第五章 結論 65第一章 前言
1.1 整合計畫背景與緣由 土石流是一個非常複雜的問題，發生原因很多。從開始有單一石頭顆粒運動到發生大 規模土石流動行為之過程中，其相關力學機制不但多而且複雜，多數力學機制到目前仍停 在學理探索階段，因此土石流之發生或災害預報必須結合多種不同資料，再配上力學與統 計才有可能得到較高百分比之預報成功率。因此將任何正在研究、已在利用或將來會研發 的設備整合成一套系統是目前極重要的工作。國外的例子如美國地調所（USGS）在聖海倫 火山(Mt. St.Hellen)地區結合地聲、紅外線、航照、水壓系統而成為火山土石流監測系統； 義大利在 Italian Atps 用地聲、攝影機與超聲波結合之系統來監測泥石流；而日本京都 大學防災研究所在澤田教授負責的觀測所（在 Gifu）更是將所有的儀器均整合並試用於一 地點，因此該所在土石流之監測技術上可說是領先全球。由這些例子可瞭解一個集中並整 合科技於現場之試驗場所對研究是非常重要的。 目前國內土石流之研究單位與學者雖然多，但方向與重點均不相同，尤其很少有完整 並公開之試驗資料可供交流討論。而理論又缺少現場資料來證明，儀器亦缺乏現場應用來 考驗，因此為了推動本土研究之進步，亟需建立一個供大家使用之平台。本群體計畫就是 要朝此方向邁出第一步，先經由儀器之整合建立一個共享之監測中心，再由此刺激土石流 研究方向之修正，最後建立所有研究成果之驗證標準。而要達成此一系統整合之目標，必 須有一群人將過去研究成果朝此方向修正應用，為最終之目標踏出第一步，本研究群因而 產生。 1.2 整合計畫總體目標 本整合計畫之總體目標預計有下列五項 1.整合國科會過去土石流監測之相關研究 2.研究並分析各類預測設備之能力與應用範圍 3.發展一套結合各種土石流預測儀器之系統 4.由土石流流變特性推估土石流之流動特性 以下就每一項分別作一說明： 1.2.1 整合國科會過去之相關研究 國科會過去的土石流相關研究均集中在一個群，其中有做基礎理論、有做應用、也有做現場，因為目標不同，因此雖有交集，但其相互間之關係並不強。但各計畫經過多年研 究，研究成果均大幅成長，本研究群即擬將過去研究群中從事監測與預警之部份成果結合， 在結合過程中，可以推動各成果間之橋樑以建立一個真正整合的系統。有了整合之成果， 才能突顯出國科會過去之研究成效。同時整合之後再加入使用者界面，來拉近理論與實務 間之鴻溝，以確實彰顯防災研究之最後目標。 1.2.2 研究並分析各類預測設備（測量元件）之能力與應用範圍 土石流之監測與預警可以分成四個步驟來談。監為監視；測為測量；預為預報；而警 為警告。許多測量儀器事實上是用來監視與記錄某個物理量之改變。例如許多壓力計是紀 錄儀器受外界作用時之電流量或變形值，再加上相關理論即可推為壓力之變化；又例如黑 白攝影機在拍攝時是記錄環境中各點之亮度值，人類再利用亮度值結合成之形狀去判斷物 理現象，如石頭形狀大小、石頭速度。因此我們在使用儀器時是在監視某一種資料，再利 用不同理論去應用監視資料時就產生了測量的行為。同樣原理下，測量數據出來後又可利 用不同理論將之應用於不同物理現象，進而預測物理現象之變化，例如用地聲探測器測得 聲音之大小變化後，用聲波理論轉為土石流之距離進而預報土石流抵達時間，或用邊坡形 變預報邊坡不穩之時間，並且運用流變特性之量測，評估土石流的運動行為與活動區域。 有了這種預報與評估再結合該地區之特性，如地形、行政支援、人口分佈、產業結構、經 濟分析、法律責任等考量後，就可以發佈警報，因此每一步驟在不同應考量之因素下都有 不同的參考數據與應用理論。 因此儀器所得數據在監測與預警之過程中，可配合不同理論而產生多種不同用途與解 釋方法。本研究之另一目的就是將各儀器所能適用之情況分清楚，在結合不同理論下，其 產生預報結果之精度也要定義清楚，尤其在預報之時間長短和空間上的精度列為研究重點。 1.2.3 發展一套結合各種土石流預測儀器之系統 對本系統而言，元件又可分為測量元件、通訊元件與控制元件三類。測量元件是用來 測量土石流相關之信號；通訊元件是用來聯繫各測量儀器與訊號傳輸；而控制元件是用來 做自動控制、自動判別。本群體計畫之總計畫偏重通訊元件與控制元件，各子計畫則負責 測量元件。而測量元件又可分為許多測量儀器。有了上節各儀器之功能界定之後，我們可 以將預報儀器區分為間接測量、遙測與接觸式測量三類儀器，而此三類儀器之結合即成為
完整測量元件，因為各類儀器是組成元件之一部份，我們可稱之為子元件。 測量元件由上述定義可分為三類，而各測量元件因測量方法各有其優缺點，相互間之 關係如圖一。 第一類：間接元件 這一類元件包括雨量測量與地下水電導度（EC），這一類資料嚴格說來與土石流並無絕 對關係，在其它條件配合下，即可對土石流之發生做機率上之預測，因此這類資料多用於 長期預報，如數小時到一天之預測時間，而這類元件在本群中由子計畫一負責（表一請見 3.1 節）。 第二類：遙測元件 這一類元件是在土石流發生前後由於土石流造成之地形、地貌改變或其造成之聲響、 影像等資料，用遙測方式來取得，這一類資料可提供土石流中、長期之預測，如數十分鐘 到數小時。 這種元件包括衛星照片、航照、攝影影像、電磁波、地聲等，本部份由第二與第四子 計畫負責。 第三類：接觸式元件 這一類是在土石流通過河道時來感測其產生之力量所得數據，包括鋼索、土壓計、水 壓計、傾斜儀、水位計、流速計、加速計等，這種子元件可依其裝設位置與預警對象位置 不同而提供約 30 分鐘以內之短期預報，這部份由第三子計畫負責。 間接測量方法一般較易獲得資料，因為非直接測得土石流之發生，因此其準確性較低。 但此較適合作長期預報。遙測式元件配合理論可作中期（如數小時）預報，且其精度高於 間接測量。而接觸式的精度最高，但一般認為只適合用來做短期預報。因此我們設計之系 統中測量元件之關係如圖一中之箭頭，而系統就靠通訊元件來串連所有測量元件與分散在 各地之資訊，最後經過經濟、法律之考量建立警告條件。
圖 1.1、系統架構 遙測方法 間接方法 接觸 式方法 通訊警報 經濟法 律考量 通報系統 而 2000 年八月在台北舉辦之「第二屆國際土石流研討會」中的一場「現場經驗與學術 研究」討論中也知道，要刺激學術研究的進步與創新，尤其需要應用者提出使用時之需求 與現象，本研究群希望在研究的過程中將我們的測量資料與應用中產生的新的研究題材不 斷提供其它研究土石流的學者，希望藉此建立土石流研究學界的一個新的、有效的互動關 係。 1.2.4 建立一個土石流相關預測科技之展示與訓練中心 結合了所有測量元件，目的就是要推廣給使用者，這其中包括公家機關、民眾與百姓。 本群擬利用國科會歷年在神木村附近累積之儀器建立一個測試區，所有過去與未來之儀器 均將架設於此，一方面測試各元件間之互動與比較，一方面方便展示。 而這個神木村測量區也就成為當然之展示區，但因該地區交通不方便，因此我們將以 虛擬實境之方式將神木村展示區在電腦上重現出來，一旦採虛擬 3D 互動式之展示後，則只 要有軟體與網路，展示中心可以搬到任何地區。 最後在推廣之過程與成果整理之過程中，我們將對系統與其中的每一個元件建立教育 與訓練之材料，使得這些方法能讓潛在使用者瞭解而自行決定要用那些方法。
三、研究群之分工架構、中心理念與互動關係 1.3 分工架構 依前一節所述之目的，本研究群之分工如圖二，而地聲與電磁波之應用子計畫因為研 究者為總計畫主持人，因此併入總計畫 子計畫五 子計畫一 子計畫二 子計畫四 子計畫三
圖 1.2、分工架構
長期預報 中期預報 短期預報 雨 量 與 地 下水 EC 監 測 方 法 之 研究 地聲與電磁 波監測方法 之應用 機械視覺 監測土石 流之研究 土石流衝擊 力與環境因 子監測技術 之研究 土 石 流 流 變 特 性 量 測 與 運 動 特 性 之 研究 總計畫—系統建立 (包括通訊、資料彙整、警報) 計 畫 項 目 主 持 人 服 務 單 位 系 所 職 稱 計 畫 名 稱 總計畫 劉格非 台灣大學土木系 教授 土石流監測與預報系統之研究 子計畫一 范正成 台灣大學農工系 教授 雨量與地下水電導度(EC)監測 方法之研究 子計畫二 劉格非 台灣大學土木系 教授 地聲與電磁波監測方法之應用 子計畫三 林炳森 中興大學土木系 教授 土石流衝擊力與環境因子監測 技術之研究 子計畫四 張守陽 台北科技大學土木系 教授 機械視覺監測土石流之研究原則上子計畫第一年均是做室內試驗，將過去之成果做一整理與驗證，同時測試與規 劃元件系統以備第二年現場工作。此部份之工作以一年來做，時間稍微不夠，因此會延續 到第二年，但在第一年期間，各元件與系統之連線規劃必須完成。第二年主要為測試現場 系統。第三、四年元件分別完成網狀佈設研究、元件功能確定與使用者界面三大工作，並 且在發生監測外引入運動分析，以利空間預警之決定，而總計畫在第三、四年之重點則在 通報元件，教育訓練中心之工作，警告部份會進行研究，但其成果較難評估。 在各測量子元件之研究過程中，本研究群將以系統發展之概念來發展各元件，期望將 來發展出來之系統可以比擬為一台電腦，而各子元件一旦發展成子系統，就好像電腦中之 晶板，可以隨插即用來增加系統之不同功能。因此所有元件之研究將遵循圖三之概念。 監視系統 理論 測量數據 測量項目與精度 理論 預測項目與能力 單一元件系統 使用手冊、使用者界面 理論 圖 1.3、元件系統研究構想 同時在各測量研究之過程中，本研究群也將配合 ISO 之觀念，定出未來任何測量元件 要發展到實用並能整合進入本系統所需要經過之標準步驟。例如我們可以訂定任何一個現 地使用之儀器應該進行（1）室內實驗（2）現地實驗（3）元件系統現地使用性（4）元件 輸出規格標準測試等四個步驟，將來儀器若能通過這四個步驟，並提出書面資料，則大家 可以接受這是一套已經由理論發展而成可應用於現場之儀器。
子計畫間之互動關係 總計畫負責結合各類測量資料並集結各地資料來做通盤研判，並將完全系統展現給使 用者，但對每一元件如何由其測量之數值轉成物理值(例如 Load cell 是測電壓值，但需轉成 位移值或壓力值才有用)，及如何由物理值轉成土石流之預報(例如雨量多少時會發生土石流 及地聲多大時代表土石流還有多遠等資訊)，均是各子計畫研究多年之成果，因此總計畫需 各子計畫提供相關之研究成果，才能做綜合研判，做到如圖一之長、中、短期預報互動之 模式。 而監測一個事件之過程中，系統之規劃是在間接資料發出初步警訊時（如雨量超過某 一程度），即將間接元件與遙測元件之監測頻率加快，若遙測元件也顯示土石流之可能，則 再把接觸元件之監測頻率也加快，如此才能在有限資源下完整監測一場土石流之發生，而 當接觸元件確認土石流通過後，則依次將遙測元件與間接元件之監測頻率降回正常。 而子計畫則靠總計畫將資料收集，並對現場做互動式遙控，同時需經過總計畫才能得 到與其它理論之比較。因此若缺了任一子計畫，這群體計畫成果將缺乏彈性且顯得不完整。
第二章 土石流監測與預報系統
2.1 前言 土石流監測是一個全世界正在進行之研究，世界上有名的監測系統如美國聖海倫火山 監測網、大陸九寨溝土石流監測系統，歐洲在阿爾卑斯山區由瑞士與義大利等國建造之監 測系統等都算是較成功且有監測成果的重要地點，不過上述系統也都只做到監測，但未到 測警之階段。本研究群即是要建立一套示範性之整合性之監測預警系統。 一般會發生土石流之地點，地質定破碎，地形也不多容易變動，因此任何有形之線路 如電源線、網路線就常常會在災害時期被破壞而造成訊號中斷，最重要的信號因此無法傳 到指揮中心，因此本系統之設計為無線系統。無線系統目前可供考慮的是大哥大(GSM, GPRS, PDA,等)、無線網域系統、無線對講機或衛星通訊系統。衛星系統之造價與維護均極昂貴， 本研究暫不考慮。無線對講機僅適合短距離使用，且不易自動化，因此本研究採大哥大為 遠距傳訊，無線網域為短距傳訊。目前通訊採先將神木測站之所有資料以無線區域網路傳 到神木國小，到了神木國小之後，平常有電、網路也通時，走一般網路，訊號可傳到任何 地方，上網即可接收。但在神木國小斷電或網路不通時，訊號將自神木國小由 GSM（大哥 大系統）傳送。因此監測人員若所在地點沒有網路可用，也可用 GSM 來傳訊，甚至人員在 移動中也可透過此系統對現場做監測，因此可說是移動式遙控中心。 2.2 監測系統 目前監測資料也可上網觀看，網址http://163.22.163.245/ ，入口密碼為Hello。監測資料 為公開資料，所有人都可上網去看，一進入網頁可以看到首頁為 只要輸入密碼「Hello」，即可進入系統觀看所有儀器的資料，系統會整合所有監測儀器之 信號為圖 2.1 之網頁圖 2.1、監控主畫面 監控畫面除展示當時所監測到的資料外，也可以按不同理論將監測資料轉換成土石流 相關資訊，如土石流發生機率、土石流可能發生距離、土石流通過時間和土石流流速等， 在監測畫面右方有轉換鈕，使用者只要點選即可。目前神木測站之相關硬體設備見表 2.1： 名稱 單位 數 量 1.雨量計 台 2 2.UPS 台 1 3.投射燈光源 台 2 4.CCD 台 2 5.Geophone(地聲探測器) 台 3 6.無線 Modem(GPRS) 台 1 7.鋼索檢知器 台 1 11.現場小屋 個 2 12.工業級中央處理器 台 1 14.電子式地下水壓計 組 1 15 GSM 組 1 16 太陽能電池 片 10 17 光遮斷器 組 1 表 2.1、神木測站之相關硬體
而且系統因為加了十片 80 瓦太陽能電池，因此即使斷電時，若 CCD 不需打光，系統可在沒 有陽光的情形下維持 72 小時的高頻率紀錄運轉，若需打光，則可支持連續 3 小時的打光。 各儀器的監測資料都可點選進入而看到，例如若想看三個地聲的資料，可點選地聲旁的 More 按鈕，即可看到以下畫面，包含三個地聲的資料與地點，畫面分別可以顯示上、中、下游 三個地聲的位置與接收之訊號，左方為原始時間序列資料，右方為左方視窗經過快速傅立 葉轉換後之結果，如果在 20-100Hz 的訊號強度到達平時雜訊強度十倍時，就會有自動警訊 「偵測到土石流訊號」出現，並啟動攝影機來錄影。 圖 2.2 地聲探測器監測畫面 又例如想看上下游的雨量資料，就可以點選雨量旁的旁的 More 按鈕，即可看到圖 2.3 圖 2.3 上下游雨量監測資料畫面 而若按分析紐，則依照子計畫理論分析出的土石流發生機率圖就會出現如圖 2.4
陳 有 蘭 溪 流 域 之 樣 本 溪 流 0 5 0 1 0 0 1 5 0 2 0 0 2 5 0 3 0 0 3 5 0 4 0 0 4 5 0 0 2 0 4 0 6 0 8 0 1 0 0 1 2 0 1 4 0 有 效 降 雨 時 間 ( h r ) 有效 累 積雨 量 (mm ) 發 生 土 石 流 9 9 . 5 0 % 9 0 % 8 0 % 7 0 % 6 0 % 5 0 % 4 3 . 6 % 臨 界 降 雨 線 3 0 % 2 0 % 1 0 % 0 . 5 % 圖 2.4 利用雨量來分析土石流發生機率 其他儀器有一樣的展示與分析功能，例如影像部分，就可以利用子計畫五來分析土石流相 關資訊，例如按了影像分析鈕並選擇「Pattern Analysis」後，分析資料與影像會顯示如 圖 Flow velocity = m/s There is 89% chance
This is (debris, water, mud) Flow Flow density
圖 2.5 利用影像來分析土石流資訊
為了能讓各儀器間有足夠的互動與相互節省資源，目前各儀器間有相互驅動的機制，還 有儀器退場機制，詳細如下圖
Critical warning Stage：All devices switched to high frequency recording mode
Preliminary Warning Stage：Geophones, Microwave sensor switched to high frequency recording mode Normal low frequency recording mode
Rain gage ERR>0.5 or > 3mm/5min; Ground water level exceeds ground surface Geophone: 40-80 Hz data diminish
Water level lowered to normal
Geophone: 40~80Hz amplitude is large ; Pore pressure changes abruptly;
Water level increases abruptly(>30cm/5min). Load cells reading becomes steady
圖 2.3 、設備互動機制
裝置好的主現場如圖 2.4，照片可以顯示整個系統安裝完畢。
圖 2.5 照片顯示太陽能面板裝置方式，十片面板形成一個長廊形。
第三章 地聲探測器
3.1 前言 台灣因地形陡峭且地形破碎，經過多年山坡地開發及九二一地震後，造成山坡地的土 石鬆動。台灣地區年平均降雨量為 2500 公厘，約為全球平均降雨量的二倍，每當颱風及豪 雨，大量的雨水與鬆動的土石混合形成土石流，常常發生土石流災害，造成人員傷亡及金 錢損失；如：桃芝颱風。 一般工法上都設置一些抑制土石流動的結構物；然而土石流是大面積的災害，在無法 有效抑制土石流的發生與流動時，預測土石流發生位置及流動方向就是一個重要的課題。 若能有效預測土石流發生位置及流動方向，便能在土石流發生或開始流動時，先行疏散下 游居民，以減少人員傷亡及金錢損失。 土石流在上游流動時，在下游會聽到低沉的聲音；此現象是在固體介質中，能量以波 的形式傳遞，此即固體波；當固體波傳遞到空氣時，固體波會產生聲波，而從固體波傳遞 到空氣中的聲波能量，只有固體波的能量的一小部分，絕大多數的能量都在固體中以固體 波傳遞。 聲波是一個定義的名詞，在空氣中傳遞的聲音屬於縱波；縱波是介質震動方向與波傳 方向平行。人類聽到的聲音是藉著空氣中分子震動，影響耳膜震動，在經由神經傳遞到大 腦的聽覺區塊，人類才聽到聲音；水中的聲納技術也是藉著縱波來傳遞能量的。在空氣及 水中，其抗剪強度非常低，都是藉著介質不停的壓縮與伸展來傳遞能量；但在固體介質中， 固體的抗剪強度增強，除了縱波之外，還有橫波；橫波是介質震動方向與波傳方向垂直， 因此我們所稱的地聲，是包含了縱波與橫波。 3.2 前人研究 固體彈性波傳理論已發展許多年。在固體波傳理論之外，Biot(1956a, 1956b)發展飽 和孔隙彈性介質波傳理論，孔隙中液體為可壓縮黏性流體。其波傳理論分為兩部分：在低 頻時，液體為 Poiseuille flow；在高頻時，由特徵頻率及四個無因次參數來描述。由 Biot 的波傳理論解得兩個縱波及一個橫波，並由 Hovem and Ingram(1979)利用實驗分別得出這 三個波傳現象。利用地聲探測器(geophone)研究土石流地聲特性方面。劉與李(1999)得實驗結果，其 訊號分析是利用快速傅利葉轉換(Fast Fourier Transform)將時間域訊號轉換到頻率域， 得出頻率與能量之間的關係。李並分析出土石滑動頻率分佈在 20~35Hz 之間，石頭滾動及 相互碰撞頻率分佈在 40~60Hz 之間。謝等人(2000)得實驗結果，利用 Gabor Transform 將
時間域訊號轉換到時間-頻率域，得出時間、頻率與能量三者之間的關係。其實驗分兩部分： 第一部分為摩擦與自由落體實驗，由實驗結果得知土石材料間相互摩擦的頻率分佈在 20~80Hz 之間，而土石材料間相互碰撞的頻率分布在 800Hz 以內，但大多集中在低頻。第 二部分為土石流渠道實驗，由實驗結果得知，礫石型土石流地聲頻率範圍分佈在 20~200Hz 之間。Arattano(2003)將一系列 ground vibration detector 設置在野外，當土石流流過 有高差的河床時，利用訊號分析土石流前鋒速度。 3.3 基本理論
圖 3.1 微小立方體受力及座標
z
x
y
Δy
Δx
Δz
zz τ zz τ yx τ yx τ zx τ (一) 推導介質受力的控制方程式 令平面 xy 為地表，如圖： zx τ 假設此介質為線性彈性、均質(homogeneous)且等向(isotropic)的彈性體。考慮介質中微 小立方體受力，其受力包含遠距力及表面力，遠距力如：重力、電磁力，表面力如：應力。 利用牛頓第二運動定律，推導介質受力的控制方程式。在 x 及 y 方向不受重力影響， 在 z 方向由於介質所受重力已經作用且平衡，介質的運動純粹是受外力影響。以彈簧位移 為例：吊在天花板的彈簧，彈簧底部綁一個鐵球，彈簧會受鐵球重力而伸長，當彈簧靜止 時，施一外力將鐵球向下拉一段距離後放手，此時考慮的是外力所造成的彈簧位移，而不 包含鐵球重力所造成的彈簧位移，因此不考慮重力。 令ρ=密度(kg/ 3) m ij τ =介質中微小立方體所受的應力(F/ 2) mvG=速度向量(m/s)，其中v1、v2、v3分別為 x、y、z 方向的速度分量 u =位移向量(G m)，其中 u 、 、 分別為 x、y、z 方向的位移分量 v w M =質量(kg) 首先考慮 x 方向受力情形，由牛頓第二運動定律：
(
Mv1)
dt d F x = ∑ ⇒[τyx(x,y+∆y,z,t)−τyx(x,y,z,t)]∆y∆z+[τyx(x,y+∆y,z,t)−τyx(x,y,z,t)]∆x∆z +[τyx(x,y+∆y,z,t)−τyx(x,y,z,t)]∆x∆y= ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ + dt dv dt d v 1 1 ρ ρ z y x∆ ∆ ∆ (3.1) 對(3-1)式作尺度分析： 令u=εU，ε
是 x 方向位移量的尺度大小，U 是無因次參數 LX x= ， 是 x 方向的尺度大小，L X 是無因次參數 T t t= ， t 是時間的尺度大小，T 是無因次參數 ρ ρ ρ = ˆ ，ρ 是密度的尺度大小， ρˆ 是無因次參數 代入(3-1)式等號右邊括弧內兩項，得： ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛ dt dv dt d v o 1 1 ρ ρ ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛ ∆ ∆ ∆∆ ∆ ∆ ∆ = ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛ = t t t t o dt u d dt d dt du o ε ρ ρ ε ρ ρ 2 2 ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∆ = ρ ρ o (3.2) 由膨脹模數 B (bulk modulus)： B ρ ρ ∆ = P P 是單位壓力，ρ是密度，∆ 是受壓後密度的變化量。 ρ 鐵的膨脹模數 B 9 2 10 160 m N × = 水的膨脹模數 B 9 2 10 2 . 2 m N × = 本研究的問題的介質是彈性體，其密度約 1.45g/cm^3，根據此密度大小找尋相當的膨脹模 數，根據此密度所對應的膨脹模數約為 B 1 109 2 m N × = ，所以(3.2)式：1 << ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∆ ρ ρ o 因此忽略(3.1)式等號右邊括弧內第一項；在(3.1)式等號左邊， 利用泰勒展開式展開，並忽略高次項，得： z y x zx yx xx ∂ ∂ + ∂ ∂ + ∂ ∂τ τ τ =ρ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ∂t du dt d 由全微分展開式將上式等號右邊括弧內展開，得： z y x zx yx xx ∂ ∂ + ∂ ∂ + ∂ ∂τ τ τ =ρ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ + •∇ ∂ ∂ u v t u dt d G (3.3) 對(2-1-2)式等號右邊括弧內兩項作尺度分析，將向量式vG、uG以張量分別表示成 、 。 因為在單位長度下的單位位移量很小，而且是從零開始變化的。所以 i v ui ) ( ) (L o L o = ∆ 、 ) ( ) (ε = o ∆ε o 且L>>ε ，得： = ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛ ∇ •∂ ∂ u v t u o G ( )>>1 ∆ ∆ = ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛ ∂ ∂ ∂ ∂ ε L o x u dt du t u o i i 因此忽略(3.3)式等號右邊括弧內第二項，得： z y x zx yx xx ∂ ∂ + ∂ ∂ + ∂ ∂τ τ τ =ρ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ∂ ∂ t u dt d (3.4) 由全微分展開式再將(3.4)式等號右邊展開，得： z y x zx yx xx ∂ ∂ + ∂ ∂ + ∂ ∂τ τ τ =ρ ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ∂ ∂ ∇ • + ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ∂ ∂ ∂ ∂ t u v t u t G (3.5) 對(3.5)式等號右邊括弧內兩項作尺度分析，得： = ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛ ∂ ∂ ∇ • ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ∂ ∂ ∂ ∂ t u v t u t G 0 ( )>>1 ∆ ∆ = ⎟⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜⎜ ⎜ ⎜ ⎜ ⎝ ⎛ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ∂ ∂ ∂ ∂ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ∂ ∂ ∂ ∂ ε L o t u x dt du t u t o i i 因此忽略(3.5)式等號右邊括弧內第二項，得： z y x zx yx xx ∂ ∂ + ∂ ∂ + ∂ ∂τ τ τ =ρ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ∂ ∂ ∂ ∂ t u t (3.6) (3.6)式為介質在 x 方向受力的控制方程式，同理可推導出 y 及 z 方向的控制方程式，並以張量式表示： j ij x ∂ ∂τ =ρ 2i 2 t u ∂ ∂ i,j=x,y,z (3.7)
令介質為等向的情況下，以張量式表示應力與應變的關係 (Timoshenko, 1951)為： j i k k ij ij x u x u ∂ ∂ + ∂ ∂ =λδ µ τ 2 (3.8) 其中 z w y v x u x u k k ∂ ∂ + ∂ ∂ + ∂ ∂ = ∂ ∂ ，u、 、v w分別為 x、y、z 方向的位移量，λ 及µ 為 Lame＇s constant。將(3.8)式代入(3.7)式，並以向量式表示，得： u u uG=(λ+µ)∇(∇•G)+µ∇2G ρ ⇒ ρuG=µ∇2uG+(λ+µ)∇(∇•uG) (3.9) 其中 = uG uiG+vGj +wkG，u為位移向量，u、v、 分別為 x、y、z 方向的位移量，(3.9)式為介質 受力的控制方程式之向量式。 G w (3.9)式是介質受力但不受阻尼力影響的控制方程式。阻尼力(damping force)是物體 在運動時，受到周圍其他不同介質的影響，使得物體運動隨時間趨緩。如：將鐵球吊在空 氣中的彈簧，當我們將鐵球向下拉一段距離後放開，鐵球開始震動，但受到空氣阻力的影 響，鐵球震動的距離會越來越小，最後會靜止，此空氣阻力即阻尼力。 此研究一開始是將介質視為單一密度、均勻的彈性體，在彈性體內是均質且等向的。 彈性體中沒有其他介質，理論上是不考慮阻尼力的影響。但在本研究中，實驗的砂石是不 飽和且孔隙中大多是空氣，因此砂石顆粒之間會有相當大摩擦，所以阻尼會很大。 考慮在(3.9)式中加入阻尼力，一般阻尼力以速度表示，阻尼力是正比於一阻尼係數乘 上物體運動的速度，在本研究中令此阻尼係數為常數。直接考慮介質中單位體積物體的阻 尼力，即： f dG t u J ∂ ∂ = G ρ =ρJ uG 其中 f dG 為阻尼力的向量式 uG=位移向量，其中u、 、v w分別為 x、y、z 方向的位移分量 J =阻尼係數 所以(3.9)式加入阻尼力，得： ) ( ) ( 2 u u d uG+ Gf =µ∇ G+ λ+µ ∇ ∇•G ρ ⇒ρuG+ρJuG=µ∇2uG+(λ+µ)∇(∇•uG) (3.10)
(二) 三維點源有限域解 本研究的實驗，砂石選用通過 4 號篩及 5 號篩的砂石，放入砂箱中並填滿砂箱但不飽 和，將地聲探測器放入砂石中，此地聲探測器所量測的資料可以轉換成速度，在實驗設備 時會描述其如何轉換。在砂石表面固定一點，鉛球離砂石表面 1 公尺的垂直高度放下，移 動地聲探測器到距離鉛球掉落處不同的水平位置，以探測距離點源不同位置的速度變化量。 配合實驗條件，建構理論的問題。
圖 3.2 三維點源受力
(1) 問題描述 考慮三維且 x、y、z 方向皆為有限邊界，在 z=0 的平面上有一垂直平面的作用力作 用在點(a,b,0)；在 x=0、H，y=0、L，z=D 處有堅硬的邊界，位移量為零且應力連續。座標 系統如圖：(H,L,D)
(a,b,0)
zz τ(H,0,0)
(0,0,0)
(0,L,0)
z
x
y
根據 Helmholtz theorem (Arfken, 1970; Phillips, 1956)，將位移向量uG表示成：
ψ φ G G =∇ +∇× u ，∇•ψG =0 (3.11) 其中 φ為 scale potential ψG為 vector potential，將其向量式表示成：ψG =ψ1iG+ψ2Gj+ψ3kG 0 = • ∇ ψG 為 gage condition φ、ψ1、ψ2、ψ3這四個變數，其中三個是獨立的，這三個獨立變數分別是φ及ψ1、ψ2、ψ3 中的其中兩個。由 gage condition 知，若已知其中兩個ψ 就可以由 gage condition 求得
令一個ψ 。 將(3.11)式代入(3.10)式，得： 0 )] ( [ )] ( 2 [ 2 2 2 2 2 2 = ∇ − ∂ ∂ + ∂ ∂ × ∇ + ∇ + − ∂ ∂ + ∂ ∂ ∇ ψ ρ µ ψ ψ φ ρ µ λ φ φ G G G t J t t J t (3.12) 由 Helmholtz theorem，(3.12)式可分解(Anchbach, 1973)成： ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∂ ∂ + ∂ ∂ = ∇ t J t cp φ φ φ 2 22 2 1 (3.13) ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∂ ∂ + ∂ ∂ = ∇ t J t cT ψ ψ ψG 2 22G G 2 1 (3.14) 其中Cp ρ µ λ+2 = ，Cp為縱波波速 ρ µ = T C ，CT為橫波波速 將(3.14)式分解成純量式： ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∂ ∂ + ∂ ∂ = ∇ t J t cT 1 2 1 2 2 1 2ψ 1 ψ ψ (3.15) ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∂ ∂ + ∂ ∂ = ∇ t J t cT 2 2 2 2 2 2 2ψ 1 ψ ψ (3.16) ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∂ ∂ + ∂ ∂ = ∇ t J t cT 3 2 3 2 2 3 2ψ 1 ψ ψ (3.17) (3.13)為縱波波 傳方程式，(3.14)、(3.15)、(3.16)、(3.17)為橫波波傳方程式。 根據基本假設，初始條件與邊界條件如下。 初始條件： t=0 x、y、z 三方向位移量為 0： 0 2 3 = ∂ ∂ − ∂ ∂ + ∂ ∂ = z y x u φ ψ ψ (3.18) 0 1 3 = ∂ ∂ + ∂ ∂ − ∂ ∂ = z x y v φ ψ ψ (3.19)
0 1 2 = ∂ ∂ − ∂ ∂ + ∂ ∂ = y x z w φ ψ ψ (3.20) x、y、z 三方向速度為 0： 0 2 3 = ∂ ∂ − ∂ ∂ + ∂ ∂ = z y x u φ ψ ψ (3.21) 0 1 3 = ∂ ∂ + ∂ ∂ − ∂ ∂ = z x y v φ ψ ψ (3.22) 0 1 2 = ∂ ∂ − ∂ ∂ + ∂ ∂ = y x z w φ ψ ψ (3.23) 邊界條件： 在 =x 0, H 處 堅硬的邊界，位移量為 0： 0 2 3 = ∂ ∂ − ∂ ∂ + ∂ ∂ = z y x u φ ψ ψ (3.24) 應力連續： 0 2 2 2 3 2 3 3 3 3 2 3 2 3 3 3 = ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ ∂ ∂ ∂ − ∂ ∂ ∂ + ∂ ∂ + ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ ∂ ∂ ∂ + ∂ ∂ ∂ + ∂ ∂ = ∂ ∂ z x y x x z x y x x x xx λ φ φ φ µ φ ψ ψ τ (3.25) 0 2 2 1 3 3 3 3 2 3 2 3 3 2 3 = ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ ∂ ∂ ∂ + ∂ ∂ − ∂ ∂ ∂ ∂ − ∂ ∂ ∂ + ∂ ∂ ∂ = ∂ ∂ z x x z y x y x y x x xy µ φ ψ ψ ψ ψ τ (3.26) 0 2 2 1 3 3 2 3 2 2 3 3 3 2 3 = ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ ∂ ∂ ∂ + ∂ ∂ − ∂ ∂ ∂ − ∂ ∂ ∂ ∂ + ∂ ∂ ∂ = ∂ ∂ y x x z x z y x z x x xz µ φ ψ ψ ψ ψ τ (3.27) 在 =y 0 , L處 堅硬的邊界，位移量為 0： 0 1 3 = ∂ ∂ + ∂ ∂ − ∂ ∂ = z x y v φ ψ ψ (3.28) 應力連續： 0 2 2 1 3 2 3 3 3 3 2 3 3 3 2 3 = ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ ∂ ∂ ∂ + ∂ ∂ ∂ − ∂ ∂ + ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ ∂ ∂ ∂ + ∂ ∂ + ∂ ∂ ∂ = ∂ ∂ z y y x y z y y y x y yy λ φ φ φ µ φ ψ ψ τ (3.29) 0 2 1 3 2 3 3 2 2 3 3 3 3 2 3 = ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ ∂ ∂ ∂ ∂ + ∂ ∂ ∂ − ∂ ∂ ∂ − ∂ ∂ + ∂ ∂ ∂ = ∂ ∂ z y x y x z y y y x y yx µ φ ψ ψ ψ ψ τ (3.30)
0 2 31 3 2 2 3 2 1 3 3 3 2 3 = ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ ∂ ∂ − ∂ ∂ ∂ + ∂ ∂ ∂ + ∂ ∂ ∂ ∂ − ∂ ∂ ∂ = ∂ ∂ y y x z y z y x z y y yz µ φ ψ ψ ψ ψ τ (3.31) 在 z=0 處 單點作用力：
鉛球從高處撞擊到(a,b,o)處的作用力是脈衝函數，用 cosine function 來代替脈衝函 數。在 x、y 方向變化為 cosine function。在隨時間的變化上，將鉛球綁上繩子，當鉛球 從高處撞擊到(a,b,o)處時馬上拉起，使鉛球不停留在(a,b,o)處，因此作用力隨 t 變化為 f(t)： ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ ∂ ∂ ∂ − ∂ ∂ ∂ + ∂ ∂ + ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ ∂ ∂ + ∂ ∂ + ∂ ∂ = z y z x z z y x zz 1 2 2 2 2 2 2 2 2 2 2 2 2µ φ ψ ψ φ φ φ λ τ ) ( ) ( ) (x f y f t f σ − = (3.32) 在此的 、 、 分別代表作用力隨 x、y、t 的變化的 function，隨著不同的實 驗方法 function 會跟著變，而本研究的實驗所對應的 function： ) (x f f( y) f(t) ) 2 cos( ) ( π c a x x f = − a−c< x<a+c =0 0<x<a−c,a+c<x<H ) 2 cos( ) ( π c b y y f ′ − = b−c′< y<b+c′ =0 0< y<b−c′,b+c′< y<L ) 2 sin( ) ( π c t t f ′′ = 0<t <c′′ =0 t >c′′ σ 是常數，c<<H，c′<<L。 自由邊界(free surface)： 0 2 1 2 2 2 2 2 2 2 3 2 2 = ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ ∂ ∂ ∂ − ∂ ∂ + ∂ ∂ − ∂ ∂ ∂ + ∂ ∂ ∂ = y x x z z y z x zx ψ ψ ψ ψ φ µ τ (3.33) 0 2 21 2 2 2 2 1 2 3 2 2 = ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ ∂ ∂ − ∂ ∂ ∂ + ∂ ∂ + ∂ ∂ ∂ − ∂ ∂ ∂ = y y x z z x z y zy ψ ψ ψ ψ φ µ τ (3.34) 在 z=D 處
堅硬的邊界，位移量為 0： 0 1 2 = ∂ ∂ − ∂ ∂ + ∂ ∂ = y x z w φ ψ ψ (3.35) 應力連續： 0 2 12 3 2 2 3 3 3 3 3 2 3 2 3 = ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ ∂ ∂ ∂ − ∂ ∂ ∂ + ∂ ∂ + ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ ∂ ∂ + ∂ ∂ ∂ + ∂ ∂ ∂ = ∂ ∂ z y z x z z z y z x z zz λ φ φ φ µ φ ψ ψ τ (3.36) 0 2 1 3 2 2 3 3 2 3 2 3 3 2 3 = ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ ∂ ∂ ∂ ∂ − ∂ ∂ ∂ + ∂ ∂ − ∂ ∂ ∂ + ∂ ∂ ∂ = ∂ ∂ z y x z x z z y z x z zx µ φ ψ ψ ψ ψ τ (3.37) 0 2 2 1 3 2 3 3 1 3 2 3 3 2 3 = ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ ∂ ∂ ∂ − ∂ ∂ ∂ ∂ + ∂ ∂ + ∂ ∂ ∂ − ∂ ∂ ∂ = ∂ ∂ z y z y x z z x z y z zy µ φ ψ ψ ψ ψ τ (3.38) (2) 解φ、ψ 、1 ψ 、2 ψ 3 根據 Helmholtz theorem，位移量u、v、w由φ、ψ 、1 ψ 、2 ψ 表示。如果解出3 φ、ψ 、1 2 ψ 、ψ ，就可以解出3 u、 、v w。由(1)中所得得控制方程式： ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∂ ∂ + ∂ ∂ = ∇ t J t cp φ φ φ 2 22 2 1 (3.39) ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∂ ∂ + ∂ ∂ = ∇ t J t cT 1 2 1 2 2 1 2ψ 1 ψ ψ (3.40) ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∂ ∂ + ∂ ∂ = ∇ t J t cT 2 2 2 2 2 2 2ψ 1 ψ ψ (3.41) ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∂ ∂ + ∂ ∂ = ∇ t J t cT 3 2 3 2 2 3 2ψ 1 ψ ψ (3.42) 現在要解上四式，利用分離變數，將φ、ψ 、1 ψ 、2 ψ 分別表示為： 3 ) ( ) ( ) ( ) (xY y Z z T t X = φ (3.43) ) ( ) ( ) ( ) ( 1 = X x Y y Z z T t ψ (3.44) ) ( ˆ ) ( ˆ ) ( ˆ ) ( ˆ 2 = X x Y y Z z T t ψ (3.45) ) ( ~ ) ( ~ ) ( ~ ) ( ~ 3 = X x Y y Z z T t ψ (3.46) 將上四式代入(3.13)、(3.15)、(3.16)、(3.17)式中，得： λ = ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ + = ′′ + ′′ + ′′ ) ( ) ( ) ( ) ( 1 ) ( ) ( ) ( ) ( ) ( ) ( 2 t T t T J t T t T c z Z z Z y Y y Y x X x X p (3.47) λ = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + = ′′ + ′′ + ′′ ) ( ) ( ) ( ) ( 1 ) ( ) ( ) ( ) ( ) ( ) ( 2 t T t T J t T t T c z Z z Z y Y y Y x X x X T (3.48)
λˆ ) ( ˆ ) ( ˆ ) ( ˆ ) ( ˆ 1 ) ( ˆ ) ( ˆ ) ( ˆ ) ( ˆ ) ( ˆ ) ( ˆ 2 ⎟= ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + = ′′ + ′′ + ′′ t T t T J t T t T c z Z z Z y Y y Y x X x X T (3.49) λ~ ) ( ~( ) ~ ) ( ~( ) ~ 1 ) ( ~( ) ~ ) ( ~( ) ~ ) ( ~( ) ~ 2 ⎟= ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + = ′′ + ′′ + ′′ t T t T J t T t T c z Z z Z y Y y Y x X x X T (3.50) 其中λ 、λ 、 、λˆ λ~與 x、y、z、t 無關 令 p x X x X ′′ = ) ( ) ( q y Y y Y′′ = ) ( ) ( r z Z z Z′′ = ) ( ) ( p x X x X = ′′ ) ( ) ( q y Y y Y = ′′ ) ( ) ( r z Z z Z = ′′ ) ( ) ( p x X x X ˆ ) ( ˆ ) ( ˆ = ′′ q y Y y Y ˆ ) ( ˆ ) ( ˆ = ′′ r z Z z Z ˆ ) ( ˆ ) ( ˆ = ′′ p x X x X ~ ) ( ~( ) ~ = ′′ q y Y y Y ~ ) ( ~( ) ~ = ′′ r z Z z Z ~ ) ( ~( ) ~ = ′′ 其中 p 、p、 、 、 q 、pˆ ~p q、 q 、 、ˆ q~ r 、r、 rˆ 、 r~ 皆為常數，共需要 24 個條件來決定。 I. 邊界條件的處理 已知φ、ψ 、1 ψ 、2 ψ 這四個變數，其中三個是獨立的，這三個獨立變數分別是3 φ及ψ 、1 2 ψ 、ψ 中的其中兩個。若令3 φ、ψ 、2 ψ 為獨立變數，將(3.39)、(3.130)、(3.131)、(3.132)3 式，代入邊界條件(3.24)式中， 在x=0處 0 2 3 = ∂ ∂ − ∂ ∂ + ∂ ∂ = z y x u φ ψ ψ (3.24) 得： 0 2 3 = ∂ ∂ − ∂ ∂ + ∂ ∂ = z y x u φ ψ ψ ⇒ X ′(0)Y(y)Z(z)T(t)+X~(0)Y~′(y)Z~(z)T~(t)-Xˆ(0)Yˆ(y)Zˆ′(z)Tˆ(t) = 0 因為φ、ψ 、2 ψ 為獨立變數，在同一點 y 或 z 時也是獨立變數，所以3 Y(y)Z(z)T(t)、 ) ( ~ ) ( ~ ) ( ~ t T z Z y Y ′ 、Yˆ(y)Zˆ′(z)Tˆ(t)彼此為獨立函數，得： ) 0 ( X ′ = X~(0) = Xˆ(0)= 0 同理 在 x=H 處
0 2 3 = ∂ ∂ − ∂ ∂ + ∂ ∂ = z y x u φ ψ ψ (3.24) ⇒ X ′(H)Y(y)Z(z)T(t)+X~(H)Y~′(y)Z~(z)T~(t)-Xˆ(H)Yˆ(y)Zˆ′(z)Tˆ(t) = 0 得： ) (H X ′ = X~(H) = Xˆ H( )= 0 在 y=0 處 0 1 3 = ∂ ∂ + ∂ ∂ − ∂ ∂ = z x y v φ ψ ψ (3.28) 由 gage condition： 0 3 2 1 = ∂ ∂ + ∂ ∂ + ∂ ∂ = • ∇ z y x ψ ψ ψ ψG ⇒ dx z y
∫
+∂∂ ∂ ∂ − = ( 2 3) 1 ψ ψ ψ 將上式代入(3.28)式，得： 0 2 3 2 2 2 3 = ∂ ∂ − ∂ ∂ ∂ − ∂ ∂ − ∂ ∂∫
∫
dx z dx z y x y ψ ψ ψ φ 將(3.39)、(3.40)、(3.41)、(3.42)式代入上式，並整理得： ) ( ) ( ) ( ) 0 ( X x Z z T t Y ′ −Y~(0)[X~′(x)Z~(z)T~(t)+∫
X~(x)Z~′′(z)T~(t)dx] ) ( ˆ ) ( ˆ ) ( ˆ ) 0 ( ˆ X x Z z T t Y′ ′ − = 0 由上式得： ) 0 ( Y ′ =Y~(0)=Y ′ˆ (0)= 0 同理由(3.28)式，得： 在 y=L 處 ) (L Y ′ =Y~(L) =Y ′ˆ L( )= 0 同理由(3.35)式，得： 在 z=D 處 0 1 2 = ∂ ∂ − ∂ ∂ + ∂ ∂ = y x z w φ ψ ψ 得：) (D
Z ′ =Zˆ D( ) =Z ′~(D)= 0
根據 Helmholtz theorem，∇•ψG =0為 gage condition，在此問題的區域內任何一點都成 立，因此 在 x=0處 0 3 2 1 = ∂ ∂ + ∂ ∂ + ∂ ∂ = • ∇ z y x ψ ψ ψ ψG ⇒ X ′(0)Y(y)Z(z)T(t)+Xˆ(0)Yˆ′(y)Zˆ(z)Tˆ(t)+X~(0)Y~(y)Z~′(z)T~(t) = 0 因為Xˆ(0)= X~(0)= 0，所以得 X ′(0)= 0 同理，在 x=H 處 0 3 2 1 = ∂ ∂ + ∂ ∂ + ∂ ∂ = • ∇ z y x ψ ψ ψ ψG 得X ′(H)= 0。在 y= 0,L 處 0 3 2 1 = ∂ ∂ + ∂ ∂ + ∂ ∂ = • ∇ z y x ψ ψ ψ ψG 得Y(0)= 0 ) (L Y = 0 在z = D處 0 3 2 1 = ∂ ∂ + ∂ ∂ + ∂ ∂ = • ∇ z y x ψ ψ ψ ψG 得 ) (D Z = 0 將以上條件整理，得： x=0 X ′(0)= X ′(0)= Xˆ(0)= X~(0)= 0 x H = X ′(H)= X ′(H)= Xˆ H( ) = X~(H)= 0 y 0 = Y ′(0)=Y(0)=Y ′ˆ (0)=Y~(0)= 0 y L = Y ′(L) =Y(L) =Y ′ˆ L( )=Y~(L) = 0 z D = Z ′(D) =Z(D)=Zˆ D( ) =Z ′~(D)= 0 從上面得知在 z=0 還少四個邊界條件。 在 z=0 處
0 2 1 2 2 2 2 2 2 2 3 2 2 = ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ ∂ ∂ ∂ − ∂ ∂ + ∂ ∂ − ∂ ∂ ∂ + ∂ ∂ ∂ = y x x z z y z x zx ψ ψ ψ ψ φ µ τ 由 gage condition： 0 3 2 1 = ∂ ∂ + ∂ ∂ + ∂ ∂ = • ∇ z y x ψ ψ ψ ψG ⇒ dx z y
∫
+∂∂ ∂ ∂ − = ( 2 3) 1 ψ ψ ψ 將上式代入(3.51)式，得： 0 2 3 2 2 2 2 2 2 2 2 2 2 3 2 2 = ∂ ∂ ∂ + ∂ ∂ + ∂ ∂ + ∂ ∂ − ∂ ∂ ∂ + ∂ ∂ ∂ z x y x z z y z x ψ ψ ψ ψ ψ φ 將(3.39)、(3.40)、(3.41)、(3.42)式代入上式，並整理得： ⇒2Z′(0)X′(x)Y(y)T(t)+Z~′(0)[X~(x)Y~′(y)T~(t)+ X~′(x)Y~(y)T~(t)] )] ( ˆ ) ( ˆ ) ( ˆ ) ( ˆ ) ( ˆ ) ( ˆ )[ 0 ( ˆ X x Y y T t X x Y y T t Z ′′ + ′′ + −Z ′′ˆ (0)Xˆ(0)Yˆ(y)Tˆ(t)= 0 (3.51) 討論(3.51)式等號左邊的第三、四項，將特徵方程式 r z Z z Z ˆ ) ( ˆ ) ( ˆ = ′′ 求解得： ) ˆ sinh( ) ˆ cosh( ) ( ˆ 2 1 sz s sz s z Z = + 其中rˆ= sˆ2 >0 z s s z Zˆ( )= 3 + 4 其中rˆ=0 ) ˆ sin( ) ˆ cos( ) ( ˆ 6 5 sz s sz s z Z = + 其中rˆ= s−ˆ2 <0 因此有三種解，不論解的形式是那一種， 與 的函數都一樣，兩者的差別是係數 不同而已，所以(3.51)式等號左邊的第三、四項可以合併，因此由(3.51)式可得： ) ( ˆ z Z Z ′′ˆ z( ) ) 0 ( Z ′ =Zˆ(0) =Z ′~(0) = 0 在 z=0 處，由 gage condition： 0 3 2 1 = ∂ ∂ + ∂ ∂ + ∂ ∂ = • ∇ z y x ψ ψ ψ ψG 可得 X(x)′Y(y)Z(0)T(t)+Xˆ(x)Yˆ′(y)Zˆ(0)Tˆ(t)+X~(x)Y~(y)Z~′(0)T~(t)= 0 因為Zˆ(0)=Z ′~(0) = 0，所以得 Z(0) = 0 在此 z=0 的邊界條件已全部找出，為 ) 0 ( Z ′ =Z(0) = Zˆ(0)= Z ′~(0) = 0II. 解出特徵值及特徵函數 由 a.所求得的邊界條件解特徵方程式。 特徵方程式： p x X x X ′′ = ) ( ) ( 邊界條件：x=0 X ′(0)= 0 x H = X ′(H)= 0 得： 特徵值： p = 2 ) ( H mπ − ⇒特徵函數： cos( x) H m c X = m π m=0,1,2... 特徵方程式： p x X x X = ′′ ) ( ) ( 邊界條件：x=0 X ′(0)= 0 x H = X ′(H)= 0 得特徵值：p = 1 2 ) ( H mπ − ⇒特徵函數： cos( 1 ) 1 x H m c X m π = m1 =0,1,2... 特徵方程式： p x X x X ˆ ) ( ˆ ) ( ˆ = ′′ 邊界條件：x=0 Xˆ(0)= 0 x H = Xˆ H( ) = 0 得特徵值： pˆ = ( 2 )2 H m π − ⇒特徵函數： ˆ sin( 2 ) 2 x H m c X = m π m2 =1,2... 特徵方程式： p x X x X ~ ) ( ~( ) ~ = ′′ 邊界條件：x=0 X~(0)= 0 x H = X~(H) = 0 得特徵值： p~ = 3 2 ) ( H m π − ⇒特徵函數： ~ sin( 3 ) 3 H x m c X = m π m3 =1,2... 特徵方程式： q y Y y Y = ′′ ) ( ) ( 邊界條件：y 0 = Y ′(0)= 0 y L = Y ′(L) = 0 得特徵值： q= 2 ) ( L nπ − ⇒特徵函數： cos( y) L n c Y = n π n=0,1,2...
特徵方程式： q y Y y Y ′′ = ) ( ) ( 邊界條件：y 0 = Y(0)= 0 y L = Y(L) = 0 得特徵值：q = 1 2 ) ( L nπ − ⇒特徵函數： sin( 1 ) 1 y L n c Y = n π n1 =1,2... 特徵方程式： q y Y y Y ˆ ) ( ˆ ) ( ˆ = ′′ 邊界條件：y 0 = Y ′ˆ (0)= 0 y L = Y ′ˆ L( ) = 0 得特徵值：qˆ = 2 2 ) ( L n π − ⇒特徵函數： ˆ cos( 2 ) 2 y L n c Y = n π n2 =0,1,2... 特徵方程式： q y Y y Y ~ ) ( ~( ) ~ = ′′ 邊界條件：y 0 = Y~(0)= 0 y L = Y~(L) = 0 得特徵值： q~ = 3 2 ) ( L nπ − ⇒特徵函數：~ sin( 3 ) 3 y L n c Y n π = n3 =1,2... 特徵方程式： r z Z z Z′′ = ) ( ) ( 邊界條件：z 0 = Z ′(0) = 0 z D = Z ′(D) = 0 得特徵值：r = ( )2 D eπ − ⇒特徵函數： cos( z) D e c Z e π = e=0,1,2... 特徵方程式： r z Z z Z ′′ = ) ( ) ( 邊界條件：z 0 = Z(0) = 0 z D = Z(D)= 0 得特徵值：r = ( 1 )2 D eπ − ⇒特徵函數： sin( 1 ) 1 D z e c Z = e π e1 =1,2... 特徵方程式： r z Z z Z ˆ ) ( ˆ ) ( ˆ = ′′ 邊界條件：z 0 = Zˆ(0)= 0
z D = Zˆ D( ) = 0 得特徵值：rˆ = 2 2 ) ( D e π − ⇒特徵函數： ˆ sin( 2 ) 2 D z e c Z = e π e2 =1,2... 特徵方程式： r z Z z Z ~ ) ( ~( ) ~ = ′′ 邊界條件：z 0 = Z ′~(0) = 0 z D = Z ′~(D) = 0 得特徵值： r~ = 3 2 ) ( D eπ − ⇒特徵函數：~ cos( 3 ) 3 z D e c Z = e π e3 =0,1,2... III. 求φ 、ψ1、ψ2、ψ 的全解 3 在此，所有的特徵值及特徵函數都求出，將所有特徵值代入下四式： λ = ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ + = ′′ + ′′ + ′′ ) ( ) ( ) ( ) ( 1 ) ( ) ( ) ( ) ( ) ( ) ( 2 t T t T J t T t T c z Z z Z y Y y Y x X x X p (3.52) λ = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + = ′′ + ′′ + ′′ ) ( ) ( ) ( ) ( 1 ) ( ) ( ) ( ) ( ) ( ) ( 2 t T t T J t T t T c z Z z Z y Y y Y x X x X T (3.53) λˆ ) ( ˆ ) ( ˆ ) ( ˆ ) ( ˆ 1 ) ( ˆ ) ( ˆ ) ( ˆ ) ( ˆ ) ( ˆ ) ( ˆ 2 ⎟= ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + = ′′ + ′′ + ′′ t T t T J t T t T c z Z z Z y Y y Y x X x X T (3.54) λ~ ) ( ~( ) ~ ) ( ~( ) ~ 1 ) ( ~( ) ~ ) ( ~( ) ~ ) ( ~( ) ~ 2 ⎟= ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + = ′′ + ′′ + ′′ t T t T J t T t T c z Z z Z y Y y Y x X x X T (3.55) 先解(3.52)式的T(t)，特徵值及特徵函數代入(3.52)式，得： 0 ) ( ) ( ) (t +JT t +c 2k2 T t = T p mne (3.56) 其中 mne k ( )2 ( )2 ( )2 D e L n H mπ + π + π = 令 t代入(3.56)式，得 e t T( )= α 0 2 2 2 + + = mne p k c Jα α 所以： 2 4 2 2 2 mne p k c J J ± − − = α 如果 2 −4 2 2 <0，得： mne p k c J
) (t T = ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ − + − − ) 4 sin( ) 4 cos( 2 2 2 2 2 2 2 1 2 C c k J t C c k J t e p mne p mne t J (3.57) 如果J2 −4cp2k2mne =0，得： ) (t T =e 2t
[
C1 C2t) J ′ + ′ −]
(3.58) 如果J2 −4cp2k2mne >0，得： ) (t T = ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ ′′ + ′′ − − − −Jt J cpkmnet J cpkmnet e C e C e 2 2 2 2 2 2 4 2 4 1 2 (3.59) 1 C 、C2、C′1、C′2、C ′′1 、C ′′2皆為常數。 同理T(t)的解也有三種情形： 如果 4 111 0 2 2 2 − < e n m T k c J ，得： ) (t T = ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ − + − − ) 4 sin( ) 4 cos( 2 2 2 4 2 2 2 3 2 1 1 1 1 1 1 t J k c C t J k c C e T mne T mne t J (3.60) 如果 4 111 0 2 2 2 − = e n m T k c J ，得： ) (t T =[
3 4 ) 2 C C t e t J ′ + ′ −]
(3.61) 如果 4 111 0 2 2 2 − > e n m T k c J ，得： ) (t T = ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ ′′ + ′′ − − − −Jt J cT k mne t J cT k mne t e C e C e 111 2 2 2 1 1 1 2 2 2 4 4 4 3 2 (3.62) 其中 1 1 1ne m k ( 1 )2 ( 1 )2 ( 1 )2 D e L n H mπ π π + + = 3 C 、C4、C′3、C′4、C ′′ 、3 C ′′4 皆為常數。 同理Tˆ t( )的解也有三種情形： 如果 4 ˆ 2 22 0，得： 2 2 2 − < e n m T k c J ) ( ˆ t T = ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ − + − − ) 4 ˆ sin( ) 4 ˆ cos( 2 2 2 6 2 2 2 5 2 2 2 2 2 2 2 t J k c C t J k c C e T mne T mne t J (3.63) 如果 4 ˆ 2 22 0，得： 2 2 2− = e n m T k c J) ( ˆ t T =e 2t
[
C5 C6t) J ′ + ′ −]
(3.64) 如果 4 ˆ 2 22 0，得： 2 2 2 − > e n m T k c J ) ( ˆ t T = ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ ′′ + ′′ − − − −Jt J cT k mne t J cT k mne t e C e C e 222 2 2 2 2 2 2 2 2 2 ˆ 4 6 ˆ 4 5 2 (3.65) 其中 2 2 2 ˆ e n m k 2 2 2 2 2 2 ) ( ) ( ) ( D e L n H m π π π + + = 5 C 、C6、C′5、C′6、C ′′5 、C ′′6 皆為常數。 同理T~(t)的解也有三種情形： 如果 4 ~ 333 0 2 2 2 − < e n m T k c J ，得： ) ( ~ t T = ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ − + − − ) 4 ~ sin( ) 4 ~ cos( 2 2 2 8 2 2 2 7 2 3 3 3 3 3 3 t J k c C t J k c C e T mne T mne t J (3.66) 如果 4 ~ 333 0，得： 2 2 2 − = e n m T k c J ) ( ~ t T =[
7 8 ) 2 C C t e t J ′ + ′ −]
(3.67) 如果 4 ~ 333 0 2 2 2 − > e n m T k c J ，得： ) ( ~ t T = ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ ′′ + ′′ − − − −Jt J cT k mne t J cT k mne t e C e C e 333 2 2 2 3 3 3 2 2 2 ~ 4 8 ~ 4 7 2 (3.68) 其中 3 3 3 ~ e n m k 3 2 3 2 3 2 ) ( ) ( ) ( D e L n H mπ π π + + = 7 C 、C8、C′7、C′8、C ′′ 、7 C ′′8 皆為常數。 由於cp、cT、 都是未定的係數，先考慮J 4 0、 2 2 2 − < mne p k c J 0 4 111 2 2 2 − < e n m T k c J 、 4 ~ 3 33 0 2 2 2 − < e n m T k c J 、 情形下的解，分別將所求得 的(3.59)、(3.62)、(3.65)、(3.68)式及特徵函數及代入(3.39)、(3.130)、(3.131)、 (3.132)，利用重疊原理，得 0 ˆ 4 2 22 2 2 2 − < e n m T k c J φ、ψ1、ψ2、ψ3的全解： φ = ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ − + − ∑ ∑ ∑ ∞ − = ∞ = ∞ = cos( 4 ) sin( 4 ) 2 2 2 2 2 2 2 0 0 0 t J k c B t J k c Ae mne p mne mne p mne t
J
e n m
) cos( x H mπ ) cos( y L nπ ) cos( z D eπ (3.69) 1 ψ = ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ − + − ∑ ∑ ∑ ∞ − = ∞ = ∞ = cos( 4 ) sin( 4 ) 2 2 2 2 2 2 2 1 1 0 t J k c B t J k c A
e mne T mne mne T mne t J e n m ) cos( x H mπ ) sin( y L nπ ) sin( z D eπ (3.70) 2 ψ = ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ − + − ∑ ∑ ∑ ∞ − = ∞ = ∞ = sin( 4 ) ˆ ) 4 cos( ˆ 2 2 2 2 2 2 2 1 0 1 t J k c B t J k c A
e mne T mne mne T mne t J e n m ) sin( x H mπ ) cos( y L nπ ) sin( z D eπ (3.71) 3 ψ = ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ − + − ∑ ∑ ∑ ∞ − = ∞ = ∞ = sin( 4 ) ~ ) 4 cos( ~ 2 2 2 2 2 2 2 0 1 1 t J k c B t J k c A
e mne T mne mne T mne t J e n m ) sin( x H mπ ) sin( y L nπ ) cos( z D eπ (3.72) 將φ、ψ1、ψ2、ψ3的全解(3.69)、(3.70)、(3.71)、(3.72)式，代入初始條件， 得： 在 t = 0 時 x 方向位移量u=0： 0 2 3 = ∂ ∂ − ∂ ∂ + ∂ ∂ = z y x u φ ψ ψ (3.18) ⇒ ( ) 0 0 0n e mne mne m B A + ∑ ∑ ∑ − ∞ = ∞ = ∞ = H mπ ) sin( x H mπ ) cos( y L nπ ) cos( z D eπ + 3 3 3 3 3 3 ~ 0 1 1n e mne m A ∞ = ∞ = ∞ = ∑ ∑ ∑ L n3π ) sin( 3 x H mπ ) cos(3 y L nπ ) cos( 3 z D eπ + 2 2 2 2 2 2 ˆ 1 0 1n e mne m A ∞ = ∞ = ∞ = ∑ ∑ ∑ D e2π ) sin( 2 x H mπ ) cos( 2 y L nπ ) cos(2 z D eπ 0 = 將m=n=e=e3 =n2 =0提出，並使 ∞ ⋅⋅ ⋅ = = =m2 m3 1,2 , m ∞ ⋅⋅ ⋅ = = =n2 n3 1,2 , n ∞ ⋅⋅ ⋅ = = =e2 e3 1,2 , e
得： mne e n m A − ∑ ∑ ∑ ∞ = ∞ = ∞ =1 1 1( H mπ mne A~ + L nπ mne Aˆ − ) D eπ ) sin( x H mπ ) cos( y L nπ ) cos( z D eπ + mn n m A − ∑ ∑ ∞ = ∞ =1 1( H mπ 0 ~ mn A + L nπ )sin( x) H mπ ) cos( y L nπ + m e e m A 0 1 1 (− ∑ ∑ ∞ = ∞ = H mπ e m Aˆ 0 − ) D eπ ) sin( x H mπ ) cos( z D eπ + 00 1 ( m m A − ∑∞ = H mπ ) sin( x) H mπ 0 = (3.73) 由上式得： mne A H m mne A~ − L n mne Aˆ + D e 0 = (3.74) 0 mn A H m 0 ~ mn A − L n 0 = (3.75) e m A 0 H m e m Aˆ 0 + D e 0 = (3.76) 00 m A =0 (3.77) 在t=0時 y 方向位移量v=0： 0 1 3 = ∂ ∂ + ∂ ∂ − ∂ ∂ = z x y v φ ψ ψ (3.19) ⇒ mne e n m A ∞ = ∞ = ∞ = ∑ ∑ ∑ − 0 0 0 L nπ ) cos( x H mπ ) cos( y L nπ ) cos( z D eπ 3 3 3 3 3 3 ~ 0 1 1n e mne m A ∞ = ∞ = ∞ = ∑ ∑ ∑ − H m3π ) cos( 3 x H mπ ) sin( 3 y L nπ ) cos( 3 z D eπ + 1 1 1 1 1 1 0 1 1 e n m e n m A ∞ = ∞ = ∞ = ∑ ∑ ∑ D e1π ) cos( 1 x H mπ ) sin( 1 y L nπ ) cos(1 z D eπ 0 = 將m=n=e=e3 =m1 =0提出，並使 ∞ ⋅⋅ ⋅ = = =m1 m3 1,2 , m ∞ ⋅⋅ ⋅ = = =n1 n3 1,2 , n ∞ ⋅⋅ ⋅ = = =e1 e3 1,2 , e 得： mne e n m A − ∑ ∑ ∑ ∞ = ∞ = ∞ =1 1 1( L nπ mne A + D eπ mne A~ − ) H mπ ) cos( x H mπ ) sin( y L nπ ) cos( z D eπ
+ 0 1 1 ( mn n m A − ∑ ∑ ∞ = ∞ = L nπ 0 ~ mn A − H mπ ) cos( x) H mπ ) sin( y L nπ + ne e n A0 1 1 (− ∑ ∑ ∞ = ∞ = L nπ ne A0 + ) D eπ ) sin( y L nπ ) cos( z D eπ + 0 0 1 ( n n A − ∑∞ = H mπ ) sin( y) L nπ 0 = (3.78) 由上式得： mne A L n mne A − D e mne A~ + H m 0 = (3.79) 0 mn A L n 0 ~ mn A + H m 0 = (3.80) ne A0 L n ne A0 − D e 0 = (3.81) 0 0n A =0 (3.82) 在t=0時 z 方向位移量w=0： 0 1 2 = ∂ ∂ − ∂ ∂ + ∂ ∂ = y x z w φ ψ ψ (3.20) ⇒ mne e n m A ∞ = ∞ = ∞ = ∑ ∑ ∑ − 0 0 0 D eπ ) cos( x H mπ ) cos( y L nπ ) sin( z D eπ 2 2 2 2 2 2 ˆ 1 0 1n e mne m A ∞ = ∞ = ∞ = ∑ ∑ ∑ + H m2π ) cos( 2 x H mπ ) cos( 2 y L n π ) sin(2 z D eπ 1 1 1 1 1 1 0 1 1 e n m e n m A ∞ = ∞ = ∞ = ∑ ∑ ∑ − L n1π ) cos( 1 x H mπ ) cos( 1 y L nπ ) sin(1 z D eπ 0 = 將m=n=e=n2 =m1 =0提出，並使 ∞ ⋅⋅ ⋅ = = =m1 m2 1,2 , m ∞ ⋅⋅ ⋅ = = =n1 n2 1,2 , n ∞ ⋅⋅ ⋅ = = =e1 e2 1,2 , e 得： mne e n m A − ∑ ∑ ∑ ∞ = ∞ = ∞ =1 1 1( D eπ mne Aˆ + H mπ mne A − ) L nπ ) cos( x H mπ ) cos( y L nπ ) sin( z D eπ + m e e m A 0 1 1 (− ∑ ∑ ∞ = ∞ = D eπ e m Aˆ 0 + H mπ ) cos( x) H mπ ) sin( z D eπ + ne e n A0 1 1 (− ∑ ∑ ∞ = ∞ = D eπ ne A0 − ) L nπ ) cos( y L nπ ) sin( z D eπ
