Bài giảng Phương pháp phần tử hữu hạn - Bài 2: Các hàm và các ma trận phần tử - Nguyễn Xuân Thành

Tóm tắt Bài giảng Phương pháp phần tử hữu hạn - Bài 2: Các hàm và các ma trận phần tử - Nguyễn Xuân Thành: ...�, 𝑧) bên trong phần tử; A(𝑥, 𝑦, 𝑧) là ma trận chứa các hàm cơ sở 𝑓𝑖(𝑥, 𝑦, 𝑧); và α là véc-tơ các hệ số 𝑎𝑖 trong các hàm chuyển vị. Hàm dạng Ma trận độ cứng Véc-tơ lực nút tương đương Một quy trình xác định hàm dạng Biểu thức trên nghiệm đúng tại các điểm, trong đó có các điểm nút c...trên. Hàm dạng Ma trận độ cứng Véc-tơ lực nút tương đương NỘI DUNG CHÍNH 1 Hàm dạng Hàm chuyển vị cho phần tử tuyến tính Một quy trình xác định hàm dạng Hàm dạng của một số phần tử thanh tuyến tính 2 Ma trận độ cứng Quy trình chung Phần tử thanh chịu kéo/nén đúng tâm Phần tử thanh chịu uốn...𝐿2 − 2𝑥 3 𝐿3 −𝑥 2 𝐿 + 𝑥3 𝐿2 ⎤⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦ 𝑇 Hàm dạng Ma trận độ cứng Véc-tơ lực nút tương đương Phần tử thanh chịu uốn ngang phẳng, nút cứng hai đầu Cụ thể: B = −𝑦 ⎡⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣ − 3 𝐿2 + 12𝑥 𝐿3 − 4 𝐿 + 6𝑥 𝐿2 6 𝐿2 − 12𝑥 𝐿3 − 2 𝐿 + 6𝑥 𝐿2 ⎤⎥⎥⎥⎥⎥⎥⎥⎥⎥...

pdf29 trang | Chia sẻ: havih72 | Lượt xem: 196 | Lượt tải: 0download
Nội dung tài liệu Bài giảng Phương pháp phần tử hữu hạn - Bài 2: Các hàm và các ma trận phần tử - Nguyễn Xuân Thành, để tải tài liệu về máy bạn click vào nút DOWNLOAD ở trên
Hàm dạng Ma trận độ cứng Véc-tơ lực nút tương đương
BÀI GIẢNG 4
(hệ thanh - phần 1/3
các hàm dạng và các ma trận phần tử)
Nguyễn Xuân Thành
tkris1004@nuce.edu.vn
Bộ môn Cơ học Kết cấu
Trường Đại học Xây dựng
Ngày 29 tháng 8 năm 2013
Hàm dạng Ma trận độ cứng Véc-tơ lực nút tương đương
NỘI DUNG CHÍNH
1 Hàm dạng
Hàm chuyển vị cho phần tử tuyến tính
Một quy trình xác định hàm dạng
Hàm dạng của một số phần tử thanh tuyến tính
2 Ma trận độ cứng
Quy trình chung
Phần tử thanh chịu kéo/nén đúng tâm
Phần tử thanh chịu uốn ngang phẳng, nút cứng hai đầu
3 Véc-tơ lực nút tương đương
Quy trình chung
Phần tử thanh chịu kéo/nén và uốn
Hàm dạng Ma trận độ cứng Véc-tơ lực nút tương đương
NỘI DUNG CHÍNH
1 Hàm dạng
Hàm chuyển vị cho phần tử tuyến tính
Một quy trình xác định hàm dạng
Hàm dạng của một số phần tử thanh tuyến tính
2 Ma trận độ cứng
Quy trình chung
Phần tử thanh chịu kéo/nén đúng tâm
Phần tử thanh chịu uốn ngang phẳng, nút cứng hai đầu
3 Véc-tơ lực nút tương đương
Quy trình chung
Phần tử thanh chịu kéo/nén và uốn
Hàm dạng Ma trận độ cứng Véc-tơ lực nút tương đương
Hàm chuyển vị cho phần tử tuyến tính
∙ ∙
∙ ∙
∙
∙ ∙
∙∙
∙
∙∙
∙∙
∙∙
∙ ∙ ∙
∙
∙
Phần tử tuyến tính
Khi chỉ có các điểm nút đặt ở đỉnh của phần tử ⇒ phần tử
"tuyến tính"*
Đơn giản cho việc thiết lập công thức các ma trận phần tử
Trong các tính toán yêu cầu độ chính xác cao, ít dùng!
Hàm dạng Ma trận độ cứng Véc-tơ lực nút tương đương
Hàm chuyển vị cho phần tử tuyến tính
Hàm chuyển vị
Thường được chọn ở dạng đa thức
Trường hợp phần tử tuyến tính, chọn các hàm chuyển vị là
các đa thức có số thành phần đúng bằng số chuyển vị tại
nút của phần tử .
Với phần tử thanh tuyến tính, bài toán 1 chiều:
𝑢(𝑥) = 𝑎1 + 𝑎2𝑥
Với phần tử 2D tam giác tuyến tính:
𝑢𝑥(𝑥, 𝑦) = 𝑎1 + 𝑎2𝑥+ 𝑎3𝑦
𝑢𝑦(𝑥, 𝑦) = 𝑎4 + 𝑎5𝑥+ 𝑎6𝑦
Hàm dạng Ma trận độ cứng Véc-tơ lực nút tương đương
Hàm chuyển vị cho phần tử tuyến tính
Hàm chuyển vị
Với phần tử 2D hình chữ nhật tuyến tính:
𝑢𝑥(𝑥, 𝑦) = 𝑎1 + 𝑎2𝑥+ 𝑎3𝑦 + 𝑎4𝑥𝑦
𝑢𝑦(𝑥, 𝑦) = 𝑎5 + 𝑎6𝑥+ 𝑎7𝑦 + 𝑎8𝑥𝑦
Với phần tử 3D hình chóp tuyến tính:
𝑢𝑥(𝑥, 𝑦, 𝑧) = 𝑎1 + 𝑎2𝑥+ 𝑎3𝑦 + 𝑎4𝑧
𝑢𝑦(𝑥, 𝑦, 𝑧) = 𝑎5 + 𝑎6𝑥+ 𝑎7𝑦 + 𝑎8𝑧
𝑢𝑧(𝑥, 𝑦, 𝑧) = 𝑎9 + 𝑎10𝑥+ 𝑎11𝑦 + 𝑎12𝑧
Hàm dạng Ma trận độ cứng Véc-tơ lực nút tương đương
Hàm chuyển vị cho phần tử tuyến tính
Hàm chuyển vị
Với phần tử 3D hình hộp tuyến tính:
𝑢𝑥(𝑥, 𝑦, 𝑧) = 𝑎1 + 𝑎2𝑥+ 𝑎3𝑦 + 𝑎4𝑧 + 𝑎5𝑥𝑦 + · · ·+ 𝑎8𝑥𝑦𝑧
𝑢𝑦(𝑥, 𝑦, 𝑧) = 𝑎9 + 𝑎10𝑥+ 𝑎11𝑦 + 𝑎12𝑧 + 𝑎13𝑥𝑦 + · · ·+ 𝑎16𝑥𝑦𝑧
𝑢𝑧(𝑥, 𝑦, 𝑧) = 𝑎17+ 𝑎18𝑥+ 𝑎19𝑦+ 𝑎20𝑧+ 𝑎21𝑥𝑦+ · · ·+ 𝑎24𝑥𝑦𝑧
Nhận xét: Điều kiện tương thích của trường chuyển vị trên
các biên chung của các phần tử hữu hạn đòi hỏi số các
thành phần của mỗi hàm chuyển vị phải bằng số nút của
của phần tử hữu hạn và ngược lại.
Chú ý! Sự có mặt của số hạng kiểu như 𝑥𝑦 hoặc 𝑥𝑦𝑧,
nhưng dạng hình học lại đặc biệt, nên ...
Hàm dạng Ma trận độ cứng Véc-tơ lực nút tương đương
Một quy trình xác định hàm dạng
Hàm dạng là hàm thể hiện mối liên hệ giữa chuyển vị tại
một điểm bất kỳ bên trong phần tử với các chuyển vị tại
nút của phần tử.
Như vậy, để xác định hàm dạng, ta xuất phát từ hàm
chuyển vị, kết hợp với điều kiện biên động học của phần tử
(tức là các chuyển vị tại nút của phần tử).
Có thể viết trường chuyển vị (tập hợp các hàm chuyển vị)
như sau:
u(𝑥, 𝑦, 𝑧) = A(𝑥, 𝑦, 𝑧)α
trong đó u là véc-tơ các chuyển vị thành phần của điểm
(𝑥, 𝑦, 𝑧) bên trong phần tử; A(𝑥, 𝑦, 𝑧) là ma trận chứa các
hàm cơ sở 𝑓𝑖(𝑥, 𝑦, 𝑧); và α là véc-tơ các hệ số 𝑎𝑖 trong các
hàm chuyển vị.
Hàm dạng Ma trận độ cứng Véc-tơ lực nút tương đương
Một quy trình xác định hàm dạng
Biểu thức trên nghiệm đúng tại các điểm, trong đó có các
điểm nút của phần tử, mà tại đó:
u(𝑥, 𝑦, 𝑧)|(𝑥,𝑦,𝑧)tại nút = q
Như vậy*:
α =
(︁
A(𝑥, 𝑦, 𝑧)|(𝑥,𝑦,𝑧)tại nút
)︁−1
q = A−1𝑒 q
Thay vào biểu thức của trường chuyển vị, có:
u(𝑥, 𝑦, 𝑧) = A(𝑥, 𝑦, 𝑧)A−1𝑒⏟ ⏞ 
N(𝑥, 𝑦, 𝑧)
q
Hàm dạng Ma trận độ cứng Véc-tơ lực nút tương đương
Hàm dạng của phần tử thanh chịu kéo/nén đúng tâm
𝑥
𝐿
𝑞4𝑞1 𝐸,𝐴∙𝑖 ∙𝑗
Hàm chuyển vị:
𝑢(𝑥) = 𝑎1 + 𝑎2𝑥 = A(𝑥)
[︂
𝑎1
𝑎2
]︂
với A(𝑥) =
[︀
1 𝑥
]︀
Xét điều kiện biên tại các đầu thanh:
q =
[︂
𝑞1
𝑞4
]︂
=
[︂
A(0)
A(𝐿)
]︂ [︂
𝑎1
𝑎2
]︂
=
[︂
1 0
1 𝐿
]︂
⏟ ⏞ 
A𝑒
[︂
𝑎1
𝑎2
]︂
Hàm dạng Ma trận độ cứng Véc-tơ lực nút tương đương
Hàm dạng của phần tử thanh chịu kéo/nén đúng tâm
Từ đó, rút ra:
N(𝑥) = A(𝑥)A−1𝑒 =
[︀
1 𝑥
]︀ [︂1 0
1 𝐿
]︂−1
=
[︁(︁
1− 𝑥
𝐿
)︁ 𝑥
𝐿
]︁
Hàm dạng Ma trận độ cứng Véc-tơ lực nút tương đương
Hàm dạng của phần tử thanh chịu uốn ngang phẳng,
nút cứng hai đầu
𝑥
𝐿
𝑞5
𝑞6
𝑞2
𝑞3
𝑦
𝐸, 𝐼∙
𝑖
∙
𝑗
Có bốn chuyển vị tại nút q =
[︀
𝑞2 𝑞3 𝑞5 𝑞6
]︀𝑇
Hàm chuyển vị:
𝑢(𝑥) = 𝑎1 + 𝑎2𝑥+ 𝑎3𝑥
2 + 𝑎4𝑥
3 = A(𝑥)
⎡⎢⎣𝑎1...
𝑎4
⎤⎥⎦
với A(𝑥) =
[︀
1 𝑥 𝑥2 𝑥3
]︀
Hàm dạng Ma trận độ cứng Véc-tơ lực nút tương đương
Hàm dạng của phần tử thanh chịu uốn ngang phẳng,
nút cứng hai đầu
Xét điều kiện biên tại các đầu thanh:
q =
⎡⎢⎢⎣
𝑞2
𝑞3
𝑞5
𝑞6
⎤⎥⎥⎦ =
⎡⎢⎢⎣
A(0)
A,𝑥 (0)
A(𝐿)
A,𝑥 (𝐿)
⎤⎥⎥⎦
⎡⎢⎢⎣
𝑎1
𝑎2
𝑎3
𝑎4
⎤⎥⎥⎦ =
⎡⎢⎢⎣
1 0 0 0
0 1 0 0
1 𝐿 𝐿2 𝐿3
0 1 2𝐿 3𝐿2
⎤⎥⎥⎦
⏟ ⏞ 
A𝑒
⎡⎢⎢⎣
𝑎1
𝑎2
𝑎3
𝑎4
⎤⎥⎥⎦
Từ đó, rút ra:
N(𝑥) = A(𝑥)A−1𝑒 =
[︀
1 𝑥 𝑥2 𝑥3
]︀ ⎡⎢⎢⎣
1 0 0 0
0 1 0 0
1 𝐿 𝐿2 𝐿3
0 1 2𝐿 3𝐿2
⎤⎥⎥⎦
−1
Hàm dạng Ma trận độ cứng Véc-tơ lực nút tương đương
Hàm dạng của phần tử thanh chịu uốn ngang phẳng,
nút cứng hai đầu
Triển khai cụ thể:
N(𝑥) =
⎡⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣
1− 3𝑥
2
𝐿2
+
2𝑥3
𝐿3
𝑥− 2𝑥
2
𝐿
+
𝑥3
𝐿2
3𝑥2
𝐿2
− 2𝑥
3
𝐿3
−𝑥
2
𝐿
+
𝑥3
𝐿2
⎤⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦
𝑇
Các trường hợp phần tử khác ...
Thực hiện theo cách tương tự trên.
Hàm dạng Ma trận độ cứng Véc-tơ lực nút tương đương
NỘI DUNG CHÍNH
1 Hàm dạng
Hàm chuyển vị cho phần tử tuyến tính
Một quy trình xác định hàm dạng
Hàm dạng của một số phần tử thanh tuyến tính
2 Ma trận độ cứng
Quy trình chung
Phần tử thanh chịu kéo/nén đúng tâm
Phần tử thanh chịu uốn ngang phẳng, nút cứng hai đầu
3 Véc-tơ lực nút tương đương
Quy trình chung
Phần tử thanh chịu kéo/nén và uốn
Hàm dạng Ma trận độ cứng Véc-tơ lực nút tương đương
Quy trình chung
Từ quan hệ chuyển vị - biến dạng và các hàm dạng đã
được xác định ở bước trên, ta đi tìm ma trận B liên hệ
giữa biến dạng với chuyển vị tại nút của phần tử
B =∇N
Từ quan hệ ứng suất - biến dạng, ta xác định ma trận E.
Trong bài toán thanh, ma trận đàn hồi E có dạng E = 𝐸
Thay vào công thức xác định ma trận độ cứng phần tử K
K =
∫︁
Ω
B𝑇EB 𝑑Ω
Hàm dạng Ma trận độ cứng Véc-tơ lực nút tương đương
Phần tử thanh chịu kéo/nén đúng tâm
𝑥
𝐿
𝑞4𝑞1 𝐸,𝐴∙𝑖 ∙𝑗
Do thành phần biến dạng duy nhất trong phần tử này chỉ
là biến dạng dọc trục 𝜀𝑥 và có liên hệ với chuyển vị dọc
trục 𝑢(𝑥) theo 𝜀𝑥 = 𝑑𝑢(𝑥)/𝑑𝑥 nên ma trận B được xác
định như sau:
B =
𝑑
𝑑𝑥
N =
𝑑
𝑑𝑥
[︁(︁
1− 𝑥
𝐿
)︁ 𝑥
𝐿
]︁
=
[︂
− 1
𝐿
1
𝐿
]︂
Hàm dạng Ma trận độ cứng Véc-tơ lực nút tương đương
Phần tử thanh chịu kéo/nén đúng tâm
Trong bài toán kéo/nén đơn, đúng tâm, không có các ứng
suất trước và biến dạng trước, thì quan hệ giữa ứng suất và
biến dạng tuân theo Định luật Hooke quen thuộc sau:
𝜎𝑥 = 𝐸𝜀𝑥
Như vậy E = 𝐸, từ đó dẫn đến:
K =
∫︁
Ω
B𝑇𝐸B 𝑑Ω = 𝐸𝐴
𝐿∫︁
0
⎡⎢⎣−
1
𝐿
1
𝐿
⎤⎥⎦[︂− 1
𝐿
1
𝐿
]︂
𝑑𝑥
hay K =
𝐸𝐴
𝐿
[︂
1 −1
−1 1
]︂
Hàm dạng Ma trận độ cứng Véc-tơ lực nút tương đương
Phần tử thanh chịu uốn ngang phẳng, nút cứng hai đầu
𝑥
𝐿
𝑞5
𝑞6
𝑞2
𝑞3
𝑦
𝐸, 𝐼∙
𝑖
∙
𝑗
Xét dầm Euler-Bernoulli, bỏ qua thành phần biến dạng
trượt, lúc này biến dạng 𝜀𝑥 ở tiết diện có tọa độ 𝑥, tại một
phân tố cách trục trung hòa đoạn bằng 𝑦 như sau:
𝜀𝑥(𝑥, 𝑦) =
𝑦
𝜌
= −𝑦𝑢,𝑥𝑥 (𝑥)
Hàm dạng Ma trận độ cứng Véc-tơ lực nút tương đương
Phần tử thanh chịu uốn ngang phẳng, nút cứng hai đầu
Như vậy, ta có ma trận B được xác định như sau:
B = −𝑦 𝑑
2
𝑑𝑥2
N = −𝑦 𝑑
2
𝑑𝑥2
⎡⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣
1− 3𝑥
2
𝐿2
+
2𝑥3
𝐿3
𝑥− 2𝑥
2
𝐿
+
𝑥3
𝐿2
3𝑥2
𝐿2
− 2𝑥
3
𝐿3
−𝑥
2
𝐿
+
𝑥3
𝐿2
⎤⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦
𝑇
Hàm dạng Ma trận độ cứng Véc-tơ lực nút tương đương
Phần tử thanh chịu uốn ngang phẳng, nút cứng hai đầu
Cụ thể:
B = −𝑦
⎡⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣
− 3
𝐿2
+
12𝑥
𝐿3
− 4
𝐿
+
6𝑥
𝐿2
6
𝐿2
− 12𝑥
𝐿3
− 2
𝐿
+
6𝑥
𝐿2
⎤⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦
𝑇
Dễ dàng xác định ma trận E = 𝐸
Hàm dạng Ma trận độ cứng Véc-tơ lực nút tương đương
Phần tử thanh chịu uốn ngang phẳng, nút cứng hai đầu
Thay vào biểu thức xác định ma trận độ cứng phần tử K
Thành phần ở hàng 𝑖 cột 𝑗 với (𝑖, 𝑗 = 1, ..., 4) của ma trận
độ cứng này được xác định như sau:
𝐾𝑖𝑗 = 𝐸
𝐿∫︁
0
⎛⎝(𝑁𝑖),𝑥𝑥 (𝑁𝑗),𝑥𝑥 ∫︁
𝐴
𝑦2𝑑𝐴
⎞⎠ 𝑑𝑥
trong đó 𝑁𝑖(𝑥) và 𝑁𝑗(𝑥) tương ứng là thành phần thứ 𝑖 và
thứ 𝑗 của véc-tơ các hàm dạng N
Hàm dạng Ma trận độ cứng Véc-tơ lực nút tương đương
Phần tử thanh chịu uốn ngang phẳng, nút cứng hai đầu
Dễ dàng nhận thấy
∫︁
𝐴
𝑦2 𝑑𝐴 là mô-men quán tính 𝐼 của
mặt cắt ngang xung quanh trục 𝑧
Thay vào trên và tiến hành các phép lấy đạo hàm và tích
phân cần thiết, ta có kết quả cuối cùng sau:
K =
𝐸𝐼
𝐿3
⎡⎢⎢⎣
12 6𝐿 −12 6𝐿
4𝐿2 −6𝐿 2𝐿2
12 −6𝐿
đx 4𝐿2
⎤⎥⎥⎦
Các trường hợp phần tử khác ...
Thực hiện theo cách tương tự trên.
Hàm dạng Ma trận độ cứng Véc-tơ lực nút tương đương
NỘI DUNG CHÍNH
1 Hàm dạng
Hàm chuyển vị cho phần tử tuyến tính
Một quy trình xác định hàm dạng
Hàm dạng của một số phần tử thanh tuyến tính
2 Ma trận độ cứng
Quy trình chung
Phần tử thanh chịu kéo/nén đúng tâm
Phần tử thanh chịu uốn ngang phẳng, nút cứng hai đầu
3 Véc-tơ lực nút tương đương
Quy trình chung
Phần tử thanh chịu kéo/nén và uốn
Hàm dạng Ma trận độ cứng Véc-tơ lực nút tương đương
Quy trình chung
Thay ma trận các hàm dạng N và ma trận B liên hệ giữa
biến dạng với chuyển vị tại nút phần tử đã tìm được ở trên,
cùng với trường tải trọng đã được cho trước vào công thức
xác định lực nút tương đương R
R𝑖 =
∫︁
B𝑇Eε0 𝑑Ω−
∫︁
B𝑇𝜎0 𝑑Ω
+
∫︁
N𝑇 f 𝑑Ω+
∫︁
N𝑇 t 𝑑𝑆 +P
Nhận xét: Không cần sử dụng đến ma trận B nếu nguyên
nhân tác dụng không bao gồm ứng suất trước và biến dạng
trước.
Hàm dạng Ma trận độ cứng Véc-tơ lực nút tương đương
Phần tử thanh chịu kéo/nén và uốn
𝑥
𝐿
𝑥
𝑅4𝑅1 𝐸,𝐴, 𝐼
𝑅5
𝑅6
𝑅2
𝑅3
𝑦
∙
𝑖
∙
𝑗
𝑝𝑥(𝑥)
𝑝𝑦(𝑥)
Do phần tử chịu đồng thời cả kéo/nén và uốn, nên ta xét
trường hợp tổng quát là mỗi đầu nút phần tử có 3 chuyển
vị tại nút. Như vậy, véc-tơ chuyển vị tại nút tổng cộng có 6
thành phần
q =
[︀
𝑞1 𝑞2 . . . 𝑞6
]︀𝑇
Hàm dạng Ma trận độ cứng Véc-tơ lực nút tương đương
Phần tử thanh chịu kéo/nén và uốn
Tương ứng với nó là véc-tơ lực nút tương đương cũng có 6
thành phần
R =
[︀
𝑅1 𝑅2 . . . 𝑅6
]︀𝑇
Chuyển vị tại một điểm bất kỳ có tọa độ 𝑥:
u =
[︂
𝑢𝑥(𝑥)
𝑢𝑦(𝑥)
]︂
Ma trận các hàm dạng sau khi kể đến tất cả các chuyển vị
tại nút:
N =
[︂
𝑁1(𝑥) 0 0 𝑁4(𝑥) 0 0
0 𝑁2(𝑥) 𝑁3(𝑥) 0 𝑁5(𝑥) 𝑁6(𝑥)
]︂
Hàm dạng Ma trận độ cứng Véc-tơ lực nút tương đương
Phần tử thanh chịu kéo/nén và uốn
Trường tải trọng tác dụng:
p =
[︂
𝑝𝑥(𝑥)
𝑝𝑦(𝑥)
]︂
Do không xét đến các nguyên nhân ứng suất ban đầu và
biến dạng ban đầu, công thức xác định lực nút tương
đương sẽ chỉ là:
R =
∫︁
Ω
N𝑇p 𝑑Ω
Hàm dạng Ma trận độ cứng Véc-tơ lực nút tương đương
Phần tử thanh chịu kéo/nén và uốn
Cụ thể các thành phần của véc-tơ tải trọng như sau:
𝑅1 =
𝐿∫︁
0
𝑁1(𝑥)𝑝𝑥(𝑥) 𝑑𝑥
𝑅2 =
𝐿∫︁
0
𝑁2(𝑥)𝑝𝑦(𝑥) 𝑑𝑥
và cứ như vậy ...
Một số trường hợp đặc biệt của tải trọng ...

File đính kèm:

  • pdfbai_giang_phuong_phap_phan_tu_huu_han_bai_2_cac_ham_va_cac_m.pdf
Ebook liên quan