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.
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]:
The 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 . The derivative terms with D refer to material derivatives. In equation (5), is the mesoscale average value of the mean curvature of the interface. The quantities 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 of the liquid-solid interaction term is given by the Ergun equation:where is the superficial liquid velocity, and are the density and viscosity of the liquid phase, is the average bed porosity, and 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) thatwhere
The solution of the above equations requires closure relations for the gas-liquid and liquid-solid force interaction terms ( and ), the momentum jump condition (), and the effective viscosity (). 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 “ 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 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.
(a)
(b)
3. Traveling Wave Analysis
Elimination of the pressure gradients in equations (3) and (4) leads to the following equation for the variables and :
The functions , , , , and are defined as follows:
In equation (9), and is the relative velocity obtained from equation (7):
In equation (10), we have assumed that in equation (7) is a constant denoted by . Equations (1), (9), and (10) provide the closed set of equations that we use in this work to solve for and . The choice of the parameter is discussed in the next section.
In the traveling wave approximation, we define z = x−ct ( is nonzero positive) and look for periodic solutions of the form and , 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:
Substitution of equation (11) into equation (9) leads to the following second-order differential equation for the liquid saturation:
In equation (12), we have
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:
Equation (14) is subsequently analyzed using the language of dynamical systems (for basic definitions, see [6]) with and as “bifurcation parameters”. It is possible to provide a detailed description of bifurcations in the 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 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 Constant
In general, it is not clear how to relate the parameter 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 :
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 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 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 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 based on equation (23) can yield saturated solutions.
3.2. The Constrained Path
The first path of interest in the plane is defined using the constraint on the average liquid volumetric flux:
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 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 above , which is slightly lower than the phase velocity of the zero wave number from linear stability.
(a)
(b)
3.3. The Path
For fixed values of L and G, and with given by equation (16), the dynamical system described in previous paragraphs can be solved provided that the value of the velocity parameter 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 . 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 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), 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, 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 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 path, the base uniform state is a fixed point for all values of . We find that the neutral velocities of linear stability (equations (17) and (19)) correspond to bifurcation points of the base uniform state along this path. As shown in Figures 2(a) and 2(b), this fixed point undergoes three bifurcations as is varied along the path (for definitions, see [6]). At the lower end, the value of equal to the phase velocity of the nonzero neutral wave number ( defined in equation (17)) corresponds to a Hopf bifurcation. In the figures, this value of is denoted by . When reaches the value denoted by 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 interval, the value of equal to the phase velocity of the zero neutral wave number ( defined in equation (19)) corresponds to a node-saddle bifurcation. In the figures, this value of is denoted by . The focus changes into a node at a slightly lower value of , denoted by 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 values defined by the bounds in equations (17) and (19), the “bubbly flow” base state solution 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 path as increases past and onto is consistent with its stability behavior in the full model as the wavenumber decreases past and onto 0. In addition, it is possible to visualize graphically the dependence of the eigenvalues of the fixed point along on those obtained from linear stability of the “full” model. In Figure 3(a), we plot the ratio of , the real part of the eigenvalue from the traveling wave analysis, to the quantity (the ratio of the growth rate from linear stability in the full model to the velocity) versus the velocity , between and . As expected, this ratio is seen to be positive and to approach 1 near the Hopf velocity. In Figure 3(b), we compare , the (positive) imaginary part of the eigenvalue from traveling wave analysis, with the wavenumber corresponding to a phase velocity 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 ), the solution of the linearized equations in the traveling wave approximation corresponds to an exponentially growing (in the negative z direction) oscillation of wavenumber and spatial growth rate . 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.
(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 velocity is associated with a wavenumber and a (positive) temporal growth rate . We define the following dimensionless quantities:
Very near the Hopf velocity, the reciprocal of the dimensionless velocity is much smaller than 1. We define the variables by first rotating the axes counterclockwise by an angle whose tangent is equal to and then scaling by the sine of the same angle. The variable is equal to the negative of , the dimensionless . We then write equations (1) and (9) in a dimensionless form and make the change of variables to . If we restrict ourselves to solutions that depend only on (the axis corresponds to an almost constant value when 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 () to those of the traveling wave approximation (). We find that
The 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 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 corresponding to the node-saddle bifurcation of the base uniform flow ().
3.4. Two Types of Constrained Paths
Figures 2(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 path and to merge gradually with a saddle connection line. This saddle connection line originates on the path, stretches towards lower values of , 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 . 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 path towards higher 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 . For 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 values (Figure 2(b)). We find that the predicted (finite-amplitude) profile is practically the same for all values of in this interval. For values beyond this interval, equation (16) is no longer linear in and the amplitude of the predicted solutions as well as their average value and wavelength steadily increase. As approaches , the constrained path of equation (16) “starts to merge” with the line of saddle connections. At , 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 is increased beyond . As mentioned in Section 2, a range of near-solitary wave solutions or “pulsing solutions” is predicted for beyond . Although the amplitude of these solutions does not seem to vary much as increases, the period does increase significantly. However, if we restrict the velocity parameter to values between and (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 , 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.”
(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 and . 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 parameter was calculated from equation (15). In contrast, saturated (large time) solutions can be obtained when is given by the following equation:
With this choice of , 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 to a little over (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 decreases (or the box size increases). Initially (near ), 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.
(a)
(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 path may be described as follows. The distance between the two eigenvalues decreases as is increased, reaches a minimum (which we will refer to as ) and then increases again. Increasing the effective viscosity increases the value of . We can extract a closure for the effective viscosity by requiring that coincides with . The rationale behind this requirement is that the 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 () of the minimum distance between the eigenvalues of the saddle point with . The distance between the eigenvalues is given by the following equation:
In equation (24), is the saddle point along the path . 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 as follows:
(a)
(b)
In equation (25), is the average interstitial velocity for the hypothetical bubbly flow corresponding to the fixed values of L and G. Figure 8(b) shows that 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.
Symbols
| : | Liquid saturation (dimensionless) |
| : | Density of liquid or gas phase (kg/m3) |
| : | Effective viscosity (Pa·s) |
| : | Liquid-solid interaction term (dimensionless) |
| : | Gas-liquid interaction term (Pa·s/m2) |
| : | Gas-liquid interfacial tension (N/m) |
| : | Average bed porosity (dimensionless) |
| : | Liquid viscosity (Pa·s) |
| : | Temporal growth rate (s−1) |
| : | Average mesoscale velocity in liquid or gas phase (m/s) |
| : | Average mesoscale pressure in liquid or gas phase (Pa) |
| : | Ergun parameter (Pa/m) |
| : | Mesoscale average of the gas-liquid interface mean curvature (m−1) |
| : | Packing diameter (m) |
| : | Liquid superficial velocity (m/s) |
| : | Ergun coefficients (dimensionless) |
| : | Integration constant (m/s) |
| : | Integration constant (m/s) |
| : | Velocity parameter (m/s) |
| : | Wavenumber (m−1). |
Data Availability
The data on average pulse properties used for comparison with the predictions made in this article may be found in the studies cited within the text as references [8–10].
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work was supported by a grant from the NASA Glenn Research Center (Grant NNX17AF91 G). The authors would like to thank V. Balakotaiah at UH and B. J. Motil and E. Ramé at GRC for their kind support over the last few years.