Development of Advanced Time-dependent Monte Carlo Neutron Transport Methods for Reactor Transient Analysis

Abstract

The analysis of the time-dependent neutron behavior is essential to understanding the kinetic properties of reactor cores in various application fields including reactor start-up analyses, reactivity measurements, accident analysis, and experiments in research reactors. Since the reactor transient analysis requires a lot of computational time compared to the steady-state analysis, it has mainly relied on deterministic or quasi-static methods. However, these methods can lead to inaccurate results due to inherent differential approximations to space, energy, or time domains. Recently, thanks to the ever-advancing computing power and the development of high-reliability methodologies, the time-dependent Monte Carlo (TDMC) neutron transport method without any approximation is being actively studied and has become an applicable alternative for reactor transient analysis in a practical time range. The development of Monte Carlo transient analysis method is very important because it can provide reference solutions for the transient analysis of various reactors such as gen-IV reactors and research reactors as well as commercial reactors. In domestic, McCARD, a Monte Carlo code developed by Seoul National University, has equipped with the capability of the TDMC simulation and performed alpha eigenvalue calculations and transient analysis. However, the existing code is limited to two-dimensional analysis and lacks the capability for an accurate and reliable reactor transient analysis such as an unbiased variance estimation. Therefore, motivated by the necessity of developing such high-reliability transient analysis code, this thesis aims to advance the TDMC algorithms for reactor transient analysis and to extend the applicability of the TDMC simulation to practical problems. The TDMC method enables sustainable neutron simulation by introducing time intervals to the conventional Monte Carlo simulation, conducting time bin-by-bin neutron simulation, and then controlling the neutron population at the end of each time interval. In McCARD, the analog MC branching method that simulates the branching of particles as it is, and the combing method, which samples the exact number of neutrons without bias, are used as neutron population control methods. The efficiency of the algorithm is improved by removing the scale factor previously used for weight normalization and introducing a dynamic weight window. The TDMC method which allows sustainable simulation of neutrons over time is applied to the PNS experiments to estimate the prompt neutron decay constant alpha. In the PNS experiment, it has been reported that the different alpha values are measured depending on the detector position and detection time due to the initial source and geometry effect. Utilizing the TDMC method which can accurately simulate the space- and time-dependent behavior of neutrons, an optimum detector position search algorithm is developed for the PNS alpha measurement. The developed method is applied to the Pb-Bi-zoned experimental benchmark at KUCA and numerical experiments at AGN-201K. In the KUCA experimental benchmark, the spallation source effect is well reflected so the relative alpha convergence at different detector positions is well predicted. As for the AGN-201K simulation, the sensitivity of the optimum detector position to the initial source location is evaluated. By comparing the detector signals at convergence times according to positions, the optimum detector position is determined where the detector signal is highest. The application results are expected to be a good reference for designing and performing actual PNS alpha measurement experiments. Along with the neutron population control method, the methodologically important parts of the reactor transient simulation are the delayed neutron simulation and the steady-state modeling. Since the prompt neutron and the delayed neutron have a large difference in generation time, the conventional method of directly sampling the delayed neutron from the fission event causes a large statistical error. In addition, a method that can simulate the distribution of the prompt neutron and delayed neutron sources in the initial steady-state is needed because most transient analysis starts from the steady-state. Therefore, the forced decay algorithm through precursor simulation is introduced for efficient delayed neutron simulation. As the initial steady-state modeling method, the TDMC steady-state simulation method is employed, which consistently uses the TDMC simulation method from the initial state to the transient state without changing the calculation mode. In this process, the algorithm for normalizing the initial number of precursors is improved, and a new feature of moving geometry treatment to deal with more realistic three-dimensional transient scenarios is developed. The developed McCARD transient analysis capability is verified for the 2D and 3D problems of the C5G7-TD reactor transient benchmark and compared with nTRACER. In the 3D problems, the axial insertion and withdrawal of the control rods are well simulated, and the trends of the core dynamic reactivity and relative fission rate show good agreement with the results of nTRACER within the stochastic errors. The Monte Carlo calculation usually provides the accuracy and reliability of the calculation result through the sample mean and the variance of the sample mean. In the course of verifying the capability of McCARD transient analysis modules, it is found that the sample variance of a TDMC tally mean obtained from the conventional stochastic processing is highly biased. This is because correlation is made between neutrons during the branching process or population control of the TDMC calculation. The biased variance gives distorted information for judging the accuracy and reliability of Monte Carlo calculation results. Moreover, there is a problem of how to allocate and statistically process the contribution of the delayed neutrons when estimating the variance of tally mean. To address these issues and to estimate accurate variance, a history-based batch method for the TDMC simulation is developed. In the history-based batch method, the neutrons and precursors are grouped in several batches to simulate separately, and the results are statistically processed batch-wisely to break the correlation between the estimates. At the same time, the allocation problem of the delayed neutron contributions can be solved naturally by assigning the contribution to the included batch tally. The developed method is verified in infinite homogeneous two-group problems and C5G7-TD benchmark problems. It gives unbiased variances for the tally means if the batch size is sufficient. In addition, error propagation is observed in the TDMC simulation. It is demonstrated that it stems from the weight normalization scheme in the population control and propagates through the weight of survival neutrons. On the other hand, in the system in which the delayed neutrons are more dominant than the survival neutrons, the error does not propagate significantly because the neutrons are not survived for a long time interval and the delayed neutrons contribution is relatively large. When establishing point kinetics model for reactor transient analysis or performing a reactivity measurement experiment, it is necessary to calculate kinetics parameters of the core. In general, since it is difficult to obtain an accurate neutron flux distribution in a transient state, kinetics parameters are calculated by assuming solutions of the steady-state transport equation and its adjoint equation as the shape function and weighting function. However, as it becomes possible to simulate the accurate time-dependent neutron flux distribution through the TDMC method, a time-dependent kinetics parameters estimation method is developed based on the exact point kinetics equations utilizing TDMC simulation. In particular, Monte Carlo algorithms that can efficiently calculate the adjoint response during TDMC forward simulation are developed instead of the conventional Contributon method, which is computationally quite burdensome. To verify the developed methods, the time-dependent kinetics parameters are evaluated in infinite homogeneous two-group problems, and the results show good agreement with the analytic solutions. In addition, it shows more than 1,000 times the computational efficiency for problems near the critical when comparing with the Contributon method. Then to check the applicability of the estimated time-dependent kinetics parameters, the point kinetics equation is established using them and the transient behavior of the system is predicted. The point kinetics analyses are performed on the C5G7-TD problem in which the initial state is critical and the beam trip simulation of the thorium-loaded ADS system at KUCA in which the initial state is subcritical. The PKEs with different kinetics parameters estimated from the k-eigenvalue calculation, fixed source calculation, and TDMC calculation are compared to each other. It is noted that the developed method can calculate the kinetics parameters reflecting the accurate neutron flux distribution regardless of the system and it provides a framework for the point kinetics analysis in the generalized time domain including the conventional steady-state and arbitrary transient states. Finally, McCARD/CUPID coupled transient analysis system is established using TCP/IP socket communication to provide transient analysis capabilities for more practical problems considering the thermal-hydraulic feedback. CUPID is a three-dimensional sub-channel code developed by KEARI, and by coupling it with McCARD, it is possible to conduct transient analysis considering the effect of coolant mixing in the sub-channels. The coupled analysis system is verified in the VERA #6 HFP assembly problem by comparing it with the existing McCARD/MATRA coupled steady-state analysis system. The steady-state results match well with the McCARD/MATRA system with a maximum difference of 0.17% in power distribution, 0.06% in coolant exit temperature, and 1.55% in fuel temperature. As for the coupled transient analysis, a simple rod ejection accident is analyzed in the modified VERA #6 problem under the HZP condition. From the temperature and power trend, the Doppler effect of the fuel temperature feedback is observed, and thus the integrity of the coupled transient analysis system is verified.

시간에 따른 중성자의 거동을 해석하는 것은 노심 기동 해석, 임계도 측정, 사고 해석, 연구용 원자로 실험 등 다양한 적용분야에서 노심의 동역학적 특성을 이해하는데 필수적이다. 원자로 과도해석은 정상상태 해석에 비해 많은 전산 시간이 필요하기 때문에 주로 결정론적 방법론이나 준정상상태 방법론에 의존해 왔으나 이 방법론들은 공간과 에너지 또는 시간에 대한 차분화 근사를 포함하기 때문에 부정확한 결과를 초래할 수 있다. 최근 지속적인 전산 성능의 향상과 고신뢰도 방법론들의 개발에 힘입어 근사가 없는 시간종속 몬테칼로(TDMC) 중성자 수송해석방법론에 대한 연구가 활발히 이루어지고 있으며 실제 시간 범위의 원자로 과도해석에 적용 가능한 대안이 되었다. 몬테칼로 방법론은 상용 원자로뿐만 아니라 차세대 원자로나 연구로 등 다양한 노심의 과도해석에 기준해를 제공할 수 있기 때문에 몬테칼로 과도해석 방법론의 개발은 매우 중요하다. 국내에서는 서울대학교에서 개발한 몬테칼로 코드인 McCARD가 TDMC 모의 기능을 갖추어 알파 고유치 계산과 과도해석을 수행한 적이 있으나 2차원 해석에 제한되고 편향된 분산을 예측하는 등 정밀하고 신뢰할만한 원자로 과도해석능을 갖추지 못하였다. 본 논문은 고신뢰도 과도해석 코드 개발의 필요성을 바탕으로 원자로 과도해석을 위한 TDMC 알고리즘을 고도화 개발하고 TDMC 모의의 적용성을 실질적인 문제에 확장하는 것을 목표로 한다. TDMC 방법론은 기존의 몬테칼로 모의계산에 시간 구간을 도입하고 시간 구간별로 중성자 수송을 모의한 다음 각 시간 구간의 끝에서 중성자수를 제어함으로써 중성자의 지속적인 모의를 가능하게 한다. McCARD에서는 입자의 분기를 그대로 모의하는 아날로그 몬테칼로 분기법과 편향 없이 정확한 수의 중성자를 추출하는 빗질 방법을 중성자수 제어 방법으로 사용하였으며, 기존에 가중치 정규화 용도로 사용하던 축적 인자를 제거하고 동적 가중치창을 도입하여 알고리즘의 효율성을 개선하였다. 시간에 따른 중성자의 지속적인 모의가 가능한 TDMC 방법은 즉발중성자붕괴상수 알파를 계측하는 펄스중성자선원(PNS) 실험에 적용하였다. PNS 실험은 초기 선원과 노심의 기하학적 영향으로 검출기의 위치와 시간에 따라 서로 다른 결과가 측정되는 문제가 보고되어 왔는데 위치와 시간에 따른 중성자의 거동을 정확히 모의할 수 있는 TDMC 방법을 활용하여 PNS 알파 측정을 위한 최적의 검출기 위치를 찾는 탐색 알고리즘을 개발하였다. 개발된 방법은 교토대 임계집합체(KUCA)의 납-비스무스 장전 실험 검증 문제와 AGN-201K의 모의 실험에 적용하였다. KUCA 실험 검증문제에서는 파쇄 중성자 선원의 영향이 잘 반영되어 서로 다른 검출기 위치에서의 알파 수렴성을 잘 예측하였으며 AGN-201K 모의 실험에서는 초기 선원 위치에 따른 계측의 민감도 평가가 수행되었다. 위치에 따른 수렴 시간에서 검출기 신호를 비교하여 신호가 가장 높은 곳을 최적의 검출기 위치로 선정하였고 이러한 적용 해석 결과는 실제 PNS 알파 계측 실험을 설계하고 수행하는데 있어 좋은 참고가 될 것으로 기대된다. 중성자수 제어 방법과 더불어 원자로 과도상태 모의에 있어 방법론적으로 중요한 부분이 지발중성자 모의법과 정상상태 모델링 방법이다. 즉발중성자와 지발중성자는 세대 시간에 큰 차이가 있기 때문에 핵분열로부터 지발중성자를 직접 추출하는 통상적인 방법은 큰 통계 오차를 일으킨다. 또한 대부분의 과도해석은 정상상태로부터 시작하기 때문에 초기 정상상태의 즉발중성자와 지발중성자 선원 분포를 모의할 수 있는 방법이 필요하다. 따라서 효율적인 지발중성자 모의를 위해 지발중성자 선행핵 모의를 통한 강제 붕괴 알고리즘을 도입하였고, 초기상태 모델링 방법으로는 계산 모드의 전환 없이 초기부터 과도상태까지 TDMC 모의법을 통해 일관적으로 초기 정상상태를 모델링할 수 있는 몬테칼로 정상상태 모의법을 사용하였다. 이 과정에서 초기 지발중성자 선행핵수를 정규화하는 알고리즘을 개선하였고 보다 실질적인 3차원 과도상태 시나리오를 모의할 수 있는 동적 구조물 처리 기능을 개발하였다. 개발된 McCARD 과도해석능은 원자로 과도상태 검증 계산 문제집인 C5G7-TD의 2차원 문제와 3차원 문제에 대해 nTRACER와 비교 검증하였다. 3차원 문제에서 제어봉의 축방향 삽입과 인출이 잘 모의되었으며 시간에 따른 노심 동적반응도 및 상대 핵분열 반응률은 통계 오차 이내에서 잘 일치하였다. 몬테칼로 계산은 일반적으로 표본 평균과 표본 평균의 분산을 통해 계산 결과의 정확도와 신뢰도를 제공하는데 McCARD 과도해석 모듈의 해석능을 검증하는 과정에서 일반적인 통계처리를 통해 얻은 TDMC 집계 평균의 표본 분산이 크게 편향되어 있는 것을 발견하였다. 이는 TDMC 계산의 분기 과정 및 중성자수 제어에서 중성자간에 상관관계가 생기기 때문인데 이러한 편향된 분산은 몬테칼로 계산 결과의 정확도 및 신뢰도를 판단하는데 왜곡된 정보를 주게 된다. 또한 집계 평균의 분산을 추정할 때 지발중성자의 기여도를 어떻게 할당하여 통계 처리해야 하는지에 대한 문제도 존재한다. 본 논문에서는 이러한 문제를 해결하고 정확한 분산을 계산하기 위해 TDMC 모의에서의 히스토리 기반 배치법을 개발하였다. 히스토리 기반 배치법은 중성자와 선행핵을 묶어 배치별로 따로 모의하고 결과를 배치별로 통계처리하여 추정치 간 상관관계를 끊는 방법으로 동시에 지발중성자의 기여도 할당 문제도 기여도를 포함된 배치의 집계에 할당함으로써 자연스럽게 해결할 수 있다. 개발된 방법은 2군 무한균질문제와 C5G7-TD 검증 계산문제를 통해 검증하였고 배치의 크기가 충분할 경우 집계 평균에 대해 편향되지 않은 정확한 분산을 계산하였다. 또한, 검증과정에서 TDMC 모의 도중 오차가 전파되는 현상을 확인하였고 이것이 중성자수 제어의 가중치 정규화로부터 파생되어 살아남은 중성자들을 통해 전파된다는 것을 입증하였다. 한편 살아남은 중성자보다 지발중성자의 영향이 지배적인 시스템에서는 중성자가 오래 살아남지 못하고 지발중성자의 기여가 크기 때문에 오차가 크게 전파되지 않았다. 원자로 과도 해석을 위해 일점운동방정식 모델을 수립하거나 반응도 계측 실험을 수행할 때는 해당 노심의 동특성 인자 계산이 필요하다. 일반적으로는 과도상태에서 정확한 중성자속 분포를 얻기 어렵기 때문에 정상상태 수송방정식과 수반방정식의 해를 형상 함수와 가중 함수로 가정하고 동특성 인자를 계산한다. 그러나 TDMC 방법을 통해 시간에 따른 정확한 중성자속 분포를 모의할 수 있게 되면서 TDMC 모의를 이용하여 정확한 일점운동방정식에 기반한 시간종속 동특성 인자 계산법을 개발하였다. 특히 수반 응답 계산에 대해서는 계산 부담이 큰 기존의 기여자 모의법(Contributon method)을 대신하여 TDMC 순방향 모의 중에 이를 효율적으로 계산할 수 있는 몬테칼로 알고리즘을 개발하였다. 개발된 방법은 검증을 위해 2군 무한균질 문제에서 시간에 따른 동특성 인자를 평가하였으며 결과는 해석해와 잘 일치하였다. 또한 기존의 기여자 모의법과 성능을 비교했을 때 임계에 가까운 문제에서는 1,000배가 넘는 계산 효율을 보였다. 다음으로 계산된 시간종속 동특성 인자의 적용성을 확인하기 위해 이를 이용하여 일점운동방정식을 세우고 시스템의 과도상태 거동을 예측하였다. 일점운동해석은 초기 상태가 임계인 C5G7-TD 문제와 초기 상태가 미임계인 KUCA의 토륨장전 가속기 구동 시스템의 빔 트립 모의 실험에 대해 수행되었으며 k 고유치 계산과 고정 선원 계산, TDMC 계산으로 추정된 동특성 인자를 이용하여 서로 다른 일점운동방정식을 비교하였다. 개발된 방법은 시스템에 상관없이 정확한 중성자속 분포를 반영하여 동특성 인자를 계산할 수 있으며 기존의 정상상태를 기준으로 수행되던 일점운동해석을 일반화하여 임의의 상태를 기준으로 일점운동해석을 수행할 수 있는 해석 체계를 제공한다는 점에서 큰 의의가 있다. 마지막으로 열수력 궤환 효과를 고려한 보다 실제적인 문제에 대한 과도해석능을 갖추기 위해 TCP/IP 소켓 통신을 이용하여 McCARD/CUPID 연계 과도해석 체계를 구축하였다. CUPID는 한국원자력연구원에서 개발한 3차원 부수로 코드로 이를 연계함으로써 부수로에서의 냉각재 섞임 효과를 고려한 과도상태 해석을 수행할 수 있다. 연계해석체계는 VERA 6번 고온전출력 집합체 문제에 대해 기존에 정상상태 해석을 위해 구축하였던 McCARD/MATRA 연계해석체계와 비교 검증하였으며, 정상상태에서 출력 분포는 0.17%, 냉각재 출구 온도는 0.06%, 핵연료 온도는 1.55%의 차이 이내에서 서로 잘 일치하였다. 과도 상태에서의 연계해석으로는 VERA 6번 문제를 수정하여 고온영출력 조건에서 간단한 제어봉 이탈사고를 가정하여 이를 해석하였다. 핵연료 온도와 출력변화로부터 핵연료 온도 궤환에 의한 도플러 효과를 확인하였고 이를 통해 연계해석체계의 건전성을 확인하였다.