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

Authors: Paul Salgi; Sanaa Krayem
Journal: International Journal of Chemical Engineering
 (2019)

---

## Abstract

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.


---

## Body

## 1. Introduction

The 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. These 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. Therefore, 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. The resulting closure leads to reasonable predictions of average pulse properties near the bubble-to-pulse transition.


## 2. The Model Equations

The continuity equations of the one-dimensional model may be written as follows [5]:(1)∂εl∂t+vl∂εl∂x+εl∂vl∂x=0liquid,(2)−∂εl∂t+vg∂εl∂x+1−εl∂vg∂x=0gas.The momentum balances are given by the following equations:(3)ρlDlvlDt=−εl∂pl∂x+μeff∂∂xεl∂vl∂x−AlS1+λB1−εl+βεlvg−vl+ρlg=0liquid,(4)1−εlρgDgvgDt−ρlDlvlDt=−1−εl∂pg∂x−βεlvg−vl+1−εlρg−ρlggas.Under typical simplifying assumptions, the momentum jump condition is given by the following equation:(5)pg−pl=2σHεl.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. The 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. The 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 Als of the liquid-solid interaction term is given by the Ergun equation:(6)Als=E11−ε2μlε3dp2Vsl+E21−ερlε3dpVsl2,where Vsl is the superficial liquid velocity, ρl and μl are the density and viscosity of the liquid phase, ε is the average bed porosity, and dp is the packing size. In this work, we use the values 150 and 1.75 for the Ergun constants E1 and E2.Lastly, it follows from equations (1) and (2) that(7)εlvl+εgvg=Ct,where(8)εg=1−εl.The 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 fluxL 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 “kl0 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. These 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.Figure 1
(a) Sample fit and extrapolation of theβ function (air/water). (b) Sample fit and extrapolation of the H′ function (air/water).
(a)
(b)


## 3. Traveling Wave Analysis

Elimination of the pressure gradients in equations (3) and (4) leads to the following equation for the variables εl and vl:(9)ρl−Δρεl2vl∂εl∂t+ρl−Δρεl2εl∂vl∂t+c1εl,vl−fεl∂εl∂x+c2εl,vl∂vl∂x−gεl∂∂xεl∂vl∂x+f2εl,vl=0.The functionsc1, c2, f, g, and f2 are defined as follows:(9a)c1εl,vl=ρl−Δρεl2vl2−ρgεl2vr2,(9b)c2εl,vl=2ρl−Δρεl2εlvl+2ρgεl3vr,(9c)fεl=2σ1−εlεl2H′εl,(9d)gεl=μeff1−εlεl,(9e)f2εl,vl=AlS1−εlεl1+λB1−εl−βεlεlvr−ρl1+εl−ρgεl1−εlεlg.In equation (9), Δρ=ρl−ρg and vr=vg−vl is the relative velocity obtained from equation (7):(10)1−εlvr+vl=Ut.In equation (10), we have assumed that Ct in equation (7) is a constant denoted by Ut. Equations (1), (9), and (10) provide the closed set of equations that we use in this work to solve for εl and vl. The choice of the parameter Ut is discussed in the next section.In the traveling wave approximation, we definez = x−ct (c is nonzero positive) and look for periodic solutions of the form εlz and vlz, at fixed 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:(11)εlvl−c=klconstant  independent  ofz.Substitution of equation (11) into equation (9) leads to the following second-order differential equation for the liquid saturation:(12)klμeff1−εld2εldz2+f1∗εl,cdεldz−klμeff1−εlεldεldz2=f2∗εl,c.In equation (12), we have(13a)f1∗εl,c=c1εl,vl−fεl−klc2εl,vlεl2−c2ρl−Δρεl2,(13b)f2∗εl,c=−AlS1−εlεl1+λB1−εl+βεl1−εlUt−cεl−kl+ρl1+εl−ρgεl1−εlεlg.The problem is reduced to a single dependent variable (the liquid saturation) which is a solution to a nonlinear second-order ordinary differential equation (equation (12)). It is then possible to define an “initial value” problem by reducing this ODE to a system of first-order equations:(14)dUdz=FU,kl,c.Equation (14) is subsequently analyzed using the language of dynamical systems (for basic definitions, see [6]) with c and kl as “bifurcation parameters”. It is possible to provide a detailed description of bifurcations in the c,kl 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,kl coordinates of a (near) solitary wave with properties adequate to describe pulsing at the given L and G values.


### 3.1. Choice of the Integration ConstantUt

In general, it is not clear how to relate the parameterUt 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. This leads to the following value for Ut:(15)Ut=1εLρl+Gρg.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 Ut 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 (“The kl0 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 ofUt 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 Ut based on equation (23) can yield saturated solutions.


### 3.2. The Constrained Path

The first path of interest in thec,kl plane is defined using the constraint on the average liquid volumetric flux:(16)kl=Lερl−c1λ∫0λεldz,forperiodicsolutionsoffinitewavelengthλ.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]. The value of kl 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 cnode, which is slightly lower than the phase velocity of the zero wave number from linear stability.Figure 2
(a) Sample bifurcation diagram very near the “transition” predicted from linear stability (Gtrans = 0.054 kg/m2·s). The air/water system under microgravity is with dp = 2 mm and L = 15 kg/m2·s. (b) Sample bifurcation diagram farther away from the “transition” predicted from linear stability (Gtrans = 0.054 kg/m2·s). The air/water system under microgravity is with dp = 2 mm and L = 15 kg/m2·s.
(a)
(b)


### 3.3. Thekl0 Path

For fixed values ofL and G, and with kl 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. This range is bounded by two “neutral” wave numbers, 0 and kN. The 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]. Therefore, 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:(17)VpkN=vl0ac1−εl02lower bound,(18)ac=vlεl,max1−εl,max2vgεl,max+1−εl,max/εl,maxvrεl,max,(19)Vp0=vg0+1−εl0εl0vr0upper bound.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:(20)kl0=Lερl−cεl0.Along thekl0 path, the base uniform state εl0 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 εl0 along this path. As shown in Figures 2(a) and 2(b), this fixed point undergoes three bifurcations as c is varied along the kl0 path (for definitions, see [6]). At the lower end, the value of c equal to the phase velocity of the nonzero neutral wave number (VpkN defined in equation (17)) corresponds to a Hopf bifurcation. In the figures, this value of c is denoted by cH. When c reaches the value denoted by chmcl in the figures, the system undergoes a homoclinic bifurcation. This value of the velocity appears to correspond with the neutral velocity at the transition gas flow 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 (Vp0 defined in equation (19)) corresponds to a node-saddle bifurcation. In the figures, this value of c is denoted by c∩. The focus changes into a node at a slightly lower value of c, denoted by cnode in the figure. Lastly, we note that linear stability of a fixed 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 infinity. For a fixed x location (lab frame location), it is clear from the definition of z that the growth of a perturbation when z goes to plus infinity will be the “opposite” of its growth when time goes to plus infinity (e.g., an eigenvalue with a positive real part will mean an unstable fixed point when z goes to plus infinity but a stable fixed point when time goes to plus infinity). For instance, over the entire range of c values defined by the bounds in equations (17) and (19), the “bubbly flow” base state solution εl0 is found to be unstable (as expected) when time goes to infinity for a fixed observer. Therefore, the stability behavior of the base uniform flow along the kl0 path as c increases past cH and onto c∩ is consistent with its stability behavior in the full model as the wavenumber decreases past kN and onto 0. In addition, it is possible to visualize graphically the dependence of the eigenvalues of the εl0 fixed point along kl0 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 cH and cnode. 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. The figure shows close agreement between these two quantities. In both Figures 3(a) and 3(b), the effective viscosity is 4.2 Pa·s, the value obtained from our closure as described in a subsequent section. As shown by the figures, 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 finite-amplitude solution that gradually displaces (by advection in the “flow” direction) an unstable uniform flow in a semi-infinite domain with “boundary” at z = 0.Figure 3
(a) Real part of the eigenvalue of theεl0 fixed point along the kl0 path as a function of the linear growth rate and phase velocity of the base uniform flow in the full model, for the same system as in Figure 2(b). (b) Imaginary part of the eigenvalue of the εl0 fixed point along the kl0 path as a function of the wavenumber with phase velocity c for the base uniform flow in the full model, for the same system as in Figure 2(b).
(a)
(b)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 velocityc is associated with a wavenumber k and a (positive) temporal growth rate ξ. We define the following dimensionless quantities:(21)t∗=tξ,x∗=kx,c∗=kcξ.Very near the Hopf velocity, the reciprocal of the dimensionless velocityc∗ is much smaller than 1. We define the variables x∗∗,t∗∗ by first 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. The 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. This same transformation provides an approximate way to relate the linear modes of the full model (eξteikx−ct) to those of the traveling wave approximation (eχz). We find that(22)e−ξ/cxeξteikx−cttransformsintoeχzwithχr=−ξcandχi=k.Thekl0 path is useful for another reason. The other fixed point along this path is a saddle point. The 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 effective 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 flow (c∩).


### 3.4. Two Types of Constrained Paths

Figures2(a) and 2(b) show typical results from the traveling wave analysis. The liquid mass flux L is 15 kg/m2·s under 0 g conditions. Figure 2(a) corresponds to a gas mass flux of 0.06 kg/m2·s, which is about 10% above the transition value predicted from linear stability (Gtrans = 0.054 kg/m2·s). In Figure 2(b), the gas mass flux is 0.07 kg/m2·s, about 30% above the transition value. In both cases, the constrained path is seen to start at the Hopf bifurcation on the kl0 path and to merge gradually with a saddle connection line. This saddle connection line originates on the kl0 path, stretches towards lower values of c, and then turns back toward higher velocities. We do not find any periodic solutions to the right of this line. There is a qualitative difference in the shape of the constrained path between the two figures. 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 cH. The constrained path, which has to go around the saddle connection line, starts toward lower velocities from the Hopf bifurcation. This 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 flow 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 kl0 path towards higher c values. The qualitative differences between Figures 2(a) and 2(b) are also observed at higher liquid flow 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).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 cnode, the constrained path of equation (16) “starts to merge” with the line of saddle connections. At cnode, 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. The distance between the two lines continues to decrease as c is increased beyond cnode. As mentioned in Section 2, a range of near-solitary wave solutions or “pulsing solutions” is predicted for c beyond cnode. 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 cnode and c∩ (consistent with linear stability results), the period will vary only by a few percents.As an example, we show in Figures4(a) and 4(b) the periodic profiles obtained at c = 0.59 m/s (slightly higher than cnode, which is about 0.54 m/s) for the liquid saturation and liquid velocity, respectively. The 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. These bands are about 2 cm wide on average, representing ten particle diameters in this case. The 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. This qualitative picture seems reasonable for a 1D model at “mesoscale.”Figure 4
(a) A pulsing solution (c = 0.59 m/s) for the liquid saturation from traveling wave analysis (same system as in Figure 2(b)). (b) A pulsing solution (c = 0.59 m/s) for the liquid velocity from traveling wave analysis (same system as in Figure 2(b)).
(a)
(b)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. The problem is initialized with the linear stability solutions. This 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 vl 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.


### 3.5. 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 parameterUt was calculated from equation (15). In contrast, saturated (large time) solutions can be obtained when Ut is given by the following equation:(23)Ut−vl0=1εLρl+Gρg.With this choice ofUt, 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π/kN to a little over 2π/kmax (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 kN), 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. The 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.Figure 5
(a) Liquid saturation profile from time integration with periodic boundary conditions (same system as in Figure2(b)). (b) Liquid velocity profile from time integration with periodic boundary conditions (same system as in Figure 2(b)).
(a)
(b)Figure 6
(a) Liquid saturation profile atc = 0.43 m/s from traveling wave analysis (same system as in Figure 2(b)). (b) Liquid velocity profile at c = 0.43 m/s from traveling wave analysis (same system as in Figure 2(b)).
(a)
(b)


## 3.1. Choice of the Integration ConstantUt

In general, it is not clear how to relate the parameterUt 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. This leads to the following value for Ut:(15)Ut=1εLρl+Gρg.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 Ut 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 (“The kl0 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 ofUt 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 Ut based on equation (23) can yield saturated solutions.


## 3.2. The Constrained Path

The first path of interest in thec,kl plane is defined using the constraint on the average liquid volumetric flux:(16)kl=Lερl−c1λ∫0λεldz,forperiodicsolutionsoffinitewavelengthλ.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]. The value of kl 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 cnode, which is slightly lower than the phase velocity of the zero wave number from linear stability.Figure 2
(a) Sample bifurcation diagram very near the “transition” predicted from linear stability (Gtrans = 0.054 kg/m2·s). The air/water system under microgravity is with dp = 2 mm and L = 15 kg/m2·s. (b) Sample bifurcation diagram farther away from the “transition” predicted from linear stability (Gtrans = 0.054 kg/m2·s). The air/water system under microgravity is with dp = 2 mm and L = 15 kg/m2·s.
(a)
(b)


## 3.3. Thekl0 Path

For fixed values ofL and G, and with kl 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. This range is bounded by two “neutral” wave numbers, 0 and kN. The 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]. Therefore, 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:(17)VpkN=vl0ac1−εl02lower bound,(18)ac=vlεl,max1−εl,max2vgεl,max+1−εl,max/εl,maxvrεl,max,(19)Vp0=vg0+1−εl0εl0vr0upper bound.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:(20)kl0=Lερl−cεl0.Along thekl0 path, the base uniform state εl0 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 εl0 along this path. As shown in Figures 2(a) and 2(b), this fixed point undergoes three bifurcations as c is varied along the kl0 path (for definitions, see [6]). At the lower end, the value of c equal to the phase velocity of the nonzero neutral wave number (VpkN defined in equation (17)) corresponds to a Hopf bifurcation. In the figures, this value of c is denoted by cH. When c reaches the value denoted by chmcl in the figures, the system undergoes a homoclinic bifurcation. This value of the velocity appears to correspond with the neutral velocity at the transition gas flow 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 (Vp0 defined in equation (19)) corresponds to a node-saddle bifurcation. In the figures, this value of c is denoted by c∩. The focus changes into a node at a slightly lower value of c, denoted by cnode in the figure. Lastly, we note that linear stability of a fixed 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 infinity. For a fixed x location (lab frame location), it is clear from the definition of z that the growth of a perturbation when z goes to plus infinity will be the “opposite” of its growth when time goes to plus infinity (e.g., an eigenvalue with a positive real part will mean an unstable fixed point when z goes to plus infinity but a stable fixed point when time goes to plus infinity). For instance, over the entire range of c values defined by the bounds in equations (17) and (19), the “bubbly flow” base state solution εl0 is found to be unstable (as expected) when time goes to infinity for a fixed observer. Therefore, the stability behavior of the base uniform flow along the kl0 path as c increases past cH and onto c∩ is consistent with its stability behavior in the full model as the wavenumber decreases past kN and onto 0. In addition, it is possible to visualize graphically the dependence of the eigenvalues of the εl0 fixed point along kl0 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 cH and cnode. 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. The figure shows close agreement between these two quantities. In both Figures 3(a) and 3(b), the effective viscosity is 4.2 Pa·s, the value obtained from our closure as described in a subsequent section. As shown by the figures, 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 finite-amplitude solution that gradually displaces (by advection in the “flow” direction) an unstable uniform flow in a semi-infinite domain with “boundary” at z = 0.Figure 3
(a) Real part of the eigenvalue of theεl0 fixed point along the kl0 path as a function of the linear growth rate and phase velocity of the base uniform flow in the full model, for the same system as in Figure 2(b). (b) Imaginary part of the eigenvalue of the εl0 fixed point along the kl0 path as a function of the wavenumber with phase velocity c for the base uniform flow in the full model, for the same system as in Figure 2(b).
(a)
(b)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 velocityc is associated with a wavenumber k and a (positive) temporal growth rate ξ. We define the following dimensionless quantities:(21)t∗=tξ,x∗=kx,c∗=kcξ.Very near the Hopf velocity, the reciprocal of the dimensionless velocityc∗ is much smaller than 1. We define the variables x∗∗,t∗∗ by first 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. The 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. This same transformation provides an approximate way to relate the linear modes of the full model (eξteikx−ct) to those of the traveling wave approximation (eχz). We find that(22)e−ξ/cxeξteikx−cttransformsintoeχzwithχr=−ξcandχi=k.Thekl0 path is useful for another reason. The other fixed point along this path is a saddle point. The 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 effective 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 flow (c∩).


## 3.4. Two Types of Constrained Paths

Figures2(a) and 2(b) show typical results from the traveling wave analysis. The liquid mass flux L is 15 kg/m2·s under 0 g conditions. Figure 2(a) corresponds to a gas mass flux of 0.06 kg/m2·s, which is about 10% above the transition value predicted from linear stability (Gtrans = 0.054 kg/m2·s). In Figure 2(b), the gas mass flux is 0.07 kg/m2·s, about 30% above the transition value. In both cases, the constrained path is seen to start at the Hopf bifurcation on the kl0 path and to merge gradually with a saddle connection line. This saddle connection line originates on the kl0 path, stretches towards lower values of c, and then turns back toward higher velocities. We do not find any periodic solutions to the right of this line. There is a qualitative difference in the shape of the constrained path between the two figures. 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 cH. The constrained path, which has to go around the saddle connection line, starts toward lower velocities from the Hopf bifurcation. This 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 flow 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 kl0 path towards higher c values. The qualitative differences between Figures 2(a) and 2(b) are also observed at higher liquid flow 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).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 cnode, the constrained path of equation (16) “starts to merge” with the line of saddle connections. At cnode, 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. The distance between the two lines continues to decrease as c is increased beyond cnode. As mentioned in Section 2, a range of near-solitary wave solutions or “pulsing solutions” is predicted for c beyond cnode. 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 cnode and c∩ (consistent with linear stability results), the period will vary only by a few percents.As an example, we show in Figures4(a) and 4(b) the periodic profiles obtained at c = 0.59 m/s (slightly higher than cnode, which is about 0.54 m/s) for the liquid saturation and liquid velocity, respectively. The 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. These bands are about 2 cm wide on average, representing ten particle diameters in this case. The 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. This qualitative picture seems reasonable for a 1D model at “mesoscale.”Figure 4
(a) A pulsing solution (c = 0.59 m/s) for the liquid saturation from traveling wave analysis (same system as in Figure 2(b)). (b) A pulsing solution (c = 0.59 m/s) for the liquid velocity from traveling wave analysis (same system as in Figure 2(b)).
(a)
(b)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. The problem is initialized with the linear stability solutions. This 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 vl 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.


## 3.5. 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 parameterUt was calculated from equation (15). In contrast, saturated (large time) solutions can be obtained when Ut is given by the following equation:(23)Ut−vl0=1εLρl+Gρg.With this choice ofUt, 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π/kN to a little over 2π/kmax (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 kN), 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. The 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.Figure 5
(a) Liquid saturation profile from time integration with periodic boundary conditions (same system as in Figure2(b)). (b) Liquid velocity profile from time integration with periodic boundary conditions (same system as in Figure 2(b)).
(a)
(b)Figure 6
(a) Liquid saturation profile atc = 0.43 m/s from traveling wave analysis (same system as in Figure 2(b)). (b) Liquid velocity profile at c = 0.43 m/s from traveling wave analysis (same system as in Figure 2(b)).
(a)
(b)


## 4. 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 kl0 path may be described as follows. The distance between the two eigenvalues decreases as c is increased, reaches a minimum (which we will refer to as cneck) and then increases again. Increasing the effective viscosity increases the value of cneck. We can extract a closure for the effective viscosity by requiring that cneck coincides with c∩. The 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.Numerically, the effective viscosity may be estimated easily by comparing the location (cneck) of the minimum distance between the eigenvalues of the saddle point with c∩. The distance between the eigenvalues is given by the following equation:(24)dλ1,λ2=−f1∗2+4kl0μeff1−εls∂f2∗/∂εlεlskl0μeff1−εls.In equation (24), εls is the saddle point along the path kl0. Figure 7 gives a graphical illustration of the approach. In Figure 8(a), we show a typical variation of the effective viscosity parameter with the gas flow rate beyond the “transition”. The drastic increase in the effective 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. This interpretation may be quantified approximately by defining the length L∗ as follows:(25)L∗=ρlvl0μeff.Figure 7
Criterion for choosing the effective viscosity (same system as in Figure2(b)).Figure 8
(a) Evolution of the effective viscosity as the gas flow rate is increased beyond the bubble-to-pulse transition. (b) Evolution of the characteristic lengthL∗ as the gas flow rate is increased beyond the bubble-to-pulse transition.
(a)
(b)In equation (25),vl0 is the average interstitial velocity for the hypothetical bubbly flow corresponding to the fixed values 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.The Newtonian stress tensor in the averaged liquid momentum balance (equation (3)) represents the lumped contribution of three different 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. The 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 cofluctuations in the bulk (a Reynolds-like stress tensor), respectively. The lumping of the three terms is a simple way to provide closure with only a single parameter, the effective viscosity. A large value of the effective 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)). The characteristic length defined in equation (25) provides a “visual” interpretation of the effective viscosity parameter. As this parameter increases, so does the characteristic length, which may be interpreted as the “size” of the “turbulent” structures or pulses.


## 5. Summary and Conclusions

In this article, we have used traveling wave analysis to provide a closure for the effective 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 flow may not be stable. Future work will focus on a more detailed comparison with available experimental data from both 0 g and 1 g down-flow data [8–10]. We will also continue exploring the feasibility of time-dependent simulations in the near-pulsing regime. This will require additional emphasis on the stability behavior of the periodic solutions predicted by the model, including those obtained just beyond the transition from bubbly flow.



---

*Source: 1635265-2019-01-08.xml*