9/3/2019
1
Khoa Khoa học ứng dụng
Bộ môn Cơ kỹ thuật
Đại học Quốc gia Thành phố Hồ Chí Minh
Trường đại học Bách Khoa
Lập trình tính toán Cơ kỹ thuật Matlab® – Matrix Laboratory 1
Giảng viên: Ths. Trần Kim Bằng
Bộ môn Cơ Kỹ Thuật, P.106B4
Khoa Khoa Học Ứng Dụng
Đại học Bách Khoa TpHCM
Khoa Khoa học ứng dụng
Bộ môn Cơ kỹ thuật
Đại học Quốc gia Thành phố Hồ Chí Minh
Trường đại học Bách Khoa
Lập trình tính toán Cơ kỹ thuật Matlab® – Matrix Laboratory 2
Chương 6
Lập trình PPPTHH cho
bài toán khung phẳng
Khoa Khoa học ứng dụng
Bộ môn Cơ kỹ thuật
Đại học Quốc gia Thành phố Hồ Chí Minh
Trường đại học Bách Khoa
Lập trình tính toán Cơ kỹ thuật
Matlab® – Matrix Laboratory 3
6.1. Tiền xử lý (Preprocessor)
Khoa Khoa học ứng dụng
Bộ môn Cơ kỹ thuật
Đại học Quốc gia Thành phố Hồ Chí Minh
Trường đại học Bách Khoa
Lập trình tính toán Cơ kỹ thuật
Matlab® – Matrix Laboratory 4
6.1.1. Mô hình cơ học của kết cấu khung phẳng
E = 2.1011 N/m2; A = 100 cm2; J = 1000 cm4;
l = 1 m; p0 = q0 = 10KN/m
q0
l
l p0 l/2
q0l2
A B
C
D
Khoa Khoa học ứng dụng
Bộ môn Cơ kỹ thuật
Đại học Quốc gia Thành phố Hồ Chí Minh
Trường đại học Bách Khoa
Lập trình tính toán Cơ kỹ thuật Matlab® – Matrix Laboratory 5
Nhập moment quán tính
Đưa kết quả số về dạng 4 số lẻ thập
phân nhân 10n
Nhập mô-đun đàn hồi
Nhập diện tích mặt cắt ngang
clc
clear all
close all
format short
E = 2e11;
A = 0.01;
J = 1e-5;
btd_nut = 3;
Nhập số bậc tự do của một nút
6.1.2. Nhập các hằng số như vật liệu, diện tích mặt cắt ngang
Khoa Khoa học ứng dụng
Bộ môn Cơ kỹ thuật
Đại học Quốc gia Thành phố Hồ Chí Minh
Trường đại học Bách Khoa
Lập trình tính toán Cơ kỹ thuật Matlab® – Matrix Laboratory 6
6.1.3. Rời rạc hóa kết cấu thành mô hình phần tử hữu hạn
i
q
O x
y
Chọn gốc
tọa độ của
hệ trục tổng
thể đặt tại
nút 1 theo
chỉ số tổng
thể.
là bậc tự do thứ i của kết
cấu theo chỉ số tổng thể.
7
q
8
q
9
q
4
q
5
q
6
q
1
q
2
q
3
q
10
q
12
q
1 2
3
4
(1)
(2 )
(3)
9/3/2019
2
Khoa Khoa học ứng dụng
Bộ môn Cơ kỹ thuật
Đại học Quốc gia Thành phố Hồ Chí Minh
Trường đại học Bách Khoa
Lập trình tính toán Cơ kỹ thuật Matlab® – Matrix Laboratory 7
e
i
q
là bậc tự do thứ i của phần tử thứ (e) theo chỉ
số địa phương khi xét trên một phần tử.
( 2 )
90
o
(2 )
( 2 )
1
q
( 2 )
3
q
( 2 )
2
q
( 2 )
4
q
( 2 )
6
q
( 2 )
5
q
,
x
,
y
(1 )
0
o
(1)
(1)
4
q
(1)
6
q
(1)
5
q
(1)
1
q
(1)
3
q
(1)
2
q
,
x
,
y
( 3 )
90
o
(3)
( 3 )
1
q
( 3 )
3
q
( 3 )
2
q
( 3 )
4
q
( 3 )
6
q
( 3 )
5
q
,
x
,
y
Khoa Khoa học ứng dụng
Bộ môn Cơ kỹ thuật
Đại học Quốc gia Thành phố Hồ Chí Minh
Trường đại học Bách Khoa
Lập trình tính toán Cơ kỹ thuật Matlab® – Matrix Laboratory 8
• Ma trận toado là ma trận có n hàng, 2 cột.
• Với n là số nút.
• Hàng 1 chứa tọa độ của nút thứ 1
• Hàng n chứa tọa độ của nút thứ n
• Cột 1 chứa tọa độ x, cột 2 chứa tọa độ y của nút.
Nút 1
Nút 4
Code Matlab:
toado(1,1)=0;
toado(1,2)=0;
toado(2,1)=1;
toado(2,2)=0;
toado(3,1)=1;
toado(3,2)=-0.5;
toado(4,1)=1;
toado(4,2)=-1;
toado =
0 0
1 0
1 -0.5
1 -1
6.1.4. Tạo ma trận chứa tọa độ các nút
Khoa Khoa học ứng dụng
Bộ môn Cơ kỹ thuật
Đại học Quốc gia Thành phố Hồ Chí Minh
Trường đại học Bách Khoa
Lập trình tính toán Cơ kỹ thuật
Matlab® – Matrix Laboratory
phantu =
1 2
2 3
3 4
9
• Ma trận phantu là ma trận có n hàng 2 cột, n là số phần tử.
• Hàng 1 chứa số thứ tự tổng thể của nút đầu và cuối thuộc phần tử thứ 1.
• Hàng n chứa số thứ tự tổng thể của nút đầu và cuối thuộc phần tử thứ n.
• Cột 1 chứa số thứ tự tổng thể của tất cả nút đầu của tất cả phần tử.
• Cột 2 chứa số thứ tự tổng thể của tất cả nút cuối của tất cả phần tử.
Phần tử 1
Phần tử 3
Nút đầu
Nút cuối
Code Matlab:
phantu(1,1)=1
phantu(1,2)=2
phantu(2,1)=2;
phantu(2,2)=3;
phantu(3,1)=3;
phantu(3,2)=4;
6.1.5. Tạo ma trận thể hiện mối quan hệ giữa phần tử và
nút tổng thể
Khoa Khoa học ứng dụng
Bộ môn Cơ kỹ thuật
Đại học Quốc gia Thành phố Hồ Chí Minh
Trường đại học Bách Khoa
Lập trình tính toán Cơ kỹ thuật
Matlab® – Matrix Laboratory 10
số phần tử
số nút tổng thể
số nút trên một phần tử
Khởi tạo vector lực tổng thể
Khởi tạo ma trận độ cứng tổng thể
Code Matlab:
sophantu = size(phantu,1);
sonut_pt = size(phantu,2);
sonut = size(toado,1);
Ktt = zeros(sonut*btd_nut,sonut*btd_nut);
Ftt = zeros(sonut*btd_nut,1);
lucphap = zeros(1,sophantu);
luctiep = zeros(1,sophantu);
6.1.6. Xác định số phần tử, số nút, khởi tạo các ma trận
Khởi tạo vector lực pháp tuyến
Khởi tạo vector lực tiếp tuyến
Khoa Khoa học ứng dụng
Bộ môn Cơ kỹ thuật
Đại học Quốc gia Thành phố Hồ Chí Minh
Trường đại học Bách Khoa
Lập trình tính toán Cơ kỹ thuật Matlab® – Matrix Laboratory
Code Matlab:
rangbuoc = [1 2 3 10 11 12];
7
q
8
q
9
q
4
q
5
q
6
q
1
q
2
q
3
q
10
q
11
q
12
q
1 2
3
4
(1)
(2 )
(3)
11
• Tạo ma trận chứa số thứ tự của bậc tự do
bị ràng buộc
• Vector rangbuoc là ma trận có 1 hàng n cột,
n là số lượng bậc tự do bị ràng buộc
• Mỗi phần tử của ma trận chứa số thứ tự
tổng thể của bậc tự do bị ràng buộc.
Phần tử thứ 4 của
ma trận rangbuoc
chứa bậc tự do
thứ 10 bị ràng
buộc
Phần tử thứ 2 của
ma trận rangbuoc
chứa bậc tự do
thứ 2 bị ràng buộc
x x
x
x
x
6.1.7. Nhập điều kiện biên ràng buộc chuyển vị
A B
C
D
x
Khoa Khoa học ứng dụng
Bộ môn Cơ kỹ thuật
Đại học Quốc gia Thành phố Hồ Chí Minh
Trường đại học Bách Khoa
Lập trình tính toán Cơ kỹ thuật Matlab® – Matrix Laboratory 12
Code Matlab:
for i = 1:length(rangbuoc)
m = mod(rangbuoc(i),3);
if m == 0
x = toado(rangbuoc(i)/3,1);
y = toado(rangbuoc(i)/3,2);
plot(x,y,'r+','markersize',20,'linewidth',3);
plot(x,y,'ro','markersize',15,'linewidth',2);
end
m = mod(rangbuoc(i)+1,3);
if m == 0
x = toado((rangbuoc(i)+1)/3,1);
y = toado((rangbuoc(i)+1)/3,2);
plot(x,y-
0.05,'^','markersize',17,'markerfacecolor','g');
end
m = mod(rangbuoc(i)+2,3);
if m == 0
x = toado((rangbuoc(i)+2)/3,1);
y = toado((rangbuoc(i)+2)/3,2);
plot(x-0.05,y,'>','markersize',17,'markerfacecolor','b');
end
end
6.1.8. Biểu diễn đồ họa điều kiện biên ràng buộc chuyển vị
9/3/2019
3
Khoa Khoa học ứng dụng
Bộ môn Cơ kỹ thuật
Đại học Quốc gia Thành phố Hồ Chí Minh
Trường đại học Bách Khoa
Lập trình tính toán Cơ kỹ thuật Matlab® – Matrix Laboratory 13
Code Matlab:
for i = 1 : sophantu
nutdau = phantu(i,1);
nutcuoi = phantu(i,2);
xi = toado(nutdau,1);
xj = toado(nutcuoi,1);
yi = toado(nutdau,2);
yj = toado(nutcuoi,2);
L(i) = sqrt((xj-xi)^2 + (yj-yi)^2);%chieu dai phan tu
C(i) = (xj-xi)/L(i);%cos(alpha) cua phan tu
S(i) = (yj-yi)/L(i);%sin(alpha) cua phan tu
plot([xi xj],[yi yj],'-','linewidth',3);
plot([xi xj],[yi yj],'o','markersize',10);
text(xi+0.03,yi,num2str(nutdau));
text(xj+0.03,yj,num2str(nutcuoi));
if xj>xi
text((xj-xi)/2+xi,(yj-yi)/2 + yi,num2str(i));
elseif xj == xi
text(xj,(yj-yi)/2 + yi,num2str(i));
elseif xj<xi
text((xi-xj)/2+xj,(yj-yi)/2 + yi,num2str(i));
end
axis equal
grid on
end
Khoa Khoa học ứng dụng
Bộ môn Cơ kỹ thuật
Đại học Quốc gia Thành phố Hồ Chí Minh
Trường đại học Bách Khoa
Lập trình tính toán Cơ kỹ thuật Matlab® – Matrix Laboratory 14
6.1.9. Nhập điều kiện biên tải trọng
1
1
1
4
4
4
1
2
3
4
5
6
7
8
9
10
11
12
0
0
0
0
0
10000
nu t
H
V
M
F
H
V
M
A B
C
D
1
H
1
V
1
M
4
H
4
V
4
M
(1)
(2 )
(3)
6.1.9.1. Nhập điều kiện biên tải trọng tại nút
Code Matlab:
Ftt(9) = 10000;
Khoa Khoa học ứng dụng
Bộ môn Cơ kỹ thuật
Đại học Quốc gia Thành phố Hồ Chí Minh
Trường đại học Bách Khoa
Lập trình tính toán Cơ kỹ thuật
Matlab® – Matrix Laboratory 15
0 0
0 0
0
e
0 0
0 0
0
p c q s
p s q c
q
6
P
p c q s
2
p s q c
q
6
qo
x'
y'
y
x
1
2
P2
P1
po
P3
P4
P5
P6
6.1.9.2. Quy đổi tải đặt trên phần tử thành tải tập trung tại nút
6.1.9.2.1. Tải đặt trên phần tử là tải phân bố pháp tuyến/tiếp tuyến
Khoa Khoa học ứng dụng
Bộ môn Cơ kỹ thuật
Đại học Quốc gia Thành phố Hồ Chí Minh
Trường đại học Bách Khoa
Lập trình tính toán Cơ kỹ thuật
Matlab® – Matrix Laboratory
Q
y
x
1
2
P2
P1
P3
P4
P5
P6
/2
/2
x'
y’
2
2
8
2
2
8
e
Q
s
Q
c
QL
PQ
s
Q
c
QL
6.1.9.2.2. Tải đặt trên phần tử là lực tập trung tại giữa phần tử
Khoa Khoa học ứng dụng
Bộ môn Cơ kỹ thuật
Đại học Quốc gia Thành phố Hồ Chí Minh
Trường đại học Bách Khoa
Lập trình tính toán Cơ kỹ thuật Matlab® – Matrix Laboratory
M
x’
y’
y
x
1
2
P2
P1
P3
P4
P5
P6
/2
/2
32
32
4
32
32
4
e
M
s
M
c
M
PM
s
M
c
M
6.1.9.2.3. Tải đặt trên phần tử là moment tập trung tại giữa phần tử
Khoa Khoa học ứng dụng
Bộ môn Cơ kỹ thuật
Đại học Quốc gia Thành phố Hồ Chí Minh
Trường đại học Bách Khoa
Lập trình tính toán Cơ kỹ thuật Matlab® – Matrix Laboratory 18
A B
C
D
1
H
1
V
1
M
4
H
4
V
4
M
(1)
(2 )
(3)
Trong mô hình đang xét, vì số
phần tử được sử dụng là 3 phần
tử nên
• Tải phân bố pháp tuyến sẽ
tác dụng lên phần tử (1)
• Tải phân bố tiếp tuyến tác
dụng lên phần tử (2) và (3)
9/3/2019
4
Khoa Khoa học ứng dụng
Bộ môn Cơ kỹ thuật
Đại học Quốc gia Thành phố Hồ Chí Minh
Trường đại học Bách Khoa
Lập trình tính toán Cơ kỹ thuật Matlab® – Matrix Laboratory 19
(1 )
0
5 0 0 0
8 3 3, 3
0
5 0 0 0
8 3 3
1
2
3
4
5
6
, 3
P
0
0
2
0
e
0
0
2
0
q
s
2
q
c
2
q
12
P' q
s
2
q
c
2
q
12
A B
C
D
1
H
1
V
1
M
4
H
4
V
4
M
(1)
(2 )
(3)
Quy đổi tải phân bố pháp tuyến thành tải tập
trung hai đầu nút cho phần tử (1)
Khoa Khoa học ứng dụng
Bộ môn Cơ kỹ thuật
Đại học Quốc gia Thành phố Hồ Chí Minh
Trường đại học Bách Khoa
Lập trình tính toán Cơ kỹ thuật Matlab® – Matrix Laboratory 20
Code Matlab:
function lucphap_tt = quydoi_lucphap(lucphap,S,C,L)
lucphap_tt(1,1) = (L/2)* (-lucphap*S);
lucphap_tt(2,1) = (L/2)* (lucphap*C);
lucphap_tt(3,1) = (L/2)* (lucphap*L/6);
lucphap_tt(4,1) = (L/2)* (-lucphap*S);
lucphap_tt(5,1) = (L/2)* (lucphap*C);
lucphap_tt(6,1) = (L/2)* (-lucphap*L/6);
Khoa Khoa học ứng dụng
Bộ môn Cơ kỹ thuật
Đại học Quốc gia Thành phố Hồ Chí Minh
Trường đại học Bách Khoa
Lập trình tính toán Cơ kỹ thuật
Matlab® – Matrix Laboratory 21
0
0
e
0
0
p c
p s
0
P '
p c
2
p s
0
(2 ) (3)
4 7
5 8
0
2500
0
0
2500
0
6 9
7 10
8 11
9 12
P P
Quy đổi tải phân bố tiếp tuyến thành tải tập
trung hai đầu nút cho phần tử (2) và (3)
A B
C
D
1
H
1
V
1
M
4
H
4
V
4
M
(1)
(2 )
(3)
Khoa Khoa học ứng dụng
Bộ môn Cơ kỹ thuật
Đại học Quốc gia Thành phố Hồ Chí Minh
Trường đại học Bách Khoa
Lập trình tính toán Cơ kỹ thuật
Matlab® – Matrix Laboratory 22
Code Matlab:
function luctiep_tt = quydoi_luctiep(luctiep,S,C,L)
luctiep_tt(1,1) = (L/2)*(luctiep*C);
luctiep_tt(2,1) = (L/2)*(luctiep*S);
luctiep_tt(3,1) = (L/2)*0;
luctiep_tt(4,1) = (L/2)*(luctiep*C);
luctiep_tt(5,1) = (L/2)*(luctiep*S);
luctiep_tt(6,1) = (L/2)*0;
Khoa Khoa học ứng dụng
Bộ môn Cơ kỹ thuật
Đại học Quốc gia Thành phố Hồ Chí Minh
Trường đại học Bách Khoa
Lập trình tính toán Cơ kỹ thuật Matlab® – Matrix Laboratory 23
6.2. Giải (Solution)
Khoa Khoa học ứng dụng
Bộ môn Cơ kỹ thuật
Đại học Quốc gia Thành phố Hồ Chí Minh
Trường đại học Bách Khoa
Lập trình tính toán Cơ kỹ thuật Matlab® – Matrix Laboratory 24
Code Matlab:
for i=1:sophantu
nutdau = phantu(i,1);
nutcuoi = phantu(i,2);
mt_bool_pt = tinh_bool([nutdau nutcuoi],btd_nut);
Kpt = tinh_Kpt(E,A,J,L(i),S(i),C(i));
Ktt = tinh_Ktt(Ktt,Kpt,mt_bool_pt);
lucphap_tt = quydoi_lucphap(lucphap(i),S(i),C(i),L(i));
luctiep_tt = quydoi_luctiep(luctiep(i),S(i),C(i),L(i));
Ftt_td = lucphap_tt + luctiep_tt;
Ftt = chuyen_Ftt_td(Ftt,Ftt_td,mt_bool_pt);
end
Giải thích:
Bắt đầu vòng lặp
• Xác định tọa độ nút đầu và nút cuối phần tử thứ i
• Tính ma trận bool cho phần tử thứ i qua hàm con tinh_bool
• Tính ma trận cứng Kpt của phần tử thứ i qua hàm con tinh_Kpt
• Từ ma trận bool, lắp ghép Kpt phần tử thứ i vào Ktt tổng thể qua hàm con
tinh_Ktt
• Quy đổi tải phân bố pháp tuyến trên phần tử về tải tập trung hai đầu nút.
• Quy đổi tải phân bố tiếp tuyến trên phần tử về tải tập trung hai đầu nút.
• Tính vector lực tổng thể Ftt qua hàm con chuyen_Ftt_td
Quay lại vòng lặp với phần tử thứ i+1
6.2.1. Xây dựng ma trận độ cứng tổng thể
9/3/2019
5
Khoa Khoa học ứng dụng
Bộ môn Cơ kỹ thuật
Đại học Quốc gia Thành phố Hồ Chí Minh
Trường đại học Bách Khoa
Lập trình tính toán Cơ kỹ thuật Matlab® – Matrix Laboratory
Code Matlab:
function [mt_bool_pt] = tinh_bool(mt_nut_pt,btd_nut)
k = 0;
sonut_pt = size(mt_nut_pt,2);
for i = 1:sonut_pt
s = (mt_nut_pt(i) - 1)*btd_nut;
for j = 1:btd_nut
k=k+1;
mt_bool_pt(k) = s + j;
end
end
25
ma trận chứa số thứ tự tổng
thể các nút của phần tử
Bậc tự do của mỗi nút
Số nút của phần tử
ma trận bool của một phần
tử thứ i
Hàm tinh_bool có nhiệm vụ tính ma trận bool của một phần tử.
6.2.1.1. Tính ma trận bool cho một phần tử
Nút i Nút j
pt (1)
pt (2)
pt (3)
( )
1
e
q
( )
3
e
q
( )
4
e
q
( )
2
e
q
3
q
4
q
1
q
2
q
5
q
6
q
7
q
8
q
7
q
1 2
q
1 1
q
8
q
( )
5
e
q
( )
6
e
q
4
q
5
q
6
q
9
q
9
q
1 0
q
Khoa Khoa học ứng dụng
Bộ môn Cơ kỹ thuật
Đại học Quốc gia Thành phố Hồ Chí Minh
Trường đại học Bách Khoa
Lập trình tính toán Cơ kỹ thuật Matlab® – Matrix Laboratory 26
6.2.1.2. Tính ma trận độ cứng K cho một phần tử
Khoa Khoa học ứng dụng
Bộ môn Cơ kỹ thuật
Đại học Quốc gia Thành phố Hồ Chí Minh
Trường đại học Bách Khoa
Lập trình tính toán Cơ kỹ thuật
Matlab® – Matrix Laboratory 27
Code Matlab:
function [Kpt] = tinh_Kpt(E,A,J,L,S,C)
B=12*J/(L^2);
Kpt = (E/L)*[A*C^2+B*S^2 (A-B)*C*S -B*L*S/2 -(A*C^2+B*S^2) -(A-B)*C*S -B*L*S/2;
(A-B)*C*S A*S^2+B*C^2 B*L*C/2 -(A-B)*C*S -(A*S^2+B*C^2) B*L*C/2;
-B*L*S/2 B*L*C/2 4*J B*L*S/2 -B*L*C/2 2*J;
-(A*C^2+B*S^2) -(A-B)*C*S B*L*S/2 A*C^2+B*S^2 (A-B)*C*S B*L*S/2;
-(A-B)*C*S -(A*S^2+B*C^2) -B*L*C/2 (A-B)*C*S A*S^2+B*C^2 -
B*L*C/2;
-B*L*S/2 B*L*C/2 2*J B*L*S/2 -B*L*C/2 4*J];
Hàm tinh_Kpt có nhiệm vụ tính ma trận cứng K của một phần tử.
Khoa Khoa học ứng dụng
Bộ môn Cơ kỹ thuật
Đại học Quốc gia Thành phố Hồ Chí Minh
Trường đại học Bách Khoa
Lập trình tính toán Cơ kỹ thuật
Matlab® – Matrix Laboratory
8
(1)
10 0 0 10 0 0
0 0,12 0,06 0 0,12 0,06
0 0,06 0,04 0 0,06 0,02
2.10 10 0 0 10 0 0
0 0,12 0,06 0 0,12 0,06
0 0,06 0,02 0 0,06 0,4
1
2
3
, /
4
5
6
K
N m
28
Ma trận bool của phần tử 1
1 2 3 4 5 6
Khoa Khoa học ứng dụng
Bộ môn Cơ kỹ thuật
Đại học Quốc gia Thành phố Hồ Chí Minh
Trường đại học Bách Khoa
Lập trình tính toán Cơ kỹ thuật Matlab® – Matrix Laboratory 29
Ma trận bool của phần tử 2
8
(2)
0,96 0 0,24 0,96 0 0,24
0 20 0 0 20 0
0,24 0 0,08 0,24 0 0,04
2.10 0,96 0 0,24 0,96 0 0,24
0 20 0 0 20 0
0,24 0 0,04 0,24 0 0,0
4
5
6
, /
7
8
9
8
NK
m
4 5 6 7 8 9
Khoa Khoa học ứng dụng
Bộ môn Cơ kỹ thuật
Đại học Quốc gia Thành phố Hồ Chí Minh
Trường đại học Bách Khoa
Lập trình tính toán Cơ kỹ thuật Matlab® – Matrix Laboratory 30
Ma trận bool của phần tử 3
8
(3)
0,96 0 0,24 0,96 0 0,24
0 20 0 0 20 0
0,24 0 0,08 0,24 0 0,04
2.10 0,96 0 0,24 0,96 0 0,24
0 20 0 0 20 0
0,24 0 0,04 0,2
7
8
9
, /
10
11
12
4 0 0,08
NK
m
7 8 9 1 0 1 1 1 2