Analyzing the influence of an uncertain geometry on fluid-structure interactions in a tube bundle

fluid-structure interactions in a tube bundle

Analyzing the influence of an uncertain geometry on

Academic year 2019-2020

Master of Science in Electromechanical Engineering

Master's dissertation submitted in order to obtain the academic degree of

Counsellor: Ir. Henri Dolfen

Supervisor: Prof. dr. ir. Joris Degroote

Student number: 01306071

Acknowledgements

I would like to thank my supervisor, prof. dr. ir Joris Degroote for all the help he has provided me. His tips and insights greatly improved this thesis work. More importantly, when I was sometimes struggling to find the right path forwards, he was always available to put me on the right track again.

Also a thank you in particular to my counsellor, ir. Henri Dolfen, whose assistance was invaluable during this thesis year. He helped me greatly, from the searching of source material, performing the simulations to the writing of the text itself. He was always available and punctual in his agreements, something which I regret to say was not always the case from my part.

Furthermore, special thanks to C´edric Bruynsteen who reviewed the thesis text, which undoubt-edly greatly improved its readability. Last but not least I would like to thank my family and friends, whose patience I have tested multiple times during this thesis year, but who never stopped supporting me.

Permission to lease-lend

The author gives permission to make this master dissertation available for consultation and to copy parts of this master dissertation for personal use. In the case of any other use, the copyright terms have to be respected, in particular with regard to the obligation to state expressly the source when quoting results from this master dissertation.

Basile Lievens

Ghent University (UGent) June 15th, 2020

Abstract

Tube bundles subjected to axial flow are commonly found in a number of applications, most notably inside the core of nuclear reactors. Although nuclear systems are subject to high levels of regulations and precision engineered systems are employed, the harsh environment is extremely demanding on the equipment. Over long periods of time, large fluid forces, high temperatures and radiation can lead to physical deformations of the fuel rods which cannot always be ade-quately modelled or predicted. The goal of this thesis is to efficiently quantify the influence of the uncertainty in different tube deformations on the hydraulic performance of a tube bundle subjected to steady axial flow. To this end polynomial chaos methods are explored, which treat the uncertainty quantification of the deformation parameters. The tube bundle is tested using Computational Fluid Dynamics (CFD) and a statistical analysis is carried out to assess the in-fluence of these uncertainties on the performance. In a first stage the deformed tube is assumed rigid, only focusing on the effects of the initial deformation. In a second stage Fluid-Structure Interaction (FSI) simulations are performed to see whether the flow itself causes a significant change in the deformation and study the consequences thereof.

Keywords

Polynomial Chaos, Fluid-Structure Interactions, Computational Fluid Dynamics, Uncertainty Quantification, Tube Bundle

Analyzing the influence of an uncertain geometry on

fluid-structure interactions in a tube bundle

Basile Lievens

Supervisors: Prof. dr. ir. J. Degroote, ir. H. Dolfen

Abstract— Tube bundles subjected to axial flow are commonly found in a number of applications, most notably inside the core of nuclear reactors. Although nuclear systems are subject to high levels of regulations and preci-sion engineered systems are employed, the harsh environment is extremely demanding on the equipment. Over long periods of time, large fluid forces, high temperatures and radiation can lead to physical deformations of the fuel rods which cannot always be adequately modelled or predicted. The goal of this thesis is to efficiently quantify the influence of the uncertainty in different tube deformations on the hydraulic performance of a tube bun-dle subjected to steady axial flow. To this end polynomial chaos methods are explored, which treat the uncertainty quantification of the deformation parameters. The tube bundle is tested using Computational Fluid Dynam-ics (CFD) and a statistical analysis is carried out to assess the influence of these uncertainties on the performance. In a first stage the deformed tube is assumed rigid, only focusing on the effects of the initial deformation. In a second stage Fluid-Structure Interaction (FSI) simulations are performed to see whether the flow itself causes a significant change in the deformation and study the consequences thereof.

Keywords— Polynomial Chaos, Fluid-Structure Interactions, Computa-tional Fluid Dynamics, Uncertainty Quantification, Tube Bundle

I. INTRODUCTION

Tube bundles can be found in countless applications, both in the field of thermodynamics and beyond. These kind of geometries have time and again proven to be an efficient and economical way to transfer heat from one fluid to another. In most appli-cations the tube bundle is placed in a cross-flow configuration, which is hence also the most researched. A specific application for axial flow tube bundles is found inside nuclear reactor cores. The tube bundle consists out of fuel rods that are arranged in a hexagonal or square lattice. It is simultaneously subjected to high temperatures, pressures and levels of radiation. This puts severe stresses on all equipment inside the core, but in particu-lar on these rods due to its slender shape. It is not uncommon for these tube bundle geometries to change over their lifespan. This can cause problems in some cases, putting additional forces on the rods and holding structures, increasing the pressure drop over the bundle or even blocking entire fluid channels if the de-formation is too large. The exact modelling of the reaction of the fuel rods on the environment conditions is challenging and not all uncertainties can be eliminated. Proper ways to deal with these uncertainties, also called Uncertainty Quantification (UQ), need to be explored.

The difficulty of access to these kind of environments makes ex-perimental testing also difficult, warranting an approach using CFD, which does not suffer from these limitations. The combi-nation of UQ and CFD methods is an active topic of research as the computational power needed to perform these calculations was not widely available until quite recently.

II. GOALSTATEMENT

The goal of this thesis is threefold. The first goal is to vali-date the proposed polynomial chaos (PC) method (see section III-B) for CFD simulation. To this end, a CFD model will be constructed where a tube inside a tube bundle is subjected to a simple deformation. The PC model will be validated with a Monte Carlo (MC) simulation. The second goal is to use the validated PC method to study other types of bundle deforma-tions that are closer to those observed actual reactors [6]. The influence of the input uncertainty on the uncertainty of certain key output variables, will be analysed. Those variables are: the normal force magnitude Fnand directionθ, and the drag force

Fz on the deformed inner tube. Fz is also directly linked to the

pressure drop∆p over the bundle. A qualitative description of the flow field phenomena will also be given.

In the last stage a FSI simulation will be set up with the same types of deformations as studied in the 1st and 2nd stage, to study the possibility of interactions that may occur between the fluid flow and the deformed tube.

III. UNCERTAINTY QUANTIFICATION

This section aims to give a theoretical overview and explanation of the different methods, relevant to this thesis, that can handle uncertainty in a group of random variables. The goal is to trans-form the continuous stochastic nature of the input variables in such a way that they can be handled by the mathematical com-putation model of the studied system. The resulting output is then analysed to obtain its probabilistic properties.

A. Monte Carlo

One of the most straightforward and consequently one of the most widely used techniques for UQ are (MC) methods [7]. Here the random input variables are sampled in a pseudo-random way to obtain the discrete sample set. Afterwards, each sample is evaluated by solving the computational model and the probabilistic properties from the output are retrieved through statistical analysis. For instance, the true meanµyand variance

σ2

y of an output variable y with any distribution can be estimated

by respectively the sample mean Y and sample variance S2if that sample set is sufficiently large

b µy=Y =_N1 N

∑

∑

where N is the total number of samples and Yi represents one

The main advantages of MC are its simple implementation and its robustness. MC also allows the system itself to be treated as a black box as the uncertain input is converted into a dis-crete set of deterministic samples. This is interesting in the case of CFD simulations as specialized deterministic solvers can be used without the need for adaptation. The main disadvantage is the slow convergence rate, requiring a large number of samples before a conclusive result can be obtained. Different mecha-nisms exist to accelerate the convergence rate to some extent [8], but faster methods are still desired, especially considering the time-consuming nature of CFD in itself.

B. Polynomial Chaos Expansions

Another technique for UQ is the use of Polynomial Chaos (PC) theory. Here, instead of sampling the input as with MC, a so-lution is explicitly constructed by a decomposition of the model using spectral expansions. The spectral expansion of a random variable rests on the principle of separating the deterministic and stochastic parts [7]. Suppose we have a random (output) variable y(xxx,ξξξ). The following then holds

y(xxx,ξξξ) =

∑

∑

where xxx(x1,x2, ...)is the vector containing all deterministic

vari-ables andξξξ(ξ1, ...,ξn)the random variable vector, with a known

distribution. yi(xxx) are deterministic coefficients andΨi(ξξξ)

rep-resent the set of polynomial basis functions which are orthogo-nal to each other w.r.t a weight function. The inner product in the domainΩ of ξξξ is defined as

Ψi(ξξξ),Ψj(ξξξ)=

Z

ΩΨi(ξξξ)Ψj(ξξξ)p(ξξξ)dξξξ (3)

where p(ξξξ) represents the weight function. The choice of the basis functionsΨi depends solely on the distribution ofξξξ.

For example, for a normal (multivariate) distribution orthogonal Hermite polynomials are used as a basis and for a uniform distri-butions Legendre polynomials are more suited [9]. For the spec-tral expansion in Eq.2 to be an exact representation of y(xxx,ξξξ) the number of expansion coefficients P should go to infinity but in practice the series is truncated.

With the aid of Eq.2 and knowledge of the basis functionsΨi(ξξξ)

the problem of finding the probabilistic properties of the output y(xxx,,,ξξξ) can be solved by looking at the coefficients yi(xxx). The

estimated mean bµyand variance bσy2can directly be derived from

these coefficients [5]. Two general approaches can be taken to retrieve yi(xxx): intrusive and non-intrusive. Non-intrusive

meth-ods have the benefit that specialized deterministic solvers can be used for the system calculation, as was the case with MC. Different non-intrusive schemes exist, like point collocation or spectral-projection with quadrature [7], in this thesis the spectral projection method will be used.

With spectral projection methods, the output y(xxx,ξξξ) is projected on the basis functionΨi(ξξξ) using the inner product (Eq.3).

Uti-lizing the orthogonality condition and following from definition, it can be shown that this projection is equal to the respective

de-terministic coefficient yi(xxx) apart for the factor hΨ2i(ξξξ)i

yi(xxx) = y(xxx,ξξξ),Ψi(ξξξ) hΨ2i(ξξξ)i for i = 0,1,...,P = 1 hΨ2i(ξξξ)i Z Ωy(xxx,ξξξ)Ψi(ξξξ)p(ξξξ)dξξξ (4)

Ideally the resulting integral is solved analytically, but in prac-tice this is impossible in almost all cases and a numerical ap-proximation is used. This can be achieved by using quadrature methods where a set of N samples ofξξξ is taken. It can be shown that the smallest number of samples needed to adequately cap-ture Eq.4 is linked to the order of chaos p and the number of random variables n by N = (p + 1)n_{[3]. These samples are}

de-noted byξξξ_j. The integral is then rewritten as yi(xxx) ≈ 1 hΨ2i(ξξξ)i N−1

∑

where the value of the sample set (ξξξ0, ...,ξξξj, ...,ξξξN−1)and the

corresponding weights wjare determined by the specific

numer-ical quadrature method that is employed. The amount of ran-dom variables as well as their distribution influence the type of quadrature method that is most efficient to use.

To summarise, the PC method works by first choosing a suitable set of basis functionsΨibased on the distribution of the given

input random variablesξξξ. Next the factors hΨ2

i(ξξξ)i can be

cal-culated using Eq.3 and a sample set of size N fromξξξ is con-structed with the chosen numerical quadrature scheme. These samples are put through the deterministic solver to obtain N out-put samples y(xxx,ξξξ_j). Then the P + 1 deterministic coefficients yi(xxx) can be determined with Eq.5, to finally be able to obtain

the probabilistic properties of the output.

The big advantage of PC methods is its speed, for a low order of chaos p (usually p ≤ 5) a good result can be obtained. For a rea-sonable amount of random variables n the resulting sample set size for PC can be one or more orders of magnitude smaller com-pared to MC. The downsides are that these methods are more cumbersome to implement and they are not as robust as MC methods, as both the choice of polynomial basis functions and the specific numerical quadrature scheme can influence the rate of convergence quite heavily. Nevertheless PC methods show great promise in CFD simulations, where all gains in computa-tional time are crucial.

IV. NON-FSI SIMULATIONS

All simulation cases performed without FSI are, despite using different UQ schemes or being subjected to different types of deformation, structured in a similar way. This structure is shown in Fig. 1 and every step will be explained in more detail below. A. CFD model

The first step is the construction of a CFD model with no defor-mations that will be used as base a for every subsequent case. Both the construction of the mesh and the CFD calculations are performed inside the ANSYS software package. The tube

Calculate N samples from the stochas-tic input variables applying UQ method

i = 1

Deform nodes of the CFD mesh according to sample i Perform steady CFD calculation until convergence

Extract output variables for sam-ple i: pressure p and forces ~Fn, ~Fz

i = i + 1

i ≤ N?

Calculate distribution (µ, σ) of the out-put variables applying UQ method

no yes

Fig. 1: Calculation Flowchart

bundle geometry is that of a hexagonally stacked bundle and is derived from the geometry described in [1]. This geometry is in turn is based on the design of MYRRHA, a prototype nuclear re-actor which is being developed at the SCK-CEN. To reduce the required computational power, a simplified mesh is constructed which consists of one inner tube that can be deformed and walls sections of the six surrounding tubes in the bundle lattice, which will be kept stationary. An overview of the geometry and mesh-ing parameters can be found in Table 1.

Notice the very slender shape of the tubes, considering the large L/D ratio. This makes the tubes more susceptible to deforma-tions and vibration in the transversal direction, perpendicular to the tube length axis. In practice that transversal movement is limited by spacer grids or wires, but modelling these falls out-side the scope of this thesis. Also note that that a Pb-Bi eutectic coolant (MYRRHA design) is used, which has fundamentally different properties compared to more traditional coolants like water. One interesting property to highlight is the density, which is about ten times that of water. This has two direct effects on the fluid dynamics: any vibrations that might occur will be damped more with Pb-Bi flow and any forces that arise from deflection of the flow due to tube deformations will be larger (if all other flow parameters remain the same).

The constructed mesh has a total of 441 000 cells and can be seen in Fig. 2. A mesh sensitivity study was performed to en-sure that the quality of the base mesh is sufficient, before any de-formation has occurred. A high Reynolds mesh was constructed (y+_{'= 30), so that the first grid point from the tube wall falls}

into the logarithmic layer. The driving factor for this decision is the reduced computational cost. A low Re mesh (y+_{=0.8) was}

also constructed for the sensitivity study as comparison and the results remain similar between the two. The deformation of the mesh is realised by manually calculating the new positions of the nodes on the inner tube wall according to the specific defor-mation function. The rest of the mesh domain nodes are adapted

with an automatic mesh smoothing function. A diffusion based smoothing method was used as it was observed that this method could maintain the mesh quality the best.

As CFD solver, ANSYS Fluent is used, which is a finite vol-ume solver. A quick summary of the relevant settings for Flu-ent is given here. For more information on the rationale behind these settings, the reader is referred to the paper of J. De Rid-der et al. [1]. The studied flow is modelled as incompressible and steady. The turbulence is modelled with the Reynolds Av-eraged Navier-Stokes (RANS) equations supplemented with the k − ω SST model, which is a good combination of the k − ω model (near-wall) and the k − ε model (free-stream). The inlet is subjected to a uniform velocity of 1.5 m/s which results in a Re number ' 60 000. The turbulence at the inlet is specified by a turbulent intensity of 5% and characteristic length scale of 0.1 mm. The outlet is modelled as a constant pressure outlet. No gravity force acts on the fluid domain. The walls of the inner and outer tubes are zero slip boundaries and the gap openings be-tween the outer tubes are modelled as periodic boundaries. The flow passing through one such opening is thus passed through to the opening that is diametrically opposed to it, creating the effect of multiple (identically) deforming tubes simultaneously. This behaviour has been observed in reality ([6]. Another op-tion was to make these boundaries symmetric, creating a very confined flow domain, representing a tube bundle with one de-forming tube isolated from the rest. However this situation is almost never encountered in real-life tube bundles and is not further investigated.

B. Deformations and UQ

Three different cases will be studied: a simple deformation where a straight tube is put at an angle (‘Angle tube’), and two more realistic deformations where the ends of the tube are clamped and the tube is deformed according to the first (‘C-shape’) and second (‘S-(‘C-shape’) mode shape of a transverse oscil-lating beam. Each deformation is characterised by two random variables, a variable describing the plane in which the deforma-tion is occurring: the direcdeforma-tion angleφ, and a variable describing the magnitude of the deformation (Fig. 3). The magnitude pa-rameter in the Angle-Tube case is defined as the angleψ that the tube centreline makes with the z-axis. The magnitude parameter in the C-shape and S-shape configuration is defined as the max-imum deflection u0of the centreline with respect to the z-axis.

For the UQ methods, a suitable distribution needs to be given for the random variables. All random parameters are given the same uniform distribution for simplicity. Other distributions or com-binations of distributions for the different parameters are also a possibility but fall outside the scope of the current investiga-tion. Legendre polynomials are used as basis functions. The distribution boundaries for the three cases are given in Table 2. The limits of the directionφ is chosen in such a way that the symmetry of the geometry is utilised to a maximum. The mag-nitude parameterψ/u0is limited by the bundle geometry, as for

larger magnitudes, the fluid domain is deformed more and more, which deteriorates the mesh quality, eventually reaching a point where the fluid solver breaks down. The practical

implemen-Fig. 2: View of the base CFD model mesh (blue) annotated with the most important boundary conditions. The CSM mesh of the inner tube (grey) for the FSI simulation (Section V) is also shown.

tation of the UQ methods is done with ChaosPy [4], a Python based statistical software package.

C. Monte Carlo vs. Polynomial Chaos

The first simulation run aims to validate the proposed PC method. To this end two cases are constructed, each with an Angle Tube deformation with deformation parameters described above. The first case is a Monte Carlo simulation with 10 000 samples. The amount of samples is high (a sample set of N =2 000 also achieves good results for this case) but was done to ensure that any error due to the random sampling is sufficiently low and a good baseline for comparison is created. The sec-ond case is the PC simulation, with order of chaos p = 4 and a Gaussian-Legendre quadrature scheme of the same order, re-sulting in a sample set size of N + 1 = 25. The output variable under consideration is the resultant force on the inner cylinder, separated into a normal component (xy-plane), with magnitude Fnand a directionθ, and an axial component Fzalong the z-axis

(Fig. 3). After the statistical validation, the results of the case itself will be discussed, focusing on the correlations between in-put and outin-put variables. With the cases completely constructed, the total calculation time per sample for each case is ca. 150 s, resulting in a total calculation time of 1.05 h for PC, 3.5 days for MC (N =2 000) and 17 days for MC (N =10 000).

The mean and variance of the output variables for both cases is shown in Table 3. a 95%-confidence interval on the MC so-lution was calculated with a bootstrap method. The data shows a good convergence of the PC solution to the values of the MC for all output variables, with one exception. The mean value of θ falls just outside the 95%-confidence interval. Additional PC simulations were performed with both lower and high order of chaos. All these results did fall into the interval of the MC solu-tion. This leads to the conclusion that the exception was just a one time error and that in general a PC simulation of order p = 4 works well for evaluating the uncertainty of the input. However, the exception does show that these methods are not that robust and can depend significantly on the different settings used. Fig.

Geometrical Data Fluid Properties

Parameter Value Parameter Value

Diameter D 6.55 mm Type Pb-Bi

Length L 700 mm Temperature T 608 K Pitch P 8.384 mm Densityρ 10 291 kg/m3

P/D 1.28 Dyn. visc.µ 0.0017 Pa.s

Table 1: Geometrical Data and Fluid Properties of the CFD Model.

(a) (b)

(c) (d)

Fig. 3: Drawing of the geometry with indicated input and output variables for an Angle Tube (b), C-shape (c) and S-shape (d) deformation. The bundle L/D ratio is reduced and deformations are amplified for a better visualization.

Case φmin φmax ψmin/u0,min ψmax/u0,max

Angle Tube 0◦ ₃₀◦ ₀◦ _0.127◦

C-shape 0◦ ₃₀◦ _{0 mm 1.45 mm}

S-shape 0◦ ₃₀◦ _{0 mm 1.45 mm} Table 2: Upper and lower limits for the uniformly distributed parametersφ and ψ/u0of all three deformation cases.

4 further substantiates the performance of the PC method, where both the PC and MC histograms show a large similarity, which is of course desired.

Furthermore, it was observed that the uniform distribution is largely maintained for the normal force magnitude Fn and

di-rectionθ (not shown) but not for the drag force Fz(Fig. 4). This

suggests that the correlation between the input and Fn,θ is of

a linear form. A more in-depth discussion of these correlations relations between the input and output is given in the next sec-tion for the C- and S-shape deformasec-tion, as these correlasec-tions are similar for all three deformation types.

D. C-shape and S-shape with Polynomial Chaos

The second simulation run again consists of two cases, one with an C-shape and one with an S-shape deformation. A PC expan-sion of order 4 with 4th order Gaussian-Legendre quadrature is used for UQ, which was validated in the previous section. The main goal here is to study the phenomena that occur in these

Mean PC MC 95%-confidence interval Fn(mN) 6.15460 6.18854 [6.11971 6.25763]

θ (◦_{) -164.933 -164.764 [-164.929 -164.596]}

Fz(mN) 957.799 957.713 [957.454 957.976]

Std. dev. PC MC 95%-confidence interval Fn(mN) 3.54239 3.52782 [3.49730 3.55939]

θ (◦_{) 8.67455 8.65322 [8.57424 8.72920]}

Fz(mN) 13.2982 13.2883 [13.1453 13.4262] Table 3: Mean and std. dev.n of Fn,θ and Fzfor the MC (N =10 000) and PC

(p =4) simulation. 930 940 950 960 970 0 300 600 900 1200 1400

(a) Monte Carlo

930 940 950 960 970 (b) Polynomial Chaos

Fig. 4: Histogram of the drag force Fz(mN) after an Angle Tube deformation

for the MC and PC method.

kind of deformations and to study the influence of the deforma-tion in quesdeforma-tion on the output forces. Fig. 5 shows the output force magnitudes (Fn, Fz) and direction (θ) plotted against

re-spectively u0andφ for both the C- and S-shape deformation.

Both Fn and Fz showed no influence of the deformation

direc-tion φ, hence only one value of φ is presented (Fig. 5a). A linear relationship is observed between Fnand the deformation

magnitude u0for the C-shape, as was also observed for the

An-gle Tube case (not shown). This relationship between Fn and

u0is lost for the S-shape, due to its more complex deformation

geometry. Fig. 5b shows that the normal force is directed in the same sense as the initial deformation, where for the Angle Tube deformation the normal force was 180◦_{opposed toφ. This is}

caused by the curvature of these deformation shapes, which is closely followed by the fluid flow (no separation observed). The fluid turning motion creates a force in the direction of the curved deformation. The Angle Tube deformation is absent of any cur-vature, here the deflection of the flow in one direction logically creates a force in the opposing direction. Note that the S-shape is point-symmetric around z =1₂L, thus the two curvatures can-cel each other out. This explains the small magnitude of Fn. The

normal force that is still being generated is probably caused by smaller, local, deflections of the tube, favouring the directionφ.

(a) C-shape (triangles) / S-shape (squares)

0 5 10 15 20 25 30 φ (◦₎ 0 5 10 15 20 25 30 θ ( ◦) θ = φ u0= 0.07 mm u0= 0.73 mm u0= 1.12 mm u0= 1.38 mm (b) S-shape

Fig. 5: Summary of the output forces. Magnitude of Fzand Fnin function of the

max. deformation u0is shown in (a) (φ = 15◦). The normal force directionθ in

function of the deformation direction angleφ, for different u0, is given in (b).

The effect of the magnitude u0on this relationship is

insignif-icant, except in the case of large deformations where a small deviation is observed. In these cases, the extremities of the in-ner tube deformation are nearly touching the surrounding tubes. This causes a displacement of the fluid in the narrowing gap, through the periodic boundaries to the opposing side. This dis-placement happens symmetrically when the deformation is lined up with a symmetry plane of the bundle geometry (φ = 0◦_,

φ = 30◦_{), however when φ is not aligned with such a plane,}

the fluid displacement is asymmetrical. The hypothesis is that this asymmetry creates an extra force on the inner tube directed away from the initial deformation direction, which causes the deviation observed in Fig. 5b. Further investigation is needed for a rigorous proof, but this falls outside the scope of this study. Fig. 6 show the planar velocity for variousφ, as clarification. The drag force magnitude is shown to decrease with increasing deformation magnitude (Fig. 5a). At first sight, this is counter-intuitive as a deflecting object usually causes more drag to be generated. To further comprehend this result, the drag force was divided into the drag generated by pressure forces and the part generated by viscous forces. It was found that the pressure drag does indeed increase with increasing u0, in the same order as

the normal force. However, the viscous drag decreases in a non-linear fashion. Due to the long and narrow channels, created by the bundle geometry, the viscous drag is completely dominant. The non-linear relationship between Fz and u also explains the

skewed histogram (Fig. 4).

V. FSISIMULATIONS

A. Setup

The possibility exists that significant interactions between the flow field and the deformations occur. To investigate this, a FSI simulation is set up. An FSI problem can be deconstructed into two smaller problems: a fluid-dynamical part (fluid domain) and a mechanical part (the deforming inner tube). These two prob-lems can be solved simultaneously with a monolithic approach, but here a partitioned approach is taken. This allows us to use a specialised solver for each part, analogue to the reasoning be-hind the non-intrusive UQ methods (Section III-B). The two solvers used are ANSYS Fluent for the CFD part, as before, and

C-shape S-shape

µ σ µ σ

Fn(mN) 0.50896 0.31386 0.13459 0.059232

θ (◦_{) 14.8775 8.68322 14.8742 8.95062}

Fz(mN) 958.236 12.8359 963.345 8.76256 Table 4: Mean and standard deviation of Fn,θ and Fzfor an C-shape and S-shape

deformation.

(a)φ = 0◦ _{(b)φ = 15}◦ _{(c)φ = 30}◦

ABAQUS for the Computational Structural Mechanics (CSM) simulations. A two-way coupling scheme is used to link these two solvers together. The coupling method used is IQN-ILS, of which more information can be found in [2]. This is a strongly coupled scheme, meaning that the solver solutions are called multiple times per time step. This is necessary as only steady simulations are performed and thus only one “time step” is cal-culated. The practical implementation of the complete FSI is performed with CoCoNut, a Python software package developed by the STFES research group at Ghent University.

The steel inner tube is modelled as a hollow cylinder with outer diameter D = 6.55mm and inner diameter Di=5.65mm, made

out of steel with densityρ = 7850kg/m3_{and Elasticity}

modu-lus E = 183GPa. The CSM mesh consists of 11 760 quadratic elements. For the CFD part the same settings were used as de-scribed in Section IV-A. For the UQ again the same PC model (p = 4) was used as in the previous parts, with a total sample set size of 25. The three deformation cases with parameters de-scribed in Table 2 were used for the FSI. The calculation time for all samples of one such case was approximately 10 hours. B. Results

A comparison between the FSI and non-FSI simulations for the mean and standard deviation of the output variables Fn,θ and Fz

is shown in Table 5 for the three cases. One can immediately see that the difference between the respective simulations is small, in all cases. The reason is that the displacements are quite low, resulting in only a small interaction between fluid and structure. In this sub-critical regime, no dynamic phenomena are expected (e.g. flutter), justifying the use of steady simulations. Other steady flow phenomena (e.g. buckling) could have appeared, resulting in a permanent deformation of the tube bundle. This would have been visible in the data by a large difference be-tween FSI and non-FSI calculations, but this was not the case.

Angle tube C-shape S-shape

µ σ µ σ µ σ

Fn +.05% +.11% +.31% +.24% +1.6% +1.5%

θ +.00% -.07% -.16% -.02% -.14% -.37% Fz +.00% -.02% +.00% +.22% +.00% -.03% Table 5: Relative difference of the mean and std. dev. of Fn,θ, Fzbetween the

FSI case and the non-FSI counterpart, given for the three deformation types.

However that does not mean that no deformations occur which were caused by the fluid flow. Figure 7 shows the mean addi-tional deflection∆u of the inner tube centreline compared to the initial deformation, for an Angle Tube and a C-shape deforma-tion (S-shape not shown). The amplitudes are small but observe the bowing of the initially straight tube due to the incoming flow that hits the tube. The C-shape deformation is less intuitive but it is observed that the areas subjected to a large initial deforma-tion, are deformed even more. This creates a more ‘accentuated’ geometry which results in more tension inside the tube.

VI. CONCLUSION

Three major simulation runs were performed and as such three conclusions are drawn. Firstly it was observed that the proposed

0 175 350 525 700 z-co (mm) −1.6 −1.2 −0.8 −0.4 0.0 ∆u (µ m)

(a) Angle Tube

0 175 350 525 700 z-co (mm) 0.0 0.3 0.6 0.9 1.2 1.5 ∆u (µ m) (b) C-shape

Fig. 7: Mean (blue) of the additional centreline deflection∆u in case of FSI simulations, shown with plus and minus one std. dev. (dashed, black).

PC method converges adequately to the correct solution, which was validated with a pseudo-random sampled MC simulation. The results of the PC method fell within the 95%-confidence interval of the MC simulation. The selection of the specific quadrature scheme for PC methods is highly important, as these methods are in general less robust.

Second, the PC method was used to study the effects of uncer-tainty in the input on two more complex deformation cases: a ‘C-shape’ and a ‘S-shape’. The obtained results showed some expected effects, such as the proportionality between deforma-tion magnitude u0and magnitude of the normal force, as well as

its direction being in-line with the deformation direction. How-ever it was also observed that as the deformation becomes more complex, these properties are not completely maintained. For the simulated S-shape, the proportionality between the defor-mation and normal force magnitude disappears. Also an effect of u0on the normal force directionθ was found, which is not

present in the other cases. The drag force was found to decrease with increasing deformation magnitude in all cases.

Lastly an FSI simulation was set up for all deformation cases. The resulting output variable statistics showed little to no differ-ence compared to their non-FSI counterparts, indicating that for the current flow regime and studied geometries no significant in-teractions between flow and structure occur. The magnitude of the additional centreline deflection by the FSI was studied. All cases showed an additional bowing caused by the fluid flow, in the same direction as the normal force acting on the tube.

REFERENCES

[1] Jeroen De Ridder, Joris Degroote, Katrien Van Tichelen, Paul Schuurmans, and Jan Vierendeels, Modal characteristics of a flexible cylinder in turbulent axial flow from numerical simulations, Journal of Fluids and Structures43 (2013), 110–123. [2] Joris Degroote, Klaus J¨urgen Bathe, and Jan Vierendeels, Performance of a new

parti-tioned procedure versus a monolithic procedure in fluid-structure interaction, Comput-ers and Structures87 (2009), no. 11-12, 793–801.

[3] Michael Eldred, Clayton Webster, and Paul Constantine, Evaluation of Non-Intrusive Approaches for Wiener-Askey Generalized Polynomial ChaosApril (2008). [4] Jonathan Feinberg and Hans Petter Langtangen, Chaospy: An open source tool for

de-signing methods of uncertainty quantification, Journal of Computational Science11 (2015nov), 46–57.

[5] Serhat Hosder, Robert W. Walters, and Michael Balch, Point-collocation nonintrusive polynomial chaos method for stochastic computational fluid dynamics, AIAA Journal 48 (2010dec), no. 12, 2721–2730.

[6] IAEA, Structural Behaviour of Fuel Assemblies for Water Cooled Reactors, Atomic EnergyJuly (2005), 324.

[7] O P Le Maˆıtre and Omar M Knio, Spectral Methods for Uncertainty Quantification: With Applications to Computational Fluid Dynamics, Springer, 2010.

[8] Robert W Walters and Luc Huyse, Uncertainty Analysis for Fluid Mechanics with Ap-plications, Office50 (2002), 50.

[9] Dongbin Xiu and George EM Karniadakis, The Wiener-Askey Polynomial Chaos For Stochastic Differential Equations, Society for Industrial and Applied Mathematics24 (2002), no. 2, 619–644.

Contents

List of Figures xii

List of Tables xvi

Nomenclature xvii

1 Introduction 1

1.1 Goal Statement . . . 2

2 Literature Study 4 2.1 Tube Bundle Deformations . . . 4

2.1.1 Deformations and Failure Mechanisms . . . 4

2.1.2 Consequences of Deformations and Failures . . . 5

2.2 Uncertainty Quantification . . . 7

2.3 Monte Carlo . . . 8

2.4 Polynomial Chaos Expansions . . . 9

2.5 Non-intrusive PC Methods . . . 12

2.5.1 Point-Collocation . . . 12

2.5.2 Quadrature Methods . . . 13

2.5.3 Sparse Grid Methods . . . 15

2.5.4 Example Polynomial Chaos . . . 16

2.6 Numerical Modelling . . . 17

2.6.1 Governing Fluid Equations . . . 18

2.6.2 Finite Volume Method . . . 19

2.6.3 Turbulence . . . 19

2.7 Flow Induced Vibrations . . . 21

2.8 FSI . . . 23

3 Methodology 25 3.1 CFD model . . . 26

CONTENTS 3.1.1 Mesh . . . 26 3.1.2 CFD Solver . . . 29 3.2 Deformations . . . 30 3.3 UQ Setup . . . 31 4 Results non-FSI 33 4.1 Polynomial Chaos vs. Monte Carlo . . . 33

4.1.1 Statistical Analysis . . . 33

4.1.2 Flow Analysis . . . 36

4.1.3 Output Correlations . . . 37

4.2 C-shape and S-shape with Polynomial Chaos . . . 40

4.2.1 Statistical Analysis . . . 41

4.2.2 Output Correlations . . . 43

4.2.3 Flow Analysis S-shape . . . 47

5 Fluid Structure Interactions 51 5.1 Setup . . . 51

5.2 Results . . . 52

5.2.1 Statistics of Fn, θ and Fz . . . 53

5.2.2 Centreline Deformation ∆u . . . 54

5.2.3 Data Visualisation . . . 60

6 Conclusion 62 6.1 Future Work . . . 63

List of Figures

1.1 Cross section of the MYRRHA reactor [1]. . . 2

2.1 Cuttrough of the MYRRHA reactor core [1]. . . 4

2.2 Diagram of a LMFR FA and fuel rod [2]. . . 4

2.3 Examples of void swelling causing fuel rod deformations [2]. . . 5

2.4 Example of fretting defect [3]. . . 5

2.5 Fuel rod failure as a result of corrosion [3]. . . 6

2.6 Observed bow deformation in two PWR reactors at the Ringhals Power Plant [4] 6 2.7 Schematic of a system with an uncertain input, containing both deterministic and random components. The listed example causes of randomness are specific for the tube bundle case. . . 7

2.8 Common sources of uncertainty and error in CFD according to [5] . . . 8

2.9 Left: First five Legendre polynomials LPc(ξ), given for a random variable ξ defined on the interval [−1, 1]. Right: Total number of PC coefficients in function of order of chaos p and the number of random variables n. n is an integer value, ranging from 1 to 5. . . 11

2.10 Total number of samples N for a Gaussian quadrature in function of order of chaos and the number of random variables n. n is an integer value, ranging from 1 to 5. . . 15

2.11 Illustration of the sample set for a two-dimensional (n = 2) random variable vector in the case of a full tensor product (left) and two isotropic Smolyak sparse grid tensorizations (middle: Clenshaw-Curtis abcissas, right: Gaussian abcissas) [6]. . . 16

2.12 Comparison of the mean ¯y(t), with indicated sample set and analytical solution, between a MC (N = 250) method and a PC Gaussian quadrature (p = 4) method. 17 2.13 Turbulent Energy cascade (left) and energy spectrum (right) with indicated mod-elling techniques [7] . . . 20

LIST OF FIGURES

2.14 Sketch indicating the different dynamic zones, with change in eigenfrequency in with increasing flow velocity (left) and maximum amplitude of the cylinder vibra-tions in function of flow velocity (right) [8]. . . 22 2.15 Axial velocity contour in a bundle cross-section [9]. . . 23 2.16 Velocity profile and flow model in the gap region [10]. . . 23 2.17 Schematic of two FSI system coupling schemes. Each system consists of a fluid

solver _{F with input x and output ˜y, and a structural solver S with input y and} output ˜x. [11] . . . 24 3.1 Calculation Flowchart . . . 25 3.2 Mesh sensitivity test results, influence of number of divisions on the pressure

drop ∆p. The number of divisions for the radial sweep was recalculated to the equivalent y+_{-wall value. . . . . 27}

3.3 View of the base CFD model mesh (blue) annotated with the most important boundary conditions. The CSM mesh of the inner tube (grey) for the FSI simu-lation (Sect. 5) is also shown. . . 28 3.4 Comparison of the mesh quality in a tube bundle cross section between a spring

based smoothing and a diffusion based smoothing method for an equal deformation. 28 3.5 Drawing of the tube bundle geometry with indicated axes, input and output

variables for the three deformation cases. The bundle L/D ratio is greatly reduced and deformations are amplified for a better visualization. . . 31 4.1 Mean (left) and standard deviation (right) of Fn, θ and Fzgiven relative to the MC

simulation with indicated 95%-confidence interval and the PC simulations ranging order p = 1→ 5. The confidence interval was constructed with a bootstrap method. 35 4.2 Histograms of Fn(left), θ (middle) and Fz (right) for the PC and MC simulations.

Each histogram contains 10 000 samples, divided in 50 equally sized bins. . . 36 4.3 Velocities and pressures in a cross section of the tube bundle at z = 0.35m = 1

2L,

for the Base Model (φ = 0.0◦_{, ψ = 0.0}◦_{), Angle Tube 1 (φ = 1.4}◦_{, ψ = 0.12}◦_{) and}

Angle Tube 2 (φ = 28.6◦_{, ψ = 0.12}◦_{) case. . . . 38}

4.4 Normal force magnitude Fn in function of deformation angle ψ. . . 39

4.5 Normal force direction θ in function of deformation direction φ. . . 40 4.6 Drag force magnitude Fz and pressure drop ∆p in function of deformation angle ψ. 40

4.7 PC coefficients yi in the 4th-order PC method for Fn, θ and Fz, given relative to

the 0th_{-order term (= mean of the distribution) and separated by corresponding}

polynomial basis function degree px. . . 41

4.8 Histograms of Fn(left), θ (middle) and Fz (right) for the C-shape and S-shape

de-formation cases, simulated with a 4th_{-order PC method. Each histogram contains}

LIST OF FIGURES

4.9 Normal force magnitude Fn, drag force magnitude Fz and pressure drop ∆p in

function of deformation magnitude u0 for an C-shape (blue) and S-shape (red)

deformation. Results are given for multiple deformation directions φ. . . 44 4.10 Sketch of the flow domain (blue) with deforming tube (grey) as seen from the flow

inlet for the reference geometry and different deformation types. . . 45 4.11 Normal force direction θ in function of deformation direction φ for a C-shape and

S-deformation. Results are given for multiple deformation magnitudes u0. Two

data points at φ = 0◦ _{and φ = 30}◦ _{were added for every amplitude u}

0. . . 45

4.12 PC coefficients yi in the 4th-order PC method for Fn, θ and Fz, given relative to

the 0th_{-order term (= mean of the distribution) and separated by corresponding}

polynomial basis function degree px. . . 46

4.13 Velocity vectors and magnitude in the fluid domain of an S-shape deformation with direction φ = 0◦ and magnitude u0 = 1.38 mm. . . 47

4.14 Planar Velocityqv2

x+ v2y in a bundle cross-section at different axial positions for

an S-shape 1 (φ = 0.0◦, u0 = 1.38mm), S-shape 2 (φ = 15.0◦, u0 = 1.38mm) and

C-shape (φ = 15.0◦_{, u}

0 = 1.38mm) case. . . 48

4.15 Decomposition of the normal force Fn into a component parallel to the

deforma-tion plane F_nk, and a perpendicular component Fn⊥. . . 49

4.16 Left: Schematic of circumference angle ϕ along the inner tube, with indicated deformation direction φ. Right: Pressure distribution along the circumference of the inner tube for the Angle Tube (φ = 15◦, ψ = 0.12◦), C-shape and S-shape deformation (φ = 15◦, u0 = 1.38mm), measured at z = 0.35m. Pressure

is presented as the difference between the pressure value and the mean pressure over the complete circumference. . . 50 5.1 Flow domain of an Angle Tube deformation case with indicated additional

deflec-tion due to FSI in blue hatched lines (exaggerated for visualisadeflec-tion). Addideflec-tional deflection of the tube centreline along the z-axis is given by ∆u with a maximum value ∆umax. . . 53

5.2 PC coefficients yi in the 4th-order PC method for Fn, θ and Fz of an S-shape

deformation case, given relative to the 0th_{-order term (= mean of the distribution)}

and separated by corresponding polynomial basis function degree px. . . 55

5.3 Mean µ of additional deflection ∆u (blue), shown with plus and minus the stan-dard deviation σ (black, dashed) along the length of the tube for the three de-formation types. The C-shape and S-shape plots also show a magnification of µ and σ for respectively the start and end portion of the tube. Red-dashed line indicates the location of the maximum (absolute) additional deflection ∆|umax|. . 56

5.4 Raw data of the additional centreline deflection ∆u for an Angle Tube deformation (φ = 1.4◦_{, ψ = 0.006}◦_{), with highlighted fluctuations. . . . 57}

LIST OF FIGURES

5.5 Histograms of |∆umax| in an FSI simulation for the three deformation types

(N =10 000). . . 59 5.6 Maximum additional deflection|∆umax| in function of φ and ψ / u0 for the Angle

Tube (black), C-shape (blue) and S-shape (red) deformations. . . 59 5.7 Additional centreline deflection ∆u, magnified with a factor 100, on top of the

initial deformation for each deformation typeommlk. Each point on the graph represents an average of 5 data points. . . 61

List of Tables

2.1 Different continuous distributions of one random variable ξ with their correspond-ing Wiener-Askey chaos polynomials and weight functions. α, β are probability parameters of the Gamma- and Beta-distribution; a, b∈ R . . . 10 3.1 Geometrical Data and Fluid Properties of the CFD Model. . . 26 3.2 Upper and lower limits for the uniformly distributed parameters φ, ψ and u0 of

all three deformation cases. . . 32 4.1 Summary of the studied MC and PC cases. Case MCb gives an indication of a

more sensible MC calculation time but was never actually calculated. . . 34 4.2 Mean and standard deviation of Fn, θ and Fz for the MC (N =10 000) with

95%-confidence interval (95%-CI) and PC (p =4). . . 34 4.3 Decomposition of the drag force Fz into its pressure and viscous components, for

a low deformation (ψ = 0.03◦) and a high deformation (ψ = 0.12◦) case. . . 40 4.4 Mean µ and standard deviation σ of Fn, θ and Fz for an C-shape and S-shape

deformation. Values for Angle Tube case (Table 4.2) are shown for comparison. . 42 4.5 Normal force decomposition into Fnk and Fn⊥ for input φ = 15◦, u0 = 1.38 mm

(θexp= φ). . . 49

5.1 Dimensions and material properties of the inner tube subjected to deformation. . 52 5.2 Mean and standard deviation of Fn, θ, Fz for the FSI simulations with

deforma-tions: Angle Tube, C-shape and S-shape. . . 54 5.3 Relative difference of the mean and standard deviation of Fn, θ, Fz between the

Nomenclature

Symbols

A [m2_{] cross-sectional area as seen by the fow}

D [mm] outer diameter

Di [mm] inner diameter

E [GPa] Young’s modulus EI [Pa m4_{] Flexural rigidity}

Fn [mN] magnitude normal force

Fnk [mN] normal force component parallel with deformation direction

Fn⊥ [mN] normal force component perpendicular to deformation direction

Fz [mN] magnitude drag force

g [m s−2] gravity vector

k [J] turbulent kinetic energy

L [mm] length tube

Lt [mm] turbulent length scale

LPc Legendre polynomials n [-] number of random variables np [-] oversampling ratio

N [-] sample set size

p [-] order of chaos

px collection of polynomials with degree x

P [mm] Pitch

Pc [-] number of PC coefficients

Re [-] Reynolds number

¯

S2 _{sample variance}

S0 [s−1] shear deformation rate tensor

Sm kg m−2s−2 shear stress source term (tensor)

t [s] time

NOMENCLATURE

u0 [mm] deformation magnitude of ‘C-shape’/‘S-shape’

huiuji [m2s−2] Reynolds stresses

U [m s−1_{] inlet velocity}

v [m s−1] velocity vector

vdim [-] dimensionless flow velocity

q v2

x+ vy2 [m s−1] planar velocity

x deterministic variable vector y(x, ξ) system output

yi polynomial chaos coefficients

¯

Y sample mean

Yi Monte Carlo sample

x, y, z Cartesian coordinates of the axis system

δij Kronecker delta

∆p [kPa] pressure drop over tube bundle

∆u [mm] additional tube centreline deflection due to FSI

∆umax [mm] maximum additional tube centreline deflection due to FSI

 [J kg−1s−1] turbulent dissipation rate

λk [-] √ωk

lDI [m] smallest Kolmogorov length scales

µ [Pa s] dynamic viscosity

µy mean

c

µy expected value of the mean

ν [-] Poisson’s ratio

ωk [-] non-dimensional natural frequency

φ [°] deformation direction

ψ [°] deformation magnitude of ‘Angle Tube’ case Ψ polynomial basis function

ω [kg−1s−1] turbulent frequency ρ [kg m−3_{] density fluid}

ρs [kg m−3] density structural material

σy standard deviation

c

σy expected value of the standard deviation

τ [Pa] shear stress tensor θ [°] direction normal force ξ stochastic variable vector

NOMENCLATURE

Subscripts

i, j, k Indices of Polynomial Chaos Expansion coefficient F SI indicates variable obtained by FSI calculation noF SI indicates variable obtained by non-FSI calculation rel relative difference between two variables

Acronyms

CANDU CANada Deuterium Uranium CFD Computational Fluid Dynamics

CI Confidence Interval

CSM Computational Structural Mechanics DNS Direct Numerical Simulation

FA Fuel Assembly

FIV Flow-Induced Vibrations

FVM Finite Volume Method

IQNI-ILS Interface Quasi-Newton Iterations with Inverse jacobian Least-Squares LES Large Eddie Simulation

LBE Lead-Bismuth Eutectic

LMFR Liquid Metal Fast Reactor

MC Monte Carlo

MYRRHA Multi-purpose hYbrid Research Reactor for High- tech Applications ODE Ordinary Differential Equation

PC Polynomial Chaos

PDF Probability Density Function PWR Pressurized Water Reactor RANS Reynolds Averaged Navier-Stokes SST Shear-Stress Transport

STFES Sustainable Thermo-Fluid Energy Systems (UGent) UQ Uncertainty Quantification

Chapter 1

Introduction

Tube bundles can be found in countless applications, both in the field of thermodynamics and beyond. These kind of geometries have time and again proven to be an efficient and economical way to transfer heat from one fluid to another. One specific application of tube bundles is inside nuclear reactor cores. Here a coolant fluid, typically water, flows axially along a bundle of fuel rods. The heat absorbed from the nuclear fission reaction can then be used for e.g. electricity generation.

A new type of nuclear reactor currently being developed at the SCK-CEN (Mol - Belgium), is MYRRHA (Multi-purpose hYbrid Research Reactor for High- tech Applications), which is shown in Fig. 1.1. This reactor is a Gen.IV reactor design intended for nuclear research. Gen.IV is a collective initiative between a host of different countries with a stake in nuclear energy. The aim is to research and develop new reactor designs, for commercial applications, based on these key features: highly economical and reliable, with improved safety and less radioactive waste. MYRRHA differs from a conventional nuclear reactor in multiple ways.

It is designed as an Liquid Metal Fast Reactor (LMFR), which means that the coolant flowing between the fuel rods is liquid metal (lead bismuth eutectic - LBE) instead of conventional water. This has several benefits, such as e.g. a higher boiling point (_{≈ 1670}◦_{C) enabling the reactor}

to be run at lower pressures, improving safety. It also imposes several challenges, requiring a sturdier support structure as LBE is ten times denser than water. There is also the increased possibility of corrosion due to the increased reactivity between LBE and the fuel rod material. These challenges combined with the harsh operating conditions, generating large amounts of radiation and heat, put serious strain on the structural components of the system. The fuel rods in particular are severely impacted, as these hollow tubes are long, slender and packed very close together. A certain amount of deformation over the lifespan of these rods is to be expected. Making an exact prediction of where and when these deformations will occur, is very difficult.

Chapter 1: Introduction

Figure 1.1: Cross section of the MYRRHA reactor [1].

Hence it is necessary to take into account a certain degree of uncertainty in the geometry. The modelling of uncertainty, also known as uncertainty quantification (UQ), is already well established and mathematical models have been around for an extensive period of time. However, the combination of UQ and Computational Fluid Dynamic (CFD) simulations is still a very active topic of research. The combination of these two techniques typically requires a lot of computational power, which was not widely available until quite recently.

The nuclear industry is highly regulated due to the inherent risks associated with nuclear fission. Furthermore, the construction of experimental setups are costly and time consuming. This makes the application of CFD for the design of nuclear reactor cores very attractive. CFD simulations can help gain insight into the nuclear reactor performance without the need for expensive experimental setups.

1.1

Goal Statement

The goals of this thesis are threefold. The first goal is to construct a suitable UQ method to use in CFD simulations. To this end, polynomial chaos (PC) expansions are explored as an efficient UQ method. A CFD model will be constructed where a single tube inside a tube bundle is subjected to a simple deformation. The model will then be treated with the PC method and the resulting output will be validated against a Monte Carlo (MC) simulation.

Chapter 1: Introduction

The second goal is to use the validated UQ method to study more complex bundle deformations, close to those found in actual reactors. The influence of geometric uncertainty will be reviewed on certain key output variables, which are: the normal force magnitude, direction and the drag force on the deformed inner tube. The total pressure drop over the length of the bundle will also be monitored. Only the mechanical and hydraulic properties of the flow will be examined, heat transfer is not considered.

In the last stage a FSI simulation will be setup with the same types of deformations as studied in the 1st and 2nd stage. Any interactions that may occur between the fluid flow and the deformed tube will be analysed.

Chapter 2

Literature Study

2.1

Tube Bundle Deformations

This section aims to give an overview of the most common sources of deformations and defects happening inside a nuclear reactor as well as the possible consequences that these might have. A reactor is constructed as a set of fuel assemblies (FA) (Fig. 2.1), which are stacked in a regular, usually square or hexagonal, pattern. Each (FA) is again made up of multiple fuel rods, which contains the nuclear fuel in the form of pellets. The FA is not completely made up of fuel rods alone, as certain positions are used to install instrumentation tubes, control rods, etc. Fig. 2.2 shows a typical FA and fuel rod with a hexagonal layout.

Figure 2.1: Cuttrough of the MYRRHA reactor core [1].

Figure 2.2: Diagram of a LMFR FA and fuel rod [2].

2.1.1 Deformations and Failure Mechanisms

As mentioned in the introduction, the environment of a nuclear reactor is very demanding. An LMFR, like MYRRHA, typically operate at a high temperature range (300_{− 700}◦_{C) and all the}

Chapter 2: Literature Study

equipment inside a reactor is constantly subjected to high doses of radiation. The most critical components are the cladding of the fuel rods and the wire spacers [2]. The wire spacer is used to prevent two adjacent tubes from hitting each other. The prolonged exposure to radiation causes changes at the molecular level of these materials (typical residence time _{∼ 2 years) and this can} cause the fuel rods to expand or warp. This phenomenon is known as void swelling (Fig. 2.3). The correct choice of material for the fuel rod cladding and wire spacer is the most important factor in mitigating the effects of void swelling.

Another important failure mechanism in fuel rods is rod fretting. This is caused by excessive vibrating of the fuel rods, incited by the fluid flow, which leads to repeated hitting of adjacent tubes. As mentioned before, spacer wires or in other cases spacer grids, are used to help prevent this. However, in some cases this only moves the problem and fretting occurs at the connection between grid and fuel rod (= grid-to-rod fretting), as shown in Fig. 2.4.

Figure 2.3: Examples of void swelling causing fuel rod deformations [2].

Figure 2.4: Example of fretting defect [3].

A last common source of failure and of particular importance in LMFRs is corrosion. The heavy metal LBE is known to be an accelerator of corrosion for the fuel rods if the choice of materials is not suitable. The parts of the fuel rod compromised by corrosion can become brittle, making them prone to failure (Fig. 2.5). The corrosion can also lead to a heavy build-up of oxide layers, which can restrict or block the flow channels [3].

2.1.2 Consequences of Deformations and Failures

The mentioned failure mechanisms affect the performance in different ways. The rod fretting (or corrosion) can cause local chips to break away from the fuel rod. These debris can clump together and completely or partially block a flow passage of the FA [12], [13]. This leads to a

Chapter 2: Literature Study

Figure 2.5: Fuel rod failure as a result of corrosion [3].

redistribution of the flow, which can cause an imbalance in the FA. The reduced flow velocity in the blocked passage also locally lowers the heat transfer capacity. This can cause the fluid to heat up excessively and lead to departure from nucleate boiling in those areas. This is less problematic in LMFR’s, due to the high boiling point of LBE. In very severe cases the fretting, or corrosion, can completely wear through the fuel rod cladding, releasing nuclear fuel into the coolant.

Apart from outright failure of the fuel rods, the deformations (commonly caused by void swelling) also have a detrimental effect on the reactor performance. An important difference between void swelling and fretting, is that the effects of void swelling are usually present on large parts or the complete FA. The fretting damage is much more localized.

An important consequence of void swelling is fuel assembly bow, where the complete FA is deformed. Different FA’s in the reactor core can experience FA bow simultaneously and these do not necessarily have to occur in the same direction [14]. Figure 2.6 shows the actual bow deformation from a FA observed in the Ringhals Vattenfal PWR plant (Sweden) [4].

(a) C-shape bow deformation (b) S-shape bow deformation

Chapter 2: Literature Study

The bow deformation has several negative effects. It can narrow fluid channel passages, which reduces the flow velocity and can cause hot-spots, as discussed above. It can also lead to an increased pressure drop over the tube bundle, which will affect the coolant circulation pump’s performance. Furthermore, the bow deformation will be subject to additional lateral forces generated by the fluid flow. This forces can amplify the bow deformation even more, creating a possible runaway-effect and increase the likelihood of fretting damage. Lastly if the FA bow is too excessive, there can be physical problems extracting the FA from the core during reloading. As the FA bow is mainly caused by void swelling, which is in turn induced by the high radiation level, one could opt to implement faster reloading cycles, limiting the exposure time of the fuel rods to the radiation. One could also lower the thermal output of the reactor core, which reduces the nuclear fuel burn-up and hence the radiation level. These measures however come at a clear economical cost.

2.2

Uncertainty Quantification

This section will focus on the different techniques used to model the uncertainty at the input of a system, handle the propagation of the uncertainty through that system, and quantify the resulting uncertainty in the output. A schematic of a general system with output y subjected to an uncertain input, which is composed of a set of deterministic components x and a set of stochastic, or random, components ξ, is given in Fig. 2.7.

Figure 2.7: Schematic of a system with an uncertain input, containing both deterministic and random components. The listed example causes of randomness are specific for the tube bundle case.

One aspect to point out is that uncertainty is not the same as error. A good definition of both error and uncertainty is given by Walters [5], “Error is a recognizable deficiency in any simulation or activity that does not arise from a lack of knowledge. Uncertainty can be defined as a recognizable deficiency in any simulation or activity that does arise from a lack of knowledge.” However it must be said that in practice, it often occurs that both terms are used somewhat

Chapter 2: Literature Study

interchangeably. A list of the most common uncertainties regarding CFD analysis is shown in Fig. 2.8.

Figure 2.8: Common sources of uncertainty and error in CFD according to [5]

It is believed that the most important uncertainty/error sources in CFD are the discretization error, the geometric uncertainty and the turbulence modelling [5]. As already mentioned, the focus will be put on the geometric uncertainty in this thesis.

Numerous mathematical methods have been developed that can quantify the uncertainty prop-agation in a set of variables, such as e.g. interval analysis or moment methods, as were used for CFD analysis by Huyse [15]. This thesis will focus on Polynomial Chaos (PC) expansions to use as an efficient UQ method. Monte Carlo (MC) simulations are also used, as a validation of the proposed PC methods. MC simulations have the characteristic of being very robust, making them the ideal candidate for this purpose.

2.3

Monte Carlo

One of the most straightforward and, as a consequence, one of the most widely used techniques for UQ are Monte Carlo (MC) methods [16]. The MC simulation’s origins lie in the Second World War, where one of its first uses was for the calculation of neutron reactions by von Neumann for the atomic bomb research [17]. The initial concept of MC is simple: a set of stochastic input variables that needs to be evaluated is sampled in a pseudo-random way to obtain a discrete

Chapter 2: Literature Study

sample set. Next, each sample is used to solve a computational model and afterwards the probabilistic properties from the output are retrieved through statistical analysis. For instance, the true mean µy and variance σy2 of an output variable y with an unknown distribution can

be estimated by respectively the sample mean Y and sample variance S2 if that sample set is sufficiently large: c µy = Y = 1 N N X i=1 Yi , cσy2 = S 2 = 1 N_{− 1} N X i=1 (Yi− Y )2 (2.1)

where N is the total number of samples and Yi represents one such sample.

The main advantages of MC are its simple implementation and its robustness. MC also allows the system itself to be treated as a black box as the uncertain input is converted into a discrete set of deterministic samples. This is interesting in the case of CFD simulations as specialized deterministic solvers can then be used without the need for adaptation. The main disadvantage is the slow convergence rate. For instance, using a pseudo-random sampling procedure, the standard deviation of the sample mean σ_Y convergences to the true mean at a rate inversely proportional to the square root of N

σ_Y = √σy

N (2.2)

A large number of samples is thus required before a conclusive result can be obtained. Different mechanisms exist to accelerate the convergence rate, such as for example the employment of more intelligent sampling schemes at the input [5]. However these mechanisms are not powerful enough to completely alleviate this problem and, combined with the already complex and time-consuming nature of CFD simulations themselves, leads to impractically long calculation times in all but the simplest of CFD cases.

2.4

Polynomial Chaos Expansions

Another technique for UQ is the use of Polynomial Chaos (PC) Theory. Here, instead of sampling the input as with MC, a solution is explicitly constructed by a decomposition of the model using spectral expansions. The spectral expansion of a random variable rests on the principle of separating the deterministic and stochastic parts from each other [16]. Suppose we have a random (output) variable y(x, ξ). The following then holds

y(x, ξ) = ∞ X i=0 yi(x)Ψi(ξ)≈ Pc X i=0 yi(x)Ψi(ξ) (2.3)

where x(x1, x2, ...) is the vector containing all deterministic variables and ξ(ξ1, ..., ξn) the

Chapter 2: Literature Study

represent the set of polynomial basis functions which are orthogonal to each other with respect to a weight function. Notice that the relationship y(x, ξ) and its decomposition is only exact when the number of expansion terms goes to infinity. In practice the series is of course truncated to a limited amount of terms. The inner product in the domain Ω of ξ is defined as

hΨi(ξ), Ψj(ξ)i =

Z

Ω

Ψi(ξ)Ψj(ξ)p(ξ)dξ (2.4)

where p(ξ) represents the weight function. The orthogonality condition for the basis functions can also be expressed as

hΨi(ξ), Ψj(ξ)i =

Ψ2_i(ξ)δij with δij = 0, if i6= j (2.5)

with δij being the Kronecker delta and hΨ2i(ξ)i being the inner product of Ψi(ξ) with itself.

The choice of the basis functions Ψi depends on the distribution of the random variables ξ.

The convergence rate of the solution is influenced heavily by the correct choice of polynomial functions [18]. The original chaos model, as developed by Wiener [19], used Hermite polynomials which are the optimal choice for a normally distributed input. For other distributions, the original Wiener polynomial chaos is expanded. The generalized Wiener-Askey polynomial chaos employs different subsets of the hypergeometric (orthogonal) polynomials to couple an input distribution with its optimal set of basis functions. The Hermite polynomials are also a subset of the hypergeometric polynomials. The classification of the different polynomial sets is achieved using the Askey scheme, hence the name.

Different probability distributions with their corresponding optimal polynomial sets and weight functions are shown in Table 2.1. Note that while only continuous distributions are shown, polynomial sets also exist that can efficiently handle discrete distributions, however these are not relevant to this thesis.

Distribution ξ Polynomial set Ψi(ξ) Weight function p(ξ) support

Normal Hermite e−ξ22 ]− ∞, ∞[

Uniform Legendre 1 [a, b]

Gamma Laguerre ξα_e−ξ _[0,∞[

Beta Jacobi (1− ξ)α_{(1 + ξ)}β _[0,∞[

Table 2.1: Different continuous distributions of one random variable ξ with their correspond-ing Wiener-Askey chaos polynomials and weight functions. α, β are probability parameters of the Gamma- and Beta-distribution; a, b_{∈ R .}

trans-Chapter 2: Literature Study

formed to or are the result from a transformation of another polynomial set [6]. For example the uniform distribution is a limit case of the beta-distribution with α = β = 0. For this rea-son the Askey scheme is commonly presented as a tree-like structure, with branches linking the polynomial subsets. The Legendre polynomials, up to an order of 4, are shown in Fig. 2.9 as illustration. The polynomial sets are, in their most basic form, only valid for a univariate random variable ξ. However it is fairly straightforward to find a multivariate expression, based on the 1D-formula, as long as the n random components of ξ all have the same probability distribution type (e.g. all normally distributed). The interested reader is referred to [16].

The amount of random variables n and the maximum order of the polynomials that one wants to use, also called the order of chaos p, determine the total size of the polynomial basis and thus the amount of coefficients yi(x) in the resulting PC model. It can be shown that the number of

coefficients Pc+ 1 is given by

Pc+ 1 =

(n + p)!

n!p! (2.6)

Note that the number of terms in the model goes up quickly with both rising order of chaos and number of random variables (in equal manner), this is shown in Fig. 2.9. This of course has a direct impact on the calculation effort required to solve the system. In practice, limits are imposed on both the maximum amount of random input variables that one can consider and the order of chaos of the model, depending on the complexity of the studied system and the available computational power. This is also referred to as the ‘curse of dimensionality’. Techniques exist that can improve these limits, with a small impact to the calculation time, which will be discussed further in Section 2.5.

−1.0 −0.5 0.0 0.5 1.0 ξ −1.0 −0.5 0.0 0.5 1.0 1.5 2.0 LPc (ξ ) L0(ξ) L1(ξ) L2(ξ) L3(ξ) L4(ξ) 1 2 3 4 5 6 7 order of chaos (-) 100 101 102 103 nr .of coef ficients (-)

n

Figure 2.9: Left: First five Legendre polynomials LPc(ξ), given for a random variable ξ de-fined on the interval [_{−1, 1].}

Right: Total number of PC coefficients in function of order of chaos p and the number of random variables n. n is an integer value, ranging from 1 to 5.