Excerpt from the Proceedings of the COMSOL Conference 2010 ParisA Model of Gas Bubble Growth by Comsol MultiphysicsB.Chinè*1,2 and M. Monno1,3Laboratorio MUSP, Macchine Utensili e Sistemi di Produzione, Piacenza, Italy; 2Instituto Tecnològicode Costa Rica, Escuela de Ciencia e Ingenieria de Materiales, Cartago, Costa Rica; 3Politecnico di Milano,Dipartimento di Meccanica, Milano, Italy.1*Corresponding author: Laboratorio MUSP, Località Le Mose, SS 10, 29100 Piacenza, Italy; firstname.lastname@example.orgAbstract: We use Comsol Multiphysics tomodel a gas bubble expansion in a viscous liquidinitially at rest, a very common system forlightweight foamed materials from metalproduction and polymer processing. The aim ofthe present work is to develop a firstcomputational model for the growth of gasbubbles under simpler conditions, modeling boththe gas and liquid flow and to verify its validityby comparing the numerical results with existinganalytical solutions. The two dimensionalisothermal model developed in ComsolMultiphysics considers a gas bubble growing dueto a pressure difference with a surroundingNewtonian liquid. Surface tension effects on thegas-liquid interface are considered. The modelequations are solved on a fixed grid, built both inthe gas than in the liquid region. In order tocapture the front between the two fluids weexploit the capability of the level set method.The numerical results of the computationalmodel compare well with analytical solutionsfrom theory and obtained for a few simple cases.This first computational work is a basis forconsidering successive more realistic foamexpansions.Keywords: Bubble growth, multiphase flow,level set.1. IntroductionGas bubble growing in a surrounding liquidmatrix is an important and complex phenomenonin many technological fields. An interesting caseis encountered in cellular metals and metal foamsproduction, when nucleated gas bubbles expandand move in a confined liquid metal before tocool and solidify. A liquid metal (e.g. Al) can befoamed directly by injecting gas (H2) or gasreleasing blowing agents (solid particles), or byproducing supersaturated metal-gas solutions .The process depends on simultaneous mass,momentum and energy transfer between threephases: solid, liquid and gas. Furthermore, otherphysical phenomena should be taken in accountin the system, like complex interface processes,bubble motion, coexistence, coalescence andcollapse of bubbles. Experimental works carriedout by observation techniques cannot besometimes applied owing to the specificproperties of liquid metals: they are hot, opaqueand very reactive with oxygen. Then, thesemechanisms can be modelled and studied byapplying computational techniques, although thecomputational work is very challenging: thephenomenon are not independent among themand many times are simultaneous. Multiphaseflow modelling is also required because of thepresence of more phases: the solid metal matrixcontaining the foaming agent not yet fullymelted, the gaseous phase made up of bubblesand the liquid phase composed of a melted metalmatrix. Mass and heat transfer should be alsotaken into account in the modelling work to givethe most accurate results. On the other handthese phenomenon have major effects on thequality of metal foams. To give an example,during mould filling the desired metal foamdensity is very strategic and this is dependent onthe ability of controlling the gas bubble growth.It is understood that the improving offoaming quality and cost-effectiveness may bealso realized by means of simplifiedcomputationalmodels. Attention will bedemanded by the presence of a dynamic interfacebetween the gas bubbles continuouslyoriginating and the surrounding liquid. Toaccurately compute the evolution of interfaces,many computational methods are been designedin the past. Lagrangian methods use a numericalgrid which follows the fluid and tracks theinterface, while by Eulerian methods theinterface is captured on a stationary grid. Oneinteresting alternative Eulerian numericalformulation is provided by the level set method,which embeds the interface as the zero level setof a function. The method was first introduced
by Osher and Sethian in 1988  and hasencountered extensive applications in multiphaseflow modelling. Level set techniques have theadditional advantage that they can easily provideaccurate values for the normal direction and thecurvature of a physical interface.As a first effort to understand metal foamingmechanisms, we consider the growth of a gasbubble embedded in a viscous liquid. We takeinto account the interface movement which isdue to a pressure difference between the twophases. To model and solve the governingequations of the problem we will use ComsolMultiphysics, a commercial simulation toolbased on finite elements method which includesthe level set method in its Chemical EngineeringModule .ypEXT(t)LiquidRΩR(t)θpG(t)xGasnFigure 1. Schematic of a gas bubble growth in aliquid region with the considered symmetry: at thebeginning R(t) R0, pG(t) pG, 0 and pEXT(t) pEXT, 0.2. Theory and model descriptionIn this section we describe the twodimensional model and the assumptions whichare been done to reasonably simplify theproblem. This is due also for managing thecomputational job in a common laptop machine.2.1 TheoryLet us to consider a gas bubble growing dueonly to a pressure difference with a surroundinglimited amount of liquid matrix. The problem isassumed two dimensional, no heat transfer andmass diffusion are taken into account. Thesystem is thus isothermal and there are notgradients of species concentration. The gas in thebubble follows the ideal gas law while the liquidis considered to be an incompressible Newtonianfluid. Furthermore, the gas and the liquid areassumed to be immiscible.Figure 1 shows a schematic of a circular gasbubble with initial radius and pressurerespectively equal to R0 and pG,0. The expansiontakes place in a liquid matrix Ω modeled as acircular region of radius RΩ. For symmetricalconsiderations it is possible to reduce themodeled region, both for the gas bubble and theliquid (see section 2.2). The gas liquid interfaceis a free surface with uniform surface tensioncoefficient σ, with κ representing the localcurvature and n the unit normal to the interface.At any time t, T is the constant absolutetemperature of the system and R(t) and pG(t) arethe bubble radius and gas pressure, respectively.At the beginning ( t 0 ), with the liquid atrest and by calling with pEXT, 0 the initial ambientpressure imposed on the boundary of Ω, theLaplace’s equation states that stress balance atthe surface of a circular bubble of radius R0 is:pG,0 pEXT, 0 (1)where: 1 / R0(2)for a circular bubble. Substituting Eq.2 into Eq.1,it is written as:R0 pG ,0 pEXT,0(3)which gives the radius of the bubble at t 0when the liquid is still at rest. Consequently,given σ, the bubble will expand or contract andthe liquid will flow if the pressure differencepG,0 pEXT, 0 changes its value. For a gas bubbleexpansion, a common case could be representedby a sudden lowering of the ambient pressure.Under this condition, a new equilibrium state
will be reached by the bubble and, with theliquid at rest again, the bubble radius atequilibrium is:Req pG pEXT(4)Here pG and pEXT are the new equilibrium valuesof the pressure, for the gas and the ambientrespectively. On the other hand, the ideal gasequation followed by the specie inside thebubble is:pV n T(5)where V represents the volume (area for a twodimensional problem) of the bubble, is theuniversal gas constant and n is the number ofmoles. We assume that during bubble expansionthe behavior of the gas is polytropic, thus for anisothermal process it follows:pG,0 A0 pG A(6)being A0 and A the equilibrium bubble areas forthe pressure pG,0 and pG , respectively.2.2 Governing equationsTo simulate numerically the isothermalgrowth of a gas bubble embedded in a viscousliquid and the fluid flows that originate, we usethe classical equations of fluid dynamics coupledto the level set method available in ComsolMultiphysics. As we said before, the method isvery well suited to describe the motion of theinterface during the gas expansion. We assumethe liquid to be an incompressible Newtonianfluid and take into account the compressibility ofthe gas in the bubble. For gas flows with lowMach numbers (approximately Ma 0.3 ), aweakly- compressible model can be used. In thismodel, the gas density ρG(t) is given by the idealgas law , after introducing the molar mass Mand mass m ( n m M ). Then, for both thefluids, the coupled partial differential equationsof the model are the following (Two Phase Flow,Level Set Application Mode, ): ( u) 0 t(7) u (u )u [ pI ( u ( u)T ) t (2 DV )( u)I] F g FST3 ] u [ (1 ) t (8)(9)In the momentum transport equation, the scalarmagnitudes ρ, η and DV are the fluid density,dynamic viscosity and bulk viscosity,respectively. Among the other terms, u is thefluid velocity, I is the identity tensor, ρg is thegravity force and F takes into accounts otherbody forces. The term FST accounts for thesurface tension force acting at the interfacebetween the two fluids (see  for more detailson FST modeling).The advection of the level set function inthe computational region is given by Eq. 9. Tocompute this field, a signed distance function att 0 is used to build the function whichcorresponds to the interface at the level set 0.5 . In the same step the values of insidethe two phases are setting as 0 0.5 for onefluid (in our model is the liquid) and 0.5 1for the other (gas in the model). The parameter γrepresents the reinitialization parameter andcontrols the re-initialization performed at somelater point in the calculation beyond t 0 , needto preserve the values of distance close to theinterface. Finally ε is the interface thicknessparameter which adds extra numerical diffusionin order to stabilize the computations of Eq. 9.To solve the governing equations we haveselect a Cartesian system of coordinates (x,y) andapplied the symmetry conditions shown inFigure 1. The computational domain Ω has beenreduced again, giving a circular surface withazimuthal angle 0 4. In the domain theexternal boundary of Ω has been set as anoutflow with a Dirichlet condition pEXT (t ) 0on the pressure and vanishing viscous stresses.At t 0 , starting with gas and liquid phasesat equilibrium, we impose a pressure valuegreater than R0 on the gas bubble. Forces ρgand F are set to zero, then we drive the interfacemotion by analytically computing the density
variation ρG(t). Using cylindrical coordinates(r,θ), the incompressible liquid motion is radiallyoriented and the velocity field u is only afunction of (r,t). The interface velocity may beobtained by the velocity field at r R, this leadsto:u (r , t ) dR(t ) R(t ), r Rdtr(10)In the absence of tension surface effects andconsidering the stress at r R, the bubble radiusmay be calculated by solving an ordinarydifferential equation, as explained in .Successively, for pEXT 0 the density variationis computed as: G (t ) G , 0(1 pG ,0 t L(11))Instead, for 0 we solve the respectiveordinary differential equation and calculate R(t),numerically, by the solution:[C AR(t )]C exp AR(t ) [C AR0 ]C exp A( R0 At)(12)2pG , 0 R0 where C and A .2 L2 LFinally, the radius R(t) from Eq.12 is substitutedinto Eq.6, modeling a density change for anisothermal process with surface tensionphenomenon at interface.3. SimulationsOur main interest is to simulate a gas bubbleexpansion in a liquid metal during foamprocessing, in order to capture the interface andanalyze fluid flows similar to those we couldhave within moulds. However, multiphasephenomenon in metal foam production givestrong property gradients at fluid interfaces,which cause computations to be carried out withsomedifficulties.Consideringthesecomplexities, we start to study simpler problems,as first examples of future more realisticsimulations. Consequently, we have set theTable 1: Fluid properties used in the simulations ofgas bubble growth in a liquid.MagnitudeUniversal gasconstantGas molar massGas densityLiquid densityGas viscosityLiquid viscosityGas bulkviscositySurface tensioncoefficientInitial bubbleradiusInitial bubblepressureAmbient pressureConstanttemperatureSymbol MρGρLηGηLκDVσR0pG,0pEXTTValue8.314 J/(mol·K)2 g/molEq.5, Eq.6, Eq.1110 kg/m310-3 Pa·s10-1 Pa·s0 Pa·s0 N/m10--2 N/m10-2 m0.2 Pa1.2 Pa; 2.2 Pa0 Pa933 KTable 2: Model parameters used in the simulations ofgas bubble growth in a liquid.MagnitudeMax element sizeof the meshTime steppingRelative ializationSymbol-Value10-4 m-set by the solver10-3 s10-4 sε10-4 mγ0.01 0.02 m/svalues of fluid properties shown in Table 1.Other parameters used in the model are given inTable 2. In this way we model flows withmoderate density and viscosity differences, aswell as surface tension. The density of the liquidis taken to be 10 kg/m3, while the initial gasdensity ρG,0, computed by Eq. 5, is 0.02613kg/m3 which gives a density ratio ofapproximately 4x102. In the same way theviscosity ratio is near 102. The surface tensioncoefficient is zero in the first simulated case andequal to 0.01 N/m in the other computations.
To solve the model equations, the globalcomputational domain Ω has been meshed by 104triangle elements approximately, correspondingto more than 8·104 degrees of freedom, numberwhich could be still managed by a commonlaptop. Finally, calculations have been carriedout with the direct solver PARDISO, ComsolMultiphysics version 3.5a. Although the partialdifferential equations of the model are non linearand time dependent, the converge obtainedduring computations was good, giving a stepsize near to 10-3 s with a solution time of 1.8hours approximately, for a laptop with 2.8 GHzIntel Core2 Duo processor and 4 GB RAM.The simulations started by initializing thelevel set function, using a pseudo-time equal to10-2 s, such that varies smoothly from zero toone across the interface. Then, the transientmode was selected and computations werecarried out in order to get bubble growth,interface motion and fluid flow, for threedifferent conditions:Figure 2. Growth of a gas bubble in a liquidcorresponding to the simulated case b) of section 3.a) σ 0 N/m, p 0.2 Pab) σ 10--2 N/m, p 1.2 Pac) σ 10--2 N/m, p 2.2 Pa4. Results and discussionFigure 2 shows the bubble growth in theliquid region when tension surface is present andthe initial gas pressure is equal to 1.2 Pa. Thevolume gas fraction is plotted 1.4 s after thegrowth starts, when the bubble has alreadyreached its maximum expansion. We observethat the final bubble radius is equal to 0.012 m,in agreement with the analytical value given byEq.4. For the same conditions, Figure 3 givesthe pressure field both for the gas and liquid,which agrees with the values of pressure for aliquid again at rest ( pL pEXT ) and with a static 0.01pressure jump equal to 0.83 PaReq 0.012Figure 4 presents the change of areas, over time,obtained by integrating the level set function in the computational region. As shown in thefigure, the gas area is growing only until to acertain time value ( t 1s ), stating that fort 1.4 s the bubble is in equilibrium with theliquid. Due to the expansion, the fluids acquire aFigure 3. Pressure field in the gas and liquidcorresponding to the simulated case b) of section 3.significant motion, with the characteristic radialpattern depicted in Figure 5. In the figure wehave plotted the pressure field superimposed onthe contour of fluid velocity at t 0.2 s , for agas bubble with an initial pressure equal to 2.2Pa. As expected, the magnitude of the gasvelocity increases with decreasing distance fromthe gas-liquid interface, where it reaches itsmaximal value, which is equal to 0.02 m/sapproximately. Continuity between gas andliquid velocity can be well noted, indicating acorrect coupling of the model at interface. Att 0.2 s the bubble is still expanding, far from
5. ConclusionsFigure 4. Values of total (black line), gas (browndashed line) and liquid (blue dashed line) area,corresponding to the case b) of section 3.A model by Comsol Multiphysics has beenpresented for the simulation of gas bubblegrowth with flow in the liquid and inside the gasregion. No mass and energy transfers has beenconsidered, bubble growth is due to only apressure difference. Flows in both the fluids arecalculated using a weakly-compressible modelcoupled to a level set transport equation whichcaptures the front between the two phases. Gasdensity has been modeled starting from the gasideal equation, also in presence of surfacetension effects. The resulting model performswell, giving results which agree with analyticalsolutions. Although the model takes intoaccounts moderate density and viscositydifference values for the fluids, simplifying inthis way the numerical computations, itrepresents a good basis for future and morerealistic simulations of foam expansions.6. ReferencesFigure 5. Pressure and velocity field during thegrowth of a gas bubble in a liquid, corresponding tothe case c) of section 3: t 0.2 s after growth starts.the static equilibrium. In fact, the radius for thislatterconditionisgivenby2pG , 0 R0Req 0.022 m , which is obtained introducing Eq.4 into Eq.6. In correspondence ofthat radius the pressure would be equal to 0.01 0.45 Pa . This value is lower thanReq 0.022the maximum pressure present at that momentwithin the bubble, close to 0.73 Pa, as we mayobserve in the same Figure 5. Banhart J., Manufacture, characterization andapplication of cellular metals and metal foams,Progress in Materials Science, 46, 559-632(2001). S. Osher and J.A. Sethian, Fronts propagatingwith curvature dependent speed: Algorithmsbased on Hamilton-Jacobi formulation, Journalof Computational Physics, 79, 12-49 (1988). Comsol AB, Comsol Multiphysics-ChemicalEngineering Module, User’s Guide, Version 3.5a, (2008). J. Bruchon, A. Fortin, M. Bousmina and K.Benmoussa, Direct 2D simulation of small gasbubble clusters: From the expansion step to theequilibrium state, Int. J. Numer. Meth. Fluids,54, 73-101 (2007).7. laboration provided by Dr. Valerio Marrawith Comsol AB, Brescia, Italy.
Abstract: We use Comsol Multiphysics to model a gas bubble expansion in a viscous liquid initially at rest, a very common system for lightweight foamed materials from metal production and polymer processing. The aim of the present work is to develop a first computational model for the growth of gas