Computational Methods for Burnup Calculations with Monte Carlo Neutronics

The work described in this thesis deals with the computational methods and algorithms used in burnup calculations, which model changes in the composition of nuclear fuel under irradiation. While only cases where the neutron transport part of the calculation is handled by the Monte Carlo method are considered, most of the results should also be applicable with deterministic neutronics. During each step of a Monte Carlo burnup calculation, changes in material compositions are solved by evaluating an explicit solution to the Bateman equations with constant coefficients. Five depletion algorithms capable of doing this while explicitly modeling all of the thousands of nuclides and reactions encountered in burnup calculations were compared. The results are quite conclusive and, together with other studies, show rational approximation based matrix exponential methods to be the best choice for Monte Carlo burnup calculations. The constant coefficients of the Bateman equations are selected by a coupling scheme that uses one or more steady state neutronics solutions to predict their time development. Because the coefficients must be constant, these predictions are further approximated with their averages. New coupling schemes that use data from the previous step to make higher order predictions are presented. Since the old values are readily available, no additional calculations are required, and the stepwise running time is not affected. The coupling is further improved by dividing the steps to substeps, which are then solved sequentially. Since each substep can use different coefficients for the Bateman equations, this allows piecewise constant, rather than constant, approximation of the predicted behavior. These new methods greatly improve the accuracy obtainable with given step lengths, thus allowing longer steps to be used. Prior studies have shown that the existing coupling schemes used in Monte Carlo burnup calculations suffer from instabilities caused by spatial xenon oscillations. The new methods are also affected, but it is shown that the simulation models used in these tests actually describe physical xenon oscillations, not a stable state. Thus it is the models, not the methods used to solve them, that are unstable. Regardless, all xenon driven oscillations can be prevented by forcing a mutual equilibrium between the neutron flux and saturated xenon distribution. The equilibrium calculation can be integrated to Monte Carlo neutronics, which provides a simple and lightweight solution that can be used with any of the existing burnup calculation algorithms. However, oscillations driven by nuclides other than xenon may still arise if step lengths are too long.Tässä väitöskirjassa esitelty tutkimus on keskittynyt laskennallisiin menetelmiin ja algoritmeihin, joita käytetään ydinpolttoaineen koostumuksen käytönaikaisia muutoksia mallintavissa palamalaskuissa. Vaikka työssä käsitelläänkin vain tilanteita, joissa laskujen neutroniikka-osio hoidetaan Monte Carlo -menetelmällä, ovat useimmat tulokset sovellettavissa myös deterministisen neutroniikan kanssa. Monte Carlo -palamalasku muodostuu askelista, joilla muutokset materiaalien koostumuksissa lasketaan evaluoimalla eksplisiittinen ratkaisu vakiokertoimisiin Bateman-yhtälöihin. Työn ensimmäisessä osassa vertailtiin viittä algoritmia, jotka pystyvät tekemään tämän huomioiden eksplisiittisesti kaikki ne tuhannet nuklidit ja reaktiot, joihin palamalaskuissa törmätään. Vertailun tulokset ovat selkeitä ja yhdessä muiden tutkimusten kanssa osoittavat rationaaliapproksimaatioihin pohjautuvat matriisieksponentiaalimenetelmät Monte Carlo -palamalaskuihin parhaiten soveltuviksi. Bateman-yhtälöiden kerrointen arvot valitaan palama-algoritmilla, joka käyttää yhtä tai useampaa aikariippumatonta neutroniikkaratkaisua ennustaakseen niiden kehityksen. Koska kertoimet halutaan vakioiksi, approksimoidaan näitä ennusteita edelleen niiden keskiarvoilla. Väitöksessä esitetään uusia palama-algoritmeja, jotka käyttävät informaatiota edellisiltä askelilta muodostaakseen korkeamman asteen ennusteita. Koska edellisten askelten arvot ovat valmiiksi saatavilla, tämä ei vaadi ylimääräisiä laskutoimituksia, eikä siten vaikuta askelkohtaiseen ajoaikaan. Tarkkuutta voidaan edelleen parantaa jakamalla askeleet aliaskeleisiin, jotka ratkaistaan peräkkäin. Koska jokainen aliaskel voi käyttää Bateman-yhtälöissä eri kertoimia, tämä mahdollistaa kerrointen ennustetun kehityksen tarkemman seuraamisen. Nämä uudet menetelmät parantavat suuresti saavutettavissa olevaa tarkkuutta mahdollistaen askelpituuksien lisäämisen. Aikaisemmat tutkimukset ovat osoittaneet vanhojen palama-algoritmien kärsivän spatiaalisista ksenon-oskillaatioista, jotka vaikuttavat myös uusiin menetelmiin. Työssä kuitenkin osoitetaan epästabiileissa testeissä käytettyjen mallien kuvaavan fysikaalisia ksenon-oskillaatioita. Oskillaatioiden taustalla on siis ratkaistavien tapausten, ei mallintamisessa käytettyjen menetelmien, epästabiilisuus. Kaikki ksenon-lähtöiset oskillaatiot voidaan estää pakottamalla saturoitunut ksenon-konsentraatio ja neutronivuo tasapainoon. Tasapainotilan laskeminen voidaan yhdistää Monte Carlo -neutroniikan laskuun, mikä tarjoaa kaikkien palama-algoritmien kanssa toimivan yksinkertaisen ja kevyen ratkaisun. Liian pitkiä askelia käytettäessä voi kuitenkin edelleen esiintyä muiden nuklidien kuin ksenonin ajamia oskillaatioita

Improving the accuracy of the Chebyshev rational approximation method using substeps

The Chebyshev rational approximation method (CRAM) for solving the decay and depletion of nuclides is shown to have a remarkable decrease in error when advancing the system with the same time step and microscopic reaction rates as the previous step. This property is exploited here to achieve high accuracy in any end-of-step solution by dividing a step into equidistant substeps. The computational cost of identical substeps can be reduced significantly below that of an equal number of regular steps, as the lower-upper decompositions for the linear solutions required in CRAM need to be formed only on the first substep. The improved accuracy provided by substeps is most relevant in decay calculations, where there have previously been concerns about the accuracy and generality of CRAM. With substeps, CRAM can solve any decay or depletion problem with constant microscopic reaction rates to an extremely high accuracy for all nuclides with concentrations above an arbitrary limit