Tổng quan luận án

Động lực học quá trình chuyển pha của hạt lưu chất, cụ thể là hiện tượng hóa rắn khi tiếp xúc bề mặt lạnh hoặc lơ lửng trong môi trường nhiệt độ thấp, là bài toán cơ học chất lỏng kết hợp truyền nhiệt có ý nghĩa thực tiễn trong kỹ thuật và tự nhiên. Luận án "Nghiên cứu tính toán động lực học hạt lưu chất đa thành phần có truyền nhiệt, chuyển pha" của nghiên cứu sinh Phạm Duy Bính tập trung giải quyết bài toán mô phỏng số đối với hạt lưu chất rỗng—cấu trúc gồm một lớp vỏ chất lỏng bao quanh một nhân bóng khí bên trong.

Tính cấp thiết của đề tài xuất phát từ hai yêu cầu ứng dụng kỹ thuật chính:

  • Chế tạo pin năng lượng mặt trời màng mỏng: Nghiên cứu của Biancardo và cộng sự [5] về hạt bán dẫn Silic dạng cầu đặt trong tấm lưới cho thấy chỉ lớp màng mỏng bên ngoài hạt đóng góp vào hiệu quả chuyển đổi quang điện (đạt hiệu suất 13,5% ở điều kiện chiếu sáng tiêu chuẩn $100\text{ mW/cm}^2$, AM1.5). Do đó, việc thay thế hạt đặc bằng hạt rỗng hóa rắn giúp tiết kiệm vật liệu bán dẫn và giảm khối lượng kết cấu mà vẫn duy trì hiệu suất.
  • Kiểm soát và phòng ngừa đóng băng trên bề mặt cánh máy bay [10] và tuabin gió [11], [12]: Hạt nước trong không khí bám vào cánh thường chứa bọt khí hòa tan hoặc bị giữ lại. Khi nhiệt độ bề mặt thấp hơn nhiệt độ đóng băng, chúng hóa rắn thành băng xốp hoặc băng rỗng, làm biến đổi biên dạng khí động học, giảm lực nâng, tăng lực cản, gây sụt giảm công suất hoặc dẫn đến tai nạn hàng không nghiêm trọng [10], [13], [14]. Ngoài ra, cấu trúc hạt đa thành phần còn liên quan đến công nghệ thực phẩm (giảm hàm lượng đường sucrose mà vẫn duy trì cảm quan độ ngọt theo Buyukkestelli và cộng sự [3]), bao bọc dược chất [1], [25]–[27], mỹ phẩm [22]–[24] và vật liệu tiêu âm [6].

Khoảng trống nghiên cứu được tác giả xác định là: các công trình lý thuyết, thực nghiệm và mô phỏng số trước đây chủ yếu tập trung vào hạt lưu chất đơn thành phần. Các nghiên cứu về hạt đa thành phần hầu như chỉ dừng lại ở bài toán biến dạng hoặc phân tách cơ học thuần túy trong dòng chảy cắt mà chưa khảo sát tương tác nhiệt - thủy động học đồng thời với quá trình hóa rắn có dịch chuyển của hai biên phân cách (biên ngoài và biên bóng khí trong).

Mục tiêu nghiên cứu của luận án gồm:

  1. Mục tiêu tổng quát: Nghiên cứu quá trình truyền nhiệt, chuyển pha (hóa rắn) của hạt lưu chất rỗng bằng phương pháp theo dấu biên (front-tracking method).
  2. Mục tiêu cụ thể:
    • 2.1. Xây dựng mô hình tính toán mô phỏng số bài toán hạt lưu chất rỗng hóa rắn trên nền tảng lưới Euler cố định và các điểm đánh dấu Lagrangian di động.
    • 2.2. Khảo sát định lượng ảnh hưởng của các thông số không thứ nguyên (Reynolds, Stefan, Bond, Prandtl, tỉ số khối lượng riêng, tỉ số bán kính) và thông số hình học (góc ướt, góc phát triển, độ lệch tâm) đến hình thái và thời gian hóa rắn của hạt.

Đối tượng và phạm vi nghiên cứu:

  • Đối tượng nghiên cứu: Các hạt lưu chất rỗng có kích thước hình học từ vài trăm micromet đến vài milimet, gồm lớp vỏ chất lỏng bao quanh một nhân khí bên trong.
  • Phạm vi nghiên cứu: Giới hạn trong mô hình lưu chất không nén được, không trộn lẫn, tuân theo định luật lưu chất Newton. Quá trình truyền nhiệt và động lực dòng chảy được mô tả qua hệ phương trình Navier-Stokes và phương trình năng lượng trong hệ tọa độ trụ hai chiều đối xứng trục $(O, r, z)$. Các đặc tính nhiệt vật lý gồm khối lượng riêng ($\rho$), độ nhớt động lực ($\mu$), hệ số dẫn nhiệt ($k$) và nhiệt dung riêng đẳng áp ($C_p$) trong từng pha được giả định là hằng số.

Tổng quan tài liệu và vị trí của luận án

Luận án tổng thuật có hệ thống các công trình nghiên cứu về hạt lưu chất theo ba nhánh chính: lý thuyết, thực nghiệm và mô phỏng số.

Về các nghiên cứu lý thuyết đối với hạt đơn:

  • Snoijer và Brunet [30] xây dựng mô hình giải tích cho hạt nước đông đặc trên mặt phẳng lạnh với giả định biên hóa rắn là đoạn thẳng phẳng, phân tích sự hình thành đỉnh chóp nhọn theo các tỉ số khối lượng riêng pha rắn/lỏng ($\rho_{sl}$) khác nhau và so sánh với giá trị tới hạn $\rho_c = 3/4$.
  • Zhang và cộng sự [31] khảo sát ảnh hưởng của trọng trường và nhiệt độ bề mặt, xác định chiều cao biên hóa rắn lý thuyết đạt 2,16 so với giá trị thực nghiệm là 2,00 (sai số thời gian hóa rắn trong khoảng 2%).
  • Tembely và Dolatabadi [32] đề xuất mô hình một chiều xem xét sự dãn nở nhiệt và thu được dạng biên hóa rắn cong lõm phù hợp với thực tế, đồng thời thiết lập phương trình dự đoán góc đỉnh hạt sau khi đông kết.
  • Anderson và cộng sự [33] phân tích bài toán hóa rắn với góc phát triển thay đổi dựa trên vận tốc trượt thẳng đứng ($V_v$), vận tốc trượt ngang ($V_s$) và góc tiếp xúc động ($\theta_a$) tại điểm chập ba pha.

Về các nghiên cứu thực nghiệm:

  • Hạt đơn: Huang và cộng sự [34] quan sát đỉnh chóp nhọn của hạt nước với góc ướt biến thiên trong dải $76,0^\circ - 154,9^\circ$. Pan và cộng sự [35] mở rộng thí nghiệm trên mặt phẳng nghiêng $30^\circ$ với góc ướt $77^\circ - 145^\circ$. Zhang và cộng sự [36] cùng Ju và cộng sự [37] khảo sát sự đông đặc trên bề mặt cong. Satunkin [38] đo đạc góc phát triển tinh thể của các vật liệu bán dẫn Silic (Si), Gecmani (Ge) và Indi Antimonua (InSb).
  • Hạt đa thành phần: Vũ và cộng sự [18] cùng Bhagat và cộng sự [19] thực nghiệm chế tạo hạt Gecmani rỗng qua vòi phun đồng trục, ghi nhận đỉnh chóp nhọn hình thành ở bề mặt ngoài nhưng mặt biên bên trong hạt bóng khí vẫn giữ độ nhẵn. Nadler và cộng sự [43] chế tạo hạt rỗng hợp kim nhôm (Al-2%Ca) cho hình cầu ổn định. Kendall và cộng sự [44] nghiên cứu dao động tắt dần của tia kim loại đa lớp để đảm bảo hạt đạt độ cầu đối xứng trước khi đông đặc.

Về các nghiên cứu mô phỏng số:

  • Hạt đơn chuyển pha: Shetabivash và cộng sự [45] sử dụng phương pháp tập mức kép (multiple level-set) phân chia biên dương và biên âm để mô phỏng đỉnh nhọn hạt nước. Vahab và cộng sự [46] kết hợp phương pháp tập mức và mô men lưu chất (moment of fluid) mô phỏng hạt nước va chạm cánh máy bay, chỉ ra rằng ở số Stefan lớn chóp nhọn không còn giữ được dạng hình học xác định. Sun và cộng sự [47] dùng phương pháp mạng Boltzmann (LBM) phân tích động học đông kết theo nhiệt độ vách. Zhou và cộng sự [48] dùng phương pháp thể tích lưu chất (VOF) khảo sát tia kim loại lỏng xâm nhập trong nước. Vũ và cộng sự [49], [50] dùng phương pháp theo dấu biên mô phỏng hạt đơn biến dạng dưới trọng trường khi chuyển pha.
  • Hạt đa thành phần (cơ học thuần túy): Hua và cộng sự [51] dùng phương pháp biên nhúng (immersed boundary) cho thấy góc nghiêng của nhân trong luôn lớn hơn hoặc bằng hạt ngoài trong dòng cắt. Vũ và cộng sự [52], [54] dùng phương pháp theo dấu biên thiết lập bản đồ chế độ tách hạt theo số Capillary ($Ca$) và Reynolds ($Re$) cho hạt một nhân và hai nhân. Chen và cộng sự [53] dùng VOF phân loại bốn dạng tách hạt. Vũ [55] và Phạm và cộng sự [56] khảo sát tương tác va chạm giữa hai hạt đa thành phần (chế độ trượt qua nhau hoặc đẩy ngược lại).

Vị trí của luận án: Luận án giải quyết khoảng trống nằm giữa hai mảng nghiên cứu nêu trên: phát triển mô hình tính toán số hoàn chỉnh có khả năng theo dõi đồng thời biên phân cách ngoài (lỏng-khí), biên phân cách trong (khí-lỏng) và biên chuyển pha (rắn-lỏng) kèm theo sự dịch chuyển của hai hệ điểm chập ba pha độc lập trong suốt quá trình truyền nhiệt và hóa rắn hạt rỗng.

Cơ sở lý thuyết và phương pháp nghiên cứu

Luận án xây dựng mô hình toán học dựa trên hệ phương trình Navier-Stokes một trường (one-field formulation) cho lưu chất không nén được kết hợp phương trình bảo toàn năng lượng có số hạng nguồn chuyển pha:

  1. Phương trình Navier-Stokes trong hệ tọa độ trụ $(r, z)$: $$\rho \left( \frac{\partial \mathbf{u}}{\partial t} + \mathbf{u} \cdot \nabla \mathbf{u} \right) = -\nabla p + \nabla \cdot \left[ \mu \left( \nabla \mathbf{u} + \nabla \mathbf{u}^T \right) \right] + \sigma \kappa \delta(\mathbf{x} - \mathbf{x}_f) \mathbf{n}_f + \rho \mathbf{g} + \mathbf{f}$$ Trong đó, $\mathbf{u} = (u, v)$ là véc tơ vận tốc, $p$ là áp suất, $\sigma$ là hệ số sức căng bề mặt, $\kappa$ là độ cong trung bình, $\mathbf{n}_f$ là véc tơ pháp tuyến đơn vị tại biên, $\delta(\mathbf{x} - \mathbf{x}_f)$ là hàm delta Dirac ba chiều định vị lực căng bề mặt tại biên phân cách $\mathbf{x}_f$, và $\mathbf{f} = \frac{\mathbf{u}_F - \hat{\mathbf{u}}}{\Delta t}$ là lực thể áp đặt điều kiện không trượt trong pha rắn ($\mathbf{u}_s = 0$).

  2. Phương trình năng lượng: $$\rho C_p \left( \frac{\partial T}{\partial t} + \mathbf{u} \cdot \nabla T \right) = \nabla \cdot (k \nabla T) + q , \delta(\mathbf{x} - \mathbf{x}_f)$$ Thông lượng nhiệt giải phóng tại biên chuyển pha tuân theo điều kiện Stefan: $$q = k_s \left( \frac{\partial T}{\partial \mathbf{n}} \right)_s - k_l \left( \frac{\partial T}{\partial \mathbf{n}} \right)_l$$

  3. Phương trình liên tục có tính đến biến thiên thể tích do chênh lệch khối lượng riêng giữa pha rắn và pha lỏng: $$\nabla \cdot \mathbf{u} = \left( \frac{1}{\rho_s} - \frac{1}{\rho_l} \right) \frac{q}{L_h} \delta(\mathbf{x} - \mathbf{x}_f)$$ với $L_h$ là nhiệt ẩn chuyển pha.

Hệ thống thông số không thứ nguyên được thiết lập thông qua các đại lượng tham chiếu ($L$, $U_t$, $t_t = L/U_t$, $\sigma_o$, $T_t$):

Đại lượng không thứ nguyên Định nghĩa toán học Ý nghĩa vật lý trong bài toán
Số Reynolds ($Re$) $\rho_l L U_t / \mu_l$ Tỉ số giữa lực quán tính và lực nhớt
Số Weber ($We$) $\rho_l L U_t^2 / \sigma_o$ Tỉ số giữa lực quán tính và lực căng bề mặt
Số Froude ($Fr$) $U_t / \sqrt{g L}$ Tỉ số giữa lực quán tính và lực trọng trường
Số Bond ($Bo$) $\rho_l g L^2 / \sigma_o$ Tỉ số giữa lực hấp dẫn và lực căng bề mặt
Số Prandtl ($Pr$) $\mu_l C_{pl} / k_l$ Tỉ số giữa khuếch tán động lượng và khuếch tán nhiệt
Số Stefan ($St$) $C_{pl} (T_m - T_c) / L_h$ Tỉ số giữa nhiệt cảm và nhiệt ẩn chuyển pha
Số Ohnesorge ($Oh$) $\sqrt{We} / Re = \mu_l / \sqrt{\rho_l \sigma_o L}$ Mối quan hệ giữa lực nhớt, lực quán tính và sức căng bề mặt
Tỉ số khối lượng riêng ($\rho_{sl}$) $\rho_s / \rho_l$ Mức độ dãn nở ($\rho_{sl} < 1$) hoặc co ngót ($\rho_{sl} > 1$) khi đông đặc
Tỉ số bán kính ($R_{io}$) $R_i / R_o$ Tỉ lệ kích thước hạt bóng khí trong so với hạt ngoài
Nhiệt độ không thứ nguyên ($\theta$) $(T - T_c) / (T_m - T_c)$ Thang nhiệt độ chuẩn hóa theo nhiệt độ bề mặt lạnh $T_c$ và nhiệt độ chuyển pha $T_m$

Phương pháp số thực tế:

  • Chương trình được lập trình bằng ngôn ngữ Fortran.
  • Phương pháp sai phân hữu hạn (finite difference) triển khai trên lưới so le MAC (Marker and Cell). Các biến vô hướng (áp suất, nhiệt độ, hàm chỉ thị) đặt tại tâm ô lưới, các thành phần vận tốc ($u, v$) đặt tại các mặt ô lưới. Phương trình áp suất Poisson được giải lặp để thỏa mãn phương trình liên tục tại mỗi bước thời gian.
  • Phương pháp theo dấu biên (Front-tracking): Các biên phân cách lỏng-khí, rắn-lỏng và khí-rắn được rời rạc hóa thành chuỗi các điểm Lagrangian liên kết di động trên nền lưới Euler cố định. Vị trí điểm biên được cập nhật theo $\mathbf{x}_f^{n+1} = \mathbf{x}_f^n + \mathbf{V}_f \Delta t$, trong đó vận tốc pháp tuyến của biên hóa rắn $\mathbf{V}_n$ được xác định trực tiếp từ điều kiện cân bằng nhiệt $q / (\rho_s L_h)$.
  • Hai hàm chỉ thị $I_s$ và $I_l$ được sử dụng để phân biệt trạng thái của 3 pha (rắn, lỏng, khí), từ đó gán chính xác các đặc tính vật lý cục bộ qua hàm bậc thang Heaviside nội suy từ vị trí các điểm mốc biên.

Nội dung chính theo từng chương

Chương 1: Tổng quan

Chương 1 hệ thống hóa các khái niệm nền tảng: định nghĩa hạt lưu chất, hạt đa thành phần, hạt rỗng (hạt đa thành phần có một nhân khí), cơ chế đối lưu cưỡng bức và hiện tượng chuyển pha hóa rắn trên bề mặt phẳng lạnh. Tác giả trình bày hai cơ chế chính hình thành hạt đa thành phần trong công nghiệp: (1) Sự tách hạt của tia lưu chất qua hệ vòi phun đồng trục thu hẹp - mở rộng theo nghiên cứu của Chen và cộng sự [16], Ji và cộng sự [17], Vũ và cộng sự [18], Bhagat và cộng sự [19]; (2) Sự co lại và đứt gãy của sợi lưu chất đa lớp theo các số Ohnesorge khác nhau theo mô phỏng của Hồ và cộng sự [20].

Chương này so sánh ưu và nhược điểm của 4 phương pháp số phổ biến trong bài toán dòng chảy đa pha:

  • Phương pháp tập mức (Level-Set): Bắt biên trơn tru nhưng khó bảo toàn khối lượng và tính toán lực căng bề mặt phức tạp.
  • Phương pháp thể tích lưu chất (VOF): Bảo toàn khối lượng tốt nhưng biên phân cách xấp xỉ dạng đoạn thẳng trong ô lưới, kém chính xác khi tính độ cong bề mặt.
  • Phương pháp mạng Boltzmann (LBM): Linh hoạt, dễ song song hóa nhưng gặp khó khăn khi mô tả biến thiên lớn về tỉ số nhiệt dung và hệ số dẫn nhiệt.
  • Phương pháp theo dấu biên (Front-Tracking): Biểu diễn biên bằng chuỗi điểm Lagrangian tường minh, cho phép xác định độ cong và tính lực căng bề mặt chính xác, dễ dàng theo dõi vị trí điểm chập ba pha. Đây là lý do tác giả chọn phương pháp theo dấu biên làm công cụ giải quyết bài toán của luận án.

Chương 2: Nghiên cứu và xây dựng chương trình tính toán cho bài toán có truyền nhiệt, chuyển pha

Chương 2 trình bày chi tiết thuật toán rời rạc hóa hệ phương trình vi phân chi phối trên lưới MAC so le, thiết lập điều kiện biên vật lý: biên đáy không trượt với nhiệt độ $T_c < T_m$, biên đối xứng trục tại $r = 0$, biên trên và biên phải áp dụng điều kiện trượt hoàn toàn.

Để kiểm chứng độ tin cậy của thuật toán, tác giả thực hiện tính toán đối chứng với dữ liệu thực nghiệm của Zhang và cộng sự [69] đối với hạt nước đơn ($R_o$) đông đặc trên bề mặt lạnh có góc ướt $\theta_0 = 85^\circ$:

  • Động học biên hóa rắn và sự thay đổi thể tích hạt theo thời gian mô phỏng số hoàn toàn trùng khớp với ảnh chụp thực nghiệm qua từng mốc thời gian không thứ nguyên ($\tau$).
  • Chương trình tiếp tục kiểm chứng sự phù hợp hình thái chóp nhọn đối với các trường hợp góc ướt nhỏ hơn $90^\circ$ ($\theta_0 = 78^\circ$) và góc ướt lớn hơn $90^\circ$ ($\theta_0 = 124^\circ$ và $\theta_0 = 155^\circ$). Trường nhiệt độ ($\theta$) và các véc tơ vận tốc không thứ nguyên hóa bằng vận tốc tham chiếu $U_t$ tái hiện chính xác hiện tượng dồn chất lỏng lên đỉnh hạt do dãn nở thể tích khi đông đặc.

Chương 3: Nghiên cứu quá trình truyền nhiệt, chuyển pha của hạt lưu chất rỗng trên một bề mặt lạnh

Chương 3 khảo sát bài toán hạt lưu chất rỗng đối xứng đặt trên bề mặt lạnh duy trì ở nhiệt độ $T_c$. Điểm chập ba pha xuất hiện tại cả biên ngoài và biên trong tiếp xúc với bề mặt đế.

Nghiên cứu độ hội tụ lưới được thực hiện qua hai độ phân giải $128 \times 256$ và $256 \times 512$ mắt lưới:

Thông số kiểm tra Sai số trung bình so với lưới $256 \times 512$
Chiều cao hạt ngoài ($H_d$) Dưới 1,2%
Chiều cao biên hóa rắn ($H_s$) Dưới 1,5%
Tọa độ trọng tâm bóng khí ($z_{c-in}$) Dưới 0,9%

Các kết quả khảo sát tham số chính trong Chương 3 gồm:

  • Số Bond ($Bo = 0,18$ và $Bo = 1,78$): Khi $Bo$ tăng, tác dụng của trọng trường lấn át sức căng bề mặt, hạt bị ép bẹp xuống, diện tích tiếp xúc với đế tăng, làm giảm nhẹ thời gian hóa rắn toàn phần ($t_s$) và làm giảm chiều cao hạt ngoài ($H_{out}$).
  • Số Stefan ($St = 0,032$ đến $0,316$): Tăng $St$ tương ứng với tăng độ chênh nhiệt $(T_m - T_c)$ hoặc giảm nhiệt ẩn $L_h$, vận tốc biên hóa rắn $V_n$ tăng mạnh, thời gian hóa rắn $t_s$ rút ngắn rõ rệt.
  • Tỉ số khối lượng riêng ($\rho_{sl} = 0,8$ và $1,1$): Với $\rho_{sl} = 0,8$ (nước dãn nở khi đông đặc), chất lỏng bị đẩy dồn lên trên tạo đỉnh chóp nhọn tại mặt ngoài và đẩy nâng bóng khí bên trong ($H_{in}$ tăng). Với $\rho_{sl} = 1,1$ (vật liệu co ngót thể tích), đỉnh hạt bị lõm xuống.
  • Tỉ số bán kính ($R_{io} = 0,3$ đến $0,7$): Hạt có nhân khí lớn ($R_{io} = 0,7$) có thể tích chất lỏng ít hơn, đường dẫn nhiệt qua màng lỏng mỏng hơn nên thời gian hóa rắn kết thúc nhanh hơn so với hạt nhân nhỏ ($R_{io} = 0,3$).
  • Góc ướt ngoài ($\theta_{0o} = 60^\circ - 130^\circ$) và góc ướt trong ($\theta_{0i} = 50^\circ - 120^\circ$): Góc ướt quyết định hình thái bán kính ướt ban đầu $R_w$, ảnh hưởng trực tiếp đến tốc độ truyền nhiệt từ bề mặt đế lạnh lên hạt.
  • Góc phát triển ($\theta_{gr} = 0^\circ - 25^\circ$): Góc phát triển kiểm soát độ nghiêng của mặt phân cách tại điểm chập ba pha, làm biến thiên độ dốc của thành hạt ngoài và biên bóng khí trong suốt quá trình kết tinh.
  • Hiện tượng điểm gián đoạn (điểm gấp khúc): Tác giả phát hiện khi biên hóa rắn tiến đến đỉnh của hạt bóng khí bên trong và khép kín xung quanh nhân khí, đường cong biến thiên vận tốc $V_n$ và chiều cao $H_d$ xuất hiện một điểm gấp khúc đặc trưng trước khi chuyển sang giai đoạn hóa rắn phần chỏm lỏng phía trên.

Chương 4: Nghiên cứu quá trình truyền nhiệt, chuyển pha của hạt lưu chất rỗng chịu ảnh hưởng của đối lưu cưỡng bức

Chương 4 mở rộng mô hình cho trường hợp hạt lưu chất rỗng lơ lửng trong không gian tự do, chịu tác động của một luồng khí lạnh chuyển động cưỡng bức từ dưới lên với vận tốc đầu vào $U_{in}$ và nhiệt độ $T_c$.

Kiểm tra độ độc lập của lưới được thực hiện trên 4 hệ lưới khác nhau so với lưới chuẩn $384 \times 1536$:

Độ phân giải lưới khảo sát Sai số trung bình so với lưới $384 \times 1536$
Lưới $96 \times 384$ 5,42%
Lưới $192 \times 768$ 1,85%
Lưới $288 \times 1152$ 0,61%

Các phát hiện nhiệt - thủy động lực học trong Chương 4 gồm:

  • Ảnh hưởng của số Reynolds ($Re$): Luồng khí cưỡng bức tạo ra sự bất đối xứng mạnh trong trường nhiệt độ giữa mặt đón gió (phía đáy) và mặt khuất gió (phía đỉnh). Số $Re$ tăng cường truyền nhiệt đối lưu ở phần đáy, làm biên hóa rắn bắt đầu từ cực dưới và tiến nhanh về cực trên, làm biến dạng tỉ số hình dạng trong ($A_{ri} = L_{is} / d_{is}$) và tỉ số hình dạng ngoài ($A_{ro} = L_{os} / d_{os}$), đồng thời làm giảm thời gian hóa rắn toàn phần $t_s$.
  • Ảnh hưởng của kích thước nhân hóa rắn ban đầu ($r_0 / R$): Tác giả áp đặt một mầm tinh thể lạnh tại cực đáy. Kích thước mầm ban đầu chi phối tốc độ giải phóng nhiệt ẩn và hướng phát triển ban đầu của mặt bao đóng băng.
  • Ảnh hưởng của độ lệch tâm ban đầu ($\varepsilon_0$): Khi tâm hạt bóng khí không trùng với tâm hạt ngoài ($\varepsilon_0 \neq 0$), gradient nhiệt độ phân bố bất đối xứng theo phương thẳng đứng. Nhân bóng khí lệch về phía đáy cản trở dòng nhiệt qua chất lỏng, làm thay đổi cục bộ vận tốc biên đông đặc và gây ra độ méo hình học đáng kể cho hạt rắn thành phẩm.
  • Ảnh hưởng của góc phát triển ($\theta_{gr}$) và tỉ số bán kính ($R_{io}$): Tác giả chỉ ra rằng dưới tác động của dòng bao cưỡng bức, sự kết hợp giữa biến dạng khí động học do lực cản áp suất và sự co dãn thể tích ($\rho_{sl}$) quyết định cấu trúc vỏ rắn rỗng có bị vỡ hay giữ được dạng hình học ổn định.

Kết quả và những đóng góp mới

Các đóng góp mới về mặt khoa học và phương pháp tính toán:

  • Đã xây dựng và hoàn thiện chương trình tính toán số viết bằng ngôn ngữ Fortran dựa trên phương pháp theo dấu biên (Front-tracking) kết hợp lưới so le MAC cho bài toán truyền nhiệt, chuyển pha của hệ lưu chất đa pha phức tạp (rắn - lỏng - khí). Mô hình giải quyết được sự tồn tại đồng thời của hai biên phân cách động và xử lý ổn định hai điểm chập ba pha (ở biên ngoài và biên trong) trong hệ tọa độ trụ đối xứng.
  • Đã thiết lập bộ số liệu định lượng mô tả chi tiết ảnh hưởng của các tổ hợp thông số không thứ nguyên ($Re, St, Bo, Pr, \rho_{sl}, R_{io}$) và các thông số tiếp xúc ($\theta_{0o}, \theta_{0i}, \theta_{gr}, \varepsilon_0$) tới hình thái biên hóa rắn, hình dạng đỉnh chóp nhọn, tỉ số hình dạng trong/ngoài ($A_{ri}, A_{ro}$) và thời gian hóa rắn hoàn toàn ($t_s$).
  • Phát hiện cơ chế hình thành điểm gấp khúc trên đường động học chiều cao hạt và vận tốc hóa rắn, tương ứng với thời điểm màng hóa rắn bao kín hoàn toàn nhân bóng khí bên trong.

Đóng góp về mặt thực tiễn và kiến nghị ứng dụng:

  • Cung cấp cơ sở dữ liệu và quy luật tương quan tham số phục vụ công nghệ chế tạo hạt bán dẫn rỗng trong pin mặt trời: Cho phép kiểm soát độ dày vỏ ngoài, bán kính nhân rỗng và độ đồng đều hình học thông qua việc điều chỉnh các số Stefan và Bond trong buồng làm nguội.
  • Đóng góp vào việc xây dựng mô hình dự báo đóng băng trên cánh tuabin gió và máy bay: Làm sáng tỏ cơ chế đông kết của các giọt nước có chứa bọt khí bám dính trên bề mặt hoặc va chạm trong luồng khí đối lưu, hỗ trợ thiết kế các giải pháp kiểm soát nhiệt và chống bám băng hiệu quả hơn.

Hạn chế và hướng nghiên cứu tiếp

Hạn chế của luận án:

  • Toàn bộ bài toán được giải quyết trong hệ tọa độ trụ giả định đối xứng trục hai chiều (2D axisymmetric), do đó chưa mô tả được các chế độ biến dạng không đối xứng ba chiều (3D) xảy ra khi có nhiễu loạn dòng chảy hoặc góc nghiêng bề mặt phức tạp.
  • Các tính chất nhiệt vật lý của lưu chất ($\rho, \mu, k, C_p$) trong từng pha được coi là hằng số độc lập với nhiệt độ; bỏ qua ảnh hưởng của bức xạ nhiệt và đối lưu tự nhiên sinh ra do chênh lệch mật độ nhiệt trong lòng chất lỏng (hiệu ứng Boussinesq).
  • Giới hạn trong mô hình lưu chất Newton thuần túy; dải thông số khảo sát tập trung ở một số giá trị đại diện do hạn chế về thời gian tính toán của phương pháp theo dấu biên.

Hướng nghiên cứu tiếp được đề xuất:

  • Mở rộng thuật toán theo dấu biên sang không gian ba chiều (3D) hoàn chỉnh để mô phỏng sự vỡ hạt, dao động bề mặt bất đối xứng hoặc bài toán dòng chảy xiên góc.
  • Nghiên cứu mô hình hóa rắn của hạt đa thành phần chứa nhiều nhân khí/lỏng (multi-core droplets) hoặc tương tác chuyển pha của cụm nhiều hạt lưu chất rỗng trong dòng khí.
  • Khảo sát các loại lưu chất phi Newton và tính toán sự phụ thuộc của các hệ số nhiệt vật lý theo nhiệt độ thực tế.

Giá trị tham khảo

Luận án là tài liệu tham khảo chuyên sâu đối với:

  • Nghiên cứu sinh, học viên cao học và giảng viên chuyên ngành Cơ học chất lỏng và chất khí, Thủy khí động lực học, Truyền nhiệt và Kỹ thuật tính toán: Cung cấp quy trình chi tiết về thuật toán rời rạc hóa phương trình Navier-Stokes một trường, kỹ thuật xử lý số hạng lực căng bề mặt CSF, phương pháp xấp xỉ biên di động Lagrangian trên lưới Euler và phương pháp lặp giải phương trình Poisson trên lưới MAC.
  • Kỹ sư và chuyên gia công nghệ vật liệu: Cung cấp phương pháp tính toán tối ưu hóa hình thái hạt rỗng phục vụ sản xuất vật liệu bán dẫn màng mỏng tiết kiệm nguyên liệu cho pin quang điện, sản xuất hạt rỗng kim loại hoặc hạt vi nang trong công nghệ thực phẩm và hóa mỹ phẩm.
  • Kỹ sư thiết kế trong ngành hàng không và năng lượng gió: Hỗ trợ dữ liệu phân tích quá trình hình thành lớp băng xốp/băng rỗng trên biên dạng cánh tuabin và cánh máy bay khi vận hành trong điều kiện sương giá khí quyển.

Phần đáng tham khảo nhất trong luận án là Chương 2 (chi tiết thiết lập thuật toán theo dấu biên cho bài toán 3 pha) và Chương 3 (phân tích định lượng sự phụ thuộc của hình thái hạt hóa rắn vào số Bond, Stefan và tỉ số khối lượng riêng).

Câu hỏi thường gặp

1. Tại sao luận án lựa chọn phương pháp theo dấu biên (Front-tracking) thay vì phương pháp VOF hay Level-set?
Phương pháp theo dấu biên sử dụng chuỗi các điểm Lagrangian liên kết di chuyển tự do trên nền lưới Euler cố định để mô tả mặt phân cách. Cách tiếp cận này giúp xác định tường minh vị trí biên, tính toán trực tiếp và chính xác độ cong bề mặt $\kappa$ cùng lực căng bề mặt $\sigma$, khắc phục nhược điểm khuếch tán số và khó tính lực mao dẫn của phương pháp thể tích lưu chất (VOF) và phương pháp tập mức (Level-set), đặc biệt khi bài toán có sự xuất hiện đồng thời của nhiều bề mặt phân cách lỏng-khí, rắn-lỏng và khí-rắn.

2. Điểm chập ba pha của hạt lưu chất rỗng trên bề mặt phẳng lạnh khác biệt như thế nào so với hạt lưu chất đơn?
Ở hạt lưu chất đơn đặt trên bề mặt phẳng, chỉ tồn tại một hệ điểm chập ba pha lỏng - khí - rắn tại biên ngoài tiếp xúc với đế. Đối với hạt lưu chất rỗng, xuất hiện thêm một hệ điểm chập ba pha thứ hai nằm ở biên bên trong, nơi mặt phân cách khí - lỏng của hạt bóng khí tiếp xúc với bề mặt đế lạnh. Do đó, chương trình tính toán phải theo dõi và xử lý đồng thời điều kiện góc tiếp xúc động và góc phát triển tại cả hai vị trí tiếp xúc độc lập này.

3. Tỉ số khối lượng riêng giữa pha rắn và pha lỏng ($\rho_{sl}$) chi phối hình thái hạt sau khi đông kết như thế nào?
Tỉ số $\rho_{sl}$ quyết định mức độ biến thiên thể tích khi chuyển pha:

  • Khi $\rho_{sl} < 1$ (ví dụ như nước, thể tích dãn nở khi đóng băng), phần chất lỏng chưa đông bị nén và đẩy dồn lên phía trên, dẫn đến việc hình thành đỉnh chóp nhọn ở bề mặt ngoài và làm tăng chiều cao hạt bóng khí bên trong ($H_{in}$).
  • Khi $\rho_{sl} > 1$ (vật liệu co ngót thể tích khi đông kết), đỉnh hạt ngoài bị kéo lõm xuống, không tạo thành chóp nhọn.
  • Khi $\rho_{sl} = 1$, hình thái hình học của hạt hầu như không biến đổi sau khi kết thúc quá trình hóa rắn.

4. Hiện tượng "điểm gấp khúc" trong Chương 3 xuất hiện vào thời điểm nào của quá trình hóa rắn?
Điểm gấp khúc xuất hiện trên đồ thị động học chiều cao hạt ($H_d$) và vận tốc biên hóa rắn ($V_n$) tại thời điểm mặt bao hóa rắn phát triển từ đáy tiến lên và khép kín hoàn toàn xung quanh đỉnh của hạt bóng khí bên trong. Sau thời điểm này, toàn bộ hạt bóng khí đã bị cô lập hoàn toàn trong pha rắn, quá trình truyền nhiệt chuyển pha chuyển sang giai đoạn hóa rắn lớp chất lỏng còn lại ở đỉnh hạt ngoài với tốc độ biến thiên diện tích biên truyền nhiệt thay đổi đột ngột.

5. Luồng khí đối lưu cưỡng bức trong Chương 4 làm thay đổi quy luật hóa rắn của hạt rỗng lơ lửng như thế nào so với hạt đặt trên đế lạnh?
Đối với hạt đặt trên bề mặt phẳng lạnh (Chương 3), nhiệt được truyền dẫn chủ yếu qua bề mặt đáy tiếp xúc phẳng, biên hóa rắn tiến đều từ đáy lên đỉnh theo phương thẳng đứng. Đối với hạt lơ lửng chịu dòng cưỡng bức (Chương 4), sự phân bố áp suất và lớp biên nhiệt xung quanh hạt bị chi phối bởi số Reynolds ($Re$). Truyền nhiệt đối lưu diễn ra mạnh nhất tại điểm dừng phía đáy (mặt đón gió) và yếu hơn ở vùng xoáy phía sau (mặt khuất gió), làm cho vỏ rắn hình thành từ cực dưới bao bọc dần lên cực trên và gây biến dạng hình học không đồng đều giữa phần dưới và phần trên của hạt.

Kết luận

Luận án của nghiên cứu sinh Phạm Duy Bính đã thiết lập thành công mô hình tính toán số dựa trên phương pháp theo dấu biên để giải quyết trọn vẹn bài toán động lực học truyền nhiệt và hóa rắn hạt lưu chất rỗng trong hai cấu hình: đặt trên bề mặt lạnh và lơ lửng dưới tác động của dòng đối lưu cưỡng bức. Kết quả mô phỏng đã làm rõ cơ chế tương tác phức tạp giữa các lực mao dẫn, lực trọng trường, lực quán tính và biến thiên thể tích do chuyển pha thông qua hệ thông số không thứ nguyên xác định. Các phát hiện về quy luật hình thái chóp nhọn, tỉ số hình dạng và thời gian hóa rắn cung cấp luận cứ khoa học có giá trị trực tiếp cho công nghệ chế tạo hạt rỗng bán dẫn màng mỏng và mô hình hóa kiểm soát đóng băng trong kỹ thuật hàng không và năng lượng.