未來研究方向

因本文採用的 DKT 板元素沒有定義元素內部的側向位移w，所以本文 推導的兩種元素幾何剛度矩陣都僅是近似的幾何剛度矩陣，這也許會影響

平衡迭代的收斂速度和偵測平衡路徑上分歧點及挫屈模態的準確性，因此 以後研究可以考慮把 DKT 板元素替換成具有側向位移場w的三角板元素，

或用不同的側向位移場，以本文採用的共旋轉 total Lagrangian 法推導元素 幾何剛度矩陣，或以共旋轉推導法來推導元素幾何剛度矩陣，並探討其對 平衡迭代和偵測平衡路徑上分歧點及挫屈模態的影響。雖然本文中決定元 素座標系統的方法是文獻上最常用的方法，但文獻[49]中提到在共旋轉推導 法中，不同的元素座標會影響分析的結果，並提出數種決定元素座標的方 法，因此以後的研究可以提出其他決定元素座標的方法，並比較各種方法 用在本文之殼元素的結果，選擇對本文之殼元素本身最適合的元素座標系 統，應可提高薄殼結構之幾何非線性分析的精度或收斂速率。

參考文獻 參考文獻參考文獻 參考文獻

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

[2] J. L. Batoz, C. L. Zheng, F. Hammadi, “Formulation and evaluation of new triangular, quadrilateral, pentagonal and hexagonal discrete Kirchhoff plate/shell elements”, International Journal for Numerical Methods in Engineering, 52, 615-630, 2001.

[3] N. Carpenter, H. Stolarski, T. Belytschko, “A flat triangular shell element with improved membrane interpolation”, Communications in Applied Numerical Methods, 1, 161-168, 1985.

[4] Y. K. Cheung, W. J. Chen, “Refined non-conforming triangular elements for analysis of shell structures”, International Journal for Numerical Methods in Engineering, 46, 433-455, 1999.

[5] J. G. Kim, J. K. Lee, Y. K. Park, “A new 3-node triangular flat shell element”, Communications in Numerical Methods in Engineering, 18, 153-159, 2002.

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

[7] K. M. Hsiao, “Nonlinear analysis of general shell structures by flat triangular shell element”, Computer & Structure, 25,665-675, 1987.

[8] R. Levy, E. Gal, “Geometrically nonlinear three-noded flat triangular shell element”, Computer & Structure, 79, 2349-2355, 2001.

[9] Y. X. Zhang, Y. K. Cheung, “A refined non-linear non-conforming triangular plate/shell element”, International Journal for Numerical Methods in Engineering, 56, 2387-408, 2003.

[10] J. M. Battini, C. Pacoste, “On the choice of the linear element for

corotational triangular shells”, Computer Methods in Applied Mechanics and Engineering, 195, 6362-6377, 2006

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

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

[13] J. M. Battini, “A modified corotational framework for triangular shell elements”, Computer Methods in Applied Mechanics and Engineering, 196, 1905-1914, 2007.

[14] 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, 69, 859-885, 2007.

[15] P. Khosravi, R. Ganesan, R. Sedaghati, “An efficient facet shell element for corotational nonlinear analysis of thin and moderately thick laminated composite structures”, Computer & Structure, 86, 850-858, 2007.

[16] Z. X. Li, L. Vu-Quoc, “An efficient co-rotational formulation for curved triangular shell element”, International Journal for Numerical Methods in Engineering, 72, 1029-1062, 2007.

[17] L. Kang, Q. Zhang, Z. Wang, “Linear and geometrically nonlinear analysis of novel flat shell elements with rotational degrees of freedom”, Finite Elements in Analysis and Design, 45, 386-392, 2009.

[18] J. M. Battini, “A non-linear corotational 4-node plane element”, Mechanics Research Communications, 35, 408-413, 2008.

[19] D. J. Allman, “A Compatible Triangular Element Including Vertex Rotations for Plane Elasticity Analysis”, Computers and Structures, 19, 1-8, 1984.

[20] P. G. Bergan, C. A. Felippa, “A triangular membrane element with rotational degrees of freedom”, Computer Methods in Applied Mechanics and Engineering, 50, 25-69, 1985.

[21] D. J. Allman, “Evaluation of the constant strain triangle with drilling rotations”, International Journal for Numerical Methods in Engineering, 26, 2645-2655, 1988.

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

[23] C. A. Felippa, “Refined finite element analysis of linear and nonlinear two-dimensional structures”, Ph.D. Dissertation, Department of Civil Engineering, University of California at Berkeley, Berkeley, CA, 1966.

[24] C. A. Felippa, “A study of optimal membrane triangles with drilling freedoms”, Computer Methods in Applied Mechanics and Engineering, 192, 2125–2168, 2003.

[25] I. Holand and P. G. Bergan, “Higher order finite element for plane stress”, Discussion, Journal of the Engineering Mechanics Division, 2, 698-702, 1968.

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

[27] G. P. Bazeley, Y. K. Cheung, B. M. Irons, O. C. Zienkiewicz,

“Triangular elements in plate bending—conforming and nonconforming solutions”, Proceeding First Conference on Matrix

in Structural Mechanics, Air Force Institute of Technology, Dayton, Ohio, 66–80, 1966.

[28] R. W. Clough, J. L. Tocher, “Finite element stiffness matrices for analysis of plate bending”, Proc. Conf. on Matrix Methods in Structural Mechanics, 515-545, 1965

[29] J. A. Stricklin, W. E. Haisler, P. R. Tisdale, R. Gunderson, “A rapidly converging triangular plate element”, AIAA J., 7, 180-181, 1969.

[30] D. J. Allman, “Triangular finite element bending with constant and linearly varying bending elements”, High Speed Comput. Elastic Struct. Tom 1, 105-107, 1971.

[31] L. M. Tang, W. J. Chen, Y. X. Liu, “The quasi-conforming element for the finite element analysis”, J. Dalian Inst. Technol., 19, 19-35, 1980.

[32] J. L. Batoz, “An explicit formulation for an efficient triangular plate-bending element”, International Journal for Numerical Methods in Engineering, 18, 1077–1089, 1982.

[33] Y. K. Cheung, W. J. Chen, “Refined triangular discrete Kirchhoff plate element for thin plate bending, vibration, and buckling analysis”, International Journal for Numerical Methods in Engineering, 41,1507-1525, 1998.

[34] C. A. Felippa, “Recent advances in finite element templates”, University of Colorado, 2000.

[35] T. Wenzel, H. Schoop, “A non-linear triangular curved shell element”, Communications in Numerical Methods in Engineering, 20, 251-264, 2004

[36] E. M. Kasparek, “An efficient triangular plate element with C1-continuity”, International Journal for Numerical Methods in

Engineering, 73, 1010-1026, 2008.

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

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

[39] H. Goldstein, Classical Mechanics, Addision-Wesley Publishing Company, 1980.

[40] G. Dhatt, G. Touzot, The Finite Element Method Displayed, John Wiley & Sons, 1984.

[41] A. P. Boresi, K. P. Chong, Elasticity In Engineering Mechanics, John Wiley & Sons, 1999.

[42] R. D. Cook, Concepts and Applications of Finite Element Analysis, John Wiley & Sons, 1981

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

[44] D. J. Dawe, Matrix and Finite Element Displacement Analysis of Structures, Clarendon Press, 1984.

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

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

[47] Y. C. Fung, Foundation of Solid Mechanics, Prentice-Hall, Englewood Cliffs, N.J., 1965.

[48] E. H. Boutyour, H. Zahrouni, M. Potier-Ferry, M. Boudi, Bifurcation points and bifurcated branches by an asymptotic numerical method and Pade approximants, International Journal for Numerical Method in Engineering, 60, 1987–2012, 2004.

[49] J. M. Battini, C. Pacoste, “On the choice of local element frame for corotational triangular shell elements”, Communications in Numerical Method in Engineering, 20, 819-825, 2004

表 4.1 圓柱殼片段受到單點集中力作用所偵測到的挫屈負荷(例題 4.3)

Element Mesh Kg 元素座標 Buckling load(N)

10×10 (1) (a) 528.513

10×10 (1) (b) 528.496

10×10 (2) (a) 528.315

Present

10×10 (2) (b) 528.296

10×10 (3) (a) 528.209

10×10 (3) (b) 528.192

表 4.2 Lateral torsional buckling 所偵測到的挫屈負荷(例題 4.4)

Element Mesh Kg 元素座標 Buckling load

(2+2+2)×14 (1) (a) 2889.74

(2+2+2)×14 (1) (b) 2888.91

(2+2+2)×14 (2) (a) 2906.06

Present

(2+2+2)×14 (2) (b) 2905.41

(2+2+2)×14 (3) (a) 2887.08

(2+2+2)×14 (3) (b) 2886.25

(4+4+4)×60 (1) (a) 2746.74

(4+4+4)×60 (1) (b) 2750.94

(4+4+4)×60 (2) (a) 2889.61

(4+4+4)×60 (2) (b) 2752.41

(4+4+4)×60 (3) (a) 2738.66

(4+4+4)×60 (3) (b) 2751.19

OPT DKT[10] (2+2+2)×15 - (b) 2618.26

(4+4+4)×60 - (b) 2529.29

表 4.3 Simply supported compressed plate 所偵測到的挫屈負荷(例題 4.5)

Element Mesh Kg 元素座標 Buckling load

Present 4×6 (1) (a) 58.4289

4×6 (1) (b) 58.4491

4×6 (2) (a) 58.2711

4×6 (2) (b) 58.2915

4×6 (3) (a) 58.2699

4×6 (3) (b) 58.2902

4×10 (1) (a) 56.5610

4×10 (1) (b) 56.5711

4×10 (2) (a) 56.4145

4×10 (2) (b) 56.4253

4×10 (3) (a) 56.4221

4×10 (3) (b) 56.4247

6×12 (1) (a) 57.1913

6×12 (1) (b) 57.1983

6×12 (2) (a) 57.0414

6×12 (2) (b) 57.0489

6×12 (3) (a) 57.0411

6×12 (3) (b) 57.0483

8×20 (1) (a) 57.2675

8×20 (1) (b) 57.2570

8×20 (2) (a) 57.1089

8×20 (2) (b) 57.1201

8×20 (3) (a) 57.1176

8×20 (3) (b) 57.1073

OPT DKT[10] 4×10 - (b) 59.1703

8×20 - (b) 58.4762

表 4.4 Transversally loaded T profile 所偵測到的挫屈負荷(例題 4.6)

Element Mesh Kg 元素座標 Buckling load

Prsent (2+2+4)×20 (1) (a) 2799.01

(2+2+4)×20 (1) (b) 2798.68

(2+2+4)×20 (2) (a) 2797.21

(2+2+4)×20 (2) (b) 2796.70

(2+2+4)×20 (3) (a) 2797.16

(2+2+4)×20 (3) (b) 2796.65

(3+3+5)×30 (1) (a) 2848.46

(3+3+5)×30 (1) (b) 2847.75

(3+3+5)×30 (2) (a) 2846.41

(3+3+5)×30 (2) (b) 2845.71

(3+3+5)×30 (3) (a) 2846.32

(3+3+5)×30 (3) (b) 2845.61

(5+5+8)×50 (1) (a) 2829.86

(5+5+8)×50 (1) (b) 2829.02

(5+5+8)×50 (2) (a) 2827.77

(5+5+8)×50 (2) (b) 2826.93

(5+5+8)×50 (3) (a) 2827.67

(5+5+8)×50 (3) (b) 2826.83

OPT DKT[10] (2+2+4)×20 - (b) 3103.86

(5+5+8)×50 - (b) 2890

表 4.5 Channel section in torsion 所偵測到的極限點(例題 4.7)

Element Mesh Kg 元素座標 J 點側向位移(-W) Limit point

(1+2+1)×20 (2) (a) 18.8209 1097.65

(1+2+1)×20 (2) (b) 18.9037 1104.65

(2+4+2)×40 (2) (a) 20.5250 1104.07

Present

(2+4+2)×40 (2) (b) 20.0766 1106.21

OPT DKT[10] (1+2+1)×20 - (b) - 1011.92

(3+8+3)×56 - (b) - 1080.86

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

,U

圖 2.2 旋轉向量

n R′

R

φ

圖 2.3 元素變形示意圖 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.4 元素節點 j中心面之x_ij^B軸受旋轉向量₋^θ_nj作用的示 意圖

θ

nj

−

B

x

₁j

B

x

₂j

B

x

₃j

B j E

x x

₃

,

₃^′

θ

nj

j

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

B

x

₁j^′ B

x

₁j^′ 0

B

x

₂j^′

B

x

₂j^′ zj 0

θ

θ

zj

E B

j B

j

x x x

_{3 3 3}

0 ^′

,

^′

,

j

X₃G

X₂G

X₁G

x₁E

x₃E

x₂E

h v³

u3

w3 3

θ

z

3

θ

y

L

3

2

1

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

s

s

s α

12

α

13

α

²³

n

12

n

13

n

23

x

₁E

x

₂E

5

1

3

6

4

2

圖 2.7 DKT 元素的節點及其三邊上的局部座標示意圖

(a)

(b) (c)

圖 4.1 半圓環受到單點集中力作用 (a)結構尺寸示意圖 (b)網格 18×1 示意 圖(c)網格 18×2 示意圖

F 2

A

=1 d

=20 R

θ V X₂^G,

U X₁^G, 25 . 0 , 10⁷ =

= v

E

A A

圖 4.2 半圓環受到單點集中力作用，不同網格下之負荷位移曲線圖。

0 5 10 15 20 25 30

0 1 2 3 4 5 6 7 8 9 10 11 12 13

(a)

(

b)

ref.[18]

20x1 30x1 40x2 80x4 100x6 L o a d F ( 1 0

3

)

Vertical displacement of point A (-V

_A

)

(a)

(b)

圖 4.3 直角構架受到端點剪力作用 (a)結構尺寸示意圖 (b)網格 M21 與網格 M22 示意圖

F A

9

1 1

9 U X₁^G, V

X₂^G, θ

1 , 3 . 0 , 10

3× ⁷ = =

= t

E ν

M21 M22

圖 4.4 直角構架受到端點剪力作用，不同網格下之負荷位移曲線圖。

0 1 2 3 4 5 6 7

0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0

(a) (b)

ref.[18]

M21 M22

L o a d F ( 1 0

4

)

Horizontal displacement of point A (U

_A

)

(a)

(b) (c)

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

L=254mm R=2540mm h=6.35mm θ=0.1rad

E=3102.75MPa ν=0.3

L

L

R

h

θ θ

E D

A

C

B

R R

θ θ

E

U X

,

W

Z

,

F

A B

C D

圖 4.6 圓柱殼片段受到單點集中力作用之負荷位移曲線圖

0 5 10 15 20 25 30

-0.4 -0.2 0.0 0.2 0.4 0.6

0.8 ref.[48]

Primary path Secondary path

L o ad F ( k N )

Vertical displacement(-W) at E (mm)

(a)

(b)

圖 4.7Lateral torsional buckling (a)結構尺寸示意圖 (b)網格(1+2+1)×2 示意圖

=150

L E=21000

=10

b ν =0.3

=10 h

=1 t

A B

C D

E F

G I

t

W X₃^G,

V X₂^G,

U X₁^G,

λ I

G

F E

D

C B

A

L

b h

(a)

(b)

圖 4.8Simply supported compressed plate (a)結構尺寸示意圖 (b)網格 4×6 示意圖

=100

L ，b =50，t =2，E =100，ν =0.3

A B

C D

L

A B

C D

λ λ

U X₁^G, V X₂^G,

b

(a)

(b)

圖 4.9Transversally loaded T profile (a)結構尺寸示意圖 (b)網格(2+2+3)×4 示意圖

=450

L E =70960

=38

b ν =0.321

=65 h

=1 t

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

h

b

L

J

(a)

t (b)

圖 4.10Channel section in torsion (a)結構尺寸示意圖 (b)網格(1+2+1)×2 示意圖

A B

C D

E F

G I

=1100

L E =2.1×10⁶

=75

b ν =0.3

=25 h

3 .

=1 t

C B

U X₁^G, V

X₂^G,

W X₃^G,

D A

E

F G

I

K J

M

λ λ

λP

λP

L

b h

圖 4.11 Channel section in torsion，不同網格下之負荷位移曲線圖。

0 5 10 15 20 25 30

0 200 400 600 800 1000 1200

ref.[13]

(1+2+1)x20 (2+4+2)x40

L o ad f ac to r

λ

Displacement(-W) at J

(a)

(b)

圖 4.12 直角梁受到單點集中力作用 (a)結構尺寸示意圖 (b)網格(2+3)×2 示意圖

=1400

L E =193300 75

.

=47

b ν =0.3 75

.

=72 h

5 .

=2 t

A

B

C

D E

F

L

U X₁^G, V

X₂^G,

W X₃^G,

A

B

C

D

F

P

t

E

b

h

圖 4.13 直角梁受到單點集中力作用，不同網格下之負荷位移曲線圖。

0 20 40 60 80 100

0 1 2 3 4 5 6 7 8 9 10

(a) (b)

OPT[10]

ALL[10]

CM FM L o ad P

(

1 0

3 )

Axial displacement of point A

(a)

(b)

圖 4.14 懸臂圓柱殼受到單點集中力作用 (a)結構尺寸示意圖 (b)網格 16×16 示意圖

R

t

F 2

F 2

L U X₁^G, V

X₂^G,

W X₃^G, A

end Free

end Clamped 016

.

=1

R E=2.0685×10⁷ 048

.

=3

L ν =0.3 03

.

=0 t

B

C D

A

B

C D

圖 4.15 懸臂圓柱殼受到單點集中力作用之負荷位移曲線圖

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

ref.[13]

Present

L o ad F

Vertical displacement at A / R

(a)

(b)

圖 4.16 半球殼受到單點集中力作用 (a)結構尺寸示意圖 (b)網格 12×12 示意圖

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

B A

C

=10 R

5 .

=0 t

104

= E

2 .

=0 ν

圖 4.17 半球殼受到單點集中力作用之負荷位移曲線圖

0 1 2 3 4 5 6 7 8 9

0 5 10 15 20 25 30 35

40 ref.[13]

Present

B A

Load F

Normal displacements at A (with reverse sign) and B

附錄 附錄 附錄

附錄 A 面積座標面積座標面積座標(area coordinates) 面積座標

A.1 面積座標的定義面積座標的定義面積座標的定義面積座標的定義

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

A A₁

λ = ^(A.1)

A A₂

ξ = ^(A.2)

A A₃

η = ^(A.3)

(

21 31 31 21

)

2

1 x y x y

A= − ^(A.4)

j i

ij x x

x = − ^(A.5)

j i

ij y y

y = − ^(A.6)

其中A為三角形 123 的面積，x_i、y_i代表三角形頂點i的x和y座標。λ^、ξ^、 η^{稱為三角形中P點的面積座標，固}

A A A

A₁ + ₂ + ₃ = (A.7)

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

=1 + +ξ η

λ ^(A.8)

因λ^、ξ ^、η之間不是互相獨立的，因此在本文中僅用ξ ^、η^{表示三角形中} 任意點的面積座標，如圖 A.2。

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

P點之面積座標與直角座標之間的關係可表示成[42]



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

圖 A.2 面積座標示意圖

P

¹

A A

2

A

3

1

2 3

x y

(0,0)

1 2

3

(1,0) (0,1)

ξ

η

附錄 附錄附錄

附錄 B 不完整三階埃爾米特元素不完整三階埃爾米特元素不完整三階埃爾米特元素的形狀函數及其微分不完整三階埃爾米特元素的形狀函數及其微分的形狀函數及其微分的形狀函數及其微分

在(2.15)式與(2.16)式中N_i

(

i =1 L^{, ,9}

)

及其對ξ^、η的偏微分，可表示為[40]：

i N_i N_i,_ξ N_i,_η

1 λ²(3−2λ)+2a ⁶λ⁽⁻¹⁺λ^{)+ 2b 6}λ⁽−¹+λ⁾+ ^2c 2 ξλ² +a/2 λ⁽λ ⁻²ξ^)+b/2 −²ξλ +^c/2 3 ηλ² +a/2 ⁻²λη ^+b/2 λ⁽λ −²η⁾+^c/2 4 ξ²(3−2ξ)+2a ⁶ξ⁽¹⁻ξ^{)+ 2b} 2c

5 ξ²(−1+ξ)−a ξ⁽⁻²⁺³ξ^)−b −^c

6 ξ²η +a/2 ²ξη ^+b/2 ξ² +c/2

7 η²(3−2η)+2a 2b 6η(1−η)+2c

8 ξη² +a/2 η² +^{b/2 2}ξη ^+c/2

9 η²(−1+η)−a −^b η(−2+3η)−c

其中

η ξ

λ =1− − (B.1)

ξηλ

=

a (B.2)

) (λ ξ η −

=

b (B.3)

) (λ η ξ −

=

c (B.4)

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

η

ξ η

ξ^{, , , ,}

,_{x x u x u}

u N N

N = + ^(B.5)

η

ξ η

ξ^{, , , ,}

,_{y y u y u}

u N N

N = + ^(B.6)

η

附錄 C DKT 元素的形狀函數

在(2.54)式裡面的H_x與H_y分別有 9 個分量，其表示式為[32]

H_x1 =1.5(a₆N₆ −a₅N₅) H_x2 =b₅N₅ +b₆N₆ H_x3 = N₁ −c₅N₅ −c₆N₆ H_x4 =1.5(a₄N₄ −a₆N₆) H_x5 =b₆N₆ +b₄N₄ H_x6 = N₂ −c₆N₆ −c₄N₄ H_x7 =1.5(a₅N₅ −a₄N₄) H_x8 =b₄N₄ +b₅N₅ H_x9 = N₃ −c₄N₄ −c₅N₅ H_y1 =1.5(d₆N₆ −d₅N₅) H_y2 =−N₁ +e₅N₅ +e₆N₆ H_y3 =−H_x2

H_y4 =1.5(d₄N₄ −d₆N₆) H_y5 =−N₂ +e₆N₆ +e₄N₄ H_y6 =−H_x5

H_y7 =1.5(d₅N₅ −d₄N₄) H_y8 =−N₃ +e₄N₄ +e₅N₅ H_y9 =−H_x8

其中

其中ξ^、η^{是元素面積座標裡的L}₂^與L₃[43]，如圖 C.1 所示。

1

圖 C.1 元素面積座標的示意圖

圖 C.2 殼元素節點的自然座標示意圖

x y

A

1

A

2

A

3

1

2 3

P

ξ η

6

5 4

3

1 2 ) 0 , 0

( ,0) (1,0)

2 (1 ) 1 , 0 (

2) ,1 0

( )

2 ,1 2 (1

附錄附錄

附錄附錄 D 元素幾何剛度矩陣元素幾何剛度矩陣元素幾何剛度矩陣元素幾何剛度矩陣

為了要推導精確的幾何剛度矩陣，必須考慮非線性的位移－應變關 係，文獻[43]中以完整的 Green’s strain 及 total Lagrangian 推導法推導出一三 維元素之切線剛度矩陣的通式，本附錄將重複文獻[43]中的推導，並加以修

 將(2.37)式、(2.38)式、(2.41)式代入(D.7)式可得：

Gq

θ= (D.14) 其中G 為u- θ 關係矩陣，q 為元素節點自由度向量。

將(D.14)式代入(D.13)式可得：

q

D.2 元素幾何剛度矩陣元素幾何剛度矩陣元素幾何剛度矩陣元素幾何剛度矩陣的推導的推導的推導的推導 將(D.20)式、(D.21)式代入(D.19)式可得：

dV matrix)，k_L為初始位移矩陣或大位移矩陣(initial displacement matrix or large

displacement matrix)，k_σ 為元素幾何剛度矩陣。

將(D.16)式代入(D.24)式可得：



將(D.26)式代入(D.25)可得：

V dV

b



將(D.31)式、(D.44)式代入(D.43)可得：

dzdA

∫ ∫

₋ 

附錄

將(E.1)至(E.4)式代入(E.11)式可得

( )

∫

⁻

= ^{t t} _{A t u n v}

P

tf q T t N t N dA

q_θ δ _{θ ξθ}

δ ^(E.12)

由(E.12)式可得

( )

∫

⁻

= ^t _{A t u n v}

P T t N t N dA

f _ξθ (E.13)

在(E.13)式中的受力面積dA為 ξ

d tL

dA= ₁₂ (E.14)

其中L₁₂為邊 12 的長度。將(E.14)式代入(E.13)式可得

( )

∫

⁻

=^tL_{12 ξθ 0}^{1 t}t u ^tn v ^dξ

t

P T N N

f (E.15)

E.2 結構分佈力之等效節點外力結構分佈力之等效節點外力結構分佈力之等效節點外力結構分佈力之等效節點外力

如圖 E.2、圖 E.3 與圖 E.4 所示，考慮邊 AB 被m+1個節點由節點 A 向節點 B 離散成等長的m段，則邊 AB 的長度與(E.15)式中的L₁₂之關係為

mL12

L_AB = ^(E.16)

將(E.16)式代入(E.15)式中，並將t_t與t_n分開寫成

= _ξθ

∫

₀^{1t d}ξ m

tL

u t AB t

t T N

f (E.17)

−

∫

= _{ξθ 0}^{1t d}ξ m

tL

v n AB t

n T N

f (E.18)

在本文例題中負荷分佈情形最高為二次變化，因此以下將就不同負荷分佈 作分類：

1.均勻分佈力(uniform distributed surface traction)

(1)邊 AB 受到平行於邊的均勻分佈力，如圖 E.2(a)所示，其合力大小為 f_ut，

則第k 段所受分佈力可表示成：

AB ut t

k

tL

t = f ^(E.19)

將(E.19)式代入(E.17)式並積分可得

F1

將(E.23)式代入(E.18)式並積分可得

F2 2.二次分佈力(quadratic distributed surface traction)

(1)邊 AB 受到平行於邊的二次分佈力，如圖 E.3(a)所示，其合力大小為

將(E.27)式代入(E.17)式並積分可得

3 T F3

將(E.35)式代入(E.18)式並積分可得

3 T F4

3.純彎分佈力(pure bending distributed surface traction)

將(E.39)式代入(E.18)式並積分可得

2 T F5

在(E.22)式、(E.26)式、(E.30)式、(E.38)式以及(E.42)式中所計算出的^kf_t與

n

kf 須經(2.57)式座標轉換至固定總體座標與節點基礎座標後，並計算k從 第1段到m段，然後將其總合成負荷外力向量P。

(a) (b)

圖 E.1 元素受力示意圖 (a)σ_t ^(b)σ_n

(a) (b)

圖 E.2 均勻分佈力示意圖 (a) f_ut (b) f_un

(a) (b)

圖 E.3 二次分佈力示意圖 (a) f_qt (b) f_qn

圖 E.4 純彎分佈力示意圖 M

1 2 L m

B A

fqt f_qn

1 2 L m 1 2 L m

A B

B A

fut f_un

1 2 L m 1 2 L m

A B

B A

1 2

E 3 x₂

x₁E

tt

1 2

E 3 x₂

x₁E

tn

附錄

令(F.5)為 0 可得：

∑

∑

=

=

+

−

= ₃

1

0 0

3

1

0 0

) (

) (

tan

j

j j j j j

j j j j

x y y x

x y y x

α (F.6)

此時當前變形位置元素座標軸轉α角使得當前變形位置的元素各節點對 應初始未變形時的元素各節點的距離平方和為最小。

圖F.1元素變形示意圖 u3

v3

Initial

Current

x₁₂B^′ 0

x₁₂B^′ 0

x₁₂B^′

2

θ

z

v x x₂^E, ₂^E,

0

u x x₁^E, ₁^E,

0

02

2

01

1 3

03

O α

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