利用有限元素波譜法建立河海三維水理水質數值模式(3/3)

(1)

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

利用有限元素波譜法建立河海三維水理水質數值模式(3/3)

計畫類別：個別型計畫計畫編號： NSC91-2211-E-002-060- 執行期間： 91 年 08 月 01 日至 92 年 07 月 31 日執行單位：國立臺灣大學土木工程學系暨研究所計畫主持人：劉格非報告類型：完整報告處理方式：本計畫可公開查詢

中華民國 92 年 12 月 4 日

(2)

利用有限元素波譜法建立河海三維水理水質數值模式3/3 Finite element spectral method applied in 3-D nearshore numerical model (III)

計畫編號：NSC91-2211-E-002-090 執行單位：國立台灣大學土木研究所研究期間：91 年 8 月 1 日至 92 年 7 月 31 日計畫主持人：劉格非國立台灣大學土木工程研究所中文摘要（關鍵詞字：波譜元素法，數值模擬，教學網站）本研究的重點在於發展波譜元素法，此方法同時包含快速計算與可求得空間中任一點的優點。本方法的基礎建立於波譜法與有限元素的原理，其為權重餘數法（MWR）的一種。在本研究中所選取的基函數為 Chebyshev 多項式，其優點在於可處理非週期性的邊界問題。並在研究中利用有限元素的概念將問題的定義域（ domain ）分割成數個次定義域（sub-domains），以加快計算的收斂速度及增加計算時的精度。再利用現場採樣來驗證程式，同時配合參數檢定程式來檢定現場水質參數如曝氣係數、延散係數。最後利用上述成果建立一個教學網站，將研究與程式結合成教學工具。

英文摘要(KEY WORDS：Finite element

spectral method, Numerical simulation, Interactive teaching Web)

The purpose of this research is to develop the Spectral Element Method. The trial functions is chosen to be Chebyshev function which can handle non-period boundary problem. We use the concept of the Finite Element Method to divide computational domain to several sub-domains to accelerate the convergence rate and the accuracy .

Field experiments are carried out in order to test the validity of the numerical program. Field measurements are used to validate the numerical program. Combined with calibration program, many in situ parameters, such as recreation coefficient and dispersion coefficient, are also obtained.

All the above then is combined to establish an interactive teaching web.

一、數值模擬與參數檢定 研究水理之數值模式比較常見的有有限差分法（王順賢 1999）、有限解析法（蘇怡中 1998）、有限元素法（簡紹欽 1998）、邊界元素法（葛家豪）、有限體積法和波譜法（郭佳虹 1997）。以上之數值方法中的有限差分、有限解析法與有限體積法皆採局部近似的觀念 (local approximation)，即處理空間微分時將其離散化。而邊界元素法、有限元素法和波譜法是採取整體近似。國內在波譜法方面的研究，有顏瑞和（ 1996 ），汪吉昌（1997），陳效吾（1996）等。皆為採用配置點法（collocation method），其理論基礎皆為選取 Chebyshev 權重函數所產生之點來作分析。本研究是利用 Chebyshev τ method 以一展開級數來求解整個空間的物理特性。免去取點的繁雜與 Fourier 轉換的麻煩，並且當展開級數之係數求得後，空間中的任一點的值都可求得。再將整個定義域依據物理變化的幅度分割次定義域，將物理變化較大的次計算域做加強，以增加

(3)

計算效率。二、波譜法理論若有一微分方程如下： L(u)=0 （1）其中 L 為一包括空間與時間的操作因子（operator），於該微分方程式中，我們將 u 以一基函數（trial function）來展開，即

∑

∞ = = 0 n n n a u φ 。而其中an為待定係數，φn 為基函數。再將微分方程式乘上一權重函數ϕ_k（weighting function）後，對空間加以積分而得一新的微分方程。

(

)

∫ ∑

L a_nφ_n ϕ_kωdx = 0 , k=0,1,2,… (2) 其中ω為一權重係數。再搭配邊界條件或起始條件求解，即可求得式中所有未知之待定係數a 。在本n 研究中，取基函數與權重函數相同，皆取 Chebyshev 多項式。任何問題其定義域轉換到 [-1,1] 後，可以 Chebyshev 多項式對變數展開 ( )

∑

( )

_∫

− ( ) ( ) ( ) ∞ = = = 1 1 0 2 ˆ ˆ u x T x w x dx c u x T u x u k k k k k k _π （3）其中權重函數

( )

2 1 1 x x w − = 。將函數 u（x）的導數以 Chebyshev 多項式展開，則函數導數根據（3）可以表示成 ( )

∑

∞ = = ′ 0 1 ˆ k k k T u u (4) 其中 ( )

_∑

∞ = + + = = odd k p k p p k k pu c u 1 1 ˆ 2 ˆ (5) ( )1 u 代表 u 的一階導數。在各子定義域之相交處其一階微分的值需相同。若為二維的定義域，由於τ scheme 並不在空間裡取點，故需將邊界條件以權重餘數法的觀念將邊界條件積分，但積分後所求得之方程式中會出現相依的方程式。為了避免這個問題，本研究於二維空間的連接面取點，以避免相依方程式所帶來之問題。所以為了符合方程式的需要，以 collocation 的觀念取點。其方式如下 N j xj π cos = （6）其中 j 為我們所選取的點的數目，N 為我們所選取的基函數所展開的數目，即為基函數的自由度。令在這些點的未知函數值相同，並且其連接點的 n-1 階導數相等。即 j j X X₁_, = ₂_, 1≤ j≤ N−1 （7） j j 2, , 1 =Γ Γ 1≤ j≤ N−1 （8）其中 Xi,j為連接面所取的點的值，而 Γi,j為連接面所取的點的導數的值。也可以說是於連接面所取的點的值必須是連續（continue）而且平滑的（smooth）。理論結果經測試波浪方程式與 Poisson 方程式後，數值程式之誤差可表示為： 504 . 4 2285 . 6 − = N error N:mode數(子定義域為 2 時) 7723 . 1 1235 . 0 − = E error E:子定義域數(mode=4) 因此收斂非常快。隨著空間分割的數目越多，其精度越高。例如一維之常微分方程式固定 8 個 mode 時其誤差與次定義域的關係式為 9306 . 1 0119 . 0 − = E error _{。而增到 16 個 mode} 時其誤差與次定義域的關係式為 196 . 1 0904 . 1 − = E error 。而 Poisson 方程式固定 3×3 個 mode 時其誤差與次定義域的關係式為 1.7723 1235 . 0 − = E error 。可發現當

(4)

固定 mode 時其誤差的遞減隨著次定義域數目的-1 至-2 次方在遞減。故當空間切割的越多其收斂也越快速，即更能增加計算時的精度。程式發展完成後，以下將其應用於淡水河的各項檢定中，即先用數值模擬水理狀況，再用參數檢定程式來檢定目標參數 1.1實際地形之模擬：已將淡水河河口自關渡橋到河口段（共10 公里）切成長50公尺、寬20公尺、深1 公尺之三維格網，此格網現已放於教學網站上，可供任何人使用 X 方向單位：百公尺 Y 方向單位：百公尺 Z 方向單位：公尺數值水理模擬即以此地形為主，配合在關渡橋與河口測出的邊界流速，來執行這計算域中的水理模擬。 1.2 延散係數 利用上述的數值程式，配合現場測量的結果，延散係數之成果也初步檢定出來（檢定的程式與方法見劉格非，蒲子超， 1998）。延散係數檢定之結果見下圖：

Ex= -143u2_{+ 187u + 14.3}

0 20 40 60 80 -0.2 0 0.2 0.4 0.6 數列1 多項式 (數列1) Ey = -15.2u2 - 100u + 58 0 10 20 30 40 50 60 70 -0.2 0 0.2 0.4 0.6 數列1 多項式 (數列1) 上圖之結果可以下二式表示： 58 100 2 . 15 3 . 14 187 143 2 2 + − − = + + − = u u E u u E y x (m2/S) Ex, Ey單位為m _S 2 , u為主流向方向之流速，單位為 S m _{，請注意資料中有}_u _為負值時，表漲潮時流速向內陸之測量值。大致上來說，延散係數與速度約成線性關係。 1.3 氨氮反應係數 我們採用四步驟機制來模擬河中之反應，即有機氮→氨氮→亞硝酸→硝酸。再分別以實驗檢定反應係數，淡水河河口之有機氮與氨氮在我們於現場取得樣品中（24組），其含量約在2ppm~4ppm之間，而亞硝酸之濃度在0~0.2ppm，而硝酸往往要在樣品培養三天之後濃度才會變大，因此在現場應可不考慮其值，而有機氮到氨氮之反應係數約為 0.4~0.8ppm/day，而氨氮到亞硝酸之反應係數約為0.1ppm/day。我們採用抑制機制型硝化模式為基礎。自試驗中得到淡水河硝化菌最大成長率與 Nitrosomanos 菌隨鹽份之成長率變化為 ) 3163 . 0 0036 . 0 ( logµ_m =µ_m₅_.₃₆ − x− 。其中x 為鹽度(g/l)，µ_m₅_.₃₆為Nitrosomanos菌在x =5.36時之最大成長率。 二、教學網頁 為了將研究成果開放給大眾，我們已

(5)

將上述成果放在「擴散污染教學網頁」網址為 140.112.12.122/diffusion/，網頁有分基礎教學與現場測量資料，網頁架構如下「污染擴散教學網」首頁 http://140.112.12.122/edu/edu01.htm 話說擴散動動腦！淡水河… 水質模擬解析模擬基本觀念線源與面源有限域點源無限域點源面源線源濃度分佈圖固定點濃度固定點濃度濃度分佈圖現地實驗水質資料淡水河地形資料水質參數檢定三維參數檢定二維參數檢定彎曲渠道直渠道河道圖及地形剖面圖歷年實驗成果第一關水質模式其中基本教學的部分又分點源、線源、面源和數值部分，學生可以學到基本理論和互動式運算機制，部分功能，我們以圖檔顯示如下：教學網頁，學生學習線源公式之應用。數值計算網頁，學生可利用本計劃發展出的數值程式來模擬，本頁介紹基本架構數值模擬結果馬上以圖形方式顯出網頁中任何程式需輸入資料時，都有一個專門的資料輸入頁，詳細解釋各項輸入需求、條件還有格式網頁中還有水質模擬與參數檢定程式，學生輸入資料後，可以自途中之到檢定結果計劃中所有現場測量的資料，也都已放到網上，若點選現場測量，並選擇時間與斷面，則可看到許多我們自現場測得的各項參數隨時間與空間的變化圖，下圖顯示的為導電度在不同斷面之測量值與流速剖面圖，

(6)

現場第一斷面在不同時間的導電度深度剖面圖，測量時間為1997/10/26 現場第一斷面在不同時間的流速深度剖面圖，測量時間為1997/10/26 三、結論本計劃經三年，執行現場全潮測量（16 小時）共計五次，並配合需求建立水理數值模式，水質模擬模式與參數檢定模式，再加上互動式教學材料，完成了「擴散污染教學網頁」，可供學生使用甚至業界參考。參考文獻 1. 楊德良（1994）：“三維淺水流場基本解與邊界積分法之研究”國立台灣大學土木工程學研究所 2. 郭佳虹（1997）：“波譜法應用於擬三維河海水理模式”，國立台灣大學土木工程學研究所 3. 簡紹欽（1998）：“大區域海岸水域非線性波浪數值模式之研究”，國立台灣大學土木工程學研究所 4. 葛家豪（1998）：“邊界元素法對三維水沖激問題之數值方法”國立台灣大學造船及海洋工程學研究所 5. 陳效吾（1996）：“寬頻元素法應用於火焰燃燒之研究”，國立台灣大學機械工程研究所 6. 汪吉昌（1996）：“交錯格點應用於寬頻元素法解二維不可壓縮 Navier-Stoke 方程式之研究”，國立台灣大學機械工程研究所 7. 劉格非，蒲子超，1998“淡水河河囗段延散係數之實測與推估模式之建立”。第九屆水利工程研討會論文集，中壢，七月，C69-C76 8. Orszag, S. A. (1970): “Transform method for calculaion of vector couple sums. Applicaton to the spectral form of the vorticity equation.”,J. Atmosph.Sci. 27,pp890-895.

12. Orszag, S. A. (1969): “Numerical methods for the simulation of turbulence.”Phys Fluids, Suppl. II, 12,pp250-257

利用有限元素波譜法建立河海三維水理水質數值模式(3/3)

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