- Email: [email protected]
Numerical study of bubble generation and transport in a serpentine channel with a T-junction
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 9 ( 2 0 1 4 ) 2 3 2 5 e2 3 3 3
Available online at www.sciencedirect.com
ScienceDirect journal homepage: www.elsevier.com/locate/he
Numerical study of bubble generation and transport in a serpentine channel with a T-junction Simo Kang, Biao Zhou* Department of Mechanical, Automotive and Materials Engineering, University of Windsor, Windsor, ON N9B 3P4, Canada
article info
abstract
Article history:
The two-phase flow problem is one of the most critical challenges for direct methanol fuel
Received 24 August 2013
cells (DMFC). With carbon dioxide generated in the anode catalyst layer, a bubble that could
Received in revised form
not be removed from the flow channel would hinder the oxidation reaction. In this paper,
21 November 2013
bubble generation and flow in a micro-serpentine channel is modeled using a multi-phase
Accepted 27 November 2013
three-dimensional NaviereStokes plus volume of fluid (VOF) method. The bubble genera-
Available online 30 December 2013
tion process, which results from the gaseliquid flow and the effects of surface tension and viscosity, are discussed and validated by a comparison with a serious of experimental
Keywords:
results. The results prove that the VOF method is an effective method to simulate flow
Direct methanol fuel cell
fields where the bubble flow phenomena exist, such as the two-phase flow fields in DMFCs.
Serpentine
Copyright ª 2013, Hydrogen Energy Publications, LLC. Published by Elsevier Ltd. All rights reserved.
Bubble Volume of fluid Viscosity Surface tension
1.
Introduction
Bubble generation and transport are always crucial subjects in the study of micro channel devices. This research is relevant to the gas e liquid two-phase flow in direct methanol fuel cells (DMFCs). A DMFC is a non-hydrogen device that can utilize methanol directly as the fuel. In the anode side of a DMFC, gas CO2 is generated by electroechemical reaction and transferred out of the cell by the methanol flow inside the anode channel. This removal process should be as efficient as possible; otherwise, the bubble of CO2 choked in the anode channel would increase the gradient of methanol concentration along the channel and hinder the continuous reaction. With the development of DMFCs, a number of studies have
verified that the bubble behaviors in the micro channel have an enormous influence on DMFC performance [1,2]. Gas management has been one of the primary issues in DMFC research. To explore the physical phenomena during bubble removal, many experiments of bubble behavior in micro channels have been conducted. Xiong et al. [3] analyzed bubble generation and transport in a microfluidic device with a high aspect ratio. The effects of liquid viscosity and surface tension were discussed accordingly. Fu et al. [4] explored the nucleation and growth of CO2 bubbles generated by chemical reactions in three different types of micro-channels. In addition to the experimental investigations of bubble or droplet behaviors, numerous numerical methods, such as front-tracking method [5e7], Lattice Boltzmann method [8],
* Corresponding author. Tel.: þ1 519 253 3000x2630. E-mail address: [email protected] (B. Zhou). 0360-3199/$ e see front matter Copyright ª 2013, Hydrogen Energy Publications, LLC. Published by Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.ijhydene.2013.11.115
2326
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 9 ( 2 0 1 4 ) 2 3 2 5 e2 3 3 3
Fig. 1 e Computational domain of the serpentine channel.
and volume of fluid method [9], have been introduced to explore the bubbles/droplet behaviors in the simple channels (either straight or U-shape) in order to obtain a reliable and fundamental understanding. In recent years, many researchers have introduced different CFD models to represent bubble behaviors in DMFC anode. For example, Kulikovsky [10,11] developed a simple numerical model based on the continuity equations that explored the relationship between gaseous bubble behavior and local methanol concentration in the anode channel. This model produced an asymptotic solution. Wang et al. [12] introduced the mixture model to study DMFC, which includes fully coupled heat transfer, methanol crossover and addition components. Liu et al. [13] advanced DMFC modeling research by developing a three-dimensioned flow model considering the uneven distribution of the reaction rate in the catalyst layer. However, in all these models, the gaseliquid interface is not tracked, because the fluid was considered to be a mixture or single phase within the flow channel. Understanding the detailed bubble removal process with an accurate simulation method that can track the gaseliquid interface is a must for improving DMFC performance.
According to the author’s previous study, employing a volume of fluid (VOF) method is advantageous to track the interface between the gas and liquid phases [14]. As stated in previous works, the VOF method has been validated for the simulation of two-phase flow in a polymer electrolyte membrane fuel cell [15]. It is possible to use the VOF methodology to solve the similar problem of two-phase flow in DMFCs. In DMFC, CO2 bubbles interact with methanol. But it’s very difficult to find a well documented DMFC experiments from the available literature in terms of CO2 and methanol interaction. As a necessary step to verify the possibility of VOF method for DMFC modeling, in this paper, a numerical simulation is conducted about water and air bubble interaction in a serpentine channel and compared with well-documented experimental results [3]. The reasons for choosing these experimental results to compare are as follows: Two-phase interface with bubble phenomena are correctly tracked with experimental method; Most necessary parameters are well documented; This experiment is artfully designed with bubble breakup occurring at the T-junction process;
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 9 ( 2 0 1 4 ) 2 3 2 5 e2 3 3 3
2327
The serpentine geometry is one type of fundamental channel design used in DMFCs; Different surface tension and viscosity have been experimentally studied and could help to validate the model. In the following sections, the numerical simulation domain, mathematical model, and initial settings are introduced. More discussion about bubble generation in T-junctions and the removal process, corner effects, the pressure and velocity fields, and the effects of surface tension and viscosity are presented.
2.
Numerical model setup
2.1. Numerical simulation domain and boundary conditions The schematic numerical simulation domain represents a micro serpentine channel with a T-junction, as shown in Figs. 1e5. The cross-section of 594 (w) 80 (h) mm2 and the high aspect ratio of 7.425 are the same in the experiment and the simulation. Because the channel length was not documented in the experiment paper, it was specified by analyzing the figures plotted in the experimental publication [3]. The gas and liquid inlets are from different channels, which are perpendicular to each other. In this simulation, 2 mm length channels for the pure gas or liquid channels are specified. Four inverse “U” channels are included in the serpentine part with a height of 2.75 mm and width of 2.6 mm. The whole dimension is 22.06 mm in the X-direction and 4.75 mm in the Y-direction. The thickness of the whole domain is 0.594 mm. A no-slip boundary condition is applied to the walls. The direction of gravity is along the negative Z-direction in the operation condition. The gas inlet is set with a flow rate of 3.5 109 kg1, which is normal to the inlet boundary. The liquid inlet has a constant flow rate of 4.31654 106 kg1, which is from a vertical direction. The pressure outlet is applied at the outlet boundary. The contact angle is specified as 35 in the side walls of channels with a hydrophilic assumption. The computational domain consisted of 20, 910 cells with a minimum grid size of 1.944 1015 m3.
2.2.
Fig. 2 e Validation of bubble generation in the T-junction.
Computational methodology
In this paper, FLUENT was used with a 3D unsteady VOF model to perform the numerical simulations and to track the gaseliquid interface inside of the computational domain. A user defined function (UDF) was applied to track the pressure and gas phase amount with time. The flow in the gas phase was assumed to be laminar. The mass conservation law that governs unsteady, laminar flow can be written as vðrÞ þ V$ðruÞ ¼ Sm vt
(1)
where sg expresses the volume fraction of the gaseous phase. The tracking of the interface between the phases is accomplished by the solution of the continuity equation for the volume fraction of gaseous phase as the following equation [16]: v sg rg vt
þ V$ sg rg u ¼ 0
(3)
The source term is specified as 0 for this simplified model. The momentum equation is shown below:
The density in the governing equation was the density of the mixture phase, which could be calculated as
v ! ðr v Þ þ Vðr! v! v Þ ¼ Vp þ V½mV! v þ Sv vt
r ¼ 1 sg rl þ sg rg
where m represents the dynamic viscosity with the units 1 kgm s1 and is calculated as
(2)
(4)
2328
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 9 ( 2 0 1 4 ) 2 3 2 5 e2 3 3 3
Fig. 3 e Gas volume fraction and velocity vector in the T-junction.
m ¼ 1 sg ml þ sg mg
b ¼ n=jnj and n ¼ Vð1 sg Þ is the surface normal that k ¼ V$ n is defined as the gradient of ð1 sg Þ.
In this paper, the surface tension was represented by using a continuum surface force (CSF) model. The source term in the anode channel is expressed in terms of a pressure rise across the interface and can be written as follows:
3. Validation with the break-up process at the T-junction
rkV 1 sg . Su ¼ slg rl þ rg 2
(5)
where s is the surface tension coefficient, r represents the density, and k represents the surface curvature, which is b: defined in terms of the divergence of the unit normal n
3.1. Comparison of the numerical simulation and experimental visualization with the break-up process at the T-junction Fig. 2 shows the entire process of the first bubble generation at the T-junction from both the experimental results (left) [3] and
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 9 ( 2 0 1 4 ) 2 3 2 5 e2 3 3 3
2329
Fig. 4 e Volume fraction of the gas and pressure distributions in the serpentine channel. (The figures in the left column show the volume fraction of gas, and the figures on the right show the pressure distribution at (a) t [ 0.0148 s, (b) t [ 0.1480 s, and (c) t [ 0.1840 s).
the simulation results (right). As discussed in the experimental study, the bubble generation (break-up) procedure is divided into three phases: gas ligament expansion, gas ligament collapse, and bubble removal phase. Fig 2 shows that this procedure is similar for the simulation result. Moreover, the modeling results are helpful for analysis of the velocity field and interface behaviors as follows. As shown in Fig 2(aeg), after gas enters from the left inlet with a constant flow rate, a bubble slowly progresses through the channel pushed by the gas supplied from the gas inlet and opposed by the resistance of the liquid already present in the channel. The head of the bubble retains a convex round shape due to the contact angle with the side wall. As the gas slug arrives, the T-junction of the two inlet channels liquid enters from the liquid inlet channel, pushing the bubble to the top side of the horizontal channel. The gas forms a neck above the perpendicular channel, where bearing from the shear force is applied. This sub-process is called gas ligament expansion. The gas ligament blocks the opening of the liquid channel, and results in an increased shear force. Thus, the neck, the only transport ligament of the gas, is reduced and finally severed by the liquid. This process is so called gas ligament collapse. Fig. 2(hel) shows the exiting bubble travel through the first inverted U shape of the serpentine channel. In Fig. 2(h) and Fig. 2(k), the bubble appears as an “S” shape attached to the top side of the lower channel and to the bottom side of the
upper channel. It appears to be attempting to take the shortest route to the outlet. This same reason leads to the bubble shapes in Fig. 2(iek) and Fig. 2(m); the bubble always attaches to the inner side of the converted U channel.
3.2.
Velocity field of the break-up process
Fig 3 shows the gaseous bubble behavior during the break-up process with the velocity field. Fig 3(a) shows that the liquid flows smoothly prior to the arrival of the bubble at the Tjunction. Due to the corner effects, the velocity vector is intense around the inner corner and unconsolidated near the external corner. The velocity field during the liquid phase near the T-junction is enlarged when the gas slug attempts to cover the exit of the liquid inlet channel, as shown in Fig 3(b). With the liquid inlet channel blocked, the liquid accumulates at the exit of the liquid inlet channel and pushes the air to the top wall and forms the “neck” shape (as presented in Section 3.1). As shown in Fig 3(f), when the gas bubbles approach, but do not touch, the wall of the channel, stress on the water will strongly increase in the gap. After the ligament has collapsed, strong liquid velocity forms a vortex between the bubble released and the gas slug in the gas inlet channel. The interaction between the airflow and the liquid water can provide strong shear stress and can fragment some of the air on the gasewater interface into small bubbles. These small bubbles
2330
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 9 ( 2 0 1 4 ) 2 3 2 5 e2 3 3 3
Fig. 5 e Volume fraction of the gas and pressure drops with the vector of the first two inverted U-shaped channels.
are rolled up by the vortex and consumed by the gasewater interface again.
3.3.
Pressure fields
Figs. 4 and 5 show the volume fraction of the gas and pressure distribution in the serpentine channel. The pressure drops gradually along the channel from inlet to outlet. Because the outlet pressure is assumed to remain relatively close to a zero pressure (atmosphere), this is used as the outlet boundary condition. The pressure around the inlet changes from 7000 Pa to 5500 Pa, as shown in the three figures of different times (t ¼ 0.0148 s, t ¼ 0.1480 s and t ¼ 0.1840 s). Fig 5 shows the pressure and velocity field of the first and second inverted U-shaped channel in the red blocks that are shown in Fig 4(b). The velocity of the gas phase is larger in the middle of the bubble. The bubble, which is located at the corner, tends to flow straight and result in the large velocity towards to the wall. Although the pressure tends to decrease along the channel, the pressure inside the gas bubbles is higher than inside the liquid and results in an interesting phenomenon, a velocity vector orientated from the lower pressure field to the high pressure field.
This phenomenon shows the effects from the capillary pressure difference across the interface between the liquid and gas with the surface tension. Because the contact angle is specified as 145 in the gas phase, the gas slug in the tube shows convex surfaces on both the head and tail sides, which indicates a higher pressure inside the gas is balanced by the surface tension.
3.4.
Pressure drop
Fig. 6 shows the pressure at the gas inlet, the liquid inlet and the channel outlet. Because the outlet is specified as the outlet pressure (atmosphere), the outlet pressure remains at 0, as shown in Fig. 6(a). Fig. 6(b) shows only the pressure at the gas inlet and the liquid inlet. As presented above, the pressure generally drops with more gas enters into the serpentine channel with strong fluctuations. With the bubble that is generated in the T-junction component, the pressure fluctuation established a rhythm. The period could be defined from a bubble breakup from the Tjunction to the next bubble generated. To see this periodical phenomenon, Fig. 6(c) is shown after a Stineman Interpolation curve was fit from Fig. 6(b).
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 9 ( 2 0 1 4 ) 2 3 2 5 e2 3 3 3
2331
T-junction section. Simultaneously, the gas slug begins to cover the liquid inlet channel. This results in the peak pressure during one period.
4. Effects of surface tension and liquid viscosity 4.1.
Effects of surface tension
The two images in Fig 7 show the results of two simulation experiments. The same boundary conditions (liquid/gas inlet velocity, contact angle, etc.) are adopted for these two cases, with the exception of surface tension coefficient. The liquid velocity is 0.273 m s1, and the gas velocity is 0.0546 m s1. The surface tension is selected as the experimental measured value of 0.072 N m1 and 0.037 N m1, separately. Also, the liquid density is slightly different between the two cases (998 kg m3 and 1000 kg m3 for the first and second cases, respectively) to maintain identical conditions with the experiments. Fig 7 shows the simulation results at the exact same time instant. With the same flow rate for the two types of fluid, the lower surface tension results in fewer bubbles. After 89.5 ms, the first two small bubbles combined together to form a tadpole-like shape as presented in the experimental paper. Because the results for the experiment are not detailed enough, the analysis below only represents the modeling results. As the surface tension could be defined as the property of the surface of a liquid that allows it to resist an external force [17], the asunder bubbles are reasonable. For the bubble generated in the T-junction, as presented in previous section, the shear force from the liquid inlet channel would be easier to cut the gas slug that has lower surface tension.
4.2.
Fig. 6 e Pressure at the gas inlet, the liquid inlet and the outlet.
At the time instant of the wave trough, approximately 0.064 s, 0.098 s, and 0.13 s, the bubble breakup at the T-junction resulted in a rapid pressure drop. At the times of 0.083 s, 0.12 s, and 0.155 s, the bubble generated previously, driven by the liquid flow, must leave the
Effects of viscosity
As shown in Table 1, the three cases are simulated to explore the effects of liquid viscosity and compared to experimental results [3]. The velocities of the gas or the liquid in the three cases are identified separately as 0.182 m s1 for the gas and 0.091 m s1 for the liquid. The liquid viscosity was increased from 0.00119 Pa s for the first case to 0.00731 Pa s for the third case. The liquid viscosity of 0.00268 Pa s was chosen for the second case. The simulation results show that the bubble generated in the case with higher liquid viscosity is shorter, as shown in the figures of the final column of Table 1, which is qualitatively similar to the schematic results obtained in the experimental study as shown in Table 1. At a fixed velocity ratio of gas and liquid phase, if the viscosity increases, then the bubble length and diameter will decrease, i.e., the bubble can change the shape from a large slug to a small one. This phenomenon has also been observed in the experimental schematic drawing from Fig. 6 in Ref. [3]. And the reason for this phenomenon has been explained in Ref. [3]: the bigger viscosity will produce a larger shear stress that causes the reduced bubble dimension.
2332
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 9 ( 2 0 1 4 ) 2 3 2 5 e2 3 3 3
Fig. 7 e Gaseous volume fraction of the two cases with different surface tensions.
5.
Conclusions
In this paper, a serpentine channel with a liquid inlet and a gas inlet was adopted as the geometry used for the VOF model application to investigate bubbles behaviors by analyzing the gaseliquid interaction. (1) The entire process of bubble generation and breakup in the T-junction resulting from the simulation is almost exactly the same as the experiment investigation. (2) By analyzing the velocity field and pressure distribution from the simulation results, some phenomena, such as
the corner effect and the surface behavior of bubbles, are explained. (3) The effect of the surface tension is studied under similar conditions with the exception of the surface tension. The results show a tadpole-like shape as in the experimental study. (4) The effect of the liquid viscosity is explored by three simulations. With the increase in the liquid viscosity, the bubble length along the channel is decreased due to the stronger shear force from the liquid. (5) The VOF algorithm could be applied to bubble effects modeling in DMFC.
Table 1 e Effect of liquid viscosity. Velocity of gas inlet (ms1)
Velocity of liquid inlet (ms1)
Liquid viscosity (Pa s)
0.182
0.091
0.00119
0.182
0.091
0.00268
0.182
0.091
0.00731
Bubble shape presented in experimental results [3]
Simulation results
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 9 ( 2 0 1 4 ) 2 3 2 5 e2 3 3 3
Acknowledgments The authors are grateful for the support of this work by the Auto21 Networks of the Centres of Excellence (Grant D303DFC), the Natural Sciences and Engineering Research Council of Canada (NSERC), the Canada Foundation for Innovation (CFI), the Ontario Innovation Trust (OIT), and the University of Windsor.
Nomenclature r u S s m s k K n R1, R2
Density (kg m3) Velocity vector (m s1) Source term Phase volume fraction Dynamic viscosity (Pa s) Surface tension (N m1) Surface curvature Permeability (m2) Normal vector Two radii in the orthogonal directions
Subscripts l Liquid phase g Gas phase
references
[1] Lu GQ, Wang CY. Electrochemical and flow characterization of a direct methanol fuel cell. J Power Sources 2004;134:33e40. [2] Wong CW, Zhao TS, Ye Q, Liu JG. Transient capillary blocking in the flow field of a micro-DMFC and its effect on cell performance. J Electrochem Soc 2005;152(8):A1600e5.
2333
[3] Xiong R, Chung JN. Bubble generation and transport in a microfluidic device with high aspect ratio. Exp Therm Fluid Sci 2009;33:1156e62. [4] Fu BR, Pan C. Bubble growth with chemical reactions in microchannels. Int J Heat Mass Transs 2009;52:767e76. [5] Quan S, Lou J. Viscosity-ratio-based scaling for the rise velocity of a Taylor drop in a vertical tube. Phys Rev E 2011;84(3):036320. [6] Quan S. Co-current flow effects on a rising Taylor bubble. Int J Multiph Flow 2011;37:888e97. [7] Kang C, Quan S, Lou J. Numerical study of a Taylor bubble rising in stagnant liquids. Phys Rev E 2010;81(6):066308. [8] Fei K, Chen TS, Hong CW. Direct methanol fuel cell bubble transport simulations via thermal lattice. J Power Sources 2010;195:1940e5. [9] Quan Peng, Zhou Biao, Sobiesiak Andrzej, Liu Zhongsheng. Water behavior in serpentine micro-channel for proton exchange membrane fuel cell cathode. J Power Sources 2005;152:131e45. [10] Kulikovsky AA. Model of the flow with bubbles in the anode channel and performance of a direct methanol fuel cell. Electrochem Commun 2005;7:237e43. [11] Kulikovsky AA. Bubbles in the anode channel and performance of a DMFC: asymptotic solutions. Electrochim Acta 2006;51:2003e11. [12] Wang ZH, Wang CY. Mathematical modeling of liquid-feed direct methanol fuel cells. J. Electrochem Soc 2003;150(4):A508e19. [13] Ge J, Liu HT. A three-dimensional two-phase flow model for a liquid-fed direct methanol fuel cell. J Power Sources 2007;163:907e15. [14] Le Anh D, Zhou Biao. A general model of proton exchange membrane fuel cell. J Power Sources 2008;182:197e222. [15] Le Anh Dinh, Zhou Biao, Shiu Huan-Ruei, Lee Chun-I, Chang Wen-Chen. Numerical simulation and experimental validation of liquid water behaviors in a proton exchange membrane fuel cell cathode with serpentine channels. J Power Sources 2010;195:7302e15. [16] Fluent 6.3 manual. Fluent Inc; 2006. [17] Wikipedia, http://en.wikipedia.org/wiki/Surface_tension.
Available online at www.sciencedirect.com
ScienceDirect journal homepage: www.elsevier.com/locate/he
Numerical study of bubble generation and transport in a serpentine channel with a T-junction Simo Kang, Biao Zhou* Department of Mechanical, Automotive and Materials Engineering, University of Windsor, Windsor, ON N9B 3P4, Canada
article info
abstract
Article history:
The two-phase flow problem is one of the most critical challenges for direct methanol fuel
Received 24 August 2013
cells (DMFC). With carbon dioxide generated in the anode catalyst layer, a bubble that could
Received in revised form
not be removed from the flow channel would hinder the oxidation reaction. In this paper,
21 November 2013
bubble generation and flow in a micro-serpentine channel is modeled using a multi-phase
Accepted 27 November 2013
three-dimensional NaviereStokes plus volume of fluid (VOF) method. The bubble genera-
Available online 30 December 2013
tion process, which results from the gaseliquid flow and the effects of surface tension and viscosity, are discussed and validated by a comparison with a serious of experimental
Keywords:
results. The results prove that the VOF method is an effective method to simulate flow
Direct methanol fuel cell
fields where the bubble flow phenomena exist, such as the two-phase flow fields in DMFCs.
Serpentine
Copyright ª 2013, Hydrogen Energy Publications, LLC. Published by Elsevier Ltd. All rights reserved.
Bubble Volume of fluid Viscosity Surface tension
1.
Introduction
Bubble generation and transport are always crucial subjects in the study of micro channel devices. This research is relevant to the gas e liquid two-phase flow in direct methanol fuel cells (DMFCs). A DMFC is a non-hydrogen device that can utilize methanol directly as the fuel. In the anode side of a DMFC, gas CO2 is generated by electroechemical reaction and transferred out of the cell by the methanol flow inside the anode channel. This removal process should be as efficient as possible; otherwise, the bubble of CO2 choked in the anode channel would increase the gradient of methanol concentration along the channel and hinder the continuous reaction. With the development of DMFCs, a number of studies have
verified that the bubble behaviors in the micro channel have an enormous influence on DMFC performance [1,2]. Gas management has been one of the primary issues in DMFC research. To explore the physical phenomena during bubble removal, many experiments of bubble behavior in micro channels have been conducted. Xiong et al. [3] analyzed bubble generation and transport in a microfluidic device with a high aspect ratio. The effects of liquid viscosity and surface tension were discussed accordingly. Fu et al. [4] explored the nucleation and growth of CO2 bubbles generated by chemical reactions in three different types of micro-channels. In addition to the experimental investigations of bubble or droplet behaviors, numerous numerical methods, such as front-tracking method [5e7], Lattice Boltzmann method [8],
* Corresponding author. Tel.: þ1 519 253 3000x2630. E-mail address: [email protected] (B. Zhou). 0360-3199/$ e see front matter Copyright ª 2013, Hydrogen Energy Publications, LLC. Published by Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.ijhydene.2013.11.115
2326
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 9 ( 2 0 1 4 ) 2 3 2 5 e2 3 3 3
Fig. 1 e Computational domain of the serpentine channel.
and volume of fluid method [9], have been introduced to explore the bubbles/droplet behaviors in the simple channels (either straight or U-shape) in order to obtain a reliable and fundamental understanding. In recent years, many researchers have introduced different CFD models to represent bubble behaviors in DMFC anode. For example, Kulikovsky [10,11] developed a simple numerical model based on the continuity equations that explored the relationship between gaseous bubble behavior and local methanol concentration in the anode channel. This model produced an asymptotic solution. Wang et al. [12] introduced the mixture model to study DMFC, which includes fully coupled heat transfer, methanol crossover and addition components. Liu et al. [13] advanced DMFC modeling research by developing a three-dimensioned flow model considering the uneven distribution of the reaction rate in the catalyst layer. However, in all these models, the gaseliquid interface is not tracked, because the fluid was considered to be a mixture or single phase within the flow channel. Understanding the detailed bubble removal process with an accurate simulation method that can track the gaseliquid interface is a must for improving DMFC performance.
According to the author’s previous study, employing a volume of fluid (VOF) method is advantageous to track the interface between the gas and liquid phases [14]. As stated in previous works, the VOF method has been validated for the simulation of two-phase flow in a polymer electrolyte membrane fuel cell [15]. It is possible to use the VOF methodology to solve the similar problem of two-phase flow in DMFCs. In DMFC, CO2 bubbles interact with methanol. But it’s very difficult to find a well documented DMFC experiments from the available literature in terms of CO2 and methanol interaction. As a necessary step to verify the possibility of VOF method for DMFC modeling, in this paper, a numerical simulation is conducted about water and air bubble interaction in a serpentine channel and compared with well-documented experimental results [3]. The reasons for choosing these experimental results to compare are as follows: Two-phase interface with bubble phenomena are correctly tracked with experimental method; Most necessary parameters are well documented; This experiment is artfully designed with bubble breakup occurring at the T-junction process;
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 9 ( 2 0 1 4 ) 2 3 2 5 e2 3 3 3
2327
The serpentine geometry is one type of fundamental channel design used in DMFCs; Different surface tension and viscosity have been experimentally studied and could help to validate the model. In the following sections, the numerical simulation domain, mathematical model, and initial settings are introduced. More discussion about bubble generation in T-junctions and the removal process, corner effects, the pressure and velocity fields, and the effects of surface tension and viscosity are presented.
2.
Numerical model setup
2.1. Numerical simulation domain and boundary conditions The schematic numerical simulation domain represents a micro serpentine channel with a T-junction, as shown in Figs. 1e5. The cross-section of 594 (w) 80 (h) mm2 and the high aspect ratio of 7.425 are the same in the experiment and the simulation. Because the channel length was not documented in the experiment paper, it was specified by analyzing the figures plotted in the experimental publication [3]. The gas and liquid inlets are from different channels, which are perpendicular to each other. In this simulation, 2 mm length channels for the pure gas or liquid channels are specified. Four inverse “U” channels are included in the serpentine part with a height of 2.75 mm and width of 2.6 mm. The whole dimension is 22.06 mm in the X-direction and 4.75 mm in the Y-direction. The thickness of the whole domain is 0.594 mm. A no-slip boundary condition is applied to the walls. The direction of gravity is along the negative Z-direction in the operation condition. The gas inlet is set with a flow rate of 3.5 109 kg1, which is normal to the inlet boundary. The liquid inlet has a constant flow rate of 4.31654 106 kg1, which is from a vertical direction. The pressure outlet is applied at the outlet boundary. The contact angle is specified as 35 in the side walls of channels with a hydrophilic assumption. The computational domain consisted of 20, 910 cells with a minimum grid size of 1.944 1015 m3.
2.2.
Fig. 2 e Validation of bubble generation in the T-junction.
Computational methodology
In this paper, FLUENT was used with a 3D unsteady VOF model to perform the numerical simulations and to track the gaseliquid interface inside of the computational domain. A user defined function (UDF) was applied to track the pressure and gas phase amount with time. The flow in the gas phase was assumed to be laminar. The mass conservation law that governs unsteady, laminar flow can be written as vðrÞ þ V$ðruÞ ¼ Sm vt
(1)
where sg expresses the volume fraction of the gaseous phase. The tracking of the interface between the phases is accomplished by the solution of the continuity equation for the volume fraction of gaseous phase as the following equation [16]: v sg rg vt
þ V$ sg rg u ¼ 0
(3)
The source term is specified as 0 for this simplified model. The momentum equation is shown below:
The density in the governing equation was the density of the mixture phase, which could be calculated as
v ! ðr v Þ þ Vðr! v! v Þ ¼ Vp þ V½mV! v þ Sv vt
r ¼ 1 sg rl þ sg rg
where m represents the dynamic viscosity with the units 1 kgm s1 and is calculated as
(2)
(4)
2328
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 9 ( 2 0 1 4 ) 2 3 2 5 e2 3 3 3
Fig. 3 e Gas volume fraction and velocity vector in the T-junction.
m ¼ 1 sg ml þ sg mg
b ¼ n=jnj and n ¼ Vð1 sg Þ is the surface normal that k ¼ V$ n is defined as the gradient of ð1 sg Þ.
In this paper, the surface tension was represented by using a continuum surface force (CSF) model. The source term in the anode channel is expressed in terms of a pressure rise across the interface and can be written as follows:
3. Validation with the break-up process at the T-junction
rkV 1 sg . Su ¼ slg rl þ rg 2
(5)
where s is the surface tension coefficient, r represents the density, and k represents the surface curvature, which is b: defined in terms of the divergence of the unit normal n
3.1. Comparison of the numerical simulation and experimental visualization with the break-up process at the T-junction Fig. 2 shows the entire process of the first bubble generation at the T-junction from both the experimental results (left) [3] and
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 9 ( 2 0 1 4 ) 2 3 2 5 e2 3 3 3
2329
Fig. 4 e Volume fraction of the gas and pressure distributions in the serpentine channel. (The figures in the left column show the volume fraction of gas, and the figures on the right show the pressure distribution at (a) t [ 0.0148 s, (b) t [ 0.1480 s, and (c) t [ 0.1840 s).
the simulation results (right). As discussed in the experimental study, the bubble generation (break-up) procedure is divided into three phases: gas ligament expansion, gas ligament collapse, and bubble removal phase. Fig 2 shows that this procedure is similar for the simulation result. Moreover, the modeling results are helpful for analysis of the velocity field and interface behaviors as follows. As shown in Fig 2(aeg), after gas enters from the left inlet with a constant flow rate, a bubble slowly progresses through the channel pushed by the gas supplied from the gas inlet and opposed by the resistance of the liquid already present in the channel. The head of the bubble retains a convex round shape due to the contact angle with the side wall. As the gas slug arrives, the T-junction of the two inlet channels liquid enters from the liquid inlet channel, pushing the bubble to the top side of the horizontal channel. The gas forms a neck above the perpendicular channel, where bearing from the shear force is applied. This sub-process is called gas ligament expansion. The gas ligament blocks the opening of the liquid channel, and results in an increased shear force. Thus, the neck, the only transport ligament of the gas, is reduced and finally severed by the liquid. This process is so called gas ligament collapse. Fig. 2(hel) shows the exiting bubble travel through the first inverted U shape of the serpentine channel. In Fig. 2(h) and Fig. 2(k), the bubble appears as an “S” shape attached to the top side of the lower channel and to the bottom side of the
upper channel. It appears to be attempting to take the shortest route to the outlet. This same reason leads to the bubble shapes in Fig. 2(iek) and Fig. 2(m); the bubble always attaches to the inner side of the converted U channel.
3.2.
Velocity field of the break-up process
Fig 3 shows the gaseous bubble behavior during the break-up process with the velocity field. Fig 3(a) shows that the liquid flows smoothly prior to the arrival of the bubble at the Tjunction. Due to the corner effects, the velocity vector is intense around the inner corner and unconsolidated near the external corner. The velocity field during the liquid phase near the T-junction is enlarged when the gas slug attempts to cover the exit of the liquid inlet channel, as shown in Fig 3(b). With the liquid inlet channel blocked, the liquid accumulates at the exit of the liquid inlet channel and pushes the air to the top wall and forms the “neck” shape (as presented in Section 3.1). As shown in Fig 3(f), when the gas bubbles approach, but do not touch, the wall of the channel, stress on the water will strongly increase in the gap. After the ligament has collapsed, strong liquid velocity forms a vortex between the bubble released and the gas slug in the gas inlet channel. The interaction between the airflow and the liquid water can provide strong shear stress and can fragment some of the air on the gasewater interface into small bubbles. These small bubbles
2330
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 9 ( 2 0 1 4 ) 2 3 2 5 e2 3 3 3
Fig. 5 e Volume fraction of the gas and pressure drops with the vector of the first two inverted U-shaped channels.
are rolled up by the vortex and consumed by the gasewater interface again.
3.3.
Pressure fields
Figs. 4 and 5 show the volume fraction of the gas and pressure distribution in the serpentine channel. The pressure drops gradually along the channel from inlet to outlet. Because the outlet pressure is assumed to remain relatively close to a zero pressure (atmosphere), this is used as the outlet boundary condition. The pressure around the inlet changes from 7000 Pa to 5500 Pa, as shown in the three figures of different times (t ¼ 0.0148 s, t ¼ 0.1480 s and t ¼ 0.1840 s). Fig 5 shows the pressure and velocity field of the first and second inverted U-shaped channel in the red blocks that are shown in Fig 4(b). The velocity of the gas phase is larger in the middle of the bubble. The bubble, which is located at the corner, tends to flow straight and result in the large velocity towards to the wall. Although the pressure tends to decrease along the channel, the pressure inside the gas bubbles is higher than inside the liquid and results in an interesting phenomenon, a velocity vector orientated from the lower pressure field to the high pressure field.
This phenomenon shows the effects from the capillary pressure difference across the interface between the liquid and gas with the surface tension. Because the contact angle is specified as 145 in the gas phase, the gas slug in the tube shows convex surfaces on both the head and tail sides, which indicates a higher pressure inside the gas is balanced by the surface tension.
3.4.
Pressure drop
Fig. 6 shows the pressure at the gas inlet, the liquid inlet and the channel outlet. Because the outlet is specified as the outlet pressure (atmosphere), the outlet pressure remains at 0, as shown in Fig. 6(a). Fig. 6(b) shows only the pressure at the gas inlet and the liquid inlet. As presented above, the pressure generally drops with more gas enters into the serpentine channel with strong fluctuations. With the bubble that is generated in the T-junction component, the pressure fluctuation established a rhythm. The period could be defined from a bubble breakup from the Tjunction to the next bubble generated. To see this periodical phenomenon, Fig. 6(c) is shown after a Stineman Interpolation curve was fit from Fig. 6(b).
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 9 ( 2 0 1 4 ) 2 3 2 5 e2 3 3 3
2331
T-junction section. Simultaneously, the gas slug begins to cover the liquid inlet channel. This results in the peak pressure during one period.
4. Effects of surface tension and liquid viscosity 4.1.
Effects of surface tension
The two images in Fig 7 show the results of two simulation experiments. The same boundary conditions (liquid/gas inlet velocity, contact angle, etc.) are adopted for these two cases, with the exception of surface tension coefficient. The liquid velocity is 0.273 m s1, and the gas velocity is 0.0546 m s1. The surface tension is selected as the experimental measured value of 0.072 N m1 and 0.037 N m1, separately. Also, the liquid density is slightly different between the two cases (998 kg m3 and 1000 kg m3 for the first and second cases, respectively) to maintain identical conditions with the experiments. Fig 7 shows the simulation results at the exact same time instant. With the same flow rate for the two types of fluid, the lower surface tension results in fewer bubbles. After 89.5 ms, the first two small bubbles combined together to form a tadpole-like shape as presented in the experimental paper. Because the results for the experiment are not detailed enough, the analysis below only represents the modeling results. As the surface tension could be defined as the property of the surface of a liquid that allows it to resist an external force [17], the asunder bubbles are reasonable. For the bubble generated in the T-junction, as presented in previous section, the shear force from the liquid inlet channel would be easier to cut the gas slug that has lower surface tension.
4.2.
Fig. 6 e Pressure at the gas inlet, the liquid inlet and the outlet.
At the time instant of the wave trough, approximately 0.064 s, 0.098 s, and 0.13 s, the bubble breakup at the T-junction resulted in a rapid pressure drop. At the times of 0.083 s, 0.12 s, and 0.155 s, the bubble generated previously, driven by the liquid flow, must leave the
Effects of viscosity
As shown in Table 1, the three cases are simulated to explore the effects of liquid viscosity and compared to experimental results [3]. The velocities of the gas or the liquid in the three cases are identified separately as 0.182 m s1 for the gas and 0.091 m s1 for the liquid. The liquid viscosity was increased from 0.00119 Pa s for the first case to 0.00731 Pa s for the third case. The liquid viscosity of 0.00268 Pa s was chosen for the second case. The simulation results show that the bubble generated in the case with higher liquid viscosity is shorter, as shown in the figures of the final column of Table 1, which is qualitatively similar to the schematic results obtained in the experimental study as shown in Table 1. At a fixed velocity ratio of gas and liquid phase, if the viscosity increases, then the bubble length and diameter will decrease, i.e., the bubble can change the shape from a large slug to a small one. This phenomenon has also been observed in the experimental schematic drawing from Fig. 6 in Ref. [3]. And the reason for this phenomenon has been explained in Ref. [3]: the bigger viscosity will produce a larger shear stress that causes the reduced bubble dimension.
2332
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 9 ( 2 0 1 4 ) 2 3 2 5 e2 3 3 3
Fig. 7 e Gaseous volume fraction of the two cases with different surface tensions.
5.
Conclusions
In this paper, a serpentine channel with a liquid inlet and a gas inlet was adopted as the geometry used for the VOF model application to investigate bubbles behaviors by analyzing the gaseliquid interaction. (1) The entire process of bubble generation and breakup in the T-junction resulting from the simulation is almost exactly the same as the experiment investigation. (2) By analyzing the velocity field and pressure distribution from the simulation results, some phenomena, such as
the corner effect and the surface behavior of bubbles, are explained. (3) The effect of the surface tension is studied under similar conditions with the exception of the surface tension. The results show a tadpole-like shape as in the experimental study. (4) The effect of the liquid viscosity is explored by three simulations. With the increase in the liquid viscosity, the bubble length along the channel is decreased due to the stronger shear force from the liquid. (5) The VOF algorithm could be applied to bubble effects modeling in DMFC.
Table 1 e Effect of liquid viscosity. Velocity of gas inlet (ms1)
Velocity of liquid inlet (ms1)
Liquid viscosity (Pa s)
0.182
0.091
0.00119
0.182
0.091
0.00268
0.182
0.091
0.00731
Bubble shape presented in experimental results [3]
Simulation results
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 3 9 ( 2 0 1 4 ) 2 3 2 5 e2 3 3 3
Acknowledgments The authors are grateful for the support of this work by the Auto21 Networks of the Centres of Excellence (Grant D303DFC), the Natural Sciences and Engineering Research Council of Canada (NSERC), the Canada Foundation for Innovation (CFI), the Ontario Innovation Trust (OIT), and the University of Windsor.
Nomenclature r u S s m s k K n R1, R2
Density (kg m3) Velocity vector (m s1) Source term Phase volume fraction Dynamic viscosity (Pa s) Surface tension (N m1) Surface curvature Permeability (m2) Normal vector Two radii in the orthogonal directions
Subscripts l Liquid phase g Gas phase
references
[1] Lu GQ, Wang CY. Electrochemical and flow characterization of a direct methanol fuel cell. J Power Sources 2004;134:33e40. [2] Wong CW, Zhao TS, Ye Q, Liu JG. Transient capillary blocking in the flow field of a micro-DMFC and its effect on cell performance. J Electrochem Soc 2005;152(8):A1600e5.
2333
[3] Xiong R, Chung JN. Bubble generation and transport in a microfluidic device with high aspect ratio. Exp Therm Fluid Sci 2009;33:1156e62. [4] Fu BR, Pan C. Bubble growth with chemical reactions in microchannels. Int J Heat Mass Transs 2009;52:767e76. [5] Quan S, Lou J. Viscosity-ratio-based scaling for the rise velocity of a Taylor drop in a vertical tube. Phys Rev E 2011;84(3):036320. [6] Quan S. Co-current flow effects on a rising Taylor bubble. Int J Multiph Flow 2011;37:888e97. [7] Kang C, Quan S, Lou J. Numerical study of a Taylor bubble rising in stagnant liquids. Phys Rev E 2010;81(6):066308. [8] Fei K, Chen TS, Hong CW. Direct methanol fuel cell bubble transport simulations via thermal lattice. J Power Sources 2010;195:1940e5. [9] Quan Peng, Zhou Biao, Sobiesiak Andrzej, Liu Zhongsheng. Water behavior in serpentine micro-channel for proton exchange membrane fuel cell cathode. J Power Sources 2005;152:131e45. [10] Kulikovsky AA. Model of the flow with bubbles in the anode channel and performance of a direct methanol fuel cell. Electrochem Commun 2005;7:237e43. [11] Kulikovsky AA. Bubbles in the anode channel and performance of a DMFC: asymptotic solutions. Electrochim Acta 2006;51:2003e11. [12] Wang ZH, Wang CY. Mathematical modeling of liquid-feed direct methanol fuel cells. J. Electrochem Soc 2003;150(4):A508e19. [13] Ge J, Liu HT. A three-dimensional two-phase flow model for a liquid-fed direct methanol fuel cell. J Power Sources 2007;163:907e15. [14] Le Anh D, Zhou Biao. A general model of proton exchange membrane fuel cell. J Power Sources 2008;182:197e222. [15] Le Anh Dinh, Zhou Biao, Shiu Huan-Ruei, Lee Chun-I, Chang Wen-Chen. Numerical simulation and experimental validation of liquid water behaviors in a proton exchange membrane fuel cell cathode with serpentine channels. J Power Sources 2010;195:7302e15. [16] Fluent 6.3 manual. Fluent Inc; 2006. [17] Wikipedia, http://en.wikipedia.org/wiki/Surface_tension.