Predictive model for the scale-out of small channel two-phase flow contactors

In this paper, double manifolds are theoretically investigated for the scale-out of two-phase incompressible flows in small channels. Statistical descriptors are proposed to characterise the maldistribution of the total flow rate and the ratio of the flow rates in the two-phase channels, based on the variances and covariance of the flow rates of the two fluids. A novel resistance network model is developed that relates the flowrates of the fluids in the two-phase channels to the hydraulic resistances of the manifold. The statistical descriptors and the resistance network model are then used to develop relationships between the maldistribution coefficients and the hydraulic resistances of the double manifold, the overall pressure drop and the pumping power requirements for different parallel channel numbers. Based on these, scaling laws are proposed that maintain a constant degree of maldistribution for a scale-up factor up to 10. Double manifolds designed using these scaling laws have a constant pressure drop as the number of channels increases, whilst the power requirements increase linearly. The power requirements are inversely proportional to the phase ratio maldistribution descriptor. Recommendations for the design of double manifolds for the scale-out of two-phase systems are proposed.


Introduction
Over the last few years there have been many studies on small-scale units aiming at the intensification of multiphase processes, including (bio)chemical syntheses [1,2], absorptions [3], extractions [4,5], and (nano)crystallisations [6]. Multiphase processes carried out in small channels are characterised by fast mixing, high heat and mass transfer rates, and narrow residence time distributions when compared to conventional equipment [7,8]. These advantages, in an industrial context, would translate to reduced costs, energy savings, diminished waste production, enhanced safety, and increased product quality [9][10][11]. The benefits of operating in small channels stem from the reduction in length scales which result in short diffusion distances, increased importance of interfacial forces over inertial, viscous, and gravitational ones, large interfacial area-to-volume ratio (beneficial for mass transfer), and large channel surface-to-volume ratio (beneficial for heat transfer). The transition of these novel and efficient devices from benchscale to the industrial application involves increasing the throughput, often by using many of the single small-scale units in parallel (scale-out or numbering-up). This contrasts with the traditional scale-up where the volume of the unit increases. However, the transition is hampered by the difficulty in reproducing the flow conditions of a single channel in many parallel ones. Well-characterised flow distribution manifolds are key to scale-out of small-scale devices via parallelisation.
There has been a significant amount of work on flow distributors for single-phase flows in small channels [12][13][14][15][16][17][18][19][20]. To bring together two fluids, double manifolds are used where both phases are distributed separately before they mix and flow together in the main contacting section (Fig. 1). The distribution of the individual phases takes place in consecutive manifolds similar to the ones used in single phase flows [12][13][14][15][16][17]. This type of two-phase flow distributor has inherent benefits, including simple arrangement, small footprint, modularity, and passive control.
The negative effects of flow maldistribution in manifolds have been well documented for single-phase flows [19][20][21]. The various channels may have different flowrates and residence times, which can affect selectivity or yield. Al-Rawashdeh et al. [22][23][24] studied the effects of fabrication tolerances to maldistribution in double manifolds for gas-liquid Taylor flows [22][23][24]. Schubert et al. [25] found that maldistribution of the phases in the flow manifold of a monolith reactor severely reduced the selectivity and conversion of glucose hydrogenation reactions. Flow maldistribution in single-phase flow distributors is usually quantified with statistical dispersion descriptors [12,15,16,19,[26][27][28][29][30]. However, to quantify maldistribution in twophase flows of immiscible fluids, the statistical descriptors must consider the bivariate nature of the system, where both the flowrate and the fraction of the phases can vary.
In this work, the design of double manifold distributors that bring together two fluid phases is studied for incompressible flows. The effects of the design variables on the flow maldistribution and the pressure drop and pumping power in the manifold are examined. In what follows, first a technique to quantify maldistribution in parallel channels for two-phase flows is proposed and appropriate statistical descriptors are suggested, based on multivariate statistical analysis. The double manifolds are then modelled using a resistance network approach which reduces the large number of design variables into two resistance ratios. The resistance model and the statistical analysis are subsequently used to quantify the maldistribution and pumping requirements for a large number of channels. Based on the relationships between the hydraulic resistances and the maldistribution descriptors scaling laws are suggested, using regression analysis. These can be used to calculate the dimensions of the double manifold as the throughput increases. Finally, the methodology for designing double manifolds for two-phase incompressible flows is summarised.

Quantification of two-phase flow maldistribution
The most common metrics of flow distribution used in single phase flow distributors are shown in Table 1, where Q j is the flow rate in the jth channel, Q is the mean flow rate, and N is the number of parallel channels. The coefficient of variation (CV), defined as the standard deviation of the flowrates relative to the mean, is the most common maldistribution metric. The mean value of the flow rate (Q) (also termed as nominal or design flow rate) plays an important role because it is the flow rate all channels should have if the flow is perfectly distributed.
Single-phase maldistribution metrics must be modified before they can be applied to two-phase cases. In addition to the overall flow rate, the ratio of the flow rates of the two phases is also important. Previous   After the main channels, a separation, collection or further processing units can be connected. works have only reported the CV either of one of the phase flow rates, of the total flow rate, of the bubble velocity in plug flows [24,34] or of both phase flow rates [23,[35][36][37]. These metrics are not sufficient in the multiphase cases because they do not consider the maldistribution of the ratio of the phases. The maldistribution in the flow of two immiscible and incompressible phases is studied here with a method derived from principal component analysis (PCA), a common statistical approach for multivariate data [38]. Following this analysis, all factors affecting maldistribution can be studied with the least number of variables. To completely characterise the dispersion of a two-variable system, in this case the maldistribution of two fluids in parallel channels, a 2 × 2 covariance matrix is required. The covariance matrix (Σ) for a bivariate distribution, shown in Eq. (1), is a symmetric matrix with the variances (σ 1 2 , σ 2 2 ) and the covariance terms (cov(1,2), cov(2,1)) in the two diagonals.
N is the number of the parallel channels of the manifold, Q j,i is the flow rate of the i-th phase in the j-th channel, and Q i is the average or nominal flow rate of the i-th phase in the channels. Σ can be calculated given the flow rate of each phase in every channel of the manifold (the Q j,i data). The correlation coefficient (ρ), defined in Eq. (4), is a dimensionless number calculated as the ratio of the covariance term and the product of the standard deviations of the two variables, σ 1 and σ 2 . It follows from the symmetry of Σ that three scalar variables are necessary to characterise the maldistribution completely. According to PCA the eigenvalues and eigenvectors of Σ are first calculated and then used to define new variables that best describe the dispersity of the data; these new variables are linked to the physical properties of the system. Eigenvalues satisfy the characteristic equation, det(Σ − λI) = 0; and eigenvectors satisfy the eigenvalue equation (Σ − λI)v = 0, where I is the identity matrix. The eigenvalues, the eigenvector slope (m), and the angle of correlation of the variables (θ) are calculated as follows: If |ρ| is small (i.e. < 0.05), Σ is approximated by a diagonal matrix (elements outside the main diagonal are insignificant) and the flow rates of each phase are uncorrelated. It follows from Eq. (5) that the eigenvalues are equal to the variances and in this case, two CVs, as defined by Eq. (7), can completely characterise the maldistribution in the flow distributor.
However, if |ρ| is not small (i.e. > 0.05), the flow rates of the two phases are correlated. In this case, λ 1 , λ 2 , and m, characterise maldistribution. In order to have a dimensionless measure of the maldistribution, the rotated coefficient of variation (RCV) is defined (Eq. (8)), which is analogue to CV in Eq. (7).
The eigenvector slope defined in Eq. (6) is related to the ratio of flow rates of both phases. The closer the value of m is to the value of the ratio, the better the fluids are distributed with respect to the flow rate ratio. Eq. (9) defines the phase ratio maldistribution descriptor (PRM).
PRM measures how different the eigenvector slope and the design phase ratio are. It tends to zero as the two values approach each other. There is an extreme case for correlated flow rates. When |ρ| is close to unity (i.e. > 0.95), the flow rate distribution of the two phases is highly correlated. In this situation λ 2 and RCV 2 are very small, and only RCV 1 , and PRM are necessary to describe the maldistribution. For the highly correlated case, the two descriptors can be interpreted as follows: • RCV 1 indicates the maldistribution of the total flow rate in each of the parallel channels, independently of the ratio of the phases.
• PRM indicates the differences in the phase ratio in the parallel channels from the required ratio. PRM is independent of the total flow rate in each channel and its maldistribution.
Three possibilities of correlation are illustrated in Fig. 2 in a map of flow rates, with Q 1 and Q 2 as coordinates. Fig. 2a shows an uncorrelated case where maldistribution is quantified with the standard deviations, non-dimensionalised as shown in Eq. (7). Fig. 2b shows a correlated case with both eigenvectors shown. Both eigenvalues and the eigenvector slope, non-dimensionalised with Eqs. (8) and (9), are used to measure maldistribution in this case. In Fig. 2c, a highly correlated case is presented and the points nearly form a line segment. Only one eigenvalue and the eigenvector slope are needed to quantify maldistribution in this case.
The analysis above describes a novel and rigorous way to measure maldistribution of two fluids in parallel channels. The three cases, dependant on the correlation coefficient, are summarised in Table 2.
In addition to the methodology above, Q 1 , Q 2 , λ 1 , λ 2 , and θ can be used to calculate confidence ellipses in flow rate maps. This two-dimensional confidence interval is centred in (Q 1 , Q 2 ) and is tilted at an angle θ. The major and minor axes lengths, a and b respectively, depend on the eigenvalues and the degree of confidence selected. The axes lengths are proportional to the square-root of the chi-squared critical value.
= a 2 5.991λ 1 (10) For a confidence of 95%, the chi-squared critical value for 2 degrees of freedom is 5.991 [39]. The confidence ellipse together with the flow rate map of the particular system can be used to define the maximum tolerable maldistribution of the system.
If the flowrates Q ji can be measured before the fluids are mixed together then the statistical analysis above can also be applied to miscible fluids. In some cases, it may be more useful to use mass flow rates rather than volumetric flow rates.

Resistance network model
Flow distributors, such as double manifolds, can be modelled using resistance network models, as reviewed by Oh et al. [40]. This type of modelling is preferred over computational fluid dynamics due to the simplicity and speed to setup and run parametric studies [12]. The resistance network model considers the analogy between fluid flow in hydraulic circuits and electrons flow in electric circuits. This modelling approach has been used extensively to design and optimise single-phase flow distributors [12][13][14]. The overall pressure drop and fabrication tolerances of the designed manifold were also considered by Amador [41].
The analogy is extended here to account for two-phase incompressible flows. The model considers volumes are additive, which means that it cannot be used when there are large volume changes due to compressibility or mass transfer.
A manifold for incompressible two-phase flows and its electricalcircuit analogue are shown in Fig. 3a and 3b respectively. Kirchhoff's current law (KCL) and Kirchhoff's voltage law (KVL) are used to solve the circuit in Fig. 3b. KCL is analogous to mass conservation, while KVL is analogous to energy conservation. Eqs. (12) and (13) are the algebraic representations of KCL and KVL respectively for flow systems.
where Q is the volumetric flow rate, Δp is the pressure drop, and R is the hydraulic resistance defined as the ratio of pressure drop over volumetric flow rate. The sum in Eq. (12) applies at every node of the circuit. The sum in Eq. (13) adds the pressure drops in a closed loop in the circuit. In a double manifold with N main channels, 2(N−1) loops can be independently defined. Between the j-th and (j+1)-th main channels, two loops are defined. These loops are highlighted in Fig. 4. Fig. 4a, shows a loop including channels labelled A1,j+1; B1,j+1; R,j +1; C,N-j+1; R,j; and B1,j, while Fig. 4b shows the other independent loop containing channels labelled A2,j+1; B2,j+1; R,j+1; C,N-j+1; R,j; and B2,j. For these two loops, Eqs. (14) and (15) are obtained by applying the KVL. In total, there are N−1 equations for each type of loop.
Applying KCL along the distribution channels (Sections A1 and A2) at each junction, Eqs. (16) and (17) The same procedure along the collection channels (Section C) yields Eq. (18).
Additionally, Eq. (19) is the result of applying KCL for the junctions where both manifolds meet before entering the main channels (Sections B1, B2, R). There are N equations like this, one per main channel.  7)). b) Correlated case where the eigenvalues and the eigenvector slope are used to measure the maldistribution (Eqs. (8) and (9)). c) Highly correlated distribution where the points nearly form a line segment; the length and slope of this line segment describe the maldistribution (Eqs. (8) and (9)).
Finally, Eqs. (20) and (21) express the mass balance at the respective manifold of each phase: Eqs. (16)- (18) are substituted in Eqs. (14) and (15). The resulting equations are divided by R R,j to yield Eqs. (22) and (23). As in the case for Eqs. (14) and (15) Eqs. (19)-(23) form a system of 3N linear equations for N main channels. 2N−2 equations are obtained with the KVL analogy (Eqs. (22) and (23)), there are N equations describing the mixing of the manifolds (Eq. (19)) and 2 equations for the overall mass balance (Eqs. (20) and (21)). The system can then be solved for 3N variables. These variables are the flow rates in the barrier channels for both phases (2N) and the main channels (N). Furthermore, with these results, and using Eqs. (16)-(18), the flow rates in the rest of the channels (distributing and collecting sections) can be calculated. Eqs. (24) and (25) represent the system to solve and the solution using matrix inversion, where R is the resistances matrix, Q is the flow rates vector and S, the solution vector, is the right-hand-side of Eqs. (19)- (23). The value of Q depends only on the values of resistance ratios in matrix R, which is dimensionless, and the inlet flow rates of both liquids, Q A1,1 and Q A2,1 in vector S. Eq. (25) can be solved for Q in a single step when all the hydraulic resistances are independent of the flow rate. When the resistance terms are not independent of the flow rates, Eq. (25) needs to be solved  iteratively. In this case, the solution for the independent resistances can be used as initial guess. The model is applicable if the volumes are conserved (i.e. the densities of both phases change by less than 5%).
In the single phase flow sections (distribution and barrier), the hydraulic resistances are independent of the flow rates if the flow is laminar and incompressible. For these conditions, the Hagen-Poiseuille equation (Eq. (26)) can be applied and the hydraulic resistance depends on the channel geometry.
The pressure drop in fittings such as reducers, expansions, elbows, and tees is not significant at low Reynolds numbers; however, these effects can also be included in the pressure drop correlation if necessary. For the pressure drop in the main channels with both phases present, two options can be considered: Either the multiphase flows have hydraulic resistances independent of the flow rate, e.g. homogeneous model [42,43]; or the resistances depend on the flow rates. In the latter case Eq. (25) is solved iteratively using a non-linear pressure drop model for the two-phase flows; Al-Rawashdeh et al. [22] showed this procedure.
In this work, the homogeneous model was selected for the twophase flow channels for two main reasons. With this model, Eq. (25) can be solved in a single step and a wide range of parameters can be studied with minimal computational costs. In addition, the effects of the design variables (the hydraulic resistance of each section) on the maldistribution and pumping requirements can be isolated from complex pressure drop correlations.
In vector S, the two inlet flow rate variables (Q A1,1 , Q A2,1 ) can be written in terms of the total flow rates in the main channels and the inlet flow rate ratio of the two phases, as follows: where Q T , Q T,1ch , and r, are the total inlet flow rate, the total flow rate in one main channel (assuming perfect distribution) and the ratio of the inlet flow rates, respectively. Both Q T,1ch and r are independent of N, which makes them more suitable for the scale-out studies of double manifolds than Q A1,1 and Q A2,1 . From the above: Both inlet flow rates, the only non-zero elements in S, have Q T (or Q T,1ch N) as a factor, therefore Eq. (25) can be non-dimensionalised dividing by Q T (or Q T,1ch N). Matrix R has as many variables as the number of channels in the manifold. In what follows, in each of the different manifold sections (distribution of the single phases, barrier, main, collection) the channels are assumed to have the same resistance. In this case the number of variables is reduced from 6N+3 (one for the number of main channels, two for the inlet flow rate of each phase and 6N for each resistor) to 7 variables (N, Q T , r, R A , R B , R R , R C ). Furthermore, if the collection section is not included (R C = 0), as is the case when the output of each main channel is collected separately, one less variable is needed. Additionally, it follows from Eqs. (22) and (23) that only the ratios of resistances are necessary, which eliminates one more variable. Thus, 5 input variables are necessary for using the model under these assumptions, namely N, Q T (or Q T,1ch ), r, R B /R R , and R A /R R (or any other combination of the resistance ratios). Of these input variables, R A and R B are the design variables of the double manifold; the rest of the variables depend on the throughput required and the requirements of the process in the single channel. The resistance network model and the flow maldistribution characterisation are studied with an in-house built Matlab code.

Results and discussion
The methodologies developed above are first used to demonstrate the effects of the hydraulic resistances on the flow maldistribution in manifolds with a certain number of parallel channels. In addition to flow distribution, the overall pressure drop and power requirements are then estimated. The effects of increasing the number of the parallel channels in the manifold on the flow distribution and on the pressure drop and power requirements are finally investigated. As an example of the resistance network model and the application of the flow maldistribution descriptors, a parametric study of a double manifold with 5 channels was carried out as shown in Table 3. The input variables are the flow rate ratio (r = [1,5]), and the resistance ratios (R A /R R = R B /R R = [0. 1,10]). For all cases, ρ > 0.05 and only for case 6 ρ < 0.95 thus requiring both RCVs and PRM to characterise the maldistribution ( Table 2).
The cases in Table 3 where the flow rate ratio is 1, show consistent values of ρ = 1 and PRM = 0. In these cases, the maldistribution is completely described by RCV 1 . The lowest RCV 1 value, thus lowest maldistribution, corresponds to case 3 where R A /R R = 0.1 and R B / R R = 10; this case is plotted in Fig. 5a to c. As can be seen in Fig. 5a and 5b the flow rate in each main channel is close to the equal distribution line. For both phases, channels 1 and 2 have more flow rate than the average and channels 3-5 have slightly lower flow rates. Fig. 5c shows the flow rates of the two phases in each main channel in a flow rate map, together with the point of perfect distribution. In the same graph, the 95% confidence ellipse is plotted. This ellipse has its centre in the point of perfect distribution and is rotated by angle, θ = 45°, calculated from Eq. (6). The major and minor axes lengths are calculated from Eqs. (10) and (11). PRM in this case is zero which indicates the ratio in all Table 3 Maldistribution descriptors for r: 1 and 5, R A /R R : 0.1 and 10, R B /R R : 0.1 and 10. The double manifold modelled has 5 channels and total flow rate of 6 × 10 −6 m 3 s −1 . the channels is equal to the inlet ratio and thus the eigenvector slope (m = tan (θ)) is equal to 1/r. Since the minor axis length is zero (from ρ = 1 and Eqs. (5) and (11)), the confidence ellipse collapses into a line segment. Case 5 in Table 3, with r = 5, and both R A and R B smaller than R R presents channeling, where one of the phases flows into a channel of the opposite phase. In this case, channeling occurs because phase 1, with a larger flow rate (r = 5) has a much larger pressure than phase two and flows into a barrier channel of phase 2. In the model this is indicated by a negative flow rate (opposite to the positive-defined direction) in the channel where channeling has occurred. The rest of the cases with r = 5, have ρ > 0.05 and thus both RCVs and PRM are used to describe maldistribution. The largest maldistribution descriptors occur when R A /R R = 10 and R B /R R = 0.1 (case 6). Case 7 with R A /R R = 0.1 and R B /R R = 10 has the lowest maldistribution with the smallest RCV 1 (ρ > 0.95 so RCV 2 is not considered). Case 7 is plotted in Fig. 5d-f. The distribution profiles are similar to the r = 1 case although the ranges for each phase are different because of the different phase ratio. The flow rates in Fig. 5f fall in a straight line (as in Fig. 5c) for r = 5 though the slope of the eigenvector (m = 0.113) is different to the inverse of r (1/ r = 0.2). This difference represents the maldistribution of the phase ratio in each main channel, as measured by the PRM descriptor.
Case 8, with high both R A and R B compared to R R has a larger RCV 1 but smaller PRM than case 7. These conflicting results show that further considerations should be taken into account when characterising maldistribution in double manifolds. For both r = 1 and 5, RCV 1 is the smallest when R A /R R = 0.1 and R B /R R = 10. These ratios compare the resistances of Sections A and B with R. However, they also indirectly compare the resistances between Sections A and B: For both flow rate ratios analysed, the lowest RCV 1 values are achieved when R A /R B is the smallest. This suggests R A /R B describes RCV 1 better than the other combinations of resistance ratios. Fig. 6 presents logarithmic contour plots of RCV 1 and PRM in R A /R B − R B /R R planes for N = 5 and r = 5. For the whole range of resistance ratios plotted, ρ is larger than 0.95 and is increasingly closer to unity as R A /R B decreases and R B /R R increases in the range plotted.  Table 3 (N = 5, Q T = 6 × 10 −6 m 3 s −1 ). a, b, d, e) Flow rates of each phase in the channels; the perfect distribution line is calculated as Q Ai,1 /N. c, f) Flow rate maps. The circle indicates the perfect distribution case, calculated as (Q A1,1 /N, Q A2,1 /N). The 95% confidence ellipses are also shown; the minor axis length, b, is zero in both cases and the ellipses collapse into line segments. Fig. 6a show RCV 1 decreases steeply as R B increases with respect to R A . This is shown by both the quasi-vertical contour lines and the arrows. By contrast, the ratio between R B and R R has very limited influence on the behaviour of RCV 1 . Similar results were found for r = 1. Furthermore, the model predicts channeling at log 10 (R B / R R ) < 0.5 for the range of R A /R B plotted, which means that the flow distribution is severely disrupted at R B < 3R R . Fig. 6b presents the log 10 PRM in the same plane as Fig. 6a. The results show that R B /R R has a stronger influence over PRM than R A /R B . However, the influence is less as R A /R B increases, where the arrows and the curvature of the contour lines show that both resistance ratios affect PRM. If the arrows of greatest decent in both Fig. 6a and 6b are superimposed, they are almost perpendicular to each other. This means that the different resistance ratios affect the maldistribution descriptors practically independently. RCV 1 , which is analogous to total flowrate and residence time distribution widening, depends mostly on R A /R B . PRM, which is related to flow rate ratio distribution, is dependent on R B /R R primarily.

Results in
The current model is compared against the results of the two-phase resistance network developed by Al-Rawashdeh et al. [22,23]. Theauthors suggested channeling occurs when the resistivity of the barrier sections is different for each phase and the flow rates of the phases are equal [23]. This contrasts with the findings of the current model which predicts channeling only for flow rate ratios different to unity and when R B /R R is below a threshold value. Regarding flow maldistribution, Al-Rawashdeh's model used a single maldistribution descriptor which decreases as R A /R R decreases and R B /R R increases. In contrast, results of the current work show maldistribution is completely characterised with at least two descriptors and these descriptors indicate that the residence time and the flow rate ratio distributions become narrower as R A /R B is reduced and R B /R R is increased.

Effect of hydraulic resistances on the pressure drop and power requirements of the double manifold
The pressure drop in the double manifold and the associated power requirements are discussed here for different channel resistance ratios. The equivalent resistance (R eq ) is the resistivity of a single resistor that has the same resistance as the whole network. Since a double manifold distributes two phases, this definition is modified and two equivalent resistances are defined by Eq. (32), one for each phase.
where Δp i is the pressure drop of the i-th fluid and is numerically equal to the head required by the i-th pump and Q i is the total flow rate of the i-th phase. To compute R eq,i for each phase, the following are considered. First, for each phase, only the sections where that phase flows affect its R eq ; this is true as long as there is no channeling. Second, and as in the previous discussion, no collection section is included in the double manifold. The equivalent resistance for each phase in double manifolds with 2 and 3 main channels can be calculated by Eqs. (33) and (34) respectively. These equations show that even for a small number of channels, the calculation of R eq,i is complex.
However, when R Ai ≪ R Bi (necessary for reducing maldistribution) the denominator in the second term of the right-hand-side depends only on R Bi + R R . An approximation can then be made to find the value of R eq,i for any number of channels, given by: Eq. (35) shows that the dimensionless equivalent resistance increases linearly with R B /R R . Fig. 7 shows plots of R eq,i /R R over R A /R B for different number of channels, N, at constant R B /R R both from the exact solutions (Eqs. (33) and (34)) and the approximation (Eq. (35)). As can be seen, the approximation improves as R A /R B decreases. This figure also shows that R eq,i /R R decreases as N increases. The overall pressure drop for each phase can then be calculated  from Eqs. (36) and (37) and the total pumping power required from Eq. (38). These equations are derived using Eqs. (29), (30), (32), and (35).
where P T is the total pumping power required, and η i is the efficiency of the pumping system for the i-th fluid. The right-hand-side of these equations depends on the design of the double manifold. Both pressure drops and the power requirement increase linearly with R B /R R , while PRM (Fig. 6b) decreases with R B /R R . The pressure drop estimates, from Eqs. (36) and (37), can be used to calculate the changes in density if one of the fluids is a gas, and assess the validity of the model in Section 3. A contour plot of P T /(Q T 2 R R ) in a log-log plane of R A /R B − R B /R R is presented in Fig. 8, for a 5-channel manifold operating at a flow rate ratio of 5. The efficiency of both pumps is set to be 80%. Arrows of steepest descent for the dimensionless power requirement are also included. Fig. 8 is plotted in the same range of resistance ratios as Fig. 6 for comparison. Figs The purpose of double manifolds is to increase the throughputs compared to single channels. The effect of the number of main channels on flow maldistribution and the changes needed in the channel resistances to reduce or keep constant the maldistribution during scaleout are considered here. In this section Q T,1ch (defined in Eq. (27)) is used instead of Q T because it is independent of N. Fig. 9 shows flow rate maps for three cases with 5, 50 and 500 channels. The same Q T,1ch , flow rate ratio, and resistance ratios are used for all cases. As can be seen, the distributions are highly correlated, with ρ close to unity. As the number of channels increases, maldistribution increases with respect to RCV 1 but not with respect to PRM. The increase in RCV 1 with N means that the residence time distribution becomes wider. RCV 1 increases by a factor of 84 and of 732 as the number of main channels increases from 5 to 50 and 500, respectively. In the 500 channels distributor, there is significant maldistribution with several channels having flow rates close to zero. PRM, which compares the flow rate ratios in the channels with the inlet flow rate ratio, shows the opposite trend to RCV 1 as N increases. It decreases by 18% and by 65% when N increases from 5 to 50 channels and then from 50 to 500 channels, respectively. In all cases, PRM is very small and the flow rate ratios are very similar in all channels despite the large maldistribution measured by RCV 1 .
The effect of resistance ratios on flow maldistribution as the number of channels increases was also analysed. Fig. 10 presents the logarithmic plots of both maldistribution descriptors as function of N for flow rate ratio equal to 5. These plots present four cases for R A /R B equal to 10 −3 or 10 −2 and R B /R R equal to 10 1 or 10 2 . The results show that the trends observed for 5 channels (Fig. 6) also apply when N is further increased. RCV 1 depends predominantly on R A /R B while the effect of R B /R R is not significant. RCV 1 , and the flow maldistribution it describes, are reduced with reducing R A /R B values. In addition, RCV 1 increases with the number of channels; the rate of increase is high at low N values and is reduced at high N values. This trend can be explained by considering the flow rates in the main channels for N = 500. In this case, the maldistribution is very high and some channels have flow rates very close to zero (Fig. 9c); as a result the increase of RCV 1 , which is statistical in nature, is not as pronounced for large channel numbers. The results for PRM in Fig. 10b are also in agreement with those of Fig. 6b. PRM depends mainly on R B /R R and less on R A /R B . PRM, and thus the distribution of flow rate ratios, are reduced when R B /R R is increased. PRM is constant at low values of N and reduces to a lower constant value at large values of N. This transition is also affected by the channels with zero flow rates at high N values, similar to RCV 1 described above. The value of N at which the PRM transition occurs depends on R A /R B .
To facilitate the use of the above observations in the design of manifolds with many parallel channels, scaling laws are suggested below, based on regression analysis. The equations will allow the design of manifolds with channel resistances that result in a given flow maldistribution or vice versa.
The variables investigated in Figs. 6 and 10 and their ranges are presented in Table 4. The first case examined is for inlet flow rate ratio equal to 1. The results for the whole range analysed have ρ larger than 0.95 and PRM equal to zero. Therefore, RCV 1 fully describes the maldistribution ( Table 2). The regression coefficients are given in Eq. (39).
R B /R R does not appear in Eq. (39) because its p-value in the regression is larger than 0.05 (p-value = 0.62) and it does not affect RCV 1 significantly. The exponents for N and R A /R B have 95% confidence intervals of ± 5.8% and ± 4.6%, respectively. Al-Rawashdeh et al. [15] Fig. 8. Contour plot of log 10 P T /(Q T 2 R R ) in a log 10 -log 10 plane of R A /R B vs. R B / R R . The double manifold modelled is symmetrical, has 5 main channels, Q A1,1 = 5 × 10 −6 m 3 s 1 , Q A2,1 = 1 × 10 −6 m 3 s −1 . The pumping efficiencies are 80%.The arrows show the direction of greatest decline for the dimensionless pumping power.
presented a similar correlation even though they measured maldistribution differently. Their correlation agrees with the results presented here in that the exponent for N is positive and larger than one.
A second regression is performed for the same variables and the inlet flow rate ratio varying between 2 and 8 ( Table 4). For the range of variables analysed, ρ is larger than 0.95. Therefore, the two descriptors, RCV 1 and PRM, contain all the statistical information needed to describe maldistribution ( Table 2). Eqs. (40) and (41) Neither R B /R R nor the inlet flow rate ratio appear in Eq. (40) because their p-values are larger than 0.05 (0.48 and 0.50, respectively). The 95% confidence intervals of the N and R A /R B exponents in Eq. (40) are ± 2.8% and ± 2.2%, respectively. Since RCV 1 is independent of r, the results for one value of r apply for other flow rate ratios within the range of r considered in the regression (2 ≤ r ≤ 8). In the PRM regression all the variables have p-values below 0.05. N, R A /R B , R B /R R , and r have 95% confidence intervals of ± 17%, ± 15%, ± 0.71%, and ± 5.8%, respectively.
With respect to RCV 1 , the regression for the cases where r > 1 show very similar results to the r = 1 case. The exponents are almost identical and the additional variable, the inlet flow rate ratio, has a very minor influence. When the flow rate ratio is larger than unity, PRM is different to zero and important for describing maldistribution. The exponent of N in Eq. (41) is negative and close to zero, which means that the maldistribution does not become worse as the number of channels increases. The negative exponents in the resistance ratios factors imply that PRM increases as these ratios decrease. The R B /R R exponent is 21 times larger than the R A /R B exponent; this confirms that PRM is dominated by R B /R R . Additionally, the exponent for R B /R R in Eq. (41) is almost −1 while, as was discussed, the pumping requirements increase linearly with R B /R R (Eqs. (35)- (38)). This indicates that maldistribution and pumping requirements are inversely associated to each other. The flow rate ratio has the largest exponent in Eq. (41), which indicates effective flow distribution is harder to achieve the larger the flow rate ratio is.
The equations provide a powerful tool to assess the changes needed in the manifold design to maintain good flow distribution when the number of channels is increased. Eqs. (39)- (41) can be used to obtain scaling laws for the resistance ratios. This is done by keeping constant the maldistribution descriptor for two different numbers of channels, N1 and N2. Eqs. (42) and (43) are the scaling laws for double manifolds. They calculate the changes in the resistance ratios required to maintain the flow distribution quality as the number of channels changes, for the same main channel resistance (R R ) and inlet flow rate ratio (r). The values of the exponents are given in Table 5 for both cases of r equal to or larger than one.
The ratio N1/N2 in Eqs. (42) and (43) is equal to the throughput ratio of manifolds with N1 and N2 main channels. R B does not change with the number of channels for r = 1. The confidence interval for Y when r > 1 is large, but since the exponent is close to zero, it only has a minor effect on the R B(N2) calculation and design. For example, when 100 times more channels (and throughput) are considered, the uncertainty on the design of R B(N2) is around 11%.
The scaling laws proposed above are applied to the scale-out of the double manifolds modelled in Fig. 9. The double manifold with 5 main channels has a very low degree of maldistribution. However, when the number of channels increased by one or two orders of magnitude, the maldistribution worsened. The scale-out in Fig. 9 was performed by keeping the resistances constant as the number of channels increased. With the scaling laws proposed in Eqs. (42) and (43), the resistances are changed as a function of the number of channels and the improved results are shown in Fig. 11. The resistance ratios for the N = 5 case are the same as in Fig. 9a, R A /R B = 10 −3 , R B /R R = 10 2 . The cases with N = 50 and 500, have the following resistance ratios: R A / R B = 2.08 × 10 −5 , and 4.33 × 10 −7 ; R B /R R = 97.5, and 95.0, respectively. In all cases, the inlet flow rate ratio is 5 and Q T,1ch is 1.2 × 10 −6 m 3 s −1 . For the three cases in Fig. 11, ρ is larger than 0.95 and the maldistribution is described by both RCV 1 and PRM (Table 2).
By comparing Figs. 9 and 11 the improvements achieved with the scaling laws are evident. In the case of N = 50, when the resistances are not scaled, RCV 1 increases 84 times from N = 5; when the resistances are scaled, RCV 1 increases by only 12%. This is not zero because the Eqs. (42) and (43) do not fit perfectly the results of the double manifold model. Similarly, when N increases from 5 to 500 without the adjustment in resistances, the flow is severely maldistributed. In contrast, when the resistances are adjusted following Eqs. (42) and (43), RCV 1 for N = 500 is only 4.4 times RCV 1 for N = 5 (instead of 732 times). The results are very similar for the PRM indicator. These findings have implications on the modularity of double manifolds. Modular double manifolds designed for large throughputs will be suitable for lower throughputs, while the opposite is not true. It should be noted that the improvements in the flow distribution are accompanied by changes in the resistances, which affect the pressure drop and power requirements; this is explored in the next section.

Resistance network model results for pressure drop and power requirements
The effects of increasing the number of channels on the pressure drop in the manifold and the pumping power requirements are investigated, using the scaling laws presented in Eqs. (42) and (43). The change in pressure drop for either fluid as the number of channels increases is given by Eq. (44). The change in total pumping power with increasing number of channels is given by Eq. (45) The X and Y exponents are given in Table 5. Both the change in pressure drop and total pumping power depend on the values of N1 and N2 and not only on their ratio. The base case double manifold modelled above in Figs. 9 and 11 (N1 = 5, r = 5, R A /R B = 10 −3 , and R B / R R = 10 2 ) is used to illustrate how the pressure drop and the pumping requirements change as the number of main channels increases using the scaling laws given by Eqs. (42) and (43). As can be seen from Fig. 12a, the pressure drop of either phase remains almost constant with increasing number of main channels. For example for 500 main channels, the pressure drop is 0.95 ± 10% times the pressure drop for 5 main channels. This shows that a double manifold scaled by keeping the The 95% confidence intervals for the coefficient and exponent in Eq. (47) are 3.8% and 2.3% respectively. Eq. (47) is not a generalised scaling law because it is based on a specific case; however, it shows that despite the complexity of Eqs. (45) and (46), the pumping power scales almost linearly with the number of channels.
The almost linear increase of the pumping power with the number of channels has cost implications for the scale-out and indicates nearabsence of economies of scale. Annual pumping costs depend mostly on the energy consumption (> 85%) [44]; also, pumps have large scaling exponents (∼0.9) on capital cost [45]. While the linearity of cost curves for the capital investment of microreactor technology had been proposed in the literature [46], the linearity of the pumping costs curve has not been shown before.

Design methodology for two-phase double manifolds
Based on the previous findings, a general methodology is proposed below for the design of double manifolds for incompressible two-phase flows or the mixing of miscible incompressible fluids. A general design methodology, including pumps and control systems selection, has been given by Zhang et al. [7].  11. Flow rate maps of scaled-out double manifolds using the scaling laws in Eqs. (42) and (43) for Q T,1ch = 1.2 × 10 −6 m 3 s −1 , r = 5. The axes do not start in the origin. The resistance ratios are R A /R B = 10 −3 R B /R R = 10 2 for N = 5, R A /R B = 2.08 × 10 −5 R B /R R = 97.5 for N = 50, and R A / R B = 4.33 × 10 −7 R B /R R = 95.0 for N = 500.  (46), respectively. The 95% confidence intervals lines in red dashed lines are plotted using the confidence intervals in Table 5. The base case has N1 = 5, r = 5, R A /R B = 10 −3 , and R B /R R = 100. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) 1) Single channel investigations. These studies will reveal key parameters, such as optimal geometry (channel diameter and length) and flow rates for each phase (Q T,1ch and r) to achieve the required conversion and yield (in the case of reactors), and a function of the pressure drop with respect to the flow rate of each phase. It will be difficult to avoid maldistribution in the scale-out configurations. For this reason, it is important that the single channel studies include a maldistribution sensitivity analysis which will show how maldistribution affects key properties, such as conversion, and what level of maldistribution parameters (i.e. RCV 1 , RCV 2 and PRM) are acceptable. 2) Design the double manifold for a required throughput. Eqs. (36)- (38). At this stage, the trade-off between power requirements and maldistribution parameters should be considered; if the power requirements are high, different, less stringent, maldistribution coefficients may need to be used in step 2.

Conclusions
In this paper, the design of double manifolds (Fig. 1) for the scaleout of incompressible two-phase processes in small channels was studied. A statistical method, derived from principal component analysis, was developed to quantify maldistribution of both phases in two-phase flow distributors. It showed that two descriptors, RCV 1 and PRM, are relevant in most cases. The former is a measure of the total flow rate and residence time distribution in the main scale-out channels, while the latter is a measure of the distribution of the flow rate ratio in the main channels. The flow distribution was modelled using a resistance network model for incompressible flows. In the cases studied, the hydraulic resistances were considered to be independent of the flow rates and only depended on the sizes of the channels. The methodology can also be used when the resistances depend on the flow rates, but in this case the network model has to be solved iteratively.
Using the statistical and the resistance network methods, relationships between the design variables and maldistribution descriptors were found (Eqs. (39)-(41)). It was shown that RCV 1 decreased as the ratio R A /R B decreased. The descriptor PRM decreased as R B /R R increased. Scaling laws of the hydraulic resistances for the distribution and barrier sections were also obtained which can be used to scale-out double manifolds, while maintaining the same degree of maldistribution. When the scaling laws were applied, the total pressure drop of a double manifold did not change significantly with increasing number of channels. The pumping power requirements, however, increased linearly with increasing number of channels and throughput. The trade-off between flow maldistribution and power requirements and the implications on economies-of-scale were discussed. Based on these results, criteria for the design of double manifolds were presented.