Ứng dụng phương pháp sai phân hữu hạn để tính toán quá trình nung cảm ứng đạt trạng thái bán lỏng của phôi hợp kim nhôm trong công nghệ tạo hình xúc biến

Trong bài báo này, phương pháp sai phân hữu hạn (FDM) được sử dụng để tính toán quá trình truyền nhiệt khi gia nhiệt bằng dòng cảm ứng nhằm đưa nhiệt độ tới vùng bán lỏng...

Application of finite differencial method in heat transfer for induction heating to semisolid state of aluminium alloy in thixoforming process

 NGUYỄN VINH DỰ ¹, LƯU PHƯƠNG MINH^{2 1} Sở Khoa học và Công nghệ TP HCM ² Trường Đại học Bách khoa TP HCM

 TÓM TẮT

Trong bài báo này, phương pháp sai phân hữu hạn (FDM) được sử dụng để tính toán quá trình truyền nhiệt khi gia nhiệt bằng dòng cảm ứng nhằm đưa nhiệt độ tới vùng bán lỏng. Đồng thời phân tích và so sánh kết quả tính toán của phương pháp số với phương pháp Comsol và phương pháp thực nghiệm. Kết quả so sánh cho thấy phương pháp sai phân hữu hạn chính xác cho quá trình truyền nhiệt với mô hình dạng trụ (đường kính D = 76 mm, cao H = 50 mm) với điều kiện biên đặc biệt vừa nung vừa làm nguội bề mặt. Bằng phương pháp sai phân hữu hạn, kết quả tính toán cho thấy việc kết hợp các phương pháp mô phỏng số, mô phỏng bằng phần mềm Comsol và thực nghiệm đã hỗ trợ đắc lực cho việc đưa ra chiến lược đồng thời vừa gia nhiệt, vừa làm nguội phôi bằng thiết bị nung cảm ứng với kết quả tốt hơn: thời gian nung được rút ngắn 50 % so với các nghiên cứu cùng mô hình vật lý và độ chênh nhiệt độ là ≤ 5 ^oC. Từ khóa: Phương pháp sai phân hữu hạn, tạo hình xúc biến, nung cảm ứng

 ABSTRACT

In this paper, the finite differential method (FDM) is applied to calculate the heat transfer process during induc- tion heating to temperature of the semisolid state. At the same time, the results of the numeric calculation method will be analyzed and compared with that of the Comsol method and experiments. Results of comparison show that the finite differential method is accurate for the heat transfer process with cylindrical form model (diameter D = 76 mm, height H = 50 mm) and particular boundary conditions - simultaneous heating and surface cooling. By FDM, calculation results show that the combination of the numeric simulation method with the Comsol simulation and the experimentation has supported the strategy of simultaneous heating and surface cooling. The used induction fur- nace gives better results with the heating time of only 50 % compared with the research of the same physical model and the difference in temperature is less than 5 ^oC. Keywords: Finite differential method (FDM), thixoforming, induction heating  MỞ ĐẦU Trong công nghệ tạo hình vật liệu ở trạng thái bán lỏng, một trong những vấn đề được các chuyên gia quan tâm trong quá trình nghiên cứu là xác định quá trình gia nhiệt đảm bảo độ đồng đều nhất về nhiệt độ trên toàn bộ thể tích của phôi. Một trong những phương án công nghệ để đạt yêu cầu đó là phương pháp nung không tiếp xúc bằng dòng cảm ứng, không chỉ dễ điều khiển mà còn tăng tỷ lệ sử dụng vật liệu, thời gian nung nhanh [1]. Vì vậy, bài toán truyền nhiệt trong công nghệ nung cảm ứng sẽ có vai trò rất quan trọng trong việc điều khiển quá trình nung này. Ở đây, đã lựa chọn Phương pháp sai phân hữu hạn (Finite Differential Method) 2 chiều để tiến hành tính toán bài toán nhiệt. Dựa trên cơ sở của các tác giả đi trước như M. Necati Ozisik [2] cho bài toán tổng quát, tác giả Ewa Majchrzak [3] - Viện Khoa học máy tính thuộc Đại học kỹ thuật Czestochowa, Ba Lan sử dụng phương pháp FDM một chiều 1D, tác giả Ozisik [4] với phương pháp một chiều trong hệ tọa độ cực; bài toán truyền nhiệt một biến của tác giả Subhash [5], bài toán truyền nhiệt hai chiều của các tác giả Ching-yu Yang [6]. Trong nghiên cứu này đã sử dụng mô hình vật lý và điều kiện biên có khác biệt lớn với việc sử dụng một nguồn nhiệt,  cụ thể là mô hình vật lý trong nghiên cứu này có xét đến điều kiện mất nhiệt do bức xạ và đối lưu tại bề mặt trên của chi tiết, đồng thời so với các mô hình vật lý nghiên cứu của các tác giả khác thì trong mô hình này có xét quá trình làm mát bề mặt trong quá trình nung, do đó điều kiện biên sẽ khác và đây là yếu tố làm cho bài toán phức tạp, đòi hỏi mô hình tính toán khó khăn và lớn hơn rất nhiều - ma trận tính toán. 2. MÔ HÌNH TRUYỀN NHIỆT VỚI BÀI TOÁN DẪN NHIỆT KHÔNG ỔN ĐỊNH HAI CHIỀU a. Mô hình vật lý của quá trình nung cảm ứng Hình 1. Mô hình vật lý sử dụng cho quá trình nghiên cứu [9] Phôi nung là hợp kim nhôm A356, hình trụ có đường kính 76 mm, chiều cao 50 mm, đặt theo phương thẳng đứng (Hình 1), nguồn nhiệt bên trong được tạo ra nhờ sự cảm ứng điện từ để đảm bảo được tỷ lệ pha rắn-pha lỏng là (45 - 50) % thì phôi được gia nhiệt đến 585 ^oC. b. Mô hình toán học Phương trình truyền nhiệt: ở đây: S(r, t) là nguồn nhiệt bên trong, ρ   là trọng lượng riêng, Cp là nhiệt dung riêng, K là độ dẫn nhiệt. Khi nhiệt độ T trong khoảng: T_S,<,T,<,TL (với TS: nhiệt độ đường đặc, TL: nhiệt độ đường lỏng) thì xuất hiện ẩn nhiệt trong kim loại nung, khi đó Cp(T) được tính theo công thức sau: với L là ẩn nhiệt của hợp kim nhôm A356. Giả sử nguồn nhiệt được biểu diễn như sau: S(r, t) = G(t) . H(r). Từ phương trình Maxwell, ta có:  ; ở đây: I0 là dòng điện xoay chiều cực đại trên bề mặt. Từ định luật Joule, có thể suy ra sự phân bố của nguồn nhiệt:   trong đó: So là nguồn nhiệt cực đại trên bề mặt tương ứng khi r = L, khi đó S_o bằng công suất cực đại của thiết bị [7]; là chiều sâu thẩm thấu. Đối với vật liệu kim loại, độ thẩm thấu của vật liệu phụ thuộc tần số f của dòng điện [7]:  Đối với hợp kim nhôm và đồng [7] thì µ_r ≅ 1, độ dẫn nhiệt χ  tỷ lệ thuận với nhiệt độ; đối với hợp kim nhôm A356 thì  .  - Điều kiện biên:  Hình 2. Tổng hợp các yếu tố nhiệt, điện từ xảy ra trong lò nung cảm ứng [7] * Tại r = 0:  (5) * Tại r = R:  (6) →  (7) → q_{làm mát}  (t): Trong mô hình vật lý này thì vật nung sẽ được làm mát bằng không khí, công suất làm mát sẽ thay đổi theo thời gian t và phụ thuộc độ chênh nhiệt độ giữa bề mặt và tâm phôi, [W/m²]. với: - Điều kiện ban đầu: * Tại z = 0:  (9)  * Tại z = H:  (10) c. Chia lưới Thời gian t được chia nhỏ thành; Δt = 21(s) từ 0 → t^f bán kính R được chia nhỏ thành; chiều cao H được chia nhỏ thành Thiết lập hệ phương trình đại số tại các điểm nút như sau:

• Nút lưới bên trong: i = 1 ÷ (m - 1); j = (l - 1)

Đây là bài toán không ổn định 2 chiều, việc giải bài toán này phức tạp hơn so với bài toán không ổn định 1 chiều do thêm chiều z và số lượng điểm nút sẽ lớn hơn rất nhiều. Đặt:  

• Nút lưới trên đường tâm đối xứng: i = 0; j = 1 ÷  (l - 1)

 (12)

• Nút lưới trên biên r = R: i = m; j = 1 ÷ (l - 1)  Với  

 (13)

• Nút lưới trên biên z = 0: i = 1 ÷ (m - 1); j = 0

• Tương tự với nút lưới trên biên z = H: i = 1 ÷ (m - 1);  j = 1

 Với  thu được:   (15) Tương tự nút lưới tại r = 0; z = 0: (i = 0; j = 0)  và   (16)

•  Nút lưới tại r = 0, z = H: (i = 0; j = 1)

• Nút lưới tại r = R, z = 0: (i = m; j = 0)

Thu được: 

• Nút lưới tại r = R, z = H: (i = m; j = 1)

Hệ phương trình của bài toán dẫn nhiệt không ổn định 2 chiều tại tất cả các điểm nút là tập hợp các phương trình (11); (12); (13); (14); (15); (16); (17); (18); (19) Từ đó hệ phương trình có dạng:  MTⁿ⁺¹ = Tⁿ + D + E₁ + E₂        (20)  3. KẾT QUẢ VÀ THẢO LUẬN Để giải bài toán này, đã sử dụng ba phương pháp: phương pháp số, phương pháp Comsol, phương pháp thực nghiệm, đồng thời so sánh đánh giá sai số của các kết quả trên. a. Phương pháp số Đã tận dụng khả năng tính toán của máy tính để hỗ trợ tính toán nhanh chóng và hiệu quả bằng cách đưa hệ (20) về dạng ma trận. b. Phương pháp Comsol Comsol là một trong những phần mềm nổi tiếng trong lĩnh vực mô phỏng đa vật lý. Với khả năng liên kết các dạng bài toán mô phỏng khác nhau, Comsol có khả năng ứng dụng trong các lĩnh vực như: điện-từ trường, cơ khí, nhiệt, lưu chất và hóa học. Kết quả so sánh hai phương pháp a và b cho trong Bảng 1. Bảng 1. Kết quả phương pháp số và mô phỏng Comsol theo tọa độ z [8, 9] [caption id="attachment_1944" align="aligncenter" width="640"] Hình 3. Sơ đồ thiết bị tiến hành thí nghiệm [10][/caption] Bảng 2. Thành phần hóa học (%) của A356 [9]

Al	Si	Mg	Cu	Fe	Mn	Ti	Cr	Fb	Zn
92	7,270	0,403	0,007	0,149	0,007	0,011	<0,002	<0,005	0,049

 Bảng 3. Tính chất vật lý của hợp kim nhôm A356 [2, 7] Nhiệt dung riêng của vật liệu (J/kg ^oC)       0,454.T (C) + 904,6 Độ dẫn nhiệt của vật liệu (W/m ^oC)             0,04.T (C) + 153,1 Khối lượng riêng của chi tiết (kg/m³)           -0,208T (C) + 2680 c. Phương pháp thực nghiệm Mô hình thực nghiệm dùng để kiểm tra chiến lược nung và làm nguội, đồng thời so sánh kết quả tính toán giữa phương pháp số, mô phỏng bằng Comsol. Nhiệt độ trên phôi được đo bằng các cặp nhiệt điện cắm vào phôi tại các vị trí khoan lỗ định trước, trong đó một lỗ khoan ở sát bề mặt phôi và một lỗi khoan ngay tâm phôi. Cảm biến nhiệt sử dụng là cảm biến K, một cực có vật liệu là Ni-Cr (Chromel), cực còn lại là Ni-Al, giới hạn nhiệt độ sử dụng là 1.150 ^oC (Hình 3). Phôi thực nghiệm: Phôi ban đầu cần phải gia nhiệt lại là phôi trụ có kích thước F76x50, làm từ hợp kim nhôm A356 có yêu cầu độ chênh nhiệt độ trong cả phôi là DT < 10 ^oC . Thành phần và tính chất vật lý của A356 nêu trong các Bảng 2 và 3. d. Thảo luận Nhiệt độ tại các điểm A và B nhận được từ các phương pháp trên cho trong Bảng 4 và các Hình 4, 5, 6 và 7. [caption id="attachment_1943" align="aligncenter" width="300"] Hình 4. Nhiệt độ tại điểm A[/caption] [caption id="attachment_1942" align="aligncenter" width="300"] Hình 5. Nhiệt độ tại điểm A[/caption] [caption id="attachment_1941" align="aligncenter" width="300"] Hình 6. Nhiệt độ tại điểm A và điểm B thể hiện bằng phương pháp thực nghiệm[/caption] [caption id="attachment_1940" align="aligncenter" width="300"] Hình 7. Nhiệt độ tại điểm B[/caption] Bảng 4. Nhiệt độ tại hai điểm A và B [9] Các đồ thị trên các Hình 8, 9 và 10 thể hiện các giá trị chênh lệch từ kết quả đo của bảng 4. [caption id="attachment_1938" align="aligncenter" width="300"] Hình 8. Chênh lệch nhiệt độ tại điểm A[/caption] [caption id="attachment_1937" align="aligncenter" width="300"] Hình 9. Chênh lệch nhiệt độ tại điểm B[/caption] [caption id="attachment_1936" align="aligncenter" width="300"] Hình 10. Chênh lệch nhiệt độ tại hai điểm A và B[/caption] [caption id="attachment_1935" align="aligncenter" width="300"] Hình 11. Nhiệt độ tại hai điểm A và B[/caption] Nhận xét: [caption id="attachment_1934" align="aligncenter" width="300"] Hình 12. Biểu đồ nung và làm nguội tại tần số f = 21kHz[/caption] Hình 11 cho biết hệ thống kiểm soát nhiệt độ và công suất tại các điểm A và B, Hình 12 cho thấy biến đổi công suất nung của thiết bị qua từng thời điểm khác nhau. Trong giai đoạn ban đầu khoảng 30 giây, thiết bị nung nung với công suất được gia tăng dần dần và được duy trì khoảng 300 giây với công suất tối đa và dần dần giảm khi gần đến cuối thời gian nung. Từ các đồ thị hình 12, nhận thấy quá trình nung hợp kim nhôm A356 đến trạng thái bán lỏng diễn ra theo hai giai đoạn chính như sau: a. Quá trình gia nhiệt: Quá trình này diễn ra trong khoảng thời gian từ 0 đến 30 giây. Lúc này, phôi nung được nâng nhiệt theo đường cong; đồng thời với quá trình này ta cũng tiến hành làm nguội bề mặt chi tiết theo đường thẳng. Do hiệu ứng bề mặt nên nhiệt độ ở tâm sẽ cao hơn ở bề mặt, tuy vậy nhiệt độ trong toàn phôi nung vẫn tăng đều. Và trong khoảng thời gian này vẫn phải duy trì quá trình làm nguội, bởi vì độ chênh nhiệt độ trong phôi nung trong giai đoạn đang tăng. b. Quá trình giữ nhiệt: quá trình này bắt đầu diễn ra từ 30 đến 300 giây phôi nung vẫn được gia nhiệt đồng thời kết hợp việc tăng cường độ làm nguội. Từ kết quả thực nghiệm ta thấy, tại điểm A, B sai số, mặc dù có sự chênh lệch sai số giữa thời gian đầu và thời gian cuối khi nung từ ≈1 % còn 0,03 % ở thời gian cuối nhưng bản chất chênh lệch nhiệt độ về cơ bản gần như nhau (do nhiệt độ này về cuối càng cao nên việc chênh lệch nhiệt độ ảnh hưởng đến kết quả sai số) nguyên nhân do chênh lệch nhiệt độ có thể là do sai số của thiết bị, độ không đồng tâm giữa phôi và cuộn dây, ngoài ra chất lượng của mức độ đồng đều phôi ban đầu cũng ảnh hưởng đến quá trình truyền nhiệt trong khi nung. Nhưng tổng thể chung, độ chính xác giữa 3 phương pháp: phương pháp số, phương pháp mô phỏng Comsol, phương pháp thực nghiệm là gần như tiệm cận. Qua đó, kết quả tính toán số và thực nghiệm gần như nhau và đảm bảo độ chênh nhiệt độ giữa tâm và bên ngoài hay nói cách khác là độ đồng đều trên toàn bộ thể tích phôi là ≤ 10 ^oC đúng theo điều kiện ban đầu đặt ra.  4. KẾT LUẬN Đã nghiên cứu bài toán dẫn nhiệt không ổn định 2 chiều, trong đó xem xét sự thay đổi nhiệt độ bên trong vật thể hình trụ tròn theo cả 2 phương bán kính R và chiều cao H. Qua quá trình tính toán cho thấy khối lượng tính toán lớn; lập hệ phương trình, đưa về ma trận tương đối phức tạp và để giải quyết bài toán này phải có sự hỗ trợ của máy tính nên cần nhiều tài nguyên bộ nhớ máy tính để chạy chương trình và điểm đặc biệt là các ma trận chỉ cần lập cho trường hợp tổng quát một lần và áp dụng tương tự cho các bài toán cùng dạng. TÀI LIỆU TRÍCH DẪN

1. V. T. Borukhov, G. M. Zayats, Identification of a time-dependent source term in nonlinear hyperbolic or par- abolic heat equation.
2. M. Necati Ozisik, Finite Difference Methods in Heat Transfer, CRC, 1994.
3. Ewa Majchrzak, Application of FDM for numberical solution of Hyperbolic heat conduction equation, Technical Report.
4. M. Necati Ozisik, On the wave theory in heat conduction, Journal Heat Transfer, 1994.
5. Subhash C. Mishra*, N. Gullipalli, A. Jain, A. Barurah, A. Mukherjee, Analysis of hyperbolic heat conduction in 1-D planar, cylindrical, and spherical geometry using the lattice Boltzmann method, 2016.
6. 6. Ching-yu Yang, Direct and inverse solutions of the two-dimensional hyperbolic heat conduction problems,
7. A. Bendada, C. Q. Zheng and N. Nardini, Investigation of temperature control parameters for inductively heat- ed semi-solid light alloys using infrared imaging and inverse heat conduction, Journal of physics D: Applied physics, Institute of physics publishing, 2004.
8. N. V. Dự, T. V. Quốc, L. P. Minh, Ứng dụng phần mềm Comsol mô phỏng quá trình nung cảm ứng đến trạng thái bán lỏng của phôi hệ hợp kim nhôm trong công nghệ tạo hình Thixoforming, Tạp chí Cơ khí Việt Nam, 2015.
9. V. D. Nguyen, P. M. Luu and T. Q. Quang, “Experimantal and Numerical Studies on the Effects of Heating Frequency in the Thixoforming Process for the 2d Aluminium Alloy Semi-solid State”, IEEE International Conference on Industrial Engineering and Engineering Management, 2016.
10. N. V. Dự, P. Q. An, P. S. Minh, N. N. Hà, Research, designed and Manufactured the equipment heating induction to aluminium alloys semisolid state in Thixoforming Process, Proceedings of the 4th National Conference on Mechanical Science and Technology, 2015.