Luận án Đánh giá khả năng chịu lửa của các cấu kiện thép chịu lực được bọc bảo vệ ứng dụng cho các công trình nhà tại Việt Nam

BỘ GIÁO DỤC VÀ ĐÀO TẠO 
TRƯỜNG ĐẠI HỌC XÂY DỰNG HÀ NỘI 
-------------&------------- 
Phạm Thị Ngọc Thu 
ĐÁNH GIÁ KHẢ NĂNG CHỊU LỬA CỦA 
CÁC CẤU KIỆN THÉP CHỊU LỰC ĐƯỢC BỌC BẢO VỆ 
ỨNG DỤNG CHO CÁC CÔNG TRÌNH NHÀ TẠI VIỆT NAM 
LUẬN ÁN TIẾN SỸ 
Hà Nội – Năm 2021 
BỘ GIÁO DỤC VÀ ĐÀO TẠO 
TRƯỜNG ĐẠI HỌC XÂY DỰNG HÀ NỘI 
-------------&------------- 
Phạm Thị Ngọc Thu 
ĐÁNH GIÁ KHẢ NĂNG CHỊU LỬA CỦA 
CÁC CẤU KIỆN THÉP CHỊU LỰC ĐƯỢC BỌC BẢO VỆ 
ỨNG DỤNG CHO CÁC CÔNG TRÌNH NHÀ TẠI VIỆT NAM 
LUẬN ÁN TIẾN SỸ 
Chuyên ngành: Kỹ thuật xây dựng 
Mã số: 9580201 
 CÁN BỘ HƯỚNG DẪN 
 GS.TS. Phạm Văn Hội 
Hà Nội – Năm 2021 
 i 
LỜI CAM ĐOAN 
Tôi xin cam đoan những nội dung trong luận án là kết quả công trình nghiên cứu 
của riêng tôi. Các số liệu và kết quả trong luận án là trung thực và chưa từng được 
công bố trong bất kỳ công trình nào khác. 
 ii 
MỤC LỤC 
 Trang 
Lời cam đoan  i 
Mục lục ..... ii 
Danh mục các chữ viết tắt, ký hiệu, thuật ngữ . vi 
Danh mục các bảng biểu .. xi 
Danh mục các hình vẽ, sơ đồ, đồ thị  xiii 
Mở đầu ..... 01 
1 Lý do lựa chọn đề tài . 01 
2 Mục đích nghiên cứu . 04 
3 Đối tượng và phạm vi nghiên cứu . 04 
4 Phương pháp nghiên cứu ... 04 
5 Cơ sở khoa học .. 05 
6 Đóng góp của luận án  05 
7 Cấu trúc và nội dung luận án . 06 
Chương 1. Tổng quan về các hình thức bọc bảo vệ KCT và các phương pháp 
tính toán kết cấu trong điều kiện chịu lửa ứng dụng cho các công trình nhà ... 
08 
1.1 Tổng quan về các giải pháp bọc bảo vệ kết cấu thép ứng dụng cho các 
công trình nhà .................... 
08 
1.1.1 Các giải pháp thiết kế phòng chống cháy cho các công trình nhà . 08 
1.1.2 Vật liệu thạch cao chống cháy 09 
1.1.2.1 Quy trình sản xuất thạch cao chống cháy .. 09 
1.1.2.2 Cơ chế làm việc của vật liệu thạch cao chống cháy .. 10 
1.1.2.3 Các hình thức cấu tạo bọc cấu kiện thép bằng tấm thạch cao chống 
cháy ... 
11 
1.1.3 Vật liệu vữa chống cháy  14 
1.2 Tổng quan về quy trình tính toán kết cấu thép trong điều kiện chịu lửa 16 
1.2.1 Bài toán tính toán kết cấu thép chịu lực khi chịu lửa  16 
 iii
1.2.2 Sự biến thiên nhiệt độ theo thời gian trong mô hình cháy .................... 16 
1.2.2.1 Mô hình cháy danh nghĩa .. 17 
1.2.2.2 Mô hình cháy tham biến  18 
1.2.2.3 Mô hình cháy thực tế . 19 
1.2.3 Sự biến thiên nhiệt độ bên trong kết cấu thép theo mô hình kết cấu ..... 20 
1.2.4 Các nguyên tắc tính toán cơ bản ............................................................ 21 
1.2.5 Tải trọng tác dụng lên kết cấu trong điều kiện chịu lửa ........................ 22 
1.2.6 Trạng thái làm việc của các cấu kiện dầm, cột thép trong điều kiện 
chịu lửa .................................................................................................. 
24 
1.3 Các nghiên cứu về tính toán kết cấu thép trong điều kiện chịu lửa ....... 26 
1.3.1 Các kết quả nghiên cứu trên thế giới ..................................................... 26 
1.3.2 Các kết quả nghiên cứu ở Việt Nam ...................................................... 29 
1.4 Các tiêu chuẩn, quy chuẩn áp dụng ....................................................... 31 
 Nhận xét chương 1 ................................................................................ 34 
Chương 2. Bài toán truyền nhiệt trong các cấu kiện thép không bọc và được 
bọc bảo vệ ......... 
36 
2.1 Bản chất của quá trình truyền nhiệt ....................................................... 36 
2.1.1 Hiện tượng bức xạ ..... 37 
2.1.2 Hiện tượng đối lưu ................................................................................ 38 
2.2 Sự thay đổi các đặc tính nhiệt của thép trong điều kiện chịu lửa .......... 39 
2.2.1 Hệ số giãn nở vì nhiệt ............................................................................ 39 
2.2.2 Nhiệt dung riêng .................................................................................... 41 
2.2.3 Hệ số dẫn nhiệt ...................................................................................... 41 
2.3 Bài toán truyền nhiệt trong không gian ba chiều .................................. 42 
2.3.1 Phương trình truyền nhiệt trong không gian ba chiều ........................... 42 
2.3.2 Thuật toán DT3D ứng dụng phương pháp PTHH trong bài toán 
truyền nhiệt ba chiều ............................................................................ 
44 
2.3.2.1 Mô hình phần tử hữu hạn ...................................................................... 44 
2.3.2.2 Xây dựng hàm dạng .............................................................................. 45 
 iv 
2.3.2.3 Xây dựng ma trận truyền nhiệt của hệ ................................................... 47 
2.3.2.4 Rời rạc hóa các bước thời gian trong quá trình khảo sát ....................... 49 
2.4 Ví dụ tính toán và kiểm chứng .............................................................. 52 
2.4.1 Ví dụ tính toán ....................................................................................... 52 
2.4.2 Ví dụ kiểm chứng thí nghiệm ................................................................ 55 
 Nhận xét chương 2 ................................................................................ 57 
Chương 3. Xác định khả nănng chịu lực của cấu kiện thép chịu lực trong điều 
kiện chịu lửa theo phương pháp đơn giản hóa ................................................. 
59 
3.1 Sự thay đổi các đặc tính cơ lý của thép trong điều kiện chịu lửa .......... 59 
3.2 Quy trình tính toán cấu kiện thép chịu lực trong điều kiện chịu lửa 
theo phương pháp tính đơn giản hóa SMD ........................................... 
63 
3.2.1 Khả năng chịu lực của các cấu kiện cơ bản ........................................... 63 
3.2.1.1 Các cấu kiện dầm chịu uốn theo điều kiện bền ..................................... 63 
3.2.1.2 Các cấu kiện cột chịu nén theo điều kiện ổn định tổng thể ................... 65 
3.2.2 Quy trình tính toán cấu kiện thép chịu lực trong điều kiện chịu lửa 
theo phương pháp tính đơn giản hóa SMD ........................................... 
69 
3.3 Ứng dụng SMD cho cấu kiện dầm thép được bọc bảo vệ chịu lực 
trong điều kiện chịu lửa ........................................................................ 
72 
3.3.1 Trạng thái làm việc của cấu kiện dầm thép chịu lực trong điều kiện 
chịu lửa .................................................................................................. 
72 
3.3.2 Dầm thép bọc bảo vệ ba mặt theo chu vi .............................................. 73 
3.3.3 Dầm thép bọc bảo vệ ba mặt dạng hình hộp ......................................... 76 
3.4 Ứng dụng SMD cho cấu kiện cột thép được bọc bảo vệ chịu lực trong 
điều kiện chịu lửa .................................................................................. 
78 
3.4.1 Trạng thái làm việc của cấu kiện cột thép chịu lực trong điều kiện 
chịu lửa .................................................................................................. 
78 
3.4.2 Cột thép bọc bảo vệ bốn mặt theo chu vi .. 79 
3.4.3 Cột thép bọc bảo vệ bốn mặt dạng hình hộp . 81 
 Nhận xét chương 3 ................................................................................ 84 
 v
Chương 4. Đánh giá khả năng chịu lực của các cấu kiện thép trong điều kiện 
chịu lửa theo phương pháp mô phỏng số ..................................... 
85 
4.1 Phân tích sự làm việc của cấu kiện thép chịu lực trong điều kiện chịu 
lửa theo phương pháp mô phỏng số ...................................................... 
85 
4.1.1 Mô hình Thermal Analysis .................................................................... 85 
4.1.1.1 Các phương trình truyền nhiệt cơ bản ................................................... 85 
4.1.1.2 Xây dựng ma trận truyền nhiệt của phần tử .......................................... 87 
4.1.1.3 Mô hình phần tử vật liệu thép trong Thermal Analysis ......................... 88 
4.1.2 Mô hình Structural Analysis .................................................................. 89 
4.1.2.1 Mối quan hệ ứng suất-biến dạng ........................................................... 89 
4.1.2.2 Xây dựng ma trận ứng suất-biến dạng của phần tử ... 92 
4.1.2.3 Mô hình phần tử vật liệu thép trong Structural Analysis .. 93 
4.1.2.4 Phân tích Eigenvalue Buckling . 95 
4.1.3 Các thí nghiệm kiểm chứng ................................................................... 95 
4.1.3.1 Thí nghiệm dầm thép không bọc bảo vệ ... 95 
4.1.3.2 Thí nghiệm cột thép chèn gạch .. 98 
4.2 Các ví dụ khảo sát ................................................................................. 101 
4.2.1 Dầm thép bọc bảo vệ ba mặt theo chu vi .............................................. 101 
4.2.2 Dầm thép bọc bảo vệ ba mặt dạng hình hộp ......................................... 109 
4.2.3 Cột thép bọc bảo vệ bốn mặt theo chu vi .. 116 
4.2.4 Cột thép bọc bảo vệ bốn mặt dạng hình hộp . 125 
 Nhận xét chương 4 ................................................................................ 132 
Kết luận chung .............. 134 
Các công trình khoa học của nghiên cứu sinh liên quan đến luận án ................... 137 
Tài liệu tham khảo ................................................................................................ 138 
Phụ lục Phần code chương trình ví dụ tính toán mục 2.4.1................................... PL1 
 vi 
DANH MỤC CÁC CHỮ VIẾT TẮT, KÝ HIỆU, THUẬT NGỮ 
Các chữ cái Latin viết hoa 
A diện tích của cấu kiện 
Ad tác động gián tiếp do lửa gây ra 
Ai diện tích phần thứ i có nhiệt độ qi 
Am chu vi của cấu kiện bị đốt nóng 
[B] ma trận biến dạng - chuyển vị 
C nhiệt dung riêng 
Ca nhiệt dung riêng của vật liệu thép 
[D] ma trận độ cứng đàn hồi 
E tiêu chuẩn về tính kín, modun đàn hồi 
Ea modun đàn hồi của vật liệu thép ở nhiệt độ thường 
Ea,q modun đàn hồi của vật liệu thép ở nhiệt độ q 
Ed nội lực trong điều kiện nhiệt độ thường 
Efi,d nội lực trong điều kiện chịu lửa 
Efi,d,t nội lực thời điểm t trong điều kiện chịu lửa 
G modun đàn hồi trượt 
Gk tải trọng thường xuyên trong điều kiện chịu lửa 
H chiều cao cột 
I tiêu chuẩn về cách nhiệt 
J ma trận Jacobien 
[K] ma trận truyền nhiệt của cả hệ 
L nhịp dầm 
[M] mômen uốn đàn hồi cho phép khi tính toán cấu kiện trong điều kiện 
 nhiệt độ thường 
Mb,fi,t,Rd mômen oằn vặn bên cho phép của tiết diện tại thời điểm t 
Mfi,t,Rd mômen uốn cho phép của tiết diện tại thời điểm t 
Mfi,q,Rd mô men uốn của tiết diện chịu sự phân bố nhiệt độ đều q 
 vii 
N hàm hình dạng 
Nfi,t,Rd lực nén (kéo) cho phép trên tiết diện tại thời điểm t 
Pr hằng số Prant 
Q lượng nhiệt phát sinh trong phân tố 
Qk hoạt tải trong điều kiện chịu lửa 
Qk,1 hoạt tải chính trong điều kiện chịu lửa 
Qk,i hoạt tải phụ thứ i trong điều kiện chịu lửa 
R tiêu chuẩn về khả năng chịu lực 
Re hằng số Raynolds 
Rfi,d,t khả năng chịu lực của kết cấu tại thời điểm t trong điều kiện chịu lửa 
Rd cườ ... K1s;K46=K1s;K47=K1s; 
g11=diff(n1,x)*N(3,1)+diff(n2,x)*N(4,1)+diff(n3,x)*N(10,1)+... 
diff(n4,x)*N(9,1)+diff(n5,x)*N(33,1)+diff(n6,x)*N(34,1)+... 
diff(n7,x)*N(40,1)+diff(n8,x)*N(39,1); 
g22=diff(n1,y)*N(3,2)+diff(n2,y)*N(4,2)+diff(n3,y)*N(10,2)+... 
diff(n4,y)*N(9,2)+diff(n5,y)*N(33,2)+diff(n6,y)*N(34,2)+... 
diff(n7,y)*N(40,2)+diff(n8,y)*N(39,2); 
PL4
g33=diff(n1,z)*N(3,3)+diff(n2,z)*N(4,3)+diff(n3,z)*N(10,3)+... 
diff(n4,z)*N(9,3)+diff(n5,z)*N(33,3)+diff(n6,z)*N(34,3)+... 
diff(n7,z)*N(40,3)+diff(n8,z)*N(39,3); 
g12=diff(n1,y)*N(3,1)+diff(n2,y)*N(4,1)+diff(n3,y)*N(10,1)+... 
diff(n4,y)*N(9,1)+diff(n5,y)*N(33,1)+diff(n6,y)*N(34,1)+... 
diff(n7,y)*N(40,1)+diff(n8,y)*N(39,1); 
g13=diff(n1,z)*N(3,1)+diff(n2,z)*N(4,1)+diff(n3,z)*N(10,1)+... 
diff(n4,z)*N(9,1)+diff(n5,z)*N(33,1)+diff(n6,z)*N(34,1)+... 
diff(n7,z)*N(40,1)+diff(n8,z)*N(39,1); 
g23=diff(n1,z)*N(3,2)+diff(n2,z)*N(4,2)+diff(n3,z)*N(10,2)+... 
diff(n4,z)*N(9,2)+diff(n5,z)*N(33,2)+diff(n6,z)*N(34,2)+... 
diff(n7,z)*N(40,2)+diff(n8,z)*N(39,2); 
g21=g12;g31=g13;g32=g23; 
J3=[g11 g12 g13;g21 g22 g23;g31 g32 g33]; 
J3nd=inv(J3); 
for e=1:snp 
 for k=1:snp 
f3(e,k)=1/det(J3)*J3nd(1,1)*diff(n(e),x)*1/det(J3)*J3nd(1,1)*diff(n(k),x)
+... 
 1/det(J3)*J3nd(2,2)*diff(n(e),y)*1/det(J3)*J3nd(2,2)*diff(n(k),y)+... 
 1/det(J3)*J3nd(3,3)*diff(n(e),z)*1/det(J3)*J3nd(3,3)*diff(n(k),z); 
A3(e,k)=det(J3)*int(int(int(f3(e,k),x,-1,1),y,-1,1),z,-1,1); 
 end; 
end; 
K3s=Lamdathep*single(A3); 
K17=K3s;K31=K3s;K45=K3s; 
l11=diff(n1,x)*N(9,1)+diff(n2,x)*N(10,1)+diff(n3,x)*N(14,1)+... 
diff(n4,x)*N(13,1)+diff(n5,x)*N(39,1)+diff(n6,x)*N(40,1)+... 
diff(n7,x)*N(44,1)+diff(n8,x)*N(43,1); 
l22=diff(n1,y)*N(9,2)+diff(n2,y)*N(10,2)+diff(n3,y)*N(14,2)+... 
diff(n4,y)*N(13,2)+diff(n5,y)*N(39,2)+diff(n6,y)*N(40,2)+... 
diff(n7,y)*N(44,2)+diff(n8,y)*N(43,2); 
l33=diff(n1,z)*N(9,3)+diff(n2,z)*N(10,3)+diff(n3,z)*N(14,3)+... 
diff(n4,z)*N(13,3)+diff(n5,z)*N(39,3)+diff(n6,z)*N(40,3)+... 
diff(n7,z)*N(44,3)+diff(n8,z)*N(43,3); 
l12=diff(n1,y)*N(9,1)+diff(n2,y)*N(10,1)+diff(n3,y)*N(14,1)+... 
diff(n4,y)*N(13,1)+diff(n5,y)*N(39,1)+diff(n6,y)*N(40,1)+... 
diff(n7,y)*N(44,1)+diff(n8,y)*N(43,1); 
l21=l12; 
l13=diff(n1,z)*N(9,1)+diff(n2,z)*N(10,1)+diff(n3,z)*N(14,1)+... 
diff(n4,z)*N(13,1)+diff(n5,z)*N(39,1)+diff(n6,z)*N(40,1)+... 
diff(n7,z)*N(44,1)+diff(n8,z)*N(43,1); 
l31=l13; 
l23=diff(n1,z)*N(9,2)+diff(n2,z)*N(10,2)+diff(n3,z)*N(14,2)+... 
diff(n4,z)*N(13,2)+diff(n5,z)*N(39,2)+diff(n6,z)*N(40,2)+... 
diff(n7,z)*N(44,2)+diff(n8,z)*N(43,2); 
l32=l23; 
J6=[l11 l12 l13;l21 l22 l23;l31 l32 l33]; 
J6nd=inv(J6); 
for e=1:snp 
 for k=1:snp 
f6(e,k)=1/det(J6)*J6nd(1,1)*diff(n(e),x)*1/det(J6)*J6nd(1,1)*diff(n(k),x)
+... 
 1/det(J6)*J6nd(2,2)*diff(n(e),y)*1/det(J6)*J6nd(2,2)*diff(n(k),y)+... 
 1/det(J6)*J6nd(3,3)*diff(n(e),z)*1/det(J6)*J6nd(3,3)*diff(n(k),z); 
A6(e,k)=det(J6)*int(int(int(f6(e,k),x,-1,1),y,-1,1),z,-1,1); 
PL5
 end; 
end; 
K6s=Lamdathep*single(A6); 
K7=K6s;K8=K6s;K9=K6s;K20=K6s;K21=K6s;K22=K6s;K23=K6s; 
K34=K6s;K35=K6s;K36=K6s;K37=K6s;K48=K6s;K49=K6s;K50=K6s;K51=K6s; 
h11=diff(n1,x)*N(19,1)+diff(n2,x)*N(20,1)+diff(n3,x)*N(26,1)+... 
diff(n4,x)*N(25,1)+diff(n5,x)*N(49,1)+diff(n6,x)*N(50,1)+... 
diff(n7,x)*N(56,1)+diff(n8,x)*N(55,1); 
h22=diff(n1,y)*N(19,2)+diff(n2,y)*N(20,2)+diff(n3,y)*N(26,2)+... 
diff(n4,y)*N(25,2)+diff(n5,y)*N(49,2)+diff(n6,y)*N(50,2)+... 
diff(n7,y)*N(56,2)+diff(n8,y)*N(55,2); 
h33=diff(n1,z)*N(19,3)+diff(n2,z)*N(20,3)+diff(n3,z)*N(26,3)+... 
diff(n4,z)*N(25,3)+diff(n5,z)*N(49,3)+diff(n6,z)*N(50,3)+... 
diff(n7,z)*N(56,3)+diff(n8,z)*N(55,3); 
h12=diff(n1,y)*N(19,1)+diff(n2,y)*N(20,1)+diff(n3,y)*N(26,1)+... 
diff(n4,y)*N(25,1)+diff(n5,y)*N(49,1)+diff(n6,y)*N(50,1)+... 
diff(n7,y)*N(56,1)+diff(n8,y)*N(55,1); 
h13=diff(n1,z)*N(19,1)+diff(n2,z)*N(20,1)+diff(n3,z)*N(26,1)+... 
diff(n4,z)*N(25,1)+diff(n5,z)*N(49,1)+diff(n6,z)*N(50,1)+... 
diff(n7,z)*N(56,1)+diff(n8,z)*N(55,1); 
h23=diff(n1,z)*N(19,2)+diff(n2,z)*N(20,2)+diff(n3,z)*N(26,2)+... 
diff(n4,z)*N(25,2)+diff(n5,z)*N(49,2)+diff(n6,z)*N(50,2)+... 
diff(n7,z)*N(56,2)+diff(n8,z)*N(55,2); 
h21=h12;h31=h13;h32=h23; 
J10=[h11 h12 h13;h21 h22 h23;h31 h32 h33]; 
J10nd=inv(J10); 
for e=1:snp 
 for k=1:snp 
f10(e,k)=1/det(J10)*J10nd(1,1)*diff(n(e),x)*1/det(J10)*J10nd(1,1)*diff(n(
k),x)+... 
1/det(J10)*J10nd(2,2)*diff(n(e),y)*1/det(J10)*J10nd(2,2)*diff(n(k),y)+... 
1/det(J10)*J10nd(3,3)*diff(n(e),z)*1/det(J10)*J10nd(3,3)*diff(n(k),z); 
A10(e,k)=det(J10)*int(int(int(f10(e,k),x,-1,1),y,-1,1),z,-1,1); 
 end; 
end; 
A10s=Lamdathep*single(A10); 
m11=diff(n1,x)*N(21,1)+diff(n2,x)*N(22,1)+diff(n3,x)*N(28,1)+... 
diff(n4,x)*N(27,1)+diff(n5,x)*N(51,1)+diff(n6,x)*N(52,1)+... 
diff(n7,x)*N(58,1)+diff(n8,x)*N(57,1); 
m22=diff(n1,y)*N(21,2)+diff(n2,y)*N(22,2)+diff(n3,y)*N(28,2)+... 
diff(n4,y)*N(27,2)+diff(n5,y)*N(51,2)+diff(n6,y)*N(52,2)+... 
diff(n7,y)*N(58,2)+diff(n8,y)*N(57,2); 
m33=diff(n1,z)*N(21,3)+diff(n2,z)*N(22,3)+diff(n3,z)*N(28,3)+... 
diff(n4,z)*N(27,3)+diff(n5,z)*N(51,3)+diff(n6,z)*N(52,3)+... 
diff(n7,z)*N(58,3)+diff(n8,z)*N(57,3); 
m12=diff(n1,y)*N(21,1)+diff(n2,y)*N(22,1)+diff(n3,y)*N(28,1)+... 
diff(n4,y)*N(27,1)+diff(n5,y)*N(51,1)+diff(n6,y)*N(52,1)+... 
diff(n7,y)*N(58,1)+diff(n8,y)*N(57,1); 
m21=m12; 
m13=diff(n1,z)*N(21,1)+diff(n2,z)*N(22,1)+diff(n3,z)*N(28,1)+... 
diff(n4,z)*N(27,1)+diff(n5,z)*N(51,1)+diff(n6,z)*N(52,1)+... 
diff(n7,z)*N(58,1)+diff(n8,z)*N(57,1); 
m31=m13; 
PL6
m23=diff(n1,z)*N(21,2)+diff(n2,z)*N(22,2)+diff(n3,z)*N(28,2)+... 
diff(n4,z)*N(27,2)+diff(n5,z)*N(51,2)+diff(n6,z)*N(52,2)+... 
diff(n7,z)*N(58,2)+diff(n8,z)*N(57,2); 
m32=m23; 
J12=[m11 m12 m13;m21 m22 m23;m31 m32 m33]; 
J12nd=inv(J12); 
for e=1:snp 
 for k=1:snp 
f12(e,k)=1/det(J12)*J12nd(1,1)*diff(n(e),x)*1/det(J12)*J12nd(1,1)*diff(n(
k),x)+... 
1/det(J12)*J12nd(2,2)*diff(n(e),y)*1/det(J12)*J12nd(2,2)*diff(n(k),y)+... 
1/det(J12)*J12nd(3,3)*diff(n(e),z)*1/det(J12)*J12nd(3,3)*diff(n(k),z); 
A12(e,k)=det(J12)*int(int(int(f12(e,k),x,-1,1),y,-1,1),z,-1,1); 
 end; 
end; 
A12s=Lamdathep*single(A12); 
for e=1:snp 
 for k=1:snp 
c(e,k)=n(e)*n(k); 
C1(e,k)=ra*ca*det(J1)*int(int(int(c(e,k),x,-1,1),y,-1,1),z,-1,1); 
 end; 
end; 
C1s=single(C1);C2=C1s;C4=C1s;C5=C1s;C15=C1s;C16=C1s;C18=C1s;C19=C1s; 
C29=C1s;C30=C1s;C32=C1s;C33=C1s; 
C43=C1s;C44=C1s;C46=C1s;C47=C1s; 
for e=1:snp 
 for k=1:snp 
c(e,k)=n(e)*n(k); 
C3(e,k)=ra*ca*det(J3)*int(int(int(c(e,k),x,-1,1),y,-1,1),z,-1,1); 
 end; 
end; 
C3s=single(C3);C17=C3s;C31=C3s;C45=C3s; 
for e=1:snp 
 for k=1:snp 
c(e,k)=n(e)*n(k); 
C6(e,k)=ra*ca*det(J6)*int(int(int(c(e,k),x,-1,1),y,-1,1),z,-1,1); 
 end; 
end; 
C6s=single(C6); 
C7=C6s;C8=C6s;C9=C6s;C20=C6s;C21=C6s;C22=C6s;C23=C6s; 
C34=C6s;C35=C6s;C36=C6s;C37=C6s;C48=C6s;C49=C6s;C50=C6s;C51=C6s; 
for e=1:snp 
 for k=1:snp 
c(e,k)=n(e)*n(k); 
C10(e,k)=ra*ca*det(J10)*int(int(int(c(e,k),x,-1,1),y,-1,1),z,-1,1); 
 end; 
end; 
C10s=single(C10); 
C11=C10s;C13=C10s;C14=C10s;C24=C10s;C25=C10s;C27=C10s;C28=C10s; 
C38=C10s;C39=C10s;C41=C10s;C42=C10s;C52=C10s;C53=C10s;C55=C10s;C56=C10s; 
PL7
for e=1:snp 
 for k=1:snp 
c(e,k)=n(e)*n(k); 
C12(e,k)=ra*ca*det(J12)*int(int(int(c(e,k),x,-1,1),y,-1,1),z,-1,1); 
 end; 
end; 
C12s=single(C12);C26=C12s;C40=C12s;C54=C12s; 
for e=1:sp 
 index(e,:)=Nut(e,:); 
end; 
Y=cat(3,C1s,C2,C3s,C4,C5,C6s,C7,C8,C9,C10s,C11,C12s,C13,C14,C15,C16,C17,C
18,C19,C20,C21,C22,C23,C24,C25,C26,C27,C28,C29,C30,C31,C32,C33,... 
C34,C35,C36,C37,C38,C39,C40,C41,C42,C43,C44,C45,C46,C47,C48,C49,C50,C51,C
52,C53,C54,C55,C56); 
C=zeros(sn,sn); 
for e=1:sp 
 for row_index=1:sn 
 for m1=1:snp 
 if index(e,m1)==row_index 
 for col_index=1:sn 
 for m2=1:snp 
if index(e,m2)==col_index 
C(row_index,col_index)=C(row_index,col_index)+Y(m
1,m2,e); 
 end; 
 end; 
 end; 
 end; 
 end; 
 end; 
end; 
n3=1/4*(1+x)*(1-y); 
n4=1/4*(1-x)*(1-y); 
n7=1/4*(1+x)*(1+y); 
n8=1/4*(1-x)*(1+y); 
n(1)=0; 
n(2)=0; 
n(3)=n3; 
n(4)=n4; 
n(5)=0; 
n(6)=0; 
n(7)=n7; 
n(8)=n8; 
o11=diff(n3,x)*N(26,1)+diff(n4,x)*N(25,1)+diff(n7,x)*N(56,1)+diff(n8,x)*N
(55,1); 
o22=diff(n3,y)*N(26,2)+diff(n4,y)*N(25,2)+diff(n7,y)*N(56,2)+diff(n8,y)*N
(55,2); 
o12=diff(n3,y)*N(26,1)+diff(n4,y)*N(25,1)+diff(n7,y)*N(56,1)+diff(n8,y)*N
(55,1); 
o21=o12; 
J10m=[o11 o12;o21 o22]; 
for e=1:snp 
 for k=1:snp 
hvt(e,k)=n(e)*n(k); 
k10(e,k)=hd*det(J10m)*int(int(hvt(e,k),x,-1,1),y,-1,1); 
end; 
PL8
 end; 
k10s=single(k10);K10s=A10s-k10s; 
K11=K10s;K13=K10s;K14=K10s;K24=K10s;K25=K10s;K27=K10s;K28=K10s; 
K38=K10s;K39=K10s;K41=K10s;K42=K10s;K52=K10s;K53=K10s;K55=K10s;K56=K10s; 
p11=diff(n3,x)*N(28,1)+diff(n4,x)*N(27,1)+diff(n7,x)*N(58,1)+diff(n8,x)*N
(57,1); 
p22=diff(n3,y)*N(28,2)+diff(n4,y)*N(27,2)+diff(n7,y)*N(58,2)+diff(n8,y)*N
(57,2); 
p12=diff(n3,y)*N(28,1)+diff(n4,y)*N(27,1)+diff(n7,y)*N(58,1)+diff(n8,y)*N
(57,1); 
p21=p12; 
J12m=[p11 p12;p21 p22]; 
for e=1:snp 
 for k=1:snp 
hvt(e,k)=n(e)*n(k); 
k12(e,k)=hd*det(J12m)*int(int(hvt(e,k),x,-1,1),y,-1,1); 
end; 
 end; 
k12s=single(k12);K12s=A12s-k12s;K26=K12s;K40=K12s;K54=K12s; 
for e=1:sp 
 index(e,:)=Nut(e,:); 
end; 
X=cat(3,K1s,K2,K3s,K4,K5,K6s,K7,K8,K9,K10s,K11,K12s,K13,K14,K15,K16,K17,K
18,K19,K20,K21,K22,K23,K24,K25,K26,K27,K28,K29,K30,K31,K32,K33,... 
K34,K35,K36,K37,K38,K39,K40,K41,K42,K43,K44,K45,K46,K47,K48,K49,K50,K51,K
52,K53,K54,K55,K56); 
K=zeros(sn,sn); 
for e=1:sp 
 for row_index=1:sn 
 for m1=1:snp 
 if index(e,m1)==row_index 
 for col_index=1:sn 
 for m2=1:snp 
 if index(e,m2)==col_index 
K(row_index,col_index)=K(row_index,col_index)+X(m1,m2,e
); 
 end; 
 end; 
 end; 
 end; 
 end; 
 end; 
end; 
ff1=zeros(8,1); 
for e=1:snp 
 ff10(e,1)=-hd*det(J10m)*int(int(n(e),x,-1,1),y,-1,1); 
 ff12(e,1)=-hd*det(J12m)*int(int(n(e),x,-1,1),y,-1,1); 
end; 
ff2=ff1;ff4=ff1;ff5=ff1;ff15=ff1;ff16=ff1;ff18=ff1;ff19=ff1; 
ff29=ff1;ff30=ff1;ff32=ff1;ff33=ff1; 
ff43=ff1;ff44=ff1;ff46=ff1;ff47=ff1; 
PL9
ff3=ff1;ff17=ff3;ff31=ff3;ff45=ff3; 
ff6=ff1;ff7=ff6;ff8=ff6;ff9=ff6;ff20=ff6;ff21=ff6;ff22=ff6;ff23=ff6; 
ff34=ff6;ff35=ff6;ff36=ff6;ff37=ff6;ff48=ff6;ff49=ff6;ff50=ff6;ff51=ff6; 
ff11=ff10;ff13=ff10;ff14=ff10;ff24=ff10;ff25=ff10;ff27=ff10;ff28=ff10; 
ff38=ff10;ff39=ff10;ff41=ff10;ff42=ff10;ff52=ff10;ff53=ff10;ff55=ff10;ff5
6=ff10; 
ff26=ff12;ff40=ff12;ff54=ff12; 
for e=1:sp 
 index(e,:)=Nut(e,:); 
end; 
Z=cat(3,ff1,ff2,ff3,ff4,ff5,ff6,ff7,ff8,ff9,ff10,ff11,ff12,ff13,ff14,ff15
,ff16,ff17,ff18,ff19,ff20,ff21,ff22,ff23,ff24,ff25,ff26,ff27,ff28,ff29,ff
30,... 
ff31,ff32,ff33,ff34,ff35,ff36,ff37,ff38,ff39,ff40,ff41,ff42,ff43,ff44,ff4
5,ff46,ff47,ff48,ff49,ff50,ff51,ff52,ff53,ff54,ff55,ff56); 
F=zeros(sn,1); 
for e=1:sp 
 for row_index=1:sn 
 for m1=1:snp 
 if index(e,m1)==row_index 
F(row_index,1)=F(row_index,1)+Z(m1,1,e); 
 end; 
 end; 
 end; 
 end; 
Ttruoc=Tdau*ones(sn,1); 
VP=(C-denta_t*0.5*K)*Ttruoc+denta_t*0.5*F*(Tc+Tdau); 
VT=-(C+denta_t*0.5*K); 
Dapso=inv(VT)*VP;