Calculating α Eigenvalues and Eigenfunctions of One-Dimensional Media

1. Calculating α Eigenvalues and Eigenfunctions of One-Dimensional Media Benjamin R. Betzler, William R. Martin, Brian C. Kiedrowski, Forrest B. Brown 19 September 2013

2. Acknowledgment • PhD research of Ben Betzler at the University of Michigan • Work started on this at Los Alamos in the summer of 2012 with collaboration with Brian Kiedrowski and Forrest Brown

3. Outline • The forward α eigenvalue problem with delayed precursors • The adjoint αeigenvalue problem with delayed precursors • The transition rate matrix (TRM) for the adjoint problem • Estimation of the adjoint TRM elements using forward keff Monte Carlo for a specified configuration with delayed neutrons • Determination of the forward TRM and its eigenfunctions and eigenvalues • Standard eigenfunction expansion to solve for arbitrary time-dependent sources for the specified configuration • Numerical results for 1-D slab problems with comparison to Green’s Function Method (GFM) and time-dependent Monte Carlo

4. Motivation to pursue the α eigenvalue problem • The time- or α-eigenvalue spectrum is necessary for calculating eigenfunction expansions for time-dependent problems • When the spectrum (and its eigenfunctions) are known, the flux response within the system may be calculated for any source • Some applications of this work are pulsed neutron experiments and accelerator-driven subcritical systems, where the flux varies spatially in time • Critical systems may have higher α eigenvalues and spatial eigenfunctions that may be excited by a change in neutron population • Objective is to show the Transition Rate Matrix Method (TRMM) is able to calculate these eigenvalues and eigenfunctions during a Monte Carlo random walk

5. Use separation of variables to obtain the α eigenvalue problem • To formulate the eigenvalue problem, separate the temporal dimension from position, energy, and direction: • The largest real eigenvalue (the fundamental) has the trend and its corresponding eigenvector must be all one sign • All other eigenvalues and eigenvectors may be real or complex (complex eigenvalue must have a conjugate)

6. The forward α-eigenvalue problem • Introducing this separation into the time-dependent neutron transport equation yields the eigenvalue problem: • We desire the eigenvalues of the combination of the left product matrix, but not the elements of this matrix can be tallied (e.g., vλ) during a Monte Carlo simulation • Formulating the time-dependent neutron importance equation (the adjoint) yields another matrix but with physical quantities that can be tallied during a forward Monte Carlo simulation

7. The adjoint α-eigenvalue problem • The adjoint α-eigenvalue problem is • The matrix elements are now reaction rates or decay rates • Discretize phase space into a set of bins and keep track of the neutron (or precursor) population in each bin • These rates describe transitions among the discrete state space bins, for example: (1) a precursor at a position decaying into a neutron with a given energy and direction or (2) a neutron with a given energy and angle scattering to another energy and angle • The transition rate matrix is the left matrix product

8. The transition rate matrix • This matrix is similar to a continuous-time Markov process with a transition rate matrix • The diagonal elements q are rates of transitions out of that state • The off-diagonal elements are rates of transition between states qij = transition rate from state j to state i (i.e., qij) • Note: the fission matrix formulation is a discrete-time Markov process where fission generations represent “time”

9. The adjoint transition matrix elements are tallied with a forward Monte Carlo simulation – with careful bookkeeping • For transition rates between energy states (top left): • For rates from energy states to precursor groups (top right): • For rates from precursor groups to energy states (bottom left): • For rates between precursor groups (bottom right):

10. Eigenfunction expansion of the solution and source • The time-dependent angular flux is assumed to be separable in time: • The time dependence of the initial source determines the function; for a pulsed neutron source, it is • The coefficients are calculated with the adjoint and forward shape eigenfunctions (using bi-orthogonality) and the given source

11. Determination of the forward transition rate matrix and its eigenvalues and eigenfunctions • The adjoint transition rate matrix is estimated from tallies during a forward k-eigenvalue Monte Carlo calculation that describes the system that will be hit with a given source • Transition rates are calculated from the tallies and the transition rate matrix is formed (careful bookkeeping is required to keep in mind these are actually adjoint transition rates) • The velocity matrix is formed and the forward matrix is obtained by swapping the speeds and transposing the resultant matrix • Eigenvalues and eigenfunctions of the forward transition matrix are determined using LAPACK (direct Schur factorization) Assertion: using arguments similar to what Forrest will use for the fission matrix, the limit of continuous bins is equivalent to using adjoint continuous energy Monte Carlo to obtain the TM

12. Test problems with comparisons to Green’s function method (GFM) and time-dependent Monte Carlo • One speed slabs – eigenvalues only • Alternating slabs with incident pulse • Loosely coupled critical system with incident pulse • 5-region problem with subcritical, critical, and supercritical configurations and several imposed time-dependent sources • Continuous energy alternating slab problem (no Monte Carlo solution for this case)

13. One group slab: the point spectrum agrees with GFM but the continuous spectrum is “pointized” • 5 mfp homogeneous slab • Has three real eigenvalues above the continuum spectrum (less than -1) • 50 equal position states and 18 equal angle states • Eigenvalues form along rings in the continuum TRMM GFM

14. Observations on the one-group slab spectrum • The number of rings of eigenvalues in the continuum is half that of the number of angular states • The diameter of the “rings” increases with the number of position states; attributed to position states being smaller resulting in faster leakage • Some eigenvalues that are part of the continuum leak into the real portion of the spectrum (greater than -1), but this decreases as the number of angles and positions is increased • These phenomena with the point spectrum in the continuum were previously seen using an Sn-method to obtain the eigenvalue spectrum

15. Results: Multi-Region Problem • Alternating materials of equal thickness • Material 1 is purely scattering and material 2 has a small absorption cross section • The speed is effectively 10 so the continuum portion of the spectrum starts at -10

16. All calculated eigenvalues match within 1% to GFM-calculated values

17. Eigenfunction expansion solution for incident pulse compared to time-dependent Monte Carlo

18. Five-region loosely coupled system • For this five-region problem, the right fuel region thickness is either 1 (symmetric) or 1.1 (asymmetric) • The fuel is varied to make subcritical, near critical, and supercritical configurations • All configurations have only two real eigenvalues

19. The calculated eigenvalues for all 5-region configurations compare well to those from the GFM

20. Eigenfunction comparison with GFM for 5-region problem – symmetric case • The first two eigenfunctions for the symmetric, near critical case agree with the GFM results • The second shape eigenfunction has rotational symmetry • These eigenfunctions have the most difficulty converging due to the proximity of the first two alpha eigenvalues (-0.00615 and -0.00644)

21. Eigenfunction comparison with GFM for 5-region problem – asymmetric critical case • The asymmetric fundamental shape eigenfunction increases nearly three orders of magnitude towards the thicker fuel region, where the second shape eigenfunction is nearly flat • These converge faster than the symmetric case because the first two eigenvalues differ greatly

22. Now impose different source conditions on several of the 5-region configurations and compare with TDMC • Symmetric with constant source from left • Asymmetric supercritical with incident pulse on left • Asymmetric near critical with incident pulse on left • Same case to determine detector response in the right fuel region

23. Subcritical symmetric configuration with a constant source incident on the left

24. Asymmetric supercritical configuration with a pulse incident on the left face

25. Eigenfunction expansion for the near critical symmetric configuration with a pulse incident on the left face

26. The estimated detector response for critical symmetric problem for a detector in the right fuel region • Shows the ability of the TRMM to accurately predict the response in a given region of the problem • The higher modes are still present throughout a large portion of the response, due to the proximity of the first two eigenvalues

27. Observations on the 5-region problems • The TRMM eigenfunction expansion solutions compare well to TDMC • The eigenfunction expansions shown use all modes in the expansion, even those belonging to the eigenvalues calculated that fall within the continuum; these modes are necessary to model the earliest times where the flux is zero in much of the problem • The smoothness of the expanded solution is not able to accurately describe sharp variations in the flux shape • Differences at the front of the pulse shapes at early times are due to the differences in the angular dependence of the initial source and the source described to the eigenfunction expansion

28. Continuous Energy Problem • This continuous energy problem is a subcritical hydrogenous medium set up with five fuel regions and an outer reflector • The total problem size is 280 centimeters • The detector response is a linear combination of the flux in each of the energy groups (mostly thermal)

29. The eigenfunction expansion solution for a pulse on the left side of the problem (no TDMC)

30. The estimated response for a detector located in the right fuel region for an incident pulse on the left • The estimated response for a detector located in the right fuel region for the pulse incident from the left side of the reactor • The prompt fundamental mode is plotted alongside the expected detector response • After 40ms, the delayed neutron modes begin to dominate the response shape

31. Summary and Future Work • Summary • Formulated the Transition Rate Matrix Method • Tallied the continuous-time Markov process transition rate matrix with Monte Carlo • Verified one-speed results with Green’s Function Method • Investigated the spectrum behavior • Compared TRMM eigenfunction expansions to TDMC • Future work • Sparse matrix storage and eigenvalue routines • Continuous energy TDMC comparisons • Calculate three-dimensional alpha modes • Benchmark to experimental data (e.g., Fort St. Vrain pulsed neutron experiments and MUSE experiments)

32. Questions?

33. The α-eigenvalue problem with delayed precursors

34. Methods: Eigenfunction expansion • Flux expansion in the time-dependent transport and precursor equations • Resulting differential equation