Tính truyền nhiệt khi cấp đông cryo bề mặt phẳng vật ẩm bằng phương pháp sai phân hữu hạn

Chia sẻ: _ _ | Ngày: | Loại File: PDF | Số trang:4

Thêm vào BST

Báo xấu

14
lượt xem 3
download

Download Vui lòng tải xuống để xem tài liệu đầy đủ

Bài viết Tính truyền nhiệt khi cấp đông cryo bề mặt phẳng vật ẩm bằng phương pháp sai phân hữu hạn trình bày tính truyền nhiệt khi cấp đông cho bề mặt phẳng vật ẩm bằng phương pháp sai phân hữu hạn, qua đó tính độ dày của lớp băng được hình thành theo thời gian cấp lạnh.

Chủ đề:

Bình luận(0) Đăng nhập để gửi bình luận!

Lưu

Nội dung Text: Tính truyền nhiệt khi cấp đông cryo bề mặt phẳng vật ẩm bằng phương pháp sai phân hữu hạn

44 Thái Ngọc Sơn, Hoàng Ngọc Đồng, Lê Minh Trí TÍNH TRUYỀN NHIỆT KHI CẤP ĐÔNG CRYO BỀ MẶT PHẲNG VẬT ẨM BẰNG PHƯƠNG PHÁP SAI PHÂN HỮU HẠN USING FINITE DIFFERENCE METHODS TO CALCULATE HEAT TRANSFER IN FREEZING THE FLAT SURFACE OF HUMID MATERIALS BY CRYOGENIC TECHNIQUE Thái Ngọc Sơn1, Hoàng Ngọc Đồng1, Lê Minh Trí2 1 Trường Đại học Bách khoa - Đại học Đà Nẵng; thaingocson@gmail.com, hndong@dut.udn.vn 2 Trường Cao đẳng Công nghiệp Huế; lmtri@hueic.edu.vn Tóm tắt - Một trong những cách để tiêu diệt tế bào ung thư là làm Abstract - One of the methods to destroy cancer cells is freezing them lạnh nó đến nhiệt độ rất thấp bằng kỹ thuật cryo. Để thực hiện liệu to very low temperatures using cryogenic technique. To perform this pháp này, cần xác định tốc độ đóng băng và phân bố nhiệt độ trong technique, we need to answer these questions: How process of khối u hóa rắn. Trong [1], các tác giả đã sử dụng phương pháp giải freezing tumor occurs, how fast the freezing and distributing of tích để giải quyết các vấn đề trên. Bài báo này trình bày tính truyền temperature in the ice is, etc. In [1], the authors use calculus and nhiệt khi cấp đông cho bề mặt phẳng vật ẩm bằng phương pháp analysis to solve the above problems. This paper presents the heat sai phân hữu hạn, [2]qua đó tính độ dày của lớp băng được hình transfer when freezing the flat surface of humid materials using finite thành theo thời gian cấp lạnh. Đây là bài toán cơ sở giúp các nhà difference method, thereby, calculating the thickness of the ice sheet khoa học xác định thời gian làm lạnh cần thiết để tiêu diệt có định formed during the freezing time. This is a fundamental problem that hướng các tế bào ung thư bằng kỹ thuật cryo. Tác giả đã thiết lập helps scientists determine the freezing time needed to destroy cancer phương trình, xây dựng phần mềm trong môi trường Excel, dùng cells with cryogenic technique. The author has set up equations, built VBA để giải bài toán. Kết quả tính toán được so sánh với kết quả software in the Excel, and used VBA to solve the problem. The results tính theo phương pháp giải tích và phương pháp số khác. Sai số are compared with the ones under analysis and other calculus kết quả tính giữa các phương pháp là nhỏ hơn 15%. methods. The calculation error between methods is less than 15%. Từ khóa - lỏng nitơ; truyền nhiệt; ung thư; kỹ thuật cryo; phương Key words - liquid nitrogen; heat transfer; cancer; cryogenic pháp sai phân hữu hạn technique; finite difference method 1. Đặt vấn đề ẩm = 75%; nhiệt chuyển pha của nước rc= 334.100J/kg. Cấp đông cryo bề mặt phẳng vật ẩm là kỹ thuật sử dụng Vật được cấp đông bằng thiết bị áp lạnh với nguồn thu nhiệt khí hóa lỏng có nhiệt độ rất thấp phun hoặc cho tiếp xúc là lỏng nitơ có nhiệt độ tf= -196°C. Hãy xác định độ dày trực tiếp với bề mặt phẳng vật ẩm. Kỹ thuật này được sử lớp băng được hình thành theo thời gian (Hình 1). dụng trong y học để điều trị bệnh sẹo lồi [3], bệnh chàm mạn tính tăng sừng [4], bệnh rối loạn da lành tính và ác tính [[5]5], bệnh mụn cóc, nốt ruồi, thẻ da, keratoses năng lượng mặt trời, u dây thần kinh Morton và ung thư da nhỏ [6] … Nhằm tiến đến phá hủy có kiểm soát các tổ chức bệnh lý và các tế bào ung thư, chúng tôi tính truyền nhiệt khi cấp đông cryo vật ẩm, qua đó xác định thời gian cấp lạnh cần thiết để làm đóng băng hóa rắn một khối tế bào có độ dày cho trước. Đây là việc làm hết sức cần thiết và có ý nghĩa trong ứng dụng kỹ thuật lạnh cryo vào lĩnh vực y tế. Hình 1. Mô hình bài toán Chúng ta có thể tính độ dày của lớp băng tạo ra bằng 2.2. Giải bài toán phương pháp giải tích với một số giả thiết [1], trong đó có giả thiết: Xem các thông số vật lý (hệ số dẫn nhiệt, nhiệt Đây có thể xem như phần mở rộng của bài toán tính dung riêng...) là không đổi trong suốt quá trình đóng băng. truyền nhiệt của thiết bị áp lạnh trong [1], nghĩa là khảo Tuy nhiên, ở đây chúng tôi giải bài toán này theo cách sát trường nhiệt độ của tấm phẳng rộng vô hạn với biên khác: sử dụng phương pháp sai phân hữu hạn, so sánh kết loại 1 với các thông số vật lý phụ thuộc nhiệt độ. Từ thời quả tính giữa 2 phương pháp. Khi giải bài toán theo phương điểm τ, trong quãng thời gian Δτ, 1 lớp vật ẩm có bề dày pháp số, trong đó có phương pháp sai phân hữu hạn, các δ sẽ hạ nhiệt độ từ t 1 xuống đến t0 rồi sau đó hóa rắn. Gọi thông số vật lý có thể tính toán theo nhiệt độ, qua đó làm quãng thời gian để mô sống hạ nhiệt độ từ t 1 xuống đến giảm sai số so với phương pháp giải tích. t0 là Δτh, ta có: Lượng nhiệt thải ra khi vật ẩm hạ nhiệt độ từ t1 xuống đến t0: 2. Kết quả nghiên cứu và khảo sát 𝑞ℎ = 𝛿𝜌𝑖 𝐶𝑖 (𝑡1 − 𝑡0 ), J/m2 2.1. Mô hình bài toán Lượng nhiệt ẩn thải ra khi vật ẩm hóa rắn: Giả sử vật ẩm có các thông số vật lý: nhiệt độ đóng băng 𝑞𝑎 = 𝛿𝜌𝑖 𝜔𝑟𝑐 , J/m2 t0= 0°C; nhiệt độ ban đầu t1= 37°C; khối lượng riêng i [kg/m3]; nhiệt dung riêng Ci [J/kgK]; hệ số dẫn nhiệt khi Phương trình dẫn nhiệt từ lớp mô qua lớp băng đã hình nhiệt độ nhỏ hơn 0°C i [W/mK]; độ ẩm toàn phần của vật thành trước thời điểm τ, có nhiệt độ t, nhiệt trở R có dạng:
ISSN 1859-1531 - TẠP CHÍ KHOA HỌC VÀ CÔNG NGHỆ ĐẠI HỌC ĐÀ NẴNG, SỐ 9(118).2017 - Quyển 2 45 𝑡1 − 𝑡 b. Quá trình hình thành lớp băng thứ N=2 từ thời Δ𝜏ℎ = 𝛿𝜌𝑖 𝐶𝑖 (𝑡1 − 𝑡0 ) 𝛿 điểm τ1=Δτ đến thời điểm τ2= τ1+Δτ=2Δτ +𝑅 𝜆𝑖 𝑡0 − 𝑡 Trong quãng thời gian từ thời điểm τ1 đến thời điểm τ2, 𝛿 Δ𝜏𝑎 = 𝛿𝜌𝑖 𝜔𝑟𝑐 trường nhiệt độ của lớp băng thứ N=1 (hình thành từ thời { 𝜆0 + 𝑅 điểm τ0 đến thời điểm τ1) thay đổi. Vì quãng thời gian Δτ Ký hiệu: rất nhỏ, nhiệt độ thay đổi từ tf đến t0 rất lớn nên có thể giả thiết trường nhiệt độ thay đổi theo quy luật tuyến tính. 𝐸1 ≜ 𝜌𝑖 𝐶𝑖 (𝑡1 − 𝑡0 ) Khi đó, sử dụng công thức (1) để tính bề dày lớp băng 𝐸2 ≜ 𝜌𝑖 𝜔𝑟𝑐 thứ N=2, nhiệt trở R là nhiệt trở dẫn nhiệt qua lớp băng Biết: chiều dày δ1. 𝛥𝜏 = 𝛥𝜏ℎ + 𝛥𝜏𝑎 𝐸1 𝐸2 𝐸= + Suy ra: 𝜆𝑖 (𝑡1 − 𝑡𝑓 ) 𝜆0 (𝑡0 − 𝑡𝑓 ) 𝛿𝐸1 𝛿 𝛿𝐸2 𝛿 𝛿1 𝐸1 𝐸2 ( + 𝑅) + ( + 𝑅) = 𝛥𝜏 𝐹= ( + ) 𝑡1 − 𝑡 𝜆𝑖 𝑡0 − 𝑡 𝜆0 𝜆1 𝑡1 − 𝑡𝑓 𝑡0 − 𝑡𝑓 𝐸1 𝐸2 𝐸1 𝑅 𝐸2 𝑅 Trong đó, λ1 là hệ số dẫn nhiệt trung bình theo nhiệt độ 𝛿2 [ + ]+𝛿[ + ] = 𝛥𝜏 𝜆𝑖 (𝑡1 − 𝑡) 𝜆0 (𝑡0 − 𝑡) 𝑡1 − 𝑡 𝑡0 − 𝑡 trung bình của lớp băng thứ N=1 (W/mK). Ký hiệu: Thay vào (1) tính được bề dày lớp băng thứ N=2, hình 𝐸1 𝐸2 thành từ thời điểm τ1 đến thời điểm τ2. 𝐸≜ + 𝜆𝑖 (𝑡1 − 𝑡) 𝜆0 (𝑡0 − 𝑡) 𝐸1 𝑅 𝐸2 𝑅 𝐸1 𝐸2 𝐹≜ + = 𝑅( + ) 𝑡1 − 𝑡 𝑡0 − 𝑡 𝑡1 − 𝑡 𝑡0 − 𝑡 Giải phương trình bậc 2, khử nghiệm âm do không thích hợp, ta có phương trình: −𝐹+√𝐹2 +4𝐸Δ𝜏 𝛿= (1) 2𝐸 Giá trị δ là bề dày của lớp băng hình thành từ thời điểm τ đến thời điểm τ +Δτ, phụ thuộc vào nhiệt trở và trường nhiệt độ của lớp băng đã hình thành trước thời điểm τ. Phương pháp sai phân hữu hạn giúp ta tìm được trường nhiệt độ trong lớp băng hình thành trước thời điểm τ. Cụ Hình 3. Trường nhiệt độ lớp băng thứ N=2 thể các bước tiến hành như sau: c. Quá trình hình thành lớp băng thứ N (N≥3) từ a. Quá trình hình thành lớp băng thứ N=1 từ thời thời điểm τN-1=(N-1)Δτ đến thời điểm τN=τN-1+Δτ=NΔτ điểm τ0=0 đến thời điểm τ1=Δτ Trong quãng thời gian từ thời điểm τN-1 đến thời điểm Trong công thức (1) dùng để tính δ, do hệ số tỏa nhiệt τN, trường nhiệt độ của lớp băng hình thành từ thời điểm τ0  khi chất lỏng đang sôi rất lớn, vách ngăn mỏng và có hệ đến thời điểm τN-1 với bề dày đã thay đổi: số dẫn nhiệt  lớn nên nhiệt trở R rất nhỏ, có thể bỏ qua. 𝑁−1 Do đó, F =0. 𝛿 = ∑ 𝛿𝑖 Thay vào (1), ta có bề dày δ1 của lớp băng thứ N=1 hình 𝑖=1 thành từ thời điểm τ0=0 đến thời điểm τ1=Δ τ: Giá trị Δτ rất nhỏ (chọn Δτ=0,1s) nên các giá trị δi cũng Δ𝜏 rất nhỏ. Chia bề dày lớp băng đó thành N-1 phần, gọi là sai 𝛿1 = √ phân theo tọa độ x: 𝐸 𝛿 𝛥𝑥𝑁 = Trường nhiệt độ tại thời điểm τ1 có dạng: 𝑁−1 Ký hiệu tọa độ của các mặt phân chia là xi, 𝑖 ∈ [0, 𝑁 − 1] ⊂ ℕ0 . Nhiệt độ tại các tọa độ tương ứng tại thời điểm τj là ti,j. Phương trình vi phân dẫn nhiệt trong lớp băng: 𝜕𝑡 𝜌𝑟 𝐶𝑟 = 𝛻(𝜆𝑟 𝛻𝑡) 𝜕𝜏 Trong đó: ρr, Cr, λr lần lượt ứng với khối lượng riêng, nhiệt dung riêng và hệ số dẫn nhiệt của băng. Trong khoảng thời gian Δτ đang xét, xem các thông số đó là không đổi và phụ thuộc theo nhiệt độ trung bình trong mỗi phần của lớp Hình 2. Trường nhiệt độ lớp băng thứ N=1 băng. Khi đó phương trình vi phân sẽ có dạng:
46 Thái Ngọc Sơn, Hoàng Ngọc Đồng, Lê Minh Trí 𝜕𝑡 2.3. Kết quả giải bằng phương pháp sai phân hữu hạn = 𝑎𝑟 𝛻 2 𝑡 Dữ liệu tính toán được thể hiện trong Bảng 1. Các thông 𝜕𝜏 Trong đó, ar là hệ số khuếch tán nhiệt của băng (m2/s). số vật lý được tính toán phụ thuộc theo nhiệt độ trung bình của mỗi lớp băng theo công thức trong tài liệu [8]. Kết quả Chuyển phương trình vi phân thành phương trình sai kích thước lớp băng tạo ra được thể hiện ở Bảng 2. phân, ta có: Bảng 1. Dữ liệu tính toán khi giải bằng phương pháp 𝛥𝑥 𝑡 𝛥𝜏 𝛥𝜏 𝑡 sai phân hữu hạn = 𝑎𝑟 ( ) 𝛥𝜏 𝛥𝑥 𝛥𝑥 Thông số Ký hiệu Trị giá Đơn vị Trong đó, các chỉ số x, τ biểu diễn biến không đổi khi Khoảng thời gian khảo sát τ 10 s tính độ chênh nhiệt độ. Phương trình sai phân theo sơ đồ ẩn cho điểm ti, j bất kỳ không nằm trên biên của vùng khảo Bước thời gian tính toán Δτ 0,1 s sát sẽ có dạng: Bước thời gian in kết quả ΔτIn 1 s 𝑡𝑖,𝑗 − 𝑡𝑖,𝑗−1 𝑡𝑖−1,𝑗 − 2𝑡𝑖,𝑗 + 𝑡𝑖+1,𝑗 Sai số tính toán nhiệt độ εΔt 0,01 °C = 𝑎𝑟 𝛥𝜏 𝛥𝑥 2 Bảng 2. Độ dày của lớp băng được hình thành theo thời gian Với ar xác định theo nhiệt độ trung bình: Thời gian [s] δs [mm] 𝑡𝑖,𝑗−1 + 𝑡𝑖,𝑗 1 0,907 𝑡𝑖 = 2 2 1,316 Ta có thể lập hệ gồm N-2 phương trình với N-2 ẩn ti, j 3 1,629 (i=1÷N-2) và giải. Tuy nhiên, do hệ số khuếch tán nhiệt ar 4 1,890 cũng phụ thuộc nhiệt độ nên giải bài toán rất phức tạp. 5 2,116 Do vậy, chúng tôi đã xây dựng phần mềm trong môi trường 6 2,317 Excel, dùng VBA để giải bài toán theo phương pháp lặp để xác định trường nhiệt độ trong lớp băng. 7 2,495 8 2,655 Sử dụng công thức (1) để tính bề dày lớp băng thứ N, với R là nhiệt trở dẫn nhiệt qua lớp băng chiều dày ΔxN, 9 2,799 nhiệt độ 2 bên bề mặt là tN-2, N, tN-1, N=t0. Khi đó: 10 2,928 𝐸1 𝐸2 So sánh kết quả tính giữa phương pháp sai phân 𝐸= + hữu hạn và phương pháp giải tích 𝜆𝑖 (𝑡1 − 𝑡𝑁−2,𝑁 ) 𝜆0 (𝑡0 − 𝑡𝑁−2,𝑁 ) Tiến hành giải bài toán bằng phương pháp giải tích theo 𝛥𝑥𝑁 𝐸1 𝐸2 𝐹= ( + ) [1], so sánh kết quả tính toán giữa hai phương pháp, ta lập 𝜆𝑁−1 𝑡1 − 𝑡𝑁−2,𝑁 𝑡0 − 𝑡𝑁−2,𝑁 được Bảng 3. Trong đó, λN-1 là hệ số dẫn nhiệt của lớp băng, tính theo Bảng 3. So sánh kết quả tính độ dày lớp băng tạo ra giữa phương nhiệt độ trung bình của lớp băng (W/mK). pháp sai phân hữu hạn và phương pháp giải tích theo tài liệu [1] 𝑡𝑁−2,𝑁 + 𝑡𝑁−1,𝑁 Sai số 𝑡𝑁−1 = Thời Giải tích Sai phân Sai số tuyệt đối 2 tương đối gian [s] x()[mm] s [mm] abs(-x) [mm]  [%] 1 1,080 0,907 0,173 16,0 2 1,528 1,316 0,212 13,9 3 1,871 1,629 0,242 12,9 4 2,161 1,890 0,271 12,5 5 2,415 2,116 0,299 12,4 6 2,646 2,317 0,329 12,4 7 2,858 2,495 0,363 12,7 8 3,055 2,655 0,400 13,1 9 3,241 2,799 0,442 13,6 10 3,416 2,928 0,488 14,3 Sai số trung bình 0,300 13,4 Nhận xét Hình 4. Trường nhiệt độ lớp băng thứ N≥3 Sai số kết quả tính toán giữa phương pháp giải tích theo Thay các giá trị E, F vào (1), tính được bề dày lớp băng [1] và phương pháp sai phân hữu hạn đã nêu trên là 13,4%. thứ N hình thành từ thời điểm τN-1 đến thời điểm τN. Trong một phạm vi nhất định đối với bài toán nhiệt, sai số Bề dày lớp băng hình thành từ thời điểm τ0=0 đến thời trên có thể chấp nhận. điểm τN là: So sánh kết quả tính giữa phương pháp sai phân 𝑁 hữu hạn và phương pháp số 𝛿𝑠 = ∑ 𝛿𝑖 So sánh kết quả xác định trường nhiệt độ theo phương 𝑖=1
ISSN 1859-1531 - TẠP CHÍ KHOA HỌC VÀ CÔNG NGHỆ ĐẠI HỌC ĐÀ NẴNG, SỐ 9(118).2017 - Quyển 2 47 pháp sai phân hữu hạn với các thông số vật lý tương đồng cấp đông cho bề mặt phẳng bằng phương pháp sai phân với phương pháp số khi phun lỏng lên bề mặt da [8]. Kết hữu hạn. Khi tính toán sử dụng các thông số vật lý phụ quả giải bằng phương pháp số khi phun lỏng lên bề mặt da thuộc theo nhiệt độ trung bình của mỗi lớp băng, khi nhiệt thể hiện trên Hình 5. độ lớp băng thay đổi thì các thông số vật lý thay đổi theo. Có thể thấy, bề mặt phun lỏng không phải là dạng mặt Trong khi theo phương pháp giải tích, tất cả các thông số phẳng, bán kính lỏng phun ra khá bé (quãng 3 mm), nên vật lý đều là hằng số trong suốt quá trình tính toán. Đây thực chất mô hình tính toán theo 2 phương pháp sai phân cũng chính là ưu điểm của các phương pháp số, trong đó hữu hạn của tác giả và phương pháp số theo [8] không có phương pháp sai phân hữu hạn đối với phương pháp tương đồng về mặt hình học. Ngoài ra, do không có các mô giải tích. hình khác thích hợp hơn nên ta lấy kích thước lớp băng tại 2. Sai số kết quả tính giữa phương pháp giải tích và vị trí trục ống phun r=0 để so sánh. Kết quả so sánh thể phương pháp sai phân hữu hạn là 13,4%. Điều đó cho thấy, hiện trong Bảng 4. các bước tiến hành cũng như phương trình thiết lập khi giải bài toán truyền nhiệt cấp đông cryo bề mặt phẳng vật ẩm là chính xác và có độ tin cậy cao. 3. Sai số kết quả tính toán giữa phương pháp số theo [8] và phương pháp sai phân hữu hạn của tác giả nhỏ hơn 15% trong phạm vi so sánh, do có 1 số điểm không tương đồng về mặt hình học, điều kiện biên của hai mô hình. Sai số trên có thể chấp nhận được trong 1 bài toán nhiệt, cho thấy bài toán đã được giải thành công theo phương pháp sai phân hữu hạn. TÀI LIỆU THAM KHẢO [1] Hoàng Ngọc Đồng, Nguyễn Bốn, Lê Minh Trí, “Quá trình truyền Hình 5. Trường nhiệt độ giải theo phương pháp số nhiệt, truyền chất không ổn định với biên di động khi hóa rắn vật ẩm theo F. Sun [8] bằng kỹ thuật cryo”, Tạp chí Năng lượng Nhiệt, năm thứ 21, Số 117, 2014, trang 15-20. Bảng 4. So sánh kết quả hình thành trường nhiệt độ theo thời [2] Исаченко В.П., Осипова В.А., Сукомел А.С., Теплопередача, gian theo phương pháp sai phân hữu hạn và tài liệu [8] Москва: Энергоиздат, 1981. Sai phân Sai số Sai số [3] Phạm Trí Dũng, Nguyễn Thị Cẩm Hường, Phạm Văn Trinh, Đánh Vùng Thời Sun HH tuyệt đối tương đối giá hiệu quả của lỏng nitơ trên bệnh nhân sẹo lồi tại Khoa Da liễu khảo gian Bệnh viện An Giang, Kỷ yếu Hội nghị Khoa học Bệnh viện An sát [s] z abs(z-x) Giang, 2015. x [mm] % [mm] [mm] [4] http://123doc.org//document/2045706-tac-dung-cua-nito-long- 1s 0,75 0,81 0,06 8,0 trong-dieu-tri-benh-cham-man-tinh-tang-sung-ppt.htm [truy cập: -10°C 3s 1,7 1,52 0,18 10,6 01/10/2016] 7s 2,7 2,35 0,35 13,0 [5] eMedicine.com – Dermatology – Cryotherapy. [6] http://www.podiatrytoday.com/closer-look-cryosurgery-neuromas 1s 0,8 0,83 0,03 3,8 [Truy cập: 01/10/2016]. -5°C 3s 1,8 1,56 0,24 13,3 [7] http://cancer.com.vn/index.php/lieu-phap/lieu-phap-phau-thuat- 7s 2,9 2,51 0,39 13,4 lanh [Truy cập: 18/4/2016]. [8] F . Sun, L. Martínez-Suástegui, G.-X. Wang, K.M. Kelly, and G. Nhận xét Aguilar, “Numerical prediction of the intracellular ice formation Sai số kết quả tính toán giữa phương pháp số theo [8] zone during cryosurgery on a nodular basal cell carcinoma using và phương pháp sai phân hữu hạn đã nêu nhỏ hơn 15%. liquid nitrogen spray”, International Journal of Spray and Combustion Dynamics, Volume 4, Number 4, 2012, pp. 341-380. Trong một phạm vi nhất định đối với bài toán nhiệt, sai số [9] Nguyễn Bốn, Các phương pháp tính truyền nhiệt: Giáo trình cao trên có thể chấp nhận. học ngành Công nghệ Nhiệt – Điện lạnh, Trường Đại học Bách khoa – Đại học Đà Nẵng, 2007. 3. Kết luận và kiến nghị [10] Đặng Quốc Phú, Trần Thế Sơn, Trần Văn Phú, Truyền nhiệt, NXB 1. Tác giả đã giải thành công bài toán truyền nhiệt khi Giáo dục, Hà Nội, 2007. (BBT nhận bài: 12/7/2017, hoàn tất thủ tục phản biện: 24/8/2017)