- Email: [email protected]
Numerical simulation of the flow in wire-wrapped pin bundles: Effect of pin-wire contact modeling
Nuclear Engineering and Design 253 (2012) 374–386
Contents lists available at SciVerse ScienceDirect
Nuclear Engineering and Design journal homepage: www.elsevier.com/locate/nucengdes
Numerical simulation of the flow in wire-wrapped pin bundles: Effect of pin-wire contact modeling夽 E. Merzari a,∗ , W.D. Pointer a , J.G. Smith a , A. Tentner a , P. Fischer b a b
Nuclear Engineering Division, Argonne National Laboratory, 9700 S. Cass Avenue, Argonne, IL 60439, USA Mathematics and Computer Science Division, Argonne National Laboratory, 9700 S. Cass Avenue, Argonne, IL 60439, USA
a r t i c l e
i n f o
Article history: Received 14 March 2011 Received in revised form 28 March 2011 Accepted 9 September 2011
a b s t r a c t The rapid advancement of numerical techniques and the availability of increasingly powerful supercomputers recently enabled scientists to use large eddy simulation (LES) to simulate numerically the flow in a full subassembly composed of wire-wrapped pins. Because of the extreme computational cost of such simulations, it was not possible to conduct a sensitivity case on the pin-wire interface modeling. Since such calculations are likely to be extended to conjugate heat transfer, however, a sensitivity study is necessary to assess the reliability of the numerical results. It is well known, for example, that conjugate heat-transfer results are often strongly influenced by near-wall modeling. The objective of the present work is to investigate the effect of pin-wire contact modeling from the point of view of both the hydraulics and the heat transfer characteristics. In particular, the focus is on the prediction of the hot spot in conjugate heat-transfer calculations. The primary test case is the simplified geometry recently proposed by Ranjan et al., which consists of a simple channel flow with a wire embedded in one of the walls. After reproducing the results using the LES code Nek5000, we examined several other choices for the wire-pin interface modeling, including the introduction of a nominal gap between the wire and wall. The results shed light on the sensitivity of CFD calculations results to the modeling of the interface region between wires and pins. © 2011 Elsevier B.V. All rights reserved.
1. Introduction The study of conjugate heat transfer is needed to support accurate simulations of the flow and temperature fields of a nuclear reactor core. With computational fluid dynamics-based techniques and the use of supercomputing platforms, one can now achieve a remarkable degree of accuracy for the most relevant turbulence factors, such as the average velocity field, the averaged temperature distribution, and the stresses inside the fluid region. The coupling of the fluid dynamics problem with the conduction equation solved inside the fuel pin is less straightforward, however. In particular, the choices made in modeling the contact between the fuel pin and other structural components (such as spacers or wires in the case of sodium-cooled fast reactors) might lead to
夽 The submitted manuscript has been created by UChicago Argonne, LLC, Operator of Argonne National Laboratory (“Argonne”). Argonne, a U.S. Department of Energy Office of Science laboratory, is operated under Contract No. DE-AC02-06CH11357. The U.S. Government retains for itself, and others acting on its behalf, a paid-up nonexclusive, irrevocable worldwide license in said article to reproduce, prepare derivative works, distribute copies to the public, and perform publicly and display publicly, by or on behalf of the Government. ∗ Corresponding author. Tel.: +1 6302526637. E-mail address: [email protected] (E. Merzari). 0029-5493/$ – see front matter © 2011 Elsevier B.V. All rights reserved. doi:10.1016/j.nucengdes.2011.09.030
sensitive perturbations on the prediction of the peak temperature in the solid (hot spot), the peak temperature in the fluid or other thermodynamic and hydrodynamic properties, since such details have effects on other geometries (Qu and Mudawar, 2002). For instance, Pointer et al. (2009) recently used large eddy simulation (LES) to simulate numerically the flow in a full subassembly composed of wire-wrapped pins. Because of the extreme computational cost of such simulations, it was not possible to conduct a sensitivity study of the pin-wire interface modeling. Since such calculations are likely to be extended to conjugate heat transfer, however, a sensitivity study is necessary in order to assess the reliability of the numerical results. In fact, the wire introduces a higher resistance to conduction in the radial direction which, combined with the high conductivity of sodium, might lead to a significant localized hot spot. The presence of local stagnation points might also affect the temperature distribution in the fluid (Tang et al., 1978). The way the interface is modeled has the potential to affect the value of the hot spot, since it would affect the hydrodynamics and the local heat transfer. Moreover, in some instances, the wire might be in contact with the pin, depending on design and operational conditions. Thus, to the uncertainty due to numerical modeling one should add uncertainty due to the variability of the physical conditions. The importance of these points is shown clearly in Fig. 1 which presents the results
E. Merzari et al. / Nuclear Engineering and Design 253 (2012) 374–386
Nomenclature y+ xi x, y, z ui u
v w T T˜ T˜max,f
␥ Rez Rex Wbulk Ubulk D h g e
r 2␣
normalized wall distance Cartesian coordinates (vector notation) Cartesian coordinates Cartesian velocity components Cartesian velocity component in direction, x Cartesian velocity component in direction, y Cartesian velocity component in direction, z temperature normalized temperature normalized maximum temperature in the fluid linear gradient coefficient density kinematic viscosity Reynolds number in direction z (defined in the text) Reynolds number in direction x (defined in the text) bulk velocity in direction, z bulk velocity in direction, x diameter of the wire half width of the channel gap size fictitious parameter representing the amount of contact. For cases with no contact it is equal to g. For cases with contact it is equal to h/2(1 − cos(˛/2)) distance from the bottom wall contact angle
of a conjugate heat transfer calculation performed for conditions similar to the Advanced Burner Test Reactor (Chang et al., 2006). Details of the calculation, performed with STAR-CCM+ 5.02 (STARCCM+ 5.02 Users Guide) are reported in Table 1. To simplify the problem periodic boundary conditions are assumed on the cross section (Fig. 1) while inlet–outlet boundary conditions are assumed in the streamwise direction (the whole core length was simulated, with uniform axial heating). The results shown are mesh converged (an example of the coarsest mesh used is reported in Fig. 1b). Fig. 2a and b shows a snapshot of the temperature distribution in an axial cross-section. The temperature peak is located downstream of the wire near the contact between wire and pin. These results are consistent with Gajapathy et al. (2007). Fig. 2c shows the difference between the temperature distribution and the linear gradient (Fig. 3) associated with uniform axial heating. The perturbations on the linear gradient are due to the presence of the wire. The main purpose of the present work is to shed some light on how pin-wire contact modeling affects the temperature prediction in a simplified geometry. Arguably, such a study would be problematic in an actual core because the hot-spot position and value are determined essentially by the axial position (due to the nonuniformity of the power distribution; see Khakim et al., 2009, and Baglietto and Ninokata, 2007) and it would be more difficult to separate the small but potentially significant effect due to the modeling of the wire or other structural components (Fig. 2c). To investigate accurately the influence of such effects it would be necessary to perform several full pin-length calculations with highfidelity models for each given power shape configuration, which would require a substantial computational effort. Another approach is to assume that the effects due to the non-uniformity of the power shape and the wire-pin contact modeling are separable. In such case a single full-length calculation could be performed assuming the correct power distribution and a pre-set contact modeling. Then a set of computationally much less expensive calculations can be performed assuming a uniform power distribution, which allows for the use of periodic boundary
375
conditions, to investigate the effect of the wire contact-modeling. By properly rescaling the temperature differences due to the contact modeling with the power distribution it is possible to determine how such temperature differences vary with the location in the reactor core (Section 3.2). In this paper we will provide an example of how such methodology would be used. The primary test case is the geometry proposed by Ranjan et al. (2010), which consists of a simple channel flow with a wire embedded in one of the walls. Ranjan et al. provided detailed and extremely accurate distributions for the velocity field and other turbulent statistics (i.e., obtained with highly resolved direct numerical simulation techniques), thus allowing us to fully verify and validate our calculations. After reproducing the results using the LES code Nek5000 (Fischer et al., 2008), we examined the conjugate heat-transfer problem (which was not studied by Ranjan et al. (2010)). We also considered an additional modeling option for the wire: the introduction of a nominal gap between the wire and wall. We then used the commercial code STAR-CCM+ in RANS mode to study the same problem, with the double purpose of verifying the results obtained and testing the capability of RANS-based methods to achieve accurate results for the present case. As a last step we investigate in greater detail, using STAR-CCM+, the influence of the contact modeling with a higher number of cases, and discuss the relevance of these results for a full pin-wire conjugate heat transfer calculation. 2. Methodology In this section we describe the wall-wire contact models used. We present the mathematical model, the relevant parameters chosen, and the methodology used for the simulations. 2.1. Problem description The direct numerical simulation (DNS) calculation of Ranjan et al. (2010) was selected for its generality, its high resolution (it was computing over 100,000,000 collocation points), and its relatively low Reynolds number. The geometry investigated in that paper corresponds to a wire immersed in a channel flow with both streamwise (direction z) and spanwise (direction x) periodic Table 1 Conditions for the preliminary case. Description
Value
Total power Number of active core assemblies Number of pin per assembly Power profile Pin diameter Pin pitch Pin material
250 MWth 54 217 Flat 8 mm 9.05 mm Metallic uranium (gap and cladding are omitted for simplicity) 1 mm 0.15 m 7.5 0.75 m Sodium 5.8 m/s 600 K k–ε realizable, two-layer 23 Inlet–outlet
Wire diameter Wire pitch Contact angle Pin length in the core region Coolant Inlet velocity Inlet temperature Turbulence model Maximum wall y+ Boundary conditions (streamwise direction) Buoyancy Discretization Number of meshes for convergence
No Second order About 8,000,000 polyhedral with prismatic layer
376
E. Merzari et al. / Nuclear Engineering and Design 253 (2012) 374–386
Fig. 1. Preliminary case. Boundary conditions (a) and coarse mesh for the fluid (b) and the entirety of the domain (c).
boundary conditions (the boundary conditions are identical for Nek5000 and STAR-CCM+). The geometry is shown in further detail in Fig. 4. In the present work, while keeping the same outline for the geometry, we explore three LES wall-wire contact models. In the first case we keep the same contact modeling used by Ranjan et al. (2010), where the wire and wall are considered in contact and merged with an angle of 7.5◦ . In the second case the wire and wall are separated by a nominal gap g. For the present case Lx = 4h, Lz = 8h, and D = h. The nominal gap is set to equal to g = 0.05 h, in order to allow a small but significant cross-flow. As a last step a simulation was performed with a gap g = 0.05 h between wire and wall with a baffle blocking the cross flow below the wire. No heat transfer calculation was performed in the work by Ranjan et al. (2010). In the present simulations a heat source is specified within the walls while it is assumed that no heat is generated in the wire. After the LES cases, a series of RANS cases was simulated using STAR-CCM+ 5.02 at different contact angles and gap sizes. A full list of the cases run can be found in Table 2. 2.2. Mathematical model The conjugate heat-transfer problem is solved by resolving the incompressible, constant-property, Navier–Stokes equation in the fluid (1)–(3) and the conduction equation in the solid (4): ∂ui =0 ∂xi
(1)
∂u 1 ∂p ∂2 ui ∂ui + uj i = − + ∂xi ∂t ∂xj ∂xj ∂xj
(2)
∂T ∂T ∂2 T = ˛f + uj ∂t ∂xj ∂xj ∂xj
(3)
∂T ∂2 T + q, = ˛s ∂xj ∂xj ∂t
(4)
where u is the velocity field, T is the temperature, is the viscosity, ˛ is the thermal diffusivity (s in the subscript refers to solid, f for the fluid), is the density of the fluid, p is the pressure, t is the time, and x is the Cartesian coordinate. Here q is not the volumetric heat source but rather the ratio between the volumetric heat source and the thermal capacity of the solid; it can be seen as a temperature source. This source is set to a constant value in the heated solid and to zero in the wire. Table 2 List of cases run for the present paper. Case
Code and method
Gap size
Contact angle
I II III IV V VI VII VIII
Nek5000/LES Nek5000/LES Nek5000/LES STAR-CCM+/RANS STAR-CCM+/RANS STAR-CCM+/RANS STAR-CCM+/RANS STAR-CCM+/RANS
0.05 D 0.05 D (with baffle) – – – – 0.025 D 0.05 D
– – 7.5 15 7.5 3.25 – –
E. Merzari et al. / Nuclear Engineering and Design 253 (2012) 374–386
377
Fig. 2. Preliminary case, Temperature in the cross section (a) with detail (b) and three-dimensional view of the wall surface temperature difference (difference between the temperature and the mean fluid temperature at each height) (c). Contour plot velocity components in the same plane of (a) velocity in direction y (d) and velocity in direction z (e).
Unlike in the case presented in Figs. 1–3, periodic boundary conditions are assumed in the streamwise and spanwise directions for the velocity, with a constant pressure gradient in both directions determined dynamically in order to maintain a constant
mass-flow rate. Similar configurations are common in the simulation of parallel channel flows. In the case of conjugate heat transfer, special attention needs to be paid to the energy equations, since the temperature field is not periodic. If the temperature source
378
E. Merzari et al. / Nuclear Engineering and Design 253 (2012) 374–386
Fig. 3. Temperature as a function of z at the center of the computational domain.
distribution is constant, however, the problem can be reduced to a periodic one for the temperature as well (Iaccarino et al., 2002). The temperature in a configuration such as the one examined here can be written as: T = T0 + z + T˜
(5)
which represents a linear gradient in direction z (streamwise direction) superposed on a periodic temperature field T˜ . The ramp of the gradient can be determined by the conservation of energy principle (since the energy is evacuated from the system only through the fluid, if higher-order axial conduction is neglected (Ghaddar et al., 1986): =q
cp,s AS 1 cp,f Af Wbulk
(6)
Here, AS /Af is the ratio between the area of the solid and the fluid (on a plane normal to z), and Wbulk is the streamwise bulk velocity. We note that depends on the ratio between the heat capacity of the solid and the heat capacity of the fluid cp,s /cp,f . In fact, is a measure of how easily the fluid is able to evacuate heat from the solid. If the heat capacity of the fluid is significantly higher than the heat capacity of the solid, the same temperature source in the solid will cause a smaller, streamwise fluid temperature increase. A high value of indicates poor cooling, whereas a low value indicates effective cooling. Substituting (5) into (3) and (4) leads to the following set of equations for the temperature: ∂T˜ ∂2 T˜ ∂T˜ + uj = ˛f − u3 ∂t ∂xj ∂xj ∂xj
(7)
∂2 T˜ ∂T˜ = ˛s +q ∂t ∂xj ∂xj
(8)
The equations are balanced in the sense that the energy produced in the solid is dissipated in the fluid. No energy is transferred at the boundaries. Therefore, periodic boundary conditions can be applied for the temperature field T˜ . 2.3. Choice of parameters The relevant parameters for the following problem are the Reynolds numbers in direction x and z; the Prandtl number; the ratio between the thermal diffusivity in the solid and the thermal diffusivity in the liquid; and, to a minor extent, the ratio between the thermal capacities. The Reynolds number in direction z can be defined as: Rez =
Wbulk h ,
(9)
which is chosen to be equal to 6000 (consistently with Ranjan et al. (2010), albeit slightly higher). Wbulk is the bulk velocity in the streamwise direction. The Reynolds number in direction x can be defined as: Rex =
Ubulk h ,
(10)
which is chosen to be equal to 1000 (consistently with Case C in Ranjan et al. (2010), albeit slightly higher). Ubulk is the bulk velocity in the spanwise direction. The Prandtl number can be defined as: Pr =
, ˛f
(11)
which is chosen to be equal to 0.01, a value representative of a liquid metal coolant. The ratio between the two thermal diffusivities is
Fig. 4. Three-dimensional representation of the case with (left) contact between wire and wall and (right) no contact between meshes.
E. Merzari et al. / Nuclear Engineering and Design 253 (2012) 374–386
379
Fig. 5. Grid structure for the nek5000 simulation: (a) very coarse grid (fourth-order polynomial) for the no-contact case; (bottom) detail (b). CCM mesh (c) near the contact between wire and wall.
assumed to be equal to 10 (the thermal diffusivity in the liquid being 10 times higher than the solid), and the ratio between the thermal capacities is assumed to be equal to 1. This corresponds loosely to the case of a liquid metal coolant and a uniform metallic fuel. We note that the source term does not appear as a relevant parameter of choice. The reason is the following. There exists a normalization of the temperature T˜ (and notably L, where L is a lenghtscale of arbitrary choice) for which the source term can be removed from the equation. Since T˜ is a fluctuation on the linear streamwise gradient, it is intuitive that the choice of q is irrelevant for the determination of its shape. An increased value of q would have a significant effect on , however, thereby influencing the distribution of T.
On the upper and lower boundaries for the solid, a zero heat flux condition is applied (the temperature derivative is set to zero on the boundary). 2.4. Numerical methods The simulations were carried out by using the spectral element code Nek5000 in LES mode on two Argonne computer systems—Cosmea (a 32-node, 128-core Linux cluster) and Blue Gene/P (an IBM system with 40,960 quad-core compute nodes)—with up to 15,000,000 collocation points. An example of structure of the macroscopic grid used is shown in Fig. 5 for the streamwise-normal cross section; the figure depicts the elements and the actual collocation points for a very coarse mesh with no
380
E. Merzari et al. / Nuclear Engineering and Design 253 (2012) 374–386
Fig. 6. Instantaneous streamwise velocity in a z-normal plane (a) and time-averaged streamwise velocity in a z-normal plane (b). Spatial units are normalized by h.
contact between wire and wall. The discretization is designed to allow for at least one point near the wall at y+ < 1 and five points within y+ < 10. Spectral elements methods are a class of higher-order methods in space. Polynomial functions of up to the 11th degree have been used to discretize the velocity field in each element. In the generalized weighted residual framework, the present spectral element method can be classified as a Galerkin method where the test functions and the basis functions for each element are Lagrange polynomials evaluated on Gauss–Lobatto–Legendre collocations points (for the velocity). In LES, large-scale turbulence is simulated while smaller scales are modeled. Since smaller scales have a nearly universal behaviour, LES is a more reliable methodology than is Reynolds-averaged Navier–Stokes (RANS), in the sense that it generally depends less on the model details. It is also preferable to use in place of direct numerical simulation because it can be performed at a reasonable computational cost. The contribution of the smaller scales to the energy cascade is in the present case modeled through a local, element-based explicit cut-off filter in wave-number space. The energy is removed from the smallest simulated scales (high wavenumbers), thus mimicking the effect of smaller eddies (Fischer and Mullen, 2001). This procedure allows for a coarser grid and therefore lower computational cost. Time advancement has been carried out through an explicit third-order, backward finite-difference scheme. An additional benefit of this method is that it does not explicitly define an eddy viscosity, leaving the definition of a turbulent Prandtl number as superfluous. The Nek5000 code has been extensively validated and is massively parallel. In the present case the time stepping was chosen to have a maximum Courant number no higher than 0.4. For comparison, some calculations were also run by using the commercial code STAR-CCM+ in the RANS framework with a realizable k–ε turbulence model and all-y+ wall modeling. The number of meshes used was considerably lower for the same cases (300,000 meshes). Polyhedral cells with a prismatic layer were used everywhere except in the contact region (Fig. 5c): local refinement was applied to this region to reach convergence at a lower cost. A
segregated flow solver and a second order scheme for the convective terms were also employed. The turbulent Prandtl number was assumed constant and equal to 0.9. While this approach presents serious limits for the simulation of the flow in liquid metal systems (Tzanos, 2011), it proved to yield accurate results for simplified geometries such as pipes and channels. 3. Results and discussion In the case of LES, several cases were run corresponding to two different geometries and different mesh resolutions. While the macroscopic mesh was not changed during the simulations, mesh convergence was proven by increasing the order of the polynomial basis from the 4th order up to the 11th. For each RANS case a mesh convergence study was conducted. Among other verification tests, each result shown here has been obtained on at least two separate meshes with a polyhedral mesh size ratio equal to 1.5. 3.1. Velocity distribution We begin by comparing the velocity field with the data of Ranjan et al. (2010). Fig. 6a shows an instantaneous contour plot for the streamwise velocity, while Fig. 6b shows a contour plot for the averaged streamwise velocity; both figures refer to the case with contact between wire and wall. Compared with similar plots published by Ranjan et al. (2010), these figures show a good degree of accuracy considering that the Reynolds number is 10% higher than that for the DNS case. Fig. 7 shows four profile comparisons of the results obtained with Nek5000 (with contact modeling) and the DNS data. The accuracy is good considering the enormous difference in resolutions (we used less than one-tenth of the collocation points) and the difference in Reynolds number. All profiles refer to time-averaged data. (For the position in x, see Fig. 7). Fig. 8 shows a comparison between the Nek5000 results for the case with no contact between wire and wall and the DNS data. Such a comparison would be meaningful, for validation purposes, only
E. Merzari et al. / Nuclear Engineering and Design 253 (2012) 374–386
381
Fig. 7. Comparisons of the streamwise velocity and the cross velocity at x/h = 0.0 (wire centerline) and x/h = 1.0. The solid line represents the DNS data. All values are normalized by Wbulk . r represents the distance from the bottom wall.
Fig. 8. Comparisons of the streamwise velocity at x/h = 6. Values are normalized by Wbulk . The solid line represents the DNS data; r represents the distance from the bottom wall.
far from the wire. We therefore did the comparison at x/h = 6.0: the agreement seems reasonable and within the scope of the present work. Time-averaging convergence was verified. Near the wire, the two cases (i.e., no contact and contact between wire and wall) behave differently. In the case of no contact, the velocity distribution appears to be altered. Fig. 9 compares the two Nek5000 results at x/h = 1.0. Near the bottom wall (r/h = 0.0, where r is the distance from the bottom wall), the streamwise velocity is lower, an effect that could be related to the presence of a higher cross-flow in this area. In fact, as it can be seen in Fig. 10, the recirculation pattern near the gap is altered when compared to the case with no contact. Moreover the cross velocity in the narrow gap is significantly higher than the surrounding areas. This effect is nontrivial, and it is a relevant from the point of view of the hydrodynamics. To determine if it the cross flow in the case with no-contact is the source of the disturbance in the streamwise velocity distribution, an additional case has been simulated. A conductive baffle was introduced between wire and wall (by fixing the velocity components to zero), thus blocking the flow passage below the wire. Fig. 11 demonstrates how the streamwise velocity distribution for this case matches well the velocity distribution for the case with contact.
382
E. Merzari et al. / Nuclear Engineering and Design 253 (2012) 374–386
Let us, however, assess in further detail the magnitude of this influence
3.2. Temperature
Fig. 9. Comparison of the streamwise velocity between the two Nek5000 cases at x/h = 1. Values are normalized by Wbulk ; r represents the distance from the bottom wall. The solid line represents the case with contact between wall and wire.
It can be concluded that the wall-contact modeling has a definite impact on the streamwise velocity distribution in the region downstream of the wire. In particular, a case with no contact would introduce a small but significant cross flow between wire and pin which in turn affects the stream-wise velocity distribution by continuity and transport of momentum. Such change in the streamwise velocity distribution would have an impact on the local heat transfer coefficient and therefore on the local temperature distribution.
Fig. 12a shows the average temperature distribution for the case with contact between wire and wall while Fig. 12b shows the average temperature distribution for the case with no-contact between wire and wall. All values are normalized by D where D is the diameter of the wire. The plots in Fig. 12a and b are similar except for the wire region. A hot spot can be observed in both cases near x = 0 and r = −h (near the lower adiabatic boundary), with the hot spot in the no-contact case shifted toward x > 0. Note that the scales are different. It is evident however, that in the case with contact between wire and pin the hot spot is slightly higher. To provide a more quantitative comparison, we are also reporting the temperature profiles for x/h = 0 in Fig. 13. While the two profiles are different and the case with contact clearly predicts a higher hot spot, the difference appears to be small. We conclude that the significant difference observed in the hydrodynamic side had little effect on the hot-spot prediction for the solid. We note, however, that the temperatures of the fluid at the wall near the contact region between wire and pin is significantly affected by the contact modeling. In the case with no contact the peak T˜ in the fluid was evaluated to be −31.5 a.u., while in the case with contact it was evaluated to be −4.77 a.u. This is a significant temperature increase (about 7.5% of the maximum temperature difference within the model). The peak fluid temperature is an important design parameter because it determines the margin to boiling and other technological considerations. Before drawing any further conclusions, however, we will proceed in examining in further details the results; Figs. 14 and 15
Fig. 10. Cross velocity: (a) vector plot for the case with no contact; (b) near-gap vector plot for the case with contact; (c) contour plot of the velocity in direction x for the case with no contact.
E. Merzari et al. / Nuclear Engineering and Design 253 (2012) 374–386
Fig. 11. Comparison of the streamwise velocity between the Nek5000 cases at x/h = 1. Values are normalized by Wbulk ; r represents the distance from the bottom wall.
provide further verification of the results obtained. Fig. 14 shows the average temperature for the STAR-CCM+ case with contact between wall and wire. Fig. 15 shows a comparison between the temperature profile at x = 0 obtained by Nek5000 and that by STARCCM+ for the case with contact between wire and wall. The results are in very good agreement. The peak fluid temperature for the same case is estimated by STAR-CCM+ at −36.2 a.u., thus with
383
a discrepancy of approximately 1.5% of the maximum temperature difference. Two very different codes provided very similar results, giving us confidence on the validity of the conclusions that can be drawn. In particular it is likely that the uncertainty due to the contact modeling is more significant than the uncertainty due to the turbulence modeling. Moreover for the present case RANS appears sufficiently accurate to predict the value of the hot spot. Let us therefore examine in greater detail the impact of the wall-wire contact modeling using the cheaper RANS methodology. A series of cases have been run in STARCCM+ changing the type of contact modeling, which can be parameterized based on the distance e (Fig. 16a). A positive e represents cases with no contact modeling (and e therefore represents the distance between wire and pin). A negative e represents cases with contact: the higher the absolute value of e the higher will be the contact angle ˛. For each different case the results were verified to be mesh converged: this was particularly difficult for small values of e, where a very strong local refinement near the wire-wall interface was required (up to 9,000,000 meshes). No fillet was applied. Fig. 16b reports the peak fluid temperature as a function of e (a fictitious parameter that accounts for the amount of contact) for all the cases run in the present work. It is clear that the peak fluid temperature is not a monotonic function of e: in fact the function has a maximum between e = −0.0016 D and e = 0.025 D based on the data collected. However, it can be assumed that the peak will actually be at e = 0.0, where the wire and the wall form a singularity since the presence of a conductive layer of fluid (e > 0) would decrease the peak temperature. In fact, the baffle case—despite the hydrodynamic in all similar to the case with contact and e = −0.016 D presented a temperature drastically lower and similar to the no-contact case with e = 0.05 D: since the fluid is more conductive than the solid the fluid layer between wire and wall helps reducing the temperature in the wire region (Fig. 13). Moreover, in low Prandtl flows, such as the one considered here, conductivity plays a much stronger role
Fig. 12. Average temperature T˜ distribution for the case with no contact (a) and average temperature T˜ distribution for the case with contact (b).
384
E. Merzari et al. / Nuclear Engineering and Design 253 (2012) 374–386
Fig. 13. Temperature T˜ profile at x/h = 0.0 for the two cases. Detail of the distribution near the lower adiabatic boundary is shown in the right box. The solid line represents the case with contact between wall and wire.
Fig. 14. Temperature T˜ for the case with contact, obtained by using STAR-CCM+.
in heat-transfer than in water and other medium to high Prandtl fluids. The fact that the peak fluid temperature has a peak near e = 0.0 D warrants caution in the estimation of peak temperature for a real bundle case. From Fig. 16a it appears evident that, in order to be conservative a simulation should assume contact, with a sufficiently small contact angle. Therefore a reasonable uncertainty due to the contact modeling should be assumed in the range of 2.5–5% of the maximum temperature difference (all of this is, of course, valid only for the present set of parameters). Let us now go back to the preliminary case examined in the introduction (Figs. 1–3, Table 1). The maximum temperature difference within the whole model at each given cross section varied between 80 K and 90 K (Fig. 17). In the fluid, the variation was considerably lower (with a variation consistently around 13 K at the wall, when the linear gradient is subtracted, Fig. 2c). Normalizing by ␥D where D is the diameter of the wire, we obtain a spread of 380 a.u. thus consistent with the values reported earlier (368 a.u. Fig. 12a). A 10 a.u. error band as suggested in Fig. 16b would be equivalent to a ±2.1 K band in the pin-wire conjugate heat transfer simulation. However, if no-contact cases are considered a larger
Fig. 15. Temperature T˜ profile at x/h = 0.0 for the case with contact, obtained by using Nek5000 and STAR-CCM+. The solid line represents the case performed using Nek5000.
E. Merzari et al. / Nuclear Engineering and Design 253 (2012) 374–386
385
Fig. 16. Peak Temperature T˜ for the cases considered in this work as a function of e the straight dashed lines represent the level of uncertainty recommended for the present case (a plus or minus 10 a.u. band around the value predicted with a reasonable contact angle).
Fig. 17. Temperature distribution in the solid for the preliminary calculation (details in Table 1).
band should be assumed (up to 8.2 K for the calculations performed here, to make sure to include the zero point in Fig. 16b). 4. Conclusions A series of LES and STAR-CCM+ simulations have been carried out on a simplified geometry in order to investigate the influence of wire-pin contact modeling on the prediction of the temperature distribution.
The results presented here suggest that the modeling of the contact between wire and wall is significant, if not crucial, from the hydrodynamic point of view. It is less so from the point of view of the hot-spot prediction, since a significant difference in hydrodynamic behaviour leads to an approximated one percentage point difference in the value of the hot spot. The peak temperature of the fluid, on the other hand, appears to be significantly affected. This warrants caution in the use of CFD for the prediction of the peak surface temperature, which should be carried out assuming a small contact angle (as small as numerically permitted). These results have then been related to a full single pin-wire conjugate heat transfer simulation with periodic boundary conditions. The simulation, albeit simplified, provided a reference and allowed us to give a quantitative value to the uncertainty in the pinwire contact modeling: for a case similar to the ABTR, depending on the simulation performed such uncertainty has been estimated between 2.1 K for a small contact angle simulation and up to 8.2 K for a no-contact simulation. While this temperature differences might appear small it is important to point out that they might increase substantially when a more realistic simulation is taken into account (the presence of cladding and a gap for instance) and constitute a significant percentage of the temperature difference in the cross section. An additional conclusion that can be drawn is that RANS-based techniques are sufficient to reach accurate results for the prediction of the hot spot and peak fluid temperature for the present geometry and parameters. This seems to be due to the conduction-dominated nature of heat transfer for low-Prandtl fluid conjugate heat transfer cases. Acknowledgments This work was supported in part by the Office of Advanced Scientific Computing Research, Office of Science, U.S. Department of Energy, under Contract DE-AC02-06CH11357.
386
E. Merzari et al. / Nuclear Engineering and Design 253 (2012) 374–386
References Baglietto, E., Ninokata, H., 2007. “Improved turbulence modeling for performance evaluation of novel fuel designs”. Nuclear Technology 158 (2), 237–248. Chang Y. I., Finck, P.J., Grandy C., 2006. “Advanced burner test reactor: preconceptual design report”, Argonne National Laboratory report ANL-ABR-1 (ANL-AFCI-173). Fischer, P.F., Lottes, J., Pointer, W.D., Siegel, A., 2008. Petascale algorithms for reactor hydrodynamics. Journal of Physics: Conference Series 125, 012076. Fischer, P.F., Mullen, J.S., 2001. Filter-based stabilization of spectral element methods. In: Comptes Rendus de l’Académie des sciences Paris, 332, Série I. Analyse numérique, pp. 265–270. Gajapathy, R., Velusamy, K., Selvaraj, P., Chellapandi, P., Chetal, S.C., 2007. CFD investigation of helical wire-wrapped 7-pin fuel bundle and the challenges in modeling full scale 217 pin bundle. Nuclear Engineering and Design 237, 2332–2342. Ghaddar, N.K., Karniadakis, G.E., Patera, A.T., 1986. A conservative isoparametric spectral element method for forced convection; application to fully developed flow in periodic geometries. Numerical Heat Transfer B: Fundamentals 9 (3), 277–300.
Iaccarino, G., Ooi, A., Durbin, P.A., Behnia, M., 2002. Conjugate heat transfer predictions in two-dimensional ribbed passages. International Jounal of Heat and Fluid Flow 23, 340–345. Khakim, A., Merzari, E., Ninokata, H., 2009. Feasibility study of the application of exotic pin for tight-lattice fuel assembly. Transactions of the American Nuclear Society 100, 747. Pointer, W.D., Smith, J.G., Fischer, P., Siegel, A., 2009. RANS simulations of turbulent diffusion in wire-wrapped sodium fast reactor fuel assemblies. In: abstract in International Conference on Fast Reactors and Related Fuel Cycles, Kyoto, December, pp. 249–250. Qu, W., Mudawar, Issam, 2002. Analysis of three-dimensional heat transfer in micro-channel heat sinks. International Journal of Heat and Mass Transfer 45, 3973–3985. Ranjan, R., Pantano, C., Fischer, P., 2010. Direct simulation of turbulent swept flow over a wire in a channel. Journal of Fluid Mechanics 651, 165–209. STAR-CCM+ 5.02 User Guide. Tang, Y.S., Coffield Jr., R.D., Markley Jr., R.A., 1978. Thermal analysis of liquid-metal fast breeder reactors. American Nuclear Society, ISBN: 0-89448-011-1. Tzanos, C., 2011. Improved simulation of heat transfer in liquid-metal flows. Nuclear Technology 174, 41.
Contents lists available at SciVerse ScienceDirect
Nuclear Engineering and Design journal homepage: www.elsevier.com/locate/nucengdes
Numerical simulation of the flow in wire-wrapped pin bundles: Effect of pin-wire contact modeling夽 E. Merzari a,∗ , W.D. Pointer a , J.G. Smith a , A. Tentner a , P. Fischer b a b
Nuclear Engineering Division, Argonne National Laboratory, 9700 S. Cass Avenue, Argonne, IL 60439, USA Mathematics and Computer Science Division, Argonne National Laboratory, 9700 S. Cass Avenue, Argonne, IL 60439, USA
a r t i c l e
i n f o
Article history: Received 14 March 2011 Received in revised form 28 March 2011 Accepted 9 September 2011
a b s t r a c t The rapid advancement of numerical techniques and the availability of increasingly powerful supercomputers recently enabled scientists to use large eddy simulation (LES) to simulate numerically the flow in a full subassembly composed of wire-wrapped pins. Because of the extreme computational cost of such simulations, it was not possible to conduct a sensitivity case on the pin-wire interface modeling. Since such calculations are likely to be extended to conjugate heat transfer, however, a sensitivity study is necessary to assess the reliability of the numerical results. It is well known, for example, that conjugate heat-transfer results are often strongly influenced by near-wall modeling. The objective of the present work is to investigate the effect of pin-wire contact modeling from the point of view of both the hydraulics and the heat transfer characteristics. In particular, the focus is on the prediction of the hot spot in conjugate heat-transfer calculations. The primary test case is the simplified geometry recently proposed by Ranjan et al., which consists of a simple channel flow with a wire embedded in one of the walls. After reproducing the results using the LES code Nek5000, we examined several other choices for the wire-pin interface modeling, including the introduction of a nominal gap between the wire and wall. The results shed light on the sensitivity of CFD calculations results to the modeling of the interface region between wires and pins. © 2011 Elsevier B.V. All rights reserved.
1. Introduction The study of conjugate heat transfer is needed to support accurate simulations of the flow and temperature fields of a nuclear reactor core. With computational fluid dynamics-based techniques and the use of supercomputing platforms, one can now achieve a remarkable degree of accuracy for the most relevant turbulence factors, such as the average velocity field, the averaged temperature distribution, and the stresses inside the fluid region. The coupling of the fluid dynamics problem with the conduction equation solved inside the fuel pin is less straightforward, however. In particular, the choices made in modeling the contact between the fuel pin and other structural components (such as spacers or wires in the case of sodium-cooled fast reactors) might lead to
夽 The submitted manuscript has been created by UChicago Argonne, LLC, Operator of Argonne National Laboratory (“Argonne”). Argonne, a U.S. Department of Energy Office of Science laboratory, is operated under Contract No. DE-AC02-06CH11357. The U.S. Government retains for itself, and others acting on its behalf, a paid-up nonexclusive, irrevocable worldwide license in said article to reproduce, prepare derivative works, distribute copies to the public, and perform publicly and display publicly, by or on behalf of the Government. ∗ Corresponding author. Tel.: +1 6302526637. E-mail address: [email protected] (E. Merzari). 0029-5493/$ – see front matter © 2011 Elsevier B.V. All rights reserved. doi:10.1016/j.nucengdes.2011.09.030
sensitive perturbations on the prediction of the peak temperature in the solid (hot spot), the peak temperature in the fluid or other thermodynamic and hydrodynamic properties, since such details have effects on other geometries (Qu and Mudawar, 2002). For instance, Pointer et al. (2009) recently used large eddy simulation (LES) to simulate numerically the flow in a full subassembly composed of wire-wrapped pins. Because of the extreme computational cost of such simulations, it was not possible to conduct a sensitivity study of the pin-wire interface modeling. Since such calculations are likely to be extended to conjugate heat transfer, however, a sensitivity study is necessary in order to assess the reliability of the numerical results. In fact, the wire introduces a higher resistance to conduction in the radial direction which, combined with the high conductivity of sodium, might lead to a significant localized hot spot. The presence of local stagnation points might also affect the temperature distribution in the fluid (Tang et al., 1978). The way the interface is modeled has the potential to affect the value of the hot spot, since it would affect the hydrodynamics and the local heat transfer. Moreover, in some instances, the wire might be in contact with the pin, depending on design and operational conditions. Thus, to the uncertainty due to numerical modeling one should add uncertainty due to the variability of the physical conditions. The importance of these points is shown clearly in Fig. 1 which presents the results
E. Merzari et al. / Nuclear Engineering and Design 253 (2012) 374–386
Nomenclature y+ xi x, y, z ui u
v w T T˜ T˜max,f
␥ Rez Rex Wbulk Ubulk D h g e
r 2␣
normalized wall distance Cartesian coordinates (vector notation) Cartesian coordinates Cartesian velocity components Cartesian velocity component in direction, x Cartesian velocity component in direction, y Cartesian velocity component in direction, z temperature normalized temperature normalized maximum temperature in the fluid linear gradient coefficient density kinematic viscosity Reynolds number in direction z (defined in the text) Reynolds number in direction x (defined in the text) bulk velocity in direction, z bulk velocity in direction, x diameter of the wire half width of the channel gap size fictitious parameter representing the amount of contact. For cases with no contact it is equal to g. For cases with contact it is equal to h/2(1 − cos(˛/2)) distance from the bottom wall contact angle
of a conjugate heat transfer calculation performed for conditions similar to the Advanced Burner Test Reactor (Chang et al., 2006). Details of the calculation, performed with STAR-CCM+ 5.02 (STARCCM+ 5.02 Users Guide) are reported in Table 1. To simplify the problem periodic boundary conditions are assumed on the cross section (Fig. 1) while inlet–outlet boundary conditions are assumed in the streamwise direction (the whole core length was simulated, with uniform axial heating). The results shown are mesh converged (an example of the coarsest mesh used is reported in Fig. 1b). Fig. 2a and b shows a snapshot of the temperature distribution in an axial cross-section. The temperature peak is located downstream of the wire near the contact between wire and pin. These results are consistent with Gajapathy et al. (2007). Fig. 2c shows the difference between the temperature distribution and the linear gradient (Fig. 3) associated with uniform axial heating. The perturbations on the linear gradient are due to the presence of the wire. The main purpose of the present work is to shed some light on how pin-wire contact modeling affects the temperature prediction in a simplified geometry. Arguably, such a study would be problematic in an actual core because the hot-spot position and value are determined essentially by the axial position (due to the nonuniformity of the power distribution; see Khakim et al., 2009, and Baglietto and Ninokata, 2007) and it would be more difficult to separate the small but potentially significant effect due to the modeling of the wire or other structural components (Fig. 2c). To investigate accurately the influence of such effects it would be necessary to perform several full pin-length calculations with highfidelity models for each given power shape configuration, which would require a substantial computational effort. Another approach is to assume that the effects due to the non-uniformity of the power shape and the wire-pin contact modeling are separable. In such case a single full-length calculation could be performed assuming the correct power distribution and a pre-set contact modeling. Then a set of computationally much less expensive calculations can be performed assuming a uniform power distribution, which allows for the use of periodic boundary
375
conditions, to investigate the effect of the wire contact-modeling. By properly rescaling the temperature differences due to the contact modeling with the power distribution it is possible to determine how such temperature differences vary with the location in the reactor core (Section 3.2). In this paper we will provide an example of how such methodology would be used. The primary test case is the geometry proposed by Ranjan et al. (2010), which consists of a simple channel flow with a wire embedded in one of the walls. Ranjan et al. provided detailed and extremely accurate distributions for the velocity field and other turbulent statistics (i.e., obtained with highly resolved direct numerical simulation techniques), thus allowing us to fully verify and validate our calculations. After reproducing the results using the LES code Nek5000 (Fischer et al., 2008), we examined the conjugate heat-transfer problem (which was not studied by Ranjan et al. (2010)). We also considered an additional modeling option for the wire: the introduction of a nominal gap between the wire and wall. We then used the commercial code STAR-CCM+ in RANS mode to study the same problem, with the double purpose of verifying the results obtained and testing the capability of RANS-based methods to achieve accurate results for the present case. As a last step we investigate in greater detail, using STAR-CCM+, the influence of the contact modeling with a higher number of cases, and discuss the relevance of these results for a full pin-wire conjugate heat transfer calculation. 2. Methodology In this section we describe the wall-wire contact models used. We present the mathematical model, the relevant parameters chosen, and the methodology used for the simulations. 2.1. Problem description The direct numerical simulation (DNS) calculation of Ranjan et al. (2010) was selected for its generality, its high resolution (it was computing over 100,000,000 collocation points), and its relatively low Reynolds number. The geometry investigated in that paper corresponds to a wire immersed in a channel flow with both streamwise (direction z) and spanwise (direction x) periodic Table 1 Conditions for the preliminary case. Description
Value
Total power Number of active core assemblies Number of pin per assembly Power profile Pin diameter Pin pitch Pin material
250 MWth 54 217 Flat 8 mm 9.05 mm Metallic uranium (gap and cladding are omitted for simplicity) 1 mm 0.15 m 7.5 0.75 m Sodium 5.8 m/s 600 K k–ε realizable, two-layer 23 Inlet–outlet
Wire diameter Wire pitch Contact angle Pin length in the core region Coolant Inlet velocity Inlet temperature Turbulence model Maximum wall y+ Boundary conditions (streamwise direction) Buoyancy Discretization Number of meshes for convergence
No Second order About 8,000,000 polyhedral with prismatic layer
376
E. Merzari et al. / Nuclear Engineering and Design 253 (2012) 374–386
Fig. 1. Preliminary case. Boundary conditions (a) and coarse mesh for the fluid (b) and the entirety of the domain (c).
boundary conditions (the boundary conditions are identical for Nek5000 and STAR-CCM+). The geometry is shown in further detail in Fig. 4. In the present work, while keeping the same outline for the geometry, we explore three LES wall-wire contact models. In the first case we keep the same contact modeling used by Ranjan et al. (2010), where the wire and wall are considered in contact and merged with an angle of 7.5◦ . In the second case the wire and wall are separated by a nominal gap g. For the present case Lx = 4h, Lz = 8h, and D = h. The nominal gap is set to equal to g = 0.05 h, in order to allow a small but significant cross-flow. As a last step a simulation was performed with a gap g = 0.05 h between wire and wall with a baffle blocking the cross flow below the wire. No heat transfer calculation was performed in the work by Ranjan et al. (2010). In the present simulations a heat source is specified within the walls while it is assumed that no heat is generated in the wire. After the LES cases, a series of RANS cases was simulated using STAR-CCM+ 5.02 at different contact angles and gap sizes. A full list of the cases run can be found in Table 2. 2.2. Mathematical model The conjugate heat-transfer problem is solved by resolving the incompressible, constant-property, Navier–Stokes equation in the fluid (1)–(3) and the conduction equation in the solid (4): ∂ui =0 ∂xi
(1)
∂u 1 ∂p ∂2 ui ∂ui + uj i = − + ∂xi ∂t ∂xj ∂xj ∂xj
(2)
∂T ∂T ∂2 T = ˛f + uj ∂t ∂xj ∂xj ∂xj
(3)
∂T ∂2 T + q, = ˛s ∂xj ∂xj ∂t
(4)
where u is the velocity field, T is the temperature, is the viscosity, ˛ is the thermal diffusivity (s in the subscript refers to solid, f for the fluid), is the density of the fluid, p is the pressure, t is the time, and x is the Cartesian coordinate. Here q is not the volumetric heat source but rather the ratio between the volumetric heat source and the thermal capacity of the solid; it can be seen as a temperature source. This source is set to a constant value in the heated solid and to zero in the wire. Table 2 List of cases run for the present paper. Case
Code and method
Gap size
Contact angle
I II III IV V VI VII VIII
Nek5000/LES Nek5000/LES Nek5000/LES STAR-CCM+/RANS STAR-CCM+/RANS STAR-CCM+/RANS STAR-CCM+/RANS STAR-CCM+/RANS
0.05 D 0.05 D (with baffle) – – – – 0.025 D 0.05 D
– – 7.5 15 7.5 3.25 – –
E. Merzari et al. / Nuclear Engineering and Design 253 (2012) 374–386
377
Fig. 2. Preliminary case, Temperature in the cross section (a) with detail (b) and three-dimensional view of the wall surface temperature difference (difference between the temperature and the mean fluid temperature at each height) (c). Contour plot velocity components in the same plane of (a) velocity in direction y (d) and velocity in direction z (e).
Unlike in the case presented in Figs. 1–3, periodic boundary conditions are assumed in the streamwise and spanwise directions for the velocity, with a constant pressure gradient in both directions determined dynamically in order to maintain a constant
mass-flow rate. Similar configurations are common in the simulation of parallel channel flows. In the case of conjugate heat transfer, special attention needs to be paid to the energy equations, since the temperature field is not periodic. If the temperature source
378
E. Merzari et al. / Nuclear Engineering and Design 253 (2012) 374–386
Fig. 3. Temperature as a function of z at the center of the computational domain.
distribution is constant, however, the problem can be reduced to a periodic one for the temperature as well (Iaccarino et al., 2002). The temperature in a configuration such as the one examined here can be written as: T = T0 + z + T˜
(5)
which represents a linear gradient in direction z (streamwise direction) superposed on a periodic temperature field T˜ . The ramp of the gradient can be determined by the conservation of energy principle (since the energy is evacuated from the system only through the fluid, if higher-order axial conduction is neglected (Ghaddar et al., 1986): =q
cp,s AS 1 cp,f Af Wbulk
(6)
Here, AS /Af is the ratio between the area of the solid and the fluid (on a plane normal to z), and Wbulk is the streamwise bulk velocity. We note that depends on the ratio between the heat capacity of the solid and the heat capacity of the fluid cp,s /cp,f . In fact, is a measure of how easily the fluid is able to evacuate heat from the solid. If the heat capacity of the fluid is significantly higher than the heat capacity of the solid, the same temperature source in the solid will cause a smaller, streamwise fluid temperature increase. A high value of indicates poor cooling, whereas a low value indicates effective cooling. Substituting (5) into (3) and (4) leads to the following set of equations for the temperature: ∂T˜ ∂2 T˜ ∂T˜ + uj = ˛f − u3 ∂t ∂xj ∂xj ∂xj
(7)
∂2 T˜ ∂T˜ = ˛s +q ∂t ∂xj ∂xj
(8)
The equations are balanced in the sense that the energy produced in the solid is dissipated in the fluid. No energy is transferred at the boundaries. Therefore, periodic boundary conditions can be applied for the temperature field T˜ . 2.3. Choice of parameters The relevant parameters for the following problem are the Reynolds numbers in direction x and z; the Prandtl number; the ratio between the thermal diffusivity in the solid and the thermal diffusivity in the liquid; and, to a minor extent, the ratio between the thermal capacities. The Reynolds number in direction z can be defined as: Rez =
Wbulk h ,
(9)
which is chosen to be equal to 6000 (consistently with Ranjan et al. (2010), albeit slightly higher). Wbulk is the bulk velocity in the streamwise direction. The Reynolds number in direction x can be defined as: Rex =
Ubulk h ,
(10)
which is chosen to be equal to 1000 (consistently with Case C in Ranjan et al. (2010), albeit slightly higher). Ubulk is the bulk velocity in the spanwise direction. The Prandtl number can be defined as: Pr =
, ˛f
(11)
which is chosen to be equal to 0.01, a value representative of a liquid metal coolant. The ratio between the two thermal diffusivities is
Fig. 4. Three-dimensional representation of the case with (left) contact between wire and wall and (right) no contact between meshes.
E. Merzari et al. / Nuclear Engineering and Design 253 (2012) 374–386
379
Fig. 5. Grid structure for the nek5000 simulation: (a) very coarse grid (fourth-order polynomial) for the no-contact case; (bottom) detail (b). CCM mesh (c) near the contact between wire and wall.
assumed to be equal to 10 (the thermal diffusivity in the liquid being 10 times higher than the solid), and the ratio between the thermal capacities is assumed to be equal to 1. This corresponds loosely to the case of a liquid metal coolant and a uniform metallic fuel. We note that the source term does not appear as a relevant parameter of choice. The reason is the following. There exists a normalization of the temperature T˜ (and notably L, where L is a lenghtscale of arbitrary choice) for which the source term can be removed from the equation. Since T˜ is a fluctuation on the linear streamwise gradient, it is intuitive that the choice of q is irrelevant for the determination of its shape. An increased value of q would have a significant effect on , however, thereby influencing the distribution of T.
On the upper and lower boundaries for the solid, a zero heat flux condition is applied (the temperature derivative is set to zero on the boundary). 2.4. Numerical methods The simulations were carried out by using the spectral element code Nek5000 in LES mode on two Argonne computer systems—Cosmea (a 32-node, 128-core Linux cluster) and Blue Gene/P (an IBM system with 40,960 quad-core compute nodes)—with up to 15,000,000 collocation points. An example of structure of the macroscopic grid used is shown in Fig. 5 for the streamwise-normal cross section; the figure depicts the elements and the actual collocation points for a very coarse mesh with no
380
E. Merzari et al. / Nuclear Engineering and Design 253 (2012) 374–386
Fig. 6. Instantaneous streamwise velocity in a z-normal plane (a) and time-averaged streamwise velocity in a z-normal plane (b). Spatial units are normalized by h.
contact between wire and wall. The discretization is designed to allow for at least one point near the wall at y+ < 1 and five points within y+ < 10. Spectral elements methods are a class of higher-order methods in space. Polynomial functions of up to the 11th degree have been used to discretize the velocity field in each element. In the generalized weighted residual framework, the present spectral element method can be classified as a Galerkin method where the test functions and the basis functions for each element are Lagrange polynomials evaluated on Gauss–Lobatto–Legendre collocations points (for the velocity). In LES, large-scale turbulence is simulated while smaller scales are modeled. Since smaller scales have a nearly universal behaviour, LES is a more reliable methodology than is Reynolds-averaged Navier–Stokes (RANS), in the sense that it generally depends less on the model details. It is also preferable to use in place of direct numerical simulation because it can be performed at a reasonable computational cost. The contribution of the smaller scales to the energy cascade is in the present case modeled through a local, element-based explicit cut-off filter in wave-number space. The energy is removed from the smallest simulated scales (high wavenumbers), thus mimicking the effect of smaller eddies (Fischer and Mullen, 2001). This procedure allows for a coarser grid and therefore lower computational cost. Time advancement has been carried out through an explicit third-order, backward finite-difference scheme. An additional benefit of this method is that it does not explicitly define an eddy viscosity, leaving the definition of a turbulent Prandtl number as superfluous. The Nek5000 code has been extensively validated and is massively parallel. In the present case the time stepping was chosen to have a maximum Courant number no higher than 0.4. For comparison, some calculations were also run by using the commercial code STAR-CCM+ in the RANS framework with a realizable k–ε turbulence model and all-y+ wall modeling. The number of meshes used was considerably lower for the same cases (300,000 meshes). Polyhedral cells with a prismatic layer were used everywhere except in the contact region (Fig. 5c): local refinement was applied to this region to reach convergence at a lower cost. A
segregated flow solver and a second order scheme for the convective terms were also employed. The turbulent Prandtl number was assumed constant and equal to 0.9. While this approach presents serious limits for the simulation of the flow in liquid metal systems (Tzanos, 2011), it proved to yield accurate results for simplified geometries such as pipes and channels. 3. Results and discussion In the case of LES, several cases were run corresponding to two different geometries and different mesh resolutions. While the macroscopic mesh was not changed during the simulations, mesh convergence was proven by increasing the order of the polynomial basis from the 4th order up to the 11th. For each RANS case a mesh convergence study was conducted. Among other verification tests, each result shown here has been obtained on at least two separate meshes with a polyhedral mesh size ratio equal to 1.5. 3.1. Velocity distribution We begin by comparing the velocity field with the data of Ranjan et al. (2010). Fig. 6a shows an instantaneous contour plot for the streamwise velocity, while Fig. 6b shows a contour plot for the averaged streamwise velocity; both figures refer to the case with contact between wire and wall. Compared with similar plots published by Ranjan et al. (2010), these figures show a good degree of accuracy considering that the Reynolds number is 10% higher than that for the DNS case. Fig. 7 shows four profile comparisons of the results obtained with Nek5000 (with contact modeling) and the DNS data. The accuracy is good considering the enormous difference in resolutions (we used less than one-tenth of the collocation points) and the difference in Reynolds number. All profiles refer to time-averaged data. (For the position in x, see Fig. 7). Fig. 8 shows a comparison between the Nek5000 results for the case with no contact between wire and wall and the DNS data. Such a comparison would be meaningful, for validation purposes, only
E. Merzari et al. / Nuclear Engineering and Design 253 (2012) 374–386
381
Fig. 7. Comparisons of the streamwise velocity and the cross velocity at x/h = 0.0 (wire centerline) and x/h = 1.0. The solid line represents the DNS data. All values are normalized by Wbulk . r represents the distance from the bottom wall.
Fig. 8. Comparisons of the streamwise velocity at x/h = 6. Values are normalized by Wbulk . The solid line represents the DNS data; r represents the distance from the bottom wall.
far from the wire. We therefore did the comparison at x/h = 6.0: the agreement seems reasonable and within the scope of the present work. Time-averaging convergence was verified. Near the wire, the two cases (i.e., no contact and contact between wire and wall) behave differently. In the case of no contact, the velocity distribution appears to be altered. Fig. 9 compares the two Nek5000 results at x/h = 1.0. Near the bottom wall (r/h = 0.0, where r is the distance from the bottom wall), the streamwise velocity is lower, an effect that could be related to the presence of a higher cross-flow in this area. In fact, as it can be seen in Fig. 10, the recirculation pattern near the gap is altered when compared to the case with no contact. Moreover the cross velocity in the narrow gap is significantly higher than the surrounding areas. This effect is nontrivial, and it is a relevant from the point of view of the hydrodynamics. To determine if it the cross flow in the case with no-contact is the source of the disturbance in the streamwise velocity distribution, an additional case has been simulated. A conductive baffle was introduced between wire and wall (by fixing the velocity components to zero), thus blocking the flow passage below the wire. Fig. 11 demonstrates how the streamwise velocity distribution for this case matches well the velocity distribution for the case with contact.
382
E. Merzari et al. / Nuclear Engineering and Design 253 (2012) 374–386
Let us, however, assess in further detail the magnitude of this influence
3.2. Temperature
Fig. 9. Comparison of the streamwise velocity between the two Nek5000 cases at x/h = 1. Values are normalized by Wbulk ; r represents the distance from the bottom wall. The solid line represents the case with contact between wall and wire.
It can be concluded that the wall-contact modeling has a definite impact on the streamwise velocity distribution in the region downstream of the wire. In particular, a case with no contact would introduce a small but significant cross flow between wire and pin which in turn affects the stream-wise velocity distribution by continuity and transport of momentum. Such change in the streamwise velocity distribution would have an impact on the local heat transfer coefficient and therefore on the local temperature distribution.
Fig. 12a shows the average temperature distribution for the case with contact between wire and wall while Fig. 12b shows the average temperature distribution for the case with no-contact between wire and wall. All values are normalized by D where D is the diameter of the wire. The plots in Fig. 12a and b are similar except for the wire region. A hot spot can be observed in both cases near x = 0 and r = −h (near the lower adiabatic boundary), with the hot spot in the no-contact case shifted toward x > 0. Note that the scales are different. It is evident however, that in the case with contact between wire and pin the hot spot is slightly higher. To provide a more quantitative comparison, we are also reporting the temperature profiles for x/h = 0 in Fig. 13. While the two profiles are different and the case with contact clearly predicts a higher hot spot, the difference appears to be small. We conclude that the significant difference observed in the hydrodynamic side had little effect on the hot-spot prediction for the solid. We note, however, that the temperatures of the fluid at the wall near the contact region between wire and pin is significantly affected by the contact modeling. In the case with no contact the peak T˜ in the fluid was evaluated to be −31.5 a.u., while in the case with contact it was evaluated to be −4.77 a.u. This is a significant temperature increase (about 7.5% of the maximum temperature difference within the model). The peak fluid temperature is an important design parameter because it determines the margin to boiling and other technological considerations. Before drawing any further conclusions, however, we will proceed in examining in further details the results; Figs. 14 and 15
Fig. 10. Cross velocity: (a) vector plot for the case with no contact; (b) near-gap vector plot for the case with contact; (c) contour plot of the velocity in direction x for the case with no contact.
E. Merzari et al. / Nuclear Engineering and Design 253 (2012) 374–386
Fig. 11. Comparison of the streamwise velocity between the Nek5000 cases at x/h = 1. Values are normalized by Wbulk ; r represents the distance from the bottom wall.
provide further verification of the results obtained. Fig. 14 shows the average temperature for the STAR-CCM+ case with contact between wall and wire. Fig. 15 shows a comparison between the temperature profile at x = 0 obtained by Nek5000 and that by STARCCM+ for the case with contact between wire and wall. The results are in very good agreement. The peak fluid temperature for the same case is estimated by STAR-CCM+ at −36.2 a.u., thus with
383
a discrepancy of approximately 1.5% of the maximum temperature difference. Two very different codes provided very similar results, giving us confidence on the validity of the conclusions that can be drawn. In particular it is likely that the uncertainty due to the contact modeling is more significant than the uncertainty due to the turbulence modeling. Moreover for the present case RANS appears sufficiently accurate to predict the value of the hot spot. Let us therefore examine in greater detail the impact of the wall-wire contact modeling using the cheaper RANS methodology. A series of cases have been run in STARCCM+ changing the type of contact modeling, which can be parameterized based on the distance e (Fig. 16a). A positive e represents cases with no contact modeling (and e therefore represents the distance between wire and pin). A negative e represents cases with contact: the higher the absolute value of e the higher will be the contact angle ˛. For each different case the results were verified to be mesh converged: this was particularly difficult for small values of e, where a very strong local refinement near the wire-wall interface was required (up to 9,000,000 meshes). No fillet was applied. Fig. 16b reports the peak fluid temperature as a function of e (a fictitious parameter that accounts for the amount of contact) for all the cases run in the present work. It is clear that the peak fluid temperature is not a monotonic function of e: in fact the function has a maximum between e = −0.0016 D and e = 0.025 D based on the data collected. However, it can be assumed that the peak will actually be at e = 0.0, where the wire and the wall form a singularity since the presence of a conductive layer of fluid (e > 0) would decrease the peak temperature. In fact, the baffle case—despite the hydrodynamic in all similar to the case with contact and e = −0.016 D presented a temperature drastically lower and similar to the no-contact case with e = 0.05 D: since the fluid is more conductive than the solid the fluid layer between wire and wall helps reducing the temperature in the wire region (Fig. 13). Moreover, in low Prandtl flows, such as the one considered here, conductivity plays a much stronger role
Fig. 12. Average temperature T˜ distribution for the case with no contact (a) and average temperature T˜ distribution for the case with contact (b).
384
E. Merzari et al. / Nuclear Engineering and Design 253 (2012) 374–386
Fig. 13. Temperature T˜ profile at x/h = 0.0 for the two cases. Detail of the distribution near the lower adiabatic boundary is shown in the right box. The solid line represents the case with contact between wall and wire.
Fig. 14. Temperature T˜ for the case with contact, obtained by using STAR-CCM+.
in heat-transfer than in water and other medium to high Prandtl fluids. The fact that the peak fluid temperature has a peak near e = 0.0 D warrants caution in the estimation of peak temperature for a real bundle case. From Fig. 16a it appears evident that, in order to be conservative a simulation should assume contact, with a sufficiently small contact angle. Therefore a reasonable uncertainty due to the contact modeling should be assumed in the range of 2.5–5% of the maximum temperature difference (all of this is, of course, valid only for the present set of parameters). Let us now go back to the preliminary case examined in the introduction (Figs. 1–3, Table 1). The maximum temperature difference within the whole model at each given cross section varied between 80 K and 90 K (Fig. 17). In the fluid, the variation was considerably lower (with a variation consistently around 13 K at the wall, when the linear gradient is subtracted, Fig. 2c). Normalizing by ␥D where D is the diameter of the wire, we obtain a spread of 380 a.u. thus consistent with the values reported earlier (368 a.u. Fig. 12a). A 10 a.u. error band as suggested in Fig. 16b would be equivalent to a ±2.1 K band in the pin-wire conjugate heat transfer simulation. However, if no-contact cases are considered a larger
Fig. 15. Temperature T˜ profile at x/h = 0.0 for the case with contact, obtained by using Nek5000 and STAR-CCM+. The solid line represents the case performed using Nek5000.
E. Merzari et al. / Nuclear Engineering and Design 253 (2012) 374–386
385
Fig. 16. Peak Temperature T˜ for the cases considered in this work as a function of e the straight dashed lines represent the level of uncertainty recommended for the present case (a plus or minus 10 a.u. band around the value predicted with a reasonable contact angle).
Fig. 17. Temperature distribution in the solid for the preliminary calculation (details in Table 1).
band should be assumed (up to 8.2 K for the calculations performed here, to make sure to include the zero point in Fig. 16b). 4. Conclusions A series of LES and STAR-CCM+ simulations have been carried out on a simplified geometry in order to investigate the influence of wire-pin contact modeling on the prediction of the temperature distribution.
The results presented here suggest that the modeling of the contact between wire and wall is significant, if not crucial, from the hydrodynamic point of view. It is less so from the point of view of the hot-spot prediction, since a significant difference in hydrodynamic behaviour leads to an approximated one percentage point difference in the value of the hot spot. The peak temperature of the fluid, on the other hand, appears to be significantly affected. This warrants caution in the use of CFD for the prediction of the peak surface temperature, which should be carried out assuming a small contact angle (as small as numerically permitted). These results have then been related to a full single pin-wire conjugate heat transfer simulation with periodic boundary conditions. The simulation, albeit simplified, provided a reference and allowed us to give a quantitative value to the uncertainty in the pinwire contact modeling: for a case similar to the ABTR, depending on the simulation performed such uncertainty has been estimated between 2.1 K for a small contact angle simulation and up to 8.2 K for a no-contact simulation. While this temperature differences might appear small it is important to point out that they might increase substantially when a more realistic simulation is taken into account (the presence of cladding and a gap for instance) and constitute a significant percentage of the temperature difference in the cross section. An additional conclusion that can be drawn is that RANS-based techniques are sufficient to reach accurate results for the prediction of the hot spot and peak fluid temperature for the present geometry and parameters. This seems to be due to the conduction-dominated nature of heat transfer for low-Prandtl fluid conjugate heat transfer cases. Acknowledgments This work was supported in part by the Office of Advanced Scientific Computing Research, Office of Science, U.S. Department of Energy, under Contract DE-AC02-06CH11357.
386
E. Merzari et al. / Nuclear Engineering and Design 253 (2012) 374–386
References Baglietto, E., Ninokata, H., 2007. “Improved turbulence modeling for performance evaluation of novel fuel designs”. Nuclear Technology 158 (2), 237–248. Chang Y. I., Finck, P.J., Grandy C., 2006. “Advanced burner test reactor: preconceptual design report”, Argonne National Laboratory report ANL-ABR-1 (ANL-AFCI-173). Fischer, P.F., Lottes, J., Pointer, W.D., Siegel, A., 2008. Petascale algorithms for reactor hydrodynamics. Journal of Physics: Conference Series 125, 012076. Fischer, P.F., Mullen, J.S., 2001. Filter-based stabilization of spectral element methods. In: Comptes Rendus de l’Académie des sciences Paris, 332, Série I. Analyse numérique, pp. 265–270. Gajapathy, R., Velusamy, K., Selvaraj, P., Chellapandi, P., Chetal, S.C., 2007. CFD investigation of helical wire-wrapped 7-pin fuel bundle and the challenges in modeling full scale 217 pin bundle. Nuclear Engineering and Design 237, 2332–2342. Ghaddar, N.K., Karniadakis, G.E., Patera, A.T., 1986. A conservative isoparametric spectral element method for forced convection; application to fully developed flow in periodic geometries. Numerical Heat Transfer B: Fundamentals 9 (3), 277–300.
Iaccarino, G., Ooi, A., Durbin, P.A., Behnia, M., 2002. Conjugate heat transfer predictions in two-dimensional ribbed passages. International Jounal of Heat and Fluid Flow 23, 340–345. Khakim, A., Merzari, E., Ninokata, H., 2009. Feasibility study of the application of exotic pin for tight-lattice fuel assembly. Transactions of the American Nuclear Society 100, 747. Pointer, W.D., Smith, J.G., Fischer, P., Siegel, A., 2009. RANS simulations of turbulent diffusion in wire-wrapped sodium fast reactor fuel assemblies. In: abstract in International Conference on Fast Reactors and Related Fuel Cycles, Kyoto, December, pp. 249–250. Qu, W., Mudawar, Issam, 2002. Analysis of three-dimensional heat transfer in micro-channel heat sinks. International Journal of Heat and Mass Transfer 45, 3973–3985. Ranjan, R., Pantano, C., Fischer, P., 2010. Direct simulation of turbulent swept flow over a wire in a channel. Journal of Fluid Mechanics 651, 165–209. STAR-CCM+ 5.02 User Guide. Tang, Y.S., Coffield Jr., R.D., Markley Jr., R.A., 1978. Thermal analysis of liquid-metal fast breeder reactors. American Nuclear Society, ISBN: 0-89448-011-1. Tzanos, C., 2011. Improved simulation of heat transfer in liquid-metal flows. Nuclear Technology 174, 41.