RESEARCH PAPER
Microfluidic analysis of CO 2bubble dynamics using thermal lattice-Boltzmann method
K.Fei ÆW.H.Chen ÆC.W.Hong
Received:17August 2007/Accepted:10October 2007/Published online:30October 2007ÓSpringer-Verlag 2007
Abstract By means of microfluidic analysis with a thermal lattice-Boltzmann method,we investigated the hydrophilic,thermal and geometric effects on the dynamics of CO 2bubbles at anode microchannels (e.g.,porous layers and flow channels)of a micro-direct methanol fuel cell.The simulation results show that a more hydrophilic wall provides an additional attractive force to the aqueous methanol in the flow direction and that moves the CO 2bubble more easily.The bubble propagates quicker in the microchannel with a positive temperature gradient impod from the inlet to the exit,mainly due to the Marangoni effect.Regarding the geometric effect of the microchannel,the bubble moves more rapidly in a divergent microchannel than in a straight or convergent channel.On the basis of the quantitative evaluation of hydrophilic,thermal and geo-metric effects,we are able to design the bubble-removal technique in micro fuel cells.
Keywords Thermal lattice-Boltzmann method ÁBubble dynamics ÁTwo-pha flow ÁMicro-direct methanol fuel cell
List of symbols e lattice velocity vector e lattice speed (cm/s)f density distribution function (g/cm 3)g thermal distribution function (g K/cm 3)G rr 0
interaction strength between the species r and the
other species r 0(cm 3/g s)
G r fluid–solid interaction potential parameter of the species r
G gravitational constant (cm/s 2)G 0non-dimensional constant T temperature (K )t time (s )
U velocity vector(u ,v )
u velocity in the x -direction (cm/s)v velocity in the y -direction (cm/s)X
position vector (x ,y )
Greek symbols b thermal expansion coefficient (1/K)q density (g/cm 3)h contact angle (°)s collision ti
me for momentum transfer s T collision time for energy transfer w function of the mass density Subscripts i lattice velocity directions ?reference state Superscripts eq equilibrium r species
1Introduction
Micro-direct methanol fuel cells (l DMFC)are considered a strong competitor of future power sources for portable equipment.The advantages are high efficiency,high power density,low operation temperature and almost zero pollu-tion;but some technical challenges remain before commercialization is practicable.A major issue is the
K.Fei ÁW.H.Chen ÁC.W.Hong (&)
Department of Power Mechanical Engineering,National Tsing Hua University,Hsinchu 30013,Taiwan
e-mail:hu.edu.tw
Microfluid Nanofluid (2008)5:119–129DOI 10.1007/s10404-007-0232-x
removal of bubbles of carbon dioxide(CO2)at the anode side,which are generated there on oxidation of methanol; the bubbles obstruct the diffusion layer and theflow channel if they are not removed e
fficiently,leading to l DMFC malfunction.The behavior of CO2bubbles at the anode has thus drawn much attention.One feasible way to investigate the bubbleflow is through experimentalflow visualization(Lu and Wang2004;Yang et al.2005;Wong et al.2005).Most rearchers focud on the bubble dynamics inside theflow channels of the fuel cell.The experimental approach is unable to provide detailed insight into the diffusion layer becau of technical difficulties in performing that visualization.Computer simulation is thus an applicable approach to study bubble dynamics from a micro-scale viewpoint.We constructed mathematical models using a thermal lattice-Boltzmann method(TLBM) to obrve the phenomenon of CO2bubbles in the anode microchannel of a l DMFC.The main focus of this rearch is on the thermal and geometric effects of the microchannelflow from a microfluidic viewpoint,the hydrophilic effect was inherited from the original lattice-Boltzmann method in our previous publications(Fei et al. 2006;Fei and Hong2007).
The thermal effect influences the bubble transport sig-nificantly becau the surface tension of a gas bubble depends on temperature.The temperature gradient around the gas bubble leads to a difference of surface tension around the bubble surface and caus a tangential force at the interface.The direction of the tangential force is opposite the temperature gradient,dragging thefluid around the bubble surface from a region of high temperature to a region of low temperature(Marang
oni effect).A drag force exerts an opposite reaction on the bubble and makes it move from a cold region to a warm region.This phenomenon is evident when the size of the bubble is on a micro-scale and the temperature gradient is large.The movement of a bubble caud by a temperature gradient is known as ther-mocapillary motion and wasfirst obrved experimentally by Young et al.(1959).The Marangoni effect is a driving force of the bubble movement and has been widely utilized in the micro-flow device(Jun and Kim1998;Takahashi et al.1999,2001).Although the lattice-Boltzmann method has undergone rapid progress in simulating multi-pha flow,little rearch has been conducted to discuss the dynamic distribution of temperature in a multi-phaflow including the thermal effect.Shan(1997)propod a pas-sive-scalar thermal lattice-Boltzmann model that treated temperature as a component of the mixture.Bad on an analogy between mass and heat transfer,the temperature field was obtained on solving an additional equation similar to the original lattice-Boltzmann equation;compression work and viscous dissipation of heat were both neglected in this model.The momentum and energy equations were coupled by an external body force through the temperature gradient in the thermal lattice-Boltzmann equation.Yuan and Schaefer(2006a,b)combined this thermal lattice-Boltzmann model with a multi-phaflow model propod by Shan and Chen(1993)to simulate the static distribution of temperature in a two-phaflow.Shi et al.(2004)derived a single-pha thermal lattice-Boltzmann model from the Boltzman
n equation and the Maxwell–Boltzmann distri-bution,and obtained the model after applying the Bhatnagar–Gross–Krook(BGK)assumption.This model was capable of simulating the incompressible thermalflow with viscous heat dissipation and a large Prandtl number. Recently,Guo et al.(2007)developed another thermal lattice-Boltzmann for low Mach numberflows from a two-relaxation-time kinetic model.They defined an additional distribution function to reprent the total energy of theflow so that the compression work and viscous dissipation can easily be included in the model.In this paper,we derived a thermal lattice-Boltzmann model for a two-phaflow in the anode microchannel of methanol fuel cells.A simpli-fied,2D microchannel was t up in this rearch to mimic the realflow channel in a porous diffusion layer.Simula-tions were conducted by alteringflow parameters one at a time to investigate the effects on the two-phaflow in the microchannel.Tho parameters include the degree of hydrophilicity,the wall temperature distribution and the geometry of the microchannel.The main purpo is to investigate the bubble removal technique in the two-pha flow of gaous CO2bubbles and the liquid of methanol–water solution.
2Thermal lattice-Boltzmann method
Details of a2D,nine-velocity(D2Q9)lattice-Boltzmann model for a two-phaflow,including the treatment of boundary conditions,were explicitly described in our previous work(Fei et al.2006;Fei an
d Hong2007).We derived the lattice-Boltzmann model from the physically-oriented concept propod by Shan and Chen(1993).The main idea of this model was to treat additional effects of the two-phaflow as source terms of the momentum equation.In addition to the momentum equation,an energy equation is required to calculate the temperaturefield in this work.The thermal lattice-Boltzmann equation for heat transfer has a form identical to that of the previous lattice-Boltzmann equation and is expresd as
g iðXþe i D t;tþD tÞÀg iðX;tÞ¼À
1
s T
g
i
ðX;tÞÀgðeqÞ
i
ðX;tÞ
h i
ð1Þin which g i(X,t)is the thermal distribution function that is a function of space X,and time,t;e i lattice velocity vector
in the i direction;s T dimensionless collision time for heat
transfer;and the superscript eq denotes an equilibrium state.The mathematical form of the thermal equilibrium distribution functions in each direction,g (eq)i
,was derived by Shi et al.(2004)as the general equation of
g ðeq Þ
i ¼w i q T 1þ3e i ÁU e þ
92ðe i ÁU Þ2
e À32ðU ÁU Þe "#
ð2Þ
in which w i =4/9for i =0,w i =1/9for i =1–4and
w i =1/36for i =5–8;U is the particle velocity vector that includes the velocity components in both x and y directions;q is the density at a certain position in the flow field.The thermal energy density is defined as the sum over the thermal distribution functions and is expressible as
q r ðX ;t ÞT r ðX ;t Þ¼X
i
g r i ðX ;t Þ
ð3Þwhere superscript r denotes the species,r .The temperature of each species T r (X ,t ),can be calculated from the thermal energy density divided by the mass density q r (X ,t ).
The fluid density can be approximated as an inverly linear function to the temperature over a small temperature range.In our ca of the l DMFC,the operation range is from 298to 333K.The additional thermal force caud by the density difference due to the temperature variation is
F r thermal ðX ;t Þ¼s T q r 1ðX ;t Þb r
G ðT r ÀT r
1ÞÂÃð4Þin which b r is the thermal expansion coefficient of species r ;G is the gravitational constant;q ?and T ?are the reference density and temperature,respectively.This additional force can be added to the original source term of the momentum equation which includes the surface tension,the fluid–solid wall interaction and the buoyancy force.The new general momentum equation is expresd as
q r ðX ;t ÞU r ðX ;t Þ¼X
i
f r i ðX ;t Þe i þF r
total ðX ;t Þð5Þ
in which f is the mass density distribution function and the total source term is
F r total ðX ;t Þ¼F r surface tension ðX ;t Þþ
F r
solid ðX ;t Þ
þF r
buoyancy ðX ;t Þþ
F r
thermal ðX ;t Þ
ð6Þ
Each external force,except the thermal force in Eq.(4),
is expressible as
F r
surface
tension ðX ;
t Þ
¼Às r
w r ðX ;t ÞX
r ‘
G rr 0X i
w r 0ðX þe i D t ;t Þe i
"
#
ð7Þ
F r
solid
ðX ;t Þ¼Àq r
ðX ;t ÞX
i
G r i s ðX þe i D t Þe i
ð8ÞF r
buoyancy ðX ;t Þ¼G 0
X
i
q r ðX þe i D t ;t Þe i
ð9Þ
In the above equations,w r (X ,t )is a function of the mass
density of species r at position vector X ;G rr 0
reprents the strength of the interaction between different species and is t constant for simplicity;s r is the collision time of the species r for momentum transfer;G r is a fluid–solid interaction potential parameter;s is a function of particle position (s =0when the particle is in the fluid,s =1when the particle is at the fluid/solid interface);G 0is a non-dimensional constant.Shan and Chen (1994)suggested that the form of w r (X ,t )can be reprented by w r ðX ;t Þ¼A exp ðÀq 1=q Þ
ð10Þ
in which A is a constant and is t to be 5.4in the following simulations.
3Boundary conditions for TLBM
To determine the unknown thermal distribution functions at the boundary nodes,thermal boundary equations were derived following the principle given by Inamuro et al.(2002).Figure 1shows the known and unknown thermal distribution functions at the inlet,exit,upper and lower walls of the microchannel.
3.1Inlet and exit boundaries
Assuming the inlet flow temperature T in ,to be known at the
boundary node,the energy equation in the TLBM
is
Fig.1Thermal distribution functions,g ,(i =0–8),at boundaries of
a microchannel;solid lines are known conditions and dashed lines are unknown ones
q T in¼
X
i
g ið11Þ
The summation of the unknown thermal distribution functions is expressible as
g1þg5þg8¼1
6
q T01þ3
u
e
þ3
u2
e
!
ð12Þ
in which T0is the unknown temperature.Substituting Eq.
(11)into Eq.(12),the unknown temperature,T0,is obtained as
T0¼
6
q1þ3u
e
ÀÁ
þ3ðu
e2
Þ
Âýq T inÀg0Àg2Àg3Àg4Àg6Àg7
ð13Þ
Tho unknown thermal distribution functions of the inlet nodes are
g1¼1
9
q T01þ3
e1ÁU
e2
þ
9
2
ðe1ÁUÞ2
e4
À
3
2
U2
e2
"#
ð14Þ
g5¼1
q T01þ3
e5ÁU
þ
9ðe5ÁUÞ2
À
3U2
"#
ð15Þ
g8¼1
36
q T01þ3
e8ÁU
e
þ
9
2
ðe8ÁUÞ2
e
À
3
2
U2
e
"#
ð16Þ
The unknown distribution functions at the exit boundary, g3,g6and g7,can also be determined with the similar procedure as above.
3.2Wall boundaries
In Fig.1,tho unknown thermal distribution functions of the boundary nodes at the bottom wall are obtainable with a similar procedure and are expresd as
g2¼1
9
q T0ðXÞ1þ3
e2ÁU
e2
þ
9
2
ðe2ÁUÞ2
e4
À
3
2
U2
e2
"#
ð17Þ
g5¼1
36
q T0ðXÞ1þ3
e5ÁU
e
þ
9
2
ðe5ÁUÞ2
e
À
3
2
U2
e
"#
ð18Þ
g6¼1
36
q T0ðXÞ1þ3
e6ÁU
e
þ
9
2
ðe6ÁUÞ2
e
À
3
2
U2
e
"#
ð19Þ
in which the unknown temperature T0(X)is
T0ðXÞ¼
6
q1þ3ðv
e
Þþ3ðv
e2
ÞÂÃ
½q T wallðXÞÀg0Àg1Àg3Àg4Àg7Àg8
ð20Þ
The wall temperature,T wall(X),is a function of the node position along the wall or it can be a constant wall temperature.The effect of wall temperature on theflow can be investigated in the following ction.The unknown distribution functions of the upper wall can also be expresd in a similar way as above,except that g4,g7 and g8are unknowns.
4Simulation results
The example simulation domain in our ca is a straight microchannel of height 1.5l m reprenting the pore diameter of a diffusion layer and length15.9l m,which is ten times the height to make the2D assumption valid(as shown in Fig.2a).Figure2b is the grid generation of the D2Q9TLBM model,in which D x and D y are both t to be 0.05l m.There are in total9,540grid cells in the example simulation ca.
4.1Effect of hydrophilicity
The hydrophilicity between thefluid and the wall surface is controlled by adjusting the value of thefluid–solid inter-action potential parameter,G r in Eq.(8),in this simulation.
A larger value of G r implies the microchannel wall to be less hydrophilic.Physically,the degree of hydrophilicity can be reprented by the contact angle in the two-pha flow.The definition of the contact angle of the bubble in a microchannel is defined as the angle between the solid wall and the line tangent to the bubble surface at the contact edge,as shown in Fig.3a.The total surface energy E s of the2D two-phaflow in the microchannel is expresd
as Fig.2a Simulation ,length L=15.9l m,height H=1.5l m)of a straight horizontal microchannel reprenting the flow path of the diffusion layer and theflow channel;b grid generation of the D2Q9model(D x=D y=0.05l m),the total number of grid cells is9,540
E s ¼x c S 1l þx c S 2l þðL Àx Þc S 1g þðL Àx Þc S 2g
ð21Þ
in which x is the length of the liquid pha,L length of the microchannel,c S i j (i =1,2;j =g ,l )is the solid–liquid or solid–gas surface energy of the upper wall S 1,and the lower wall S 2(e Fig.3a).Young’s equation (Young 1805)describes the relation between the surface energy and the contact angle h ,c sg ¼c sl þc gl cos h ;
ð22Þ
where c sg ,c sl and c gl are the solid–gas,solid–liquid and gas–liquid surface energies,respectively.After substituting Eq.(22)into Eq.(21)and differentiating with respect to the variable x ,the forward pressure drop across the gas–liquid interface acting on the rear of the bubble D P gl ,is obtained as
D P gl ¼c gl
H ðcos h 1þcos h 2Þð23Þ
in which H is the height of the horizontal microchannel;h 1is the upper contact angle and h 2is the lower contact angle.The pressure drop across the gas–liquid interface evidently depends on the contact angles (or the hydrophilicity at the upper and lower walls),the gas–liquid surface tension and the channel width.When a bubble moves in a micro-channel,as illustrated in Fig.3b,an advancing contact angle is formed at the front contact edge between the bubble surface and the channel wall,while a receding contact angle is formed at the rear contact edge.The advancing contact angle may be different from the receding contact angle depending on flow conditions and surface property at the solid wall.According to the above analysis,when the advancing contact angle of the moving bubble is smaller than the receding contact angle,the bubble sustains
a larger backward pressure drop than the forward pressure drop;the transport of the bubble is thus obstructed.
To study the effect of the hydrophilicity,we cho a t of G r values from –0.007to 0.007in the step of 0.001,which are related to the contact angle as shown in Fig.4.The G r values and other flow pa
rameters ud in the simulations are listed in Table 1.In Fig.4,the contact angle of the bubble increas from 10.0°to 51.0°with G r increasing from –0.007to 0.007,indicating that the larger the value of G r ,the greater the contact angle,the less the hydrophilicity.The velocity of the bubble movement in the microchannel also varies with hydrophilicity.The bubble velocity is calculated on recording the traveling distance of the bubble center during a certain interval.When the channel wall is hydrophilic,the wall provides an attractive force on the liquid,making the liquid flow more rapidly than in a less hydrophilic microchannel.The bubble movement in the liquid is likewi affected;it transports with greater velocity if the liquid flows more rapidly.In addition,the pressure drop across the gas–liquid interface,D P gl ,increas if the microchannel is more hydrophilic (Eq.23).The CO 2bubble in a hydrophilic microchannel is thus able to transport more rapidly than in a hydrophobic microchannel,as indicated by Fig.4.The first conclusion to be drawn here is that a hydrophilic microchannel is preferable for the purpo of bubble removal at the anode of the l DMFC.
Further analysis regarding the effect of hydrophilicity is conducted on imposing a linear gradient of hydrophilicity on both walls of the microchannel.From the previous results,the bubble velocity is 280.2l m/s when the contact angle is 10.0°(G r =–0.007),and the velocity is 266.9l
m/
Fig.3a Diagram of the contact angle definition of bubble flow in a microchannel (D P gl forward pressure drop across the gas–liquid interface;L length of the microchannel;H height of the horizontal microchannel;S 1upper surface of the microchannel;S 2lower surface of the microchannel;x length of the liquid pha;h 1upper contact angle;h 2lower contact angle);b illustration of a moving bubble in the microchannel (h rec the receding contact angle;h adv the advancing contact
angle)
Fig.4Contact angle and bubble velocity versus the varied value of the fluid–solid interaction strength G r .The diagram indicates that the larger the value of G r ,the greater the contact angle,the less the hydrophilicity is,and the lower the bubble transport velocity
