有限元素法

在數值模擬中，除了有限差分法之外，有限元素法亦是相當常用的一 種技術。其藉由將實際的物理或數學問題模型簡化為許多元素

(element)

的組 合，先針對每個離散後的元素進行計算後，再將所有的元素整合而得到整 體的表現。在操作上，其可以概略的分為六個步驟：

(1)

將問題進行離散

(Discritization)

；

(2)

選擇每個元素中的內插

(interpolation)

方式；

(3)

導用有 限元素法計算式；

(4)

組結

(assembly)

各個離散的元素；

(5)

邊界條件給定

(constrain)

；最後

(6)

求解。

模型離散化 模型離散化 模型離散化

模型離散化：在使用有限元素法進行數值解時的第一個步驟便是將物理或 數學模型離散化，將原本在空間中連續的模型，分割為許多無重疊元素的 組合

(

如圖

2.25

所示

)

。很直觀的可以知道，當我們使用更多更小的元素去 組合成原本連續的模型時，我們將能有更接近真實的解。在圖

2.25

的例子 中，原本是弧狀的邊界，經由採用多邊形元素將其離散化後組成一多邊形 形體，因而有誤差的產生，然而，在有限元素法中是容許具有弧形的元素，

使得邊界的定義上可以完全符合，但這需採用更為複雜的元素作為組合元 件，也需要更高技巧的解決方式。

常被用於有限元素法中的元素，對二維空間問題來說有三角形與四邊 形，又因為在這一元素中所使用的節點數目不同而可分為線性元件

(linear member)

、二次元件

(quadratic member)

及三次元件

(cubic member)

等不同元素

種類

(

如圖

2.26)

。元素種類的採用間接影響了內插方式的選用，亦會對解決

的問題產生不同程度的誤差。例如對於解決二階微分的問題，若以線性元

圖 圖圖

圖 2.2.2.2.25252525 二維二維二維模型離散化示意圖二維模型離散化示意圖模型離散化示意圖 模型離散化示意圖

圖 圖 圖

圖 2.2.2.2.26262626 二維空間問二維空間問二維空間問題常用元素二維空間問題常用元素題常用元素題常用元素

件的元素組成具有弧度的邊界，在邊界周圍所得到的近似值常會向錯誤的 值收斂；而若是同樣的離散方式要解決四階微分的問題，則其值必不收斂

(Kaliakin, 2002)

。因此，在不同的問題情況下，必須選用不同的元素，以期

能達到最佳的結果。

要對於每個元素的特徵進行了解，除了其形狀以及節點數外，還須了 解節點的屬性、在每個節點上的自由度以及所將採用的內插方法。如圖

2.26

中三次元件的三角形元素，在此元素中的節點便有不同的屬性存在，一種 是屬於圍外節點

(exterior)

，一種是內節點

(interior)

。圍外節點

(

圖

2.26

中填

實的位置

)

與內節點

(

圖

2.26

中中心未填實的圓圈處

)

的主要差異在於，圍外 節點將會與其他元素的圍外節點共用或是在邊界上，而內節點則與鄰近元 素皆無關。其次是在節點上的自由度，亦即在節點處所須計算的獨立未知 數，以二維空間問題來說，節點處除了可以有水平方向及垂直方向的移動，

亦可以有轉動方向的移動量

(

如圖

2.27)

，每個節點上的自由度會因為問題的 不同而有所不同。

圖圖

圖圖 2.2.2.272.272727 節點自節點自節點自節點自由度示意圖由度示意圖由度示意圖由度示意圖 內插方式

內插方式 內插方式

內插方式：與節點相關的特徵中，還有一個極需注意的是內插方式。內插 方式說的是一個元素內沒有節點的位置的值是如何以節點的值做表示。對 應於不同的節點數的元素使用，內插的方式可以有線性

(linear)

或二次式

(quadratic)

等高次的內插方式，使得

e

元素中的

u

物理量可表示為

) ( )

(^e

( X ) Nu

_n^e

u = (2. 113)

其中，X 是一空間向量，表示在ｅ元素中任一位置ｕ物理量的座標；N 為 一權重係數矩陣，用來表示元素中不同節點對此一位置物理量的不同權 重；u_n^(e)為一向量，表示節點處物理量ｕ的值。舉例來說，在一個二維空間 中，以單自由度三個節點的三角形元素做離散後

(

如圖

2.28

所示

)

，在ｅ元 素中的ｕ物理量可表示為：

) ( )

( )

( )

(^e

( x , y ) N

_i

( x , y ) u

_i^e

N

_j

( x , y ) u

_j^e

N

_k

( x , y ) u

_k^e

u = + + (2. 114)

其中，N_i,N_j,N_k分別為節點i, j,k的權重係數，與元素中所欲求的物理量座 標位置

(x,y)

有關，以不同的內插方式會有不同的權重係數函數；u_i^(e),u⁽_j^e),u_k^(e)

為在元素

e

中節點處的物理量

u

的值。

圖 圖 圖

圖 2.2.2.282.282828 元素與內插方式示意圖元素與內插方式示意圖元素與內插方式示意圖元素與內插方式示意圖

內插的權重係數函數主要是在其相對的節點上為

1

，其他節點為

0

，而 在這中間便是由內插方式計算。為使這樣的概念得以表達，將問題降至一 維空間來看。圖

2.29 a)

中

i

，

j

表示一個有兩個節點的元素，

b)

中

i

，

j

，

k

表 示一個有三個節點的元素。首先看圖

2.29 a)

，對節點

i

而言，其內插權重在

i

的地方為

1

，在節點

j

的地方為

0

，而在

i

與

j

中間的地方便是線性的內插，

反之，對於節點

j

亦然。再看圖

2.29 b)

，其以二次式為內插方式，但同樣的，

對節點

i

來說，其內插權重在

i

為

1

，在點

j

，

k

的地方為

0

，而在這之間的 便是以二次式進行內插。反之，對其他節點亦然。從這也可知前述提及元 素的選用會間接的影響內插方式的選用原因：欲使用二次式為內插方式在 一個一維元素中至少需要有三個節點，而若使用線性內插方式在一個一維 元素中只需兩個節點。

圖圖

圖圖 2.2.2.292.292929 一維空一維空一維空一維空間問題內插權重係數示意圖間問題內插權重係數示意圖間問題內插權重係數示意圖間問題內插權重係數示意圖 收斂性

收斂性 收斂性

收斂性、、、、協和協和協和協和性與完整性性與完整性性與完整性：對有限元素法的分析來說，其主要產生的誤差性與完整性 一部分是來自於離散化的過程，而有一大部分是來自於內插方式的選用。

內差方式的選用並非是任一選取的，為了讓有限元素法得到的解有單一收

斂

(monotonic convergence)

，亦即當離散時的元素使用越多時，其解的準確

度 會 越 高 ， 內 插 方 式 的 選 用 必 須 符 合 協 和 性

(compatibility)

及 完 整 性

(completeness)

的要求。協和性所要求的可以分為兩項，一是平滑

(smoothness)

要求，亦即在元素內部內插的結果必須是連續的；另一是連續性

(continuity)

要求，亦即在每個元素邊界處

(

或圍外節點

)

對內插函數進行與控制方程式相 同階數微分後，其不可為

0

。由此兩項要求便可知，協和性的要求是要確保 當我們將元素與其周圍元素接合時不致會有不連續的間隙

(gap)

產生。而完 整性所欲保障的是，確保即使當離散過程中所使用的元素數目增加，其將 含有內插函數的值微分至常數的值都不會有所改變，確保這些微分後的常 數在不同的元素數目下都有再現性。當一內插方式能達到此兩點要求便能 使得此有限元素法達到單一收斂的目的。

有限元素法計算式導用 有限元素法計算式導用 有限元素法計算式導用

有限元素法計算式導用：內插方式確認可以使有限元素法的解具有單一收 斂結果後，便可將此內插式代入針對不同問題所得到的控制方程式，再因 應不同的問題採用戈林肯有限元素模型

(Galerkin finite element model)

或變 異有限元素模型

(Variational finite element model)

導出有限元素運算式

(

詳細 過程可參考

Kaliakin, 2002

或

Bickford, 1994)

。導出的有限元素運算式以矩 陣方式表達，其一般式為：

) ( )

( )

( e e

n

e

U f

K = ^{(2. 115)}

其中，K^(e)是第ｅ個元素的常數矩陣；Un^(e)是第ｅ個元素在節點上的未知 物理量向量； f ^(e)是第ｅ元素節點處的＂力＂向量

(

此處力向量的意義並非 是真的力，而是表示已知的作用物理量，只是因為有限元素法起初是應用 於結構的計算，因此習慣上以力矩陣稱之

)

。

組結 組結 組結

組結：以上述的式子我們將可得到每個元素處的區域性

(local)

有限元素運算 式，當每個元素都計算完成後，我們需要將這些區域性的式子依照當初離 散的相互關係組結起來，而得到全域性

(global)

有限元素運算關係式，其一 般是可以表示為：

f

KU

_n

= (2. 116)

其中，K是模型的常數矩陣；U_n 是模型中每個節點上的未知物理量向量，

其未知數的個數與整個模型的全部自由度相同； f _{是整個模型的＂力＂向} 量。這意思也就是說，將所有的未知數整合為Un _{矩陣，依造}U_{n 矩陣中未} 知數擺放的位置將得到的K^(e)適當擺放至相對應的K矩陣位置，而對於相 同位置的係數只需將其相加便可，而 f _{亦是將相對應的} f ^(e)_{係數置入便完} 成組結的步驟。

邊界條件給定 邊界條件給定 邊界條件給定

邊界條件給定：待全域性的有限元素關係式給定後，依照問題不同的條件 給予自然邊界條件

(Neumann boundary condition)

、強制邊界條件

(Dirichlet

boundary condition)

或是兩者皆有的混合邊界條件，代入全域性的有限元素 關係式中，而後以矩陣關係去除多餘的陣列，得到：

r r

n

r

U f

K = ^{(2. 117)}

其中，

r

表示為去除多餘陣列後的全域有限元素關係式。而後求解便可求得 未知數矩陣：

( )

^{r r}

r

n

K f

U =

⁻¹

^{(2. 118)}

而獲得節點上未知數的解。

有限元素法震波模擬 有限元素法震波模擬 有限元素法震波模擬

有限元素法震波模擬：相較於在震波模擬領域文獻豐富的有限差分法，有 發表出的有限元素法相關討論文章極為缺乏，這除了有限元素法的使用較 有限差分法困難外，也因為有限元素法的基本架構已大致發展齊全

(Bathe and Wilson, 1976

；

Zienkiewicz, 1977)

，在震波模擬的使用上只是處理邊界條 件以及導用有限元素計算式的不同

(Kelly and Marfurt, 1990)

。震波模擬在有 限元素法上最大的突破便是為了解決對地層反應的龐大計算量而發展出的 數值計算方式

(Marfurt et al., 1987)

。

早期以有限元素法進行震波模擬的文章有

Lysmer

團隊在

1970

年代以 在頻 率域解得 表面波 波傳行為 的一系 列文章

(Lysmer and Drake, 1971, 1972

；

Drake, 1972a, b)

，這之後又有

Smith

對於有限元素法對實體波波傳以 及地表自由震盪

(Free Oscillations)

情形所做的模擬

(Smith, 1975, 1980)

，而

Marfurt

以半無限域空間為例，探討以有限差分法及有限元素法進行波傳模

擬的準確性

(Marfurt, 1984)

。這之後在震波模擬領域上有限元素法的討論便 失之闕如。直至近年，因為對於複雜的地形效應以及具有高度側向變化的 非均質非均向甚至具有黏滯性地層的探討需求，使得有限元素法逐漸恢復

擬的準確性

(Marfurt, 1984)

。這之後在震波模擬領域上有限元素法的討論便 失之闕如。直至近年，因為對於複雜的地形效應以及具有高度側向變化的 非均質非均向甚至具有黏滯性地層的探討需求，使得有限元素法逐漸恢復

在文檔中 Pseudo-section概念於表面波震測應用之數值模擬探討 (頁 71-78)