Silva et al. (2008) proposed a two-dimensional test case to evaluate the performance of the direct quadrature method of moments (DQMOM) (Marchisio and Fox 2005) implemented in CFD codes, ANSYS CFX, and OpenFoam. They defined a fictitious problem of water in oil emulsion laminar flow through a backward-facing step. The same case was adopted for this study.

Figure 4.5 shows the backward-facing step geometry and dimensions. Analysis of the results is shown on the two dashed lines shown in Fig. 4.5, referred to as horizontal and vertical lines. Water in oil emulsion (i.e., dispersed water droplets in

Fig. 4.4 Particle size distribution at t = 450 s computed using 12 moments (This figure was originally published in Chemical Engineering Science 102, 2013 and has been reused with permission)

Fig. 4.5 Backward-facing step geometry (This figure was originally published in Chemical Engineering Science 102, 2013 and has been reused with permission)

continuous oil flow) enters the domain with a predefined initial droplet size distribution. The droplet size distribution evolves inside the domain as a result of convection and coalescence of the droplets.

Laminar flow in backward-facing step possesses unique features including presence of circulation zones with strong gradients in the laminar regime that makes it suitable for our purpose without bringing other uncertainties such as turbulence modeling into the picture. It is well known that there is a primary circulation zone behind the step in which its length (reattachment length) increases by increasing the Reynolds number. As shown in Fig. 4.6, at higher Reynolds numbers (~300), a secondary circulation zone will form at the channel upper wall due to the adverse pressure gradient coming from sudden expansion at the edge of the step (Biswas et al. 2004).

Existence of circulation zones means longer residence time of the disperse phase in those regions; therefore, a higher degree of coalescence is expected in these zones.

Fig. 4.6 Circulation zones in a laminar backward-facing step flow at different Reynolds number (This figure was originally published in Chemical Engineering Science 102, 2013 and has been reused with permission)

Table 4.1 Inlet condition and emulsion properties

Dimensionless initial size distribution of water droplets was defined by Eq. (4.33), and it is assumed that coalescence is the only process that changes the droplet size distribution according to the sum aggregation kernel (Eq. 4.34). A set of governing transport equations for the moments of the droplet size distribution was derived based on Eq. (4.23). In this equation, it is assumed that all of the droplets share the same velocity regardless of their size. Mean droplet diameter, d_{10}, was defined as the ratio of the first-order moment (ц_{1}) to the zeroth-order moment (p_{0}).

The quadrature method of moments, QMOM (Marchisio et al. 2003a, b), which is a built-in population balance solution method in ANSYS Fluent code and FCMOM, was used to generate results using two different numerical techniques. This could consider a verification of one numerical technique versus another. Baseline simulations were performed using 4, 6, and 8 moments based on QMOM to study the sensitivity of the solution to the number of moments and to compare the results obtained by FCMOM and QMOM.

Moreover, the effect of the mesh resolution was studied using both the FCMOM and QMOM methods having 4 moments and using three different uniform hexahedral meshes, referred to as coarse, medium, and fine meshes with 8400, 33,600 and 134,400 cells, respectively. Table 4.1 shows the inlet conditions and flow properties used in the simulations.

Spatial discretization of the governing equations was performed using a first- order upwind scheme along with a first-order and implicit temporal discretization. Pressure-velocity coupling was accomplished using the Phase Coupled SIMPLE algorithm. The time step of 10^{—3} s ensures convergence within 50 iterations per time step. In addition, the higher-order spatial discretization schemes were employed for discretization of the transport equations of both QMOM (Mazzei et al. 2012) and FCMOM and both resulted in moment corruption. Moment corruption means that a certain set of moments are non-realizable; hence, they no longer represent a valid distribution function. For FCMOM, the moments are calculated in the [—1, 1] interval, which results in negative odd and positive even moments. Therefore, moment corruption was observed when moments changed their sign (i.e., an odd moment became positive or an even moment became negative) (Abbasi and Arastoopour 2013).

The same behavior was reported by many other authors. Vikas et al. (2011) and Desjardins et al. (2008) showed that realizability is guaranteed only with the first- order finite volume schemes in spite of the numerical diffusion associated with these schemes. Mazzei et al. (2012) reported that the second-order upwind scheme for a volume-density-function-based QMOM significantly affects the stability of the solution and eventually causes moment corruption. Therefore, the first-order upwind scheme was used for the method of moments techniques. Figures 4.7 and 4.8 show the mean droplet size along the vertical line for different mesh resolutions, using 4 moments in QMOM and FCMOM, respectively. Both graphs show similar

Fig. 4.7 Mean droplet size along the vertical line for different mesh resolutions with QMOM and four moments at t = 1 s (This figure was originally published in Chemical Engineering Science 102, 2013 and has been reused with permission)

Fig. 4.8 Mean droplet size along the vertical line for different mesh resolutions with FCMOM and four moments at t = 1 s (This figure was originally published in Chemical Engineering Science 102, 2013 and has been reused with permission)

behaviors, having larger droplet size in the circulation region because of the longer residence time in this region; however, FCMOM appears to be slightly less sensitive to mesh resolution. Based on the grid independence study, the mesh with medium resolution was chosen for further simulations.

Figure 4.9 shows the mean droplet size along the vertical line for different number of moments calculated by QMOM. In the absence of exact solution and knowing that the error of the QMOM method decreases by increasing the number of moments, QMOM-8 was used as the basis for comparison of all results. The average error for d_{10} calculated by QMOM-4 and QMOM-6 with respect to QMOM-8 was 0.3 % and 0.03 %, respectively.

Figure 4.10 shows the mean droplet size along the vertical line calculated by FCMOM using 4 and 6 moments. Although d_{10} shows the same behavior along the vertical line, it has larger values of average error that reveal greater sensitivity of FCMOM to the number of moments compared to QMOM. The average error for d_{10 }calculated by FCMOM-4 and FCMOM-6 with respect to QMOM-8 was 1 % and 0.9%, respectively. Comparison between the calculated mean droplet sizes using FCMOM-6 and QMOM-8 along the horizontal line is shown in Fig. 4.11.

The mean droplet diameter contours calculated by FCMOM-6 at t = 1 s are presented in Fig. 4.12, which are essentially very close to those of QMOM. The disperse phase enters the domain with a 12.5 ^m mean droplet diameter. The values of d_{10} continuously evolve because of the coalescence of droplets. As expected, larger droplets form in the lower and upper circulation zones mainly because of the longer residence time of the dispersed phase in those regions. Therefore, the d_{10}

Fig. 4.9 Mean droplet size along the vertical line for different number of moments with QMOM att = 1 s (This figure was originally published in Chemical Engineering Science 102, 2013 and has been reused with permission)

Fig. 4.10 Mean droplet size along the vertical line for different number of moments with FCMOM compared to QMOM-8 at t = 1 s (This figure was originally published in Chemical Engineering Science 102, 2013 and has been reused with permission)

profile directly depends on the path-line and velocity field calculated by the code. The velocity field calculations are the same for both the QMOM and FCMOM method; however, the d_{10} contours are slightly different because of the possible inaccuracy in the calculation of moments.

Fig. 4.11 Mean droplet size along the horizontal line for FCMOM-6 compared to QMOM-8 at t = 1s (This figure was originally published in Chemical Engineering Science 102, 2013 and has been reused with permission)

Fig. 4.12 Mean droplet diameter d_{10} contours calculated by FCMOM-6 (contour levels are in cm) att = 1 s (This figure was originally published in Chemical Engineering Science 102, 2013 and has been reused with permission)

A difference between FCMOM and QMOM techniques is that FCMOM converges faster in comparison to QMOM. Figure 4.13 compares the CPU time requirement as a function of the number of grid cells for FCMOM and QMOM in an inhomogeneous aggregation problem using 4 moments. Simulations were performed using a 3.33 GHz Intel Xenon CPU on 6 parallel cores for 50 time steps with the time step size of 10~^{3} s. The comparison of the results clearly shows the better performance of FCMOM over QMOM, especially for a higher number of grid cells. For the finest grid with 134,000 grid cells, the CPU time of FCMOM is about half of that for QMOM (Abbasi and Arastoopour 2013).

Strumendo and Arastoopour (2008) have shown that increasing the number of moments increases the accuracy of FCMOM in homogeneous cases. However, Abbasi and Arastoopour (2013) were not able to reach a converged solution by increasing the number of moments to eight. In fact, tracking higher-order moments has the same effect as using a second-order discretization scheme. In other words, higher-order moments are less stable and will form non-realizable moment sets much faster than lower-order moments. In the earlier stages of the simulation, when

Fig. 4.13 Comparison of the CPU time as a function of the number of grid cells for QMOM and FCMOM (This figure was originally published in Chemical Engineering Science 102, 2013 and has been reused with permission)

the disperse phase reached the circulation zone where sharp gradients exist, moment corruption started and the error propagated into the domain that eventually leads to instabilities and divergence of the solution. To make sure that this phenomenon is not due to other factors, Abbasi and Arastoopour (2013) performed various simulations utilizing four times finer mesh, three orders of magnitude smaller time steps, and significantly smaller convergence criteria. In all cases, divergence occurred at the same flow time. Mazzei et al. (2012) also reported that quadrature nodes of the distribution become negative and corrupt the moment set when the number of moments is increased to eight or, as mentioned before, when second-order schemes were used for spatial discretization of moment transport equations, which confirms our observations. McGraw (2006), Wright (2007), and Kah et al. (2012) described this issue as an inherent problem of finite volume convective schemes as each moment is convected correctly but with independent equations. Therefore, the overall transport algorithm does not preserve the moment space and corrupts the relations among the moments.

A remedy proposed in the literature (McGraw 2006; Mazzei et al. 2012) is using moment-correction algorithms to force the moment set to always be valid (i.e., nonnegative Hankel-Hadamard determinants). The algorithm makes the least possible change to the corrupted (invalid) moments in a way that the new moment set meets the validity criteria (i.e., Hankel-Hadamard determinants become nonnegative). An alternative approach can be the use of high-order finite volume schemes similar to those proposed by Vikas et al. (2011) and Kah et al. (2012) that naturally preserve the moment space with very limited numerical diffusion. Such methods guarantee nonnegative distribution functions through an advective transport. However, before moving toward using higher-order finite-volume methods, it should be noted again that there is an inherent problem associated with these techniques such as FCMOM which approximates the distribution function with a finite number of terms in series expansion of orthogonal functions (e.g., Legendre polynomials). In this case, even with a valid set of moments, particle size distribution (PSD) gains negative values at some points on the internal coordinate, known as the Gibbs phenomenon (Gottlieb and Shu 1997; Arias-Zugasti 2012). Therefore, it seems that the first step toward mitigation is to eliminate, or at least minimize, the effect of the Gibbs phenomenon. There is a rich history of using Gegenbauer (i.e., Legendre, Chebyshev, etc.) polynomials in other engineering fields such as image processing in addition to strategies to overcome the Gibbs phenomenon (Silver et al. 1996; Shu et al. 2010; Yap et al. 2001). For instance, the use of well-developed noise-filtering techniques such as the kernel polynomial method (KPM) (Silver et al. 1996) seems promising in our application of interest.