行政院國家科學委員會專題研究計畫 成果報告
含裂紋之功能梯度複合材料受熱載重之邊界元素法分析
計畫類別: 個別型計畫
計畫編號: NSC92-2212-E-011-046-
執行期間: 92 年 08 月 01 日至 93 年 07 月 31 日 執行單位: 國立臺灣科技大學營建工程系
計畫主持人: 張燕玲
計畫參與人員: 陳韋廷、廖建璋
報告類型: 精簡報告
處理方式: 本計畫可公開查詢
中 華 民 國 93 年 9 月 20 日
行政院國家科學委員會專題研究計畫成果報告
含裂紋之功能梯度複合材料受熱載重之邊界元素法分析
Boundary Element analysis of the Thermal Stressed Cracks in a Composite with Functionally Graded Materials
計畫編號:NSC 92-2212-E-011-046
執行期限:2003 年 8 月 1 日至 2004 年 7 月 31 日 主持人:張燕玲 國立台灣科技大學 營建系教授
計畫參與人員: 陳韋廷、廖建璋 國立台灣科技大學 營 建系研究生
一、摘要
本研究以邊界元素法,探討含邊裂縫之 單一 FGM 受溫度變化對裂縫尖端熱應力 強度因子(以下簡稱 TSIF)之影響,以及 TSIF 和裂縫長度之關係。本研究所探討之 FGM,其彈力係數
E
、υ
為常數,熱膨脹係數呈函數變化,稱為熱膨脹係數功能梯 度材料(簡稱
α −
FGM)。對於功能梯度材料承受溫度變化之問 題,其全解可分為特解與齊次解兩部份處 理;其中特解部份為
α −
FGM 受溫度變化 所產生,可由 Fourier 級數展開而求得,經 由特解可推知齊次解之邊界條件,再由α −
FGM 之多區域邊界元素法求得齊次 解,特解與齊次解兩者之合即為全解。數值分析之結果,首先在相同的徹體 力下與文獻之結果相比較,經比較後其結 果很接近。進而推廣本程式之使用範圍,
探討
α −
FGM 含邊界裂縫受溫度變化時,其 TSIF 與材料常數及裂縫長度間的關係。
關鍵詞:邊界元素法、熱應力強度因子、
功能梯度材料 Abstract
This thesis applies the method of boundary element method to analyze the functionally graded material(FGM) subjected to constant temperature loadings. The Young’s modulus and Poisson’s ratio of the FGM using in the thesis are constant, while the thermal expansion coefficient
is presented as functionally distributions, called
α − FGM. The thesis is to inspect the effects of the thermal stress intensity factors(TSIF) at the crack tip and the relation among the TSIF, material characteristics, and the crack lengths when the α − FGM containing edge-crack subjected to constant temperature loadings.
For the problem of α − FGM subjected to constant temperature loadings, first at all, the problem can be separated into homogeneous and particular problems. The particular solution is obtained by expanding the first order derivative of thermal expansion coefficient into Fourier series. Then by the inserting the obtained particular solution, the boundary conditions of homogeneous problem can be specified. Then the homogeneous solution is evaluated by the
α − FGM multi-region boundary element method. Finally, the complete solutions can be obtained by summing the homogeneous and particular solutions.
Results in this thesis agree well with those of references that have the same body-force with this thesis. When the proposed program is applied to the various one-layer α − FGM
containing edge-crack subjected to constant temperature loadings, the relations of TSIF with the crack lengths are obtained.
Keywords: Boundary element method, Thermal stress intensity factor, Functionally graded material
二、緣由與目的
複合材料在工業上已被廣為使用,然而 複合材料在交界處由於相鄰材料之差異
性,容易出現應力不連續和應力集中之現 象,而造成材料之破壞。若能降低相鄰材 料間的差異性,將可減少材料之破壞,遂 有 功 能 梯 度 材 料 (Functionally Graded Material,簡稱 FGM)產生。
使用 FGM 的目的在降低複合材料相鄰 兩材料間之差異,以降低熱殘留應力,使 其有效發揮材料之特性。但 FGM 或複合材 料之製造,或使用過程,往往需承受巨大 的溫度變化,而產生很高之熱殘留應力,
因此材料內可能產生裂縫,故使用 FGM 時,評估其損傷或破壞程度是很重要的。
本文探討含邊裂縫之 FGM 受溫度變化 後,對裂縫尖端 TSIF 之影響,以及熱應力 強度因子和裂縫長度之關係,藉由對 FGM 含邊裂縫之研究,深入了解 FGM 受溫度變 化後之行為,以增加 FGM 之了解,並增進 FGM 之經濟效益。
一般功能梯度材料之彈性係數、熱膨 脹係數及熱傳導係數均呈連續之函數曲線 分佈,但亦有楊模數為常數的 FGM,如 Zirconia/nickel FGM【1】 其楊氏模數在 FGM 內幾乎不變,此乃因為鎳合金(nickel alloy)及鋅(Zirconia)有相似之楊氏模數。此 類 FGM 尚 有 MoSi / Al O2 2 3 FGM, TiC/SiC FGM, Zirconia/Steel FGM 等。因 此本文欲考慮:楊氏模數為常數但熱膨脹 係數為一連續函數之功能梯度材料(簡稱
α −
FGM),受到熱載重之龜裂行為。
三、研究方法與成果 1.研究方法:
邊界元素法應用於非均質等向性材 料之 FGM,所須用到之基本解不昜獲得,
且考慮溫度載重作用,相當於物體受到徹 體力作用,因此,若使用邊界元素法分析 時,須修改其基本解,或另外處理徹體力。
而 BEM 的優點是將二維或三維的問題簡 化為一維或二維的問題來處理,可以減少 計算時間,故本文擬以 BEM 分析功能梯度 材料受熱載重作用之龜裂行為。
功能梯度材料承受溫度載重時,於 Navier 之平衡式中表示出溫度之效應與徹 體力之關係為:
b
i*= − b
ic
,i∆ T − c T ∆
,i。因此,若溫度之變化與位置有關,則溫度 之影響可視同徹體力之影響。利用 BEM 處
理徹體力問題時,本文將以宋等【2】之方 法,將含徹體力之問題視為齊性問題及特 解來處理。特解部分,利用 Tang【3】之方 法將徹體力展為 Fourier 級數,先對級數之 各代表項求出對應之無限域之特解,然後 將各項特解疊加即得徹體力作用於無限域 問題之特解。將原為有限域問題之邊界條 件加以變換,原為非齊性方程之問題即可 化為齊性方程問題,然後再運用邊界元素 法處理此齊性問題。
假 設 徹 體 力 *
b
j 在 區 域Ω 內 能 滿 足 Dirichlet 條件【4】;換言之,已知熱膨脹 係 數 函 數 的 一 階 微 分 在 區 域 內 滿 足 Dirichlet 條件。如圖 1,徹體力 *b
j作用在區域
Ω
內,假設區域Ω
包含於一[2a
1×
2a
2]之Ω
f區域內,則徹體力b
*j在Ω
f 中可用雙層傅利葉級數展開【3】:
+
∑∑ ∑+
∑∑=
= = = = =4
1 1 1
4
1 1
0 0
*
l N
n M
m
l nm l jmn l
n l
N
n l jn j
j
K f K f K f
b
並分為(A)常數項;(B)含正弦或餘弦之項;
(C)多調和項(multiple harmonic)等三部分。
y
x
Ω
Ω
f 22a
2a1圖 1 在徹體力作用下之區域Ω示意圖
(A) 常數項
由於徹體力常數
b
j*必須滿足平衡方程 式σ
ij,j+ b
*j=
0,故將應力分量假設為一次 函數,以求得位移場:
− −
= + C x Cxy uxp E 1 2 2
2 1
1
υ υ υ
− −
= + C y Cxy
uyp E 1
2
2 2
1
1
υ υ υ
及曳引力向量:
[ K x n K x K y n ]
t
xp2
10 1(
20 10)
22
1 − + +
= − υ
υ
[ K x K y n K y n ]
t
yp(
20 10)
12
20 22
1 + −
= − υ
υ
(B) 含正弦或餘弦之項徹體力含正弦或餘弦項即
= ∑∑
= = 4
1 1
* l
N
n l n l
j
K
jnf
b
之部分,位移向量及曳引力向量分別為:
2 ) ( 1 2 2 3 3 3 4
1 1 1
n n
n N
n
n p
x f
n f d n f d n f d n u d
+ ′′
− ′′
+ ′
− ′
=
∑
=
2 )
(2 3 1 2 4 2 3
1 3 3
n n n N
n
n p
y f
n f d n f d n f d n u d
− ′′
+ ′′
− ′
=
∑
′=
∑
= ′ − + + ′ −= N
n
n n
n n
p
x k f n k f n
n n f k n f n k t
1
2 2 2 1 1 1 2 1 2 1 2
1 1( )
) 1(
3 4 4 3
2 1 1 2 2 1 1 2
1 1
( ) ( )
1 k f nn k f nn 1 k f nn k f nn
n n
υ υ
υ υ
+ ′′ − − + ′′ − − +
∑
= −′ − +
− +
′ −
= N
n
n n
n n
p
y k f n k f n
n n f k n f n k
t
1
1 2 2 2 1 1 1
1 2 2 2
1 )
(1 ) 1 ( 1
1
υ υ υ
υ
3 4 4 3
2 2 1 1 2 2 1 1
1 1
(k f nn k f nn ) ( k f nn k f nn )
n n
+ ′′ − + ′′ − +
(C) 含多調和項 即
= ∑∑ ∑
= = =
4
1 1 1
* l
N
n M
m
l nm l
j
K
jmnf
b
之部分位移向量及曳引力向量分別為:
)
( 1 3 1 1 1 4
1 1 2 1
nm nm nm N
n
nm M
m p
x f
n f d n f d n f d n u d
+ ′
− ′ + ′
− ′
=
∑∑
= =
)
( 2 1 2 4 2 1
1 2 3 1
nm nm n N
n
nm M
m p
y f
m f d m f d m f d m u d
− ′ + ′ + ′
− ′
=
∑∑
= =
∑∑ [
= = + + +
= N
n
nm nm nm
nm M
m p
x c f n c f n cf n c f n
t
1
2 4 3 1 1 1 2 1 3 1 4 1
1 2)
]
3 3 1 2 1 2 2 3 1 3
1f n c f n c f n c f n
c nm + nm + nm + nm
+
∑∑ [
= =
+ +
+
= N
n
nm nm
nm nm
M
m p
y c f n c f n c f n c f n
t
1
2 1 2 1 4 3 2 4 2 1 1 3
+c3fnm12n1+c2fnm3n2+c3fnm3n1+c2fnm2n2)
]
將(A)、(B)、(C)三部份解疊加,得到 徹體力作用下區域內任一點的位移向量與 曳引力向量,此即問題之特解部份,再將 邊界上由徹體力造成之位移與曳引力場代 入已知邊界條件中,求得齊次解之邊界條 件,並利用邊界元素法求得齊次解。最後,
將齊次解與特解相加,求得問題之全解。
2.研究成果:
考慮一E=2.1×106kg cm2
υ
=0.3且熱 膨脹係數為α ( x ) = ( x w )
nα
m+ ( 1 − ( x w )
n) α
c 之α
−FGM,如圖 2 所示,α (0) = α
c;( ) w
mα = α
;其間0≤ x ≤ w
,則膨脹係數呈 梯度變化。此α
−FGM 含有長度為a
之邊裂 縫,且受到溫度變化T = 1
oC
之作用。(A)
α ( ) x
為線性變化之分析設
α
c= ×1 10−5 oC ,α
m =1.38×10−5 oC, n=1次 方 , 改 變 裂 縫 長 度 使
a / w =
0.01、0.03、0.05、0.07、0.1、0.15、0.2、0.3、0.4、0.5、0.6、0.7、0.8、0.9 其中
cm
w =
10 ,裂縫寬度a
為變數,以平面應 變分析,求得無因次 TSIF 之K
Ισ
Tπ w
與w
a /
關係繪於圖 3, 其中σT =EαcT (1−2υ)。 圖 3 顯示,在裂縫很小的時候,由於α
(x)為線性且隨
x
增加而增加,故在微小裂縫0
/ w ≈
a
時裂縫端點右邊會有很大面積其熱c n m
n
w x w
x x α α
α( )=( ) +(1−( ) ) C
T =1o T
E , υ , w
w 2 y
a x
) α( x
圖 2
α −
FGM 含邊裂縫受溫度載重示意圖膨脹係數比裂縫端點左邊大,溫度上升之 結果,會使裂縫尖端右側膨脹比裂縫尖端 左側大很多,在邊界條件的限制下,相對 於裂縫尖端左側會產生拉力效應,把裂縫 拉開,因此縱使裂縫長度很小,其 TSIF 不 會很小,且隨裂縫延伸而增大;在
a / w
為 0.18 時有最大值。當a / w ≥
0.18 時,隨著 裂縫長度增加,裂縫尖端右側的面積逐漸 變小,同時左側面積增大,使裂縫尖端右 側對左側的拉力效應減小,故應力強度因 子隨裂縫增長而減小;直到a/w ≥0.7時,因裂縫已經很大,裂縫尖端右側所產生的 拉力效應不足使裂縫張開, 故溫度上升裂 縫產生閉合現象,TSIF 開始進入負值。若 與張傳煜論文【5】中之單一均質材料受線 性溫度載重分析做比較,兩者之無因次整 體 TSIF 趨勢相同,因為 n=1 之
α
−FGM 與 均質材料受線性溫度載重,只影響兩者之 徹體力,由徹體力公式b
i*= − c
,iT − cT
,i計算,可知兩者之徹體力均為常數,故分析 得到相同之結果。
圖 3 線性
α −
FGM 受溫度載重之無因次 TSIF(B) 線性
α −
FGM,改變α
c及α
m的差值之分析
考慮圖 2 的問題,熱膨脹係數
α ( ) x
亦 為 線 性 函 數 , 即 n=1 ;α
c 固 定 為圖 4線性
α −
FGM 改變α
c與α
m差值之無因次 TSIFoC 10 5
0 .
1 × − ,而
α
m 則改變其值,分別為oC
m
5 1=1.2×10−
α , αm2=1.38×10−5 oC ,
oC
m
5 3=1.5×10−
α , αm4=1.6×10−5 oC ,
oC
m
5 5=1.7×10−
α ;將無因次 TSIF 與
a w
關 係繪於圖 4。由圖 4 可發現,當 n=1 之
α −
FGM 其( ) (
x x w)
αm(
x w)
αcα = + 1− ;當αm與αc差
距愈大時,經由徹體力公式可知,溫度載 重下其所產生的徹體力也愈大,故 TSIF 亦 愈大;在αc為定值下,αm愈大會有愈大的 TSIF,而在a w ≥ 0.7 時,其閉合壓應力 愈大,所以 TSIF 為較大之負值。注意的是,
不論αm如何改變,其 TSIF 為零的臨界點均 在a/w=0.7處,這表示在線性
α
−FGM 下,αc與αm差距之值並不影響裂縫閉合時之裂 縫長度,但會影響裂縫閉合壓應力大小。
(C) 固定
α
c及α
m的差值,改變n
值分析 圖 2 之問題,假設其熱膨脹係數函數 為 冪 次 方 , 分 析 時 將α
c 與α
m 固 定 為oC 10 5
0 .
1 × − 和 1.7×10 −5 oC ,而改變
n
為 1、2、3、4、15。將無因次 TSIF KΙ σT πw 與a w關係繪於圖 5。由圖 5 可發現,無因次 TSIF 與a w的 關係在不同之 n 值下曲線幾乎重合。對於 本節所探討之
α −
FGM,在x =0其α
c 為oC 10 5
0 .
1 × − ;在x = w其αm為 1.7×10−5 oC; 不論
n
之值為何,x =
0及x = w
兩端的α
c 與α
m都為固定,所差別只是物體沿x
方向的熱膨脹係數梯度分佈。另外,為了 解材料受力的情況,由徹體力公式bi* =−c,iT 計算此函數梯度材料之徹體力,其結果繪圖 5 單一冪次方
α −
FGM 改變 n 值之無因次 TSIF於圖 6。由圖 6 發現,縱使材料在
x =
0及w
x =
兩端的α
c 與α
m 為定值,在不同的n 次方下,其徹體力分佈明顯不同,而在
圖 5 的 TSIF 曲線則幾乎沒有變化。
為了能較明顯比較改變
n 值後 TSIF 的
微小差異,將n=1和n=15單獨繪於圖 7。由圖 7 發現,n=1和n=15時其 TSIF 為零的 裂縫長度會有微小移動,這和圖 4 的結果 有所不同,這表示不同的
n
次梯度函數,會微小地影響材料其 TSIF 為零的裂縫長度 位置。為了解
α
c和α
m的差值,和n
次方的 改變對材料 TSIF 的影響有何差別,我們再 分析一次相同的模型,並將 n 次方固定為 3 次,觀察在不同的α
c與α
m差值下之結果。0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Dimensionless Crack Length:a/w -0.5
-0.25 0 0.25 0.5
Dimensionless TSIF:K1*
n=1 n=2 n=3 n=5 n=4
) ( ,
, x
Eυα C T=1o
w
w 2 y
a x
c n m n
w x w
x x α α
α( )=( ) +(1−( ) )
C C
o m
o c
5 5
10 7 . 1
10 0 . 1
−
−
×
=
×
= α α
) 2 1
0 (
1
* 1
υ α σ
π σ
−
=
= T E
w K K
c T
T
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Dimensionless Crack Length:a/w -0.2
-0.1 0 0.1 0.2 0.3
Dimensionless TSIF:K1*
) ( ,
, x
Eυα C T=1o
w
w 2 y
a x
c n m n
w x w
x x α α
α( )=( ) +(1−( ) )
) 2 1
0 (
1
* 1
υ α
σ
π σ
−
=
= T E
w K K
c T
T
C C C C C C n
o m
o m
o m
o m
o m
o c
5 5
5 4
5 3
5 2
5 1
5
10 7 . 1
10 6 . 1
10 5 . 1
10 38 . 1
10 2 . 1
10 0 . 1 1
−
−
−
−
−
−
×
=
×
=
×
=
×
=
×
=
×
=
=
α α α α α α
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Dimensionless Crack Length:a/w -0.15
-0.1 -0.05 0 0.05 0.1 0.15 0.2
Dimensionless TSIF:K* 1
) ( ,
, x
Eυα C T=1o
w
w 2 y
a x
c n m n
w x w
x x α α
α( )=( ) +(1−( ) )
C C n
o m
o c
5 5
10 38 . 1
10 0 . 1 1
−
−
×
=
×
=
= α α
) 2 1
0 (
1
* 1
υ α σ
π σ
−
=
= T E
w K K
c T
T
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Dimensionless Location:x/w 0
2 4 6 8 10 12 14 16
Body-Force bi*(kg/cm3)
n=1 n=2 n=3 n=4
c n m n
w x w
x x α α
α( )=( ) +(1−( ) )
C C
o m
o c
5 5
10 7 . 1
10 0 . 1
−
−
×
=
×
= α α
圖 6 冪次方
α −
FGM 不同 n 值下之徹體力分佈令n=3,αc固定為1.0×10 −5 oC,分別 輸入αm1=1.2×10−5 oC, αm2 =1.38×10−5 oC,
oC
m
5 3 =1.5×10−
α , αm4 =1.6×10−5 oC ,
oC
m
5 5 =1.7×10−
α ; 將 所 得 無因次 TSIF
w
KΙ σT π 與a w關係繪於圖 8,比較發現
圖 7 冪次方
α −
FGM 在 n=1 與 n=15 之無因次 TSIF圖 8 一冪次方
α −
FGM 在 n=3 時改變α
c與α
m差值之無因次 TSIF
圖 8 和圖 4 很相似,TSIF 為零同樣交於一
點,
α
c與α
m的差值愈大其 TSIF 也愈大;並將對應於圖 8 的徹體力分佈繪於圖 9。
四、結果與討論
在本節將以徹體力的圖形解釋圖 5,在圖 6 中不同 n 次方下,徹體力曲線明顯不同,
徹體力單位為kg cm3,且為
x
之函數,若將 徹體力沿 x 軸積分,即求∫
0wbidx* ,所得即為 徹體力在 x 方向之合力。將圖 6 中不同徹 體力曲線沿
x
軸積分,發現結果其值均相 同,表示在不同的n
次方下,徹體力在x
方 向之合力都相同,這可以解釋為何在圖 5 中n
次方增加時,其 TSIF 仍保持不變,因 為徹體力作用之x
方向合力一樣。圖 8 中,n
固定為 3,α
c與α
m的差值增大時,圖 9中徹體力曲線下面積明顯增加, 亦即徹體 力 對 材 料 作 用 之 合 力 增 大 , 所 以
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 Dimensionless Length:a/w
0 1 2 3 4 5 6 7 8 9 10 11 12
Body-Force bi*(kg/cm3)
αm1
αm2
αm3
αm4
αm5
n=3
αm1=1.2×10-5/oC αm2=1.38×10-5/oC αm3=1.5×10-5/oC αm4=1.6×10-5/oC αm5=1.7×10-5/oC αc=1.0×10-5/oC
c n m n
w x w
x x α α
α( )=( ) +(1−( ))
TSIF 也因此增大。至於 TSIF 為零之點,
可以觀察在圖 6 中,
n
次方增加時其徹體 力作用之合力大小雖然不變,但徹體力曲 線下面積的重心位置會隨n
增大而向右 移,表示徹體力合力作用點之位置會隨著n
增大而在x
軸上向右移動,不同n
之下徹 體力合力之作用點在x
軸的位置可見表 1;徹體力之合力位置向右移,造成徹體力 對裂縫尖端左側的影響減小,使裂縫更易 於閉合,故 TSIF 為零的點會向左移動,所 以在圖 5 中當n
次方增加時,K
Ι= 0
之裂 縫長度減小。反觀圖9
,由徹體力合力位置圖 9 單一冪次方
α −
FGM 在 n=3 時改變 αc與αm差值之徹體力分佈 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Dimensionless Crack Length:a/w -0.5
-0.25 0 0.25 0.5
Dimensionless TSIF:K1*
n=1
n=15 )
( ,
, x
Eυα C T=1o
w
w 2 y
a x
c n m n
w x w
x x α α
α( )=( ) +(1−( ) )
C C
o m
o c
5 5
10 7 . 1
10 0 . 1
−
−
×
=
×
= α α
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Dimensionless Crack Length:a/w -0.2
0 0.2 0.4
Dimensionless TSIF:K1*
) ( ,
, x
Eυα C T=1o
w
w 2 y
a x
c n m n
w x w
x x α α
α( )=( ) +(1−( ) )
C C C C C C n
o m
o m
o m
o m
o m
o c
5 5
5 4
5 3
5 2
5 1
5
10 7 . 1
10 6 . 1
10 5 . 1
10 38 . 1
10 2 . 1
10 0 . 1 3
−
−
−
−
−
−
×
=
×
=
×
=
×
=
×
=
×
=
=
α α α α α α
之觀念可知,在
n
固定下改變α
c與α
m的差 值時,其徹體力圖形下面積的重心在x
軸 位置不會改變,故在圖8
中徹體力作用合 力於x
軸上不變,其TSIF
為零之點不變。由圖
4
,圖5
與圖8
整理得以下結論:( Ι )
對α − FGM
受定溫載重之問題,會直 接影響TSIF
的是徹體力作用之合力的大 小,當徹體力之合力愈大,TSIF
也愈大;而若徹體力之合力相同,但徹體力曲線下
x
軸之重心位置不同亦會影響TSIF
為零時 的裂縫長度;而只要徹體力合力在x
軸重 心位置相同,則TSIF
圖為零處會交於一點。( ΙΙ )
不論熱膨脹係數α (x )
如何分布,在7 .
≥0 w
a 時裂縫都會產生閉合的現象,而在
w
a / 靠近
0.7
時,TSIF
都趨近於零,並開 始進入負值,可見α − FGM
中其裂縫本身 的長短,會直接影響受溫度載重後之行為。五、參考文獻
【
1
】Miyamoto, Y., “Application of FGMs in Japan,” Ceramic Transactions, 76:
Functionally Graded Materials (ed. A.
Ghosh, Y. Miyamoto, I. Reimanis and J. J. Lannutti), pp. 171-189, American Ceramic Society, Westerville, Ohio (1997).
【
2
】宋見春,譚建國,蕭任智,涂聰明,「裂縫受徹體力作用之邊界元素法 解析
,
」中國土木水利工程學刊,Vol.
2, No. 2, June (1991)
。【
3
】Tang, W., Transformation Domain into Boundary Integrals in BEM , Spring-Verlag, New York (1988).
【
4
】Cherepanov, G. P., Mechanics of Brittle , McGraw-Hill, New York (1979).
【