Static and Dynamic Modelling of 1D Heterogeneous Molten Salt Reactors Development of a Numerical Solver for Stationary and Fre- quency Domain Dynamic Two-Group Neutron Diffusion Prob- lems Master’s thesis in Physics ERIK ANDERSSON DEPARTMENT OF PHYSICS CHALMERS UNIVERSITY OF TECHNOLOGY Gothenburg, Sweden 2021 www.chalmers.se www.chalmers.se Master’s thesis 2021 Static and Dynamic Modelling of 1D Heterogeneous Molten Salt Reactors Development of a Numerical Solver for Stationary and Frequency Domain Dynamic Two-Group Neutron Diffusion Problems ERIK ANDERSSON Department of Physics Division of Subatomic, High Energy and Plasma Physics Task Force on Deterministic REActor Modelling (DREAM) Chalmers University of Technology Gothenburg, Sweden 2021 Static and Dynamic Modelling of 1D Heterogeneous Molten Salt Reactors Development of a Numerical Solver for Stationary and Frequency Domain Dynamic Two-Group Neutron Diffusion Problems ERIK ANDERSSON © ERIK ANDERSSON, 2021. Examiner: Christophe Demazière, Department of Physics Master’s Thesis 2021 Department of Physics Division of Subatomic, High Energy and Plasma Physics Task Force on Deterministic REActor Modelling (DREAM) Chalmers University of Technology SE-412 96 Gothenburg Telephone +46 31 772 1000 Typeset in LATEX, template by Magnus Gustaver Printed by Chalmers Reproservice Gothenburg, Sweden 2021 iv Static and Dynamic Modelling of 1D Heterogeneous Molten Salt Reactors Development of a Numerical Solver for Stationary and Frequency Domain Dynamic Two-Group Neutron Diffusion Problems ERIK ANDERSSON Department of Physics Chalmers University of Technology Abstract Molten salt reactors represent one of the so called Generation IV concepts. In this design the fuel is dissolved in a molten salt, which also acts as a coolant. The moving fuel introduces a new physical phenomenon not present in reactor designs with stationary fuel. Delayed neutrons, emitted through the decay of certain fission fragment nuclei, called precursors, will in this case be emitted down-stream from where the original fission event occurred. Since the precursor nuclei have mean lives comparable to the temporal scale of the velocity, this will affect the neutron flux in the reactor. There is also the possibility of precursors reentering the core since the molten salt is circulating in the system. This thesis describes the process of developing a novel numerical simulation tool, using the diffusion approximation for neutron flux, for the calculation of the neutron flux and precursor concentration of a one dimensional heterogeneous reactor core. The tool solves the diffusion equa- tions with two energy groups and one group of precursors in the stationary case and for a monochromatic perturbation in the frequency domain, respectively. The eigenvalue problem arising in the stationary case with no external sources is solved using the power iteration method with a Wielandt’s shift technique and the inver- sion of matrices is done using an LUPQ factorization. For the purpose of testing the implementation, a scaled down solver for one energy group and one group of precursors was developed in the same way as the full solver. When tested, the full solver reproduces the results of a previously developed similar solver for the case of zero fuel velocity and the scaled down version reproduces the results of a pre- viously developed one-group solver as well as the analytical solution to the special case of infinite fuel velocity. A condition on the mesh resolution was found that greatly impacted the dynamic solution for the precursor concentration as its spatial distribution exhibits an oscillatory behaviour. Keywords: Molten Salt Reactors, two-group neutron diffusion theory, neutron noise, numerical reactor modelling v Acknowledgements I would like to convey my sincerest gratitude to all the people who have helped and supported me during this project, as well as during my time at Chalmers Institute of Technology. I would like to begin by thanking my supervisor and examiner on this Master’s thesis work, Christophe Demazière, for giving me the opportunity to do a project in nuclear reactor modelling and for his invaluable help and support. I would also like to extend my gratitude to Antonios’ Mylonakis, for the access to his MatLab code, and for the brief discussions we had, helping me to understand some key deeper aspects of the numerical procedures. Finally, I would like to thank my friends at and outside Chalmers as well as my family, for continuously challenging both my understanding and communication of physics by asking questions, intriguing and ridiculous ones alike. Erik Andersson, Gothenburg, June 2021 vii Contents List of Figures xi List of Tables xiii 1 Introduction 1 2 Molten Salt Reactors 3 2.1 Fission Reactors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 2.1.1 Criticality . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 2.2 Molten Salt Reactors . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 3 Neutron transport and diffusion 7 3.1 Neutron Transport . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 3.2 From Transport to Diffusion . . . . . . . . . . . . . . . . . . . . . . . 8 3.3 The Diffusion equation . . . . . . . . . . . . . . . . . . . . . . . . . . 11 3.3.1 Boundary conditions . . . . . . . . . . . . . . . . . . . . . . . 11 3.3.2 External sources . . . . . . . . . . . . . . . . . . . . . . . . . 12 3.3.3 Stationary equation . . . . . . . . . . . . . . . . . . . . . . . . 12 3.3.4 Dynamic equations in the frequency domain . . . . . . . . . . 13 4 Solving the Diffusion Equations 15 4.1 Discretization of the equations . . . . . . . . . . . . . . . . . . . . . . 15 4.1.1 Node average . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 4.1.2 Mesh definition . . . . . . . . . . . . . . . . . . . . . . . . . . 16 4.1.3 Matrix form of the discretized problem . . . . . . . . . . . . . 16 4.1.3.1 Source problem . . . . . . . . . . . . . . . . . . . . . 18 4.1.3.2 Dynamic matrix equation . . . . . . . . . . . . . . . 19 4.2 Numerical methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 4.3 Steps of verification . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 5 Performance of the Solver 23 5.1 Verification test results . . . . . . . . . . . . . . . . . . . . . . . . . . 23 5.1.1 Benchmark case 1 . . . . . . . . . . . . . . . . . . . . . . . . . 23 5.1.2 Benchmark case 2 . . . . . . . . . . . . . . . . . . . . . . . . . 24 5.1.3 Benchmark case 3 . . . . . . . . . . . . . . . . . . . . . . . . . 26 5.2 Mesh dependence in the precursor noise . . . . . . . . . . . . . . . . . 27 ix Contents 6 Conclusions and outlook 33 Bibliography 35 A Matrix coefficients I A.1 Static problem matrices . . . . . . . . . . . . . . . . . . . . . . . . . I A.2 Dynamic problem matrices . . . . . . . . . . . . . . . . . . . . . . . . II B Input and Output V B.1 Input . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . V B.2 Output . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . VI x List of Figures 4.1 Schematic of the mesh used to discretize the diffusion problem. . . . . 16 5.1 Static solution for the neutron fluxes (left) and the precursor concen- tration (right) for a heterogeneous system with u0 = 0 cm/s. Com- parison between the developed MSR solver and the 1D CORE SIM stationary solver. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 5.2 Noise solution for a heterogeneous system with u0 = 0 cm/s. The noise source is a perturbation in the removal cross-section, positioned at z ≈ 189 cm, with an amplitude of 1% of the static cross-section at a frequency of 1Hz. Comparison between the developed MSR solver and the 1D CORE SIM solver. . . . . . . . . . . . . . . . . . . . . . . 25 5.3 One-group static solution for the neutron flux (left) and precursor concentration (right) of a thermal reactor with fuel velocity u0 = 50 cm/s. Comparison to the solution developed by A. Mylonakis. . . . 26 5.4 One-group noise solution, at 1Hz, of a thermal reactor with fuel veloc- ity u0 = 50 cm/s. The noise source is a perturbation in the absorption cross-section, positioned at z = 150 cm, with an amplitude of 1% of the static cross-section. Comparison to the solver developed by A. Mylonakis. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 5.5 One-group static solution of a critical thermal reactor with very high fuel velocity. The neutron flux (left) compares the MSR solver at u0 = 100 000 cm/s to the (normalized) analytical solution for u0 =∞. The convergence of the precursor concentration (right) is shown for a range of velocities. . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 5.6 One-group noise solution of a thermal reactor with very high fuel velocity. The noise source is a perturbation in the absorption cross- section, positioned at z = 150 cm, with an amplitude of 1% of the static cross-section at a frequency of 1Hz. The absolute value (top left) and argument (bottom left) of the neutron flux noise compare the MSR solver at u0 = 100000 cm/s to the (normalized) analytical solutions for u0 =∞. The convergence of the precursor concentration absolute value (top right) and argument (bottom right) are shown for a range of velocities. . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 xi List of Figures 5.7 Absolute value of the one-group neutron noise, considering a per- turbation in the absorption cross section at different frequencies and location. The source is located at z = 150 cm for the top left and right and the bottom left figures. In the bottom right, the source is located at z ≈ 183 cm. In all cases the MSR solver at u0 = 100000 cm/s is compared to the analytical solution for u0 =∞. . . . . . . . . . . . . 30 5.8 The precursor concentration noise for a test problem with a source of frequency 1Hz, calculated at different numbers of cells N . . . . . . . 31 xii List of Tables 5.1 Values of α using f = 1Hz, u0 = 50 cm/s and L = 300 cm for varying numbers of cells N . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 B.1 Specification of the input matrix file XS_data_LIQ.mat. . . . . . . . . V B.2 Specification of the input matrix file S_data_LIQ.mat. One source is sufficient to run the program. . . . . . . . . . . . . . . . . . . . . . . VI B.3 Specification of the input matrix file DYN_data_LIQ.mat. . . . . . . . VI B.4 Specification of the input matrix file dS_data_LIQ.mat. One source is sufficient to run the noise calculations. . . . . . . . . . . . . . . . . VII B.5 Specification of the output matrix file RESULTS_LIQ.mat. Note that the units of all these quantities are only applicable in the case of an external static source, in which case keff will not be calculated. Without static source all units are arbitrary due to normalisation. . . VII xiii List of Tables xiv 1 Introduction The climate crisis is one of the biggest issues of our times and hence the importance of cleaner energy production cannot be understated. Electrifying the transport sec- tor, as well as reducing the use of fossil fuels in other areas of society, will require a substantial increase in low-emission electricity production. The total energy con- sumption is estimated to increase by about a third until 2060, and today more than half of the electricity production in the world (and about 80% of the total energy consumption) comes from coal, oil and gas [1]. Replacing this high-emission produc- tion while also expanding the production of energy is a very challenging task. Fission technology could have a huge impact on speeding up the process of replacing the fossil fuels in energy production. While currently utilized technology carries with it some very serious and well known problems, such as the risk of meltdowns and the storage of the nuclear waste for hundreds of thousands of years, new techniques can dramatically improve the situation. New fission reactor concepts, generally referred to as generation IV reactors, are being developed around the world, that could be cheaper to build and run, be safer and produce less waste than current running reactors. In the beginning of this century, the Generation IV International Forum (GIF), today consisting of experts from a wide range of countries including the US, the UK, Japan, China and Euratom (the collaboration for atomic energy research within the European Union), chose six reactor concepts that were deemed promising and should be investigated further. These concepts are all referred to as Generation IV fission reactors but vary greatly in function and design. From the Supercritical-Water- Cooled Reactors (SCWR) that resemble the water-cooled reactors in use today, to more novel concepts such as the Lead- and Sodium-cooled Fast Reactors (LFR and SFR), designs with helium cooled reactors (GFR and VHTR) and Molten Salt Re- actors (MSR) with liquid fuel circulating in the primary system. They were however chosen with a common set of goals in mind, concerning sustainability, economics, safety and reliability as well as a goal to be an unattractive route to extract weapons- usable material [2]. Common to all the chosen concepts is the higher efficiency in fuel usage and some of them could even be used as burner reactors, which could be loaded with nuclear waste from presently running reactors, decreasing the amount and storage time of high-level waste. The Molten Salt Reactor concept with fluid fuel circulating through the core offers many advantages. The salt acts as both a solvent to the fuel as well as a coolant, offering passive cooling and containment of the radioactive fuel and fission products. It is also one of the concepts offering the possibility of burning nuclear waste, reducing the amount and storage time [3]. Molten Salt Reactors have one 1 1. Introduction distinctive feature, setting them apart from the other concepts. In these reactors the nuclear fuel is dissolved in a molten salt which is circulating in the primary system of the reactor. Two important upsides of this are that the fuel is always cooled passively and that the reactor may be constructed to turn itself off in case of overheating [4]. The molten salt mixture will then simply cool and solidify, trapping the radioactive material and preventing far spreading of toxic material. The velocity of the fuel impacts the modelling of the neutron flux, which governs the behaviour of the reactor, due to the mechanism of delayed neutrons. In a fission reaction induced by the absorption of a neutron by the target nucleus, two to three new neutrons are released. However these are not always released instantly. A small fraction of them, called delayed neutrons, are released through the decay of certain fission fragments, called precursors, created at the moment of fission [5, p. 345]. The precursors then flow with the molten salt and decay at a later time, and thus, at a different position in the reactor. This fact is the key difference in modelling the neutronic properties of a Molten Salt Reactor as opposed to a reactor with solid stationary fuel. Taking this feature into account is the goal of this thesis. A tool, CORE SIM [6], for modelling the neutron flux in reactors with stationary fuel was previously developed by the task force on Deterministic REActor Modelling (DREAM group) at Chalmers University of Technology. This tool solves the system for the neutron flux, using two-group diffusion theory, meaning that the neutrons are divided into two energy groups, and the properties of each group are weighted averages. The system is solved both in a stationary case as well as for a perturbation in the frequency domain. An extended version of the tool, called CORE SIM+, is under development, with the objective to model reactors with greater spatial resolution [7]. The aim of this thesis work was to develop a similar tool to CORE SIM, that can handle moving precursors, for a 1-dimensional heterogeneous system. This solver will, given the input of the spatial distribution of material properties as well as properties of the system, solve the two-group diffusion equations with one averaged group of precursors. Since many of the properties, such as cross-sections, depend on the thermo- and fluid dynamics of the system, it is assumed that these properties are given. In order to fully model an MSR system one would have to couple this tool to other ones that handle the modelling of these properties given the neutron flux. The ambition has been that the neutron solver should be able to handle a very general input, with spatially heterogeneous material properties and fuel velocity, so as to allow for the most general systems to be investigated. It should also be able to handle external point-wise noise sources at any position in the reactor, as well as a possible external source of precursors at the inlet. The user is able to define properties of the system, such as reactor size and recirculation properties. In Chap. 2, the basics of fission reactors and the concept of MSRs is described in more depth. Then, in Chap. 3, the equations governing the neutron behaviour in the reactor core are presented as well as a derivation of a simpler set of equations implemented in the simulation tool. In Chap. 4, the implementation of the equations into a numerical solver is described, together with an analytical solution to a special case problem. Finally the solver is presented with some key features of the solutions as well as comparisons to the reference solutions in Chap. 5. 2 2 Molten Salt Reactors Fission reactors are complex multi-scale and multi-physics systems. From the small scale of the fission reactions to the large scale of the fluid dynamics. In nuclear fission reactors, heat is produced through the splitting of atoms. The total mass of the particles produced is smaller than the mass of the atom that was split and the mass difference is released primarily as kinetic energy of the fission products. These fission fragments then collide with other atoms in the fuel, heating it up. The heat generated in the fuel in turn heats up the coolant, which is pumped through the core, transporting the heat away. This heat is then used to boil water, creating steam at high pressure that drives a turbine, generating electricity. Simple as this may sound, there are certain complications. The material quantities that determine the neutron behaviour are all strongly dependent on the temperature of the reactor core. This means that the neutron transport, the thermodynamics and the fluid dynamics are all interdependent, and have to be solved simultaneously. In this chapter, some basic principles of fission reactors are presented, with emphasis on the neutronic behaviour in the reactor core, especially in the case of molten salt reactors. 2.1 Fission Reactors The basis of fission reactors is that a neutron induced reaction produces 2 to 3 new neutrons which in turn can induce fission in other nuclei, causing a chain reaction. In each such event, about 200 MeV of energy is released, some as radiation, but mostly as kinetic energy of the fission fragments, i.e. the nuclei resulting from the fission event. The fission fragments then collide with other atoms in the fuel, heating it up. The more neutrons moving around in the reactor core, the more fission events can occur, leading to more heat being produced. A very important part of modelling nuclear reactors is therefore the modelling of the neutron flux, which is strongly interdependent on the temperature distribution in the reactor [8]. The heat created in the fuel is transferred to some coolant, transporting the heat away to be used to generate electricity. In Light Water Reactors (LWR), which are the most common type of commercial reactors in use today, the fuel is in the form of solid pellets and the coolant is ordinary water which is either liquid under high pressure, or allowed to boil in the core, at lower (but still high) pressure. The neutrons produced in fission reactions are almost exclusively fast neutrons, with a mean energy of about 2 MeV, but the highest fission cross-sections are en- countered for the slower thermal neutrons with energies of a few tenths of an eV and below [5, pp. 344-345]. Therefore, neutrons created in fission reactions gener- 3 2. Molten Salt Reactors ally have to slow down before they induce new fission reactions. The creation of neutrons in the high end of the energy spectrum, the neutrons slowing down and the absorption of the thermal neutrons are together referred to as the neutron cycle, which is very important for the behaviour of the reactor. 2.1.1 Criticality In nuclear reactors it is of great interest to look at the ratio of neutrons in one generation of this cycle compared to the previous. This number, denoted keff , is the effective multiplication factor and describes the evolution of the reactor [5, p. 346]. If keff < 1 the reactor is said to be sub-critical. The neutron density in the reactor will then decrease exponentially and the reaction will die out. If keff > 1, the reactor is said to be super-critical. In this case the neutron density will increase exponentially and the reaction will quickly diverge unless the growth is quickly countered by introducing some neutron absorption, e.g. by inserting control rods of some highly absorptive material. During normal operation the reactor will be critical with keff = 1. In this case the neutron density is constant over time. The neutrons are not all produced at once though. Most of them, called prompt neutrons, are produced more or less instantaneously at the moment of fission, but there are also delayed neutrons [5, p. 345]. These are produced through the decay of certain fission fragments, known as precursors, created at the moment of fission. The half life of these may vary but is on the order of tenths of a second to tens of seconds. This is of great importance for the stability of the reactor. If there was only prompt neutrons keeping the reactor critical, any deviation from criticality would cause a very fast increase in neutron density [5, pp. 349-350]. Much faster than the control rods could be operated to safely counter the deviation, even if some other feedback mechanisms (such as the Doppler effect) would prevent any accidental situation. In practice however, the reactor is constructed to be sub-critical if disregarding the contribution of delayed neutrons. It is the delayed neutrons that push the system to be critical and due to the mean-lives on the order of seconds, deviations do not grow as fast and the system can be safely monitored and controlled. 2.2 Molten Salt Reactors Molten Salt Reactors (MSR) are a Gen. IV concept in which the fuel is in liquid form, dissolved in a molten salt. This sets them apart from other concepts, where the fuel is solid and stationary, introducing physical effects not present in other types of reactors. However, MSR systems offer several advantages over traditional reactors. The molten salt acts as a coolant and is circulated through the core region in a primary loop, carrying the fissile fuel particles with it. Inside the core, fission takes place and criticality is achieved, heating up the molten salt. This heat is then extracted outside the core through the use of heat-exchangers and used to generate electricity. It is also possible to process the fuel on-line, e.g. removing fission products or continuously refueling, while the reactor is running. The fuel is also always cooled passively by the molten salt and since it is already in a liquid state, which the system is built to handle, it cannot melt down [9]. Instead, it is 4 2. Molten Salt Reactors possible to construct the reactor to turn itself off in the case of overheating, without the need for external interaction or even power. This can be done by inserting a plug of frozen salt at the bottom of the reactor, which is cooled externally to keep it solid during normal operation [9]. In the case of overheating, the cooling will not be enough and the plug will melt, draining the molten salt and dumping it into cooling tanks with sub-critical conditions. In these tanks the salt will simply cool and solidify. This would trap the radioactive material, preventing far spread of toxic material. Another advantage of molten salt reactors is that they operate at, or close to, atmospheric pressure. This means that there is no risk of a large explosion happening as a consequence of any malfunction. The molten salt also allows for operation at higher temperatures, possibly offering higher efficiency in the electricity production. While the velocity of the fuel is very important for the thermo- and fluid dynamics of these types of reactors, this project is only concerned with the modelling of the neutron flux, which is also affected by this property. Since the precursors move with the fuel, delayed neutrons will be emitted downstream from where the fission event took place. Due to the comparable time scales of the decay and the fuel movement this becomes important when modelling the neutron distribution of these reactors. There is also the possibility of precursors exiting and reentering the system as well as the possibility of some on-line processing during operation. 5 2. Molten Salt Reactors 6 3 Neutron transport and diffusion The distribution of neutrons in a nuclear reactor is governed by the neutron transport equation. This is a balance equation depending on position, direction of movement, energy and time. These seven degrees of freedom, three for position, two for direction and one each for energy and time, make calculating general solutions a cumbersome task. For this reason, the simpler diffusion equation is often used. The basis of the diffusion equation is integrating the directional properties of individual neutrons, effectively replacing them with a net movement of neutrons through the concept of diffusion. The energy dependence can also be handled by simply integrating over the energy in certain intervals, called energy groups, yielding a set of equations for the neutron flux which now only depend on position and time, leaving four degrees of freedom. This chapter describes the derivation of the neutron diffusion equation from the transport equation. The important quantities concerning the modelling of neu- tron distribution in nuclear reactors are explained. From the diffusion equation, a stationary equation is shown, as well as the equation in phase space determining the response of the system to a monochromatic local perturbation. The necessary boundary conditions for solving these equations are also discussed. 3.1 Neutron Transport The neutron transport equations for a system with one group of delayed neutrons are derived as a balance equations per unit solid angle and energy in a short time interval dt and are given by [10, eq. 2:29] 1 v(E) ∂ ∂t ψ(r,Ω, E, t) + Ω · ∇ψ(r,Ω, E, t) + Σt(r, E, t)ψ(r,Ω, E, t) (3.1) = ∫ (4π) ∫ ∞ 0 Σs(r,Ω′ → Ω, E ′ → E, t)ψ(r,Ω′, E ′, t)dE ′dΩ′ + χprompt(r, E) 4π [1− β0(r)] ∫ ∞ 0 νΣf (r, E ′, t)φ(r, E ′, t)dE ′ + 1 4πχdelayed(r, E)λ0(r)C(r, t) and ∂C ∂t (r, t) = β0(r) ∫ ∞ 0 νΣf (r, E, t)φ(r, E, t)dE − λ0(r)C(r, t). (3.2) 7 3. Neutron transport and diffusion These equations determine the angular neutron flux ψ(r,Ω, E, t), which is the den- sity of neutrons at point r moving in the direction Ω, having energy E at time t multiplied with the neutron speed v(E) = √ 2E/m. They also contain the scalar neutron flux φ(r, E, t) = ∫ (4π) ψ(r,Ω, E, t)dΩ. (3.3) Together with the precursor concentration C(r, t), these quantities are the unknowns of these equations. Additionally there are several other quantities that need to be explained. The macroscopic total cross-section Σt accounts for all ways that neutrons can disappear from the region around (r,Ω, E) (e.g. absorption and scattering out of the region) except for traveling out of it spatially which is handled by the second term of eq. (3.1). The macroscopic scattering cross-section Σs accounts for the neutrons that are scattered into the region, i.e. it determines the number of neutrons scattered from direction Ω′ and energy E ′ to direction Ω and energy E. The third cross-section is the macroscopic fission cross-section Σf accounting for the neutron production due to fission. This quantity is combined with the average neutron production per fission event ν. The fission production term in eq. (3.1) is also multiplied with a factor (1 − β0(r)) and a factor β0(r) in eq. (3.2), where β0 is the fraction of delayed neutrons produced through fission. In addition, λ0 is the average decay constant for the precursor nuclei. Lastly χprompt and χdelayed are the energy spectra of produced neutrons coming from prompt and delayed neutrons respectively. For the purpose of MSRs where the fuel is moving, we need, however, to change eq. (3.2) to account for the fact that the precursors flow with the fuel. This means that an advection term ∇ · [u(r, t)C(r, t)] is needed and we can restate eq. (3.2) as ∂C ∂t (r, t) +∇ · [u(r, t)C(r, t)] = β0(r) ∫ ∞ 0 νΣf (r, E, t)φ(r, E, t)dE − λ0(r)C(r, t), (3.4) where u(r) is the velocity of the fuel. The above stated transport equation determines the distribution of neutrons in three space coordinates, two directional coordinates, energy and time. Altogether this makes seven degrees of freedom, making it a very difficult task to solve these equations. A simpler approximation, the diffusion equation, can be derived by making some assumptions about the system. This approach entails the loss of some information, such as the angular flux. It is however possible to extract the neutron current density vector J and from this approximate the angular flux. This is not as complete as the actual angular flux though. 3.2 From Transport to Diffusion The first step in reaching the diffusion equation is to get rid of the angular depen- dency by integrating over the solid angle Ω in eq. (3.1). Using eq. (3.3), in all terms whose angular dependency comes only from ψ, this quantity can just be exchanged for φ directly. Terms with no angular dependency, i.e. the two last terms on the r.h.s. get multiplied by 4π. This leaves the second term on the l.h.s. and the first 8 3. Neutron transport and diffusion term on the r.h.s. to be handled separately. The former can be handled by noting that Ω · ∇ψ = ∇ · [Ωψ] (since Ω is independent of position) giving∫ (4π) Ω · ∇ψ(r,Ω, E, t)dΩ = ∇ · ∫ (4π) Ωψ(r,Ω, E, t)dΩ ≡ ∇ · J(r, E, t), where the last uses the definition of the neutron current density vector J(r, E, t) = ∫ (4π) Ωψ(r,Ω, E, t)dΩ. For the scattering term it is assumed that the medium is isotropic (i.e. behaves the same in all directions). This means that the macroscopic cross-section only depends on the angle between Ω and Ω′. If, additionally, it is assumed that the scattering is isotropic in the laboratory system, the cross-section can be rewritten as Σs(r,Ω′ → Ω, E ′ → E, t) = 1 4πΣs,0(r, E ′ → E, t), allowing the simplification of the scattering term as∫ (4π) ∫ (4π) ∫ ∞ 0 Σs(r,Ω′ → Ω, E ′ → E, t)ψ(r,Ω′, E ′, t)dE ′dΩdΩ′ = ∫ ∞ 0 Σs,0(r, E ′ → E, t)φ(r, E ′, t)dE ′ where again eq. (3.3) was used to evaluate the Ω′ integral. Applying Fick’s law [11], J = −D∇φ, the full equation now reads 1 v(E) ∂ ∂t φ(r, E, t) +∇ · [−D(r, E)∇φ(r, E, t)] + Σt(r, E, t)φ(r, E, t) = ∫ ∞ 0 Σs,0(r, E ′ → E, t)φ(r, E ′, t)dE ′ + χprompt(r, E)[1− β0(r)] ∫ ∞ 0 νΣf (r, E ′, t)φ(r, E ′, t)dE ′ + χdelayed(r, E)λ0(r)C(r, t), where D(r, E) is the diffusion coefficient, assumed to be time independent. The validity of Fick’s law is based on the assumption of a medium that is not strongly absorbing, source free, infinite and homogeneous. In practice though, the approximation holds as long as these conditions are met on the scale of a few mean free paths. The equation has already become much simpler, but there is still the energy dependence to be handled. This can be done through a multi-group formalism. The idea is to integrate the equation over a discrete set of energy intervals [Eg, Eg−1], also called energy groups. Doing this allows for defining a set of group averaged quantities: • The scalar neutron flux integrated over group g φg(r, t) = ∫ Eg−1 Eg φ(r, E, t)dE. 9 3. Neutron transport and diffusion • The group averaged total macroscopic cross-section Σt,g(r, t) = 1 φg(r, t) ∫ Eg−1 Eg Σt(r, E, t)φ(r, E, t)dE, which can also be written as Σt,g(r, t) = Σa,g(r, t) + ∑ g′ Σs0,g→g′(r, t), where Σa,g is the macroscopic absorption cross-section, describing all neutron absorption in the material. • The group averaged scattering cross-section from group g′ to group g Σs0,g′→g(r, t) = 1 φg′(r, t) ∫ Eg−1 Eg ∫ Eg′−1 Eg′ Σs,0(r, E ′ → E, t)φ(r, E ′, t)dE ′dE. • The group averaged macroscopic fission cross-section for fission events caused by neutrons in group g′ νΣf,g′(r, t) = 1 φg′(r, t) ∫ Eg′−1 Eg′ νΣf (r, E ′, t)φ(r, E ′, t)dE ′. (3.5) • The group averaged diffusion coefficient Dg(r, t) = 1 ‖∇φg′(r, t)‖ ∫ Eg′−1 Eg′ D(r, E ′, t) ‖∇φ(r, E ′, t)‖ dE ′. • The integrated fission spectrum for delayed and prompt neutrons over group g χg = ∫ Eg−1 Eg χ(r, E)dE. Using two energy groups, fast and thermal (meaning E0 =∞, E2 = 0 and E1 is the boundary between the two groups), neutrons will be emitted as fast neutrons, i.e. χ1 = 1 and χ2 = 0 for both prompt and delayed neutrons. Furthermore it is convenient to define a removal cross section, representing the net scattering of neutrons from the fast to the thermal group, as Σr(r, t) = Σs0,1→2(r, t)− Σs0,2→1(r, t)φ2(r, t) φ1(r, t) . 10 3. Neutron transport and diffusion The diffusion equations can then be written as 1 v1,0(r) ∂ ∂t φ1(r, t) = ∇ · [D1(r, t)∇φ1(r, t)] + [(1− β0(r)) νΣf,1(r, t)− Σa,1(r, t)− Σr(r, t)]φ1(r, t) + (1− β0(r)) νΣf,2(r, t)φ2(r, t) + λ0(r)C(r, t) 1 v2,0(r) ∂ ∂t φ2(r, t) = ∇ · [D2(r, t)∇φ2(r, t)] + Σr(r, t)φ1(r, t) −Σa,2(r, t)φ2(r, t) ∂ ∂t C(r, t) +∇ · [u0(r)C(r, t)] = β0(r) [νΣf,1(r, t)φ1(r, t) +νΣf,2(r, t)φ2(r, t)]− λ0(r)C(r, t) (3.6) where eq. (3.5) was used to replace the energy integral in eq. (3.4). 3.3 The Diffusion equation In two-group diffusion theory with one group of delayed neutrons, the equations become a system of three equations. One for each of the fast and thermal neutron fluxes, respectively, and one for the concentration of precursor nuclei. The system of equations to be solved is then eq. (3.6). 3.3.1 Boundary conditions Finding distinct solutions to eq. (3.6) requires spatial boundary conditions for the fast and thermal neutron fluxes φg as well as for the precursor concentration C. For the neutron flux, the Marshak boundary conditions can be used, which are given by Jg(rb, t) · nb = 1 2φg(rb, t), ∀t (3.7) where rb is a coordinate on the boundary and nb is the outward unit normal to the boundary. This condition states that, in general, the amount of neutrons that leave the system (the L.H.S.) is equal to half the scalar neutron flux at the boundary (the R.H.S.). Another possible boundary condition is to set φg(rb, t) = 0 which will be used for most of the verification, due to its simplicity. For the precursor concentration the boundary condition C(rin, t) = Cin(rin, t), ∀t (3.8) is used. Here Cin(rin, t) is some function that may or may not depend on the precursor concentration C(rout, t′ < t) at the outlet of the system. This freedom of choice is necessary because certain models consider an on-line reprocessing of the fuel, filtering out certain components and adding others. Mathematically, the 11 3. Neutron transport and diffusion simplest case would be that of a closed primary system without fuel reprocessing. This would mean that the precursor concentration is reduced through decay while the fuel circulates outside the core. In that case, we could use the simple formula Cin(rin, t) = C(rout, t − τ) exp (−λ0τ), where τ > 0 is the recirculation time, i.e. the time the fuel spends outside the core each circulation. In the general case, the function takes the form C(rin, t) = Cin,ext(t) + ηC(rout, t− τ) exp (−λ0τ) (3.9) where Cin,ext is a source of precursors added at the inlet and η ∈ [0, 1] is a fraction of precursor nuclei filtered out in a potential on-line process. 3.3.2 External sources In certain cases, it is important to model a system with an external neutron source present. Adding a source term in the form of a monochromatic point neutron beam Sext(r,Ω, E, t) = Sextδ(r− rs)δ(Ω −Ωs)δ(E − Es) to eq. (3.1) and performing the same steps as for the rest of the equation, it is easy to see that the corresponding term in the diffusion equation becomes Sext(r) = Sg,extδ(r− rs), where g is the energy group containing the beam energy Es. In molten salt reactors it is also possible to introduce a source of precursors. Since this source would have to be introduced outside the core region, it can be modelled as an extra term in the boundary condition as shown in eq. (3.8). 3.3.3 Stationary equation To find the stationary solution to eq. (3.6), we simply set the time derivatives to zero and remove the time dependencies from all quantities. In the source-free case, a stationary solution is only possible for a critical reactor (with keff = 1) in which the number of neutrons is constant. If the reactor is not critical, the process is either dying or escalating exponentially. We can however rescale the neutron production terms, i.e. the ones containing νΣf,g, with a factor 1/keff , the reciprocal of the effective multiplication factor. This will, mathematically, make the system critical and thus a stationary solution can be found, giving valuable information about the properties of the equivalent critical system. The system of equations for the stationary problem without sources is thus 0 = ∇ · [D1,0(r)∇φ1,0(r)] + [ (1− β0(r)) νΣf,1(r) keff − Σa,1(r)− Σr(r) ] φ1,0(r) + (1− β0(r)) νΣf,2(r) keff φ2,0(r) + λ0(r)C0(r) 0 = ∇ · [D2,0(r)∇φ2,0(r)] + Σr(r)φ1,0(r)− Σa,2(r)φ2,0(r) ∇ · [u0(r)C0(r)] = β0(r) keff [νΣf,1(r)φ1,0(r) + νΣf,2(r)φ2,0(r)]− λ0(r)C0(r) (3.10) 12 3. Neutron transport and diffusion which is linear in the unknown quantities φ1, φ2 and C, allowing the problem to be stated in the form of a matrix eigenvalue problem as will be described in the following chapter. In the case of a source, the rescaling with keff is not needed or even possible. The appropriate source terms are simply added to the r.h.s. of the equations and the problem can be stated in the form of an ordinary matrix equation. However, in this case, a solution is only possible if the system is sub-critical in the absence of sources, since otherwise the solution would blow up exponentially. 3.3.4 Dynamic equations in the frequency domain For solving the dynamic case we assume small variations around a stationary mean, that is, for all time-dependent quantities we make the substitution X(r, t) = X0(r) + δX(r, t) and insert into eq. (3.6). Eliminating terms using eq. (3.10), disregarding higher order terms in the small perturbations and transforming to the frequency domain yields the set of inhomogeneous equations  iω v1,0(r)δφ1(r, ω) = ∇ · [D1,0(r)∇δφ1(r, ω)] + (1− β(r))νΣf,2,0(r) keff δφ2(r, ω) + [ (1− β(r))νΣf,1,0(r) keff − Σa,1,0(r)− Σr,0(r) ] δφ1(r, ω) + λ(r)δC(r, ω) + [ (1− β(r)) δνΣf,1(r,ω) keff − δΣa,1(r, ω)− δΣr(r, ω) ] φ1,0(r) +(1− β(r)) δνΣf,2(r,ω) keff φ2,0 + δS1(r) iω v2,0(r)δφ2(r, ω) = ∇ · [D2,0(r)∇δφ2(r, ω)] + δΣr(r, ω)φ1,0(r) −δΣa,2(r, ω)φ2,0(r) + Σr,0(r)δφ1(r, ω)− Σa,2,0(r)δφ2(r, ω) + δS2(r) iωδC(r, ω) +∇ · [u0(r)δC(r, ω)] + λ(r)δC(r, ω) = β keff [δνΣf,1(r, ω)φ1,0(r) + δνΣf,2(r, ω)φ2,0(r)] + β keff [νΣf,1(r)δφ1(r, ω) + νΣf,2(r)δφ2(r, ω)] (3.11) where keff is the same eigenvalue obtained from solving the stationary equations. If the stationary problem contained a source, this rescaling of course disappears here too. Note that the stationary source does not appear here since it cancels out when eliminating the stationary equations. With the exception of the source terms (those independent of δφ1, δφ2 and δC) and the terms containing iω, the structure of these equations is the same as for the stationary problem. 13 3. Neutron transport and diffusion 14 4 Solving the Diffusion Equations The system of equations derived in the previous chapter can be readily solved by dividing the reactor core into cells, defining node-averaged quantities for these cells and constructing matrices to solve one of two types of matrix equations. In the case of no external source driving the system the neutron contribution from fission is scaled by the effective multiplication factor keff to yield an eigenvalue equation. If, on the other hand, the system is driven by an external neutron source or added precursor nuclei the problem is an ordinary source equation. In both cases the solu- tion is a vector containing the discretized values for the neutron flux and precursor concentration. This chapter describes the process of deriving these matrices and how the matrix equations are solved numerically. Starting by defining the node averaged quantities and the node averaged equations, the coefficients of the matrices involved are de- rived. The Power Iteration method for solving eigenvalue matrix equations is then presented. Lastly, the verification process is described. 4.1 Discretization of the equations As both eqs. (3.10) and (3.11) are systems of linear equations they can be converted to matrix equations. The stationary problem becomes an eigenvalue equation in the case of no external sources, while in the presence of such sources it is an ordinary source equation. The dynamic system always has source terms by definition and thus always becomes an ordinary source equation. However, in order to convert the systems into matrix equations they have to be discretized on some mesh. One way to do this is to use node-averaged quantities. 4.1.1 Node average The node average of a certain quantity is defined as a volume average An = 1 Vn ∫ Vn A(r)dV (4.1) where Vn is the volume of the node. Applying node averages to eq. (3.10) or eq. (3.11) will for most terms just exchange the continuous quantities for their node averaged counterparts, resulting in one equation per node. The special cases are the terms containing ∇. From now on, the generality of three spatial dimensions is dropped and the systems considered will be one-dimensional and in the z-direction. In this 15 4. Solving the Diffusion Equations case, a node is defined by its boundaries, i.e. node n is the region z ∈ [zn−1, zn], where zn = zn−1 + ∆z. Starting with the diffusion term, applying the node average gives 1 ∆z ∫ zn zn−1 ∂ ∂z [ Dg,0(z) ∂ ∂z φg,0(z) ] dz = − 1 ∆z [Jg,0(zn)− Jg,0(zn−1)] , where the fundamental theorem of calculus was used together with Fick’s law J = −D∇φ = −D ∂ ∂z φ (in this 1D case, J is the net current in the positive z-direction). The advection term for the precursor concentration can be handled similarly as 1 ∆z ∫ zn zn−1 ∂ ∂z [u0(z)C(z)]dz = 1 ∆z [u0(zn)C(zn)− u0(zn−1)C(zn−1)] . 4.1.2 Mesh definition The mesh is defined by dividing the system into N cells of size ∆z. A schematic of the mesh and definition of the values used in computations is shown in Fig. 4.1. The node averaged quantities over these cells are then defined in the center. However, since the boundary condition in eq. (3.8) is defined by the outlet and inlet surface values and since the advection term contains the node surface values, the most natural points to use for the precursor concentration is at the edges of the cells. These are also the most natural points at which to define the velocity distribution u0,n, which also needs to be set at the inlet. Since the precursor concentration at the inlet, Cin, is given directly by the boundary condition, there is no need to explicitly include this value in the computations giving N points for C just as for φ. Using the surface values for C requires an approximation to be made for the node value, which also enters into the equations. The approximation made is simply taking the average of the surface values of each cell to get the cell volume average, i.e. Cn,node = Cn + Cn−1 2 . (4.2) The perturbations δφg,n and δCn are defined analogously to their stationary coun- terparts. Cin C1 Cn−1 Cn CN−1 CN Inlet Outlet φg,1 φg,n φg,N zn∆z∆z ∆z Figure 4.1: Schematic of the mesh used to discretize the diffusion problem. 4.1.3 Matrix form of the discretized problem Constructing matrix equations from the node averaged equations requires writing the neutron current Jg,0(zn) in terms of the node averaged quantities φn. This is done 16 4. Solving the Diffusion Equations by defining a temporary boundary value φb and then approximating J using Fick’s law and an ordinary finite difference scheme. We then demand that the estimation is the same from above as from below and get −Dg,n φb − φg,n ∆z/2 = Jg,0(zn) = −Dg,n+1 φn+1 − φb ∆z/2 . Eliminating φb gives Jg,0(zn) = − Dg,nDg,n+1 Dg,n +Dg,n+1 φg,n+1 − φg,n ∆z/2 . Letting n→ n− 1 gives Jg,0(zn−1) and the full term can then be written as − 1 ∆z [Jg,0(zn)− Jg,0(zn−1)] = 1 ∆z2 (anφg,n + bnφg,n+1 + cnφg,n−1), where an = − 2Dg,n−1Dg,n Dg,n−1 +Dg,n − 2Dg,nDg,n+1 Dg,n +Dg,n+1 , bn = 2Dg,nDg,n+1 Dg,n +Dg,n+1 and cn = 2Dg,n−1Dg,n Dg,n−1 +Dg,n . However, at the top and bottom of the system (node N and node 1 respectively) one of the nodes, n+ 1 or n− 1, respectively, is missing. In this case the boundary conditions for the neutron flux has to be employed. The Marshak condition in eq. (3.7) gives for the bottom node − φg(0) 2 = Jg,0(0) = −Dg,1 φg,1 − φg(0) ∆z/2 and again, eliminating φg(0) and combining with the expression for Jg(z1) yields a1 = − 2Dg,1Dg,2 Dg,1 +Dg,2 − 1/2 1 ∆z + 1 4Dg,1 , b1 = 2Dg,1Dg,2 Dg,1 +Dg,2 , and c1 = 0. A similar derivation for the top node gives aN = − 2Dg,N−1Dg,N Dg,N−1 +Dg,N − 1/2 1 ∆z + 1 4Dg,N , bN = 0, 17 4. Solving the Diffusion Equations and cN = 2Dg,N−1Dg,N Dg,N−1 +Dg,N . The discretized stationary source-free equations can now be stated as 0 = 1 ∆z2 [anφ1,n + bnφ1,n+1 + cnφ1,n−1] + [ (1− βn)νΣf,1,n keff − Σa,1,n − Σr,n ] φ1,n +(1− βn)νΣf,2,n keff φ2,n + λn 2 (Cn−1 + Cn) 0 = 1 ∆z2 [anφ2,n + bnφ2,n+1 + cnφ2,n−1] + Σr,nφ1,n − Σa,2,nφ2,n u0,n ∆z Cn − u0,n−1 ∆z Cn−1 = βn keff [νΣf,1,nφ1,n + νΣf,2,nφ2,n]− λn 2 (Cn−1 + Cn) and by defining a solution vector φ = φ1 φ2 C  , where φ1, φ2 and C are column vectors of length N (the number of nodes), the problem can be reshaped into a matrix equation. Since keff is an unknown eigenvalue it is convenient to define two different ma- trices. The matrix F contains all terms rescaled by keff (but does not contain the factor 1/keff) and the matrix M represents all other terms. The matrix equation is then written as Mφ = 1 keff Fφ (4.3) where the matricesM and F have coefficients as presented in App. A.1. Multiplying with keffM −1 from the left gives the ordinary eigenvalue problem M−1Fφ = keffφ. This equation has many solutions, in fact 3N solutions, but the most interesting is the solution corresponding to the highest eigenvalue, which is the true multiplication factor of the system. 4.1.3.1 Source problem In the case of an external source, the problem is no longer an eigenvalue problem. Instead a source vector S is added to the equations, giving Mφ+ S = Fφ which has the solution φ = (M − F )−1(−S). The source vector is S = Sφ1 Sφ2 SC  , 18 4. Solving the Diffusion Equations where Sφ1 = Sφ1,neutron + [λ1Cin,ext/2, 0, 0, ..., 0]T, Sφ2 = Sφ2,neutron and SC = [ Cin,ext [ −u0 ∆z + λ1 2 ] , 0, 0, ..., 0 ]T , where u0 in this case is the fuel velocity at the inlet. It is thus possible to have neutron sources in any of the two energy groups as well as a source of precursors at the inlet. 4.1.3.2 Dynamic matrix equation In the dynamic case there is always a source term in the equations, caused by one or more of the possible perturbations. Due to the similar structure of the equations it is possible to use the same matrices used in the stationary case, only having to add the terms containing −iω on the diagonal and multiplying terms affected by the precursor boundary condition with exp[−iωτ ] due to the time lag between outlet and inlet. The details can be found in App. A.2 but the equation becomes (M ′ − F )δφ = δS = δSφ1 δSφ2 δSC  , where δS is given by δSφ1,n = − [(1− β0)δνΣf,1,n − δΣa,1,n − δΣr,n]φ0,1,n − (1− β0)δνΣf,2,nφ0,2,n + δS1, δSφ2,n = −δΣr,nφ0,1,n + δΣa,2,nφ0,2,n + δS2 and δSC,n = β0[δνΣf,1,nφ0,1,n + δνΣf,2,nφ0,2,n]. We also have to add the term λ0δCin,const/2 to δSφ1,1 and λ0δCin,const to δSC,1 if we want to account for noise in the external precursor concentration source. To run the solver, only one of these terms is required to be non-zero, but combining several terms is possible. The equation has the solution δφ = (M ′ − F )−1δS. 4.2 Numerical methods In order to solve the matrix equations, numerical methods need to be employed to invert matrices and solve eigenvalue equations. To preserve the sparsity of the matrices, LUPQ decomposition can be used to handle the inversion of matrices. If we want to calculate A−1 for some N ×N matrix A we begin by employing LUPQ decomposition to write PAQ = LU ⇒ A = P−1LUQ−1, where L and U are lower and upper triangular matrices, respectively. The permu- tation matrices P and Q correspond to permutations of the rows and columns, and 19 4. Solving the Diffusion Equations are employed to make use of the sparsity of the equations. This decomposition is handled by the matlab function lu(A). We can then write A−1 = QU−1L−1P which can be easily calculated. To solve eigenvalue equations we can use the power iteration method. This is a very robust method for finding the largest eigenvalue of a matrix A and its corre- sponding eigenvector if the largest eigenvalue is relatively well separated from the second largest. The method consists of repeatedly applying the matrix A to an initial vector x and normalizing with the euclidean norm in every step, i.e. xn = Axn−1 ‖Axn−1‖ . As long as the initial vector x0, which we will choose to have all elements equal to 1/ √ N , has a non-zero component in the direction of the eigenvector we seek (i.e. the one with the largest eigenvalue) this vector will be amplified most in every step and for a large number of steps xn will approach the solution. If, however the highest eigenvalues are closely spaced, the convergence of this method can be very slow. This can be remedied by employing the Wielandt’s shift technique [6], through which convergence can be sped up by guessing the eigenvalue we seek. By estimating an eigenvalue kest and subtracting 1 kest F from both sides of eq. (4.3) we can rewrite the problem as( M − 1 kest F ) φ = ( 1 keff − 1 kest ) Fφ, which can be stated in the form of an ordinary eigenvalue problem as( M − 1 kest F )−1 Fφ = 1( 1 keff − 1 kest )φ = keffkest kest − keff φ. If the estimation is close enough to the value we seek this will generate a very large eigenvalue for this equation, making the power iteration method much more stable. The only requirement is that the estimation is closer to the eigenvalue we want than any other eigenvalue. This also makes it possible to use the power iteration method to extract all the eigenvalue-eigenvector pairs of the problem as long as no eigenvalues are degenerate. In order to get a good initial guess for the eigenvalues, the Arnoldi method is used [6]. It is based on building a Krylov subspace from the repeated applications of the matrix whose eigenvalues are to be computed, i.e. the space is spanned by the set of vectors {Anx0}Nn=0, where N is generally much smaller than the dimension of A. We can then project the matrix A onto the Krylov subspace, and compute the eigenvalues of the projection. From these we can then estimate the N largest eigenvalues of the original matrix. This method is very powerful, and by an iterative scheme it would be possible to close in on the true eigenvalues. However, the method does not guarantee conver- gence, which is one reason it is only used to get estimates in this case. Another is that, in the one dimensional case studied, the matrices do not generally become suf- ficiently large to warrant using this method, as may be the case in three dimensions, and thus the large speedup offered by the Arnoldi method is not as significant. 20 4. Solving the Diffusion Equations 4.3 Steps of verification In order to check the validity of the solver, three steps of verification have been carried out. Firstly, the full two-group solver was tested using u0 = 0 against the already existing CORE SIM [6] developed for reactor cores with stationary fuel. This requires changing the boundary conditions for the precursor concentration to be proportional to the fission neutron production at the inlet (in accordance with the Marshak boundary condition) instead of the recirculation condition of eq. (3.9). Secondly, the solver was simplified to only handle one group of neutrons. This simpler version was then tested against a one-group solver for a homogeneous MSR, developed and verified [12] by A. Mylonakis of the Chalmers DREAM group. The method of converting the full two-group solver into a one-group solver amounts to simply removing the matrix blocks concerned with the second group, ending up with 2N × 2N matrices. Thirdly the one-group version of the solver was tested against an analytical so- lution to the special case of u0 = ∞. The reference solution presented by Pázsit and Jonsson [13] amounts to analytical expressions for the neutron flux and precur- sor concentration of a stationary critical system (keff = 1) as well as an analytical expression for the dynamic problem in the frequency domain. The noise solution is given in the form of a Green’s function solving the equations in the case where the perturbation of the absorption cross-section is a delta function in position at a given frequency ω. Since the noise source term in general is of the form δS = δΣδ(z − z0)φ0(z), multiplying the Green’s function by δΣφ0(z0), where φ0 is the stationary solution, gives the correct solution for the noise. In both the stationary and dynamic cases, due to the infinite velocity, the solutions for the precursor concentration will be a constant given by integrating the solutions for the neutron flux. 21 4. Solving the Diffusion Equations 22 5 Performance of the Solver The MSR solver developed has very general capabilities in terms of the inputs. It is able to handle heterogeneous distributions in all material properties, including cross-sections, diffusion coefficients, decay constant λ, delayed fraction β and fuel velocity u0. The solver is also able to handle several types of hypothetical sources. Neutron sources in either energy group can be introduced as point sources in any node of the mesh and a source of precursor nuclei can be introduced at the inlet of the system. For the noise problem, the calculation is carried out at a set frequency for a perturbation to any cross-section and/or source. The full details of the format and options regarding the input are enclosed in App. B.1. It is worth noting that the solver is not capable of handling varying mesh size ∆z and thus it is not possible to have a more refined mesh in one part of the reactor. The solver outputs the stationary distribution for neutron fluxes and precursor concentration, together with the effective multiplication factor keff if no source was given. If a perturbation is supplied, the output also contains the noise solutions for these quantities. The format of the output is enclosed in App. B.2. 5.1 Verification test results To check that the solver performs as expected, three steps of verification have been carried out. 5.1.1 Benchmark case 1 First, the solver was tested on a heterogeneous two-group system with stationary fuel (u0 = 0 cm/s) and of height about 366 cm. The solution was then compared to the result of the previously developed CORE SIM and the static solution is shown in Fig. 5.1. The neutron flux solution can be seen to agree well with the stationary CORE SIM solver. Since the precursor concentration is not explicitly calculated in CORE SIM, the reference is evaluated by inserting the neutron fluxes into the equation for the precursor concentration in eq. (3.10) and eq. (3.11). In this case, due to instability in the MSR solver for u0 = 0 cm/s, the node averaged value Cnode,0 (given by eq. (4.2)) is used. This instability is probably due to the boundary condition given by eq. (3.9), which states that for u0 = 0 cm/s (implying τ = ∞), the precursor concentration at z = 0 cm is zero at all times. This is not true in this case. Instead, the concentration should be proportional to the neutron fluxes at the boundary, as seen from setting u0 = 0 cm/s in eq. (3.10). In conjunction with the 23 5. Performance of the Solver Figure 5.1: Static solution for the neutron fluxes (left) and the precursor concen- tration (right) for a heterogeneous system with u0 = 0 cm/s. Comparison between the developed MSR solver and the 1D CORE SIM stationary solver. fact that only Cnode,n enters into the equations, the face values Cn are not determined correctly. To compensate for the value used at the inlet being too low (zero instead of proportional to the fluxes), the value C1 is determined to be slightly higher to compensate, causing C2 to be too low and so on, creating a discontinuous oscillating solution for the face values Cn. Implementing the correct boundary conditions to fix this issue is not trivial and since the purpose of the solver is not to handle stationary fuel, this has not been done. However the node averages are still determined well, as can be seen in Fig. 5.1. The noise solution to the same problem is shown in Fig. 5.2. Here, the noise source was a perturbation to the removal cross-section at a frequency of 1Hz, located at z ≈ 189 cm. The amplitude of the perturbation was set to 1% of the static removal cross-section at that point. The results show good agreement between the solvers, but again, the δCnode values have been used. This is due to the same reasons, related to the boundary conditions, as in the static case. 5.1.2 Benchmark case 2 For the second step, the one-group version of the solver was tested on a homogeneous system with a fuel velocity u0 = 50 cm/s. The parameters chosen were those of a typical, slightly sub-critical (keff ≈ 0.9998), thermal reactor. In this case, a zero boundary condition was used for the neutron flux. The solution was compared to a one-group homogeneous solution by A. Mylonakis [12] and the result for the static solution is shown in Fig. 5.3. Here, the provided reference solution uses a 300 cell mesh, while the MSR solver has been set to 901 cells. The two solutions show good agreement in both neutron flux and precursor concentration. The noise solution for a perturbation to the absorption cross-section with fre- quency 1Hz, located at z = 150 cm, is shown in Fig. 5.4. Here the amplitude of the perturbation is 1% of the static absorption cross-section. For the neutron flux noise, the reference solver was run with a 300 cell mesh, while the precursor concentration 24 5. Performance of the Solver Figure 5.2: Noise solution for a heterogeneous system with u0 = 0 cm/s. The noise source is a perturbation in the removal cross-section, positioned at z ≈ 189 cm, with an amplitude of 1% of the static cross-section at a frequency of 1Hz. Comparison between the developed MSR solver and the 1D CORE SIM solver. noise is shown for both 300 and 4800 cells. The reference solution for 4800 cells was obtained by supplying the reference noise solver with the static MSR solver solution at this resolution. This more refined solution still agreed well with the lower resolution solution provided with the reference solver. In all cases, the MSR solver was run at 901 cells. The agreement in the neutron noise is very good, but in the precursor concentration, a clear mesh dependence is present in the reference solver. This dependence is also present in the MSR solver, however the effect is not as prominent at a few hundred cells and above. This will be investigated further in Sec. 5.2, but for now it suffices to say that both solvers seem to agree at high mesh resolution. Tests were also carried out at different velocities and for different frequencies. All of these gave similar agreement, with the same mesh dependence being visible for quickly varying solutions for δC. 25 5. Performance of the Solver Figure 5.3: One-group static solution for the neutron flux (left) and precursor con- centration (right) of a thermal reactor with fuel velocity u0 = 50 cm/s. Comparison to the solution developed by A. Mylonakis. 5.1.3 Benchmark case 3 Thirdly, the solver was tested against the analytical solution [13] for u0 = ∞. The system chosen was again the same typical thermal reactor, but the analytical solution involves a criticality equation that has to be fulfilled. This was achieved by adjusting the fission cross-section to fulfill the criticality equation. The static solutions for neutron flux and precursor concentration are shown in Fig. 5.5. Due to the nature of the numerical methods, it is not possible to set the fuel velocity to infinity. This has very little effect on the neutron flux solutions for sufficiently large velocity, which in this case was set to 100 000 cm/s, but the precursor concentration, which analytically is constant at infinite velocity, will never become constant in the MSR solver. For this reason, the precursor concentration is shown for a range of increasing velocities. The noise solutions, again calculated for a perturbation of the absorption cross- section at 1 Hz and located at z = 150 cm, are shown in Fig. 5.6. The amplitude of the perturbation was again set to 1% of the static cross-section value. As for the static solution, the neutron flux solutions was calculated at 100 000 cm/s, while the precursor concentration noise is shown for a range of velocities. Additionally, the absolute value of the neutron noise, calculated for a few other perturbations, are shown in Fig. 5.7. Here, the solutions for a perturbation at z = 150 cm are shown for different frequencies, 0.1Hz, 5Hz and 10Hz, respectively, as well as a solution at 1 Hz but where the position has been shifted to z ≈ 183 cm. In all cases the MSR solver was run at 100 000 cm/s and the amplitude of the perturbations was, as before, set to 1% of the static cross-section value. In all cases, the MSR solver reproduces the analytical solutions very well, showing that the solver performs well at high velocities. 26 5. Performance of the Solver Figure 5.4: One-group noise solution, at 1Hz, of a thermal reactor with fuel velocity u0 = 50 cm/s. The noise source is a perturbation in the absorption cross- section, positioned at z = 150 cm, with an amplitude of 1% of the static cross-section. Comparison to the solver developed by A. Mylonakis. 5.2 Mesh dependence in the precursor noise Returning to the mesh dependence in the precursor concentration noise, the result of the MSR solver for the problem corresponding to the one considered in Fig. 5.4 are shown in Fig. 5.8 for varying numbers of cells. A clear dependence on the mesh resolution is visible, and it is also clear that this is not just a problem with the so- lution being unresolved. Instead, at low numbers of cells, N , the peaks are spatially shifted towards the center. The most probable reason for this mesh dependence is the approximations made in the equation for the precursor concentration. In the MSR solver, the node values are estimated from the surface values using eq. (4.2), which is a linear approximation. The error in this equation can be estimated by applying the node average eq. (4.1) to the second order term in the expansion of δC(z). We get ∆δC ∝ ∆z2δC ′′(z). 27 5. Performance of the Solver Figure 5.5: One-group static solution of a critical thermal reactor with very high fuel velocity. The neutron flux (left) compares the MSR solver at u0 = 100 000 cm/s to the (normalized) analytical solution for u0 =∞. The convergence of the precursor concentration (right) is shown for a range of velocities. This explains why the static solution is not affected since it does not have the oscillatory behavior with larger second-order derivatives that the noise solution has. To get an estimate of what number of cells are required we can assume that the largest contribution to the second derivative comes from a sinusoidal function of wavelength u0/(2πf). That the wavelength should have this form can be seen from the fact that the oscillating source will cause the precursor concentration at the source to vary with time, and these variations are then carried away with velocity u0. Disregarding the amplitude of this sinusoidal (it will generally differ by at most O(1) relative to the overall magnitude of the solution) we can estimate the maximum relative error ∆δC ∼ ( L N )2 (2πf u0 )2 and formulate the requirement α = ( 2πfL u0N )2 � 1. Values of α for the solutions shown in Fig. 5.8 are given in Tab. 5.1. It seems that α < 0.01 would give good results, however, this should be seen as an indication of the size of N , rather than an exact condition, due to the approximations made. Checking the convergence will still be necessary. Table 5.1: Values of α using f = 1Hz, u0 = 50 cm/s and L = 300 cm for varying numbers of cells N . N 31 51 75 101 151 301 901 α 1.48 0.55 0.25 0.14 0.062 0.016 0.0018 28 5. Performance of the Solver Figure 5.6: One-group noise solution of a thermal reactor with very high fuel ve- locity. The noise source is a perturbation in the absorption cross-section, positioned at z = 150 cm, with an amplitude of 1% of the static cross-section at a frequency of 1Hz. The absolute value (top left) and argument (bottom left) of the neutron flux noise compare the MSR solver at u0 = 100000 cm/s to the (normalized) analytical solutions for u0 =∞. The convergence of the precursor concentration absolute value (top right) and argument (bottom right) are shown for a range of velocities. 29 5. Performance of the Solver Figure 5.7: Absolute value of the one-group neutron noise, considering a perturba- tion in the absorption cross section at different frequencies and location. The source is located at z = 150 cm for the top left and right and the bottom left figures. In the bottom right, the source is located at z ≈ 183 cm. In all cases the MSR solver at u0 = 100000 cm/s is compared to the analytical solution for u0 =∞. 30 5. Performance of the Solver Figure 5.8: The precursor concentration noise for a test problem with a source of frequency 1Hz, calculated at different numbers of cells N . 31 5. Performance of the Solver 32 6 Conclusions and outlook The next generation of fission reactors, currently under development, offers great improvements over, and solutions to several problems associated with, conventional reactor types. As a very promising Gen. IV concept, molten salt reactors may prove very important in the energy production of the future. The modelling of these reactors is thus of great importance. In this work, a simulation tool for calculating the neutron distribution in the core of MSRs, focused on simplicity and flexibility, has been developed. This solver was built on the foundation of the previously developed CORE SIM and the goal has been to handle very flexible inputs. The solver is also able to handle sub- and super-critical systems, with or without sources, and can also calculate the response of the system to perturbations in the frequency domain. The developed solver has been tested for a few special cases and performs as expected. For the precursor concentration noise solution, a dependence of the mesh size has been found and given a sufficiently high-resolved mesh, the solver gives good results. In this area, however, there are additional tests that can be performed to further validate the solver. Firstly, a comparison to analytical solutions for homo- geneous one-group systems can be made at finite fuel velocities [13], which could further verify the solver for a large range of fuel velocities. This would be similar to the test at infinite velocity done in this work, but the analytical solution would be more complicated. Secondly, it is possible to test the full heterogeneous two-group solver by a method that also involves solving the adjoint equations for the system [14]. This would require implementing the adjoint problem in the solver. Having carried out these additional tests, the next natural step would be to extrapolate the one-dimensional solver to the full three dimensions of a physical system. The theory presented in Chap. 4 naturally extends to three dimensions and while the implementation may prove more complicated, it should not be much more so than the implementations made in the existing three dimensional CORE SIM solver. The solver could also be expanded to handle several groups of precursors. This may be of great interest in MSRs, since the distribution of precursors may be very different for different groups, due to the transport mechanism. It may also be valuable to, in the context of online processing of the fuel, be able to set different removal fractions η for the different groups of precursors. 33 6. Conclusions and outlook 34 Bibliography [1] T. Kober, H.-W. Schiffer, M. Densing, and E. Panos, “Global energy perspectives to 2060 – wec’s world energy scenarios 2019,” Energy Strategy Reviews, vol. 31, p. 100523, 2020. [2] “A technology roadmap for generation iv nuclear energy systems,” tech. rep., 2002. [3] R. Moir, T. Dolan, S. McDeavitt, D. Williams, C. Forsberg, E. Greenspan, and J. Ahn, “Deep-burn molten-salt reactors,” 2008. [Accessed: 29/01/2021]. [4] W. N. Association, “Molten salt reactors,” 2020. [5] B. R. Martin and G. Shaw, Nuclear and particle physics : an introduction, 3rd ed. Wiley, 2019. [6] C. Demazière, “Core sim: A multi-purpose neutronic tool for research and education,” Annals of Nuclear Energy, vol. 38, no. 12, pp. 2698 – 2718, 2011. [7] A. Mylonakis, P. Vinai, and C. Demazière, “Core sim+: A flexible diffusion-based solver for neutron noise simulations,” Annals of Nuclear Energy, vol. 155, p. 108149, 2021. [8] V. Dykin and C. Demaziere, “Description of the models and algorithms used in the coupled core sim neutronic and thermo-hydraulic tool.,” 2014. [9] B. M. Elsheikh, “Safety assessment of molten salt reactors in comparison with light water reactors,” Journal of Radiation Research and Applied Sciences, vol. 6, no. 2, pp. 63 – 70, 2013. [10] C. Demazière, Modelling of nuclear reactor multi-physics : from local balance equations to macroscopic models in neutronics and thermal-hydraulics. Academic Press, 2020. [11] J. Lamarsch and A. Baratta, Introduction to Nuclear Engineering. Prentice-Hall, 2001. [12] C. Demazière, A. Mylonakis, and P. Vinai, “Development and test of a novel neutronic ver- ification scheme for molten salt reactors,” Proceedings of the 2021 American Nuclear Society Virtual Annual Meeting, June 14-16, 2021. [13] I. Pázsit and A. Jonsson, “Reactor kinetics, dynamic response, and neutron noise in molten salt reactors,” Nuclear Science and Engineering, vol. 167, pp. 61–76, 01 2011. [14] C. Demazière, A. Mylonakis, and P. Vinai, “Development and test of a novel verification scheme applied to the neutronci modelling of molten salt reactors,” (as of yet unpublished), May 2021. 35 Bibliography 36 A Matrix coefficients A.1 Static problem matrices To solve Mφ(+S) = [1/keff ]Fφ, choosing either the () or the [] to be true depending on the input (S being an external source vector), we have M = φ1φ1 0 φ1C φ2φ1 φ2φ2 0 0 0 CC  , consisting of nine N ×N blocks, and where the block AB is the matrix determining how quantity A depends on quantity B. The blocks are then constructed according to eq. (3.10) as follows: φgφg =  ǎ1 b1 c2 a2 b2 c3 a3 b3 . . . . . . cN−1 aN−1 bN−1 cN âN  , where ǎ1 = − 2Dg,1Dg,2 (∆z)2(Dg,1 +Dg,2) − 1/2 ∆z + (∆z)2 4Dg,1 − Σa,g + (g − 2)Σr, âN = − 2Dg,N−1Dg,N (∆z)2(Dg,N−1 +Dg,N) − 1/2 ∆z + (∆z)2 4Dg,N − Σa,g + (g − 2)Σr, an = − 2Dg,n−1Dg,n (∆z)2(Dg,n−1 +Dg,n) − 2Dg,nDg,n+1 (∆z)2(Dg,n +Dg,n+1) − Σa,g + (g − 2)Σr, bn = 2Dg,nDg,n+1 (∆z)2(Dg,n +Dg,n+1) and cn = 2Dg,n−1Dg,n (∆z)2(Dg,n−1 +Dg,n) . I A. Matrix coefficients Further we have φ2φ1 = diag(Σr,n) and CC =  d1 η exp[−λ0τ ]e1 e2 d2 . . . . eN dN  , where dn = [ un ∆z + λ0 2 ] , en = [ −un−1 ∆z + λ0 2 ] and η exp[−λ0τ ] ∈ [0, 1] is the fraction of the concentration at the outlet left when reaching the inlet. Lastly φ1C =  λ0 2 η exp[−λ0τ ]λ0 2 λ0 2 λ0 2 . . . . λ0 2 λ0 2  . For the fission matrix F we have F = diag[−(1− β0)νΣf,1,n] diag[−(1− β0)νΣf,2,n] 0 0 0 0 diag[β0νΣf,1,n] diag[β0νΣf,2,n] 0  . A.2 Dynamic problem matrices Writing the dynamic problem as (M ′ − F )δφ = δS = δSφ1 δSφ2 δSC  , we can use the same matrix F as for the static problem and in a similar manner we can write M ′ = δφ1δφ1 0 δφ1δC δφ2δφ1 δφ2δφ2 0 0 0 δCδC  . Making use of the similarities between the stationary and dynamic equations we can then write M ′ = M(η → η exp[−iωτ ]) +  diag ( −iω v1 ) 0 0 0 diag ( −iω v2 ) 0 0 0 diag(iω)  , II A. Matrix coefficients where we have added a phase to the terms related to the precursor boundary con- dition because of the time-lag τ . The source vectors is then δSφ1,n = − [(1− β0)δνΣf,1,n − δΣa,1,n − δΣr,n]φ0,1,n − (1− β0)δνΣf,2,nφ0,2,n, δSφ2,n = −δΣr,nφ0,1,n + δΣa,2,nφ0,2,n and δSC,n = β0[δνΣf,1,nφ0,1,n + δνΣf,2,nφ0,2,n]. We also have to add λ0δCin,const/2 to δSφ1,1 and λ0δCin,const to δSC,1 if we want to account for noise in the constant precursor concentration input. III A. Matrix coefficients IV B Input and Output B.1 Input The input to the program is divided over several matrices. Firstly the node size ∆z is given separately, in cm, with the name DZ, in a matlab matrix file called GEOM_data_LIQ.mat. The rest of the most basic input, required to run the program, is given in the file XS\data\LIQ.mat and the contents is specified in Tab. B.1. Table B.1: Specification of the input matrix file XS_data_LIQ.mat. Variable Name Dimension Description β0 [1] BETA (1× 1) Fraction of delayed neutrons. λ0 [s−1] LAMBDA (1× 1) Average decay constant of the precursors. Σa,1 [cm−1] ABS1 (N × 1) Absorption cross-section for group 1. Σa,2 [cm−1] ABS2 (N × 1) Absorption cross-section for group 2. Σr [cm−1] REM (N × 1) Removal cross-section. νΣf,1 [cm−1] NUFIS1 (N × 1) Fission cross-section for group 1. νΣf,2 [cm−1] NUFIS2 (N × 1) Fission cross-section for group 2. D1 [cm] D1 (N × 1) Diffusion coefficient for group 1. D2 [cm] D2 (N × 1) Diffusion coefficient for group 2. η [1] ETA (1× 1) The fraction of precursors removed through on-line processing. τ [s] tau (1× 1) The recirculation time. u0 [cm/s] u0 (N + 1× 1) The fuel velocity, defined at the interfaces of the nodes, including the inlet. In order to run the program with external stationary sources, these need to be specified in the matrix file S_data_LIQ.mat. The possible sources are shown in Tab. B.2 and several can be used at once. However, it is sufficient to specify one of these sources to run the program. V B. Input and Output Table B.2: Specification of the input matrix file S_data_LIQ.mat. One source is sufficient to run the program. Variable Name Dimension Description Cin [cm−3] C_in (1× 1) External source of precursors at the inlet. S1 [cm−3s−1] S1 (N × 1) External source of neutrons in group 1. S2 [cm−3s−1] S2 (N × 1) External source of neutrons in group 2. For the noise calculations, the neutron speeds and the frequency at which the cal- culations are performed are given in the matrix file DYN_data_LIQ.mat as specified in Tab. B.3. Table B.3: Specification of the input matrix file DYN_data_LIQ.mat. Variable Name Dimension Description v1 [cm/s] v1 (1× 1) Average speed of the neutrons in group 1. v2 [cm/s] v2 (1× 1) Average speed of the neutrons in group 2. f [Hz] f (1× 1) Frequency at which the noise calculations are done. For these calculations, a noise source needs to be specified in the file dS_data_LIQ.mat. The options for different noise sources are shown in Tab. B.4. Note that while it is possible to define several sources, one is sufficient to run the program. B.2 Output The output of the program is saved to the matrix file RESULTS_LIQ.mat. The description of these quantities are shown in Tab. B.5. The output files also contain information about the numerical methods, e.g. convergence criteria and number of iterations. VI B. Input and Output Table B.4: Specification of the input matrix file dS_data_LIQ.mat. One source is sufficient to run the noise calculations. Variable Name Dimension Description δΣa,1 [cm−1] dABS1 (N × 1) Perturbation of the absorption cross-section for group 1. δΣa,2 [cm−1] dABS2 (N × 1) Perturbation of the absorption cross-section for group 2. δΣr [cm−1] dREM (N × 1) Perturbation of the removal cross-section. δΣf,1 [cm−1] dNUFIS1 (N × 1) Perturbation of the fission cross-section for group 1. δΣf,2 [cm−1] dNUFIS2 (N × 1) Perturbation of the fission cross-section for group 2. δS1 [cm−3s−1] dS1 (N × 1) Perturbation of the neutron source for group 1. δS2 [cm−3s−1] dS2 (N × 1) Perturbation of the neutron source for group 2. δCin [cm−3] dC_in (1× 1) Perturbation of the external precursor source. Table B.5: Specification of the output matrix file RESULTS_LIQ.mat. Note that the units of all these quantities are only applicable in the case of an external static source, in which case keff will not be calculated. Without static source all units are arbitrary due to normalisation. Variable Name Dimension Description φ1,0 [cm−2s−1] FLX1 (N × 1) Stationary neutron flux for group 1. φ2,0 [cm−2s−1] FLX2 (N × 1) Stationary neutron flux for group 2. Cnode [cm−3] C_node (N × 1) Stationary node average precursor concentration. C [cm−3] C (N × 1) Stationary precursor concentration. keff [1] k_eff (1× 1) The effective multiplication factor. δφ1 [cm−2s−1] dFLX1 (N × 1) Neutron noise solution for group 1. δφ2 [cm−2s−1] dFLX2 (N × 1) Neutron noise solution for group 2. δC [cm−3] dC (N × 1) Precursor concentration noise solution. δCnode [cm−3] dC_node N × 1 Node averaged precursor noise. X (3N × 1) The full solution vector (only included for debugging purposes) VII DEPARTMENT OF PHYSICS CHALMERS UNIVERSITY OF TECHNOLOGY Gothenburg, Sweden www.chalmers.se www.chalmers.se List of Figures List of Tables Introduction Molten Salt Reactors Fission Reactors Criticality Molten Salt Reactors Neutron transport and diffusion Neutron Transport From Transport to Diffusion The Diffusion equation Boundary conditions External sources Stationary equation Dynamic equations in the frequency domain Solving the Diffusion Equations Discretization of the equations Node average Mesh definition Matrix form of the discretized problem Source problem Dynamic matrix equation Numerical methods Steps of verification Performance of the Solver Verification test results Benchmark case 1 Benchmark case 2 Benchmark case 3 Mesh dependence in the precursor noise Conclusions and outlook Bibliography Matrix coefficients Static problem matrices Dynamic problem matrices Input and Output Input Output