Phân tích ứng xử tấm composite bằng phương pháp phần tử hữu hạn tam giác 3 nút MITC3, sử dụng thuyết biến dạng cắt bậc cao

Nội dung text: Phân tích ứng xử tấm composite bằng phương pháp phần tử hữu hạn tam giác 3 nút MITC3, sử dụng thuyết biến dạng cắt bậc cao

1. PHÂN TÍCH ỨNG XỬ TẤM COMPOSITE BẰNG PHƯƠNG PHÁP PHẦN TỬ HỮU HẠN TAM GIÁC 3 NÚT MITC3, SỬ DỤNG THUYẾT BIẾN DẠNG CẮT BẬC CAO. ANALYSIC OF LAMINATED COMPOSITE PLATES USING THE HIGH-ORDER SHEAR DEFORMATION THEORY BY AN MITC3 ELEMENT (1) (1) Nguyễn Hòa , Châu Đình Thành (1) Khoa xây dựng và cơ học ứng dụng, trường đại học Sư phạm kỹ thuật thành phố Hồ Chí Minh. Từ khóa: Finite element Tóm tắt method (FEM), Edge-based Trong bài báo này, phương pháp phần tử hữu hạn trơn dựa trên cạnh smoothed FEM (ES-FEM), ES-FEM (the edge-based smoothed finite element method) được phát MITC3, C0-type HSDT, triển cho bài toán phân tích tĩnh và dao động tự do của kết cấu tấm Laminated composite plates. composite lớp. Biến dạng trượt do lực cắt sẽ được kể tới theo lý thuyết biến dạng cắt bậc cao dựa trên hàm bậc 3 của Reddy. Miền hình học được rời rạc thành lưới phần tử tam giác ba nút với bảy bậc tự do cho mỗi nút. Trong phương pháp ES-FEM, ma trận độ cứng được tính toán bởi kỹ thuật trơn hóa biến dạng trên miền trơn (smoothing domains) dựa trên cạnh của phần tử. Để giải quyết hiện tượng khóa cắt “shear locking” khi tấm có chiều dày mỏng dần, kỹ thuật nội suy các thành phần ten xơ (mixed interpolation tensorial components viết tắt là MITC) được sử dụng kết hợp với kỹ thuật trơn hóa biến dạng và được gọi là phương pháp phần tử hữu hạn trơn dựa trên cạnh ES-MITC3. Tính hiệu quả và độ chính xác của phương pháp ES-MITC3 được kiểm chứng thông qua các ví dụ số phân tích cho bài toán tĩnh và dao động tự do của kết cấu tấm composite lớp. 1. Giới thiệu: Trong những năm gần đây, kết cấu tấm nhiều lớp luôn là một nhu cầu thiết yếu. Trong vài thập kỷ qua làm bằng vật liệu composite đã được sử dụng rộng rãi nhiều lý thuyết cho việc giải quyết ứng xử của tấm và chuyên sâu trong nhiều ứng dụng kỹ thuật như composite nhiều lớp đã được các nhà khoa học nghiên hàng không, hàng hải và cơ sở hạ tầng dân dụng, vv cứu và giới thiệu. Những lý thuyết này có thể được bởi vì chúng có nhiều thuận lợi về tính chất cơ học phân loại thành: lớp tương đương (ESL), zig-zag như độ cứng cao so với trọng lượng. Điều này đặc biệt (ZZ), và lớp thông minh (LW) [7]. Lần đầu tiên có ý nghĩa trong kết cấu không gian vũ trụ, tàu ngầm Noor [5, 6] đưa ra một lý thuyết đàn hồi ba chiều (3D) và trong lĩnh vực xây dựng các kết cấu cao tầng. Tuy để cải thiện tính chính xác của ứng suất cắt ngang. nhiên, để có được các điều kiện thuận lới này thường Trong lý thuyết đàn hồi 3D, mỗi lớp được mô phỏng đi kèm với sự phức tạp về phân tích, mô hình tính như một chất rắn 3D, và độ chính xác của ứng suất cắt toán, vv Để sử dụng các tấm ghép nhiều lớp có hiệu ngang được cải thiện một cách đáng kể. Tuy nhiên quả trong thực tiễn thì việc cần thiết là phải phát triển việc sử dụng lý thuyết này làm chi phí tính toán tăng các lý thuyết phân tích thích hợp [1, 2] nhằm dự đoán lên đáng kể [9]. Các lý thuyết tấm nhiều lớp cổ điển chính xác ứng xử của chúng dưới các dạng tải trọng (CLPT) dựa trên giả thuyết Kirchhoff cung cấp kết khác nhau. Từ nhiều thập kỷ qua phương pháp phần quả hợp lý cho tấm mỏng. Tuy nhiên, CLPT bỏ qua tử hữu hạn đã được xem như là phương pháp hiệu quả các hiệu ứng cắt ngang và do đó không thích hợp cho và chiếm ưu thế trong việc phân tích các kết cấu tấm, mô hình tấm dày vừa phải hoặc tấm dày trong đó tác vỏ nói chung và tấm composite nói riêng. Tuy nhiên dụng cắt ngang rõ rệt hơn. Để khắc phục những hạn hiệu quả của việc phân tích tính toán còn phụ thuộc chế của CLPT và kết hợp chính xác các hiệu ứng cắt vào nhiều yếu tố như mô hình toán học, lưới phần ngang, nhiều giả thuyết biến dạng cắt đã được phát tử , bên cạnh đó tốc độ hội tụ của bài toán nên được triển. Reissner và Mindlin đề xuất lý thuyết cắt biến tối ưu [3]. Vì vậy việc tìm kiếm những phương pháp dạng cắt bậc nhất (FSDT) trong đó ứng suất cắt ngang tính toán hiệu quả với độ tin cậy cao trong phân tích được giả định không đổi và do đó nó đòi hỏi việc sử 1
2. dụng hệ số điều chỉnh biến dạng cắt để đáp ứng các các phần tử tam giác chiếm ưu thế trong việc rời rạc điều kiện biên tự do tại các bề mặt trên và dưới của hóa hình học của các kết cấu phước tạp. tấm. Những hạn chế của FSDT có thể được khắc Do liên tục trong hình học (độ cong hay độ dày), phục bằng cách giới thiệu các lý thuyết biến dạng cắt không tương thích của điều kiện biên, và sự bất bậc cao (HSDT). Các HSDT đã được phát triển bởi thường của tải trọng, các ứng suất, chuyển vị của cấu Reddy [8], Matsunaga [11], Kant và Swaminathan trúc tấm vỏ thường thay đổi rất nhanh và gây ra sự tập [12] và Liu et al. [13], vv. Những mô hình này có thể trung năng lượng biến dạng. Nên để có được một giải bỏ qua yếu tố điều chỉnh biến dạng và cho ra ứng suất pháp mang lại độ chính xác như mong muốn như vậy, cắt ngang chính xác và ổn định hơn. Từ đó đến nay một hệ lưới tốt hoặc là chức năng nội suy bậc cao phải nhiều nghiên cứu về lý thuyết biến dạng cắt bậc cao được sử dụng. Gần đây Kim và Bathe [22] đã phát đã ra đời và có những cải tiến đáng kể [9], [23-29]. triển và nghiên cứu một phương pháp phần tử hữu hạn Năm 1970 Ahmad, Irons và Zienkiewicz [14, 15] đã hàm nội suy bậc cao, nhằm làm tăng độ chính xác của giới thiệu một phần tử với tham số C0 độc lập nội suy các phần tử mà không có bất cứ sự thay đổi nào trong cho chuyển vị và góc xoay. Phương diện phù hợp hệ lưới. Phương pháp này không những sử dụng các nhất của phần tử này là các hàm nội suy chỉ cần thỏa gradient cao hơn mà còn làm giảm ứng suất nhảy giữa điều kiện C0 và giới thiệu về sự ảnh hưởng của các các phần tử. [4] [19-22] [30-34]. biến dạng cắt. Phần tử này được biết đến như các Trong nỗ lực phát triển xa hơn phương pháp PTHH, phần tử vỏ Reissner / Mindlin [16, 17]. Mặc dù bao tác giả Liu và cộng sự đã ứng dụng kỹ thuật làm trơn gồm các biến dạng cắt để phân tích vỏ dày nhưng khó hóa biến dạng [36] để thiết lập công thức phương khăn chính của phần tử là hiện tượng “khóa cắt” pháp PTHH trơn dựa trên phần tử con (cell) còn được (shear locking) khi chiều dày tấm giảm dần. Vào gọi là SFEM hoặc CS-FEM cho bài toán 2D cơ vật những năm 1970 hầu hết các nghiên cứu trên lĩnh vực rắn và sau đó CS-FEM được phát triển cho tấm và vỏ. vỏ đều dựa trên phương pháp mà Ahmad, Irons và Bằng cách sử dụng các phần tử con (cell) trong mỗi Zienkiewicz đã xây dựng và tìm các biện pháp để phần tử (element) (ví dụ phần tử con 4 nút), CS-FEM khắc phục hiện tượng “khóa cắt”. Tuy nhiên với việc đã làm tăng mức độ chính xác của lời giải. Kỹ thuật đưa ra giả thuyết các dạng năng lượng bằng không đã làm trơn hóa biến dạng còn được kết hợp với phương làm giảm đi độ tin cậy trong các kết quả nghiên cứu pháp PTHH mở rộng (XFEM) để giải quyết bài toán của họ [18]. rạn nứt 2 chiều trong cơ học vật rắn liên tục và kết cấu Năm 1980 Bathe và Dvorkin đề xuất phương pháp tấm. Với cách tiếp cận khác của phương pháp PTHH nội suy hỗn hợp các thành phần ten xơ (mixed trơn là kết hợp với nút (Node-based Smoothed Finite interpolation tensorial components viết tắt là MITC) Element Method – NS-FEM) được áp dụng vào phân đã giải quyết được các vấn đề về khóa cắt. Các tích thích nghi (adaptive analysis). Sau đó bằng cách phương pháp MITC rất hiệu quả trong giải quyết các kết hợp NS-FEM và FEM với hệ số tỷ lệ, một phương bài toán tấm vỏ và cho kết quả tin cậy. Từ khi ra đời pháp mới được đặt tên là phương pháp PTHH alpha phương pháp MITC đã được sử dụng rất thành công (αFEM) được đề xuất và kết quả thu được năng lượng trong việc loại bỏ các hiện tượng “khóa cắt” cho bài biến dạng gần với lời giải chính xác, sử dụng phần tử toán tấm/vỏ với phần tử tứ giác. Kỹ thuật này được tam giác và phần tử tứ diện. Bathe và Dvorkin đề xuất với phần tử 4 nút và phần tử Năm 2008, tác giả Liu và cộng sự [37] đã đề xuất 8 nút (các MITC4 và MITC8) [19, 20]. Sau đó, Bathe phương pháp PTHH trơn với cách tiếp cận dựa trên và Bucalem đã phát triển nó lên 9 nút và 16 nút cạnh (the Edge-based Smoothed Finite Element (MITC9 và MITC16) [21]. Bên cạnh thành công của Method –ES-FEM) cho bài toán phân tích tĩnh, dao các phần tử tứ giác thì việc xây dựng các phần tử tam động tự do và dao động cưỡng bức áp dụng cho bài giác cũng đã và đang được nghiên cứu ứng dụng do toán phân tích cơ học vật rắn hai chiều. Kết quả số đã 2
3. chứng minh rằng phương pháp PTHH trơn dựa trên 44zz33 u(,,) x y z u0 z 22xx cạnh ES-FEM [37] đã đạt được những thành tựu ấn 33hh tượng: 44zz33 v(,,) x y z v0 z 22yy ; (1) -FEM cho kết quả hội tụ nhanh và 33hh chính xác hơn phương pháp PTHH truyền thống sử w(,) x y w0 dụng phần tử tứ giác 4 nút với cùng số nút. Trường chuyển vị (1) chứa 7 ẩn số độc lập uvw000 xyxy  cần xác định. (spurious mode) và do đó phương pháp này cho lời T Trong đó: d000 uv  là các chuyển vị màng, w là giải ổn định nhất là với bài toán phân tích dao động tự 0 T do. độ võng,    xy là các góc xoay quanh trục y, và T trục x,  xy,   là các hàm độ cong. sử dụng thông số phạt, và hiệu quả tính toán tốt hơn 2.2. Trường biến dạng: phương pháp PTHH truyền thống với cùng số nút Từ trường chuyển vị (1) các biến dạng được xác định khảo sát. như sau: Phương pháp phần tử hữu hạn trơn gần đây đã được  Biến dạng trong mặt phẳng áp dụng nghiên cứu trên nhiều lĩnh vực như phân tích T 3 vết nứt, phân tích kết cấu dạng tấm vỏ trong môi ε pxxyyxy  εκκ02zz1 (2) trường đa vật lý, phân tích phi tuyến hình học kết cấu u  [38~42]. 0 x x x 2. Lý thuyết biến dạng cắt bậc cao:  v0 1 y ε0 (3); κ1  (4) y 2 y uv   00 x y yx yx  xx xx c yy κ 2  (5) 6 yy yy xx  yxyx  Biến dạng ngoài mặt phẳng w Hình 1: Hình học ban đầu và hình học biến dạng trên một x  xz x cạnh của tấm với các lý thuyết tấm cổ điển (CLPT), biến γ ε z2 κ (6) ; ε (7) ss s w dạng cắt bậc nhất (FSDT), và biến dạng cắt bậc 3 (TSDT).  yz  y y 2.1. Trường chuyển vị: Reddy [8] đã xây dựng trường chuyển vị của lý thuyết  4 tấm biến dạng cắt bậc cao dựa trên hàm xấp xỉ bậc 3, κ c xx (8) Với c (9) s  2 và từ điều kiện ứng suất cắt trên biên bằng 0 trường yy t chuyển vị được xác định như sau: 3
4.  N xx  xx 2.3. Trường ứng suất: h/2 N;  N yy  yy dz (15) h/2 Quan hệ ứnng suất biến dạng trong một lớp  Nxy  xy composite của vật liệu trực hướng  Mômen tổng ngoài mặt phẳng:  Ứng suất trong mặt phẳng M xx   xx   QQQ h/2 xx 111216 xx (16) MM yy   yy zdz;  QQQ (10) h/2  yy 212226 yy M   QQQ xy  xy  xy  616266 xy  Lực tổng ngoài mặt phẳng:  Ứng suất ngoài mặt phẳng  Qx h/2  xz  QQ  Q   dz; (17) yz 44 45 yz Q h/2    (11)  y yz xz  QQ54 55 xz  3. Phần tử tam giác 3 nút MITC3. Trong đó: y s E1 122E E2 d (3) (3) Q11 ; Q12 ; Q22 ; 1 1 1221 1 1221 1 1221 c QG ; QG ; QG . (12) 6 6 1 2 4 4 2 3 5 5 1 3 s (2) b r (2) 0 (1) a x 0 (1) r (a) (a) 1 Hình 3: Phần tử tam giác trong hệ tọa độ tự nhiên (a) và hệ tọa độ quy chiếu (b). Dạng hình học của phần tử có thể viết lại dưới dạng: 3 3 xNx  ii; yNy  ii. (18) 1 1 Hình 2.: Tấm composite gia cường sợi một phường với hệ Với: Nr1 ; Ns2 ; Nrs3 1 . (19) tọa độ tổng thể (x,y,z) và hệ tọa độ địa phương (x1,x2,x3). Vectơ chuyển vị nút phần tử: T k k k  deiiixiyixiyi u v w   i 1,3 (20) kkQQQ11 12 16 00 xx  xx  k k k QQQ21 22 26 00 Quan hệ giữa các biến dạng và chuyển vị của phần tử yy yy k k k (13) n n n xy  QQQ 00 xy  61 62 63 εBd m ; κ Bb1 d ; κBd b2  kk 0  Ii 1  Ii 2  Ii yz yz i 1 i 1 i 1 0 0 0 QQ44 45  xz  k xz  n n 0 0 0 QQ54 55 s0 s1 εs  BI d i ; κsI  B di (21) k k 4 k (k ) 2 2 k 4 i 1 i 1 QQ11 11 cos 2 QQ12 2 66 sin cos Q22 sin k k k k 2 2 k 4 4 Trong đó: QQQQ12 11 22 4 66 sin cos Q12 sin cos (14) k k 4 k k 2 2 k 4 QQ sin 2QQ 2 sin cos Q cos Ni 22 11 12 66 22 0 0 0 0 0 0 k k k k 33 k k k x QQQQ16 11 12 2 66 sin cos QQQ12 22 2 66 sin cos m Ni (22) k k k k 33 k k k BI 0 0 0 0 0 0 Q26 QQQ 11 12 2 66 sin cos QQQ12 22 2 66 sin cos y k k k k k 2 2 k 4 4 QQQQQ66 11 22 2 12 2 66 sin cos Q66 sin cos NN ii00000 yx 2.4. Nội lực trong tấm:  Lực tổng trong mặt phẳng: 4
5. N hiện tượng này cần xấp xỉ lại các biến dạng cắt bằng 000000 i x hàm mới thông qua những “điểm buộc” (tying point) b1 Ni (23) BI 000000 y s q s NNii 000 00 1 yx est eqt r 135° ert 0 Ni 1 r 00 Ni 000 x s0 (24) BI Ni Hình 4: Cách xác định biến dạng trượt ngang eqt 00000 Ni y Theo Lee và Bathe [4] 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ư 0 0 0NN 0 0 sau: s1 ii (25) BI 0 0 0 0NNii 0 1 e e e () (32) qt 2 strt NN ert a b r c s111 (33) 00000 ii xx est a b r c s222 (34) b2 c NNii (26) BI 00000 3 yy 1 eeeqt strt NNNN 2 000 iiii (35) yxyx 1 ab222111 rc sab rc s  2 Ma trận độ cứng phần tử s s TT KB  D BS  D S dd (27) 1 e ij isj 1 (3) (2) Trong đó: st eqt e mbb 12 ss01 BBBBiIII ; SBBiII (28) 0 r 0 1 (1) 1 r ert ABE ABSS Hình 5: Vị trí các điểm “tying point” cho phần tử tam * * giác 3 nút. Với D = B D F ; D=s (29) BDSS EFH Các điều kiện ràng buộc: h/2 (1) (1) 2346 (k ) ee(0,0) ; ee(1,0) (,,,,,)(1,A B D ,,,,),1,2,6 E F Hz z z z zQdzi j th rt rt rt rt ijijijijijij ij h/2 (2) (2) h/2 eest (0,0) st ; eest (0,1) st (36) (As , B s , Dz s )(1, zQ , dz ),4,524 i j (kth ) (30) ij ij ij ij h/2 (3) 1 (3) (3) eqt (1,0) eqt est e rt Vectơ tải phần tử: 2 1 FN pd (31) (3) (3) (3) e  eqt (0,1) eqt est e rt 2 Giải phương trình cân bằng: Kd = F để tìm các Từ sáu điều kiện trên ta có: chuyển vị. ae (1) ; b 0 ; c e(2) e (1) e (3) e (3) . Nếu chỉ dừng lại ở việc dùng phương trình (31) để 1 rt 1 1 st rt st rt (2) tìm các chuyển vị thì bài toán chỉ đúng với những tấm ae2 rt ; bc21 ; c2 0 . (37) có bề dày lớn hoặc vừa phải, còn đối với tấm mỏng sẽ Biến dạng trượt có thể được viết lại: dẫn đến hiện tượng “khóa cắt” (Shear locking) dẫn (1) (2) ert e rt cs ; est e st cr (38) đến kết quả không còn chính xác nữa. Để khắc phục 5
6. (2)(1)(3)(3) Trong đó: ceeee strtstrt (39) Trong đó (x) là một hàm có chức năng làm trơn và Quan hệ giữa các biến dạng cắt và các chuyển vị được thỏa mãn các điều kiện sau [37]: viết lại  0 và  1 (43) x k (x) nij  hrsks(,)| Bd=Bd s . (40) Để đơn giản hàm (x) được chọn là một hàm hằng ij  ijij kk 0iMITC (,s,rt ) k 1 ijij số theo từng mảng: 4. Phương pháp phần tử hữu hạn trơn ES-FEM 1 ( k) với phần tử MITC3 x   (x) A( k) (44) ( k) Trong phương pháp phần tử hữu hạn trơn, các biến 0 x  dạng được làm trơn trên các miền trơn địa phương, và k 1 Ne dĩ nhiên việc tính toán ma trận độ cứng không còn Với AdA( k)  là diện tích của miền  3  i phụ thuộc vào phần tử mà dựa trên các miền trơn này. i 1 k Những miền trơn này được xây dựng dựa trên cạnh trơn  của phần tử sao cho k Ne 1 đối với các cạnh biên. N ed k k ij N 2   và    với ij (41) e đối với các cạnh bên trong. k 1 Ai là diện tích của phần tử thứ i xung quanh cạnh k. Với Ned là tổng số cạnh của tất cả các phần tử trên Từ phương trình (51) và (53) biến dạng của phần tử miền khảo sát. trở thành: Đối với phần tử tam giác 3 nút miền trơn k liên 1 quan đến cạnh thứ k được tạo ra bằng cách kết nối hai  h d k ()k (45) A ()k nút cuối của cạnh với trọng tâm của phần tử liền kề.  Từ việc sử dụng các miền trơn dựa trên cạnh phần tử, Ma trận chuyển vị và biến dạng trơn được xác định bởi các biến dạng trơn k có thể xác định được từ sự 1 B()B()dxx  (46) tính toán các biến dạng tương ứng h kết hợp với Ik ()k I A ()k một thao tác làm trơn trên miền k liên quan đến  Đối với phẩn tử tam giác 3 nút MITC3, hàm dạng là cạnh thứ k: tuyến tính. Do đó ma trận biến dạng chuyển vị là hằng số trên phần tử. Kết hợp các phương trình từ   h (x)(x)d k (42) (21)~(26) và (46) ta có: ; ()k  k 1 Ne B()Bbb11x A (47) Ik ()k  ij A i 1 k 1 Ne Bbb22 ()Bx A Ik ()k  ij (48) A i 1 k 1 Ne Bss00 ()Bx A Ik ()k  ij MITC 3 (49) A i 1 kk NNee ss111 c * k B()xx B()dA B (50) I k ()kk A I () i Hình 6 Miền trơn  liên kết với cạnh trong ES FEM AAii 11i 3 6
7. m b1 b2 s0 Để khảo sát khả năng hội tụ của phần tử ESMITC3, B j ; B j ; B j ; B j M I T C 3 là các ma trận biến dạng xét trường hợp tấm có tỷ lệ t/L = 0,01 và số phần tử của phần tử thứ j xung quanh cạnh k. trên mỗi cạnh của tấm thay đổi NxN = 2x2, 4x4, 8x8, m b1 b2 B j ; B j ; B j được xác định từ (22) (23) và(26), 16x16, 32x32. Kết quả độ võng tương đối tại tâm tấm s0 w w/ ( qL4 /100 D ) hội tụ đến lời giải chính xác B j M I T C 3 được xác định từ (40) và c khi số phần tử tăng lên được thể hiện trong hình 8 đối * 0001010 B với tấm biên ngàm và hình 9 đối với biên tựa đơn. 0000101 Chú ý, wc là độ võng tại tâm tấm và D là độ cứng trụ Ma trận độ cứng phần tử ES-MITC3 y của tấm. Ned KK k e  (51) L/2 L/2 k 1 x x T T T KBABBBBBEBkmmmbmbk 12A IIIIII L/2 L/2 T T T BBBBDBBFBbmbbbbk1 1112 A IIIIII (NxN) (NxN) L L T T T bmbbbbk2 2122 BEBBFBBHB() A 52 q IIIII I q TT t t s0 ss 0 s0 ss 1 BAIMITC BBB3 B IMITCIMITC 3 3 I (a) (b) Ak TT B s1 Bss BBD0 B sss1 1 I IMITCI 3 I Hình 7 Tấm liên kết tựa đơn (a) và tấm liên kết ngàm (b) chịu tải trọng phân bố đều Với bài toán phân tích tĩnh giải phương trình: Từ khảo sát độ hội tụ của tấm, sử dụng lưới phần tử Kd = F (53) NxN = 32x32 để khảo sát khả năng khử hiện tượng Với bài toán phân tích động giải phương trình: khóa cắt của phần tử ES-MITC3 khi chiều dày tấm t 2 ()0KMd  (54) giảm dần (10-1 ~ 10-4) lần chiều dài tấm (tức là t/L = Trong đó M là ma trận khối lượng tổng thể được ghép 10-1 ~ 10-4). từ các ma trận khối lượng phần tử m. Kết quả độ võng tương đối tại tâm tấm w ứng với I1 0 0I2 0 c / 3 I4 0 chiều dày tấm giảm dần cho trường hợp tấm biên II0 0 0c / 3 I 1 2 4 ngàm và biên tựa đơn được lần lượt trình bày trong I1 0 0 0 0 (55) bảng 1 và bảng 2 m I3 0 c / 3 Is 0 ; I3 0c / 3 Is 2 cI/ 97 0 2 sym c/9 I7 Với h/2 (,IIIIII , , , , ) (1,z,z2 ,z, 3 z 4 ,z) 6 dz (56) 1 2 3 4 5 7 h/2 5. Các ví dụ số 5.1. Phân tích tấm đẳng hướng Cho tấm vuông đẳng hướng có kích thước cạnh a = b = L, chiều dày t, mô-đun đàn hồi E = 1092000 và hệ số Possion v = 0.3. Tấm liên kết tựa đơn (hình 7a) Hình 8 Độ võng tại tâm tấm (t/L=0.01) của các phần hoặc liên kết ngàm (hình 7b) trên bốn cạnh biên và tử theo số phần tử trên biên. ( Tấm chịu liên kết tựa chịu tải phân bố đều q = 1. đơn trên bốn cạnh biên) 7
8. Hình 10 Biểu đồ độ võng tại tâm tấm theo log(L/t) với Hình 9 Độ võng tại tâm tấm (t/L=0.01) của các phần tấm liên kết ngàm trên 4 cạnh biên. tử theo số phần tử trên biên. ( Tấm chịu liên kết ngàm trên bốn cạnh biên) Bảng 1: Độ võng tại tâm tấm với tấm liên kết ngàm 4 cạnh biên theo t/L Tỷ số t/L Phần tử 10-1 10-2 10-3 10-4 MITC3[4] 0.1504 0.1267 0.1263 0.1265 ESMITC3 0.1505 0.1268 0.1265 0.1265 Exact [35] 0.1499 0.1267 0.1265 0.1265 Bảng 2: Độ võng tại tâm tấm với tấm liên kết tựa đơn Hình 11 Biểu đồ độ võng tại tâm tấm theo log(L/t) với tấm liên kết ngàm trên 4 cạnh biên. cạnh biên theo t/L 5.2. Phân tích tấm composite 4 lớp (0/90/90/0) Tỷ số t/L Phân tích một tấm composite trực hướng vuông cạnh Phần tử a có chiều dày h, liên kết tựa đơn trên bốn cạnh biên 10-1 10-2 10-3 10-4 với 4 lớp vật liệu (0/90/90/0). Các mô đun đàn hồi của vật liệu như sau: E = 25E , G = G = 0.5E , G = MITC3[4] 0.4271 0.4062 0.4059 0.4059 1 2 12 13 2 23 0.2E2, ν12 = 0.25. Giả thiết các lớp có chiều dày bằng ESMITC3 0.4272 0.4064 0.4062 0.4062 nhau, tấm chịu tải trọng phân bố đều (hình 12) và tải Exact [35] 0.4273 0.4064 0.4062 0.4062 hình sin (hình 13). Kết quả phân tích ở bảng 3 và bảng 4, các biểu đồ ứng suất tiếp và ứng suất cắt xem Để so sánh khả năng khử “hiện tượng khóa cắt” của ở hình 14 đến hình 20 phần tử ESMITC3 so với những phần tử có khả năng khử “hiện tượng khóa cắt” khác như MIN3, DSG3, ESDSG3 có thể quan sát đồ thị ở hình 10 và hình 11 8
9. q 0 0 o 90 90 x b 0 h/2 h/2 y z a Hình 12 Tấm composite bốn lớp liên kết tựa đơn với tải Hình 15 Biểu đồ ứng suất tiếp  yy theo tọa độ z tại vị trọng phân bố đều. trí (a/2;b/2)(Trường hợp tải phân bố đều hình 12) 0 o 90 90 x b 0 q 0 h/2 h/2 y z a Hình 13 Tấm composite bốn lớp liên kết tựa đơn với tải xy trọng hình sin pqz 0 sinsin aa Các giá trị không thứ nguyên được xác định như sau: Hình 16 Biểu đồ ứng suất cắt  theo tọa độ z tại vị 100Eh3 aa haah2 xz 2 trí (0;b/2)(Trường hợp tải phân bố đều hình 12) ww 4 ;;0 xx 2 xx ;; qa0 22 qa0 222 h2 a a h hh2 ;yy 2 yy ;; xy 22 xy 0;0; ; qa0 2 2 4 qa0 2 h b h ha xz xz 0; ; ; yz yz ;0;0 qa0 22 qa0 2 Hình 17 Biểu đồ ứng suất cắt  yz theo tọa độ z tại vị trí (0;a/2) (Trường hợp tải phân bố đều hình 12) Hình 14 Biểu đồ ứng suất tiếp  xx theo tọa độ z tại vị trí (a/2;b/2)(Trường hợp tải phân bố đều hình 12) 9
10. E1 = 40E2, G12 = G13 = 0.6E2, G23 = 0.5E2, ν12 = 0.25. Tần số khảo sát được xác định bởi: 222 1/2  = bhD//0 với 3 DEh02 /12(1) 12 Kết quả thu được trong bảng 5 được so sánh với một số lời giải khác đã được công bố của những tác giả như sau: Liew (p-Ritz [43], 1996), Zhen và Wanji (Global-local higher-order theory [44], 2006), Ferreira, A. J. M., & Fasshauer [45] và Liew, K. M., Huang, Y. Q., & Reddy, J. N [46]. Bốn dạng dao động Hình 18: Biểu đồ ứng suất tiếp  theo tọa độ z tại vị xx (mode shapes) đầu tiên của kết cấu tấm vuông trí (a/2;0)(Trường hợp tải hình sin hình 13) composite 3 lớp [0/90/0] điều kiện biên ngàm 4 cạnh cũng được biểu diễn thông qua hình 21~25. Hình 21 Dạng dao động riêng của tấm ứng với mode 1 Hình 19: Biểu đồ ứng suất cắt  xz theo tọa độ z tại vị trí (0;b/2)(Trường hợp tải hinh sin hình 13) Hình 22 Dạng dao động riêng của tấm ứng với mode 2. Hình 20: Biểu đồ ứng suất cắt  yz theo tọa độ z tại vị trí (a/2;0)(Trường hợp tải hình sin hình 13) 5.3. Phân tích dao động tấm. Khảo sát tần số dao động riêng của tấm vuông composite ba lớp (0, 90, 0) kích thước a=b=L chiều dày t liên kết ngàm trên bốn cạnh biên với các đặc Hình 23 Dạng dao động riêng của tấm ứng với mode 3. trưng vật liệu như sau: 10
11. Hình 24 Dạng dao động riêng của tấm ứng với mode 4. Hình 25 Dạng dao động riêng của tấm ứng với mode 5. 6. Kết luận Việc xây dựng mô hình tính toán bằng phương pháp phần tử hữu hạn đã mang lại hiệu quả trong việc phân tích các bài toán ở các điều kiện biên, điều kiện tải trọng đa dạng và phức tạp. Khả năng khử “shear locking” của phần tử một lần nữa đã được kiểm tra và cho kết quả phân tích tốt. Phần tử tam giác 3 nút MITC3 được xây dựng lại kết hợp với lý thuyết biến dạng trơn đem lại một phần tử mới ES-MITC3 có tốc độ hội tụ nhanh và cho kết quả chính xác hơn một số phần tử tam giác 3 nút khác đã được công bố trên các tạp chí khoa học. Việc áp dụng phần tử ESMITC3 kết hợp với lý thuyết biến dạng cắt bậc cao mà cụ thể là hàm bậc 3 của Reddy cũng đã đem lại những kết quả như ứng suất cắt trên hai biên tấm bằng 0, và ứng suất phân bố theo chiều dày trong từng lớp của tấm là một đường cong. Điều này chứng minh lý thuyết biến dạng cắt bậc cao mô tả các ứng suất cắt thực tế hơn lý thuyết biến dạng cắt bậc nhất. Bảng 3:: Kết quả khảo sát độ võng không thứ nguyên và các ứng suất tấm vuông Composite chịu tải phân bố đều.(Trường hợp ứng với hình 12) Phương pháp w  xx  yy  xy  xz  yz [9] 3D-FEM 0.7794 0.8207 0.4870 0.0444 0.6040 0.4738 FSDT-MITC3 0.7689 0.8044 0.397 0.042 0.6710 HSDT-MITC3 0.7881 0.8136 0.4093 0.0435 0.5530 0.3486 HSDT-ESMITC3 0.7926 0.8179 0.4108 0.0433 0.5623 0.3612 [9] FSDT (NS-DSG3) 0.7745 0.80469 0.4010 0.0441 0.5797 [9] HSDT (NS-DSG3) 0.7966 0.8251 0.4144 0.0469 0.6940 0.3934 tL/0.05 [9] 3D-FEM 1.0576 0.8201 0.5605 0.0577 0.5781 0.498 FSDT-MITC3 1.0246 0.7577 0.5008 0.047 0.6425 HSDT-MITC3 1.1053 0.8206 0.5428 0.0544 0.5101 0.3875 HSDT-ESMITC3 1.1123 0.8252 0.5452 0.0538 0.5178 0.3991 [9] FSDT (NS-DSG3) 1.0304 0.7578 0.5040 0.0534 0.5312 [9] HSDT (NS-DSG3) 1.1184 0.8331 0.5495 0.0615 0.6534 0.4274 tL/0.1 [9] 3D-FEM 2.1044 0.8995 0.7386 0.0991 0.534 0.4487 FSDT-MITC3 1.9372 0.6504 0.7398 0.0572 0.4297 HSDT-MITC3 2.1645 0.9049 0.8189 0.0813 0.4173 0.4736 HSDT-ESMITC3 2.180 0.9104 0.8232 0.0803 0.4234 0.4843 [9] FSDT (NS-DSG3) 1.9454 0.6507 0.7419 0.0704 0.4316 [9] HSDT (NS-DSG3) 2.1936 0.9202 0.8287 0.0900 0.5846 0.5086 tL/ 0.2 11
12. Bảng 4: Kết quả khảo sát độ võng không thứ nguyên và các ứng suất tấm vuông Composite chịu tải hình sin.(Trường hợp ứng với hình 13) Phương pháp w  xx  yy  xy  xz  yz [9] RPF-PS-Layerwise 1.9023 0.7057 0.6309 0.0441 0.2138 [9] Elasticity 1.9540 0.7200 0.6660 0.0467 0.2700 [9] FSDT (NS-DSG3) 1.7316 0.4058 0.5765 0.0351 0.1378 [9] HSDT (NS-DSG3) 1.9266 0.7076 0.6303 0.0475 0.2084 0.2404 HSDT-MITC3 1.8863 0.6945 0.6211 0.0466 0.2101 0.2420 HSDT-ESMITC3 1.9203 0.7006 0.6251 0.0454 0.2122 0.2447 tL/ 0 . 2 5 [9] RPF-PS-Layerwise 0.7204 0.5609 0.3911 0.0273 0.2843 [9] Elasticity 0.7430 0.5590 0.4030 0.0276 0.3010 [9] FSDT (NS-DSG3) 0.6688 0.4983 0.3623 0.0257 0.1663 [9] HSDT (NS-DSG3) 0.7246 0.5609 0.3909 0.0288 0.2812 0.1566 HSDT-MITC3 0.7132 0.5522 0.3852 0.0275 0.2794 0.1562 HSDT-ESMITC3 0.7183 0.5558 0.3875 0.0270 0.2821 0.1558 t / 0L . 1 [9] RPF-PS-Layerwise 0.5078 0.5436 0.3052 0.0230 0.3066 [9] Elasticity 0.5170 0.5430 0.3090 0.0230 0.3280 [9] FSDT(NS-DSG3) 0.2968 0.0225 0.1784 [9]HSDT(NS-DSG3) 0.5089 0.5433 0.3050 0.0234 0.3051 0.1266 HSDT-MITC3 0.5024 0.5352 0.3005 0.0229 0.3013 0.1249 HSDT-ESMITC3 0.5055 0.5385 0.3023 0.0227 0.3043 0.1275 tL/0.05 Bảng 3.5: Kết quả phân tích tần số dao động riêng của tấm. Phương pháp Modes 1 2 3 4 5 [43]Liew 4.447 6.642 7.700 9.185 9.738 [44]Zhen and Wanzi 4.45 6.524 8.178 9.473 9.492 [45]RBF-PS 4.5141 6.508 8.0361 9.3468 9.3929 [9]HSDT(NS-DSG3) 4.4802 6.4046 7.9291 9.1682 9.1763 HSDT(ES-MITC3) 4.5092 6.4854 8.0354 9.3238 9.3731 Lt/5 [43]Liew 7.411 10.393 13.913 15.429 15.806 [44]Zhen and Wanzi 7.484 10.207 14.340 14.863 16.070 [45]RBF-PS 7.4727 10.2544 14.244 14.9363 15.9807 [46] MLSDQ 7.4320 10.3990 13.9580 15.4670 15.8380 [9]HSDT(NS-DSG3) 7.4224 10.1266 14.0657 14.6205 14.6205 HSDT(ES-MITC3) 7.4723 10.2540 14.2562 14.9547 15.9948 Lt/ 10 [43]Liew 10.953 14.028 20.388 23.196 24.978 [44]Zhen and Wanzi 11.003 14.064 20.321 23.498 25.350 [45]RBF-PS 10.968 13.9636 20.0983 23.3572 25.0859 [9]HSDT(NS-DSG3) 10.9042 13.8634 19.8429 23.0955 24.8382 HSDT(ES-MITC3) 10.9821 14. 0114 20.2343 23.4191 25.2004 Lt/ 20 12
13. 7. Tài liệu tham khảo higher-order plate theory, Compos. Struct. 48 (2000) [1] A.A. Khdeir, L. Librescu, Analysis of symmetric 231–244. cross-ply laminated elastic plates using a higher-order [12] T. Kant, K. Swaminathan, Analytical solutions theory, Compos. Struct. 9 (1988) 189–213. for free vibration of laminated composite and [2] N.J. Pagano, Exact solutions for rectangular sandwich plates based on a higher order refined bidirectional composites and sandwich plates, J. theory, Compos. Struct. 53 (2001)73–85. Compos. Mater. 4 (1970) 20–34. [13] L. Liu, L.P. Chua, D.N. Ghista, Mesh-free radial [3] Hyeong-Min Jeon, Phill-Seung Lee, Klaus-Jürgen basis function method for static, free vibration and Bathe. The MITC3 shell finite element enriched by buckling analysis of shear deformable interpolation covers.Computers and Structures 134 compositelaminates, Compos. Struct. 78 (2007) 58– (2014) 128–142 69. [4] Phill-Seung Lee, Klaus-Jurgen Bathe. [14] Ahmad, S., Irons, B.M. and Zienkiewicz, O.C. Development of MITC isotropic triangular shell finite (1970), “Analysis of thick and thin shell structures by elements. Computers and Structures 82 (2004) 945– curved finite elements”, Int. J. Num. Meth. Engrg, 2, 962. pp. 419–451. [5] A.K. Noor, Free vibration of multilayered [15] Zienkiewicz, O.C. and Taylor, R.L. (1989), “The composite plates, AIAA J.11 (1973)1038–1039. Finite Element Method”, (Fourth Edition), McGraw [6] A.K. Noor, Stability of multilayered composite Hill. plates, Fibre Sci. Technol. 8 (1975) 81–89. [16] Reissner, E. (1945), “The effect of transverse [7] Neeraj Grover, D.K. Maiti, B.N. Singh. A new shear deformation on the bending of elastic plates”, J. inverse hyperbolic shear deformation theory for static Appl. Mech., 12, pp. 69–76. and buckling analysis of laminated composite and [17] Mindlin, R.D. (1951), “Influence of rotatory sandwich plates. Composite Structures 95 (2013) inertia and shear in flexural motions of isotropic 667–675 elastic plates”, J. Appl. Mech., 18, pp. 31–38. [8] J.N. Reddy, Mechanics of Laminated Composite [18] Eduardo N. Dvorkin, “Nonlinear Analysis of Plates – Theory and Analysis, CRC Press, New York, Shells Using the MITC Formulation”. Archives of 1997. Computational Methods in Engineering, Vol. 2, 2, 1– [9] Chien H. Thai, Loc V. Tran, Dung T. Tran, T. 50 (1995). Nguyen-Thoi, H. Nguyen-Xuan. Analysis of [19] Dvorkin EN, Bathe KJ. A continuum mechanics laminated composite plates using higher-order shear based four-node shell element for general nonlinear deformation plate theory and node-based smoothed analysis. Eng Comput 1984;1: 77–88. discrete shear gap method. [20] Bathe KJ, Dvorkin EN. A formulation of general [10] An Edge-based smoothed discrete shear gap shell elements – the use of mixed interpolation of method (ES-DSG) using the C0-type Higher-order tensorial components. Int J Numer Meth Eng 1986; shear deformation Theory for Analysis of Laminated 22: 697–722. Composite Plates. Loc V. Tran , T. Nguyen-Thoi, [21] Bucalem ML, Bathe KJ. Higher-order MITC Chien H. Thai, H. Nguyen-Xuan. Mechanics of general shell elements. Int J Numer Meth Eng Advanced Materials and Structures. 1993;36:3729–54. [11] H. Matsunaga, Vibration and stability of cross- [22] Kim J, Bathe KJ. The finite element method ply laminated composite plates according to a global enriched by interpolation covers. Comput Struct 2013;116:35–49. 13
14. [23] Y.X. Zhang, K.S. Kim. A simple displacement- [29] A cell-based smoothed three-node Mindlin plate based 3-node triangular element for linear and element (CS-FEM-MIN3) based on the C0-type geometrically nonlinear analysis of laminated higher-order shear deformationfor geometrically composite plates. Comput. Methods Appl. Mech. nonlinear analysis of laminated composite plates. Engrg. 194 (2005) 4607–4632. [30] Computers and Structures 110–111 (2012) 93– [24] R.K. Kapania, P. Mohan, Static, free vibration 106.Youngyu Lee, Kyungho Yoon, Phill-Seung and thermal analysis of composite plates and shells Lee.Improving the MITC3 shell finite element by using a flat triangular shell element, Comput. Mech. using the Hellinger–Reissner principle. 17 (1996) 343–357. [31] Lee PS, Bathe KJ. Insight into finite element [25] E. Madenci, A. Barut, A free-formulation-based shell discretizations by use of the ‘‘basic shell flat shell element for nonlinear analysis of thin mathematical model’’. Comput Struct 2005;83:69–90. composite structures, Int. J.Numer. Methods Engrg. [32] Dvorkin EN, Bathe KJ. A continuum mechanics 37 (1994) 3825–3842. based four-node shell element for general nonlinear [26] O. Polit, M. Touratier, New laminated triangular analysis. Eng Comput 1984;1:77–88. finite element assuring interface continuity for [33] Bathe KJ, Dvorkin EN. A formulation of general displacements and stresses, Compos. Struct. 38 (1997) shell elements – the use of mixed interpolation of 37–44. tensorial components. Int J Numer Meth Eng [27] O. Polit, M. Touratier, High-order triangular 1986;22:697–722. sandwich plate finite element for linear and non-linear [34] Bucalem ML, Bathe KJ. Higher-order MITC analyses, Comput. Methods Appl. Mech. Engrg. 185 general shell elements. Int J Numer Meth Eng (2000) 305–324. 1993;36:3729–54. [28] Trung-Kien Nguyen, Thuc P. Vo, Huu-Tai Thai, [35] H. Nguyen-Xuan, T. Rabczuk, St´ephane Bordas Vibration and buckling analysis of functionally J.F.Debongnie. A smoothed finite element method for graded sandwich plates with improved transverse plate analysis. shear stiffness based on the first-order shear [36] Liu, G.R., Dai, K.Y., Nguyen-Thoi T., A deformation theory. Institution of Mechanical smoothed finite element method for mechanics Engineers, Part C: Journal of Mechanical Engineering problems, Computational Mechanics 39(6) (2007) Science, 228 (12), 2110-2131, 2014 (SCI). 859–877. [37] Liu, G.R., Nguyen-Thoi T., K.Y., Lam. An edge- [39] Nguyen-Van Hieu, Vo Anh Vu, Nguyen Hoai based smoothed finite element method (ES-FEM) for Nam, Chau Dinh Thanh, Nguyen Ngoc Duong static and dynamic problems of solid mechanics, J. (2012). “Analysis of shell structures via a smoothed Sound.Vib. 320(2008) 1100–1130. four-node flat element”. Proceeding of the [38] Nguyen Hoai Nam, Nguyen Van Hieu, Luong International Conference On Advances In Van Hai, Chau Dinh Thanh, Nguyen Ngoc Duong Computational Mechanics (ACOME), August 14-16, (2013). "A smoothed strain based element for 2012, Ho Chi Minh City, Viet Nam, pp.219-233. geometrically nonlinear analysis of plate/shell [40] Nguyen-Xuan H, Liu GR, Nguyen-Thoi T, and structures". The 4th International Conference of Euro Nguyen-Tran C. 2009. An edge-based smoothed finite Asia Civil Engineering Forum 2013 (EACEF): element method (ES-FEM) for analysis of two Innovations in Civil Engineering for Society and the dimensional piezoelectric structures. Smart Materials Environment, National University of Singapore, and Structures. SINGAPORE, June 26-28, 2013, pp. S81-S86. [41] Liu GR, Nguyen-Thoi T, Nguyen-Xuan H, and Lam KY. 2009. A node-based smoothed finite element 14
15. method (NS-FEM) for upper bound solution to solid [45] Ferreira, A. J. M., & Fasshauer, G. E. (2007). mechanics problems. Computers and Structures; 87: Analysis of natural frequencies of composite plates by 14 –26. an RBF-pseudospectral method. Composite [42] Cheng L and Liu GR. 2008. Adaptive analysis Structures,79(2),202-210.doi:DOI:10.1016/j. using the edge-based smoothed finite element method; compstruct. 2005.12.004 (submitted). [46] Liew, K. M., Huang, Y. Q., & Reddy, J. N. [43] Liew, K. M. (1996). Solving the vibration of (2003). Vibration analysis of symmetrically laminated thick symmetric laminates by Reissner/Mindlin plate plates based on FSDT using the moving least squares theory and the p-Ritz method. Journal of Sound and differential quadrature method. Computer Methods in Vibration, 198(3), 343-360. doi:DOI: Applied Mechanics and Engineering, 192(19), 2203- 10.1006/jsvi.1996.0574. 2222. doi:DOI: 10.1016/S0045-7825(03)00238 [44] Zhen, W., & Wanji, C. (2006). Free vibration of laminated composite and sandwich plates using global–local higher-order theory. Journal of Sound 15
16. 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 2016-2017 của Thư viện Trường Đại học Sư phạm Kỹ thuật Tp. Hồ Chí Minh.