- Email: [email protected]
Numerical investigation of turbulent flow in an annular sector channel with staggered semi-circular ribs using large eddy simulation
International Journal of Heat and Mass Transfer 123 (2018) 705–717
Contents lists available at ScienceDirect
International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt
Numerical investigation of turbulent flow in an annular sector channel with staggered semi-circular ribs using large eddy simulation Samuel J. Lee ⇑, Saya Lee, Yassin A. Hassan Department of Nuclear Engineering, Texas A&M University, 3133 TAMU, College Station, TX 77845, United States
a r t i c l e
i n f o
Article history: Received 13 December 2017 Received in revised form 2 March 2018 Accepted 6 March 2018
a b s t r a c t Computational fluid dynamics (CFD) simulations are conducted to explore flow inside a channel with staggered semi-circular ribs attached to both the inner and outer annular walls. Large eddy simulation (LES) with the wall-adapting local eddy-viscosity (WALE) sub-grid scale (SGS) model is used to simulate turbulent flow in the computational domain. The flow channel is modeled on the test section of a simplified helical coil steam generator (HCSG) design. Owing to their compactness and higher heat transfer coefficient in comparison to straight tube steam generators, HCSGs are employed in recent designs of advanced light water reactors in the nuclear industry. A Reynolds number of 13,900 is considered based on the mean inlet velocity, channel height, and water properties. Results from LES are compared with both an unsteady Reynolds-Averaged Navier-Stokes (URANS) simulation and experimental data acquired from the particle image velocimetry (PIV) method. Flow features observed inside the channel are presented by visualizing contours of velocities, vorticity, and Reynolds stresses on a monitoring plane. Time-averaged profiles of velocities and Reynolds stresses at five monitoring lines are also presented to validate the flow field resolved by LES with PIV. Additionally, vortical structures of the flow are captured using the Q-criterion. This study aims at providing flow characteristics of a specifically designed flow channel and acquiring the scalability of the numerical experiment. Ó 2018 Elsevier Ltd. All rights reserved.
1. Introduction A helical coil heat exchanger (HCHX) is a type of shell-and-tube heat exchanger characterized by its compact structure and high heat transfer coefficient with relatively low pressure drop [1]. Thus, it has been employed in various types of engineering applications such as chemical process plants, cryogenic systems, heat recovery processes, and other energy intensive industries [2–4]. Notably, a recent design of advanced light water reactor (LWR) adopted a HCHX as a steam generator [5]. This type of nuclear reactor is emerging as a reflection of demands for increased scalability, lower initial investment, and greater operational safety than conventional large reactors [6]. They are small to medium sized integral type pressurized water reactors (PWRs) in which the helical coil steam generator (HCSG) is placed inside the same vessel as the reactor core and immersed in the coolant which is circulated by natural convection [7]. To ensure the safe and reliable operation of the reactor, investigation of the characteristics of the flow field around the helical coil bundle is essential. Understanding how the flow field around coil bundle affects the pressure distribution ⇑ Corresponding author. E-mail addresses: [email protected] (S.J. Lee), [email protected] (S. Lee), [email protected] (Y.A. Hassan). https://doi.org/10.1016/j.ijheatmasstransfer.2018.03.026 0017-9310/Ó 2018 Elsevier Ltd. All rights reserved.
on the helical coil surface is crucial to determine the structural and thermal safety of the reactor pressure vessel (RPV) [8]. Most of the numerous investigations on HCHX are concerned with the convective heat transfer and pressure drop in a coil tube for the tube-side and outside of a coil tube, which is simplified as a cylindrical tube array for the shell-side [9–16]. However, it is rare to find studies of detailed flow around a realistic HCHX model. An experimental test section simplifying the practical helical coil steam generator (HCSG) model [17] was constructed by Texas A&M University to examine the fluid phenomena around the shellside helical coil bundle [18]. Although the test facility was designed to simulate the natural circulation of primary coolant inside the reactor vessel, the driving force in the test facility is provided by a hydraulic pump. A honeycomb structure was installed prior to the inlet of the test section to make the flow as uniform as possible under the assumption of isothermal and uniform downward inlet velocity. To be specific, the HCSG has multiple concentric layers of helical coils in the radial direction, and each layer has a spiral direction different from the adjacent layers. To simplify the geometry, two adjacent layers were sliced in the radial direction. Further simplification is achieved by cutting the annular section into 1/25 sectors. The acrylic test section was fabricated based on this geometry, as illustrated in Fig. 1. Consequently, the flow channel
706
S.J. Lee et al. / International Journal of Heat and Mass Transfer 123 (2018) 705–717
Fig. 1. Simplified HCSG model for the experimental test section.
inside the test section has concentric annular-shaped walls and semi-circular ribs on each wall surface staggered with those on the other wall. The rib represents the coil surface of the HCSG. Previously, an unsteady Reynolds-Averaged Navier-Stokes (URANS) simulation was carried out using the flow channel extracted from the simplified test section as the computational domain for comparison with the experimental results [19]. In the URANS simulation, the time-averaged static pressure drop between the inlet and outlet of the test section showed good agreement with the experimental result. Typically, URANS simulation produces reasonable averaged statistics, but it is not appropriate to discern detailed instantaneous fields. To accurately resolve the instantaneous field and detailed flow structures, a small-scale resolving turbulence model is required. The most accurate method for simulating turbulent flow is direct numerical simulation (DNS). However, DNS is still restricted to low Reynolds numbers, which differ considerably from most practical engineering applications. Although the computing power has dramatically improved in the last few decades owing to high performance computer clusters, DNS remains as a computationally voracious method. For this circumstance, the LES approach is a realistic alternative with which investigates instantaneous turbulent flow fields. In this paper, the LES results of the simplified HCSG test section are found to provide reliable numerical data. The flow field resolved by LES was validated by the experimental PIV measurement. The flow inside the channel was studied by analyzing the velocity field. The vortex identification method was also employed to identify vortical structures that affect the pressure distribution of the channel. In addition, the static pressure drop of the test section was calculated and compared with the experimental data. 2. Numerical details The numerical experiment was performed on the flow channel inside the simplified HCSG test section using the finite volume method (FVM) based commercial CFD package, STAR-CCM+ ver.11.06-R8. To resolve the turbulent flow field of the computational domain, LES models the flow field at a scale smaller than the grid size. 2.1. Computational domain and boundary conditions The computational domain extracted from the test section is shown in Fig. 2. This geometry has semi-circular ribs attached to both concavely and convexly curved walls, each at an angle, and the ribs on one wall are staggered with those on the other wall. The diameter of the ribs is D = 15.9 mm. The height of the channel
at the inlet is h = 23.6 mm, and the gap between the ribs on both walls is d = 7.7 mm. A uniform velocity inlet and a constant pressure outlet were applied as boundary conditions of the computational domain. inlet = 0.525 m/s, corresponds to The value imposed at the inlet, u a volume flow rate of 3.15 103 m3/s (50 gpm). The resulting stream = mean streamwise velocity at the ribbed region was u 1.601 m/s. The non-uniform time varying inlet condition, discussed in Section 4.6, caused no significant inlet effect at the monitoring position of the test section; thus, a uniform velocity profile was applied to the inlet. For the outlet condition, a static pressure of 0 Pa was imposed. The density and dynamic viscosity of the working fluid are those of water at room temperature, which are q = 997.561 kg/ m3 and l = 8.887 104 kg/m s, respectively. The Reynolds number based on the mean inlet velocity, channel height at the inlet, inlet/l 13,900. and properties of water is Reh = qhu 2.2. Computational setup The unsteady isothermal incompressible fluid flow field inside the simplified HCSG test section was calculated adopting LES as a turbulence model. The SIMPLE algorithm was used to solve the discretized Navier-Stokes equations. The second-order implicit scheme was applied to discretize the governing equations temporally. The integral time scale and the time step were T = 1 s and Dt = 104 s, respectively. The bounded-central difference scheme [20] was employed to discretize the convection term with an upwind blending factor of r = 0.0 to reduce numerical diffusion. The hybrid Gauss-least squares (LSQ) method was used to reconstruct the gradient terms. The number of iterations per time step was 20, and the convergence of the solution at each time step had residuals below 104. The computational domain was equally divided into three sub-parts, cutting the domain using the x-y plane for vortical structure monitoring. Neighboring boundaries of sub-parts are nonconformal mesh connected with the interface. To ensure that a fully developed field inside the computational domain was used in the analyses, the initial 0.2 s flow fields were excluded. 2.3. Sub-grid scale model for LES The wall-adapting local eddy-viscosity (WALE) sub-grid scale (SGS) model was used to model the scale smaller than the grid size j , where the overbar i u [21]. The SGS stress tensor, sij ¼ ui uj u means filtered quantity, can be modeled based on an eddyviscosity assumption as follows:
1 3
sij skk dij ¼ 2mt Sij
ð1Þ
S.J. Lee et al. / International Journal of Heat and Mass Transfer 123 (2018) 705–717
707
Fig. 2. Computational domain of the HCSG test section.
i =@xj þ @ u j =@xi Þ is the filtered strain rate tensor where Sij ¼ ð1=2Þð@ u and mt is the SGS eddy-viscosity. To acquire the local eddy-viscosity, the square of the velocity gradient tensor which considers both the strain and the rotation rate was used to formulate the WALE model. In the WALE model, SGS eddy viscosity is defined as mt = (CwD)2Sw where Cw is a model coefficient, D the grid filter width, and Sw a deformation parameter defined as
3=2 Sdij Sdij Sw ¼ 5=4 Sij Sij 5=2 þ Sd Sd ij ij
ð2Þ
where Sdij is the traceless symmetric part of the square of the velocity gradient tensor.
Sdij ¼
i @ u k @ u j @ u k k @ u l 1 @u 1 @u dij þ 2 @xk @xj @xk @xi 3 @xl @xk
lent flow: energy-containing range, inertial subrange, and dissipation range. It is expected that LES can resolve the inertial subrange in the computational domain. That being so, the velocity component of the y-directional streamwise was monitored at the midpoint, represented by the yellow1 dot in Fig. 3, which is 70.9 mm in the y-direction from the top of the first semi-circular rib on the convex wall. The instantaneous value of the streamwise velocity component at this monitoring point was acquired at each time step in the 0.5 s period. The fast Fourier transform (FFT) was then taken from this velocity fluctuation and is illustrated in Fig. 4. It is observed that the power density spectrum covered the distinguishing 5/3 slope which represents the inertial subrange. Therefore, it can be assumed that the grid system has sufficient resolution to resolve the inertial subrange.
ð3Þ
The salient features of the WALE model are that it yields proper near-wall scaling for eddy-viscosity without requiring the Van Driest damping function, and it is less sensitive to the model coefficient Cw [21]. 2.4. Grid system A nested Cartesian grid system with ten layers of prism mesh was constructed on the extracted computational domain, as illustrated in Fig. 3. Since the nested grid system uses predominantly regular hexahedral meshes, high quality cells with minimum cell skewness fill the computational domain. The base size of the nested grid system was 1.5 mm. The grid on the ribbed region was refined with 25% of the base size as 0.375 mm to resolve detailed flow structures. The first off-wall grid distance, Dy = 103 mm, was selected to make Dy+ 1 after performing preliminary simulations. There were approximately 37 million cells inside the domain. Kolmogorov’s 5/3 law with Taylor’s hypothesis was invoked to check the validity of this grid system for the present numerical study [22]. There are three ranges in the energy spectrum of turbu-
2.5. URANS simulation In this study, the Reynolds stress transport (RST) model with the linear pressure strain two-layer model was selected as a turbulence model for the URANS simulation. The same grid system as that in LES was used in the URANS simulation for comparison. For the temporal and spatial discretization schemes, the secondorder implicit scheme with 104 s time-step, and the MUSCL/CD scheme were applied, respectively. The SIMPLE algorithm was used as a pressure correction. In addition, all boundary conditions applied to the computational domain were the same as in LES. 3. Experimental data The experimental data presented in this paper are the result of PIV analysis [23] on the acrylic test section shown in Fig. 2. The schematic of the simplified HCSG test facility for the PIV measurement is illustrated in Fig. 5. The PIV system consisted of a highspeed camera and 532 nm continuous laser with 10 W power. 1 For interpretation of color in Figs. 3, 6 and 9, the reader is referred to the web version of this article.
708
S.J. Lee et al. / International Journal of Heat and Mass Transfer 123 (2018) 705–717
concave wall side Ǹ
convex wall side Ǹ
Fig. 3. Section view of the grid system on the computational domain.
Fig. 4. FFT of streamwise velocity fluctuation of a point inside the channel.
Fig. 6. Monitoring plane and lines for data acquisition.
Fig. 5. Schematic of the simplified HCSG test section facility for the PIV measurement.
The frame rate of the PIV picture was 3500 fps. The size of the image was 332 1024 pixels with a magnification factor of 70.9 lm/pixel. Silver-coated trace particles with a 13 lm mean cutoff diameter and density of 1.6 g/cm3 were seeded to water. In the post processing of PIV analysis, the adaptive FFT-correlation [24] with three-step interrogation window sizes of 32 16, 16 8, and 8 4 pixels was utilized. The uncertainty of the PIV measurement in accordance with the ITTC guide [25] was uc = 6.6% of the mean streamwise velocity. In addition to the PIV measurement, the static pressure drop across the test section was obtained using pressure transducers installed at two red dot positions as shown in Fig. 6.
S.J. Lee et al. / International Journal of Heat and Mass Transfer 123 (2018) 705–717
PIV
URANS
709
LES
v [m/s] Fig. 7. Instantaneous streamwise velocity, v, contour.
PIV
URANS
LES
[m/s] Fig. 8. Mean streamwise velocity, hvi, contour.
4. Results and discussion In the CFD simulation, the same monitoring positions were set inside the computational domain to compare the simulation result with the experimental data, as shown in Fig. 6. The monitoring plane was made by slicing the x-y plane of the fluid domain 62.7 mm from the right-side wall. The velocity field in the boxed plane was resolved, and various properties were calculated to illustrate each contour. The velocity components are defined as ui = huii + ui0 , where the angled bracket and prime indicate the mean and fluctuating components of velocity, respectively. In this study, over
Fig. 9. Mean streamwise velocity, hvi, profiles on lines C1 to C5.
10,000 instantaneous fields were used to obtain a mean value. There are five monitoring lines around the sixth semi-circular rib on the convex wall from the upstream position used to plot various
710
S.J. Lee et al. / International Journal of Heat and Mass Transfer 123 (2018) 705–717
PIV
URANS
LES
u [m/s] Fig. 10. Instantaneous transverse velocity, u, contour.
PIV
URANS
LES
[m/s] Fig. 11. Mean transverse velocity, hui, contour.
mean properties. In addition, the center section of the computational domain was utilized to show the three-dimensional vortical structures resolved by LES. 4.1. Streamwise and transverse velocity fields Fig. 7 illustrates the instantaneous streamwise (negative ydirection) velocity, v; and the contours resolved by PIV, URANS, and LES. As seen in the figure, both PIV and LES exhibit the characteristic features of the instantaneous velocity field, whereas the contour of the URANS simulation shows a velocity field similar to its mean field. Although the URANS simulation solves the Navier-
Fig. 12. Mean transverse velocity, hui, profiles on lines C1 to C5.
Stokes equation using a turbulence closure model with an unsteady manner, the RANS turbulence model resolves the instantaneous velocity field time-averaged by its very nature.
S.J. Lee et al. / International Journal of Heat and Mass Transfer 123 (2018) 705–717
PIV
URANS
711
LES
[m/s] Fig. 13. Mean Reynolds normal stress, hv0 v0 i, contour.
The mean, or time-averaged, streamwise velocity, hvi, contours are presented in Fig. 8. Contours for all three cases show the characteristic feature of the mean velocity field. In other words, there are no abrupt changes in the streamwise velocity in the observed area. Regions of relatively high magnitude of downstream velocity appear near the channel core. Specifically, regions of the highest streamwise velocity magnitude are captured right next to the rib surfaces. The cavity between two consecutive ribs of the same wall-side exhibits counter streamwise velocity due to flow recirculation. The mean streamwise velocity profiles on lines C1 to C5 depicted in Fig. 6 are plotted in Fig. 9. As indicated in the legend box on the first plot, the black solid circle symbol, black dashed line, and blue solid line are the PIV, URANS, and LES data, respectively. The 0 mm position in each plot indicates the leftmost wall of each corresponding line. The results of the LES show good agreement with the PIV near the wall and cavities and are underestimated near the channel core region. For the case of URANS, the simulation results underestimate the mean streamwise velocity compared to PIV on the whole monitoring lines. Figs. 10 and 11 illustrate the instantaneous and mean transverse x-directional velocity, u and hui, contours, respectively. Likewise, in the case of streamwise velocity, the contours for the instantaneous transverse velocity fields of PIV and LES show the characteristic feature of instantaneous field; however, URANS shows that the instantaneous field looks similar to the mean field in Fig. 11. The notable feature observed in the mean transverse velocity contour is the stagnation point, which can be distinguished by finding the point where the sign change occurs on the semi-circular rib surface facing the upstream direction. The mean transverse velocity profiles on lines C1 to C5 are shown in Fig. 12. Profiles for the LES show better agreement with the PIV results than the URANS simulation. However, for two lines, C1 and C5, the LES overestimated the mean transverse velocity compared to the PIV in close vicinity of the right region near the concave ribbed wall. For the cases of C2 and C3, the profiles of LES skewed toward the right wall.
Fig. 14. Mean Reynolds normal stress, hv0 v0 i, profiles on lines C1 to C5.
The mean transverse velocity profiles of the URANS simulation on C1, C2, and C3 deviate considerably from the PIV results. The trend of the URANS profile on C4 follows the PIV results, but its
712
S.J. Lee et al. / International Journal of Heat and Mass Transfer 123 (2018) 705–717
PIV
URANS
LES
[m/s] Fig. 15. Mean Reynolds normal stress, hu0 u0 i, contour.
magnitude is substantially underestimated. For the C5 profile, the URANS simulation matched well with the PIV data in the 5 < x < 10 mm range but deviated over other ranges. 4.2. Mean Reynolds stresses The contours of the mean Reynolds normal stresses, hv0 v0 i and hu0 u0 i, are shown in Figs. 13 and 15, respectively. In addition, the corresponding profiles for the mean Reynolds normal stresses captured on lines C1 to C5 are plotted in Figs. 14 and 16, respectively. As seen in the contours, the LES results have better agreement with those of the PIV than the URANS simulation. The regions of local maximum values of Reynolds normal stresses of the PIV and LES contours are similar. However, in the case of the URANS simulation, regions of local maximum values deviate from the others. Observation of the profiles of these mean Reynolds normal stresses on lines C1 to C5 shows obvious differences between LES and URANS. The LES profiles agree well with the PIV data, whereas the profiles of the URANS simulation are underestimated compared to PIV. Fig. 17 shows the contours of the mean Reynolds shear stress, hu0 v0 i. In addition, the corresponding profiles on lines C1 to C5 are plotted in Fig. 18. Regions of local maximum values are observed behind the ribs on both the PIV and LES contours, whereas the local maximum regions on the URANS contour show different aspects to the other cases. The profiles of the mean Reynolds shear stress for LES follow the PIV trends well, and the magnitudes of the values are underestimated at several locations. For the URANS simulation, the profiles partially follow the trends of PIV; however, it is observed that the profiles deviate from the PIV data. 4.3. Vorticity The instantaneous and mean vorticity, xz and hxzi, contours are illustrated in Figs. 19 and 20, respectively. As seen in the previous contours of velocities, the vorticity contours of URANS show no difference between the instantaneous and mean fields, and the
Fig. 16. Mean Reynolds shear stress, hu0 u0 i, profiles on lines C1 to C5.
instantaneous vorticity contours of PIV and LES show the distinctive features of the instantaneous field. In the case of mean vorticity contours, the three contours show similar patterns of mean
S.J. Lee et al. / International Journal of Heat and Mass Transfer 123 (2018) 705–717
PIV
URANS
713
LES
[m/s] Fig. 17. Mean Reynolds shear stress, hu0 v0 i, contour.
vorticity. Relatively high vorticity magnitude regions are observed behind the top surface of the ribs in the downstream direction. Fig. 21 shows the profiles of mean vorticity on lines C1 to C5. Both URANS and LES have good agreement with PIV, other than in proximity to the wall. For all five lines, LES has a high magnitude of mean vorticity near the wall, because LES can resolve the velocity gradient better in close vicinity to the wall than PIV. 4.4. Vortical structures High magnitude of vorticity is not a sufficient condition for the existence of a strong vortex. Therefore, Q-criterion, which is a vortex identification method suggested by Hunt et al. [26], was utilized to identify vortical structures inside the computational domain. Indeed, Q is the second invariant of the velocity gradient tensor, @ui/@xj. Assuming an incompressible flow, Q is defined as
Q
1 2
@ui @ui @uj 1 @ui @uj 1
Xij 2 Sij 2 ¼ ¼ 2 @xj @xi 2 @xi @xj @xi
ð4Þ
where Xij and Sij are the rotation and strain rate tensors, respectively, and |||| is the Frobenius norm. The Q-criterion utilizes the positive value of Q to represent the vortex in a fluid region. In other words, it identifies the vortical structure, finding a fluid region in which the rotation rate is larger than the strain rate of the velocity field. Fig. 22 shows vortical structures resolved by LES inside the computational domain using the Q-criterion. Velocity magnitude contours on the iso-surface of each Q value are presented. Observation of the figure indicates that vortical structures move down in the streamwise direction after they are generated next to the rib surface near the channel core region. This phenomenon is clearly observed on the contour of Q on the center plane, as shown in Fig. 23. Some vortices are captured inside the cavity between two consecutive ribs, and the others travel downstream.
Fig. 18. Mean Reynolds shear stress, hu0 v0 i, profiles on lines C1 to C5.
4.5. Static pressure drop To measure the static pressure drop across the test section, pressure transducers were installed at both the inlet and outlet
pressure taps of the test section shown in Fig. 6. The pressure drop, Dp = pin pout, was 2855 ± 8.4 Pa which was obtained from the ensemble average of eighteen data sets. Each data set has over
714
S.J. Lee et al. / International Journal of Heat and Mass Transfer 123 (2018) 705–717
URANS
PIV
z
LES
[s-1] Fig. 19. Instantaneous vorticity, xz, contour.
PIV
URANS
LES
< z> [s-1] Fig. 20. Mean vorticity, hxzi, contour.
10,000 samples. In the CFD simulations, fluctuating pressure drops were averaged from each T = 1 s simulation. URANS and LES calculated the pressure drop as 2343 Pa and 2648 Pa, respectively. Both CFD simulations yielded lower pressure drops than the experiment, which are 17.9% and 7.3%, for URANS and LES, respectively. LES predicted the pressure drop measured from the experiment better than the URANS simulation. Fig. 21. Mean vorticity, hxzi, profiles on lines C1 to C5.
4.6. Discussion In this study, LES resolved the flow field better than the URANS simulation, as seen in both the instantaneous and mean flow fields.
Since the geometry of the flow channel is not symmetric and the Reynolds number, Reh = 13,900, is in the turbulent regime, the instantaneous flow field inside the domain showed strong
S.J. Lee et al. / International Journal of Heat and Mass Transfer 123 (2018) 705–717
715
metric difference between the CAD model used to extract the computational domain and the actual flow channel of the test section can be considered a source of uncertainty. The CFD simulations in this study used the uniform inlet boundary condition. However, the experiment showed that the upper
Q = 4.0 105 s -2
Q = 2.0 105 s -2
Q = 1.5 105 s -2
Q = 1.0 105 s -2
Q [s −2 ] ur | V | [m / s]
Fig. 23. Contour of Q value resolved by LES on the center plane of the middle section.
Fig. 22. Vortical structure resolved by LES using the iso-surface of Q-criterion with various values.
unsteadiness. Therefore, to examine the flow phenomena for this flow properly, selection of an appropriate turbulence model that can resolve small scale unsteady features in the flow is essential. Relatively large discrepancies between the PIV and the LES profiles are observed on line C3. In raw pictures of the PIV measurement, a shadow line induced by optical distortion due to the complex geometry of the acrylic test section on the laser sheet appears right below line C3. Consequently, it is assumed that the shadow has an effect on the PIV analysis near line C3, and the resolved field deviates from the field resolved by unaffected PIV pictures. This shadow line effect is clearly observed in Figs. 8 and 13. In addition, PIV is not able to capture the high velocity gradient very close to the wall, as seen in Fig. 21, since this method has limited resolution capabilities in close proximity to the wall. The geo-
LES
LES w/ SEM
|V| [m/s] Fig. 24. Velocity magnitude contours near the inlet boundary for the uniform velocity and synthetic-eddy-method cases.
716
S.J. Lee et al. / International Journal of Heat and Mass Transfer 123 (2018) 705–717
Fig. 25. Mean streamwise and transverse velocities, hv i and hui; mean Reynolds normal and shear stresses, hv0 v0 i, hu0 u0 i, and hu0 v0 i; and mean vorticity hxzi profiles on line C5 for LES with uniform inlet (blue solid line) and for LES with SEM (red dashed line). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
part of the test section generated a non-uniform inlet velocity profile with a turbulent intensity (TI) of 0.107. The monitoring positions are located after the third rib in the streamwise direction. It is necessary to ascertain whether the inlet condition affects the flow field at the monitoring positions. There is a numerical method to generate realistic inlet flow conditions for LES; it is called the synthetic-eddy-method (SEM) and was suggested by Jarrin et al. [27]. The turbulent intensity measured at the inlet position of the experimental test section was TI = 0.107, and the honeycomb mesh size can be considered as the eddy turbulent length scale, le = 3.97 mm. The SEM is applied to the inlet boundary condition of the computational domain using these two values. Fig. 24 illustrates the velocity magnitude contours on the center plane of the computational domain for both the uniform velocity and the SEM. Unlike in the uniform velocity inlet case, velocity fluctuations are observed in the SEM case owing to the existence of synthetic eddies. Comparisons of various mean properties between LES with a uniform inlet and LES under the SEM condition, are plotted on line C5 in Fig. 25. As seen in the figure, even if the inlet conditions for both cases are different, results for the mean profiles show no significant difference at the region of interest. Therefore, it can be assumed that the flow fields resolved by the experimental PIV method with the unsteady non-uniform inlet and LES with a uniform inlet condition are comparable and reliable. According to the review paper of Païdoussis [28] on vortex shedding cylinder arrays in cross flow, the flow pattern after the third row remains almost unchanged owing to the influence of the turbulence gener-
ated by the cylinder placed at an upstream position. Since this study deals with channel flow, which can be seen as a part of cross flow of the helical coil tube bundle, a similar phenomenon is expected.
5. Summary and conclusion A numerical investigation of the turbulent flow field inside a channel with a specific configuration of semi-circular ribs that simulates the simplified HCSG geometry was conducted using LES. The results of the LES were validated by experimental data. In particular, the WALE SGS model was used to resolve the flow in the computational domain with a Cartesian nested grid system. The grid resolution was sufficient to resolve the inertial subrange of turbulence. In addition, a URANS simulation with an RST model was performed on the same grid system to provide comparison data. Instantaneous and mean flow fields resolved by PIV, URANS, and LES were compared at the region of interest. Owing to the nature of each turbulence model, LES resolved the unsteady flow features better than the URANS simulation. Moreover, the profiles of the mean velocities and Reynolds stresses showed that LES correlated better with the PIV data. Vortical structures resolved by LES were visualized using the Qcriterion. Vortices were induced next to the rib surface near the channel core region, and then they detached from the rib wall and were carried away by the bulk flow downstream. Contours
S.J. Lee et al. / International Journal of Heat and Mass Transfer 123 (2018) 705–717
of the Q value on a monitoring plane enabled to locate the vortex cores. Additionally, some vortical structures, which were resolved in the cavity between two consecutive ribs, were observed. The static pressure drop calculated by both numerical and experimental studies was presented. LES predicted the pressure drop of the test section with higher accuracy than the URANS simulation when compared to the experimental data. Several sources of discrepancy between the experiment and the CFD simulation were discussed. On the experimental side, the uncertainty of the raw pictures of PIV and the geometric difference were expected to affect this comparison. On the CFD side, the nonrealistic uniform inlet boundary condition applied to the computational domain was examined as a possible cause of discrepancy. However, the profiles of the mean flow fields on the monitoring position showed no significant difference to the case of the more realistic inlet boundary condition generated by the SEM. This leads to the conclusion that the inlet condition does not significantly affect the results of the CFD simulation at the region of interest. This is also applicable to the experimental facility; thus, the inlet condition is not considered a valid cause of discrepancy between the experimental and CFD results. For future research, a more realistic test section geometry will be constructed, and the extracted computational domain will capture the actual HCSG flow phenomena more accurately. Conflict of interest There are no known conflicts of interest associated with this publication. Acknowledgements This study is supported by the Integrated University Program (IUP) graduate fellowship, the U.S. Department of Energy Office of Nuclear Energy (DOE-NE). The authors acknowledge Texas A&M University High Performance Research Computing for providing computing power to perform the CFD simulations of this study. References [1] R.K. Patil, B.W. Shende, P.K. Ghosh, Designing a helical-coil heat exchanger, Chem. Eng. 13 (12) (1982) 85–88. [2] M.E. Ali, Experimental investigation of natural convection from vertical helical coiled tubes, Int. J. Heat Mass Transf. 37 (4) (1994) 665–671. [3] R.C. Xin, A. Awwad, Z.F. Dong, M.A. Ebadian, H.M. Soliman, An investigation and comparative study of the pressure drop in air-water two-phase flow in vertical helicoidal pipes, Int. J. Heat Mass Transf. 39 (4) (1996) 735–743. [4] D.G. Prabhanjan, G.S.V. Raghavan, T.J. Rennie, Comparison of heat transfer rates between a straight tube heat exchanger and a helically coiled heat exchanger, Int. Commun. Heat Mass Transf. 29 (2) (2002) 185–191.
717
[5] International Atomic Energy Agency, Status of advanced light water reactor designs, IAEA-TECDOC-1391, Vienna, Austria, 2004. [6] International Atomic Energy Agency, Instrumentation and control systems for advanced small modular reactors, IAEA Nuclear Energy Series No. NP-T-3.19, Vienna, Austria, 2017. [7] International Atomic Energy Agency, Natural circulation in water cooled nuclear power plants: Phenomena, models, and methodology for system reliability assessments, IAEA-TECDOC-1474, Vienna, Austria, 2005. [8] International Atomic Energy Agency, Design safety considerations for water cooled small modular reactors incorporating lessons learned from the Fukushima Daiichi accident, IAEA-TECDOC-1785, Vienna, Austria, 2016. [9] A.N. Dravid, K.A. Smith, E.W. Merrill, P.L.T. Brian, Effect of secondary fluid motion on laminar flow heat transfer in helically coiled tubes, AICHE J. 17 (5) (1971) 1114–1122. [10] S.V. Patankar, V.S. Pratap, D.B. Spalding, Prediction of laminar flow and heat transfer in helically coiled pipes, J. Fluid Mech. 62 (3) (1974) 539–551. [11] S.A. Berger, L. Talbot, L.S. Yao, Flow in curved pipes, Annu. Rev. Fluid Mech. 15 (1983) 461–512. [12] K. Futagami, Y. Aoyama, Laminar heat transfer in a helically coiled tube, Int. J. Heat Mass Transf. 31 (2) (1988) 387–396. [13] P. Naphon, S. Wongwises, A review of flow and heat transfer characteristics in curved tubes, Renew. Sust. Energy Rev. 10 (5) (2006) 463–490. [14] S. Kakaç, H. Liu, A. Pramuanjaroenkij, Heat Exchangers: Selection, Rating, and Thermal Design, CRC Press, 2012. [15] P.V. Gilli, Heat transfer and pressure drop for cross flow through banks of multistart helical tubes with uniform inclinations and uniform longitudinal pitches, Nucl. Sci. Eng. 22 (3) (1965) 298–314. [16] S. Majumdar, A. Moisseytsev, K. Natesan, Assessment of the Next Generation Nuclear Plant Intermediate Heat Exchanger Design, Argonne National Laboratory EXT-08/32, 2008. [17] D.K. Williams, D.P Fassett, B.J. Webb, W.J. Bees, A.S. Kruskamp, Helical Coil Steam Generator, U.S. Patent 20120111287 A1, Issued Apr. 12, 2013. [18] M. Delgado, S.J. Lee, S. Lee, Y.A. Hassan, Pressure drop measurements across a helical coil steam generator configuration, in: Proceedings of International Topical Meeting on High Temperature Reactor Technology (HTR 2016), Las Vegas, Nevada, USA, November 6–10, 2016. [19] S.J. Lee, S. Lee, Y.A. Hassan, Numerical investigation of the flow inside the helical coil steam generator test section, in: Proceedings of International Topical Meeting on High Temperature Reactor Technology (HTR 2016), Las Vegas, Nevada, USA, November 6–10, 2016. [20] B.P. Leonard, The ULTIMATE conservative difference scheme applied to unsteady one-dimensional advection, Comput. Meth. Appl. Mech. Eng. 88 (1) (1991) 17–74. [21] F. Nicoud, F. Ducros, Subgrid-scale stress modelling based on the square of the velocity gradient tensor, Flow Turbul. Combust. 62 (3) (1999) 183–200. [22] S.B. Pope, Turbulent Flows, Cambridge University Press, 2000. [23] S. Lee, M. Delgado, C. Guilbault, Y.A. Hassan, Visualization of flow interaction in a channel flow produced by a pair of helically coiled rods, Transactions of ANS Annual Meeting, San Francisco, California, USA, June 11–15, 2017. [24] W. Thielicke, E.J. Stamhuis, PIVlab - towards user-friendly, affordable and accurate digital particle image velocimetry in MATLAB, J. Open Res Software 2 (1) (2014) e30, https://doi.org/10.5334/jors.bl. [25] ITTC, ITTC recommended procedures and guidelines - uncertainty analysis particle image velocimetry, Technical Report, 2008. [26] J.C.R. Hunt, A.A. Wray, P. Moin, Eddies, streams, and convergence zones in turbulent flows, in: Proceedings of the 1988 Center for Turbulence Research Summer Program (CTR-S88), Stanford, California, USA, 1988. [27] N. Jarrin, S. Benhamadouche, D. Laurence, R. Prosser, A synthetic-eddy-method for generating inflow conditions for large-eddy simulations, Int. J. Heat Fluid Flow 27 (4) (2006) 585–593. [28] Païdoussis, A review of flow-induced vibrations in reactors and reactor components, Nucl. Eng. Des. 74 (1) (1983) 31–60.
Contents lists available at ScienceDirect
International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt
Numerical investigation of turbulent flow in an annular sector channel with staggered semi-circular ribs using large eddy simulation Samuel J. Lee ⇑, Saya Lee, Yassin A. Hassan Department of Nuclear Engineering, Texas A&M University, 3133 TAMU, College Station, TX 77845, United States
a r t i c l e
i n f o
Article history: Received 13 December 2017 Received in revised form 2 March 2018 Accepted 6 March 2018
a b s t r a c t Computational fluid dynamics (CFD) simulations are conducted to explore flow inside a channel with staggered semi-circular ribs attached to both the inner and outer annular walls. Large eddy simulation (LES) with the wall-adapting local eddy-viscosity (WALE) sub-grid scale (SGS) model is used to simulate turbulent flow in the computational domain. The flow channel is modeled on the test section of a simplified helical coil steam generator (HCSG) design. Owing to their compactness and higher heat transfer coefficient in comparison to straight tube steam generators, HCSGs are employed in recent designs of advanced light water reactors in the nuclear industry. A Reynolds number of 13,900 is considered based on the mean inlet velocity, channel height, and water properties. Results from LES are compared with both an unsteady Reynolds-Averaged Navier-Stokes (URANS) simulation and experimental data acquired from the particle image velocimetry (PIV) method. Flow features observed inside the channel are presented by visualizing contours of velocities, vorticity, and Reynolds stresses on a monitoring plane. Time-averaged profiles of velocities and Reynolds stresses at five monitoring lines are also presented to validate the flow field resolved by LES with PIV. Additionally, vortical structures of the flow are captured using the Q-criterion. This study aims at providing flow characteristics of a specifically designed flow channel and acquiring the scalability of the numerical experiment. Ó 2018 Elsevier Ltd. All rights reserved.
1. Introduction A helical coil heat exchanger (HCHX) is a type of shell-and-tube heat exchanger characterized by its compact structure and high heat transfer coefficient with relatively low pressure drop [1]. Thus, it has been employed in various types of engineering applications such as chemical process plants, cryogenic systems, heat recovery processes, and other energy intensive industries [2–4]. Notably, a recent design of advanced light water reactor (LWR) adopted a HCHX as a steam generator [5]. This type of nuclear reactor is emerging as a reflection of demands for increased scalability, lower initial investment, and greater operational safety than conventional large reactors [6]. They are small to medium sized integral type pressurized water reactors (PWRs) in which the helical coil steam generator (HCSG) is placed inside the same vessel as the reactor core and immersed in the coolant which is circulated by natural convection [7]. To ensure the safe and reliable operation of the reactor, investigation of the characteristics of the flow field around the helical coil bundle is essential. Understanding how the flow field around coil bundle affects the pressure distribution ⇑ Corresponding author. E-mail addresses: [email protected] (S.J. Lee), [email protected] (S. Lee), [email protected] (Y.A. Hassan). https://doi.org/10.1016/j.ijheatmasstransfer.2018.03.026 0017-9310/Ó 2018 Elsevier Ltd. All rights reserved.
on the helical coil surface is crucial to determine the structural and thermal safety of the reactor pressure vessel (RPV) [8]. Most of the numerous investigations on HCHX are concerned with the convective heat transfer and pressure drop in a coil tube for the tube-side and outside of a coil tube, which is simplified as a cylindrical tube array for the shell-side [9–16]. However, it is rare to find studies of detailed flow around a realistic HCHX model. An experimental test section simplifying the practical helical coil steam generator (HCSG) model [17] was constructed by Texas A&M University to examine the fluid phenomena around the shellside helical coil bundle [18]. Although the test facility was designed to simulate the natural circulation of primary coolant inside the reactor vessel, the driving force in the test facility is provided by a hydraulic pump. A honeycomb structure was installed prior to the inlet of the test section to make the flow as uniform as possible under the assumption of isothermal and uniform downward inlet velocity. To be specific, the HCSG has multiple concentric layers of helical coils in the radial direction, and each layer has a spiral direction different from the adjacent layers. To simplify the geometry, two adjacent layers were sliced in the radial direction. Further simplification is achieved by cutting the annular section into 1/25 sectors. The acrylic test section was fabricated based on this geometry, as illustrated in Fig. 1. Consequently, the flow channel
706
S.J. Lee et al. / International Journal of Heat and Mass Transfer 123 (2018) 705–717
Fig. 1. Simplified HCSG model for the experimental test section.
inside the test section has concentric annular-shaped walls and semi-circular ribs on each wall surface staggered with those on the other wall. The rib represents the coil surface of the HCSG. Previously, an unsteady Reynolds-Averaged Navier-Stokes (URANS) simulation was carried out using the flow channel extracted from the simplified test section as the computational domain for comparison with the experimental results [19]. In the URANS simulation, the time-averaged static pressure drop between the inlet and outlet of the test section showed good agreement with the experimental result. Typically, URANS simulation produces reasonable averaged statistics, but it is not appropriate to discern detailed instantaneous fields. To accurately resolve the instantaneous field and detailed flow structures, a small-scale resolving turbulence model is required. The most accurate method for simulating turbulent flow is direct numerical simulation (DNS). However, DNS is still restricted to low Reynolds numbers, which differ considerably from most practical engineering applications. Although the computing power has dramatically improved in the last few decades owing to high performance computer clusters, DNS remains as a computationally voracious method. For this circumstance, the LES approach is a realistic alternative with which investigates instantaneous turbulent flow fields. In this paper, the LES results of the simplified HCSG test section are found to provide reliable numerical data. The flow field resolved by LES was validated by the experimental PIV measurement. The flow inside the channel was studied by analyzing the velocity field. The vortex identification method was also employed to identify vortical structures that affect the pressure distribution of the channel. In addition, the static pressure drop of the test section was calculated and compared with the experimental data. 2. Numerical details The numerical experiment was performed on the flow channel inside the simplified HCSG test section using the finite volume method (FVM) based commercial CFD package, STAR-CCM+ ver.11.06-R8. To resolve the turbulent flow field of the computational domain, LES models the flow field at a scale smaller than the grid size. 2.1. Computational domain and boundary conditions The computational domain extracted from the test section is shown in Fig. 2. This geometry has semi-circular ribs attached to both concavely and convexly curved walls, each at an angle, and the ribs on one wall are staggered with those on the other wall. The diameter of the ribs is D = 15.9 mm. The height of the channel
at the inlet is h = 23.6 mm, and the gap between the ribs on both walls is d = 7.7 mm. A uniform velocity inlet and a constant pressure outlet were applied as boundary conditions of the computational domain. inlet = 0.525 m/s, corresponds to The value imposed at the inlet, u a volume flow rate of 3.15 103 m3/s (50 gpm). The resulting stream = mean streamwise velocity at the ribbed region was u 1.601 m/s. The non-uniform time varying inlet condition, discussed in Section 4.6, caused no significant inlet effect at the monitoring position of the test section; thus, a uniform velocity profile was applied to the inlet. For the outlet condition, a static pressure of 0 Pa was imposed. The density and dynamic viscosity of the working fluid are those of water at room temperature, which are q = 997.561 kg/ m3 and l = 8.887 104 kg/m s, respectively. The Reynolds number based on the mean inlet velocity, channel height at the inlet, inlet/l 13,900. and properties of water is Reh = qhu 2.2. Computational setup The unsteady isothermal incompressible fluid flow field inside the simplified HCSG test section was calculated adopting LES as a turbulence model. The SIMPLE algorithm was used to solve the discretized Navier-Stokes equations. The second-order implicit scheme was applied to discretize the governing equations temporally. The integral time scale and the time step were T = 1 s and Dt = 104 s, respectively. The bounded-central difference scheme [20] was employed to discretize the convection term with an upwind blending factor of r = 0.0 to reduce numerical diffusion. The hybrid Gauss-least squares (LSQ) method was used to reconstruct the gradient terms. The number of iterations per time step was 20, and the convergence of the solution at each time step had residuals below 104. The computational domain was equally divided into three sub-parts, cutting the domain using the x-y plane for vortical structure monitoring. Neighboring boundaries of sub-parts are nonconformal mesh connected with the interface. To ensure that a fully developed field inside the computational domain was used in the analyses, the initial 0.2 s flow fields were excluded. 2.3. Sub-grid scale model for LES The wall-adapting local eddy-viscosity (WALE) sub-grid scale (SGS) model was used to model the scale smaller than the grid size j , where the overbar i u [21]. The SGS stress tensor, sij ¼ ui uj u means filtered quantity, can be modeled based on an eddyviscosity assumption as follows:
1 3
sij skk dij ¼ 2mt Sij
ð1Þ
S.J. Lee et al. / International Journal of Heat and Mass Transfer 123 (2018) 705–717
707
Fig. 2. Computational domain of the HCSG test section.
i =@xj þ @ u j =@xi Þ is the filtered strain rate tensor where Sij ¼ ð1=2Þð@ u and mt is the SGS eddy-viscosity. To acquire the local eddy-viscosity, the square of the velocity gradient tensor which considers both the strain and the rotation rate was used to formulate the WALE model. In the WALE model, SGS eddy viscosity is defined as mt = (CwD)2Sw where Cw is a model coefficient, D the grid filter width, and Sw a deformation parameter defined as
3=2 Sdij Sdij Sw ¼ 5=4 Sij Sij 5=2 þ Sd Sd ij ij
ð2Þ
where Sdij is the traceless symmetric part of the square of the velocity gradient tensor.
Sdij ¼
i @ u k @ u j @ u k k @ u l 1 @u 1 @u dij þ 2 @xk @xj @xk @xi 3 @xl @xk
lent flow: energy-containing range, inertial subrange, and dissipation range. It is expected that LES can resolve the inertial subrange in the computational domain. That being so, the velocity component of the y-directional streamwise was monitored at the midpoint, represented by the yellow1 dot in Fig. 3, which is 70.9 mm in the y-direction from the top of the first semi-circular rib on the convex wall. The instantaneous value of the streamwise velocity component at this monitoring point was acquired at each time step in the 0.5 s period. The fast Fourier transform (FFT) was then taken from this velocity fluctuation and is illustrated in Fig. 4. It is observed that the power density spectrum covered the distinguishing 5/3 slope which represents the inertial subrange. Therefore, it can be assumed that the grid system has sufficient resolution to resolve the inertial subrange.
ð3Þ
The salient features of the WALE model are that it yields proper near-wall scaling for eddy-viscosity without requiring the Van Driest damping function, and it is less sensitive to the model coefficient Cw [21]. 2.4. Grid system A nested Cartesian grid system with ten layers of prism mesh was constructed on the extracted computational domain, as illustrated in Fig. 3. Since the nested grid system uses predominantly regular hexahedral meshes, high quality cells with minimum cell skewness fill the computational domain. The base size of the nested grid system was 1.5 mm. The grid on the ribbed region was refined with 25% of the base size as 0.375 mm to resolve detailed flow structures. The first off-wall grid distance, Dy = 103 mm, was selected to make Dy+ 1 after performing preliminary simulations. There were approximately 37 million cells inside the domain. Kolmogorov’s 5/3 law with Taylor’s hypothesis was invoked to check the validity of this grid system for the present numerical study [22]. There are three ranges in the energy spectrum of turbu-
2.5. URANS simulation In this study, the Reynolds stress transport (RST) model with the linear pressure strain two-layer model was selected as a turbulence model for the URANS simulation. The same grid system as that in LES was used in the URANS simulation for comparison. For the temporal and spatial discretization schemes, the secondorder implicit scheme with 104 s time-step, and the MUSCL/CD scheme were applied, respectively. The SIMPLE algorithm was used as a pressure correction. In addition, all boundary conditions applied to the computational domain were the same as in LES. 3. Experimental data The experimental data presented in this paper are the result of PIV analysis [23] on the acrylic test section shown in Fig. 2. The schematic of the simplified HCSG test facility for the PIV measurement is illustrated in Fig. 5. The PIV system consisted of a highspeed camera and 532 nm continuous laser with 10 W power. 1 For interpretation of color in Figs. 3, 6 and 9, the reader is referred to the web version of this article.
708
S.J. Lee et al. / International Journal of Heat and Mass Transfer 123 (2018) 705–717
concave wall side Ǹ
convex wall side Ǹ
Fig. 3. Section view of the grid system on the computational domain.
Fig. 4. FFT of streamwise velocity fluctuation of a point inside the channel.
Fig. 6. Monitoring plane and lines for data acquisition.
Fig. 5. Schematic of the simplified HCSG test section facility for the PIV measurement.
The frame rate of the PIV picture was 3500 fps. The size of the image was 332 1024 pixels with a magnification factor of 70.9 lm/pixel. Silver-coated trace particles with a 13 lm mean cutoff diameter and density of 1.6 g/cm3 were seeded to water. In the post processing of PIV analysis, the adaptive FFT-correlation [24] with three-step interrogation window sizes of 32 16, 16 8, and 8 4 pixels was utilized. The uncertainty of the PIV measurement in accordance with the ITTC guide [25] was uc = 6.6% of the mean streamwise velocity. In addition to the PIV measurement, the static pressure drop across the test section was obtained using pressure transducers installed at two red dot positions as shown in Fig. 6.
S.J. Lee et al. / International Journal of Heat and Mass Transfer 123 (2018) 705–717
PIV
URANS
709
LES
v [m/s] Fig. 7. Instantaneous streamwise velocity, v, contour.
PIV
URANS
LES
4. Results and discussion In the CFD simulation, the same monitoring positions were set inside the computational domain to compare the simulation result with the experimental data, as shown in Fig. 6. The monitoring plane was made by slicing the x-y plane of the fluid domain 62.7 mm from the right-side wall. The velocity field in the boxed plane was resolved, and various properties were calculated to illustrate each contour. The velocity components are defined as ui = huii + ui0 , where the angled bracket and prime indicate the mean and fluctuating components of velocity, respectively. In this study, over
Fig. 9. Mean streamwise velocity, hvi, profiles on lines C1 to C5.
10,000 instantaneous fields were used to obtain a mean value. There are five monitoring lines around the sixth semi-circular rib on the convex wall from the upstream position used to plot various
710
S.J. Lee et al. / International Journal of Heat and Mass Transfer 123 (2018) 705–717
PIV
URANS
LES
u [m/s] Fig. 10. Instantaneous transverse velocity, u, contour.
PIV
URANS
LES
[m/s] Fig. 11. Mean transverse velocity, hui, contour.
mean properties. In addition, the center section of the computational domain was utilized to show the three-dimensional vortical structures resolved by LES. 4.1. Streamwise and transverse velocity fields Fig. 7 illustrates the instantaneous streamwise (negative ydirection) velocity, v; and the contours resolved by PIV, URANS, and LES. As seen in the figure, both PIV and LES exhibit the characteristic features of the instantaneous velocity field, whereas the contour of the URANS simulation shows a velocity field similar to its mean field. Although the URANS simulation solves the Navier-
Fig. 12. Mean transverse velocity, hui, profiles on lines C1 to C5.
Stokes equation using a turbulence closure model with an unsteady manner, the RANS turbulence model resolves the instantaneous velocity field time-averaged by its very nature.
S.J. Lee et al. / International Journal of Heat and Mass Transfer 123 (2018) 705–717
PIV
URANS
711
LES
The mean, or time-averaged, streamwise velocity, hvi, contours are presented in Fig. 8. Contours for all three cases show the characteristic feature of the mean velocity field. In other words, there are no abrupt changes in the streamwise velocity in the observed area. Regions of relatively high magnitude of downstream velocity appear near the channel core. Specifically, regions of the highest streamwise velocity magnitude are captured right next to the rib surfaces. The cavity between two consecutive ribs of the same wall-side exhibits counter streamwise velocity due to flow recirculation. The mean streamwise velocity profiles on lines C1 to C5 depicted in Fig. 6 are plotted in Fig. 9. As indicated in the legend box on the first plot, the black solid circle symbol, black dashed line, and blue solid line are the PIV, URANS, and LES data, respectively. The 0 mm position in each plot indicates the leftmost wall of each corresponding line. The results of the LES show good agreement with the PIV near the wall and cavities and are underestimated near the channel core region. For the case of URANS, the simulation results underestimate the mean streamwise velocity compared to PIV on the whole monitoring lines. Figs. 10 and 11 illustrate the instantaneous and mean transverse x-directional velocity, u and hui, contours, respectively. Likewise, in the case of streamwise velocity, the contours for the instantaneous transverse velocity fields of PIV and LES show the characteristic feature of instantaneous field; however, URANS shows that the instantaneous field looks similar to the mean field in Fig. 11. The notable feature observed in the mean transverse velocity contour is the stagnation point, which can be distinguished by finding the point where the sign change occurs on the semi-circular rib surface facing the upstream direction. The mean transverse velocity profiles on lines C1 to C5 are shown in Fig. 12. Profiles for the LES show better agreement with the PIV results than the URANS simulation. However, for two lines, C1 and C5, the LES overestimated the mean transverse velocity compared to the PIV in close vicinity of the right region near the concave ribbed wall. For the cases of C2 and C3, the profiles of LES skewed toward the right wall.
Fig. 14. Mean Reynolds normal stress, hv0 v0 i, profiles on lines C1 to C5.
The mean transverse velocity profiles of the URANS simulation on C1, C2, and C3 deviate considerably from the PIV results. The trend of the URANS profile on C4 follows the PIV results, but its
712
S.J. Lee et al. / International Journal of Heat and Mass Transfer 123 (2018) 705–717
PIV
URANS
LES
[m/s] Fig. 15. Mean Reynolds normal stress, hu0 u0 i, contour.
magnitude is substantially underestimated. For the C5 profile, the URANS simulation matched well with the PIV data in the 5 < x < 10 mm range but deviated over other ranges. 4.2. Mean Reynolds stresses The contours of the mean Reynolds normal stresses, hv0 v0 i and hu0 u0 i, are shown in Figs. 13 and 15, respectively. In addition, the corresponding profiles for the mean Reynolds normal stresses captured on lines C1 to C5 are plotted in Figs. 14 and 16, respectively. As seen in the contours, the LES results have better agreement with those of the PIV than the URANS simulation. The regions of local maximum values of Reynolds normal stresses of the PIV and LES contours are similar. However, in the case of the URANS simulation, regions of local maximum values deviate from the others. Observation of the profiles of these mean Reynolds normal stresses on lines C1 to C5 shows obvious differences between LES and URANS. The LES profiles agree well with the PIV data, whereas the profiles of the URANS simulation are underestimated compared to PIV. Fig. 17 shows the contours of the mean Reynolds shear stress, hu0 v0 i. In addition, the corresponding profiles on lines C1 to C5 are plotted in Fig. 18. Regions of local maximum values are observed behind the ribs on both the PIV and LES contours, whereas the local maximum regions on the URANS contour show different aspects to the other cases. The profiles of the mean Reynolds shear stress for LES follow the PIV trends well, and the magnitudes of the values are underestimated at several locations. For the URANS simulation, the profiles partially follow the trends of PIV; however, it is observed that the profiles deviate from the PIV data. 4.3. Vorticity The instantaneous and mean vorticity, xz and hxzi, contours are illustrated in Figs. 19 and 20, respectively. As seen in the previous contours of velocities, the vorticity contours of URANS show no difference between the instantaneous and mean fields, and the
Fig. 16. Mean Reynolds shear stress, hu0 u0 i, profiles on lines C1 to C5.
instantaneous vorticity contours of PIV and LES show the distinctive features of the instantaneous field. In the case of mean vorticity contours, the three contours show similar patterns of mean
S.J. Lee et al. / International Journal of Heat and Mass Transfer 123 (2018) 705–717
PIV
URANS
713
LES
[m/s] Fig. 17. Mean Reynolds shear stress, hu0 v0 i, contour.
vorticity. Relatively high vorticity magnitude regions are observed behind the top surface of the ribs in the downstream direction. Fig. 21 shows the profiles of mean vorticity on lines C1 to C5. Both URANS and LES have good agreement with PIV, other than in proximity to the wall. For all five lines, LES has a high magnitude of mean vorticity near the wall, because LES can resolve the velocity gradient better in close vicinity to the wall than PIV. 4.4. Vortical structures High magnitude of vorticity is not a sufficient condition for the existence of a strong vortex. Therefore, Q-criterion, which is a vortex identification method suggested by Hunt et al. [26], was utilized to identify vortical structures inside the computational domain. Indeed, Q is the second invariant of the velocity gradient tensor, @ui/@xj. Assuming an incompressible flow, Q is defined as
Q
1 2
@ui @ui @uj 1 @ui @uj 1
Xij 2 Sij 2 ¼ ¼ 2 @xj @xi 2 @xi @xj @xi
ð4Þ
where Xij and Sij are the rotation and strain rate tensors, respectively, and |||| is the Frobenius norm. The Q-criterion utilizes the positive value of Q to represent the vortex in a fluid region. In other words, it identifies the vortical structure, finding a fluid region in which the rotation rate is larger than the strain rate of the velocity field. Fig. 22 shows vortical structures resolved by LES inside the computational domain using the Q-criterion. Velocity magnitude contours on the iso-surface of each Q value are presented. Observation of the figure indicates that vortical structures move down in the streamwise direction after they are generated next to the rib surface near the channel core region. This phenomenon is clearly observed on the contour of Q on the center plane, as shown in Fig. 23. Some vortices are captured inside the cavity between two consecutive ribs, and the others travel downstream.
Fig. 18. Mean Reynolds shear stress, hu0 v0 i, profiles on lines C1 to C5.
4.5. Static pressure drop To measure the static pressure drop across the test section, pressure transducers were installed at both the inlet and outlet
pressure taps of the test section shown in Fig. 6. The pressure drop, Dp = pin pout, was 2855 ± 8.4 Pa which was obtained from the ensemble average of eighteen data sets. Each data set has over
714
S.J. Lee et al. / International Journal of Heat and Mass Transfer 123 (2018) 705–717
URANS
PIV
z
LES
[s-1] Fig. 19. Instantaneous vorticity, xz, contour.
PIV
URANS
LES
< z> [s-1] Fig. 20. Mean vorticity, hxzi, contour.
10,000 samples. In the CFD simulations, fluctuating pressure drops were averaged from each T = 1 s simulation. URANS and LES calculated the pressure drop as 2343 Pa and 2648 Pa, respectively. Both CFD simulations yielded lower pressure drops than the experiment, which are 17.9% and 7.3%, for URANS and LES, respectively. LES predicted the pressure drop measured from the experiment better than the URANS simulation. Fig. 21. Mean vorticity, hxzi, profiles on lines C1 to C5.
4.6. Discussion In this study, LES resolved the flow field better than the URANS simulation, as seen in both the instantaneous and mean flow fields.
Since the geometry of the flow channel is not symmetric and the Reynolds number, Reh = 13,900, is in the turbulent regime, the instantaneous flow field inside the domain showed strong
S.J. Lee et al. / International Journal of Heat and Mass Transfer 123 (2018) 705–717
715
metric difference between the CAD model used to extract the computational domain and the actual flow channel of the test section can be considered a source of uncertainty. The CFD simulations in this study used the uniform inlet boundary condition. However, the experiment showed that the upper
Q = 4.0 105 s -2
Q = 2.0 105 s -2
Q = 1.5 105 s -2
Q = 1.0 105 s -2
Q [s −2 ] ur | V | [m / s]
Fig. 23. Contour of Q value resolved by LES on the center plane of the middle section.
Fig. 22. Vortical structure resolved by LES using the iso-surface of Q-criterion with various values.
unsteadiness. Therefore, to examine the flow phenomena for this flow properly, selection of an appropriate turbulence model that can resolve small scale unsteady features in the flow is essential. Relatively large discrepancies between the PIV and the LES profiles are observed on line C3. In raw pictures of the PIV measurement, a shadow line induced by optical distortion due to the complex geometry of the acrylic test section on the laser sheet appears right below line C3. Consequently, it is assumed that the shadow has an effect on the PIV analysis near line C3, and the resolved field deviates from the field resolved by unaffected PIV pictures. This shadow line effect is clearly observed in Figs. 8 and 13. In addition, PIV is not able to capture the high velocity gradient very close to the wall, as seen in Fig. 21, since this method has limited resolution capabilities in close proximity to the wall. The geo-
LES
LES w/ SEM
|V| [m/s] Fig. 24. Velocity magnitude contours near the inlet boundary for the uniform velocity and synthetic-eddy-method cases.
716
S.J. Lee et al. / International Journal of Heat and Mass Transfer 123 (2018) 705–717
Fig. 25. Mean streamwise and transverse velocities, hv i and hui; mean Reynolds normal and shear stresses, hv0 v0 i, hu0 u0 i, and hu0 v0 i; and mean vorticity hxzi profiles on line C5 for LES with uniform inlet (blue solid line) and for LES with SEM (red dashed line). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
part of the test section generated a non-uniform inlet velocity profile with a turbulent intensity (TI) of 0.107. The monitoring positions are located after the third rib in the streamwise direction. It is necessary to ascertain whether the inlet condition affects the flow field at the monitoring positions. There is a numerical method to generate realistic inlet flow conditions for LES; it is called the synthetic-eddy-method (SEM) and was suggested by Jarrin et al. [27]. The turbulent intensity measured at the inlet position of the experimental test section was TI = 0.107, and the honeycomb mesh size can be considered as the eddy turbulent length scale, le = 3.97 mm. The SEM is applied to the inlet boundary condition of the computational domain using these two values. Fig. 24 illustrates the velocity magnitude contours on the center plane of the computational domain for both the uniform velocity and the SEM. Unlike in the uniform velocity inlet case, velocity fluctuations are observed in the SEM case owing to the existence of synthetic eddies. Comparisons of various mean properties between LES with a uniform inlet and LES under the SEM condition, are plotted on line C5 in Fig. 25. As seen in the figure, even if the inlet conditions for both cases are different, results for the mean profiles show no significant difference at the region of interest. Therefore, it can be assumed that the flow fields resolved by the experimental PIV method with the unsteady non-uniform inlet and LES with a uniform inlet condition are comparable and reliable. According to the review paper of Païdoussis [28] on vortex shedding cylinder arrays in cross flow, the flow pattern after the third row remains almost unchanged owing to the influence of the turbulence gener-
ated by the cylinder placed at an upstream position. Since this study deals with channel flow, which can be seen as a part of cross flow of the helical coil tube bundle, a similar phenomenon is expected.
5. Summary and conclusion A numerical investigation of the turbulent flow field inside a channel with a specific configuration of semi-circular ribs that simulates the simplified HCSG geometry was conducted using LES. The results of the LES were validated by experimental data. In particular, the WALE SGS model was used to resolve the flow in the computational domain with a Cartesian nested grid system. The grid resolution was sufficient to resolve the inertial subrange of turbulence. In addition, a URANS simulation with an RST model was performed on the same grid system to provide comparison data. Instantaneous and mean flow fields resolved by PIV, URANS, and LES were compared at the region of interest. Owing to the nature of each turbulence model, LES resolved the unsteady flow features better than the URANS simulation. Moreover, the profiles of the mean velocities and Reynolds stresses showed that LES correlated better with the PIV data. Vortical structures resolved by LES were visualized using the Qcriterion. Vortices were induced next to the rib surface near the channel core region, and then they detached from the rib wall and were carried away by the bulk flow downstream. Contours
S.J. Lee et al. / International Journal of Heat and Mass Transfer 123 (2018) 705–717
of the Q value on a monitoring plane enabled to locate the vortex cores. Additionally, some vortical structures, which were resolved in the cavity between two consecutive ribs, were observed. The static pressure drop calculated by both numerical and experimental studies was presented. LES predicted the pressure drop of the test section with higher accuracy than the URANS simulation when compared to the experimental data. Several sources of discrepancy between the experiment and the CFD simulation were discussed. On the experimental side, the uncertainty of the raw pictures of PIV and the geometric difference were expected to affect this comparison. On the CFD side, the nonrealistic uniform inlet boundary condition applied to the computational domain was examined as a possible cause of discrepancy. However, the profiles of the mean flow fields on the monitoring position showed no significant difference to the case of the more realistic inlet boundary condition generated by the SEM. This leads to the conclusion that the inlet condition does not significantly affect the results of the CFD simulation at the region of interest. This is also applicable to the experimental facility; thus, the inlet condition is not considered a valid cause of discrepancy between the experimental and CFD results. For future research, a more realistic test section geometry will be constructed, and the extracted computational domain will capture the actual HCSG flow phenomena more accurately. Conflict of interest There are no known conflicts of interest associated with this publication. Acknowledgements This study is supported by the Integrated University Program (IUP) graduate fellowship, the U.S. Department of Energy Office of Nuclear Energy (DOE-NE). The authors acknowledge Texas A&M University High Performance Research Computing for providing computing power to perform the CFD simulations of this study. References [1] R.K. Patil, B.W. Shende, P.K. Ghosh, Designing a helical-coil heat exchanger, Chem. Eng. 13 (12) (1982) 85–88. [2] M.E. Ali, Experimental investigation of natural convection from vertical helical coiled tubes, Int. J. Heat Mass Transf. 37 (4) (1994) 665–671. [3] R.C. Xin, A. Awwad, Z.F. Dong, M.A. Ebadian, H.M. Soliman, An investigation and comparative study of the pressure drop in air-water two-phase flow in vertical helicoidal pipes, Int. J. Heat Mass Transf. 39 (4) (1996) 735–743. [4] D.G. Prabhanjan, G.S.V. Raghavan, T.J. Rennie, Comparison of heat transfer rates between a straight tube heat exchanger and a helically coiled heat exchanger, Int. Commun. Heat Mass Transf. 29 (2) (2002) 185–191.
717
[5] International Atomic Energy Agency, Status of advanced light water reactor designs, IAEA-TECDOC-1391, Vienna, Austria, 2004. [6] International Atomic Energy Agency, Instrumentation and control systems for advanced small modular reactors, IAEA Nuclear Energy Series No. NP-T-3.19, Vienna, Austria, 2017. [7] International Atomic Energy Agency, Natural circulation in water cooled nuclear power plants: Phenomena, models, and methodology for system reliability assessments, IAEA-TECDOC-1474, Vienna, Austria, 2005. [8] International Atomic Energy Agency, Design safety considerations for water cooled small modular reactors incorporating lessons learned from the Fukushima Daiichi accident, IAEA-TECDOC-1785, Vienna, Austria, 2016. [9] A.N. Dravid, K.A. Smith, E.W. Merrill, P.L.T. Brian, Effect of secondary fluid motion on laminar flow heat transfer in helically coiled tubes, AICHE J. 17 (5) (1971) 1114–1122. [10] S.V. Patankar, V.S. Pratap, D.B. Spalding, Prediction of laminar flow and heat transfer in helically coiled pipes, J. Fluid Mech. 62 (3) (1974) 539–551. [11] S.A. Berger, L. Talbot, L.S. Yao, Flow in curved pipes, Annu. Rev. Fluid Mech. 15 (1983) 461–512. [12] K. Futagami, Y. Aoyama, Laminar heat transfer in a helically coiled tube, Int. J. Heat Mass Transf. 31 (2) (1988) 387–396. [13] P. Naphon, S. Wongwises, A review of flow and heat transfer characteristics in curved tubes, Renew. Sust. Energy Rev. 10 (5) (2006) 463–490. [14] S. Kakaç, H. Liu, A. Pramuanjaroenkij, Heat Exchangers: Selection, Rating, and Thermal Design, CRC Press, 2012. [15] P.V. Gilli, Heat transfer and pressure drop for cross flow through banks of multistart helical tubes with uniform inclinations and uniform longitudinal pitches, Nucl. Sci. Eng. 22 (3) (1965) 298–314. [16] S. Majumdar, A. Moisseytsev, K. Natesan, Assessment of the Next Generation Nuclear Plant Intermediate Heat Exchanger Design, Argonne National Laboratory EXT-08/32, 2008. [17] D.K. Williams, D.P Fassett, B.J. Webb, W.J. Bees, A.S. Kruskamp, Helical Coil Steam Generator, U.S. Patent 20120111287 A1, Issued Apr. 12, 2013. [18] M. Delgado, S.J. Lee, S. Lee, Y.A. Hassan, Pressure drop measurements across a helical coil steam generator configuration, in: Proceedings of International Topical Meeting on High Temperature Reactor Technology (HTR 2016), Las Vegas, Nevada, USA, November 6–10, 2016. [19] S.J. Lee, S. Lee, Y.A. Hassan, Numerical investigation of the flow inside the helical coil steam generator test section, in: Proceedings of International Topical Meeting on High Temperature Reactor Technology (HTR 2016), Las Vegas, Nevada, USA, November 6–10, 2016. [20] B.P. Leonard, The ULTIMATE conservative difference scheme applied to unsteady one-dimensional advection, Comput. Meth. Appl. Mech. Eng. 88 (1) (1991) 17–74. [21] F. Nicoud, F. Ducros, Subgrid-scale stress modelling based on the square of the velocity gradient tensor, Flow Turbul. Combust. 62 (3) (1999) 183–200. [22] S.B. Pope, Turbulent Flows, Cambridge University Press, 2000. [23] S. Lee, M. Delgado, C. Guilbault, Y.A. Hassan, Visualization of flow interaction in a channel flow produced by a pair of helically coiled rods, Transactions of ANS Annual Meeting, San Francisco, California, USA, June 11–15, 2017. [24] W. Thielicke, E.J. Stamhuis, PIVlab - towards user-friendly, affordable and accurate digital particle image velocimetry in MATLAB, J. Open Res Software 2 (1) (2014) e30, https://doi.org/10.5334/jors.bl. [25] ITTC, ITTC recommended procedures and guidelines - uncertainty analysis particle image velocimetry, Technical Report, 2008. [26] J.C.R. Hunt, A.A. Wray, P. Moin, Eddies, streams, and convergence zones in turbulent flows, in: Proceedings of the 1988 Center for Turbulence Research Summer Program (CTR-S88), Stanford, California, USA, 1988. [27] N. Jarrin, S. Benhamadouche, D. Laurence, R. Prosser, A synthetic-eddy-method for generating inflow conditions for large-eddy simulations, Int. J. Heat Fluid Flow 27 (4) (2006) 585–593. [28] Païdoussis, A review of flow-induced vibrations in reactors and reactor components, Nucl. Eng. Des. 74 (1) (1983) 31–60.