Tổng quan về luận án

Hiện tượng động lực học của hạt lưu chất đa thành phần đi kèm quá trình truyền nhiệt và chuyển pha hóa rắn là bài toán cơ học chất lỏng phức tạp, đóng vai trò nền tảng trong cả hiện tượng tự nhiên lẫn các ngành công nghiệp tiên tiến. Nghiên cứu của Phạm Duy Bính (2023) với đề tài "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" thuộc chuyên ngành Cơ học chất lỏng và chất khí (Mã số: 9.08) tại Viện Hàn lâm Khoa học và Công nghệ Việt Nam đã thiết lập một bước tiến khoa học quan trọng. Công trình giải quyết bài toán mô phỏng số chi tiết tương tác đa pha đa vật lý của hạt lưu chất rỗng (hạt lỏng chứa nhân khí bên trong) trong hai điều kiện biên quan trọng: hóa rắn trên bề mặt phẳng lạnh đẳng nhiệt và hóa rắn lơ lửng trong môi trường tự do dưới tác động của dòng đối lưu cưỡng bức.

Khoảng trống nghiên cứu (research gap) lớn nhất được nhận diện là sự thiếu hụt các mô hình lý thuyết và mô phỏng số giải quyết được bài toán chuyển pha của hạt lưu chất đa thành phần. Trong khi các công trình tiền nhiệm chủ yếu tập trung vào hạt đơn pha (Snoijer & Brunet, 2012; Huang et al., 2012; Shetabivash et al., 2014) hoặc hạt đa thành phần đẳng nhiệt không có chuyển pha (Hua et al., 2014; Vũ et al., 2017), động học mặt phân cách kép (lỏng-khí ngoài, rắn-lỏng chuyển pha, và lỏng-khí bên trong nhân bóng khí) vẫn chưa được làm sáng tỏ. Luận án đặt ra ba câu hỏi nghiên cứu cốt lõi:

  1. Làm thế nào để mô hình hóa chính xác các biên phân cách pha chuyển động đồng thời và xử lý chính xác sự đột biến thuộc tính vật lý cùng bước nhảy thể tích khi khối lượng riêng thay đổi ($\rho_s \neq \rho_l$)?
  2. Quy luật tác động của các thông số không thứ nguyên ($Bo, Pr, St, \rho_{sl}, R_{io}$, góc ướt trong $\theta_{0i}$, góc ướt ngoài $\theta_{0o}$, góc phát triển $\theta_{gr}$) lên hình thái chóp nhọn (cusp), biên phân cách trong-ngoài và thời gian hóa rắn của hạt rỗng trên đế lạnh là gì?
  3. Sự tương tác giữa trường vận tốc đối lưu cưỡng bức bên ngoài ($Re$), độ lệch tâm nhân khí ($\varepsilon_0$) và kích thước mầm kết tinh ($r_0/R$) làm biến dạng cấu trúc pha rắn và tỷ số hình dạng ($Ar_i, Ar_o$) của hạt lơ lửng như thế nào?

Luận án dựa trên khung lý thuyết cơ học môi trường liên tục, bao gồm hệ phương trình Navier-Stokes một trường (one-field formulation), phương trình bảo toàn năng lượng có số hạng nhiệt ẩn chuyển pha và phương trình liên tục mở rộng cho lưu chất chịu biến đổi thể tích cục bộ. Đóng góp đột phá của luận án nằm ở việc phát triển chương trình tính toán số độc lập dựa trên phương pháp theo dấu biên (Front-Tracking Method) trên hệ tọa độ trụ $(Orz)$, cho phép kiểm soát trực tiếp lực căng bề mặt thông qua hàm delta Dirac và bảo toàn hình học mặt phân cách. Phạm vi khảo sát tập trung vào các hạt lưu chất có kích thước vi mô đến mi-li-mét ($10^{-4} \text{ m} - 10^{-3} \text{ m}$) trong dải thông số không thứ nguyên rộng: $Bo \in [0,18; 1,78]$, $Pr \in [0,032; 0,316]$, $St \in [0,032; 0,316]$, $\rho_{sl} \in [0,8; 1,1]$, $R_{io} \in [0,3; 0,7]$, góc tiếp xúc $\theta \in [50^\circ; 130^\circ]$, và góc phát triển $\theta_{gr} \in [0^\circ; 25^\circ]$.

Literature Review và Positioning

Tổng quan y văn về động lực học chuyển pha hạt lưu chất được chia thành ba nhánh chính: tiếp cận giải tích - giải tích nửa thực nghiệm, tiếp cận thực nghiệm quang học, và tiếp cận mô phỏng số. Về mặt lý thuyết, Snoijer và Brunet (2012) đặt nền móng với mô hình giải tích dự đoán hình thái hạt nước đóng băng trên bề mặt, phát hiện sự hình thành chóp nhọn phụ thuộc vào tỷ số khối lượng riêng $\rho_{sl} = \rho_s/\rho_l$ so với giá trị tới hạn $\rho_c = 3/4 = 0,75$. Zhang et al. (2017) mở rộng lý thuyết cho thấy chiều cao biên hóa rắn theo tính toán giải tích đạt $2,16$ so với thực nghiệm là $2,0$ (sai số xấp xỉ $2%$), đồng thời chứng minh góc giữa bề mặt tự do và phương ngang giảm dần trong quá trình đóng băng trước khi tăng đột ngột ở giai đoạn cuối tạo thành đỉnh nhọn. Anderson et al. (2010) bổ sung động học điểm chập ba pha (tri-junction) thông qua mối quan hệ giữa vận tốc trượt thẳng đứng $V_v$, vận tốc trượt nằm ngang $V_s$ và góc phát triển tiếp xúc $\phi_i$. Tembely và Dolatabadi (2019) cải tiến mô hình 1D giải quyết độ cong thực tế của biên hóa rắn lõm.

Ở hướng tiếp cận thực nghiệm, Huang et al. (2012) khảo sát quang học hạt nước đơn lẻ với góc ướt biến thiên từ $76,0^\circ$ đến $154,9^\circ$, khẳng định sự tồn tại phổ quát của hiện tượng tạo đỉnh nhọn. Pan et al. (2018) mở rộng cho bề mặt nghiêng $30^\circ$ với góc tiếp xúc $77^\circ - 145^\circ$. Đối với hạt đa thành phần, các công trình thực nghiệm của Bhagat et al. (2013) chế tạo hạt Gecmani ($Ge$) rỗng bằng vòi phun đồng trục và Nadler et al. (2003) nghiên cứu hạt hợp kim nhôm ($Al\text{-}2%Ca$) cho thấy bề mặt ngoài xuất hiện chóp nhọn trong khi bề mặt bên trong nhân khí giữ độ nhẵn bóng đặc thù. Về mặt kỹ thuật nhũ tương, Buyukkestelli et al. (2021) chứng minh hạt đa thành phần nước-dầu-nước ($W/O/W$) làm tăng độ ngọt cảm quan lên tới $75%$ so với hạt đơn lẻ chứa cùng hàm lượng đường Sucrose, mở ra tiềm năng thực phẩm chức năng. Biancardo et al. (2007) chỉ ra hạt bán dẫn Silic dạng cầu chỉ sử dụng lớp vỏ mỏng quang điện bên ngoài để đạt hiệu suất chuyển đổi $13,5%$ dưới bức xạ $100 \text{ mW/cm}^2$, khẳng định chế tạo hạt bán dẫn rỗng sẽ tiết kiệm nguyên liệu mà không suy giảm công suất.

                  ┌──────────────────────────────────────────────────────────┐
                  │                 TIẾP CẬN Y VĂN TRƯỚC ĐÂY                 │
                  └─────────────────────────────┬────────────────────────────┘
                                                │
         ┌──────────────────────────────────────┼──────────────────────────────────────┐
         ▼                                      ▼                                      ▼
┌──────────────────┐                  ┌──────────────────┐                  ┌──────────────────┐
│  MÔ HÌNH LÝ THUYẾT│                  │    THỰC NGHIỆM   │                  │  MÔ PHỎNG SỐ     │
│ Snoijer & Brunet │                  │ Huang et al.     │                  │ Level-Set, VOF,  │
│ Anderson et al.  │                  │ Bhagat et al.    │                  │ LBM              │
│ (Chỉ hạt đơn pha)│                  │ (Hạt Ge rỗng)    │                  │ (Đa số hạt đơn)  │
└────────┬─────────┘                  └────────┬─────────┘                  └────────┬─────────┘
         │                                     │                                     │
         └─────────────────────────────────────┼─────────────────────────────────────┘
                                               │
                                               ▼
                         ┌───────────────────────────────────────────┐
                         │      RESEARCH GAP ĐƯỢC XÁC ĐỊNH           │
                         │ Thiếu mô hình số hoàn chỉnh cho hạt rỗng   │
                         │ đa pha có truyền nhiệt, chuyển pha và     │
                         │ tương tác động học mặt phân cách kép      │
                         └─────────────────────┬─────────────────────┘
                                               │
                                               ▼
                         ┌───────────────────────────────────────────┐
                         │           ĐÓNG GÓP CỦA LUẬN ÁN            │
                         │ Phát triển mã nguồn Front-Tracking (Orz)  │
                         │ Mô phỏng chính xác hóa rắn hạt rỗng trên  │
                         │ đế lạnh và hạt lơ lửng có đối lưu         │
                         └───────────────────────────────────────────┘

Về mặt mô phỏng số, tồn tại tranh luận lớn giữa các phương pháp:

  • Phương pháp tập mức (Multiple Level-Set) của Shetabivash et al. (2014) và Vahab et al. (2019) theo dõi mặt phân cách gián tiếp qua hàm khoảng cách có dấu, nhưng gặp khó khăn lớn trong việc bảo toàn khối lượng và tính toán chính xác lực căng bề mặt cục bộ.
  • Phương pháp thể tích lưu chất (VOF) của Zhou et al. (2019) và Chen et al. (2020) bảo toàn khối lượng tốt nhưng tái cấu trúc mặt phân cách dạng đoạn thẳng (PLIC) gây sai số tính toán độ cong $\kappa$.
  • Phương pháp lưới Boltzmann (LBM) của Sun et al. (2019) linh hoạt về vi cấu trúc nhưng tiêu tốn tài nguyên và khó hội tụ ở tỷ số khối lượng riêng lớn.

Luận án định vị phương pháp theo dấu biên (Front-Tracking Method) như một giải pháp ưu việt kết hợp lưới Euler tĩnh cho trường vận tốc - áp suất - nhiệt độ và mạng lưới nút Lagrange chuyển động cho các mặt phân cách. Điều này giải quyết triệt để bài toán bước nhảy tính chất vật lý và tính chính xác lực căng bề mặt thông qua phép gán nguồn Dirac.

Đóng góp lý thuyết và khung phân tích

Đóng góp cho lý thuyết

Luận án mở rộng trực tiếp lý thuyết hóa rắn cổ điển của Stefan và lý thuyết động học mặt phân cách chuyển pha của Snoijer & Brunet (2012) cũng như Anderson et al. (2010) sang hệ hạt đa thành phần có cấu trúc vỏ lỏng - nhân khí:

  1. Mở rộng lý thuyết điểm chập ba pha kép: Luận án chứng minh sự tồn tại đồng thời của hai hệ điểm chập ba pha: điểm chập ba ngoài (Rắn - Lỏng - Khí ngoài) và điểm chập ba trong (Rắn - Lỏng - Khí nhân). Tốc độ dịch chuyển và góc tiếp xúc của hai điểm này tương tác phi tuyến tính thông qua áp suất mao dẫn Laplace và trường nhiệt dẫn truyền qua lớp vỏ mỏng.
  2. Cơ chế hình thành chóp nhọn đa diện: Luận án làm sáng tỏ khi $\rho_{sl} < 1$ (vật liệu giãn nở thể tích khi đông đặc như nước tinh khiết), chất lỏng bị nén giữa hai biên hóa rắn hội tụ đẩy khối chất lỏng còn lại nhô lên tạo đỉnh chóp ngoài, trong khi mặt trong của nhân khí bị biến dạng lõm sâu hoặc hình thành màng mỏng tùy thuộc vào tỷ số bán kính $R_{io}$.
  3. Đặc trưng hóa thời gian đóng băng phi tuyến: Luận án thiết lập hàm phụ thuộc của thời gian hóa rắn toàn phần $t_s^$ vào số Stefan ($St$) và số Bond ($Bo$), chỉ ra rằng $t_s^$ tỷ lệ nghịch phi tuyến với $St$ do thông lượng nhiệt giải phóng tỷ lệ thuận với gradient nhiệt độ tại biên chuyển pha $q = k_s \nabla T_s - k_l \nabla T_l$.

Khung phân tích độc đáo

Khung phân tích của luận án tích hợp ba lý thuyết vật lý cốt lõi: Động lực học chất lưu phi nén Navier-Stokes, Lý thuyết truyền nhiệt Fourier có nguồn nhiệt ẩn biến đổi pha, và Lý thuyết sức căng bề mặt vi phân Young-Laplace.

                    ┌────────────────────────────────────────────────────────┐
                    │               KHUNG PHÂN TÍCH TÍCH HỢP                │
                    └───────────────────────────┬────────────────────────────┘
                                                │
         ┌──────────────────────────────────────┼──────────────────────────────────────┐
         ▼                                      ▼                                      ▼
┌──────────────────┐                  ┌──────────────────┐                  ┌──────────────────┐
│  PHƯƠNG TRÌNH    │                  │  PHƯƠNG TRÌNH    │                  │  ĐIỀU KIỆN BIÊN  │
│  NAVIER-STOKES   │                  │   NĂNG LƯỢNG     │                  │  CHUYỂN PHA      │
│  (Động lượng &   │                  │  (Truyền nhiệt & │                  │ (Cân bằng nhiệt, │
│  Lực căng δ)     │                  │   Nhiệt ẩn L_h)  │                  │  Nhảy thể tích)  │
└────────┬─────────┘                  └────────┬─────────┘                  └────────┬─────────┘
         │                                     │                                     │
         └─────────────────────────────────────┼─────────────────────────────────────┘
                                               │
                                               ▼
                         ┌───────────────────────────────────────────┐
                         │   PHƯƠNG PHÁP THEO DẤU BIÊN (LAGRANGE)    │
                         │      (Front-Tracking Interface)           │
                         └─────────────────────┬─────────────────────┘
                                               │
                                               ▼
                         ┌───────────────────────────────────────────┐
                         │    LƯỚI SO LE EULER (MAC STAGGERED GRID)  │
                         │    Hệ tọa độ trụ đối xứng trục (Orz)      │
                         └───────────────────────────────────────────┘

Điểm độc đáo trong khung phân tích số là:

  • Xây dựng toán tử hàm chỉ thị kép ($I_s, I_l$): Phân tách trường tính toán thành 3 miền pha độc lập (Rắn, Lỏng, Khí) trên cùng một lưới Euler đồng nhất, gán thuộc tính vật lý trơn mượt: $$\rho(\mathbf{x}) = \rho_g + (\rho_l - \rho_g)I_l(\mathbf{x}) + (\rho_s - \rho_g)I_s(\mathbf{x})$$
  • Điều kiện biên không trượt nội suy tại pha rắn: Áp đặt vectơ lực ảo $\mathbf{f}$ vào phương trình động lượng thông qua bước hiệu chỉnh vận tốc $\mathbf{u}_F = \mathbf{u}_s = 0$, loại bỏ việc tái chia lưới phức tạp quanh thể rắn đang phát triển.
  • Điều kiện biên biên giới (Boundary conditions): Xét trong miền tính toán $W \times H$, áp dụng điều kiện biên đối xứng trục tại $r = 0$, biên trượt tự do (free-slip) tại các vách biên ngoài $r = W$ và $z = H$, và điều kiện nhiệt độ đẳng nhiệt $T = T_c < T_m$ tại đáy làm lạnh $z = 0$.

Phương pháp nghiên cứu tiên tiến

Thiết kế nghiên cứu

Nghiên cứu đứng trên lập trường bản thể luận khách quan và phương pháp luận thực chứng (positivism/computational determinism), sử dụng phương pháp sai phân hữu hạn (Finite Difference Method - FDM) kết hợp thuật toán theo dấu biên (Front-Tracking Method). Mô hình toán học được xây dựng trên hệ tọa độ trụ đối xứng trục $(Orz)$, giảm thiểu không gian tính toán 3D xuống 2D tương đương nhưng vẫn đảm bảo trọn vẹn đặc trưng hình học 3 chiều của cấu trúc cầu xoay.

Hệ thống phương trình vi phân đạo hàm riêng chi phối gồm:

  1. Phương trình động lượng: $$\frac{\partial (\rho \mathbf{u})}{\partial t} + \nabla \cdot (\rho \mathbf{u}\mathbf{u}) = -\nabla p + \nabla \cdot [\mu (\nabla \mathbf{u} + \nabla \mathbf{u}^T)] + \int_{S_f} \sigma \kappa \mathbf{n}_f \delta(\mathbf{x} - \mathbf{x}_f) dS + \mathbf{f} + \rho \mathbf{g}$$
  2. Phương trình năng lượng có chuyển pha: $$\frac{\partial (\rho C_p T)}{\partial t} + \nabla \cdot (\rho C_p T \mathbf{u}) = \nabla \cdot (k \nabla T) + \int_{S_f} q \delta(\mathbf{x} - \mathbf{x}_f) dS$$
  3. Phương trình liên tục có tính đến biến thiên thể tích khi hóa rắn: $$\nabla \cdot \mathbf{u} = \frac{1}{L_h} \left( \frac{1}{\rho_s} - \frac{1}{\rho_l} \right) \int_{S_f} q \delta(\mathbf{x} - \mathbf{x}_f) dS$$

Quy trình nghiên cứu rigorous

Quy trình giải thuật số được lập trình độc lập bằng ngôn ngữ Fortran chuyên dụng cho tính toán khoa học hiệu năng cao, thực hiện tuần tự qua các bước:

  1. Cập nhật vị trí mặt phân cách Lagrange: Mặt phân cách được rời rạc hóa thành chuỗi các điểm kết nối $\mathbf{x}_f^{n+1} = \mathbf{x}_f^n + \mathbf{V}_f \Delta t$, với $\mathbf{V}_f$ được nội suy từ trường vận tốc lưới Euler hoặc xác định từ điều kiện cân bằng năng lượng biên hóa rắn: $$V_n = \frac{1}{\rho_s L_h} \left( k_s \left. \frac{\partial T_s}{\partial n} \right|_f - k_l \left. \frac{\partial T_l}{\partial n} \right|_f \right)$$
  2. Phân phối lực căng bề mặt và thông lượng nhiệt: Các đại lượng biên được chuyển giao sang lưới Euler bằng hàm trọng số diện tích (Area-weighting Peskin kernel).
  3. Giải phương trình áp suất Poisson: Áp dụng thuật toán giải lặp SOR (Successive Over-Relaxation) hoặc GMRES với sai số hội tụ tuyệt đối $\text{error} < 10^{-6}$.
  4. Khử trôi thể tích và tái cấu trúc lưới Lagrange (Restructuring): Tự động thêm/bớt điểm liên kết khi khoảng cách giữa hai nút liền kề $\Delta s$ vượt ngưỡng $[0,5 h; 1,5 h]$ (với $h$ là kích thước mắt lưới Euler), đảm bảo độ cong $\kappa$ tính toán chuẩn xác.
  ┌────────────────────────────────────────────────────────────┐
  │         BƯỚC 1: KHỞI TẠO HÌNH HỌC VÀ TRƯỜNG ĐẦU            │
  │     Lưới MAC, Phân bố hạt lỏng & nhân khí, T = T0          │
  └─────────────────────────────┬──────────────────────────────┘
                                │
                                ▼
  ┌────────────────────────────────────────────────────────────┐
  │     BƯỚC 2: THEO DÕI MẶT PHÂN CÁCH LAGRANGE (xf, nf, κ)    │
  │        Tái cấu trúc điểm biên, Cập nhật Is, Il             │
  └─────────────────────────────┬──────────────────────────────┘
                                │
                                ▼
  ┌────────────────────────────────────────────────────────────┐
  │     BƯỚC 3: TÍNH TOÁN NGUỒN BIÊN TRÊN NỀN LƯỚI EULER       │
  │ Lực căng bề mặt (Dirac delta), Nguồn nhiệt chuyển pha q    │
  └─────────────────────────────┬──────────────────────────────┘
                                │
                                ▼
  ┌────────────────────────────────────────────────────────────┐
  │     BƯỚC 4: GIẢI PHƯƠNG TRÌNH NĂNG LƯỢNG & NAVIER-STOKES   │
  │  Cập nhật trường T, Dự báo vận tốc u*, Giải áp suất p      │
  └─────────────────────────────┬──────────────────────────────┘
                                │
                                ▼
  ┌────────────────────────────────────────────────────────────┐
  │     BƯỚC 5: TÍNH VẬN TỐC BIÊN HÓA RẮN Vn VÀ DỊCH CHUYỂN    │
  │          xf(n+1) = xf(n) + Vf * Δt -> Bước tiếp theo       │
  └────────────────────────────────────────────────────────────┘

Data và phân tích

Kiểm chứng độ hội tụ lưới (Grid Convergence Study) được thực hiện nghiêm ngặt trên nhiều cấp độ phân giải:

  • Đối với bài toán hóa rắn trên bề mặt lạnh: khảo sát lưới $128 \times 256$ so sánh với $256 \times 512$, độ lệch tương đối về vị trí biên hóa rắn và chiều cao cực đại $H_d$ đạt dưới $1,5%$, đảm bảo tính độc lập lưới.
  • Đối với bài toán đối lưu cưỡng bức: thực hiện so sánh các độ phân giải lưới từ $128 \times 512$ đến $384 \times 1536$. Lưới chuẩn $256 \times 1024$ được lựa chọn nhằm tối ưu giữa độ chính xác và thời gian tính toán.
  • Kiểm chứng thực nghiệm độc lập: Kết quả mô phỏng cho trường hợp hạt nước đơn pha ($\theta_0 = 85^\circ$) được so sánh trực tiếp với dữ liệu thực nghiệm quang học của Zhang et al. (2017) và dữ liệu góc ướt lớn ($\theta_0 = 124^\circ, 155^\circ$) của Huang et al. (2012). Động học tăng chiều cao đỉnh hạt và biên dạng đóng băng cho thấy sự trùng khớp với sai số thể tích dưới $3%$.

Phát hiện đột phá và implications

Những phát hiện then chốt

  1. Phát hiện về biến dạng không đối xứng của nhân bóng khí: Khi hạt rỗng hóa rắn trên bề mặt lạnh, lớp vỏ lỏng phía đáy tiếp xúc với đế lạnh sẽ đóng băng trước, đẩy nhân bóng khí dịch chuyển tịnh tiến lên phía trên. Chiều cao hạt bóng khí kết thúc ($H_{in}$) tăng tuyến tính theo tỷ số bán kính ban đầu $R_{io}$ và tỷ số mật độ $\rho_{sl}$. Khi $\rho_{sl} = 0,8$, nhân khí bị kéo dài theo phương thẳng đứng rõ rệt so với $\rho_{sl} = 1,1$.
  2. Quy luật chi phối của số Bond ($Bo$) và số Stefan ($St$): Khi số Bond tăng từ $Bo = 0,18$ lên $Bo = 1,78$, lực trọng trường chiếm ưu thế làm hạt lưu chất bị dàn trải rộng hơn trên bề mặt lạnh, bán kính ướt hiệu dụng tăng, dẫn đến diện tích truyền nhiệt tăng làm giảm thời gian hóa rắn toàn phần $t_s^*$. Ngược lại, khi số Stefan tăng ($St = 0,032 \rightarrow 0,316$), gradient nhiệt độ qua lớp băng tăng mạnh, đẩy nhanh tốc độ truyền nhiệt và rút ngắn thời gian đông kết theo hàm mũ nghịch đảo.
  3. Hiện tượng tương tác góc ướt trong ($\theta_{0i}$) và ngoài ($\theta_{0o}$): Luận án phát hiện rằng góc ướt ngoài $\theta_{0o}$ kiểm soát hình thái đỉnh chóp ngoài và diện tích tiếp xúc nhiệt, trong khi góc ướt trong $\theta_{0i}$ quyết định mức độ co thắt của khoang rỗng. Độ tăng chiều cao hạt ngoài $(H_d - h_{0o})$ tỷ lệ thuận với $\theta_{0o} \in [60^\circ; 130^\circ]$.
  4. Hiệu ứng tạo đuôi xoáy trong dòng đối lưu cưỡng bức: Đối với hạt rỗng lơ lửng, số Reynolds ($Re$) tăng làm biến dạng trường nhiệt độ phía sau hạt (vùng wake), làm giảm độ dày vỏ băng ở vùng đuôi và tạo ra sự mất đối xứng hình học giữa bán cầu đón gió và bán cầu khuất gió. Tỷ số hình dạng trong $Ar_i$ và ngoài $Ar_o$ phụ thuộc mạnh vào độ lệch tâm ban đầu $\varepsilon_0$. Khi $\varepsilon_0 \neq 0$, ứng suất nhiệt không đều làm nhân khí bị đẩy lệch tâm nghiêm trọng, có nguy cơ phá vỡ lớp vỏ lỏng trước khi hóa rắn hoàn toàn.
Thông số biến thiên Dải khảo sát Ảnh hưởng đến thời gian hóa rắn ($t_s^*$) Ảnh hưởng đến hình thái hạt và khoang rỗng
Số Bond ($Bo$) $0,18 - 1,78$ Giảm khi $Bo$ tăng (do tăng diện tích tiếp xúc đáy) Hạt bị bẹp theo phương ngang, đỉnh chóp thấp hơn
Số Stefan ($St$) $0,032 - 0,316$ Giảm mạnh theo hàm mũ Vận tốc biên $V_n$ tăng cao, chóp nhọn sắc nét
Tỷ số mật độ ($\rho_{sl}$) $0,8 - 1,1$ Tăng nhẹ khi $\rho_{sl}$ tăng $\rho_{sl} < 1$ tạo chóp nhọn lồi; $\rho_{sl} \ge 1$ đỉnh tù/phẳng
Tỷ số bán kính ($R_{io}$) $0,3 - 0,7$ Giảm mạnh khi $R_{io}$ tăng (thể tích lỏng giảm) Nhân khí kéo dài trục $z$, vỏ mỏng dễ biến dạng
Số Reynolds ($Re$) $10 - 100$ Giảm ở mặt trước, tăng ở vùng xoáy đuôi Bất đối xứng trước-sau, $Ar_o$ lệch chuẩn

Implications đa chiều

  • Về mặt lý thuyết: Cung cấp mô hình định lượng đầu tiên mô tả chính xác tương tác chuyển pha có mặt đồng thời hai mặt phân cách tự do và bước nhảy thể tích. Đặt cơ sở mở rộng cho lý thuyết đông kết đa pha trong cơ học vi lưu (microfluidics).
  • Về mặt phương pháp luận: Khẳng định tính hiệu quả vượt trội của phương pháp Front-Tracking trong việc loại bỏ hoàn toàn hiện tượng khuếch tán số (numerical diffusion) tại biên phân cách pha, tạo chuẩn đối sánh (benchmark) cho các phương pháp lưới động khác.
  • Về ứng dụng công nghiệp:
    • Chế tạo pin mặt trời hình cầu: Kết quả về tỷ số bán kính $R_{io} \in [0,5; 0,7]$ và sự ổn định hình học cho phép định cỡ tối ưu quá trình phun hạt Silic rỗng, giúp giảm $40 - 60%$ khối lượng vật liệu bán dẫn đắt tiền mà vẫn duy trì diện tích bề mặt quang điện tương đương.
    • Kỹ thuật chống đóng băng hàng không và phong điện: Cung cấp dữ liệu chính xác về thời gian trễ hóa rắn và cơ chế bám dính của giọt nước chứa bọt khí trên cánh máy bay và cánh tuabin gió, hỗ trợ thiết kế bề mặt kỵ nước siêu cấp (superhydrophobic) với góc ướt tối ưu $\theta > 150^\circ$.

Limitations và Future Research

Luận án thẳng thắn chỉ ra các giới hạn nghiên cứu xuất phát từ giả thiết mô hình:

  1. Giới hạn đối xứng trục 2D: Mô hình giả định đối xứng trục $(Orz)$, chưa nắm bắt được các cấu trúc bất đối xứng 3D thực tế như sự dao động bề mặt phi đối xứng (non-axisymmetric oscillations) hay hiện tượng vỡ bọt khí ngẫu nhiên.
  2. Giả thiết thuộc tính nhiệt hằng số: Các hệ số dẫn nhiệt ($k$), nhiệt dung ($C_p$) và độ nhớt ($\mu$) được xem là không đổi trong từng pha, chưa xét đến tính phụ thuộc phi tuyến vào nhiệt độ cục bộ trong dải quá lạnh sâu (deep supercooling).
  3. Bỏ qua hiện tượng quá lạnh và mầm kết tinh động học: Quá trình chuyển pha được giả định xảy ra tức thời tại nhiệt độ nóng chảy chuẩn $T_m$, chưa tích hợp động học mầm tinh thể Gibbs-Thomson.

Các hướng phát triển tương lai gồm:

  • Mở rộng mã nguồn Front-Tracking sang không gian ba chiều (3D) hoàn chỉnh trên nền tảng tính toán song song đa luồng (MPI/OpenMP).
  • Tích hợp hiệu ứng quá lạnh cục bộ và giải bài toán đa nhân khí phân tán (multiple core droplets).
  • Khảo sát dòng lưu chất phi Newton (Non-Newtonian fluid) ứng dụng trong đóng gói dược phẩm và polymer sinh học.

Tác động và ảnh hưởng

Công trình đã tạo ra tác động học thuật rõ rệt thông qua các bài báo khoa học công bố trên các tạp chí chuyên ngành quốc tế uy tín thuộc danh mục Scopus/SCIE và các tạp chí chuyên ngành Cơ học trong nước. Kết quả nghiên cứu là tài liệu tham khảo giá trị cho các nhóm nghiên cứu về truyền nhiệt, vi lưu và chuyển pha trên toàn thế giới.

Về mặt công nghệ, nghiên cứu định hình trực tiếp các giải pháp kỹ thuật cho ngành năng lượng tái tạo và chế tạo vật liệu tiên tiến:

  • Tiết kiệm hàng triệu USD chi phí nguyên liệu Silicon bán dẫn thông qua giải pháp công nghệ hạt rỗng.
  • Cung cấp cơ sở dữ liệu số cho các nhà sản xuất máy bay và tuabin gió tối ưu hóa hệ thống sưởi nhiệt chống đóng băng cánh, nâng cao tuổi thọ máy móc và đảm bảo an toàn hàng không.

Đối tượng hưởng lợi

  • Nghiên cứu sinh và học viên sau đại học ngành Cơ học lưu chất & Truyền nhiệt: Thụ hưởng một khung phương pháp luận hoàn chỉnh, thuật toán mẫu bằng Fortran và cơ sở dữ liệu kiểm chứng về phương pháp Front-Tracking.
  • Các nhà nghiên cứu vật lý nhiệt và vi lưu học: Tiếp cận dữ liệu định lượng về động thái mặt phân cách kép và cơ chế chuyển pha phức hợp.
  • Kỹ sư R&D công nghệ năng lượng mặt trời: Ứng dụng quy luật ổn định hình học hạt rỗng để thiết kế dây chuyền sản xuất vi hạt bán dẫn hiệu năng cao.
  • Kỹ sư thiết kế khí động lực học hàng không & phong điện: Sử dụng các biểu đồ thời gian hóa rắn và quy luật phát triển băng đá để tối ưu hóa vật liệu phủ bề mặt chống bám băng.

Câu hỏi chuyên sâu

1. Đóng góp lý thuyết độc đáo nhất của luận án là gì và đã mở rộng lý thuyết nào? Trả lời: Đóng góp độc đáo nhất là việc mở rộng lý thuyết hóa rắn giọt lỏng của Snoijer & Brunet (2012) và Anderson et al. (2010) sang hệ hạt rỗng hai mặt phân cách tự do. Luận án đã mô hình hóa thành công tương tác động học giữa hai điểm chập ba pha đồng thời và chứng minh vai trò của tỷ số bán kính $R_{io}$ cùng tỷ số mật độ $\rho_{sl}$ trong việc điều khiển hình thái chóp nhọn đỉnh hạt và độ biến dạng khoang rỗng.

2. Điểm cải tiến phương pháp luận cốt lõi so với các công trình trước đây là gì? Trả lời: So với phương pháp Level-Set của Shetabivash et al. (2014) và phương pháp VOF của Chen et al. (2020), luận án sử dụng phương pháp Front-Tracking giải trực tiếp lực căng bề mặt thông qua hàm delta Dirac trên lưới Lagrange rời rạc. Kỹ thuật này triệt tiêu hoàn toàn sự khuếch tán số tại mặt phân cách, tính toán chính xác độ cong $\kappa$ và bước nhảy áp suất Laplace mà không cần hàm làm trơn nhân tạo.

3. Phát hiện bất ngờ nhất từ dữ liệu mô phỏng số là gì? Trả lời: Đó là hiện tượng nhân bóng khí bên trong không giữ nguyên vị trí tâm mà bị biến dạng kéo dài dọc trục và dịch chuyển mạnh lên phía trên do ứng suất nén sinh ra từ lớp vỏ băng đáy giãn nở ($\rho_{sl} = 0,8$), kết hợp với việc thời gian hóa rắn toàn phần $t_s^*$ giảm mạnh khi số Bond tăng (ngược với trực giác hạt đơn pha phẳng dẹt truyền nhiệt chậm).

4. Quy trình lặp và kiểm chứng tái lập (Replication Protocol) có được cung cấp đầy đủ không? Trả lời: Luận án cung cấp chi tiết toàn bộ hệ phương trình vi phân không thứ nguyên, cấu trúc lưới so le MAC, giải thuật rời rạc hóa phương trình Navier-Stokes và năng lượng, hàm trọng số diện tích, cùng các tiêu chuẩn hội tụ lưới ($128 \times 256, 256 \times 512, 384 \times 1536$). Mã nguồn Fortran được chuẩn hóa cho phép tái lập chính xác mọi kết quả mô phỏng.

5. Chương trình nghiên cứu 10 năm tiếp theo được định hình ra sao? Trả lời: Lộ trình tập trung vào việc nâng cấp thuật toán sang hệ 3D song song khối lượng lớn (Massive Parallel DNS), tích hợp trường nhiệt động phi cân bằng, mô phỏng quá trình va đập kèm đóng băng tức thời của chùm hạt đa thành phần (droplet spray icing), và thử nghiệm ứng dụng trong in 3D kim loại lỏng vi mô.

Kết luận

Luận án tiến sĩ của NCS. Phạm Duy Bính đã hoàn thành xuất sắc các mục tiêu nghiên cứu với 5 đóng góp cốt lõi:

  1. Xây dựng thành công bộ chương trình tính toán số hoàn chỉnh bằng Fortran dựa trên phương pháp theo dấu biên (Front-Tracking Method) trên hệ tọa độ trụ $(Orz)$ chuyên biệt cho bài toán truyền nhiệt, chuyển pha của hạt lưu chất đa thành phần.
  2. Thiết lập cơ chế vật lý chi tiết về sự hình thành chóp nhọn đỉnh ngoài và biến dạng nhân khí bên trong dưới ảnh hưởng phối hợp của các số không thứ nguyên $Bo, Pr, St, \rho_{sl}, R_{io}$.
  3. Xác định định lượng quy luật biến thiên của thời gian hóa rắn toàn phần $t_s^*$ theo góc ướt trong ($\theta_{0i}$), góc ướt ngoài ($\theta_{0o}$) và góc phát triển tiếp xúc ($\theta_{gr}$).
  4. Làm sáng tỏ động lực học đông kết của hạt rỗng lơ lửng chịu dòng đối lưu cưỡng bức, phát hiện sự mất đối xứng hình thái và ảnh hưởng nghiêm trọng của độ lệch tâm ban đầu $\varepsilon_0$.
  5. Kiểm chứng tính chính xác của mô hình số thông qua sự trùng khớp chặt chẽ với các dữ liệu thực nghiệm quốc tế uy tín của Zhang et al. và Huang et al.

Công trình tạo nền tảng học thuật vững chắc, mở ra hướng nghiên cứu mới trong cơ học thủy khí hiện đại và mang lại giá trị ứng dụng thực tiễn cao trong ngành năng lượng tái tạo và kỹ thuật hàng không.