Bài giảng Chuyên đề Phương pháp tính - Chương 8: Phương pháp phần tử hữu hạn (Phần 2)

Sở Kỹ Thuật 
theo công thức (3.16) 
Điểm tích phân 
iξ 
Số điểm tích 
phân r 
 Trọng số wi 
 0.0000000000 Một điểm 0000000000.2
5773502692.0± Hai điểm 1 0000000000.
 0000000000.0 Ba điểm 8888888889.0
7745966692.0± 5555555555.0
3399810.0± 435 Bốn điểm 6521451548.0
8611363116.0± 0.3478548451 
 0000000000.0 0.5688888889 
5384693101.0± Năm điểm 0.4786286705 
9061798459.0± 0.2369268850 
2386191861.0± 0.4679139346 
6612093865.0± Sáu điểm 0.3607615730 
9324695142.0± 0.1713244924 
8.4 Các bước tính toán cơ bản và kỹ thuật lập trình cho máy tính số 
theo phương pháp phần tử hữu hạn 
Để áp dụng cách giải bài toán theo phương pháp phần tử hữu hạn 
người ta thực hiện các bước sau: 
- Bước 1: Rời rạc hoá miền khảo sát 
 Chia miền khảo sát V thành ne miền con V(e) hay các phần tử có 
dạng hình học nhất định. 
Ta có: V (3.18) ,V
en
1e
)e(∑
=
=
Với cách chia miền tính toán V bằng tổng các miền con V(e) , mô hình thực 
tế được thay bằng mô hình tính toán với ne phần tử hữu hạn được liên kết 
với nhau bởi các điểm nút và tại mỗi điểm nút tồn tại các đại lượng thể hiện 
sự tác động qua lại của các phần tử kề nhau, như vậy bài toán hệ liên tục có 
bậc tự do vô hạn được thay bằng bài toán tính hệ có bậc tự do hữu hạn đơn 
giản hơn nhiều. 
Ví dụ với các bài toán thấm thường có các dạng sơ đồ sau: 
Bài Giảng Chuyên Đề Phương Pháp Tính Trang 65 
Khoa Xây Dựng Thủy Lợi Thủy Điện Bộ môn Cơ Sở Kỹ Thuật 
- Một chiều: 
Nút 
Phần tử 
Lớp không thấm 
Mưa 
- Hai chiều: 
Phần tử 
Mực nước
ầ
Mặt đất 
hần Phần tửP
- Ba chiều: 
- Bước 2: Chọn hàm xấp xỉ thích hợp 
Phương pháp phần tử hữu hạn áp dụng ở đây thường là phương pháp 
Galerkin- gọi tắt là phương pháp phần tử hữu hạn Galerkin. 
Để tìm được nghiệm trên các miền con điều quan trọng là phải chọn 
hàm toạ độ Np(e) ( hay còn gọi là hàm nội suy, hàm dạng) đảm bảo sự liên 
tục của các đại lượng cần tìm giữa các phần tử trong miền D. 
 -Bước 3: Xây dựng phương trình phần tử 
Miền V được chia thành ne phần tử (miền con V(e) ) bởi R điểm nút. Tại một 
nút có s bậc tự do, thì số bậc tự do của cả hệ là: n = R.s 
Gọi { q } là véc-tơ ẩn của toàn hệ, { q }e là véc-tơ ẩn của mỗi phần tử; giả sử 
mỗi phần tử có r nút, thì số bậc tự do của mỗi phần tử là: r. s 
Ta có liên hệ { q }e = [L]e . { q } (3.19) 
Bài Giảng Chuyên Đề Phương Pháp Tính Trang 66 
Khoa Xây Dựng Thủy Lợi Thủy Điện Bộ môn Cơ Sở Kỹ Thuật 
 (ne.1) = (ne.n) x (n .1) 
Với [L]e được gọi là ma trận định vị. 
Ứng với mỗi phần tử, ta có phương trình ma trận: 
 [K]e{ }q e = {C}e (3.20) 
 [K]e ma trận phần tử , {C]e vectơ vế phải phần tử 
{q}e là tập hợp các giá trị cần tìm tại các nút của phần tử 
 -Bước 4 : Ghép nối các phần tử 
 Tập hợp cho tất cả các phần tử trong miền V, ta có: 
∑
=
ne
e 1
[K]e { }q e = {C}∑
=
ne
e 1
e 
Viết lại: [K ].{ q } = {C } (3.21) 
Trong đó: [K ] = [K]∑
=
ne
e 1
e = [L]∑
=
ne
e 1
e
T[K]e[L]e 
 { } { }∑
=
=
ne
e
eCC
1
 = ∑ [L]
=
ne
e 1
e
T{C}e 
{ }K - Ma trận tổng thể 
{ }q - Vectơ tập hợp tổng các ẩn cần tìm tại các nút (tổng bậc tự do 
tại các nút) 
 {C } Vectơ các số hạng tổng thể ở vế phải 
Như vậy việc sử dụng ma trận định vị [L]e để tính [K ] và {C }, thực chất 
là sắp xếp các phần tử [K]e , {C}e vào vị trí của nó ở trong [K ] và {C }. Tuy 
nhiên trong thực hành người ta không dùng cách nầy. 
Sau đây, sẽ giới thiệu một cách ghép nối trực tiếp để thiết lập ma trận 
tổng thể và vectơ vế phải tổng thể mà không cần xử dụng ma trận định vị 
[L]e . 
Giả sử xét bài toán thấm có áp trong miền Ω (A B C D E F), miền 
được chia thành 8 phần tử tam giác (ne =8), có 9 điểm nút (R =9), tại mỗi 
điểm nút có s bậc tự do (số ẩn số tại nút ), ở đây s =1 là cột nước thấm, mỗi 
phần tử tam giác có 3 điểm nút (r = 3); thì số bậc tự do của mỗi phần tử là: 
r ×s = 3×1 = 3 (xem Hình 3.5). 
Bài Giảng Chuyên Đề Phương Pháp Tính Trang 67 
Khoa Xây Dựng Thủy Lợi Thủy Điện Bộ môn Cơ Sở Kỹ Thuật 
X(m)
Ω 
j 
k 
i
i
Vn = 0 
 Y(m) 
8 
7
6
5
4 3 
2 
1 
9
8
7
6 
5 
4 
2 
1 
F 
3 
E D
A 
Vn = 0 
B C
 Hình 3.5: Ví dụ bài toán thấm có áp miền tính toán (ABCDEF) 
Nếu cũng với phần tử tam giác có ba điểm nút nầy r = 3, tại mỗi nút có ba ẩn 
h, u, v như bài toán dòng chảy hở hai chiều ngang s = 3, thì số bậc tự do của 
mỗi phần tử là r. s = 3x3 = 9, ta sẽ được ma trận phần tử (9,9). Để đơn giản 
ta xét phần tử tam giác tại mỗi nút có một bậc tự do. Mỗi phần tử (ở đây là 
tam giác) được đánh số các nút (i,j, k), theo chiều được qui ước (chẳng hạn 
ngược chiều kim đồng hồ), nút i được qui ước là nút ở bên trái và thấp 
nhất. Với mỗi phần tử bất kỳ ne ta có ma trận phần tử [K]e và vectơ vế phải { }C e như sau: 
 , {C}[ ]








=
e
kk
e
kj
e
ki
e
jk
e
jj
e
ji
e
ik
e
ij
e
ii
e
KKK
KKK
KKK
K e = 





e
k
e
j
e
i
c
c
c
Với cách đánh số nút và phần tử như trên ta có 8 phần tử với các nút tương 
ứng (i,j,k) như sau: e1(1,4,5), e2(1,5,2), e3(2,5,6), e4(2,6,3), e5(4,7,8), 
e6(4,8,5), e7(5,8,9), e8(5,9,6) 
Ví dụ phần tử: e4(i,j,k) e≡ 4 (2,6,3) 
 [K]e=4 = , và {C}








4
33
4
36
4
32
4
63
4
66
4
62
4
23
4
26
4
22
KKK
KKK
KKK
e=4 = 





4
3
4
6
4
2
c
c
c
Mỗi hệ số Kije chử e chỉ số trên, chỉ hệ số nầy thuộc ma trận phần tử 
nào; i là hàng nào trong ma trận tổng thể, j là cột nào trong ma trận tổng thể. 
Ví dụ K624 đây là hệ số của ma trận phần tử e = 4, nằm trong hàng 6 cột 
2 của ma trận tổng thể.và ma trận tổng thể: 
Bài Giảng Chuyên Đề Phương Pháp Tính Trang 68 
Khoa Xây Dựng Thủy Lợi Thủy Điện Bộ môn Cơ Sở Kỹ Thuật 
 [K ] = = [ ]∑
=
ne
e
eK
1
[ ]∑
=
8
1e
K e = [X] 
[X] 
11 11 12 14 15 15
21 22 22 22 23 25 25 26 26
32 33 36
41 44 44 44 45 45 47 48 48
51 51 52 52 54 54 55 55 55 55 55 55 56 56 5
1 2 2 1 1 2
2 2 3 4 4 2 3 3 4
4 4 4
1 1 5 6 1 6 5 5 6
1 2 2 3 1 6 1 2 3 6 7 8 3 8
1 2 3 4 5 6 7 8 9
1 k +k k k k +k
2 k k +k +k k k +k k +k
3 k k k
4 k k +k +k k +k k k +k
=
5 k +k k +k k +k k +k +k +k +k +k k +k k
                   
8 58 59 59
62 62 63 65 65 66 66 66 69
74 77 78
84 84 85 85 87 88 88 88 89
95 95 96 98 99 99
6 7 7 8
3 4 4 3 8 3 4 8 8
5 5
5 6 6 7 5 5 6 7 7
7 8 8 7 7 8
+k k +k
6 k +k k k +k k +k +k k
7 k k
8 k +k k +k k k +k +k k
9 k +k k k k +k
5k
..(3.22) 
Cọng một cách tương tự cho vectơ vế phải {C }, với chú ý phép cọng nầy 
giống cọng các số hạng trên đường chéo chính của ma trận tổng thể [K ]: 
{C } = {C}∑
=
ne
e 1
e (3.23) 
 Ta thấy ở ma trận tổng thể các phần tử khác không có dạng đường chéo 
(hay còn gọi là dạng Band). Để tiết kiệm bộ nhớ và thời gian tính của máy 
tính, người ta chỉ lưu trữ các phần tử khác không nầy và thuật toán cũng chỉ 
tính toán với các phần tử khác không. 
Người ta phải lưu trữ cả ma trận dạng band nầy khi ma trận band có 
chiều rộng Band hẹp (liên quan đến cách đánh số nút của các phần tử), 
không đối xứng (Hình 3.6). Chỉ cần lưu trữ một nữa band khi ma trận đối 
xứng. Khi chiều rộng Band lớn và trong các hàng của Band còn nhiều phần 
tử bằng không, người ta có thể dồn ma trận lại thành ma trận Band hẹp hơn, 
như vậy sẽ cần thêm ma trận định vị nữa. Tuy nhiên với cách lưu trữ ma trận 
Band dù theo kiểu nào, thì trong Band vẫn còn một số hệ số phần tử bằng 
không; do đó để loại bỏ các phần tử bằng không ở trong Band, người ta còn 
có cách lưu trữ các phần tử khác không nầy ở dạng vectơ gọi là kỷ thuật 
frontal method. 
Thiết lập ma trận tổng thể của bài toán ở dạng ma trận Band 
Bài Giảng Chuyên Đề Phương Pháp Tính Trang 69 
Khoa Xây Dựng Thủy Lợi Thủy Điện Bộ môn Cơ Sở Kỹ Thuật 
Ở đây ma trận tổng thể được lưu trữ ở dạng Band, ví dụ ma trận tổng thể 
không đối xứng, nên lưu trữ cả hai Band (KIJ ≠ K J I ) 
 b 
KIIKII 
2b+1
n [VK] [K]= n
n 
b 
 Hình 3.6: Cách lưu trữ ma trận dạng Band 
 Ta có KIJ = V Ki j (3.24) 
 Với i = I 
 j = J - I + 1 + b 
( Nếu ma trận đối xứng chỉ cần lưu trữ một Band, lúc đó j = J - I + 1 ) 
Sau đây là thuật toán theo phương pháp khử Gauss, viết cho ma trận Band 
đối xứng, chỉ lưu trữ một Band có chiều rộng b: 
 Ước lượng thuận Thế ngược 
- 
Bài Giảng
 Bước 5: áp đặt các điều kiện biên của bài toán ta sẽ nhận được hệ 
phương trình để giải như sau: 
 Chuyên Đề Phương Pháp Tính Trang 70 
Khoa Xây Dựng Thủy Lợi Thủy Điện Bộ môn Cơ Sở Kỹ Thuật 
Cách áp đặt điều kiện biên 
Sau khi có được ma trận hệ thống ở dạng Band, để việc lập chương trình 
được đơn giản, kích thước ma trận thổng thể của bài toán được cố định khi 
có số điều kiện biên là bất kì. 
Cách làm như sau: 
 Dạng phương trình [ K ].{ q } = { c } (3.25) 
 Nếu ẩn số thứ i = r được biết là αi , tức là: 
 qr = αi 
thì các hệ số của ma trận hệ thống được biến đổi như sau: 
 b
1
2b+1 
n [VK] = 
[K]= 
n 
n 
1 b
Hình 3.7: Cách áp đặt điều kiện biên 
 Krj = 0 nếu j ≠ r 
 Kir = 0 nếu i ≠ r (3.26) 
 Krr = 1 
 Vec-tơ vế phải của hệ thống sẽ là: 
 { = (3.27) }C














−
−
−
•
•
•
irnn
i
ir
ir
kc
kc
kc
α
α
α
α
M
M
22
11
Cũng có thể đưa điều kiện biên vào bằng cách nhân hệ số trên đường chéo 
chính của ma trận [VK] với một số rất lớn (từ 108 - 1015), khi ma trận [K] 
có tính chất trội hoặc không xấu (các hệ số kii là không quá bé so với các hệ 
số khác). 
Bài Giảng Chuyên Đề Phương Pháp Tính Trang 71 
Khoa Xây Dựng Thủy Lợi Thủy Điện Bộ môn Cơ Sở Kỹ Thuật 
- Bước 6: Giải hệ phương trình đại số 
 { }{ } { }*** CqK = (3.28) 
Cách giải hệ phương trình ở dạng ma trận (3.28) nầy tuỳ theo từng loại bài 
toán (dừng hoặc không dừng), tính chất của ma trận lưu trữ, cách lưu trữ ma 
trận tổng thể mà chọn cách giải thích hợp; chẳng hạn khử Gauss trực tiếp, 
phép tách LU hay Cholexski hoặc giải lặp Gauss-seidel có hệ số giảm dư 
hay lặp theo phương pháp gradient liên hợp, (xem N.T. Hùng, 2000) 
Bài Giảng Chuyên Đề Phương Pháp Tính Trang 72