Abstract
A core catcher has been developed to maintain the integrity of nuclear reactor containment from molten corium during a severe accident. It uses a two-phase natural circulation for cooling molten corium. Flow in a typical core catcher is unique because (i) it has an inclined cooling channel with downwards-facing heating surface, of which flow processes are not fully exploited, (ii) it is usually exposed to a low-pressure condition, where phase change causes dramatic changes in the flow, and (iii) the effects of a multidimensional flow are very large in the upper part of the core catcher. These features make computational analysis more difficult. In this study, the MARS code is assessed using the two-phase natural circulation experiments that had been conducted at the CE-PECS facility to verify the cooling performance of a core catcher. The code is a system-scale thermal-hydraulic (TH) code and has a multidimensional TH component. The facility was modeled by using both one- and three-dimensional components. Six experiments at the facility were selected to investigate the parametric effects of heat flux, pressure, and form loss. The results show that MARS can predict the two-phase flow at the facility reasonably well. However, some limitations are obviously revealed.
1. Introduction
The Fukushima nuclear power plant accident in 2011 was caused by the damage to the power supply system of the safety system due to the tsunami following the great earthquake. The unprecedented long-term station blackout evoked the reactor core cooling function to be completely lost, resulting in a core melting and subsequent reactor vessel failure. Since the accident, the importance of a passive cooling system has become more prominent. In particular, the so-called core catcher has been actively developed as one of the major accident coping facilities against a severe accident.
The core catcher aims at protecting the integrity of a reactor containment building from corium when a core is molten and the molten corium leaks out of the reactor vessel. Figure 1 shows a conceptual diagram of a core catcher [1, 2]. It consists of a carbon steel body, a sacrificial material layer, and an inclined cooling channel. The carbon steel body is a place where discharged corium from the reactor vessel is rearranged and stabilized. It prevents molten core from interacting with concrete at the bottom of the reactor containment building. The sacrificial material placed in the surface of the body controls the melt properties and modifies melt conditions favorable to corium retention [3]. Beneath the body, there is an inclined cooling channel to cool down the molten corium. This cooling channel is connected to a downcomer to establish a passive closed-loop cooling system.
Decay heat from the corium is removed by conduction and convection through the inclined cooling channel. Convective heat transfer also occurs from the corium to the coolant above it. Coolant in the cooling channel starts boiling due to the heat transfer from the corium and moves upwards by buoyancy. The coolant flowing out from the channel is mixed with the water inside the reactor cavity and thus cooled. The cooled water again flows into the downcomer and the cooling channel, which closes a natural circulation loop.
The researches for core catcher design have been focused mainly on the verification of its performance to reinforce the safety of advanced light water reactors (LWRs). AREVA carried out critical heat flux experiment for the core catcher of European PWR and verified that the measured maximum heat flux at the heating surface under the condition of low mass flow rate of 0.1~0.2 kg/s did not exceed critical heat flux [4]. Toshiba performed a natural circulation test using a fan-shaped test facility and confirmed a cooling performance of a round-shaped core catcher for an advanced boiling water reactor [5]. They also showed that the results of computational analysis are in a good agreement with the experiments. Korea Maritime University conducted a flow boiling test in an inclined channel with downwards-facing heating wall. They measured wall temperature and bubble size to investigate heat transfer and two-phase flow phenomena in the inclined channel [6]. KAERI has conducted the conceptual design of a core catcher for an advanced PWR and constructed a test facility named the CE-PECS [7]. They demonstrated that the prototypic core catcher has a sufficient heat removal capacity. The facility, which is similar to the right half of Figure 1, consists of an inclined cooling channel and a downcomer, a cooling water tank to simulate the water in the reactor cavity, and heating blocks to simulate the decay heat from corium. The heating blocks are installed on the upper wall of the inclined cooling channel.
While many experiments have been conducted for core catcher design and verification, few computational studies have been found in the literatures. The flow in a core catcher, including that in the CE-PECS facility, is unique because it has an inclined cooling channel with downwards-facing heating surface, of which flow processes are not fully exploited. Furthermore the core catcher is usually exposed to a low-pressure condition, where the phase change causes dramatic changes in a two-phase flow, and the effects of a multidimensional two-phase flow are very large in the cooling water tank. These features make computational analysis more difficult.
Computational analysis capability of two-phase flows in a core catcher is essential for the design and safety analysis of advanced nuclear reactors. The primary objective of this study is to assess the MARS code against the two-phase natural circulation experiments using the CE-PECS facility. The MARS code is a system-scale thermal-hydraulic (TH) code [8, 9]. It has been developed from the consolidated version of the RELAP5/MOD3 and COBRA-TF codes [10] and has a multidimensional TH component as well as one-dimensional components [11]. This code assessment is focused on the prediction of a two-phase natural circulation flow at low-pressure conditions at the facility, of which key phenomena involve a three-dimensional two-phase flow in the cooling water tank, a two-phase flow and heat transfer from the downwards-facing heating surface in an inclined channel, and the balance between the buoyance and frictional pressure drop along the circulation loop. Six experiments at the facility were selected to investigate the parametric effects of heat flux, pressure, and form loss on natural circulation flow. For accurate simulation of a multidimensional two-phase flow, the three-dimensional TH component was used for the cooling water tank. The results of the MARS code are presented and its limitations are discussed. Some results are also compared with those of a component-scale three-dimensional code to find out an efficient way of a multiscale analysis using the two different scale codes.
2. CE-PECS Test Facility and Its Modeling
In this section, the test facility and selected experiments for code assessment are briefly introduced first. Then, the MARS input model for the facility are presented, which is refined by reflecting the results of preliminary calculations.
2.1. CE-PECS Facility and Selected Experiments
The CE-PECS is an experimental facility to study the cooling performance of a core catcher designed for an advanced 1400 MWe PWR [7]. It consists of a coolant circulation system to simulate natural circulation of cooling water, a water tank to simulate the water inside the reactor cavity which surrounds the reactor vessel, and heating blocks to simulate decay heat of corium as shown in Figure 2. A coolant circulation system is made up of an inclined cooling channel and a downcomer. The inclined cooling channel has a square cross section with 0.3 m in width and 0.1 m in gap. It is inclined 10 degrees horizontally. As illustrated in Figure 2, nine studs acting as supports for the core catcher body are located in the channel, which evoke additional flow resistance. The downcomer, which is a circular pipe with 0.1 m diameter, is connected to the inlet of the cooling channel from the lower part of the water tank so that coolant can be continuously circulated. The height of the downcomer is 2.35 m from inlet of the cooling channel. An orifice is placed in the downcomer to match the form loss of the CE-PECS loop with that of the real core catcher system. The water tank is 0.5 m × 0.5 m in cross-sectional area and ~8 m in total height. This height was determined to simulate the cooling water which is supplied to the core catcher from in-containment refueling water storage tank during a severe accident. When water in the tank is saturated to the top and the steam is discharged to its free surface, the steam is condensed by a heat exchanger and the condensate flows back to the tank.
The heating blocks are installed on the upper wall of the cooling channel. They are partitioned into seven zones from A to G as shown in Figure 3. The heating blocks A through E are downwards-facing heating walls. The heating blocks F and G are installed in the vertical part of the channel, forming one-side heating walls. Different heat fluxes can be imposed to each zone. The heating blocks are 0.3 m in width and 0.07 m in thickness. A total of 114 K-type thermocouples were installed to measure the temperature in the heating blocks.
(a)
(b)
The measuring parameters are the power supplied to the heating blocks, mass flow rate at the downcomer, temperature at the cooling channel inlet, and outlet, pressure, and temperature in the heating blocks. The measurement gauge errors of the total power, temperature in the coolant, temperatures in the heating blocks, pressure, and mass flow rate are bounded within ±1.0 kW, ±1.0 K, ±2.2 K, ±0.001 bar, and ±0.03 kg/sec, respectively [7].
Initially the cooling channel is filled with water at a predetermined temperature up to a specified level at the water tank. Then, the heater power is gradually increased up to a preset level as a function of time for each test. When the natural circulation reaches a nearly quasi-steady state, the power supply is terminated and the test ends.
Natural circulation flow rate in the loop is dependent on the heat flux, pressure, subcooling, friction and form loss, and so forth. To confirm whether the MARS code can represent these effects appropriately, we selected 6 cases among the series of the CE-PECS experiments. Table 1 shows the initial and boundary conditions for the selected tests:(i)Cases and 4-6 are used to compare the effect of the pressure. As listed in Table 1, the water levels are different, leading to different pressures at the inlet of the cooling channel.(ii)Cases 4-6, 4-7, and 4-8 are selected to examine the effect of heat flux imposed to heating blocks.(iii)Cases 3-1 and 4-1 are used to compare the effect of the loss coefficient of the orifice. Both the loss coefficient and pressure are different in the two cases. Because the two parameters can affect the natural circulation flow, sensitivity calculations for the loss coefficient are additionally performed.
2.2. MARS Input Model
The MARS code is a best-estimate system code and has a multidimensional TH component. To prepare a MARS input for the CE-PECS experiments, preliminary calculations were carried out first and the results were reflected into the input data. In this study, the MARS-KS 1.3 version is used.
2.2.1. System Nodalization
Figure 4 shows the MARS nodalization for the CE-PECS facility. It is modeled by using a combination of one-dimensional (1D) components and a three-dimensional (3D) component: the inclined cooling channel and the downcomer are modeled by a 1D pipe component for computational efficiency. However, the cooling water tank is modeled by a 3D component.
If a 1D pipe component is used for the cooling water tank, the 3D flow behavior in the tank cannot be simulated, which has a strong effect on the natural circulation behavior in the system. Figure 5 schematically shows the flow behavior in the cooling water tank in the case of the 1D and 3D models. If the 1D model is used for the tank, the flow in the water tank just can move up or down vertically and, thus, it is impossible to simulate the natural circulation and flow mixing in the water tank. In order to see the effect, preliminary calculations for Case 4-1 were performed with 1D and 3D component for the water tank. The resulting coolant temperatures at the inlet of the cooling channel are compared in Figure 6(a). It shows the 3D model predicts the temperature behavior more realistically. In the 1D calculation, the coolant in the two bottom cells of the water tank only is involved in the natural circulation, which results in an unrealistic rapid temperature increase as shown in Figure 6(a). Because relatively hot water is introduced into the cooling channel in the 1D calculation very early, boiling occurs too early. This results in a sudden increase in the mass flow rate at the channel at ~1700 s as shown in Figure 6(b). Thus, the 3D component is used for the water tank. To select an optimum 3D grid, mesh sensitivity calculations were carried out and a grid with 5 × 3 × 24 cells was adopted for the tank as shown in Figure 5(b).
(a) 1D modeling
(b) 3D modeling
(a) Temperature
(b) Mass flow rate
The heating blocks are modeled by the heat structure components. For each heat structure component, the heat flux is given as a function of time at the left boundary and the convective boundary condition is specified at the right boundary. Figure 7, for example, shows total power behavior of the heating blocks for Case 4-1. The heat fluxes on each heating surface were slightly different at each test. For the convective boundary type, the option “flat plate above fluid” was specified. The environmental heat loss during the experiment was not modeled in the MARS calculation.
2.2.2. Turbulent Mixing Model
As the 3D components of the MARS code are used, the coefficient for the turbulent mixing length model must be specified as a user input. The turbulent viscosity in the MARS code [11] is given bywhere , , and are fluid density, bulk deformation tensor, and the mixing length, respectively. Subscript means liquid or vapor phases. The default value of the mixing length is 0.1 m.
Figure 8 presents the temperature measurement locations in the water tank. Figure 9 compares the calculated coolant temperatures using the default mixing length and a reduced one with the experimental data. In the experiment, the temperature difference between the upper (TTC-WT-04) and lower (TTC-WT-01) water tank is approximately 10 K as shown in Figure 9(a). However, in the calculation with the default value, the temperature difference is approximately 3 K. This might be due to both excessive turbulent mixing and numerical diffusion. To realistically simulate the temperature difference, should be modified appropriately. Figure 9(b) compares the temperature behaviors in the water tank with a reduced mixing length of 0.01 m, which is close to a theoretical minimum, with the experimental data. This yields better agreements with the experimental data. In the subsequent calculations, the mixing length of 0.01 m is used.
(a) = 0.1 m (default)
(b) = 0.01 m (reduced)
2.2.3. Interfacial Drag
To evaluate the cooling performance of a core catcher, it is important to accurately predict the two-phase natural circulation flow rate. The interfacial drag has a direct effect on the slip between the vapor and liquid phase. If the interfacial drag is smaller, the slip increases and, in turn, the void fraction decreases. This leads to the reduction of a driving force for natural circulation, resulting in a low natural circulation flow in the loop. Preliminary calculations showed that, generally, the natural circulation flow rates are overpredicted and the void fractions at the exit of the cooling channel are underpredicted by the MARS code.
Sensitivity calculations with a reduced interfacial drag are carried out, using the multipliers of 0.9 and 0.8 for the cooling channel zone, where a two-phase flow mainly occurs. Figure 10(a) shows the flow rate with different interfacial drag multipliers. There is no change in the flow rate before ~2600 s because the void fraction is small as shown in Figure 10(b). After 2600 s, the smaller multiplier yields the lower flow rates at the inlet of the cooling channel. From these results, it is clear that the drag has a direct influence on the flow rate when the void fraction reaches a certain value. However, the underprediction of the void fraction became worse by decreasing the interfacial drag multiplier. Meanwhile, there is a negligible effect on the inlet coolant temperature according to the interfacial drag multipliers. Considering these results, the default drag model, that is, a multiplier of 1.0, was adopted.
(a) Mass flow rate
(b) Void fraction at the exit
(c) Temperature at the inlet and exit
3. The Results of Calculations and Discussions
Using the MARS input data, the six cases of the CE-PECS experiments were simulated. In the calculations, the natural circulation flow rate is a key parameter. In a single-phase natural circulation loop with a heat source and a sink, a steady-state mass flow rate can be approximated using the Boussinesq approximation as follows [12]: where and are the frictional factor and form loss factor along the closed loop, respectively. The sum symbol in the denominator means to sum up along the natural circulation loop. is the power of the heat source. For a two-phase loop, a similar but more complicated formula can be derived by using some assumptions. General features would be similar in single- and two-phase flow loops. It is clear that natural circulation flow rate depends on the driving force and resistance; heat flux is related to the former, whereas the friction and form loss are related to the latter.
3.1. The Effect of Pressure
If the pressure at the inlet of the cooling channel is reduced, the boiling point decreases and, thus, boiling occurs earlier. At the same time, because the steam density at a lower pressure is lower, buoyance force increases, resulting in an increased natural circulation. Therefore we can expect an earlier and greater flow rate by reducing the inlet pressure. For Cases 4-2 and 4-6, the water levels of 6.4 m and 8.9 m from the inlet of the cooling channel are different, which lead to different pressures of ~1.6 bar and ~1.9 bar, respectively. Figure 11 shows the natural circulation flow rates for the two cases. Both the experimental data and the MARS results show the trend reasonably well. However, the MARS code overpredicts the flow rate. This is to be discussed later.
(a) Case 4-2: 1.6 bar
(b) Case 4-6: 1.9 bar
3.2. Heat Flux
The flow rate is proportional to the 1/3 power of the power for a single-phase natural circulation flow. This trend is similar in a two-phase case. Cases 4-6, 4-7, and 4-8 are carried out with the heat flux of 25%, 50%, and 75% of the rated power, respectively. The resulting natural circulation flow behaviors are compared in Figure 12. In the case of 25% power, the MARS code overpredicts the flow rate and, for the case of 50% and 75% power, the code predictions are better.
(a) Case 4-6: 25% power
(b) Case 4-7: 50% power
(c) Case 4-8: 75% power
In a low power, the driving force is relatively small and this leads to greater uncertainties in the calculated mass flow rates. It is shown in Figure 12 that the flow rate gradually increases in the experiments but discontinuously increases in the calculations. Also, after the electricity was shut off at each test, the flow rate decreases rapidly in the experiments and, however, tends to decrease more or less slowly in the calculations. This might be related to the environmental heat loss, which was not considered in the calculations. It seems that thermal inertia was greatly evaluated because there was no heat loss in the calculations.
3.3. Form Loss of the Orifice
To confirm the effect of the form loss of the orifice, the experiments of Cases 3-1 and 4-1 were carried out. However, because both the form loss and the pressure were different in the two experiments, it is not easy to separately check the influence of the form loss on the natural circulation flow. Figures 13(a) and 13(b) show the results of calculations for the two cases. The lower pressure in Case 4-1 leads to an increase in the mass flow and, meanwhile, the greater form loss results in a decrease. The effect of the greater form loss is dominant in this case. Sensitivity calculations of Case 4-1 for the form loss were performed additionally. Figure 13(c) shows that the flow rate decreases as the loss coefficient increases and its effect is greater in the two-phase flow period. This is consistent with analytical study. The results in Figure 13 also imply that the mass flow rate can be better predicted by increasing the friction and/or form loss.
(a) Case 3-1: = 1.90 (1.9 bar)
(b) Case 4-1: = 20.3 (1.6 bar)
(c) Calculation with different ’s
3.4. The Mass Flow Rate versus the Inlet Subcooling
In Figure 14, the transient calculation result for Case 4-1 is represented in terms of the mass flow rate versus the inlet subcooling. The greater the subcooling degree, the smaller the bubble generation and the smaller the natural circulation flow rate. However, when the subcooling is smaller than a certain value (e.g., smaller than 12 K in Figure 14), the flow rate tends to decrease because some bubbles from the exit of the cooling channel are introduced into the downcomer and this induces a smaller hydrostatic head and a greater two-phase frictional loss along the cooling channel. This feature is well captured by the MARS code. However, the overestimation of the mass flow rates is clearly shown again.
3.5. Discussions on the Problems of Calculation Results
It was shown that the MARS code can capture the key features of the two-phase natural circulation flow at the CE-PECS facility. However, some qualitative and quantitative problems were shown clearly. Among the various ones, three problems are discussed in the remainder of this subsection.
The natural circulation flow rate fluctuated severely but showed a continuous pattern in the experiments. However, in the calculations, the flow rate showed a discontinuous pattern. It is shown in Figure 15(a) that the calculated flow rate suddenly changes at ~2600 s, ~3050 s, ~3100 s, ~3200 s, and ~3700 s. When the flow rate suddenly increases, the void fraction also sharply rises as shown in Figure 15(b). These discontinuities coincided with the change in the two-phase flow regime in the vertical part of the cooling channel. That is, at each discontinuity, the two-phase flow regime changed from a bubbly flow to a slug flow. It seems that interfacial drag and wall friction suddenly changed at the moment of flow regime changes. This leads to discontinuity of the flow rate in the calculations, which needs improvement.
(a) Mass flow rate and flow regime
(b) Void fraction at the exit
Depending on the degree of flow mixing in the water tank, the temperature of the water entering the downcomer from the water tank is determined. As mentioned in Section 2.2.1, the flow inside the water tank is analyzed in three dimensions with a small flow mixing length, and good results could be obtained. However, as shown in Figure 16, the flow mixing inside the water tank is still overestimated by the MARS code. As a result, relatively cold water is introduced into the downcomer as compared to the experimental results. Flow mixing occurs due to the turbulence model and numerical diffusion in the MARS code. Excessive flow mixing appears to be due to the effect of the numerical diffusion that results from the coarse computational cells for the water tanks and the first-order upwind scheme. This is considered to be one of the limitations of a system-scale 3D analysis.
(a) Experiment
(b) MARS calculation
Generally, the MARS code tended to overestimate the natural circulation flow rate and underestimate the void fraction at the exit of cooling channel. This is a consistent result in terms of energy balance. An overestimation of the flow rate appears to be associated with underestimating the pressure drop along the cooling channel, especially in the upper vertical part of the cooling channel. For more detailed analysis, Case 4-1 was analyzed again using the CUPID code [13], which is a 3D two-phase flow analysis code [14, 15]. The results showed a recirculation flow in the vertical part of the cooling channel at a quasi-steady state. Figure 17(a) shows the 3D calculation results in the vertical cross section between the stud and the side wall at the cold channel exit. As can be seen in Figure 17(a), cold water enters into the exit of the cooling channel from the bottom of the water tank and the recirculation flow is generated. The length of the recirculation region, where the downward flow occurs in the cooling channel, depends on interfacial heat transfer and interfacial drag. This recirculation flow (i.e., a countercurrent flow), which is shown in Figure 17(b), did not appear in the 1D analysis using the MARS code. The experimental results showed that the void fraction at the exit of the cooling channel is significantly oscillating. It is noted that the gap size of the rectangular cooling channel is 0.1 m and one side (i.e., left wall) of the cooling channel is heated, resulting in vapor generation from the left side only. If the countercurrent flow is true, the total pressure drop along the cooling channel should increase. This leads to a decrease in the natural circulation flow rate. In this regard, it is no wonder that the MARS code, which does not predict the countercurrent flow in the cooling channel, overpredicts the natural circulation flow rate. Unfortunately, the flow velocities near the exit were not measured in the experiment, so it is not possible to confirm the fact but this seems to be possible. 3D measurement and a multiscale analysis [16] are necessary to resolve this issue. The lessons from a fine-scale calculation, thereafter, can be used for improving the coarse scale code.
(a) Velocity vector and void fraction near the exit of the cooling channel
(b) Velocity vector and temperature in the cooling channel
4. Summary and Conclusions
In this study, the system-scale thermal-hydraulic code, MARS-KS1.3, has been assessed using the two-phase natural circulation experiments at the CE-PECS facility, which had been conducted to verify the cooling performance of a core catcher. The facility was modeled by using a combination of 1D and 3D components: the inclined cooling channel and the downcomer were modeled by 1D pipe components and the cooling water tank was modeled by a 3D component. Because natural circulation flow in the loop depends on heat flux, pressure, subcooling, friction and form losses, and so forth, we selected six cases among the series of the CE-PECS experiments that can show parametric effects.
The results of assessment show that the MARS code can predict the two-phase natural circulation flow at the facility reasonably well by using the 3D component for the water tank. Also, the code captured parametric effects of heat flux, pressure, subcooling, and form loss very well. However, the following limitations of the MARS code were revealed:(i)Although the natural circulation flow rate fluctuated severely, it showed a continuous pattern in the experiments. However, in the calculations, the flow rate showed a discontinuous pattern. These discontinuities coincided with the flow regime changes in the vertical part of the cooling channel. It means that when the two-phase flow regime at the vertical cells changes from a bubbly flow to a slug flow, interfacial drag and wall friction also suddenly changed. This led to discontinuities of the flow rate in the calculations. This needs improvement.(ii)Depending on the degree of flow mixing in the water tank, the temperature of the water entering the downcomer is determined. The MARS code overestimated the flow mixing inside the water tank. As a result, relatively cold water is introduced into the downcomer as compared to the experimental results. Flow mixing occurs due to the turbulence model and numerical diffusion in the MARS code. Excessive flow mixing appears to be due to the numerical diffusion that results from the coarse computational cells and the first-order upwind scheme. This is considered to be one of the limitations of a system-scale 3D analysis.(iii)Generally the MARS code tended to overestimate the natural circulation flow rate and underestimate the void fraction at the exit of cooling channel. An overestimation of the flow appears to be related to underestimating the pressure drop along the one-side-heated cooling channel. The results of a refined 3D analysis using the CUPID code showed a recirculation flow in the vertical part of the cooling channel; cold water flows down at the right side of the vertical part of the cooling channel. If this is true, the total pressure drop should increase. This will lead to a decrease in the natural circulation flow rate. It is no wonder that the MARS code, which does not predict the recirculation at the vertical part of the cooling channel, overpredicts the natural circulation flow rate. In the experiments, the void fraction at the exit of the cooling channel is vigorously oscillated and, unfortunately, the flow velocities at the exit were not measured. So it is not possible to confirm the fact but this seems to be possible.
Despite these limitations, it can be said that the MARS code can predict the two-phase natural circulation experiments at the CE-PECS facility reasonably well and can be used for the design and safety analysis of advanced water-cooled reactors equipped with a core catcher. For a more accurate calculation, multiscale analysis combined with fine-scale 3D calculations is required, and discontinuities in the MARS results due to flow pattern changes must be improved. These conclusions are valid for the RELAP5/MOD3 code as well as the MARS code.
Nomenclature
| : | Specific heat capacity (J/kg·K) |
| : | Pipe diameter (m) |
| : | Gravitational acceleration (m/s2) |
| : | Form loss (—) |
| : | Length (m) |
| : | Mass flow rate (kg/s). |
| : | Thermal expansion coefficient (1/K) |
| : | Density (kg/m3) |
| : | Elevation difference between heat source and sink (m). |
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This research was supported by the Nuclear Safety Research Center Program of the Korea Foundation of Nuclear Safety (KOFONS) grant funded by Nuclear Safety and Security Commission of the Korean government (Grant code 1305011) and the National Research Foundation (NRF) grant funded by the Ministry of Science, ICT and Future Planning of the Korean government (2012M2A8A4025647).