第三章 有限元素數值模式與套裝軟體 ANSYS 分析
3.1 二維軸對稱有限元素模式
本論文研究之現地案例與數值模擬對象主要為圓柱形的單樁,因 此整個基礎結構系統是以軸對稱的型式對稱於基樁的中心軸線,如圖 3.1 所示即為一基樁埋入土壤中的情況。根據這種假設,則可使用二 維軸對稱之有限元素模式來模擬[4],圖 3.2 即是對應於圖 3.1 之結構 系統所建立之軸對稱有限元素模式的右半邊前視平面圖。在對稱中心 線上的節點是以受輥(roller)支承束制的方式來模擬實際基樁在該線 上之質點僅能作垂直方向上之運動,在此模式中所使用的有限元素為
m m 0 . 1 1 .
0 ×
之4 個節點的線性元素,則令N
1、N
2、N
3及N
4分別代表元 素中節點1 到節點 4 的形狀函數,如圖 3.3 所示。則這些形狀數的方 程式為:) )(
(
) )(
) ( , (
1 2 1 2
2 2
1
r r z z
z z r z r
r
N − −
−
= −
………(3.1)) )(
(
) )(
) ( , (
1 2 1 2
2 1
2
r r z z
z z r z r
r
N − −
−
− −
=
………(3.2)) )(
(
) )(
) ( , (
1 2 1 2
1 1
3
r r z z
z z r z r
r
N − −
−
= −
………(3.3)) )(
(
) )(
) ( , (
1 2 1 2
1 2
4
r r z z
z z r z r
r
N − −
−
− −
=
………(3.4)基 樁 長 度
L
基樁
基樁半徑
r
基樁的中心線
P(t)
土壤 地表
圖3.1 二維軸對稱有限元素法中基樁與土壤關係及施力狀況示意圖。
N
v 元 素圖3.3 有限元素數值模型所使用之軸對稱 4 節點元素。
3.2 動力積分方程式
將控制前述基樁與土壤結構系統行為之微分方程式以有限元素 法來離散化後,可得如下在任一瞬間的動力平衡矩陣式:
M U &&
t+ C U &
t+ KU
t= P
t ………(3.5)式中令
U
t、
U& 及t U&& 分別代表質點在時刻tt
時之位移、速度與加速度,且 M
、
C 與 K 分別為該結構系統之質量、阻尼和勁度矩陣,P
t則是在 該時刻 t 時於各有限元素節點上之外加作用力的大小。為了簡化計算 起見,本研究不考慮阻尼效應,即 C = 0。上式中,各矩陣亦可表示 為:z
2z
1r
1r
2r z
4
3
1 2
M = ∫ ρ N
TN dv
………(3.6)
K = ∫ B
TD B dv
… ………(3.7)
P = ∫ N
TP ( t ) dA
………(3.8)上式中
ρ
為密度;∫ ⋅ dv
為對定義域之體積積分;∫ ⋅ dA
為受力區域之面積積分;B 為 應 變 - 位 移 關 係 矩 陣 ( strain-displacement relationship
matrix);
D 為材料性質矩陣(material property matrix);
N 為有限元素模式中之形狀(或內插)函數矩陣。
令應力向量為
T r z
r }
{σ σ σθ τθ σ=
而應變向量為
ε={εr εz εθ γrθ}T 則
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂
∂ ∂
∂
∂
∂
∂
∂
∂
∂
=
N N
N N
N N
N
N r
N r
N r
N r
N
z N z
N z
N z
N
r N r
N r
N r
N
B
4 4
3 3
2 2
1 1
4 3
2 1
4 3
2 1
4 3
2 1
0 0
0 0
0 0
0 0
0 0
0
0
+
− +
−
− +
− +
− +
− +
−
− +
− +
− +
− +
−
=
) 1 ( 0 2
0 0
) 0 2 1 )(
1 (
) 1 ( )
2 1 )(
1 ( ) 2 1 )(
1 (
) 0 2 1 )(
1 ( ) 2 1 )(
1 (
) 1 ( )
2 1 )(
1 (
) 0 2 1 )(
1 ( ) 2 1 )(
1 ( ) 2 1 )(
1 (
) 1 (
ν ν
ν ν ν
ν ν ν
ν
ν ν ν
ν ν
ν ν ν
ν
ν ν ν
ν ν
ν ν ν
ν ν
E E
E E
E E
E
E E
E
D
令
∆ t
表示動力積分時所使用之時間步幅[16],在進行分析時,本研究 所用的時間步幅∆t = 2 . 225 × 10
−5秒,則利用中央差分法可得在時刻t
之速度及加速度方程式為:(
t t 1)
t
U U
U = −
−t 1
&
∆ …….……….………(3.9)(
t t t 1)
t
U U U
U =
+− 2 +
−t
1
2 1
&&
∆ ….……….………(3.10)式中:Ut−1 = U
( t − ∆ t )
( ) t
t
U
U
=( t t )
1
t+ =
U
+∆
U
代入方程式(3.5)後,整理可得下述反應矩陣方程式
(
t t)
t t 1 12 1
t+ =
t M
−P
−KU
+2 U
−U
−U ∆
.….……….…(3.11)由上式可利用某瞬間與之前的位移及外力情況而求出下一瞬間之 反應[22-24]。
3.3 商用有限元素分析軟體 ANSYS 之三維實體元素模式
為了驗證上述軸對稱有限元素模式的正確性,也為了更加逼近真
實結構系統的反應,本研究亦另外使用了商用套裝軟體ANSYS 來進 行單樁基礎的暫態動力反應模擬。其中,模擬基樁所使用之元素為 SOLID92 四面體立體元素(Tetrahedral Solid Element),如圖 3.4 所 示。令
N
1、N
2、N
3及N
4分別代表 SOLID92 元素中節點 1 到節點 4 的形狀函數(shape functions),則在自然座標系統(natural coordinate system) {ξ , η , ζ
} 下之形狀函數可表示為[25]:
2
) 1
(
1
ξ η ζ − −
−
= −
N
………(3.12)
2 ) 1 (
2
ξ
= +
N
………(3.13)
2 ) 1 (
3
η
= +
N
………(3.14)
2 ) 1 (
4
ζ
= +
N
………(3.15)以ANSYS 軟體來進行模擬分析時,基樁及土壤依四面體立體元素做 網格切割,該軟體本身提供了自動網格生成的功能。針對本基樁問 題,該軟體所產生之有限元素網格如圖 3.5 所示。而圖 3.6 所示即為 此網格之局部放大圖。中間密集之柱狀網格部份即為基樁。周圍之大 型網格則為土壤。
圖3.4 商用軟體 ANSYS 中之 SOLID92 四面體立體元素示意圖。
圖3.5 以 ANSYS 之四面體立體元素來模擬單樁與土壤系統之網格示 意圖。
1
2 4
3
圖3.6 以 ANSYS 之四面體元素來模擬基樁系統之網格局部放大圖。
3.4 衝擊荷重之定義
在有限元素數值模擬中,模擬基樁之音波回音或敲擊反應檢測
時,其衝擊錘敲擊基樁的外加動力荷重
P(t)為一作用於樁頂中心的垂
直向均佈荷重,P(t)常以半個正弦波函數來模擬。其形狀如圖 3.7 所 示,其定義為:t P
t
P ( ) =
0sin ω
當0 ≤ t ≤ T
D ….…….…………... (3.16)0 ) ( t =
P
當T
D≤ t
….…….…………... (3.17)其中, T
ω = π ………
(3.18)0.00E+00 2.00E+03 4.00E+03 6.00E+03 8.00E+03 1.00E+04 1.20E+04
0.00E+00 2.00E-04 4.00E-04 6.00E-04 8.00E-04 1.00E-03 1.20E-03 1.40E-03 1.60E-03 1.80E-03 2.00E-03
時間 T
DP0
作 用 力
圖 3.7 模擬衝擊外力所使用之半個正弦波荷重的示意圖。
且
T
D為外加動力荷重P(t)之作用時間或延時(duration),P
0 為此動 力荷重之尖峰值。此衝擊荷重的位置,是在基樁頂面且在中心軸上;其作用寬度為
W
L,作用方向為垂直向下,尖峰值為P
0且均勻的分佈 在受力區域內。本研究初步採用的數值如下述:
P
0 = 9000N / m
2W
L = 0.0508m
T
D =1 . 5 × 10
−3sec
(或 1.5ms)
3.5 驗證用之基樁的幾何組態
為了能較完全地模擬現場基樁之檢測反應,數值模式中之參數應 儘量符合現場之土壤環境及材料性質。本節即定義與實際案例相似長 度之兩種基樁,其長度分別為14 公尺及 16 公尺,其直徑皆為 1.5 公 尺。而土壤面則位於樁頂面以下0.5 公尺。假設土壤與基樁皆為均質 材料。在進行衝擊反應檢測時,衝擊錘敲擊的位置在於基樁頂面之中 心點處(施力之模擬詳見前一節);接收器則置放在樁頂面上介於基 樁中心點與基樁邊緣的二分之一距離處,衝擊錘施力之接觸樁頂時間 為
1 . 5 × 10
−3秒(即1.5 ms)。混凝土基樁與土壤之材料性質則如下所 述:混凝土的材料性質:
楊氏係數
E
c : 3.31 ×1010 N/m 2
浦松比ν
c : 0.2密度
ρ
c : 2300kg/m 3
土壤的材料性質:(模擬現地土壤性質)楊氏係數
E
s : 4.20 ×108 N/m 2
浦松比ν
s : 0.4密度
ρ
s : 2070kg/m 3
上述混凝土所對應之桿件縱波波速為
V
bar = 3800 m/s,剪力波波速 為2450 m/s,表面波波速為 2100 m/s。本研究有限元素模式中所使用之網格大小取為0.1m × 0.1m,此大小是進行收斂研究後之結果。本 研究曾將網格細分至0.05m × 0.05m 及 0.01m × 0.01m,所得之曲線與 網格為0.1m × 0.1m 所得者幾近相同,因此可知原來 0.1m × 0.1m 之 網之網格大小的精度,已收斂至滿意之結果。
3.6 兩種有限元素模擬驗證結果之比較
為了驗證本研究中所使用的二維軸對稱有限元素模式與商用套 裝軟體ANSYS 之正確性,本研究首先對上述之兩種基樁進行模擬之 研究,此兩種案例之基樁長度分別為 14m 及 16m。以下即提出這兩 種模式在時間領域以及頻率領域下之初步分析結果與比較。
3.6.1 時間域結果之比較
首先考慮長 14m 之無樁帽單樁。模擬現地檢測之衝擊力的歷時 曲線,如圖3.8 所示。其衝擊力之接觸時間為 1.5
ms,施力大小幅度
之尖峰值為 9000。由圖中可清楚看出當時間為 0 秒時,衝擊錘開始 接觸到樁頂面,直至達到尖峰值P
0 後開始回彈,衝擊力逐漸遞減,最終回歸至零點。
圖 3.9 所示即為使用二維軸對稱模式與 ANSYS 實體元素模式所 得之位移反應歷時曲線之比較圖。圖3.10 則為圖 3.9 之部份放大圖。
從圖中可以較清楚地標出樁底之反射波的抵達時間為
PD = 7.5 ms
由此可反推出樁長,其計算如下:
t V s m s m
L
p14 . 3
2
3800 10
5 . 7 2
3
× =
= ×
×
= ∆
−
誤差百分比為
1 . 8 14
14 3 .
14 − =
%由圖上亦可看出兩種有限元素模式的結果相當的接近。由此可互相驗 證此有限元素模式之正確性。
接著考慮長 16m 之無樁帽單樁的案例。圖 3.11 所示即為使用軸 對稱元素與ANSYS 之實體元素模擬所得之位移比較圖。為了更清楚 展現反射波的結果,本文將圖 3.11 作局部放大而得圖 3.12。由圖中 可清楚地標出樁底反射波的抵達時間為
PD = 8.66 ms,藉此可推估出
樁長如下:
t V m
L
p16 . 5
2
3800 10
66 . 8 2
3
× =
= ×
×
= ∆
−誤差百比為
2 . 8 16
16 5 .
16 − =
%針對圖3.10 所示之位移反應圖,圖 3.13 所示即為此 14m 長之無 樁帽單樁的速度反應曲線。將此曲線經過低通濾波的技術處理,可得 圖 3.14 所示之清楚曲線。再從該圖中可清楚地標示出樁底反射波的 抵達時間為
PV=7.4 ms,則樁長可推估為:
L 14 . 1 m 2
3800 10
4 .
7 ×
3× =
=
−14
1
.
14 −
同理,可將對應於圖 3.12 所示位移曲線做處理,即得圖 3.15 所 示之 16m 長之無樁帽單樁的速度反應曲線。再將圖 3.15 經過低通濾 波之技術處理可得圖 3.16 所示之清楚曲線。從該曲線中可以很清楚 地標示出樁底反射波的抵達時間為
PV = 8.46 ms,則樁長可推估為:
L 16 . 1 m 2
3800 10
46 .
8
3× =
= ×
−誤差百分比為
0 . 4 16
16 1 .
16 − =
%0.00E+00 2.00E+03 4.00E+03 6.00E+03 8.00E+03 1.00E+04 1.20E+04
0.00E+00 5.00E-04 1.00E-03 1.50E-03 2.00E-03 2.50E-03 3.00E-03 3.50E-03 4.00E-03
時間(秒) 二維軸對稱
P
0作 用 力 振 幅 大 小
圖3.8 模擬 14m 長之單樁所受之衝擊力歷時曲線圖。
-1.6E-09 -1.4E-09 -1.2E-09 -1E-09 -8E-10 -6E-10 -4E-10 -2E-10 0 2E-10
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1
時間 (秒)
ANSYS 二維軸對稱 位
移︵ m︶
圖3.9 以二維軸對稱元素與 ANSYS 實體元素模擬 14m 長之單樁的位 移歷時曲線。
-1.6E-09 -1.4E-09 -1.2E-09 -1E-09 -8E-10 -6E-10 -4E-10 -2E-10 0 2E-10
0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 0.018 0.02
時間 ( 秒 )
ANSYS 二維軸對稱
PD
位 移︵
m︶
圖3.10 以二維軸對稱元素與 ANSYS 實體元素模擬 14m 長之單樁的 位移曲線局部放大圖。
-1.60E-09 -1.40E-09 -1.20E-09 -1.00E-09 -8.00E-10 -6.00E-10 -4.00E-10 -2.00E-10 0.00E+00 2.00E-10
0.00E+00 1.00E-02 2.00E-02 3.00E-02 4.00E-02 5.00E-02 6.00E-02 7.00E-02 8.00E-02 9.00E-02 1.00E-01
時間 ( 秒 )
位移 ( m )
ANSYS 二維軸對稱
圖3.11 以二維軸對稱元素與 ANSYS 實體元素模擬 16m 長之單樁的 位移歷時曲線。
-1.60E-09 -1.40E-09 -1.20E-09 -1.00E-09 -8.00E-10 -6.00E-10 -4.00E-10 -2.00E-10 0.00E+00 2.00E-10
0.00E+00 2.00E-03 4.00E-03 6.00E-03 8.00E-03 1.00E-02 1.20E-02 1.40E-02 1.60E-02 1.80E-02 2.00E-02
時間 ( 秒 )
位移 (m )
ANSYS 二維軸對稱
PD
圖3.12 以二維軸對稱元素與 ANSYS 實體元素模擬 16m 長之單樁的 位移曲線局部放大圖。
-1.50E-06 -1.00E-06 -5.00E-07 0.00E+00 5.00E-07 1.00E-06 1.50E-06 2.00E-06
0.00E+00 2.00E-03 4.00E-03 6.00E-03 8.00E-03 1.00E-02 1.20E-02 1.40E-02 1.60E-02 1.80E-02 2.00E-02
時間 ( 秒 )
速度 (m / sec)
ANSYS 二維軸對稱
圖3.13 以二維軸對稱元素與 ANSYS 實體元素模擬 14m 長之單樁的 速度歷時曲線。
-1.00E-06 -5.00E-07 0.00E+00 5.00E-07 1.00E-06 1.50E-06 2.00E-06
0.00E+00 2.00E-03 4.00E-03 6.00E-03 8.00E-03 1.00E-02 1.20E-02 1.40E-02 1.60E-02 1.80E-02 2.00E-02
時間 ( 秒 )
速度 ( m / sec )
ANSYS 二維軸對稱
PV
圖3.14 以二維軸對稱元素與 ANSYS 實體元素模擬 14m 長之單樁經 低通濾波處理後的速度曲線圖。
-1.50E-06 -1.00E-06 -5.00E-07 0.00E+00 5.00E-07 1.00E-06 1.50E-06 2.00E-06 2.50E-06
0.00E+00 2.00E-03 4.00E-03 6.00E-03 8.00E-03 1.00E-02 1.20E-02 1.40E-02 1.60E-02 1.80E-02 2.00E-02
時間 ( 秒 )
速度 ( m / sec )
ANSYS 二維軸對稱
圖3.15 以二維軸對稱元素與 ANSYS 實體元素模擬 16m 長之單樁的 速度歷時曲線。
-1.00E-06 -5.00E-07 0.00E+00 5.00E-07 1.00E-06 1.50E-06 2.00E-06
0.00E+00 2.00E-03 4.00E-03 6.00E-03 8.00E-03 1.00E-02 1.20E-02 1.40E-02 1.60E-02 1.80E-02 2.00E-02
時間 ( 秒 )
速度 ( m / sec )
ANSYS 二維軸對稱
PV
圖3.16 以二維軸對稱元素與 ANSYS 實體元素模擬 16m 長之單樁經 低通濾波處理後的速度曲線圖。
3.6.2 頻率域結果之比較
首先考慮針對 14m 長之單樁的檢測反應。如圖 3.14 所示為其速 度反應歷時曲線。將其訊號以快速傅立葉轉換送至頻率域再將如圖 3.8 所示之施力歷時曲線亦送至頻率域。將二者相除,即得對應於此 基樁之力學導納曲線,其結果將如圖3.17 所示。該圖中之曲線即為 針對14m 長之基樁,使用二維軸對稱有限元素模式與商用 ANSYS 軟 體之實體元素模式所得之力學導納曲線比較圖。由圖3.17 可知,在 低頻率區域(即0~1000 Hz),可以很容易的找出大小級數相當且重 複出現的週期性波峰,如此可依第二章所列之公式推估出基樁之長 度。相關之計算如下述:
( )
5 . 2 137
274 549 − =
=
∆f
Hz
m
Hz m f
L V
c13 . 8
5 . 137 2
3800 sec
2 =
= ×
= ∆
誤差百分比為
1 . 3 14
14 8 .
13 − =
%同理,針對 16m 之基樁可進行類似之處理,而得到頻率域的曲 線。圖3.18 所示即為長度 16m 之無樁帽單樁,使用二維軸對稱有限 元素模式與商用ANSYS 軟體之實體元素模式所得之力學導納曲線比 較圖。由圖3.18 可知,在低頻率區域(即 0~1000 Hz),可以很容易
的找出大小級數相當且重複出現的週期性波峰,如此可依上述作法推 估出基樁之長度。其相關之計算為:
121 2
362 604 − =
=
∆f
Hz
m f
L V
c15 . 7 121 2
3800
2 =
= ×
= ∆
誤差百比為
1 . 9 16
16 7 .
15 − =
%0.00E+00 5.00E-11 1.00E-10 1.50E-10 2.00E-10 2.50E-10 3.00E-10 3.50E-10 4.00E-10 4.50E-10 5.00E-10
0.00E+00 2.00E+02 4.00E+02 6.00E+02 8.00E+02 1.00E+03 1.20E+03 1.40E+03
頻率 ( Hz )
力學導納值 ( V / F )
ANSYS 二維軸對稱
2 ∆ f=549-274=275
圖3.17 以二維軸對稱元素與 ANSYS 實體元素模擬 14m 長之單樁的 力學導納曲線比較圖。
圖3.18 以二維軸對稱元素與 ANSYS 實體元素模擬 16m 長之單樁的 力學導納曲線比較圖。
3.7 小結
針對本章使用二維軸對稱有限元素模式及商用 ANSYS 軟體之實 體元素模式,進行基樁動力反應之模擬驗證比較中發現:
(1)在時間域上,由位移圖來標示出樁底反射波之波抵時間,再反算出 基樁長度之誤差百分比都很小,顯示出此方法對於檢測十幾公尺 之單樁的精準度相當高。
(2) 同樣在時間域上,若使用速度曲線,亦可輕易標出樁底反射波之 波抵時間,且推算出的基樁長度亦有很好的可靠性與精度。
0.00E+00 5.00E-11 1.00E-10 1.50E-10 2.00E-10 2.50E-10 3.00E-10 3.50E-10 4.00E-10 4.50E-10 5.00E-10
0.00E+00 2.00E+02 4.00E+02 6.00E+02 8.00E+02 1.00E+03 1.20E+03 1.40E+03
頻率 ( Hz )
力學導納值 ( V / F )ANSYS 二維軸對稱
2
Δf=604-362=242
範圍內)可以輕易地尋得週期性的一些波峰。而使用這些週期性 波峰間的頻率差
∆ f
,可反算出基樁的長度。此方式所得之樁長資 料與實際長度的誤差很小。(4) 綜合上述,二維軸對稱有限元素模式與 ANSYS 之三維實體有限 元素模式所得之結果,無論是在時間域或頻率域上,皆非常相近,
標示出之樁底反射波的波抵時間及頻率上之頻率差亦都相當接 近。所以,此兩種有限元素模式皆能有效地模擬無樁帽單樁之反 應。