Lattice Boltzmann Method Applied to Nuclear Reactors—A Systematic Literature Review

Nuclear engineering requires computationally efficient methods to simulate different components and systems of plants. The Lattice Boltzmann Method (LBM), a numerical method with a mesoscopic approach to Computational Fluid Dynamic (CFD) derived from the Boltzmann equation and the Maxwell–Boltzmann distribution, can be an adequate option. The purpose of this paper is to present a review of the recent applications of the Lattice Boltzmann Method in nuclear engineering research. A systematic literature review using three databases (Web of Science, Scopus, and ScienceDirect) was done, and the items found were categorized by the main research topics into computational fluid dynamics and neutronic applications. The features of the problem addressed, the characteristics of the numerical method, and some relevant conclusions of each study are resumed and presented. A total of 45 items (25 for computational fluid dynamics applications and 20 for neutronics) was found on a wide range of nuclear engineering problems, including thermal flow, turbulence mixing of coolant, sedimentation of impurities, neutron transport, criticality problem, and other relevant issues. The LBM results in being a flexible numerical method capable of integrating multiphysics and hybrid schemes, and is efficient for the inner parallelization of the algorithm that brings a widely applicable tool in nuclear engineering problems. Interest in the LBM applications in this field has been increasing and evolving from early stages to a mature form, as this review shows.


Introduction
The need for achieving energy sustainability and reliability, as well as cutting air pollutants, greenhouse gases, and ozone-depleting substances, could require nuclear energy to be part of the energy mix, since renewable energies (such solar, wind, etc.) are inherently intermittent [1][2][3]. The penetration of the intermittent renewable energies in the power system demands an increased flexibility, a characteristic typical of nuclear energy that can help to achieve a higher renewable energy penetration [4]. Moreover, with the advent of modern and foreseeable technologies, such as fast reactors and seawater uranium extraction, nuclear energy could also be considered as substantially "renewable" and potentially inexhaustible [5,6]. It has even been suggested that "The support of a thermal reactor fleet in the mix will in all cases be needed until the end of the present century an even beyond, independently of the to predict the behavior of the systems in a nuclear power plant using highly accurate simulations can help to increase confidence concerning safety from stakeholders.
The validation of computer codes and calculation methods in nuclear power plants must fulfill a series of criteria required by a framework of international organizations (e.g., Nuclear Energy Agency, International Atomic Energy Agency, etc.) [17,21].
As will be shown in this paper, the Lattice Boltzmann Method (LBM) has great potential, and its formal validation will require more and more studies and research. In fact, in the last years, fluid dynamic simulations can rely on the LBM, a relatively new simulation technique that is based on streaming and collision processes of probability density functions on a lattice (a graphical representation in Figure 1, see the Appendix A for details of the method). Instead of directly solving the Navier-Stokes (NS) equations, this method solves the Boltzmann equation by a discretization procedure and using an approximation for the collision operator. Some peculiar characteristics of the LBM, in comparison with the most common computational fluid dynamics (CFD) Navier-Stokes based techniques, are presented in Table 1 [22].  Table 1. Comparison between the LBM and Navier-Stokes (NS) based methods (adapted from [22]). Presently, there is no literature review that summarizes the research works that apply the Lattice Boltzmann Method in the nuclear engineering field. An example of a review on simulation methods coupling neutronic and thermal-hydraulic analysis of nuclear reactors can be found in [23], with no mention of the use of the LBM, which represents a novel perspective, as will be duly explained within this work.

Feature NS Based Methods LBM
Identifying the applicability of the Lattice Boltzmann Method in the nuclear reactor engineering field is important not only for the original approach that proposes this numerical method to the various problems, but also for the very high capabilities of the method to handle complex geometries, and for being inherently parallelizable, making the LBM an optimal multiphysics tool to be applied to nuclear reactors simulations.
The purpose of this paper is to give an overview of the framework and the application of the Lattice Boltzmann Method in nuclear engineering, to identify which problems in the nuclear engineering field have been addressed with the LBM, how the LBM has been adapted to be applied in nuclear engineering, how results were validated, and to explore some trends in the research of LBM applied to nuclear engineering. To this aim, a systematic literature review (SLR) was carried out, identifying the relevant literature items in the scientific databases and extracting the relevant information. The specific nuclear engineering problems were identified and the included items were classified into two categories (CFD applications and neutronic applications). The adaptation of the LBM to solve nuclear engineering problems was then identified (such as grid type, collision model, turbulence model) and the validation techniques were subdivided in experimental, numerical, or analytic. The literature found was sorted by date to provide an overview of the research trends. Finally some interesting and specific features and/or conclusions emerging from each study were described.
After the introduction, in the Section 2, the methodology adopted for the systematic literature review is described. In the Section 3, the relevant research items are analyzed, in the Section 4, a discussion concerning the results found is presented. Finally, in the last section, a summary of the conclusions is carried out. In Appendix A, a brief description of the LBM can be found. Results of the review show that the LBM is a flexible numerical method capable to integrate multiphysics phenomena and hybrid schemes, and it is found to be very efficient for its intrinsic parallelization of the algorithm, which makes it a widely applicable tool in different nuclear engineering problems, ranging from CFD to neutronic applications.

Methodology
The systematic literature review is a particular kind of review commonly used in the medical field, but useful also in other research fields, such as engineering, data sciences, or environmental studies (see for example [24,25]). Following the general recommendations from [26,27], the systematic approach to the recent literature includes the definition of research questions (RQ), according to the research objectives, inclusion criteria (IC), exclusion criteria (EC), and a search strategy (as a part of the methodology), data extraction (to present coherently the results) and, finally, the conclusions. According to this methodology, the following RQ, IC, EC, search strategy, and data extraction were defined: To present the search strategy, it is common to include the Preferred Reporting Items for Systematic Reviews and Meta-Analyses (PRISMA) [28] flow chart presented in Figure 2. This flow diagram resumes the identification, screening, and eligibility steps followed to select the research items included in this SLR. Three databases were consulted: Web of Science, ScienceDirect, and Scopus. Subsequently, more items hand-searched were added (including some pertinent conference papers).
The adopted search string in the accessed databases was derived from the logical string: ("nuclear" OR "nuclear reactor" OR "reactor" OR "nuclear energy" OR "Neutron") AND ("LBM" OR "Lattice Boltzmann" OR "Lattice Boltzmann Method"): in each database, the search string was written in a congruent form (sub-step identification I), and then the results where filtered applying the inclusion criteria (IC) and exclusion criteria (EC) (sub-step Identification II).
To accomplish the data extraction, the included items were manually filtered and classified into two categories, i.e., CFD applications and neutronic applications, then the specific nuclear engineering problem was identified (RQ1). The adaptation of the LBM to solve the nuclear engineering problem (such as grid type, collision model, turbulence model, etc.) was then identified (RQ2), and the validation techniques (RQ3) were subdivided in: experimental data comparison, comparison with other numerical methods, and performance analysis between LBM variants. The literature found was sorted by date to provide an overview of the research trends (RQ 4).

Results: Research Items on LBM Applied to Nuclear Reactors
Firstly, the main description of the relevant literature found will be presented, and then the analysis will focus on the two general application fields, as previously described, i.e., computational fluid dynamics (including isothermal and thermal flow) and neutronics (including transport and diffusion), each one in a dedicated subsection. It has to be underlined that all the figures included in this review are original artworks obtained by our implementation of the LBM, replicating the results reported by the cited authors (using the open source Palabos library [29,30] for CFD, and in-house codes developed in MATLAB for neutronics).
In total, 45 records are included in this review and categorized regarding the main topic: (25) records for CFD applications and (20) records for neutronic applications. Figure 3 represents the dates of publication for all the records included.

Computational Fluid Dynamics
Before presenting the general results of the literature, a general review of the basic concepts about CFD will be presented. The use of the LBM as a CFD simulation tool is the straightforward use of this numerical method, calculating the time evolution of the density, and the mean velocity field in each grid point from the first and second moments of the distribution function as referred in the Annex A. In this way, LBM has been used to study, for example, the coolant mixing in the reactor vessel. The simple Bhatnagar-Gross-Krook (BGK) scheme described in the Appendix A is isothermal, and for this reason, to apply the method with a realistic approach to the reactor coolant circulation problem, and simulate how the particular mixing conditions and flow regime can enhance the heat transfer coefficients, the traditional isothermal LBM must be modified. This is usually done by choosing one of the following four different approaches that improve the LBM to a "thermal LBM": • Multi-speed: is an extension of the classical LBM model with a higher-order expansion for the equilibrium distribution, this method can be unstable and can have a lack of precision in heat transfer calculations. • Passive scalar transport: a numerical scheme can be adapted to the LBM for the advection of the thermal field (a comparison of some schemes can be found in [31]).

•
Multi-field (secondary distribution): the transport of a secondary energy distribution function is calculated with the Boltzmann equation, and the temperature is obtained from the integration of this secondary distribution (first integral moment). • Hybrid: refers to using a Finite Element Method (FEM) or Finite Volume Method (FVM) scheme to calculate the thermal field from the energy equation and coupling the results with the LBM fluid simulation; a comparison of hybrid and multi-field schemes in Natural Convection can be found in [32]. The Hybrid Finite Volume Lattice Boltzmann Method (FV-LBM) scheme is as accurate as a classical FVM numerical method [33].
The multi-field (or secondary distribution) approach implies the inclusion of a new energy distribution g i (r, t) overlapped to the density distribution f i (r, t), with a correlated new collision operator and a second relaxation time τ 2 . This gives the advantage of having a direct way to calculate the thermal field. The discretized Boltzmann equation for this secondary energy distribution is: The equilibrium function for the thermal (energy) field depends on the temperature T; the mean velocity u adopts the following form: The multi-field approach and the passive scalar transport approach can be extended, not only for the thermal field, but also for other scalar fields, such as chemical species. The second distribution g i (r, t) can be overlapped to the density distribution and transported with the fluid, allowing the simulation of deposition processes and some other multi-physics problems. The thermal field in each point T(r, t) can be obtained by the first moment of the secondary distribution as: It is interesting to note that the secondary distribution can be expressed in another grid quadrature DxQz different from the DxQy used for the density field. From the Chapman-Enskog expansion analysis (used to recovery the macroscopic equations from the discretized Boltzmann equation), it can be found that the secondary relaxation time τ 2 is related to the thermal diffusivity α of the fluid as: For thermal analyses, in addition to the modified LBM, it is commonly assumed that the Boussinesq hypothesis for the convective motion of the fluid (in particular natural convection), assumes that in natural circulation phenomena, the variations in density (governed by the volumetric expansion coefficient of the fluid β) are only related to temperature changes in the fluid and affect the buoyancy force term. This force term is added to the density distribution and can be written as: In this way, the thermal LBM can be applied to simulate heat transfer processes as the natural convection. In [34], a study of the Rayleigh-Bénard convection, including thermal radiation, is presented; here, the developed formulation for the Radiative Transfer Equation is applied in transient (time-dependent) convection, in the presence of volumetric radiation in a rectangular cavity containing an absorbing, emitting, and scattering medium. This new algorithm can be extended to complex geometries with a possible application in nuclear engineering problems.
In diverse studies, turbulence models are used together with the LBM scheme, such as the Large Eddy Simulation (LES), which allows filtering the turbulence effects, selecting a threshold for the vortices resolved. These kinds of models are easily adapted to the LBM and enhance the possibility to simulate fluxes with high Reynolds numbers that are of practical interest in nuclear technology applications. In [35], a simulation of a turbulent flux into a pipe with several geometries (circular and square cross-section) with two grid models is presented: to validate the scheme, the wall fluid interaction was analyzed and compared with experimental data. A relevant result points out that the grid model D3Q27 matches better than the D3Q19 in both geometries. Other interesting results are found in [36], where the LBM framework is adapted to simulate a high Rayleigh number natural convection in two coaxial cylinders, in a turbulent regime using a static LES-Smagorinsky turbulence model, which can be applied to a cavity with the internal volumetric heat source as a simplified representation of a nuclear reactor.
In the following, a description of the main literature items regarding the direct use of LBM as a CFD simulator applied to nuclear engineering problems can be found, and in Table 2, these findings are summarized.  Rod bundle flow in a Pressurized Water Reactor (PWR)-VVER-440 type (with a triangular array of fuel rods)-was simulated by Házi and Mayer in 2005 [37] using an isothermal LBM as a Direct Numerical Simulation (DNS), and the results are compared with an LES-LBM turbulent model. The objective in this kind of simulation is to know the mixing process of the coolant around the fuel rods, finding, with both methodologies, LBM and LES-LBM, pulsation in the main flux, and a secondary flow pattern that can be strongly dependent on the exact geometry and the Reynolds number value.
In subsequent work, Mayer and Házi, in 2006 [38], simulate the flux around the fuel rods in a similar triangular symmetry, using periodic boundaries in a hexagonal three-dimensional (3D) primitive cell. They compare the performance in laminar and turbulent regimes for two grids (D3Q19 and D3Q27) using the BGK operator, finding that, in a turbulent regime, only D3Q27 reproduces the expected flux characteristics. Mayer et al. in 2007 [39] have simulated a similar geometry in the turbulent regime, using only D3Q27 and an LES turbulence model with an isothermal LBM framework to test the grid size independence of the solution; obtaining periodic oscillations (flow pulsation) in the axial and azimuthal velocities that agree with experimental data (but with a lower periodicity). Specifically, the mean axial velocity agrees with the reference values, and the secondary flow is obtained with some cells recirculation, they found some systematic deviations in the normal Reynolds stresses in the near-wall region. In conclusion, they validated the use of LES-LBM D3Q27 to study the turbulent regime in bundle rods geometries.
H. Chen et al. in 2007 [40] used the LBM to study corrosion mechanisms in structural materials of nuclear plants that use lead-bismuth eutectic as coolant. The active control of oxygen in the LBE coolant can be a strategy to create a thin protective oxide layer that reduces the corrosion process, which is a relevant topic related to safety and lifespan of the plant. In this study, the LBM is modified using the secondary distribution approach to simulate the heat and mass transfer: the authors use a secondary field for temperature calculations and a third field for the oxygen mass transfer and another one for iron mass transfer. They found that through the control of oxygen concentration in the coolant, the corrosive effect diminishes, and a thin protective oxide layer is formed; the natural convection can be implemented as a way to obtain a better mix of the oxygen; in particular, this simulation compares corrosion process under three different thermal boundary conditions (BC) and, consequently, three types of convective motion (regarding the number of vortices).
Y. Chen et al. in 2007 [41] expanded the previous research taking into account forced convection instead of natural convection and tested some inlet-outlet positions for the lead-bismuth eutectic coolant flow, finding that all the strategies enhance the oxygen transfer, and that the best inlet position is at the top of the cavity. Fan et al. in 2009 [42] simulated extreme conditions in a PWR vessel, when a cold water injection from a High-Pressure Injection (HPI) system flows from the cold leg tubes to the upper side of the reaction chamber, and a pressurized thermal shock with fracture risks is possible: these extreme conditions are of remarkable interest because they can represent a risk after a Small Break-Loss Of Coolant Accident (SB-LOCA) when an emergency High-Pressure Injection System (HPIS) works in absence of coolant circulation as an emergency system. The cold plums of water injected into the vessel by the HPIS represent a risk due to the high-temperature gradients; in this study, the cold plumes are simulated by a 3D model of the reactor using LBM, and include thermal effects. Details on LBM implementation on a graphics processing unit (GPU, hereafter) are provided in [43], where the strategies to allocate the relevant cells in the GPU cores (when a complex geometry is simulated) are highlighted.
Park in 2011 [44] developed a three-dimensional LBM simulation of a Very High-Temperature Reactor (VHTR) lower plenum. The simulation was isothermal and represented the turbulent flow in a fuel assembly. The turbulent flow phenomena in a scale model of the lower plenum of a typical prismatic VHTR was performed by two turbulent models: wall adapting local eddy viscosity (WALE) and a Smagorinsky-based LES. The velocity contours for turbulence flows and streamlines were compared, qualitatively, with benchmark data, and with results from a commercial NS-based CFD tool (Ansys Fluent), symmetric solution. Moreover, flow dynamics with the same vortex characteristics were found, as expected.
The thesis of Li (2011) [45] is a contribution concerning the crud formation, a corrosion phenomenon that affects the performance of PWR by deposition of several compounds changing its geometry and heat transfer characteristics (this problem is also known as Axial Offset Anomaly and can affect the fuel rods and other elements in the reactor). A thermal LBM framework using a secondary distribution for thermal field and species transport was proposed; the model simulates, effectively, the crud formation and estimates the effects of the crud deposition on the clad temperature and the ion concentration.
Park et al. in 2012 [46,47] presented results from a deposition process in a Reactor Building Floor Weir as an option for particle deposition after a LOCA. The deposition of several materials in this safety device was simulated for different geometries using 2D models, and the particle sedimentation process around the weirs (barriers) was studied analyzing the velocity fields of the fluid near the weir as an indicator of sedimentation probability, testing different parameters of the weirs as height, inclination, and the number of weirs; the study shows the potential for emergency control of the Reactor Building Floor Weir using relatively tall, double weirs separated by "one or two heights between them". Tamura et al. in 2012 [48] applied the Finite Differences LBM (FD-LBM) to simulate the flow acoustic resonance present in the steam dyer pipes of a Boiling Water Reactor (BWR), which can generate fatigue and malfunction due to the sound-generated vibrations during extended power uprate operation. The authors found that the LBM clarify the resonance mechanism, starting with vortex generation in the safety relief valves, and generating an acoustical wave that travels in the main dryer steam lines. The model was validated by experimental data and the proposed scheme was optimized to reduce computational time using D3Q13 instead of D3Q15, and a modified probability function instead of the distribution function, which considers the distribution function and the equilibrium function at the same time.
Tiftikçi et al. in 2013 [49] developed a study concerning a 3D simulation of flux in a 2 × 2 rod bundle array of a PWR, in laminar regime, with no turbulent or thermal model, showing the capability of the LBM to handle velocity and pressure calculations in complex geometries (including some geometrical details as mixing vanes) and opening a way to further development as a turbulence model implementation.
Carrasco in 2013 [50] applied the LBM to a coolant flow in the lower plenum of a PWR core. A simplified 3D coaxial geometry was simulated and a successive grid refinement was applied to search complex flow phenomena and instabilities. To overcome the computational effort for the calculations of high Reynolds number flows, a multiscale approach was developed to increase resolution by interpolating only in the interesting zones, while a coarser grid was used in the relatively less-complex geometry.
Park, in 2015 [51], presented the simulation of flow on rod bundles in turbulent regime (2D geometry, using primitive cells and periodic boundary conditions) applying a Multiple Relaxation time (MRT) LBM, with an LES-Smagorinsky turbulence model. The comparison of its results with experimental and CFD benchmark data showed a good agreement in velocity profiles, and reproduced some phenomena as flow oscillations; the need for a grid refinement near the walls to increase the precision of the results was highlighted, as well as the computational time reduction if compared to other CFD techniques. In Figure 4, a similar 2D simulation is reproduced, with the primitive cell representing the VVER-440 reactor geometry (with a triangular array of fuel rods, represented as one central circle and four quarter circles in the corners) for a high Reynolds number (equal to 20,000). The effect of the LES model to filter the turbulence and present smoother results is evident: the simulation shows how the main flux is positive in the x-direction, but also the y-velocity field is linked to recirculation and vortex presence in each cell. Tiftikçi and Kocar in 2015 [52] simulated a three-dimensional isothermal flux (in particular the vortices) around the fuel rods, using hexagonal primitive domains with seven rods, and comparing two turbulence models, i.e., LES (Smagorinsky) and Very-Large Eddies Simulation (VLES) (k-ε). They developed the Palabos implementation of the VLES (k-ε) turbulence model, and, to validate the results, they compared the calculated friction factor with experimental values and found that the VLES model drives to better results than the LES model. This work shows how the LBM can simulate turbulent regimes around complex geometries: in this case, not only the bare rods were simulated, but also the rods with a helical wire bundle separator enhancing the turbulence and the mixing process.
Tiftikçi and Kocar in 2016 [53] continued their research comparing thermal LBM with FVM in a VVER-440 triangular rod array geometry: using a reduced triangular primitive cell (instead of a hexagonal one) they simulated the thermal fluid dynamic flow around the rods with and without rod spacers; by this approach, it is possible to obtain heat transfer coefficients and pressure loss estimations. The LBM showed a good agreement with the FVM results, and it highlighted the need for a grid refinement around the rods to capture the circular geometry. The same authors presented in [54] the application of the isothermal regularized LBM to calculate the turbulent flux around a triangular array of rods. Here, an extensive set of turbulence models was compared with FVM simulations and with experimental results; authors conducted a sensitivity analysis of the results: testing grid size and turbulence model parameters (Smagorinsky constant, filter size), and testing the performance of the grid types D3Q15, D3Q19, and D3Q27. They found that the simplest LES turbulent model does not reproduce near-wall shear stresses results as well as VLES, with a better and consistent performance in D3Q19 and D3Q27, while D3Q15 fails to reproduce the main turbulent features. In 2016, Tiftikçi and Kocar [55] presented a thermal turbulent flow simulation in a prototype Fast Breeder Reactor (FBR) geometry, with uniform heat flux conditions and a hexagonal 7-wire rod primitive cell; the use of metal-liquid as coolant (sodium) implies a low Prandtl number and a relevant role of conduction in the heat transfer process. The study compares two turbulence models, i.e., LES-WALE and VLES, validating the obtained results by comparison with experimental data. Results showed that several thermal-fluid-dynamic variables can be obtained by this kind of simulation and demonstrated the great LBM capability to simulate, at the same time, thermal effects and fluid dynamics at high Reynolds numbers. Additionally, this study made evident a better performance of the VLES turbulence model.
In order to have a visual representation of the application of the secondary field thermal model in nuclear engineering problems, in Figure 5 the results from a thermal fluid-dynamic simulation of a coolant (at 293 K) around a single hot rod (at 600 K) into a channel with a constant flow velocity at the inlet in the x-direction are represented. The Palabos library was used in our simulations with a D3Q19 lattice for CFD and a secondary thermal field approach with a D3Q7 lattice and LES-Smagorinsky turbulence model. It is interesting to see how the flow characteristics are linked with the thermal field, and that the heat is transferred and mixed following the characteristic vorticity of the turbulence. An interesting simulation of flux around and inside a semi-porous thermal-active square medium was presented by Vijaybabu et al. in 2017 [56]. This could be considered as a parallel approach to simulate and understand the control rod action in the reactor vessel. In this study, they used the secondary distribution thermal LBM for a flow in a 2D domain, testing low Reynolds number fluxes and several permeability Darcy numbers for the active zone (low permeability levels mean a solid behavior). Results showed the effects of the permeability in the heat transfer performance; the porous approach can simplify computational effort for thermal simulations of rod bundles arrays with control rods.
Wahba, in 2017 [57], developed a study using the LBM to simulate natural convection in simple geometries with a low Prandtl number liquid medium (e.g., mercury, sodium, potassium), which could be used in liquid metal fast nuclear reactors as a coolant given its low viscosity and high heat capacity. This study validates the use of a secondary distribution thermal LBM (D2Q9 for fluid, D2Q5 for thermal field) in the simulation of a low-Prandtl number fluid, and finds that for the same Rayleigh number, lower Prandtl fluids weaken the effect of convection.
In the work of Tiftikçi and Kocar (2018) [58], a simulation of a Light Water Reactor (LWR) unit cell subchannel, a typical rod bundle square array of a PWR, can be found. In this study they compared three turbulence models, i.e., LES-WALE, VLES (k-ε), and VLES (k-ω), implementing all models in the open-source library Palabos. A comparison was also made for the grid resolution influence on the three models. Following previous results, they only used the D3Q27 lattice for isothermal turbulent simulations, finding that the VLES models present a better approximation of experimental results with a relatively coarse grid size, and, subsequently, a reduction in computational costs if compared to WALE.
Tamura et al. in 2017 [59] and 2018 [60] presented a simulation of a cold trap in a FBR to control hydrogen and oxygen concentration when liquid sodium is used as a coolant. The cold trap is simulated in a three-dimensional model using a modified Low Reynolds number LBM for the Impurity Deposition of particles (LR-LBM-IDM); this model replaces the fluid cells with solid cells after a deposition condition is reached. The study confirms the capability to simulate the deposition of oxygen controlled by mass transfer at low Reynolds number, but remarks that hydrogen deposition follows a different mechanism and needs a different approach. Specifically, in [59], the authors experimentally validate the model using an analogous setup for a precipitation experiment in water instead of sodium, comparing the increment in pressure after the deposition process, and in [60], a theoretical validation is presented for LR-LBM via Chapman-Enskog expansion recovering NS equations; here, the mesh resolution effect is also studied, finding good results regardless of the tested mesh sizes.
Gui et al. in 2019 [61] analyzed the characteristics of the flux in a pebble type High-Temperature Gas Reactor (HTGR). The LBM was integrated with the Boundary Immersed Method (BIM) and Discrete Element Method (DEM) to simulate the fluid-solid interaction between the gas and the pebbles: the proposed hybrid model was validated by pebble sedimentation experimental data. This model was able to simulate the recirculation process, an intermittent periodic flow pattern of pebbles, which implies a corresponding intermittent pattern in the gas circulation.

Neutronics: Transport and Diffusion
The second group of applications of LBM in nuclear engineering problems regards the neutronics; both neutrons transport (NT) and neutrons diffusion (ND). In this section, a very general recall of the simulation of the neutrons transport using the LBM is presented, and subsequently, the research items published regarding this topic.
The neutronics is the basis to understand the energy generation in the cores of nuclear reactors. A standard nuclear plant uses thermal neutrons and/or high-speed neutrons produced by the fission process to obtain new fissions in the fuel. In nuclear engineering, the neutron radiation/transport problem inside the reactor is a crucial factor, since the neutrons (and also the gamma radiations) emitted in the fission process travel inside different media from the fuel domain to the surroundings, composed by different materials, such as light-water or solid walls that can act as shields and/or as reflectors. The radiation can interact with these media through different reactions, such as scattering or absorption. The probabilities to have the different types of interactions are normally taken into account by the material properties, of which the most important one is the total macroscopic cross section Σ t = 1/λ, inversely proportional to the mean free path value λ, and including the absorption cross-section and scattering cross section.
The linear Boltzmann Neutron Transport Equation (NTE) substantially describes the processes by a conservation law. The first two terms on the left side of the next equation represent the effects of the advection on the neutron flux Φ(r, E, Ω, t), and the scattering and absorption are represented by Σ t (r,E), while the right side represents the neutron source S(r,E,Ω,t).
The NTE is a seven-dimensional equation (3 for space r, 1 for energy E, 2 for angular space Ω, and 1 for time t) for a group of neutrons with velocity v. A first way to solve this equation could be a stochastic approach such as a Monte Carlo (MC) based method: it can be applied using ad-hoc random functions; generally, single neutrons are simulated during all their time-life, and the process is repeated several times for each neutron to derive statistically the Neutron Flux Φ(r, E, Ω, t). Alternatively, deterministic methods or different statistical methods can be applied (by its statistical nature, LBM cannot be consider as a deterministic method or an ad-hoc random method). In this case, due to the difficulty in obtaining analytical solutions for complex geometries, some approximations are commonly used to simplify the partial differential equations before solving them by a numerical method:

I
The energy is discretized, normally using multigroup approximations. II The angular dimension is discretized, commonly using the Sn Discrete Ordinate Method (DOM), and/or Pn Spherical Harmonics Method (SHM).
The multigroup simplification discretizes the continuous energy space of the neutrons into a series of non-overlapping energy groups, assuming an average and constant energy for all of the neutrons belonging to each group, a classical (but rather simplified) approach uses only two groups: fast and thermal neutrons.
The Sn method discretizes the continuous angular variable into a series of non-overlapping solid angles, and subsequently replaces the continuous direction variables by a discrete set of direction vectors: this discretization implies the use of Gaussian quadrature, i.e., for 1D, is common to find Q2, Q4, Q8, and higher discretization quadrature up to Q20. The use of higher discretization quadrature can supply more accurate results as proved by the reviewed articles cited in this section. In this way, a multigroup-Sn-NTE is obtained, and the related partial differential equations can be solved by numerical methods, such as FEM or FVM, as well as the LBM, which also emerges as an option to discretize and solve the NTE, working like a DOM, but with a simpler angular discretization implementation. This makes the LBM more adaptable to different meshes (structured or arbitrary grid frameworks) and can provide a possible unique framework to neutronics and CFD calculations.
Methods such as DOM and SHM are similar in regards to their formulation, using linear variants of the Boltzmann equation to calculate the evolution of the density distribution of neutrons.
In the LBM based approach, collisions among neutrons are neglected, and neutron scattering is only related to the interaction between the neutrons and the nuclei of the propagation media. The discretized transport equation for the monochromatic neutron group g expresses the balance between collision and streaming of the neutron density distribution functions Ψ α,g (r, t) in each of the discretized angular directions, α, and is the basis for the basic Neutron Transport Lattice Boltzmann Method (NT-LBM); in this equation, S represents the source terms (truncation to first order is common and truncation to second order is called high order LBM or H-LBM scheme): To ensure the neutrons conservation, the neutron angular flux Φ g (r, t) must be equal to the summation of the local distribution functions Ψ α,g (r, t) over all the angular directions α; the same condition must be fulfilled by the equilibrium distribution function Ψ α,g eq (r, t). Some models consider the first right term of Equation (7) (BGK collision term) redundant for not-interacting particles; in this way, no equilibrium function is calculated, and the source term includes all the relevant effects, such as the scattering and absorption. In this review, this approach is called no BGK term.
In neutronic applications of the LBM, the common notation for the grid type DxQ M differs from the DxQy notations used in the CFD-LBM, because here M represents the number of possible angular directions in the scattering process, i.e., the lattice D1Q20 is spatially similar to the CFD lattice D1Q3, but here, Q20 represents the discretization scheme adopted for the angular space for scattering. In Figure 6, the simplest D1Q2 lattice, two different neutrons density distributions propagating to the left and right neighbor cells in each time step, is represented. If the first integral moment of the probability function Ψ α,g (r, t) is calculated, the scalar neutrons flux Φ g (r, t) in each grid point can be obtained from a simple summation over the α index for all the M angular directions as: Applying a Chapman-Enskog expansion, it is possible to recovery the NTE (6) from the discrete Neutron Transport Lattice Boltzmann Equation (7), relating the relaxation time τ N for the neutrons probability distribution with the transport process: the value of τ N is physically related to the probability of a collision of the neutrons with the background nuclei of the media, consequently, with the mean free path λ, and, consequently, with the macroscopic cross-section. In particular, the relationship for a group of neutrons with velocity v is: The boundary conditions in NT problems include volumetric sources or inlet flux conditions, reflective boundaries, vacuum boundaries, and bare boundaries, each one with an implementation on the LBM. Specific descriptions of these BC can be found in the related research items included in this review. There are a great number of benchmarks in one and two dimensions that are solved by validated numerical methods based on DOM or MC, and used by the researchers to test the LBM approach to neutronics; many of them are related to the calculation of the scalar neutron flux in different configurations with two or more regions, including an active zone that emits the neutrons and represents the radioactive nuclei in the nuclear reactor core. Of great importance is the calculation of the K eff eigenvalue, known as the criticality problem. It consists of finding eigenvalue solutions to the NTE because, at least in principle, a reactor might work in critical condition (i.e., with K eff = 1).
After early works in the 1990s, such as [62], the viability of using cellular automata to simulate neutron transport is, nowadays, established, and from this point, the LBM evolved taking profit of the advantages of the cellular automata models as an alternative to traditional stochastic models. However, in more recent work [63], a cellular automata model for neutron transport is applied to two benchmark problems, i.e., a neutrons pulse that propagates from the center of the computational domain (this particular problem is not easy to be solved by traditional harmonic expansion due to the difficulty in representing the initial single event using function expansions), and a second benchmark showing the evolution of a distributed source in a scattering medium, which puts in evidence the isotropic properties of the model.
The subsequent part of this section describes the main research papers regarding the use of the LBM in neutronic problems, which are summarized in Table 3.
Erasmus, in 2012 [64], confirmed the applicability of the LBM for neutron transport and proposed a modified collision scheme named "first collision source method" to overcome blind spots in the domain, which are drawbacks in the numerical method linked to the fact that the propagation of neutral particles in LBM (neglecting the scattering) are constrained to move on the grid, and some zones of the domain cannot be reached by a neutron propagating in a straight line. With the proposed scheme and the refinement in the angular and spatial discretization, it is possible to overcome this drawback, and LBM becomes as accurate as other numerical methods, although it can require huge computational resources. Another problem emerges in the calculation of the neutrons flux in a volume because the LBM calculates the flux in grid points and then is sensitive to the interpolating term in the flux integration over a volume (the effect of constant, linear, and exponential interpolations was calculated in the study). In the same way, Erasmus and Heerden in 2013 [65] present an angular refinement scheme to mitigate the ray effects (anomalous propagation in some diagonal grid directions) and compares their simulations with benchmark solutions using classical stochastic methods based on the Monte Carlo approach.
Bindra and Patil in 2012 [66] used the NT-LBM (no BGK term) to simulate NT problems in one-dimensional and two-dimensional arrays, searching a steady-state solution for the neutrons flux solving the Linear Boltzmann Equation (inter-neutrons collision term is neglected and the equation remains linear). The proposed scheme applies to anisotropic scattering terms, and the scattering integral is replaced by a summation on the angular Gaussian quadrature. They tested the NT-LBM with the Heaviside benchmark problems, with and without scattering terms. As an example of the implementation, they calculated a source problem that represents a nuclear system composed of a neutron source in a moderator-fuel matrix. The numerical simulations of this study showed that D2Q4 (LBM) is equivalent to S2 (DOM), but a direct comparison cannot be performed between higher-order schemes as D2Q8 and S4. Figure 7 shows a reproduction of the one-dimensional results based on our MATLAB implementation. The basis of [66] is the previous work of Ma et al., 2011, for one-dimensional radiation transfer [67].   In Bindra, 2013 [68], the paper [66] is extended using an NT-LBM (no BGK term) framework to solve NT problems. Here a one-dimensional problem, including two media, source and reflector, is simulated, and the criticality constant (eigenvalue problem) for a simple benchmark was found validating the application of the NT-LBM in criticality problems. Further benchmark validations are suggested to develop a multiphysics LBM framework for nuclear engineering applications.
Starting from the Boltzmann equation, McCulloch in 2013 [69] analyzed some possible solutions for heat radiation and neutron transport problems implementing the weighted summations for the scattering term proposed in the previous models. For heat radiation, the radiative heat transfer was calculated, including conduction and convection heat transfer mechanisms. The solutions for the radiation problem in 1D, and the conductive-convective-radiative heat transfer in 2D were presented, showing the NT-LBM (no BGK term) capability to solve multi-physics problems in a single framework.
For NT, the steady-state criticality problems in one-dimensional geometries were calculated, with a multigroup approach (fast and thermal neutrons groups considered). Different slab geometries were tested and a good agreement with the benchmark values for fuel domains bigger than 3 mean free paths (without reflectors) was found. One relevant issue highlighted in the LBM is the simple parallelization scheme and the easy code-development needed to implement the neutron transport or the thermal field. Agarwal et al. in 2015 [70] used the NT-LBM (no BGK term) for criticality problems in one-dimensional geometry (Isaa cell benchmark) and 2D arrays (BWR matrix benchmark with fuel and moderator) composed of two materials and using a multigroup approach (fast and thermal neutrons). They compared the performance of D1Q2 vs. D1Q4 and D2Q4 vs. D2Q8 angular quadrature and found that the LBM has a good degree of accuracy for solving eigenvalue problems in NT and a high quadrature model improves the accuracy in the determination of the K eff eigenvalue; moreover, the convergence analysis concerning the effects of the grid resolution on the eigenvalue shows a quasi-linear trend. Figure 8 shows a reproduction of one-dimensional benchmark problems used for criticality problems (obtained by an in-house developed MATLAB routine). Here, the power iteration method was used to determine the value of the constant K eff . This iterative method starts guessing the source term and the value of K eff , and then these two values are updated using the neutrons flux calculated in each time step. The process is repeated until a convergence criterion is reached (further details on the implementation can be found in [69]). Gairola et al. in 2016 [71] implemented the NT-LBM (no BGK term) framework to solve time-dependent non-equilibrium neutron transport and criticality problems for a reactor through a model that assumes isotropic scattering. Two grids (D2Q8 and D2Q16) were proposed for the time-dependent problem to solve the neutrons flux in a heterogeneous medium composed of several materials. Two benchmarks were solved: a 2D nuclear checkboard with an array of pure scattering and absorption regions around a source and a BWR primitive cell with two different regions. The used model allows the simulation of source, diffusion zones, and absorption zones. The results calculated using the two grids were compared and agreed with DOM reference values: from these results, it is evident that from both lattices quadrature D2Q8 and D2Q16 coherent solutions in two-dimensional NT problems can be obtained, but these two lattices have drawbacks due to the spurious calculations of some ray effects.
Using the NT-LBM (no BGK term) framework in a time-dependent one-dimensional problem, and then in a 2D NT checkboard benchmark problem with a monoenergetic neutron source at the center of the domain, and pure scattering and pure absorptive array of blocks, Gairola and Bindra, in 2017 [72], compared the results from the D2Q8 LBM with the DOM S6 solutions: a good agreement near the source was found, but numerical ray effects appear in both methods far from the source (near the boundaries of the domain). These effects do not match between them because of the difference in the discretization schemes. In Figure 9, the bi-dimensional checkboard benchmark is represented (obtained by an in-house developed MATLAB routine) using D2Q8 and D2Q4 lattices, showing the diffusive patterns and the numerical ray effects for the steady-state solution.  [68]; (b) bare nuclear fuel as presented in [68] (without scattering) and [69] (isotropic scattering); (c) color map for Isaa cell benchmark [70]; (d) scalar neutron flux for Isaa cell benchmark [70].  [74] developed the finite LBM scheme for transient problems of radiation transport (RT) and NT that handles multidimensional scattering terms in transient problems: the continuous problem is spatially and angular discretized independently. While the works previously discussed presented a low angular resolution, the method proposed in this work is an integration of the Sn technique into LBM. The angle discretization and the spatial discretization are explicitly unlinked in this model. Each angular direction is independent, following a Lattice Boltzmann equation for each spatial direction, i, and each angle, α. In this way, the distribution functions must be written using an additional subindex, i.e., Ψ α,i,g (r, t) instead of Ψ α,g (r, t) (used in the previous models). This is the clue for having local dynamic rules for a high angular resolution and being highly parallelizable.
A hybrid framework FV-LBM was proposed by Wang et al. in 2017 [75] to simulate NT and RT. This method showed flexibility in the geometric representation of the domains using unstructured meshes (tetrahedral or triangular) with good computational efficiency. The physical meaning of the collision term in the Lattice Boltzmann equation was highlighted by the authors. This term represents the collision between the neutrons and the background nuclei instead of the collisions among the neutrons themselves. The tested benchmarks were a one-dimensional array composed of two material domains, a two-dimensional square array with two materials, and a two-dimensional representation of a spherical source-shell array. These simulations showed the possibility of using the hybrid scheme in transient and steady-state neutron flux problems with complex curved geometries.
In the work of Ma et al., (2017) [76], the optimization of the NT-LBM is proposed, based on a refinement technique called Structured Adaptive Mesh Refinement (SAMR). This technique allows for making a progressive mesh (grid) refinement, subdividing the domain into blocks, and applying a successive grid refinement in the different blocks until some convergence criterion is reached, allowing, in this way, the improvement of the resolution in a specific space region. The proposal overcomes the classic Adaptive Mesh Refinement drawbacks, such as data storage problems, and subsequently reduces computational demand, allowing accurate results in the simulations of NT problems. This study describes a technique to optimize the inter-block communication and overcome the jumps in the fine-coarse grid. The SAMR model was tested with two-dimensional benchmarks at steady-state and time-dependent conditions using the multigroup approach. Wang et al., in 2018 [77] improved the grid refinement using a technique called Streaming-based block Structured Adaptive Mesh Refinement (SSAMR) to overcome discontinuities in the flux calculation, using the streaming process in the LBM to optimize the communication between the blocks. Moreover, this study presented the theoretical deduction of the Neutrons Diffusion (ND) equation from the NT-LBM to show the physical meaning of the scheme. The dimensionless relaxation time describes the characteristic time needed to evolve from a non-equilibrium local distribution to a local equilibrium distribution, and depends on the collisions between the neutrons and the background nuclei in the media (mean free path). Specifically, through the Chapman-Enskog expansion, the neutrons diffusion equation was derived as a first-order approximation of the NT-LBM, then, it was theoretically showed that NT-LBM is more accurate than the ND equation to model neutronic problems. The model was tested with source-driven three-dimensional and two-dimensional benchmarks and in criticality problems.
Ma et al., in 2019 [78] developed a Multi-Block Adaptive Mesh Refinement (MB-AMR) scheme for NT in the LBM framework, using the finite Boltzmann equation and a multigroup Sn approach for the discretization problem. The multi-block approach divides the domain into fixed subdomains and the AMR dynamically defines finer or coarser grids to solve neutrons fluxes. Since the data structure is complicated, the communication between blocks is not a simple task; the combination MB-AMR simplifies this issue, using predefined blocks. The proposed method was tested with several benchmarks in one-and two-dimensions, such as a multilayer shield problem, Reed's cell short pulse transient, a square homogenous source-driven problem, and a mixed-oxide (MOX)-type reactor, i.e., a complex geometric configuration with UO 2 and mixed-oxide (MOX) fuels surrounded by light water reflectors.
Another hybrid model (neutron finite volume) NT-FV-LBM to solve the Neutron Discrete Velocity Boltzmann Equation was presented by Wang et al. in 2019 [79]. The effect of the discretization of the angular direction was studied, and it was found that it is possible to converge to a solution independent from the velocity discretization scheme. The need to include the polar angle in the scattering, i.e., implementing a D3QM scheme, even if the simulation represents a bidimensional domain, was highlighted. The neutron diffusion solution converges to NT-FV-LBM in relatively low scattering media, but is more accurate for higher scattering media.
In 2019, Wang et al. [80] developed the High-order Lattice Boltzmann Method (HLBM) for multigroup neutrons diffusion. The idea was to derive a method with a higher-order expansion from the Chapman-Enskog method, which could lead to more accurate calculations of the source term. This ND-HLBM solves the multigroup neutron diffusion in the transient regime and criticality problems. Comparing HLBM vs. LBM in 2D and 3D source driving and eigenvalue benchmarks shows how the modified model drives to more accurate simulations.
It is interesting to see how new technologies allow simulating problems with complex geometries thanks to the increasing computational power and the more efficient schemes implemented for the numerical methods. The use of GPU hardware architectures allows the LBM to improve its capabilities, thanks to the intrinsic parallelism of the method. In Wang et al., (2019) [81], the implementation of the LBM (GPU-ND-LBM) to solve stationary and transient problems in neutrons diffusion (using single and multiple GPUs) was presented. The code used a Predictor-Corrector Quasi-Static Method to compute the accelerated neutrons diffusion. This method factorizes the neutrons flux in two parts: shape (quasi-stationary) and amplitude (time-dependent only). In this way, the computational time is reduced. The validation was done using benchmarks for steady-state eigenvalue problem, and one space-time kinetic problem (a step perturbation transient system). The method can be implemented in a multi-group approximation using, e.g., two neutrons groups: fast and thermal. It was highlighted that "for the LBM approach, the non-linear term is local, and the non-local term is linear"; this implies an advantage for parallel computation. The diffusion equation is a low-order isotropic approximation of the neutron transport equation, and for this reason, Wang et al., in 2019 [82], improved the previous work with a GPU Neutron Transport LBM (GPU-NT-LBM). The method was tested with a mono-dimensional transient benchmark, a two-dimensional benchmark of a typical PWR, and a two-dimensional iron-water benchmark, making a comparison with the DOM solutions.
Agarwal et al., in 2020 [83] implemented a multigroup NT-LBM (no BGK term) for criticality problems in optically thin regions (a mean free path close to the system's dimension) and highly anisotropic scattering systems. In this paper, one-dimensional benchmarks with a single group and multigroup approaches were simulated. The proposed scheme was tested with isotropic and anisotropic scattering media using a high number of angular directions in slab geometries. Findings confirmed that, for optically thin slab regions, when mesh spacing decreases and angular discretization increases, the calculation of the multiplication factor improves, overcoming the drawbacks that can emerge when the leakage from the boundary is not captured. In conclusion, this work improved the NT-LBM and validated the use of the method to solve criticality problems with multigroup and non-homogenous scattering, as well as for optically thin slabs.
Wang et al. in 2020 [84] presented an analysis of the use of LBM in multiphysics simulations of nuclear reactors, coupling neutronics, and heat transfer physics in a single model. In particular, the proposed framework uses the same evolution equation and same mesh to solve the different fields (neutrons flux field and temperature field) and only needs to change some variable definitions for each field. In this way, the computational scheme remains simple and the interaction between the physical fields can be implemented directly. The proposed Neutron Transport Heat Transfer Lattice Boltzmann Method NT-HT-LBM was validated in transient and steady-state simulations (compared with hybrid FVM-LBM results) and shows a simple unified framework that can be used in large-scale parallel computations of the nuclear reactor physics with the integration of classical neutronic approximations as LBM-Sn or LBM-Pn.
The wide field of applications of the LBM in nuclear engineering problems includes: • Mixing of the coolant in the vessel of the reactor, which helps to understand if a thermocouple in the top head of the vessel can measure the mean temperature of the coolant at the exit, knowing the mixing conditions. The isothermal flux around fuel rods is presented in [37][38][39]44,[49][50][51][52]54,58]. A thermal approach to flow is explored in [53,55,56]. The investigation of the flux in the turbulent regime around the fuel rods is of great interest, this includes geometric details, such as the grid separators.

•
Possible fractures in the reactor after a Small Break-Loss Of Coolant Accident (SB-LOCA), by a high-temperature gradient following the emergency cold injection by a high-pressure injection system [42,43]. • How, after a LOCA, the debris (contaminant particulate material) can be trapped in weir traps under a wire floor [46,47].

•
How flow acoustic resonances in steam dyer pipes are generated and can cause problems related to fatigue [48].
Some works in the CFD field related with nuclear engineering do not need a relevant modification of the traditional LBM scheme [46,47,49]; however, to enhance the applications in CFD problems, it is common to adapt the LBM with high discretization schemes in 3D, such as D3Q27, and with turbulence models, LES being the most basic turbulence model [37][38][39]48,51,53]. More sophisticated turbulence models, such as VLES or LES-WALE [44,52,54,55,58], are also used, leading to better results than LES, and that of the use of the LBM (only without taking into account turbulence). The inclusion of thermal field models is relevant [42,43,53,55,56] to upgrade the classical isothermal LBM to a numerical method capable of simulating the heat transfer mechanisms present in the reactor (in particular, the convective heat transfer). The most commonly used is based on the secondary distribution approach. The inclusion of species transport is used to simulate chemical reactions, such as corrosion and deposition [40,41,45,59,60].
The LBM treats the complicated neutrons transport phenomena as simple linear calculations, and the use of angular discretization schemes leads to convergence. In early works, the solutions were obtained for isotropic scattering medium and with domain dimensions significantly larger than the mean free path. The multigroup solution was obtained and tested with some 1D, 2D, and 3D simulations, for both simple and complex geometries. The use of primitive cells is a commonly adopted strategy in both CFD and neutronic problems. It is possible to enhance the efficiency of the model for complex geometries using FV-LBM or mesh refinement schemes; thus, reducing the computational cost.
In all of the research items related to neutronic problems, the LBM is adapted to solve NT (or ND) equations with the inclusion of an angular discretization scheme DxQm that considers scattering effects. Hybrid schemes using the FV approach to handle different meshes and represent complex domains are also found [75,79]. Mesh refinement strategies are adopted to optimize resolution, handle curved geometries, and to reduce computational cost [76][77][78].
It is possible, at first sight, to consider the LBM as a sort of Sn-like method, where angular directions are tied to the spatial mesh, and to consider that the comparative advantages for parallelization of the LBM over Sn could be lost under anisotropic flux due to the use of distributions with high angular quadratures that could induce coupling between not-neighboring lattice points. However, the finite LBM scheme or Sn-LBM (introduced in [74,84]) decouples the angular and spatial discretization (avoiding the coupling between not-neighboring lattice points) and shows that LBM is a flexible solver for the Sn approximation that can hold the parallelization properties by applying the spatial-angular parallel mode proposed in [78], which splits the angular discretization as an additional dimension on the computer memory array.
It must be noted that, in both CFD and neutronic applications, the validation of the numerical LBM is usually done by comparison with benchmark solutions obtained by other well-established numerical methods, such as FVM in CFD or DOM and SHM in neutronics. In a few cases, experimental data are compared with the simulations; the validation of the LBM comes from analytic or experimental data of simple benchmarks (as the velocity profile in a pipe) and not necessarily from an experimental setup related with the simulated NE problem.
Recent research not only suggests the possibility to use the LBM in the nuclear engineering field, by the direct comparison of benchmark results obtained by other numerical methods, but also by the recovering of the macroscopic equations, such the Navier-Stokes equations (including the energy equation), the neutron transport equation and the neutron diffusion equation from the basic Lattice Boltzmann equation, thanks to a strong theoretical basis derived mostly from the Chapman-Enskog expansion. The LBM can, thus, be considered more than a computational automata, involving some physical principles that emerge from the statistical mechanics.
It is noted that each evolutionary step in the implementation of the LBM drove (and will drive) to validated engineering tools with flexible features as a single framework to solve fluid dynamics, convective, conductive, and radiative heat transfer, neutron transport, and more complex microphysics phenomena, such as chemical reactions and species transport. This shows that the LBM has an intrinsic capability to integrate powerful multi-physics simulation frameworks. This capability is enhanced with the integration of hybrid models, such as the Finite Volume method to adjust unstructured meshes [75,79], or the Boundary Immersed Method to calculate fluid-solid interactions [61].
To date, a complete multi-physics simulation of a nuclear reactor, including NT, heat transfer (with the integration of radiative, convective, and conductive mechanisms), CFD and transport-deposition mechanisms for a secondary chemical species has not yet been done, but the LBM is a candidate for a possible comprehensive framework to study the nuclear reactor as a whole. A great advance is presented in [84], coupling neutron transport and heat transfer in a single framework. Using a single framework, the LBM could overcome, partially, the limitation related to the use of several interfaces "Multi-physics assessment is subject to the use of several methods and, just as importantly, several interfaces between physical domains. For this reason, only limited credit may be claimed in safety analysis" [21].
Some commercial software implementing LBM are available, as well as some open-source libraries (see Appendix A for further details), but up to now, the general trend of research is to use and test self-developed codes, proving relatively simple implementation, but also the customization possibilities of this numerical method. To access all of the potential of the method, it is worth mentioning the possibility of implementing the LBM in parallel computing platforms (GPUs and clusters [43,81,82]) that not only enhance the spatial and temporal resolution of the results, but also optimize the computational cost with the inclusion of advanced mesh refinement schemes [76][77][78].
The application of the LBM in NE problems could have many more applications, such as the inclusion of multi-phase (liquid-gas) CFD simulations of BWR, or of the jet breakup (relevant after a core disruptive accident of a sodium-cooled fast reactor) [85,86], or some other topics (such as in [87]) concerning the development of a hybrid framework that couples microscale models and the LBM to simulate flows around reactive blocks (both chemical and nuclear) or a reactive flow: here, two different hybrid models were compared, i.e., the mean-field deterministic equation and the reaction kinetics coupled to LBM for the transport process, and the microscale kMC (kinetic Monte Carlo) coupled to LBM. "Apparently, the coupled kMC-LBM framework can fulfill the requirements of a unified multi-scale reactive flow simulator" [87].
Another possible application of the LBM regards the simulation of nanofluids to enhance the heat transfer properties (e.g., [88][89][90][91]), with the great advantage being that the simulation uses a distribution function for the fluid, a second one for the nanoparticles, and a third one for the thermal field. All three lattices are overlapped, with simple rules interacting between them, making it unnecessary to define complex geometrical conditions or interfaces between the nanoparticles and the main fluid.
Moreover, there are some characteristics of the Lattice Boltzmann Method that justify the use of the LBM over other numerical techniques, or even work together in a hybrid scheme:

1.
The algorithmic implementation is simple.

2.
The meshing technique is time-efficient, generating a structured lattice with relatively easy refinement schemes. 3.
LBM is highly parallelizable; the collision term is local and then the calculations can be easily distributed between different cores or GPUs. The communication for the streaming process only uses neighbor nodes, and this leads to a highly parallelizable interface.

4.
The multiphysics coupling is reduced to a single framework. The communication between a series of overlapping lattices, fluid lattice, thermal lattice, and neutronic lattice is done by a simple local scheme before the collision step, then the parallelization property is conserved.

5.
Directly linked with the previous point is the possible inclusion of a multiphase flow using an additional distribution function, also applicable for nanofluids. 6.
The thermal-hydraulic and neutronic results are validated by comparison with benchmarks. 7.
The results are accurate for the standard methods. 8.
The scheme is highly flexible and adaptable with hybrid schemes, using standard methods and turbulence models. 9.
The use of the common framework for neutronics proposed in [84] enables the use of the LBM as a solver for the Sn, Sp3, or P1 approximations, maintaining the parallelization properties of the method. 10. The use of Sn-LBM, the spatial-angular parallelization mode, and the multi-GPU technique can speed up the NT-LBM, solving the neutron transport problem accurately and effectively without affecting the convergence characteristics [82].

Conclusions
In this paper, a systematic literature review concerning the use of LBM in nuclear engineering problems was presented, considering data from 45  In CFD problems, it is possible to adapt the LBM with turbulence models, high discretization schemes in three-dimensional lattices, thermal models, and mesh refinement strategies.
In neutronic problems, the LBM is adapted to solve the NT and/or ND equations with the inclusion of an angular discretization scheme, such as the DOM technique, which takes into account scattering effects. A common framework to implement Sn-LBM, Sp3-LBM, and P1-LBM models was developed.
The validation of the LBM results is commonly made by the comparison with experimental measurements of derived macroscopic variables (as friction factor, or velocity field for CFD) or with a direct comparison of the solutions of well-established benchmarks with other numerical methods, such as FVM (for CFD) or DOM (for neutronics).
The interest in applying the LBM to NE problems is increasing, and the related research items show this tendency, with an increase in the number of publications in the last years. Furthermore, the perspective to develop a multi-physics platform based on the LBM is foreseeable in the literature, and may include: • Two-dimensional (2D) and 3D high discretization schemes. LBM is rapidly approaching to a mature development, and the potential of this numerical method is a strong reason to conduct a better-directed effort to validate its use as a highly accurate multiphysics simulation tool for the nuclear engineering field. It is fundamental to point out that theoretical research in LBM is going on, and is aiming (and succeeding) at resolving issues related to stability, speed, and model expansions.
Finally, it is worthy to remember that to access all of the potential of this numerical method, it is possible to exploit the intrinsic parallelization features of the LBM, using, for example, computing platforms, such as GPUs and clusters. This will not only enhance the spatial and temporal resolution of the results, but will also improve the computational effort. It is important to consider that the flexibility of the LBM opens the opportunity to develop hybrid methods to study complex processes, including microphysics (as deposition, oxidation, and chemical reactions).

Appendix A Overview of the Lattice Boltzmann Method
In this section, a basic description of the LBM is presented. It is not meant to be a rigorous deduction of the method, but just to serve as an introduction to it. A more detailed overview of the method can be found in [93,94]. If interested, the reader can find parallel computing optimized LBM open-source libraries, such as Palabos [29,30] and OpenLB [95,96].
A novel way to solve the mechanics of a fluid is with the Lattice Gas Automata (LGA) method, which solves the fluid dynamics by applying a local interaction rule between a discrete number of molecules that move with a discrete set of velocities in a grid (or lattice). In 1974, a first collision rule for the LGA on a square grid was published by Hardy, Pomeau, and de Pazzis [97].
After that, in 1986, a new isotropic rule for a hexagonal grid was published by Frisch, Hasslacher, and Pomeau [98], including the recovery of the Navier-Stokes equations. These two models are considered the base of the LGA method. Lattice Gas Automata has some drawbacks, such as statistical molecular noise, the non-physical meaning of some relations (as pressure dependence to velocity), and lack of Galilean invariance. The LBM overcomes the LGA drawbacks, taking also a restricted set of velocities for the molecules in a grid, but using a linearized form of the collision operator known as Bhatnagar-Gross-Krook (BGK) 1954 [99]), and using a probabilistic distribution of molecules instead of single molecules in the grid [100].
The first paper regarding the LBM is from 1988, by McNamara and Zanetti [101], and the first viable computational implementation is from 1989 by Higuera et al. [102]. The idea is that macroscopic fluid properties emerge from microscopic interactions of many particles, but the details of the interaction of a single molecule do not affect the statistical macroscopic behavior of the fluid. The statistical distribution of particles (or density distribution) is represented by f (r, c, t); for a given time, t, the probability to find a particle in a position, r, with a velocity, c, is the quantity f (r, c, t) that evolves following the classic continuum Boltzmann equation, representing an equilibrium between the transport and collision process.
To solve the Boltzmann equation, in addition to the BGK operator, the discretization of the velocity space and the position space in the grid DxQy is used, where x is the number of dimensions D (1, 2 or 3) and y is the number of possible velocities or velocity quadrature Q (see Figure 1). In each position of the grid, some y degrees of freedom allow the streaming and the interaction between the density functions. Each one of the y degrees of freedom has a discrete velocity vector c i and a distribution function f i (r, t) is directly associated. In this way, the collision and transport between neighbor distribution functions are calculated individually with the following discrete lattice Boltzmann equation for a given time step ∆t: In the LBM algorithm, this discretized Boltzmann equation is divided into two steps, called stream-step and collision-step. In the collision-step, the first term in the right side of the equation represents the BGK collision operator and is used to calculate from the incoming population f i (r, t) an outcoming f i (r, t) population relaxed to the equilibrium f eq i (r, t) (τ represents the Single Relaxation Time (SRT) of this process). In the streaming-step, the left side of the equation is used for a single step of time ∆t to transport the value of a particular density function f i (r, t) to the next grid position f i (r + c i ∆t, t + ∆t). The external source term S i (r, t) describes the coupling of the fluid with the external environment, and is the key to include external forces or potential energy interactions.
The equilibrium distribution function f eq i (r, t) is approximated in the LBM as a Taylor expansion of the Maxwell-Boltzmann distribution. This equilibrium distribution can be calculated from the following equation that only depends on the macroscopic quantities density ρ and velocity u of the fluid: The sound speed c s is a propagation constant in the computational domain that depends on the particular quadrature model. In the simplest case of D2Q9 c s = 1/ √ 3, the discrete set of velocities c i and the correlated weights w i are also fixed in each quadrature model.
The macroscopic physics quantities can be calculated from the populations f i (r, t) in each grid point, using the different moment integrals of f i (r, t) (integration in the quadrature space). The integration process leads to macroscopic variables as density ρ(r, t), velocity u(r, t), momentum flux, or stresses. The calculation is relatively simple and only requires a simple summation over the y degrees of freedom in each position, i.e., for local density ρ(r, t) the first moment integral: and for the fluid velocity the second moment integral: The pressure P(r, t) can be obtained by the simple state equation: The Chapman-Enskog expansion is an analytical method that allows the recovering of the NS equations from the Lattice Boltzmann Equation (A1); [103] presents a comparison of the macroscopic equations recovered by the Chapman-Enskog expansion or the Taylor expansion. Following the Chapman-Enskog expansion, it is possible to obtain the value of the relaxation time τ and its dependence on some physical properties. In the described BGK-SRT model, the relaxation time is related to the kinematic viscosity γ of the fluid by: The basic LBM-BGK(SRT) approximation is valid for relative low Mach numbers (Ma < 0.3) to be within the incompressible limit, but the model has evolved to be applicable to compressible flow and enhance the Mach number limit, some relevant modifications are called the regularized model or the entropic compressible LBM, a review of the progress in the development of the LBM for a high Mach number can be found in [104]. The Boltzmann equation has no limit for Kn and, for this reason, the LBM can be applied to study microfluidics problems. Regarding the Re number, LBM can handle high values in the turbulent regime and can converge to a solution if a high dimensional discretization scheme is applied, and some turbulence model is adopted.
To enhance the stability of the BGK model, d'Humières et al. [105] introduced a Multiple Relaxation Time (MRT) approach, writing the collision operator, using the moments space instead of the velocity space. This drives to multiple equilibrium equations and multiple relaxation times, allowing the model to remain stable at high Reynolds numbers but increasing the computational time. The entropic model also enhances the stability of the LBM. Modified relaxation time is calculated after each time step, and if the solution began to diverge from a stable solution, the modified relaxation time should control the divergence. In order to know which modified relaxation time will correct the instability, an entropy function is determined. The first entropic model was proposed in 1998 by Karlin and Succi [106], and in 1999 by Karlin et al. [107] using a discretized form of the continuum Boltzmann entropy function.
In 2001, Boghosian et al. [108] used instead the Tsallis entropy function. Entropic LBM is stable at very high velocities without the computational drawback of the MRT model.
A third attempt to give more stability to the LBM was developed in 2006 by Latt and Chopard [109], i.e., the regularized LBM. In this model, the density distribution function is separated into an equilibrium part and a non-equilibrium part. The collision operator was modified to operate on the non-equilibrium part. The enhancement of the stability is lower than in the entropic method, but has no impact on the computational time, and can be combined with the MRT. A common framework and a comprehensive analysis of the different collision models can be found in [110].
Implementation of the Boundary Conditions in LBM is an open field of research. Some early works, as the one by Zou and He in 1997 [111], expand the bounce back boundary conditions to include conservation of mass and momentum. In 2002, Zhao et al. [112] proposed an improvement of the BC based on extrapolation schemes. In 2003, Yu et al. [113] proposed a treatment for the curved boundaries based on a single interpolation scheme. An overview of the boundary conditions for the LBM, with special attention to heat and mass transfer can be found in [114].
Some other interesting developments are related to the classic 1993 Shan and Chen [115] model for multiple fluids, based on adding an interaction force to the equilibrium function.
The current status of the Lattice Boltzmann Method applied to aerodynamics, aeroacoustics, and thermal flows, was presented in [116].
In this additional section, some basic concepts of the Lattice Boltzmann Method are presented. The collected literature gives an idea of the continuous development of the method in the last decades (it is not an exhaustive timeline, but some key papers are cited). The amount of literature concerning LBM is increasing; the theoretical development does not stop, and future perspectives of the method can be found in [117].