Modeling Pulse Properties near the Bubble-to-Pulse Transition in Randomly Packed Beds

Traveling wave analysis of a recently developed two-fluid model for bubbly flow in lab-size packed beds is used to propose a constitutive closure for the effective viscosity, a nonzero parameter that is needed in the liquid momentum balance to avoid the prediction of disturbances with an infinite growth rate. Near-solitary wave profiles are predicted over a range of velocity parameters consistent with linear stability analysis. Centimeter-scale periodic disturbances are predicted in the near-pulsing regime. Preliminary estimates of average pulse properties compare well with typically reported experimental values. Initial comparison with time integration subject to periodic boundary conditions shows agreement of the liquid saturation profiles but differences in the liquid velocity profiles.


Introduction
e volume averaged two-fluid model [1,2] has been successfully applied in the literature to predict key engineering quantities (average liquid saturation, pressure drop, and flow regime transitions) that are needed for the design of packed beds with cocurrent and downward gas-liquid flow. Grosser et al. [3], who were the first to use such an approach in packed beds, were able to predict with some success the observed transition from trickle to pulse flow. In a subsequent paper, Dankworth et al. [4] provided improved closure relations for the momentum transfer terms and included a Newtonian-like stress tensor with an effective viscosity in the liquid momentum balance. ese authors confirmed that the location of the transition predicted from linear stability is insensitive to the value of the effective viscosity. However, they noted that when a zero effective viscosity is used, all the modes become unstable at once, with growth rates that increase monotonically to infinity. erefore, a nonzero effective viscosity is required to formulate a working dynamic model in the near-pulsing regime. As pointed out in a recent paper [5], the scaling argument used by Dankworth et al. to obtain a closure for the effective viscosity is based on the assumption that typical variations in the interstitial liquid velocity, a mesoscale quantity, occur over distances that are of the order of the packing diameter. In this paper, we show that it is possible to extract a closure for the effective viscosity without making any a priori choices for the characteristic length. First, we use a traveling wave approximation to translate the problem into a 2D dynamical system; then, we impose a requirement on the eigenvalues of one of the fixed points along a special path on the bifurcation diagram. e resulting closure leads to reasonable predictions of average pulse properties near the bubble-to-pulse transition.

The Model Equations
e continuity equations of the one-dimensional model may be written as follows [5]: − zε l zt + v g zε l zx + 1 − ε l zv g zx � 0 (gas). (2) e momentum balances are given by the following equations: Under typical simplifying assumptions, the momentum jump condition is given by the following equation: In equations (3) and (4), the gas-liquid interaction term is written as the product of the relative velocity and a function of liquid saturation, denoted by β(ε l ). e derivative terms with D refer to material derivatives. In equation (5), H(ε l ) is the mesoscale average value of the mean curvature of the interface. e quantities ρ g and σ in equations (4) and (5) refer to the average gas density in the column and the surface tension at the gas-liquid interface, respectively. In equation (3), the parameter A ls of the liquid-solid interaction term is given by the Ergun equation: where V sl is the superficial liquid velocity, ρ l and μ l are the density and viscosity of the liquid phase, ε is the average bed porosity, and d p is the packing size. In this work, we use the values 150 and 1.75 for the Ergun constants E 1 and E 2 . Lastly, it follows from equations (1) and (2) that where e solution of the above equations requires closure relations for the gas-liquid and liquid-solid force interaction terms (λ B and β), the momentum jump condition (H ′ ), and the effective viscosity (μ eff ). A detailed discussion of the first three closures may be found in Salgi and Balakotaiah [5]. In this work, we show that the remaining closure for the effective viscosity may be obtained from the bifurcation diagram of the dynamical problem formulated in the so-called traveling wave approximation.
To estimate the average pulse properties under given flow conditions (gas and liquid flow rates), we fix the liquid mass flux L and choose the gas mass flux G to be greater than the "transition" value predicted from our linear stability analysis [5]. We then consider a hypothetical bubbly flow (as described by our model) under these conditions. As long as we do not go far into the pulsing regime, the average properties (pressure gradient and liquid holdup) of this hypothetical bubbly flow will still be in good agreement with those of the experimental pulsing flow [5]. As discussed by Dankworth et al. [4], visible pulses may be described as solitary waves (i.e., waves where the width of the disturbance is much smaller than the wavelength). In the traveling wave approximation, which translates the model equations into a 2D dynamical problem, solitary waves are typically modeled as "homoclinic" or "heteroclinic" orbits, which consist of saddle points and the paths connecting these points [6]. In our case, the 2D dynamical problem has a single saddle point (along with another fixed point that is a focus in the relevant range of parameters). Pulsing solutions may be sought as the merging point between two paths: (i) the path along which the periodic solutions have an average liquid flux equal to the liquid flux entering the column (the "constrained path") and (ii) a line of (homoclinic) saddle connections. We find that two such paths emerge from a third special path (the "k 0 l Path" defined in Section 3.3), which provides a connection with the linear stability results of the full model. As already pointed out by Dankworth et al. [4], the merging of the constrained path with the line of saddle connections is gradual and near-solitary wave profiles ("pulsing solutions") may be generated over a range of velocity parameters. However, the average pulse properties extracted from these solutions are rather similar when the velocity parameter is restricted to values consistent with linear stability analysis. ese predictions are also found to compare favorably with typically reported experimental values for average pulse velocity, amplitude, and frequency.
When generating our "pulsing solutions," it was necessary to evaluate the functions β and H ′ at liquid saturations that lie outside the range over which these closures can be obtained. Typically, we first generate a best polynomial fit over the range where the function is given by the closure and then use the polynomial fit to extrapolate (see Figures 1(a) and 1(b) for an example). In addition to its simplicity, the justification for such fitting and extrapolation procedure lies ultimately in the reasonable agreement of the resulting predictions with experimental values of average pulse properties.

Traveling Wave Analysis
Elimination of the pressure gradients in equations (3) and (4) leads to the following equation for the variables ε l and v l : e functions c 1 , c 2 , f, g, and f 2 are de ned as follows: In equation (9), Δρ (ρ l − ρ g ) and v r (v g − v l ) is the relative velocity obtained from equation (7): In equation (10), we have assumed that C(t) in equation (7) is a constant denoted by U t . Equations (1), (9), and (10) provide the closed set of equations that we use in this work to solve for ε l and v l . e choice of the parameter U t is discussed in the next section.
In the traveling wave approximation, we de ne z x−ct (c is nonzero positive) and look for periodic solutions of the form ε l (z) and v l (z), at xed values of L and G, with G being above its "transition" value (linear instability limit). For solutions that depend only on the traveling wave coordinate z, equation (1) gives the following equation: Substitution of equation (11) into equation (9) leads to the following second-order di erential equation for the liquid saturation: In equation (12), we have e problem is reduced to a single dependent variable (the liquid saturation) which is a solution to a nonlinear second-order ordinary di erential equation (equation (12)). It is then possible to de ne an "initial value" problem by reducing this ODE to a system of rst-order equations: Equation (14) is subsequently analyzed using the language of dynamical systems (for basic de nitions, see [6])

International Journal of Chemical Engineering
with c and k l as "bifurcation parameters". It is possible to provide a detailed description of bifurcations in the (c, k l ) plane [4]. Here, we will only describe those features that are needed to (i) define a closure for the effective viscosity parameter and (ii) find the (c, k l ) coordinates of a (near) solitary wave with properties adequate to describe pulsing at the given L and G values.

Choice of the Integration Constant U t .
In general, it is not clear how to relate the parameter U t to the mass fluxes L and G in the case of periodic solutions. One possibility is to average equation (10) over a period and to require the average (total) volumetric flux to be equal to the total volumetric flux entering the actual column. is leads to the following value for U t : For wave profiles that do not depend explicitly on time, such as those we seek in the traveling wave approximation, time averaging over a period across a cross section translates easily into spatial averaging and equation (15) simply states that the total flow over a time period is equal to the flow obtained from the externally imposed L and G values. Equation (15) is used in all our traveling wave calculations. Since the value of U t in equation (15) is equal to the value corresponding to the base uniform flow, we expect local bifurcations of the base uniform flow in the traveling wave approximation to be consistent with the results of linear stability analysis. We find that this is indeed the case as described in a subsequent section (" e k 0 l Path"). However, we should note that time integration of the partial differential equations with periodic boundary conditions could not be carried out with the value of U t given in equation (15) (which fails before saturation as we move far enough from the initial linear stability solutions). Instead, we find that a value of U t based on equation (23) can yield saturated solutions.

e Constrained Path.
e first path of interest in the (c, k l ) plane is defined using the constraint on the average liquid volumetric flux: for periodic solutions of finite wavelength λ. (16) In the remainder of this paper, we will refer to the path defined in equation (16) as the "constrained path". We should note that similar paths were examined by Dankworth et al. [4]. e value of k l along this path is consistent with the externally imposed flow conditions and, consequently, the "physical pulsing solution" should belong to this path. When integration with periodic boundary conditions that could be carried out to follow the time evolution of a particular unstable linear mode, the liquid saturation profile evolved to a solution that is in good agreement with the constrained traveling wave profile corresponding to a velocity parameter equal to the phase velocity of the unstable linear mode. However, significant differences were found in the liquid velocity profiles. We also find that the constrained path merges gradually with a line of saddle connections as the value of the velocity parameter is increased (see Figures 2(a) and 2(b)). Near-solitary wave profiles are obtained for values of c above c node , which is slightly lower than the phase velocity of the zero wave number from linear stability.

e k 0
l Path. For fixed values of L and G, and with k l given by equation (16), the dynamical system described in previous paragraphs can be solved provided that the value of the velocity parameter c is known. For a given L, when G is above the "transition" value (as predicted by linear stability analysis of the base flow), a range of unstable modes emerges. is range is bounded by two "neutral" wave numbers, 0 and k N . e phase velocities of the emerging unstable modes are also bounded by the values corresponding to the two neutral wave numbers. In our model, these limiting phase velocities were assigned specific values [5]. erefore, we expect the range of relevant c values to correspond with the range defined by these two limits. For convenience, we reproduce here the velocity values assigned in these limits: In equation (18), ε l,max is the liquid saturation value corresponding to the maximum of the β profile and the quantities with subscript 0 are those of the base bubbly flow solution (for more details, see [5]). For fixed L, when G is equal to the "transition" value, the two velocities in equations (17) and (19) are equal and, in this case, c can only take this one value. As G increases beyond its transition value, the size of the interval enclosed by the two bounds increases and c can take an increasingly wider range of values. For this range of velocities to emerge from the traveling wave analysis, we define the following path: Along the k 0 l path, the base uniform state ε 0 l is a fixed point for all values of c. We find that the neutral velocities of linear stability (equations (17) and (19)) correspond to bifurcation points of the base uniform state ε 0 l along this path. As shown in Figures 2(a) and 2(b), this fixed point undergoes three bifurcations as c is varied along the k 0 l path (for definitions, see [6]). At the lower end, the value of c equal to the phase velocity of the nonzero neutral wave number (V p (k N ) de ned in equation (17)) corresponds to a Hopf bifurcation. In the gures, this value of c is denoted by c H . When c reaches the value denoted by c hmcl in the gures, the system undergoes a homoclinic bifurcation. is value of the velocity appears to correspond with the neutral velocity at the transition gas ow rate (as obtained from linear stability). At the higher end of the c interval, the value of c equal to the phase velocity of the zero neutral wave number (V p (0) de ned in equation (19)) corresponds to a node-saddle bifurcation. In the gures, this value of c is denoted by c ∩ . e focus changes into a node at a slightly lower value of c, denoted by c node in the gure. Lastly, we note that linear stability of a xed point in the traveling wave analysis is determined by the observed behavior of a perturbation as z, the "time variable" of the dynamical system, goes to plus in nity. For a xed x location (lab frame location), it is clear from the de nition of z that the growth of a perturbation when z goes to plus in nity will be the "opposite" of its growth when time goes to plus in nity (e.g., an eigenvalue with a positive real part will mean an unstable xed point when z goes to plus in nity but a stable xed point when time goes to plus in nity). For instance, over the entire range of c values de ned by the bounds in equations (17) and (19), the "bubbly ow" base state solution ε 0 l is found to be unstable (as expected) when time goes to in nity for a xed observer. erefore, the stability behavior of the base uniform ow along the k 0 l path as c increases past c H and onto c ∩ is consistent with its stability behavior in the full model as the wavenumber decreases past k N and onto 0. In addition, it is possible to visualize graphically the dependence of the eigenvalues of the ε 0 l xed point along k 0 l on those obtained from linear stability of the "full" model. In Figure 3(a), we plot the ratio of χ r , the real part of the eigenvalue from the traveling wave analysis, to the quantity −(ξ/c) (the ratio of the growth rate from linear stability in the full model to the velocity) versus the velocity c, between c H and c node . As expected, this ratio is seen to be positive and to approach 1 near the Hopf velocity. In Figure 3(b), we compare χ i , the (positive) imaginary part of the eigenvalue from traveling wave analysis, with the wavenumber k corresponding to a phase velocity c in the linear stability analysis of the full model. e gure shows close agreement between these two quantities. In both Figures 3(a) and 3(b), the e ective viscosity is 4.2 Pa·s, the value obtained from our closure as described in a subsequent section. As shown by the gures, very near the Hopf velocity (when 1/c * (ξ/kc) ≪ 1), the solution of the linearized equations in the traveling wave approximation corresponds to an exponentially growing (in the negative z direction) oscillation of wavenumber k and spatial growth rate −(ξ/c). When the nonlinear traveling wave equations are solved, the exponential growth gets saturated, which leads to a nite-amplitude solution that gradually displaces (by advection in the " ow" direction) an unstable uniform ow in a semi-in nite domain with "boundary" at z 0.
It is possible to provide a more formal comparison of linear stability in the full model with linear stability in the traveling wave approximation very near the Hopf velocity. In the full model, each phase velocity c is associated with a wavenumber k and a (positive) temporal growth rate ξ. We de ne the following dimensionless quantities: Very near the Hopf velocity, the reciprocal of the dimensionless velocity c * is much smaller than 1. We de ne the variables (x * * , t * * ) by rst rotating the (x * , t * ) axes counterclockwise by an angle whose tangent is equal to 1/c * and then scaling by the sine of the same angle. e variable t * * is equal to the negative of z * , the dimensionless z. We then write equations (1) and (9) in a dimensionless form and make the change of variables to (x * * , t * * ). If we restrict ourselves to solutions that depend only on t * * (the x * * axis corresponds to an almost constant t * value when 1/c * is very small), the resulting (dimensional) equations will match equations (11) and (12) of the traveling wave approximation.
is same transformation provides an approximate way to relate the linear modes of the full model (e ξt e ik(x−ct) ) to those of the traveling wave approximation (e χz ). We nd that e −(ξ/c)x e ξt e ik(x−ct) transforms into e χz with χ r − ξ c and χ i k.
e k 0 l path is useful for another reason. e other xed point along this path is a saddle point. e distance between the trajectories described by the eigenvalues of this point as c is varied exhibits a minimum. As described later, a closure for the e ective viscosity parameter can be obtained by requiring the location of this minimum to coincide with the value of c corresponding to the node-saddle bifurcation of the base uniform ow (c ∩ ).  which is about 10% above the transition value predicted from linear stability (G trans 0.054 kg/m 2 ·s). In Figure 2(b), the gas mass ux is 0.07 kg/m 2 ·s, about 30% above the transition value. In both cases, the constrained path is seen to start at the Hopf bifurcation on the k 0 l path and to merge gradually with a saddle connection line. is saddle connection line originates on the k 0 l path, stretches towards lower values of c, and then turns back toward higher velocities. We do not nd any periodic solutions to the right of this line. ere is a qualitative di erence in the shape of the constrained path between the two gures. For G close to the transition (Figure 2(a)), the "turning point" of the saddle connection line is located at a velocity that is lower than the Hopf velocity c H . e constrained path, which has to go around the saddle connection line, starts toward lower velocities from the Hopf bifurcation.

Two Types of Constrained
is results in two constrained solutions for every velocity at or below the Hopf velocity, indicating that at least one has to be unstable. At the higher gas ow rate (Figure 2(b)), the "turning point" of the saddle connection line takes place at a velocity that is higher than the Hopf velocity. In this case, the constrained path starts from the Hopf bifurcation on the k 0 l path towards higher c values. e qualitative di erences between Figures 2(a) and 2(b) are also observed at higher liquid ow rates. Further analysis of solutions of the same type as those in Figure 2(a) is beyond our present scope and will be addressed in future work. However, these solutions suggest that pulses formed just beyond the bubble-to-pulse transition may be unstable. In the remainder of the paper, we focus exclusively on solutions of the same type as those in Figure 2(b).  International Journal of Chemical Engineering First, we note that equation (20) is linear in c. For c values slightly above the Hopf bifurcation (see Figure 2(b)), equation (16) is found to be equivalent to equation (20). Beyond this initial range of velocities, equation (16) is still found to be linear over a short interval of c values (Figure 2(b)). We find that the predicted (finite-amplitude) ε l profile is practically the same for all values of c in this interval. For c values beyond this interval, equation (16) is no longer linear in c and the amplitude of the predicted solutions as well as their average value and wavelength steadily increase. As c approaches c node , the constrained path of equation (16) "starts to merge" with the line of saddle connections. At c node , the constrained path is within less than one percent of the line of saddle connections (Figure 2(b)). However, the consistency with the externally imposed flow rates (equation (16)) is grossly violated on the latter. At this velocity, the solutions on both lines have the same amplitude and width of disturbance but the period is about 4 times larger on the line of saddle connections. e distance between the two lines continues to decrease as c is increased beyond c node . As mentioned in Section 2, a range of near-solitary wave solutions or "pulsing solutions" is predicted for c beyond c node . Although the amplitude of these solutions does not seem to vary much as c increases, the period does increase significantly. However, if we restrict the velocity parameter to values between c node and c ∩ (consistent with linear stability results), the period will vary only by a few percents.
As an example, we show in Figures 4(a) and 4(b) the periodic profiles obtained at c � 0.59 m/s (slightly higher than c node , which is about 0.54 m/s) for the liquid saturation and liquid velocity, respectively. e liquid saturation profile has a period of about 8.5 cm (or a frequency of about 6.9 Hz). Narrow gas-rich bands representing "pulses" (as expected near the bubble-to-pulse "transition") are predicted. ese bands are about 2 cm wide on average, representing ten particle diameters in this case. e profile for the liquid velocity suggests that as a gas-rich "pulse" moves by, liquid is "pumped" from the front and "ejected" to the back of the pulse.
is qualitative picture seems reasonable for a 1D model at "mesoscale." Lastly, we ask if it is possible to integrate equations (1), (9), and (10) in time, starting from the linear stability solutions, and asymptotically (after saturation at large time) reach the finite-amplitude periodic profiles predicted from the traveling wave analysis. Another available method for generating finite-amplitude periodic solutions consists of solving equations (1), (9), and (10) with periodic boundary conditions in "boxes" that have dimensions equal to the wavelengths of the linear modes. We will refer to this method as PBC in short. e problem is initialized with the linear stability solutions. is method presents an important restriction in that the evolution from linear-like initial conditions to finite-amplitude solutions takes place at a fixed wavelength. However, in this method, no specific relationship is assumed between v l and ε l . Although both methods are based on different assumptions/restrictions, it would be interesting to compare their predictions, at least in a preliminary fashion in this work.

Some Finite-Amplitude Periodic Solutions from the PBC Method.
As discussed previously, time integration with periodic boundary conditions could not be carried out all the way up to a saturated solution when the parameter U t was calculated from equation (15). In contrast, saturated (large time) solutions can be obtained when U t is given by the following equation: With this choice of U t , and for the same system properties and flow conditions as those in Figure 2(b), we were able to generate solutions for box sizes ranging from about 2π/k N to a little over 2π/k max (this corresponds to box sizes ranging from 2.3 cm to 5 cm). Figures 5(a) and 5(b) show the saturated profiles for liquid saturation and liquid velocity, respectively, in the case of a 3.2 cm box. We find a wave velocity of 0.43 m/s, which is very slightly lower than the (linear) phase velocity of the corresponding wave number. For better visual clarity, the solutions are repeated over three periods. In Figures 6(a) and 6(b), we show corresponding profiles from traveling wave analysis for a velocity of 0.43 m/s. We see rather good agreement in the case of liquid saturation, but the liquid velocity profiles are significantly different. More work is needed to further analyze these discrepancies. For the case examined here, we were able to generate PBC solutions for velocities ranging from about 0.41 m/s, which is close to the Hopf velocity, to 0.47 m/s. Similar to the traveling wave approximation, the amplitude of the solutions increases when k decreases (or the box size increases). Initially (near k N ), the amplitude increases rapidly and then reaches a "plateau" where the rate of increase is extremely slow, similar to the previously described behavior of the traveling wave solutions over the two linear portions of the constrained path. e PBC solutions reach this "saturation" for a box size of about 2.6 cm, and the same practically saturated solutions can be generated for box sizes up to about 5 cm. Smaller amplitude solutions are obtained for box sizes between 2.3 and 2.6 cm. We were not able to generate near-"solitary wave" solutions in the PBC. However, efforts are still currently under way to explore such a possibility.

Closure Relation for the Effective Viscosity
In the near-pulsing regime (G above the "transition"), the impact of the effective viscosity parameter on the behavior of the eigenvalues of the saddle point along the k 0 l path may be described as follows. e distance between the two eigenvalues decreases as c is increased, reaches a minimum (which we will refer to as c neck ) and then increases again. Increasing the effective viscosity increases the value of c neck . We can extract a closure for the effective viscosity by requiring that c neck coincides with c ∩ . e rationale behind this requirement is that the c location of "nearest approach" of the eigenvalues may be expected to be the same location where the two actually "collide" at the transition.

International Journal of Chemical Engineering
Numerically, the e ective viscosity may be estimated easily by comparing the location (c neck ) of the minimum distance between the eigenvalues of the saddle point with c ∩ . e distance between the eigenvalues is given by the following equation: In equation (24), ε s l is the saddle point along the path k 0 l . Figure 7 gives a graphical illustration of the approach. In Figure 8(a), we show a typical variation of the e ective viscosity parameter with the gas ow rate beyond the "transition". e drastic increase in the e ective viscosity may be interpreted as an increase in the "average length scale" of disturbances from near pore scale at the transition to centimeter scale (or mesoscale) in the near-pulsing regime. is interpretation may be quanti ed approximately by de ning the length L * as follows: In equation (25), v 0 l is the average interstitial velocity for the hypothetical bubbly ow corresponding to the xed values  Figure 5: (a) Liquid saturation pro le from time integration with periodic boundary conditions (same system as in Figure 2(b)). (b) Liquid velocity pro le from time integration with periodic boundary conditions (same system as in Figure 2(b)).
of L and G. Figure 8(b) shows that L * is predicted to increase rapidly from a very low value at the "transition" to a "mesoscale" size (a few centimeters) in the near-pulsing regime. e Newtonian stress tensor in the averaged liquid momentum balance (equation (3)) represents the lumped contribution of three di erent tensors that arise formally from the averaging procedure (see, for instance, [7]). For a Newtonian liquid, only one of these tensors (the so-called bulk stress tensor) is Newtonian in form and proportional to the viscosity of the liquid. e other two tensors involve the deviation of local velocities near the interface from the averaged "bulk" velocity (the so-called extra-deformation stress tensor) and velocity co uctuations in the bulk (a Reynolds-like stress tensor), respectively. e lumping of the three terms is a simple way to provide closure with only a single parameter, the e ective viscosity. A large value of the e ective viscosity would imply that, in the pulsing regime, the contribution of the bulk stress tensor is negligible compared to the "turbulence" source term. It is reasonable to expect this last contribution to increase when we go further into the pulsing regime, as predicted by our model (Figure 8(a)). e characteristic length de ned in equation (25) provides a "visual" interpretation of the e ective viscosity parameter. As this parameter increases, so does the characteristic length, which may be interpreted as the "size" of the "turbulent" structures or pulses.

Summary and Conclusions
In this article, we have used traveling wave analysis to provide a closure for the e ective viscosity parameter, which is needed to obtain predictions in the near-pulsing regime. Preliminary results show that the closure predicts centimeter-scale periodic disturbances in this regime. Such disturbances are obtained as near-solitary wave solutions in the traveling wave approximation, over a range of velocities that is consistent with the predictions of linear stability analysis. In addition, realistic estimates are made for routinely measured average pulse properties (such as speed, amplitude, and frequency). Lastly, we predict that pulses formed just beyond the transition from bubbly ow may not be stable. Future work will focus on a more detailed comparison with available experimental data from both 0 g and   Figure 8: (a) Evolution of the e ective viscosity as the gas ow rate is increased beyond the bubble-to-pulse transition. (b) Evolution of the characteristic length L * as the gas ow rate is increased beyond the bubble-to-pulse transition.