殼元素的變形描述

如圖2.1 所示，殼元素中心面上有三個節點，每個節點的自由度分別是

、 、 軸方向的位移 、 、 ( j=1,2,3)，繞 、 、 軸方 向的位移轉角

x₁E x₂^E x₃^E u_j v_j w_j x₁^E x₂^E x₃^E

xj、_yj、_zj (j 1,2,3)，應變分量_xj、_yj、_xyj (j 1,2,3) 以及側向位移的二次微分w,_xxj、w,_xyj、w,_yyj (j 1,2,3)。

在圖 2.1 中的殼元素節點自由度u_j、v_j及_zj (j1,2,3)是 QST 平面元 素[4]節點自由度，而w_j、_xj、_yj (j 1,2,3)為 RQT 板元素[16-18]的節點

自由度。本文以下敘述的元素變形、元素內力及元素剛度都是定義在元素 座標下。

2.3.1 QST 平面元素的變形描述

本文採用的三角平面元素其位移場為三次變化，應變場則為二次變化。

此元素有 3 個節點，每個節點有 6 個自由度，此元素可使用彼此間能互相 轉換的不同節點參數。本文中採用的節點參數為節點 j在元素座標 、 軸 的位移分量 、 ，應變分量

x₁E x₂^E

uj v_j _xj、_yj、_xyj及逆時鐘方向的旋轉_zj。但為 了方便推導，本文在元素推導時使用的節點參數是節點 j在 、 軸的位 移 分 量 、 以 及 、 軸 的 位 移 分 量 分 別 對 面 積 座 標(area

x₁E x₂^E

uj v_j x₁^E x₂^E

coordinates)、的微分u,_j、u _j、 、 ，面積座標的介紹詳見附錄B。

)

(2.21)

將(2.20)式與(2.23)式代入(2.1)式與(2.2)式可得：

(2.26)



向 微 分w,_n為 三 次 函 數 的 三 個 限 制 條 件 ， 可 以 把₁₉ ~₂₁各 自 表 示 成



將(2.39)式代入(2.31)式推得：

b

將(2.35)、(2.36)式是代入(2.46)至(2.48)式推得：

2

)

} 其中E 是楊氏模數(Young’s module)， 是蒲松比(Poisson ratio)。

2.4.1 QST 平面元素之節點內力與剛度矩陣 將(2.26)式、(2.27)式代入(2.7)至(2.8)式可得：

(2.59)

將(2.59)式代入(2.54)式可得：

(2.62)

(2.65)

將(2.59)式、(2.62)式代入(2.63)式可得

(2.66) 將(2.52)式代入(2.54)式可得：

b

其中 為對應於 的元素節點內力， 為對應於 的節點力， 、 為

將(2.52)、(2.70)式代入(2.71)式推得：

dV

中以完整的Green’s strain 及 total Lagrangian 推導法推導出三維元素之切線 剛度矩陣的通式。文獻[6]修改、簡化文獻[26]推導元素幾何剛度矩陣的方 法，使其適用於推導殼元素的幾何剛度。本文採用文獻[6]推導殼元素幾何 剛度的方法，元素幾何剛度矩陣可表示式成

V dV

移和節點旋轉之關係式為：

(2.93)

殼元素在元素座標下的節點內力可以表示成： 由反梯度法則(contragradient law)及(2.96)式可得：

(2.103)

2.6 元素節點變形參數的決定方法

由上述的變形過程可知節點 當前變形位置的固定總體座標 可以由

(3) 側向位移二次微分κ_j

則⁰x₃^E軸和⁰x₃^B_j軸平行，所以⁰θ_nj 0。

fB疊加得出，Q 為系統位移向量，為負荷參數，P為參考外力負荷向量。

若外力為與變形位置相關(configuration dependent)的外力，則 在每一個變 形位置都須更新。本文以不平衡力 的weighted Euclidean norm 做為平衡迭 代時的誤差度量，且收歛準則表示為

P Ψ

etol

e ^ΨN 

 P (2.122)

其中N 代表離散系統的自由度總數，e_tol為一給定的容許誤差值。

第三章 數值計算方法與程序

本文解(2.121)式的非線性平衡方程式所使用平衡迭代的數值計算方法 是採用文獻[27]中所提出基於牛頓-拉福森(Newton-Raphson)法和弧長控制 (arc length control)法的增量迭代法。本文中為了求得分歧點，以系統切線剛 度矩陣之行列式值為零來判斷。本文採用二分法決定增量位移向量的長 度，以求得系統切線剛度矩陣之行列式值為零的平衡位置。為了求得次要 平衡路徑，本文在平衡路徑的第一個挫屈負荷分歧點加入一個與第一挫屈 模態向量方向一致的擾動位移。為了本文的完整性，以下將簡單介紹文獻[27]

中提出的增量迭代數值計算方法與程序。

3.1 增量迭代法

若第 個增量的平衡位置為已知，令其位移向量為I Q_I、負荷參數為_I， 則第I 1個增量的初始增量位移向量 ，可以利用尤拉預測值(Euler Q predictor)求得[28]

rT

Q 

 (3.1)

(3.2)

 

^{K P}

r_T  _T^{I 1}

其中 為初始增量負荷參數， 為第 個平衡位置的系統切線剛度矩陣。

(3.1)式中的

I

KT I

 可利用下式求出[27] 

 

T ^1/2 t

l rTr







 (3.3)

其中正負符號決定方法為，當第I 1與 個增量收歛時，其系統切線剛度 矩陣之行列式值同號，則

I

 的正負符號和第 個增量時相同；若異號則符 號相反。

I

 表示第l I 1個增量的增量弧長，其值可以如下決定[27]



_{D I}



¹²

I J J

l l 

 (3.4)

其中J_D為給定的期望迭代次數，J_I為第 個增量迭代至平衡所使用的迭代I

次數， 為第 個增量的增量弧長。本文第一個增量的增量弧長 是由下 作用下的系統線性解之Euclidean norm， 為給定之最大增量次數，

P I_max r_c

其中r  Q為增量位移修正量，整理(3.10)式，可得到位移修正量

(3.11)



_I



_T

T Ψ P r r

K

r  ^1 _1  ₀ 

其中r₀ K_T^1Ψ_I1為標準牛頓法的位移修正量，r_T在(3.2)式中已定義。

3.2 弧長控制法

(3.11)式中的負荷參數修正量可利用文獻[27]中所提出的定弧長控制

法決定，其方法在每一個增量中固定其增量位移向量的Euclidean norm為一 定值Δl，由新的增量位移向量



ΔQr



可以得到

(3.12)



^Q^r

 

^Q^r

 Δ Δ

Δl^{2 t}





將(3.11)式代入(3.12)式可以得到

(3.13)



_T

 

^t _T

l²  ΔQr₀ r ΔQr₀ r Δ

(3.13)式經過整理後可以得到的二次方程式

0 (3.14)

3 2 2

1 a a 

a 

其中

(3.15)

T t

a₁rTr

(3.16)

 

^t _T

a₂ 2 ΔQr₀ r

(3.17)



₀

 

₀



²

3 Δ Δ Δl

a  Qr ^t Qr 

將(3.11)式求得之增量位移修正量 以及(3.14)式求得之負荷參數修正量r 

加入上次迭代之ΔQ與中，可以得到新的增量位移向量與增量負荷參 數，再進行下一次迭代，迭代過程將一直重複至(3.6)式中的系統不平衡力 向量滿足(2.122)式的收斂準則為止。

3.3 二分法

利用 3.1 節的增量迭代法可以求得結構之主要平衡路徑。在每個增量的 迭代收歛時，可以得到該增量在其平衡位置的負荷參數及結構剛度矩陣的

行列式值D()。令_I及D(_I)分別表示第 I 個增量在其平衡位置的及 )

(

D 值。_I1及D(_I1)分別表示第I 1個增量在其平衡位置的及D() 值。l_I1表示第I 1個增量的增量位移向量之弧長。若D(_I)大於零且

1) ( _I

D  小於零則可利用以下二分法求得挫屈負荷參數_cr：

(1) 令l_L 0，l_R l_I1，_L _I，_R _I1。其中下標 L 及 R 表示左 界及右界。

(2) 取

2



1  l

 L

lI

 l_R

， 重 做 第 I 1個 增 量 迭 代 ， 並 求 得 新 的 _I1 及 )

( _I1 D  。

(3) 若D(_I1)大於零，則令_L _I1，l_L l_I1。若D(_I1)小於零，則 令_R _I1，l_R l_I1。

(4) 若下列二式挫屈誤差準則同時滿足，則_I1即為系統挫屈負荷，否則回 到步驟(2)，重新展開下一次二分法迭代。

^I e ) 0 (

1)

 D D(

D

_





 e

I L

R  

1

其中 及 為給定的容許誤差值，本文例題之計算給定 及

。 eD

5

e e_D 10^9

10 e

3.4 數值程序

本文使用的增量迭代法之數值程序可以分為三個主要部分：

1. 輸入並計算開始分析所需的資料

(a)輸入結構資料、邊界條件與負荷參數。

(b)選擇一個參考自由度，並給定期望此自由度應達到的位移。

(c)給定最大增量次數、最大迭代次數、期望迭代次數與容許誤差值。

(d)形成系統剛度矩陣並求得(3.5)式中的R₀。

(e)利用(3.1)式、(3.3)式與(3.5)式計算初始增量位移向量、初始增量負荷 參數與第一次增量弧長。

2. 使用迭代法求增量的收歛解

(a)由 2.6 節的方法，利用已知的增量位移求得當前元素的變形向量 和

，並計算(2.100)式中元素節點內力 ，接著由2.5 節座標轉換將 利 用(2.103)式轉換至 ，然後將元素節點內力 組合成結構系統節點內力

。

q

fE

ub

F

fE

fB f_B

(b)計算(2.121)式的不平衡力 Ψ 。

(c)檢查(2.122)式的收斂準則，若滿足則進入第 3 部分；否則檢查迭代數，

如果小於給定之最大迭代次數，則進行步驟(d)；否則減少增量弧長並以 (3.1)式與(3.3)式計算新的增量位移向量與增量負荷參數，回到步驟(a)重 新計算。

(d)利用(3.11)式與(3.14)式計算增量位移修正量與增量負荷參數修正量，

回到步驟(a)重新計算。

3. 計算下一次增量所需的資料

(a)檢查參考自由度的位移及進行的增量次數是否已達給定值，若已達到 給定值則停止分析工作；否則進行步驟(b)。

(b)計算(3.11)式中的切線剛度矩陣 。本文中 的計算方法是利用當 前變形位置的元素座標重新計算元素剛度矩陣

KT K_T

k ，再加上元素幾何剛度 矩陣 ，本文中忽略元素其他的幾何變形與內部應力對元素剛度矩陣造 成的影響，再利用 2.52 節座標轉換將 利用(2.107)式轉換至 ，然後 將元素剛度矩陣 疊加成結構系統剛度矩陣 。

k

kE k_B

kB K_T

(c)利用(3.1)式、(3.3)式與(3.4)式計算下一次增量的增量位移向量、增量 負荷參數與增量弧長。

(d)回到第 2 部分進行迭代工作。

第四章 數值分析與結果

本章節為了測試本文 RQT+QST 殼元素的性能，先使用一個線性例題測試 RQT 元素的準確性與收斂性，而 QST 元素的線性分析在文獻[11]有詳細的 測試，故本文不再多加描述 QST 元素的性能。本章節亦使用多個幾何非線 性例題和挫屈例題測試本文 RQT+QST 元素的性能，並且與文獻

[6]DKT+QST 元素和文獻[4]、[24]DKT+OPT 元素的結果相互比較。本章節 最後測試了文獻[13]中聚酯圓柱薄殼受兩階段負荷作用的例題，並且跟文獻 [13]的實驗結果和文獻[5]、[13]的數值分析的結果相互比較。

4.1 RQT 元素線性分析

本例題考慮一個平面方形薄板分別受到中點集中力 P 與均勻分布力 作 用。文獻[16]中曾用此 RQT 分析此例題，本例題採用與文獻[16]相同形式的 網格及邊界條件來驗證本文程式的正確性及確認 RQT 元素的準確性。圖 4.1(a)為本例題的結構示意圖，圖 4.1(b)為本例題的網格示意圖，根據不同 網格形式，網格可以分為 Q type 和 T type。本例題考慮了兩種邊界條件：

簡支邊界(Simply supported boundary)與固定邊界(Clamped boundary)。因對 稱之故，與文獻[16]一樣本文分析時僅考慮左下角

q0

4

1 薄板，分析時的邊界條

件如下

簡支邊界的邊界條件：

0

w X 0,Y 0 0

,_x

w Y 0,X  L 0

,_y

w X 0,Y  L 0

,_xx 

w Y 0,X  (不包含L (L,L)) 0

,_xy 

w X  ,L Y  L 0

,_yy 

w X 0,Y  (不包含L (L,L)) 固定邊界的邊界條件：

0

w X 0,Y 0

0 ,_x

w X 0,Y 0,X  L 0

,_y

w X 0,Y 0,Y  L 0

,_xx 

w X 0,Y 0,X  (不包含L (L,L)) 0

,_xy 

w X 0,Y 0,X  ,L Y  L 0

,_yy 

w X 0,Y 0,Y  (不包含L (L,L))

本例題分析的結果文獻[16]的結果完全一樣，所以本文的程式應是正確的，

另外本文還用文獻[6]的 DKT 元素分析了本例題。表 4.1 至 4.4 為本文的結

果及文獻[16]、[22]的解析解。表 4.1 至 4.1 中

) 1 ( 12 ²

3



 Eh

D ，W_C為中點C的 側向位移，表中的百分誤差是相對於解析解的誤差。表 4.1、4.2 分別是簡 支板受到中點集中力及分佈力的無因次化C點側向位移，表 4.3、4.4 分別

4.4 可以看出在本例題，不管是 RQT 元素或是 DKT 元素，使用 Q type 的網 格在收斂性上是優於 P type 的網格，而 RQT 元素的收斂速度也是遠大於 DKT 元素。

4.2 RQT+QST 殼元素的幾何非線性與挫屈分析 例題一：圓柱殼片段受集中力作用

本例題為一個圓柱殼片段受到單點集中力作用。圖 4.2(a)及圖 4.1(b)為圓柱 殼片段的幾何形狀及受到集中力之示意圖。圓柱殼片段 AD 和BC兩邊為鉸 接， AB 和DC兩邊為自由端，結構 E 點受到集中力負荷，其邊界條件設定 為：AD和BC兩邊U V W _x _y 0。本例題使用的網格為10 ，10 容許誤差值取 。圖4.3是本例題的平衡路徑圖，其中包含主要路徑及次 要路徑。在主要路徑分析過程中使用了 個增量，每個增量平均迭代次數 為 4 次。在次要路徑分析過程中使用了 個增量，每個增量平均迭代次數 約為

4

10

27 21 5

~

4 次。在主要路徑上偵測到的挫屈負荷F_cr 530.803，此時 E 點向 下位移9 mm。由圖4.3 可以得知，本文的結果與文獻[23]相當吻合，故 本文的數值程序可以準確的找出主要平衡路徑、分歧點及次要路徑。

783 .

例題二：槽型梁受扭矩(Channel section in torsion)

圖4.4(a)、4.4(b)為槽型梁受力和網格示意圖，結構線段 AB、BC、CD和 EF 、

FG、GI 為鉸接，結構J、K 點受到集中力負荷，線段 AB、BC、CD和 EF 、 FG、GI 受到均勻力負荷，其邊界條件為：在 M 點的U 0， B 、 點的F

0

V ，線段BC及FG的W 0。本例題使用(121)20和(242)40兩 種網格，取用的容許誤差為 。表4.1為本文和文獻[4]、[6] 使用相同網

V ，線段BC及FG的W 0。本例題使用(121)20和(242)40兩 種網格，取用的容許誤差為 。表4.1為本文和文獻[4]、[6] 使用相同網

在文檔中 高階平面三角形殼元素之研究 (頁 21-0)