Lattice Boltzmann method for simulation of time-dependent neutral particle transport

NUCLEAR ENERGY SCIENCE AND ENGINEERING

Lattice Boltzmann method for simulation of time-dependent neutral particle transport

Ya-Hui Wang，

Li-Ming Yan，

Bang-Yang Xia，

Yu Ma

Nuclear Science and Techniques

Vol.28, No.3

Article number 36

Published in print 01 Mar 2017

Available online 03 Feb 2017

DOI：10.1007/s41365-017-0185-z

35000

In this paper, a novel model is proposed to investigate the neutron transport in scattering and absorbing medium. This solution to the linear Boltzmann equation is expanded from the idea of lattice Boltzmann method (LBM) with the collision and streaming process. The theoretical derivation of lattice Boltzmann model for transient neutron transport problem is proposed for the first time. The fully implicit backward difference scheme is used to ensure the numerical stability, and relaxation time and equilibrium particle distribution function are obtained. To validate the new lattice Boltzmann model, the LBM formulation is tested for a homogenous media with different sources, and both transient and steady-state LBM results get a good agreement with the benchmark solutions.

Transient neutron transportLattice Boltzmann methodLinear Boltzmann equation

1. Introduction

The neutral particle transport [1] has attracted great attention in different fields of science [2], such as neutron capture therapy (mainly in boron neutron capture therapy) [3], transport of charged carriers in semiconductor devices [4], nuclear reactor design [5], radiation protection, radiotherapy, astrophysics, traffic flow, medical imaging, shielding design, oil well- logging and nuclear fuel transport [6].

It’s well known that the particle density can be obtained by solving the Boltzmann transport equation (BTE) (more popularly known as transport equation) [6]. Analogously, the distribution of neutral particles for a determinant system and equipment can be express with the neutron transport equation (NTE) [7]. It has 7 independent variables, 3 for space, 1 for time, 2 for transfer direction and 1 for energy. However, it is a tough task to solve the BTE, and researchers have changed their attentions to numerical solution and computational techniques [8].

Several mathematical and numerical techniques have been taken to solve the NTE [9]. For complex models, the Monte Carlo (MC) method is widely used to simulate particle transport process [10-15]. The deterministic finite element analysis has been applied to solve BTE for time-dependent and arbitrary geometry problem [16], and the finite element method (FEM) has been used in nuclear power reactor simulation [17]. The method of characteristics (MOC), a numerical integration method for partial differential equations, has been widely used for nuclear reactor physics lattice calculation [18, 19]. As the expansion of MOC, the OpenMOC code has been developed to research the speedup computation technique and parallel algorithms of MOC [20]. To improve the MOC efficiency, an acceleration method named CMFD is developed. Its implementations in 3D MOC simulation show that it can reduce the iteration number and computation time [21]. The finite difference method (FDM) emerged as the simplest technique to solve the BTE, but this requires fully implicit formulation for desirable numerical stability [22].

The discrete-ordinate method (S_N) and spherical harmonics method (P_N) are widely used to discretize the BTE in angle because of its simple construction, but their major disadvantage is the unphysical oscillation [23]. Thus efforts have been made in developing new methods, which include the cellar neural networks (CNN) coupled with P_N method to solve the time-independent NTE in x-y geometry [24]; the Random Variable Transformation technique (RVT) to solve the stochastic transport problem with the continuous stochastic media [25]; the GTNEUT code, based on transmission and escape probabilities method, to replace the classical Monte Carlo methods [26]; the PARAFISH to solve the second-order even-parity scheme of the time-independent BTE with an algebraic domain-decomposition method [27]; and the quasi-diffusion (QD) method to solve the neutral particle transport problems in Cartesian XY geometry with unstructured quadrilateral meshes [1]. Similarly, the ADAFENT (adaptive finite elements for neutron transport) code was developed to solve the second order even-parity BTE based on the principle of K+ variation [6]. A variable order spherical harmonics method was used to solve the BTE [28]. Two new wavelets on basis of Sweden’s second generation wavelets [29] were applied to discretize the angular dimension of the first-order BTE [23]. Then, a second generation spherical wavelet method [30], was improved into the 3-D wavelet-based method to solve the first-order BTE [31]. A novel hybrid technique for solving the time-dependent NTE was proposed [32], and an expansion of the hybrid transport point kinetic model appeared for initially-critical multiplying structures [33].

Lattice Boltzmann method (LBM), developed for computational fluid dynamics (CFD) originally, is used increasingly in the field of heat and mass transfer [34-36]. It assumes that the macroscopic appearance of a liquid is the result of collective behavior of many microscopic particles in flow field and it is not susceptible to individual particles. Numerical calculation of LBM is popular, with advantages of clear physical model and background, easy treatment of boundary condition and convenient parallel algorithms (because of the locality of evolution equation) [37]. LBM has been used in nuclear reactor engineering, such as large eddy simulation of a sub-channel of a rod bundle [38], and a two-phase flow simulation [39].

In this paper, we propose a novel lattice Boltzmann model to solve the neutron transport problem. After the detailed derivation of the LBM model, the neutron flux densities in one-dimensional homogenous media with different sources are calculated. The time-dependent and steady-state LBM results agree well with benchmarks and existing results.

2. 1-D neutron transport lattice Boltzmann model

The time-dependent linear Boltzmann equation for mono energetic neutral particles with an isotropic source in homogeneous medium can be described as [9]

(\frac{1}{v} \frac{\partial}{\partial t} + \hat{Ω} \cdot \vec{\nabla} + \sum_{t}) ψ (\vec{r}, \hat{Ω}, t) = \frac{Q (\vec{r}, t)}{4 π} .

(1)

where $\hat{Ω}$ is the angular direction, and $Q (\vec{r}, t)$ is the total source term in position $\vec{r}$ and time t.

For one-dimensional problem, Eq. (1) can be rewritten as [2]

\begin{array}{l} [\frac{1}{v} \frac{\partial}{\partial t} + μ \frac{\partial}{\partial x} + \sum_{t} (x)] ψ (x, μ, t) \\ = \frac{\sum_{s} (x)}{2} \int_{- 1}^{+ 1} f (μ' \to μ) ψ (x, μ', t) d μ' + q (x, μ, t) . \end{array}

(2)

where, v is the velocity of neutrons emitted from the medium; μ is the direction cosine of the polar angle (−1≤ μ ≤1); ψ(x, μ, t) is the transient neutron flux corresponding to position x and direction μ; and Σ_t(x) = Σ_a(x)+Σ_s(x) is the macroscopic cross section, with Σ_a(x) being macroscopic absorption cross section and Σ_s(x) being macroscopic scattering cross section; $f (μ' \to μ)$ is the anisotropic scattering phase function; and q is the external particle source.

To homogeneous medium, the source in each direction is isotropic, and the neutron transport can be described as

[\frac{1}{v} \frac{\partial}{\partial t} + μ \frac{\partial}{\partial x} + \sum_{t} (x)] ψ (x, μ, t) = S_{t} (x, t),

(3)

where S_t is the total source term given as

S_{t} (x, t) = \frac{Σ_{s} (x)}{2} \int_{- 1}^{1} f (μ' \to μ) ψ (x, μ', t) d μ' + q (x, t) .

(4)

The anisotropic scattering can be express by a linear anisotropic phase function $f (μ' \to μ) = 1 + a μ μ'$ , where a (−1≤ a ≤1) denoting the anisotropy scattering coefficient. Eq. (4) of the source term can be written as

S_{t} (x, t) = \frac{\sum_{s} (x)}{2} \int_{- 1}^{+ 1} (1 + a μ μ') ψ (x, μ', t) d μ' + q (x, t) .

(5)

The anisotropy scattering coefficient in isotropic medium equals to zero (i.e. a=0), while a>0 for forward scattering and a<0 for backward scattering. In this paper, we just study the isotropic condition (a=0). To ensure the numerical stability, fully implicit backward difference scheme in time is used [40]. Then, Eq. (3) can be written as

\frac{ψ (x, μ, t) - ψ (x, μ, t - Δ t)}{v Δ t} + μ \frac{\partial ψ (x, μ, t)}{\partial x} + \sum_{t} (x) ψ (x, μ, t) = S_{t} (x, t),

(6)

where Ψ(x, μ, t−Δt) is the neutron flux of the last time step and Δt is the time step.

According to the common terminology referred in LBM, D1QM lattice structure (Fig.1) is used, where 1 represents one dimensional geometric space, and M refers to the number of lattice velocity, i.e. the number of linkages of a lattice node. We note that this D1QM model is not the standard D1Qq model in terms of conventional LBM, in which the base vectors in different lengths are parallel to a dedicated direction [41]. Eq. (1) depends on not only time and spatial coordinates, but also solid angle. This differs greatly from flow governing equations. In order to describe the direction characteristic of the neutron transport process, the multi-direction definition is usually used to discretize the angle of neutron transport equation, as S_N and P_N methods. But the spatial correlation of the proposed D1QM model in this paper is similar to the standard D1Q2 model, in which the forward and backward vectors of node x point to node x+Δx and node x−Δx, respectively. In the LBM calculation, particles collide at the mesh nodes and travel along the mesh line between nodes. All physical information is stored at the nodes [42].

Fig. 1.

Schematics of the lattice mesh used in the D1QM of lattice Boltzmann method. (a) The one-sided D1QM lattice used in the lattice Boltzmann method with M directions. (b) The D1QM model with two directions (M=2) for one-dimensional isotropic medium. (c) The direction cosine of the polar angle and the interval of two approach direction.

The NTE along any lattice linkage designated by the index m can be written as

μ \frac{\partial ψ_{m} (x, t)}{\partial x} + (\sum_{t} (x) + B) ψ_{m} (x, t) = S_{t} (x, t) + B ψ_{m} (x, t - Δ t),

(7)

where Ψ_m is the distribution function of neutron flux in the m^th direction, and B = 1/(vΔt). The lattice speed along the m^th lattice direction corresponding to the D1QM lattice structure (Fig.1)is defined as e_m = (Δx/μ_m)/Δt. Discrete Eq. (7) along the characteristic directions and keep the right-hand side constant, the evolution equation to Eq. (7) can be given as

\begin{array}{l} ψ_{m} (x + e_{m} Δ t, t + Δ t) \\ = ψ_{m} (x, t) + Δ t \cdot e_{m} (\sum_{t} (x) + B) [\begin{array}{l} \frac{1}{1 + \sum_{t} (x) v Δ t} ψ_{m} (x, t - Δ t) \\ - ψ_{m} (x, t) \end{array}] + Δ t \cdot e_{m} S_{t} (x, t) . \end{array}

(8)

Corresponding to standard lattice Boltzmann-BGK equation [37], the last two terms of right side of Eq. (8) can be treated as the collision and source terms, respectively. So Eq. (8) can be converted to the standard lattice Boltzmann equation as

ψ_{m} (x + e_{m} t, t + Δ t) = ψ_{m} (x, t) + \frac{Δ t}{{\tilde{τ}}_{m}} [ψ_{m}^{eq} (x, t) - ψ_{m} (x, t)] + Δ t \cdot e_{m} S_{t} (x, t),

(9)

where ${\tilde{τ}}_{m}$ = 1/[Σ_t(x)+B] is relaxation time of the collision process, and $ψ_{m}^{eq} (x, t)$ is the equilibrium particle distribution function.

In the LBM, information is exchanged through the collision processes firstly. After that, the particles relax toward the equilibrium distribution and carry information to neighboring lattice nodes in the directions of their characteristics. From the initial condition, the evolution process to the steady-state happens by multiple collisions and streaming, and the collision step is highly influenced by the relaxation time, hence an important parameter to attain the equilibrium parameters.

In the relaxation time, the macroscopic cross-section Σ_t is the reciprocal of the mean free path (mfp for short, which means the average distance traveled by a moving neutron between successive impacts [43, 44]. In this process, 1 mfp=1/Σ_t) of neutron transport in the medium and B is a parameter related to the time step and the velocity of neutron released from medium. Both two parameters represent the effect of diffusion process and relaxation time is a measurement to the strength of diffusion process.

Comparing Eqs. (8) and (9), the equilibrium distribution function is

ψ_{m}^{eq} (x, t) = \frac{1}{1 + \sum_{t} (x) v Δ t} ψ_{m} (x, t - Δ t) .

(10)

This function shows that the equilibrium state is influenced by the distribution of the last time and, like the relaxation time, by macroscopic cross-section and the velocity of neutrons released from medium.

Eq. (9) is the final equation we need to solve in numerical simulation and it is a linear equation. Because of its simple form, Eq. (9) can be used to solve neutron transport problem easily, and to treat more complex neutron transport process by changing its source term and initial/boundary condition.

To calculate the neutron flux, a weighted calculation method is considered. The neutron flux is define as

Φ (x, t) = \sum_{m = 1}^{M} ψ_{m} (x, t) w_{m},

(11)

where w_m is the weight function, in one-dimensional problem, weight function is define as

w_{m} = Δ μ = μ_{m + 1 / 2} - μ_{m - 1 / 2} .

(12)

The proposed LBM model is not restricted to 1-D neutron transport problem. The multi-dimensional and multi-group neutron transport processes can be simulated by the similar method with more coordinate directions and more distribution functions. Owing to the localized computational characteristics of LBM method, the parallel computing can be easily realized and the efficiency can be enhanced especially by GPU acceleration.

Then the neutron flux of each nodes Φ can be numerical calculated by the following steps:

Step 1: Initialize the parameters; choose the number of grid to mesh the computational domain;

Step 2: Define the time step and mesh size;

Step 3: Calculate the relaxation time;

Step 4: Loop at each time step:

(1) Loop for each direction iterations;

a) Calculate the equilibrium distribute functions for each direction with Eq. (10).

b) Implement the collision and streaming processes for each direction with Eq. (9).

c) Impose boundary condition on the boundary nodes.

d) Calculation the neutron flux with Eq. (11). Terminate the global iteration process if the maximum relative error of neutron flux Φ is not bigger than a very small value. Else, go back to step (a).

(2) If the difference of neutron flux between previous and present is smaller than the maximum relative error, terminate the iteration process of time loop. Else, update the neutron flux and go back to Step 1.

Step 5: Save the flux data for the next work.

3. Numerical results

In this section, similar example problems solved to test hybrid method in Ref. [2] are used to validate the new LBM method. Four neutron transport problem in finite homogenous medium are considered. The first problem is a one-dimensional slab with uniformly distributed source, and the next one is a slab with localized uniformly distributed source in the central position. The third is for a localized asymmetric uniformly distributed source, and with the last condition, the localized asymmetric uniformly distributed source that switch-off at one mft (in this process, 1 mft =1 mfp/v).

To test the applicability of different medium, two types of medium are considered. The macroscopic cross section is Σ_t =1 cm⁻¹ for Configurations A and B. The macroscopic scattering cross section is Σ_s = 0.9 cm⁻¹ and 0.5 cm⁻¹ for Configuration A and B, respectively.

For these simulations, the source turns on at t=0⁺, and for the first three problems it keeps for the whole simulation. For the last issue, the source switches off at 1 mft. The source remains at a constant strength level of 1.0 n/cm³ s throughout the turned on situation. The velocity of particles released from the source is 1.0 cm/s. All processes have the vacuum boundary (i.e., the medium is surrounded by a vacuum or nonreflecting region, so that the neutron current density entered from outside equals to zero). Eq. (13) shows the homogeneous initial/boundary conditions used in all of the following examples. And in each simulation, the mesh step is 0.01 mfp and the time step equals to 0.025 mft. The number of discrete directions is M =20. The other non-trivial numerical parameters used in simulations are Δx =0.001 cm, Δt =0.025 s, q =1 cm⁻¹·s⁻¹, and the convergence criterion is set to that the relative error of two iterations is no more than 1×10⁻⁶, for Configurations A and B; while H=10 cm for Configuration A, and H=5 cm for Configuration B.

\begin{array}{l} ψ (- H / 2, μ > 0, t) = 0, \\ ψ (+ H / 2, μ < 0, t) = 0, \\ ψ (x, μ, 0) = 0. \end{array}

(13)

3.1 Uniform volume source

This section refers to a uniform volume source problem with homogeneous conditions, and the geometry of the uniform source is shown in Fig. 2. For Configurations A and B, the slab thickness is 10 mfp and 5 mfp, respectively, and the number of uniform grid is 1000 and 500, respectively.

Fig. 2.

Schematics for 1-D uniform source in H thick slab. Shadow area represents the source region. The left and right side are vacuum boundaries. H=10 mfp for configuration A and H=5 mfp for Configuration B.

The time-dependent flux driven by the uniform volume source for Configuration A is shown in Fig. 3, where the particle densities from the LBM at 1, 2.5, 5, 7.5 and 10 mft are compared with the benchmark [9] results. From Fig.3, our model agrees well with the benchmark, and the maximum relative error between them is below 0.5%. To express the process of neutron transport clearly, the transient particle densities of the first 10 mft with 0.5 mft interval are shown in Fig. 4. At the beginning, the particle density is close to the uniform distribution except the regions nearby the boundary, because in the finite time interval, the volume source and boundary leakage can affect the finite geometric scope, i.e. 1 mft corresponding to 1 mfp. In process of time, a parabola-like particle density profile is formed, because in the central region, more particles are contributed from adjacent volume sources; and in the boundary regions, more particles leak out through the boundaries.

Fig. 3.

Particle densities at various instants of mean free time within a one-dimensional slab of Configuration A with a uniformly distributed source. Solid points are the data of benchmark solutions from Ref. [9] and the lines are the LBM results at 1, 2.5, 5, 7.5 and 10 mft.

Fig. 4.

Transient particle densities of the first 10 mft within a one-dimensional slab of Configuration A with a uniform volume source. The solid lines are calculated with 0.5 mft intervals from bottom to top.

The particle density distribution in the medium converges to steady state, as shown in Fig. 5. The lines are the LBM results, and the symbols are the benchmark solutions of Ref. [9]. Fig. 5(a) shows the convergence trend at 10, 20, 30 and 40 mft. The LBM results and benchmark solutions are consistent with each other. The steady-state particle distribution is shown in Fig. 5(b), where the LBM results are consistent with both benchmark [9] and diffusion approximation solutions [2]. Because the albedo (Σ_s/Σ_t) is large for Configuration A, the diffusion approximation is reasonable.

Fig.5.

Convergence to steady-state of the particle densities of Configuration A with a source of uniform distribution. (a) Converged time-dependent results at 10, 20, 30 and 40 mft. Lines: LBM results, symbols: benchmark solutions [9]. (b) Steady-state comparisons between LBM (line), benchmark (solid square) [9] and diffusion approximation (solid triangle) [2] solutions.

To verify the generality of the LB model, the particle densities driven by the uniform volume source of Configuration B are studied with the same boundary and initial conditions as Configuration A. As shown in Fig. 6(a), the LBM results (lines) at 1, 2 and 5 mft agree well with the transport solutions (symbols) [2]. The steady-state solutions of LBM (solid line), transport solutions (solid circle) [2], and diffusion approximation (solid triangle) [2] are shown in Fig. 6(b). The LBM results are consistent to the transport solutions, with some deviations from the diffusion approximation, which is not reasonable for smaller albedo like Configuration B. From Figs. 5 and 6, it can be seen that the transient and steady-state particle densities of Configuration B are smaller than those of Configuration A, because more particles leak out through the boundaries when the scattering cross-section is smaller. All the relative errors between LBM solutions and benchmark or transport values are less than 0.5% but deviate from diffusion approximation more especially for Configuration B, because the diffusion approximation is not accuracy any more in high scattering medium.

Fig. 6.

Particle densities of Configuration B with a uniform source at 1, 2 and 5 mft. (a) Comparisons of time-dependent particle densities between the LBM results (lines) and the transport solutions (symbols) [2]. (b) Steady-state comparisons between the LBM results (line), the transport solutions (solid circle) [2] and diffusion approximation solutions (solid triangle) [2].

3.2 Symmetric localized source

Transient and steady-state neutron particles transport problem with symmetric localized source are studied in this section. The 1-D slab in thickness of H is shown in Fig. 7 (H = 10 mfp for Configuration A and H = 5 mfp for Configuration B). The symmetric source in width of H/2 is at the center of medium. The vacuum boundary condition is used and the initial particle density is zero distribution. As shown in Fig. 8, the LBM results at first 10 mft agree well with the benchmark solutions.

Fig. 7.

Schematics for 1-D localized symmetric source in H thick slab. H=10 mfp for Configuration A and H=5 mfp for Configuration B.

Fig. 8.

Particle densities at various instants of mean free time in a 1-D slab of Configuration A with a localized symmetric source. Solid symbols are the data of benchmark solutions from Ref. [9] and the lines are the LBM results at 1, 2.5, 5, 7.5 and 10 mft.

Trends of the transient particle densities at the first 10 mft with 0.5 mft interval are shown in Fig. 9. For the several lines in the bottom region, i.e. the initial mean free times, there is a platform in the middle part with localized source, and there are two zero-density zones in the non-source parts nearby the boundaries, because of the effects of source volume and boundary leakage on the local regions. Thereafter, the particles released from the source regions have enough time to get the boundaries, hence the smoothness of density profiles.

Fig. 9.

Transient particle densities of the first 10 mft in a 1-D slab of Configuration A with a symmetric volume source. The solid lines from bottom to top are calculated with 0.5 mft intervals.

The convergence trend of steady-state in symmetric localized source problem with Configuration A is shown in Fig.10. The particle densities in first 40 mft are shown in Fig. 10(a), with each line corresponding to a passage of 10 mft. The LBM results are consistent with the benchmark solutions of Ref. [9]. And the steady-state solutions of LBM, the benchmark values [9] and the diffusion approximations [2] are consistent with others (Fig. 10b).

Fig.10.

Convergence to steady-state of the particle densities of Configuration A with a symmetric distributed source. (a) Converged time-dependent results at 10, 20, 30 and 40 mft. Lines: LBM results, symbols: benchmark solutions [9]. (b) Steady-state comparisons between LBM (line), benchmark (solid circle) [9] and diffusion approximation (solid triangle) [2] solutions.

The same problem of Configuration B with symmetric localized source is simulated at 1, 2 and 5 mft. As shown in Fig.11(a), the LBM results (lines) agree well with the transport solutions (symbols) [2]. In Fig. 11(b) the steady-state solutions of LBM (solid line), transport solutions (solid square) [2], and diffusion approximations (solid triangle) [2] are compared. The LBM results are consistent to the transport solutions, but are lower than the diffusion approximation and higher than those in the region far from the source region. The albedo of Configuration B is smaller, so the diffusion approximation is not exact anymore.

Fig. 11.

Particle densities with a symmetric source of Configuration B. (a) Time-dependent particle densities comparisons between the LBM results (lines) and the transport solutions (symbols) [2] at 1, 2 and 5 mft. (b) Steady-state comparisons between the LBM results (line), the transport solutions (solid square) [2] and diffusion approximation solutions (solid triangle) [2].

For these solutions, the relative errors are less than 0.4% between LBM results and transport values in either transient processes or steady-state conditions. And for steady-state condition, the results of diffusion approximations are different with LBM results slightly for Configuration A but obviously for Configuration B because of the strong scattering effect.

3.3 Asymmetric localized source

In this section, a non-symmetric source-driven problem is considered, with the 1-D slab in thickness of H =10 mfp for Configuration A and H =5 mfp for Configuration B. As shown in Fig. 12, the source region in 1 mfp width is located in the left part of the medium (except with special instructions). Fig.13 shows the time-dependent particle densities at 1, 2 and 5 mft. The LBM results (lines) agree well with the transport solutions (symbols) [2].

Fig. 12.

Schematics for 1-D asymmetric localized source problem (1 mfp width, centered at x=H/4) in H thick slab. Shadow area represents the source region, and left and right boundaries are vacuum boundary. For configuration A, H=10 mfp; for configuration B, H=5 mfp.

Fig. 13.

Particle densities of Configuration A at 1, 2 and 5 mft in a 1-D slab with a localized asymmetric source (width 1 mfp, centered at x=H/4). Solid symbols are the data of transport solutions from reference [2] and the lines are the results of LBM..

The transient particle densities of the first 10 mft with 0.5 mft intervals from LBM are shown in Fig.14. It is interesting to note that the evolution of density profiles is similar to Fig. 9, but not symmetrical to the center of localized source anymore, i.e., the particle densities in the right side are a little larger than that of the symmetric positions in the left side, when the time is long enough. As the right boundary is farther from the source center, more particles are scattered back to the source direction; whereas the left boundary is closer to the source, more particles are leaked out. Also, the particle density is small in the regions adjacent to the right boundary, because optical thickness between the source and right boundary is large, and just a small portion of the particles can reach the boundary. Fig. 15 compares the steady-state particle density of the LBM results (solid line) of Configurations A and B, transport solutions (solid square) and diffusion approximation (solid triangle) [2]. In Fig. 15(a), the LBM results agree well with the transport solutions, but both of them are a little larger than the diffusion approximation in the localized source region; while in Fig. 15(b), the LBM results of Configuration B agree well with the transport solutions [2], but are obviously larger than the diffusion approximation results in the localized source region, and lower than diffusion approximation in the region near the source region.

Fig. 14.

Transient particle densities of the first 10 mft in a 1-D slab of Configuration A with an asymmetric localized source. The solid lines are calculated with 0.5 mft intervals from bottom to top.

Fig. 15.

Comparisons of steady-state particle density of Configuration A (a) and Configuration B (b) with the asymmetric source. Lines, the LBM results; Solid square, the transport solutions; Solid triangle, the diffusion approximation solutions [2].

Finally, the asymmetric localized source switch-off at 1 mft is studied using the geometry similar to Fig. 12, i.e. the 1-D slab thickness of H = 10 mfp for Configuration A and H =5 mfp for configuration B, but the asymmetric localized source is 1mfp in width for both Configurations A and B. The transient particle densities of LBM results (lines) agree well with the transport solutions (symbols) [2] in all domains (Fig. 16).

Fig. 16.

Transient particle densities of Configuration A (a) and Configuration B (b) at various instants of time with asymmetric pulse source (1 mfp width, centered at x=H/4) switch-off at 1 mft. Solid symbols: the transport solutions [2]; Lines: LBM results.

Like the previous cases, the relative errors are less than 0.5% between LBM results and transport values, except the last two comparisons, with the maximum error of <1%, since the neutron flux is very small.

4. Conclusion

A novel lattice Boltzmann model for one dimension time-dependent neutron particle transport problem by solving linear Boltzmann equation is proposed, and the detailed derivation process is presented in this paper. Typical examples of neutron particle transport in homogeneous medium driven by symmetric and asymmetric source are studied. Both time-dependent and steady-state LBM results agree well with the benchmark values in the references. The relative error between LBM solutions and transport values of most cases are less than 0.5% and the maximum error is less than 1%. Comparisons show that the proposed LBM model is precise for neutron transport. The LBM method will be further developed as a powerful computational tool for neutron transport problems. This method can be developed to solve 2-D and 3-D neutron transport problems, and parallel computation of the models can be easily realized since the strong locality of the relevant variables.

References

[1]

W.A. Wieselquist, D.Y. Anistratov, J.E. Morel,

A cell-local finite difference discretization of the low-order quasidiffusion equations for neutral particle transport on unstructured quadrilateral meshes

. J. Comput. Phys. 273, 343-357 (2014). doi: 10.1016/j.jcp.2014.05.011