WHOLE CORE COUPLING METHODOLOGIES WITHIN WIMS

The ANSWERS® WIMS reactor physics code is being developed for whole core multiphysics modelling. The established neutronics capability for lattice calculations has recently been extended to be suitable for whole core modelling of Small Modular Reactors (SMRs). A whole core transport, SP3 or diffusion flux solution is combined with fuel assembly resonance shielding and pin-by-pin differential depletion. An integrated thermal hydraulic solver permits differential temperature and density variations to feedback to the neutronics calculation. This paper presents new methodology developed in WIMS to couple the core neutronics to the integrated core thermal hydraulics solver. Two coupling routes are presented and compared using a challenging PWR SMR benchmark. The first route, called GEOM, dynamically calculates the resonance shielding and homogenisation with the whole core flux solution. The second coupling route, called CAMELOT, separates the resonance shielding and pincell homogenisation from the whole core solution via generating tabulated cross sections. Both routes can use the MERLIN homogenised pin-by-pin whole core flux solver and couple to the same integrated thermal hydraulic solver, called ARTHUR. Heterogeneous differences between the neutronics and thermal hydraulics are mapped via thermal identifiers for neutronics materials and thermal regions. The ability for the integrated thermal hydraulic solver to call an external code via a FortranC-Python (FCP) interface is also summarised. This flexible external coupling permits one way coupling to an external fuel performance code or two way coupling to an external thermal hydraulic code.


INTRODUCTION
Whole core multi-physics modelling using higher fidelity solutions is becoming increasingly more viable due to improvements in algorithms and hardware. This approach aims to gain improved accuracy for existing and innovative reactor types. The ANSWERS® reactor physics code WIMS [1] is actively being developed for whole core multi-physics modelling including neutronics (resonance shielding, depletion, homogenisation and flux solve), thermal hydraulics and fuel performance.
In this paper, two new coupling routes between the neutronics and thermal hydraulics currently being developed in WIMS are presented. The KAIST PWR SMR benchmark [2] is used to analyse and compare these new routes. Differences between the two routes are investigated and their computational resources compared.
The capability for a user to define their own Python functions that use and manipulate WIMS state data at run time is also shown. This is being developed as a route to link to external fuel performance or thermal hydraulic codes.
These developments and analysis are an important step towards understanding the current models' capability and the computational resources required for this type of high fidelity whole core multi-physics modelling using WIMS.

WHOLE CORE COUPLING METHODOLOGIES
The industry standard for deterministic whole core modelling is the two-step approach. Initially a reactor physics lattice code generates tabulated cross sections for each fuel assembly type to be loaded into the core. The lattice code applies resonance shielding methods at different thermal conditions for a set number of irradiation intervals. Fuel assembly data is homogenised and collapsed to a reduced group structure. Assumptions are required regarding the environment and history during the lattice calculation, which can then be included as dependencies in the tabulation (e.g. buckling and historic coolant density).
Following the lattice calculation, a whole core nodal code solves for the global flux using the tabulated cross sections and then reconstructs pin powers using a shape factor from the lattice code. Typically, the whole core nodal code permits microscopic depletion and couples to a thermal hydraulic model that solves for an average pin and channel per fuel assembly.
Like all models, the two-step approach is bound in accuracy by any assumptions and approximations in the formulation. This includes tabulating fuel assembly cross sections using approximate environmental boundary conditions (e.g. infinite lattice) and with approximate histories (e.g. coolant density). Pin power reconstruction is another inherent approximation associated with this approach.
Advances in numerical methods and the availability of High Performance Computing (HPC) are permitting the development of higher fidelity whole core coupled models. For example, the Virtual Environment for Reactor Analysis (VERA) model [3] combines whole core transport flux solves (using MPACT) with dynamic resonance shielding and is coupled to the established COBRA-TF subchannel code. Another example, the NURISM framework, uses the SALOME open source platform to integrate established physics codes (such as the ANSWERS Monte Carlo code MONK® [4]) together to permit high fidelity simulations. While the aforementioned frameworks focus on coupling existing codes, the MOOSE framework [5] focuses on building each physics component within a common mathematical and software framework. This more cohesive approach offers direct implicit coupling and more effective utilisation of massive HPC resources.
ANSWERS is actively developing WIMS for whole core coupled modelling, initially targeting multiprocessor workstations and modest sized HPC resources. WIMS has multiple flux solvers suitable for whole core modelling for a range of reactor types. This includes the CACTUSOT Method of Characteristics (MoC) solver that uses modular ray tracing by axial slice, permits parallelisation by track and is suitable for reference calculations [6]. The MERLIN SPN flux solver uses a Finite Volume/Element Method spatial discretization and is suitable for homogenised pin-by-pin calculations. The WIMS module GEOM defines the whole core problem, controls the cross section preparation, permits control rod movement and performs critical searches [7]. Recently the ARTHUR thermal hydraulic solver has been included in WIMS and applied to transients for a typical individual PWR pincell [8]. This used an initial coupling algorithm that exchanges data via the WIMS interface data structure and relied on the WIMS controller to define the algorithm. While flexible, this initial approach is considered impractical for a user to create a whole core model.
For WIMS to offer a useful whole core coupled model, two new coupling routes have been developed. The first new coupling route, called GEOM, performs dynamic re-shielding and homogenisation of cross sections using the latest whole core state. The whole core state includes the latest thermal hydraulic solution, control rod positions, boron concentration and depleted isotopic compositions. The GEOM route can negate approximations (e.g. environmental and historic) inherent within the conventional two-step tabulation approach. This route is currently suitable for depletion analysis and reference solutions. The second new coupling route, called CAMELOT, is based on the conventional approach of tabulated cross sections. This route is currently suitable for design studies and transient analysis. Both routes can use the MERLIN homogenised pin-by-pin whole core flux solver. Therefore, an inter-pin power reconstruction is not necessary as in the conventional two-step approach. The following sections describe the two coupling routes in more detail.

GEOM coupling route
The GEOM coupling route starts with the WIMS module GEOM [9] to define the whole core problem using an engineering description of the geometry, materials and their mappings. Calculation meshes are overlaid on the core geometry to define the resolution of the resonance shielding, burnup, core flux and differential thermal feedback. Each material is tagged with a thermal identifier, which can be either fuel, can, coolant, moderator, absorber, structure or null. The thermal identifier is used to map between different heterogeneous geometric representations in the core neutronics and core thermal hydraulic solvers.
The GEOM module initiates a cross section preparation calculation for the whole core. Typically 2D lattice calculations for each shielding element (being a fuel assembly) using either Equivalence or Subgroup theory are performed. A collision probability flux solution in the WIMS 172 library group structure is then used to condense cross sections for the core solve (for example to 22 groups). A transport flux solution using the CACTUS module is used to form homogenised pincell cross sections via the Super-homogenisation (SPH) method. SPH factors are calculated for each pin type or for each pin position within the fuel assembly. The SPH factors are generated within an infinite lattice environment using the same homogeneous pin-by-pin flux solution method as will be used in the whole core solve that follows. This can either be a diffusion or SP3 calculation using the WIMS MERLIN module.
Following the cross section preparation calculation initiated by the GEOM module, the MERLIN module solves for the whole core flux using the homogenised pin-by-pin data. The MERLIN whole core flux solve can be nested within a critical search (e.g. concentration or control rod) that is controlled by the GEOM module. A coupled iteration can be defined to include thermal feedback using the WIMS ARTHUR subchannel module. Within a coupled iteration, the GEOM module can unsmear the whole core flux from MERLIN. This synthesizes the whole core flux with the heterogeneous flux for each shielding element that has been retained from the GEOM cross section preparation calculation. The heterogeneous flux is then used to determine the power in each thermal element for each thermal identifier (i.e. fuel, can, coolant etc.). A relaxation value in the range 0.4-0.6 is applied to this power to stabilise the coupled solver. The ARTHUR subchannel module uses the generated power to determine updated thermal feedback data within the active core regions. The heterogeneous power from the neutronics solution is mapped to the ARTHUR domain via assigning the thermal identifiers to each pin and channel region. This mapping is also used to average the feedback temperatures and coolant density from the ARTHUR domain back to the GEOM thermal mesh. The latest temperatures and coolant density defined on the thermal mesh are used to update the neutronics heterogeneous materials properties (microscopic cross section temperatures and isotopic number densities). The GEOM module then initiates the cross section preparation calculation using the updated neutronics materials properties. The coupled iteration continues until the change in core power between iterations is below a specified tolerance.
The GEOM route negates any conventional parameterised cross section error, including historic effects, and has the potential to capture more accurately environmental effects during the cross section calculation. The GEOM route is only applicable for rectangular lattice cores. The GEOM route also permits the CACTUSOT MoC [6] or MDLTRAN Sn [9] modules to be used for the whole core flux solve rather than MERLIN.

CAMELOT coupling route
The CAMELOT coupling route starts with the generation of tabulated material data for each fuel assembly type that will be loaded into the core. The preparation of the cross sections can use any established method within WIMS. Typically, this involves 2D lattice calculations using the same methods as discussed for the GEOM route (Equivalence theory for resonance shielding, collision probability solution for condensing spectrum, MoC solution for pincell homogenisation combined with SPH). Data can be generated by pin type or pin position for each fuel assembly type. Permitted tabulated data dependencies include temperature by thermal identifier (i.e. fuel, can, coolant etc.), coolant density and irradiation.
The CAMELOT module in WIMS uses the tabulated cross sections combined with the MERLIN and ARTHUR modules to solve for the coupled whole core state. MERLIN solves for the homogenised pinby-pin flux using either diffusion or SP3. Using the MERLIN flux solution and tabulated cross sections, CAMELOT forms a power for ARTHUR. ARTHUR can either solve for every individual pin and channel or more typically for an average pin and channel within each fuel assembly. CAMELOT maps data between MERLIN and ARTHUR using thermal identifiers similar to the GEOM coupling route. The feedback temperatures and coolant density within the active core regions determined by ARTHUR are used in subsequent coupled iterations to define the cross sections used for each region in MERLIN. Identical to the GEOM route, coupled iterations continue until the change in power falls below a specified tolerance. To stabilise the coupled iterations, the flux used to form the power can be relaxed. For a transient, either a tightly or loosely coupled approach is permitted. The tightly coupled approach iterates to convergence the coupled state for each time step. The loosely coupled approach performs a restricted number of coupled iterations, typically one or two.
To reduce the homogenisation error in the MERLIN pin-by-pin core solution, the SPH method generates correction factors during the tabulation generation step. The SPH factors can be stored in the tabulated data and interpolated between in two ways. Either they are folded directly into the cross section data or they are formed explicitly in CAMELOT from the corrected and uncorrected data. The latter approach permits the SPH factors to be applied correctly for the SP3 second moment equation (identical to the GEOM route). The latter approach also means the SPH factors and cross section data are interpolated between in CAMELOT separately rather than implicitly together, which has been found to produce more accurate coupled core solutions.
The CAMELOT route relies on the conventional parameterised cross section approach. By separating the cross section preparation and whole core solve, less computational resources are required to obtain coupled solutions for multiple states. This currently makes the CAMELOT route suitable for design studies and transient simulations. Separating the cross section preparation from the whole core solve permits a wider range of geometries in the latter. The CAMELOT route is applicable for rectangular lattice, hexagonal lattice and cylindrical coordinate represented (RZ or RTZ) cores. It also has potential for an unstructured mesh representation of the core.

PWR SMR BENCHMARK ANALYSIS
The GEOM and CAMELOT coupling routes are applied to cases based on the KAIST core benchmark [2]. The KAIST benchmark problem is a quarter core SMR containing 13 fuel assemblies. Surrounding the active core is an explicit steel baffle and a water reflector. Loaded into the active core are UOX assemblies at different enrichments and MOX assemblies. One of the UOX and one of the MOX assemblies contain Gd burnable poisons. Each fuel assembly contains 17x17 pins (of which 25 are guide tubes) on a standard 1.26cm pitch. The accuracy of the MERLIN homogeneous pin-by-pin model using different variations of the SPH factors is presented initially for a 2D model. Following this, 3D coupled solutions are represented for each fuel assembly type and then for the whole core exploiting quarter core symmetry. All results presented use Equivalence theory for the resonance shielding and 22 energy groups for the core solve.

2D core results
Both the GEOM and CAMELOT routes use the MERLIN homogeneous pin-by-pin solver for the whole core flux solution. The accuracy of the coupled models therefore depends on an accurate method to solve for the homogeneous pin fluxes, which are used to form the power. Both routes use SPH factors calculated by pin type or pin position. The GEOM route makes an approximation by combining pin position SPH factors with pin type smeared cross sections. WIMS inherently smears microscopic cross sections, which are mixed to form macroscopic cross sections. Microscopic cross sections are processed in the resonance shielding calculation by pin type, whilst differential irradiation depletes pin compositions by position. This methodology reduces the memory requirements associated with the whole core material data (primarily the microscopic cross sections). The CAMELOT route currently requires pin position SPH factors to be combined with pin position smeared cross sections. MERLIN diffusion or SP3 solutions can be used with both routes with pin type or pin position SPH factors.  Table I for a beginning of cycle core (BOC). The K-effective and pin power error measures are noticeably reduced when using SP3 rather than diffusion in MERLIN, and when using SPH factors by position rather than by type. The approximation the GEOM route uses when applying SPH factors by position incurs an insignificant error. All results are compared to a reference heterogeneous CACTUS MoC solution in the same 22 group structure.  For a BOC core WIMS is able to calculate pin powers using the MERLIN homogenised pin-by-pin solver, combined with SPH factors, to within an acceptable error compared to a reference heterogeneous transport solution.

3D coupled results
For the 3D coupled cases the original KAIST benchmark specification is extended to include a top and bottom reflector with a thickness of 21.42cm. The thermal-hydraulic data used in ARTHUR is characteristic of a typical PWR. Coupled GEOM and CAMELOT solutions for each fuel assembly type alone within an infinite lattice are compared in Table II. Close agreement between the GEOM and CAMELOT routes for these coupled single assembly cases is expected. The sole difference between the two routes for these cases is that the GEOM route determines cross sections by plane at exactly the calculated thermal conditions, whilst the CAMELOT route interpolates between tabulated cross sections defined at fixed interval thermal conditions (being 100 K for the fuel temperature and 0.1 g/cc for the coolant density). For fuel assemblies without poisons, close agreement exists irrespective of how the SPH factors are handled in CAMELOT. For the poisoned fuel assemblies, close agreement only exists if the SPH factors are handled explicitly in CAMELOT.  Both coupling routes have been applied to the BOC core solution, exploiting the geometric symmetry, to compare their computational requirements. Figure 2 shows the core geometry and pin powers and Figure  3 shows the thermal solution calculated as an average for each fuel assembly, all using the GEOM route. The coupled solution exhibits the expected physical behavior. Due to a lower coolant density at the top of the core, the power peak shifts below the midpoint causing the peak fuel center line temperature to shift also.
The total run-times for both coupling routes are found to be comparable (within 6%). For the CAMELOT route, this includes the cross section preparation which, for the current setup, takes about 90% of total run-time. Similarly, the cross section preparation constitutes 90% of the total run-time for the GEOM route. Both of these are dominated by the CACTUS MoC solve. This initial comparison indicates that the GEOM route with dynamic cross section preparation is a viable higher fidelity methodology, whilst the CAMELOT route permits multiple subsequent core solves at less cost if the same cross section tabulation is used.

MULTIPHYSICS EXTERNAL CODE COUPLING
A flexible multi-physics coupling between WIMS and an external fuel performance or thermal hydraulic code has been implemented. Python is embedded within WIMS using a general Fortran-C-Python (FCP) library developed by ANSWERS. WIMS state data at runtime is passed through an FCP interface to the Python interpreter. During the simulation, the Python interpreter executes Python code that can manipulate data exposed from WIMS. Python is ideally suited, being widely available, easy to learn and with specific features for scientific computing (e.g. NumPy). A bespoke Python adapter between WIMS and a specific external code handles the exchange of data. Because Python is interpreted rather than compiled, a user can couple an external code to WIMS without requiring direct access to the WIMS source code (assuming the necessary data from WIMS is accessible in the Python code). The Python adaptor can either manipulate the coupling data directly, pass the coupling data to a compiled library or exchange data with an external code via reading/writing to disk. The latter approach would often be required if there is no access to the external codes source or no desire to alter it for practicality or quality assurance reasons.
An FCP coupling interface exists in the ARTHUR module for one way coupling to a fuel performance code. An initial coupling between WIMS and the ANSWERS TRAFIC fast reactor fuel performance code has been created. The Python adapter generates a TRAFIC input file for each required pin using a template file combined with the power and grid passed from WIMS. A Python subprocess then executes TRAFIC. The TRAFIC output can then be analysed as if run alone. During a WIMS irradiation calculation, TRAFIC is executed multiple times using its ability to restart from a previous simulation. A second FCP coupling interface exists within the ARTHUR module for two way coupling to an external thermal hydraulic code. This permits an external code to be used instead of the ARTHUR subchannel solver. As a demonstration, this has been used to develop an initial coupling between WIMS and the external subchannel code COBRA-EN [10]. The flow of data between the two codes is shown in Figure 4. During the WIMS calculation the ARTHUR module is called and processes the neutronics data to obtain the core power and feedback fields (temperatures and coolant density) on the coupling grid. Rather than performing a subchannel calculation, ARTHUR uses the FCP interface to pass the core power and coupling fields to a Python function that exists within a Python module that is external to WIMS (i.e. wims2cobraen.py). The FCP library imports the external Python module into an internal host module. A Python dictionary is created in the host module such that its lifetime can be controlled by WIMS. This dictionary is populated with the coupling data from WIMS. NumPy arrays are used to efficiently pass array data from WIMS. The FCP library then executes the required Python function that exists in wims2cobraen.py, passing the dictionary populated with the coupling data as the only argument (which keeps the function signature simple). The Python adapter code reads a template input file, generates a valid input file via writing the coupling data from WIMS (i.e. power), executes COBRA-EN as a subprocess, parses the COBRA-EN output file to extract temperatures and coolant densities and then updates the coupling data held within the Python dictionary. The FCP library then transfers the coupling data back to the ARTHUR module which itself feeds this into the whole core coupling algorithm in WIMS.
This initial two way coupling to an external code has been verified by comparing WIMS using the ARTHUR subchannel solver and WIMS coupled to COBRA-EN for a single steady state PWR UOX fuel assembly using the GEOM route. GEOM generates homogenised pin-by-pin data for MERLIN which solves for the pincell homogenised flux. GEOM recovers the heterogeneous flux via unsmearing and forms a power on the thermal mesh. One thermal mesh element in each of the axial planes is used to define the thermal feedback. ARTHUR and COBRA-EN model one average pin and channel for the fuel assembly. ARTHUR and COBRA-EN apply the same physical models for a single pin/channel and have previously been shown to produce comparable solutions for this scenario [8]. The differences between the two coupled solutions is 1.8pcm for the K-effective and 0.12K for the maximum pin temperature. This confirms that the coupling between WIMS and COBRA-EN via FCP functions correctly.
The two-way FCP interface to an external code has the potential to couple to a variety of thermal hydraulic models including system, subchannel and CFD based. This approach is currently being applied to couple WIMS to the ANSYS CFD code CFX for high fidelity HTR modelling.
A third FCP interface exists as a base module within WIMS. This exposes all the WIMS interface data that is passed internally between modules. This includes the cross sections, isotopic number densities and the current flux solution. This low-level FCP interface has the potential to permit coupling to any applicable external code, not just thermal hydraulic models. For example this approach is currently being investigated to train and drive a surrogate Reduced Order Model.
More generally, the FCP library could be used to permit any WIMS state data to be exposed in Python. An FCP interface exists within ARTHUR to define the Nusselt number used in the heat transfer between pins and channels. This permits users to apply their own closure relationship without exposing the specific details. This approach will be applied to all closure relationships within ARTHUR in the future.

CONCLUSIONS
Two new coupling routes between core neutronics and core thermal hydraulics for the WIMS reactor physics code have been presented. The CAMELOT route uses tabulated cross sections which separates the resonance shielding and pincell homogenization from the coupled core solution. The GEOM route performs resonance shielding and homogenization dynamically with the coupled core solution, removing key approximations involved in the two-step approach. Both routes use homogenised pin-by-pin flux solutions, directly removing the requirement for reconstruction. The SPH method was shown to reduce the pin power error to acceptable values compared to a transport reference. The approximation in the GEOM route for pin position SPH factors was shown to not incur any significant error. The capability of both approaches was demonstrated via an extension of the KAIST core benchmark. Analysis of the runtime for the GEOM route showed that around 90% is in the dynamic generation of the cross sections, which is dominated by the run-time for the fuel assembly MoC solves. The dynamic shielding implementation within the GEOM route is currently being parallelised to reduce the effective run-time.
Developments in WIMS to communicate with external codes via the FCP interface was presented. The data flow between WIMS and an external code was demonstrated via an initial coupling to COBRA-EN. The FCP interface offers the potential for coupling with a variety of codes including fuel performance, subchannel, system and CFD codes. It also allows user controlled Python functions to manipulate key WIMS parameters directly, such as to define heat transfer correlations within the ARTHUR module.
The multi-physics capability presented is currently being applied to a range of reactor types and scenarios to extend the verification and validation. This includes depletion and transient analysis for conventional large core PWR benchmarks.