Numerical prediction of flow induced vibrations in nuclear reactor applications

Nuclear Engineering and Design 319 (2017) 81–90

Contents lists available at ScienceDirect

Nuclear Engineering and Design journal homepage: www.elsevier.com/locate/nucengdes

Numerical prediction of flow induced vibrations in nuclear reactor applications E. ter Hofstede a,b, S. Kottapalli a,b, A. Shams a,⇑ a b

Nuclear Research and Consultancy Group NRG, Petten, The Netherlands Delft University of Technology, Delft, The Netherlands

h i g h l i g h t s
Numerical simulations of flow induced vibration of nuclear fuel rods in axial turbulent flows.
Fluid-structure interaction simulations of strongly coupled problems.
Assessment of partitioned coupling algorithms and discretization schemes suitable for strongly coupled FSI problems.
Checking the effect of the use of different URANS models to simulate Turbulence Induced Vibrations.
Validating the results obtained from the URANS models to experimental test cases.

a r t i c l e

i n f o

Article history: Received 1 December 2016 Received in revised form 10 March 2017 Accepted 17 April 2017

Keywords: Flow induced vibrations Fluid-structure interaction URANS Nuclear reactor Strong-coupling

a b s t r a c t Flow induced vibration (FIV) plays an important role in nuclear industry. In nuclear reactors, the FIV are caused by a strong interaction between the fuel rods and the turbulent coolant flow around these rods. Numerical prediction of these strongly coupled fluid structure interaction (FSI) has been a challenge and is the main focus of this work. In this article, three aspects of FSI problems are numerically studied. In this first part, two different coupling schemes namely IQN-ILS and Guess-Seidel are thoroughly assessed for a laminar flow single rod experiment performed by Vattenfall (Lillberg et al. 2015). As a next step, a turbulent flow single rod experiment is selected to assess the capabilities of two different RANS models, i.e. a linear k  x SST and a non-linear Reynolds Stress Model. Lastly, an application based experiment of Chen and Wambsganns (1972) is selected to validate the combined effect of the selected RANS model along with the coupling method. Ó 2017 Elsevier B.V. All rights reserved.

1. Introduction In nuclear power plants, flow induced vibration (FIV) may cause fatigue problems, stress corrosion cracking, possible failure modes and fretting wear (Luk, 1993). This may eventually leads to nuclear safety issues and substantial standstill costs due to unplanned outage. Reports of flow-excited failures of heat exchanger tubes began appearing in the 1950s (Weaver et al., 2000). As was the case for the San Onofre Nuclear Generation Station, where FIV led to premature wear in over 3000 tubes, causing a leak of radioactive coolant in recently renewed steam generators in Units 2 and 3. The increase in power density of nuclear plants often results in an increase of coolant flow, or a change of cooling liquid, or a change in the component material or dimensions. These changes may alter the flow and structural behavior, and cause flow-induced vibrations to become more prominent (Weaver et al., 2000). It is ⇑ Corresponding author. E-mail address: [email protected] (A. Shams). http://dx.doi.org/10.1016/j.nucengdes.2017.04.026 0029-5493/Ó 2017 Elsevier B.V. All rights reserved.

therefore important to asses this phenomenon early in the design process. Because of this, the field of FIV is becoming an increasingly important area of research in the nuclear field. To correctly predict the FIV, a number of analytical models have been developed. However, these models are developed for slender bodies in axial flow by decomposing the fluid forces into an inviscid and viscous part. The contributions of these inviscid forces were first derived by Lighthill (1960), while the viscous forces were based on empirical relations, obtained from experiments on specific cases. The downside being the lack in accuracy when applied to other or less simplified cases. In a nuclear power plant, the use of high-density moderators leads to a low density ratio between the densities of the solid (qs ) and the fluid (qf ). Due to the high density fluid, the inertial forces that interact with solid bodies are of a higher order of magnitude. The high inertial forces have a direct effect on the dynamic behaviour of the solids, resulting in a strong coupling between fluid and the in-core structural elements. This phenomenon is also known as the Added-Mass Effect. When the

82

E. ter Hofstede et al. / Nuclear Engineering and Design 319 (2017) 81–90

added-mass effects are high, the coupling gets stronger, and hence, it is more viable to solve the problems numerically. Numerically simulating strongly coupled problems require a coupling between such fluid and structure solvers. Fluid-Structure Interaction (FSI) problems can be numerically solved either by monolithic or the partitioned approach. The Monolithic approach solves the fluid and solid equations simultaneously as one set of equations, which can solve fully coupled problems. These solvers are problem specific and implementation of a new code could be required depending on the simulation. On the other hand, partitioned approach makes use of existing fluid and solid solvers which are coupled through an external coupling algorithm. Thus, this provides more flexibility and reliability in terms of the used solvers. The coupling algorithm for a partitioned solver acts as a post-processing step for each solver (solid and fluid), it iterates between the data inputs and outputs from the solvers till the required conditions are satisfied. Since the fluid and solid equations are solved separately, the coupling algorithm introduces a coupling error within the solution. For loosely coupled problems partitioned approach is very efficient, since the coupling errors involved are very low. However, most FSI problems faced in nuclear applications are strongly coupled. When dealing with strongly coupled problems, it is shown that partitioned solvers suffer from poor convergence or even instability, therefore a strongly coupled method is often necessary (Causin et al., 2005; Yang et al., 2008; Banks et al., 2014; Degroote, 2013). With the progress in numerical methods, there have been several methods that have been introduced in the past, which can solve strongly coupled problems. In this study, two of these coupling algorithms, i.e. IQN-ILS and Gauss Seidel, are assessed for strongly coupled cases. The Gauss Seidel method is one of the earliest and a popular coupling scheme used in FSI solvers. Although, this method is simple and efficient in most of the cases, it has its own drawbacks, which would be discussed further in this article. The IQN-ILS is a Quasi-Newton method, which estimates the inverse Jacobian to attain the faster convergence. This is a state of the art coupling algorithm, which has shown better performance in some earlier studies (Degroote, 2010; Degroote, 2013; Banks et al., 2014). One of the aims of the present study is therefore to validate the capabilities of both the coupling methods to solve the strongly coupled problems. As a next step, one of these tested methods in combination with a turbulence model is used to predict flow induced vibrations. The description of these used coupling methods is given in Section 2. This is followed by the results and discussions related to the selected validation and the application cases in Section 3. The conclusions drawn from these test cases are summarized in Section 4. 2. Numerical methods In a partitioned approach, the fluid and solid equations are solved with different numerical methods. The fluid domain is solved using a computational fluid dynamic (CFD) method and the solid using a computational solid mechanics (CSM) method. To be able to model the interaction of these two models, the two methods are coupled at the fluid-solid interface. A stable and efficient numerical technique is essential for the study of FSI. In the case of strongly coupled interaction, the solvers are called multiple times during a time step until both the kinematic and the dynamic equilibrium conditions are satisfied. 2.1. Governing equations of fluid The Navier-Stokes (N-S) equations, govern the flow of fluid. The two main equations for an incompressible fluid in the N-S equations

are the conservation of mass and momentum. The mass and momentum conservation equations are given by Liu et al. (2013):

@ui ¼0 @xi

qf

ð1Þ


 @ui @ui uj ¼ r  rf þ f f þ @t @xj

ð2Þ

with, u being the fluid velocity, qf the fluid density, t the time, rf the fluid stress tensor and f is the body force. For Newtonian fluids, the stress tensor can be written as:



r ¼ pI þ l

@ui @uj þ @xj @xi



with, p the pressure and

ð3Þ

l the fluid dynamic viscosity.

2.2. Governing equations of solid The deformation of an elastic incompressible solid is governed by the conservation of momentum:

qs

@ 2 ds ¼rrþf @t 2

ð4Þ

with, ds being the displacement of the structure, qs the solid density, r the solid stress tensor and f the body forces. 2.3. Coupling between the fluid and solid equations At the interface between the solid and the fluid, the kinematic condition requires the velocity of the fluid to be equal to the time derivative of the displacement of the solid interface:

@ds ¼ uf @t

ð5Þ

The dynamic condition requires the stress on the interface due to the fluid and solid surface normal to be equal, also called as equality of traction:

rs  ns ¼ rf  nf

ð6Þ

The subscript s and f depicts the structural solver and fluid respectively. When using a Dirichlet-Neumann decomposition, the flow equations can be solved for a given velocity (or displacement) at the fluid-structure interface (Dirichlet boundary condition), and the solid equations are solved for a given traction distribution on the interface (Neumann boundary condition). Assuming, the displacement vector as, x ¼ ds and traction as y ¼ rf  nf , the response of the structure solver can be therefore written as:

x ¼ Sð y Þ

ð7Þ

where SðyÞ is the function resolving the structure equations (Eq. (4)). Similarly, for fluid solver:

y ¼ F ðxÞ

ð8Þ

where FðxÞ is the function resolving the fluid equations (Eqs. (1) and (2)). At each time step, the fixed point equation must be satisfied

x ¼ S  F ðxÞ

ð9Þ

where S  FðxÞ ¼ SðFðxÞÞ. Therefore, the residual function, RðxÞ is calculated as

RðxÞ ¼ S  F ðxÞ  x

ð10Þ

To reach convergence, this residual should be minimized, which can be done using an optimization algorithm. Well known approaches in FSI are the Gauss-Seidel method, fixed underrelaxation, Aitken under-relaxation, and the IQN-ILS method (Degroote, 2013).

E. ter Hofstede et al. / Nuclear Engineering and Design 319 (2017) 81–90

The IQN-ILS method is an efficient algorithm for strongly coupled problems (Degroote, 2010; de Ridder et al., 2013). The IQNILS method stands for Interface Quasi-Newton with an Inverse of the Jacobian by Least-Squares approximation. Unlike, the Newton-Raphson method, IQN-ILS estimates an approximate Jacobian. Hence, the convergence is not quadratic. However, it is significantly better than zeroth order methods. The IQN-ILS coupling scheme can be explained as follows:

~i ¼ xi1 þ Dxi x

ð11Þ

83

STAR-CCM+ is a commercial CFD package, widely used in industries. STAR-CCM+ v10.06 release introduced the fluid structure interaction solver, which couples FEM and FVM solvers (Bloxom, 2014). FEM is well suited for Solid Mechanics problems. The introduction of FEM-FVM coupling makes the code very appealing when simulating complex structural entities. Presently, only GaussSeidel coupling algorithm has been implemented in the code.

3. Results and discussion

!

where, Dxi is an estimate of the correction in the value of x for the iteration i. Therefore, the residual is given by,

 i ~ x ~i Ri ¼ S  F x

ð12Þ

Therefore for the next iteration, the correction factor Dxiþ1 is estimated by,

¼ Ri  1 ¼ Ri dR dx

dR Dxiþ1 dx

Dxiþ1

ð13Þ

 1 Here, dR is known as the inverse Jacobian, which is estidx mated using the method of least squares explained in detail by Degroote (2013). The Gauss-Seidel method is another coupling method used to couple the fluid and solid solvers. The Gauss-Seidel method (Degroote, 2010; Degroote, 2013) is a zeroth order method which does not involve the estimation of the Jacobian. The information from previous iterations is used to estimate the value of the variable for the next iteration. Unlike the previous method the Gauss-Seidel method approximates its Inverse Jacobian to be an identity matrix (I). Since, the Jacobian need not be calculated; the Gauss- Seidel method takes less time per iteration.


1 dR ¼I dx

ð14Þ

Gauss-Seidel method can be unstable for strongly coupled problems. Therefore, an under-relaxation factor is used to increase the stability of convergence. A fixed under-relaxation factor (x), or an adaptive under-relaxation factor (e.g. Aitkens under-relaxation) is used to stabilize the convergence.

dR1 dx

¼ xI; x  constant under-relaxation actor

dx

¼ xk I; xk  Aitkens under-relaxation factor

dR1

ð15Þ

Since, IQN-ILS uses the approximate Inverse-Jacobian to converge to a solution, the convergence rates of IQN-ILS are better than Gauss-Seidel method, which uses the Identity matrix. Moreover, the Inverse-Jacobian is updated every iteration which improves the convergence rates. In the current study both the IQN-ILS method (Degroote et al., 2009; Degroote, 2013) and the Gauss-Seidel method (Degroote, 2013) are applied. 2.4. Numerical code The IQN-ILS and the Gauss Seidel coupling methods are validated using two software packages, OpenFOAM Extend (OpenFOAM) and STAR-CCM+ (STAR-CCM+, 2015) respectively. OpenFOAM Extend is a fork of OpenFOAM and an open source CFD software package also with some capabilities to solve CSM. The fluid and the structure equations uses FVM discretization in OpenFOAM. The IQN-ILS coupling algorithm has been implemented in OpenFOAM for solving FSI problems (Campbell, 2010). The considered version of OpenFOAM extend lacks FE capabilities, therefore, the FV structure solver has been used to perform FSI simulations with OpenFOAM.

The study is set to attain two main objectives. The first objective is to compare and study the accuracy and efficiency of using FVMFEM coupling to the FVM-FVM coupling. To achieve the first objective, Vattenfall-I FSI test case (Lillberg et al., 2015) is performed with the aforementioned methods. The second objective of the study is to validate the use of URANS turbulence models to simulate flow induced vibrations, and to predict the dynamic behaviour of structure. To achieve the second objective, two test cases are performed, the Vattenfall-II test case by Lillberg (2015), and the experimental case by Chen and Wambsganns (1972). These test cases focus on a simplistic models of fuel rods in nuclear reactors. Nuclear reactor applications are the main driving force behind the selection of the test cases. 3.1. Vattenfall-I FSI case Vattenfall Research and Development (VRD) in Sweden has performed an experiment where a vertical slender rod/beam with a roller boundary condition and a clamped side is displaced and then released, after which its vibrations and damping under the influence of the surrounding fluid are investigated (Lillberg et al., 2015; Vu and Truc, 2014). The flow of water is in longitudinal direction of the beam through the channel. The center of the rod is given an initial displacement of 10 mm after which it is released to freely vibrate, shown in Fig. 1. The displacement of the rod is measured using a laser. This case is well suited to test the influence of the added mass on the structure. The stainless steel rod has a length of L ¼ 1:5 m, and a density of qs ¼ 8000 kg=m3 , the modulus of elasticity of E ¼ 193 GPa and the moment of inertia, I ¼ 8:53  1010 m4 . For a fixed-pinned beam, the analytical model to calculate the natural frequency of the first mode in vacuum is give by Timoshenko (1953)

f ¼

 sffiffiffiffiffiffiffiffiffiffi 1 15:418 EI 2p M=L L2

where E is the Young’s Modulus, I is the area moment of inertia, M is the mass and L is the length of the beam. The experimental natural frequency of the beam vibration without water was found to be 12:3 Hz. The analytical approximation of the beams frequency results in 12:37 Hz, which is agreement with the experiment result. For numerical computations, a hexahedral mesh is used for the fluid and the solid domains. With water as the surrounding fluid, the density ratio (qs =qf ), of the experimental case is  8. For the FSI simulation, the IQN-ILS coupling algorithm with FVM-FVM coupling in OpenFOAM extended. Whereas, the Gauss-Seidel coupling method was employed using FEM-FVM coupling in STAR-CCM+. A mesh and time convergence study is performed for the solid solvers (both FVM and FEM) to check their influence on the final solution. It is worth mentioning that for the fluid solver, a second order upwind scheme was used for spatial discretization and second order backward difference scheme was used for the temporal discretization in OpenFOAM. Whereas, a central differencing scheme with skew correction is used for the spatial discretization, whereas, second

84

E. ter Hofstede et al. / Nuclear Engineering and Design 319 (2017) 81–90

Fig. 1. Vattenfall-I case domain (Lillberg et al., 2015; Vu and Truc, 2014). Here, L ¼ 1:5 m; H ¼ 80 mm; h ¼ 8 mm; and l ¼ 20 mm. The initial displacement d ¼ 10 mm.

order backward scheme is used as the temporal discretization in STAR-CCM+. 3.1.1. Mesh convergence study for the solid steel beam The mesh convergence study for the solid steel beam was carried out for both the codes to check the effects of discretization schemes on the solution. To perform the study, the beam is subjected to an initial displacement of 10 mm in the vertical direction as depicted in Fig. 1 and is then allowed to vibrate freely in vacuum. A fixed-pinned boundary condition has been applied to the beam at the inlet side and outlet respectively. The solid beam mesh study was performed by varying the number of divisions in the longitudinal direction. This study was performed in the 2D as well as the 3D case. Initial study was performed with FVM discretization for structural cases using OpenFOAM extend. Fig. 2a, b depicts the effect dependence of the solution with respect to number of

divisions in the longitudinal direction for a 2D and 3D beam respectively. The term DSRS160 in the figures denotes a 2-dimensional mesh with 160 divisions in the longitudinal direction for the beam. Similarly, the term DSRS160_3D denotes a 3-dimensional mesh with 160 divisions in longitudinal direction. Table 1 gives a complete overview of the study performed for FVM discretization study. The modal frequencies shown in the table are calculated by performing FFT (Fast Fourier Transforms) (Heideman et al., 1984) on the vertical displacement of the beam center patch obtained form the simulation. The FFT was performed the complete displacement signal for a simulation time of 1 s to compute the natural frequency. From Table 1 it can be concluded that the FV discretization methods for solving structural equations is highly sensitive to mesh refinement. The mesh sensitivity in longitudinal direction

Fig. 2. The amplitude of displacement magnitude of the beam center has been plotted with respect to time. The beam is vibrating in vacuum. The figure shows the effect of mesh refinement on the computed natural frequency of the solid beam.

85

E. ter Hofstede et al. / Nuclear Engineering and Design 319 (2017) 81–90 Table 1 Effect of mesh refinement on the calculated natural frequency of the beam through FV discretization method.

Table 2 Effect of mesh refinement on the calculated natural frequency of the beam through FE discretization method.

Case

Number of Cells

Modal Frequency [Hz]

Error percentage

Case

Number of divisions

Frequency [Hz]

Experiment DSRS80 DSRS160 DSRS160_3D DSRS320 DSRS640 DSRS640_3D DSRS1280

– 486 966 2898 1926 3846 11538 7686

12.30 20.16 15.96 15.13 13.45 12.71 12.32 12.71

– 63.9 29.8 23.0 9.35 3.33 1.33 3.33

disp_100e disp_150e disp_300e disp_650e

100 150 300 650

12.34 12.34 12.34 12.34

is very high. Whereas, the mesh sensitivity from 2D to 3D conversion is lower. In case of FE discretization for structures, the simulation was performed in STAR-CCM+. Similar to the previous study the mesh refinement in the longitudinal direction was varied, by changing the number of divisions. The number of mesh divisions are varied from 100 (disp_100e) to 650 (disp_650e), this effect is noticed in Fig. 3. It can be observed that the mesh sensitivity of FE discretization method is negligible when compared to the FV method. The results virtually are not affected by the change in the number of mesh cells. This can be observed in the Table 2.

3.1.2. Time convergence study of the solid steel beam Time convergence study for the beam is performed for temporal discretization with FVM and FEM discretization. The boundary and initial conditions are similar to the study for mesh convergence. A mesh containing 160 cells in longitudinal direction (SRS160), was used. The beam solved with FVM discretization was subjected to time steps from 102 s (SRS160dt02) to 105 s (SRS160dt05). The effect of each of the cases is depicted in Fig. 4a. Since the beam is a freely vibrating beam, it should not damp (neglecting mass damping). However, it can be observed from the graph that time steps P 103 s lead to excessive numerical damping. A time step of dt 6 105 s is required to keep lower the effect of numerical damping. This explains that FV discretization is sensitive to the time step used for the unsteady simulations. For FE discretization, a similar study was performed. It was found that for FE discretization, the use of time step dt 6 5  103 s gives results closer to the experimental result. Moreover, it was also observed that FEM solvers do not suffer from the problem of numerical damping. Table 3 compiles the results that were obtained during the study. The results can also be observed in Fig. 4b.

From the above study it becomes very clear that, in comparison to FVM solvers, FEM solvers have advantage in terms of mesh sensitivity and also in terms of sensitivity to time steps. Therefore, it is more advantageous to make use of FEM rather than FVM to solve structure equations. This conclusion would be put to use in the test cases that would discussed in the later sections of this paper.

3.1.3. Results of the FSI simulation The FSI simulation was carried out in 3-D for both discretization schemes to get a better comparison of results, and to have a better co-relation with the experimental set up. Three cases with water inlet velocities 1 m/s and 3 m/s were simulated with both the codes. The bulk Reynolds number based on the bulk velocity and half channel height is Reb ¼ 3600 at an inlet velocity of 1 m/s, which is laminar in nature. The inlet of the fluid mesh is provided with a velocity inlet boundary condition, and the outlet is given a pressure outlet boundary condition. The outer walls are given no slip boundary conditions. At the interface of fluid and structure, a zero relative velocity is provided with respect to the interface itself. A fixed-pinned boundary condition was provided at the ends of the beam. The solid beam was provided an initial displacement of 10 mm at the center of the beam. This was carried out by applying an external force and carrying out a steady state simulation for the beam. Later, the force was removed, to allow the beam to vibrate freely. The simulation is then coupled with the fluid domain to replicate the conditions observed in the experiment. The simulation was performed with FEM-FVM and FVM-FVM coupling. The FVM-FVM coupling using the IQN-ILS algorithm and the FVMFEM coupling using the Gauss-Seidel coupling algorithm. The results of the simulation are shown in Fig. 5a & b and the comparison of the frequencies are shown in Table 4. In the literature study (Degroote et al., 2009; Degroote and Vierendeels, 2012; Blom et al., 2015) it has been clearly seen that for a given time-step, IQN-ILS coupling method has a higher convergence rate (lesser number of iterations). Therefore, one might opt to choose IQN-ILS over Gauss-Seidel. However, in the current scenario due to limitation of the structure solver, the time step

Fig. 3. Effect of mesh refinement on the computed natural frequency of the solid beam by Finite Element (FE) discretization. The amplitude of displacement magnitude at the beam center has been plotted in the figure. All the response curves of the beam overlap each other even with varying grid refinement.

86

E. ter Hofstede et al. / Nuclear Engineering and Design 319 (2017) 81–90

Fig. 4. Effect of time step on the computed natural frequency of the solid beam.

Table 3 Effect of mesh refinement on the calculated natural frequency of the beam through FE discretization method. Case

Natural Frequency [Hz]

Error percentage

Analytical dt_0.1 ms dt_1 ms dt_5 ms dt_10 ms

12.37 12.49 12.26 12.15 11.76

– 0.97 0.89 1.78 4.93

used for the FSI case in FVM-FVM coupling was 0.1 ms compared to a time-step of 1 ms in the case of FEM-FVM coupling. Thus, providing a major advantage in terms of the computational costs. The FEM-FVM coupling lead to the simulation being at least 10 times more cheaper, with the added advantages of a coarser mesh for the structure solver in case of FEM (150 cells in longitudinal direction) compared to FVM solid solver (640 cells in longitudinal direction). Hence, from the results two major things can be concluded; the FEM-FVM coupling is a more reliable and cheaper way of performing FSI simulations compared to FVM-FVM coupling. The second conclusion being, that the use of a better discretization scheme in an FSI simulation far outweighs the advantages of the use of a better coupling algorithm. 3.2. Vattenfall-II case study (fuel rod in a confined space) In the past, simulation of FSI cases on turbulence induced vibrations (TIV) have been performed for fuel rods with the high fidelity fluid solvers like LES (Liu et al., 2013). These cases mainly involved relatively low Reynolds number flows which cannot represent an industrial application scenario. Therefore, to reduce the costs of

FSI simulations, the use of URANS models have to be explored. To validate the use of URANS models, an FSI simulation of an experimental test case by the Vattenfall Research and Development (VRD), Sweden has been performed. Vattenfall-II case is an application based case which involves simulation of vibrations induced on a steel rod due to turbulent flow of water around it. The experimental set-up is shown in the Fig. 6. The bounded white are shown in the figure is cavity of fluid flow. The central rod is shown by the circle within the white cavity. Note that there is no flow of fluid within the circle where the rod is situated. The case simulates the vibrations on a fuel rod in a nuclear reactor surrounded by flow of turbulent water. The length of the steel rod is 1.486 m with a cylindrical crosssection. The outer diameter of the tube, Do ¼ 8 mm with a thickness of, t ¼ 0:6 mm. The inlet side of the tube is fixed (Zero displacement and no rotation) and the outlet side of the tube is pinned (zero displacement, but free to rotate in transverse direction). The domain consists of a wide inlet followed by sudden constriction of the cross-sectional area. The sudden decrease in the area of the fluid flow leads to a turbulent flow in the water. The constricted area again opens up into a wide area. The mass flow rate of water is 10 kg/s. The Bulk Reynolds number for the flow is, Reb  1:2  106 , based on the hydraulic diameter of the tube. The pressure fluctuations resulted due to the turbulent flow around the tube results in vibrations of the steel tube. This simulation aims to capture the dynamic behavior of the steel tube in the given conditions and compare it to the experimental results. It can be observed that for large domains like the Vattenfall-II, to perform an FSI simulation with a high-fidelity fluid solver such as LES (Large Eddy Simulation), the number of cells required is in the order of 107 cells and Courant numbers 6 0:2 (ter Hofstede,

E. ter Hofstede et al. / Nuclear Engineering and Design 319 (2017) 81–90

87

Fig. 5. Comparison of the results for water inlet velocities 1 m=s and 3 m=s. The amplitudes are measured as displacement magnitudes at the center of the beam in longitudinal direction.

Table 4 Comparison of results of the Vattenfall case-I FSI simulation. Result

Natural Frequency [Hz]

Error %

Experiment 1 m=s FVM-FVM coupling FEM-FVM coupling

10.41 11.21 10.22

– 7.68 1.83

Experiment 3 m=s FVM-FVM coupling FEM-FVM coupling

10.26 11.17 10.16

– 8.86 0.55

2015). This increases the computational costs of an FSI simulation by several folds. Therefore, the use of URANS models is being validated in this particular case study. The main objective of this current study is to test the effect of turbulence models on the dynamic behaviour of the structure. 3.2.1. Simulation setup In the framework of this test case, an assessment with two URANS models is performed. These models include a linear k  x SST (Menter, 1994) and a non-linear RSM (Reynolds Shearstress Model) with Elliptic Blending (Launder et al., 1975; Manceau and Hanjalic´, 2002). The two models present two different assumptions for turbulence, the former assumes isotropic condition, and the latter assumes anisotropic behaviour of turbulence. This simulation would also compare the variation in the results for the two models with respect to experimental measurements by VRD, Sweden (Lillberg, 2015). The density ratio of the present configuration is qs =qf ¼ 8. A time step of 1  103 s was used to simulate the case, with a hexahedral trimmed mesh for the fluid and a hexahedral mesh for the structure. The number if cells in the fluid

region is  4  106 cells. The first cell thickness of the mesh has a wall yþ  5. The structure region is a hexahedral mesh with number of cells  20; 000. The simulation considers a FE discretization for the structure solver and FV for the fluid solver. Similar to the previous case, the Gauss Seidel coupling algorithm was used to couple the two solvers. Table 5 gives the exact specifications of the material properties which were used in the experiment. An initial perturbation is provided in form of a uniform body force, which is applied initially and then removed. The sudden removal of the force triggers the vibrations. Due to the uniform loading, mainly the first mode of the tube is excited. The simulation set up for both the turbulence models are given in Table 6. 3.2.2. Results and discussion The vibrational frequencies of the rod are calculated by performing FFT on the obtained displacements. The results of which have been tabulated in Table 7. It can be noted here that although the frequencies in both the cases (models) have been observed to be close to the experimental values. The frequencies of each mode can be obtained by providing the initial perturbation accordingly. It can be observed that the frequencies in the x-direction and the y-direction are similar to each other. Moreover, the xdirection frequency deviates from the experimental results by approximately 10%, in the y-direction however, the frequency prediction is more accurate. There can be several reasons behind the deviation. In the experiment, the vibrations induced onto the structure are the result of random pressure fluctuations that are convected throughout the domain, due to the constriction of flow near the inlet. The LES of the fluid domain performed by ter Hofstede (2015), show the absence of large structures present

88

E. ter Hofstede et al. / Nuclear Engineering and Design 319 (2017) 81–90

Fig. 6. Vattenfall case-II domain specifications.

Table 5 Material properties in the Vattenfall case-II experiment.

X-dir

Y-dir

The second observation to be noted in the results, is that the modal frequencies obtained from the structure solver are not directly dependent on the type of URANS model used. The k  x SST model is based on the Boussinesq approximation of isotropic turbulence, whereas, the RSM model is based on calculating the turbulent quantities using transport equations for the Reynolds stress. Both the models, eventually produced similar results. This fact leads to a conclusion that the calculations of turbulent quantities do not have a direct role to play in the considered FSI simulations. The main quantities which are used are the time averaged quantities of pressure (p), calculated from RANS equations. High fidelity solvers such as LES are capable of simulating the high frequency pressure fluctuations which lead to TIV (Liu et al., 2013). URANS models can capture pressure oscillations which occur due to large structures (low frequencies) de Ridder et al. (2016) but cannot capture high frequency fluctuations which are generated mainly due the turbulent pressure and velocity fields. It can be concluded that URANS can be potentially used to calculate the eigen-frequencies of the structures when subjected to vibrations due to pressure fluctuations resulting from turbulent flows. Such vibrations are also known as Turbulence Induced Vibrations (TIV). However, the main drawback of the method is the requirement of an initial perturbation to the structures, which externally excites the structure for a brief period of time, unlike observed in turbulent induced vibrations.

8.29 9.67 9.93

9.52 9.96 9.96

3.3. Flexible brass beam in turbulent water

Fluid Properties (Water) Density (qf )

1000 kg=m3

Kinematic Viscosity (mf )

1  106 m2 =s

Structural Properties (Steel) Density (qs ) Elasticity Modulus Poisson Ratio

8055 kg=m3 193 GPa 0.285

Table 6 Turbulence Properties at Inlet of the fluid domain for k  x SST model and the Elliptic-Blending RSM. Turbulence Properties Turbulent Intensity Turbulent Length scale

1% 2.5 mm

Table 7 Compiled results for experimental results and the simulated results. Results

Experiment k  x SST Model RSM Model

Natural Frequency ðHzÞ

within the flow. Therefore, URANS models fail to capture the effects of the flow on the tube. The vibrations captured are hence a result of a time averaged flow which are uniform in radial directions and since the geometry of the domain is symmetric about the X and Y axis, the modal frequency of structure in both the directions is expected to be the same. Another reason for the errors may be due to non conformal meshes. The non conformal mesh lead to errors in mapping of forces from the fluid to the structure domain therefore causing the deviations. Generation of a mesh finer of fluid domain does not necessarily address this problem as the fluid points at the interface will become even more dense, as compared to structural mesh points, leading to greater errors.

A Flexible beam in turbulent water flow is an experimental case by Chen and Wambsganns (1972). This involves an annular flow domain across a cylindrical brass beam. The hydraulic diameter of the fluid domain is, Dh ¼ 0:0127 m. The outer-wall is a rigid steel cylinder, the vibrational frequency and the damping ratios are then measured, for three constant inlet velocity conditions, 10 m/s, 20 m/s and 30 m/s. The beam is clamped (fixed) on both the ends. The beam has a core diameter of Dc ¼ 0:0127 m and the annulus has an outer diameter of, Do ¼ 0:0254 m. The length of the cylinder is L ¼ 1:19 m. The geometry of the domain can be visualised in Fig. 7. The main objective of the test case is to investigate the effect of flow rate on the natural frequency of the structure. For this case, three inlet velocities are used, 10 m/s, 20 m/s and 30 m/s.

E. ter Hofstede et al. / Nuclear Engineering and Design 319 (2017) 81–90

89

Fig. 7. Geometry of the domain for the flexible beam in water experiment.

Table 8 Material properties used in flexible beam in water experiment. Fluid Properties (Water) Density (qf )

1000 kg=m3

Kinematic Viscosity (mf )

1  106 m2 =s

Structural Properties (Brass) Density (qs ) Elasticity Modulus Poisson Ratio

8400 kg=m3 107 GPa 0.331

3.3.1. Simulation Set-up The FSI simulation was performed with water inlet velocities 10 m/s, 20 m/s and 30 m/s. The k  x SST turbulence model was used for the three inlet conditions. The turbulence intensity at the inlet of the domain is set to 5% with a turbulent length scale of 0.1 cm. The material properties of fluid and the structure are shown in Table 8. An impulse force of 0:2 m=s2 was applied for a period of 0.1 s, at a central patch of the beam. The impulse force has a Gaussian distribution curve as a function of time. The vibrations imposed due to the impulse force were analysed using FFT. The damping ratios of the beam were calculated using a least square method. The simulation was performed till a steady state was attained of the solid.

3.3.2. Results and discussion Initial simulation of a freely vibrating beam in vacuum was performed to calculate its natural frequency. The natural frequency of the beam in vacuum was calculated to be 28.49 Hz analytically and 28.57 Hz through simulation. The experimental measurements (Chen and Wambsganns, 1972) show that the beam had a natural frequency of 28 Hz in quiescent water. A simulation in quiescent water with similar boundary conditions was performed, and the frequency was calculated to be 26 Hz. This suggests that the boundary conditions of the beam in simulation do not completely match the boundary conditions of the beam in the experiment. This can be due to a possibility of pre-tension on the beam during the experiment, which has not been recorded. A higher frequency of the beam is possible if the structure is pre-stressed. However the paper by Chen and Wambsganns (1972) does not mention any changes related to boundary conditions, therefore, a deviation in the results of about 10% can be observed. Fig. 8 shows a comparison and the trend of the modal frequencies obtained as a function of velocity. It can be observed that the natural frequency of vibration decreases by 0.4 Hz, as the average velocity increases by 10 m/s. The similar trend has also been observed in the experimental results for natural frequencies. On the other hand, the damping ratios in both the cases are a linear function of the flow velocity.

Fig. 8. Plots comparing the results obtained from experimental and simulation values.

90

E. ter Hofstede et al. / Nuclear Engineering and Design 319 (2017) 81–90

4. Conclusion

Acknowledgements

An extensive validation of two different coupling methods, namely the FVM-FVM and the FEM-FVM are performed for strong coupled FSI problems. The IQN-ILS coupling method is used in combination with an FVM-FVM solver in OpenFOAM. Whereas, the Gauss Seidel coupling scheme is used in combination with an FEM-FVM solver in STAR-CCM+. Three experimental test cases, the Vattenfall-I, Vattenfall-II and the Flexible brass beam in turbulent water experimental case by Chen and Wambsganns (1972) were selected to perform the validation study. The discretization schemes were validated for Vattenfall-I case. A convergence study was performed for the solid solver to check the influence of mesh and time step on the results, as well as the effect of using a FEM or FVM solid discretization method. It was found that FEM is less sensitive to the mesh size of the solid compared to the FVM solver. In addition, the FEM-FVM coupling is 10 times more efficient and leads to a faster convergence compared to an FVM-FVM coupling. With the respect to the coupling algorithms, i.e. the IQN-ILS and the Gauss Seidel method, it can be concluded that both methods are capable of solving the FSI problems with strong coupling. The use of the FVM-FEM combination resulted in a computationally cheaper solution; it was preferred over using the FVMFVM combination. To benefit from this, the FVM-FEM coupling with Gauss Seidel scheme was used to simulate the application based problem, i.e. Vattenfall-II and the Flexible brass beam in turbulent water case. The simulation of turbulence induced vibrations in the Vattenfall-II case was performed to validate the effect of URANS models on prediction of the natural frequencies of the structural components. Results conclude that the modal frequencies of the structures can be predicted accurately. Two different turbulence models were used to simulate TIV. The k  x SST and the RSM model with Elliptic Blending. Both the models represented different assumptions for the modeling of turbulence. The test was performed to validate whether a change in the assumptions of the flow model could effect the results of the FSI simulation. The conclusion was that a change in the URANS models do not have direct influence on the selected FSI simulation for simulating TIV. The traction on the structure is mainly calculated with the use of time averaged quantities. The URANS models in an FSI simulation of TIV can be used for calculation of modal frequencies and damping ratios. The turbulence models are incapable of simulating external excitation on the structure solver. This resulted in the damping out of the vibrations that are initiated through the external perturbation. The simulation of the Flexible brass beam in turbulent water case was performed to validate the effects of the varying inlet velocities on the natural frequency of the structure. Three inlet velocities 10 m/s, 20 m/s, and 30 m/s were selected to validate the test case. Although a deviation of 10% was observed in the values of natural frequencies, the trend of the natural frequency dependence on the flow velocities were similar. The natural frequencies are observed to be a linear function of the flow velocity with a negative slope. Whereas, the damping ratios was found to be an increasing linear function of the inlet velocity. Hence, it can be concluded that the behaviour of the beam is accurately predicted in comparison to the experiments through the use of URANS models.

The work described in the paper is funded by the Dutch Ministry of Economic Affairs. The Authors would like to acknowledge the contribution of Eric Lillberg (VRD) and David Blom (TUD) for their support with the experimental data and with the coupling solver respectively. References Banks, J.W., Henshaw, W.D., Schwenddeman, D.W., 2014. An analysis of a new stable partitioned algorithm for fsi problems. part i: Incompressible flow and elastic solids. J. Comput. Phys. 269, 108–137. Blom, D.S., van Zuijlen, A.H., Bijl, H., 2015. Multi-level acceleration with manifold mapping of strongly coupled partitioned fluid-structure interaction. Comput. Methods Appl. Mech. Eng. 296, 211–231. Bloxom, A.L., 2014. Numerical Simulation of the Fluid-structure Interaction of a Surface Effect Ship Bow Seal (Ph.D. thesis). Virginia Polythechnic Institute and State University Blackburg, VA. Campbell, R.L., 2010. Fluid-structure Interaction and Inverse Design Simulations for Flexible Turbomachinery (Ph.D. thesis College of Engineering). Pennsylvania State University. Causin, P., Gerbeau, J.F., Nobile, F., 2005. Added-mass effect in the design of partitioned algorithm for fluid-structure problems. Comput. Methods Appl. Mech. Eng. 194, 4506–4527. Chen, S.S., Wambsganns, M.W., 1972. Parallel flow induced vibration of fuel rods. Nucl. Eng. Des. 18, 253–278. Degroote, J., 2010. Partitioned simulation of interaction between an elastic structure and free surface flow. Comput. Methods Appl. Mech. Eng. 199, 2085–2098. Degroote, J., 2013. Archives of Computational Methods in Engineering. Chapter Partitioned Simulation of Fluid-Structure Interaction, Coupling BlackBox Solvers with Quasi-Newton Techniques. Springer International Publishing. Degroote, J., Vierendeels, J., 2012. Multi-level quasi-newton coupling algorithms for the partitioned simulation of fluid-structure interaction. Comput. Methods App. Mech. Eng. Degroote, J., Bathe, K.J., Vierendeels, J., 2009. Performance of a new partitioned procedure versus a monolithic procedure in fluid-structure interaction. Comput. Struct., 793–801 de Ridder, J., Degroote, J., Tichelen, K., Schuurmans, P., 2013. Modal characteristics of a flexible cylinder in turbulent axial flow from numerical simulations. J. Fluids Struct. 43, 110–123. de Ridder, J., van Tichelen, K., Degroote, J., Vierendeels, J., 2016. Vortex-induced Vibrations by Axial Flow in a Bundle of Cylinders. In: FIV. Heideman, M.T., Johnson, D.H., Burrus, C.S., 1984. Gauss and the history of the fast fourier transform. IEEE ASSP Mag. 1, 14–21. http://dx.doi.org/10.1109/ MASSP.1984.1162257. Launder, B.E., Reece, G.J., Rodi, W., 1975. Progress in the development of a reynoldsstress turbulent closure. J. Fluid Mech. 68 (3). Lighthill, M., 1960. Note on the swimming of slender fish. J. Fluid Mech. 9, 305–317. Lillberg, E., 2015. Energiforsk FIV Experiment, Data for Computational Analysis. Technical Report Vattenfall Stockholm, Sweden. Lillberg, E., Angele, K., Lundqvist, G., Edh, N., 2015. Tailored experiments for validation of cfd with fsi for nuclear reactor applications. In: NURETH16. Chicago, IL. Liu, Z.G., Liu, Y., Lu, J., 2013. Numerical simulation fluid-structure interaction for two simple fuel assemblies. Nucl. Eng. Des. 258, 1–12. Luk, K., 1993. Pressurised Water Reactor Internals Aging Degradation Study. Technical Report Oak Ridge National Laboratory. Manceau, R., Hanjalic´, K., 2002. Elliptic blending model: a new near-wall reynoldsstress turbulence closure. Phys. Fluids 14, 744–754. http://dx.doi.org/10.1063/ 1.1432693. Menter, F.R., 1994. Two-equation eddy viscosity turbulence models for engineering applications. AIAA J. 32, 1568–1605. OpenFOAM (). User Guide, Version 2.2.0. OpenFOAM Foundation. STAR-CCM+, 2015. User Manual. CD ADAPCO London. ter Hofstede, E., 2015. Numerical Study of Fluid Structure Interaction in Nuclear Reactor Applications (Masters thesis TU Delft). Timoshenko, S., 1953. History of strength of materials. McGraw-Hill New York. Vu.,Truc, T., 2014. Testing and optimization of unicorn fluid-structure Interaction solver for simulating an industrial problem (Master’s thesis KTH Computer Science and Communication). Weaver, D.S., Ziada, S., Au-Yang, M.K., Chen, S.S., Païdoussis, M.P., Pettigrew, M.J., 2000. Flow-induced vibrations in power and process plant components progress and prospects. J. Pressure Vessel Technol. 122, 339–348. Yang, J., Preidikman, S., Balraras, E., 2008. A strongly coupled embedded boundary method for fluid structure interaction of elastic mounted rigid bodies. J. Fluids Struct. 24, 167–182.