第二章 文獻回顧
2.2 數值方法
2.2.2 波譜元素法
一般除了有限元素法為常見的數值模擬方法之外，Patera 於 1984 年結 合有限元素法 (Finite element method, FEM) 和擬頻譜法 (Pseudo-spectral method)，提出波譜元素法 (Spectral element method, SEM)，此技術基於運 動方程式之弱式(Weak formulation)，並結合有限元素法的彈性與擬頻譜法 的精度。網格主要以六面體元素組成，元素中波場的離散化透過高階拉格 朗日內插 (High-degree Lagrange interpolants)，積分的計算則經由 GLL 積分 法(Gauss-Lobatto-Legendre quadrature)，結合上述方法可以形成對角質量矩
f=− M∙∇δ(x − xs)S(t) (2.16) xs為點作用力的位置，δ(x − xs)為位於 xs之 Dirac delta distribution，S(t)為震 源時間函數(source time function)。
運動方程式(2.14)必須求解於自由應力邊界條件之表面∂Ω：
T∙n=0 (2.17)
n是垂直於自由表面之單位向量。
圖 2.15 有限地球模型示意圖，震源 xs可置於模型內任意位置(Komatitsch and Tromp, 1999)
然而震波能量傳遞至人造邊界 Γ 需作吸收處理，為此才可模擬一個半 無限域空間，所以經常採用平行軸方程式來降低波場邊緣震波的能量 (Clayton and Engquist, 1977)：
T∙n=ρ[vn(n∙∂ts)n+v1(t̂1∙∂ts)t̂1+v2(t̂2∙∂ts)t̂2] (2.18) t̂1, t̂2為相互垂直的單位向量，且與人造邊界 Γ 相切；vn為n方向傳遞之準壓 縮波波速；v1, v2分別為t̂1, t̂2方向偏振之準剪力波速度。此邊界條件能完美 地吸收單一方向垂直邊界之衝擊，但是對於擦過邊界傳遞的波卻無法有效 降低影響。
於 SEM 只考慮積分項表示的運動方程式，(2.14)須與一任意的向量 w 作內積，於 FEM 稱為測試函數，再對整個模型體積 Ω 作分部積分，最後加 入上述提到的邊界條件，得到：
ρw∙∂t2sd3x
可以發現 f 項已利用 Dirac delta distribution 之特性直接積分；而使用弱式還 有另一好處是分部積分式中於自由應力表面(2.17)之積分可以抵銷。
2. 網格設計
SEM 所要模擬的空間域 Ω 將被分為數量 ne相互不重疊之元素 Ωe， e=1,…,ne使得Ω = ⋃ne=1e Ωe (如圖 2.16)，同樣的吸收邊界 Γ 也被分為 nb數量 之元素 Γb，b=1,…,nb使得Γ = ⋃nb=1b Γb。
圖 2.16 網格分割示意圖(Komatitsch and Tromp, 1999)
一般來說，FEM 可以使用四面體或六面體作為體積元素，但 SEM 限制 只能使用六面體，同樣的定義吸收邊界其表面元素形狀為四邊形，接下來 的部分將會說明如何將這些元素參數化。
3. 映射函數(mapping function)：轉換物理座標至標準區間
經由座標轉換，每一六面體元素 Ωe之物理座標將轉換成參考座標 ξ = (ξ,
圖 2.17 六面體元素之幾何形狀示意圖，空心代表可以省略(Schuberth, 2005) Na為形狀函數，由一階或二階 Lagrange 多項式之乘積；若多項式的階次為 nL，將有 nL+1 個控制點，-1≤ξi≤1, i=1,…,nL定義如下：
linL(ξ)= ξ− ξj ξi− ξj
nL
j=0 j≠i
(2.21)
若是 nL=1 即可產生直線之形狀函數，nL=2 即可產生曲線的形狀函數(如圖 2.18)
(a) (b)
圖 2.18 一階與二階形狀函數示意圖 (a) 兩個節點 (b) 三個節點 (Schuberth, 2005)
所以，一個節點數為 8 的三維元素，其形狀函數可以用 nL=1 之 Lagrange 多 為了計算 Jacobian 首先須對式(2.20)偏微分，得到 Jacobi-matrix(J)，再取行 列式值：
4. 表示元素格點之內插函數 (interpolation function)
在模型網格化與使用弱式表示運動方程式後，下一步便是定義基函數 來表示元素內未知的位移向量。因此解析式(2.19)，現在必須藉由離散函數，
使用適當的插值方法估計離散點的值。傳統的 FEM 使用低階多項式表示元 素之幾何與向量，而 SEM 也用低階多項式表示元素之幾何，但所有未知的 函數則採用高階 Lagrange 多項式來內插，典型的 SEM 常使用 4~10 階 Lagrange 多項式作內插，這是 SEM 與 FEM 兩者主要的差別。
再仔細觀察式(2.21)，多項式li的值在座標 ξj剛好是 1 或 0(圖 2.19)，而 在這些離散點以外多項式會有任何可能的值，因為只對這些離散點的值有 興趣，所以並不違反下列的數學推導式。
可以用 Kronecker-delta 來表示這個特性：
linL(ξj)=δij (2.30) 而定義於式(2.21)控制點 ξi, i=0…nL的選擇則是透過 Gauss-Lobatto-Legendre points，也是下列式子之根：
(1 − ξ2)Pn'L(ξ)=0 (2.31) Pn'L就是 nL階 Lagrange 多項式的一次微分(Canuto et al., 1988)。這樣的選擇 是基於 Lagrange 內插與 GLL 求積的結合可使質量矩陣成為對角矩陣，因此 一個二維元素將有(nL+1)2個格點(圖 2.20)，三維元素將有(nL+1)3個格點。
圖 2.19 Lagrange 多項式於[-1, 1]內插，nL=4，有 nL+1=5 個控制點(Komatitsch et al., 2005)
圖 2.20 二維波場離散化示意圖(Komatitsch et al., 2005)
∇g(x(ξ, η, ζ))≈
傳統的 FEM 使用高斯求積分(Gauss quadrature)，SEM 則是採用 GLL 積分 法，結合 GLL 點可形成對角質量矩陣。體積元素 Ωe之積分可估計為： 定義全域(global)與地域(local)的映射關係，可透過 FEM 或分類演算法去搜 尋之，並指派獨一無二的全域座標。一旦映射定義完成，分別在元素間計
M：質量矩陣 C：吸收邊界矩陣 K：勁度矩陣 F：震源項