Next Article in Journal
An Optimized-Parameter Spectral Clustering Approach to Coherent Structure Detection in Geophysical Flows
Next Article in Special Issue
Numerical Prediction of Turbulent Spray Flame Characteristics Using the Filtered Eulerian Stochastic Field Approach Coupled to Tabulated Chemistry
Previous Article in Journal
Comparison of Rheological Behaviour of Bio-Based and Synthetic Epoxy Resins for Making Ecocomposites
Previous Article in Special Issue
Effects of Mean Inflow Velocity and Droplet Diameter on the Propagation of Turbulent V-Shaped Flames in Droplet-Laden Mixtures
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Turbulent Bubble-Laden Channel Flow of Power-Law Fluids: A Direct Numerical Simulation Study

1
Institute of Applied Mathematics and Scientific Computing, Bundeswehr University Munich, Werner-Heisenberg-Weg 39, 85577 Neubiberg, Germany
2
Multiscale Modeling and Simulation, Faculty EEMCS, University of Twente, P.O Box 217, 7500 AE Enschede, The Netherlands
*
Author to whom correspondence should be addressed.
Submission received: 25 November 2020 / Revised: 16 December 2020 / Accepted: 18 December 2020 / Published: 12 January 2021
(This article belongs to the Special Issue Modelling of Reactive and Non-reactive Multiphase Flows)

Abstract

:
The influence of non-Newtonian fluid behavior on the flow statistics of turbulent bubble-laden downflow in a vertical channel is investigated. A Direct Numerical Simulation (DNS) study is conducted for power-law fluids with power-law indexes of 0.7 (shear-thinning), 1 (Newtonian) and 1.3 (shear-thickening) in the liquid phase at a gas volume fraction of 6%. The flow is driven downward by a constant volumetric flow rate corresponding to a friction Reynolds number of R e τ 127.3 . The Eötvös number is varied between E o = 0.3125 and E o = 3.75 in order to investigate the influence of quasi-spherical as well as wobbling bubbles and thus the interplay of the bubble deformability with the power-law behavior of the liquid bulk. The resulting first- and second-order fluid statistics, i.e., the gas fraction, mean velocity and velocity fluctuation profiles across the channel, show clear trends in reply to varying power-law indexes. In addition, it was observed that the bubble oscillations increase with decreasing power-law index. In the channel core, the bubbles significantly increase the dissipation rate, which, in contrast to its behavior at the wall, shows similar orders of magnitude for all power-law indexes.

1. Introduction

Turbulent bubble-laden channel flow plays an important role in a wide range of industrial applications, such as chemical reactors in the process industry, heat exchangers in power plants or boiling water reactors in the nuclear industry. Although many processes of practical interest can be considered as two-phase flows with Newtonian behavior in both phases, there is a large number of relevant applications where the continuous liquid phase exhibits non-Newtonian flow characteristics. Numerous examples can be found in the biomedical and biochemical, the metallurgical as well as the food-processing industries [1]. Common applications with bubbles moving in a non-Newtonian continuous phase are the degassing and devolatilization of polymer melts [2], the gas-assisted displacement of liquids in channels [3,4,5,6] as well as the utilization of bubble columns [7,8,9,10]. An extensive list of examples and previous investigations is provided by Chhabra [1], including some analytical, experimental and numerical studies of rise velocities and detailed velocity profiles for single bubbles and bubble chains in power-law fluids [11,12]. Many applications, such as the operation of bubble column reactors and purification processes using purge gases, are based on the heat and mass exchange between both phases. For these processes, it is beneficial if the bubbles make the maximum possible contact with the surrounding liquid. The bubble distribution and interface area are therefore key factors influencing the process efficiency. However, turbulence and the flow velocities also have a significant impact on the contact duration during which both phases can react or exchange mass and heat. Moreover, in processes with unwanted bubble generation, such as oil production, an accurate prediction of bubbly flow can help to improve the efficiency. It is therefore of great interest to analyze the dependency of the flow characteristics on the rheological behavior of the continuous liquid phase.
Since experimental investigations of bubbly flows remain challenging [13] and are often limited to dilute void fractions [14], DNS plays a major role in improving the fundamental understanding of the interaction between the continuous liquid phase and the dispersed bubbles. Due to the increased amount of computational resources available by now, DNS has become feasible for reduced-complexity two-phase flows at moderate Reynolds numbers. The pioneering numerical studies of Lu and Tryggvason [15,16,17] investigated the influence of different parameters on flows with up to 72 bubbles in the so-called minimum channel (see Section 3.2). While the computational cost and limited scalability restricted the applicability of DNS to rather simple setups for many years, recently developed massively parallel bubbly-flow solvers, such as the one used for this study [18,19], can efficiently handle up to O 10 4 fully deformable bubbles, while resolving all scales in DNS detail. Highly scalable DNS solvers for bubbly flow, which can simulate a large number of bubbles embedded in full-scale channels, offer the possibility to perform advanced evaluations, such as detailed investigations of flow topologies in dispersed bubbly flows [20,21].
Compared to the large number of studies investigating single-phase channel flow of Newtonian fluids, the effects of power-law based viscosity have rarely been investigated in this setup. Ohta and Miyashita [22] performed both single-phase DNS and Large Eddy Simulation (LES) to investigate wall turbulence for a channel setup similar to the one used in this study. To mimic the viscosity of general polymer solutions, both power-law as well as Casson fluids were used. The power-law indexes used in this study were 0.85 for the shear-thinning, and 1.15 for the shear-thickening fluid. The LES sub-grid scale stress was predicted by applying an extended Smagorinsky model.
This study extends the numerical investigation of turbulent channel flow of power-law fluids to two-phase downflow scenarios by including either quasi-spherical ( E o = 0.3125 ) or wobbling ( E o = 3.75 ) bubbles. The number of bubbles is chosen such that the total gas fraction is at a technically relevant level of 6%. The shear-thinning and shear-thickening fluids are mimicked using power-law indexes of 0.7 and 1.3, respectively. The results are compared to those of single-phase simulations with identical power-law indexes as well as corresponding bubbly flow setups with a Newtonian liquid phase. This allows for studying both the effect the presence of a second phase exerts on the flow characteristics of the non-Newtonian continuous phase, as well as the influence of the rheological behavior in the continuous phase on the bubbles. For this purpose, the resulting viscosity profiles, average wall shear stresses, and first and second order fluid statistics are compared in the statistically steady state. In addition, the wall-normal profiles of the dissipation rate are provided for all investigated cases. Finally, the bubble deformation is quantified by investigating the average bubble surface area and its fluctuation.
The mathematical formulation of the power-law model determining the viscosity in the non-Newtonian continuous phase is provided in Section 2, the numerical method and the flow configuration are described in Section 3, and the results are presented and discussed in Section 4. Section 5 summarizes the findings and draws conclusions.

2. Mathematical Description of the Power-Law Model

In this study, the non-Newtonian continuous liquid phase obeys a power-law model for viscosity. The viscous stress tensor τ i j is defined as [23]
τ i j = μ a e i j = K γ ˙ n 1 e i j
with the apparent viscosity μ a , the consistency index K, the shear rate γ ˙ = e i j e i j / 2 , the rate of strain tensor e i j = u i / x j + u j / x i and the power-law index n. The apparent viscosity of the liquid phase, which is used to implement the power-law model in the solver, can therefore be expressed as
μ a = K γ ˙ n 1 .
For shear-thinning fluids ( n < 1 ), the apparent viscosity μ a decreases with increasing shear rate, while the opposite effect occurs for shear-thickening fluids ( n > 1 ). In this study, μ a is not limited by a minimum or a maximum value. For Newtonian fluids ( n = 1 ), the shear rate does not influence the apparent viscosity, which leads to μ a = μ = K . Figure 1 illustrates the dependency of the normalized apparent viscosity μ a / μ on the shear rate for the power-law indexes 0.7, 1, and 1.3, which are used in this study.
With only one free parameter, the power-law model provides the simplest possible relation between shear stress and shear rate for non-Newtonian fluids. Although the model has some limitations and does not represent all effects present in many real-world non-Newtonian fluids perfectly [1], its simplicity allows a straightforward analysis regarding the influence of shear rate dependent viscosity on the flow characteristics. The model implementation has been validated by comparing the wall-normal profile of the stream-wise velocity component resulting from a laminar single-phase test case to an analytical solution [24] for laminar two-dimensional channel flow of power-law fluids.

3. Numerical Method and Flow Configuration

3.1. Governing Equations and Numerical Framework

The mathematical model implemented in the highly scalable (hybrid openMP-MPI) open-source code “TBFsolver” [18] uses the one-fluid formulation [25] of the governing equations. It describes the entire flow field with a single set of equations and uses a marker function to account for interfacial terms and differences in material properties, namely density and viscosity. Using the usual notation (velocity components u i , pressure p, density ρ , dynamic viscosity μ , components of the gravitational acceleration g i ), the incompressible Navier–Stokes and continuity equation are given as
ρ u i t + u i u j x j = p x i + ρ ρ 0 g i + x j μ u i x j + u j x i + σ n i κ δ S ,
u i x i = 0 .
The weight of the gas-liquid mixture is balanced by adding the additional force ρ 0 g i , with the average density of the domain volume ρ 0 = 1 V V ρ d V . The surface tension force is computed by the Continuous-Surface-Force (CSF) approach [26], using the surface tension coefficient σ , the interface normal n i = f x i / | f | and the interface curvature κ = n i / x i . The latter is accurately calculated with the Generalized Height Function method [27]. While the surface tension force is singular in theory, the interface indicator function δ S is numerically approximated as δ S = | f | , where f is the marker function determining the local volume fraction of the gas phase in the context of the geometrical Volume of Fluid (VOF) method [28]. The local density and dynamic viscosity are derived from the value of the marker function:
ρ = f · ρ g + 1 f · ρ l
μ = f · μ g + 1 f · μ l
The subscripts l and g denote the liquid and gas phase, respectively. For the cases with a non-Newtonian liquid phase, μ l in Equation (6) is calculated using the power-law model, see Equation (2).
In order to handle collisions between the bubbles and avoid the unphysical merging of bubbles in grid cells containing a part of multiple bubbles, a phenomenon referred to as numerical coalescence, the so-called multiple-marker formulation [29] is used. For each bubble b = 1 , , N present in the domain, a respective individual marker function f b is advected:
f b t + u i f b x i = 0
A geometrical reconstruction algorithm [30] performing the transport of fb ensures the sharpness of the volume fraction field, and its value on the global grid is computed from the single markers f b as f = max f b . In addition, the total surface tension contribution of multiple interface segments present in one computational cell is determined as the sum of the individual surface tension terms [18].
The Poisson equation for pressure in the framework of the projection method is transformed into a constant–coefficient Poisson equation, which allows the utilization of fast Poisson solvers, as originally shown by Dodd and Ferrante [31]. In particular, an improved projection method, designed to retain accurate numerical solutions at high density ratios, is used as detailed in [32].
Time integration is performed with a fully-explicit second-order Adams–Bashforth (AB2) method. The convective term in the momentum equation is discretized with the third-order “Quadratic Upstream Interpolation for Convective Kinematics” (QUICK) scheme, and central finite differences are used for the diffusive terms. The velocity components are arranged on a staggered grid, while the pressure and the volume fraction, as well as the density and viscosity resulting from it, are computed in the cell centers in the framework of the finite-volume method.

3.2. Case Definition

The setup investigated in this study, shown in Figure 2, is the so-called minimum channel [33], a vertical channel of size L x = π H , L y = 2 H and L z = π / 2 H . The channel half-width H is taken to be unity, and x, y, z are the stream-wise, the wall normal and the span-wise direction, respectively. A downflow configuration is investigated, with the bubbles rising relative to the surrounding liquid. In main flow direction x as well as in the second homogeneous direction z, periodic boundary conditions are imposed, whereas the domain is bounded by no-slip walls in y-direction. The flow is driven downward by a constant volumetric flow rate V ˙ / L z = y = 0 2 H u y d y , which is identical for all setups investigated in this study. The prescribed constant flow rate is the average steady-state flow rate of a single-phase channel flow simulation with a Newtonian fluid at a friction Reynolds number of R e τ = τ w / ρ l H / ν l = 127.3 , where τ w is the average wall shear stress and ν l is the kinematic viscosity of the liquid phase.
The density of the liquid phase ρ l is taken to be unity, and the density ratio of both phases is ρ l / ρ g = 10 . For the Newtonian fluid ( n = 1 and thus μ l = K ), the dynamic viscosity ratio is specified as μ l / μ g = 1 . This choice of the dynamic viscosity ratio is mainly made to avoid further complexity in addition to the variable viscosity of the liquid phase. The bulk Reynolds number, based on the total channel width 2 H , is R e = 3790 .
While the three single-phase setups, one for each power-law index, are primarily considered for comparison purposes, the main interest lies in the six two-phase setups. These are investigated at a total gas volume fraction of 6%, which is achieved by including 72 bubbles with an identical initial diameter of d b = 0.25 H . The simulation is initialized by inserting them into the velocity field of a fully developed single-phase channel flow using a uniformly distributed (x,y,z)-array of (6 × 3 × 4). To study the interplay of bubble deformability and power-law behavior of the liquid bulk, two different Eötvös numbers, E o = ρ l g d b 2 / σ , are considered. For the three setups with quasi-spherical bubbles, the surface tension coefficient is set in a way that the Eötvös number is E o = 0.3125 . Furthermore, the gravitational acceleration g is chosen such that, for the two-phase setups with the Newtonian fluid, the Galilei number G a = ρ l g d b d b / μ l has a value of 118.59 . This completes the set of material properties, which is selected to make the setup comparable to the pioneering study of Lu and Tryggvason [15]. Apart from the fact that the setups in their work are controlled by a constant pressure gradient instead of a constant volumetric flow rate, their case with the highest gas fraction is identical to one of the cases presented in this work ( n = 1 and E o = 0.3125 ). The higher Eötvös number of E o = 3.75 is achieved by decreasing the surface tension coefficient. Together with the other dimensionless numbers provided, this fully defines the investigated setup. For completeness, the values of the single parameters are given in Table A1 (Appendix A).
A joint analysis of E o and the bubble Reynolds number R e b allows an estimation of the bubble deformation. The well-known Grace diagram [34], shown in Figure 3, divides the E o - R e b -space into different regimes representing a multitude of bubble shapes ranging from simple spheres to complex, skirted geometries. In this context, R e b is defined as
R e b = u r e l d b ν a l = u l u g d b ν a l .
For the characterization of the bubble dynamics, the phase-conditional averaging · l , g is limited to the bulk part of the channel 0.5 < y / H < 1.5 . The wall boundary layers are excluded, as it is observed that the bubbles stay away from the walls (see Figure 6). To account for the variable viscosity of the non-Newtonian fluids and ensure consistency between the definitions of the relative velocity in stream-wise direction and the power-law based viscosity, the average apparent kinematic viscosity ν a l of the liquid phase is used in Equation (8). However, as will be shown later, the differences in R e b are mainly due to the different relative velocities for this setup, while ν a l has a negligible influence on this parameter. The resulting values of R e b are provided in Table 1, which gives an overview of all nine cases investigated in this study. All other conditions remaining identical, R e b increases for the shear-thinning setup while the opposite effect occurs for the shear-thickening case. As R e b just deviates slightly from the Newtonian case for n = 0.7 and n = 1.3 , only the two setups assuming n = 1 are marked in Figure 3.
All setups investigated in this study are computed using a uniform grid of 256 × 208 × 128 cells. Lu and Tryggvason [15] used 192 × 160 × 96 computational points, uniformly spaced in the homogeneous directions, but unevenly spaced and clustered near the wall in the wall-normal direction. Moreover, it was reported that a finer grid of 320 × 256 × 160 produced essentially the same statistical properties as the coarser grid. One of the setups simulated in this study (“tph10_loEo”) is identical to the setup investigated in [15], apart from the fact that the flow is driven downward by a constant volumetric flow rate instead of a constant pressure gradient. There are several reasons why a higher resolution of 256 × 208 × 128 cells is convenient for the cases investigated here. First of all, as can be expected from Figure 3, the bubble deformations increase significantly for E o = 3.75 . A finer mesh more accurately captures bubbles exhibiting more complex shapes and oscillations, i.e., ellipsoidal or wobbling instead of quasi-spherical. In [17,35], a mesh of 256 × 192 × 128 was used to simulate turbulent bubbly upflow for different Eötvös numbers, including E o = 3.75 . Moreover, the constant volumetric flow rate for the two-phase setups is taken to be that of the corresponding statistically converged single-phase flow at R e τ = 127.3 . As shown in [15], the volumetric flow rate decreases with increasing gas fraction when a constant pressure gradient is prescribed. Since the flow rate for the two-phase setups is determined using the initial conditions resulting from a single-phase setup, and kept constant, the observed flow rate is approximately 35% higher as the one reported in [15], which results from prescribing a constant pressure gradient. As a consequence, the bulk Reynolds number R e = 3790 defined above also holds true for the two-phase setups in this study. Finally, as will be shown later, the wall-normal profile of the stream-wise velocity component close to the walls is steepest for the shear-thinning fluid, i.e., u / y at the walls is, on average, observed to be the highest for n = 0.7 and the lowest for n = 1.3 . In order to accurately resolve the boundary layers at the walls in the shear-thinning case, it is therefore also advisable to increase the number of grid points from 192 × 160 × 96 to 256 × 208 × 128 6.8 × 10 6 . Using the latter resolution, the bubbles are resolved with approximately 20 cells in x- and z-directions, and 26 cells in the y-direction. The y + values of the first cell center at the wall, determined from y + = y τ w / ρ l / ν a with the apparent kinematic viscosity ν a of the liquid phase, can be found in Table 1. For all investigated setups, the values are close to one. They are determined using the viscosity in the first cell at the wall. For this purpose, the viscosity fields have been averaged over the time as well as over the homogeneous directions x and z, see Section 4.1. The variable time step is controlled by setting C F L = 0.2 for all setups.

4. Results and Discussion

This study aims to investigate both the effect the bubbles exert on the flow characteristics of the surrounding flow, as well as the influence the power-law behavior in the continuous phase has on the dispersed phase. To gain the respective insights, a number of quantities is evaluated after reaching the statistically steady state. Statistics are collected for a period corresponding to approximately 180 flow-through times ( 4.24 < t / T < 42.40 ) of the selected channel, starting after approximately 20 flow-through times ( t / T = 4.24 ) to exclude the initial transient. In this context, the reference time T is defined as T = H / u τ i n i t , with the friction velocity u τ i n i t = τ w / ρ l of the initialization setup.
All statistics are computed using the definitions introduced in [18,19]. This allows for either collecting statistics for the whole domain, or for the liquid and the gas phase separately. The respective definitions are
q c = c · q ¯ c ¯
( q 2 c ) 1 2 = c q 2 ¯ c ¯ c q ¯ c ¯ 2
where q is the quantity of interest, q its fluctuation and · ¯ denotes averaging over the time and the homogeneous directions x and z. Depending on which phase the statistics are conditioned on, c is either 1 (whole domain), f (only gas phase) or 1 f (only liquid phase). Unless indicated otherwise, the statistics are calculated for the whole domain (thus c = 1 ) and averaged over the two halves of the channel. However, some statistics are not yet perfectly symmetric after the large number of flow-through times simulated. Averaging these statistics over the two channel halves could be misleading, since the lack of convergence can be confused with physical effects that are not really present. Therefore, these specific profiles have not been averaged with respect to the symmetry of the channel. Consequently, all wall-normal profiles, including the symmetric ones, are consistently shown for 0 y / H 2 .

4.1. Average Wall Shear Stress and Viscosity Profiles

Typically, the convergence of bubble-laden channel flows towards a statistically steady state is measured using the average wall shear stress τ w [15,19]. By normalizing it with μ l u τ i n i t / H , τ w can be expressed in terms of the friction Reynolds number R e τ = τ w / ρ l H / ν l . Figure 4 shows R e τ , based on the wall-averaged shear stress τ w , for all setups as a function of the time. The starting time of statistical averaging is illustrated by the dotted line at t / T = 4.24 . As explained in Section 3.2, all cases are initialized from a single-phase channel flow of a Newtonian fluid controlled by a constant pressure gradient corresponding to R e τ = 127.3 . This value can be observed at t / T = 0 for the three cases corresponding to n = 1 . While R e τ of the single-phase case continues to fluctuate around this theoretical value when the flow control is changed to a constant volumetric flow rate, the presence of the dispersed gas fraction increases the average wall shear stress by approximately 20%. This effect was already reported in [19], where the authors attributed it to the stirring of the bubbles, which affects the near wall region above a certain gas fraction. The observed increase is 19.5% for E o = 0.3125 and 21.7% for E o = 3.75 . Although the time-averaged wall shear stress of the single-phase setup for n = 1 deviates from the theoretical value by less than 2%, the respective fluctuations are found to be stronger compared to the two-phase cases. The presence of the bubbles has a clearly visible dampening effect on the fluctuations of R e τ , irrespective of the power-law index.
As Figure 4 shows, the wall shear stress increases when the continuous liquid phase exhibits shear-thickening behavior ( n = 1.3 ), while the opposite effect can be observed for the shear-thinning scenario ( n = 0.7 ). This behavior is dominated by the viscosity values close to the walls. Due to the high shear rates γ ˙ in the near-wall region, the viscosity of the shear-thickening fluid exceeds the viscosity of the Newtonian fluid, while the viscosity of the shear-thinning fluid follows an opposite trend, see Figure 5. Although the mean velocity gradient at the wall is slightly higher (lower) for the shear-thinning (shear-thickening) fluid (see Figure 7), the wall shear stress still reveals the lowest (highest) values due to the viscosity behavior. The findings observed for the Newtonian fluid, i.e., smoother curves of R e τ with fluctuations around an increased average value for two-phase setups, also hold true for the non-Newtonian fluids.
Since the three cases in each setup category, i.e., single-phase or two-phase with either quasi-spherical or deformable bubbles, are identical apart from the different power-law indexes, it is expedient to initially investigate the viscosity statistics of each case. The viscosity values do not only explain the above-mentioned observations regarding the average wall shear stress, but they are also responsible for all other differences occurring in the results of one setup category. Figure 5 shows the average apparent dynamic viscosity of the liquid phase μ a l , normalized with the viscosity for n = 1 , as a function of the wall-normal coordinate y.
For the single-phase cases with a non-Newtonian continuous phase, the dynamic liquid viscosity obeys a monotonic trend from the wall to channel center. As the shear rate is high close to the walls and low in the channel center, the shear-thinning fluid produces the minimum viscosity at the walls, while the maximum can be found in the center. The opposite behavior can be observed for the shear-thickening fluid. While the profiles for the two-phase cases are relatively similar to those of the respective single-phase setups close to the walls, they clearly show the range where the bubbles are present, which spans from y / H 0.25 to y / H 1.75 (see Section 4.2). Due to the bubbles accumulating in this part of the domain, there are elevated local velocity differences causing the shear rate to exceed the values from the single-phase setups. As a consequence, the viscosity increases (decreases) for the shear-thickening (shear-thinning) fluid and the profiles for different power-law indexes are very close to each other in the section around the channel center. Moreover, due to the homogeneous bubbly flow in the core part of the channel (see Figures 7 and 9), e i j e i j , the average shear rate (see Section 2) as well as the average apparent viscosity are almost independent of y and constant in this region. It is also interesting to note that the Eötvös-number, and therefore the bubble deformability, seem to have a negligible influence on the viscosity profiles. Table 2 provides the normalized average viscosities of the liquid phase, which result from averaging over all spatial directions and over the time.
When investigating the average viscosity profiles for the bubble-laden flows, the quasi-identical values in the channel center lead to the question of whether there is some underlying physical cause for this behavior. As discussed in Section 2, the apparent liquid viscosity for the power-law fluids is determined as μ a = K γ ˙ n 1 . In order to observe identical average liquid viscosities across all investigated power-law indexes in the channel center, the average shear rate in this region would have to be equal to one. Since there is no obvious reason why this should generally be the case, it is assumed that this is just a characteristic of this specific setup and not a general observation. To verify this, a set of additional simulations has been performed at an increased void fraction of 10%. Indeed, the results support this assumption. For 120 instead of 72 bubbles, the average shear rate in the channel center increases further. This leads to a situation where, in contrast to the results illustrated in Figure 5, the viscosity profiles, including the peaks, for the shear-thinning fluid fall below the constant viscosity μ l = K of the Newtonian fluid for all values of y. The same effect, yet in the opposite direction, holds for the shear-thickening fluid: the minima are still larger than μ l = K . While, in principle, the shape of the profiles remains very similar, this shows that increased shear rates in the channel center due to higher gas fractions shift the viscosity profiles in the center towards the values at the wall. The single-phase setups in Figure 5 demonstrate this observation for the opposite extreme case. Another change to the setup that should, in principle, evoke similar consequences on the viscosity profiles would be to choose the power-law indexes for the shear-thinning and thickening fluid to be further away from 1. However, additionally verifying this assumption was beyond the scope of this work, as, with regard to Figure 5, it should only be clarified that the quasi-identical values in the channel center are indeed a characteristic of this specific setup and not a general observation.

4.2. First- and Second-Order Fluid Statistics

As was already reported in earlier studies investigating bubble-laden downflow channels [15,16,19], the bubbles show a clear tendency to accumulate in the channel center, see Figure 6. According to Lu and Tryggvason [15], this can be explained by the lift force acting on the bubbles. While all profiles exhibit the characteristics known from earlier studies, namely a relatively constant distribution in the channel center together with peaks in the transition to the bubble-free regions near the walls, there are clear tendencies arising for the shear-thinning and the shear-thickening fluid. It is important to note that the profiles shown in Figure 6 have not been averaged over the two channel halves. As can be recognized, the profiles are not perfectly symmetric. This indicates that the volume fraction field tends to converge slower than the other statistics provided in the following sections, which achieved a high degree of symmetry. Nevertheless, the profiles illustrated in Figure 6 are very similar to what is known from literature, with only one exception: the distribution is less uniform in the channel center. In order to make sure that this behavior only occurs due to a lack of convergence, identical setups have been simulated for a void fraction of 10%. As was observed in [19], statistics collected for dense cases tend to converge faster than those gathered at dilute gas fractions. For the increased gas fraction, it was possible to observe profiles that exhibit a high degree of symmetry together with flat profiles in the bubble-rich core. Consequently, it can be estimated that the profiles presented in Figure 6 would only flatten in the center for significantly longer simulation periods, while otherwise retaining their characteristics. Since the gathered statistics have already been computed over 180 flow-through times it is likely that, in order to reach void fraction profiles exhibiting perfect symmetry, a significantly longer simulation is required. As this study primarily aims to show tendencies arising for a non-Newtonian continuous liquid phase, it was decided to not continue the simulations.
When comparing the void fraction profiles for the cases with n = 1 to the profiles with a non-Newtonian liquid phase, clear tendencies can be observed. First of all, the bubble-free zone near the walls is wider for the shear-thickening and narrower for the shear-thinning setup. As a consequence, the profiles for n = 1.3 are shifted towards the channel center. Since the overall gas fraction is identical for all cases, this results in a stronger accumulation of bubbles in the core of the channel. Therefore, both the peaks in the transition between core and walls as well as the plateau-like section in the center are shifted to higher values. The opposite holds for n = 0.7 , where the bubbles are distributed more uniformly along the wall-normal direction. The shear lift force occurs due to the mean flow gradient near the wall. As the boundary layer becomes thinner (thicker) for shear-thinning (shear-thickening) fluids (see Figure 7), the bubbles move closer to (stay away further from) the wall compared to the Newtonian case. However, the peaks do not only increase for higher values of n. A comparison of the profiles for E o = 0.3125 and E o = 3.75 reveals that, independent of the power-law index, the profiles are also shifted towards the channel center for increased E o . This impacts the peaks in the same way as described for increased values of n.
The profiles for the average stream-wise velocity u , normalized with the friction velocity u τ i n i t to yield u + = u / u τ i n i t , are shown in Figure 7. Since the actual friction velocity u τ = τ w / ρ l is different for each of the nine setups (see Figure 4), the value u τ i n i t is taken from the initialization setup. The theoretical value for the respective case is provided in [15]. All profiles are similar to the ones shown in the literature, i.e., the single-phase setups show a monotonous increase towards their maximum value in the channel center, while the two-phase setups exhibit an almost uniform profile in the region of the bubble-rich core. In principle, the effect of increasing or decreasing n is the same as shown for the void fraction profiles in Figure 6. Once again, the profiles for the shear-thickening setups are shifted towards the channel core and, due to the constant volumetric flow rate, exhibit higher average velocities there. On the other hand, the profiles for the cases considering the shear-thinning fluid are wider and reach lower maximum flow velocities. The thickness of the wall boundary layer increases with n. This observation is consistent with the viscosity profiles illustrated in Figure 5, which demonstrate that the viscosity in the near-wall region is high for the shear-thickening fluid and low for the shear-thinning fluid. It is also interesting to note that a plateau-like section is formed in the middle of the two-phase velocity profiles, between the locations of the volume fraction peaks. In this regard, the profiles resemble the one shown in [18], although the plateaus are more pronounced in the present study, especially for n = 1.3 . This can be partly attributed to the fact that, instead of excluding the velocities of the gas phase, the velocities of both phases are taken into consideration here. Since the bubbles rise relative to the continuous liquid phase, including the gas-phase velocity would result in lower flow velocities in the bubble-rich channel core. For the two-phase setups, the velocity profiles exhibit only minor differences between E o = 0.3125 and E o = 3.75 . As can be seen, the peaks and plateau-like sections of the profiles are slightly more pronounced for the higher Eötvös number.
Figure 8 shows a semilogarithmic plot of the stream-wise velocity u , normalized by the friction velocity u τ to yield the dimensionless velocity U + , over the the dimensionless wall distance y + . It is important to note that, in this context, u τ = τ w / ρ l is determined using the actual value of the wall shear stress τ w for each setup (see Figure 4) instead of the theoretical value corresponding to the initialization setup, which was used for the profiles in Figure 7. This is indicated by using U + instead of u + in the following. For the calculation of the friction velocity, τ w is averaged over the same time span as all other statistics. The dimensionless wall distance is defined as y + = y u τ / ν . However, the kinematic viscosity ν is not constant for the non-Newtonian fluids. Since the aim is to investigate the flow behavior close to the wall, the viscosity values for the power-law fluids are taken from the first cells at the wall, i.e., ν corresponds to the average apparent viscosity values ν a of the liquid phase at y / H = 0 and y / H = 2 , see profiles shown in Figure 5. Since the near-wall region is free of bubbles, the kinematic viscosity of the gas phase does not play a role for the determination of y + . The resulting profiles are compared to the theoretical approximations suggested in the literature. For the viscous sublayer ( y + < 5 ), the predicted behavior is U + = y + . In the logarithmic layer ( y + > 30 ), the figure compares the normalized velocity to the logarithmic “law of the wall”, or “log law”, which is defined as U + = 1 / κ l n y + + B . Here, κ = 0.41 is the von Karman constant and B is taken to be 5.95 [15].
As Figure 8 shows, all nine profiles agree very well with U + = y + in the viscous sublayer. While the velocity of the single-phase flow at n = 1 also matches the “log law” very accurately, the profiles are shifted to higher values for the shear-thinning and to lower values for the shear-thickening fluid. Due to the absence of bubbles close to the walls, the two-phase profiles are relatively similar to those of the single-phase flow up to the position where the velocity peaks are formed. At larger distances from the wall, the stream-wise velocity profiles for the bubble-laden flows increasingly deviate and flatten towards the bubble-rich core region. As has also been reported in [15,18], the thickness of the wall boundary layer is reduced due to the presence of the bubbles and the formation of a logarithmic layer can therefore not be observed.
Figure 9, Figure 10 and Figure 11 show the root mean square of the stream-wise, the wall-normal and the span-wise velocity fluctuations, normalized by the friction velocity of the initialization case u τ i n i t . Overall, the fluctuations for n = 1 strongly resemble the profiles reported in the literature [15,18,19]. The introduction of a dispersed phase at a significant void fraction leads to clear deviations of all profiles from their single-phase counterparts. The stream-wise velocity fluctuations for the single-phase setups, illustrated in Figure 9, reach their maximum close to the wall, followed by a decrease towards the channel center. Independent of the power-law index, the opposite behavior can be observed for the two-phase setups. Figure 9 shows that, although the first peaks of the bubble-laden flows are located at the same values of y, they are significantly lower. This indicates a strong turbulent intensity reduction in the near-wall region. As discussed in [18], the reason for this observation might be that, due to their lower velocity compared to the liquid bulk, the bubbles act as obstacles and intermit the otherwise present turbulence structures. On the other hand, compared to the single-phase setups, the two-phase flows reveal increasing fluctuations of the stream-wise velocity component in the channel core, which are generated by the high number of bubbles present in this region [15]. As can be seen, the stream-wise fluctuations are virtually independent of y in the bubble-rich core, where the bubbles are distributed evenly (see Figure 6).
The latter observation also holds true for the velocity fluctuations in the other spatial directions. While the two-phase profiles for the wall-normal fluctuations, illustrated in Figure 10, show a slight increase from the locations of the void fraction peaks towards the channel center for E o = 0.3125 , the profiles for E o = 3.75 and those for the span-wise velocity fluctuations, shown in Figure 11, are relatively independent of y in the middle section. As shown in [18,19], all fluctuation profiles generally tend to flatten completely in the core region, if the void fraction is increased to values of 10% or more.
In the next step, the influence of the power-law index on the velocity fluctuations is investigated. For the single-phase setups, the difference in the velocity fluctuations close to the wall is insignificant. However, this changes towards the channel center, where the highest fluctuations can be observed for n = 1.3 , while the single-phase flow at n = 0.7 reveals the lowest fluctuations. These findings hold true for the stream-wise, the wall-normal and the span-wise direction. The observed profiles can be explained by the differences in the shear rate. For the single-phase channel, the shear rates are high in the near-wall region and low in the channel center. As a consequence, the shear-thinning (shear-thickening) fluid exhibits a low (high) dynamic viscosity close to the wall and a high (low) dynamic viscosity in the channel center, see Figure 5. As the strength of the fluctuations is expected to increase with decreasing dynamic viscosity, this explains the trends observed for the single-phase setups. However, this effect is partially countered by the increasing gradients close to the wall for decreasing power-law indexes, such that the peaks of the stream-wise velocity fluctuations remain more or less comparable for all cases.
As shown in Figure 9, Figure 10 and Figure 11, the two-phase profiles for the fluctuations of all three velocity components reveal remarkable differences compared to their single-phase counterparts. It is especially interesting to note that the order regarding the fluctuation strength is reversed in the bubble-rich core once a significant amount of bubbles is present. The two-phase setups with a shear-thinning continuous phase not only exhibit the strongest fluctuations close to the walls, but also in the bubble-rich core and for almost all values of y / H in the transition between these regions. While the three profiles of the stream-wise fluctuations overlap marginally for some small sections of y / H , the wall-normal and span-wise velocity fluctuations for n = 0.7 are higher than those for n = 1.3 independent of the wall-normal coordinate. The profiles for n = 1 are generally located between those of the non-Newtonian fluids apart from small overlaps on each side. The introduction of a significant number of bubbles leads to higher local velocity differences and increasing shear rates in the bubble-rich core. As a result, the apparent viscosity of the shear-thinning fluid decreases compared to the single-phase setup, while the opposite effect holds for the shear-thickening fluid. However, as shown in Figure 5, the profiles of the average apparent viscosities are quasi identical in the channel core and therefore cannot fully explain the results observed for the fluctuations. The findings can, however, be reasoned by the different viscosities in the vicinity of the phase interface, where the velocity fluctuations are generated. For the single-phase setups, the fluctuations for y / H -positions very close to the walls are strongest (weakest) for n = 0.7 ( n = 1.3 ), as the shear rate is high close to the walls. As an additional analysis revealed, high shear rates occur in the vicinity of the phase interface as well. This can also be recognized in Figure 15, where strong velocity gradients are clearly visible around the bubbles. As a consequence, the apparent liquid viscosity for the shear-thinning (shear-thickening) fluid is significantly lower (higher) than the one of the Newtonian fluid in the vicinity of the phase interface. Although, partly due to the migration of the bubbles, this effect is not visible in the averaged profiles (see Figure 5), it is responsible for the stronger (weaker) velocity fluctuations (compared to n = 1 ) observed for the shear-thinning (shear-thickening) fluid in the core region. Moreover, the fluctuations in the core occur mainly due to the vortex shedding in the bubble wake and the velocity differences between the bubbles and the surrounding liquid. As discussed in Section 3.2, all other conditions remaining identical, the relative velocity of the bubbles, and therefore also the bubble Reynolds number, increases with decreasing power-law index. This also explains the trends observed for the velocity fluctuations in the channel core.
Compared to the strong influence of the power-law index, the changes observed when increasing the Eötvös number from 0.3125 to 3.75 are relatively small. For the stream-wise velocity fluctuations (see Figure 9), the peaks in the near-wall region are slightly more pronounced. Moreover, a minimal displacement of the profiles away from each other in the transition to the bubble-rich core can be observed. The same applies to the span-wise (Figure 11) and the wall-normal (Figure 10) fluctuations. For the latter, E o = 3.75 additionally leads to a flattening in the center. However, it is interesting to note that the fluctuation of all velocity components is generally stronger for the lower Eötvös number. Although the differences are small when comparing the velocity fluctuation profiles, this effect is clearly visible when investigating the turbulent kinetic energy k, see Figure 13. Further analysis regarding this observation, as well as the impact on the dissipation rate ε , can be found in Section 4.3.

4.3. Dissipation Rates

In this section, the dissipation rate ε of the turbulent kinetic energy k = 1 2 ( u i u i ) is investigated. It can be defined as
ε = ν a u i x k u i x k
with the apparent kinematic viscosity ν a and the velocity fluctuations u i . The averaging over time and the homogeneous directions are consistent with the procedure used for the first and second order fluid statistics (see Section 4.2). The apparent kinematic viscosity is ν a = μ a / ρ . Figure 12 shows the dissipation rate, normalized with the reference value ε r e f = ( μ l / ρ l ) ( u τ i n i t / H ) 2 for n = 1 and evaluated without phase-dependent conditioning, as a function of the wall-normal coordinate. The profiles have not been averaged over the two channel halves.
For the three single-phase setups, the dissipation rate reaches its maximum at the walls, while the minimum is located in the channel center. In the near-wall region, the dissipation rate decreases with increasing power-law index. However, the order reverses at a distance of approximately H / 8 from the walls. This behavior can be deduced from the fluctuations of the three velocity components (see Figure 9, Figure 10 and Figure 11), where similar trends can be observed. It is worth noting that viscosity plays a role as well, but the velocity fluctuation gradients have a quadratic contribution to the dissipation rate. The two-phase setups are dominated by the presence of the bubbles in the core of the channel. Since ρ l / ρ g = 10 (see Section 3.2), the kinematic viscosity of the gas phase is significantly higher than the apparent kinematic viscosity of the liquid phase. For the Newtonian fluid, ν g / ν l is 10. As a consequence, there is a direct relation between the void fraction (see Figure 6) and the dissipation rate profiles (see Figure 12) and the corresponding profiles exhibit a very similar shape. The two-phase profiles do not only deviate from their single-phase counterparts in the center of the channel, but also close to the walls. This reflects the significant deviations of the two-phase velocity fluctuation profiles from their single-phase counterparts, occurring for practically all values of y (see Figure 9, Figure 10 and Figure 11).
The Eötvös number has a negligible influence on the dissipation rate in the near-wall region. In the channel center, however, the dissipation rates for E o = 0.3125 are significantly higher than those for E o = 3.75 . This can be reasoned by the behavior of the average turbulent kinetic energy k , illustrated in Figure 13, and the integral turbulent length scale . The relation between the dissipation rate and the turbulent kinetic energy is ε = k 3 / 2 / . As shown, the profiles of k for the setups at E o = 0.3125 are shifted to higher values compared to their counterparts at E o = 3.75 . As the difference in k is less pronounced than the difference in ε , the behavior of must also be considered. As will also be discussed in Section 4.4, the bubbles are deformed to an ellipsoidal shape for E o = 3.75 , while they remain quasi-spherical for E o = 0.3125 . For spherical bubbles, it is possible to assume d b . For the ellipsoidal bubbles, on the other hand, the integral length scale approximately corresponds to the lateral expansion of the ellipsoid perpendicular to the main flow direction. Therefore, it is reasonable to assume that > d b for the setups at E o = 3.75 . Combined, the differences in k and explain the observations for the dissipation rate. It is, however, not trivial to identify the reason for the dependency of k on the the bubble deformability. As discussed in Section 3.2, the bubble Reynolds numbers R e b are higher for the setups at E o = 0.3125 (see Table 1), which is indicative of stronger velocity fluctuations in the bubble wake. These higher values of R e b are caused by higher stream-wise velocity differences between the bubbles and the liquid phase. Nearly spherical bubbles encounter less resistance while rising through the liquid, which allows them to rise faster than the deformable bubbles. Moreover, it is worth noting that there is a constant interplay between interfacial surface energy and turbulent kinetic energy, as discussed by Dodd and Ferrante [36].
In order to exclude the effect the high kinematic viscosity of the bubbles has on the dissipation rate, Figure 14 provides the average normalized dissipation rate conditioned on the liquid phase ε l . The profiles have not been averaged over the two channel halves. The single-phase profiles, given for comparison, are identical to those shown in Figure 12. As observed for ε , the presence of the bubbles completely changes the profiles for all power-law indexes. In the near-wall region, the two-phase cases for n = 1 and n = 1.3 exhibit significantly lower dissipation rates compared to the single-phase flow. For n = 0.7 , this observation only holds true directly at the wall. At larger distances from the wall and especially in the core of the channel, the bubbles significantly increase the liquid phase dissipation rate irrespective of the power-law index. For n = 1 , the dissipation rate in the channel center is approximately as high as its peak value in the near-wall region. It is especially interesting to note that, for the shear-thickening fluid, the presence of the bubbles completely reverses the trends resulting for the single-phase setups. The minimum dissipation rate can be observed at the walls. The highest dissipation rates occur from y / H 0.5 to y / H 1.5 in the form of a plateau-like section. This behavior can partly be attributed to the apparent viscosity in the channel core, which, compared to the corresponding single-phase flow, reaches higher values for the bubbly flow setups at n = 1.3 (see Figure 5). Moreover, we refer to the relation between k and ε explained above. Compared to the single-phase flows, including the bubbles significantly lowers the values of k close to the walls (see Figure 13), in particular for n = 1.3 . Together with the higher turbulent kinetic energy in the channel core, which can be observed irrespective of n, this explains the two-phase profiles of the dissipation rate at n = 1.3 , shown in Figure 14.

4.4. Bubble Deformation

As explained in Section 3.2, the dynamic bubble behavior can be estimated from the Eötvös number E o and the bubble Reynolds number R e b using the Grace diagram, see Figure 3. A straightforward way to investigate the bubble shapes occurring for the different setups is to visualize the isocontours corresponding to f = 0.5 , i.e., the interface between gas and liquid phase. Snapshots for the bubbly-flow cases with a Newtonian fluid in the continuous phase are shown in Figure 15. Since there are hardly any visible differences for the non-Newtonian fluids, the respective plots are omitted. The isosurfaces confirm the prediction of the Grace diagram. As expected, the bubbles are quasi-spherical for E o = 0.3125 . However, when E o is increased to 3.75 , the bubbles are wobbling and assume an ellipsoidal shape.
While this allows a first qualitative assessment of the bubble deformations, snapshots at a single time step do not yield further insights into the transient deformation behavior as well as the effects of different power-law indexes. In order to provide a meaningful quantitative measure for the bubble deformation, it is expedient to analyze the temporal behavior of the interface area for the single bubbles. For this purpose, the interface area of a bubble b is numerically approximated as the volume integral of δ S = | f b | , which can be interpreted as the surface density, over all cells belonging to the bubble’s individual VOF marker field. It is necessary to perform this analysis for all bubbles separately. Otherwise, the oscillations of the individual bubbles would annihilate each other due to the different instantaneous deformation behavior. In the following, A b t denotes the the surface area of bubble b at the time t, normalized by the surface area of a sphere of the same volume V. The time-averaged, normalized surface area is given as
A b ¯ = 1 t e t s t s t e A b t d t
with the start and end times t s and t e for gathering statistics. The fluctuation of the normalized surface area around this value is then evaluated using the root mean square:
A b = 1 t e t s t s t e A b t A b ¯ 2 d t
Finally, the time-averaged, normalized surface areas A b ¯ and the fluctuations A b are averaged over all 72 bubbles. The results are presented in Table 3.
Table 3 shows that the normalized surface area averaged over time and all bubbles grows by 6.40% for the shear-thinning, by 5.84% for the Newtonian and by 5.08% for the shear-thickening fluid when the Eötvös number is increased from 0.3125 to 3.75. Moreover, it is clearly visible that the bubbles are not only deformed to ellipsoids for E o = 3.75 . Measured in terms of their surface area fluctuations, they are also subject to much stronger oscillations. Compared to E o = 0.3125 , the values of A increase by a factor of 40.42 for n = 0.7 , by a factor of 48.25 for n = 1 and by a factor of 55.73 for n = 1.3 . Finally, A ¯ and A show clear tendencies across the investigated power-law indexes. Both quantities increase for decreasing values of n. This behavior can be explained by the values of the bubble Reynolds number R e b given in Table 1, which also increase with decreasing power-law index. As shown in the Grace diagram (see Figure 3), the bubble deformability is governed by the combination of E o and R e b . At identical E o , higher values of R e b are indicative of stronger deformations and oscillations. It was found that the differences in the values of R e b are almost exclusively caused by the different relative stream-wise velocities u r e l of both phases, see Equation (8). On the other hand, the apparent liquid dynamic viscosities in the center part of the channel ( 0.5 < y / H < 1.5 ) seem to only have a negligible influence, as can be seen in Figure 5.

5. Conclusions

In this paper, turbulent bubble-laden channel flow of power-law fluids at a bulk Reynolds number equal to 3790 was studied by means of Direct Numerical Simulation (DNS). The power-law index n was varied between 0.7 (shear-thinning), 1 (Newtonian) and 1.3 (shear-thickening). The downflow channel was loaded with 72 either quasi-spherical ( E o = 0.3125 ) or ellipsoidal/wobbling ( E o = 3.75 ) bubbles, which resulted in a total gas fraction of 6%. In addition, corresponding single-phase setups were investigated for comparison. The results demonstrate the intricate interaction of non-Newtonian and bubbly flows:
For the single-phase flow at n = 0.7 ( n = 1.3 ), the apparent viscosity of the liquid phase reaches its maximum (minimum) directly in the center of the channel, due to the low shear rates away from the walls. It was observed that the presence of the dispersed phase significantly increases the shear rate in the channel core, where the bubbles accumulate. Therefore, compared to the single-phase flow at identical power-law index, the apparent viscosity of the liquid phase in this region decreases (increases) for the shear-thinning (shear-thickening) fluid. In comparison to the single-phase flow, a flattening of the viscosity profiles, irrespective of E o , can be observed in the channel core due to the homogeneous bubbly flow in this region. In the bubble-free near-wall region, the behavior of the viscosity remains nearly unchanged when including the bubbles. The bubble-free zone near the walls is wider for n = 1.3 , which leads to a stronger accumulation of bubbles in the core. The opposite effect is observed for n = 0.7 , i.e., the bubbles, on average, move closer to the walls. The profiles of the average stream-wise velocity demonstrate that the thickness of the wall boundary layer increases with n. In the near-wall region, the shear-thickening fluid exhibits the lowest flow velocities, while the highest velocities are observed for the shear-thinning fluid. Since the constant volumetric flow rate is identical for all setups, this order reverses in the core of the channel. As reported in other studies, the velocity profiles of the two-phase setups flatten and form a plateau-like section in the bubble-rich ore.
The presence of the bubbles significantly changes the fluctuation profiles of all three velocity components. For the single-phase flows, the strength of the fluctuations in the channel core increases with n. For the two-phase setups, this order is reversed. It is interesting to note that, overall, the velocity fluctuations and therefore the turbulent kinetic energy for E o = 0.3125 are higher compared to E o = 3.75 . For the less deformed bubbles at E o = 0.3125 , this difference, together with the lower integral turbulent length scale due to the quasi-spherical bubble shape, leads to a higher turbulence dissipation rate in the channel core. It was additionally shown that, compared to the single-phase setups, the presence of the dispersed phase significantly increases the dissipation rate in the region where the bubbles accumulate. For n = 1.3 , the trend observed for the single-phase setup reverses entirely, i.e., the dissipation rate of the liquid phase reaches its maximum in the bubble-rich core, while the minimum was found to be located directly at the walls.
Finally, an analysis of the mean value and the fluctuation of the bubble surface areas revealed that, irrespective of the Eötvös number, the average deformation and oscillation strength increases with decreasing power-law index. As expected, the average deformation and oscillation strength were observed to be more pronounced for E o = 3.75 .
The present analysis has been carried out with the power-law model. To take into account additional effects present in real-world non-Newtonian fluids, more sophisticated models, such as the Bird–Carreau or Cross-Power-Law model, can be considered for further investigations.

Author Contributions

Conceptualization, M.K. and J.H.; formal analysis, F.B. and E.T.; Investigation, F.B.; Software, F.B., E.T., and P.C.; Supervision, M.K., E.T., and J.H.; Visualization, E.T.; writing—original draft, E.T.; writing—review and editing, M.K., J.H., and P.C. All authors have read and agreed to the published version of the manuscript.

Funding

Support by the German Research Foundation (Deutsche Forschungsgemeinschaft (DFG), GS: KL1456/4-1) is gratefully acknowledged.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A. Parameters for the Case Definition

To allow a simple reproduction of the investigated setup, Table A1 provides the parameters used for the case definition. These are selected to make the setup comparable to the study presented in [15]. The dimensionless numbers given in Section 3.2 fully define the case.
Table A1. Parameters for the case definition.
Table A1. Parameters for the case definition.
SymbolParameterValue
HChannel half width1
ρ l Density of the liquid phase1
ρ g Density of the gas phase 0.1
μ l Dynamic viscosity of the liquid phase for n = 1 1 / 3000
μ g Dynamic viscosity of the gas phase 1 / 3000
KConsistency index 1 / 3000
d b Bubble diameter 0.25
gGravitational acceleration 0.1
σ Surface tension coefficient 0.02 ( E o = 0.3125 ), 1.667 × 10 3 ( E o = 3.75 )
NNumber of bubbles72
V ˙ ( and V ˙ / L z ) Constant volumetric flow rate 1.2636 ( 0.8044 )
u τ i n i t Friction velocity of the initialization setup 0.0424

References

  1. Chhabra, R.P. Bubbles, Drops, and Particles in Non-Newtonian Fluids, 2nd ed.; CRC Press: Boca Raton, FL, USA, 2006. [Google Scholar]
  2. Gestring, I.; Mewes, D. Degassing of molten polymers. Chem. Eng. Sci. 2002, 57, 3415–3426. [Google Scholar] [CrossRef]
  3. Kamisli, F.; Ryan, M.E. Perturbation method in gas-assisted power-law fluid displacement in a circular tube and a rectangular channel. Chem. Eng. J. 1999, 75, 167–176. [Google Scholar] [CrossRef]
  4. Kamisli, F.; Ryan, M.E. Gas-assisted non-newtonian fluid displacement in circular tubes and non-circular channels. Chem. Eng. Sci. 2001, 56, 4913–4928. [Google Scholar] [CrossRef]
  5. Yamamoto, T.; Suga, T.; Nakamura, K.; Mori, N. The Gas Penetration Through Viscoelastic Fluids With Shear-Thinning Viscosity in a Tube. J. Fluids Eng. 2004, 126, 148–152. [Google Scholar] [CrossRef]
  6. Delgado, M.A.; Franco, J.M.; Partal, P.; Gallegos, C. Experimental study of grease flow in pipelines: Wall slip and air entrainment effects. Chem. Eng. Process. 2005, 44, 805–817. [Google Scholar] [CrossRef]
  7. Hecht, V.; Voigt, J.; Schugerl, K. Absorption of oxygen in countercurrent multistage bubble columns—III Viscoelastic liquids. Comparison of systems with high viscosity. Chem. Eng. Sci. 1980, 35, 1325–1330. [Google Scholar] [CrossRef]
  8. Sada, E.; Kumazawa, H.; Lee, C.H. Chemical absorption in a bubble column loading concentrated slurry. Chem. Eng. Sci. 1983, 38, 2047–2051. [Google Scholar] [CrossRef]
  9. Chhabra, R.P. Wall effects on terminal velocity of non-spherical particles in non-Newtonian polymer solutions. Powder Technol. 1996, 1, 39–44. [Google Scholar] [CrossRef]
  10. Vandu, C.O.; Koop, K.; Krishna, R. Large Bubble Sizes and Rise Velocities in a Bubble Column Slurry Reactor. Chem. Eng. Technol. 2004, 27, 1195–1199. [Google Scholar] [CrossRef]
  11. Frank, X.; Li, H.Z.; Funfschilling, D.; Burdin, F.; Ma, Y. Bubble Motion in Non-Newtonian Fluids and Suspensions. Can. J. Chem. Eng. 2003, 81, 483–490. [Google Scholar] [CrossRef]
  12. Frank, X.; Li, H.Z.; Funfschilling, D. Bubble Motion in Non-Newtonian Fluids and Suspensions. Eur. Phys. J. E 2005, 16, 29–35. [Google Scholar] [CrossRef]
  13. Lau, Y.M.; Deen, N.G.; Kuipers, J.A.M. Development of an image measurement technique for size distribution in dense bubbly flows. Chem. Eng. Sci. 2013, 94, 20–29. [Google Scholar] [CrossRef]
  14. Mercado, J.M.; Gomez, D.C.; Van Gils, D.; Sun, C.; Lohse, D. On bubble-clustering and energy spectra in pseudo-turbulence. J. Fluid Mech. 2010, 650, 287–306. [Google Scholar] [CrossRef] [Green Version]
  15. Lu, J.; Tryggvason, G. Numerical study of turbulent bubbly downflows in a vertical channel. Phys. Fluids 2006, 18, 103302. [Google Scholar] [CrossRef]
  16. Lu, J.; Tryggvason, G. Effect of bubble size in turbulent bubbly downflow in a vertical channel. Chem. Eng. Sci. 2007, 62, 3008–3018. [Google Scholar] [CrossRef]
  17. Lu, J.; Tryggvason, G. Effect of bubble deformability in turbulent bubbly upflow in a vertical channel. Phys. Fluids 2008, 20, 040701. [Google Scholar] [CrossRef] [Green Version]
  18. Cifani, P.; Kuerten, J.G.M.; Geurts, B.J. Highly scalable DNS solver for turbulent bubble-laden channel flow. Comput. Fluids 2018, 172, 67–83. [Google Scholar] [CrossRef] [Green Version]
  19. Cifani, P.; Kuerten, J.G.M.; Geurts, B.J. Flow and bubble statistics of turbulent bubble-laden downflow channel. Int. J. Multiph. Flow 2020, 126, 103244. [Google Scholar] [CrossRef]
  20. Chu, X.; Liu, Y.; Wang, W.; Yang, G.; Weigand, B.; Nemati, H. Turbulence, pseudo-turbulence, and local flow topology in dispersed bubbly flow. Phys. Fluids 2020, 32, 083310. [Google Scholar] [CrossRef]
  21. Hasslberger, J.; Cifani, P.; Chakraborty, N.; Klein, M. A direct numerical simulation analysis of coherent structures in bubble-laden channel flows. J. Fluid Mech. 2020, 905, A37. [Google Scholar] [CrossRef]
  22. Ohta, T.; Miyashita, M. DNS and LES with an extended Smagorinsky model for wall turbulence in non-Newtonian viscous fluids. J. Nonnewton. Fluid Mech. 2014, 206, 29–39. [Google Scholar] [CrossRef]
  23. Ozoe, H.; Churchill, S.W. Hydrodynamic stability and natural convection in Ostwald-de Waele and Ellis fluids: The development of a numerical solution. AIChE J. 1972, 18, 1196–1207. [Google Scholar] [CrossRef]
  24. Bird, R.B.; Stewart, W.E.; Lightfoot, E.N. Transport Phenomena, 2nd ed.; John Wiley & Sons: New York, NY, USA, 2007. [Google Scholar]
  25. Prosperetti, A.; Tryggvason, G. Computational Methods for Multiphase Flow; Cambridge University Press: Cambridge, UK, 2007. [Google Scholar]
  26. Brackbill, J.U.; Kothe, D.B.; Zemach, C. A continuum method for modeling surface tension. J. Comput. Phys. 1992, 100, 335–354. [Google Scholar] [CrossRef]
  27. Popinet, S. An accurate adaptive solver for surface-tension driven interfacial flows. J. Comput. Phys. 2009, 228, 5838–5866. [Google Scholar] [CrossRef] [Green Version]
  28. Hirt, C.W.; Nichols, B.D. Volume of Fluid (VOF) method for the dynamics of free boundaries. J. Comput. Phys. 1981, 39, 201–225. [Google Scholar] [CrossRef]
  29. Coyajee, E.; Boersma, B.J. Numerical simulation of drop impact on a liquid-liquid interface with a multiple marker front-capturing method. J. Comput. Phys. 2009, 228, 4444–4467. [Google Scholar] [CrossRef]
  30. Cifani, P.; Michalek, W.R.; Priems, G.J.M.; Kuerten, J.G.M.; van der Geld, C.W.M.; Geurts, B.J. A comparison between the surface compression method and an interface reconstruction method for the VOF approach. Comput. Fluids 2016, 136, 421–435. [Google Scholar] [CrossRef] [Green Version]
  31. Dodd, M.S.; Ferrante, A. A fast pressure-correction method for incompressible two-fluid flows. J. Comput. Phys. 2014, 273, 416–434. [Google Scholar] [CrossRef]
  32. Cifani, P. Analysis of a constant-coefficient pressure equation method for fast computations of two-phase flows at high density ratios. J. Comput. Phys. 2019, 398, 108904. [Google Scholar] [CrossRef]
  33. Jiménez, J.; Moin, P. The minimal flow unit in near-wall turbulence. J. Fluid Mech. 1991, 225, 213–240. [Google Scholar] [CrossRef] [Green Version]
  34. Grace, J.R.; Wairegi, T.; Nguyen, T.H. Shapes and Velocities of Single Drops and Bubbles Moving Freely through Immiscible Liquids. Trans. Inst. Chem. Eng. 1976, 54, 167–173. [Google Scholar]
  35. Dabiri, S.; Lu, J.; Tryggvason, G. Transition between regimes of a vertical channel bubbly upflow due to bubble deformability. Phys. Fluids 2013, 25, 102110. [Google Scholar] [CrossRef]
  36. Dodd, M.S.; Ferrante, A. On the interaction of Taylor length scale size droplets and isotropic turbulence. J. Fluid Mech. 2016, 806, 356–412. [Google Scholar] [CrossRef] [Green Version]
Figure 1. Apparent dynamic viscosity μ a , normalized with the viscosity of the Newtonian fluid μ , as a function of the shear rate γ ˙ .
Figure 1. Apparent dynamic viscosity μ a , normalized with the viscosity of the Newtonian fluid μ , as a function of the shear rate γ ˙ .
Fluids 06 00040 g001
Figure 2. Downflow channel configuration. The gray areas represent the no-slip walls, the x- and z-direction are treated with periodic boundaries.
Figure 2. Downflow channel configuration. The gray areas represent the no-slip walls, the x- and z-direction are treated with periodic boundaries.
Fluids 06 00040 g002
Figure 3. Grace diagram for bubble shapes [34]. Purple point and lines specify case “tph10_loEo”, blue point and lines specify case “tph10_hiEo”.
Figure 3. Grace diagram for bubble shapes [34]. Purple point and lines specify case “tph10_loEo”, blue point and lines specify case “tph10_hiEo”.
Fluids 06 00040 g003
Figure 4. Total instantaneous wall-averaged shear stress τ w , normalized by μ l ( u τ i n i t / H ) to yield R e τ , as a function of the dimensionless time t / T for the single-phase (a) and the two-phase setups with E o = 0.3125 (b) and E o = 3.75 (c). Statistics are collected starting at t / T = 4.24 , which is indicated by the vertical dotted line.
Figure 4. Total instantaneous wall-averaged shear stress τ w , normalized by μ l ( u τ i n i t / H ) to yield R e τ , as a function of the dimensionless time t / T for the single-phase (a) and the two-phase setups with E o = 0.3125 (b) and E o = 3.75 (c). Statistics are collected starting at t / T = 4.24 , which is indicated by the vertical dotted line.
Fluids 06 00040 g004
Figure 5. Average apparent dynamic viscosity of the liquid phase μ a l , normalized by the dynamic viscosity of the Newtonian fluid μ l , as a function of the wall-normal coordinate for the single-phase (a) and the two-phase setups with E o = 0.3125 (b) and E o = 3.75 (c).
Figure 5. Average apparent dynamic viscosity of the liquid phase μ a l , normalized by the dynamic viscosity of the Newtonian fluid μ l , as a function of the wall-normal coordinate for the single-phase (a) and the two-phase setups with E o = 0.3125 (b) and E o = 3.75 (c).
Fluids 06 00040 g005
Figure 6. Average gas volume fraction field f ¯ as a function of the wall-normal coordinate for the two-phase setups with E o = 0.3125 (a) and E o = 3.75 (b).
Figure 6. Average gas volume fraction field f ¯ as a function of the wall-normal coordinate for the two-phase setups with E o = 0.3125 (a) and E o = 3.75 (b).
Fluids 06 00040 g006
Figure 7. Average stream-wise velocity, normalized by the friction velocity of the initialization setup to yield u + , as a function of the wall-normal coordinate for the single-phase (a) and the two-phase setups with E o = 0.3125 (b) and E o = 3.75 (c).
Figure 7. Average stream-wise velocity, normalized by the friction velocity of the initialization setup to yield u + , as a function of the wall-normal coordinate for the single-phase (a) and the two-phase setups with E o = 0.3125 (b) and E o = 3.75 (c).
Fluids 06 00040 g007
Figure 8. Semilogarithmic plot of the average stream-wise velocity, normalized by the actual friction velocity to yield U + , as a function of the dimensionless distance from the wall for the single-phase (a) and the two-phase setups with E o = 0.3125 (b) and E o = 3.75 (c).
Figure 8. Semilogarithmic plot of the average stream-wise velocity, normalized by the actual friction velocity to yield U + , as a function of the dimensionless distance from the wall for the single-phase (a) and the two-phase setups with E o = 0.3125 (b) and E o = 3.75 (c).
Fluids 06 00040 g008
Figure 9. Root mean square of the stream-wise velocity fluctuations u , normalized by the friction velocity of the initialization setup, as a function of the wall-normal coordinate for the single-phase (a) and the two-phase setups with E o = 0.3125 (b) and E o = 3.75 (c).
Figure 9. Root mean square of the stream-wise velocity fluctuations u , normalized by the friction velocity of the initialization setup, as a function of the wall-normal coordinate for the single-phase (a) and the two-phase setups with E o = 0.3125 (b) and E o = 3.75 (c).
Fluids 06 00040 g009
Figure 10. Root mean square of the wall-normal velocity fluctuations v , normalized by the friction velocity of the initialization setup, as a function of the wall-normal coordinate for the single-phase (a) and the two-phase setups with E o = 0.3125 (b) and E o = 3.75 (c).
Figure 10. Root mean square of the wall-normal velocity fluctuations v , normalized by the friction velocity of the initialization setup, as a function of the wall-normal coordinate for the single-phase (a) and the two-phase setups with E o = 0.3125 (b) and E o = 3.75 (c).
Fluids 06 00040 g010
Figure 11. Root mean square of the span-wise velocity fluctuations w , normalized by the friction velocity of the initialization setup, as a function of the wall-normal coordinate for the single-phase (a) and the two-phase setups with E o = 0.3125 (b) and E o = 3.75 (c).
Figure 11. Root mean square of the span-wise velocity fluctuations w , normalized by the friction velocity of the initialization setup, as a function of the wall-normal coordinate for the single-phase (a) and the two-phase setups with E o = 0.3125 (b) and E o = 3.75 (c).
Fluids 06 00040 g011
Figure 12. Average dissipation rate ε , normalized by ε r e f , as a function of the wall-normal coordinate for the single-phase (a) and the two-phase setups with E o = 0.3125 (b) and E o = 3.75 (c).
Figure 12. Average dissipation rate ε , normalized by ε r e f , as a function of the wall-normal coordinate for the single-phase (a) and the two-phase setups with E o = 0.3125 (b) and E o = 3.75 (c).
Fluids 06 00040 g012
Figure 13. Average turbulent kinetic energy k , normalized by ( u τ i n i t ) 2 , as a function of the wall-normal coordinate for the single-phase (a) and the two-phase setups with E o = 0.3125 (b) and E o = 3.75 (c).
Figure 13. Average turbulent kinetic energy k , normalized by ( u τ i n i t ) 2 , as a function of the wall-normal coordinate for the single-phase (a) and the two-phase setups with E o = 0.3125 (b) and E o = 3.75 (c).
Fluids 06 00040 g013
Figure 14. Average dissipation rate of the liquid phase ε l , normalized by ε r e f , as a function of the wall-normal coordinate for the single-phase (a) and the two-phase setups with E o = 0.3125 (b) and E o = 3.75 (c).
Figure 14. Average dissipation rate of the liquid phase ε l , normalized by ε r e f , as a function of the wall-normal coordinate for the single-phase (a) and the two-phase setups with E o = 0.3125 (b) and E o = 3.75 (c).
Fluids 06 00040 g014
Figure 15. Isosurfaces for volume fraction f = 0.5 (gray) and wall-normal slice of the velocity magnitude on the mid-plane. (a) setup with n = 1 and E o = 0.3125 ; (b) setup with n = 1 and E o = 3.75 . It is worth noting that some bubbles intersect with the periodic boundaries, which should not be confused with a concave surface topology.
Figure 15. Isosurfaces for volume fraction f = 0.5 (gray) and wall-normal slice of the velocity magnitude on the mid-plane. (a) setup with n = 1 and E o = 0.3125 ; (b) setup with n = 1 and E o = 3.75 . It is worth noting that some bubbles intersect with the periodic boundaries, which should not be confused with a concave surface topology.
Fluids 06 00040 g015
Table 1. DNS case overview in terms of power-law index n, bubble count N b , Eötvös number E o , bubble Reynolds number R e b , and dimensionless distance from the wall y + .
Table 1. DNS case overview in terms of power-law index n, bubble count N b , Eötvös number E o , bubble Reynolds number R e b , and dimensionless distance from the wall y + .
Casen N b Eo Re b y +
sph070.70--0.9574
sph1010--0.6178
sph131.30--0.4428
tph07_loEo0.7720.3125132.301.1743
tph10_loEo1720.3125123.350.6687
tph13_loEo1.3720.3125117.580.4559
tph07_hiEo0.7723.75109.121.1860
tph10_hiEo1723.75105.850.6750
tph13_hiEo1.3723.75100.510.4582
Table 2. Average normalized liquid dynamic viscosities for all nine setups. The viscosity fields have been averaged over all spatial directions and over time, which is indicated by · x , y , z , t .
Table 2. Average normalized liquid dynamic viscosities for all nine setups. The viscosity fields have been averaged over all spatial directions and over time, which is indicated by · x , y , z , t .
μ a l x , y , z , t / μ l
Single Phase Eo = 0.3125 Eo = 3.75
n = 0.71.230.970.97
n = 1111
n = 1.30.991.081.09
Table 3. Mean A ¯ and fluctuation A of the normalized surface area, averaged over all bubbles.
Table 3. Mean A ¯ and fluctuation A of the normalized surface area, averaged over all bubbles.
Eo = 0.3125 Eo = 3.75
n = 0.7 n = 1 n = 1.3 n = 0.7 n = 1 n = 1.3
A ¯ 1.0190 1.0188 1.0187 1.0842 1.0783 1.0705
A 6.8385 × 10 4 4.8121 × 10 4 3.2722 × 10 4 2.7642 × 10 2 2.3220 × 10 2 1.8235 × 10 2
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Bräuer, F.; Trautner, E.; Hasslberger, J.; Cifani, P.; Klein, M. Turbulent Bubble-Laden Channel Flow of Power-Law Fluids: A Direct Numerical Simulation Study. Fluids 2021, 6, 40. https://0-doi-org.brum.beds.ac.uk/10.3390/fluids6010040

AMA Style

Bräuer F, Trautner E, Hasslberger J, Cifani P, Klein M. Turbulent Bubble-Laden Channel Flow of Power-Law Fluids: A Direct Numerical Simulation Study. Fluids. 2021; 6(1):40. https://0-doi-org.brum.beds.ac.uk/10.3390/fluids6010040

Chicago/Turabian Style

Bräuer, Felix, Elias Trautner, Josef Hasslberger, Paolo Cifani, and Markus Klein. 2021. "Turbulent Bubble-Laden Channel Flow of Power-Law Fluids: A Direct Numerical Simulation Study" Fluids 6, no. 1: 40. https://0-doi-org.brum.beds.ac.uk/10.3390/fluids6010040

Article Metrics

Back to TopTop