未來研究方向

本文使用線性位移場推導幾何剛度，因此本文的切線剛度皆為近似剛度，

以後的研究可以使用元素內部真實的位移場推導幾何剛度，求得準確的切 線剛度，並且探討使用不同的切線剛度進行幾何非線性分析，會對迭代的 過程和偵測分歧點、極限點的結果造成什麼影響。

參 考 文 獻

[1] K. J. Bathe, L. W. Ho, “A simple and effective element for analysis of general shell structures”, Computer and Structure, vol. 13, pp. 673-681, 1981.

[2] K. M. Hsiao, “Nonlinear analysis of general shell structures by flat triangular shell element”, Computer and Structure, vol. 25, pp. 665-674, 1987.

[3] 楊禮龍, “薄殼結構在位移負荷作用之下的幾何非線性分析＂, 交通大 學機械工程學系碩士論文, 台灣, 新竹, 2006.

[4] J. M. Battini, C. Pacoste, “On the choice of the linear element for corotational triangular shells”, Computer Methods in Applied Mechanics and Engineering, vol. 192, pp. 2125-2168, 2003.

[5] 楊水勝, “拘限彈性薄板受側向位移負荷的行為研究”, 交通大學機械工 程學系碩士論文, 台灣, 新竹, 2007.

[6] 林寬政, “平面三角形殼元素之改善研究＂, 交通大學機械工程學系碩 士論文, 台灣, 新竹, 2010.

[7] P. Khosravi, R. Ganesan, R. Sedaghati, “Corotational non-linear analysis of thin plates and shells using a new shell element”, International Journal for Numerical Methods in Engineering, vol. 69, pp. 859-885, 2007.

[8] D. J. Allman, “A compatible triangular element including vertex rotations for plane elasticity analysis”, Computers and Structures, vol. 19, pp. 1-8, 1984.

[9] L. Damkilde, M. Gronne, “An improved triangular element with drilling rotations”, Proceedings of the 15th Nordic Seminar on Computational Mechanics, pp. 135-138, 2002.

[10] C. A. Felippa, “A study of optimal membrane triangles with drilling freedoms”, Computer Methods in Applied Mechanics and Engineering , vol.

192, pp. 2125-2168, 2003.

[11] 林育丞, “具旋轉自由度之三角平面元素的共旋轉推導法＂, 交通大學 機械工程學系碩士論文, 台灣, 新竹, 2008.

[12] J. L. Batoz, K. J. Bathe, L. W. Ho, “A study of three-node triangular plate bending elements”, International Journal for Numerical Methods in Engineer, vol. 15, pp. 1771-1812, 1980.

[13] A. Boudaoud, P. Patrício, Y. Couder, M. B. Amar, “Dynamics of singularities in a constrained elastic plate”, Nature, vol. 407, pp. 718–720, 2000.

[14] P. Patrício, W. Krauth, “Numerical solutions of the Von Karman equations for a thin plate”, Journal of Modern Physics, vol. 8, pp. 427-434, 1997.

[15] J. H. Argyris, I. Fried, D. W. Scharpf, “The tuba family of plate elements for the matrix displacement method”, The Aeronautical Journal of the Aeronautical Society, Vol. 72, pp. 701-709, 1968.

[16] G. R. Cowper, E. Kosko, G. M. Lindberg, M. D. Olson, “Static and dynamic applications of a high-precision triangular plate bending element”, AIAA Journal, Vol. 7, pp. 1957-1965, 1969.

[17] K. Bell, “A refined triangular plate bending finite element”, International Journal for Numerical Methods in Engineering, Vol. 1, pp. 101-122, 1969.

[18] Suman Dasgupta, Dipak Segupta,“A higher-order triangular plate bending element revisited”, International Journal for Numerical Methods in Engineering, Vol. 30, pp. 419-430, 1990.R. W. Clough, J. L. Tocher,

“Finite element stiffness matrices for analysis of plate bending”, Proc. Conf.

on Matrix Methods in Structural Mechanics, WPAFB, Ohio, pp. 515-545, 1965.

[19] K. T. Joseph, K. Singa Rao, “A fast algorithm for triangular plate bending element”, International Journal for Numerical Methods in Engineering, vol.

14, pp. 1100-1104, 1979.

[20] S. C. Jardin, “A triangular finite element with first-derivative continuity applied to fusion MHD applications”, Journal of Computational Physics, Vol. 200, pp. 133-152, 2004.

[21] C. Militello, C. A. Felippa, “The first ANDES elements:9-DOF plate bending triangles”, Computer Methods in Applied Mechanics and Engineering, vol. 93, pp. 217-246, 1991.

[22] Rober L. Taylor, Sanjay Govindjee, “Solution of clamped rectangular plate problems”, Technical Report:UCB/SEMM, 2002.

[23] E. H. Boutyour, H. Zahrouni, M. Potier-Ferry and M. Boudi, “Bifurcation points and bifurcated branches by an asymptotic numerical method and Pade approximants”, International Journal for Numerical Methods in Engineering, vol. 60, pp. 1987-2012, 2004.

[24] J. M. Battini, “A modified corotational framework for triangular shell elements”, Computer Methods in Applied Mechanics and Engineering, vol.

196, pp. 1905-1914, 2007.

[25] J. M. Battini, C. Pacoste, “On the choice of local element frame for corotational triangular shell elements”, Communication in Numerical Methods in Engineering, vol. 20, pp. 819-825, 2004.

[26] O. C. Zienkiewicz, The Finite Element Method in Engineering Science, CENTRAL BOOK CO., 1971.

[27] M. A. Crisfield, “A fast incremental/iterative solution procedure that handles 'snap-through'”, Computers and Structures, vol. 13, pp. 55-62, 1981.

[28] M. A. Crisfield, Non-linear Finite Element Analysis of Solids and Structures, John Wiley & Sons, 1991.

表 4.1 簡支板受中點集中力 P 之無因次化中點側向位移

1.13635 (-2.046%)

1.24840 (7.613%)

1.40314 (20.952%) N=4 1.15742

(-0.229%)

1.15333 (-0.583%)

1.16880 (0.751%)

1.28389 (10.672%) N=6 1.15890

(-0.102%)

1.15706 (-0.261%)

1.16470 (0.398%)

1.22499 (5.595%) N=8 1.15942

(-0.057%)

1.15838 (-0.147%)

1.16340 (0.286%)

1.19995 (3.436%)

N=2 4.0609374 (-0.03484%)

4.06848 (0.15095%)

- (31.73%)

- (4.55%) N=4 4.0623473

(-0.00013%)

4.06273 (0.00920%)

- (4.49%)

- (5.37%) N=6 4.0623517

(-0.00002%)

4.06239 (0.00091%)

- -

- - N=8 4.0623524

(-0.00001%)

4.06236 (0.00018%)

- (1.01%)

- (1.56%)

4.06235

表4.3 固定板受中點集中力P之無因次化中點側向位移

5.18580 (-7.479%)

5.66900 (1.142%)

6.21891 (10.953%) N=4 5.58233

(-0.404%)

5.54273 (-1.111%)

5.85600 (4.478%)

6.35977 (13.466%) N=6 5.59990

(-0.091%)

5.58137 (-0.422%)

5.70800 (1.838%)

6.07079 (8.310%) N=8 5.60530

(0.005%)

5.59485 (-0.181%)

5.64100 (0.642%)

5.91075 (5.455%)

表4.5 槽型梁受扭矩的極限點

Mesh Present DKT+QST[6] DKT+OPT[4]

6

Mesh Present DKT+QST[6] DKT+OPT[4]

14

Mesh Present DKT+QST[6] DKT+OPT[4]

20

X₃G

X₂G

X₁G

x₁E

x₃E

x₂E

h

v3

u3

w3 3



z

3



y

3

2

1

3

x

圖 2.1 三角元素的示意圖及節點自由

n

1

n

2

n

3

t

1

t

2

t

3

1

2 y 3

x

圖 2.2 邊上切線和法線方向示意圖

v , x , x₂^E ₂^E

0

u3 3

v3

03 Current

Initial

1

01,

x₁₂B^ 0

x₁₂B^

x₁₂B^ 2 0



z

02

u2 2 ⁰x₁^E,x₁^E,u

圖 2.3 元素變形示意圖

θ

nj



B

x

₁j

x

B

x

₃j

B j E

x x

₃

,

₃^

B j 2

θ

nj

j

圖 2.4 元素節點 j中心面之x_ij^B軸受旋轉向量θ_nj作用的示意圖

B j

B

x

₁j^ B

x

₁j^ 0

x

₂^

x

₂^B_j

zj 0





zj

E B

j B

j

x x x

_{3 3 3}

0 ^

,

^

,

j

^

圖 2.5 元素節點 中心面之j _{zj 3}e^E為將⁰x₁^B_j^軸旋轉到x₁^B_j^軸旋轉向量的示意圖

n

R

R



圖 2.6 旋轉向量

(a)

O F

C G H

X₁G

X₂G

L L

L L

(b)

O

C

F H G

Q type N 4

O

C

F H G

P type N 4

圖 4.1 方形薄板 (a)結構示意圖 (b)網格種類及切割示意圖

(a)

254mm L

2540mm R 

6.35mm h 

0.1rad

 

3102.75MPa E

0.3

 L

L

R

h

 

E D

A

C

B

(b) (c)

` F

R R

 

X₁G

X₃G

U W

E

A B

C D

圖 4.2 圓柱殼片段受到單點集中力作用 (a)結構尺寸示意圖 (b)力負荷圖 (c)網格 10×10 示意圖

0 5 10 15 20 25 30 -400

-200 0 200 400 600 800

Load F(N )

Vertical displacement at E(mm)

present ref.[23]

ref.[6]

圖 4.3 圓柱殼片段受集中力作用的位移-負荷曲線圖

(a)

(b)

圖 4.4 槽型梁受扭矩之 (a)結構尺寸示意圖 (b)網格(1+2+1)×2 示意圖

1100

L E 2.110⁶

75

b  0.3

25

h P125

3 .

1 t

  _I

E

K ^h

J P

G

C B

U X₁^G, V

X₂^G,

W X₃^G,

D t _F

A b

M

L

P

A B

C D

E F

G I

0 5 10 15 20 25 30 35 40 45 50 0

150 300 450 600 750 900 1050 1200

L oad in g p arameter



Displacement -W

_J

present (1+2+1)x20 ref. [4] (1+2+1)x20 present (2+4+2)x40

圖4.5槽型梁受扭矩的位移-負荷曲線圖

(a) L100，b50，t 2，E 100， 0.3

E (b)

圖4.6 受均勻壓縮位移之簡支板的 (a)結構尺寸示意圖 (b)網格 4×6示意圖

A B

C D

U X₁^G,

V X₂^G,

L

A

B C D

 

b 4

L

2 b

E

-2 0 2 4 6 8 10 12 0

50 100 150 200 250 300 350

Reaction Force

Displacement

U_AD() W_E I

J

K

圖4.7 受均勻壓縮位移負荷之簡支板的位移-負荷曲線圖

(a)

(a)

t

W X₃^G,

V X₂^G,

U X₁^G,

 I

G

F E

D

C B

A

L

b h

150

L E 21000

10

b  0.3

10 h

1 t

(b)

A B

C D

E G F

I

圖4.9 槽型梁 (a)結構尺寸示意圖

(a)

450

L E 70960

38

b  0.321

65 h

1 t

(b)

圖4.10T形斷面梁 (a)結構尺寸示意圖 (b)網格(2+2+3)×4示意圖

A C B

E G F

D I

C G P

W X₃^G,

V X₂^G,

U X₁^G,

t I

G F

E

D

C B

A

b

J

h L

-0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 0.0

0.5 1.0 1.5 2.0 2.5 3.0 3.5

Load P(10

3

)

Displacement

W_J U_J V_J

圖 4.11 T形斷面梁受側力之位移-負荷曲線圖

(a)

L

U X₁^G, V

X₂^G,

W X₃^G,

A

B

C

D

F

P

t

E

b

h

1400

L E 193300

75 .

47

b  0.3

75 .

72 h

5 .

2 t

(b)

A

B

C

D E

F

圖4.12角型梁 (a)結構尺寸示意圖 (b)網格(2+3)×2 示意圖

0 20 40 60 80 100 0

1 2 3 4 5 6 7 8 9 10

Load P

Displacement -U

_A

RQT+QST (2+2) x 25 RQT+QST (4+6) x 60 DKT+ALL[4] (4+6) x 60 DKT+OPT[4] (4+6) x 60 DKT+QST[6] (2+2) x 25 DKT+QST[6] (4+6) x 60

圖 4.13 角型梁受軸向集中力之位移-負荷曲線圖

(a) R1.016 E2.068510⁷ 048

.

3

L  0.3 03

.

0 t

(b)

圖4.14 懸臂圓柱殼 (a)結構尺寸示意圖

(b)網格 16×16示意圖 R

F 2

U X₁^G, V

X₂^G,

W X₃^G, A

t

end L Free

end Clamped B

C D

F 2

C D

A

B

0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 0

100 200 300 400 500 600 700 800 900

Load F

Displacement W

_A

/R

RQT+QST DKT+OPT[24]

DKT+QST[6]

圖4.15 懸臂圓柱殼受一對端點集中力之位移負-荷曲線圖

(a)

B

x₃j x₂^Bj

B

x₁j

R

W X₃^G,

V X₂^G, U

X₁^G,

F

B A 2

t

F F 2

2

F 2

C

10 R

5 .

0 t

104

 E

2 .

0



(b)

圖4.16半球殼 (a)結構尺寸示意圖 (b)網格 12×12示意圖

B A

C

0 1 2 3 4 5 6 7 8 9 10 0

5 10 15 20 25 30 35 40 45

V_A DKT+QST[6]

L oad in g Parameter



Displacement

U_B RQT+QST V_A RQT+QST DKT+OPT[24]

U_B DKT+QST[6]

圖4.17 半球殼受到二對集中力之位移-負荷曲線圖

(a)

 

d

L W

h

35

L

cm

5 .

17

W

cm

5 .

16

d

cm

35 .

0

h

mm

20

 

2 9Nm 10 8 .

3  ^

 E

4 .

0

 (b)

E D

A B

C G

F

I

H X

Y

(c)

Z

圖4.18 圓柱薄殼 (a)結構示意圖 (b)俯視圖 (c)前視圖

X

X Z

Z

 

F H

X

0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0 5

10 15 20 25 30

Displacement W

E

(mm)

Loading parameter



(rad)

present ref.[5]

圖 4.19 圓柱薄殼第一階段受力E點之負荷參數-位移曲線圖

0 5 10 15 20 0

5 10 15

20

ref.[13] experiment ref.[13] numerical ref.[5] 30x60

L oad in g p arameter



Displacement -W

_E

(mm)

present 30x60

圖4.20 圓柱薄殼第二階段受力E點之位移-負荷曲線圖

O

0 5 10 15 20

0 5 10 15 20

secendary

L o adin g p arameter



Displacement -W

_E

primary

I J

K

M N

圖4.21 圓柱薄殼第二階段受力E點位移-負荷曲線之主要路徑與次要路徑

(a)

50 100 150

0 50 100 150 200 250 300 350

-2 0 2 4 6 8 10 12 14 16 18 20 22 24 26

(b)

50 100 150

0 50 100 150 200 250 300 350

-2 0 2 4 6 8 10 12 14 16 18 20 22 24 26

圖4.22 (a)圖4.21中O點(W_E 0  0) (b)圖4.21中I點 (W_E 7.5309 6.5386) 側向位移場的等高線圖

(a)

50 100 150

0 50 100 150 200 250 300 350

-2 0 2 4 6 8 10 12 14 16 18 20 22 24 26

(b)

50 100 150

0 50 100 150 200 250 300 350

-2 0 2 4 6 8 10 12 14 16 18 20 22 24 26

圖4.23 (a)圖4.21中J點(W_E 9.1726  9.1452) (b)圖4.21中K 點 (W_E 11 5455.  14.5136) 側向位移場的等高線圖

(a)

50 100 150

0 50 100 150 200 250 300 350

-2 0 2 4 6 8 10 12 14 16 18 20 22 24 26

(b)

50 100 150

0 50 100 150 200 250 300 350

-2 0 2 4 6 8 10 12 14 16 18 20 22 24 26

圖 4.24 (a)圖4.21中M 點(W_E 9.4736  9.4156) (b)圖 4.21中N點 (W_E 14 5441.  9.14933) 側向位移場的等高線圖

附錄 A

文獻[13]的實驗是探討一聚酯圓柱薄殼受位移負荷作用後的非線性行 為，文中將一薄片的兩個長邊以夾鉗固定，再將其彎成柱狀殼，兩邊夾持 端相距一固定距離 ，並與水平面維持一固定角度d  (如圖 A.1)，然後在薄 殼中心施加一向下的位移負荷。在其實驗中隨著位移負荷的增加，結構連 續產生四個特殊的幾何變形，如圖 A.2 所示。第一個變形是中心向下位移 時，薄殼中心附近出現兩個對稱 X、Y 軸的 d-cone (developable cone) (圖 A.2a)。第二個變形是中心向下位移 時，出現兩個新的 d-cone，而四 個 d-cone 圍成一個對薄殼中心轉了一個角度的菱形(圖 A.2b)。第三個變形 是中心向下位移12 時，四個 d-cone 的連線形成一個梯形(圖 A.2c)。第 四個變形是中心向下位移 時，梯形底邊兩個 d-cone 移到薄殼自由端的 邊界時，產生一個不連續的變化，使薄殼變成波浪狀的圓柱殼(圖 A.2d)。文 獻[13]提到殼的長度若不夠長，則無法觀察到菱形及後來的梯形圖形這兩種 型態，會提早產生波浪狀的柱狀殼。薄片的尺寸和材料參數如下：

mm 4

mm 5 . 11

mm 5 .

mm 15

長度L35cm 寬度W 17.5cm 厚度h 0.35mm

楊氏係數E 3.810⁹N/m²

20

 

cm d 16.5

 

X d Y Z

mm d 165

20





^



圖 A.1 柱狀薄板結構示意圖

圖 A.2 文獻[13]實驗所觀察到四種變形轉換(a-d)及 對應於 a-c 圖結構的上視圖(e-g)

附錄 B

B.1 面積座標的定義

如圖 B.1 所示，x、y為三角形中任一點 在直角座標中的座標值，將 點與三角形的三個頂點作連線，就形成了三個小三角形，三個頂點1、2、3 相對應的三個小三角形的面積分別為 、 、 ，令面積座標

P

A2

P

A1 A₃ A

L₁  A¹

  (B.1)

A L₂  A²

  (B.2)

A L₃  A³

  (B.3)



_{21 31 31 21}

2

1 x y x y

A 



(B.4)

(B.5)

j i

ij x x

x  

(B.6)

j i

ij y y

y  

其中A為三角形 123 的面積， 、 代表三角形頂點 的x_i y_i i x和 座標。y 、、

稱為三角形中P點的面積座標，固

(B.7) A

A A

A₁  ₂  ₃ 

由(B.1)式至(B.3)式與(B.7)式可以得出

1



 

 (B.8)

本文在 QST 平面元素推導中， 使用、、為面積座標表示方式，且、

 、之間不是互相獨立，因此任意函數可以表示為 

 

,

L1

。在 RQT 板元素推導中， 使用 、 、 為面積座標表示方式，且 、 、 之 間視為互相獨立，因此任意函數可以表示為

L1 L₂ L₃ L₂ L₃



L₁,L₂, L₃





  。

B.2 面積座標與直角座標的關係

(B.11)式可以利用(B.10)式微分得出



i3，則 j 1、k 2

圖 B.1 面積座標表示方法

P

¹

A A

2

A

3

1

2 y 3

x

附錄 C

在(2.1)式與(2.2)式中 及其對

C.1 QST 元素的形狀函數及其微分

Ni



i1,_,9



 、的偏微分，可表示為：

i N_i N_i,_ N_i,_

1 ²(32)2a ⁶⁽¹⁾^2b 6(1)2c 2 ² a/2 ⁽ ²⁾^b/2 2 c/2 3 ² a/2 ^2 ^b/2 ( 2)c/2 4 ²(32) 2a ⁶⁽¹⁾ ^2b _2c

5 ²(1)a (23)^b ^c

6 ² a/2 ² ^b/2 ² c/2

7 ²(32) 2a 2b 6(1)2c

8 ² a/2 ² b/2 ² ^c/2

a



1 )

2( 

 b (23)c

9 中 其





  1  (B.1)





a (B.2)

) ( 

 



b (B.3)

( )

 



c (B.4)

在(2.73)式中的N_u,_x、N_u,_y、N_v,_x、N_v,_y為：

形狀函數 的表示式為：

下標替換的順序(cyclic permutation)：

1 3

3 2

2 1







附錄 D

此元素座標系統為文獻[25]所提出。元素座標系統為原點在形心 O ，並且當 前變形位置元素座標軸轉 角使得當前變形位置的元素各節點對應初始未

座標軸轉

變形時的元素各節點的距離平方和為最小，如圖 D.1 所示。

當前變形位置元素 角後，角逆時鐘方向為正，元素各節點 j 座 標值的關係式：

j j

j x y

x cos sin (D.1)

j j

j x y

y sin cos (D.2) 當前變形位置元素各節點 j 座標值表示 y_j)，

為(x_j, 當前變形位置元素座標 軸轉角後元素各節點 j 座標值表示為(x_j,y_j)。

當前變形位置元素座標軸轉角後的當前變形位置的元素各節點對應初始

未變形時的元素各節點的距離平方和：

] ) (

) [(

3

1 j

j j j

j ^{0 2 0 2}



^{x  x  y  y} (D.3)

標值表示為 。

1)和(D.2)代入(D.3)可得：

(D.4)

對

) , (⁰x_j ⁰y_j 初始未變形時的元素各節點 j 座

(D.

1 j

j sin ) ( sin cos ) ]

[(cos

3 0 2 0 2



 x ^  yj^ xj ^{ }  xj ^  yj^ yj

(D.4)  微分可得：

(D.5)

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