Large-Eddy Simulation of convective wall-boiling flow along an idealized PWR rod bundle

Nuclear Engineering and Design xxx (2016) xxx–xxx

Contents lists available at ScienceDirect

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

Large-Eddy Simulation of convective wall-boiling flow along an idealized PWR rod bundle Djamel Lakehal ⇑, Daniel Caviezel ASCOMP AG, Switzerland

a r t i c l e

i n f o

Article history: Received 31 May 2016 Received in revised form 1 November 2016 Accepted 19 December 2016 Available online xxxx

a b s t r a c t A LES analysis of turbulent convective boiling flow along the rods of a PWR sub-channel was performed using the code TransAT. The campaign is an extension of a former attempt in which the flow was predicted by means of DNS (or super resolved LES). The extension refers to accounting for wall-boiling heat transfer using the triple flux decomposition model, employed within the LES approach constructed on the basis of the filtered mixture equations. The results in the wall boiling configuration differ from the singleflow case in many ways, but more particularly: the flow seems to be portioned in three zones, (i) preboiling, (ii) transitional boiling, and (iii) post-boiling, where turbulence activity increases with distance from the onset of nucleate boiling location, affecting heavily the dynamic and thermal boundary layers. The boundary layers thin quite substantially in the post-boiling section, and the flow accelerates due to higher volumetric flow rate induced by phase change, notably in the gap region. The vorticity distribution between the single flow and boiling flow cases suggests that boiling causes a higher rate of vorticity generation. These results corroborate with the other finding, that is: the average frictional velocity in the boiling case increases by a factor of 1.5. Scrutinizing the power density spectra reveals that the boiling flow embodies about twice more energy than the single phase flow, with marked large-scale energetic eddies at the upper end of the rods, i.e. z/L > 0.8. A clear transitional behavior is observed in the boiling flow PSD, in that the amount of energy embodied in the flow increases with vapor production. Ó 2016 Elsevier B.V. All rights reserved.

1. Introduction The onset of nucleate boiling on the rod surface occurs when the temperature of the wall slightly exceeds saturation (Incropera and DeWitt, 1996). The small vapor bubbles form and stay attached to the solid wall. Past the point of net vapor generation, the bubbles detach and remain trapped within a layer relatively close to the wall, beyond which – towards the core flow – they condense. This is the simplest picture one can draw, but the reality is way more complex than that. Under turbulent conditions indeed, the flow is affected by large and small eddies, impacting in turn the rate of wall-to-core-flow heat transfer and thus phase change, both near the wall (boiling) and far in the core flow (condensation) when the liquid is subcooled. Without listing all the features associated with turbulent flow in narrow gaps of sub-channels, it is perhaps useful to evoke the most important ones which the authors believe as key issues in modelling using

⇑ Corresponding author. E-mail address: [email protected] (D. Lakehal).

mainstream CFD. The reader may refer to the review paper of Rehme (1992) compiling early research findings on the subject. Early experiments (Rowe et al., 1974) revealed the existence of macroscopic pulsating flow structures (not necessarily turbulence) in the regions adjacent to the gaps, with strong implications on the mixing between adjacent sub-channels. Other important features were identified experimentally as marked characteristics of flow along rod bundles (Seale, 1982; Vonka, 1988), including secondary flow motion and large-scale turbulent motions enhancing heat extraction from the heated wall. Clearly, there are incentives to resort to 3D CFD/CMFD for the prediction of the detailed fluid flow and temperature distribution in rod bundles for safety issues and operational reliability of the fuel elements. Subchannel analysis ignores the fine structures of velocity and temperature distributions in the flow passages, and can thus not account for mixing effects caused by the presence of spacers or other geometrical disturbances. Large-scale turbulent motion and larger periodic pulsating structures responsible for mixing between sub-channels are out of reach of steady state approaches resorting to statistical averaged models. Further, secondary flow motion can be predicted using anisotropic models or

http://dx.doi.org/10.1016/j.nucengdes.2016.12.024 0029-5493/Ó 2016 Elsevier B.V. All rights reserved.

Please cite this article in press as: Lakehal, D., Caviezel, D. Large-Eddy Simulation of convective wall-boiling flow along an idealized PWR rod bundle. Nucl. Eng. Des. (2016), http://dx.doi.org/10.1016/j.nucengdes.2016.12.024

2

D. Lakehal, D. Caviezel / Nuclear Engineering and Design xxx (2016) xxx–xxx

full Reynolds stress models only. While the velocity magnitude associated with secondary flow moving within the elementary cells of the sub-channels is small compared to the axial one, its implication on heat transfer is important, in particular near the wall, which appeals for the use of wall-resolved strategies. We present the results of a moderately-resolved LES of turbulent convective boiling flow upward along the heated rods of an idealized PWR sub-channel. The campaign is an extension of a former attempt in which the non-boiling flow was first predicted by means of DNS (Caviezel et al., 2013) using 6.2 million cells. The main objective is to address the feasibility of the boiling model employed here in the LES context, and whether the model combination reveals any important information about the flow that statistical approaches are not capable to reach. As a first step though, the study assesses two different subgrid-scale models for the prediction of the simple convective turbulent flow along the rod, prior to their application to the boiling case. The problem is inspired by the PSBT OECD single sub-channel benchmark (Rubin et al., 2010), in terms of radial dimensions in particular, albeit the deliverables are different from the actual OECD case. The operating conditions selected here are made on purpose different from PSBT, namely the power, which has been adjusted according to the reduced length (1 m instead of 3.65 m). We use for the purpose LES based on the filtered mixture model equations, combined with the triple heat-flux decomposition (the so-called RPI model) for wall-boiling. To avoid having to deal with issues associated with filtering the mixture model equations, the model has been employed under its homogeneous form, ignoring the interphase slip. The boiling model has been used as any other approach to triggering bubble creation and departure from the wall. 2. Simulation and modelling strategies 2.1. On TransAT The CFD/CMFD code TransATÓ developed at ASCOMP is a multiscale, multi-physics, conservative finite-volume solver for singlephase and multi-fluid Navier-Stokes equations (TransAT, 2015). The grid arrangement is collocated, and the solver is pressure based (Projection Type). High-order schemes are employed; up to 3rd-order schemes in space and 2nd-order in time for implicit time marching, respectively. Turbulent flows can be treated in two ways: RANS statistical models and Scale Resolving Approaches including LES and its V-LES variant for very high-Re flows. LES is built in TransAT, with dedicated routines for pressure coupling, boundary conditions, diffusive fluxes and near-wall stress integration. The code proved very efficient for this class of simulation strategy, including for DNS as well (Chatzikyriakou et al., 2015). Various schemes are available for pressure-velocity coupling, including Simple and Simple-C used in the LES context. Pressure solvers range from the SIP scheme to more sophisticated variants, with or without pre-conditioners, including a fully integrated AMG solver with various cycle possibilities. 2.2. The filtered Navier-Stokes (NS) equations In LES the motion of the super-grid turbulent eddies is directly captured whereas the effect of the smaller scale eddies is modelled or represented statistically by means of simple models, the same way as in Reynolds-averaged models (RANS). The usual practice is to model the sub-grid stress tensor by an Eddy viscosity model. In the present work, use was made two different subgrid-scale (SGS) model both in the single-phase and boiling flow case. LES involves the use of a spatial filtering operation


tÞ ¼ Fðx;

Z

1

Fðx0 ; tÞGðx  x0 Þdx0 ; with GðxÞ Z 1 GðxÞdx ¼ 1; ¼ GðxÞ; and 1

1

where the fluctuation of any variable F from its filtered value is 0 denoted by f ¼ F
 F . Filter function Gðx  x0 Þ is invariant in time and space, and is localized. Applying the filtering operation to the instantaneous Navier-Stokes equations under incompressible flow conditions leads to the system of filtered transport equations for turbulent convective flow:


¼0 ru @
Þ þ r  ðqu
u
Þ ¼ r  ðp
 sÞ þ F b þ F c ð q u @t

ð1Þ

where u is the velocity vector, q is the density, p ¼ ðP:I þ rÞ is the Cauchy stress embedding pressure and viscous terms P and r. The source terms in the RHS of the momentum equation represent the body force, Fb, and the convolution-induced residual, Fc. Further, fil
u
Þ, of tering introduces the SGS stress tensor defined as s ¼ qðuu  u which the deviatoric part is to be statistically modelled. This way, turbulent scales larger than the filter width imposed by function Gðx  x0 Þ are directly solved, whereas the diffusive effects of the SGS scales are modelled. 2.3. The filtered dispersed gas-liquid flow NS equations LES of dispersed gas-liquid flows (sometimes known as LESS, short for Large-Eddy & Structures Simulation) invokes filtering of the phase averaged gas and liquid equations. This combination is expected to return the interaction of the dispersed phase with turbulence without resorting to statistical modelling. LESS is well suited for bubbly flows in particular. The derivation of the filtered mixture equations follows from the filtered phase-average twofluid equations introduced first by Lakehal et al. (2002). These equations were derived from the instantaneous phasic equations, marked by the phase-indicator (distribution) function, ðCðxÞ ¼ 1 if x 2 L; and ¼ 0 if x 2 GÞ; and using the Component Weighted Volume Averaging, where filtered variables are defined as follows (under incompressible flow conditions):

C k ðx; tÞF k ðx; tÞ 00 e Fe k ðx; tÞ ¼ ; f ¼ F k  Fk ak ðx; tÞ

ð2Þ

where ak ðx; tÞ ¼ C k ðx; tÞ is the filtered void fraction. The application of the filtering operation to the instantaneous gas and liquid equations yields:

qk ak Þ þ r  ðqk ak ue k Þ ¼ Rm;k qk ak ue k Þ þ r  ðqk ak ue k u~ k Þ ¼ r  ak ðpk  sk Þ þ Mk þ qk ak g þ Rc;k

@ ð @t @ ð @t

ð3Þ where sk the subgrid stress tensor, Mk ¼ pk :nk dðx  xI Þ is the filtered interfacial force marked by the Dirac function dðx  xI Þ, and the normal vector to the interface nk , and Rm;k and Rc;k are the filtering-induced residual terms, which can be neglected for now. To proceed, we start from the concept of mixture phase averaged introduced by Manninen et al. (1996). The filtered mixture equations follow from summing up the filtered two-fluid mass and momentum equations (3) over phases k. Assuming that the bubbles are small and well mixed with the carrier phase, the drift/slip between the phases (uG ¼ uL ¼ um ) can be neglected and the final form of the filtered mixture equations under homogeneous form read:

Please cite this article in press as: Lakehal, D., Caviezel, D. Large-Eddy Simulation of convective wall-boiling flow along an idealized PWR rod bundle. Nucl. Eng. Des. (2016), http://dx.doi.org/10.1016/j.nucengdes.2016.12.024

3

D. Lakehal, D. Caviezel / Nuclear Engineering and Design xxx (2016) xxx–xxx

qm Þ þ r  ðqm ue m Þ ¼ 0 qk ak Þ þ r  ðqk ak ue m Þ ¼ 0 X @ e m Þ þ r  ðqm u emu e m Þ ¼ r  ðpm  sm Þ þ ðqm u M k þ qm g @t @ ð @t @ ð @t

2.5. Wall-boiling model

ð4Þ

k



where we define the filtered mixture quantities as follows:

em ¼ u

sm ¼

1 X

qm

ak qk ue k ; qm ¼

X

X

k

k

ak sk ¼ 

X

X

k

k

ak qk ; pm ¼

ak qk ðu00 u00 Þk

qf ¼ A1 St p ql C Pl ulp ðT w  T l Þ

ak pk ; ð5Þ

Note that in this simplified homogeneous formulation, all the terms involving the diffusion interphase relative velocity do not appear. Further, in case of phase change, the filtered mass conservation equation and energy equation should involve a corresponding term, which needs to be modelled as well. The interfacial term P k M k should rigorously represent the surface tension effects on the mixture, but since the dispersed phase cannot be resolved in the phase-average context, the term should be neglected. The SGS stress tensor of the mixture sm can thus be defined assuming that each phasic field is a superposition of a filtered and a nonresolved part, the effect of which on momentum diffusion needs to be modelled.

2.4. Subgrid-scale modelling In LES the super-grid part of the flow with its turbulent fluctuating content is directly predicted whereas the sub-grid scale (SGS) part is modelled, assuming that these scales are more homogeneous and universal in behavior. For turbulent flows featuring a clear inertial subrange the modelling of the SGS terms in the statistical sense could thus safely borrow ideas from the RANS context, in particular use of the zero-equation model to mimic the momentum diffusive effects. Use is generally made of the Eddy Viscosity Concept, linking linearly the SGS Eddy viscosity and thermal diffusivity to the filtered rate of strain tensor Sij and filtered temperature gradients, respectively:

sij ¼ 2lsgs Sij þ 13 dij sll ; lsgs ¼ ðCsDÞ2 qjSj2 l @ T
q00j ¼ ah @x ; ah ¼ Prsgst j

According to the standard wall boiling model, the wall heat flux qw is partitioned into three parts: qw ¼ qf þ qq þ qe . The first part is the single-phase heat transfer (convective heat flux):

ð6Þ

where Prt is the turbulent Prandtl number. The model constant (Cs) is either fixed or made dependent on the flow; this latter option is precisely the spirit of the Dynamic Smagorinsky model, short for DSM (Moin et al., 1991). A damping function is often introduced for the model constant Cs to accommodate the asymptotic behavior of near-wall turbulence, i.e. lsgs / yþ3 . Similarly, the same strategy could be used to close the turbulent SGS heat flux, where the thermal diffusivity could be determined either based on the resolved thermal-flow field, or alternatively based on the Eddy viscosity (defined dynamically) and a fixed Prt. In the present simulations, we have tested two variants based on fixed and variable model coefficient Cs, namely the dynamic DSM model with the variable Eddy diffusivity variant and the WALE SGS model (Nicoud and Ducros, 1999), the latter has been shown to behave very well in wall-bounded flows, and to be less dissipative, capable to capture the thin-shear layer accurately. The DSM SGS model was already validated and tested in the context of turbulent shear bubbly flow by Lakehal et al. (2002). Note that the subgrid bubble-induced diffusion part introduced by Lakehal et al. (2002) does not apply in the homogeneous mixture context, since it is tied to the relative velocity between the phases.

 W ; A1 ¼ 1  A2 m2 K

ð7Þ

where A1 and A2 represent the fraction of the wall surface influenced by liquid and by vapor bubbles formed on the wall, respectively, Tl and ulp are the liquid temperature and velocity at the cell adjacent to the wall, CPl is the liquid heat capacity, and Stp is the Stanton number Stp ¼ Nu=ðRel Prl Þ where Nu is the Nusselt number (Kurul and Podowski, 1990). The second part of the wall heat flux is attributed to the quenching mechanism and is defined by:

    W 2 qq ¼ A2 aq ðT w  T l Þ 2 ; A2 ¼ min 0:25  p  dw  n; 1 m K

ð8Þ

where aq is the associated heat transfer coefficient, dw is the bubble departure diameter calculated using Ünal’s (1976) correlation, and n is the active nucleation site density determined using the correlation of Eddington and Kenning (1979). The portion of heat flux due _ e  hfg ), in which hfg denotes the latent to liquid evaporation (qe ¼ m heat of phase change, involves one additional unknown besides the active nucleation site density, namely the bubble detachment frequency ‘f’, which needs to be modelled as well:

me ¼

p  d3w 6



qv  f  n

kg m2 s

 ð9Þ

The quenching heat transfer coefficient aq is dependent on the waiting or relaxation time tw between the bubble departure and the next bubble formation period (tw = 1/f):

aq ¼ 2  kl  f 

rffiffiffiffiffiffiffiffiffiffiffi  tw W p  al m2 K

ð10Þ

The cubic dependence in Eq. (9) suggests that small uncertainties on the bubble departure diameter are greatly magnified in the heat transfer model, thus deteriorating the accuracy of the overall CFD simulation. Therefore, using more accurate bubble departure models are key to the successful prediction of subcooled flow boiling heat transfer. The model as built in TransAT has been validated for different 2D and 3D problems (Thomas et al., 2013, 2015) used in the RANS context. Further, our model accounts for wall superheat required for the onset of subcooled boiling by reference to Bergles and Rohsenow (1964),

ðT Wall  T SAT ÞONB

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ! u u 8rq00 T SAT t ¼ kL hfg qg

ð11Þ

increasing the effective saturation (or the onset of nucleate boiling) temperature from 618° K to 625° K in the problem treated here. 3. Problem description The problem is inspired by the PSBT OECD sub-channel benchmark (Rubin et al., 2010) in terms of radial dimensions. The PWR fuel assembly consists of a rod bundle with water coolant flowing upward along the rods at a high Reynolds number. The rod diameter is 9.5 mm, the rod pitch is 12.6 mm and the active fuel length is typically 3.7 m. The hydraulic diameter for a unit cell is De = 11.8 mm. The coolant pressure is 15.5 MPa with temperature ranging from 290 °C to 340 °C. Representative values of fluid prop-

Please cite this article in press as: Lakehal, D., Caviezel, D. Large-Eddy Simulation of convective wall-boiling flow along an idealized PWR rod bundle. Nucl. Eng. Des. (2016), http://dx.doi.org/10.1016/j.nucengdes.2016.12.024

4

D. Lakehal, D. Caviezel / Nuclear Engineering and Design xxx (2016) xxx–xxx

erties are density: q = 710 kg/m3; surface tension coefficient: r = 0.0048 N/m; dynamic viscosity: l = 9.105 Pas; thermal conductivity: k = 0.54 W/m-K; heat capacity: Cp = 5.9 kJ/kg-K. The mass flux is G  3700 kg/m2s, corresponding to Re = GDe/l  4.8.105 and Res = us De/2m  104, where De is the hydraulic diameter. Several simplifications were adopted in the non-boiling case: First, the temperature of nucleation is assumed to be equal to the saturation temperature; second, a configuration consisting of the flow along a single sub-channel with periodic conditions is assumed; third, since the shear Reynolds number for the PWR channel is too high, performing DNS for such flow conditions will be computationally prohibitive. As such, the problem has been scaled down to more reasonable conditions, i.e. Res = 300. The length of the domain was reduced from the PWR case to relax the meshing requirements in the axial direction (1 m instead of 3.65 m), and thus the operating conditions such as the power and mass flow rate, which have been adjusted according to the reduced length. Since the distance to the onset of nucleate boiling depends on the integrated power (heat flux times rod surface area) supplied to the fluid, the heat flux was scaled accordingly. The PSBT reference and reduced (used here) operating conditions are summarized in Table 1. Fluid properties assumed independent of temperature and are the same as in the actual case. The heat flux is to be applied at the rod outer surface. Prior to this work, a similar large-scale simulation campaign of this flow (without phase change) was undertaken by the group within the DOEsponsored CASL project, where the flow was predicted using DNS (or super resolved LES) using up to 6.2 million cells. The results of this study can be found in Caviezel and Lakehal (2013) and Caviezel et al. (2013), who compared their results with the DNS data of pipe flow of Eggels et al. (1994). The simulations presented here were obtained on a coarser grid though for the same flow conditions, including the turbulent Reynolds number.

use the wall functions, though without a need for near-wall damping. The grid was distributed over 32 parallel cores.

4.2. Initial conditions & simulation parameters The initial flow conditions for the single-phase flow case were generated from an earlier flow solution obtained using cyclic inflow-outflow in a shorter 2pDe long domain. During this process, to speed up the turbulent flow generation, various grids of different refinements (coarse, medium and refined) were employed in a sequential way: the solution obtained on the coarse grid is mapped into the medium one, the solution of which is then passed to the fine mesh, which at the end is transferred to the final run as an initial/inflow/outflow condition. This process allows generating the fluctuating field to the finest mesh rather fast. In the boilingflow case, the initial flow conditions were generated from the earlier ‘fully converged’ single-phase flow solution. The wall-boiling model described above was activated. The DSM and WALE SGS models were used for LES. Flow averaging started around 20,000 time steps, with 100,000 additional steps, depending on which simulation and SGS model, to infer ergodic conditions, achieving up to ten flow-through times (t/tL).

5. Single-phase flow simulations results We first present the LES results of the convective single-phase flow using the reduced-model setup and flow conditions as discussed above, without wall boiling model, and using the two SGS models, separately. We provide the results of macroscopic flow structures, energy spectra and 1st order time-averaged quantities to compare later on with the convective boiling case.

5.1. SGS model comparison

4. Computational procedure 4.1. Domain, BC’s and mesh The dimensions of the simulation domain, boundary conditions and BFC mesh are indicated in Fig. 1. The case feature crosssectional dimensions similar to the PSBT OECD case; rod diameter D = 10 mm and pitch P = 13 mm. The length is reduced from 3.6 m to 1.0 m. A novel approach was used to generate proper boundary conditions for this flow. First, periodic conditions were applied in the radial and circumferential directions to mimic the effect of the neighboring rods. In the axial direction, however, a hybrid ‘Developed & Developing Flow Hybrid Approach’ has been developed, consisting first of generating turbulence in a periodic domain of length 2pDe, then the resulting fluctuating (scaled to maintain the mass flow rate) field is imposed in the entire domain, recycled periodically: temperature is updated using inflow-outflow BC’s. The BFC grid employed consists of 658,000 cells providing a near-wall resolution of an average y+  5, which avoids having to Table 1 Reference and reduced operating conditions for PSBT tests. Characteristics

Reference (Rubin et al.)

Reduced model

Pressure Saturation temperature Inlet temperature Mass flux Res = us De/2m Heat Flux Power

15.5 MPa 344.6 °C 290 °C 3333 kg/m2s 10,000 581 kW/m2 MW

15.5 MPa 344.6 °C 290 °C 74.1 kg/m2s 300 50 kW/m2 1.57 kW

The most important indicator of turbulence being developed is the wall shear velocity; us, which is shown in Fig. 2 to converge towards 0.0061 m/s for both SGS model simulations, against 0.0065 m/s for the earlier fine-mesh simulation of the group. The 5% difference resulting from the use of two different SGS models is clearly due to meshing of the near-wall region, and more precisely in the narrow gap, which in the present case is less-well resolved (y+  1–3) than in the DNS (y+  0.4–1.0). Fig. 3 below shows the contours of turbulent Eddy viscosity normalized by molecular viscosity at a later time step, obtained using both SGS models. In both simulations, the ratio is within the limit of LES requirements that is of the order of 5 and less than 10, meaning that most of the important flow scales are resolved already with the mesh. The WALE model provides a somewhat more patchy structure of this quantity, while the DSM returns a smoother and larger turbulent diffusion of momentum in the core region. Fig. 4 compares the time-averaged streamwise velocity and normal stress profiles obtained by the DSM and WALE SGS models. The results plotted at both 0 and 45 Deg. segments are also compared to the former DNS results of Caviezel et al. (2013). The velocity profiles compare rather well with the DNS regardless which SGS model is employed. The normal-stress profiles show differences between LES and DNS, reaching about 20% as to peak in particular. Overall, while the core flow region (45°) seems to be well predicted, the narrow-gap zone where the flow accelerates the most seems to be rather less well resolved. Finally, it looks like the DSM behaves slightly better than the WALE; both can be used for the rest of the investigation.

Please cite this article in press as: Lakehal, D., Caviezel, D. Large-Eddy Simulation of convective wall-boiling flow along an idealized PWR rod bundle. Nucl. Eng. Des. (2016), http://dx.doi.org/10.1016/j.nucengdes.2016.12.024

D. Lakehal, D. Caviezel / Nuclear Engineering and Design xxx (2016) xxx–xxx

5

P=13.0 mm L=1000 mm

Cross flow section

R =5 mm

θ=450 θ=00 Heat flux q” Fig. 1. Computational domain: Dimensions & BC’s. BFC grid used shown in the (x-y) cross section.

5.2. Flow structures results

Fig. 2. Space-averaged frictional velocity Comparison of the SGS models with DNS.

evolution

(single-phase

flow).

The flow structure at one instance near the central area obtained with the WALE model is shown in Fig. 5, depicting the cross-sectional velocity and temperature contours on the rod surface, with the black line indicating the location of saturation temperature (618°K). The WALE and DSM model provide virtually the same picture with marginal differences to be commented. The cross-sectional view shows flow structures of different size developing on the thermal boundary layer. Temperature contours on the surface of the rod show intermittent patches related to turbulence eddies transporting heat from the wall to the core flow. Important point to note is that the flow exhibits a strong secondary flow motion (intensity not to scale), the effect of which is on the wall temperature. These secondary cells were already predicted in the former refined LES and DNS studies of the group. This is clearly indicated by the time averaged temperature contours plotted in Fig. 6. Statistical convergence of the simulation is particularly well translated by the temperature contours plotted in the lower left panel.

Fig. 3. Instantaneous cross-sectional contours of SGS to molecular diffusivities (single-phase flow).

Please cite this article in press as: Lakehal, D., Caviezel, D. Large-Eddy Simulation of convective wall-boiling flow along an idealized PWR rod bundle. Nucl. Eng. Des. (2016), http://dx.doi.org/10.1016/j.nucengdes.2016.12.024

6

D. Lakehal, D. Caviezel / Nuclear Engineering and Design xxx (2016) xxx–xxx

Fig. 4. Mean W-velocity and w0 w0 -stress profiles at 0 and 45° segments. Results of simulations with WALE and DSM SGS models compared to DNS data.

Fig. 5. Instantaneous cross-sectional velocities and temperature contours along the rod (size not to scale).

As was to be expected from the flow statistics that were discussed in detail by Caviezel and Lakehal (2013), energy carried out by the largest scales (f < 1 Hz) is essentially concentrated in the axial velocity component W, although the other two components carry similar amounts in the high-frequency range (f > 1 Hz), with a decay behavior smoother in the flow direction than the other directions. The only difference between the two models is that the DSM preserves more energy than the WALE, which indicates more dissipation in the flow direction in particular (lower panels). The two panel of Fig. 8 is interesting in that it depicts other phenomena occurring at the highest-frequency range: The decay in the smallest scales range seems indeed to be rather peculiar, exhibiting a sharp drop in the range 10 < f < 25–30 Hz, before adjusting to the 5/3 slope in the interval 30 < f < 700 Hz. This peculiar behavior was actually observed in the experiments conducted by Rehme (1992), who pointed out to the effect of largescale structures. The absence of a clear peak in the power spectrum translating the coherent pulsations as was found by Rehme (1992) in the low frequency range (<10 Hz) is due to the limitation of size of the domain and periodic conditions applied in this LES context.

5.3. Energy spectra results 6. Convective boiling flow results DSM and WALE SGS model results are compared in this section (Figs. 7 and 8). Overall the same trend is observed, albeit some subtle differences do pertain as to the decay rate and dissipation rate.

We now turn to the results of the convective boiling flow using the reduced-model setup and flow conditions as discussed above,

Please cite this article in press as: Lakehal, D., Caviezel, D. Large-Eddy Simulation of convective wall-boiling flow along an idealized PWR rod bundle. Nucl. Eng. Des. (2016), http://dx.doi.org/10.1016/j.nucengdes.2016.12.024

D. Lakehal, D. Caviezel / Nuclear Engineering and Design xxx (2016) xxx–xxx

7

Fig. 6. Mean and instantaneous cross-sectional contours of W-velocity and temperature.

using the two SGS models separately. Again, we provide the results of figures-of-merit, macroscopic flow structures, energy spectra and 1st order time-averaged quantities to compare with the previous reference case. 6.1. Figure-of-merits results As stated previously, the frictional velocity us is a very important quantity in this context, the evolution of which is shown in Fig. 9 against flow through time. While both SGS models predicted it around 0.0061 m/s, the boiling-flow case shows a way higher value of around 0.009 m/s, almost 50% higher than the singlephase flow finding. This result cannot be a turbulence modelling or meshing issue since both SGS models return the same order of magnitude, and the limited grid effect has already been shown in the non-boiling flow. The result can clearly be interpreted as an increase in turbulence activity in the flow due to vapor generation, as it will be discussed in more detail next. Fig. 10 shows the contours of SGS Eddy viscosity normalized by molecular viscosity at a later time step, obtained using both SGS models. Overall the ratio is within the limit of LES requirements that is below ten, justifying the use of this mesh, however, there is a marked increase of this quantity beyond the boiling region,

and in particular at the later section, at z/L = 0.9. This suggests that a more refined grid would be helpful in the boiling section of the rod. The rate of Eddy viscosity generated in the transitional boiling region is comparable to the single-phase flow results discussed previously (with a maximum ratio of 4), and increases with the distance thus with vapor generation up to two folds far downstream (a maximum ratio of 8). Again, while WALE model provides a somewhat more patchy structure of this quantity, DSM returns larger values with a smoothed distribution reflecting momentum diffusion of additional boiling-induced scales. 6.2. Vapor generation results Fig. 11 shows the resulting mean and instantaneous steam concentration iso-contour (marked by 1%) on the rod resulting from water boiling at the wall, delimited at XONB. Note that these and next figures are not to scale; the dimensions were deliberately reduced to present the full picture of the flow. The threshold of a = 1% was selected by purpose such as to infer a better insight of what might happen near the wall. The main difference compared to the single-flow case discussed in the context of Fig. 5 is that the onset of nucleation occurs earlier, as was to be expected, since the mechanism is now entirely controlled by the wall boiling

Please cite this article in press as: Lakehal, D., Caviezel, D. Large-Eddy Simulation of convective wall-boiling flow along an idealized PWR rod bundle. Nucl. Eng. Des. (2016), http://dx.doi.org/10.1016/j.nucengdes.2016.12.024

8

D. Lakehal, D. Caviezel / Nuclear Engineering and Design xxx (2016) xxx–xxx

Fig. 7. Low-frequency PSD of U, V and W at various locations. (left) DSM results, (right) WALE results.

Fig. 8. Low- & high-frequency PSD of W velocity at various locations. Same nomenclature as above.

model. The results also indicate that the boiling activity increases with distance from the onset of bubble nucleation location. The boiling flow structure obtained with the WALE model only is illustrated in Fig. 12, depicting instantaneous cross-sectional velocity and temperature contours along the rod; the surface of the rod is colored by the vapor void fraction. While there are no major differences to be noted compared to single-phase flow simulations as to the global picture, the boiling case reveals stronger secondary flow motion than in the reference case. The rod surface colored by the vapor volume fraction shows a patchy structure due to intermittent wall streaks too. We could notice that the flow

seems to be portioned in three zones, (i) pre-boiling section at z/L < 0.35–0.4, (ii) transitional boiling in the range 0.4 < z/l < 0.55–0.6, and (iii) post-boiling section at z/L > 0.65. The boiling transition zone seems to be populated by isolated nucleation sites (note that the Logarithmic scale exaggerates their area) which should affect the flow in one way or the other; the post-boiling section of the rod is colored in green. Fig. 13 compares the time-average void-fraction profiles at three section; both SGS model results are reported. In both segments, although nucleation takes place already at z/L = 0.4, most R
dr) is actually generated in the post boiling of it (in the sense of a

Please cite this article in press as: Lakehal, D., Caviezel, D. Large-Eddy Simulation of convective wall-boiling flow along an idealized PWR rod bundle. Nucl. Eng. Des. (2016), http://dx.doi.org/10.1016/j.nucengdes.2016.12.024

D. Lakehal, D. Caviezel / Nuclear Engineering and Design xxx (2016) xxx–xxx

Fig. 9. Space-averaged frictional velocity evolution. Single-phase vs. & boiling-flow results using two SGS models.

section. The transitional boiling section is an interesting part, because the flow there seems not that much affected by the vapor spots, as we will be discussed later. 6.3. Flow structures results Fig. 14 depicts the time averaged cross-sectional contours of streamwise velocity (left panels) and its fluctuating component (right panels) at two stations, transitional boiling and postboiling. The results are obtained with the WALE model only;

9

results of the DSM model are quite similar. At the frontal section where little boiling has occurred only, the secondary flow is almost identical to the single-flow result discussed earlier, albeit with a certain skewness in the boiling case. The downward cross sections reveal a totally different picture of the boiling case: the boundary layer seems to be affected by the generated vapor, becoming clearly thinner than in the transitional boiling section. The flow in the core region is decelerated too. Looking now at the lower right panels reveals that the fluctuating turbulence intensity in the flow direction is also increased by the presence of the second phase, almost doubling in magnitude. The average-temperature contours plotted in Fig. 15 in the transitional boiling and postboiling sections corroborate with the previous finding, whereby the scalar is more diffuse in the core flow region of the postboiling section than in the transitional boiling flow portion of it. The real effect on the boundary layer can be inspected when looking the mean profiles plotted in Figs. 16 and 17. Fig. 16 compares the mean velocity in the axial direction: the plots include single and boiling flows at z/L = 0.25, 0.55 and 0.9, at both 0 and 45° segments. The flow is virtually the same in the pre-boiling section (z/L = 0.25) as the single phase flow, and marginally affected in the transitional boiling section (z/L = 0.55), it accelerates in the narrow-gap zone. The boiling section (z/L = 0.9) exhibits another picture, that is: the flow accelerates substantially near the wall because of the higher volumetric flow rate due to phase change, and the boundary layer thins as if the flow has increased in Reynolds number. The profiles reflect the increase of the frictional velocity discussed previously. In the core flow zone, the stream has considerably decelerated away from the wall. The profiles of the mean temperature shown in Fig. 17 are also affected by boiling. In the transitional boiling section (z/L = 0.55), the temperature remains almost identical to the single-phase flow case, where it develops in the axial direction quite substantially.

Fig. 10. Instantaneous cross-sectional contours of SGS to molecular diffusivities.

Please cite this article in press as: Lakehal, D., Caviezel, D. Large-Eddy Simulation of convective wall-boiling flow along an idealized PWR rod bundle. Nucl. Eng. Des. (2016), http://dx.doi.org/10.1016/j.nucengdes.2016.12.024

10

D. Lakehal, D. Caviezel / Nuclear Engineering and Design xxx (2016) xxx–xxx

Fig. 11. Mean and instantaneous vapor void-fraction contours on the rod. WALE model results (size not to scale).

The normal-stress profiles shown in Fig. 18 translate best the effect of boiling on turbulence activity, in the narrow gap zone in particular. The figure shows that in the transitional boiling section, turbulence has in fact diminished in intensity compared to the single-phase flow, which is perfectly replicated in the pre-boiling section. In the post-boiling, however, turbulence has substantially increased compared to the preceding sections. The mechanisms behind turbulence reduction (towards relaminarization) in the transitional-boiling region are unclear and require further refined studies. Be it as it may, one can say that while nucleation is triggered, the source(s) of turbulence has(have) somewhat vanished. Note finally that the DNS profiles were compared already to existing DNS of a pipe flow, showing a perfect agreement. A better insight into the effects of boiling on turbulence could be gained by inspecting the vorticity distribution in the singlephase flow and boiling flow cases (Fig. 19). Both SGS model results are plotted and compared to single-phase flow results (left panels). Interesting to note first is that both SGS models provide virtually the same picture: the pre-boiling zone, the transition zone, and the post-boiling section are all clearly visible with various degrees of vorticity population. While the left panels (single-phase flow) show a monotonic vorticity generation near the wall, the right panels suggest that boiling causes a decrease in vorticity in the transitional boiling zone followed by a higher rate of vorticity generation in the post-boiling section, which further accentuates beyond the ONB location. In the pre-boiling section, the rate of vorticity is similar to the single-phase flow, but reduces in the transitional boiling zone before picking up abruptly in the post-boiling section. This important result corroborates well with the previous observations concerning the modulation of turbulence activity by boiling. 6.4. Energy spectra results

Fig. 12. Instantaneous cross-sectional velocities and temperature contours (linear scale) along the rod, colored by the steam void fraction (Log. scale) on the surface. WALE model results (size not to scale).

Fig. 20 presents the energy spectra for W-velocity component (the most important one) obtained in the single-phase and boiling flow cases using both DSM and WALE SGS models. The figure clearly indicates that independent of the SGS model employed, the boiling flow always embodies about twice more energy than the single phase flow, with marked large-scale energetic eddies (f < 1 Hz) at the very high end of the rods, i.e. z/L = 0.8. The figure also shows an interesting flow transition result, in that at z/L = 0.5 where XONB hasn’t been reached yet, the energy spectrum falls within the single-phase flow range (purple line). This flow transition behavior is well illustrated in Fig. 21 depicting the boiling flow PSD’s of U and

Fig. 13. Time-average void-fraction profiles at 0 and 45° segments in the pre-boiling, transition and post-boiling flow sections. Simulations using WALE and DSM model.

Please cite this article in press as: Lakehal, D., Caviezel, D. Large-Eddy Simulation of convective wall-boiling flow along an idealized PWR rod bundle. Nucl. Eng. Des. (2016), http://dx.doi.org/10.1016/j.nucengdes.2016.12.024

D. Lakehal, D. Caviezel / Nuclear Engineering and Design xxx (2016) xxx–xxx

11

Fig. 14. Time averaged cross-sectional contours of stream wise mean velocity (left panels) and fluctuating component (right panels) along the rod: transitional boiling section; (top) post-boiling flow section (below).

Fig. 15. Time averaged cross-sectional contours of temperature (left panels) and fluctuating component (right panels) along the rod: transitional boiling; (top) post-boiling flow section (below).

Please cite this article in press as: Lakehal, D., Caviezel, D. Large-Eddy Simulation of convective wall-boiling flow along an idealized PWR rod bundle. Nucl. Eng. Des. (2016), http://dx.doi.org/10.1016/j.nucengdes.2016.12.024

12

D. Lakehal, D. Caviezel / Nuclear Engineering and Design xxx (2016) xxx–xxx

Fig. 16. Mean W-velocity profiles at 0 and 45° segments. Single-phase vs. & boiling-flow results at pre-boiling flow section, transitional boiling, and post-boiling flow section. Simulations using WALE model.

Fig. 17. Mean temperature profiles at 0 and 45° segments. Single-phase vs. & boiling-flow results at transitional boiling section and post-boiling flow section. Simulations using WALE model.

Fig. 18. w0 w0 -stress profiles at 0 and 45° segments. Single-phase vs. & boiling-flow results at pre-boiling section; transitional boiling (z/L = 0.55), and post-boiling flow section (z/L = 0.9). Simulations using WALE model.

Please cite this article in press as: Lakehal, D., Caviezel, D. Large-Eddy Simulation of convective wall-boiling flow along an idealized PWR rod bundle. Nucl. Eng. Des. (2016), http://dx.doi.org/10.1016/j.nucengdes.2016.12.024

D. Lakehal, D. Caviezel / Nuclear Engineering and Design xxx (2016) xxx–xxx

13

Fig. 19. Vorticity distribution along the rod bundle for single-phase and boiling flow cases. Upper: DSM model; lower: WALE model.

Fig. 20. Single-phase flow vs. boiling flow case PSD for the stream wise velocity component at various locations along the rod. (left) DSM results, (right) WALE results.

V velocity components taken at center of the channel at various cross sections, in that the amount of energy embodied in the flow systematically increases with distance along the rod, or in better words, with vapor production. These results are obtained with WALE model and are similar to the DSM model ones.

7. Conclusions The paper reports a detailed LES-based analysis of turbulent convective flow (with and without vapor-bubble nucleation) upward along the heated rods of a PWR sub-channel. Use is made of the CFD/CMFD code TransAT. The single-phase flow LES results have been compared to DNS data of the same flow, showing no major differences between the two SGS models employed, though quite different in nature and concept. The comparison has rather pointed out to importance of meshing of the gap region; a more suitable grid would consist in 1.2 million cells, the results of which were found close to the DNS data of the non-boiling case by

Caviezel et al. (2013). The simulation campaign has been extended to boiling heat transfer, using the triple flux decomposition model, coupled within the LES approach, where the main difference compared to the single-phase flow is that the onset of nucleation is now entirely controlled by the wall boiling model. The results in the wall boiling configuration differ from the single-flow case in many ways, but more particularly: the flow seems to be portioned in three zones, (i) pre-boiling where the flow is virtually the same as the single-phase flow, (ii) transitional boiling where turbulence production mechanisms diminish, and (iii) post-boiling, where turbulence activity increases with distance from the onset of nucleate boiling location along the rod, or with vapor production, marked by higher magnitude of the axial normal-stress component. The boundary layers thin quite substantially in the post-boiling section, and the flow accelerates due to higher volumetric flow rate induced by phase change. The wall frictional velocity increases considerably due to boiling. Finally, boiling causes a higher rate of vorticity generation, which further accentuates with the distance from XONB.

Please cite this article in press as: Lakehal, D., Caviezel, D. Large-Eddy Simulation of convective wall-boiling flow along an idealized PWR rod bundle. Nucl. Eng. Des. (2016), http://dx.doi.org/10.1016/j.nucengdes.2016.12.024

14

D. Lakehal, D. Caviezel / Nuclear Engineering and Design xxx (2016) xxx–xxx

Fig. 21. Convective boiling flow energy spectra for all velocity components at various locations. WALE model simulations.

Acknowledgements This work is supported by the NURESAFE (Nuclear Reactor Safety Simulation Platform) funded from EU resources under FP7-EURATOM-FISSION program with contract number 323263. References Bergles, A.E., Rohsenow, W.M., 1964. The determination of forced-convection surface-boiling heat transfer. J. Heat Transfer 86 (3), 365–372. Caviezel, D., Lakehal, D., Large Eddy & Interface Simulation (LEIS) of turbulent flow and convective boiling in a PWR rod bundle. In: Proc. 8th International Conference on Multiphase Flow, ICMF 2013, Jeju, Korea, May 26–31, 2013. Caviezel, D., Narayanan C., Lakehal, D., Highly-resolved LES (or DNS) of turbulent convective flow along a PWR rod bundle. In: Procs. 15th International Topical Meeting on Nuclear Reactor Thermal-hydraulics (NURETH-15), Pisa, Italy, May 12–17, 2013. Chatzikyriakou, D., Buongiorno, J., Caviezel, D., Lakehal, D., 2015. DNS and LES or turbulent flow in a closed channel featuring a pattern of hemispherical roughness elements. Int. J. Heat Fluid Flow 53, 29–43. Eddington, R.I., Kenning, D.B.R., 1979. The effect of contact angle on bubble nucleation. Int. J. Heat Mass Transfer 22, 1231–1236. Eggels, J.G.M., Unger, F., Weiss, M.H., Wetserweel, J., Adrian, R.J., Friedrich, R., Nieuwstadt, F.T.M., 1994. Fully developed turbulent pipe flow: a comparison between direct numerical simulation and experiment. J. Fluid Mech. 268, 175– 209. Incropera, H., DeWitt, H.J., 1996. Intr. Heat Transfer. Wiley & Sons. Kurul, N., Podowski, M.Z., Multidimensional effects in forced convection subcooled boiling. In: Proc. 9th Int. Heat Transfer Conference, Jerusalem, Israel, 21–26 August 1990.

Lakehal, D., Smith, B.L., Milelli, M., 2002. Large-eddy simulation of bubbly turbulent shear flows. J. Turbul. Manninen, M., Taivassalo, V., Kallio, S., 1996. On the Mixture Model for Multiphase Flow. Technical Research Centre of Finland, VTT Publications, Espoo. 288–354. Moin, P., Squires, K., Cabot, W., Lee, S., 1991. A dynamic subgrid-scale model for compressible turbulence and scalar transport. Phys. Fluids 3 (11), 2746–2757. Nicoud, F., Ducros, F., 1999. Subgrid-scale stress modelling based on the square of the velocity gradient tensor. Flow Turbul. Combust. 62, 183–200. Rehme, K., 1992. The structure of turbulence in rod bundles and the implications on natural mixing between the subchannels. Int. J. Heat Mass Transfer 35 (2), 567– 581. Rowe, D.S., Johnson, B.M., Knudsen, J.G., 1974. Implications concerning rod bundle crossflow mixing based on measurements of turbulent flow structure. Int. J. Heat Mass Transfer 17, 407–419. A. Rubin, M. Acramoca, H. Utsuno, OECD/NRC Benchmark based on NUPEC PWR subchannel and Bundle Tests (PSBT), NEA/NSC/DOC (2010). Seale, W.J., 1982. Measurements and predictions of fully developed turbulent flow in a simulated rod bundle. J. Fluid Mech. 123, 399–423. Thomas, S., Narayanan, C., Lakehal, D., Progress in the simulations of convective subcooled boiling in PWR sub-channels using TransAT. In: Procs. 15th International Topical Meeting on Nuclear Reactor Thermalhydraulics (NURETH-15), Pisa, Italy, May 12–17, 2013. Thomas, S., Narayanan, C., Lakehal, D., Validation of 3D simulations of the NUPEC PSBT in sub-channels using TransAT. In: Procs. 23rd International Conference on Nuclear Engineering, ICONE23, May 17–21, Makuhari Messe, Chiba, Japan, 2015. TransAT. ASCOMP software package. http://www.ascomp.ch/transat, (2015). Ünal, H.C., 1976. Maximum bubble diameter, maximum bubble growth time and bubble growth rate during subcooled nucleate flow boiling of water up to 17.7 MW/m2. IJHMT 19, 643–649. Vonka, V., 1988. Measurement of secondary flow cortices in a rod bundle. Nucl. Eng. Des. 106, 191–207.

Please cite this article in press as: Lakehal, D., Caviezel, D. Large-Eddy Simulation of convective wall-boiling flow along an idealized PWR rod bundle. Nucl. Eng. Des. (2016), http://dx.doi.org/10.1016/j.nucengdes.2016.12.024