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

(

sin )

)(

cos -(1

cos R n R n n R

R′= φ + φ ⋅ + φ × (2.1)

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

本文中採用增量迭代法解非線性平衡方程式，假設第I 個位置已知，此 處的第I 個位置，是指第I 個增量迭代收斂後的平衡位置。x_j 為元素節點 j (j =1,2,3)在當前元素座標上的位置向量，^IX_j、^Ix_ij^B、∆U_j、∆ 以及ε^B_j ∆φφφφj

分別為元素節點 j在固定總體座標中第I 個位置的位置向量、元素節點 j在 第I 個位置的節點基礎座標、元素節點 j的增量位移向量、元素節點 j的增 量應變向量以及元素節點 j的增量旋轉。本文中假設元素節點 j受∆U_j^、

B

εj

∆ 以及∆φφφφj作用後的變形過程如下：

1. 節點 j受∆U_j^{的作用由第}I 個位置^IX_j平移到節點當前變形位置X_j，在 移動過程中，節點 j無剛體旋轉，即元素各邊在節點的切線方向維持不 變，剛接在其上的節點基礎座標的方位亦不變。

2. 節點 j及剛接在其上的節點基礎座標軸^Ix^B受到_∆φφφφj作用，旋轉到當前變

形位置的x_ij^B。

3. 節點 j在當前變形位置的節點基礎座標上之應變分量增加一增量∆ε^B_j ^。

由上述的變形過程可知節點 j當前變形位置的固定總體座標X_j可以由^IX_j 加上∆U_j^得到，由X_j可以利用元素座標定義求出當前變形位置的元素座標

E

xi 。如圖 2.3 所示，將元素變形後當前變形位置的元素座標x_i^E與初始未變 形時的元素座標⁰x_i^E重疊，則元素節點變形參數可決定如下：

(1) 節點位移u_j

節點位移u_j可由元素節點 j在當前變形位置的元素座標上的位置向量x_j與 元素節點 j在初始未變形時的元素座標上位置向量的⁰x_j的差得出：

j j

j x x

u = −⁰ (2.2)

{

j j j

}

j = u v w

u (2.3)

{

^{j j j}

}

j = x y z

x (2.4)

{

^{j j j}

}

j ⁰x ⁰y ⁰z

0x = (2.5)

由元素座標定義的方式可知⁰z_j =z_j =w_j =0。 (2) 節點應變ε_j

利用第I 個位置在節點基礎座標上的節點應變向量^Iε^B_j 加上增量應變向量

B

εj

∆ 可以求出當前變形位置的節點基礎座標上的節點應變向量：

B j B

j I B

j ε ε

ε = +∆ (2.6)

{

xyj^B

}

B yj B

xj B

j ∆ε ∆ε ∆γ

∆ε = (2.7)

其中ε^B_j為當前變形位置的基礎節點座標上的節點應變向量。

(3) 節點變形角θ_ij

之x₁^B_j^′軸的單位向量。

令⁰θ_nj^E及θ 為將^E_nj ⁰θ_nj與θ 分別以初始未變形時的元素座標及當前變形位置_nj 的元素座標為基底，節點變形角θ_ij(i =x,y; j=1,2,3)^{可表示為：}

E nj E nj yj

xj

θ θ ⁰ 0

−

=











 θ θ

(2.11)

2.5 殼元素的變形殼元素的變形殼元素的變形描述殼元素的變形描述描述描述

如圖 2.6 所示之殼元素中心面上有三個節點，每個節點的自由度至少分

別是x₁^E、x₂^E、x₃^E軸方向的位移u_j、v_j、w ( _j j=1,2,3)，繞x₁^E、x₂^E、x₃^E軸方 向的位移轉角θ_xj^、θ_yj^、θ_zj^{( j=1,2,3)，以及應變分量}ε_xj^、ε_yj^、γ_xyj^。此殼元 素 的 變 形 可 分 為 薄 膜 變 形(membrane deformation)與 彎 曲 變 形(bending

deformation)兩部份，薄膜變形是來自平面元素，彎曲變形則是來自板元素，

由基本假設(1)可知此元素的變形可由薄膜變形與彎曲變形疊加而成。

在圖 2.6中的元素節點位移u 、_j v 及_j θ_zj ^{( j=1,2,3)是QST平面元素[26]}

節點位移，而θ_xj^、θ_yj^及w_j^{( j=1,2,3)為DKT 板元素[32]的節點位移。}

因本文使用共旋轉推導法，且假設在元素座標上，元素的變形位移及 旋轉為小位移及小旋轉，所以本章中除了推導幾何剛度矩陣外，都僅考慮 線性的位移－應變關係。

2.6 QST 平面平面平面元素平面元素元素的變形描述元素的變形描述的變形描述的變形描述

本文採用的三角平面元素為文獻[26]中所提出的 QST元素，其位移場為

三次變化，應變場則為二次變化。此元素有 3 個節點，每個節點有 6 個自 由度，此元素可使用彼此間能互相轉換的不同節點參數。本文中採用的節 點參數為節點 j在元素座標x₁^E、x₂^E軸的位移分量u_j、v_j，應變分量ε_xj^、εyj^、

γxyj^{及逆時鐘方向的旋轉}θ_zj。但為了方便推導，本文在元素推導時使用的節 點參數是節點 j在x₁^E、x₂^E軸的位移分量u_j、v_j以及x₁^E、x₂^E軸的位移分量 分別對面積座標(area coordinates)ξ ^、η^的微分u,_ξj、u,_ηj、v,_ξj、v,_ηj，面積 座標的介紹詳見附錄 A。

此元素的位移場可表示為[40]：

qξ

N^t_u

u= (2.12)

qξ

N^t_v

v= (2.13)

} { _{ξ1 ξ2 ξ3}

ξ q q q

q = (2.14)

0 0 0

{N₁ N₂ N₃

u = N

} 0 0 0

0 0

0 _{5 6 7 8 9}

4 N N N N N

N (2.15)

3 2

1 0 0

0

{ N N N

v = N

} 0

0 0

0 0

0 N₄ N₅ N₆ N₇ N₈ N₉ (2.16)

} , , , ,

{ _{j j j j j j}

j u v u_ξ u_η v_ξ v_η

ξ =

q (2.17)

其中^u =^u(ξ,η)、^v =^v(ξ,η)分別為在x₁^E、x₂^E軸的位移分量，N_u為對應於u 的 形 狀 函 數 ， N_v 為 對 應 於 v 的 形 狀 函 數 ， q_ξ 為 節 點 參 數 向 量 ，

Ni (i=1,2,3,L,9)為面積座標的三次函數，其表示式詳見附錄 B，另外本文 中

{ }

^{表行矩陣。}

由小變形的假設，元素內任意點的正應變、剪應變及逆時鐘方向的剛體旋

}

將(2.31)式與(2.34)式代入(2.12)式與(2.13)式可得：

θ

nj

nj θ

θ −⁰ <<1 時，DKT 元素的位移場可表示成：

u =zθ_y(x,y) v=−zθ_x(x,y) w= w(x,y) (2.41) 其中x、y、z 為元素上任一點分別在x₁^E、x₂^E、x₃^E軸的座標值，θ_y^是(θn−⁰θn⁾

在x₁^E軸方向的分量，θ_x^是(θ_n−^0θ_n⁾在x₂^E軸方向的分量，θ_n為當前變形位置 的元素變形角，⁰θ_n為初始未變形時的元素變形角，u是在x₁^E軸方向的位 移，v是在x₂^E軸方向的位移，w是在x₃^E軸方向位移。

DKT 元素的應變變形包含面內(in plane)正應變ε_x^、εy^與剪應變γxy^以

及橫向剪應變(transverse shear strain)γ_yz^、γxz^。

因本文假設元素的變形為小變形，ε_x^、εy^和γxy^{可表示成(2.18)至(2.20)}

式，γ_yz^、γxz^{可表示成：}

z x xz =w_, +u_,

γ γ_yz =^w_,y +^v_,z (2.42)

將(2.41)式代入(2.18)至(2.20)式可得：

κ

ε_b =z (2.43) ε_b ={ε_x ,ε_y ,γ_xy} (2.44)

x x y y y

x x

y, , , , , ,

{θ −θ θ −θ

=

κ (2.45)

將(2.41)式代入(2.42)式可得：

γ={γ_xz ,γ_yz}={w_,x +θ_y , w_,y −θ_x} (2.46)

由(2.41)式可知 w、θ_y^、θ_x^與x₃^E無關，所以可由(2.45)式知橫向剪應變在厚

度方向為常數。

本文中稱圖2.7中沿著元素邊緣方向s為切線方向，而垂直於元素邊緣

方向n為法線方向。

在文獻[32]中對於其所提出的 DKT 元素做了下列的假設：

(1) θ_y^、θx在元素內為二次變化，也就是：

∑

=

= ⁶

1 i

yi i

y Nθ

θ ^；

∑

=

−

= ⁶

1 i

xi i

x N θ

θ (2.47)

其中_θyi、θ_xi^是θ_y^、θ_x^{在圖 2.7 中節點i的節點值，}N_i⁽i =^1,L^,6)為形 狀函數，其表示式詳見附錄 C。

(2) 元素的三個頂點以及三個邊的中點滿足克希霍夫板理論(Kirchhoff plate

theory)的假設，即

(a) 在三個頂點

γ_xzi =w_,xi +θ_yi =0 i =1,2,3 (2.48a) γ_yzi =w_,yi −θ_xi =0 i =1,2,3 (2.48b)

其中w_xi、w_yi、θ_yi^、θ_xi^{分別是 ( )} x w_x w

∂

= ∂ 、 ( )

y w_y w

∂

= ∂ 、θ_y^、θ_x^在

節點i的值。

(b) 在三個邊的中點

, =0 +

−θ_nk w_sk k =4,5,6 (2.49a)

, =0 + _nk

sk w

θ k =4,5,6 (2.49b)

其中θ_nk^、θ_sk^分別是θ_n^、θ_s^{在節點k的值，}θ_n^與θ_s^{分別是θ在 n 與 s} 方向的分量，w_,sk、w_,nk分別是 _, ( )

s w_s w

∂

= ∂ 、 _, ( ) n w_n w

∂

= ∂ 在節點 k 的值。

(3) w 在元素邊緣的方向上是呈現三次變化，也就是：

sj

由(2.48)至(2.53)式可以把(2.47)式表示成[32]：

b

其中u_b為 DKT 元素的節點位移，H_x與H_y是對應於元素節點位移的新形狀 函數，其表示式詳見附錄 C，ξ^與η是元素內任一點在元素自然座標[]的座 標值，其中1≤ξ ≤0^、1≤η ≤^0。

將(2.54)式代入(2.45)式可以得到：

κ =B_bu_b (2.56)



將(2.37)式、(2.38)式代入(2.18)至(2.20)式可得：

θ

力，m_εxj、m_εyj、m_γj分別為對應於ε_xj^、ε_yj^、γ_xyj^{的廣義節點力矩，m}_zj^為對 應於θ_zj^{的節點傳統力矩，V} 為元素體積。由於應變的因次(dimension)與旋 轉相同，因此應變所對應的節點力m_εxj、m_εyj、m_γj應與力矩有相同的因次為 一廣義節點力矩，而非應力。將(2.63)式、(2.66)式代入(2.67)式可得

θ ξθ ξθ

θ θ

θ δ

δ^{qtf = qtTt}

∫

_V^{BtEBdVT q (2.70)}

由(2.70)式可得：

θ θ

θ k q

f = (2.71)

ξθ ξ ξθ

θ T k T

k = ^T (2.72)

=

∫

_VB^TEBdV

k_ξ (2.73)

其中k_ξ、k_θ 分別為應於_qξ、q_θ的 QST 元素剛度矩陣。

2.10 DKT 板元素之節點內力與剛度矩陣板元素之節點內力與剛度矩陣板元素之節點內力與剛度矩陣 板元素之節點內力與剛度矩陣 將(2.43)式、(2.56)式代入(2.58)式可得：

b b

b zEB u

σ = (2.74) 由虛功原理可得：

V bdV

t b b

t

b^{f =}

∫

^{ε σ}

u δ

δ (2.75) 其中f_b是板元素對應於δu_b^{的節點內力，V為元素的體積。}

將(2.43)式、(2.56)式、(2.74)式代入(2.75)式可得：

= ^t_b

∫∫

_A ^t_{b b b}

b t

bf u zB zEB u dzdA

u δ

δ ⁼

∫

_A b b b

t b t

b B D B u dA

δu (2.76) 其中

 獻[43]中以完整的 Green’s strain 及 total Lagrangian 推導法推導出一三維元素 之切線剛度矩陣的通式，本節中將採用文獻[43]中推導元素幾何剛度矩陣的



 



=

I I

I M I

y xy

xy x

σ τ

τ

σ I為 3×3 單位矩陣 (2.82)

b

m σ

σ

σ = + σ ={σ_x ,σ_y ,τ_xy} (2.83) 其中 j = 1,2,3為第 j 個節點，N_uθj、N_vθj為 (2.37)式之N_uθ、(2.38)式之N_vθ對 應於節點 j 的部份，H_xj、H_yj為(2.54)式之H_x、H_y對應於節點 j 的部份，I 為3×3的單位矩陣、0為 1×6的零矩陣，σ_m、σ_b為(2.66)式、(2.74)式之QST 平面元素及DKT板元素的應力場。

因為DKT元素的內部沒有定義側向位移場w，所以本文假設w為線性位移 場w= N_wu_b，其中N_w與CST元素[12]的形狀函數相同，本文除了使用QST 元素本身的位移場u=N_uθ,xq_θ、v=N_vθ,xq_θ外，還用u、v為線性位移場推導

(2.81)式之元素幾何剛度矩陣，本文推導的兩種元素幾何剛度矩陣都僅是近

似的幾何剛度矩陣，但由本文的數值結果發現這兩種元素幾何剛度矩陣在 非線性分析的平衡迭代時都能提高收斂速度且性相能近，在挫屈分析時也 都能偵測到相當接近的挫屈負荷。

2.12 座標系統轉換座標系統轉換座標系統轉換座標系統轉換

為了建立結構的平衡方程式，必須將元素節點參數q 中的u_j、v_j、w_j和

θxj^、θyj^、θzj轉換成固定總體座標的分量，ε_xj^、ε_yj^、γ_xyj^{轉換成對應於節點} 基礎座標的分量，才能將各元素節點內力組合成結構系統節點內力以及將 元素剛度矩陣疊加成結構系統剛度矩陣。在元素座標的節點位移u_j、v_j、w_j 和節點旋轉θ ^、θ ^、θ 與在固定總體座標的節點位移和節點旋轉之關係式

為：

} { _{j j j xj yj zj xj yj xyj}

Ej = u v w θ θ θ ε ε γ

q (2.92) 將(2.80)式、(2.81)式、(2.86)式代入(2.90)式可得：

B

由(2.71)和(2.78)可得：

b

由反梯度法則(contragradient law)[44]及(2.93)式可得：

E

式代入(2.101)式可得：

B B

B k q

f = (2.104)

EB E t EB

B T k T

k = (2.105)

其中k_B是對應於q_B的元素剛度矩陣。

2.13 系統平衡方程式與收斂準則系統平衡方程式與收斂準則系統平衡方程式與收斂準則系統平衡方程式與收斂準則

結構系統受外力負荷時，其平衡方程式可表示為

(

^Q

) ( )

^{=F Q − P = 0}

Ψ ^,λ λ ^(2.106)

其中Ψ為系統不平衡力向量，系統節點內力F可由(2.102)式的元素節點內力

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

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

etol

e= N ≤

P Ψ

λ ^(2.107)

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

第三章 第三章 第三章

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

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

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

3.1 增量迭代法增量迭代法增量迭代法增量迭代法

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

rT

Q= ∆λ

∆ (3.1)

( )

^{K P}

rT = T^{I −1} (3.2)

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

(3.1)式中的∆λ可利用下式求出[45]

( )

^{Tt T 1/2}

l r r

∆

±

=

∆λ (3.3)

其中正負符號決定方法為，當第I −1與I 個增量收歛時，其系統切線剛度 矩陣之行列式值同號，則∆λ的正負符號和第I 個增量時相同；若異號則符

其中正負符號決定方法為，當第I −1與I 個增量收歛時，其系統切線剛度 矩陣之行列式值同號，則∆λ的正負符號和第I 個增量時相同；若異號則符

在文檔中 平面三角形殼元素之改善研究 (頁 22-0)