The effect of liquid co-flow on gas fractions, bubble velocities and chord lengths in bubbly flows. Part I: Uniform gas sparging and liquid co-flow

Abstract Unique experiments were performed in a homogeneously sparged rectangular 400 × 200 × 2630 mm ( W × D × H ) bubble column with and without liquid co-flow. Bubbles in the range 4-7 mm were produced by needle spargers, which resulted in a very uniform bubble size. Dual-tip optical fibre probes were used to measure horizontal profiles of gas fractions, bubble velocities and bubble chord lengths for superficial gas velocities U s g in the range 0.63-6.25 cm/s and superficial liquid velocities U s l up to 20 cm/s. Images of the bubble column were captured and a Bubble Image Velocimetry technique was adopted to calculate bubble (parcel) velocities. For low gas fractions, when a homogeneous flow regime occurred, both methods agreed very well and the optical fibre probes were found to be rather accurate for our bubbles. A liquid co-flow was found to have a calming effect and to stabilize a homogeneous bubbly flow regime, with less spatial variation in gas fractions and bubble velocities. Bubble chord lengths were almost normally distributed and do not exhibit the theoretical triangular probability density functions. The mean cord lengths were in the range 1.9-3.5 mm and found to increase with U s g and to decrease slightly with increasing U s l , while a liquid co-flow significantly reduced the standard deviation of the chord length distribution.


Introduction
Due to the continuous increase in computational power and demand for more accurate multiphase CFD simulations, there is an obvious need for more precise and detailed experimental data on bubbly flows for development and validation purposes. Euler-Euler CFD simulations, where both liquid and gas phases are modeled as interpenetrating fluids, require proper modeling of two-phase turbulence and of the interfacial forces such as drag, virtual mass, lift, wall lubrication and turbulent dispersion ( Van den Akker, 1998; Van den Akker, 1998;Dhotre et al., 2013;Liao et al., 2015;Van Den Akker, 2015 ). These sub-models dealing with interfacial momentum transfer rates and bubble induced turbulence, are a strong function of (local) bubble size, slip velocity and void fraction. Many CFD models (see e.g. Dhotre et al., 2013 ) assume a constant, single bubble size to keep the computational burden of the simulation Most recent developments of multi-phase CFD codes deal with breakup and coalescence kernels to more precisely model interfacial momentum transfer rates and bubble induced turbulence based on a local bubble size distribution (modeled by a population balance and a limited amount of bubble size classes), at the cost of increased complexity, longer simulation times and convergence issues. Huang et al. (2018) studied the impact of bubble size modeling in CFD simulations of bubble columns and found that: (1) the single bubble size models gave surprisingly good agreement with experimental data for symmetrical bubble columns, but less accurate agreement for a-symmetrically sparged configurations (which will be dealt with in Part II); and (2) there is no agreement on accurate models describing bubble coalescence and breakup rates (such as h-and i-MUSIG models). According to McClure et al. (2017) average bubble size and bubble size distribution have an effect on the swarm drag correction factor. Besagni et al. (2017) explored various approaches to represent the bubble size distribution and also concluded that the correct simulation of the fluid dynamics in the bubble column requires (more) accurate coalescence and break-up closures.
The rationale behind the research reported both in this Part I and in Part II is to acquire experimental data with bubble size as uniform as possible. Eventually, such data may be instrumental for distinguishing between effects caused by a swarm of uniform bubbles and the more complicated interactions between non-uniformly sized bubbles. In the latter case, more sophisticated models for coalescence and break-up are needed, see e.g. Mukherjee et al. (2019) . In the ideal case of a single unique and constant bubble size, accounting for effects of a bubble size distribution and for breakup/coalescence is no longer necessary, such that a model for bubble induced turbulence may be validated independently of the enactment of models for interfacial momentum transfer.
It is then obvious that there is a coexisting requirement for bubble size measurements along with local bubble velocities and gas fractions. In this paper, we report new experimental data, in terms of local void fractions, bubble velocities and bubble chord lengths, on bubbly flows in a homogeneously sparged bubble column with a very uniform bubble size. For diluted bubbly flows , or in shallow (pseudo 2D) bubble columns such as in Lau et al. (2013) , an image analysis approach can be performed to obtain gas fractions and bubble size and shape measurements. For denser bubbly flows and/or larger bubble columns, image analysis becomes increasingly difficult due to overlapping bubbles. Measurement methods are then limited to X-ray densiometry ( Hernandez-Alvarado et al., 2018;Mandalahalli et al., 2020 ), electrical resistance tomography ( Singh et al., 2017 ) or intrusive measurement methods such as a borescope ( Hernandez-Alvarado et al., 2018 ), wire-mesh sensors ( Prasser et al., 1998;Hampel et al., 2009;Hernandez-Alvarado et al., 2018 ), shadowgraphic optical probes ( Lichti and Bart, 2018 ), or (multi-point) electrical resistance ( Buwa and Ranade, 2005;Singh et al., 2017 ) or optical fibre probes ( Frijlink, 1987;Bakker, 1992;Harteveld, 2005 ).
Optical fibre or electrical resistance probes are regularly used in single tip, dual tip, or four-tip configurations.Probe signals are typically sampled at a rate in the order of O (10 4 ) Hz ( Chaumat et al., 2005;Le Corre et al., 2003;Shen and Nakamura, 2014 ), or sometimes even lower as long as the phase indicator function (PIF) can be registered with sufficient temporal resolution ( Kiambi et al., 2003 ). Fibre probes in a single tip configuration are used to obtain the local phase indicator function (PIF) for determining the local gas fraction ( Enrique Juliá et al., 2005 ) and the power spectral density of the PIF ( Singh et al., 2017;Tyagi and Buwa, 2017 ). Dual-tip probes (with a vertical spacing y ), can measure, in addition to the local PIF, bubble velocities (in the y -direction) based on the flying time and bubble diameters and chord lengths by assuming aligned uni-directional flow Frijlink, 1987;Bakker, 1992;Harteveld, 2005;Tyagi and Buwa, 2017;Groen, 2004;Simonnet et al., 2007;Dias et al., 20 0 0;Barrau et al., 1999;Chaumat et al., 2005;Kiambi et al., 2003;Murzyn et al., 2005 ). Four-point probes (in a triangular pyramid ( Guet et al., 20 03;20 05;Ojha and Dahhan, 2018;Bai et al., 2008;Lucas and Mishra, 2005;Xue et al., 2003;Shen and Nakamura, 2014 )) are adopted to also determine bubble velocity directions and shapes at the cost of a more complex algorithm for the probe signal analysis.
All these types of optical fibre probes are inherent to certain measurement inaccuracies and sampling bias caused by: (1) the blinding effect due to improper (de-)wetting of the probe; (2) the crawling effect as a result of deformation and/or deceleration of a pierced bubble; and (3) the drifting effect as bubble trajectories are altered due to the presence of an intrusive probe. The latter effect causes challenges for calculating correct flying times for dual-tip or four-point bubble probes as bubbles are deflected by the first (lower) probe tip and not measured by the upper probe tip(s).
While Cartellier (1992) correlated the rise time (or signal derivative as in Mizushima et al. (2013) ) of a signal from a single optical fibre with interface velocity measurements from a digital camera, Cartellier and Barrau (1998a) used conical fiber tips and Frijlink (1987) ; Groen (2004) ; Cartellier and Barrau (1998b) introduced improved fibre tip shapes (as adopted later by Pjontek et al. (2014) ) to determine the interface velocity. However, these authors still recommend calibration of the velocity measurements (using piercing experiments) for each manufactured tip or fluid.
More recent developments on velocity measurements using a single fibre are based on resolving the coherent beat frequency between the Fresnel reflection (fibre-fluid interface) and the reflections of an approaching interface by using a very high sampling rate (10 MHz) as in Chang et al. (2003) ; Lim et al. (2008) . As this technique resolves the velocity of an approaching interface before bubble piercing takes place (when the distance between the fibre tip and bubble surface is 10 0-30 0 μm), effects of blinding, crawling and drifting are eliminated. However, as the intensity of the scattered light is limiting, velocity realizations are only found possible if the angle of attack is almost normal. Consequently, only as little as 2% of the detected bubbles contain velocity information (when pierced in the center of a bubble). While the velocity measurements may be still very representative, chord length measurements are largely biased for both single tip and dual tip optical fibres.
We now report new bubbly flow data acquired in the "Lim-BuRig" test facility, which is meticulously described in our previous paper ( Muilwijk and Van den Akker, 2019 ). This setup, with its rectangular geometry for visualization purposes, has been developed with the view of generating experimental data for CFD validation and development. We opted for a rectangular configuration to facilitate visual observations and optical measurements since anyhow the test facility is intended to provide data for CFD validation rather than to mimic cases of direct industrial relevance. With the two separate inlet compartments, configurations with (a-)symmetric liquid and/or gas flow distributions can be studied and a wide range of operating conditions can be investigated; see also our Part II paper. A needle sparger was carefully constructed, such that a maximally uniform initial bubble size distribution of large bubbles in the range 4-7 mm (with an almost constant terminal rise velocity ( Clift et al., 1978 )) was achieved for U sg up to 3.1 cm/s without liquid co-flow and beyond for higher co-flow velocities. In our previous study ( Muilwijk and Van den Akker, 2019 ), a correlation was developed to describe the initial (volume equivalent) bubble size d b,eq (at gas sparger level) as a function of the superficial gas and liquid velocities. This correlation is used to calcu- Table 1 Operating conditions. OFP: Optical Fibre Probe; BIV: Bubble Image Velocimetry. x : lateral position; y : vertical position (from trailing edge of the splitter plate). Note that the trailing edge of the splitter plate is located 17 cm above the sparger. Photographs, see Fig. 5 , bubble velocity and chord length distributions (OPF measurements in the center of the column at x = 0 cm), see Figs. 9 , and 14 respectively, are shown for four cases near the limits of the operating conditions ( U sg = 1.25, 6.25 cm/s and U sl = 0, 20 cm/s).

Method
U sg cm/s U sl cm/s x cm y cm Fig. OFPs 0. 10,63 Figs. 7,8,10 ,11,12 ,15 ,16 ( α,v b ,c) (in steps of 2.5 cm) 10, 20 −20 . . . 20 −15 . . . 1250 Fig. 10 b late the bubble sizes in the present study. For bubbles in this size range, the terminal rise velocity (in contaminated water) is estimated as ≈ 24 cm/s according to the parameterization developed by Park et al. (2017) . We developed our in-house dual-tip optical fibre probes which were used to measure gas fractions, bubble velocities and bubble chord lengths. Bubble velocities measured with dual-tip optical fibre probes were then compared with bubble velocities measured using a Bubble Image Velocimetry ( Cheng et al., 2005 ) approach at low gas fractions. Unlike other reported data on bubbly flow in the "pseudo-homogeneous" (poly-dispersed homogeneous) regime ( Besagni et al., 2018 ), we report chord length distribution measurements for conditions with known ( Muilwijk and Van den Akker, 2019 ) uniform single sized bubbles for various U sg and U sl .
The main objectives of this paper are to show how: (1) the (uniformity of) the gas fraction α is a function of the superficial liquid and gas velocities and to validate the previously proposed correlation for the gas hold-up; (2) the bubble velocity v b is influenced by liquid co-flow, where the optical fibre probe (OFP) and Bubble Image Velocimetry (BIV) methods are compared; and (3) the distribution of the bubble chord length c is a function of the superficial gas velocity and liquid co-flow velocity and how they relate to the mean bubble diameter. The structure of this paper is then as follows. Section 2 briefly describes the setup and adopted measurement methods; Section 3 shows results on the measured void fractions, bubble velocities and chord lengths respectively; and conclusions and recommendations for future work are presented in Section 4 . Fig. 1 shows a schematic of the apparatus used for this study. Two parallel streams of bubbly flows, for this work with equal superficial gas, U sg , and superficial liquid velocities, U sl , start interacting downstream of the trailing edge of a splitter plate. The depth of the channel D is 0.2 m and the width W of the channel is 2 D = 0 . 4 m. The superficial velocities U sg and U sl are the independent variables, defined as volumetric flow rates (at P 0 = 1 . 013 bar and T = 20 • C) divided by the cross-sectional area of the column. U sg was varied in the range 0.63-6.25 cm/s, and U sl in the range 0-20 cm/s. Bubbles are formed by 2 × 14 2 = 392 ∅ 1 . 55 mm i.d. needles to establish formation of a uniform homogeneously distributed bubble size. More details on the design of the test facility, overall gas fractions (using a bed expansion technique) and bubble formation rates/diameters can be found in our previous paper ( Muilwijk and Van den Akker, 2019 ). Optical Fibre Probes and a Bubble Image Velocimetry technique are used to study the flow. These techniques are discussed below and the experimental settings are summarized in Table 1 .

Bubble image velocimetry
Images of the bubble column (of the area annotated with BIV in Fig. 1 ) were captured at a rate of 100 and 120 Hz (Jai Go 2400 M camera, Kowa LMVZ166HC 16-64 mm varifocal lens) and corrected for lens distortion, see Fig. 2 a. The camera was calibrated using both the width and height of the column and resulted in a spatial resolution of 0.70 mm/pix. A contrast limited adaptive histogram equalization (CLAHE) algorithm ( Zuiderveld, 1994 ) was used to improve the contrast of the images as shown in Fig. 2 b. A direct image cross-correlation technique as in Cheng et al. (2005) was adopted to obtain the velocity of bubble parcels.
The image cross-correlation coefficient between f (time t) and g (time t + t) was calculated according to where n, m are sub-ranges (interrogation windows) of the full image and i, j are the pixel shifts in vertical and horizontal direction. Fig. 3 shows a surface plot of R i, j as function of i and j. A window size of 32 × 32 pixels (120 Hz image acquisition rate) or 40 × 40 pixels (100 Hz) was found as an optimal compromise between spatial resolution and correlation intensity. A window overlap of 50% was used, which resulted in a spatial resolution (half window size) of ≈ 1 . 25 -1.5 cm. Quadratic interpolation was used to obtain a sub-pixel displacement resolution and spurious velocities were removed using a mean ±3 σ outlier detection algorithm.

Dual-tip optical fibre probes
Horizontal profiles of gas fraction α, bubble velocity v b , and bubble chord length c were measured in the center of the column (between front and rear wall, see Fig. 1 ) at 23 and 63 cm downstream of the trailing edge of the splitter plate (40 and 80 cm above the sparger level) by using two in-house developed doublepoint optical fiber probes.
The probe response signals from each probe tip were sampled with a NI USB-6001 I/O device at a rate of 5 kHz per channel, and the signals were normalized using the span of a signal (see black and grey curve in Fig. 4 ). It was found that the adopted acquisition rate is sufficient to record the phase indicator functions, as multiple points were sampled on the steep parts of the signal (when tip (de-)wetting takes place). In order to determine the threshold intersections with a higher temporal resolution, the signals were upsampled at a rate of 20 kHz using quadratic interpolation and binarized to obtain the phase indicator functions, using a threshold value of 3 × the standard deviation of the baseline noise (red dashed line).
Void fractions were then calculated using the average value of the phase indicator functions for the lower probe tips, as those are less subjected to the drifting effect. A measurement duration of 100 s was found to be sufficient to obtain an accurate estimate of the local void fraction.  Table 1 . The direct cross correlation between the phase indicator functions from the upper and lower fibre tips was computed to obtain the time lag τ max between the signals. Then, for each bubble detected on the lower fibre tip at time t i ( ), a matching crenel is searched on the upper fibre tip signal within a time interval [ t i + aτ max ; t i + bτ max ] as in Chaumat et al. (2005) . The value of a was set as 0.5 and b was chosen in the range 2.5-6, dependent of the gas fraction.
The velocity of a bubble i traveling through both lower and upper tips, is then calculated according to where the flying time t i f , is measured as the time the i th bubble front interface takes to rise from the lower to the upper fibre tip as indicated in Fig. 4 a.
The vertical distance between the two tips of a probe, y, was 3.75 mm for probe No. 1 and 2.45 mm for Probe No. 2. As we expect uniform large ( d eq > 4 mm) oblate/wobbling bubbles rising mostly vertically, the probe tip separations y and sampling rate were optimized to deliver accurate results and no significant difference was found between the two dual-tip probes. For our application, bubble velocities > 1 . 0 m/s are highly unlikely and a sufficient resolution in the velocity domain is achieved with our optical fibre probes.
The dwelling time t i d is the duration a probe tip spends inside a bubble (time between the piercing of the front and rear interface of a bubble, see Fig. 4 ). Fig. 4 where the subscripts l, u refer to the lower and upper fibre tip, respectively. Under ideal probe conditions (rigid ellipsoidal bubbles, no crawling, no drifting), it is assumed that this parameter exhibits a uni-modal distribution for bubbles traveling with a lateral velocity component. Hence, a value of 1.2 was adopted for the pairing criterion (see dashed lines in Fig. 4 b as the probability density function estimate of the pairing parameter was found to have a minimum at approximately 1.2. Bubble pairs not passing the pairing criterion, with an absolute relative dwell time difference exceeding the value of 1.2, were ascribed to extreme drifting and crawling. Barrau et al. (1999) studied the effect of relaxing the pairing criterion on the resulting phase averaged gas velocity and found no significant difference when the pairing rate (percentage of bubbles detected on the lower tip matched to a crenel on the upper tip signal) exceeded 20%. Fig. B.19 shows the pairing rates for our experiments for further reference.
The phase weighed mean bubble velocity (phasic gas velocity) is then calculated as: where N is the validated amount of bubbles measured on both fibre tips passing the pairing criterion.
For each valid bubble velocity measurement, the chord length of a bubble is then calculated according to: The mean and standard deviation of the chord length was then calculated from the chord length distribution and denoted as c and Stdev (c) The mean chord length can also be calculated according to Chaumat et al. (2005) ; Lim et al. (2008) where v mp is the most probable velocity, calculated as y/τ max and f b is the bubble detection rate on the lower fiber tip, which can be determined with high accuracy. Both methods for determining c will be compared in Section 3 . For spherical and (oblate) ellipsoidal bubbles, the mean vertical diameter of the bubbles d can be calculated according to Simonnet et al. (2007) ; Colombet et al. (2015) : Assuming uniformly sized bubbles, the volume equivalent bubble diameter then relates to the aspect ratio and vertical diameter of the bubbles according to  ; Simonnet et al. (2007) ; Colombet et al. (2015) : where the aspect ratio of the bubble ϕ = d /d ⊥ , with d ⊥ and d being the major and minor axes of an oblate ellipsoidal bubble. At low aeration rates (left), the gas fraction is low (5.5% (a) and 3% (c)) and the column is sufficiently optically accessible to obtain a good contrast for image analysis. Individual bubbles can be distinguished and a very uniform bubble size distribution is observed. A BIV technique may be used in this regime due to the justified assumption of uniformity of the flow (the absence of wall/center peaking in velocity and/or void fraction), while the effects of front and rear wall on the measured velocity profiles are ignored.

Visual observations
At intermediate to high gas fractions, a boundary layer develops at both sides of the splitter plate, thereby creating a wake region in the center of the bubble column, leading to horizontal gradients of the void fraction and bubble velocities in the center of the column.
At high aeration rates (right) the fluid is opaque, α = 24% (b) and 14% (d), and only bubbles in the front wall region can be recorded on a camera. Bubbles do not follow organized rectilinear bubble paths (see also Fig. 6 ), and a three dimensional turmoil develops as clearly visible in the videos in the Supplementary Material. Back-mixing occurs, while a center peaking void fraction/bubble velocity and flow reversal at the column walls is emerging. Under these conditions, a BIV technique cannot be used anymore as the assumption of quasi-2D flow (no gradient in the collinear direction) no longer holds.
Without liquid co-flow and high aeration rates ( Fig. 5 b), the uniformity of the bubble size distribution suffers due to coalescence occurring during bubble formation at the needles (see sparger region). With increasing liquid co-flow ( Fig. 5 d), coalescence (e.g. at the needle) was prevented (liquid co-flow reduces bubble-bubble interactions ( Muilwijk and Van den Akker, 2019a )) 1 Note that the factor 3 2 in Eq. (6)   and a uniform bubble size distribution was reestablished. Parallel trains of bubbles were formed in the sparger region, which started to interact some 5-10 cm above the needle outlets, also leading to the development turbulent structures and backmixing. Fig. 6 shows simulated bubble streaks by averaging the 500 photographs of a 5 s image series. Under typical conditions our column is free from lateral dispersion, at least in the lower part of the column as the rise velocities of our bubbles are very uniform thanks to the already limited variation in bubble size. At low superficial gas velocities ( Fig. 6 a), bubble streaks originating from a single needle can be distinguished up to a height of ≈ 20 cm above the needle outlets, while a wake effect of the splitter plate is clearly visible in the center of the column. For somewhat higher U sg ( Fig. 6 b, see also Fig. 5 a), lateral dispersion is emerging and separate bubble streaks cannot be observed beyond ≈ 5 -10 cm above the needle outlets. A coflow ( Fig. 6 c) was then found to reduce lateral dispersion due to a strong advection as straight bubble streaks can be detected up to ≈ 40 above the needle outlets, see also    Fig. 7 a shows lateral profiles of the void fraction for U sg in the range 0.63-6.25 cm/s without co-flow at a height of 80 cm above the sparger level. For low gas flow rates, the gas fraction profile was very uniform, whereas for increasing gas flow rates, steep gradients of α were found in the close vicinity of the column wall.

Horizontal gas fraction profiles
Due to wall effects (including splitter plate), bubbles migrated to the center of the bubble column, thereby creating local maxima in the void fraction profile and liquid downflow in the direct vicinity of the column walls. Fig. 7 b shows lateral profiles of the void fraction for U sl in the range 0-0.2 m/s for U sg = 1 . 25 cm/s (note the difference in scale between (a) and (b)). While a slightly uneven void fraction profile was observed for U sl = 0 m/s (grey ), a liquid co-flow flattens the void fraction profile and the effect of the splitter plate and wall effects on the lateral bubble migration were reduced.   Fig. 7 . The filled markers denote measurements at a height of 80 cm above the sparger level, while the open markers show measurements at a height of 40 cm. Measurements at U sg = 1.25, 3.13 and 6.25 were repeated thrice and were found be reproducible. For higher U sg and for U sl = 0 m/s, lower void fractions were measured at y = 40 cm compared to the measurements at y = 80 cm. This can be explained by a more homogeneous flow pattern closer to the sparger, in the developing region of the bubble column, as a center-peaking void fraction profile with liquid downflow at the walls is emerging.

Gas hold-up curves
It should be noted that U sg on the x -axis is given as the setting at standard conditions while the actual superficial gas velocity is a function of the height in the column. More data is required in terms of local pressure measurements to derive the local gas flux as a function of the height in the column.
The dotted lines in Fig. 8 show the overall gas hold-up as predicted by the correlation developed in our previous study ( Muilwijk and Van den Akker, 2019 ), which was obtained through bed expansion experiments. Good agreement was achieved with void fraction measurements using optical fibre probes (markers) at a height of 80 cm above the sparger. A comprehensive comparison of our gas hold-up curves with data found in the literature is given in our previous paper ( Muilwijk and Van den Akker, 2019 ). It was found that gas hold-up curves are highly case-specific, which is also clearly illustrated by a thorough overview of gas hold-up curves given in Gandhi and Joshi (2010) . As such, it can be concluded that gas hold-up curves are a strong function of the bubble size and sparging method (initial bubble size distribution). Many drift-flux correlations, for estimating the gas hold-up, ignore the effect of the bubble size distribution and may thus give inaccurate predictions. Fig. 9 shows velocity histograms for four operating conditions in the range studied, e.g. U sg = 1 . 25 (left) and 6.25 cm/s (right) and co-flow velocities U sl = 0 (top) and 0.2 m/s (bottom). For all cases, the velocity histograms exhibited an almost Gaussian distribution. As bubbles were almost uniform in size and had an almost equal rise velocity in this size range ( Clift et al., 1978 ), the spreading of the measured velocities can be explained by: (1) swarm effects, as bubbles may be hindered or accelerated due to the presence of other bubbles, thereby developing large circulation patterns; (2) interface oscillations, as bubbles of this size behave as wobbling bubbles (instead of rigid ellipsoids); (3) crawling, as the velocity of an interface may be influenced by the presence of the probe and may also depend on the bubble velocity itself, interface curvature (phase of oscillation) and piercing position; (4) coalescence at gas sparger level, when the uniformity of the initial bubble size is compromised (only for high U sg and no liquid co-flow, see Fig. 9 b). For increasing liquid co-flow the spread in bubble velocities nar- rows as bubbles carry more momentum (in their added mass) and are less likely to suffer from crawling and swarming. for U sg in the range 0.63-6.25 cm/s without co-flow at a height of 80 cm above the sparger level. Low gas flow rates resulted in slightly wavy bubble velocity profiles with bubble velocities in the range 19-24 cm/s. For increasing gas flow rates, bubbles traversed away from the column side walls (and splitter plate), thereby creating two maxima in the bubble velocity profiles. Steep gradients of v b show up in the vicinity of the column wall with bubbles even moving downward close to the column side walls, which cannot be properly detected using optical fibre probes in the current configuration.

Horizontal profiles of the bubble velocity
Standard deviations of the bubble velocity (see Fig. 10 c) are almost constant across the lateral direction and increase with increasing U sg . Slight increases of Stdev (v b ) are noticed close to the column wall for high U sg due to the emerging, unsteady, down flux at the column walls. Fig. 10 b shows lateral profiles of v b obtained using OFPs (large markers + dotted line) and our BIV technique (small markers + solid line) for U sl in the range 0 − 0 . 2 m/s and U sg = 1 . 25 cm/s. The shape of the velocity profiles, obtained by using the two methods, is comparable for all cases, while for U sl = 0 . 1 m/s the agreement is almost perfect.
Due to wake effects of the splitter plate, a double peaking bubble velocity profile develops for U sl = 0 m/s (grey ). A liquid coflow then has a uniforming effect on the bubble velocity profile as (1) the effect of the splitter plate and wall effects on the lateral bubble migration is reduced as bubbles are entrained by a (uniform) liquid flow and (2) rather lower standard deviations (see Fig. 10  For U sl = 0 m/s (higher α), BIV showed lower bubble velocities than the OFPs, which may be due to a gentle down flux reducing the bubble velocities in the vicinity of the column walls. For U sl = 0 . 2 m/s, velocities from BIV were ≈9% higher than mean velocities obtained by using the dual-tip optical fibre probes. BIV experiments at U sg = 1 . 25 cm/s were repeated twice. Whereas images were initially acquired at a frame rate of 100 Hz, a different set of experiments was carried out with a framerate of 120 Hz. No significant differences were found between the experiments for U sl > 0 m/s, as the flow patterns without co-flow were found to be very sensitive to small inaccuracies of the gas flow rates in the left and right inlet, which is further studied in our Part II paper. The discrepancies between OFP and BIV may be due to 3-D effects (wall peaking of the gas fraction and bubble velocity profiles). Shen et al. (2016Shen et al. ( , 2017 ; Sun et al. (2014) observed wall peaking of the void fraction in the near vicinity ( < 10 − 20 mm distance) of the column walls as a function of the gas flux. This near-wall regions are outside the reach of the bubble probes. A future investigation, with a different probe design, should be dedicated to zoom-in on the phenomena happening in the near vicinity of the side or front walls. A wall peaking void fraction may explain the higher bubble velocities from the BIV captured in the near-wall region. Experiments with a smaller depth of view, higher frame rate and improved homogeneous illumination are recommended to optimize the BIV technique.

Bubble velocities as a function of U sg
We now average the mean bubble velocities and standard deviations measured on the 15 lateral positions as shown in Fig. 1 b. Fig. 11 a shows plots of this center-line averaged mean bubble velocity v b W as a function of superficial gas velocity. The error bars in Fig. 11 a are a measure of the uniformity of the profiles of v b (in Fig. 10 a,b) as a result of center peaking velocity profiles. Fig. 11  Without liquid co-flow ( •, •), average mean bubble velocities increase with increasing superficial gas velocity up to U sg ≈ 3 . 0 cm/s. An almost constant average mean bubble velocity was measured in the range U sg = 3 . 0 − 4 . 5 cm/s, whereafter the bubble velocities rapidly increase with increasing U sg .
The standard deviation of the bubble velocities continuously increases with increasing superficial gas velocity for the whole range of U sg . Higher bubble velocities and standard deviations were measured at a height of y = 80 cm ( •), compared to y = 40 cm ( •), as the flow regime gradually departed from a homogeneous bubbly flow and center peaking void fraction and bubble velocity profiles started to develop. While at the lowest aeration rates, bubbles in the range 4 − 8 mm essentially have a constant rise velocity of ≈ 24 cm/s ( Clift et al., 1978 ), a standard deviation of ≈ 5 cm/s was Also, co-flow leads to lower standard deviations of v b (see Fig. 11

Gas flux measurements
The local gas flux, as measured by the optical fibre probes, can be calculated as αv b , where the phasic quantity v b is the phase averaged bubble velocity as per Eq.
(3) . The centerline averaged gas flux J g W , along the accessible positions along the column center line, is then given by: and shown in Fig. 12 . The brackets denote spatial averaging of the locations shown in Fig. 1 b. The closest position the probes can travel to the wall is 2.5 cm, hence, x L = −17 . 5 and x R = 17 . 5 cm. Under ideal circumstances, in a uniform homogeneous bubbly flow, the measured J g W equals the U sg as calculated from the gas flow rate settings. Bai et al. (2010) ; Shen et al. (2016) used J g to assess the accuracy of the bubble probes, whereas ( Colombet et al., 2015 ) adopted J g to show the degree of homogeneity of the bubbly flow. While Bai et al. (2010) found their J g being some 30% smaller than the applied gas flow rate (for 2 < U sg < 10 cm/s), we found quite good agreement between both values for U sg up to about 2 cm/s without co-flow. At higher U sg , the measured gas flux departs from the applied U sg which can be explained by the center peaking void fraction and bubble velocity profiles (with down flow emerging at the column walls). This divergence starts at U sg = 2 . 0 cm/s and becomes increasingly distinct for U sg > 4 . 5 cm/s, while this effect is less prominent for the region closer to the sparger where the flow pattern is still developing (open •).
The change in hydrostatic pressure between the two measurement planes, at y = 40 and y = 80 cm ( H = 0 . 4 m), can be approximated by ρ w g H(1 − α) . This was found to be less than 4% of the atmospheric pressure, hence the difference in local superficial gas velocity between y = 40 and y = 80 cm does not explain the difference in j g W we measure between y = 40 and y = 80 cm for U sg > 3 cm/s. No significant difference of the pairing rate was found between the measurements taken at y = 40 cm and y = 80 cm, see also Fig. B.19 , hence the bubble probes at both heights work equally well.
Our measured gas fluxes J g W , for U sl = 0 m/s at y = 80 cm, are higher as (1) more measurements were performed in the bulk of the column (the plateau in the void fraction and velocity profiles, see Figs. 7 a and 10 a respectively); and (2) very low/negative bubble velocities in the vicinity of the column wall cannot be measured due to the limitations of the probe.
The measurements at y = 40 cm were taken closer to the sparger, still in the developing region of the bubble column, and follow the parity line due to a still rather uniform bubbly flow. In the developed region of the bubble column, liquid downflow at the side walls emerged, leading to a center-peaking gas fraction/bubble velocity profile, hence, J g W departs from this parity line. Considering data valuable for CFD validation, centreline profiles would give systematic and valuable information on the bubble flow behavior.
It should be noted that U sg on the x −axis is at standard conditions, while the measured J g is at a higher pressure. Therefore, a somewhat lower measured gas flux J g W can be expected due to the higher gas density. Measurements of the local pressure or gas fraction distribution along the height of the column are required to properly account for the change in gas fraction difference. In case a liquid co-flow is applied, the gas flux is almost equal to the applied U sg for all cases considered. It may then be concluded that (1) the dual-tip bubble probes work very well for determining α and v b ; (2) the flow regime gradually departs from homogeneous bubbly flow at U sg > 2 . 0 cm/s and U sl = 0 m/s; and (3) a liquid co-flow stabilizes the homogeneous bubbly flow regime. Fig. 13 shows the relative bubble (slip) velocity as a function of the gas fraction, where U s is calculated assuming uniform (co-)flow by using a drift-flux model according to Muilwijk and Van den Akker (2019) :

Literature comparison
The black markers represent the data as shown in Figs. 8 and 11 a, whereas the colored markers and dashed line show relative gas velocities as a function of α as reported in a selection of relevant literature.
While the relative gas velocities reported by Colombet et al. (2015) ; Alméras et al. (2018) ;  (no co-flow) and Garnier et al. (2002) ; Simonnet et al. (2007) (with co-flow) show decreasing gas velocities with increasing gas fraction at low gas fractions, we found U s to increase with increasing α for all conditions ( α > 3%) investigated. As the estimated terminal rise velocity U t for isolated bubbles in the range 4-8 mm is ≈ 24 cm/s, U s at low α was slightly smaller than U t . Hindered rise (as in Besagni et al., 2016 up to α ≈5% andSimonnet et al. (2007) up to α ≈ 15%) may occur to a some extent for lower gas fractions, but no data could be obtained in that regime due to the risk of weeping. The early onset of swarming (increasing U s with increasing α) can thus be explained by the large bubble size in agreement with Simonnet et al. (2007) , as small (uniform) bubbles, which stabilize homogeneous bubbly flow, show hindered bubble rise behavior as in Garnier et al. (2002) ; Alméras et al. (2018) ; Colombet et al. (2015) .
Relative velocities measured with the optical probes are higher than derived from global gas hold-up measurements by using a drift-flux model as in Muilwijk and Van den Akker (2019) (see cases B1, C2 in Fig. 13 ). At low α (uniform bubbly flow), good agreement was found, but, due to a center peaking gas fraction/bubble velocity profile at higher void fractions, the measured gas flux is overestimated (see Fig. 12 ). U s is then also overestimated as the probes are biased to high (upward) velocities and low (downward) velocities in the vicinity of the column wall cannot be measured by the current probe configuration.
It was observed that when a liquid co-flow was applied, it reduces the emerging down-flow at the wall as the liquid entrained by the bubble wakes can overflow from the top of the column, thereby resulting in a more homogeneous flow. A critical value of U sl may exist when the liquid flow rate ( U sl ) equals the liquid entrainment rate in the bubble wakes, which is roughly estimated at C AM U sg , with C AM being the added mass coefficient. More experiments, at lower co-flow rates will be required to investigate this hypothesis.

Histograms of bubble chord length distributions
Fig. 14 shows chord length distributions for four operating conditions in the studied range, e.g. U sg = 1 . 25 (left) and 6.25 cm/s (right) and co-flow velocities U sl = 0 (top) and 0.2 m/s (bottom). Bubble chord lengths clearly increase with increasing gas flow rates (left to right), while mean chord lengths become slightly smaller and distributions more narrow with increasing co-flow velocity (top to bottom). A small shoulder shows up at the left tail for U sl = 0 m/s, which is not present for chord length distributions obtained for (higher) co-flow velocities, which agrees well with the bubble contours presented in Muilwijk and Van den Akker (2019a) , indicating that liquid co-flow tends to make the bubble more spherical. Roig et al. (1998) found a similar chord length distribution as shown in Fig. 14 a, while bubbles were produced by porous tube spargers, which typically results in a polydisperse bubble size distribution. Constant bubble chord length distributions were found across the lateral and streamwise position and they concluded that "either break-up and coalescence did not occur in the flows, or they were in mutual equilibrium." While their chord length distribution was found to be independent of the superficial liquid velocity, where c = 2 . 1 ± 0 . 6 mm for α = 1.9% and U sl = 0 . 2 -0.5 m/s, our results presented in Fig. 14 show a clear dependence of c on U sl and U sg as per the design of the needle sparger ( Muilwijk and Van den Akker, 2019 ).
A theoretically derived triangular left skewed chord length distribution for rigid ellipsoidal (uniform) bubbles ( Liu and Clark, 1995;Clark and Turton, 1988 ), with the largest measured chord length being the maximum (vertical) diameter d , was not observed in our cases. Also Chaumat et al. (2005) reported that 15% of the measured chord length values were larger than the maximum diameter obtained from image analysis. Although our previous paper ( Muilwijk and Van den Akker, 2019 ) shows that very uniform bubbles are formed at gas sparger level for U sg up to 3 cm/s (and beyond when a liquid co-flow is applied), theoretically predicted triangular chord length histograms were not recovered. This discrepancy is explained by: (1) (2019) ) of non-rigid bubbles, which also explains the longer right tail for larger bubbles (see Fig. 14 c).
Theoretical chord length distributions were simulated (see A.1 ) for mono-disperse rigid and wobbling bubbles, as well as for a poly-disperse bubble mixture and we found very similar chord length distributions for mono-disperse wobbling bubbles and polydisperse rigid (oblate) bubbles. Therefore, we refrain from performing a transformation from a chord length distribution to a bubble size distribution as in Hoang et al. (2015) as assuming a constant ellipsoidal shape renders such a transformation invalid. A more complex method should be developed to account for a variable bubble shape when transforming chord length distributions into a bubble size distributions as a constant oblate ellipsoidal shape is proven invalid for large bubbles. Fig. 15 shows lateral profiles of the mean chord lengths c (top) and standard deviations Stdev( c) of the chord length distribution at a height of 80 cm above the sparger level for (a,c) U sg in the range 0.63-6.25 cm/s without co-flow; and (b,d) U sl in the range 0-0.2 m/s for U sg = 1 . 25 cm/s (note the difference in scale between the left and right column). While the void fraction and bubble velocity profiles are clearly non-uniform for some of these cases (see Figs. 7 and 10 ), a very constant mean chord length c was observed for almost all U sg and U sl , which indicates that the bubble size (distribution) is very constant over the cross section of the column. A slight spanwise spreading was found for lower gas fractions, as smaller bubbles are more subjective to drifting and the number of valid bubble measurements decreased. Fig. 16 a plots the centerline averaged mean bubble chord length c W as a function of the applied superficial gas velocity and No significant difference was found between chord lengths measured at y = 80 cm and y = 40 cm (latter not shown), which indicates that the bubble size (distribution) is independent of the vertical position. The averaged mean chord lengths c (a) and the standard deviation thereof (b) increase with increasing U sg for cases with and without liquid co-flow. Whereas the lateral spread (error bars) of both c W and Stdev (c) W is almost independent of U sg , these distributions become narrower with increasing co-flow velocity, which confirms that liquid co-flow has an organizing effect on the (homogeneous) bubbly flow and the bubble formation process.

Bubble chord lengths as a function of U sg
The open markers ( Fig. 16 a) show the average chord length c as calculated by Eq. (5) . Compared to the mean of the chord length distributions (solid markers, Eq. (4) ), these chord lengths are somewhat overestimated as they depend on the direct cross correlation between the lower and upper fibre tip signal, which are biased to the larger (vertical) velocities. While α and f b can be obtained with a high degree of accuracy ( Chaumat et al., 2005 ), it should be noted that (mean) chord length measurements are subject to a similar relative error as bubble velocity measurements, see Eqs.
A clear bend in the slope of c W is observed at U sg ≈ 2 . 5 -3 cm/s (without liquid co-flow), which coincides with the apparent transition from single separate bubble formation to bubble formation with coalescence at the needle. As smaller satellite bubbles may split off from bubbles formed at the needle sparger and a uniform single bubble size cannot be assumed beyond U sg ≈ 3 . 1 cm/s ( Muilwijk and Van den Akker, 2019 ), a smaller chord length is expected, which agrees well with the trend in Fig. 16 a and the chord length histogram in Fig. 14 b. While still assuming a uniform volume equivalent bubble size and an average oblate ellipsoidal shape, the mean aspect ratio of a bubble may be calculated using Eqs. (6) and (7) . d b,eq was then calculated as a function of U sg and U sl using a correlation developed in our previous study ( Muilwijk and Van den Akker, 2019 ) and the volume equivalent bubble diameter was found to be in the range of 4-8 mm ( 2 . 2 < Eo < 8 . 7 ).
The calculated bubble aspect ratios ϕ were in the range 0.4-0.55 (see A.2 ). This agrees quite well with the values ( ≈ 0 . 5 ) re-ported in Garnier et al. (2002) ; Riboux et al. (2010) for the smaller bubble size. The aspect ratios for the larger bubbles obtained in the present study are slightly lower than reported in Colombet et al. (2015) (smaller bubbles) and Ziegenhein and Lucas (2017) ( 0 . 5 < ϕ < 0 . 6 for Eo > 2 ), while ϕ was found to decrease slightly with increasing bubble diameter for Eo > 2 . More experiments using close-up photographs of the bubble column are needed to further study bubble sizes and shapes, which also requires a very complex image analysis algorithm due to the high density of bubbles in the column.

Conclusions
New experiments were performed in a homogeneously sparged rectangular bubble column with large, almost uniformly sized bubbles operated with and without liquid co-flow.
Very uniform void fraction profiles were obtained for U sg up to 3.0 cm/s at void fractions up to ≈14%. A gradual transition to inhomogeneous bubbly flow was observed for larger U sg . Average gas fraction measurements, with and without liquid co-flow, agree very well with the correlation for the overall gas hold-up developed in our previous paper ( Muilwijk and Van den Akker, 2019 ).
Bubble velocities were measured using dual-tip optical fibre probes and compared with parcel velocities obtained using a Bubble Image Velocimetry approach. Very good agreement between both methods was found for gas fractions < 5% and U sl up to 0.1 m/s, whereas a < 10% higher velocity was found from BIV measurements compared to OFP measurements for higher U sl . While the biases of (dual-tip) optical fibres are well-known, further research is required to calibrate the BIV method. As the optical fibres are centered in the column between the front and rear wall, while the images for the BIV method are captured of the front column wall, gradients of v b in co-linear direction may impede a fair comparison of the two methods and further study is required to investigate the possibility of wall peaking gas fractions and bubble velocities.
A liquid co-flow was found to reduce spatial variations of v b and mitigate the wake effect of the splitter plate, therefore stabilizing a homogeneous bubbly flow. Bubble rise velocities as a function of gas fraction were compared to similar studies and hindered bubble rise was not observed in any case due to the large(r) size of the bubbles.
Very similar bubble chord length distributions were measured along the horizontal position in the column which confirms the uniformity of the bubble size and homogeneous sparging. Mean bubble chord lengths were found to increase with increasing U sg and to decrease with increasing U sl in a similar way as d eq depends on U sg and U sl .
As bubbles do not behave as rigid ellipsoids and intrusive measurements are subjected to measurement biases, an improved method should be developed to distinguish between bubble size (chord) and shape more accurately.
For this Part I paper, we reported accurate data for the effect of liquid co-flow on gas hold-up, bubble velocities and chord lengths, for bubbly flows characterized by large, uniformly sized, bubbles in the range 4-7 mm. Such data may be useful for validating CFD simulations, specifically the Euler-Euler (two-fluid) type, as models for bubble coalescence/breakup and segregated size classes are no longer required. Hence, models for interfacial momentum transfer mechanisms can be validated for a truly uniform bubble size. Future work then may include to study the effect of a bi-modal bubble size distribution, such that terminal rise velocities are no longer similar. Models for lateral dispersion can be properly validated based on empirical knowledge of the bubble size distribution, and depending on the interfacial tension (water type), coalescence could start to play a role.
Our Part II paper will deal with flow patterns in an asymmetrically sparged bubble column with an uneven co-flow by a parametric study using the same dual-tip optical fibre probes and Bubble Image Velocimetry.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. where E, D denote the distributions of aspect ratio ϕ and the bubble equivalent diameter d b,eq respectively. The variance in E (No. 3-5) was kept constant.
Then, for each bubble shape in E (No. 1,3-5) or size in D b,eq (No. 2), two large sets of random numbers (for the two horizontal directions R x , R y ) proportional to the horizontal diameter ( a 2 i ) were drawn from a uniform distribution, U (−a i , a i ) . The chord length distribution for each subset i was then calculated according to: where R = R 2 x + R 2 y . a i , b i are calculated using d b,eq and ϕ = b i /a i , see Eq. (7) . Fig. A.17 shows the simulated probability density function estimates for the five cases. For uniformly sized rigid bubbles with an  aspect ratio of ϕ = 0 . 5 , a triangular chord length distribution was recovered (black solid line) with a modal chord length (almost) equal to the maximum chord length ( Clark and Turton, 1988 ). The maximum chord length observed is 3.15 mm, which agrees with d b,eq = 5 mm and E = 0 . 5 , see Eq. (7) .
The gray solid line shows the simulated chord length histogram for univariate rigid (geometrically similar, ϕ = 0 . 5 ) bubbles in D b,eq of an average diameter d b,eq = 5 mm and a Coefficient of Variation of 5% (which is even larger than obtained experimentally in Muilwijk and Van den Akker (2019) ). Due to the dispersity in the (vertical) diameter, the right tail of the distribution is elongated, while the sharp peak at the modal chord length flattened. For equally sized bubbles with a uniformly distributed aspect ratio in E, a triangular chord length distribution was observed (dotted line) with a modal chord length at ≈ 0 . 8 × the maximum chord length. Almost no difference was found between a uniformly and arcsinusoidal (dash-dotted line) distributed aspect ratio.
For a normally distributed aspect ratio (dashed line), the bubble chord length histogram approaches a chord length distribution for a normally distributed bubble size distribution, while the bubbles are uniform in terms of volume equivalent diameter.
The deviation from a triangular chord length distribution, as given by the solid black line in Fig. A.17 , can be due to both (1) the non-uniformity of the bubble size; and (2) the non-constant shape in the wobbling regime. We therefore conclude that a chord length distribution cannot be transformed into a bubble size distribution when bubbles exhibit a non-constant shape.

A2. Aspect ratio as a function of the Eötvös number
The bubble aspect ratios are calculated, using Eq. (7) , from the mean of the chord length distributions and the volume equivalent bubble diameter, which is given by the correlation developed in our previous paper ( Muilwijk and Van den Akker, 2019 ). The bubble aspect ratio as a function of the Eötvös number ( Eo = ρ w gd 2 b,eq /σ ) is given in Fig. A.18 . It should be noted that, for U sl = 0 m/s and U sg > 3 . 5 cm/s, the bubble size can no longer be assumed as strictly uniform due to coalescence (and breakup) happening at the needle outlets.

Appendix B. Pairing rates
The pairing rates as a function of x, U sg and U sl are given in Fig. B.19 . The error bars (c) denote the spreading in lateral ( x ) direction, the 15 measurement locations as shown in Fig. 1 . The pairing rate was almost independent of U sg , but increased from 50% to 65% for U sl in the range 0-0.2 m/s due to (1) increased alignment of the bubble velocity direction and the probe; and (2) the higher inertia of the bubbles' added mass reducing the effect of drifting.

Supplementary material
Supplementary material associated with this article can be found, in the online version, at doi: 10.1016/j.ijmultiphaseflow.2020. 103498 .