Giáo trình Phương pháp phần tử hữu hạn - Chương 1: Phương pháp phần tử hữu hạn (Phần 3)

h phần tử hữu hạn, kể cả mô hình siêu phần tử (superelement) dùng nghiên cứu các phương tiện giao thông vận tải được giới thiệu tiếp theo như tài liệu tham khảo. Hình 7. Mô hình superelement dùng cho máy bay 35 Hình 8. Mô hình cánh máy bay Hình 9a. Mô hình FEM tính độ bền vỏ tàu vận tải Hình 9b . Mô hình tính độ bền tàu hàng rời 36 Hình 10 Kích cỡ phần tử dùng trong mô hình hóa Mô hình hóa những vật thể có hình dáng dạng điều hòa chúng ta có quyền chọn ít hoặc rất ít kiểu phần tử. Kích thước các phần tử chênh lệch nhau không nhiều. Có thể lấy dầm thành mỏng nêu tại hình 11 làm ví dụ mô hình hóa khi tính. Yêu cầu đặt ra với người tính là xác định ứng suất do uốn và do cả cắt khi dầm chịu tác động lực tập đặt tại P1 và P2. Mô hình thích hợp, dễ chấp nhận trong trường hợp cụ thể là sử dụng phần tử tấm chữ nhật, chịu uốn, kích thước các phần tử gần như nhau, hình 11b. Số lượng phần tử dùng khi tính cũng không nhất thiết quá lớn, điều này dẫn đến kích cỡ phần tử nên ở mức “thô”. Hình 11 37 Cũng bài toán tương tự, dầm côn xôn tại hình 12 cần chia thành nhiều phần tử, song kích thước các phần tử nên gần nhau. Chúng ta có thể chọn và nên chọn phần tử khối sáu mặt mô hình hóa dầm. Hình 12 Một trong những vấn đề mà bài toán trạng thái ứng suất phẳng quan tâm là tập trung ứng suất tại những vùng có kết cấu bất bình thường, ví dụ quanh mép lỗ kích thước hạn chế. Hình 13 giới thiệu hình ảnh tấm sàn chịu kéo, có lỗ khóet hình tròn. Có thể dự đoán, ứng suất tính tại vùng sát mép lỗ sẽ lớn hơn so với ứng suất tại vùng xa lỗ khoét. Mô hình hóa kết cấu này có thể nhờ các phần tử 2D trong trạng thái ứng suất phẳng. Kích thước các phần tử dùng cho trường hợp này thay đổi tùy thuộc vị trí phần tử đang chiếm giữ. Nguyên tắc chung là các phần tử nằm sát miền tập trung ứng suất có kích thước nhỏ hơn phần tử nằm xa. Điều này có thể quan sát trên hình 13b. Hình 13 Độ chính xác đạt được trong tính toán theo phương pháp PTHH không chỉ phụ thuộc vào việc chọn đúng kiểu phần tử, phù hợp mô hình toán mà còn chọn hợp lý kích cỡ phần tử. Điều chung có thể nói rằng, kích cỡ phần tử nhỏ, đồng nghĩa tăng số lượng phần tử khi mô hình hóa vật thể có thể dẫn đến độ chính xác cao hơn. Có thể nêu một ví dụ tính toán rất xa xưa nhằm tính hệ số tập trung ứng suất, minh họa tại 38 hình 14. Độ chính xác kết quả tính phụ thuộc trực tiếp số phần tử tam giác được sử dụng, được tổng kết tại bảng. Mô hình giá trị p xσ tại mép lỗ A 1,902 B 2,585 C 2,903 D 3,049 Lời giải từ phép giải tích (chính xác) 3,181 Hình 14 Đối xứng hình học, đối xứng tải Những vật thể về mặt hình học đối xứng qua trục mang tính chất đặc trưng, kết cấu hai phần, tính qua trục giống nhau song hướng ngược nhau. Trong rất nhiều trường hợp tải tác động lên kết cấu cũng mang tính đối xứng, xem hình 13. Sử dụng tính đối xứng hình học và đối xứng hoặc phản đối xứng tải trọng chúng ta có thể tháo dỡ phần đối xứng, giữ lại một phần kết cấu khi tính toán. Trường hợp tấm chịu ứng suất theo trục dọc hình 13 là một ví dụ tiêu biểu của kết cấu đối xứng. Để giảm bớt số lượng phép tính chúng ta nên tháo dỡ ½ tấm khỏi mô hình. Cách làm này được thể hiện tại hình 14. Điều kiện động học cần đàm bảo là chuyển vị theo hướng ngang các điểm vật chất trên trục đối xứng phải bằng 0. Trường hợp chung với tấm chịu uốn, ngoài điều kiện hạn chế chuyển vị tuyến tính còn đảm bảo không quay theo trục dọc của các điểm vật chất đang đề cập. Trường hợp kết cấu 2D đối xứng qua hai trục, chúng ta có quyền giữa lại chỉ ¼ kết cấu cho việc tính toán. Mô hình này chỉ phù hợp khi tải trọng tác động trực tiếp mang tính đối xứng hoặc chí ít phản đối xứng. 39 Hình 15 Hình 13 chúng ta đã tiếp xúc tấm 2D cùng lỗ khoét đối xứng qua cả hai trục Ox và Oy, tải mang tính đối xứng qua Ox. Giai đoạn đầu chúng ta đã thành lập mô hình tính dựa vào ½ tấm nguyên. Tuy nhiên nếu sử dụng đầy đủ tính đối xứng hình học và tải, chỉ cần đưa ¼ tấm vào tính toán, với điều kiện điều kiên biên trên mép trùng trục đối xứng phải đảm bảo. Theo cách làm này, khi tính ứng suất xuất hiện trong tấm dưới tac động ngọai lực, chúng ta có quyền sử dụng mô hình vừa nêu. Bài toán uốn tấm hình chữ nhật, chịu tác động lực phân bố tác động pháp tuyến được xử lý theo cách tương tự. Lợi dụng tính đối xứng hình học và tải, ¼ tấm được đưa vào tính toán nhằm xác định độ võng điểm giữa tầm. Độ võng w xác định bằng quan hệ D qLw 4 α= , trong đó α mang giá trị khác nhau, tùy thuộc số phần tử được dùng trên ¼ tấm đang đề cập. Số phần tử hệ số α 1 1,1485 4 1,2643 9 1,2653 lời giải giải tích 1,26 Dưới đây trình bày mô hình dầm công xôn thành mỏng được mô hình hóa từ các phần tử tấm chữ nhật. 40 Hình 16 4.2 Tập họp ma trận cứng của hệ thống Công thức chung trong phương pháp phần tử hữu hạn có dạng: [K}{Δ} = {P} (a) trong đó [ ] ; { } ; { } ; (b) ∑ = = E e ekK 1 ∑ = =Δ E e e 1 δ ∑ = = E e epP 1 Ví dụ 1 : Viết ma trận cứng, vector tải trọng cho dầm bị ngàm hai đầu. Dầm dài L =2l, độ cứng dầm EJ, chịu tác động lực tập trung P đặt tại giữa dầm. Mô hình hoá dầm. Chia dầm làm hai phần tử, đánh số PT 1 cho phần tử nằm bên trái, PT 2 cho phần tử còn lại. Với hai phần tử trên dầm có 3 nút, đánh theo thứ tự : nút 1 tại ngàm trái, nút 2 giữa dầm, và nút 3 tại ngàm phải. Theo cách làm này, phần tử đầu có hai nút 1 và 2; nút 2 và 3 thuộc phần tử bên phải. Tải trọng tập trung tại nút sẽ chia hai cho hai phần tử, mỗi phần tử chịu ½ P. 1 2 3 Hình 17 1 2 1 41 Trong phần mở đầu này chúng ta xác định ma trận cứng trong hệ tọa độ cục bộ song trùng hướng với hệ tọa độ chung của bài toán 1 chiều. Trong các phép tính không xuất hiện ma trận chuyển. Có thể tính các thành phần ma trận [k]e như sau. ⎥⎥ ⎥⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎣ ⎡ − −−− − − = 22 22 3 4626 612612 2646 612612 ][ llll ll llll ll l EJk e = (c) ekk kk ⎥⎦ ⎤⎢⎣ ⎡ 2221 1211 Vector lực xác định cho mỗi phần tử có dạng: { } e e f f p p p p p ⎭⎬ ⎫ ⎩⎨ ⎧= ⎪⎪⎭ ⎪⎪⎬ ⎫ ⎪⎪⎩ ⎪⎪⎨ ⎧ = 2 1 4 3 2 1 (d) Nếu thay các gá trị EJ = 1; l = 1, có thể tính: 1 2 3 4 ⎥⎦ ⎤⎢⎣ ⎡= ⎥⎥ ⎥⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎣ ⎡ − −−− − − = )1( 22 )1( 21 )1( 12 )1( 11 1 4 3 2 1 4626 612612 2646 612612 ][ kk kk k ; (e) 3 4 5 6 ⎥⎦ ⎤⎢⎣ ⎡= ⎥⎥ ⎥⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎣ ⎡ − −−− − − = )2( 22 )2( 21 )2( 12 )2( 11 2 6 5 4 3 4626 612612 2646 612612 ][ kk kk k ; (f) Trong hệ tọa độ chung tổng của hai ma trận cứng 1 và 2 sắp xếp như sau. [K] = [ hay là: ] [ ] [ ] ( ) ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢ ⎣ ⎡ +=+= )2( 22 )2( 21 )2( 12 )2( 11 )1( 22 )1( 21 )1( 12 )1( 11 21 kk kkkk kk kkK 42 1 2 3 4 5 6 [ ] [ ] [ ] 6 5 4 3 2 1 4626 612612 26446626 612661212612 2646 612612 21 ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎣ ⎡ − −−− −++− −+−+−− − − =+= kkK (g) Kết quả nhận được dạng: 1 2 3 4 5 6 [ ] 6 5 4 3 2 1 4626 612612 26)8(026 6120)24(612 2646 612612 ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎣ ⎡ − −−− − −−− − − =K (h) { } { } { } 6 5 4 3 2 1 0 0 0 0 0 6 5 4 3 0 0 0 2/ 0 2/ 0 0 2 1 2 1 21 ⎪⎪ ⎪⎪ ⎭ ⎪⎪ ⎪⎪ ⎬ ⎫ ⎪⎪ ⎪⎪ ⎩ ⎪⎪ ⎪⎪ ⎨ ⎧ += ⎪⎪⎭ ⎪⎪⎬ ⎫ ⎪⎪⎩ ⎪⎪⎨ ⎧ + ⎪⎪⎭ ⎪⎪⎬ ⎫ ⎪⎪⎩ ⎪⎪⎨ ⎧ =+= PP P P ppP 4 3 2 1 (i) Phương trình chung có dạng: ⎪⎪ ⎪⎪ ⎭ ⎪⎪ ⎪⎪ ⎬ ⎫ ⎪⎪ ⎪⎪ ⎩ ⎪⎪ ⎪⎪ ⎨ ⎧ = ⎪⎪ ⎪⎪ ⎭ ⎪⎪ ⎪⎪ ⎬ ⎫ ⎪⎪ ⎪⎪ ⎩ ⎪⎪ ⎪⎪ ⎨ ⎧ × ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎣ ⎡ − −−− − −−− − − 0 0 0 0 0 4626 612612 26)8(026 6120)24(612 2646 612612 3 3 2 2 1 1 P w w w θ θ θ (j) Điều kiện biên bài toán được thể hiện: w1 = θ1 = θ2 = w3 = θ3 = 0; Sau xử lý điều kiện biên, hệ phương trình mang dạng: 24.w2 = P (k) 43 Ví dụ 2: Lập ma trận cứng, vector tải cho giàn hình chữ L, cạnh đều, ngàm tại hai đầu. Tải trọng phân bố đều q. Độ cứng chịu uốn dầm EJ, độ cứng chịu xoắn GJT, diện tích mặt cắt ngang dầm A. Phần tử 1 ngàm tại 3, gặp phần tử thứ hai, vuông góc với phần tử đầu tại 1. Phần tử thứ hai ngàm tại 2, xem hình. Hệ toạ độ chung bố trí để trục Oy trùng với chiều dọc của phần tử thứ hai, trục Oz hướng lên trên. Hệ tọa độ cục bộ gắn liền với phần tử, trong đó trục O’x’ trùng tâm dầm. Chuyển vị mỗi nút thể hiện bằng vector {δ} = [ wj θXj θYj]T, j = 1, 2, 3. z y x 3 1 2 q Hình 18 Ma trận cứng phần tử 3-1, và vector lực tại phần tử tính trong hệ tọa độ cục bộ trùng với hệ tọa độ chung: [ ] ⎥⎦ ⎤⎢⎣ ⎡= )1( 22 )1( 21 )1( 12 )1( 11 1 kk kk k với { } 6 5 4 0 1 2 6 5 4 406 00 6012 6 2 2 23 )1( 22 ⎪⎭ ⎪⎬ ⎫ ⎪⎩ ⎪⎨ ⎧ = ⎥⎥ ⎥⎥ ⎥⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎢⎢ ⎣ ⎡ = L T qLp L EJ L EJ L GJ L EJ L EJ k (a) Với phần tử 1-2, các phép tính được thực hiện đầu tiên trong hệ tọa độ cục bộ: [ ] ⎥⎦ ⎤⎢⎣ ⎡= )2( 22 )2( 21 )2( 12 )2( 11 2 kk kk k với { } ⎪⎭ ⎪⎬ ⎫ ⎪⎩ ⎪⎨ ⎧ = ⎥⎥ ⎥⎥ ⎥⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎢⎢ ⎣ ⎡ = 6 )2( 2 2 23 )2( 11 0 1 2 406 00 6012 L T qLp L EJ L EJ L GJ L EJ L EJ k (b) Ma trận chuyển, tính cho trường hợp dầm trực giao, góc chuyển 90°. [ ] ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢ ⎣ ⎡ − = 010 100 001 R (c) 44 Tính chuyển sang hệ tọa độ chung, ma trận cứng và vector tải phần tử thứ hai tính theo công thức [k]chung = [R]{k]cục_bộ[R]T còn {p}chung = [R]T{p}cục_bộ : [ ] ⎥⎦ ⎤⎢⎣ ⎡= )2( 22 )2( 21 )2( 12 )2( 11 2 kk kk k với { } 6 5 4 0 1 2 6 5 4 00 046 0612 6 )2( 22 23 )2( 11 ⎪⎭ ⎪⎬ ⎫ ⎪⎩ ⎪⎨ ⎧ −= ⎥⎥ ⎥⎥ ⎥⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎢⎢ ⎣ ⎡ − − = L T qLp L GJ L EJ L EJ L EJ L EJ k Sau tập họp ma trận cứng và vector tải sẽ nhận được: 4 5 6 [ ] [ ] [ ] 6 5 4 0 0 )2( 22 )2( 21 )2( 12 )2( 11 )1( 22 )1( 21 )1( 12 )1( 11 21 ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢ ⎣ ⎡ +=+= kk kkkk kk kkK ch và (d) { } ⎪⎭ ⎪⎬ ⎫ ⎪⎩ ⎪⎨ ⎧ −+ ⎪⎪⎭ ⎪⎪⎬ ⎫ ⎪⎪⎩ ⎪⎪⎨ ⎧ = 0 6 1 1 2 6 1 0 1 2 pLpLp ch (e) Từ đó có thể thấy: 4 5 6 ⎪⎪ ⎪ ⎭ ⎪⎪ ⎪ ⎬ ⎫ ⎪⎪ ⎪ ⎩ ⎪⎪ ⎪ ⎨ ⎧ −= ⎪⎭ ⎪⎬ ⎫ ⎪⎩ ⎪⎨ ⎧ ⎥⎥ ⎥⎥ ⎥⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎢⎢ ⎣ ⎡ + +− − 6 6 2 2 406 046 6624 1 1 1 2 2 223 L LpL w L EJ L GJ L EJ L EJ L GJ L EJ L EJ L EJ L EJ Y X T T θ θ (f) Để ý đến quan hệ θX1 = θY1 có thể viết lại hệ phương trình trên dưới dạng: ⎪⎭ ⎪⎬ ⎫ ⎪⎩ ⎪⎨ ⎧ −=⎭⎬ ⎫ ⎩⎨ ⎧ ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢ ⎣ ⎡ +− − 3 2 28212 1224 1 1 2 23 LpLw L EJ L GJ L EJ L EJ L EJ XT θ (g) Chúng ta cùng quay lại câu chuyện tập họp ma trận cứng trong hệ tọa độ chính, trước khi giải hệ phương trình. Hình 19 trình bày kết cấu trong thực tế gồm tấm vuông, chiều dài cạnh xác định, làm việc trong trạng thái ứng suất phẳng, nới với hai thanh dạng truss, chịu lực kéo hoặc nén dọc trục, hai thanh này được nối cứng với dầm phẳng chịu uốn. 45 Hình 19 Khi thành lập ma trận cứng toàn hệ thống, động tác đầu tiên vẫn không khác cách chúng ta đã làm trong các chương trước, xác lập ma trận cứng từng phần tử trong hệ tọa độ riêng, gắn liền phần tử và chuyển ma trận cứng đó sanh hệ tọa độ chung. Có thể nhận chung khi xử lý bài toán phẳng này, hệ tọa độ riêng dùng cho các phần tử tấm, phần tử nằm ngang, ký hiệu B trên hình, của thanh chịu lực dọc và phần tử chịu uốn D trùng với hướng hệ tọa độ chung. Ma trận cứng và vector lực thanh chịu lực kéo C phải được chuyển đổi trước khi tập họp. Hệ thống đang trình bày tuy tập họp ba kiểu phần tử song cỡ bài toán rất nhỏ, ngay từ đầu bạn đọc có thể hình dung các ma trận cứng của hệ như sau. Ma trận phần tử B: 5 4 7 6 2 2 1 1 44 3433 242322 14131211 U U U U v u v u bDX bb bbb bbbb ⎥⎥ ⎥⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎣ ⎡ Ma trận phần tử C sau chuyển hệ tọa độ sang hệ tọa độ chính [cij], cỡ bằng ma trận của phần tử B. Ma trận phần tử tấm A: U2 U3 U1 U4 U5 u1 v1 u2 v2 u3 v3 u4 v4 5 4 1 3 2 4 4 3 3 2 2 1 1 88878281 78777271 68676261 2827262221 1817161211 U U U U U u u u u aaaa aaaa aaaa aaaaa aaaaa v v v v ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎢ ⎣ ⎡ L L L MM L MM L L 46 Ma trận phần tử chịu uốn D: 8 7 6 2 2 2 1 1 1 666564 565554 464544 U U Uu u ddd ddd ddd θ θ v v ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎣ ⎡ LL LL LL LL LL LL Tập họp ba nhóm ma trận cứng trong hệ thống chung [K] sẽ có dạng: 8 7 6 5 4 2 2 1 66 56552222542121 4645121244111114 424144884387 32313478337772 4241282744224321 3231181734123311 6867626166 0 0 0 0 000 U U U U U U U U d ddcbdcb ddcbdcbb bbbaba bbbabaa ccaacaca ccaacaca aaaaa ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎢ ⎣ ⎡ ++++ ++++ ++ ++ ++ ++ Cách làm vừa nêu đưa lại kết quả, ma trận cứng toàn hệ [K] được sắp xếp thành ma trận vuông, đối xứng qua đường chéo chính. Ưu điểm hàng đầu của phương pháp chuyển vị này chính ở chỗ đó. Tuy nhiên cần nói thêm, các thành phầân mang giá trị 0 trong ma trận này thường rất nhiều, tùy thuộc vào cách tổ chức đánh dấu số thứ tự bậc tự do bài toán. Thông lệ đây là ma trận vuông bậc bằng số bậc tự do. Tùy thuộc mật độ chiếm chỗ trong ma trận, điều này phụ thuộc vào đánh số thứ tự nút, có thể gặp ma trận có rất ít ô mang giá trị khác 0. Trường hợp này chúng ta gặp ma trận loãng. Trường hợp đánh số nút theo sơ đồ thỏa đáng người ta thu được ma trận [K] dạng ma trận giải hẹp. Ma trận dạng này có cấu trúc đặc biệt, đường chéo chính điền đầy, giá trị các thành phần trên đường chéo chính luôn lớn. Một số đường chéo nằm gần đường chéo chính, trên và dưới đường này điền gần đầy dữ liệu thu được từ phép cọng ma trận phần tử. Chiều rộng giải ma trận cứng tính bằng số đường chéo chứa dữ liệu khác 0. Ví dụ chứng minh cách sắp xếp này tìm thấy tại ma trận ghi ở (j) ví dụ 1 và (f) ví dụ 2 đang nêu. Từ ma trận giải hẹp, sử dụng tính đối xứng của ma trận cứng toàn hệ người ta chỉ sử dụng một nửa giải cộng với đường chéo chính khi tính toán. Cách sắp xếp ma trận cứng dưới dạng ma trận giải hẹp thuận lợi khi giải quyết những bài toán cỡ bé. 47 Thịnh hành ở USA và các nước châu Âu cách sắp xếp theo đường chân trời (sky line)1. Các thành phần có nghĩa của nửa ma trận đối xứng cùng các thành phần trên đường chéo chính ma trận [K] được chuyển địa chỉ và sắp vào vector {A}, độ lớn bằng tổng các thành phần đang nêu. Hình 20. Sắp xếp ma trận cứng trong skyline. 4.3 Xử lý điều kiện biên Xử lý điều kiện biên là điều bắt buộc trước khi giải hệ phương trình tìm chuyển vị. Những cách làm thực tế được xem xét dưới đây. 1) Thay thế điều kiện biên vào vị trí của nó trong vector chuyển vị Ma trận [K] xác lập cho hệ thống, về nguyên tắc sẽ là bậc đúng bằng bậc của độ tự do các nút. Thực tế trong tổng số bậc tự do đó, không phải tất cả đều là ẩn số. Giả sử rằng vector chuyển vị {U}, tức vector chứa ẩn theo nghĩa ban đầu gồm hai nhóm, nhóm thỏa mãn điều kiện biên, có nghĩa không còn là ẩn và nhóm thứ hai chỉ chứa ẩn của bài toán: { } ⎭⎬ ⎫ ⎩⎨ ⎧= 2 1 U U U , trong đó U2 - giá trị biết trước, U1 - ẩn cần tìm. Bài toán {K}{U} = {P} giờ có thể viết lại: (a) [ ] [ ] [ ] [ ] ⎭⎬ ⎫ ⎩⎨ ⎧= ⎭⎬ ⎫ ⎩⎨ ⎧⎥⎦ ⎤⎢⎣ ⎡ 2 1 2 1 2221 1211 P P U U KK KK Từ (a) tiến hành giải theo cách sau: [ ] [ ] { }1212111 }{}{ PUKUK =+ 1 Xem Klaus-Jürgen Bathe,1996, Finite Element Procedures, Prentice Hall Int. Editions, trang 985, và Gouri Dhatt, Gilbert Touzot, 1984, Une Présentation de la méthode des él1ments finis, Paris, trang 258-262 48 [ ] { } [ ] }{}{ 2121111 UKPUK −= (b) và [ ] [ ] { }2222112 }{}{ PUKUK T =+ (c) Trường hợp {K11] không suy biến, nghiệm U1 sẽ được tìm theo cách thức: { } [ ] { } [ ]{ }( 21211111 UKPKU −= − ) (d) Trường hợp {U1} được biết trước, hay còn nói theo cách thỏa mãn điều kiện biên, nghiệm {P2} xác định từ biểu thức (c ). Trường hợp {U2} chỉ chứa các giá trị 0, cần tiến hành loại bỏ tất cả dòng và cột của ma trận [K] tương ứng với dòng của {U2} . Cách tính tiếp tục chỉ còn là: [K11] {U1} = {P1} Các ví dụ bạn đọc tại phần trên được xử lý điều kiện biên theo cách làm này. 2) Phương pháp tăng thành phần ma trận cứng nằm trên đường chéo chính, tương ứng với chuyển vị biết trước, làm thay đổi đặc trưng ma trận, thỏa mãn điều kiện biên. Giả sử {K} đã đượcc lập theo cách vừa trình bày trên, vector chuyển vị nút cùng bậc với hệ thống phương trình, trong đó có thành phần thứ j đã xác định Uj = U*. Thay Kjj bằng Kjj + α, với α rất lớn. Thay giá trị tương ứng tại vector tải bằng giá trị αU*. Có thể viết rằng: * 1 UUKU n k kjkj αα =⎟⎠ ⎞⎜⎝ ⎛+ ∑ = Nghiệm Uj ≈ U* nếu ∑ = >> n k kjkUKU 1 *α Ví dụ minh họa cách làm trình bày tiếp theo: Hệ phương trình đại số tuyến tính thành lập từ phép tính có dạng sau: ⎪⎪⎭ ⎪⎪⎬ ⎫ ⎪⎪⎩ ⎪⎪⎨ ⎧ = ⎪⎪⎭ ⎪⎪⎬ ⎫ ⎪⎪⎩ ⎪⎪⎨ ⎧ ⎥⎥ ⎥⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎣ ⎡ 4 3 2 1 4 3 2 1 442414 3323 242322 1411 0 00 0 00 P P P P U U U U KKK KK KKK KK Điều kiện biên đặt ra U1 = U*. Các bước tiếp theo trong việc xử lý điều kiện biên: - Thay đổi thành phần K11 của ma trận cứng 49 ⎪⎪⎭ ⎪⎪⎬ ⎫ ⎪⎪⎩ ⎪⎪⎨ ⎧ × = ⎪⎪⎭ ⎪⎪⎬ ⎫ ⎪⎪⎩ ⎪⎪⎨ ⎧ ⎥⎥ ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎢ ⎣ ⎡ + 4 3 2 *15 4 3 2 1 442414 3323 242322 14 15 11 10 0 00 0 0010 P P P U U U U U KKK KK KKK KK - Chuyển thành phần thứ 1-1 tại đường chéo chính về đơn vị, có nghĩa chia tất cả thành phần trong dòng và cột chứa đại lượng vô cùng lớn cho chính đại lượng đó. ⎪⎪⎭ ⎪⎪⎬ ⎫ ⎪⎪⎩ ⎪⎪⎨ ⎧ − = ⎪⎪⎭ ⎪⎪⎬ ⎫ ⎪⎪⎩ ⎪⎪⎨ ⎧ ⎥⎥ ⎥⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎣ ⎡ * 144 3 2 * 4 3 2 1 4424 3323 242322 00 00 0 0001 UKP P P U U U U U KK KK KKK - Hệ phương trình giờ mang dạng: *; 0 0 1 * 144 3 2 4 3 2 4424 3323 242322 UU UKP P P U U U KK KK KKK = ⎪⎭ ⎪⎬ ⎫ ⎪⎩ ⎪⎨ ⎧ + = ⎪⎭ ⎪⎬ ⎫ ⎪⎩ ⎪⎨ ⎧ ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢ ⎣ ⎡ 4.4 Giải hệ phương trình đại số tuyến tính Phần việc quan trọng hàng đầu trong giải những bài toán cơ học theo phương pháp PTHH là giải hệ phương trình đại số tuyến tính vừa lập dưới dạng [K}{D} = {P}. Vì rằng ma trận {K] có những đặc tính đáng để ý nên các cách xử lý ma trận được chú ý riêng. Ma trận [K] thành lập trong khuôn khổ phương pháp PTHH có tính đối xứng, và nếu có cách tổ chức hợp lý lúc đánh số thứ tự cho các bậc tự do các nút ma trận này có giải hẹp. Các phương pháp thuận tiện dùng xử lý ma trận cứng trong quá trình giải hệ phương trình đại số tuyến tính đang được dùng rộng rãi. Phương pháp loại trừ của Gauss Dùng phương pháp loại trừ phải qua hai giai đoạn: a) chuyển hóa ma trận {K] về ma trận tam giác, gọi là quá trình “tam giác hóa”, b) giải hệ phương trình với ma trận tam giác trên vừa lập hay là quá trình lùi. Tam giác hóa ma trận K Ma trận K trong thành phần phương trình cân bằng có dạng như sau: ⎪⎪⎭ ⎪⎪⎬ ⎫ ⎪⎪⎩ ⎪⎪⎨ ⎧ = ⎪⎪⎭ ⎪⎪⎬ ⎫ ⎪⎪⎩ ⎪⎪⎨ ⎧ ⎥⎥ ⎥⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎣ ⎡ nnnnnn n n P P P D D D KKK KKK KKK MM L MLMM L L 2 1 2 1 21 22221 11211 (a) Tam giác hóa tiến hành loại trừ Dj, j = 1, 2, , n –1, tính từ phương trình j + 1 đến n. Trong quá trình này các cột, lần lượt từ j = 1 đến n – 1 bị thay bằng 0. 50 Lần thứ nhất loại trừ tiến hành theo cách cụ thể: ( )nnDKDKPKD 12121111 1 −−= L Thay D1 vào các phương trình sau 1 của (a) chúng ta có dạng sau đây của ma trận cứng sau lần loại trừ đầu tiên. ⎪⎪ ⎪ ⎭ ⎪⎪ ⎪ ⎬ ⎫ ⎪⎪ ⎪ ⎩ ⎪⎪ ⎪ ⎨ ⎧ − − = ⎪⎪⎭ ⎪⎪⎬ ⎫ ⎪⎪⎩ ⎪⎪⎨ ⎧ ⎥⎥ ⎥⎥ ⎥⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎢⎢ ⎣ ⎡ −− −− 1 11 21 1 11 21 2 1 2 1 1 11 1 12 11 1 2 1 11 21 212 11 21 22 11211 0 0 P K KP P K KP P D D D K K KKK K KK K K KKK K KK KKK nnn n nn n n nn n MM L MLMM L L Công thức chung xác định các thành phần ma trận K và vector P cho lần tính kế tiếp được viết theo cách sau: nji K PKPP K K KKK iii j iijij ,,3,2, )0( 11 )0( 1)0( 1 )0()1( )0( 11 )0( 1)0( 1 )0()1( L= ⎪⎪⎭ ⎪⎪⎬ ⎫ −= −= Sau khi loại trừ D1 D2 . . . Dn-1 sẽ nhận được công thứcsau: ⎪⎪ ⎪⎪ ⎪ ⎭ ⎪⎪ ⎪⎪ ⎪ ⎬ ⎫ ⎪⎪ ⎪⎪ ⎪ ⎩ ⎪⎪ ⎪⎪ ⎪ ⎨ ⎧ = ⎪⎪ ⎪⎪ ⎪ ⎭ ⎪⎪ ⎪⎪ ⎪ ⎬ ⎫ ⎪⎪ ⎪⎪ ⎪ ⎩ ⎪⎪ ⎪⎪ ⎪ ⎨ ⎧ ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎢ ⎣ ⎡ − − − −− 1 1 1 3 1 2 1 3 2 1 1 11 2 3 2 33 1 2 1 22 111 00000 0 00 0 00 00 0 n n s s n s n nn s sn s ss n n n P P P P P D D D D D K KK KK KK KK M M M M L MMM MOMM LM MOM LLL LLLLL Quá trình ngược bắt đầu từ phương trình cuối: 11. −− = nnnnnn PDK từ đây 1 1 − − = n nn n n n K P D Các nghiệm tìm lần lượt: ii n ij jiji i K DKP D ∑ += − = 1 51 Ví dụ 1: Giải hệ phương trình sau theo phương pháp loại trừ Gauss: ⎪⎭ ⎪⎬ ⎫ ⎪⎩ ⎪⎨ ⎧ = ⎪⎭ ⎪⎬ ⎫ ⎪⎩ ⎪⎨ ⎧ ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢ ⎣ ⎡ 180 101 34 46186 25114 842 3 2 1 D D D Bước 1: loại trừ D1 ⎪⎪ ⎪ ⎭ ⎪⎪ ⎪ ⎬ ⎫ ⎪⎪ ⎪ ⎩ ⎪⎪ ⎪ ⎨ ⎧ =− =−= ⎪⎭ ⎪⎬ ⎫ ⎪⎩ ⎪⎨ ⎧ ⎥⎥ ⎥⎥ ⎥⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎢⎢ ⎣ ⎡ =−=− =−=− 7834 2 6180 3334 2 4101 34 228 2 64664 2 6180 98 2 42534 2 4110 842 3 2 1 D D D Bước 2: loại D1 và D2 ⎪⎪⎭ ⎪⎪⎬ ⎫ ⎪⎪⎩ ⎪⎪⎨ ⎧ =− = ⎪⎭ ⎪⎬ ⎫ ⎪⎩ ⎪⎨ ⎧ ⎥⎥ ⎥⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎣ ⎡ =− 1233 3 678 34 34 49 3 62200 930 842 3 2 1 D D D Hay là: ⎪⎭ ⎪⎬ ⎫ ⎪⎩ ⎪⎨ ⎧ = ⎪⎭ ⎪⎬ ⎫ ⎪⎩ ⎪⎨ ⎧ ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢ ⎣ ⎡ 12 33 34 400 930 842 3 2 1 D D D Tìm nghiệm theo quá trình ngược sau: Từ hệ phương trình: ⎪⎭ ⎪⎬ ⎫ ⎪⎩ ⎪⎨ ⎧ = ⎪⎭ ⎪⎬ ⎫ ⎪⎩ ⎪⎨ ⎧ ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢ ⎣ ⎡ 12 33 34 400 930 842 3 2 1 D D D 3 4 12 3 ==D ( ) 23933 3 1 2 =×−=D ( ) 1382434 2 1 1 =×−×−=D Hàm viết bằng ngôn ngữ C sau đây có thể giúp bạn đọc hình dung các bước thực hiện khi thực hiện giải hệ phương trình đại số tuyến tính bằng phương pháp loại trừ lần lượt. 52 #include #include #include void gauss( double a[], double b[], int N ) { double x; int i,j,k,max,r; r = N +1; /* triangulation */ for ( i=0; i < N; i++) { max = i; for ( j = i+1; j < N; j++) if ( fabs( a[j*r+i] ) > fabs( a[max*r +i]) ) max = i; for ( k =i; k < N+1; k++) { x = a[i*r +k]; a[i*r+k] = a[max*r+k]; a[max*r+k] = x; } for ( j=i+1; j <N; j++) { for ( k=N; k >=i; k--) { a[j*r+k] = a[j*r+k] - a[i*r+k] * a[j*r+i] / a[i*r+i]; } } } /* back */ for ( j=N-1; j >=0; j--) { x = 0; for ( k=j+1; k<N; k++) x += a[j*r+k] * b[k]; b[j] = ( a[j*r+N] - x) / a[j*r +j]; } } main() { double aa[4][5] = { 1.116, 0.125, 0.14, 0.149, 1.547, 0.158, 1.167, 0.177, 0.187, 1.647, 0.197, 0.207, 1.217, 0.227, 1.747, 0.237, 0.247, 0.257, 1.267, 1.847 } ;