Phân tích truyền nhiệt 3D kết cấu bê tông non tuổi bằng phần tử tứ diện nội suy kép

87 Tạp chí Khoa học Lạc Hồng Số Đặc Biệt
Journal of Science of Lac Hong University
Special issue (11/2017), pp. 87-92
Tạp chí Khoa học Lạc Hồng
Số đặc biệt (11/2017), tr. 87-92
PHÂN TÍCH TRUYỀN NHIỆT 3D KẾT CẤU BÊ TÔNG NON TUỔI BẰNG 
PHẦN TỬ TỨ DIỆN NỘI SUY KÉP
3D heat transfer analysis of early age concrete by a new enhanced gradient 
finite element
Nguyễn Đình Dư1, Nguyễn Bá Ngọc Thảo1, Bùi Quốc Tính2
1dinhdu85@gmail.com
1Khoa Kỹ Thuật Công Trình Trường Đại học Lạc Hồng, Đồng Nai, Việt Nam
2Department of Mechanical and Environmental Informatics, Tokyo Institute of Technology,
2-12-1-W8-22, Ookayama, Meguro-ku, Tokyo, Japan, Việt Nam
Đến tòa soạn: 08/06/2017; Chấp nhận đăng: 14/06/2017
Tóm tắt. Vết nứt hoặc những khuyết tật thường xuất hiện trong kết cấu giòn như bê tông. Những khiếm khuyết đó có thể do nhiều 
nguyên nhân khác nhau như tải cơ học, sự tác động từ môi tường. Ngoài ra, nhiều vết nứt trong bê tông có thể bắt nguồn từ sự
thay đổi thể tích nội tại hoặc phản ứng hóa học có hại làm thay đổi độ ẩm, nhiệt độ bê tông. Trong nghiên cứu này, một phần tử
hữu hạn nội suy kép (CTH4) dựa trên các thủ tục nội suy liên tiếp (CIP) với sự liên tục các Gradient nhiệt tại các nút cũng như 
trên cạnh phần tử cho bài toán truyền nhiệt 3D. Kết quả thu được là chính xác cho trường nhiệt độ, từ đó hạn chế được những 
phá hoại do nhiệt.
Từ khoá: Bê tông non tuổi; Phần tử TH4; Nội suy kép; Truyền nhiệt
Abstract. Defects or cracks are commonly observed in quasi-brittle materials like concrete structures. It is well-known to 
understand that cracks may be caused by different reasons due to, for instance, mechanical loading, some deleterious reactions or 
environmental loadings. Also, many of cracks in concrete may be traced to intrinsic volumetric changes or the deleterious chemical 
reactions, resulting in response to moisture, chemical, and thermal effects in concrete. In this work, we develop a new finite element 
based on consecutive-interpolation procedure (CTH4) for heat transfer analysis of early age concrete in three-dimension. 
Temperature distribution and its gradient over time in early age concrete will be analyzed through some representative numerical 
examples. The present numerical results are also validated against the conventional finite element (TH4).
Keywords: Early Age Concrete; Tetrahedral element; Consecutive-interpolation element; Heat transfer
1. GIỚI THIỆU
Kết cấu chỉ làm việc hiệu quả và đạt mong muốn của nhà 
thiết kế khi ở một nhiệt độ nhất định ban đầu hoặc có nhiệt 
độ có thể thay đổi nhưng trong phạm vi hẹp. Tuy nhiên, trong 
quá trình thủy phân của xi măng, một lượng nhiệt lớn được 
sinh ra dẫn đến sự chênh lệch nhiệt độ rất lớn giữa bờ mặt 
kết cấu và bên trong kết cấu. Đặc biệt là các kết cấu bê tông 
khối lớn. Sự chênh lệch nhiệt độ này là nguyên nhân chủ yếu 
gây ra các vết nứt cho kết cấu, là mầm mống gây phá hoại và 
giảm tuổi thọ công trình. Do đó, việc tính toán mô phỏng 
trường nhiệt độ cũng như trường biến dạng, ứng suất là cần 
thiết. 
Với hình thể và điều kiện biên phức tạp của bài toán thì 
việc tìm lời giải chính xác bằng phương pháp giải tích là 
không thể. Do đó, phương pháp số ra đời như một đáp ứng 
cần thiết, đặc biệt là phương pháp PTHH. Tuy nhiên, PTHH 
vẫn còn tồn tại một số vấn đề cần khắc phục như sự bất liên 
tục về Gradient tại nút và cạnh giữa các phần tử. Hiện nay, 
đã có nhiều phương pháp dựa trên nền tảng PTHH nhưng có 
sự cải tiến như XFEM, SFEMv.v. Trong nghiên cứu này, 
một đề xuất có tính kế thừa từ phương pháp PTHH nhưng 
được cải tiến hàm nội suy sao cho các Gradient nhiệt có tính 
liên tục tại nút và trên cạnh biên giữa các phần tử nhằm mô 
phỏng chính xác trường biến dạng nhiệt cũng như ứng suất 
do nhiệt. Sự cải tiến này được gọi là thủ tục nội suy kép 
(consecutive-interpolation), được giới thiệu đầu tiên bởi 
Zheng và các cộng sự (2011) cho phần tam giác (CT3)[1]. 
Tiếp đến là một loạt các công bố cho phần tử tứ giác (CQ4) 
được phát triển bởi Bui Quoc Tinh và các cộng sự, dễ dàng 
tìm thấy trong [1][2][3]. Tiếp theo sự thành công đó, trong 
phạm vi bài viết này, một nghiên cứu mới cho bài toán truyền 
nhiệt 3D được áp dụng cho móng bê tông khối lớn. Bài toán 
được mô phỏng bằng phần tử CTH4 và được lập trình trên 
phần mềm Matlab. Kết quả thu được từ nghiên cứu được so 
sánh với tài liệu tham khảo [5] nhằm kiểm chứng hiệu suất 
của phương pháp.
2. XÂY DỰNG PHẦN TỬ CTH4
2.1 Tổng quan nội suy kép
Trong bài toán 3D, xét một miền Î R3 và chịu ràng buộc 
bởi biên . Hàm số u(x) được xấp xỉ gần đúng theo thủ tục 
nội suy kép (CIP) như [1, 2, 3, 4]
(1)
Trong đó n là chỉ số của nút phần tử và u[I] là giá trị của hàm 
u(x) tại nút I được nội suy theo phương pháp PTHH truyền 
thống.
(2)
Trị số, , và là giá trị đạo hàm trung bình của 
hàm u(x) tại nút I. Giá trị đạo hàm đầu tiên của nút 
I chính là đạo hàm của hàm u(x) tại nút I thuộc phần tử e và 
được viết bởi PTHH chuẩn như bên dưới:
å
=
==
ne
l
xlxlI
e
x uNxu
1
,,
][
,
ˆˆ)( uN (3)
với ne là nút bắt đầu thuộc phần tử e. Sau khi lần lượt tính 
đạo hàm tại nút I cho tất cả các phần tử e Î SI có 
Nguyễn Đình Dư, Nguyễn Bá Ngọc Thảo, Bùi Quốc Tính
 88 Tạp chí Khoa học Lạc Hồng Số Đặc Biệt
chứa nút I, giá trị đạo hàm trung bình có thể được tính 
bằng trung bình hàm trọng số như bên dưới.
(4)
Trong phương trình (4), e là hàm trọng số được định 
nghĩa là tỷ số diện tích của phần tử e với tổng diện tích miền 
SI. Giá trị , được tính tương tự như .
Trong phương trình (1), các hàm I, Ix, Iy, Iz được gọi 
là hàm hổ trợ và phụ thuộc vào hình dáng phần tử nhưng 
chúng phải thỏa mãn điều kiện bên dưới [1, 2]:
I (xJ) = IJ , I,x (xJ) = 0 , I,y (xJ) = 0
Ix (xJ) = 0 , Ix,x (xJ) = IJ , Ix,y (xJ) = 0 (5)
Iy (xJ) = 0 , Iy,x (xJ) = 0 , Iy,y (xJ) = IJ
Những điều kiện ở phương trình (5) được diễn dãi chi tiết 
trong [2]. Phương trình (1) có thể được viết lại như sau
(6)
Trong CFEM, hàm dạng RI có liên quan đến nút I được 
viết như sau
(7)
Để đơn giản hơn khi nghiên cứu thủ tục nội suy kép (CIP) 
trong 3D, một ví dụ minh họa thủ tục CIP được áp đặt vào 
miền bài toán 2D như hình 1[1]. Một điểm cần nội suy x(x,y) 
trong phần tử tứ giác như hình 1, theo phương pháp PTHH 
truyền thống thì sẽ được nội suy từ 4 nút i, j, k, m; Tuy nhiên, 
khi áp dụng CFEM thông qua thủ tục CIP thì miền nội suy 
sẽ rộng hơn và được minh họa trong hình 1. Hàm dạng của 
phần tử CQ4 và Q4 được minh họa trong hình 2, dễ dàng 
nhận thấy CQ4 trơn và liên tục tại nút cũng như trên cạnh 
phần tử.
2.2 Phần tử CTH4
Xây dựng mới phần tử 3D tứ diện nội suy kép (CTH4)
được trình bày trong mục này. Từ phần tử tứ diện trong 
PTHH truyền thống, chúng tôi áp dụng thủ tục CIP để xây 
dựng lại hàm dạng có bậc cao hơn nhưng không làm tăng bậc 
tự do của bài toán. Mô hình phần tử tứ diện trong không gian 
vật lý và không gian tự nhiên được thể hiện như Hình 3. Miền 
hổ trợ của phần tử CTH4 được minh họa như trong Hình 4.
(a)
(b)
Hình 2. Hàm dạng của Q4 (a), CQ4 (b) trong 2D
Hình 1. Hình minh họa CQ4 trong 2D
89 
Nguyễn Đình Dư, Nguyễn Bá Ngọc Thảo, Bùi Quốc Tính
Tạp chí Khoa học Lạc Hồng Số Đặc Biệt
Hình 3. Phần tử TH4 trong không gian vật lý và tự nhiên
Chúng ta định nghĩa bốn nút của phần tử tứ diện là i, j, k, 
m và hàm dạng của 4 nút được viết như sau:
(8)
(9)
(10)
(11)
các đạo hàm từng phần được tính như sau:
(12)
trong đó J là ma trận Jacobian được tính toán như sau:
ú
ú
ú
ú
û
ù
ê
ê
ê
ê
ë
é
´
ú
ú
ú
ú
ú
ú
ú
û
ù
ê
ê
ê
ê
ê
ê
ê
ë
é
¶
¶
¶
¶
¶
¶
¶
¶
¶
¶
¶
¶
¶
¶
¶
¶
¶
¶
¶
¶
¶
¶
¶
¶
=
ú
ú
ú
ú
ú
ú
ú
û
ù
ê
ê
ê
ê
ê
ê
ê
ë
é
¶
¶
¶
¶
¶
¶
¶
¶
¶
¶
¶
¶
¶
¶
¶
¶
¶
¶
=
mmm
kkk
jjj
iii
mkji
mkji
mkji
zyx
zyx
zyx
zyx
LLLL
LLLL
LLLL
zyx
zyx
zyx
J
(13)
Các hàm hỗ trợ như “(5)” được viết bởi
(14)
(15)
(16)
(17)
với p = 0,5. Các hàm số j , jx , jy , jz ; k , kx , ky , kz ;
m , mx , my , mz được viết tương tự bằng cách xoay vòng 
các chỉ số i , j , k , m.
3. PHƯƠNG TRÌNH TRUYỀN NHIỆT VÀ MÔ 
HÌNH THỰC TẾ
3.1 Phương trình chủ đạo quá trình truyền nhiệt
Phương trình vi phân chủ đạo của quá trính truyền nhiệt 
3D ở dạng tổng quát như [5]
(18)
với các điều kiện biên được cho như bên dưới
, trên : Biên nhiệt độ đã biết (19)
, trên : Biên truyền nhiệt (20)
, trên 3G : Biên đối lưu (21)
Trong phương trình (18), k = diag(kxx , kyy , kzz) là ma trận 
hệ số dẫn nhiệt, T là trường nhiệt độ, Q là dòng nhiệt, là 
khối lượng riêng, c là nhiệt dung riêng. Trong các phương 
trình (19-21), T0 là nhiệt độ đã biết; q0 là dòng nhiệt đã biết; 
T là nhiệt độ môi trường; n(nx, ny, nz) là cosin chỉ phương.
Các dạng yếu của vấn đề truyền nhiệt thu được bằng cách 
nhân hai vế phương trình (8) với một hàm thử T và tích 
phân trên toàn miền
(22)
Áp dụng lý thuyết Gaussian, áp đặt các điều kiện biên và 
thực hiện các phép biến đổi ta có:
(23)
Sau khi trường nhiệt độ được xác định, năng lượng nhiệt 
trên toàn miền được tính toán thông qua công thức sau:
(24)
trong đó là nhiệt độ tại nút và B là mà ma trận đạo hàm 
của hàm dạng
Nguyễn Đình Dư, Nguyễn Bá Ngọc Thảo, Bùi Quốc Tính
 90 Tạp chí Khoa học Lạc Hồng Số Đặc Biệt
Hình 4. Miền hổ trợ của phần tử CTH4
(25)
Trong phương trình (25), RI là hàm dạng phần tử CTH4 
được thể hiện trong “(7)”.
3.2 Mô hình thực tế
Trong quá trình ninh kết bê tông, một lượng nhiệt rất lớn 
tỏa ra do quá trình thủy phân của xi măng. Đặc tính của bê 
tông là dẫn nhiệt kém và diện tích trao đổi nhiệt với môi 
trường là nhỏ cùng với kết cấu là bê tông khối lớn nên lượng 
nhiệt này tích tụ bên trong kết cấu bê tông và tạo nên sự
chênh lệch giữa tâm kết cấu và bờ mặt bên ngoài. Các thông 
số bài toán như nhiệt sinh ra do quá trình thủy phân, nhiệt 
độ môi trường, nhiệt tại các biên, nhiệt độ bê tông khi đổ
(nhiệt độ ban đầu) dễ dàng tìm thấy trong [6].
Bảng 1. Các thông số vật liệu trong mô hình
Bê Tông Đài 
Móng B40
Nền 
Đất
Tỷ nhiệt C, (kcal/kg0C) 0.27 0.2
Khối lượng thể tích, 
(kG/m3)
2400 1800
Hệ số dẫn nhiệt, 
(kcal/m.h.0C)
2.5 1.7
Hệ số trao đổi nhiệt, 
(kcal/m2.h.0C)
12 12
Nhiệt độ bê tông khi đổ (0C) 30
Hằng số hàm tăng nhiệt độ
đoạn nhiệt (0C)
K= 59.60C, a
= 1.113
Nhiệt độ môi trường (0C) 28.7
4. KẾT QUẢ SỐ
Một khối móng bê tông có kích thước 4600x4600x4000
đặt trên nền đất như Hình 3 là mô hình được phân tích trong 
nghiên cứu này. Do tính đối xứng về vật liệu và các điều 
kiện biên cũng như để giảm thời gian phân tích nên ¼ mô 
hình được phân tích. Khối hình học bê tông và đất được chia 
lưới có quy tắc như Hình 4. Tổng thời gian phân tích của bài 
toán là 700h, đây là khoảng thời gian mà nhiệt độ trong khối 
bê tông giảm gần ngang với nhiệt độ môi trường. Bước thời 
gian trong mỗi lần phân tích càng nhỏ thì sai số bài toán 
càng nhỏ đồng nghĩa với vấn đề sẽ tốn tài nguyên cũng như 
thời gian máy tính, do đó trong nghiên cứu này chọn bước 
thời gian phân tích là 1,75h. Các số liệu của bài toán được 
cho như trong Bảng 1.
Hình 
5. Mô hình thực tế và kích thước hình học 
Hình 6. Chia lưới và vị trí điểm cần khảo sát
91 
Nguyễn Đình Dư, Nguyễn Bá Ngọc Thảo, Bùi Quốc Tính
Tạp chí Khoa học Lạc Hồng Số Đặc Biệt
Hình 7. Biểu đồ minh họa nhiệt độ tại A theo thời gian
Hình 8. Biểu đồ minh họa sự chênh lệch nhiệt độ tại A và B theo 
thời gian
Hình 9. Biểu đồ minh họa sự chênh lệch nhiệt độ tại A và C theo 
thời gian
Nhằm thấy rõ được sự chênh lệch nhiệt độ theo thời gian
giữa vùng biên và tâm của khối móng nên nhiệt độ tại ba 
điểm A, B, C được khảo sát như hình 6. Hình 7 minh họa 
nhiệt độ tại điểm A theo thời gian với bước thời gian phân 
tích là 1,75h. Trong khoảng từ lúc bắt đầu đổ đến 98h thì 
nhiệt độ tăng một cách nhanh chóng và đạt cực đại sau thời 
gian khoảng 98h. Đây là thời điểm tương ứng với quá trình 
thủy hóa đang diễn ra rất mạnh, lượng nhiệt sinh ra lớn hơn 
rất nhiều so với lượng nhiệt tỏa ra môi trường. Sau đó nhiệt 
độ bắt đầu giảm dần theo thời gian do phản ứng thủy phân 
đi đến giai đoạn kết thúc và sẽ tiến sát đến nhiệt độ môi 
trường sau khoảng thời gian 700h.
Hình 10. Biểu đồ nhiệt độ tại tâm và biên khối móng, giữa phân tích số và thí nghiệm thực tế, tham khảo [5]
Hình 8 và Hình 9 minh họa sự khác biệt giữa nhiệt độ tại 
tâm khối móng và bờ mặt khối móng. Một trong những 
nguyên nhân gây ra nứt là do sự chênh lệch nhiệt độ giữa 
bên trong và bên ngoài khối móng. Trong khoảng 98h đầu, 
Phân tích truyền nhiệt 3D kết cấu bê tông non tuổi bằng phần tử tứ diện nội suy kép
 92 Tạp chí Khoa học Lạc Hồng Số Đặc Biệt
nhiệt độ chênh lệch rất lớn, giữa điểm A và điểm B là 38oC 
còn giữa điểm A và điểm C là 42oC. Có sự khác biệt này là 
do khoảng cách giữa hai điểm A và C là lớn hơn giữa A và 
B. Nhiệt độ tại bờ mặt có sự tăng nhẹ và đạt cực đại trễ pha
so với nhiệt độ tại tâm khối móng do khoảng cách về không 
gian tiếp xúc với môi trường.
Trong Hình 10, được tham khảo tại [5], nhiệt độ tại tâm 
khối móng và mặt biên khối móng được thể hiện bằng phân 
tích số (Phần mềm Midas) và mô hình thí nghiệm thực tế.
Có thể nhận thấy quy luật thay đổi nhiệt độ giữa [5] và 
nghiên cứu này là giống nhau. Tuy nhiên có sự khác nhau 
tương đối về nhiệt độ cực đại là do bước thời gian phân tích 
cũng như quan niệm về các điều kiện biên khi phân tích.
5. KẾT LUẬN
Kết quả phân tích trường nhiệt độ bằng phương pháp 
nghiên cứu mô tả đúng được qui luật nhiệt độ trong khối bê 
tông do quá trình thủy phân. Từ đó đề xuất được các biện 
pháp thi công hợp lý để hạn chế các vết nứt trong bê tông 
nhằm tránh phá hoại diễn ra trong lòng kết cấu. Đặc biệt là 
bê tông khối lớn, thường là kết cấu của những công trình 
quan trọng như móng nhà cao tầng, mố trụ cầu, đường băng 
sân bay...v.v. Tuy nhiên, để rõ ràng hơn thì cần phải phân 
tích thêm trường ứng suất theo thời gian để xác định được 
chính xác thời điểm xảy ra vết nứt và đó là hướng cho những 
nghiên cứu tiếp theo và phân tích những kết cấu phức tạp 
hơn.
6. TÀI LIỆU THAM KHẢO
[1] Bui QT, Vo QD, Zhang Ch, Nguyen DD, “A consecutive-
interpolation quadrilateral element (CQ4): Formulation and 
applications”, Finite elements in analysis and design, 84, 14-
31, 2014.
[2] Zuoyi Kang, Tinh Quoc Bui, Du Dinh Nguyen, Takahiro 
Saitoh, Sohichi Hirose, “An extended consecutive-
interpolation quadrilateral element (XCQ4) applied to linear 
elastic fracture mechanics”, Acta Mech, Springer-Verlag 
Wien, Volume 226, issue 12, pp 3991–4015, 2015.
[3] Tinh Quoc Bui, Du Dinh Nguyen, Xiaodong Zhang, Sohichi 
Hirose, Romesh C. Batra, “Analysis of 2-dimensional transient 
problems for linear elastic and piezoelectric structures using 
the consecutive-interpolation quadrilateral element (CQ4)”,
European Journal of Mechanics A/Solids, 58, pp 1-19, 2016.
[4] Bùi Đức Vinh, Nguyễn Văn Chánh, Chu Quốc Thắng “Phân 
tích phân bốnhiệt hydrat và ứng suất trong cấu trúc bê tông để 
kiểm soát sựgây nứt của công trình bê tông cốt thép”, Tuyển 
tập công trình khoa học - Hội nghị khoa học sự cố và hư hỏng 
công trình xây dựng, 2012.
[5] Hồ Ngọc Khoa, Vũ Chí Công, “Phân tích trường nhiệt độvà 
ứng suất nhiệt trong bê tông khối lớn bằng phương pháp phần 
tử hữu hạn”, Tạp chí KHCN Xây Dựng, Số 14/12, 2012.
TIỂU SỬ TÁC GIẢ
Nguyễn Đình Dư
Năm sinh 1985, Tuy Phước, Bình Định. Tốt nghiệp Đại học và Thạc sĩ tại Trường Đại học
Bách khoa Thành phố Hồ Chí Minh năm 2009 và 2013. Giảng viên cơ hữu khoa Kỹ thuật
Công trình Trường Đại học Lạc Hồng. Lĩnh vực nghiên cứu: Cơ học vật rắn biến dạng, phương
pháp số, động lực học kêt cấu.
Nguyễn Bá Ngọc Thảo
Năm sinh 1986, Phan Thiết, Bình Thuận. Tốt nghiệp Đại học và Thạc sĩ tại trường đại học Lạc
Hồng, Đồng Nai năm 2009 và 2013. Giảng viên cơ hữu khoa Kỹ thuật Công trình Trường Đại
học Lạc Hồng. Lĩnh vực nghiên cứu: Betong cốt thép cường độ cao.