Phân tích kết cấu vỏ bằng phần tử MITC3 được làm trơn trên cạnh (ES-MITC3)

Nội dung text: Phân tích kết cấu vỏ bằng phần tử MITC3 được làm trơn trên cạnh (ES-MITC3)

1. PHÂN TÍCH KẾT CẤU VỎ BẰNG PHẦN TỬ MITC3 ĐƯỢC LÀM TRƠN TRÊN CẠNH (ES-MITC3) ANALYSIS OF SHELL USING AN EDGE-BASED SMOOTHED FINITE ELEMENT METHOD (ES-MITC3) Quách Văn Nghiệp. Học viên cao học trường ĐHSPKT TPHCM. TÓM TẮT Trong đề tài này, phương pháp phần tử hữu hạn trơn trên cạnh phần tử ES-FEM (the edge smoothed finite element method) được phát triển cho phần tử MITC3 dùng để phân tích tĩnh kết cấu vỏ phẳng theo lý thuyết biến dạng cắt bậc nhất. Kết cấu được mô phỏng bằng các phần tử tam giác ba nút với năm bậc tự do cho mỗi nút. Nhờ vào kỹ thuật làm trơn ES-FEM, ma trận độ cứng phần tử được tính toán dựa vào tích phân trên cạnh của phần tử. Từ khóa: Vỏ, phần tử MITC3, lý thuyết biến dạng cắt bậc nhất. ABSTRACT In this topic, the edge smoothed finite element method ES-FEM was developed for elements MITC3 used to analyze static structural of shell under the first order shear deformation theory. The structural is modeled by triangular elements of three nodes with five degrees of freedom for each node. Thanks to a smoothing technique ES-FEM, element stiffness matrix is calculated based on the integral on the edge of the element. Keywords: Shell, element MITC3, the first order shear deformation theory. 1. Giới thiệu. Từ nhiều thập kỷ qua, cùng với sự phát triển vượt bậc của khoa học máy tính, các phương pháp số được nghiên cứu và đạt được những thành tựu quan trọng. Phương pháp phần tử hữu hạn (FEM) là một trong những phương pháp số được sử dụng rộng rãi với những ưu điểm nổi trội như kết quả tính toán có độ chính xác cao, mạnh mẽ khi phân tích các bài toán phức tạp, kết hợp tốt các lý thuyết khác nhau trong quá trình phân tích. Một số lý thuyết và phương pháp đã từng được giới thiệu như: lý thuyết tấm/vỏ cổ điển (vỏ mỏng) của Kirchhoff – Love (CPT) và lý thuyết tấm/vỏ dày Reissner và Mindlin [1], phần tử với tham số C0 Ahmad, Irons và Zienkiewicz [2,3] , các phần tử vỏ Reissner / Mindlin [4,5]. Mặc dù bao gồm các biến dạng cắt để phân tích vỏ dày nhưng khó khăn chính của phần tử là hiện tượng “khóa cắt” (shear locking) khi chiều dày tấm giảm dần. Năm 1980 Bathe và Dvorkin [6] đề xuất phương pháp nội suy hỗn hợp các thành phần ten xơ (mixed interpolation tensorial components viết tắt là MITC) đã giải quyết được các vấn đề về khóa cắt. Các phương pháp MITC rất hiệu quả trong giải quyết các bài toán tấm vỏ và cho kết quả tin cậy. Kỹ thuật này được Bathe và Dvorkin đề xuất với phần tử tứ giác 4 nút và phần tử tứ giác 8 nút (các MITC4 và MITC8) [6,7]. Sau đó, Bathe và Bucalem đã phát triển nó lên tứ giác 9 nút và 16 nút (MITC9 và MITC 16) [8]. Kỹ thuật MITC đã được phát triển cho các phần tử vỏ tam giác 3 nút MITC3 và tam giác 6 nút 1
2. MITC6 [9, 10, 11]. Tác giả Liu và cộng sự đã ứng dụng kỹ thuật làm trơn hóa biến dạng [12] để thiết lập công thức phương pháp PTHH trơn dựa trên phần tử con (cell) còn được gọi là SFEM hoặc CS-FEM cho bài toán 2D cơ vật rắn và sau đó CS-FEM được phát triển cho tấm và vỏ [13]. Năm 2008, tác giả Liu và cộng sự [14] đã đề xuất phương pháp PTHH trơn với cách tiếp cận dựa trên cạnh (the Edge-based Smoothed Finite Element Method – ES-FEM) cho bài toán phân tích tĩnh, dao động tự do và dao động cưỡng bức áp dụng cho bài toán phân tích cơ học vật rắn hai chiều. Tác giả Cui và cộng sự [15] với phương pháp PTHH trơn trên cạnh, tác giả Sơn và cộng sự [16] kết hợp PTHH trơn trên cạnh và trên nút đã giải quyết khá tốt các bài toán kết cấu vỏ. Nghiên cứu này sẽ xây dựng cơ sở lý thuyết công thức PTHH trơn cho phần tử MITC3, gọi là phần tử ES-MITC3, dùng để phân tích kết cấu vỏ phẳng dựa trên lý thuyết biến dạng cắt bậc nhất. So sánh và đánh giá hiệu quả của phần tử đề xuất với các công thức PTHH khác. 2. Cơ sở lý thuyết. 2.1 Lý thuyết biến dạng cắt bậc nhất: Lý thuyết tấm cổ điển của Kirchoff giả thiết là: bỏ qua ứng suất pháp vuông góc với mặt phẳng tấm dạng ( z 0). Đoạn thẳng vuông góc với mặt trung bình (mặt phẳng chia đôi chiều cao tấm) vẫn thẳng và vuông góc với mặt trung bình sau khi biến dạng. Hệ quả của giả thiết này là ta đã bỏ qua các thành phần biến dạng cắt ngang xz  yz 0. Tuy nhiên, khi tính tấm dày hoặc tấm nhiều lớp thì việc bỏ qua các biến cắt ngang sẽ không hợp lý. Trong nhiều năm qua, nhiều giả thuyết biến dạng cắt được đưa ra và phát triển để khắc phục những hạn chế của lý thuyết tấm cổ điển, trong đó có Reissner và Mindlin [1] với đề xuất về lý thuyết biến dạng cắt bậc nhất. Theo lý thuyết tấm biến dạng cắt bậc nhất Reissner và Mindlin xem rằng các đoạn thẳng pháp tuyến vuông góc với mặt trung gian trước khi biến dạng, sau khi biến dạng vẫn thẳng nhưng không vuông góc với mặt phẳng biến dạng (xz 0,  yz 0) tức là góc xoay trung bình của mặt cắt ngang có thể xem như góc xoay của pháp tuyến cộng thêm góc xoay do biến dạng cắt gây ra. Mặt khác, ứng suất pháp vẫn xem như bỏ qua ( z 0) như giả thiết Kirchoff. 2.2 Trường chuyển vị: Lý thuyết biến dạng cắt bậc nhất được áp dụng cho luận văn này. Trong tọa độ vuông góc, 3 thành phần chuyển vị được thể hiện như sau: u u0 zy (2.1) v v0 zx w w0 Trong đó: u0,, v 0 w 0 là thành phần chuyển vị của điểm trên mặt phẳng trung bình của vỏ (z = 0) theo phương x, y và z ; ,  là góc xoay pháp tuyến của mặt phẳng trung x y bình quanh các trục x, y. 2
3. 2.3 Trường biến dạng: Từ trường chuyển vị (2.1) các vec-tơ biến dạng được xác định như sau: u   0 z y   x  x xx v0 x ε  yy  z  ε m z ε b (2.2a) y  y  xy  u  v   0 0 z()y x y  x  y  x  w0   y  xz  x εs   (2.2b)  w yz  0  x y  Trong đó εm là biến dạng màng,εb là biến dạng uốn và εs là biến dạng cắt được xác định bởi:  u   y 0 x x w  0  y v0 x x εm  , εb  ,ε (2.3) y s w  y 0  x u  v   y  0 0 y x y  x  y  x  2.4 Trường ứng suất trong tọa độ cục bộ Trạng thái ứng suất phẳng được sử dụng để phân tích kết cấu tấm, vỏ như hình 2.1. x xz x  y xy yz  yx y Hình 2.1 Sự phân bố các thành phần ứng suất theo bề dày tấm 3
4. Theo định luật Hooke, mối quan hệ giữa ứng suất và biến dạng được biểu diễn như sau:     x x y  D m  y  (2.4) xy   xy  xz   xz   Ds  (2.5) yz   yz  Trong đó: DDm, s là các ma trận độ cứng vật liệu tương ứng được xác định theo công thức (2.6) và (2.7) 1v 0 E D v 1 0 m 1 v2 1 v 0 0 2 (2.6) 1 0 E 1 0 Ds G (2.7) 0 1 2(1 v ) 0 1 Trong đó: v, E là hệ số possion và mô đun đàn hồi của vật liệu. Công thức (2.4) và (2.5) có thể được viết lại dưới dạng như sau : σ Dm ε D m() ε m z ε b D m m z D m ε b (2.8) τ Ds ε s (2.9) 2.5 Nội lực trong vỏ trong tọa độ cục bộ Mx Mxy N My x Nxy Qx Myx Nxy Ny Q y Hình 2.2. Các thành phần nội lực trong tấm Nội lực chính là hợp lực của ứng suất phân bố theo bề dày tấm trên một đơn vị chiều như hình 2.2 , được xác định như sau: 4
5. - Lực màng: N h   x 2 x N  dz (2.10) y  y  h Nxy 2  xy  Hay: h 2 N σdz (2.11) h 2 Với: T N NNNx y xy (2.12) T σ x  y  xy (2.13) Thay công thức (2.8) vào công thức (2.11) ta được: h h h h h 2 2 2 2 2 NDD () zε dz D  dzz D ε dz  D dz ε zdz D mmmb mm mb mm bm h h h h h 2 2 2 2 2 (2.14) Đặt: h 1v 0 2 Eh AD dz v 1 0 (2.15) m 2 h 1 v 1 v 2 0 0 2 h 2 BD z dz 0 m (2.16) h 2 Thay công thức (2.15) và (2.16) vào công thức (2.14) ta được : NA ε (2.17) m - Moment uốn: M h   x 2 x M  zdz (2.18) y  y  h M  xy 2 xy  Hay: h 2 M σzdz (2.19) h 2 Với: T M MMMx y xy (2.20) 5
6. T σ x  y  xy (2.21) Thay công thức (2.8) vào công thức (2.19) ta được: h h h h h 2 2 2 2 2 MDD () zε zdz D  zdz z2 D ε dz  D zdz ε z 2 D dz mmmb mm mbmm b m h h h h h 2 2 2 2 2 (2.22) Đặt: h 2 CD z dz 0 (2.23) m h 2 h 1v 0 2 Eh3 DD z2 dz v 1 0 (2.24) m 2 h 12(1 v ) 1 v 2 0 0 2 Thay công thức (2.23) và (2.24) vào công thức (2.22) ta được : MD ε (2.25) b - Lực cắt ngang: h 2 Qx   xz   k  dz (2.26) Qy h  yz   2 Với là k 5 / 6 hệ số điều chỉnh cắt cho vật liệu đẳng hướng. Hay h 2 Q k τ dz (2.27) h 2 Với: T Q QQx y (2.28) T τ xz  yz (2.29) Thay công thức (2.9) vào công thức (2.27) ta được: h h 2 2 QD kε dz k ε D dz (2.30) s s s s h h 2 2 Đặt : h 2 khE 1 0 ED k dz (2.31) s h 2(1 v ) 0 1 2 Thay công thức (2.31) vào công thức (2.30) ta có: 6
7. QE εs (2.32) Từ các công thức (2.17), (2.25) và (2.32) ta có phương trình thể hiện mối quan hệ giữa nội lực và bi ế n d ạ ng mặt trung bình thông qua đ ịnh luậ t Hooke như sau : NA  0 0 εm  MD 0 0 ε  b  (2.33) QE  0 0 εs  3. Công thức phần từ hữu hạn trơn ES-MITC3. 3.1 Công thức PTHH vỏ phẳng MITC3 trong hệ trục tọa độ cục bộ Xét phần tử tam giác 3 nút trong hệ tọa độ cục bộ, trường chuyển vị tổng quát u của phần tử được nội suy bằng cách sử dụng các điểm chuyển vị tại các nút của phần tử với hàm dạng tuyến tính như sau : d1  T u u,,,, v w x  y  N1 N 2 N 3 d 2  (3.1) d3  T Trong đó, di u0 i,,,, v 0 i w 0 i xi  yi là chuyển vị tổng quát của điểm tại nút thứ i, Ni ()x là ma trận chéo của hàm dạng được xác định bởi: Ni diagNx i()()()()() Nx i Nx i Nx i Nx i  (3.2) Trong đó, N ()x là hàm dạng; trong hệ trục tọa độ tự nhiên ( r, s): i N1 1 r s N2 r (3.3) N3 s Thay biểu thức (3.1) vào biểu thức (2.3), biến dạng màng ε và biến dạng uốn ε và m b biến dạng cắt ε được viết lại : s d1 ε B d B,, B B d (3.4) m m m1 m 2 m 3 2 d3 d1 ε B d B,, B B d (3.5) b b b1 b 2 b 3 2 d3 d1 ε B d B,, B B d (3.6) s s s1 s 2 s 3 2 d3 Với : Ni, x 0 0 0 0 (3.7) Bmi 0N i,y 0 0 0 NNi,y i , x 0 0 0 7
8. 0 0 0 0 Ni, x (3.8) Bbi 0 0 0 N i,y 0 0 0 0 NNi, x i ,y 0 0NNi, x i 0 Bsi (3.9) 0 0NNi,y 0 i Với dấu , là phép lấy đạo hàm. Nếu chỉ dừng lại ở việc dùng phương trình (3.6) để tìm các chuyển vị thì bài toán chỉ đúng với những tấm có bề dày lớn hoặc vừa phải, còn đối với tấm mỏng sẽ dẫn đến hiện tượng “khóa cắt” (Shear locking) dẫn đến kết quả không còn chính xác nữa. Để khắc phục hiện tượng này cần xấp xỉ lại các biến dạng cắt bằng hàm mới thông qua những “điểm buộc” (tying point) [23]. Theo Lee và Bathe [23] các biến dạng trượt ngoài mặt phẳng trong hệ tọa độ tự nhiên được xấp xỉ lại như sau: 1 e () e  e  (3.10) qt2 st rt ert a1 b 1 r c 1 s (3.11) est a2 b 2 r c 2 s (3.12) 1 1 eqt ee  st  rt abrcs2 2 2 abrcs 1 1 1  (3.13) 2 2 Hình 3.1 Vị trí các điểm “tying point” cho phần tử tam giác 3 nút. Các điều kiện ràng buộc: (1) (1) ert(0,0) e rt ; ert(1,0) e rt (2) (2) est(0,0) e st ; est(0,1) e st (3)1 (3) (3) eqt(1,0) e qt e st e rt (3.14) 2 (3)1 (3) (3) eqt(0,1) e qt e st e rt 2 Từ sáu điều kiện trên ta có: (1) (2) (1) (3) (3) a1 ert ; b1 0; c1 est e rt e st e rt . 8
9. (2) a2 ert ; b2 c 1 ; c2 0 . (3.15) Biến dạng trượt có thể được viết lại: (1) ert e rt cs (2) est e st cr (3.16) (2) (1) (3) (3) Trong đó: c est e rt e st e rt Do đó, biến dạng cắt ngoài mặt phẳng được biểu diễn lại theo chuyển vị nút trong hệ tọa độ cục bộ như sau: 3 MITC3 MITC 3 εs  B sI d I (3.17) I 1 Với tọa độ nút phần tử trong hệ trục tọa độ cục bộ lần lượt là (,)x1 y 1 , (,)x2 y 2 , (,)x3 y 3 và định nghĩa a () x2 x 1 ,b () y2 y 1 , c () y3 y 1 , d () x3 x 1 các ma trận quan hệ giữa biến dạng và chuyển vị nút có thể viết lại dưới dạng tường minh như sau: b c 0 0 0 0 0 1 B 0d a 0 0 0 0 m1 2Ae d a b c 0 0 0 0 c 0 0 0 0 0 1 B 0 d 0 0 0 0 (3.18) m2 2Ae d c 0 0 0 0 b 0 0 0 0 0 1 B 0a 0 0 0 0 m3 2Ae a b 0 0 0 0 0 0 0b c 0 0 1 B 0 0 0 0d a 0 b1 2Ae 0 0 0d a b c 0 0 0 0c 0 0 1 B 0 0 0 0 d 0 (3.19) b2 2Ae 0 0 0 d c 0 0 0 0 b 0 0 1 B 0 0 0 0a 0 b3 2Ae 0 0 0a b 0 9
10. MITC3 1 0 0bc bcbc 6 Adabce 6 0 Bs1 e 2A 0 0daAbcad e 6 daad 6 0 MITC3 1 0 0c bc 2 cbc 6 ac 2 dbc 6 0 Bs2 e (3.20) 2A 0 0 d bd 2 c a d 6 ad 2 d a d 6 0 MITC3 1 0 0 bbc 2 bbc 6 bd 2 abc 6 0 Bs3 e 2A 0 0a ac 2 b a d 6 ad 2 a a d 6 0 Thế các quan hệ giữa biến dạng và chuyển vị, ứng suất và biến dạng vào phương trình công ảo và thực hiện các bước của phương pháp PTHH ta xác định được ma trận độ cứng của phần tử trong hệ trục tọa độ cục bộ K e, loc với thành phần độ cứng liên quan các bậc tự do của nút I và J. eloc, T e T e MITC 3T MITC 3 e KBDBBDBBDBIJ mI m mJAAA bI b bJ sI s sJ (3.21) Trong đó: D cho bởi (2.24), D không sử dụng hệ số ổn định và có công thức: b s 1 0 E 1 0 Eh D và D  1 0 (3.22) s m 2 2 1  0 1 1  0 0 1  2 Tương tự phần tử tấm MITC3, các ma trận quan hệ giữa biến dạng và chuyển vị của phần tử vỏ phẳng tam giác 3 nút MITC3 là hằng số. Do đó, ma trận độ cứng phần tử (3.21) được viết dưới dạng tường minh mà không cần sử dụng phép tính tích phân Gauss. 3.2 Công thức PTHH vỏ phẳng MITC3 trong hệ trục tọa độ toàn cục Gọi dI UVW I,,,,, I I XI  YI  ZI  là chuyển vị của nút I trong hệ tọa độ toàn cục OXYZ. Quan hệ giữa chuyển vị nút dloc trong hệ trục tọa độ cục bộ oxyz và chuyển vị nút I dI trong hệ trục tọa độ toàn cục OXYZ cho bởi: dloc Td (3.23) II Trong đó: nxX n xY n xZ 0 0 0 n n n 0 0 0 yX yY yX nzX n zY n zZ 0 0 0 T (3.24) 0 0 0 nxX n xY n xZ 0 0 0 n n n yX yY yZ 0 0 0 nzX n zY n zZ với (,,)n n n , (,,)n n n và (,,)n n n lần lượt là cosin chỉ phương của xX xY xZ yX yY yZ xZ zY zZ phương x,, y z trong hệ trục tọa độ toàn cục OXYZ. Thế (3.23) vào (3.17) ta được quan hệ giữa biến dạng và chuyển vị nút phần tử trong hệ trục tọa độ toàn cục: 3 3 3 ε B Td;; ε B Td εMITC3 B MITC 3 Td (3.25) m mI I b  bI I s  sI I III 1 1 1 10
11. Do đó, ma trận độ cứng của phần tử K e trong hệ trục tọa độ toàn cục với thành phần J độ cứng liên quan các bậc tự do của nút I và e TT eTT eTMITC3T * MITC 3 e KTBDBTTBDBTTBDBTIJ I mI m mJ JAAA I bI b bJ J I sI s sJ J (3.26) với TI và TJ là ma trận chuyển trục tọa độ (2.58) tính tại nút I và J . Từ (3.18), (3.19) và (3.20), ta nhận thấy rằng thành phần của ma trận độ cứng phần e tử K liên quan bậc tự do  ZI bằng không. Vì vậy, trong ma trận độ cứng kết cấu, thành phần liên quan đến bậc tự do  ZI tại nút kết nối giữa các phần tử đồng phẳng sẽ bằng không, tức là ma trận độ cứng kết cấu bị suy biến. Để khắc phục hiện tượng này, tại vị trí 3 liên quan bậc tự do  ZI của ma trận độ cứng phần tử, một giá trị bằng 10 lần giá trị lớn nhất của các thành phần trên đường chéo chính của ma trận K e được thêm vào [13]. 3.3 Công thức phần tử hữu hạn trơn ES-FEM với phần tử MITC3 Xét kết cấu vỏ có mặt trung bình  được rời rạc bằng các phần tử vỏ phẳng 3 nút. Tương tự phương pháp làm trơn trên cạnh cho kết cấu tấm, miền làm trơn k trên cạnh k của kết cấu vỏ cũng được xác định bằng cách nối 2 nút của cạnh k với trọng tâm của các phần tử chung cạnh k như Hình 3.2 [12, 15, 16]. Các miền làm trơn này phải thỏa mãn điều Ned k i j kiện k 1  ,     , với Ned là số cạnh của tất cả phần tử thuộc miền kết cấu . (a) (b) Hình 3.2 (a) Kết cấu vỏ  được rời rạc bằng các phần tử vỏ tam giác 3 nút và các miền làm trơn k trên cạnh phần tử; (b) Hệ trục tọa độ cục bộ oxyz được định nghĩa trên mặt phẳng tiếp xúc miền làm trơn k tại cạnh Các trường biến dạng trơn được xác định trên miền làm trơn k 1 1 1 εk εd  ; ε k ε d  ; ε  k, MITC 3 ε MITC 3 d  (3.27) mAAAk k m b k  k b s k  k s với Ak là diện tích của miền làm trơn k. 11
12. Do miền làm trơn k của cạnh k không nằm trên biên kết cấu gồm 2 phần tử vỏ không đồng phẳng nên các trường biến dạng được định nghĩa trong hệ trục tọa độ cục bộ oxyz của mỗi phần tử được biến đổi thành trường biến dạng trong hệ trục tọa độ cục bộ oxyz oxyz ox chung của 2 phần tử. Hệ trục tọa độ được định nghĩa bởi trục theo phương cạnh chung của 2 phần tử, trục oz có véc-tơ chỉ phương là trung bình của 2 véc-tơ pháp tuyến của 2 phần tử, và trục oy tạo với trục ox và trục oz thành 1 tam diện thuận như Hình 3.2b. Sử dụng ma trận biến đổi quan hệ giữa ứng suất và biến dạng từ hệ trục cục bộ sang hệ trục tọa độ toàn cục cho bởi [23], ta xác định được quan hệ giữa biến dạng e e e, MITC 3 k k k, MITC 3 εm, ε b , ε s trong hệ trục cục bộ phần tử oxyz và biến dạng εm, ε b , ε s trong hệ oxyz k trục tọa độ cục bộ của miền làm trơn  như sau εk Q e ε e;; ε  k Q e ε e ε k, MITC 3 Q e ε e , MITC 3 (3.28) m m m b b b s s s Trong đó, 2 2 nxx n yx 2 n xx  n yx  2nzx n xx  2nyx n zx  2 2 n n2 n n 2n n 2 n n xy yy  xy  yy  zy xy  yy  zy  2 2 n n2 n n 2nzz n xz  2 n yz  n zz  QQQe e xz yz  xz  yz  ; e (3.29) m b s n n n n n n n n nxx n xy  n yx  n yy  n xx  n yy  n yx  n xy  zx xy  xx  zy  yx  zy  zx  yy  n n n n n n n n n n n n n n n n xz xx yz  yx  xz  yx  yz  xx  zy xz xy  zz  yy  zz  zy  yz  n n n n n n n n nxy n xz n yy  n yz  n xy  n yz  n yy  n xz  zz xx xz  zx  yz  zx  zz  yx  với nxx,,,,,,,, n xy  n xz n yx  n yy  n yz  n zx  n zy  n zz  lần lượt là cosin chỉ phương của phương x, y, z trong hệ trục tọa độ oxyz . Áp dụng (3.28), các trường biến dạng trơn (3.27) có thể được viết lại k k k 1NNNeAAAe 1 e e 1 e e 3 k  e e e e e εm k ε m k  Q m ε m k  Q m  B mI T I d I AAAe 13 e 1 3 e 1 3 I 1 k k k 1NNNeAAAe 1 e e 1 e e 3 εk ε e Q e ε e Q e B e B e T d (3.30) bk b k  b b k  b b  bI I I AAAe 13 e 1 3 e 1 3 I 1 k k k 1Ne AAe 1 NNe e 1 e Ae 3 k, MITC 3  e , MITC 3 e e, MITC 3 e e , MITC 3 εs k ε s k Qsε s k  Q s  B sI T I d I AAe 1 3e 1 3 A e 13 I 1 Trong công thức (3.30), chỉ số e tương ứng với giá trị tính toán của phần tử e thuộc k k k miền làm trơn  ; N e 1 nếu cạnh k thuộc biên kết cấu và Ne 2 nếu cạnh k không thuộc biên kết cấu (xem Hình 3.2a). Để áp dụng nguyên lý công ảo xác định phương trình cân bằng theo PPPTHH, các trường biến dạng trơn định nghĩa trong hệ trục tọa độ cục bộ oxyz được biến đổi thành các trường biến dạng trong hệ trục tọa độ toàn cục OXYZ [15, 16] như sau: ε Q εk;; ε Q  ε k ε MITC3 Q  ε k , MITC 3 (3.31) m m m b b b s s s Trong đó, 12
13. 2 2 2 nxX nyX  n zX2 n xX  nyX  2 nzX  nxX  2 nyX  nzX    2 2 2 QQm b n xY nyY  n zY2 nxY  nyY  2 nzY  nxY  2 n yY nzY  nnxX xY nn yX yY nn zXzY  nn xX  yY nn yX xY nnzX xY  nn xXzY  nn yX zY nnzX  yY  (3.32) n n n n n n nn nn nn nn nn nn  xZ xX yZ yX  zZ zX xZyX yZ xX zZ xX xZ zX yZ zX  zZ yX  Qs nnxY xZ nn yY yZ nn zYzZ  nnxY  yZ nn yY xZ nnzY  xZ  nn xYzZ  nnyY zZ  nn zY yZ với nxX,,,,,,,, n xY n xZ n yX n yY n yZ n zX n zY n zZ lần lượt là cosin chỉ phương của phương x,, y z trong hệ trục tọa độ toàn cục OXYZ. Từ (3.30) và (3.31), quan hệ giữa biến dạng và chuyển vị trong hệ tọa độ toàn cục được biểu diễn như sau k k 1 Ne Ae 3 N ε Q ε k Q  Q e B e T d B d m m mk  m m  mI I I  mI I A e 13 I 1 I 1 k k 1 Ne Ae 3 N ε Q ε k Q  Q e B e T d B d (3.33) b b bk  b b  bI I I  bI I A e 13 I 1 I 1 k k 1 Ne Ae 3 N εMITC3 Q ε k , MITC 3 Q  Q e B e , MITC 3 T d B MITC 3d s s sk  s s  sI I I  sI I A e 13 I 1 I 1 k k    MITC 3 Ở đây, N là số nút của các phần tử thuộc miền làm trơn  và BBBmI,, bI sI là các ma trận quan hệ giữa biến dạng và chuyển vị tại nút I của phần tử e thuộc miền làm trơn k k 1 Ne Ae BQQBT  e e mIk  m m mI I A e 1 3 k 1 Ne Ae BQQBT  e e (3.34) bIk  b b bI I A e 1 3 k 1 Ne Ae BQQBT MITC3  e e , MITC 3 sIk  s s sI A e 1 3 Do đó, ma trận độ cứng K của kết cấu vỏ được lắp ghép từ các ma trận độ cứng K k k  k xây dựng trên miền làm trơn  có thành phần độ cứng KIJ liên quan bậc tự do của nút I, J k xác định bởi: T KBDBBDBBDBk  T AAA k  T  k  MITC3 *  MITC 3 k (3.35) IJ mI m mJ bI b bJ sI s sJ 4. Các ví dụ số 4.1 Kiểm tra điều kiện patch test Trong ví dụ này, khả năng tính toán được nội lực không đổi của phần tử tấm ES- MITC3 được kiểm tra thông qua bài toán tấm hình chữ nhật có kích thước và chia lưới phần tử như Hình 4.1. Tấm có chiều dày h 0,01m và làm bằng vật liệu có mô-đun đàn hồi E 105 KN/m2, hệ số Poisson  0,25. Tấm chịu chuyển vị cưỡng bức bởi phương trình độ võng w (1 x 2 y x2 xy y 2 ) / 2. Kết quả trong Bảng 4.1 cho thấy phần tử 13
14. ES-MITC3 có khả năng tính toán được mô-men uốn MMx, y và mô-men xoắn M xy không đổi, tức là vượt qua được điều kiện patch test. Hình 4.1 Hình học và tọa độ nút của lưới phần tử trong bài toán patch test. Bảng 4.1: Kết quả chuyển vị và nội lực tại nút 5 của bài toán patch test. w   M M M Phần tử 5 x5 y5 x5 y5 xy5 ES-MITC3 0,6422 1,1300 -0,6400 -0,0111 -0,0111 -0,0033 Chính xác 0,6422 1,1300 -0,6400 -0,0111 -0,0111 -0,0033 4.2 Kiểm tra điều kiện không phụ thuộc vào thứ tự đánh số nút phần tử Cho vỏ hình tam giác dày h 0, 02 m và có tọa độ các đỉnh trong không gian lần lượt là (0, 0, 0) , (9,10, 4) và (5,8, 2) . Vỏ làm bằng vật liệu có mô-đun đàn hồi E 3 x 10 6 KN/m2, hệ số Poisson  0,3 . Vỏ chịu tải trọng và mô-men tập trung T -2 T []PPPMMMx y z x y z = 10 x[10 30 20 5 8 7] tại đỉnh có tọa độ (0, 0, 0) và ngàm tại 2 đỉnh còn lại. Để kiểm tra khả năng tính toán của phần tử không phụ thuộc vào thứ tự đánh số nút, kết cấu vỏ đã cho được rời rạc bằng 1 phần tử vỏ ES-MITC3 được đánh số nút theo 3 sơ đồ khác nhau như Hình 4.2. Kết quả chuyển vị như nhau tại điểm (0, 0, 0) trong 3 sơ đồ đánh số nút phần tử được cho Bảng 4.2. Vì vậy, ứng xử của phần tử vỏ ES-MITC3 không phụ thuộc vào thứ tự đánh số nút phần tử. (a) (b) (c) Hình 4.2 Hình học, điều kiện biên và 3 sơ đồ đánh số nút phần tử vỏ ES-MITC3 14
15. Bảng 4.2: Kết quả chuyển vị tại tọa độ (0, 0, 0) trong các sơ đồ đánh số nút khác nhau Sơ đồ u u u    đánh số nút X Y Z X Y Z (a) -0,5959 0,4354 -0,2521 0,3893 -0,1462 0,0007 (b) -0,5959 0,4354 -0,2521 0,3893 -0,1462 0,0007 (c) -0,5959 0,4354 -0,2521 0,3893 -0,1462 0,0007 4.3 Phân tích tĩnh vỏ trụ liên kết ngàm chịu tải trọng tập trung Vỏ trụ với hai đầu liên kết ngàm phẳng (u 0, w 0) chịu tải trọng tập trung P 1 0 0 N như Hình 4.3. Các thông số hình học và vật liệu như sau như sau: R 300in, L 600 in, t 3in, E 3 x 10 6 N/in2. Vì tính chất đối xứng qua 2 trục của bài toán, chỉ 1/8 vỏ được mô hình với hệ lưới 8 x 8, 12 x 12, 16 x 16, 20 x 20 và 24 x 24 phần tử. Với các loại lưới NxN phần tử khác nhau, Bảng 4.3 trình bày kết quả độ võng tại vị trí đặt lực tập trung khi phân tích bằng các phần tử vỏ phẳng tam giác 3 nút ES-MITC3, MITC3+, CS-DGS3, ES-DSG3, phần tử vỏ tứ giác 4 nút MITC4 và lời giải tham khảo [17] và được so sánh ở biểu đồ trong Hình 4.5. Hình 4.3 Vỏ trụ liên kết ngàm ở hai đầu chịu tải trọng tập trung. Hình 4.4 Chia lưới 8x8 phần tử tam giác. 15
16. Bảng 4.3: So sánh kết quả chuyển vị w (10-5) tại điểm đặt lực của vỏ. Kiểu Số phần tử Kết quả phần tử 8x8 12x12 16x16 20x20 24x24 chính xác MITC3+ 1,3945 1,6097 1,6981 1,7433 1,7697 ES-DSG3 1,1132 1,4689 1,6386 1,7295 1,7825 CS-DSG3 1,1358 1,4299 1,5768 1,6586 1,7081 1,8248 MITC4 1,3546 1,5922 1,6922 1,7425 1,7713 ES-MITC3 1,3535 1,6318 1,7486 1,8043 1,8338 Hình 4.5 Biểu đồ so sánh chuyển vị tại điểm đặt lực của vỏ với các kiểu phần tử khác nhau. Từ Bảng 4.3 và biểu đồ Hình 4.5 nhận thấy kết quả về chuyển vị tại điểm đặt lực của vỏ trụ đối với phần tử ES-MITC3 cho kết quả hội tụ đến lời giải tham khảo khá tốt, thậm chí hơn cả phần tử tứ giác. Từ đó cho thấy phần tử ES-MITC3 cho kết quả tốt khi phân tích kết cấu vỏ trụ chịu tải trọng tập trung. 16
17. 4.4 Phân tích tĩnh mái vòm ngàm hai cạnh tải phân bố đều Mái vòm liên kết ngàm phẳng hai đầu u w 0 chịu trọng lượng bản thân q 90 KN như Hình 4.6. Các thông số hình học và vật liệu như sau: R 25m, L 50m, 8 2 t 0, 25 m, E 4,32 x 10 KN/m ,  0. Vì tính chất đối xứng qua 2 trục của bài toán, chỉ 1/4 vỏ được mô hình với hệ lưới 8 x 8, 12 x 12, 16 x 16, 20 x 20 và 24 x 24 phần tử. Kết quả phân tích chuyển vị tại điểm A của các phần tử được trình bày trong Bảng 4.4 và được so sánh ở biểu đồ trong Hình 4.8. Kết quả chính xác của chuyển vị tại điểm A được tham khảo từ [18] là 0,3024. Hình 4.6 Mái vòm liên kết ngàm ở hai đầu chịu tải trọng bản thân. Hình 4.7 Chia lưới 8x8 phần tử tam giác. 17
18. Bảng 4.4: So sánh kết quả chuyển vị lớn tại điểm A của vỏ. Kiểu Số phần tử Kết quả phần tử 8x8 12x12 16x16 20x20 24x24 chính xác MITC3+ 0,2633 0,2822 0,2901 0,2941 0,2963 ES-DSG3 0,2943 0,2976 0,2981 0,2983 0,2984 CS-DSG3 0,2578 0,2793 0,2885 0,2933 0,2960 0,3024 MITC4 0,2947 0,2988 0,3005 0,3014 0,3019 ES-MITC3 0,2935 0,2966 0,2974 0,2976 0,2978 Hình 4.8 Biểu đồ so sánh chuyển vị tại điểm A của vỏ với các kiểu phần tử khác nhau. Từ Bảng 4.4 và biểu đồ Hình 4.8 nhận thấy kết quả về chuyển vị tại điểm A của mái vòm so sánh với phần tử MITC3+, phương pháp làm trơn trên cạnh đã giúp phần tử vỏ phẳng ES-MITC3 cho kết quả chuyển vị tại A tốt hơn. Cùng kỹ thuật làm trơn trên cạnh, phần tử vỏ phẳng ES-MITC3 với phương pháp khử “khóa cắt” MITC3 cho kết quả tiệm cận với lời giải tham khảo nhanh hơn phần tử vỏ phẳng ES-DGS3 dùng phương pháp khử “khóa cắt” DSG3. Tuy nhiên, kết quả tính toán bởi phần tử ES-MITC3 không tốt bằng phần tử vỏ phẳng dùng phương pháp làm trơn trên phần tử CS-DSG3 và đặc biệt là phần tử vỏ phẳng tứ giác 4 nút MITC4. Từ đó cho thấy phần tử ES-MITC3 cho kết quả khá tốt khi phân tích kết cấu mái vòm chịu tải trọng bản thân. 18
19. 4.5 Phân tích tĩnh vỏ panel cầu chịu tải trọng tập trung Vỏ panel cầu liên kết tựa đơn chịu tải trọng tập trung P 45,4 kgf như Hình 4.3. Các thông số hình học và vật liệu như sau như sau: a b 400 mm, R 2400mm, h 2,54 mm, E 7037kgf và  0.3 . Vì tính chất đối xứng qua 2 trục của bài toán, chỉ 1/4 vỏ được mô hình với hệ lưới 8 x 8, 12 x 12, 16 x 16, 20 x 20 và 24 x 24 phần tử. Kết quả phân tích chuyển vị tại điểm đặt lực của các phần tử được trình bày trong Bảng 4.5 và được so sánh ở biểu đồ trong Hình 4.11. Kết quả chính xác của chuyển vị được tham khảo [19] là 1. Hình 4.9 Vỏ panel cầu liên kết tựa đơn chịu tải trọng tập trung. Hình 4.10 Chia lưới 8x8 phần tử tam giác. 19
20. Bảng 4.5: So sánh kết quả chuyển vị w tại điểm đặt lực của vỏ. Kiểu Số phần tử Kết quả phần tử 8x8 12x12 16x16 20x20 24x24 chính xác MITC3+ 0,8942 0,9406 0,9609 0,9716 0,9781 ES-DSG3 0,5501 0,7295 0,8224 0,8695 0,8947 CS-DSG3 0,6433 0,8134 0,8875 0,9237 0,9439 1,0000 MITC4 0,8785 0,9382 0,9616 0,9733 0,9799 ES-MITC3 0,7359 0,8478 0,8881 0,907 0,9174 Hình 4.11 Biểu đồ so sánh chuyển vị tại điểm đặt lực của vỏ panel cầu với các kiểu phần tử khác nhau. Từ Bảng 4.5 và biểu đồ Hình 4.11 nhận thấy kết quả về chuyển vị tại điểm đặt lực của vỏ trụ đối với phần tử ES-MITC3 so sánh với lời giải tham khảo, kết quả chuyển vị tính bởi phần tử vỏ phẳng ES-MITC3 chỉ đạt 90% khi dùng 24x24 phần tử, kém xa kết quả tính toán bởi phần tử vỏ phẳng không dùng phương pháp làm trơn MITC3+, phần tử được làm trơn trên miền phần tử CS-DSG3 hoặc phần tử vỏ phẳng tứ giác 4 nút MITC4. Đây là kết cấu có độ cong theo 2 phương nên việc làm trơn trên cạnh có thể dẫn tới sai số do trung bình biến dạng trên miền làm trơn là mặt tiếp xúc của 2 phần tử vỏ không đồng phẳng tại cạnh chung. Tuy nhiên, so với phần tử cùng sử dụng phương pháp làm trơn trên cạnh, phần tử ES-MITC3 vẫn cho kết quả tốt hơn phần tử ES-DSG3, trong cả hai trường hợp có và không có sử dụng hệ số ổn định, tuy hội tụ kém hơn phần tử MITC4 nhưng lại hội tụ nhanh hơn các phần tử khác và kết quả thu được tiệm cận với kết quả chính xác. 20
21. 5. Kết luận Qua kết quả phân tích và so sánh trong các ví dụ số của chương 4, ta có thể rút ra một số kết luận cho luận văn như sau: - Trên cơ sở lý thuyết biến dạng cắt bậc nhất kết hợp kỹ thuật khử khóa cắt của phần tử MITC3, cùng kỹ thuật làm trơn trên cạnh ES-FEM luận văn đã xây dựng phần tử ES- MITC3 dùng phân tích tĩnh kết cấu vỏ. Kết quả số cho thấy tùy bài toán kết quả có thể tốt hoặc không tốt hơn các phần tử khác nhưng kết quả cũng tiệm cận với lời giải chính xác. - Đối với phần lớn các bài toán thì kết quả sẽ hội tụ nhanh khi tăng hệ lưới. - Phần tử ES-MITC3 giải quyết khá tốt cho kết cấu cong 1 phương, tuy nhiên lại cho kết không tốt bằng đối với kết cấu có độ cong theo 2 phương do việc làm trơn trên cạnh có thể dẫn tới sai số do trung bình biến dạng trên miền làm trơn là mặt tiếp xúc của 2 phần tử vỏ không đồng phẳng tại cạnh chung. Lời cảm ơn. Tôi xin chân thành cảm ơn thầy hướng dẫn khoa học là PGS.TS Nguyễn Văn Hiếu đã đưa ra những gợi ý đầu tiên để hình thành nên ý tưởng của đề tài và hướng dẫn tôi một cách sâu sát để giúp tôi nhận định đúng đắn trong những vấn đề nghiên cứu mà quan trọng nhất là sự trung thực trong làm nghiên cứu khoa học. Thầy cũng đã hướng dẫn tôi cách tiếp cận nghiên cứu hiệu quả cũng như những nguồn tài liệu quý báu. Và với sự hướng dẫn khoa học, nghiêm túc, tận tình đó của thầy đã giúp tôi đạt đến kết quả nghiên cứu cuối cùng. Tài liệu tham khảo [1] E. Rank, R. Krause, K. Preusch. On the accuracy of p-version elements for the Reissner–Mindlin plate problem. Int. J. Numer. Methods Eng, 43, pp.51–67, 1988. [2] Ahmad et al. Analysis of thick and thin shell structures by curved finite elements. Int. J. Num. Meth. Engrg, 2, pp.419–451, 1970. [3] Zienkiewicz, Olgierd Cecil, et al. The finite element method. Vol. 3. London: McGraw- hill, 1977. [4] Reissner, E. The effect of transverse shear deformation on the bending of elastic plates. J. Appl. Mech., 12, pp.69–76, 1945. [5] Mindlin, R.D. Influence of rotatory inertia and shear in flexural motions of isotropic elastic plates. J. Appl. Mech., 18, pp.31–38, 1951. [6] Dvorkin EN, Bathe KJ. A continuum mechanics based four-node shell element for general nonlinear analysis. Engineering Computations, 1(1), pp.77–88, December 1984. [7] Bathe KJ, Dvorkin EN. A formulation of general shell elements – the use of mixed interpolation of tensorial components. Int J Numer Meth Eng, 22, pp.697–722, 1986. [8] Bucalem ML, Bathe KJ. Higher-order MITC general shell elements. Int J Numer Meth Eng, 36:37, 29–54, 1993. [9] Phill-Seung Lee, Klaus-Jurgen Bathe. Development of MITC isotropic triangular shell finite elements. Computers & Structures, 2014. [10] Do-Nyun Kim, Klaus-Jurgen Bathe. A triangular six-node shell element. Computers & Structures, 2019. [11] Youngyu Lee, Phill-Seung Lee, Klaus-Jurgen Bathe. The MITC3+ shell element and its performance. Computers & Structures, 2014. [12] Liu, G.R., Dai, K.Y., Nguyen-Thoi T A smoothed finite element method for mechanics problems. Computational Mechanics, 39, pp.859–877, 2007. 21
22. [13] Nguyen Thoi T, Phung Van P, Nguyen Xuan H. A cell - based smoothed discrete shear gap method (CS - DSG3) using triangular elements for static and free vibration analyses of shell structures. Computers & Structures, 2013. [14] Cheng L and Liu GR. Adaptive analysis using the edge-based smoothed finite element method; (submitted), 2008. [15] Xiangyang Cui, Guo-Rong Liu, Guang-yao Li, GuiYong Zhang, Gang Zheng. Analysis of plates and shells using an edge-based smoothed finite element method. Comput Mech, 2010. [16] Son Nguyen-Hoang, Phuc Phung-Van, Sundararajan Natarajan, Hyun-Gyu Kim. A combined scheme of edge-based and node-based smoothed fiite element methods for Reissner–Mindlin flt shells. Engineering with Computers, 2015. [17] Fluge, W Stress in shells. Berlin: Springer,1960, pp.2-34. [18] Macneal, R.H., and Harder, R.L. . A proposed standard set of problems to test finite element accuracy. Finite elements in analysis and design, 1, pp.3-20,1985. [19] Mousa, A.I., and El Naggar, M.H Shallow Spherical Shell Rectangular Finite Element for Analysis of Cross Shape Shell Roof Electronic, Journal of Structure Engineering, 7, pp.41-51, 2007. Thông tin liên hệ tác giả: Họ tên: Quách Văn Nghiệp Đơn vị: Học viên cao học trường ĐHSPKT TPHCM Điện thoại: (+84) 0983777076 Email: vannghiepquach@gmail.com 22
23. BÀI BÁO KHOA HỌC THỰC HIỆN CÔNG BỐ THEO QUY CHẾ ĐÀO TẠO THẠC SỸ Bài báo khoa học của học viên có xác nhận và đề xuất cho đăng của Giảng viên hướng dẫn Bản tiếng Việt ©, TRƯỜNG ĐẠI HỌC SƯ PHẠM KỸ THUẬT TP. HỒ CHÍ MINH và TÁC GIẢ Bản quyền tác phẩm đã được bảo hộ bởi Luật xuất bản và Luật Sở hữu trí tuệ Việt Nam. Nghiêm cấm mọi hình thức xuất bản, sao chụp, phát tán nội dung khi chưa có sự đồng ý của tác giả và Trường Đại học Sư phạm Kỹ thuật TP. Hồ Chí Minh. ĐỂ CÓ BÀI BÁO KHOA HỌC TỐT, CẦN CHUNG TAY BẢO VỆ TÁC QUYỀN! Thực hiện theo MTCL & KHTHMTCL Năm học 2017-2018 của Thư viện Trường Đại học Sư phạm Kỹ thuật Tp. Hồ Chí Minh.