An OpenFOAM solver for multiphysics modeling of fusion reactor design: The nemoFoam code

This paper introduces a multiphysics code, named nemoFoam, designed within the OpenFOAM environment to address the coupling of neutron and photon transport with thermal-hydraulics in nuclear reactor simulations. The code, conceived both for fusion and fission applications, employs a modular approach where the neutronic module is currently based on multi-group neutron diffusion equations, including a mono-kinetic treatment for photons. The thermal-hydraulic module is built on the standard OpenFOAM solver. The paper focuses on presenting the key features of nemoFoam and discussing the structure of the code, the nuclear and thermal-hydraulic models, and the coupling strategy between them. In order to assess the performance of the neutronic module, the code is applied to a fusion case study, indeed a benchmark against the Serpent Monte Carlo code for neutron and photon transport is performed, applying nemoFoam to the Affordable, Robust and Compact (ARC) fusion reactor design. Simulation results demonstrate good agreement with Monte Carlo benchmarks, emphasizing the potential of the code for multiphysics simulations. The modular design of nemoFoam allows users to extend the implemented model equations to study additional phenomena, focus only on selected aspects independently and, potentially, to add new solvers in each module.


Introduction
The design and safety analysis of nuclear reactors requires modeling different physical phenomena, mainly thermal-hydraulics and neutronics, which are typically addressed by several different codes.Usually, the high-fidelity codes used to address the 3D thermal-hydraulics of the system are CFD software (e.g.StarCCM+ [1], Ansys [2], COM-SOL [3], OpenFOAM [4]), while neutronics and photonics are solved by means of Monte Carlo codes such as MCNP [5], Serpent [6], or OpenMC [7].Thus, the modeling approach required for nuclear reactor design is naturally multiphysics [8,9].The coupling of different codes in multiphysics simulations often encounters challenges related to data exchange interfaces and code modifications.
Several attempts have been made to solve this class of problems, either by coupling single-physics models and codes with proper implicit or explicit algorithms for data exchange, or addressing the whole set of equations describing the multiphysics problem within a single solver.For example, the FERMI environment [10], developed at Lawrence Livermore National Laboratory, couples three different softwares, MCNP (neutron and gamma transport), OpenFOAM (thermal-hydraulics), and Diablo [11] (structural mechanics) through the preCICE coupling library [12], to allow the high-fidelity modeling of a fusion reactor blanket.On the other hand, the code GenFoam [13] is a multiphysics solver for fission reactor analysis, which couples thermal-hydraulics, neutronics, and thermal-mechanics thanks to the use of the OpenFOAM library.An alternative approach is represented by MOOSE [14], developed at Idaho National Laboratory (INL), which is an open-source platform for the development of multiphysics codes based on FEM discretization.It allows focusing on the development of physics-based modules without dealing with the interfaces between the various modules.
All the codes mentioned earlier focus on the detailed resolution of the fundamental equations on a 3D grid, offering a precise representation of the physical phenomena at the cost of significant computational resources.In contrast, system-level codes serve a different purpose, focusing on the simulation of the entire system rather than on the detailed modeling of individual components or groups of them.These codes tend to use simplified models and approximations in order to study the operation of complete reactors.For example, the RELAP5-3D code [15], a tool developed at INL for the analysis of transients and accidents in nuclear power plants and the analysis of advanced reactor designs, or, the GETTHEM code [16], developed ad Politecnico di Torino, which is a system level code based on Modelica language for the simulation of the Primary Heat Transfer System and Balance of Plant of tokamak reactors.https://doi.org/10.1016/j.nme.2024 This work aims at developing a multiphysics code, nemoFoam, where both the neutronic and thermal-hydraulic models are solved in the same environment, thanks to the adoption of the C++-based library OpenFOAM.The nemoFoam code exploits the OpenFOAM capabilities in discretization, parallel solution of partial differential equations, and finite volume methods on unstructured meshes, allowing to develop a computational tool which solves the multiphysics problem with a consistent level of accuracy for all the physical models included in the solver.The approach utilized in the nemoFoam multiphysics solver consists in the development of a multi-group neutron and mono-energetic photon diffusion module, which is coupled with a thermal-hydraulics module based on the standard multi-region solver of OpenFOAM.Designed to model fusion reactors, nemoFoam allows to define an external source of neutrons and/or photons outside of the discretized domain.The nemoFoam code is a dual-scale detail tool: ''high-fidelity'' for thermal-hydraulics modeling (CFD) and ''low-order'' for the neutronic part (multi-group diffusion).Compared to a dual high-fidelity approach (CFD + Monte Carlo), like the FERMI project, the advantage of such an approach lies in the reduction of computational cost, providing a lighter and faster tool for conducting blanket design analyses.The simulations are, in principle, 3D, but the solver can perform 2D simulations as well, provided that suitable geometries and boundary conditions (BCs) are given.The nemoFoam code is fully modular so, if needed, single physics analysis can be performed.
The primary objective of this paper is to present the key features of the code nemoFoam, using the ARC reactor design [17,18] as a case study.We have utilized this reactor design to conduct a code benchmark of the neutron and photon module against the Serpent Monte Carlo code.The paper is structured as follows: first, the single modules employed in nemoFoam are described in Section 2; then, a sample case study is set up and presented in Section 3, and the results are presented in Section 4.

Structure of the nemoFoam code
The nemoFoam code has been developed in the framework of the OpenFOAM environment, given its open-source nature and flexibility.It is available at the GitHub repository: https://github.com/NEMO-Group/nemoFoam.
The code has been designed to model fusion systems.For the purposes of this work, it has been tested with the ARC reactor.The main contribution of the nuclear module consists in the evaluation of the volumetric power deposition field due to the interactions of neutrons and photons with the surrounding medium.The thermal-hydraulics module, based on the standard thermal-hydraulic multi-region OpenFOAM solver, implements an additional power source term in the energy equation, which corresponds to the previously mentioned volumetric power deposition field.
Thanks to the modularity of the code, it is expected that it could be used also for fission systems.Future works are foreseen to apply the nemoFoam code to Molten Salt Reactors, where the fuel, as component of the molten salt, circulates within the system, for which many experimental data for validation exist [19].

Model for neutrons and photons
The nemoFoam code implements the transient neutron diffusion multi-group equations: where  indicates the th energy group, and  represents the overall number of energy groups,  is the neutron scalar flux,  is the diffusion coefficient, the ensemble of the   values are collected in the scattering matrix,   is the absorption cross section and  represents an additional source of neutrons,  is the fission spectrum,  is the average number of neutrons emitted per fission,    is the fission cross section,   is the total precursor fraction,  is the precursor half-life and  is the precursor concentration, which evolution equation is implemented in nemoFoam but not shown here, in order to focus on the aspects of the model that corresponds to the physics of the fusion case study.
In fact, in fusion systems, it is possible to neglect the fission and the decay concentration terms in Eq. (1).Therefore, in the proceeding of the article, the equations are reported as they are implemented in the code, even if they not correspond to the physics of the steady-state fusion case study analyzed in this work.
In a multiphysics tool developed for the modeling of nuclear fusion reactors, like nemoFoam, it is mandatory to take into account photons too.In fact, photons contribute for ∼ 30% of the total power deposition [20].Neglecting them will substantially distort the thermalhydraulic results, both in terms of power deposition magnitude and of power deposition shape, since the shape of the photon flux is related to, but different from, the neutron one.
Regarding the photon modeling, the mono-kinetic diffusion equation has been implemented (Eq.( 2)).The adoption of one single energy group is an approximation due to the fact that, at the moment, the nuclear properties used in nemoFoam are calculated by the Serpent Monte Carlo code, which currently does not provide the scattering matrix for photons.A possible solution to this limitation is to evaluate the neutron and photon spectra of the reactor under study using Serpent and then to use them as an input to the NJOY processing code [21,22], which allows to evaluate photon scattering matrices too; however, for the time being, the mono-kinetic approximation is adopted in the code.The resulting equation for the photon flux   is the following: where the subscripts  and  indicate that the term is related to photons or neutrons, respectively.The photon diffusion coefficient   is calculated consistently with the P 1 approximation of transport assuming isotropic scattering: where    is the absorption attenuation coefficient, consisting in the sum of the attenuation coefficients for pair production and photoelectric effect,    is the sum of the absorption and scattering (i.e.Compton plus Rayleigh scattering) attenuation coefficients, and   represents an external source of photons.In this work,   has been used to model the contribution of the so-called secondary photons generated by Bremsstrahlung, annihilation and atomic relaxation, as will be explained in Section 3.2.1.The  produced by neutrons (e.g., after radiative capture and inelastic scattering) are considered in the term        .
Concerning the BCs for the closure of the neutronic and the photonic problem, the albedo BCs n is implemented: where n is the outward normal to the boundary and  is the albedo coefficient, defined as the ratio of the currents entering and leaving the system from the boundary under analysis.In the case in which  = 0 (i.e.no incoming current) the surface is treated as free surface and the condition imposed becomes the vacuum BCs.For nuclear fusion systems, the source of neutrons is located outside of the domain, for this reason the incoming partial current condition has been implemented: where   is the incoming current expected from the plasma, computed, as shown in Section 3.2.2, by means of CHERAB open-source tool [23].
For what concerns the interface condition, the continuity of the fluxes and of the normal component of the currents must be satisfied.Finally, after having computed the neutronic and photonic distributions, the power deposition field is computed by means of the KERMA (Kinetic Energy Released in Matter) coefficients (denoted as ), expressed in eV barn: Being a multiphysics tool, nemoFoam offers the possibility to consider thermal feedbacks on the nuclear properties thanks to the Doppler effects and to the variation of the materials density due to the variation of the temperature (which impacts the macroscopic cross sections).In fact, it is possible to feed the module with several sets of properties for neutrons and photons (each for a certain reference temperature and for the corresponding densities computed at this temperature), enabling the code to perform the required 2D interpolation, as customarily done in multiphysics calculations of fission reactors.

Thermal-hydraulic model
The thermal-hydraulic module of nemoFoam solves the standard transient thermal-hydraulic problem based on the chtMultiRegionFoam solver available in any OpenFOAM distribution, which, for fluid regions, consists of the Reynolds Averaged continuity, Navier-Stokes (RANS) and the energy conservation equations: where  is the fluid density,  is the mean fluid velocity,  is the mean pressure,  is the gravitational acceleration,  is the dynamic viscosity,   is the turbulent viscosity, ⃗   is an additional mean source of momentum,  is the mean internal energy,  is the thermal conductivity,  is the turbulent thermal conductivity, π is the mean viscous stress tensor and   is the volumetric power deposition computed from the nuclear model.Here again the equation are reported as implemented in the standard OpenFOAM solver, even if the conducted analysis is steady-state.The turbulent model selected is the standard k- model [24].
For solid regions, the energy conservation equation is solved: where ℎ is the specific enthalpy,  = ∕  is the thermal diffusivity, defined as the ratio between the solid thermal conductivity () and the solid specific heat capacity (  ).The conjugate heat transfer (solid-liquid) exploits the conservation of the heat flux at the interface: where n1,2 is the normal vector in the direction of the wall of the two regions in contact.

Coupling strategy
In principle, the nemoFoam code can solve the thermal-hydraulic (TH) and the neutrons and photons (NP) equations on two different meshes.The rationale behind this choice is to reduce the computational cost of the simulations and to increase the flexibility of the tool.In fact, the different problems require different mesh refinements: an effective neutron distribution can be obtained on a coarser mesh, at the periphery of the domain, without boundary layers, while it is mandatory to resolve accurately the boundary layer in thermalhydraulics simulations.The necessary mapping of the relevant fields (e.g., power deposition and temperature) from one mesh to the other one is based on the post-processing utility mapFields available in Open-FOAM.Firstly, the code solves the thermal-hydraulic problem with no nuclear power source, then computes the multi-group constants (if the thermal feedback is enabled) and proceeds with the nuclear calculations.
At every ''call'', each module requires the updated field computed from the other one: the NP module requires the updated temperature field and TH the power deposition field.Iterations between NP and TH are performed until convergence is reached for the steady-state calculations (i.e.residuals falls below the value set as limit in the fvSolution dictionary of the case as shown in Section 3.2), respectively, on the neutron fluxes and photon flux, and on the thermal-hydraulic fields.A sketch of the procedure is reported in Fig. 1.

Code-to-code benchmark
The nemoFoam code has been used for the simulation of the steady state case of the ARC fusion reactor, assuming nominal operation.With the aim of assessing the accuracy of the nuclear results, a Monte Carlo simulation has been carried by means of the Serpent code (version 2.2.1) for the same problem.Serpent can be used to analyze fusion systems thanks to its capability to easily import CAD geometries in .stlformat [25] and to perform neutron-photon coupled simulations [26].

Geometry and materials
The ARC reactor is currently in the design phase.The structure of the reactor comprises multiple layers of different materials (e.g., In-conel718 and Beryllium), housing coolant channels and an external blanket tank containing FLiBe.An internal sheet of 1 mm of tungsten is foreseen as plasma facing material, but in the present work it has been neglected, to reduce the computational cost.The geometry used in the simulations represents a 10 • section of the full reactor (Fig. 2), in order to exploit the toroidal symmetry of the problem.In principle, the nuclear problem can be assumed to be quasi two-dimensional, while the thermal-hydraulics is fully three dimensional.To reduce the amount of computational cells, and given the hydraulic arrangement of the inlet and outlet manifolds, the minimum section to be discretized has been found to be equivalent to the 10 • sector considered in the work.A comparison of the volumes of the CAD model used in this work and the literature ones is reported in [20].

Simulation setup
The nemoFoam code has been developed following the structure of chtMultiRegionFoam, the only multi-region solver within the OpenFOAM distribution.For this reason, the setup of a nemoFoam simulation follows the rationale behind the setup of a generic chtMultiRegionFoam case.The setup consists in feeding the module with information on group constants, geometry, boundary and initial conditions.The first two points are addressed in the case sub-folder ''constant '', while the other two in the sub-folder ''0''.The last sub-folder required by the module is ''system'', where information on numerical schemes, and other parameters associated with the solution procedure, can be set.For instance, inside the directory constant (Fig. 3) it is required to have a file with the list of the regions (regionProperties) and for each region a sub-directory (named as the region) containing, the ''polyMesh'' sub-directory for the mesh, a file with the properties of the neutrons groups (nuclearProperties) and a file with the properties of photons (photonProperties).
The boundary and the initial conditions are addressed in the subfolder ''0'' (Fig. 4), where for each region it is required to have a sub-directory (named as the region) containing a file for the temperature field, a file for the photon flux and a file for the external photon source and a set of  files for the neutron fluxes and for eventual external neutron sources, where  is the number of neutron groups.
In the folder ''system'' the standard dictionaries for the control of the simulation, the choice of the numerical schemes and tolerances of the iterative solution must be provided.For the case study, a relative tolerance of 10 −5 has been set on the neutron and thermal-hydraulics variables.

Nuclear properties
Nuclear properties have been calculated using the Serpent Monte Carlo code (for all the details about the Serpent methodology for the generation of the multi-group cross sections, refer to [27]).Multiple Serpent simulations were performed using the same geometry files as in the previous section but at different temperatures (and densities), determined from an initial nemoFoam simulation.The Serpent results have been obtained with a statistical error smaller than 1% (1 sigma) for all the cross sections, despite of some elements of the scattering matrix with higher errors, but which are negligible.This allows to obtain a library of nuclear properties at different temperatures that can be used in multiphysics simulations.Moreover, in the ARC reactor, where neutron multiplication is expected (due to the presence of beryllium), it is necessary to consider the neutron multiplication reactions (n,xn).For  In this case, it is assumed to have two regions and two neutron groups ( = 2).
this reason, the-so called ''reduced absorption cross section'' has been used, which is defined as: Note that, in the case of strongly scattering multiplying media,   can become negative.The reduced absorption cross section can be used in place of the standard absorption cross section, with no changes to the diffusion Eq. ( 1) [28].
The Serpent code generated the group constants for six energy groups, with an energy grid (Table 1) compliant with the fusion neutron energy spectrum: more refined for high energy values and less refined for thermal energy ranges.The choice of the energy grid has been made with the aim of reducing the computational time, while remaining consistent with the fusion spectrum.However, it is not guaranteed that this grid is the best choice, thus, in the future, it will be useful to test also other energy grids and to try to find the optimal one using different techniques, like genetic algorithms [29].Serpent allows to automatically estimate multi-group neutron cross sections, but this same feature is not available for multi-group photon cross sections.Thus, the photon properties have to be evaluated defining specific detectors in the reactor model.This is also the reason why currently it is not possible to evaluate the scattering matrix for photons, since Serpent detectors cannot provide multi-group scattering matrices.The strategy employed in this work to obtain the single-group photon cross sections is to define a detector for each chemical element and for each reaction of interest (i.e.photo-electric effect, Compton scattering, Rayleigh scattering and pair production) in each reactor region for which photon properties are required.The procedure is the one proposed in [30], where the attenuation coefficient for a certain reaction x can be expressed as: The sum is performed over all the elements in each region (  is the density of the th element in atoms/cm 3 ).A similar approach has been used also for the KERMA coefficients   .Particular attention has been dedicated to the definition of the   term in Eq. ( 2).As already stated, this term accounts for the generation of secondary photons from Bremsstrahlung, annihilation and atomic relaxation in a simplified form.This term can be estimated performing two Serpent simulations: the first one is a neutron-photon coupled simulation that allows to obtain the number of photons generated by neutron interactions (e.g.radiative capture and inelastic scattering) and to generate a file containing the photon source points, the second one is a simulation with only photons emitted according to the above mentioned file, which allows to obtain the total number of photons generated in the reactor.Finally,   is computed as the difference between the total number of photons (from the second simulation) and the photons generated by neutron interactions (from the first simulation).Results characteristic of ARC are shown in Fig. 5, underlining the importance of secondary photons to the total balance.However, it is important to underline that using the external source term to model secondary photons is a simplifying assumption, since the generation of secondary photons actually depends on the photon flux   in a specific point of the reactor.To the best of the authors' knowledge, even though this is a simplifying assumption, there is no deterministic code that can be used for the simulation of coupled neutron and photon transport able to take into account secondary photons, even in the nuclear fission community (see [31][32][33]).In the future, a more realistic photonic model will be implemented in nemoFoam, but at the moment this simplification can be considered satisfactory in order to obtain preliminary, yet sufficiently accurate results.Taking all the previous considerations into account, a set of nuclear data has been evaluated for each component of ARC (i.e., the inner vacuum vessel, the cooling channel, the neutron multiplier, the outer vacuum vessel and the breeding blanket), using the ENDF/B-VIII.0nuclear data library [34].Then, the nemoFoam files containing the nuclear properties employed by its NP module have been generated thanks to a Python script and exploiting the serpentTools Python package [35], developed for the post-processing of the Serpent outputs.

Boundary conditions
BCs are established depending on the type of the surface under consideration: free surface BCs for the outer blanket tank (external surface of the geometry), interface BCs for internal patches within different regions, and incoming current BCs for plasma-facing patches (inner surface of the geometry).Particular attention is given to the definition of incoming neutron currents, which is challenging due to the complex geometry of the ARC blanket.Neutron currents were computed by means of the CHERAB as shown in Fig. 6.CHERAB is  based on the 3D Monte-Carlo inverse ray-tracing code Raysect [36] and was used in the past specifically for fusion applications, like the evaluation of the radiative heat load distribution on the EU-DEMO first wall [37].In a similar way, in the present work CHERAB evaluates the neutron current impacting on the first wall, which then can be used to evaluate the BCs of the problem in nemoFoam.To this aim, it is sufficient to define the mesh geometry of the first wall and the plasma source (in this case, a simplified annular source emitting isotropically mono-energetic neutrons at 14.1 MeV).
However, CHERAB can be used only to evaluate the incoming current in the highest neutron energy group, since the neutrons generated by the plasma source belong to this group.For the other groups, one can assume unitary values for the albedo coefficients, since it is reasonable to expect that the number of neutrons belonging to the slower groups that enter the first wall is almost equal to the number of neutrons belonging to the same group entering the plasma chamber (i.e., from the first wall) at every location.Actually, the albedo coefficients are exactly equal to 1 in the case of symmetric plasma chambers (e.g. with a circular shape), while in the case of asymmetric geometries (like ARC) this equivalence is not fully satisfied.An alternative solution is to directly evaluate incoming currents and/or albedo coefficients using Serpent, thanks to the definition of a series of surfaces along the plasma chamber profile.This solution obviously reduces the geometric detail of the BCs, but it does not require the use of another software like CHERAB.The albedo coefficients evaluated using this method are similar to 1, showing that the previous assumption is reasonable.Moreover, the results obtained with nemoFoam using the two different sets of BCs are similar, suggesting that the user can arbitrarily select the preferred one.The latter method has been used to evaluate the BCs for the photon diffusion equation at the plasma-facing interface, since it was not possible to define a photon source in CHERAB due to a lack of knowledge.In fact, the evaluation of the number of photons emitted inside the plasma is much more complicated than the neutrons one and it is out of the scope of this work.
For the fluid domain, FLiBe enters from the inlet patch at 2.0 m/s and at a fixed temperature of 800 K, at this location the Reynolds number at the inlet is in order of  = 6 ⋅ 10 4 .The outlet patch uses the inletOutlet type, allowing outflow with zero gradient condition.This choice is based on the absence of circuit information and on the elongated position of the outlet away from the tank exit.Interface patches are set as turbulentTemperatureCoupledBaffleMixed for thermal treatment and as walls for hydraulic cases.The other patches, representing the external tank surface, are adiabatic walls.Patches resulting from domain cutting are designated as symmetry boundaries.

Mesh generation
The nemoFoam code requires the definition of meshes in the poly-Mesh format, so meshes must be generated or converted from the meshing tools available for OpenFOAM.The meshes used in this study have been generated by means of cfMesh [38], primarily designed for meshing single domains, producing separated meshes of individual regions.While this approach offers flexibility in parameter selection for each mesh, it requires a careful redefinition of how patches communicate between themselves before the simulations.To address this aspect, contiguous patches were redefined through the OpenFOAM standard utility createPatch, based on the information from the createPatchDict file found in each region sub-directory within the simulation case system directory.The type of patch selected has been mappedWall, with nearestPatchFaceAMI as sample mode in order to account for non conformal meshes at the interfaces.The choice of the sample mode was influenced by the nature of the sampling, enabling OpenFOAM to extract values from the nearest cell face on the communicating patch as defined in the dictionary.This approach accommodates non-conformal patches (where the grid on the two patches may differ) and employs Arbitrary Mesh Interface (AMI) interpolation, allowing sampling across disconnected but adjacent mesh regions.In figures Figs. 7(a) and 7(b) two zooms of the first 15 cm of the blanket with the plasma standing on the right side of the mesh are reported, where the major differences in the meshes are present.The advantage of having the possibility to use different meshes for the different problems becomes clearer with Table 2, where the total number of cells is reduced by ∼ 22% in the neutronic mesh.Note that in the thermal-hydraulic mesh the two regions ''Breeding Blanket'' and ''Coolant Channel'' are merged in one single domain, since they are the same under the hydraulic point of view.

Results
The nuclear performance of the code has been assessed against Monte Carlo simulations.The resulting neutron flux, photon flux and power deposition are compared.Moreover, the neutronic spectrum is calculated and compared to the Monte Carlo reference, in order to assess the accuracy of the diffusion approximation.Regarding the thermal-hydraulic module, the flow field and the temperature map are presented at the end of the section, although no benchmark data is available for them.

Neutronics
The neutron and the photon fluxes are the fundamental quantities of the analysis, since the power deposition is computed from those.The nemoFoam simulations for the NP part were run in parallel with 24 computational cores, imposing as tolerance criteria the value of 10 −5 for both the photon flux and the group-wise neutron fluxes, requiring 2.5 h (60 core-hours) in total, for a total of 1100 iterations mainly required by the first call of the module, which had a constant initial condition of 10 13 n∕cm 2 ∕s.The qualitative agreement of nemoFoam (see Fig. 8) with the neutron and photon flux maps obtained with Serpent 2.2.1 (Figs. 9 and 10) is fairly good.The Serpent simulation has been run in the coupled neutron-photon transport mode with a total number of 10 8 neutron active histories divided in 100 batches using 40 computational cores, requiring 14 h (560 core-hours), meaning that the nuclear module of nemoFoam is 9× times faster.It is important to observe that, despite the large number of neutron histories, the statistical error obtained with Serpent is still quite large in the peripheral regions.It is important to notice that Serpent, being a Monte Carlo code, is able to simulate the transport of neutrons and photons in the vacuum (i.e., the plasma chamber) too.This is the reason why in Figs. 9 and 10 the flux is non-zero in the plasma chamber.On the other hand, nemoFoam, being a diffusion-based code, is not able to do this, but can fairly reproduce Serpent results in the remaining (and more interesting) part of the domain thanks to the definition of proper BCs (using either CHERAB or Serpent).The total amount of power deposited in each region was calculated by nemoFoam, integrating the volumetric power field computed by the code.The relative error between the nemoFoam results and the Serpent results are below 10% (except for the photon power deposition inside the tank, see Table 3).These results confirm that the impact of photons is non-negligible, justifying the implementation of the photonic module in the nemoFoam code.The region in which the worst results are obtained is the tank, where the photon power deposition is strongly underestimated.This can be explained by the mono-kinetic treatment of the photon diffusion and to the fact that only a single homogeneous value of   (i.e., the secondary photons production term) has been evaluated with Serpent for the whole tank region.This means that, according to nemoFoam, the density of secondary photons generated at the periphery of the tank and near the outer vacuum vessel (where it is actually much higher, due to the higher and faster photon flux), is considered equal, leading to an underestimation of the photon power deposition in the tank.Another interesting response is the neutron flux spectrum, which describes the energy distribution of neutrons.Neutrons are initially generated with high energy inside the plasma, at 14.1 MeV, and subsequently slow down interacting with the nuclei of the media through which they propagate.The specific characteristics of this slowing-down process depends on various factors, including the properties of the medium, the geometry and the arrangement of the system.Defining the spectrum of neutrons per unit lethargy (  ) as: where  0 is the upper limit of each energetic group and  the lower one, it is possible to compare the solutions obtained by the nemoFoam code and the Monte Carlo Serpent code from a qualitative point of view.
As shown in Fig. 11, the agreement is very good for the less energetic groups, while the major differences regard the two most energetic groups.However, the final impact on the power deposition is limited, as shown in Table 3.
Finally, the distribution of the total power in the radial direction, on the equatorial section of the ARC reactor, is reported in Fig. 12, showing again that the results computed by nemoFoam are consistent with the ones obtained with Serpent.

TH results
In this section the computed flow field and the temperature map are presented on a 2D section of the ARC reactor.The cross-sectional area    of the channel is significantly smaller than that of the tank and of the inlet patch.Therefore, with an inlet velocity of 2.0 m/s, the channel reaches a maximum velocity of 15 m/s, creating a highly turbulent flow with effective heat removal from surrounding solid materials.The limit of the legend of Fig. 13  the flow bifurcates into two streams after the inlet, with comparable mass flow rates.Interestingly, the combined effects of different crosssectional flow areas and localized losses at elbows produce similar global hydraulic resistance in two seemingly distinct pathways.The velocity of the FLiBe in the inner and outer branches depends on the channel's cross-sectional flow area, which is influenced by the distance from the torus' rotation axis.Particular attention is required at the elbows, where the flow changes direction and accelerates, potentially creating hot spots in low-velocity regions and corrosion concerns in high-velocity locations.Downstream of the tank connection port, due to a sudden change in cross-sectional area, the flow rapidly transforms into a jet-type flow, leading to a decrease in velocity.In the blanket tank, the flow field is characterized by low velocities, impacting heat transfer capabilities.However, this is compensated by a volumetric power deposition that is one order of magnitude smaller than the one in the channel, even just 20 cm from the Inconel718 surface.Regions of relatively higher velocities occur near divertors and the outlet pipe.The FLiBe exiting the port bifurcates into two streams -one through the outboard region and one through the inboard.The region after the port, where the fluid expands, exhibits a complex structure requiring further investigation, likely involving specific turbulence modeling to keep into account the high-Prandtl nature of the molten salt [39], which will be addressed in future code releases.Analyzing the temperature map of the coolant channel (Fig. 13(b)), the inboard flow experiences a lower temperature gain with respect to the outboard flow.This is attributed to the longer distance and lower velocity of the outboard flow, resulting in a hot spot near the last elbow on the outboard segment.Examining the temperature distribution in solid layers, the outer Vacuum Vessel is the hottest, reaching a peak temperature of 890 K, within the thermal limits of Inconel718.The thermal-hydraulic results presented serve as a demonstration of the tool capabilities.However, it is important to note that key aspects such as MHD and high Prandtl number turbulence have not been taken into account.Therefore, it is mandatory not to consider them as final, as they cannot represent the operating scenario of the reactor.

Conclusions and perspective
Given the complex nature of nuclear reactors design, addressing multiple physical phenomena such as thermal-hydraulics and neutronics often involves using various specialized codes.This paper introduces nemoFoam, a multiphysics code developed within the OpenFOAM environment, that deals with both neutronics and thermal-hydraulics for the modeling of ARC-type fusion reactors, using two distinct modules.The nuclear module implements the neutron diffusion multi-group equations.It also takes into account the significant contribution of photons, which can affect thermal-hydraulic results, by means of the mono-kinetic diffusion equation, due to current limitations in generating photon scattering matrices.The thermal-hydraulic module, based on the standard thermal-hydraulic multi-region OpenFOAM solver, incorporates a power source term derived from the volumetric power deposition field obtained from the nuclear model.The coupling strategy involves solving the nuclear and thermal-hydraulic problems on separate meshes, with data exchange and convergence achieved through iterations between the two modules.To demonstrate the capabilities of nemoFoam, the paper presents a code-to-code benchmark simulation using the ARC fusion reactor as a test-case.The results are compared against Monte Carlo simulations using the Serpent code, showing good agreement in terms of neutron and photon fluxes, power deposition, and neutron flux spectrum.The code demonstrates its capability to provide reliable results with relative errors below 10% (except for the FLiBe tank).The neutron flux spectrum comparison between nemo-Foam and Serpent indicates good agreement, particularly for lower energy groups.Concerning the thermal-hydraulic analysis, nemoFoam exhibits the ability to capture the complex flow patterns and temperature distribution in the ARC reactor.In conclusion, nemoFoam has demonstrated to be a flexible and modular code capable of addressing both steady-state neutronics and thermal-hydraulics in fusion reactors.The benchmark results for the ARC reactor show the code's reliability and its potential.Future developments will include more refined photon models and an MHD solver for the thermal-hydraulic part, including dedicated turbulence modeling for high-Prandtl fluids.

Fig. 2 .
Fig. 2. Scheme of the CAD produced, with a focus on the inlet, outlet and connection port from the coolant channel to the tank.The different material layers are reported on the right.

Fig. 3 .
Fig. 3. Structure of the constant folder for a nemoFoam simulation of the nuclear module.In this case it is assumed to have two regions.

Fig. 4 .
Fig. 4. Structure of the 0 folder for a nemoFoam simulation of the nuclear module.

Fig. 6 .
Fig. 6.Incoming current distribution on the plasma facing surface.Currents are computed with CHERAB and subsequently mapped onto the OpenFOAM mesh.

Fig. 11 .
Fig. 11.Ratio of the neutronic spectrum computed by Serpent with respect to the results obtained with nemoFoam.
(a) is set to 4.0 m/s to allow to see the actual flow patterns inside the blanket.Examining the coolant channel,

Fig. 12 .
Fig. 12.Total volumetric power deposition on the mid-plane section of ARC reactor in radial direction at the equatorial location.The error bars of the Serpent results are not visible because the relative error is too low.

Fig. 13 .
Fig. 13.Thermal-hydraulics results showed on the mid-plane section of ARC reactor.

Table 1
Energy grid employed for the group constants generation in the ARC model.

Table 3
Comparison of the computed power deposition between the nemoFoam code and the Monte Carlo Serpent code.