Xây dựng chương trình lập ma trận độ cứng tổng thể [k] của kết cấu giải bằng phương pháp phần tử hữu hạn

1 XÂY DỰNG CHƯƠNG TRÌNH LẬP MA TRẬN ĐỘ CỨNG TỔNG THỂ [K] CỦA KẾT CẤU GIẢI BẰNG PHƯƠNG PHÁP PHẦN TỬ HỮU HẠN KS.PHẠM VĂN ĐOAN Email: doan281@gmail.com K19-UCT(60.58.30) Tóm tắt: Trong bài viết này, tôi sẽ giới thiệu với các bạn cách thành lập ma trận độ cứng tổng thể [K] trong bài toán giải kết cấu bằng phương pháp phần tử hữu hạn. Abstract: In this article, I will show you how to set up the overall stiffness matrix [K] in the problem resolution structures of the finite element method. I. ĐẶT VẤN ĐỀ Khi giải bài toán kết cấu bằng phương pháp phần tử hữu hạn, ta thường phải lập ma trận độ cứng [K] tổng thể của kết cấu từ các ma trận độ cứng phần tử. Nếu chỉ có một hoặc hai, thậm chí là năm ẩn thì việc thiết lập ma trận độ cứng tổng thể có thể coi là đơn giản. Nhưng khi số ẩn tính toán lên đến 20, thậm chí là 30 ẩn (ví dụ trong bài toán tấm phẳng chẳng hạn) thì việc thiết lập ma trận độ cứng tổng thể trở lên khó khăn hơn, mất thời gian hơn mà độ chính xác chắc chắn khó mà đảm bảo được khi làm bằng phương pháp thủ công. Chính vì vậy mà việc nghiên cứu và xây dựng chương trình (mà cụ thể là một hàm) có khả năng tổng hợp và thiết lập ma trận độ cứng tổng thể trở nên cần thiết cho người tính toán. II. GIẢI QUYẾT VẤN ĐỀ 1. Lựa chọn môi trường phát triển MS.Excel của MS là một trong những công cụ hàng đầu trong việc xử lý số liệu. Đặc biệt là số liệu được trình bày dưới dạng bảng, dạng cột, dạng ma trận… Chính vì vậy mà tác giả lựa chọn MS.Excel và VBA trong MS.Excel để xây dựng chương trình (hay hàm mới) có khả năng thiết lập ma trận độ cứng K tổng thể của kết cấu. 2. Cơ sở của việc lập ma trận độ cứng tổng thể Khi thiết lập ma trận độ cứng tổng thể của kết cấu, thông thường có 2 phương pháp đó là: phương pháp ma trận chỉ số và phương pháp ma trận định vị. Tuy nhiên trong bài viết này, tác giả nghiên cứu phương pháp ma trận chỉ số. Cơ sở 2 của phương pháp này như sau: ta cần phân biệt hai khái niệm chỉ số cục bộ và chỉ số tổng thể. Ø Chỉ số cục bộ là chỉ số gắn với phần tử. Ø Chỉ số tổng thể là chỉ số gắn với ẩn. Nghĩa là trong ma trận độ cứng tổng thể K sẽ có các phần tử của ma trận là: K11, K12, K13, K21, K22, K23, K31, K32 và K33. Trong ma trận độ cứng của phần tử 1 ta đã có: K111, K112, K121 và K122. Trong ma trận độ cứng của phần tử 2 ta đã có: K222, K223, K232 và K233. Như vậy ta có: Ø K11 = K111; K12 = K112; K13 = 0 Ø K21 = K121; K22 = K122 + K222; K23 = K223 Ø K31 = 0; K32 = K232; K33 = K233 Hay: [ ] ú ú ú û ù ê ê ê ë é += ú ú ú û ù ê ê ê ë é = 33 2 23 2 32 2 22 2 22 1 12 1 21 1 11 1 332313 322212 312111 0 0K KK KKKK K KKK KKK KKK K Kết luận: Như vậy, việc thành lập ma trận độ cứng K tổng thể của kết cấu có bản chất là ta phải đi tính toán từng phần tử Kij trong ma trận đó (với 1≤ i, j ≤ N: tổng số ẩn chuyển vị cần tìm). Do đó, ta chỉ việc đi xây dựng hàm tính toán từng Kij với ba tham số đầu vào là i, j và ma trận độ cứng của từng phần tử. 3 3. Thiết kế cấu trúc dữ liệu cho ma trận độ cứng của phần tử Ke và của kết cấu K a. Đối với ma trận độ cứng của phần tử Ke Để thuận tiện trong việc thiết kế thuật toán xử lý thành lập ma trận độ cứng tổng thể của kết cấu ta phải thiết kế cấu trúc dữ liệu các ma trận độ cứng của phần tử có dạng như mô hình dưới đây. Trong đó, các ẩn chuyển vị nút i và nút j của phần tử được ghi ngay lên trên ma trận độ cứng phần tử Ke theo thứ tự i, j từ trái qua phải. Phần tử 1 ẩn 1 ẩn 2 ẩn 3 … ẩn n Phần tử 2 … … … … … … … … … … … Phần tử m … … … … … Phần tử 1 Phần tử 2 … Phần tử m K11 K12 K13 … K1n ẩn 1 … K21 K22 K23 … K2n ẩn 2 … [Ke]= K31 K32 K33 … K3n ẩn 3 … … … … … … … … Kn1 Kn2 Kn3 … Knn ẩn n … Trong đó, ma trận độ cứng phần tử Ke là ma trận vuông kích thước nn ´ . Các ẩn: ẩn 1, ẩn 2, ẩn 3… ẩn n của các nút của phần tử được ghi ngay lên trên ma trận Ke. Tương tự ta sẽ đặt các ẩn này theo phương thẳng đứng ở bên phải ma trận Ke. Giả sử có m phần tử hữu hạn có cùng ma trận độ cứng Ke (phần tử 1, phần tử 2… phần tử m) nằm phía trên ma trận Ke. Khi đó, ma trận liệt kê các ẩn của các phần tử nằm phía trên ma trận Ke có kích thước là nm´ . Suy ra, ma trận Kpt bao gồm ma trận Ke và ma trận liệt kê các ẩn của các phần tử nằm trên ma trận Ke có kích thước là nnm ´+ )( . Phần tử 1 ẩn 1 ẩn 2 ẩn 3 … ẩn n Phần tử 2 … … … … … … … … … … … Phần tử m … … … … … [Kpt]= K11 K12 K13 … K1n K21 K22 K23 … K2n K31 K32 K33 … K3n … … … … … Kn1 Kn2 Kn3 … Knn Việc đặt các phần tử từ 1 đến m có thể theo trình tự từ trên xuống dưới hoặc ngược lại hoặc lộn xộn cũng được nhưng tên ẩn và vị trí đặt ẩn theo nút i, j phải chính xác. Ví dụ 1: ma trận độ cứng của phần tử thanh số 6 trong giàn phẳng được lập với mô hình trên như hình dưới đây. trong đó: (3), (4) là ẩn của nút i và (7), (8) là ẩn của nút j của phần tử hữu hạn số 6. Ma trận 4 phần tử Ke có kích thước 44´ . Ma trận Kpt có kích thước 45´ (vì thêm một dòng phần tử số 6 chứa các ẩn nên nó là 45´ ). Ví dụ 2: ma trận độ cứng của phần tử tam giác trong bài toán tấm phẳng được lập với mô hình trên như hình dưới đây. trong đó: các phần tử (1), (2), (3), (4), (5), (6) và (7) với các ẩn tương ứng được ghi lên trên ma trận Ke. Suy ra, m = 7. Vậy, ma trận phần tử Ke có kích thước 66´ . Ma trận Kpt có kích thước 613´ (vì thêm 7 dòng của 7 phần tử (từ số (1) đến số (7)) chứa các ẩn nên nó là 6136)67( ´=´+ ). b. Đối với ma trận độ cứng của kết cấu K ẩn 1 ẩn 2 ẩn 3 … ẩn N K11 K12 K13 … K1N ẩn 1 K21 K22 K23 … K2N ẩn 2 [K]= K31 K32 K33 … K3N ẩn 3 … … … … … … KN1 KN2 KN3 … KNN ẩn N Khi thành lập ma trận độ cứng của kết cấu ta thiết kế như trên. Trong đó, ma trận K là ma trận vuông kích thước NN ´ (N là tổng số ẩn chuyển vị cần tìm). Các ẩn chuyển vị: ẩn 1, ẩn 2, ẩn 3… ẩn N (gọi là các An_X, nằm ngang) của các nút được ghi ngay lên trên ma trận K. Tương tự ta cũng có các ẩn của các nút đặt thẳng đứng mà ta gọi là An_Y được ghi ở bên phải ma trận K. C ác ẩ n nà y gọ i l à: A n_ Y Các ẩn này gọi là: An_X 5 Ví dụ 3: Ma trận tổng thể của kết cấu giàn phẳng với 12 ẩn tương ứng với ma trận 1212´=´ NN như hình dưới đây. 4. Thiết kế chương trình a. Mô tả lại cách thành lập ma trận độ cứng của kết cấu K Bước 1: Xác định xem đang cần tính Kij nào của K. Ví dụ đang cần tính K37 chẳng hạn. Tương ứng với K37 là giá trị nằm ở giao của cột có chỉ số (3) và hàng có chỉ số (7) trong ma trận độ cứng của phần tử. Xem lại ví dụ 1 ta sẽ thấy K37 = -90.000. Bước 2: Dựa vào (i) và (j) của Kij ta sẽ biết được ẩn này thuộc nút nào từ đó suy ra có các phần tử nào liên quan. Trong ví dụ này chỉ có 1 phần tử số 6 là có K37. Bước 3: Khi đã xác định được các phần tử nào liên quan rồi thì ta cứ tìm Kij trong các phần tử ấy rồi cộng các giá trị đó lại thì ta được Kij của ma trận độ cứng tổng thể K của kết cấu. Trong ví dụ này, chỉ có phần tử số 6 là có K37 nên trong ma trận K ta sẽ thấy K37 = - 90.000 như trong ví dụ 3 (giao của hàng (3) và cột (7)). b. Trình tự thuật toán tính Kij Bước 1: Đọc vào 3 tham số cơ bản: An_X (chỉ số i của Kij), An_Y (chỉ số j của Kij) và ma trận Kpt. Bước 2: Tính số hàng (mRow) và số cột (mCol) của ma trận Kpt. Từ đó suy ra số phần tử hữu hạn nằm phía trên ma trận Ke là mAn = mRow – mCol. Bước 3: Nếu mAn > 0 thì làm việc sau: + Duyệt từng phần tử hữu hạn (duyệt từ 1 đến m) sau đó: +) Việc thứ nhất: Kiểm tra nếu thấy có An_X có trong danh sách các ẩn của phần tử hữu hạn đang duyệt thì đánh dấu vị trí của nó vào iX. Đồng thời, đánh dấu xác nhận là có mặt bằng biến ktX và thoát khỏi việc kiểm duyệt này vì nếu có thì chỉ có 1 ẩn giống An_X 6 thôi. +) Việc thứ hai: Kiểm tra nếu thấy có An_Y có trong danh sách các ẩn của phần tử hữu hạn đang duyệt thì đánh dấu vị trí của nó vào iY. Đồng thời, đánh dấu xác nhận là có mặt bằng biến ktY và thoát khỏi việc kiểm duyệt này vì nếu có thì chỉ có 1 ẩn giống An_Y thôi. +) Việc thứ ba: Nếu tồn tại đồng thời ktX và ktY thì tính S = S + Kij (Kij = giá trị tại giao của hàng (iY+mAn) và cột (iX) của ma trận Kpt). Ví dụ 4: Tính K37: + 3 tham số đầu vào sẽ là An_X = (3), An_Y = (7), Kpt = Ma trận 45´ . Kết hợp xem lại ví dụ 1 ta suy ra: mRow = 5, mCol = 4 và mAn = 1. + Sau khi duyệt: ta có iX = 1 (vì ẩn (3) đứng ở ngay vị trí đầu tiên) và iY=3 (vì ẩn (7) đứng ở vị trí số 3). Suy ra: iY+mAn = 3+1 = 4 và iX = 1. + Vậy, K37 sẽ là giá trị là giao của hàng 4 (tính từ (3)) với cột 1 trong ma trận 45´ của phần tử số 6. Mà giá trị này bằng -90.000 (xem ví dụ 1). c. Đoạn mã nguồn chi tiết viết bằng VBA trong MS.Excel2003 7 5. Kết quả đạt được Chỉ bằng việc viết một hàm tính các Kij là ta có thể lập được ma trận độ cứng tổng thể của kết cấu [K], xem hình như trong ví dụ 3. Trong đó, bài toán giàn phẳng có 16 phần tử với 12 ẩn. Khi sử dụng, ta chỉ vào các ma trận độ cứng của phần tử có chứa các An_X và An_Y. Tuy nhiên, để tránh sai sót ta nên dùng hàm trên với tất cả các ma trận độ cứng của các phần tử rồi cộng chúng lại. Chú ý: tên hàm là K_KC với 3 tham số đầu vào là An_X, An_Y và ma trận độ cứng của phần tử Kpt. Tải hàm này theo địa chỉ sau: III. KẾT LUẬN VÀ KIẾN NGHỊ 1. Kết luận Việc tạo ra hàm này có ý nghĩa to lớn trong việc lập ma trận độ cứng của kết cấu K. Nó vừa tính chính xác các Kij mà thời gian thực hiện lại được giảm bớt đáng kể so với việc cộng thủ công. Là một tài liệu hữu ích cho môn học Phương pháp phần tử hữu hạn giải trên MS.Excel. 2. Kiến nghị Từ kết quả này, có thể phát triển lên thành chương trình tính kết cấu bằng phương pháp phần tử hữu hạn trên MS.Excel. TÀI LIỆU THAM KHẢO [1] Nguyễn Xuân Lựu. Phương pháp phần tử hữu hạn. NXB Giao thông vận tải, 2010. [2] Chu Quốc Thắng. Phương pháp phần tử hữu hạn. NXB Khoa học và kỹ thuật, 1997. [3] Lương Xuân Bính. Bài giảng môn học Phương pháp phần tử hữu hạn ứng dụng (dành cho Cao học khóa 19). Bộ môn SBVL, Khoa Công trình, Trường ĐH.GTVT, 2011.