i
PHƯƠNG PHÁP PHẦN
TỬ HỮU HẠN
Lý thuyết
Bài tập
Chương trình MATLAB
THÁI NGUYÊN 2011
NGÔ NHƯ KHOA
Ngô Như Khoa
PHƯƠNG PHÁP PHẦN TỬ HỮU HẠN

THÁI NGUYÊN 2011
MỞ ĐẦU
Giáo trình Phương pháp Phần tử hữu hạn (PP PTHH) được biên soạn dựa trên
cuốn: Giáo trình Phương pháp Phần tử hữa hạn – Lý thuyết, bài tập và chương trình
Matlab. GS.TS. Trần Ích Thịnh, TS. Ngô Như Khoa. NXB Khoa học Kỹ thuật 2007.
Và kinh nghiệm giảng dạy môn học cùng tên trong những năm gần đây cho sinh viên
khoa Cơ khí, trường Đại học Bách khoa Hà Nội và học viên cao học ngành Cơ học Kỹ
thuật, trường Đại học Kỹ thuật Công nghiệp - Đại học Thái Nguyên. Nội dung giáo
trình có mục đích trang bị cho sinh viên các ngành kỹ thuật: Công nghệ chế tạo máy,
Kỹ thuật cơ khí, v.v. Với các nội dung:
- Những kiến thức cơ bản nhất của PP PTHH ứng dụng,
- Áp dụng phương pháp để giải quyết một số bài toán kỹ thuật khác nhau,
- Nâng cao kỹ năng lập trình Matlab trên cơ sở thuật toán PTHH.
Giáo trình biên soạn gồm 11 chương.
Sau phần giới thiệu phương pháp PTHH, một số loại phần tử thực và phần tử qui
chiếu hay gặp (Chương 1), giáo trình đề cập đến một số phép tính ma trận, phương
pháp khử Gauss (Chương 2) và thuật toán xây dựng ma trận độ cứng và véctơ lực nút
chung cho kết cấu (Chương 3). Phương pháp Phần tử hữu hạn trong bài toán một
chiều chịu kéo (nén) được giới thiệu trong Chương 4 và ứng dụng vào tính toán hệ
thanh phẳng (Chương 5). Tiếp theo, giáo trình tập trung vào mô tả phần tử hữu hạn

1.4. Nhân ma trận với hằng số 10
1.5. Nhân hai ma trận 10
1.6. Chuyển vị ma trận 11
1.7. Đạo hàm và tích phân ma trận 11
1.8. Định thức của ma trận 11
1.9. Nghịch đảo ma trận 12
1.10. Ma trận đường chéo 12
1.11. Ma trận đối xứng 12
1.12. Ma trận tam giác 13
2. PHÉP KHỬ GAUSS 13
2.1. Mô tả 13
2.2. Giải thuật khử Gauss tổng quát 14
Chương 3 17
THUẬT TOÁN XÂY DỰNG MA TRẬN ĐỘ CỨNG CHUNG 17
VÀ VÉCTƠ LỰC NÚT CHUNG 17
1. CÁC VÍ DỤ 17
1.1. Ví dụ 1 17
1.2. Ví dụ 2 19
2. THUẬT TOÁN GHÉP K VÀ F 21
2.1. Nguyên tắc chung 21
2.2. Thuật toán ghép nối phần tử: 23
Chương 4 24
PHẦN TỬ HỮU HẠN TRONG BÀI TOÁN MỘT CHIỀU 24
1. MỞ ĐẦU 24
2. MÔ HÌNH PHẦN TỬ HỮU HẠN 24
3. CÁC HỆ TRỤC TOẠ ĐỘ VÀ HÀM DẠNG 25
4. THẾ NĂNG TOÀN PHẦN 27
5. MA TRẬN ĐỘ CỨNG PHẦN TỬ 28
6. QUI ĐỔI LỰC VỀ NÚT 28
7. ĐIỀU KIỆN BIÊN, HỆ PHƯƠNG TRÌNH PHẦN TỬ HỮU HẠN 30

Chương 8 74
PHẦN TỬ TỨ GIÁC 74
1. MỞ ĐẦU 74
5.1. PHẦN TỬ TỨ GIÁC 74
5.2. HÀM DẠNG 74
5.3. MA TRẬN ĐỘ CỨNG CỦA PHẦN TỬ 76
5.4. QUI ĐỔI LỰC VỀ NÚT 78
5.7. TÍCH PHÂN SỐ 78
5.8. TÍNH ỨNG SUẤT 81
5.9. VÍ DỤ 82
BÀI TẬP CHƯƠNG 8 83
Chương 9 84
PHẦN TỬ HỮU HẠN TRONG TÍNH TOÁN KẾT CẤU DẦM VÀ KHUNG 84
1. GIỚI THIỆU 84
5.10. THẾ NĂNG 84
5.11. HÀM DẠNG HERMITE 85
5.12. MA TRẬN ĐỘ CỨNG CỦA PHẦN TỬ DẦM 86
5.13. QUY ĐỔI LỰC NÚT 88
5.14. TÍNH MÔMEN UỐN VÀ LỰC CẮT 89
5.15. KHUNG PHẲNG 89
5.16. VÍ DỤ 91
BÀI TẬP CHƯƠNG 9 95
Chương 10 97
PHẦN TỬ HỮU HẠN TRONG BÀI TOÁN DẪN NHIỆT 97
1. GIỚI THIỆU 97
2. BÀI TOÁN DẪN NHIỆT MỘT CHIỀU 97
iii
2.1. Mô tả bài toán 97
2.2. Phần tử một chiều 97
2.3. Ví dụ 98

Để có thể khai thác hiệu quả những phần mềm PTHH hiện có hoặc tự xây dựng
lấy một chương trình tính toán bằng PTHH, ta cần phải nắm được cơ sở lý thuyết, kỹ
thuật mô hình hoá cũng như các bước tính cơ bản của phương pháp.
5.2.XẤP XỈ BẰNG PHẦN TỬ HỮU HẠN
Giả sử V là miền xác định của một đại lượng cần khảo sát nào đó (chuyển vị, ứng
suất, biến dạng, nhiệt độ, v.v.). Ta chia V ra nhiều miền con v
e
có kích thước và bậc tự
do hữu hạn. Đại lượng xấp xỉ của đại lượng trên sẽ được tính trong tập hợp các miền
v
e
.
Phương pháp xấp xỉ nhờ các miền con v
e
được gọi là phương pháp xấp xỉ bằng các
phần tử hữu hạn, nó có một số đặc điểm sau:
- Xấp xỉ nút trên mỗi miền con v
e
chỉ liên quan đến những biến nút gắn vào nút
của v
e
và biên của nó,
- Các hàm xấp xỉ trong mỗi miền con v
e
được xây dựng sao cho chúng liên tục
trên v
e
và phải thoả mãn các điều kiện liên tục giữa các miền con khác nhau.
- Các miền con v
e

hạn hay gặp.
Phần tử một chiều
Phần tử hai chiều
Phần tử ba chiều
Phần tử tứ diện
2
biên giới
biên giới
v
2
v
1
biên giới
v
2
v
1
v
1
v
2
Hình 1.1. Các dạng biên chung giữa các phần tử
Phần tử bậc nhất
Phần tử bậc hai
Phần tử bậc ba
Phần tử bậc nhất
Phần tử bậc hai
Phần tử bậc ba
5.5. PHẦN TỬ QUY CHIẾU, PHẦN TỬ THỰC
Với mục đích đơn giản hoá việc xác định giải tích các phần tử có dạng phức tạp,

3
v
2
v
1
1,00,0
y
ξ
x
(1)
(2)
(3)
(4)
(5)
η
r
3
r
2
r
1
0,1
Hình 1.2. Phần tử quy chiếu và các phần tử thực tam giác
- Một phần tử qui chiếu v
r
được biến đổi thành tất cả các phần tử thực v
e
cùng loại nhờ các phép biến đổi khác nhau. Vì vậy, phần tử qui chiếu còn được
gọi là phần tử bố-mẹ.
- Có thể coi phép biến đổi hình học nói trên như một phép đổi biến đơn giản.

v
r
10,0
1
ξ
v
r
10,0
1
η
η η
1
/
2
,1
/
2
1
/
2
1
/
2
1
/
3
,
2
/
3

r
0,1,0
0,0,1
v
r
ζ
η
η
1,0,0
ζ
1,0,0
ξ
η
ζ
0,1,0
1,0,0
0,0,1
v
r
ξ
Phần tử bậc nhất Phần tử bậc hai Phần tử bậc ba
0,1,1
ξ
v
r
ζ
η
η
1,1,0
ζ

T
- Lực tập trung P
i
:

P
i
= P
i
[ P
x
, P
y
, P
z
]
T
Chuyển vị của một điểm thuộc vật được ký hiệu bởi:
u = [u, v, w]
T
(1.1)
Các thành phần của tenxơ biến dạng được ký hiệu bởi ma trận cột:
ε
= [
ε
x
,
ε
y
,

y
v
x
u






∂
∂
+
∂
∂
∂
∂
+
∂
∂
∂
∂
+
∂
∂
∂
∂
∂
∂
∂

= D
ε
(1.5)
Trong đó:
( )( )




















−
−
−
−
−

Với vật thể đàn hồi tuyến tính thì năng lượng biến dạng trên một đơn vị thể tích
được xác định bởi:
εσ
T
2
1

Do đó năng lượng biến dạng toàn phần:
∫
=
V
T
dvU
εσ
2
1
(1.7)
6
Công của ngoại lực được xác định bởi:
∑
∫∫
=
−−−=
n
i
i
T
i
S
T

i
Áp dụng nguyên lý cực tiểu thế năng: Đối với một hệ bảo toàn, trong tất cả các di
chuyển khả dĩ, di chuyển thực ứng với trạng thái cân bằng sẽ làm cho thế năng đạt
cực trị. Khi thế năng đạt giá trị cực tiểu thì vật (hệ) ở trạng thái cân bằng ổn định.
5.9.SƠ ĐỒ TÍNH TOÁN BẰNG PHƯƠNG PHÁP PHẦN TỬ HỮU HẠN
Một chương trình tính bằng PTHH thường gồm các khối chính sau:
Khối 1: Đọc các dữ liệu đầu vào: Các dữ liệu này bao gồm các thông tin mô tả nút và
phần tử (lưới phần tử), các thông số cơ học của vật liệu (môđun đàn hồi, hệ
số dẫn nhiệt ), các thông tin về tải trọng tác dụng và thông tin về liên kết
của kết cấu (điều kiện biên);
Khối 2: Tính toán ma trận độ cứng phần tử k và véctơ lực nút phần tử f của mỗi phần
tử;
Khối 3: Xây dựng ma trận độ cứng tổng thể K và véctơ lực nút F chung cho cả hệ
(ghép nối phần tử);
Khối 4: Áp đặt các điều kiện liên kết trên biên kết cấu, bằng cách biến đổi ma trận độ
cứng K và vec tơ lực nút tổng thể F;
Khối 5: Giải phương trình PTHH, xác định nghiệm của hệ là véctơ chuyển vị chung
Q;
Khối 6: Tính toán các đại lượng khác (ứng suất, biến dạng, gradiên nhiệt độ, v.v.) ;
Khối 7: Tổ chức lưu trữ kết quả và in kết quả, vẽ các biểu đồ, đồ thị của các đại lượng
theo yêu cầu.
Sơ đồ tính toán với các khối trên được biểu diễn như hình sau (Hình 1.3);
7
8
Tính toán ma trận độ cứng phần tử k
Tính toán véctơ lực nút phần tử f
Giải hệ phương trình KQ = F
(Xác định véctơ chuyển vị nút tổng thể Q)
Đọc dữ liệu đầu vào
- Các thông số cơ học của vật liệu

=++
=++




2211
22222121
11212111
(2.1)
trong đó, x
1
, x
2
, …, x
n
là các nghiệm cần tìm. Hệ phương trình (2.1) có thể được biểu
diễn ở dạng thu gọn:
Ax = b (2.2)
trong đó, A là ma trận vuông có kích thước (n
×
n), và x và b là các véctơ (n
×
1), được
biển diễn như sau:












=
n
x
x
x
x

2
1














=

=
34
2
11
c
1.2. Ma trận đơn vị
Ma trận đơn vị là ma trận đường chéo với các phần tử trên đường chéo chính bằng
1, ví dụ:
9










=
100
010
001
I
1.3. Phép cộng và phép trừ ma trận.
Cho 2 ma trận A và B, cùng có kích thước là (m
×
n). Tổng của chúng là 1 ma trận
C = A + B và được định nghĩa như sau:
c


−
34
75
21
58
15
23
phép trừ được định nghĩa tương tự.
1.4. Nhân ma trận với hằng số
Nhân 1 ma trận A với hằng số c được định nghĩa như sau:
cA=[ca
ij
] (2.4)
Ví dụ:






−
=






−

=
n
k
kjikij
bac
1
(2.6)
Ví dụ:






=










×





] kích thước (m
×
n) là 1 ma trận, ký hiệu là A
T
có
kích thước là (n
×
m), được tạo từ ma trận A bằng cách chuyển hàng của ma trận A
thành cột của ma trận A
T
. Khi đó, (A
T
)
T
= A.
Chuyển vị của một tích các ma trận là tích các chuyển vị của ma trận thành phần
theo thứ tự đảo ngược, có nghĩa là:
(A
×
B
×
C)
T
=C
T
×
B
T
×
A







=
dx
xda
xA
dx
d
ij
)(
)(
(2.8)
[ ]
∫ ∫
=
dxdyaAdxdy
ij
(2.9)
Chúng ta sẽ sử dụng thường xuyên biểu thức (2.8) để xây dựng hệ phương trình
PTHH trong các chương sau. Xét ma trận vuông A, kích thước (n
×
n) với các hệ số
hằng, véctơ cột x = {x
1
x
2

n
n i j
n n ij ij
j
A a A a A a A a A
+ +
=
= − + − = −
∑
L
(2.11)
trong đó, A
ij
là ma trận kích thước (n-1
×
n-1) thu được bằng cách loại đi hàng i cột j
của ma trận A.
11
Ví dụ:



















32
33332
22322
11
21
22221
11211
Công thức (2.11) là công thức tổng quát. Theo công thức này, định thức của ma trận
vuông có kích thước (n
×
n) được xác định theo phương pháp truy hồi từ định thức các
ma trận có kích thước (n-1
×
n-1). Trong đó, ma trận chỉ có 1 phần tử (1
×
1) có:
det(a
pq
) = a
pq
(2.12)
1.9. Nghịch đảo ma trận
Cho ma trận vuông A, nếu det(A) ≠ 0, thì A có ma trận nghịch đảo và ký hiệu là A

Aa
+
−=
và A
ji
là ma trận
thu được từ A bằng cách loại đi hàng thứ j và cột thứ i.
Ví dụ:
Nghịch đảo của ma trận A kích thước (2
×
2) là:






−
−
=






=
−
−
1121

D
1.11. Ma trận đối xứng
Ma trận đối xứng là một ma trận vuông có các phần tử thoả mãn điều kiện:
a
ij
= a
ji
hay: A = A
T
(2.15)
Như vậy, ma trận đối xứng là ma trận có các phần tử đối xứng qua đường chéo chính.
12
1.12. Ma trận tam giác
Ma trận được gọi là ma trận tam giác trên hay ma trận tam giác dưới, tương ứng là
các ma trận có tất cả các phần tử nằm dưới hay nằm trên đường chéo chính bằng
không.
Ví dụ, các ma trận được minh hoạ dưới đây tương ứng là ma trận tam giác trên A
và ma trận tam giác dưới B:










−
−

-1
và nhận được nghiệm: x =
A
-1
b. Tuy nhiên, trong hầu hết các bài toán kỹ thuật, kích thước của ma trận A là rất
lớn và các phần tử của A thường là số thực với miền xác định rất rộng; do đó, việc tính
toán ma trận nghịch đảo của A là rất phức tạp và dễ gặp phải sai số do việc làm tròn
trong các phép tính. Vì vậy, phương pháp khử Gauss là một công cụ rất hữu ích cho
việc giải hệ phương trình đại số tuyến tính.
2.1. Mô tả
Chúng ta sẽ bắt đầu mô tả phương pháp khử Gauss thông qua một ví dụ minh hoạ sau
đây; sau đó tìm hiểu giải thuật khử Gauss tổng quát.
Xét hệ phương trình:
152
321
=++
xxx
(1)
2352
321
−=++
xxx
(2)
415
321
=+−−
xxx
(3)
Bước 1: bằng các phép biến đổi tương đương để khử x
1

(1)
470
321
=+−
xxx
(2
1
)
92700
321
=++
xxx
(3
2
)
13
Ở đây, ta nhận được hệ phương trình mà ma trận các hệ số lập thành ma trận tam
giác trên. Từ phương trình cuối cùng (3
2
), ta tìm được nghiệm x
3
, lần lượt thế các
nghiệm tìm được vào phương trình trên nó, (2
1
) và (1). Sẽ nhận được các ẩn số cần tìm
như sau:
3
8
;
3


−⇒










−−
−
92700
4710
1521
52010
4710
1521
41511
2352
1521
bằng phương pháp thế ngược, cuối cùng ta nhận được các nghiệm:
3
8
;
3
5
;

=






























x
x
x
x
aaaaa
aaaaa
aaaaa
aaaaa
aaaaa











3
2
1
3
2
1
321
321
33333231
22232221

nj
aaaaa
aaaaa
aaaaa
aaaaa
aaaaa







321
321
33333231
22232221
11131211













14
( )
( )







=−=
−=
njib
a
a
bb
a
a
a
aa
i
ii
j
i
ijij
, ,2,;
1
11
1
1







111
3
1
2
111
3
1
2
1
3
1
3
1
33
1
32
1
2
1
2
1
23
1
22




















1
1
1
3
1
2
1
n
i
b
b












−−−
+
−−−
+
−
+
−
+
−
++
1
,
1
,
1
1,
1
,
1

nn
k
jn
k
kn
k
ni
k
ji
k
ki
k
nk
k
jk
k
kk
nj
nj
nj
aaa
aaa
aaa
aaa
aaaa
aaaaa
















−
−
−
+
1
1
1
1
3
3
1
2
1
k
n
k
i
k
k

−
−
−
−
−
−
−
−
nkjib
a
a
bb
nkjia
a
a
aa
k
k
k
kk
k
ik
k
i
k
i
k
kj
k
kk











=




















1
4
3
2
1
)1(
)3(
4
)3(
44
)2(
3
)2(
34
)2(
33
)1(
2
)1(
24
)1(
23
)1(
22
114131211
0
n
nn
n
nn

1
,,n,ni;
a
xab
x,;
a
b
x
ii
n
ij
jiji
i
nn
n
n

−−=
−
==
∑
+=
(2.23)
16
Chương 3
THUẬT TOÁN XÂY DỰNG MA TRẬN ĐỘ CỨNG CHUNG
VÀ VÉCTƠ LỰC NÚT CHUNG
Việc ghép các ma trận độ cứng k và các véctơ lực f của các phần tử để tạo ra ma
trận độ cứng K và véctơ lực nút F chung cho cả hệ, từ đó thiết lập hệ phương trình
PTHH là một vấn đề quan trọng.

k
;










=
432
371
218
2
k
;










=

1 2 3
1 1 2 4
2 4 2 5
3 2 3 5
2. Xét từng phần tử
Với phần tử 1, các dòng và cột được nhận dạng như sau:
4
2
1
521
263
137
421
1










=k
Ma trận này được cộng vào ma trận độ cứng chung và ta sẽ được:









=K
Ma trận độ cứng của phần tử 2 được gán số bởi:
5
2
4
432
371
218
524
2










=k
Các số hạng của ma trận k
2
được cộng thêm vào ma trận chung, cho ta









++
++
=K
Với phần tử 3:
18
5
3
2
501
064
149
532
3










=k












+++
+
++
=K
Việc cộng các véctơ lực phần tử vào véctơ lực chung được tiến hành hoàn toàn
tương tự.
1.2. Ví dụ 2
Giả sử có hai phần tử 1 và 4 trong bài toán hai chiều; mỗi phần tử có 3 nút, mỗi
nút có 2 bậc tự do (Hình 3.2). Hãy mô tả quá trình ghép nối ma trận độ cứng chung K
và véctơ lực nút chung, theo các ma trận độ cứng phần tử và véctơ lực nút phần tử k
1
,
k
4
, f
1
và f
4
cho trước như sau:

5363097
7199293
2647322
1
k
;






















=
5

−−−−−
−−−−−
−−−−−
=
2874755
7272873
4225768
7873026
5742191
5386123
4
k
;
















