Phân tích ổn định tĩnh vỏ mỏng bằng phương pháp phần tử hữu hạn - Vũ Tùng Lâm

C¬ kü thuËt & Kü thuËt c¬ khÝ ®éng lùc V.T. Lâm, N.V.Chúc, T.N. Thanh, “Phân tích ổn định tĩnh phần tử hữu hạn.” 158 ph©n tÝch æn ®Þnh tÜnh vá máng b»ng ph­¬ng ph¸p phÇn tö h÷u h¹n VŨ TÙNG LÂM, NGUYỄN VĂN CHÚC, TRẦN NGỌC THANH Tóm tắt: Đánh giá sự ổn định của vỏ có vai trò quan trọng trong việc thiết kế các kết cấu vỏ mỏng. Mục đích chính của bài toán nhằm kiểm tra sự ổn định của vỏ trong các điều kiện chịu lực khác nhau. Trong báo đề cập tới vấn đề phân tích ổn định của vỏ mỏng bằng phương pháp phân tử hữu hạn. Phần tử được sử dụng để mô hình hóa vỏ là phần tử vỏ cong hai chiều đẳng tham số. Mô hình toán học, giải thuật và chương trình tính toán đã được xây dựng và áp dụng tính cho vỏ trụ mỏng đơn giản. Các kết quả được so sánh với kết quả tính trên phần mềm ANSYS. Từ khóa: Ổn định tĩnh, Phần tử hữu hạn, Kết cấu vỏ mỏng. 1. ĐẶT VẤN ĐỀ Trong thiết kế các kết cấu thành mỏng nếu chỉ tính đến độ bền và độ cứng thì chưa đủ đánh giá khả năng làm việc của kết cấu. Trong nhiều trường hợp đặc biệt là đối với các kết cấu thành mỏng chịu tải trọng nén hoặc nén và uốn tuy tải trọng chưa đạt tới giá trị giới hạn theo độ bền hoặc độ cứng nhưng kết cấu đã mất khả năng làm việc. Do đó, bắt buộc phải đánh giá được khả năng ổn định của kết cấu. Bài toán ổn định tuyến tính của vỏ mỏng đã được nghiên cứu rộng rãi từ rất lâu. Trong tài liệu [2], phương pháp tính toán đánh giá khả năng ổn định của vỏ mỏng thường sử dụng các công thức kinh nghiệm. Khi áp dụng các công thức này có ưu điểm là kết quả tính toán rất nhanh tuy nhiên độ chính xác không cao. Phương pháp phần tử hữu hạn (PTHH) cho phép đánh giá chính xác khả năng ổn định của kết cấu và có thể áp dụng để giải quyết nhiều dạng kết cấu khác nhau. Khảo sát bài toán ổn định của vỏ bằng phương pháp PTHH trên cơ sở phương giải bài toán trị riêng [4, 6, 7]:      ,i i G iK K   (1) trong đó, [K] là ma trận độ cứng tuyến tính của kết cấu;[KG] – ma trận độ cứng hình học (độ cứng ứng suất) của kết cấu;{ψi} – véc tơ riêng; λi – trị riêng, trong bài toán ổn định tuyến tính được gọi là hệ số tải trọng. Giải phương trình trị riêng trên ta xác định được các trị riêng 1 2 ... n     , để đánh giá độ ổn định của hệ thì nhân tử tải nhỏ nhất λ1 là có ý nghĩa. Dựa vào giá trị riêng tìm được ta có thể đánh giá sự ổn định của kết cấu như sau nếu 1>1 thì hệ ổn định; nếu 1=1 thì hệ ở trạng thái tới hạn còn 1<1 hệ không ổn định. Như vậy để thực hiện một phân tích ổn định bằng PTHH trước tiên ta phải rời rạc hóa kết cấu thành những phần tử sau đó xác định các ma trận độ cứng tuyến tính, giải bài toán tĩnh để xác định trạng thái ứng suất, tính ma trận độ cứng hình học và cuối cùng là giải bài toán trị riêng. Trong phân tích vỏ, cách tiếp cận đầu tiên là kết hợp phần tử màng với phần tử tấm uốn trên cơ sở các giả thiết của Kirchhoff. Lớp phần tử này có nhược điểm là mô phỏng hình học vỏ chỉ xấp xỉ với bề mặt thật, tính chất của vỏ được mô tả bởi hệ phương trình vi phân không hoàn toàn đồng nhất với tấm, sự không liên tục của các đạo hàm qua biên phần tử gây ra ứng suất uốn mà trên thực tế là không có. Để tăng độ chính xác người ta sử dụng phần tử khối 3 chiều tổng quát để phân tích vỏ tuy nhiên chúng làm tăng khối lượng tính toán lên nhiều lần. Để giảm khối lượng tính toán ta có thể sử dụng phần tử giảm bậc tổng Nghiên cứu khoa học công nghệ Tạp chí Nghiên cứu KH&CN quân sự Số 32, 08 - 2014 159 quát của phần tử bằng cách áp dụng giả thiết của Mindlin, tuy nhiên chúng dễ gặp phải hiện tượng “nghẽn cắt” khi tính toán vỏ mỏng[1]. Để khắc phục những nhược điểm người ta đã phát triển phần tử vỏ cong để phân tích vỏ. Bài báo đề cập đến vấn đề sử dụng phần tử vỏ cong kép đẳng tham số được phát triển bởi Saetta và Vitaliani [5] để mô hình hóa vỏ mỏng. Trên cơ sở đó tiến hành khảo sát sự ổn định của vỏ mỏng khi chịu tác dụng của các ngoại tải. Các kết quả tính toán được so sánh, đánh giá với các kết quả mô phỏng trên phần mềm ANSYS. 2. MÔ HÌNH PHẦN TỬ HỮU HẠN CỦA VỎ MỎNG 2.1. Rời rạc hóa kết cấu vỏ Phần tử vỏ cong kép là là phát triển của phần tử tấm Mindlin-Reissner. Dạng hình học của phần tử được mô tả trên hình 1. Trong hệ tọa độ cục bộ phần tử có dạng tứ giác cong, trong hệ tọa độ tự nhiên () phần tử có dạng hình vuông trong mặt phẳng () với các đỉnh là các điểm có tọa độ (,) là đơn vị.   Hình 1. Dạng hình học, các hệ tọa độ và mô hình phần tử vỏ. Hàm dạng của phần tử là hàm dạng hai chiều thuộc lớp C0 [3]. Biểu thức tiệm cận cho hình học của vỏ có dạng:     8 1 . T T i i i i i x y z N x y z   (2) Đạo hàm của các tọa độ tổng thể theo các tọa độ tự nhiên được tính như sau: 8 8 8 1 1 1 8 8 8 1 1 1 ; ; ; ; ; . i i i i i i i i i i i i i i i i i i N N Nx y z x y z N N Nx y z x y z                                                     (3) Ta nhận được ma trận chuyển tọa độ Jacobian:   . x y z J x y z                          (4) Định thức của ma trận chuyển hệ tọa độ không tính được theo cách thông thường mà được tính như sau: 2 2 2 1 2 3 1 2 3; ; ; . y z x z x y J J J J J J J y z x z x y                                           (5) C¬ kü thuËt & Kü thuËt c¬ khÝ ®éng lùc V.T. Lâm, N.V.Chúc, T.N. Thanh, “Phân tích ổn định tĩnh phần tử hữu hạn.” 160 Trên mặt trung hòa của vỏ ta định nghĩa hệ tọa độ cong cục bộ (sr) với:     2 2 2 2 2 2 0 0 ; . x y z x y z s d r d                                                                  (6) Độ cong của hệ tọa độ cong theo Lame: 3/2 3/22 2 2 2 2 2 1 2 1 2 1 1 ; , x y z x y z C C R R                                                                      (7) trong đó: 1/22 2 22 2 2 2 2 2 1 2 2 2 2 2 2 1/2 2 22 2 2 2 2 2 2 2 2 2 2 2 2 ; y z z y z x x z x y y x C y z z y z x x z x y y x C                                                                                                                           2 . (8) Đạo hàm của hàm dạng theo s và r: 1/22 2 2 1/22 2 2 ; . i i i i N N x y z s N N x y z r                                                                         (9) 2.2. Nội suy các biến trường Chuyển vị tại mỗi nút của phần tử bao gồm 5 thành phần: chuyển vị theo các trục s, r và theo phương pháp tuyến với mặt rs cùng với hai góc xoay quanh các trục này    T i i i i si rid u v w   . Chuyển vị tại điểm bất kỳ trong phần tử cũng bao gồm 5 thành phần và được tính theo quan hệ:     ,s r eu v w N d   (10) với:    1 1 1 1 1 8 8... T e s r s rd u v w     - véc tơ chuyển vị nút phần tử; [N] – ma trận hàm dạng:            1 2 8 5... ; ,i iN N N N N N I    (11) [I5] là ma trận đơn vị cấp 5. Quan hệ ứng suất, biến dạng vỏ được xây dựng theo lý thuyết biến dạng trượt của vỏ cong tính đến bán kính cong R1, R2và giả thuyết của Mindlin có thể được biểu diễn như sau [4], [5]: Nghiên cứu khoa học công nghệ Tạp chí Nghiên cứu KH&CN quân sự Số 32, 08 - 2014 161 1 2 2 2 0 2 1 1 2 1 2 ( ) ( ) s r rsr s sr sr s s r r r s u w s R v w w r R su v wr s rz k ws k szk r v u z C s r s r w v s R w u r R                                                                                             , 0 0 0 0 0 T L w r                                (12) trong đó, {εT}, {εL} là véc tơ biến dạng đặc trưng cho trạng thái biến dạng tuyến tính và phi tuyến tương ứng;  1 20,5 1/ 1/ .oC R R  (13) Đối với thành phần biến dạng tuyến tính sau khi tính toán các thành phần biến dạng xếp chúng vào dạng ma trận ta nhận được quan hệ giữa biến dạng và chuyển vị nút sau:     ,T eB d  (14) trong đó ma trận [B] là ma trận chuyển vị nút biến dạng có kích thước 40x40. Với thành phần biến dạng phi tuyến ta có thể biểu diễn như sau:      1 . 2 L A  (15) Ma trận [A] và véc tơ {θ} có thể được biểu diễn dưới dạng:         0 0 0 0 0 0 ; ,0 0 0 0 0 0T s r s er s rA G d             (16) trong đó:          1 2 8 0 0 0 0 ... ; . 0 0 0 0 i i i N G G G G G N         (17) Trường ứng suất trong phần tử được biểu diễn như sau:       , T s r sr s r sr s r TN N N M M M Q Q C   (18) trong đó [C] là ma trận vật liệu của phần tử vỏ phụ thuộc vào mô đun đàn hồi E,hệ số poát xông µ,và độ dày vỏ h. Ma trận vật liệu là ma trận vuông có kích thước 8x8 đã được trình bày cụ thể trong[4]. 2.3. Tính toán các ma trận độ cứng Ma trận độ cứng tuyến tính của phần tử được cho theo công thức tổng quát sau [1, 3]:        . T tev V k B C B dV  (19) Chuyển về hệ tọa độ tự nhiên ta có: C¬ kü thuËt & Kü thuËt c¬ khÝ ®éng lùc V.T. Lâm, N.V.Chúc, T.N. Thanh, “Phân tích ổn định tĩnh phần tử hữu hạn.” 162        1 1 1 1 det . T ek B C B Jd d       (20) Ma trận độ cứng tuyến tính tổng thể được xây dựng trên cơ sở ghép nối các ma trận độ cứng phần tử vào ma trận độ cứng tổng thể theo đúng thứ tự các bậc tự do [1, 3]. Để xác định được ma trận độ cứng hình học trước tiên ta phải xác định được trạng thái ứng suất trong kết cấu. Để làm được điều này ta thực hiện một phân tích trạng thái ứng suất biến dạng tĩnh. Phương trình cân bằng có dạng sau [1, 3]:     ,K d Q (21) trong đó {Q} là véc tơ tải tổng thể có được bằng cách ghép nối các véc tơ tải phần tử {Qe} theo đúng thứ tự. Véc tơ tải phần tử được cho theo công thức tổng quát [1, 3]:           , T T e S V S Q N X dV N p dS   (22) trong đó: {X} – véc tơ lực khối;{p} – véc tơ lực diện;NS – hàm dạng trên biên S;[NS] – ma trận hàm dạng trên biên S được xác định tương tự như ma trận hàm dạng [N]. Chuyển về hệ tọa độ tự nhiên ta có:           1 1 1 1 1 1 1 1 det det . T T e SQ h N X Jd d N p Jd d             (23) Sau khi xác định được trạng thái ứng suất ma trận độ cứng hình học của phần tử được tính như sau [3, 4, 6]:           1 1 1 1 det , T T ge V k G S G dV G S G Jd d            (24) với [S] – là ma trận ứng suất:   .s sr sr r N N S N N      (25) Ma trận độ cứng hình học tổng thể được xây dựng tương tự ma trận độ cứng tuyến tính. 2.4. Giải bài toán trị riêng Thông thường người ta tiến hành giải bài toán giá trị riêng ở dạng bài toán giá trị riêng chuẩn. Bằng cách nhân trái phương trình (1) với [K]-1 chúng ta nhận được:       1 0,i i I D          (26) trong đó [I] là ma trận đơn vị, còn ma trận [D] xác định bởi:       1 GD K K   . Bài toán trị riêng của phương trình (26) được gọi là bài toán giá trị riêng chuẩn. Để nghiệm {ψ} là không tầm thường, định thức đặc trưng cần phải bằng không, nghĩa là:         1 1 det 0. i i I D I D            (27) Sau khi khai triển, phương trình (27) cho một phương trình đa thức bậc n. Viêc khai triển định thức đặc trưng, sau đó giải phương trình đa thức bậc n thu được để tìm ra nhân tử tải (hệ số tải trọng) có thể sẽ rất mất thời gian một khi n có giá trị lớn. Người ta đã phát Nghiên cứu khoa học công nghệ Tạp chí Nghiên cứu KH&CN quân sự Số 32, 08 - 2014 163 triển một vài phương pháp số để tính toán nhân tử tải của hệ nhiều bậc tự do, thường được sử dụng nhất là phương pháp lặp ma trận. Phương pháp lặp ma trận thừa nhận rằng các nhân tử tải của hệ là khác nhau và phân bố tách biệt theo trật tự 1 2 ... n     . Phép lặp được bắt đầu với việc chọn một vectơ làm giá trị thử cho vectơ {ψ1}, rồi đem nhân trái {ψ1} với ma trận [D]. Ma trận cột thu được tiếp đó được chuẩn hóa. Đem vectơ cột đã được chuẩn hóa nhân trái với [D] để nhận được vectơ cột thứ ba, vectơ cột này lại được chuẩn hóa theo cách tương tự và đi đến một vectơ cột thử khác. Lặp lại quá trình trên cho đến khi các vectơ cột chuẩn hóa kế tiếp nhau hội tụ đến một vectơ cột chung, đó là vectơ riêng cơ bản. Hệ số chuẩn hóa cho giá trị lớn nhất của 1/λ, hay cho giá trị nhỏ nhất của λ. 3. CHƯƠNG TRÌNH VÀ VÍ DỤ TÍNH TOÁN SỐ 3.1. Thuật toán và chương trình tính Từ mô hình phân tích phần tử hữu hạn trên các tác giả đã xây dựng chương trình phân tích ổn định kết cấu vỏ mỏng. Các tích phân tính ma trận độ cứng và véc tơ tải nút được tính theo phương pháp tích phân số cầu phương Gauss theo sơ đồ 2x2 [1, 3]. Sơ đồ giải thuật cho bài toán phân tích ổn định tĩnh vỏ mỏng như hình 2: Hình 2. Sơ đồ giải thuật bài toán bài toán phân tích ổn định tĩnh vỏ mỏng. Để giải bài toán trị riêng ta sử dụng thuật toán lặp không gian con, thuật toán này đã được tích hợp sẵn trong một số chương trình lập trình như Matlab.Kết quả của bài toán là giá trị nhân tử tải nhỏ nhất. Các chương trình được xây dựng bằng Matlab là ngôn ngữ phổ biến và có nhiều tiện ích thích hợp cho việc phân tích PTHH. 3.2. Ví dụ tính toán số Một kết cấu vỏ trụ mỏng chịu tải nén dọc trục được phân tích bằng chương trình đã xây dựng. Kết cấu cũng được khảo sát bằng ANSYS, các kết quả của ANSYS được dùng như là một sự đối chứng để đánh giá tính chính xác và ổn định của phương pháp cũng như chương trình đã thiết lập. Bắt đầu Nhập các thông số đầu vào Xây dựng mô hình hình học Tính ma trận độ cứng tuyến tính Tính véc tơ tải trọng nút Giải phương trình tổng quát xác định ứng suất Tính ma trận độ cứng hình học Giải bài toán trị riêng xác định nhân tử tải Kết thúc C¬ kü thuËt & Kü thuËt c¬ khÝ ®éng lùc V.T. Lâm, N.V.Chúc, T.N. Thanh, “Phân tích ổn định tĩnh phần tử hữu hạn.” 164 Mô hình hình học và các tham số kích thước được mô tả trên hình 3. Kết cấu là một vỏ hình trụ trơn có đường kính ngoài chiều dài và độ dày thành cho trước, chịu tải trọng nén dọc trục. Một đầu được ngàm cứng.Vật liệu sử dụng cho kết cấu được xác định là hợp kim nhôm. Hình 3. Mô hình hình học, tải trọng và điều kiện biên của bài toán. Giá trị của các tham số kích thước, tải trọng, vật liệu được cho trong bảng 1. Bảng 1. Giá trị tham số đầu vào. TT Tham số Đơn vị Giá trị 1 Đường kính ngoài D mm 420 2 Độ dàyδ mm 0,5;1;1,5;2 3 Chiều dài L mm 720 4 Lựcphânbốtheochu vi p N/mm 10 5 Mô đun Young E MPa 7.1x104 6 Hệ số poát xông µ 0,33 Kết quả tính toán bằng phần mềm ANSYS và chương trình tự xây dựng được cho trong bảng 2. Bảng 2. Các kết quả tính toán. Độ dày vỏ  tính theo ANSYS  tính theo chương trình Sai lệch (%) 0,5 0,5957 0,6125 2,8 1 2,6966 2,5438 3,8 1,5 6,8920 6,6079 4,1 2 12,9696 12,6768 2,3 Dựa trên các kết quả tính toán ta xây dựng được đồ thị biểu diễn mối quan hệ giữa độ dày vỏ và hệ số tải trọng (hình 4). Ta nhận thấy các kết quả tính bằng phần mềm tự xây dựng là ổn định và phù hợp tốt với các kết quả tính toán nhận được bằng phần mềm ANSYS. Sai số lớn nhất giữa 2 chương trình tính không vượt quá 5%. Mô hình PTHH và chương trình có thể sử dụng trong tính toán thiết kế thân vỏ khoang KCB với thời gian tính toán thấp hơn sử dụng ANSYS. Hình 3. Quan hệ giữa hệ số tải trọng và độ dày vỏ. Nghiên cứu khoa học công nghệ Tạp chí Nghiên cứu KH&CN quân sự Số 32, 08 - 2014 165 4. KẾT LUẬN Trong bài báo đã xây dựng mô hình PTHH phân tích ổn định tuyến tính cho vỏ mỏng trên cơ sở phần tử vỏ cong kép đẳng tham số 8 nút. Mô hình này có thể sử dụng để phân tích PTHH cho vỏ có hình dạng bất kỳ. Thuật toán và chương trình tính toán đánh giá khả năng ổn định của vỏ mỏng cũng đã được xây dựng. Một dạng kết cấu vỏ trụ với độ dày vỏ khác nhau đã được khảo sát. Kết quả nhận được đã được so sánh với các kết quả tính toán mô phỏng bởi chương trình ANSYS cho thấy mô hình toán và chương trình đã xây dựng là tin cậy và phù hợp tốt.Các nghiên cứu trên là cơ sở ban đầu giải quyết bài toán ổn định trong các trường hợp vỏ tổng quát. Dựa trên mô hình cơ sở này hoàn toàn có thể phát triển, mở rộng để khảo sát đối với các vỏ có gân gia cường và vỏ compsite... TÀI LIỆU THAM KHẢO [1]. Nguyễn Quốc Bảo, Trần Nhất Dũng (2003). Phương pháp phần tử hữu hạn lý thuyết và lập trình. NXB Khoa học và Kỹ thuật, Hà nội. [2]. Nguyễn Hoa Thịnh, Hoàng Xuân Lượng, Nguyễn Đức Cương, Trịnh Hồng Anh, Nguyễn Minh Tuấn (2005). Kết cấu và tính toán độ bền khí cụ bay. NXBKH&KT. [3]. Zienkiewicz, O. C. and R. L. Taylor (2000). The Finite Element Method. Butterworth Heinemann, Oxford, fifth edition. [4]. Jordanka Ivanova, Franco Pastrone(2003). Geometric method for stability of non- linear elastic thin shells. KluwerAcademicpublishers. [5]. Saetta, A.V., Vitaliani, R.V., (1990). A finite element formulation for shells of arbitrary geometry. Computers and Structures 37, pp. 781–793 [6]. А. И. Голованов, О. Н. Тюленева, А. Ф. Шигабудиновю (2006). Метод конечных элементов в статике и динамике тонкостеных конструкций. Физматлит, Москва. [7]. Образцов И.Ф., Савельев Л.М., Хазанов Х.С (1985). Метод конечных элементов в задачах строительной механики летательных аппаратов. Высшая школа, Москва. ABSTRACT LINEAR BUCKLING ANALYSIS OF THIN SHELL BY FINITE ELEMENT METHOD Linear buckling analysis of thin shell plays an important role in the design of thin shell structures. The main purpose of the problem is checking the stability of the shell under various loading. This report refers to the linear buckling analysis of thin shellby finite element method. Doubly curved isoparametric shell element has been used to model thin shell. Mathematical model, algorithm flowcharts and some Matlab code were built and applied to solve simple cylinder shell structure. Results are discussed and compared with those obtained by ANSYS software. Keywords: Linear buckling, Finite element method, Thin shell. Nhận bài ngày 14 tháng 04 năm 2014 Hoàn thiện ngày 18 tháng 5 năm 2014 Chấp nhận đăng ngày 12 tháng 7 năm 2014 Địa chỉ: Viện Tên lửa – Viện Khoa học và Công nghệ Quân sự