- Email: [email protected]
Numerical simulation of oscillatory falling liquid film flows
Fusion Engineering and Design 65 (2003) 387 /392 www.elsevier.com/locate/fusengdes
Numerical simulation of oscillatory falling liquid film flows Tomoaki Kunugi *, Chiaki Kino Department of Nuclear Engineering, Kyoto University, Yoshida, Sakyo, Kyoto 606-8501, Japan
Abstract Two-dimensional numerical simulations of falling liquid film on a vertical wall with and without some artificial velocity oscillations have been performed in order to investigate the relationship between the interfacial wave behavior and heat transfer. Fully developed two-dimensional wave consists of a solitary and a capillary wave. Numerical results indicate that the heat transfer of falling film flow is enhanced by small vortices between the solitary and capillary waves in addition to the well-known flow recirculation in the solitary wave. Resulting from these findings, the heat transfer of the falling liquid film flow can be controlled by the artificial oscillation of the film flow. # 2003 Elsevier Science B.V. All rights reserved. Keywords: Liquid film; Heat transfer; Numerical simulation; Oscillatory falling liquid film flow
1. Introduction The liquid wall concept holds a central position in the advanced power extraction (APEX) study [1]. The APEX design idea includes both thin liquid films flowing with very high velocity over the first wall solid surface, and thick liquid layers acting as both the first wall and blanket flow. The working fluid is lithium-containing liquid metal or molten salt, such as Flibe (2LiF /BeF2). In general, a high Prandtl number (Pr /n/a, here n is viscosity and a is thermal diffusivity) fluid like Flibe (Pr / /30) has lower heat transport capability than that of liquid metal (Pr / /0.03 for
* Corresponding author. Tel./fax: /81-75-753-5823. E-mail address: [email protected] (T. Kunugi).
Li), because of very low thermal diffusivity and low thermal mixing due to thinner thermal boundary layer at the interfacial region. As for applications of a free surface cooling concept and/or a wall protection scheme to fusion reactors, it is necessary to understand the instabilities of the liquid film flow on the plasma facing components because the liquid droplets caused by the wave breaking and the vapor from the free surface increase the plasma density and could cause too much plasma contamination. The falling film flow is very unstable [6] and includes a wide range of wave numbers. However, applying a frequency like a specified sound wave, the free surface behavior can be controlled according to existing experiments [2]. In the present study, two-dimensional numerical simulations of falling film flow with and without
0920-3796/03/$ - see front matter # 2003 Elsevier Science B.V. All rights reserved. doi:10.1016/S0920-3796(03)00007-3
T. Kunugi, C. Kino / Fusion Engineering and Design 65 (2003) 387 /392
388
Fig. 1. Schematic computational model.
an artificial velocity oscillation have been carried out by means of the multi-interface advection and reconstruction solver (MARS) method [3]. Fig. 2. Film thickness behavior of developing film flow.
2. Numerical procedures
U(t; y) [1osin(2pft)]U0 (y)
The governing equations are the continuity for multi-phase flows, momentum equation based on the one-field model and the energy equation as follows:
here, U0 (y) is initial velocity profile, o amplitude of velocity fluctuation, t time and f is a specified frequency. In the present calculation, o is set to 0.03 as same as reference [4], f is 0, 13, 20 and 45 Hz for water and 20 Hz for Flibe. The film Reynolds number (Re /Umh0/n) is 75 for water and Re /30 and 60 for Flibe, here Um is mean velocity, h0 is inlet film thickness and n is kinetic viscosity. Especially, Re /75 is a critical Re [6]. Pressure gradient of outflow boundary is assumed to be zero, pressure of the gas side boundary surrounding the solution domain is assumed to be constant. At the wall, the non-slip velocity condition is applied. Temperature at wall
@Fm @t
9(Fm V ) Fm 9V 0
@V 1 9(VV ) G (9P FV ) 9t @t r
(1) (2)
@ hrCv iM T 9(hrCv iM TV ) @t 9(hliM 9T)Q
(3)
Here, F is the volume fraction of fluid, l is the thermal conductivity, Cv is the specific heat at constant volume and the suffix m denotes the m th fluid or phase, denotes an average of material properties and FV is body force due to the surface tension based on the continuum surface force model [5]. The interface tracking technique is based on the MARS [3]. The mesh sizes are (Dx , Dy ) /(0.3h0, 0.075h0). Here, h0 is an equilibrium liquid film thickness, for instance h0 /8.0 /10 4 m in case of Re /30. The schematic computational model is as shown in Fig. 1. The forcing function for the inlet velocity U (t, y ) is as follows:
(4)
Fig. 3. Comparison of wave peak height between simulation and experimental data [2].
T. Kunugi, C. Kino / Fusion Engineering and Design 65 (2003) 387 /392
389
Fig. 4. Film thickness distribution of the falling water film flow along the wall: (a) f /13 Hz, (b) f /20 Hz, (c) f /45 Hz.
is constant. Constant radiation heat flux at free surface is q /0.1 MW/m2. The wall is made of stainless steel.
3. Results and discussions
Fig. 5. (a) Streamlines observed from moving coordinates with wave celerity of water film flow for Re /75 and f /20 Hz. (b) Vector of velocity for Re /75 and f /20 Hz.
At first, the validation test of this numerical program are carried out for f /0 and 20 Hz with Re /75 for water as a reference fluid and the numerical results are compared with the experimental data by Nosoko [2]. Fig. 2(a and b) show
390
T. Kunugi, C. Kino / Fusion Engineering and Design 65 (2003) 387 /392
Fig. 6. Streamlines observed from moving coordinates with wave celerity of Flibe film flow for Re /30 and 60 and f /20 Hz.
Fig. 7. Temperature contours of Flibe film flow for Re /30 and 60.
the film thickness distribution for f /0 and 20 Hz, respectively. According to the stability theory, the vertical falling film flows are unconditionally unstable [6]. The film thickness shown in Fig. 2(a) is very stable. This discrepancy often appears in numerical simulation problems in a low parameter region because of some indirect restrictions such as a strong two-dimensionality of the flow, limitation of solution domain and number of grids, etc. However, as for the oscillatory case as shown in Fig. 2(b), the liquid film becomes unstable in the downstream region, and eventually the solitary wave and the capillary waves in front of the solitary one are developed. The numerical results of the relationship between the peak height
of solitary wave and the frequency show very good agreement with the experiments [2] as shown in Fig. 3. Fig. 4 shows the film thickness h distribution of the falling water film flow along the wall for f / 13, 20 and 45 Hz. The dimensionless film thickness and distance from the inlet are d /h/h0, and X / x/h0, respectively. In case of lower frequency, the film thickness becomes higher than that of higher frequency case. Fig. 5(a) shows the streamlines on the moving coordinate with the wave celerity of water film flow for Re /75 and f /20 Hz. A large recirculation flow can be seen in the solitary wave. Fig. 5(b) shows the velocity vector and the vortices can be
T. Kunugi, C. Kino / Fusion Engineering and Design 65 (2003) 387 /392
Fig. 8. Surface temperature variation along the flow direction.
seen between the waves. There is no report on these vortices in the previous experimental and numerical studies except our ongoing study [7]. Fig. 6 shows the streamlines observed from moving coordinates with wave celerity of Flibe film flow for Re /30 and 60 and f /20 Hz. The recirculation region appears at the upper part of solitary wave for Re /30, but it appears at the downstream front roll for Re /60. From these flow field results, the effect of streamline pattern on the heat transfer characteristics of Flibe film flow could be very large. Fig. 7 shows the temperature contours of Flibe film flow for Re /30 and 60. The temperature contours appears in parallel to the streamlines in case of Re /30. As for Re /60, a high temperature region appears at the downstream front of solitary wave and the heat flows across the streamlines there. Fig. 8 shows the local surface temperature variation along the flow direction for Re /30 and 60. The local surface temperature varies with
391
Fig. 9. Local heat transfer coefficent variation at wall surface for Flibe.
Fig. 10. Temporal variation of mean surface temperature.
the wave height, but the temperature level seems to be saturated soon for both cases. Fig. 9 shows the local heat transfer coefficient distribution along the wall for Re /30 and 60: hlocal /[(Twall /T )/(Twall /Tb)]l/Dy , here Twall is the wall temperature, Dy is the distance from wall and Tb is the bulk temperature, l is the thermal
392
T. Kunugi, C. Kino / Fusion Engineering and Design 65 (2003) 387 /392
into the falling film flow in order to enhance the heat transfer.
4. Summary
Fig. 11. Temporal variation of mean heat transfer coefficient.
conductivity of Flibe. The local heat transfer coefficient increases towards the downstream for both cases. The heat transfer for Re /30 is much higher than that of Re /60. This tendency is the same as the well-known Nusselt theory for heat transfer of falling liquid film. Fig. 10 shows the temporal variation of mean surface temperature of Flibe for Re /30 and 60. The mean surface temperature for Re /60 is relatively larger than that of Re /30. This means the thermal mixing would be enhanced in the solitary wave as shown in Fig. 7 (b). The trend of temporal mean surface temperature seems to be saturated. Since the local surface temperature reaches steady state quickly as shown in Fig. 8, we expect steady state temperature in the downstream region. Fig. 11 shows the temporal variation of mean heat transfer coefficient for Re /30 and 60. The heat transfer coefficient for Re /30 case is almost twice of that for Re /60. This means laminar film flow shows the limitation of heat transfer performance even as the Re of the flow increases. This tendency is very good agreement with the wellknown facts and the theory. Therefore, some turbulent flow characteristics must be introduced
Two-dimensional numerical simulations of the falling liquid film on a vertical wall with and without some artificial velocity oscillations have been performed in order to investigate the relationship between the interfacial wave behavior and heat transfer. Numerical simulation shows the feasibility of controlling the liquid film flow on the vertical plate stably by means of artificial velocity oscillation. It was found that the relation between wave structure and heat transfer characteristics is very important and some turbulent flow characteristics are necessary to reduce the surface temperature and enhance the heat transfer at the wall.
References [1] M.A. Abdou, et al., Exploring novel high power density concepts for attractive fusion systems, Fusion Eng. Des. 45 (1999) 145 /167. [2] T. Nosoko, Characteristics of two-dimensional waves on a falling liquid film, Chem. Eng. Soc. 51 (1996) 725 /732. [3] T. Kunugi, MARS for multiphase calculation, Comput. Fluid Dynam. J. 9 (2001) 563 /571. [4] A. Miyara, Numerical simulation of wavy liquid film flowing down on a vertical wall and an inclined wall, Int. J. Therm. Sci. 39 (2000) 1015 /1027. [5] J.U. Brackbill, D.B. Kothe, C. Zemach, A continuum method for modeling surface tension, J. Comput. Phys. 100 (1992) 335 /354. [6] S.V. Alekseenko, V.E. Nakoryakov, B.G. Pokusaev, Wave Flow of Liquid Films, Begell House, New York, 1994. [7] T. Kunugi, C. Kino, A. Serizawa, Surface wave structure and heat transfer of vertical liquid film flow to be presented at Thirty-ninth National Heat Transfer Conference, Sapporo, Japan, June 4 /7, 2002 and also German /Japanese Workshop on Multiphase Flow, Karlsruhe, August 25 /27, 2002.
Numerical simulation of oscillatory falling liquid film flows Tomoaki Kunugi *, Chiaki Kino Department of Nuclear Engineering, Kyoto University, Yoshida, Sakyo, Kyoto 606-8501, Japan
Abstract Two-dimensional numerical simulations of falling liquid film on a vertical wall with and without some artificial velocity oscillations have been performed in order to investigate the relationship between the interfacial wave behavior and heat transfer. Fully developed two-dimensional wave consists of a solitary and a capillary wave. Numerical results indicate that the heat transfer of falling film flow is enhanced by small vortices between the solitary and capillary waves in addition to the well-known flow recirculation in the solitary wave. Resulting from these findings, the heat transfer of the falling liquid film flow can be controlled by the artificial oscillation of the film flow. # 2003 Elsevier Science B.V. All rights reserved. Keywords: Liquid film; Heat transfer; Numerical simulation; Oscillatory falling liquid film flow
1. Introduction The liquid wall concept holds a central position in the advanced power extraction (APEX) study [1]. The APEX design idea includes both thin liquid films flowing with very high velocity over the first wall solid surface, and thick liquid layers acting as both the first wall and blanket flow. The working fluid is lithium-containing liquid metal or molten salt, such as Flibe (2LiF /BeF2). In general, a high Prandtl number (Pr /n/a, here n is viscosity and a is thermal diffusivity) fluid like Flibe (Pr / /30) has lower heat transport capability than that of liquid metal (Pr / /0.03 for
* Corresponding author. Tel./fax: /81-75-753-5823. E-mail address: [email protected] (T. Kunugi).
Li), because of very low thermal diffusivity and low thermal mixing due to thinner thermal boundary layer at the interfacial region. As for applications of a free surface cooling concept and/or a wall protection scheme to fusion reactors, it is necessary to understand the instabilities of the liquid film flow on the plasma facing components because the liquid droplets caused by the wave breaking and the vapor from the free surface increase the plasma density and could cause too much plasma contamination. The falling film flow is very unstable [6] and includes a wide range of wave numbers. However, applying a frequency like a specified sound wave, the free surface behavior can be controlled according to existing experiments [2]. In the present study, two-dimensional numerical simulations of falling film flow with and without
0920-3796/03/$ - see front matter # 2003 Elsevier Science B.V. All rights reserved. doi:10.1016/S0920-3796(03)00007-3
T. Kunugi, C. Kino / Fusion Engineering and Design 65 (2003) 387 /392
388
Fig. 1. Schematic computational model.
an artificial velocity oscillation have been carried out by means of the multi-interface advection and reconstruction solver (MARS) method [3]. Fig. 2. Film thickness behavior of developing film flow.
2. Numerical procedures
U(t; y) [1osin(2pft)]U0 (y)
The governing equations are the continuity for multi-phase flows, momentum equation based on the one-field model and the energy equation as follows:
here, U0 (y) is initial velocity profile, o amplitude of velocity fluctuation, t time and f is a specified frequency. In the present calculation, o is set to 0.03 as same as reference [4], f is 0, 13, 20 and 45 Hz for water and 20 Hz for Flibe. The film Reynolds number (Re /Umh0/n) is 75 for water and Re /30 and 60 for Flibe, here Um is mean velocity, h0 is inlet film thickness and n is kinetic viscosity. Especially, Re /75 is a critical Re [6]. Pressure gradient of outflow boundary is assumed to be zero, pressure of the gas side boundary surrounding the solution domain is assumed to be constant. At the wall, the non-slip velocity condition is applied. Temperature at wall
@Fm @t
9(Fm V ) Fm 9V 0
@V 1 9(VV ) G (9P FV ) 9t @t r
(1) (2)
@ hrCv iM T 9(hrCv iM TV ) @t 9(hliM 9T)Q
(3)
Here, F is the volume fraction of fluid, l is the thermal conductivity, Cv is the specific heat at constant volume and the suffix m denotes the m th fluid or phase, denotes an average of material properties and FV is body force due to the surface tension based on the continuum surface force model [5]. The interface tracking technique is based on the MARS [3]. The mesh sizes are (Dx , Dy ) /(0.3h0, 0.075h0). Here, h0 is an equilibrium liquid film thickness, for instance h0 /8.0 /10 4 m in case of Re /30. The schematic computational model is as shown in Fig. 1. The forcing function for the inlet velocity U (t, y ) is as follows:
(4)
Fig. 3. Comparison of wave peak height between simulation and experimental data [2].
T. Kunugi, C. Kino / Fusion Engineering and Design 65 (2003) 387 /392
389
Fig. 4. Film thickness distribution of the falling water film flow along the wall: (a) f /13 Hz, (b) f /20 Hz, (c) f /45 Hz.
is constant. Constant radiation heat flux at free surface is q /0.1 MW/m2. The wall is made of stainless steel.
3. Results and discussions
Fig. 5. (a) Streamlines observed from moving coordinates with wave celerity of water film flow for Re /75 and f /20 Hz. (b) Vector of velocity for Re /75 and f /20 Hz.
At first, the validation test of this numerical program are carried out for f /0 and 20 Hz with Re /75 for water as a reference fluid and the numerical results are compared with the experimental data by Nosoko [2]. Fig. 2(a and b) show
390
T. Kunugi, C. Kino / Fusion Engineering and Design 65 (2003) 387 /392
Fig. 6. Streamlines observed from moving coordinates with wave celerity of Flibe film flow for Re /30 and 60 and f /20 Hz.
Fig. 7. Temperature contours of Flibe film flow for Re /30 and 60.
the film thickness distribution for f /0 and 20 Hz, respectively. According to the stability theory, the vertical falling film flows are unconditionally unstable [6]. The film thickness shown in Fig. 2(a) is very stable. This discrepancy often appears in numerical simulation problems in a low parameter region because of some indirect restrictions such as a strong two-dimensionality of the flow, limitation of solution domain and number of grids, etc. However, as for the oscillatory case as shown in Fig. 2(b), the liquid film becomes unstable in the downstream region, and eventually the solitary wave and the capillary waves in front of the solitary one are developed. The numerical results of the relationship between the peak height
of solitary wave and the frequency show very good agreement with the experiments [2] as shown in Fig. 3. Fig. 4 shows the film thickness h distribution of the falling water film flow along the wall for f / 13, 20 and 45 Hz. The dimensionless film thickness and distance from the inlet are d /h/h0, and X / x/h0, respectively. In case of lower frequency, the film thickness becomes higher than that of higher frequency case. Fig. 5(a) shows the streamlines on the moving coordinate with the wave celerity of water film flow for Re /75 and f /20 Hz. A large recirculation flow can be seen in the solitary wave. Fig. 5(b) shows the velocity vector and the vortices can be
T. Kunugi, C. Kino / Fusion Engineering and Design 65 (2003) 387 /392
Fig. 8. Surface temperature variation along the flow direction.
seen between the waves. There is no report on these vortices in the previous experimental and numerical studies except our ongoing study [7]. Fig. 6 shows the streamlines observed from moving coordinates with wave celerity of Flibe film flow for Re /30 and 60 and f /20 Hz. The recirculation region appears at the upper part of solitary wave for Re /30, but it appears at the downstream front roll for Re /60. From these flow field results, the effect of streamline pattern on the heat transfer characteristics of Flibe film flow could be very large. Fig. 7 shows the temperature contours of Flibe film flow for Re /30 and 60. The temperature contours appears in parallel to the streamlines in case of Re /30. As for Re /60, a high temperature region appears at the downstream front of solitary wave and the heat flows across the streamlines there. Fig. 8 shows the local surface temperature variation along the flow direction for Re /30 and 60. The local surface temperature varies with
391
Fig. 9. Local heat transfer coefficent variation at wall surface for Flibe.
Fig. 10. Temporal variation of mean surface temperature.
the wave height, but the temperature level seems to be saturated soon for both cases. Fig. 9 shows the local heat transfer coefficient distribution along the wall for Re /30 and 60: hlocal /[(Twall /T )/(Twall /Tb)]l/Dy , here Twall is the wall temperature, Dy is the distance from wall and Tb is the bulk temperature, l is the thermal
392
T. Kunugi, C. Kino / Fusion Engineering and Design 65 (2003) 387 /392
into the falling film flow in order to enhance the heat transfer.
4. Summary
Fig. 11. Temporal variation of mean heat transfer coefficient.
conductivity of Flibe. The local heat transfer coefficient increases towards the downstream for both cases. The heat transfer for Re /30 is much higher than that of Re /60. This tendency is the same as the well-known Nusselt theory for heat transfer of falling liquid film. Fig. 10 shows the temporal variation of mean surface temperature of Flibe for Re /30 and 60. The mean surface temperature for Re /60 is relatively larger than that of Re /30. This means the thermal mixing would be enhanced in the solitary wave as shown in Fig. 7 (b). The trend of temporal mean surface temperature seems to be saturated. Since the local surface temperature reaches steady state quickly as shown in Fig. 8, we expect steady state temperature in the downstream region. Fig. 11 shows the temporal variation of mean heat transfer coefficient for Re /30 and 60. The heat transfer coefficient for Re /30 case is almost twice of that for Re /60. This means laminar film flow shows the limitation of heat transfer performance even as the Re of the flow increases. This tendency is very good agreement with the wellknown facts and the theory. Therefore, some turbulent flow characteristics must be introduced
Two-dimensional numerical simulations of the falling liquid film on a vertical wall with and without some artificial velocity oscillations have been performed in order to investigate the relationship between the interfacial wave behavior and heat transfer. Numerical simulation shows the feasibility of controlling the liquid film flow on the vertical plate stably by means of artificial velocity oscillation. It was found that the relation between wave structure and heat transfer characteristics is very important and some turbulent flow characteristics are necessary to reduce the surface temperature and enhance the heat transfer at the wall.
References [1] M.A. Abdou, et al., Exploring novel high power density concepts for attractive fusion systems, Fusion Eng. Des. 45 (1999) 145 /167. [2] T. Nosoko, Characteristics of two-dimensional waves on a falling liquid film, Chem. Eng. Soc. 51 (1996) 725 /732. [3] T. Kunugi, MARS for multiphase calculation, Comput. Fluid Dynam. J. 9 (2001) 563 /571. [4] A. Miyara, Numerical simulation of wavy liquid film flowing down on a vertical wall and an inclined wall, Int. J. Therm. Sci. 39 (2000) 1015 /1027. [5] J.U. Brackbill, D.B. Kothe, C. Zemach, A continuum method for modeling surface tension, J. Comput. Phys. 100 (1992) 335 /354. [6] S.V. Alekseenko, V.E. Nakoryakov, B.G. Pokusaev, Wave Flow of Liquid Films, Begell House, New York, 1994. [7] T. Kunugi, C. Kino, A. Serizawa, Surface wave structure and heat transfer of vertical liquid film flow to be presented at Thirty-ninth National Heat Transfer Conference, Sapporo, Japan, June 4 /7, 2002 and also German /Japanese Workshop on Multiphase Flow, Karlsruhe, August 25 /27, 2002.