Large Eddy Simulation of the Inlet Cross-Flow in the CiADS Heat Exchanger Using the Lattice Boltzmann Method

: The liquid lead-bismuth eutectic (LBE) is the coolant of the heat exchanger in China initiative Accelerator Driven System, which may have a risk of structural failure due to the washout of the coolant in the inlet of the heat exchanger. It is necessary to study the mechanical properties of the heat exchanger bundles of CiADS, especially the fatigue life of the bundle structure in the transverse ﬂow of the LBE. The numerical simulation is the Lattice Boltzmann method combined with the large eddy simulation by Python codes. The velocity distribution of the ﬂow ﬁeld and the time domain characteristics of the heat exchanger bundles’ force are calculated, and the frequency domain characteristics of the heat exchanger bundles’ vibration are obtained by Fourier transform. The bundles vibrate at high cycle fatigue in turbulent ﬂow at high Reynolds number. The transverse ﬂow of LBE does not affect the fatigue life of the bundle. No structural failure occurs in the CiADS heat exchanger due to the transverse ﬂow of LBE.


Introduction
Nuclear energy has significant advantages in stability, reliability, and zero carbon dioxide emissions; see Chen [1].It is considered a strategic choice to solve the future energy supply and ensure sustainable economic and social development; see Clulow [2].According to Drace [3] and Yang [4], The lack of public awareness of nuclear energy has resulted in a low level of acceptance of nuclear energy.The accelerator-driven subcritical system (ADS) drives subcritical reactors with high energy provided by the accelerator to produce high-throughput and broad-spectrum spallation neutrons, which will not create criticality accidents in the common critical reactors.With the advantages of ADS reactors being subcritical systems with inherent safety, they can contribute to increasing public acceptance of nuclear energy.
The Accelerator Driven System is an advanced technological approach to the sustainable development of nuclear energy and an effective solution for achieving carbon peak and neutrality targets; see Zhan [5].The Nobel Prize owner, Professor C. Rubbia [6], has proposed the European EUROTRANS program since 1986 to promote ADS technology in the future.The Belgian Nuclear Research Centre has been committed to the long-term development of MYRRHA (Multi-purpose hybrid research reactor for high-tech applications), which is the world's first large-scale ADS [7]; see Abderrahim [7].From 2011 to 2016 in China, the strategic pilot project "Future Advanced Fission Energy-ADS Transmutation Systems" was launched by the Chinese Academy of Sciences, which is also led by the Institute of Modern Physics of the Chinese Academy of Sciences (IMPCAS), focusing on solving individual vital technical problems in the ADS system.The Chinese government has decided to support IMPCAS in leading the design and construction of the integrated ADS project, the so-called CiADS, backed by the national primary science and technology infrastructure project of the "12th Five-Year Plan"; see Su [8] and Gu [9].The China initiative Accelerator Driven System (CiADS) is also devoted to establishing an experimental and verification ADS facility globally [10]; see Xunchao [10].
Liquid metals are widely considered coolants with good thermal-hydraulic characteristics for many energy systems, such as fast reactors and subcritical reactors; see Shi [11].According to Yanyi [12], The lead-cooled fast reactor (LFR) is applied by the China Initiative Accelerator Driven System due to its good nuclear waste transmutation and nuclear fuel breading ability.The China initiative Accelerator Driven System aims to be the world's first megawatt ADS device.Benefiting from the advantages of the LBE with a good thermalhydraulic capability, high boiling point, and low chemical inertness, an LBE-cooled fast reactor is chosen as the subcritical reactor of the CiADS; see Fan [13].The main systems include a reactor core system, reactor coolant system, reactor auxiliary system, dedicated safety facilities, nuclear island common support system, etc. Thermal and hydraulic characteristics of the LBE will affect the safety performance of the CiADS Lead-bismuth-cooled reactor.Although the accuracy and confidence of experimental measurements are higher than the numerical calculation method, due to the corrosive and opaque nature of the LBE, a number of studies of the LBE are based on CFD methods, according to Qing [14], Pavel [15] and Zhan [16].
In recent years, there have been several accidents in which the heat exchanger tubes were damaged by vibration, mostly due to the vibration of the tubes caused by the coolant transverse flow, see Jo [17].The lateral flow of LBE coolant at the heat exchanger inlet can cause flow separation and create vortices on the back side of the heat exchanger tubes.This causes continuous vibration of the heat exchanger bundle, resulting in vibration fatigue and reduced heat exchanger bundle life.The flow characteristics of the CiADS heat exchanger will also differ from conventional heat exchangers due to the physical characteristics of LBE, which are very different from those of water, e.g., the density of LBE is approximately ten times that of water, see Su [18].The flow characteristics of LBE will directly affect the force characteristics of the heat exchanger tubes, so it is important to study the force characteristics of the heat exchanger tubes.

Materials and Methods
The structural strength of the CiADS heat exchanger must meet the design requirements to ensure the structural integrity of the CiADS heat exchanger, ensuring the integrity of the radioactive material containment boundary and making the CiADS system safe.
The flow of lead-bismuth coolant in the CiADS heat exchanger is shown in Figure 1; after flowing from the hot legs, the coolant enters the channel between the casing and the tube bundle through the coolant inlet receiver at the top of the casing, flows down through the outer heat exchanger tube, transfers the heat to the second circuit cooling water, cools and then flows out through the outlet receiver to the core, where it is pumped by the main pump.After cooling, the water flows out of the outlet receiver and into the core cold pool via the main pump.The low-temperature water on the second circuit side enters the upper inlet of the CiADS heat exchanger from the second circuit main pipe, where it is heated by the first circuit coolant and then flows back to the second circuit.

Computational Methods
The Lattice Boltzmann Method (LBM) is a mesoscopic numerical method that solves discrete distribution functions by collision and streaming steps to obtain macroscopic physical parameters, according to Marie [19], Hartono [20], and Kruger [21].With the advantages of the explicit evolutionary process and simple algorithmic, LBM has been a vibrant research topic in a wide range of applications in porous media seeing Farshad [22], particle flow seeing Afra [23], Delouei [24], non-Newtonian fluids seeing Kefayati [25], and turbulent flows currently seeing Kuwata [26].

Original Lattice Boltzmann Method
In order to facilitate the understanding of the algorithms of the overall lattice Boltzmann method, the single-relaxation time (SRT) collision model is used as a common collision model describing the framework of LBM; see An [27].In what follows, a brief description of the original LBM is presented; this introduction will be later used to implement the MRT-LBM and LES-LBM methodologies.
The lattice Boltzmann equation with the SRT collision operator is where  is the unit velocities vector along discrete direction α,  is the spatial position,  is the steaming time step,  is the relaxation time,  (, ) is the distribution function in the direction of α at time of  in the position of , and  (, ) is the equilibrium distribution function in the direction of α at time of , in the position of .Notice that all parameters presented in this equation, as well as the ones introduced in the rest of the equations presented, are non-dimensional.
A class of widely used discrete models is DnQb, which indicates that there are b discrete velocities at each grid point on the n-dimensional space.The D2Q9 model is a typical DnQb model with nine discrete velocities in two-dimensional space.Figure 2 shows the discrete velocities of the D2Q9 model employed in all simulations presented in this paper.

Computational Methods
The Lattice Boltzmann Method (LBM) is a mesoscopic numerical method that solves discrete distribution functions by collision and streaming steps to obtain macroscopic physical parameters, according to Marie [19], Hartono [20], and Kruger [21].With the advantages of the explicit evolutionary process and simple algorithmic, LBM has been a vibrant research topic in a wide range of applications in porous media seeing Farshad [22], particle flow seeing Afra [23], Delouei [24], non-Newtonian fluids seeing Kefayati [25], and turbulent flows currently seeing Kuwata [26].

Original Lattice Boltzmann Method
In order to facilitate the understanding of the algorithms of the overall lattice Boltzmann method, the single-relaxation time (SRT) collision model is used as a common collision model describing the framework of LBM; see An [27].In what follows, a brief description of the original LBM is presented; this introduction will be later used to implement the MRT-LBM and LES-LBM methodologies.
The lattice Boltzmann equation with the SRT collision operator is where e α is the unit velocities vector along discrete direction α, x is the spatial position, δ t is the steaming time step, τ 0 is the relaxation time, f α (x, t) is the distribution function in the direction of α at time of t in the position of x, and f eq α (x, t) is the equilibrium distribution function in the direction of α at time of t, in the position of x.Notice that all parameters presented in this equation, as well as the ones introduced in the rest of the equations presented, are non-dimensional.
A class of widely used discrete models is DnQb, which indicates that there are b discrete velocities at each grid point on the n-dimensional space.The D2Q9 model is a typical DnQb model with nine discrete velocities in two-dimensional space.Figure 2 shows the discrete velocities of the D2Q9 model employed in all simulations presented in this paper.
where  =  / ,  is the grid spacing, and the equilibrium distribution function can be expressed as where ρ and  are the macroscopic density and velocity of the fluid, and  is the weighting factor, which is expressed as  = 0,  = 0 4 9 ⁄ ,  = 1,2,3,4 1 9 ⁄ ,  = 5,6,7,8 Macroscopic densities and velocities are calculated by the equilibrium distribution function: The original lattice Boltzmann method is usually applied to simulate incompressible flow, with the Reynolds number commonly limited to low values.

Multiple-Relaxation-Time Lattice Boltzmann Method
The multi-relaxation-time (MRT) collision model demonstrates superiority in numerical stability and accuracy compared to the SRT model at high Reynolds flows, according to Jahanshaloo [28] and Luo [29].The MRT collision process is accomplished using the operation of matrices and vectors.
The evolution equation of LBM with the MRT model is where  = diag( ,  ,  ,  ,  ,  ,  ,  ,  , ) denotes the diagonal collision relaxation matrix with the expression. is the orthogonal transformation matrix, velocity space and moment space are achieved using linear transformation: D2Q9 model discrete speed expressions: where c = δ x /δ t , δ x is the grid spacing, and the equilibrium distribution function can be expressed as where ρ and u are the macroscopic density and velocity of the fluid, and ω α is the weighting factor, which is expressed as Macroscopic densities and velocities are calculated by the equilibrium distribution function: The original lattice Boltzmann method is usually applied to simulate incompressible flow, with the Reynolds number commonly limited to low values.

Multiple-Relaxation-Time Lattice Boltzmann Method
The multi-relaxation-time (MRT) collision model demonstrates superiority in numerical stability and accuracy compared to the SRT model at high Reynolds flows, according to Jahanshaloo [28] and Luo [29].The MRT collision process is accomplished using the operation of matrices and vectors.
The evolution equation of LBM with the MRT model is where S = diag(s 0 , s 1 , s 2 , s 3 , s 4 , s 5 , s 6 , s 7 , s 8 ,) denotes the diagonal collision relaxation matrix with the expression.M is the orthogonal transformation matrix, velocity space and moment space are achieved using linear transformation: Converting the distribution function and its equilibrium state to a distribution function in moment space: where ρ is the density, e is the kinetic energy, ε is the quadratic of the kinetic energy, j x and q x are the momentum and energy fluxes in the x-direction, j y and q y are the momentum and energy fluxes in the y-direction, and p xx and p yy are the diagonal and diagonal viscous stress tensor, respectively.

Large-Eddy-Simulation Lattice Boltzmann Method
The large-eddy-simulation lattice Boltzmann method is another way to improve the numerical stability from the original lattice Boltzmann method.For high Reynolds number turbulence flows, it is important to choose a reasonable turbulence model of numerical simulation, according to Krafczyk [30], Dong [31], Rahmati [32], and Chengwu [33].Three turbulent models can be coupled with LBM, which are Direct Numerical Simulation (DNS), Large Eddy Simulation (LES), and Reynolds Average Navier-Stokes (RANS), respectively.The DNS can obtain the highest numerical accuracy, but it pays numerous computing costs.The RANS can save computing costs, but it cannot capture the transient fluctuation of the turbulence.A turbulence model needs to be introduced into the computational model, and the MRT-LBM has better numerical stability for turbulence simulations.It is well understood that the numerical stability of LBGK models is substantially inferior to the corresponding multiple relaxation time (MRT) LBE models, see Haussmann [34].The significant improvement in numerical stability by using the MRT-LBE models can directly result in a drastic gain in computational efficiency.Thus, it is natural to incorporate MRT-LBE models with the LES to improve the capability of turbulent flows.
The large eddy simulation is processed using a filter function to obtain the variables, where τ ij is the subgrid stress, with the expression: Introducing the Smagorinsky model for the solution, the subgrid stress is calculated as: S is the modal length of the large-scale strain rate tensor, ν t is the vortex viscosity coefficient and C s is the Smagorinsky constant.
where S ij is the large-scale strain rate tensor, expressed as: Introducing the Smagorinsky sublattice stress model into the LBM model, the total relaxation time varies, and the LBM equation evolves in the form of:

Calculation Models
It is common in engineering practice to see a cluster of cylinders arranged in a cluster, such as a nuclear reactor core assembly, a tube array system in a heat exchanger, a cluster of high-rise buildings in a wind farm, a diagonal cable in a large span bridge, an overhead transmission cable line, a cluster of cooling towers, etc.The heat exchanger tube bundles consist of multiple cylinders arranged in a column group.When the fluid flows through, due to the mutual interference between multiple cylinders and the strong nonlinearity of the interference, the composition of the cylinders group of single cylinders subject to the fluid force characteristics and flow structure and single-column perturbation is significantly different.As the basic constituent unit in the column group, the three-column is also a common structural form in engineering, and the wake structure of the three-cylinder disturbance is much more complicated compared with that of the single-cylinder.The simulation is carried out by a bundle of three cylinders in a two-dimensional summary.The location of cylinders #1, 2, and 3 is shown in Figure 3.
of high-rise buildings in a wind farm, a diagonal cable in a large span bridge, an overhead transmission cable line, a cluster of cooling towers, etc.The heat exchanger tube bundles consist of multiple cylinders arranged in a column group.When the fluid flows through, due to the mutual interference between multiple cylinders and the strong nonlinearity of the interference, the composition of the cylinders group of single cylinders subject to the fluid force characteristics and flow structure and single-column perturbation is significantly different.As the basic constituent unit in the column group, the three-column is also a common structural form in engineering, and the wake structure of the three-cylinder disturbance is much more complicated compared with that of the single-cylinder.The simulation is carried out by a bundle of three cylinders in a two-dimensional summary.The location of cylinders #1, 2, and 3 is shown in Figure 3.
Compared to the axial flow of the fluid, the cross-flow effect is strongest at the inlet, and the effect of the disturbance is most pronounced.The transverse scour of the fluid has a greater effect on the oscillations of the tube bundle; therefore, in this paper, the inlet section of the heat exchanger where the transverse flow is strongest is calculated, and the model of the flow is considered to be reduced to the transverse flow of a two-dimensional triangularly arranged cylindrical tube bundle.Both 2D and 3D simulations could obtain satisfactory results of the St number in agreement with the experiments.The St number is the key parameter, which is conclusive for the calculating of the frequency of the drag and lift forces.The basic flow parameters of the vortex shedding are acceptable and obtained in both 2D and 3D simulations.The drag coefficients calculated in 2D are larger, which means that the drag coefficients obtained from 2D have a numerical margin compared with 3D; see Hongtao [35].Li et al. [36] compared a three-dimensional model of the heat exchanger with a two-dimensional model.The results show that the maximum errors in pressure drop and flow velocity for different axial sections are 0.1% and 14%, respectively, which are negligible.The two-dimensional simplification is reasonable and feasible, and the flow problem of the heat exchanger is calculated using the two-dimensional model.Considering the simple structure of the two-dimensional model and the small number of meshes required for the calculation, the two-dimensional model will be used for the simulation of the heat exchanger in the study of this paper.Additionally, the model of the flow is simplified to the transverse flow of a two-dimensional triangularly arranged cylindrical tube bundle.
The numerical simulation in this article is non-dimensional.Reynolds number is an important non-dimensional parameter that represents the ratio of the inertial force of the fluid to the viscous force, which is a key parameter for studying the subject of the flow Compared to the axial flow of the fluid, the cross-flow effect is strongest at the inlet, and the effect of the disturbance is most pronounced.The transverse scour of the fluid has a greater effect on the oscillations of the tube bundle; therefore, in this paper, the inlet section of the heat exchanger where the transverse flow is strongest is calculated, and the model of the flow is considered to be reduced to the transverse flow of a two-dimensional triangularly arranged cylindrical tube bundle.Both 2D and 3D simulations could obtain satisfactory results of the St number in agreement with the experiments.The St number is the key parameter, which is conclusive for the calculating of the frequency of the drag and lift forces.The basic flow parameters of the vortex shedding are acceptable and obtained in both 2D and 3D simulations.The drag coefficients calculated in 2D are larger, which means that the drag coefficients obtained from 2D have a numerical margin compared with 3D; see Hongtao [35].Li et al. [36] compared a three-dimensional model of the heat exchanger with a two-dimensional model.The results show that the maximum errors in pressure drop and flow velocity for different axial sections are 0.1% and 14%, respectively, which are negligible.The two-dimensional simplification is reasonable and feasible, and the flow problem of the heat exchanger is calculated using the two-dimensional model.
Considering the simple structure of the two-dimensional model and the small number of meshes required for the calculation, the two-dimensional model will be used for the simulation of the heat exchanger in the study of this paper.Additionally, the model of the flow is simplified to the transverse flow of a two-dimensional triangularly arranged cylindrical tube bundle.
The numerical simulation in this article is non-dimensional.Reynolds number is an important non-dimensional parameter that represents the ratio of the inertial force of the fluid to the viscous force, which is a key parameter for studying the subject of the flow past cylinders.The flow pattern is mainly affected by the Reynolds number.Within a specific range of Reynolds numbers, different flow patterns can be distinguished by Reynolds numbers.Animasaun et al. [37] disclosed that laminar flow patterns degrade into turbulent eddies and vortices at high Reynolds numbers due to the dominance of nonlinear interactions between fluid particles.A wide variety of fluid systems, from airflow over an airplane wing to water movement in rivers and seas, exhibit this transition to turbulence as a universal feature.
The physical properties of LBE can be obtained from Table 1, as the temperature of LBE in the inlet of the CiADS heat exchanger is 380 K.The key non-dimensional parameter Reynolds number is 43,936.55.The unit of physical quantities in the LBM model is lu (lattice unit) after the dimensionless process.The time can be obtained by T = υ /υ D/D 2 in which parameters without prime are practical quantities, while parameters with prime are lattice non-dimensional parameters.

Calculation Results When Re = 200
The velocity distribution clouds for a single cylindrical disturbance with Reynolds number Re = 200 are shown in Figure 4, which represents the velocity distribution clouds at different moments in an oscillation period.A periodic vortex shedding "Carmen vortex street" is formed downstream of the cylinder, with two rows of vortices regularly and equidistantly arranged in the wake area, each vortex in the middle of the other vortex row, with the vortices in the same vortex row rotating in the same direction, but with the opposite vortex row rotating.

Calculation Results When Re = 200
The velocity distribution clouds for a single cylindrical disturbance with Reynolds number Re = 200 are shown in Figure 4, which represents the velocity distribution clouds at different moments in an oscillation period.A periodic vortex shedding "Carmen vortex street" is formed downstream of the cylinder, with two rows of vortices regularly and equidistantly arranged in the wake area, each vortex in the middle of the other vortex row, with the vortices in the same vortex row rotating in the same direction, but with the opposite vortex row rotating.The drag and lift coefficients of the cylinder are shown in Figure 5.The cylinder is sustained with lift and drag forces periodically.The results of the numerical simulation in this article are compared with the results available in the references.The drag coefficients are 1.585 and 1.575 for this article and the reference, respectively.The frequency of the drag force is twice as much as the lift force due to the influence of the Carmen vortices.The lift coefficients ∆C l are 1.70 and 1.67 for the manuscript and the reference, respectively.The calculations are in agreement with the article by Chen [38].

CiADS Heat Exchanger Calculation Results
The flow simulation of the inlet cross-flow in the CiADS Heat Exchanger us Large-eddy-simulation lattice Boltzmann method at Reynolds number Re = 43936 complished using Python code.The velocity cloud diagram of the three cylinders is in Figure 6.The flow is no longer a simple laminar flow at high Reynolds numbe and the flow field is turbulent.As seen in Figure 7, vortices are formed on the cy surface, but the wake structure becomes irregular, the vortices grow downstream w main flow, and the cylinders vibrate, subjecting to the fluid excitation force.The v behind cylinder #1 and cylinder #2 deviate from the direction of the inlet center-l to the high velocity and low pressure, which deflects the wake outwards.Due to th ding and merging of the vortices, the flow field produces periodic changes in pres that its wake takes on the form of a jet oscillation.The vorticity diagram shows t surface of the cylinder generates small vortices, which grow in size as they flow stream, with some of the vortices growing and breaking up and others merging wi other.

CiADS Heat Exchanger Calculation Results
The flow simulation of the inlet cross-flow in the CiADS Heat Exchanger using the Large-eddy-simulation lattice Boltzmann method at Reynolds number Re = 43,936.5 is accomplished using Python code.The velocity cloud diagram of the three cylinders is shown in Figure 6.The flow is no longer a simple laminar flow at high Reynolds number flow, and the flow field is turbulent.As seen in Figure 7, vortices are formed on the cylinders' surface, but the wake structure becomes irregular, the vortices grow downstream with the main flow, and the cylinders vibrate, subjecting to the fluid excitation force.The vortices behind cylinder #1 and cylinder #2 deviate from the direction of the inlet centerline due to the high velocity and low pressure, which deflects the wake outwards.Due to the shedding and merging of the vortices, the flow field produces periodic changes in pressure so that its wake takes on the form of a jet oscillation.The vorticity diagram shows that the surface of the cylinder generates small vortices, which grow in size as they flow downstream, with some of the vortices growing and breaking up and others merging with each other.

CiADS Heat Exchanger Calculation Results
The flow simulation of the inlet cross-flow in the CiADS Heat Exchanger using the Large-eddy-simulation lattice Boltzmann method at Reynolds number Re = 43936.5 is accomplished using Python code.The velocity cloud diagram of the three cylinders is shown in Figure 6.The flow is no longer a simple laminar flow at high Reynolds number flow, and the flow field is turbulent.As seen in Figure 7, vortices are formed on the cylinders' surface, but the wake structure becomes irregular, the vortices grow downstream with the main flow, and the cylinders vibrate, subjecting to the fluid excitation force.The vortices behind cylinder #1 and cylinder #2 deviate from the direction of the inlet center-line due to the high velocity and low pressure, which deflects the wake outwards.Due to the shedding and merging of the vortices, the flow field produces periodic changes in pressure so that its wake takes on the form of a jet oscillation.The vorticity diagram shows that the surface of the cylinder generates small vortices, which grow in size as they flow downstream, with some of the vortices growing and breaking up and others merging with each other.The interaction of the underdeveloped wake of the upstream cylinders, the turbulence in the gap among cylinders, and the shedding vortex of the three cylinders make for a very complex flow field, which also makes the forces of the cylinders complicated.The force analysis of the three cylinders in the time domain is shown in Figure 6.The force of cylinder #3 downstream is more complex than the two cylinders upstream due to the interaction of the wake of the upstream cylinder #1 and cylinder #2.The turbulence in the cylinder gap and the shedding of the vortex so that the vibrations occurring in the lift direction do not have a distinct vibration period and are not of equal amplitude.It can be seen that the drag coefficients of cylinder #1 and cylinder #2 upstream are almost the same, which are higher than that of cylinder #3 downstream.The reason for this is that the cylinders upstream play a role in shielding the cylinders downstream.
The interaction of the underdeveloped wake of the upstream cylinders, the turbulence in the gap among cylinders, and the shedding vortex of the three cylinders make for a very complex flow field, which also makes the forces of the cylinders complicated.The force analysis of the three cylinders in the time domain is shown in Figure 8.The force of cylinder #3 downstream is more complex than the two cylinders upstream due to the interaction of the wake of the upstream cylinder #1 and cylinder #2.The turbulence in the cylinder gap and the shedding of the vortex so that the vibrations occurring in the lift direction do not have a distinct vibration period and are not of equal amplitude.It can be seen that the drag coefficients of cylinder #1 and cylinder #2 upstream are almost the same, which are higher than that of cylinder #3 downstream.The reason for this is that the cylinders upstream play a role in shielding the cylinders downstream.
The vibration frequencies of the three cylinders are obtained using Fourier analysis, as shown in Figure 9.The drag coefficients of the cylinders at different positions are not the same.The maximum amplitude of the drag coefficients of cylinder #1 corresponds to a frequency of 398 Hz, the maximum amplitude of the drag coefficients of cylinder #2 corresponds to a frequency of 398 Hz, and the maximum amplitude of the drag coefficients of cylinder #3 corresponds to a frequency of 100 Hz.The turbulent wake of the cylinders upstream is less effective.The drag coefficients of the cylinders upstream are lower due to the shielding effect on the downstream cylinder, and the turbulent wake of the cylinders upstream has a reduced drag on the cylinder downstream.The vortex shedding of cylinder #1 and cylinder #2 upstream gradually develops fully, and the wake vortices from the upstream cylinder shedding impact the downstream cylinder #3, resulting in a gradual increase in the lift and drag forces of cylinder #3 downstream.The direct crossflow impact on the upstream cylinder and the heat transfer tubes in the back row influence the development of its vortex.The superposition of the two parts of the vibration causes The interaction of the underdeveloped wake of the upstream cylinders, the turbulence in the gap among cylinders, and the shedding vortex of the three cylinders make for a very complex flow field, which also makes the forces of the cylinders complicated.The force analysis of the three cylinders in the time domain is shown in Figure 6.The force of cylinder #3 downstream is more complex than the two cylinders upstream due to the interaction of the wake of the upstream cylinder #1 and cylinder #2.The turbulence in the cylinder gap and the shedding of the vortex so that the vibrations occurring in the lift direction do not have a distinct vibration period and are not of equal amplitude.It can be seen that the drag coefficients of cylinder #1 and cylinder #2 upstream are almost the same, which are higher than that of cylinder #3 downstream.The reason for this is that the cylinders upstream play a role in shielding the cylinders downstream.
The interaction of the underdeveloped wake of the upstream cylinders, the turbulence in the gap among cylinders, and the shedding vortex of the three cylinders make for a very complex flow field, which also makes the forces of the cylinders complicated.The force analysis of the three cylinders in the time domain is shown in Figure 8.The force of cylinder #3 downstream is more complex than the two cylinders upstream due to the interaction of the wake of the upstream cylinder #1 and cylinder #2.The turbulence in the cylinder gap and the shedding of the vortex so that the vibrations occurring in the lift direction do not have a distinct vibration period and are not of equal amplitude.It can be seen that the drag coefficients of cylinder #1 and cylinder #2 upstream are almost the same, which are higher than that of cylinder #3 downstream.The reason for this is that the cylinders upstream play a role in shielding the cylinders downstream.
The vibration frequencies of the three cylinders are obtained using Fourier analysis, as shown in Figure 9.The drag coefficients of the cylinders at different positions are not the same.The maximum amplitude of the drag coefficients of cylinder #1 corresponds to a frequency of 398 Hz, the maximum amplitude of the drag coefficients of cylinder #2 corresponds to a frequency of 398 Hz, and the maximum amplitude of the drag coefficients of cylinder #3 corresponds to a frequency of 100 Hz.The turbulent wake of the cylinders upstream is less effective.The drag coefficients of the cylinders upstream are lower due to the shielding effect on the downstream cylinder, and the turbulent wake of the cylinders upstream has a reduced drag on the cylinder downstream.The vortex shedding of cylinder #1 and cylinder #2 upstream gradually develops fully, and the wake vortices from the upstream cylinder shedding impact the downstream cylinder #3, resulting in a gradual increase in the lift and drag forces of cylinder #3 downstream.The direct cross-flow impact on the upstream cylinder and the heat transfer tubes in the back row influence the development of its vortex.The superposition of the two parts of the vibration causes the amplitude to decay over a longer period.The frequency of the lift of downstream cylinder #3 is lower compared to the frequency of the cylinders upstream.The strong interference effect between the three cylinders makes the disturbed vortex shedding appear nonlinear.The pressure distribution of the cylinders upstream is symmetrical due to the symmetrical geometric arrangement and the flow.The more the upstream cylinder is affected by the cross-flow impact, the greater the drag coefficient in comparison to the downstream cylinder.The vibration balance of the downstream cylinder is positioned close to the concenter of the flow; this is due to the presence of the upstream heat transfer tube weakening the impact of the cross-flow on the downstream tube, making the effect of the cross-flow impact on the vibration in the direction of resistance weaker.The vortex development behind the downstream heat transfer tube is not suppressed, and the lift effect is enhanced.
Sustainability 2023, 15, x FOR PEER REVIEW 10 of 14 the amplitude to decay over a longer period.The frequency of the lift of downstream cylinder #3 is lower compared to the frequency of the cylinders upstream.The strong interference effect between the three cylinders makes the disturbed vortex shedding appear nonlinear.The pressure distribution of the cylinders upstream is symmetrical due to the symmetrical geometric arrangement and the flow.The more the upstream cylinder is affected by the cross-flow impact, the greater the drag coefficient in comparison to the downstream cylinder.The vibration balance of the downstream cylinder is positioned close to the concenter of the flow; this is due to the presence of the upstream heat transfer tube weakening the impact of the cross-flow on the downstream tube, making the effect of the cross-flow impact on the vibration in the direction of resistance weaker.The vortex development behind the downstream heat transfer tube is not suppressed, and the lift effect is enhanced.

Fatigue Life Analysis
Engineering needs fatigue life analysis because it can forecast a component's lifespan under cyclic loads, assuring the dependability and safety of goods and buildings.Fatigue the amplitude to decay over a longer period.The frequency of the lift of downstream cylinder #3 is lower compared to the frequency of the cylinders upstream.The strong interference effect between the three cylinders makes the disturbed vortex shedding appear nonlinear.The pressure distribution of the cylinders upstream is symmetrical due to the symmetrical geometric arrangement and the flow.The more the upstream cylinder is affected by the cross-flow impact, the greater the drag coefficient in comparison to the downstream cylinder.The vibration balance of the downstream cylinder is positioned close to the concenter of the flow; this is due to the presence of the upstream heat transfer tube weakening the impact of the cross-flow on the downstream tube, making the effect of the cross-flow impact on the vibration in the direction of resistance weaker.The vortex development behind the downstream heat transfer tube is not suppressed, and the lift effect is enhanced.

Fatigue Life Analysis
Engineering needs fatigue life analysis because it can forecast a component's lifespan under cyclic loads, assuring the dependability and safety of goods and buildings.Fatigue

Fatigue Life Analysis
Engineering needs fatigue life analysis because it forecast a component's lifespan under cyclic loads, assuring the dependability and safety of goods and buildings.Fatigue analysis is usually based on the S-N curve of the material, the cumulative damage relationship, and the damage criterion to determine the fatigue life of the structure.Vibratory fatigue analysis is a method of fatigue life determination that uses the Miner linear accumulation theory to assess the fatigue life based on structural response analysis and finding the applicable S-N curve using the Fourier transform to find the basic frequency of the vibration, according to the Shun [39] and Tridello [40].This paper uses the widely used Miner theory to analyze the fatigue life of the heat transfer tubes of CiADS heat exchangers.
According to the damage rule of linear accumulation, when the parameter for measuring destruction D t = 1, the structural strength failure occurs.Based on the tube material of CiADS heat transfer, according to the S-N curve, the constant c = 6.4 × 10 −8 and b = 4.
∆ε is the maximum strain difference occurring in one cycle, and the expression for ∆ε is where A is the amplitude and the corresponding number of cycle cycles n(∆ε i ) can be expressed as n(∆ε i ) = f i t i n(∆ε i ) is the number of periods in which the alternating strain occurs for a range of variations of ∆ε i .According to the damage rule of linear accumulation, when when D t reaches 1, damage to the structure will occur.where f i is the frequency corresponding to the i-th amplitude and t i is the corresponding time.The value of D t in one year is D t = ∑ i f i t i A i 4 (D/L) 8   6.745 × 10 −12 (21) where A i is the i-th amplitude in the pipe, let T i be the number of hours per day that the pipe vibrates at the i-th amplitude, and the expression for T i is Then, the fatigue life of the heat transfer tube in years is obtained as The fatigue life of the tube bundle was assessed based on the spectral response values of the lift resistance of the three cylinders.The calculated result of the force on the cylinders is a maximum of 1.45 MPa, which is much smaller than the threshold value of 150 MPa for fatigue life stress calculation, while the calculated threshold value of 100 MPa for the ultrahigh circumference is relatively small, the tube bundle will not suffer from life depreciation due to fatigue of the structural material under the washout of the lead-bismuth alloy, which will affect the service life of the heat exchanger.Therefore, in accordance with the fatigue life prediction method, the calculations show that the fatigue life of the heat exchanger tube bundle under the lateral flushing of the lead-bismuth alloy meets the design life of the heat exchanger, and the heat exchanger will not suffer structural failure due to lateral flow.

Conclusions
The fatigue life of the bundle structure of the CiADS heat exchanger in the transverse flow of the LBE is studied in this paper, which is the criterion to assess the mechanical properties of the bundle meeting the operation requirements.The flow characteristics at the inlet of the CiADS heat exchanger at a high Reynolds number are investigated using the Large-eddy-simulation lattice Boltzmann method.The fatigue life of the bundles is analyzed based on the frequency domain characteristics of the bundles obtained using Fourier analysis, and the following conclusions are obtained: The maximum stress on the bundle is 1.45 MPa, which is smaller than the threshold value of fatigue life stress calculation.The transverse flow of LBE does not affect the fatigue life of the bundle.No structural failure occurs in the CiADS heat exchanger due to the transverse flow of LBE.

Figure 2 .
Figure 2. Discrete velocities of the D2Q9 model.D2Q9 model discrete speed expressions:

Figure 3 .
Figure 3.The location of the three cylinders.

Figure 3 .
Figure 3.The location of the three cylinders.

Figure 4 .
Figure 4. Velocity profile of the flow past a cylinder at Re = 200.(A) at time of 1/4 T; (B) at time of 1/2 T; (C) at time of 3/4 T; (D) at time of 4/4 T. The drag and lift coefficients of the cylinder are shown in Figure 5.The cylinder is sustained with lift and drag forces periodically.The results of the numerical simulation in this article are compared with the results available in the references.The drag coefficients are 1.585 and 1.575 for this article and the reference, respectively.The frequency of the

Figure 4 .
Figure 4. Velocity profile of the flow past a cylinder at Re = 200.(A) at time of 1/4 T; (B) at time of 1/2 T; (C) at time of 3/4 T; (D) at time of 4/4 T.

Sustainability 2023 ,
15,  x FOR PEER REVIEW drag force is twice as much as the lift force due to the influence of the Carmen v The lift coefficients ∆C are 1.70 and 1.67 for the manuscript and the reference, tively.The calculations are in agreement with the article by Chen[38].

Figure 5 .
Figure 5. Drag coefficient and lift coefficient of the cylinder at Re = 200.(A) Drag coefficien cylinder; (B) lift coefficient of the cylinder.

Figure 6 .
Figure 6.Velocity profile of flow around three cylinders.

Figure 5 .
Figure 5. Drag coefficient and lift coefficient of the cylinder at Re = 200.(A) Drag coefficient of the cylinder; (B) lift coefficient of the cylinder.

Sustainability 2023 ,
15,  x FOR PEER REVIEW 8 of 14 drag force is twice as much as the lift force due to the influence of the Carmen vortices.The lift coefficients ∆C are 1.70 and 1.67 for the manuscript and the reference, respectively.The calculations are in agreement with the article by Chen[38].

Figure 5 .
Figure 5. Drag coefficient and lift coefficient of the cylinder at Re = 200.(A) Drag coefficient of the cylinder; (B) lift coefficient of the cylinder.

Figure 6 .
Figure 6.Velocity profile of flow around three cylinders.Figure 6. Velocity profile of flow around three cylinders.

Figure 6 .
Figure 6.Velocity profile of flow around three cylinders.Figure 6. Velocity profile of flow around three cylinders.

Figure 7 .
Figure 7.The vorticity of flow around three cylinders.

Figure 7 .
Figure 7.The vorticity of flow around three cylinders.

Figure 8 .
Figure 8.Time domain response of the force coefficient of the flow around three cylinders.(A) Drag coefficient; (B) lift coefficient.

Figure 9 .
Figure 9. Frequency domain response of the force coefficient of the flow around three cylinders.(A) Drag coefficient; (B) lift coefficient.

Figure 8 .
Figure 8.Time domain response of the force coefficient of the flow around three cylinders.(A) Drag coefficient; (B) lift coefficient.

Figure 8 .
Figure 8.Time domain response of the force coefficient of the flow around three cylinders.(A) Drag coefficient; (B) lift coefficient.

Figure 9 .
Figure 9. Frequency domain response of the force coefficient of the flow around three cylinders.(A) Drag coefficient; (B) lift coefficient.

Figure 9 .
Figure 9. Frequency domain response of the force coefficient of the flow around three cylinders.(A) Drag coefficient; (B) lift coefficient.

Table 1 .
Physical properties of LBE.

Table 1 .
[37]flow pattern is mainly affected by the Reynolds number.Within a specific range of Reynolds numbers, different flow patterns can be distinguished by Reynolds numbers.Animasaun et al.[37]disclosed that laminar flow patterns degrade into turbulent eddies and vortices at high Reynolds numbers due to the dominance of nonlinear interactions between fluid particles.A wide variety of fluid systems, from airflow over an airplane wing to water movement in rivers and seas, exhibit this transition to turbulence as a universal feature.The physical properties of LBE can be obtained from Table1, as the temperature of LBE in the inlet of the CiADS heat exchanger is 380 K.The key non-dimensional parameter Reynolds number is 43936.55.The unit of physical quantities in the LBM model is lu (lattice unit) after the dimensionless process.The time can be obtained by  = υ /υ (D/D ) in which parameters without prime are practical quantities, while parameters with prime are lattice non-dimensional parameters.Physical properties of LBE.