Tổng quan luận án
Luận án tiến sĩ vật lý với tiêu đề "Nghiên cứu tối ưu thay đảo nhiên liệu lò phản ứng hạt nhân VVER-1000" (tên tiếng Anh: "Study on fuel loading pattern optimization for VVER-1000 nuclear reactor") được tác giả Trần Việt Phú thực hiện tại Trung tâm Đào tạo hạt nhân thuộc Viện Năng lượng nguyên tử Việt Nam, dưới sự hướng dẫn của PGS. Trần Hoài Nam và GS. Yamamoto Akio, bảo vệ năm 2022. Luận án thuộc chuyên ngành Vật lý Nguyên tử và Hạt nhân, mã số chuyên ngành 9.06 (theo phân loại trong văn bản).
Về tính cấp thiết và khoảng trống nghiên cứu, tác giả trình bày bối cảnh nhu cầu sử dụng điện năng ngày càng tăng trên phạm vi toàn cầu, đòi hỏi bảo đảm an ninh năng lượng và giảm thiểu phát thải khí nhà kính gây biến đổi khí hậu. Trong cơ cấu kinh tế của nhà máy điện hạt nhân, chi phí đầu tư xây dựng chiếm khoảng 70% tổng chi phí và thường cố định ngay từ giai đoạn chuẩn bị xây dựng, trong khi chi phí nhiên liệu chiếm khoảng 15–25% tổng chi phí và có thể tối ưu hóa được trong suốt quá trình vận hành thông qua công tác quản lý nhiên liệu trong vùng hoạt (In-core Fuel Management - ICFM). Đối với các lò phản ứng nước nhẹ (LWR) nói chung và lò nước áp lực VVER nói riêng, sau mỗi chu kỳ vận hành thông thường từ 12 đến 18 tháng, khoảng 1/3 số bó nhiên liệu đã cháy kiệt được đưa ra bể chứa, 2/3 số bó nhiên liệu còn lại kết hợp cùng các bó nhiên liệu mới được sắp xếp lại để tạo thành cấu hình nạp tải (Loading Pattern - LP) cho chu kỳ mới.
Bài toán tối ưu hóa cấu hình nạp nhiên liệu là một bài toán tối ưu hóa tổ hợp đa mục tiêu rất phức tạp với không gian tìm kiếm cực lớn và chứa đựng các mục tiêu có xu hướng triệt tiêu lẫn nhau:
- Tối đa hóa độ dài chu kỳ nhiên liệu nhằm nâng cao hiệu quả kinh tế sử dụng nhiên liệu, thường được quy đổi thành bài toán tối đa hóa hệ số nhân hiệu dụng ($k_{eff}$) tại thời điểm bắt đầu chu kỳ (Beginning of Cycle - BOC) nhằm giảm rò rỉ nơtron ra vùng biên và tiết kiệm thời gian tính toán độ cháy;
- Duy trì các giới hạn an toàn nhiệt thủy lực của vùng hoạt trong suốt quá trình vận hành, tiêu biểu là việc giảm thiểu hệ số đỉnh công suất hướng tâm (radial Power Peaking Factor - PPF) nhằm giảm thiểu xác suất hư hỏng thanh nhiên liệu và tăng cường biên độ an toàn.
Do vùng hoạt lò VVER-1000 gồm 163 bó nhiên liệu lục giác với nhiều mức làm giàu và độ cháy khác nhau, số lượng cấu hình khả dĩ là vô cùng lớn (lên tới $4,5 \times 10^{10}$ cấu hình ngay cả khi đã áp dụng tính đối xứng quay 1/6 và xét tính tương đồng của các bó nhiên liệu). Các thuật toán cổ điển dễ bị rơi vào cực trị địa phương, trong khi các thuật toán ngẫu nhiên đơn lẻ như mô phỏng tôi kim truyền thống (SA) lại có tốc độ hội tụ chậm. Khoảng trống nghiên cứu then chốt mà luận án xác định là việc thiếu các phương pháp tối ưu hóa ngẫu nhiên cải tiến có tốc độ hội tụ cao, đồng thời thiếu một công cụ mô phỏng vật lý nơtron chuyên dụng giải nhanh và chính xác phương trình khuếch tán trên lưới tam giác cho lò VVER-1000 để phục vụ quá trình tìm kiếm tự động.
Mục tiêu nghiên cứu của luận án được xác định rõ ràng thành hai mục tiêu cụ thể:
- Nghiên cứu và áp dụng các phương pháp tối ưu hóa tiên tiến (gồm phương pháp Mô phỏng tôi kim tiến hóa - ESA và phương pháp Tiến hóa vi phân thích ứng dựa trên lịch sử thành công dạng rời rạc - discrete SHADE) vào bài toán tối ưu hóa cấu hình nạp tải nhiên liệu của lò phản ứng hạt nhân nhằm tiết kiệm chi phí tính toán, nâng cao tốc độ hội tụ và tăng xác suất tiếp cận nghiệm tối ưu toàn cục;
- Xây dựng và phát triển một công cụ tính toán mô phỏng vật lý vùng hoạt chuyên dụng (chương trình LPO-V) cho lò phản ứng VVER với cấu trúc hình học lưới tam giác, kết hợp trực tiếp với các thuật toán tìm kiếm tối ưu hóa để thực hiện bài toán thay đảo nhiên liệu.
Đối tượng và phạm vi nghiên cứu:
- Đối tượng nghiên cứu: Cấu hình sắp xếp các bó nhiên liệu trong vùng hoạt của lò phản ứng hạt nhân nước áp lực VVER-1000 nạp 30% nhiên liệu oxit hỗn hợp U-Pu (MOX) và 70% nhiên liệu oxit urani ($UO_2$);
- Phạm vi không gian: Mô hình hai chiều (2D) của vùng hoạt lò VVER-1000 có tính đối xứng quay góc 60 độ (mô hình 1/6 vùng hoạt gồm 28 vị trí bó nhiên liệu), bao gồm vùng hoạt chứa các bó nhiên liệu, khe nước, vách ngăn thép có lỗ nước, thùng đỡ, khoang hạ lưu và thùng áp lực lò phản ứng, với điều kiện biên ngoài là biên chân không;
- Phạm vi thời gian và trạng thái: Khảo sát tại thời điểm bắt đầu chu kỳ nạp tải (BOC) ở các trạng thái vận hành chuẩn theo bài toán benchmark của OECD/NEA, tập trung vào trạng thái làm việc chuẩn (State S1) và không tính toán diễn biến suy giảm độ cháy theo thời gian trong quá trình tìm kiếm tối ưu vì hệ số nhân vô hạn $k_{\infty}$ của toàn bộ các bó nhiên liệu đều suy giảm đơn điệu theo độ cháy.
| Thông số cấu hình vùng hoạt VVER-1000 MOX |
Giá trị số lượng cấu hình khả dĩ |
Thời gian tính toán ước tính (năm)* |
| Không xét đối xứng, không xét tính tương đồng của bó nhiên liệu |
$1,9 \times 10^{280}$ |
$6,0 \times 10^{272}$ |
| Không xét đối xứng, có xét tính tương đồng của bó nhiên liệu |
$3,2 \times 10^{120}$ |
$1,0 \times 10^{113}$ |
| Đối xứng 1/6 vùng hoạt, không xét tính tương đồng của bó nhiên liệu |
$1,5 \times 10^{19}$ |
$4,8 \times 10^{10}$ |
| Đối xứng 1/6 vùng hoạt, có xét tính tương đồng của bó nhiên liệu |
$4,5 \times 10^{10}$ |
$1,4 \times 10^{2}$ |
*Ghi chú: Giả định tốc độ tính toán trung bình đạt 10 cấu hình/giây.
Tổng quan tài liệu và vị trí của luận án
Trong phần tổng quan tài liệu, tác giả điểm lại quá trình phát triển của các phương pháp tối ưu hóa áp dụng cho bài toán quản lý nhiên liệu trong vùng hoạt (ICFM) từ những giai đoạn đầu của công nghệ lò phản ứng hạt nhân đến nay. Các phương pháp được phân loại thành hai nhóm chính:
- Các thuật toán tìm kiếm dựa trên vùng lân cận (Neighborhood-based algorithms):
- Các phương pháp cổ điển như Tìm kiếm trực tiếp (Direct Search - DS) và Hoán vị nhị phân (Binary Exchange - BE) được tác giả trích dẫn theo công trình của Parks (1992) [6]. Các thuật toán này đánh giá số lượng cấu hình nhỏ và hội tụ nhanh về cực trị địa phương nhưng hoàn toàn không có khả năng thoát khỏi cực trị địa phương.
- Phương pháp Mô phỏng tôi kim (Simulated Annealing - SA) được trích dẫn qua các công trình nghiên cứu [6–15]. SA có khả năng thoát khỏi cực trị địa phương nhờ xác suất chấp nhận các nghiệm kém hơn theo phân bố Boltzmann, song nhược điểm lớn là tốc độ hội tụ rất chậm, đòi hỏi phải đánh giá số lượng cấu hình thử nghiệm rất lớn.
- Phương pháp Mô phỏng tôi kim thích ứng (Adaptive Simulated Annealing - ASA) được nghiên cứu nhằm tăng tốc độ hội tụ thông qua cải tiến biểu đồ giảm nhiệt độ và cơ chế sinh nghiệm thử nghiệm [9, 10, 70–72]. Trong đó, hai chiến lược chính được áp dụng là chiến lược "quay lại nghiệm tốt nhất" (return to the best) [70, 71] và chiến lược sử dụng ma trận xác suất chuyển trạng thái hoặc danh sách cấu hình hạn chế [72].
- Phương pháp Tìm kiếm cấm (Tabu Search - TS) được áp dụng để tránh việc lặp lại các nghiệm đã khảo sát trong quá khứ [16–24].
- Các thuật toán tìm kiếm dựa trên quần thể (Population-based algorithms):
- Giải thuật di truyền (Genetic Algorithms - GA) [25–34] và các phương pháp tiến hóa (Evolution Methods - EM) [35–43] mô phỏng quá trình chọn lọc tự nhiên.
- Thuật toán Tối ưu hóa bầy đàn (Particle Swarm Optimization - PSO) [44–53] và Tối ưu hóa đàn kiến (Ant Colony Optimization - ACO) [54–60].
- Thuật toán Tiến hóa vi phân (Differential Evolution - DE) do Storn và Price đề xuất năm 1997 [73], được áp dụng rộng rãi trong các bài toán kỹ thuật phức tạp nhờ khả năng hội tụ và tính bền vững cao [74–76], có khả năng thăm dò không gian tìm kiếm vượt trội so với GA [77, 78].
- Biến thể Tiến hóa vi phân thích ứng dựa trên lịch sử thành công (Success-History based Adaptive Differential Evolution - SHADE) do Tanabe và Fukunaga đề xuất năm 2013 [79], sử dụng cơ chế thích nghi tham số điều khiển dựa trên lịch sử các tham số thành công trong quá khứ, chứng minh được ưu thế so với các biến thể DE khác như JADE, CoDE và EPSDE.
- Các công trình áp dụng DE cho bài toán thay đảo nhiên liệu hạt nhân: Do DE nguyên bản được thiết kế cho biến liên tục, việc ứng dụng DE vào bài toán tổ hợp biến số nguyên đòi hỏi các chiến lược rời rạc hóa [62, 63, 80–82]. Nghiên cứu trước đây của chính nhóm tác giả [63] đã phát triển thuật toán DE rời rạc cho lò phản ứng nghiên cứu dạng bể 500 kW của Viện Nghiên cứu Hạt nhân Đà Lạt [83, 84] và chỉ ra hiệu năng vượt trội của DE so với GA.
Về các nghiên cứu tối ưu hóa cấu hình nhiên liệu dành riêng cho lò VVER:
- Các chương trình OPTIMAL (sử dụng phương pháp trao đổi nhị phân) và ATHENA (sử dụng phương pháp ngẫu nhiên) được phát triển từ năm 1994 cho lò VVER-440 Dukovany tại Cộng hòa Séc [95, 96].
- Năm 2013, nghiên cứu kết hợp phương pháp trọng số, GA và ACO được thực hiện cho lò VVER-1000 Bushehr của Iran [97]; công cụ LONSA kết hợp mạng nơ-ron nhân tạo và SA cũng được áp dụng cho lò Bushehr [98].
- Các công cụ SCAM-W, CIGARO kết hợp giải thuật di truyền và chương trình tính toán vật lý nơtron được ứng dụng cho lò VVER-440 Dukovany và VVER-1000 Temelin [99].
- Tại Việt Nam, nhiều nghiên cứu tính toán cho lò VVER đã được thực hiện tại Trung tâm Năng lượng hạt nhân [102–104], cùng đề tài cấp Bộ do nghiên cứu sinh chủ trì thực hiện giai đoạn 2015–2016 [105].
Khoảng trống mà luận án lựa chọn giải quyết: Các công trình nghiên cứu trên lò VVER vẫn đang tiếp tục vì chưa có phương pháp nào khẳng định tìm được nghiệm tối ưu toàn cục tuyệt đối trong không gian tìm kiếm khổng lồ. Luận án lựa chọn hướng đi cải tiến thuật toán SA thành thuật toán lai Evolutionary Simulated Annealing (ESA) kết hợp toán tử trao đổi chéo kiểu di truyền, đồng thời xây dựng biến thể rời rạc cho thuật toán thích ứng SHADE thông qua kỹ thuật đánh chỉ số vị trí tương đối (Relative Position Indexing - RPI). Đi đôi với đó là việc xây dựng chương trình mô phỏng vật lý nơtron LPO-V độc lập trên hệ tọa độ tam giác để phục vụ trực tiếp bài toán tối ưu hóa.
Cơ sở lý thuyết và phương pháp nghiên cứu
Cơ sở lý thuyết của mô hình vật lý nơtron trong luận án dựa trên hệ phương trình khuếch tán nơtron đa nhóm ở trạng thái dừng, được thiết lập từ phương trình liên tục bảo toàn nơtron và định luật thực nghiệm Fick ($J = -D \nabla \phi$). Phương trình khuếch tán đa nhóm đối với nhóm năng lượng $g$ ($g = 1, 2, ..., G$) có dạng tổng quát như sau:
$$-\nabla \cdot (D_g \nabla \phi_g) + \Sigma_{a,g} \phi_g + \sum_{h=1, h \neq g}^G \Sigma_{g \to h} \phi_g = \sum_{h=1, h \neq g}^G \Sigma_{h \to g} \phi_h + \frac{\chi_g}{k_{eff}} \sum_{h=1}^G \nu \Sigma_{f,h} \phi_h$$
Trong đó:
- $D_g$ là hệ số khuếch tán nơtron của nhóm $g$;
- $\phi_g$ là thông lượng nơtron nhóm $g$;
- $\Sigma_{a,g}$ là tiết diện vĩ mô hấp thụ nơtron của nhóm $g$;
- $\Sigma_{g \to h}$ và $\Sigma_{h \to g}$ là các tiết diện vĩ mô tán xạ chuyển nhóm nơtron;
- $\nu$ là số nơtron trung bình sinh ra trên mỗi phân hạch;
- $\Sigma_{f,h}$ là tiết diện vĩ mô phân hạch của nhóm $h$;
- $\chi_g$ là phổ phân hạch nơtron sinh ra trong nhóm năng lượng $g$;
- $k_{eff}$ là hệ số nhân hiệu dụng của vùng hoạt.
Về phương pháp thực tế, tác giả áp dụng phương pháp định lượng và mô phỏng số thông qua các bước phát triển cụ thể:
- Rời rạc hóa không gian bằng phương pháp sai phân hữu hạn (Finite Difference Method - FDM) trên lưới tam giác:
Do các bó nhiên liệu của lò VVER có dạng hình lăng trụ lục giác, vùng hoạt được chia thành các lưới tam giác hai chiều (2D triangular mesh). Trong mô hình của chương trình LPO-V, mỗi bó nhiên liệu lục giác được chia thành 24 lưới tam giác bằng nhau. Tích phân phương trình khuếch tán trên thể tích $V_i$ của lưới tam giác thứ $i$, số hạng rò rỉ nơtron được chuyển đổi từ tích phân thể tích sang tích phân mặt của dòng nơtron $J_{g,i,j}$ qua các mặt giao diện $A_{i,j}$ với các lưới láng giềng. Dựa trên điều kiện liên tục dòng nơtron và thông lượng tại mặt tiếp xúc giữa hai lưới tam giác $i$ và $k$, dòng nơtron trung bình qua bề mặt được xác định bởi:
$$J_{g,i,k} = \frac{3 D_{g,i} D_{g,k}}{h (D_{g,i} + D_{g,k})} (\phi_{g,i} - \phi_{g,k}) = d_{g,i,k} (\phi_{g,i} - \phi_{g,k})$$
với $h$ là chiều cao của tam giác đều, diện tích bề mặt giao diện là $A_{i,k} = \frac{2h}{\sqrt{3}} = L$ (độ dài cạnh tam giác), và thể tích lưới là $V_i = 0,5 h L$.
-
Phương pháp giải hệ phương trình đại số:
Hệ phương trình đại số tuyến tính của các thông lượng nơtron trong toàn bộ lưới tam giác được giải bằng phương pháp quá dãn liên tiếp (Successive Over-Relaxation - SOR) để xác định giá trị riêng $k_{eff}$ và véctơ riêng phân bố thông lượng nơtron.
-
Nguồn tư liệu và chuẩn bị dữ liệu tiết diện nơtron:
Dữ liệu tiết diện nơtron được chuẩn bị từ bài toán chuẩn VVER-1000 MOX benchmark do OECD/NEA ban hành [5]. Phân tích tế bào mạng (lattice calculation) được thực hiện bằng mô-đun PIJ của hệ chương trình SRAC (Cơ quan Năng lượng Nguyên tử Nhật Bản - JAEA) sử dụng phương pháp xác suất va chạm (collision probability method). Thư viện dữ liệu hạt nhân sử dụng là ENDF/B-7.0 với cấu trúc 107 nhóm năng lượng (gồm 62 nhóm nhanh và 45 nhóm nhiệt). Bó nhiên liệu được mô phỏng hình học 1/6 đối xứng quay với 400 vùng không gian, mỗi thanh nhiên liệu chia làm 5 vùng đồng tâm. Bộ tiết diện 107 nhóm sau đó được thu gọn thành bộ hằng số vĩ mô 4 nhóm năng lượng (3 nhóm nơtron nhanh và 1 nhóm nơtron nhiệt) để nạp vào mã LPO-V.
-
Phát triển các phương pháp tối ưu hóa:
- Phương pháp Mô phỏng tôi kim tiến hóa (Evolutionary Simulated Annealing - ESA): Được cải tiến từ thuật toán SA nguyên bản bằng cách sử dụng các toán tử lai ghép (crossover) kết hợp đột biến giữa hai cấu hình cơ sở để tạo ra cấu hình thử nghiệm mới, thay vì chỉ hoán đổi vị trí nhị phân hoặc tam phân cục bộ. Hai toán tử lai ghép được phát triển gồm: Toán tử C1 (hoán đổi ngẫu nhiên 2 bó nhiên liệu giữa hai cấu hình cha mẹ) và Toán tử C2 (hoán đổi ngẫu nhiên một khối gồm nhiều bó nhiên liệu liền kề giữa hai cấu hình cha mẹ).
- Phương pháp Tiến hóa vi phân thích ứng dựa trên lịch sử thành công rời rạc (discrete SHADE): Thuật toán SHADE sử dụng một bộ nhớ lịch sử $H$ để tự động điều chỉnh hai tham số điều khiển là hệ số đột biến $F$ (theo phân bố Cauchy) và tỷ lệ trao đổi chéo $CR$ (theo phân bố chuẩn Gaussian). Nhờ cơ chế này, người dùng chỉ cần thiết lập hai tham số ban đầu là kích thước quần thể $NP$ và kích thước bộ nhớ $H$. Để áp dụng cho biến số nguyên của bài toán nạp nhiên liệu, tác giả tích hợp kỹ thuật đánh chỉ số vị trí tương đối (Relative Position Indexing - RPI) để chuyển đổi véctơ biến thực liên tục thành véctơ số nguyên chỉ định vị trí các bó nhiên liệu mà không làm thay đổi số lượng các loại bó nhiên liệu.
- Thiết lập hàm mục tiêu và kiểm định thống kê:
Hai hàm mục tiêu (Fitness Function - FF) được xây dựng:
- Hàm mục tiêu $F_1$: Dùng để kiểm tra và so sánh năng lực tái lập một cấu hình tham chiếu định trước giữa các thuật toán;
- Hàm mục tiêu $F_2$: Dùng cho bài toán tối ưu hóa thực tế, kết hợp tối đa hóa $k_{eff}$ và làm phẳng phân bố công suất (giảm $PPF$ và độ phân tán công suất $Flatness$):
$$F_2 = k_{eff} - w_p \cdot \max(0, PPF - PPF_{target}) - w_f \cdot Flatness$$
với các trọng số thành phần $w_p = 2,5$ và $w_f$ được xác định qua các tính toán khảo sát. Kiểm định phi tham số Mann-Whitney U Test ở mức ý nghĩa thống kê $\alpha = 0,05$ được thực hiện trên 50 lần chạy độc lập để đánh giá tính vượt trội giữa các thuật toán.
Nội dung chính theo từng chương
Chương 1: Giới thiệu chung và bài toán tối ưu hóa nạp tải nhiên liệu
Chương 1 trình bày tổng quan về lịch sử phát triển của công nghệ năng lượng hạt nhân thế giới từ những năm 1960 với các nhà máy thương mại đầu tiên như Dresden-1 (BWR), Yankee Rowe (PWR) tại Mỹ và Beloyarsk, Novovoronezh-1 (VVER) tại Nga. Tác giả phân tích các yếu tố kinh tế, sự cố hạt nhân lịch sử (Three Mile Island 1979, Chernobyl 1986, Fukushima 2011) và vai trò của các dòng lò phản ứng thế hệ III+ và thế hệ IV.
Chương này đi sâu vào cấu trúc và đặc tính kỹ thuật của các dòng lò VVER do Nga thiết kế với các đặc trưng khác biệt so với PWR của phương Tây: sử dụng bó nhiên liệu hình lục giác với các thanh nhiên liệu xếp theo ô mạng tam giác; sử dụng bình sinh hơi đặt nằm ngang; bình điều áp dung tích lớn; và không có các lỗ đâm xuyên qua đáy thùng áp lực lò. Bảng so sánh thông số kỹ thuật chính giữa các thế hệ VVER được trình bày chi tiết trong văn bản:
| Thông số kỹ thuật |
VVER-440 (V213) |
VVER-1000 (V-320) |
VVER-1200 (V-392M) |
| Công suất nhiệt (MWt) |
1375 |
3000 |
3212 |
| Công suất điện (MWe) |
440 |
1000 |
1198 |
| Mật độ công suất ($MW/m^3$) |
83 |
111 |
108,5 |
| Số vòng tuần hoàn sơ cấp |
6 |
4 |
4 |
| Áp suất mạch sơ cấp (MPa) |
12,27 |
15,7 |
16,2 |
| Lưu lượng chất tải nhiệt ($m^3/h$) |
41000 |
84800 |
86000 |
| Nhiệt độ lối vào / lối ra ($^\circ C$) |
270 / 300 |
290 / 320 |
298,2 / 328,9 |
| Số lượng bó nhiên liệu |
312 |
163 |
163 |
| Số thanh nhiên liệu trên mỗi bó |
126 |
312 |
312 |
| Vật liệu thanh hấp thụ |
Thép chứa Boron |
$B_4C$ |
$B_4C + Dy_2O_3 \cdot TiO_2$ |
| Độ làm giàu nhiên liệu ($wt%$) |
3,6 |
4,4 |
4,95 |
| Thời gian chu kỳ nhiên liệu (tháng) |
12 |
12 |
12 (18) |
| Chiều cao vùng hoạt (m) |
2,5 |
3,55 |
3,73 |
| Tuổi thọ thiết kế thùng lò (năm) |
30 |
40 |
60 |
Chương 1 cũng mô tả chi tiết bài toán ICFM và chứng minh bản chất đối nghịch giữa hai mục tiêu: để tăng chu kỳ nhiên liệu cần giảm rò rỉ nơtron ra biên ngoài, song điều này làm công suất dồn vào vùng trung tâm, khiến hệ số $PPF$ tăng vọt và vi phạm giới hạn an toàn ($PPF \le 1,35$ đối với lò PWR thông thường). Dữ liệu tính toán trên 3 cấu hình mẫu của vùng hoạt VVER-1000 MOX (Bảng 1.1) minh họa rõ ràng: Cấu hình (a) có $k_{eff} = 1,06197$ nhưng $PPF = 1,446$ (vi phạm an toàn); cấu hình (b) có $k_{eff} = 1,04859$ và $PPF = 1,316$; cấu hình (c) có $k_{eff} = 1,05763$ và $PPF = 1,384$.
Chương 2: Phương pháp nghiên cứu và phát triển công cụ
Chương 2 mô tả chi tiết bài toán chuẩn VVER-1000 MOX core benchmark của OECD/NEA [5]. Vùng hoạt gồm 163 bó nhiên liệu lục giác. Mô hình 1/6 đối xứng góc 60 độ chứa 28 vị trí bó nhiên liệu, gồm 19 bó $UO_2$ và 9 bó MOX. Bó $UO_2$ gồm 312 thanh nhiên liệu, 18 ống dẫn thanh điều khiển và 1 ống đo trung tâm; gồm 66 thanh có độ làm giàu 3,7 wt% $^{235}U$, 240 thanh có độ làm giàu 4,2 wt% $^{235}U$, và 6 thanh hấp thụ $UO_2-Gd_2O_3$ (UGD) chứa 3,6 wt% $^{235}U$ cùng 5 wt% $Gd_2O_3$. Bó MOX nạp 3 vùng nồng độ Plutoni với 234 thanh chứa 2,6 wt% Pu phân hạch và 18 thanh UGD. Các mức độ cháy và giá trị $k_{\infty}$ tương ứng của các loại bó nhiên liệu trong mô hình 1/6 được tổng hợp như sau:
| Loại bó |
Ký hiệu |
Độ cháy (MWd/kg) |
Số lượng bó trong 1/6 vùng hoạt |
Hệ số nhân vô hạn $k_{\infty}$ |
| $UO_2$ |
U0 |
0 |
5 |
1,30907 |
| $UO_2$ |
U15 |
15 |
5 |
1,13437 |
| $UO_2$ |
U32 |
32 |
5 |
1,01897 |
| $UO_2$ |
U40 |
40 |
4 |
0,97607 |
| MOX |
M0 |
0 |
3 |
1,18956 |
| MOX |
M17 |
17 |
3 |
1,07727 |
| MOX |
M33 |
33 |
3 |
0,99307 |
Tác giả trình bày thuật toán xây dựng chương trình mô phỏng vùng hoạt LPO-V bằng ngôn ngữ Fortran, giải hệ phương trình khuếch tán 4 nhóm năng lượng trên 24 lưới tam giác cho mỗi bó lục giác. Kiểm chứng chương trình LPO-V được thực hiện bằng cách so sánh với nghiệm tham chiếu Monte Carlo MCNP4c trên 5 trạng thái vận hành của bài toán chuẩn (S1: làm việc có boron; S2: nhiệt độ không đổi có boron; S3: trạng thái nguội boron cao; S4: làm việc không boron; S5: nhiệt độ không đổi không boron). Kết quả kiểm chứng chứng minh mã LPO-V đạt độ chính xác cao, độ sai lệch $k_{eff}$ rất nhỏ và phân bố công suất hoàn toàn tương đồng với MCNP4c, đồng thời tốc độ tính toán nhanh hơn mô-đun CITATION.
Chương 2 cũng trình bày chi tiết quy trình phát triển thuật toán ESA với 10 phương án tiếp cận kết hợp toán tử C1, C2 và đột biến, cùng thuật toán discrete SHADE với kỹ thuật biến đổi vector thực - nguyên RPI.
Chương 3: Tối ưu hóa cấu hình nạp nhiên liệu lò phản ứng VVER-1000
Chương 3 trình bày các kết quả tính toán số mô phỏng bài toán tối ưu hóa cấu hình nạp nhiên liệu cho vùng hoạt VVER-1000 nạp 30% MOX bằng hai phương pháp phát triển mới (ESA và discrete SHADE), đặt trong sự so sánh trực tiếp với SA, ASA và DE kinh điển:
- Khảo sát và tối ưu hóa bằng phương pháp ESA:
- Thực hiện khảo sát 10 biến thể ESA trên hai hàm mục tiêu $F_1$ và $F_2$ qua 50 lần chạy độc lập. Biến thể ESA kết hợp toán tử lai C1 và tỷ lệ đột biến thích hợp (ký hiệu ESA-C1A4) cho hiệu năng tốt nhất.
- Khi so sánh trên hàm mục tiêu tái lập cấu hình tham chiếu $F_1$, phương pháp ESA đạt độ hội tụ trung bình sau khoảng 3800 lượt đánh giá, nhanh hơn đáng kể so với ASA (khoảng 7200 lượt) và SA (vượt quá 10000 lượt).
- Đánh giá bằng kiểm định phi tham số Mann-Whitney U test xác nhận sự vượt trội về mặt thống kê của ESA so với SA và ASA ($p < 0,05$).
- Khảo sát và tối ưu hóa bằng phương pháp discrete SHADE:
- Khảo sát các tham số điều khiển ban đầu của SHADE: Sau khi quét các giá trị kích thước quần thể $NP$ từ 10 đến 40 và kích thước bộ nhớ lịch sử $H$ từ 4 đến 12, bộ tham số thích hợp nhất được xác định là $NP = 20$ và $H = 6$.
- Phân tích sự tiến hóa của hai tham số $F$ và $CR$ qua các thế hệ cho thấy cơ chế thích ứng lịch sử đã tự động điều chỉnh $F$ dao động quanh khoảng 0,5–0,7 và $CR$ tăng dần từ 0,5 lên xấp xỉ 0,8–0,9 ở các thế hệ sau nhằm tăng tốc độ hội tụ cục bộ.
- So sánh hiệu năng giữa SHADE, DE kinh điển và ESA qua 50 lần chạy độc lập cho thấy giá trị hàm mục tiêu tối ưu trung bình của SHADE đạt mức tương đương DE và ESA, vượt trội rõ rệt so với SA và ASA.
- Đánh giá cấu hình nạp nhiên liệu tối ưu tìm được:
Hai phương pháp ESA và discrete SHADE khi chạy độc lập đều tìm được một cấu hình nạp nhiên liệu tối ưu giống hệt nhau. Bảng so sánh giữa cấu hình tối ưu và cấu hình tham chiếu của bài toán chuẩn OECD/NEA được tổng kết như sau:
| Đại lượng vật lý |
Cấu hình tham chiếu (Reference LP) |
Cấu hình tối ưu (Optimal LP) |
Mức độ cải thiện |
| Hệ số nhân hiệu dụng ($k_{eff}$) tại BOC |
1,04071 |
1,05651 |
Tăng 1580 pcm |
| Hệ số đỉnh công suất hướng tâm ($PPF$) |
1,326 |
1,294 |
Giảm khoảng 2,4% |
| Độ không đồng đều công suất ($Flatness$) |
0,165 |
0,154 |
Giảm 6,7% |
Cấu hình tối ưu bố trí các bó nhiên liệu mới (độ cháy 0 MWd/kg) vào các vị trí xen kẽ hợp lý ở vùng giữa và đưa các bó nhiên liệu đã cháy sâu (như U40, M33) ra vùng biên ngoài để giảm rò rỉ nơtron, đồng thời tránh việc tập trung công suất quá mức tại tâm lò, giúp đường cong phân bố công suất xuyên tâm phẳng hơn đáng kể.
Chương 4: Kết luận và hướng nghiên cứu tiếp theo
Chương 4 đúc kết toàn bộ các kết quả đạt được của luận án về mặt phát triển phần mềm mô phỏng vật lý nơtron và hai thuật toán tối ưu hóa tổ hợp; đồng thời phân tích các giới hạn còn tồn tại và đề xuất các hướng phát triển kỹ thuật trong tương lai.
Kết quả và những đóng góp mới
Luận án tự xác định các đóng góp mới về mặt khoa học lý luận và thực tiễn kỹ thuật như sau:
- Đóng góp mới về mặt phương pháp và lý luận:
- Đã phát triển thành công phương pháp Mô phỏng tôi kim tiến hóa (ESA) kết hợp toán tử trao đổi chéo (crossover C1, C2) và đột biến kế thừa từ giải thuật di truyền vào thuật toán tôi kim, khắc phục hạn chế hội tụ chậm của thuật toán SA và ASA truyền thống;
- Đã xây dựng thành công biến thể rời rạc của thuật toán Tiến hóa vi phân thích ứng dựa trên lịch sử thành công (discrete SHADE) cho bài toán thay đảo nhiên liệu hạt nhân bằng cách tích hợp kỹ thuật đánh chỉ số vị trí tương đối (RPI) để ánh xạ không gian nghiệm biến thực liên tục sang biến nguyên hoán vị;
- Ứng dụng phương pháp kiểm định thống kê phi tham số Mann-Whitney U Test để đánh giá định lượng và khẳng định sự khác biệt có ý nghĩa thống kê về hiệu năng giữa các thuật toán tối ưu hóa ngẫu nhiên trong bài toán vật lý lò phản ứng.
- Đóng góp mới về mặt thực tiễn và ứng dụng:
- Xây dựng hoàn chỉnh chương trình tính toán vật lý nơtron LPO-V bằng ngôn ngữ Fortran, giải hệ phương trình khuếch tán đa nhóm dừng trên lưới tam giác dựa trên phương pháp FDM và SOR. Chương trình có độ chính xác tương đương nghiệm MCNP4c trên bài toán chuẩn VVER-1000 MOX và có hiệu năng tính toán nhanh hơn mô-đun CITATION;
- Tích hợp thành công mã LPO-V với hai thuật toán ESA và discrete SHADE tạo thành bộ công cụ tự động tìm kiếm cấu hình nạp đảo nhiên liệu cho lò VVER;
- Xác định được cấu hình nạp đảo nhiên liệu tối ưu cho vùng hoạt lò VVER-1000 nạp 30% MOX: giá trị $k_{eff}$ tại thời điểm đầu chu kỳ tăng khoảng 1580 pcm so với cấu hình tham chiếu (cho phép kéo dài chu kỳ nhiên liệu), đồng thời hệ số đỉnh công suất hướng tâm $PPF$ giảm khoảng 2,4% (tăng cường biên độ an toàn nhiệt).
Hạn chế và hướng nghiên cứu tiếp
Các hạn chế được văn bản luận án chỉ rõ bao gồm:
- Mô hình tính toán vật lý vùng hoạt trong luận án mới chỉ dừng lại ở mô hình hình học hai chiều (2D) đối xứng 1/6 góc quay 60 độ, chưa tính đến sự phân bố thông lượng nơtron và công suất theo phương thẳng đứng dọc theo chiều cao vùng hoạt (mô hình 3D);
- Quá trình tối ưu hóa cấu hình nạp tải được thực hiện tại thời điểm đầu chu kỳ (BOC) mà chưa thực hiện quá trình tính toán suy giảm độ cháy (burnup calculation) theo từng bước thời gian trong suốt chu kỳ vận hành;
- Ràng buộc an toàn nhiệt thủy lực mới chỉ xem xét đơn lẻ thông qua hệ số đỉnh công suất hướng tâm $PPF$ tại BOC, chưa xét đến các ràng buộc về hệ số phản ứng nhiệt độ hoặc độ sâu chèn của các chùm thanh điều khiển.
Hướng nghiên cứu phát triển tiếp theo được tác giả đề xuất:
- Nâng cấp chương trình tính toán LPO-V lên mô hình không gian ba chiều (3D) giải phương trình khuếch tán đa nhóm cho lò VVER;
- Tích hợp mô-đun tính cháy nhiên liệu phụ thuộc thời gian vào quy trình tìm kiếm tối ưu hóa tự động để theo dõi liên tục các chỉ số an toàn xuyên suốt chu kỳ vận hành;
- Mở rộng phạm vi ứng dụng của phương pháp ESA và discrete SHADE sang các bài toán nạp đảo nhiên liệu cho các lò phản ứng hạt nhân nghiên cứu có cấu trúc bất đối xứng hoặc các dòng lò thế hệ mới như lò nước siêu tới hạn VVER-SCWR.
Giá trị tham khảo
Luận án có giá trị tham khảo chuyên môn đối với các nhóm đối tượng cụ thể:
- Các nghiên cứu sinh, học viên cao học và giảng viên chuyên ngành Vật lý nguyên tử, Kỹ thuật hạt nhân: Bản mô tả chi tiết phương pháp sai phân hữu hạn giải phương trình khuếch tán trên lưới tam giác hai chiều (Mục 2.4) và quy trình chuẩn bị dữ liệu tiết diện vĩ mô 4 nhóm bằng hệ chương trình SRAC/PIJ kết hợp thư viện ENDF/B-7.0 (Mục 2.3) là tài liệu tham khảo có giá trị cao trong nghiên cứu tính toán vùng hoạt;
- Các kỹ sư thiết kế và vận hành nhà máy điện hạt nhân: Cấu hình nạp nhiên liệu tối ưu, phân bố công suất xuyên tâm và các thông số đánh đổi giữa $k_{eff}$ và $PPF$ đối với vùng hoạt VVER-1000 nạp 30% MOX (Mục 3.3.4) cung cấp cơ sở dữ liệu thực nghiệm phục vụ việc hoạch định chiến lược quản lý nhiên liệu trong lò;
- Các nhà nghiên cứu trong lĩnh vực tối ưu hóa ứng dụng: Kỹ thuật Relative Position Indexing (RPI) để chuyển đổi thuật toán SHADE từ miền biến liên tục sang miền biến rời rạc tổ hợp (Mục 2.6.4), cùng phương pháp luận đánh giá thuật toán ngẫu nhiên bằng kiểm định Mann-Whitney U Test (Mục 2.8), là những giải pháp kỹ thuật có thể tham khảo và áp dụng cho các bài toán tối ưu hóa tổ hợp phức tạp khác.
Câu hỏi thường gặp
1. Vùng hoạt chuẩn VVER-1000 MOX benchmark được cấu hình như thế nào về tỷ lệ nhiên liệu và các mức độ cháy?
Vùng hoạt bài toán chuẩn của OECD/NEA gồm 163 bó nhiên liệu lục giác, nạp 30% nhiên liệu MOX và 70% nhiên liệu $UO_2$. Trong mô hình 1/6 vùng hoạt (gồm 28 bó), có 19 bó $UO_2$ phân chia theo 4 mức độ cháy: 0 MWd/kg (5 bó), 15 MWd/kg (5 bó), 32 MWd/kg (5 bó), 40 MWd/kg (4 bó); và 9 bó MOX phân chia theo 3 mức độ cháy: 0 MWd/kg (3 bó), 17 MWd/kg (3 bó), 33 MWd/kg (3 bó).
2. Chương trình LPO-V sử dụng phương pháp số nào để giải phương trình khuếch tán nơtron trong vùng hoạt lò VVER?
Chương trình LPO-V sử dụng phương pháp sai phân hữu hạn (FDM) trên lưới tam giác hai chiều, chia mỗi bó nhiên liệu lục giác thành 24 lưới tam giác bằng nhau để giải hệ phương trình khuếch tán đa nhóm dừng (4 nhóm năng lượng), và sử dụng phương pháp quá dãn liên tiếp (SOR) để giải hệ phương trình đại số tuyến tính tìm trị riêng $k_{eff}$ cùng phân bố thông lượng nơtron.
3. Kỹ thuật Relative Position Indexing (RPI) trong phương pháp discrete SHADE giải quyết bài toán biến số nguyên như thế nào?
Kỹ thuật RPI chuyển đổi một véctơ số thực liên tục thành một véctơ số nguyên thông qua việc sắp xếp thứ tự độ lớn của các phần tử thực để tạo ra thứ tự chỉ số vị trí (vị trí tương đối). Thứ tự chỉ số này được dùng để gán các loại bó nhiên liệu vào các vị trí trong vùng hoạt, giúp thuật toán SHADE tối ưu hóa bình thường trên miền thực mà luôn ánh xạ chính xác về một cấu hình bó nhiên liệu rời rạc hợp lệ.
4. Phương pháp ESA cải tiến cơ chế sinh nghiệm thử nghiệm như thế nào so với SA và ASA?
Thay vì chỉ sử dụng phép hoán vị nhị phân hoặc tam phân ngẫu nhiên giữa các bó nhiên liệu trong một cấu hình đơn lẻ như SA và ASA, phương pháp ESA sử dụng hai toán tử lai ghép trao đổi chéo giữa hai cấu hình cha mẹ: Toán tử C1 (hoán đổi ngẫu nhiên 2 bó nhiên liệu) và Toán tử C2 (hoán đổi ngẫu nhiên một khối bó nhiên liệu), kết hợp toán tử đột biến ngẫu nhiên, giúp mở rộng không gian tìm kiếm và tăng tốc độ hội tụ.
5. Cấu hình nạp nhiên liệu tối ưu tìm được mang lại những cải thiện cụ thể nào so với cấu hình tham chiếu?
Cấu hình tối ưu được tìm thấy đồng thời bởi cả hai phương pháp ESA và discrete SHADE giúp giá trị $k_{eff}$ tại thời điểm đầu chu kỳ tăng khoảng 1580 pcm (từ 1,04071 lên 1,05651), đồng thời giảm hệ số đỉnh công suất hướng tâm $PPF$ khoảng 2,4% (từ 1,326 xuống 1,294) và giảm độ không đồng đều công suất $Flatness$ khoảng 6,7% (từ 0,165 xuống 0,154).
Kết luận
Luận án tiến sĩ của tác giả Trần Việt Phú đã giải quyết bài toán tối ưu hóa cấu hình nạp đảo nhiên liệu cho lò phản ứng hạt nhân VVER-1000 nạp 30% nhiên liệu MOX thông qua việc xây dựng chương trình mô phỏng vật lý nơtron LPO-V và phát triển hai phương pháp tối ưu hóa hiện đại là ESA và discrete SHADE. Kết quả tính toán số khẳng định chương trình LPO-V đạt độ chính xác tương đương MCNP4c với tốc độ giải nhanh, đồng thời hai thuật toán tối ưu hóa đề xuất đã chứng minh được tính ưu việt thống kê so với SA và ASA qua kiểm định Mann-Whitney U test. Cấu hình tối ưu tìm được giúp tăng $k_{eff}$ thêm 1580 pcm và giảm hệ số đỉnh công suất hướng tâm 2,4%, góp phần hài hòa giữa hiệu quả kinh tế chu kỳ nhiên liệu và giới hạn an toàn nhiệt của vùng hoạt lò phản ứng hạt nhân.
